diff --git a/qiskit/providers/aer/__init__.py b/qiskit/providers/aer/__init__.py
index b76fdb31db..0090e198d2 100644
--- a/qiskit/providers/aer/__init__.py
+++ b/qiskit/providers/aer/__init__.py
@@ -34,6 +34,7 @@
     QasmSimulator
     StatevectorSimulator
     UnitarySimulator
+    PulseSimulator
 
 Job Class
 =========
@@ -42,6 +43,15 @@
 
    AerJob
 
+OpenPulse
+=========
+
+.. autosummary::
+    :toctree: ../stubs/
+
+    PulseSystemModel
+    duffing_system_model
+
 Exceptions
 ==========
 .. autosummary::
@@ -64,6 +74,7 @@
 from .aerjob import AerJob
 from .aererror import AerError
 from .backends import *
+from .openpulse import *
 from . import noise
 from . import utils
 from .version import __version__
diff --git a/qiskit/providers/aer/backends/__init__.py b/qiskit/providers/aer/backends/__init__.py
index 79fccf46e0..669afa8a78 100644
--- a/qiskit/providers/aer/backends/__init__.py
+++ b/qiskit/providers/aer/backends/__init__.py
@@ -17,3 +17,4 @@
 from .qasm_simulator import QasmSimulator
 from .statevector_simulator import StatevectorSimulator
 from .unitary_simulator import UnitarySimulator
+from .pulse_simulator import PulseSimulator
diff --git a/qiskit/providers/aer/backends/pulse_simulator.py b/qiskit/providers/aer/backends/pulse_simulator.py
index a41d32223c..fa61e08e97 100644
--- a/qiskit/providers/aer/backends/pulse_simulator.py
+++ b/qiskit/providers/aer/backends/pulse_simulator.py
@@ -1,6 +1,6 @@
 # This code is part of Qiskit.
 #
-# (C) Copyright IBM 2018, 2019.
+# (C) Copyright IBM 2018, 2019, 2020.
 #
 # This code is licensed under the Apache License, Version 2.0. You may
 # obtain a copy of this license in the LICENSE.txt file in the root directory
@@ -32,74 +32,57 @@
 
 
 class PulseSimulator(AerBackend):
-    """
-    Aer OpenPulse simulator
-
-    The `PulseSimulator` simulates pulse `Schedules` on a model of a quantum system, where a model
-    is specified by a `PulseSystemModel` object, which stores Hamiltonian and control channel
-    information. Simulation is performed in the rotating frame of the drift Hamiltonian in the
-    `PulseSystemModel`.
-
-    `PulseSystemModel` objects can be constructed from backends with a hamiltonian description, or
-    from the `transmon_system_model` function.
+    r"""Aer OpenPulse simulator.
 
-    Results are returned in the same format as when jobs are submitted to actual devices.
+    The ``PulseSimulator`` simulates continuous time Hamiltonian dynamics of a quantum system,
+    with controls specified by pulse :class:`Schedule` objects, and the model of the physical
+    system specified by :class:`PulseSystemModel` objects. Results are returned in the same format
+    as when jobs are submitted to actual devices.
 
-    **Example of usage**
+    **Example**
 
-    To use the simulator, `assemble` a `PulseQobj` object from a list of pulse `Schedules`, using
-    `backend=PulseSimulator()`.
-
-    In the following, `schedules` is a list of pulse `Schedule` objects, and `system_model` is a
-    `PulseSystemModel` object.
+    To use the simulator, first :meth:`assemble` a :class:`PulseQobj` object
+    from a list of pulse :class:`Schedule` objects, using ``backend=PulseSimulator()``.
+    Call the simulator with the :class:`PulseQobj` and a :class:`PulseSystemModel`
+    object representing the physical system.
 
     .. code-block:: python
 
-        backend_sim = qiskit.Aer.get_backend('pulse_simulator')
+        backend_sim = qiskit.providers.aer.PulseSimulator()
 
-        # assemble pulse_qobj with backend=backend_sim
+        # Assemble schedules using PulseSimulator as the backend
         pulse_qobj = assemble(schedules, backend=backend_sim)
 
-        # Run simulation
+        # Run simulation on a PulseSystemModel object
         results = backend_sim.run(pulse_qobj, system_model)
 
-    **Important parameters**
-
-    * `qubit_lo_freq`: The local oscillator frequency for each `DriveChannel`. This can be drawn
-                       from several places, listed in order of importance:
-           * passed as an argument to `assemble`
-           * determined from PulseSystemModel attribute `_qubit_freq_est`, if not passed to assemble
-           * computed from the dressed energy gaps of the drift Hamiltonian
-
-    **Measurement and output**
+    **Supported PulseQobj parameters**
 
-    The measurement results are from projections of the state vector in dressed energy basis of
-    the drift Hamiltonian.
+    * ``qubit_lo_freq``: Local oscillator frequencies for each :class:`DriveChannel`.
+      Defaults to either the value given in the :class:`PulseSystemModel`, or
+      is calculated directly from the Hamiltonian.
+    * ``meas_level``: Type of desired measurement output, in ``[1, 2]``.
+      ``1`` gives complex numbers (IQ values), and ``2`` gives discriminated states ``|0>`` and
+      ``|1>``. Defaults to ``2``.
+    * ``meas_return``: Measurement type, ``'single'`` or ``'avg'``. Defaults to ``'avg'``.
+    * ``shots``: Number of shots per experiment. Defaults to ``1024``.
 
-    There are three measurement levels that return the data, specified when using assemble.
-    Measurement level `0` gives the raw data.
-    Measurement level `1` gives complex numbers (IQ values).
-    Measurement level `2` gives the discriminated states, `|0>` and `|1>`.
 
-    **Simulation method**
+    **Simulation details**
 
-    The simulator uses the `zvode` differential equation solver method through `scipy`.
+    The simulator uses the ``zvode`` differential equation solver method through ``scipy``.
+    Simulation is performed in the rotating frame of the diagonal of the drift Hamiltonian
+    contained in the :class:`PulseSystemModel`. Measurements are performed in the `dressed basis`
+    of the drift Hamiltonian.
 
     **Other options**
 
-    The `run` function additionally takes an argument `backend_options` for additional
-    customization. It accepts keys:
+    :meth:`PulseSimulator.run` takes an additional ``dict`` argument ``backend_options`` for
+    customization. Accepted keys:
 
-    * `'ode_options'`: a dictionary containing options to pass to the `zvode` solver. Accepted keys
-      for this option are `'atol'`, `'rtol'`, `'nsteps'`, `'max_step'`, `'num_cpus'`, `'norm_tol'`,
-      and `'norm_steps'`
-
-    **Default behaviors**
-
-    Defaults filled in for `assemble` parameters if not specified:
-    * `meas_level`: `2`
-    * `meas_return`: `'avg'`
-    * `shots`: `1024`
+    * ``'ode_options'``: A ``dict`` for ``zvode`` solver options. Accepted keys
+      are ``'atol'``, ``'rtol'``, ``'nsteps'``, ``'max_step'``, ``'num_cpus'``, ``'norm_tol'``,
+      and ``'norm_steps'``.
     """
 
     DEFAULT_CONFIGURATION = {
@@ -133,7 +116,17 @@ def __init__(self, configuration=None, provider=None):
                          provider=provider)
 
     def run(self, qobj, system_model, backend_options=None, validate=False):
-        """Run a qobj on the backend."""
+        """Run a qobj on system_model.
+
+        Args:
+            qobj (PulseQobj): Qobj for pulse Schedules to run
+            system_model (PulseSystemModel): Physical model to run simulation on
+            backend_options (dict): Other options
+            validate (bool): Flag for validation checks
+
+        Returns:
+            Result: results of simulation
+        """
         # Submit job
         job_id = str(uuid.uuid4())
         aer_job = AerJob(self, job_id, self._run_job, qobj, system_model,
diff --git a/qiskit/providers/aer/openpulse/__init__.py b/qiskit/providers/aer/openpulse/__init__.py
index 74d5c7b286..fa21b1638a 100644
--- a/qiskit/providers/aer/openpulse/__init__.py
+++ b/qiskit/providers/aer/openpulse/__init__.py
@@ -18,6 +18,9 @@
 import numpy as np
 from .qutip_lite.cy import pyxbuilder as pbldr
 
+from .duffing_model_generators import duffing_system_model
+from .pulse_system_model import PulseSystemModel
+
 # Remove -Wstrict-prototypes from cflags
 CFG_VARS = distutils.sysconfig.get_config_vars()
 if "CFLAGS" in CFG_VARS:
diff --git a/qiskit/providers/aer/openpulse/duffing_model_generators.py b/qiskit/providers/aer/openpulse/duffing_model_generators.py
index 3e0e4dced6..1d3aa4cf14 100644
--- a/qiskit/providers/aer/openpulse/duffing_model_generators.py
+++ b/qiskit/providers/aer/openpulse/duffing_model_generators.py
@@ -27,50 +27,92 @@ def duffing_system_model(dim_oscillators,
                          drive_strengths,
                          coupling_dict,
                          dt):
-    """Returns a PulseSystemModel for a specified Duffing oscillator system.
-
-    Note: Frequencies are assumed to be in frequency units (as opposed to radial).
-
-    Single oscillators are specified by three parameters: frequency v, anharmonicity alpha, and
-    drive strength r, which enter into the Hamiltonian model via the terms:
-        pi*(2*v-alpha)*O + pi*alpha*O*O + 2*pi*r*X*D(t),
-    where O is the number operator, X is the drive operator, and D(t) is the signal of the
-    corresponding DriveChannel.
-
-    Couplings between oscillators are specified in the argument coupling_dict, with keys being
-    edges, and values being the coupling strenghts. Couplings enter the Hamiltonian model via an
-    exchange coupling term:
-        2*pi*j*(A0*C1+C0*A1),
-    where j is the coupling strength, and A0*C1+C0*A1 is the exchange coupling operator between the
-    two oscillators, where A is the annihiliation operator, and C is the creation operator.
-
-    For each pair of coupled oscillators, pulse ControlChannels are prepared for doing cross
-    resonance (CR) drives between oscillators. A CR drive on oscillators 0 with target oscillators
-    1 is modeled with the Hamiltonian term:
-        2*pi*r0*X0*U(t),
-    where r0 is the drive strength for oscillators 0, X0 is the drive operator for oscillators 0,
-    and U(t) is the signal for the CR drive. U(t) is set to have carrier frequency equal to that of
-    the target oscillators.
-
-    Indices for ControlChannels corresponding to CR drives can be retrieved from the returned
-    PulseSystemModel using the control_channel_index() method, with keys being tuples. E.g. the
-    index for the ControlChannel corresponding to a CR drive on oscillator 0 with target 1 is
-    given by system_model.control_channel_index((0,1)).
+    r"""Returns a :class:`PulseSystemModel` representing a physical model for a
+    collection of Duffing oscillators.
 
-    Args:
-        dim_oscillators (int): Dimension of truncation for each oscillator
-        oscillator_freqs (list): Oscillator frequencies in frequency units
-        anharm_freqs (list): Anharmonicity values in frequency units
-        drive_strengths (list): Drive strength values in frequency units
-        coupling_dict (dict): Specification of the coupling graph with keys being edges, and values
-                              the coupling strengths in frequency units.
+    In the model, each individual oscillator is specified by the parameters:
+
+        * Frequency: :math:`\nu`, specified in the list ``oscillator_freqs``
+        * Anharmonicity: :math:`\alpha`, specified in the list ``anharm_freqs``, and
+        * Drive strength: :math:`r`, specified in the list ``drive_strengths``.
+
+    For each oscillator, the above parameters enter into the Hamiltonian via the terms:
+
+    .. math::
+        \pi(2 \nu - \alpha)a^\dagger a +
+        \pi \alpha (a^\dagger a)^2 + 2 \pi r D(t) (a + a^\dagger),
+
+    where :math:`a^\dagger` and :math:`a` are, respectively, the creation and annihilation
+    operators for the oscillator, and :math:`D(t)` is the drive signal for the oscillator.
+
+    Each coupling term between a pair of oscillators is specified by:
+
+        * Oscillator pair: :math:`(i,k)`, and
+        * Coupling strength: :math:`j`,
+
+    which are passed in the argument ``coupling_dict``, which is a ``dict`` with keys
+    being the ``tuple`` ``(i,k)``, and values the strength ``j``. Specifying a coupling
+    results in the Hamiltonian term:
+
+    .. math::
+        2 \pi j (a_i^\dagger a_k + a_i a_k^\dagger).
+
+    Finally, the returned :class:`PulseSystemModel` is setup for performing cross-resonance
+    drives between coupled qubits. The index for the :class:`ControlChannel` corresponding
+    to a particular cross-resonance drive channel is retreived by calling
+    :meth:`PulseSystemModel.control_channel_index` with the tuple ``(drive_idx, target_idx)``,
+    where ``drive_idx`` is the index of the oscillator being driven, and ``target_idx`` is
+    the target oscillator (see example below).
+
+    Note: In this model, all frequencies are in frequency units (as opposed to radial).
+
+    **Example**
+
+    Constructing a three Duffing Oscillator :class:``PulseSystemModel``.
 
-                              For example:
+    .. code-block:: python
 
-                              * ``{(0,1): 0.02, (1,3): 0.01}`` specifies a system in which
-                                oscillators (0,1) are coupled with strength 0.02, and (1,3) are
-                                coupled with strength 0.01
-        dt (float): Sample width for pulse instructions
+        # cutoff dimensions
+        dim_oscillators = 3
+
+        # single oscillator drift parameters
+        oscillator_freqs = [5.0e9, 5.1e9, 5.2e9]
+        anharm_freqs = [-0.33e9, -0.33e9, -0.33e9]
+
+        # drive strengths
+        drive_strengths = [0.02e9, 0.02e9, 0.02e9]
+
+        # specify coupling as a dictionary; here the qubit pair (0,1) is coupled with
+        # strength 0.002e9, and the qubit pair (1,2) is coupled with strength 0.001e9
+        coupling_dict = {(0,1): 0.002e9, (1,2): 0.001e9}
+
+        # time
+        dt = 1e-9
+
+        # create the model
+        three_qubit_model = duffing_system_model(dim_oscillators=dim_oscillators,
+                                                 oscillator_freqs=oscillator_freqs,
+                                                 anharm_freqs=anharm_freqs,
+                                                 drive_strengths=drive_strengths,
+                                                 coupling_dict=coupling_dict,
+                                                 dt=dt)
+
+    In the above model, qubit pairs (0,1) and (1,2) are coupled. To perform a
+    cross-resonance drive on qubit 1 with target 0, use the :class:`ControlChannel`
+    with index:
+
+    .. code-block:: python
+
+        three_qubit_model.control_channel_index((1,0))
+
+    Args:
+        dim_oscillators (int): Dimension of truncation for each oscillator.
+        oscillator_freqs (list): Oscillator frequencies in frequency units.
+        anharm_freqs (list): Anharmonicity values in frequency units.
+        drive_strengths (list): Drive strength values in frequency units.
+        coupling_dict (dict): Coupling graph with keys being edges, and values
+                              the coupling strengths in frequency units.
+        dt (float): Sample width for pulse instructions.
 
     Returns:
         PulseSystemModel: The generated Duffing system model
@@ -127,7 +169,7 @@ def duffing_system_model(dim_oscillators,
     return PulseSystemModel(hamiltonian=hamiltonian_model,
                             u_channel_lo=u_channel_lo,
                             control_channel_labels=coupling_graph.sorted_two_way_graph,
-                            qubit_list=oscillators,
+                            subsystem_list=oscillators,
                             dt=dt)
 
 
diff --git a/qiskit/providers/aer/openpulse/hamiltonian_model.py b/qiskit/providers/aer/openpulse/hamiltonian_model.py
index 54a9742aea..294a80aa40 100644
--- a/qiskit/providers/aer/openpulse/hamiltonian_model.py
+++ b/qiskit/providers/aer/openpulse/hamiltonian_model.py
@@ -15,7 +15,6 @@
 
 "HamiltonianModel class for system specification for the PulseSimulator"
 
-from warnings import warn
 from collections import OrderedDict
 import numpy as np
 import numpy.linalg as la
@@ -29,15 +28,13 @@ class HamiltonianModel():
     def __init__(self,
                  system=None,
                  variables=None,
-                 qubit_dims=None,
-                 oscillator_dims=None):
+                 subsystem_dims=None):
         """Initialize a Hamiltonian model.
 
         Args:
             system (list): List of Qobj objects representing operator form of the Hamiltonian.
             variables (OrderedDict): Ordered dict for parameter values in Hamiltonian.
-            qubit_dims (dict): dict of qubit dimensions.
-            oscillator_dims (dict): dict of oscillator dimensions.
+            subsystem_dims (dict): dict of subsystem dimensions.
 
         Raises:
             ValueError: if arguments are invalid.
@@ -50,9 +47,7 @@ def __init__(self,
         self._variables = variables
         # Channels in the Hamiltonian string
         # Qubit subspace dimensinos
-        self._qubit_dims = qubit_dims or {}
-        # Oscillator subspace dimensions
-        self._oscillator_dims = oscillator_dims or {}
+        self._subsystem_dims = subsystem_dims or {}
 
         # The rest are computed from the previous
 
@@ -73,12 +68,12 @@ def __init__(self,
         self._compute_drift_data()
 
     @classmethod
-    def from_dict(cls, hamiltonian, qubit_list=None):
+    def from_dict(cls, hamiltonian, subsystem_list=None):
         """Initialize from a Hamiltonian string specification.
 
         Args:
             hamiltonian (dict): dictionary representing Hamiltonian in string specification.
-            qubit_list (list or None): List of qubits to extract from the hamiltonian.
+            subsystem_list (list or None): List of subsystems to extract from the hamiltonian.
 
         Returns:
             HamiltonianModel: instantiated from hamiltonian dictionary
@@ -94,15 +89,15 @@ def from_dict(cls, hamiltonian, qubit_list=None):
 
         # Get qubit subspace dimensions
         if 'qub' in hamiltonian:
-            if qubit_list is None:
-                qubit_list = [int(qubit) for qubit in hamiltonian['qub']]
+            if subsystem_list is None:
+                subsystem_list = [int(qubit) for qubit in hamiltonian['qub']]
 
-            qubit_dims = {
+            subsystem_dims = {
                 int(key): val
                 for key, val in hamiltonian['qub'].items()
             }
         else:
-            qubit_dims = {}
+            subsystem_dims = {}
 
         # Get oscillator subspace dimensions
         if 'osc' in hamiltonian:
@@ -116,11 +111,11 @@ def from_dict(cls, hamiltonian, qubit_list=None):
         # Parse the Hamiltonian
         system = HamiltonianParser(h_str=hamiltonian['h_str'],
                                    dim_osc=oscillator_dims,
-                                   dim_qub=qubit_dims)
-        system.parse(qubit_list)
+                                   dim_qub=subsystem_dims)
+        system.parse(subsystem_list)
         system = system.compiled
 
-        return cls(system, variables, qubit_dims, oscillator_dims)
+        return cls(system, variables, subsystem_dims)
 
     def get_qubit_lo_from_drift(self):
         """ Computes a list of qubit frequencies corresponding to the exact energy
@@ -129,13 +124,13 @@ def get_qubit_lo_from_drift(self):
         Returns:
             qubit_lo_freq (list): the list of frequencies
         """
-        qubit_lo_freq = [0] * len(self._qubit_dims)
+        qubit_lo_freq = [0] * len(self._subsystem_dims)
 
         # compute difference between first excited state of each qubit and
         # the ground energy
         min_eval = np.min(self._evals)
-        for q_idx in range(len(self._qubit_dims)):
-            single_excite = _first_excited_state(q_idx, self._qubit_dims)
+        for q_idx in range(len(self._subsystem_dims)):
+            single_excite = _first_excited_state(q_idx, self._subsystem_dims)
             dressed_eval = _eval_for_max_espace_overlap(
                 single_excite, self._evals, self._estates)
             qubit_lo_freq[q_idx] = (dressed_eval - min_eval) / (2 * np.pi)
@@ -223,12 +218,6 @@ def _compute_drift_data(self):
             evals_mapped[pos] = evals[i]
             estates_mapped[:, pos] = estate
 
-        overlap_threshold = 0.6
-        if min_overlap < overlap_threshold:
-            warn('Warning: The minimum overlap of an eigenstate of the drift to an element of '
-                 'the computational basis is below ' +
-                 '{0}, and may result in unexpected behavior.'.format(str(overlap_threshold)))
-
         self._evals = evals_mapped
         self._estates = estates_mapped
         self._h_diag = np.ascontiguousarray(np.diag(ham_full).real)
@@ -247,7 +236,7 @@ def _hamiltonian_parse_exceptions(hamiltonian):
         raise AerError('Oscillator-type systems are not supported.')
 
 
-def _first_excited_state(qubit_idx, qubit_dims):
+def _first_excited_state(qubit_idx, subsystem_dims):
     """
     Returns the vector corresponding to all qubits in the 0 state, except for
     qubit_idx in the 1 state.
@@ -259,18 +248,18 @@ def _first_excited_state(qubit_idx, qubit_dims):
     Parameters:
         qubit_idx (int): the qubit to be in the 1 state
 
-        qubit_dims (dict): a dictionary with keys being qubit index, and
-                        value being the dimension of the qubit
+        subsystem_dims (dict): a dictionary with keys being subsystem index, and
+                        value being the dimension of the subsystem
 
     Returns:
         vector: the state with qubit_idx in state 1, and the rest in state 0
     """
     vector = np.array([1.])
     # iterate through qubits, tensoring on the state
-    qubit_indices = [int(qubit) for qubit in qubit_dims]
+    qubit_indices = [int(qubit) for qubit in subsystem_dims]
     qubit_indices.sort()
     for idx in qubit_indices:
-        new_vec = np.zeros(qubit_dims[idx])
+        new_vec = np.zeros(subsystem_dims[idx])
         if idx == qubit_idx:
             new_vec[1] = 1
         else:
diff --git a/qiskit/providers/aer/openpulse/pulse_system_model.py b/qiskit/providers/aer/openpulse/pulse_system_model.py
index c0b32b5697..66d37d6af6 100644
--- a/qiskit/providers/aer/openpulse/pulse_system_model.py
+++ b/qiskit/providers/aer/openpulse/pulse_system_model.py
@@ -23,7 +23,35 @@
 
 
 class PulseSystemModel():
-    """PulseSystemModel containing all model parameters necessary for simulation.
+    r"""Physical model object for use in :class:`PulseSimulator`.
+
+    This class contains model information required by :class:`PulseSimulator`. It contains:
+
+        * ``"hamiltonian"``: a :class:`HamiltonianModel` object representing the
+          Hamiltonian of the system.
+        * ``"qubit_freq_est"`` and ``"meas_freq_est"``: optional default values for
+          qubit and measurement frequencies.
+        * ``"u_channel_lo"``: A description of :class:`ControlChannel` local oscillator
+          frequencies in terms of qubit local oscillator frequencies.
+        * ``"control_channel_labels"``: Optional list of identifying information for
+          each :class:`ControlChannel` that the model supports.
+        * ``"subsystem_list"``: List of subsystems in the model.
+        * ``"dt"``: Sample width size for OpenPulse instructions.
+
+    A :class:`PulseSystemModel` object can be instantiated from the
+    helper function :meth:`duffing_system_model`, or using the
+    :meth:`PulseSystemModel.from_backend` constructor.
+
+    **Example**
+
+    Constructing from a backend:
+
+    .. code-block: python
+
+        provider = IBMQ.load_account()
+        armonk_backend = provider.get_backend('ibmq_armonk')
+
+        system_model = PulseSystemModel.from_backend(armonk_backend)
     """
     def __init__(self,
                  hamiltonian=None,
@@ -31,7 +59,7 @@ def __init__(self,
                  meas_freq_est=None,
                  u_channel_lo=None,
                  control_channel_labels=None,
-                 qubit_list=None,
+                 subsystem_list=None,
                  dt=None):
         """Initialize a PulseSystemModel.
 
@@ -44,7 +72,7 @@ def __init__(self,
             u_channel_lo (list): list of ControlChannel frequency specifications.
             control_channel_labels (list): list of labels for control channels, which can be of
                                            any type.
-            qubit_list (list): list of valid qubit indicies for the model.
+            subsystem_list (list): list of valid qubit indicies for the model.
             dt (float): pixel size for pulse Instructions.
         Raises:
             AerError: if hamiltonian is not None or a HamiltonianModel
@@ -60,16 +88,16 @@ def __init__(self,
         self.hamiltonian = hamiltonian
         self.u_channel_lo = u_channel_lo
         self.control_channel_labels = control_channel_labels or []
-        self.qubit_list = qubit_list
+        self.subsystem_list = subsystem_list
         self.dt = dt
 
     @classmethod
-    def from_backend(cls, backend, qubit_list=None):
-        """Returns a PulseSystemModel constructed from a backend object.
+    def from_backend(cls, backend, subsystem_list=None):
+        """Returns a PulseSystemModel constructed from an OpenPulse enabled backend object.
 
         Args:
             backend (Backend): backend object to draw information from.
-            qubit_list (list): a list of ints for which qubits to include in the model.
+            subsystem_list (list): a list of ints for which qubits to include in the model.
 
         Returns:
             PulseSystemModel: the PulseSystemModel constructed from the backend.
@@ -93,10 +121,10 @@ def from_backend(cls, backend, qubit_list=None):
         meas_freq_est = defaults.get('meas_freq_est', None)
 
         # draw from configuration
-        # if no qubit_list, use all for device
-        qubit_list = qubit_list or list(range(config['n_qubits']))
+        # if no subsystem_list, use all for device
+        subsystem_list = subsystem_list or list(range(config['n_qubits']))
         ham_string = config['hamiltonian']
-        hamiltonian = HamiltonianModel.from_dict(ham_string, qubit_list)
+        hamiltonian = HamiltonianModel.from_dict(ham_string, subsystem_list)
         u_channel_lo = config.get('u_channel_lo', None)
         dt = config.get('dt', None)
 
@@ -136,7 +164,7 @@ def from_backend(cls, backend, qubit_list=None):
                    meas_freq_est=meas_freq_est,
                    u_channel_lo=u_channel_lo,
                    control_channel_labels=control_channel_labels,
-                   qubit_list=qubit_list,
+                   subsystem_list=subsystem_list,
                    dt=dt)
 
     def control_channel_index(self, label):
@@ -155,7 +183,7 @@ def control_channel_index(self, label):
             return self.control_channel_labels.index(label)
 
     def calculate_channel_frequencies(self, qubit_lo_freq=None):
-        """Calculate frequencies for each channel.
+        """Calculate frequencies for each channel given qubit_lo_freq.
 
         Args:
             qubit_lo_freq (list or None): list of qubit linear
diff --git a/qiskit/providers/aer/openpulse/qobj/digest.py b/qiskit/providers/aer/openpulse/qobj/digest.py
index 69542df63f..4d4ff60c08 100644
--- a/qiskit/providers/aer/openpulse/qobj/digest.py
+++ b/qiskit/providers/aer/openpulse/qobj/digest.py
@@ -106,7 +106,7 @@ def digest_pulse_obj(qobj, system_model, backend_options=None):
     # ###############################
 
     # Get qubit list and number
-    qubit_list = system_model.qubit_list
+    qubit_list = system_model.subsystem_list
     if qubit_list is None:
         raise ValueError('Model must have a qubit list to simulate.')
     n_qubits = len(qubit_list)
@@ -123,8 +123,8 @@ def digest_pulse_obj(qobj, system_model, backend_options=None):
     out.h_diag = ham_model._h_diag
     out.evals = ham_model._evals
     out.estates = ham_model._estates
-    dim_qub = ham_model._qubit_dims
-    dim_osc = ham_model._oscillator_dims
+    dim_qub = ham_model._subsystem_dims
+    dim_osc = {}
     out.freqs = system_model.calculate_channel_frequencies(qubit_lo_freq=qubit_lo_freq)
     # convert estates into a Qutip qobj
     estates = [op.state(state) for state in ham_model._estates.T[:]]
diff --git a/test/terra/backends/test_pulse_simulator.py b/test/terra/backends/test_pulse_simulator.py
index c5cb7c1da1..8167af91d3 100644
--- a/test/terra/backends/test_pulse_simulator.py
+++ b/test/terra/backends/test_pulse_simulator.py
@@ -641,12 +641,12 @@ def _system_model_1Q(self, omega_0, omega_a, qubit_dim=2):
         ham_model = HamiltonianModel.from_dict(hamiltonian)
 
         u_channel_lo = []
-        qubit_list = [0]
+        subsystem_list = [0]
         dt = 1.
 
         return PulseSystemModel(hamiltonian=ham_model,
                                 u_channel_lo=u_channel_lo,
-                                qubit_list=qubit_list,
+                                subsystem_list=subsystem_list,
                                 dt=dt)
 
     def _system_model_2Q(self, omega_0, omega_a, omega_i, qubit_dim=2):
@@ -678,12 +678,12 @@ def _system_model_2Q(self, omega_0, omega_a, omega_i, qubit_dim=2):
 
         u_channel_lo = [[{'q': 0, 'scale': [1.0, 0.0]}],
                         [{'q': 0, 'scale': [-1.0, 0.0]}, {'q': 1, 'scale': [1.0, 0.0]}]]
-        qubit_list = [0, 1]
+        subsystem_list = [0, 1]
         dt = 1.
 
         return PulseSystemModel(hamiltonian=ham_model,
                                 u_channel_lo=u_channel_lo,
-                                qubit_list=qubit_list,
+                                subsystem_list=subsystem_list,
                                 dt=dt)
 
     def _simple_1Q_schedule(self, phi, total_samples, shape="square", gauss_sigma=0):
diff --git a/test/terra/openpulse/test_duffing_model_generators.py b/test/terra/openpulse/test_duffing_model_generators.py
index 09de8100d0..3bb9aebb39 100644
--- a/test/terra/openpulse/test_duffing_model_generators.py
+++ b/test/terra/openpulse/test_duffing_model_generators.py
@@ -46,7 +46,7 @@ def test_duffing_system_model1(self):
         cr_idx_dict = {label: idx for idx, label in enumerate(system_model.control_channel_labels)}
 
         # check basic parameters
-        self.assertEqual(system_model.qubit_list, [0, 1])
+        self.assertEqual(system_model.subsystem_list, [0, 1])
         self.assertEqual(system_model.dt, 1.3)
 
         # check that cr_idx_dict is correct
@@ -70,7 +70,7 @@ def test_duffing_system_model1(self):
                          'r0': 1.1, 'r1': 1.2,
                          'j01': 0.02}
         self.assertEqual(ham_model._variables, expected_vars)
-        self.assertEqual(ham_model._qubit_dims, {0: 2, 1: 2})
+        self.assertEqual(ham_model._subsystem_dims, {0: 2, 1: 2})
         self._compare_str_lists(list(ham_model._channels), ['D0', 'D1', 'U0', 'U1'])
 
         # check that Hamiltonian terms have been imported correctly
@@ -127,7 +127,7 @@ def test_duffing_system_model2(self):
         cr_idx_dict = {label: idx for idx, label in enumerate(system_model.control_channel_labels)}
 
         # check basic parameters
-        self.assertEqual(system_model.qubit_list, [0, 1, 2])
+        self.assertEqual(system_model.subsystem_list, [0, 1, 2])
         self.assertEqual(system_model.dt, 1.3)
 
         # check that cr_idx_dict is correct
@@ -154,7 +154,7 @@ def test_duffing_system_model2(self):
                          'r0': 1.1, 'r1': 1.2, 'r2': 1.3,
                          'j01': 0.02, 'j12': 0.03}
         self.assertEqual(ham_model._variables, expected_vars)
-        self.assertEqual(ham_model._qubit_dims, {0: 3, 1: 3, 2: 3})
+        self.assertEqual(ham_model._subsystem_dims, {0: 3, 1: 3, 2: 3})
         self._compare_str_lists(list(ham_model._channels), ['D0', 'D1', 'D3', 'U0', 'U1', 'U2', 'U3'])
 
         # check that Hamiltonian terms have been imported correctly
@@ -221,7 +221,7 @@ def test_duffing_system_model3(self):
         cr_idx_dict = {label: idx for idx, label in enumerate(system_model.control_channel_labels)}
 
         # check basic parameters
-        self.assertEqual(system_model.qubit_list, [0, 1, 2, 3])
+        self.assertEqual(system_model.subsystem_list, [0, 1, 2, 3])
         self.assertEqual(system_model.dt, 1.3)
 
         # check that cr_idx_dict is correct
@@ -253,7 +253,7 @@ def test_duffing_system_model3(self):
                          'r0': 1.1, 'r1': 1.2, 'r2': 1.3, 'r3': 1.4,
                          'j01': 0.02, 'j02': 0.03, 'j03': 0.14, 'j12': 0.33, 'j13': 0.18}
         self.assertEqual(ham_model._variables, expected_vars)
-        self.assertEqual(ham_model._qubit_dims, {0: 2, 1: 2, 2: 2, 3: 2})
+        self.assertEqual(ham_model._subsystem_dims, {0: 2, 1: 2, 2: 2, 3: 2})
         self._compare_str_lists(list(ham_model._channels), ['D0', 'D1', 'D3', 'D4',
                                                             'U0', 'U1', 'U2', 'U3', 'U4',
                                                             'U5', 'U6', 'U7', 'U8', 'U9'])
diff --git a/test/terra/openpulse/test_pulse_digest.py b/test/terra/openpulse/test_pulse_digest.py
index e22573180f..e07d6278fb 100644
--- a/test/terra/openpulse/test_pulse_digest.py
+++ b/test/terra/openpulse/test_pulse_digest.py
@@ -134,12 +134,12 @@ def _system_model_2Q(self,
         # set up u channel freqs,
         u_channel_lo = [[{'q': 1, 'scale': [1.0, 0.0]}],
                         [{'q': 0, 'scale': [1.0, 0.0]}]]
-        qubit_list = [0, 1]
+        subsystem_list = [0, 1]
         dt = 1.
 
         return PulseSystemModel(hamiltonian=ham_model,
                                 u_channel_lo=u_channel_lo,
-                                qubit_list=qubit_list,
+                                subsystem_list=subsystem_list,
                                 dt=dt)
 
 
diff --git a/test/terra/openpulse/test_system_models.py b/test/terra/openpulse/test_system_models.py
index b98d4543fe..ed47168222 100644
--- a/test/terra/openpulse/test_system_models.py
+++ b/test/terra/openpulse/test_system_models.py
@@ -55,13 +55,13 @@ def _simple_system_model(self, v0=5.0, v1=5.1, j=0.01, r=0.02, alpha0=-0.33, alp
                                'alpha1': alpha1}
         ham_model = HamiltonianModel.from_dict(hamiltonian)
 
-        qubit_list =[0, 1]
+        subsystem_list =[0, 1]
         dt = 1.
 
         return PulseSystemModel(hamiltonian=ham_model,
                                 qubit_freq_est=self._default_qubit_lo_freq,
                                 u_channel_lo=self._u_channel_lo,
-                                qubit_list=qubit_list,
+                                subsystem_list=subsystem_list,
                                 dt=dt)
 
 
@@ -195,16 +195,7 @@ def test_eigen_sorting(self):
                       'vars': {'a': 100, 'b': 32.1, 'c' : 0.12},
                       'qub': {'0': 2}}
 
-        # construc the model, verify the off-diagonal warning is thrown
-        # test that it gives a warning when a key has no corresponding control channel
-        with warnings.catch_warnings(record=True) as w:
-            # Cause all warnings to always be triggered.
-            warnings.simplefilter("always")
-
-            ham_model = HamiltonianModel.from_dict(simple_ham)
-
-            self.assertEqual(len(w), 1)
-            self.assertTrue('minimum overlap' in str(w[-1].message))
+        ham_model = HamiltonianModel.from_dict(simple_ham)
 
         # check norm
         for estate in ham_model._estates: