diff --git a/qiskit/providers/aer/openpulse/qobj/operators.py b/qiskit/providers/aer/openpulse/qobj/operators.py
index 5546e7f795..b180128c45 100644
--- a/qiskit/providers/aer/openpulse/qobj/operators.py
+++ b/qiskit/providers/aer/openpulse/qobj/operators.py
@@ -23,8 +23,8 @@ def gen_oper(opname, index, h_osc, h_qub, states=None):
     ----------
     opname (str): Name of the operator to be returned.
     index (int): Index of operator.
-    h_osc (list): Dimension of oscillator subspace
-    h_qub (list): Dimension of qubit subspace
+    h_osc (dict): Dimension of oscillator subspace
+    h_qub (dict): Dimension of qubit subspace
     states (tuple): State indices of projection operator.
 
     Returns
@@ -35,32 +35,34 @@ def gen_oper(opname, index, h_osc, h_qub, states=None):
     # get number of levels in Hilbert space
     if opname in ['X', 'Y', 'Z', 'Sp', 'Sm', 'I', 'O', 'P']:
         is_qubit = True
-        dim = h_qub[index]
+        dim = h_qub.get(index, 2)
     else:
         is_qubit = False
-        dim = h_osc[index]
+        dim = h_osc.get(index, 5)
 
     if opname == 'P':
         opr_tmp = op.get_oper(opname, dim, states)
     else:
         opr_tmp = op.get_oper(opname, dim)
 
-    # qubit_0 * … * qubit_n * osc_0 * … * osc_n
+    # reverse sort by index
+    rev_h_osc = sorted(h_osc.items(), key=lambda x: x[0])[::-1]
+    rev_h_qub = sorted(h_qub.items(), key=lambda x: x[0])[::-1]
+
+    # osc_n * … * osc_0 * qubit_n * … * qubit_0
     opers = []
-    for ii, dd in enumerate(h_qub):
-        if ii == index and is_qubit:
+    for ii, dd in rev_h_osc:
+        if ii == index and not is_qubit:
             opers.append(opr_tmp)
         else:
             opers.append(op.qeye(dd))
-    for ii, dd in enumerate(h_osc):
-        if ii == index and not is_qubit:
+    for ii, dd in rev_h_qub:
+        if ii == index and is_qubit:
             opers.append(opr_tmp)
         else:
             opers.append(op.qeye(dd))
 
-    # return in reverse order
-    # osc_n * … * osc_0 * qubit_n * … * qubit_0
-    return op.tensor(opers[::-1])
+    return op.tensor(opers)
 
 
 def qubit_occ_oper(target_qubit, h_osc, h_qub, level=0):
@@ -71,28 +73,29 @@ def qubit_occ_oper(target_qubit, h_osc, h_qub, level=0):
     Parameters
     ----------
     target_qubit (int): Qubit for which operator is built.
-    h_osc (list): Dimension of oscillator subspace
-    h_qub (list): Dimension of qubit subspace
+    h_osc (dict): Dict of number of levels in each oscillator.
+    h_qub (dict): Dict of number of levels in each qubit system.
     level (int): Level of qubit system to be measured.
 
     Returns
     -------
     out_oper (qutip.Qobj): Occupation number operator for target qubit.
     """
+    # reverse sort by index
+    rev_h_osc = sorted(h_osc.items(), key=lambda x: x[0])[::-1]
+    rev_h_qub = sorted(h_qub.items(), key=lambda x: x[0])[::-1]
 
-    # qubit_0 * … * qubit_n * osc_0 * … * osc_n
+    # osc_n * … * osc_0 * qubit_n * … * qubit_0
     opers = []
-    for ii, dd in enumerate(h_qub):
+    for ii, dd in rev_h_osc:
+        opers.append(op.qeye(dd))
+    for ii, dd in rev_h_qub:
         if ii == target_qubit:
-            opers.append(op.fock_dm(dd, level))
+            opers.append(op.fock_dm(h_qub.get(target_qubit, 2), level))
         else:
             opers.append(op.qeye(dd))
-    for ii, dd in enumerate(h_osc):
-        opers.append(op.qeye(dd))
 
-    # return in reverse order
-    # osc_n * … * osc_0 * qubit_n * … * qubit_0
-    return op.tensor(opers[::-1])
+    return op.tensor(opers)
 
 
 def measure_outcomes(measured_qubits, state_vector, measure_ops,
@@ -136,8 +139,8 @@ def apply_projector(measured_qubits, results, h_qub, h_osc, state_vector):
     ----------
     measured_qubits (list): measured qubit indices.
     results (list): results of qubit measurements.
-    h_osc (list): Dimension of oscillator subspace
-    h_qub (list): Dimension of qubit subspace
+    h_qub (dict): Dict of number of levels in each qubit system.
+    h_osc (dict): Dict of number of levels in each oscillator.
     state_vector (ndarray): State vector.
 
     Returns:
@@ -145,20 +148,21 @@ def apply_projector(measured_qubits, results, h_qub, h_osc, state_vector):
     proj_state (qutip.Qobj): State vector after projector applied, and normalized.
     """
 
-    # qubit_0 * … * qubit_n * osc_0 * … * osc_n
+    # reverse sort by index
+    rev_h_osc = sorted(h_osc.items(), key=lambda x: x[0])[::-1]
+    rev_h_qub = sorted(h_qub.items(), key=lambda x: x[0])[::-1]
+
+    # osc_n * … * osc_0 * qubit_n * … * qubit_0
     opers = []
-    for ii, dd in enumerate(h_qub):
+    for ii, dd in rev_h_osc:
+        opers.append(op.qeye(dd))
+    for ii, dd in rev_h_qub:
         if ii in measured_qubits:
             opers.append(op.fock_dm(dd, results[ii]))
         else:
             opers.append(op.qeye(dd))
-    for ii, dd in enumerate(h_osc):
-        opers.append(op.qeye(dd))
-
-    # return in reverse order
-    # osc_n * … * osc_0 * qubit_n * … * qubit_0
-    proj_oper = op.tensor(opers[::-1])
 
+    proj_oper = op.tensor(opers)
     psi = op.opr_apply(proj_oper, state_vector)
     psi /= la.norm(psi)
 
@@ -177,12 +181,13 @@ def init_fock_state(h_osc, h_qub, noise_dict={}):
     Returns:
         qutip.Qobj: State vector
     """
+    # reverse sort by index
+    rev_h_osc = sorted(h_osc.items(), key=lambda x: x[0])[::-1]
+    rev_h_qub = sorted(h_qub.items(), key=lambda x: x[0])[::-1]
 
-    # qubit_0 * … * qubit_n * osc_0 * … * osc_n
+    # osc_n * … * osc_0 * qubit_n * … * qubit_0
     sub_state_vecs = []
-    for ii, dd in enumerate(h_qub):
-        sub_state_vecs.append(op.basis(dd, 0))
-    for ii, dd in enumerate(h_osc):
+    for ii, dd in rev_h_osc:
         n_thermal = noise_dict['oscillator']['n_th'].get(str(ii), 0)
         if n_thermal == 0:
             # no thermal particles
@@ -196,7 +201,7 @@ def init_fock_state(h_osc, h_qub, noise_dict={}):
             cum_sum = np.cumsum(diags)
             idx = np.where(np.random.random() < cum_sum)[0][0]
         sub_state_vecs.append(op.basis(dd, idx))
+    for ii, dd in rev_h_qub:
+        sub_state_vecs.append(op.basis(dd, 0))
 
-    # return in reverse order
-    # osc_n * … * osc_0 * qubit_n * … * qubit_0
-    return op.tensor(sub_state_vecs[::-1])
+    return op.tensor(sub_state_vecs)
diff --git a/qiskit/providers/aer/openpulse/qobj/opparse.py b/qiskit/providers/aer/openpulse/qobj/opparse.py
index 9541c99010..e8b79f353d 100644
--- a/qiskit/providers/aer/openpulse/qobj/opparse.py
+++ b/qiskit/providers/aer/openpulse/qobj/opparse.py
@@ -41,15 +41,15 @@
 class HamiltonianParser:
     """ Generate QuTip hamiltonian object from string
     """
-    def __init__(self, hamiltonian, dim_osc, dim_qub):
+    def __init__(self, h_str, dim_osc, dim_qub):
         """ Create new quantum operator generator
 
         Parameters:
-            hamiltonian (list): list of Hamiltonian string
-            dim_osc (list): dimension of oscillator subspace
-            dim_qub (list): dimension of qubit subspace
+            h_str (list): list of Hamiltonian string
+            dim_osc (dict): dimension of oscillator subspace
+            dim_qub (dict): dimension of qubit subspace
         """
-        self.hamiltonian = hamiltonian
+        self.h_str = h_str
         self.dim_osc = dim_osc
         self.dim_qub = dim_qub
         self.__td_hams = []
@@ -72,7 +72,7 @@ def parse(self):
         self._expand_sum()
 
         # convert to reverse Polish notation
-        for ham in self.hamiltonian:
+        for ham in self.h_str:
             if len(re.findall(r"\|\|", ham)) > 1:
                 raise Exception("Multiple time-dependent terms in %s" % ham)
             p_td = re.search(r"(?P<opr>[\S]+)\|\|(?P<ch>[\S]+)", ham)
@@ -102,7 +102,7 @@ def _expand_sum(self):
         sum_str = re.compile(r"_SUM\[(?P<itr>[a-z]),(?P<l>[a-z\d{}+-]+),(?P<u>[a-z\d{}+-]+),")
         brk_str = re.compile(r"]")
 
-        ham_list = copy.copy(self.hamiltonian)
+        ham_list = copy.copy(self.h_str)
         ham_out = []
 
         while any(ham_list):
@@ -144,7 +144,7 @@ def _expand_sum(self):
                                           trg_s, ham[p_brks[ii].end():]]))
                 ham_list.extend(_temp)
 
-        self.hamiltonian = ham_out
+        self.h_str = ham_out
 
         return ham_out
 
@@ -279,18 +279,21 @@ def _token2qobj(self, tokens):
 
 class NoiseParser:
     """ Generate QuTip noise object from dictionary
-    Qubit noise is given in the format of list of dictionary:
-
-        "qubit": [
-            {"Sm": 0.006}
-        ]
-
-    and oscillator noise is given in the nested list of (n_th, coupling):
-
-        "oscillator": [
-            [0.001, 0.05]
-        ]
-
+    Qubit noise is given in the format of nested dictionary:
+        "qubit": {
+            "0": {
+                "Sm": 0.006
+            }
+        }
+    and oscillator noise is given in the format of nested dictionary:
+        "oscillator": {
+            "n_th": {
+                "0": 0.001
+            },
+            "coupling": {
+                "0": 0.05
+            }
+        }
     these configurations are combined in the same dictionary
     """
     def __init__(self, noise_dict, dim_osc, dim_qub):
@@ -298,11 +301,11 @@ def __init__(self, noise_dict, dim_osc, dim_qub):
 
         Parameters:
             noise_dict (dict): dictionary of noise configuration
-            dim_osc (list): dimension of oscillator subspace
-            dim_qub (list): dimension of qubit subspace
+            dim_osc (dict): dimension of oscillator subspace
+            dim_qub (dict): dimension of qubit subspace
         """
-        self.noise_osc = noise_dict.get('oscillator', [])
-        self.noise_qub = noise_dict.get('qubit', [])
+        self.noise_osc = noise_dict.get('oscillator', {'n_th': {}, 'coupling': {}})
+        self.noise_qub = noise_dict.get('qubit', {})
         self.dim_osc = dim_osc
         self.dim_qub = dim_qub
         self.__c_list = []
@@ -317,7 +320,7 @@ def parse(self):
         """ Parse and generate quantum class object
         """
         # Qubit noise
-        for index, config in enumerate(self.noise_qub):
+        for index, config in self.noise_qub.items():
             for opname, coef in config.items():
                 # TODO: support noise in multi-dimensional system
                 # TODO: support noise with math operation
@@ -327,16 +330,19 @@ def parse(self):
                     raise Exception('Unsupported noise operator %s is given' % opname)
                 self.__c_list.append(np.sqrt(coef) * opr)
         # Oscillator noise
-        ndic = [nconf[0] for nconf in self.noise_osc]
-        cdic = [nconf[1] for nconf in self.noise_osc]
-        for ii, (n_coef, c_coef) in enumerate(zip(ndic, cdic)):
-            if c_coef > 0:
-                opr = gen_oper('A', int(ii), self.dim_osc, self.dim_qub)
-                if n_coef:
-                    self.__c_list.append(np.sqrt(c_coef * (1 + n_coef)) * opr)
-                    self.__c_list.append(np.sqrt(c_coef * n_coef) * opr.dag())
-                else:
-                    self.__c_list.append(np.sqrt(c_coef) * opr)
+        ndic = self.noise_osc['n_th']
+        cdic = self.noise_osc['coupling']
+        for (n_ii, n_coef), (c_ii, c_coef) in zip(ndic.items(), cdic.items()):
+            if n_ii == c_ii:
+                if c_coef > 0:
+                    opr = gen_oper('A', int(n_ii), self.dim_osc, self.dim_qub)
+                    if n_coef:
+                        self.__c_list.append(np.sqrt(c_coef * (1 + n_coef)) * opr)
+                        self.__c_list.append(np.sqrt(c_coef * n_coef) * opr.dag())
+                    else:
+                        self.__c_list.append(np.sqrt(c_coef) * opr)
+            else:
+                raise Exception('Invalid oscillator index in noise dictionary.')
 
 
 def math_priority(o1, o2):