Repository: qutip/qutip-qip (language: python)

Issue:
Add more code examples in the docstring



Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index 606c4e0e..cb6b8e59 100755
--- a/README.md
+++ b/README.md
@@ -36,7 +36,7 @@ Documentation and tutorials
 The documentation of `qutip-qip` updated to the latest development version is hosted at [qutip-qip.readthedocs.io/](https://qutip-qip.readthedocs.io/en/stable/).
 Tutorials related to using quantum gates and circuits in `qutip-qip` can be found [*here*](https://qutip.org/tutorials#quantum-information-processing) and those related to using noise simulators areavailable at [*this link*](https://qutip.org/tutorials#nisq). 
 
-Code examples used in the preprint [*Pulse-level noisy quantum circuits with QuTiP*](https://arxiv.org/abs/2105.09902), updated for the latest code version, are hosted in [this folder](https://github.com/qutip/qutip-qip/tree/master/doc/pulse-paper).
+Code examples used in the article [*Pulse-level noisy quantum circuits with QuTiP*](https://quantum-journal.org/papers/q-2022-01-24-630), updated for the latest code version, are hosted in [this folder](https://github.com/qutip/qutip-qip/tree/master/doc/pulse-paper).
 
 Installation from source
 ------------------------
@@ -65,22 +65,28 @@ to test the code in the documentation.
 
 Testing
 ------------
-To test the installation from a download of the source code, run from the `qutip-qip` directory
+To test the installation, choose the correct branch that matches with the version, e.g., `qutip-qip-0.2.X` for version 0.2. Then download the source code and run from the `qutip-qip` directory
 ```
 pytest tests
 ```
+
 Citing `qutip-qip`
 ------------
-If you use `qutip-qip` in your research, the  [*preprint*](https://arxiv.org/abs/2105.09902) can be cited as
+If you use `qutip-qip` in your research, please cite the [article](https://quantum-journal.org/papers/q-2022-01-24-630) as
 
 ```bibtex
-@misc{li2021pulselevel,
-      title={Pulse-level noisy quantum circuits with QuTiP},
-      author={Boxi Li and Shahnawaz Ahmed and Sidhant Saraogi and Neill Lambert and Franco Nori and Alexander Pitchford and Nathan Shammah},
-      year={2021},
-      eprint={2105.09902},
-      archivePrefix={arXiv},
-      primaryClass={quant-ph}
+@article{Li2022pulselevelnoisy,
+  doi = {10.22331/q-2022-01-24-630},
+  url = {https://doi.org/10.22331/q-2022-01-24-630},
+  title = {Pulse-level noisy quantum circuits with {Q}u{T}i{P}},
+  author = {Li, Boxi and Ahmed, Shahnawaz and Saraogi, Sidhant and Lambert, Neill and Nori, Franco and Pitchford, Alexander and Shammah, Nathan},
+  journal = {{Quantum}},
+  issn = {2521-327X},
+  publisher = {{Verein zur F{\"{o}}rderung des Open Access Publizierens in den Quantenwissenschaften}},
+  volume = {6},
+  pages = {630},
+  month = jan,
+  year = {2022}
 }
 ```
 Support
diff --git a/doc/pulse-paper/README.md b/doc/pulse-paper/README.md
index 6245e0b4..d3fc4cdd 100644
--- a/doc/pulse-paper/README.md
+++ b/doc/pulse-paper/README.md
@@ -1,4 +1,4 @@
-This folder contains code examples used in the publication [*Pulse-level noisy quantum circuits with QuTiP*](https://arxiv.org/abs/2105.09902). To run the examples, please first install the software package qutip-qip
+This folder contains code examples used in the publication [*Pulse-level noisy quantum circuits with QuTiP*](https://quantum-journal.org/papers/q-2022-01-24-630). To run the examples, please first install the software package qutip-qip
 ```
 pip install qutip_qip[full] joblib
 ```
diff --git a/doc/pulse-paper/customize.py b/doc/pulse-paper/customize.py
index dce51c41..40ee51a6 100644
--- a/doc/pulse-paper/customize.py
+++ b/doc/pulse-paper/customize.py
@@ -38,7 +38,7 @@ def get_control(self, label):
         if label[:2] == "sx":
             return 2 * np.pi * sigmax() / 2, [targets]
         elif label[:2] == "sy":
-            return 2 * np.pi * sigmax() / 2, [targets]
+            return 2 * np.pi * sigmay() / 2, [targets]
         else:
             raise NotImplementError("Unknown control.")
 
@@ -181,7 +181,7 @@ def single_crosstalk_simulation(num_gates):
     """ A single simulation, with num_gates representing the number of rotations.
 
     Args:
-        num_gates (int): The random number of gates to add in the simulation.
+        num_gates (int): The number of random gates to add in the simulation.
 
     Returns:
         result (qutip.solver.Result): A qutip Result object obtained from any of the
diff --git a/doc/pulse-paper/dj_algorithm.py b/doc/pulse-paper/dj_algorithm.py
index 17810a68..5cee511f 100644
--- a/doc/pulse-paper/dj_algorithm.py
+++ b/doc/pulse-paper/dj_algorithm.py
@@ -134,13 +134,15 @@
 fig3.show()
 
 
-plt.rcParams.update({"text.usetex": False})
-
+plt.rcParams.update({"text.usetex": False, "font.size": 9})
 # Plot hinton
 from qutip import hinton
 fig4, ax4 = plt.subplots(figsize=(LINEWIDTH*0.9, LINEWIDTH*0.7), dpi=200)
 first_two_qubits = result1.states[-1].ptrace([0,1])
-_, ax4 = hinton(first_two_qubits, ax=ax4)
+_, ax4 = hinton(
+    first_two_qubits, ax=ax4,
+    xlabels=['$|00\\rangle$', '$|01\\rangle$', '$|10\\rangle$', '$|11\\rangle$'],
+    ylabels=['$|00\\rangle$', '$|01\\rangle$', '$|10\\rangle$', '$|11\\rangle$'])
 fig4.savefig("hinton.pdf")
 fig4.show()
 
@@ -150,7 +152,7 @@
     tmp = state.ptrace([0,1])
     tmp = basis([2,2], [0,0]).dag() * tmp * basis([2,2], [0,0])
     expect.append(np.real(tmp[0, 0]))
-# fig5, ax5 = plt.subplots(figsize=(LINEWIDTH*0.7, LINEWIDTH*0.7*0.6), dpi=200)
+
 fig5, ax5 = plt.subplots(figsize=(LINEWIDTH, LINEWIDTH*0.7), dpi=200)
 ax5.plot(t_record, expect, color="slategrey")
 ax5.set_ylabel(r"Population of $|00\rangle$")
diff --git a/doc/pulse-paper/qft.py b/doc/pulse-paper/qft.py
index e7ea0fa5..49e50e0d 100644
--- a/doc/pulse-paper/qft.py
+++ b/doc/pulse-paper/qft.py
@@ -97,12 +97,13 @@ def get_control_latex(model):
  11.539598500728607,
  31.34040207862854]
 
-fig, ax = plt.subplots(figsize=(LINEWIDTH, LINEWIDTH/3*2), dpi=200)
-ax.plot(range(1, 11), compiling_time, "-s", markersize=4, label="solver")
-ax.plot(range(1, 11), simulation_time, "-D", markersize=4, label="compiling")
+fig, ax = plt.subplots(figsize=(TEXTWIDTH, LINEWIDTH/3*2), dpi=200)
+ax.plot(range(1, 11), compiling_time, "-s", markersize=4, label=r"Compiler (\texttt{Processor.load\_circuit)}")
+ax.plot(range(1, 11), simulation_time, "-D", markersize=4, label=r"Solver (\texttt{Processor.run\_state})")
 ax.set_ylabel("Simulation time [s]")
 ax.set_xlabel("Number of qubits")
 ax.set_yscale("log")
+ax.set_ylim((ax.get_ylim()[0], 150))
 ax.legend()
 fig.tight_layout()
 fig.savefig("runtime.pdf")
diff --git a/doc/requirements.txt b/doc/requirements.txt
index 6603c3c6..098c1bb7 100644
--- a/doc/requirements.txt
+++ b/doc/requirements.txt
@@ -1,10 +1,11 @@
-numpy==1.19
-scipy==1.5.3
-qutip==4.6.0
-cython==0.29.21
-sphinx==3.5.2
+numpy==1.21
+scipy==1.5
+qutip==4.6
+cython==0.29
+sphinx==3.5
 sphinx_rtd_theme==0.5.1
 readthedocs-sphinx-search==0.1.0rc3
-numpydoc==1.1.0
-matplotlib==3.3.4
+numpydoc==1.1
+matplotlib==3.3
 docutils==0.16
+sphinxcontrib-bibtex==2.4
diff --git a/doc/source/apidoc.rst b/doc/source/apidoc.rst
index 695f6639..9c9de1cd 100644
--- a/doc/source/apidoc.rst
+++ b/doc/source/apidoc.rst
@@ -6,7 +6,7 @@ qutip\_qip package
    :maxdepth: 1
 
 
-Operator-level simulation
+Gate-level simulation
 -------------------------
 Simulation based on operator-state multiplication.
 
diff --git a/doc/source/bibliography.rst b/doc/source/bibliography.rst
new file mode 100644
index 00000000..86024a37
--- /dev/null
+++ b/doc/source/bibliography.rst
@@ -0,0 +1,7 @@
+************
+Bibliography
+************
+
+.. bibliography::
+    :cited:
+    :style: unsrt
diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst
index f280f2ef..89d33d59 100644
--- a/doc/source/changelog.rst
+++ b/doc/source/changelog.rst
@@ -1,8 +1,22 @@
-.. _changelog:
-
-**********
+*********
 Changelog
-**********
+*********
+
+Version 0.2.1 (Feb 20, 2021)
+++++++++++++++++++++++++++++
+
+This release includes a revision of the documentation and adds more code examples in the API doc strings.
+
+Bug Fixes
+---------
+- Remove the check on the initial state dimension in ``QubitCircuit.run()`` (#127 <https://github.com/qutip/qutip-qip/pull/127>`_)
+
+Improvements
+------------
+-  Rewrite the documentation for the pulse-level simulation. (#121 <https://github.com/qutip/qutip-qip/pull/121>`_)
+-  Add more code examples in the doc strings. (#126 <https://github.com/qutip/qutip-qip/pull/126>`_)
+
+
 Version 0.2.0 (Nov 26, 2021)
 ++++++++++++++++++++++++++++
 This release adds a few new features to the pulse-level simulator.
diff --git a/doc/source/conf.py b/doc/source/conf.py
index 65aaa551..a28cf9a5 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -57,6 +57,8 @@ def qutip_qip_version():
               'sphinx.ext.viewcode',
               'sphinx.ext.ifconfig',
               'sphinx.ext.napoleon',
+              'sphinxcontrib.bibtex',
+              'sphinx.ext.intersphinx',
 ]
 
 # Add any paths that contain templates here, relative to this directory.
@@ -171,4 +173,15 @@ def qutip_qip_version():
 # attention that some api docs files are adjusted manually for better illustration
 # and should not be overwritten.
 # autosummary_generate = True
-autosummary_imported_members = True
\ No newline at end of file
+autosummary_imported_members = True
+
+# -- Options for biblatex ---------------------------------------
+
+bibtex_bibfiles = ['references.bib']
+bibtex_default_style = 'unsrt'
+
+# -- Options for intersphinx ---------------------------------------
+
+intersphinx_mapping = {
+    'qutip': ('https://qutip.org/docs/latest/', None),
+}
\ No newline at end of file
diff --git a/doc/source/figures/illustration.jpg b/doc/source/figures/illustration.jpg
new file mode 100644
index 00000000..659d6a5a
Binary files /dev/null and b/doc/source/figures/illustration.jpg differ
diff --git a/doc/source/figures/illustration.png b/doc/source/figures/illustration.png
deleted file mode 100644
index c92fc7c9..00000000
Binary files a/doc/source/figures/illustration.png and /dev/null differ
diff --git a/doc/source/figures/structure.jpg b/doc/source/figures/structure.jpg
new file mode 100644
index 00000000..677750c7
Binary files /dev/null and b/doc/source/figures/structure.jpg differ
diff --git a/doc/source/figures/workflow.png b/doc/source/figures/workflow.png
deleted file mode 100644
index 0f43f637..00000000
Binary files a/doc/source/figures/workflow.png and /dev/null differ
diff --git a/doc/source/index.rst b/doc/source/index.rst
index f285691a..f8176fdb 100644
--- a/doc/source/index.rst
+++ b/doc/source/index.rst
@@ -4,9 +4,9 @@
     :align: center
     :width: 5in
 
-***********************************************
-qutip-qip: QuTiP quantum information processing
-***********************************************
+*********************************************************
+qutip-qip |version|: QuTiP quantum information processing
+*********************************************************
 
 .. toctree::
     :maxdepth: 2
@@ -49,6 +49,11 @@ qutip-qip: QuTiP quantum information processing
 
     apidoc.rst
 
+.. toctree::
+    :caption: Bibliography
+
+    bibliography.rst
+
 Indices and tables
 ==================
 
diff --git a/doc/source/introduction.rst b/doc/source/introduction.rst
index ff6e4ecf..6abf53c8 100644
--- a/doc/source/introduction.rst
+++ b/doc/source/introduction.rst
@@ -19,23 +19,28 @@ It offers two different approaches for simulating quantum circuits, one with :cl
 Citing
 ===========
 
-If you use ``qutip.qip`` in your research, cite the `preprint <https://arxiv.org/abs/2105.09902>`_
+If you use ``qutip-qip`` in your research, please cite the `article <https://quantum-journal.org/papers/q-2022-01-24-630>`_
 as
 
 .. code-block:: text
 
-  Boxi Li, Shahnawaz Ahmed, Sidhant Saraogi, Neill Lambert, Franco Nori, Alexander Pitchford, & Nathan Shammah. (2021). Pulse-level noisy quantum circuits with QuTiP.
-
+    Boxi Li, Shahnawaz Ahmed, Sidhant Saraogi, Neill Lambert, Franco Nori, Alexander Pitchford, & Nathan Shammah. (2022). Pulse-level noisy quantum circuits with QuTiP. Quantum, 6, 630.
 
 The bibtex can be downloaded directly :download:`here<qutip_qip.bib>` or
 copy-pasted using :
 
 .. code-block:: text
 
-  @misc{li2021pulselevel,
-        title={Pulse-level noisy quantum circuits with QuTiP},
-        author={Boxi Li and Shahnawaz Ahmed and Sidhant Saraogi and Neill Lambert and Franco Nori and Alexander Pitchford and Nathan Shammah},
-        year={2021},
-        eprint={2105.09902},
-        archivePrefix={arXiv},
-        primaryClass={quant-ph}}
+    @article{Li2022pulselevelnoisy,
+        doi = {10.22331/q-2022-01-24-630},
+        url = {https://doi.org/10.22331/q-2022-01-24-630},
+        title = {Pulse-level noisy quantum circuits with {Q}u{T}i{P}},
+        author = {Li, Boxi and Ahmed, Shahnawaz and Saraogi, Sidhant and Lambert, Neill and Nori, Franco and Pitchford, Alexander and Shammah, Nathan},
+        journal = {{Quantum}},
+        issn = {2521-327X},
+        publisher = {{Verein zur F{\"{o}}rderung des Open Access Publizierens in den Quantenwissenschaften}},
+        volume = {6},
+        pages = {630},
+        month = jan,
+        year = {2022}
+    }
diff --git a/doc/source/qip-processor.rst b/doc/source/qip-processor.rst
index 2a4ccfb5..8c374568 100644
--- a/doc/source/qip-processor.rst
+++ b/doc/source/qip-processor.rst
@@ -4,51 +4,122 @@
 Pulse-level circuit simulation
 ******************************
 
-Modelling quantum hardware with Processor
-=========================================
+This documentation page is adapted from the publication :cite:`Li2022pulselevelnoisy` under `CC BY 4.0 <https://creativecommons.org/licenses/by/4.0/>`_.
 
-Based on the open system solver, :class:`.Processor` in the :mod:`qutip_qip` module simulates quantum circuits at the level of time evolution. One can consider the processor as a simulator of a quantum device, on which the quantum circuit is to be implemented.
+Overview
+========
+
+Based on the open system solver, :class:`.Processor` in the :mod:`qutip_qip` module simulates quantum circuits at the level of time evolution. One can consider the processor as an emulator of a quantum device, on which the quantum circuit is to be implemented.
 
 The procedure is illustrated in the figure below.
-It first compiles circuit into a Hamiltonian model, adds noisy dynamics and then uses the QuTiP open time evolution solvers to simulation the evolution.
+It first compiles the circuit into a Hamiltonian model, adds noisy dynamics and then uses the QuTiP open time evolution solvers to simulate the evolution.
 
-.. image:: /figures//illustration.png
+.. figure:: /figures//illustration.jpg
+    :scale: 40%
+    :align: center
 
-Like a real quantum device, the processor is determined by a list of Hamiltonians, i.e. the control pulses driving the evolution. Given the intensity of the control pulses and the corresponding time slices for each pulse, the evolution is then computed. A control pulse is characterized by :class:`.pulse.Pulse`, consisting of the control Hamiltonian, the targets qubit, the pulse coefficients and the time sequence. We can either use the coefficients as a step function or with cubic spline. For step function, ``tlist`` specifies the start and the end of each pulse and thus is one element longer than the ``coeffs``. One example of defining the control pulse coefficients and the time array is as follows:
+    Illustration of the workflow :cite:`Li2022pulselevelnoisy`.
 
+In the following, we illustrate our framework with an example simulating a 3-qubit Deutsch-Jozsa algorithm on a chain of spin qubits.
+We will work through this example and explain briefly the workflow and all the main modules.
 .. testcode::
 
     import numpy as np
-    from qutip import sigmaz
-    from qutip_qip.device import Processor
+    from qutip import basis
+    from qutip_qip.circuit import QubitCircuit
+    from qutip_qip.device import LinearSpinChain
 
-    processor = Processor(2)
-    processor.add_control(sigmaz(), cyclic_permutation=True)  # sigmaz for all qubits
-    processor.set_coeffs = np.array([[1.0, 1.5, 2.0], [1.8, 1.3, 0.8]])
-    processor.set_tlist = np.array([0.1, 0.2, 0.4, 0.5])
+    # Define a circuit
+    qc = QubitCircuit(3)
+    qc.add_gate("X", targets=2)
+    qc.add_gate("SNOT", targets=0)
+    qc.add_gate("SNOT", targets=1)
+    qc.add_gate("SNOT", targets=2)
 
-It defines a :math:`\sigma_z` operator on both qubits and a pulse that acts on the first qubit.
-An equivalent approach is using the :meth:`.Processor.add_pulse` method.
+    # Oracle function f(x)
+    qc.add_gate("CNOT", controls=0, targets=2)
+    qc.add_gate("CNOT", controls=1, targets=2)
 
-.. testcode::
+    qc.add_gate("SNOT", targets=0)
+    qc.add_gate("SNOT", targets=1)
 
-    from qutip_qip.pulse import Pulse
+    # Run gate-level simulation
+    init_state = basis([2,2,2], [0,0,0])
+    ideal_result = qc.run(init_state)
 
-    processor = Processor(2)
-    coeff=np.array([0.1, 0.2, 0.4, 0.5])
-    tlist=np.array([[1.0, 1.5, 2.0], [1.8, 1.3, 0.8]])
-    pulse = Pulse(sigmaz(), targets=0, coeff=coeff, tlist=tlist)
-    processor.add_pulse(pulse)
+    # Run pulse-level simulation
+    processor = LinearSpinChain(num_qubits=3, sx=0.25, t2=30)
+    processor.load_circuit(qc)
+    tlist = np.linspace(0, 20, 300)
+    result = processor.run_state(init_state, tlist=tlist)
 
-One can also use choose the ``pulse_mode`` attribute of :class:`.Processor`
-between ``"discrete"`` and ``"continuous"``.
+In the above example, we first define a Deutsch-Jozsa circuit, and then run the simulation first at the gate-level.
+We then choose the spin chain model for the underlying physical system, which is a subclass of :obj:`.Processor`.
+We provide the number of qubits and the :math:`\sigma_x` drive strength 0.25MHz.
+The other parameters, such as the interaction strength, are set to be the default value.
+The decoherence noise can also be added by specifying the coherence times (:math:`T_1` and :math:`T_2`) which we discuss hereafter.
 
-.. note::
+By initializing this processor with the hardware parameters, a Hamiltonian model for a spin chain system is generated, including the drift and control Hamiltonians.
+The Hamiltonian model is represented by the :obj:`.Model` class and is saved as an attribute of the initialized processor. In addition, the :obj:`.Processor` can also hold simulation configurations such as whether to use a cubic spline interpolation for the pulse coefficients. Such configurations are not directly part of the model but nevertheless could be important for the pulse-level simulation.
 
-   If the coefficients represent discrete pulse, the length of each array is 1 element shorter than ``tlist``. If it is supposed to be a continuous function, the length should be the same as ``tlist``.
+Next, we provide the circuit to the processor through the method :obj:`.LinearSpinChain.load_circuit`.
+The processor will first decompose the gates in the circuit into native gates that can be implemented directly on the specified hardware model.
+Each gate in the circuit is then mapped to the control coefficients and driving Hamiltonians according to the :obj:`.GateCompiler` defined for a specific model.
+A :obj:`.Scheduler` is used to explore the possibility of executing several pulses in parallel.
+
+With a pulse-level description of the circuit generated and saved in the processor, we can now run the simulation by
+The :obj:`.Processor.run_state` method first builds a Lindblad model including all the defined noise models (none in this example, but options are discussed below) and then calls a QuTiP solver to simulate the time evolution.
+One can pass solver parameters as keyword arguments to the method, e.g., ``tlist`` (time sequence for intermediate results), ``e_ops`` (measurement observables) and ``options`` (solver options).
+In the example above, we record the intermediate state at the time steps given by ``tlist``.
+The returned result is a :class:`~qutip.solver.Result` object, which, depending on the solver options, contains the final state, intermediate states and the expectation value.
+This allows one to extract all information that the solvers in QuTiP provide.
+
+Theory
+======
+
+Down to the physical level, quantum hardware, on which a circuit is executed, is described by quantum theory.
+The dynamics of the system that realizes a unitary gate in quantum circuit is characterized by the time evolution of the quantum system.
+For isolated or open quantum systems, we consider both unitary time evolution and open quantum dynamics.
+The latter can be simulated either by solving the master equation or sampling Monte Carlo trajectories.
+Here, we briefly describe those methods as well as the corresponding solvers available in QuTiP.
+
+For a closed quantum system, the dynamics is determined by the Hamiltonian and the initial state.
+From the perspective of controlling a quantum system, the Hamiltonian is divided into the non-controllable drift :math:`H_{\rm{d}}` (which may be time dependent) and controllable terms combined as :math:`H_{\rm{c}}` to give the full system Hamiltonian
+
+.. _the control Hamiltonian:
+
+.. math::
+
+    H(t) = H_{\rm{d}}(t) + H_{\rm{c}}(t) = H_{\rm{d}}(t) + \sum_j c_j(t) H_j,
+
+where the :math:`H_j` describe the effects of available physical controls on the system that can be modulated by the time-dependent control coefficients :math:`c_j(t)`, by which one drives the system to realize the desired unitary gates.
+
+The unitary :math:`U` that is applied to the quantum system driven by the Hamiltonian :math:`H(t)` is a solution to the Schrödinger operator equation
+
+.. math::
+
+    i \hbar \frac{\partial U(t)}{\partial t}
+    = H(t) U(t)
+    .
+
+By choosing :math:`H(t)` that implements the desired unitaries (the quantum circuit) we obtain a pulse-level description of the circuit in the form of `the control Hamiltonian`_. 
+The choice of the solver depends on the parametrization of the control coefficients :math:`c_j(t)`. 
+The parameters of :math:`c_j(t)` may be determined through theoretical models or automated through control optimisation, as introduced later.
+
+
+Processor and Model
+===================
+
+The :obj:`.Processor` class handles the routine of a pulse-level simulation.
+It connects different modules and works as the main API interface, as the figure below illustrates.
 
+.. figure:: /figures//structure.jpg
+    :scale: 40%
+    :align: center
 
-The above example shows the framework and the most essential part of the simulator's API. So far, it looks like just a wrapper for the open system solvers. However, based on this, we can implement different physical realizations. They differ mainly in how to find the control pulse for a quantum circuit, which gives birth to different sub-classes:
+    Code structure :cite:`Li2022pulselevelnoisy`.
+
+We provide a few predefined processors with Hamiltonian models and compilers routines. They differ mainly in how to find the control pulse for a quantum circuit, which gives birth to different sub-classes:
 
 | Processor
 | ├── ModelProcessor
@@ -57,8 +128,7 @@ The above example shows the framework and the most essential part of the simulat
 | │   └── SpinChain
 | └── OptPulseProcessor
 
-In general, there are two ways to find the control pulses. The first one, :class:`.ModelProcessor`, is more experiment-oriented and based on physical models. A universal set of
-gates is defined in the processor as well as the pulse implementing them in this particular physical model. This is usually the case where control pulses realizing those gates are well known and can be concatenated to realize the whole quantum circuits. Two realizations have already been implemented: the spin chain and the Cavity QED model for quantum computing. In those models, the driving Hamiltonians are predefined. Another approach, based on the optimal control module in QuTiP (http://qutip.org/docs/latest/guide/guide-control.html), is called :class:`.OptPulseProcessor`. In this subclass, one only defines the available Hamiltonians in their system. The processor then uses algorithms to find the optimal control pulses that realize the desired unitary evolution.
+In general, there are two ways to find the control pulses. The first one, :class:`.ModelProcessor`, is more experiment-oriented and based on physical models. An initialized processor has a :obj:`Processor.model` attributes that save the control Hamiltonians. This is usually the case where control pulses realising those gates are well known and can be concatenated to realize the whole quantum circuits. Three realizations have already been implemented: the spin chain, the Cavity QED and the circuit QED model. In those models, the driving Hamiltonians are predefined, saved in the corresponding :obj:`.Model` object. Another approach, based on the `optimal control module in QuTiP <http://qutip.org/docs/latest/guide/guide-control.html>`_, is called :class:`.OptPulseProcessor`. In this subclass, one only defines the available Hamiltonians in their system. The processor then uses algorithms to find the optimal control pulses that realize the desired unitary evolution.
 
 Despite this difference, the logic behind all processors is the same:
 
@@ -66,74 +136,68 @@ Despite this difference, the logic behind all processors is the same:
 
 * The control pulse coefficients and time slices are either specified by the user or calculated by the method :meth:`.Processor.load_circuit`, which takes a :class:`.QubitCircuit` and find the control pulse for this evolution.
 
-* The processor calculates the evolution using the QuTiP solvers. Collapse operators can be added to simulate decoherence. The method :meth:`.Processor.run_state` returns a object :class:`qutip.solver.Result`.
+* The processor calculates the evolution using the QuTiP solvers. Collapse operators can be added to simulate decoherence. The method :meth:`.Processor.run_state` returns a object :class:`~qutip.solver.Result`.
 
 It is also possible to calculate the evolution analytically with matrix exponentiation by setting ``analytical=True``. A list of the matrices representing the gates is returned just like for :meth:`.QubitCircuit.propagators`. However, this does not consider the collapse operators or other noise. As the system size gets larger, this approach will become very inefficient.
 
 In the following, we describe the predefined subclasses for :class:`.Processor`:
 
-**SpinChain**
+SpinChain
+---------
 
 :class:`.LinearSpinChain` and :class:`.CircularSpinChain` are quantum computing models based on the spin chain realization. The control Hamiltonians are :math:`\sigma_x`, :math:`\sigma_z` and :math:`\sigma_x \sigma_x + \sigma_y \sigma_y`. This processor will first decompose the gate into the universal gate set with ISWAP or SQRTISWAP as two-qubit gates, resolve them into quantum gates of adjacent qubits and then calculate the pulse coefficients.
 
-An example of simulating a simple circuit is shown below:
+In the following example we plot the compiled Deutsche Jozsa algorithm:
 
-.. testcode::
+.. plot::
+    :context: reset
 
-    from qutip import basis
+    # Deutsch-Jozsa algorithm
     from qutip_qip.circuit import QubitCircuit
-    from qutip_qip.device import LinearSpinChain
+    qc = QubitCircuit(3)
+    qc.add_gate("X", targets=2)
+    qc.add_gate("SNOT", targets=0)
+    qc.add_gate("SNOT", targets=1)
+    qc.add_gate("SNOT", targets=2)
 
-    qc = QubitCircuit(2)
-    qc.add_gate("X", targets=0)
-    qc.add_gate("X", targets=1)
-    processor = LinearSpinChain(2)
-    processor.load_circuit(qc)
-    result = processor.run_state(basis([2,2], [0,0]))
-    print(result.states[-1].tidyup(1.0e-6))
+    # Oracle function f(x)
+    qc.add_gate("CNOT", controls=0, targets=2)
+    qc.add_gate("CNOT", controls=1, targets=2)
 
-.. testoutput::
-    :options: +NORMALIZE_WHITESPACE
+    qc.add_gate("SNOT", targets=0)
+    qc.add_gate("SNOT", targets=1)
 
-    Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket
-    Qobj data =
-    [[ 0.]
-    [ 0.]
-    [ 0.]
-    [-1.]]
-
-We can also visualize the pulses implementing this circuit:
-
-.. plot::
-
-    from qutip import basis
-    from qutip_qip.circuit import QubitCircuit
     from qutip_qip.device import LinearSpinChain
-
-    qc = QubitCircuit(2)
-    qc.add_gate("X", targets=0)
-    qc.add_gate("X", targets=1)
-    processor = LinearSpinChain(2)
-    processor.load_circuit(qc)
-    fig, axis = processor.plot_pulses()
+    spinchain_processor = LinearSpinChain(num_qubits=3, t2=30)  # T2 = 30
+    spinchain_processor.load_circuit(qc)
+    fig, ax = spinchain_processor.plot_pulses(figsize=(8, 5))
     fig.show()
 
-**Superconducting qubits**
+Superconducting qubits
+----------------------
+
+.. _superconducting qubits:
 
 For the :class:`.SCQubits` model, the qubit is simulated by a three-level system, where the qubit subspace is defined as the ground state and the first excited state.
 The three-level representation will capture the leakage of the population out of the qubit subspace during single-qubit gates.
 The single-qubit control is generated by two orthogonal quadratures :math:`a + a^{\dagger}` and :math:`i(a - a^{\dagger})`, truncated to a three-level operator.
 Same as the Spin Chain model, the superconducting qubits are aligned in a 1 D structure and the interaction is only possible between adjacent qubits.
-As an example, the default interaction is implemented as a Cross Resonant pulses.
-Parameters for the interaction strength are taken from [1]_ [2]_.
+As an example, the default interaction is implemented as a Cross Resonant pulse.
+Parameters for the interaction strength are taken from :cite:`magesan2020effective`:cite:`blais2021circuit`.
 
-.. [1] Easwar Magesan and Jay M. Gambetta.Effective hamiltonian models of the cross-resonance gate. *Phys. Rev. A*, 101:052308, 2020.
+.. plot::
+    :context: close-figs
 
-.. [2] Blais A, Grimsmo A L, Girvin S M, et al. Circuit quantum electrodynamics[J]. *arXiv preprint arXiv:2005.12667*, 2020.
+    from qutip_qip.device import SCQubits
+    scqubits_processor = SCQubits(num_qubits=3)
+    scqubits_processor.load_circuit(qc)
+    fig, ax = scqubits_processor.plot_pulses(figsize=(8, 5))
+    fig.show()
 
-**DispersiveCavityQED**
+DispersiveCavityQED
+-------------------
 
-Same as above, :class:`.DispersiveCavityQED` is a simulator based on Cavity Quantum Electrodynamics. The workflow is similar to the one for the spin chain, except that the component systems are a multi-level cavity and a qubits system. The control Hamiltonians are the single-qubit rotation together with the qubits-cavity interaction :math:`a^{\dagger} \sigma^{-} + a \sigma^{+}`. The device parameters including the cavity frequency, qubits frequency, detuning and interaction strength etc.
+Same as above, :class:`.DispersiveCavityQED` is a simulator based on Cavity Quantum Electrodynamics. The workflow is similar to the one for the spin chain, except that the component systems are a multi-level cavity and a qubits system. The control Hamiltonians are the single-qubit rotation together with the qubits-cavity interaction :math:`a^{\dagger} \sigma^{-} + a \sigma^{+}`. The device parameters include the cavity frequency, qubits frequency, detuning and interaction strength etc.
 
 .. note::
 
@@ -142,7 +206,8 @@ Same as above, :class:`.DispersiveCavityQED` is a simulator based on Cavity Quan
    hence including the cavity.
    To obtain the circuit result, one needs to first trace out the cavity state.
 
-**OptPulseProcessor**
+OptPulseProcessor
+-----------------
 
 The :class:`.OptPulseProcessor` uses the function in :func:`~qutip.control.pulseoptim.optimize_pulse_unitary` in the optimal control module to find the control pulses. The Hamiltonian includes a drift part and a control part and only the control part will be optimized. The unitary evolution follows
 
@@ -152,41 +217,20 @@ The :class:`.OptPulseProcessor` uses the function in :func:`~qutip.control.pulse
 
 To let it find the optimal pulses, we need to give the parameters for :func:`~qutip.control.pulseoptim.optimize_pulse_unitary` as keyword arguments to :meth:`.OptPulseProcessor.load_circuit`. Usually, the minimal requirements are the evolution time ``evo_time`` and the number of time slices ``num_tslots`` for each gate. Other parameters can also be given in the keyword arguments. For available choices, see :func:`~qutip.control.pulseoptim.optimize_pulse_unitary`. It is also possible to specify different parameters for different gates, as shown in the following example:
 
-.. testcode::
-
-      from qutip_qip.device import OptPulseProcessor
-      from qutip import sigmaz, sigmax, sigmay, tensor
-
-
-      # Same parameter for all the gates
-      qc = QubitCircuit(1)
-      qc.add_gate("SNOT", 0)
-
-      num_tslots = 10
-      evo_time = 10
-      processor = OptPulseProcessor(1, drift=sigmaz())
-      processor.add_control(sigmax())
-      # num_tslots and evo_time are two keyword arguments
-      tlist, coeffs = processor.load_circuit(
-      qc, num_tslots=num_tslots, evo_time=evo_time)
-
-      # Different parameters for different gates
-      qc = QubitCircuit(2)
-      qc.add_gate("SNOT", 0)
-      qc.add_gate("SWAP", targets=[0, 1])
-      qc.add_gate('CNOT', controls=1, targets=[0])
-
-      processor = OptPulseProcessor(N=2, drift=tensor([sigmaz()]*2))
-      processor.add_control(sigmax(), cyclic_permutation=True)
-      processor.add_control(sigmay(), cyclic_permutation=True)
-      processor.add_control(tensor([sigmay(), sigmay()]))
-
-      setting_args = {"SNOT": {"num_tslots": 10, "evo_time": 1},
-                      "SWAP": {"num_tslots": 30, "evo_time": 3},
-                      "CNOT": {"num_tslots": 30, "evo_time": 3}}
-
-      tlist, coeffs = processor.load_circuit(
-                      qc, setting_args=setting_args, merge_gates=False)
+.. plot::
+    :context: close-figs
+
+    from qutip_qip.device import OptPulseProcessor, SpinChainModel
+    setting_args = {"SNOT": {"num_tslots": 6, "evo_time": 2},
+                    "X": {"num_tslots": 1, "evo_time": 0.5},
+                    "CNOT": {"num_tslots": 12, "evo_time": 5}}
+    opt_processor = OptPulseProcessor(
+        num_qubits=3, model=SpinChainModel(3, setup="linear"))
+    opt_processor.load_circuit(  # Provide parameters for the algorithm
+        qc, setting_args=setting_args, merge_gates=False,
+        verbose=True, amp_ubound=5, amp_lbound=0)
+    fig, ax = opt_processor.plot_pulses(figsize=(8, 5))
+    fig.show()
 
 Compiler and scheduler
 ======================
@@ -194,103 +238,43 @@ Compiler and scheduler
 Compiler
 --------
 
-In order to simulate quantum circuits at the level of time evolution.
-We need to first compile the circuit into the Hamiltonian model, i.e.
-the control pulses.
-Hence each :class:`.Processor` has a corresponding
-:class:`.compiler.GateCompiler` class.
-The compiler takes a :class:`.QubitCircuit`
-and returns the compiled ``tlist`` and ``coeffs``.
-It is called implicitly when calling the method
-:class:`.Processor.run_state`.
-
-.. testcode::
-
-    from qutip_qip.compiler import SpinChainCompiler
-    qc = QubitCircuit(2)
-    qc.add_gate("X", targets=0)
-    qc.add_gate("X", targets=1)
-
-    processor = LinearSpinChain(2)
-    compiler = SpinChainCompiler(2, params=processor.params)
-    resolved_qc = qc.resolve_gates(["RX", "RZ", "ISWAP"])
-    tlists, coeffs = compiler.compile(resolved_qc)
-    print(tlists)
-    print(coeffs)
+A compiler converts the quantum circuit to the corresponding pulse-level controls :math:`c_j(t)H_j` on the quantum hardware.
+In the framework, it is defined as an instance of the :obj:`.GateCompiler` class.
+The compilation procedure is achieved through the following steps.
 
-**Output**
+First, each quantum gate is decomposed into the native gates (e.g., rotation over :math:`x`, :math:`y` axes and the CNOT gate), using the existing decomposition scheme in QuTiP.
+If a gate acts on two qubits that are not physically connected, like in the chain model and superconducting qubit model, SWAP gates are added to match the topology before the decomposition. Currently, only 1-dimensional chain structures are supported.
 
-.. testoutput::
-    :options: +NORMALIZE_WHITESPACE
+Next, the compiler maps each quantum gate to a pulse-level control description.
+It takes the hardware parameter defined in the Hamiltonian model and computes the pulse duration and strength to implement the gate.
+A pulse scheduler is then used to explore the possibility of executing multiple quantum gates in parallel.
 
-    {'sx0': array([0., 1.]), 'sx1': array([0., 1., 2.])}   
-    {'sx0': array([0.25]), 'sx1': array([0.  , 0.25])} 
+In the end, the compiler returns a time-dependent pulse coefficient :math:`c_j(t)` for each control Hamiltonian :math:`H_j`.
+They contain the full information to implement the circuit and are saved in the processor.
+The coefficient :math:`c_j(t)` is represented by two NumPy arrays, one for the control amplitude and the other for the time sequence.
+For a continuous pulse, a cubic spline is used to approximate the coefficient.
+This allows the use of compiled Cython code in QuTiP to achieve better performance.
 
-Here we first use :meth:`.QubitCircuit.resolve_gates`
-to decompose the X gate to its natural gate on Spin Chain model,
-the rotation over X-axis.
-We pass the hardware parameters of the :class:`.SpinChain` model, ``processor.params``, as well as a map between the pulse name and pulse index ``pulse_dict`` to the compiler.
-The latter one allows one to address the pulse more conveniently in the compiler.
+For the predefined physical models described in the previous subsection, the corresponding compilers are also included and they will be used when calling the method :obj:`.Processor.load_circuit`.
 
-The compiler returns a list of ``tlist`` and ``coeff``, corresponding to each pulse.
-The first pulse starts from ``t=0`` and ends at ``t=1``, with the strengh :math:`\pi/2`.
-The second one is turned on from ``t=1`` to ``t=2`` with the same strength.
-The compiled pulse here is different from what is shown in the plot
-in the previous subsection because the scheduler is turned off by default.
+.. Note::
+    It is expected that the output of the compiler will change after the official release of qutip-v5.
+    The use of the new ``Coefficient`` class will allow more flexibility and improve performance.
+    Therefore, we recommend to not parse the compiled result themselves but use :obj:`.Processor.set_tlist` and :obj:`.Processor.set_coeffs` instead.
 
 Scheduler
 ---------
 
-The scheduler is implemented in the class :class:`.compiler.Scheduler`,
-based on the idea of https://doi.org/10.1117/12.666419.
-It schedules the order of quantum gates and instructions for the
-shortest execution time.
-It works not only for quantum gates but also for pulse implementation of gates
-(:class:`.compiler.Instruction`) with varying pulse duration.
-
-The scheduler first generates a quantum gates dependency graph,
-containing information about which gates have to be executed before some other gates.
-The graph preserves the mobility of the gates,
-i.e. commuting gates are not dependent on each other, even if they use the same qubits.
-Next, it computes the longest distance of each node to the start and end nodes.
-The distance for each dependency arrow is defined by the execution time of the instruction
-(By default, it is 1 for all gates).
-This is used as a priority measure in the next step.
-The gate with a longer distance to the end node and a shorter distance to the start node has higher priority.
-In the last step, it uses a list-schedule algorithm with hardware constraint and
-priority and returns a list of cycles for gates/instructions.
-Since the algorithm is heuristics, sometimes it does not find the optimal solution.
-Hence, we offer an option that randomly shuffles the commuting gates and
-repeats the scheduling a few times to get a better result.
-
-.. testcode::
-
-    from qutip_qip.circuit import QubitCircuit
-    from qutip_qip.compiler import Scheduler
-    circuit = QubitCircuit(7)
-    circuit.add_gate("SNOT", 3)  # gate0
-    circuit.add_gate("CZ", 5, 3)  # gate1
-    circuit.add_gate("CZ", 4, 3)  # gate2
-    circuit.add_gate("CZ", 2, 3)  # gate3
-    circuit.add_gate("CZ", 6, 5)  # gate4
-    circuit.add_gate("CZ", 2, 6)  # gate5
-    circuit.add_gate("ISWAP", [0, 2])  # gate6
-    scheduler = Scheduler("ASAP")
-    result = scheduler.schedule(circuit, gates_schedule=True)
-    print(result)
-
-**Output**
-
-.. testoutput::
-
-    [0, 1, 3, 2, 2, 3, 4]
+The scheduling of a circuit consists of an important part of the compilation.
+Without it, the gates will be executed one by one and many qubits will be idling during the circuit execution, which increases the execution time and reduces the fidelity.
+In the framework, the scheduler is used after the control coefficient of each gate is computed.
+It runs a scheduling algorithm to determine the starting time of each gate while keeping the result correct.
 
-The result shows the scheduling order of each gate in the original circuit.
+The heuristic scheduling algorithm we provide offers two different modes: ASAP (as soon as possible) and ALAP (as late as possible).
+In addition, one can choose whether permutation among commuting gates is allowed to achieve a shorter execution time.
+The scheduler implemented here does not take the hardware architecture into consideration and assumes that the connectivity in the provided circuit matches with the hardware at this step.
 
-For pulse schedule, or scheduling gates with different duration,
-one will need to wrap the :class:`.Gate` object with :class:`.compiler.instruction` object,
-with a parameter `duration`.
-The result will then be the start time of each instruction.
+In predefined processors, the scheduler runs automatically when loading a circuit and hence there is no action necessary from the side of the user.
 
 Pulse shape
 -----------
@@ -310,73 +294,42 @@ The method :obj:`.GateCompiler.generate_pulse_shape` allows one to generate puls
     ax.set_xlabel("Time")
     fig.show()
 
-For predefined compilers, the compiled pulse shape can also be configured by the key word ``"shape"`` and ``"num_samples"`` in the dictionary attribute :attr:`.GateCompiler.args`
+For predefined compilers, the compiled pulse shape can also be configured by the keyword ``"shape"`` and ``"num_samples"`` in the dictionary attribute :attr:`.GateCompiler.args`
 or the ``args`` parameter of :obj:`.GateCompiler.compile`.
 
 Noise Simulation
 ================
 
-In the common way of QIP simulation, where evolution is carried out by gate matrix product, the noise is usually simulated with bit flipping and sign flipping errors.
-The typical approaches are either applying bit/sign flipping gate probabilistically
-or applying Kraus operators representing different noisy channels (e.g. amplitude damping, dephasing) after each unitary gate evolution. In the case of a single qubit, they have the same effect and the parameters in the Kraus operators are exactly the probability of a flipping error happens during the gate operation time.
+.. _noise:
 
-Since the processor simulates the state evolution at the level of the driving Hamiltonian, there is no way to apply an error operator to the continuous-time evolution. Instead, the error is added to the pulses (coherent control error) or the collapse operators (Lindblad error) contributing to the evolution. Mathematically, this is no different from adding an error channel probabilistically (it is actually how :func:`qutip.mcsolve` works internally). The collapse operator for single-qubit amplitude damping and dephasing are exactly the destroying operator and the sign-flipping operator. One just needs to choose the correct coefficients for them to simulate the noise, e.g. the relaxation time T1 and dephasing time T2. Because it is based on the open system evolution instead of abstract operators, this simulation is closer to the physical implementation and requires less pre-analysis of the system.
+The :obj:`Noise` module allows one to add control and decoherence noise following the Lindblad description of open quantum systems.
+Compared to the gate-based simulator, this provides a more practical and straightforward way to describe the noise.
+In the current framework, noise can be added at different layers of the simulation, allowing one to focus on the dynamics of the dominant noise, while representing other noise, such as single-qubit relaxation, as collapse operators for efficiency.
+Depending on the problem studied, one can devote the computing resources to the most relevant type of noise.
 
-Compared to the approach of Kraus operators, this way of simulating noise is more computationally expensive. If you only want to simulate the decoherence of single-qubit relaxation and the relaxation time is much longer than the gate duration, there is no need to go through all the calculations. However, this simulator is closer to the real experiment and, therefore, more convenient in some cases, such as when coherent noise or correlated noise exist. For instance, a pulse on one qubit might affect the neighbouring qubits, the evolution is still unitary but the gate fidelity will decrease. It is not always easy or even possible to define a noisy gate matrix. In our simulator, it can be done by defining a :class:`.noise.ControlAmpNoise` (Control Amplitude Noise).
+Apart from imperfections in the Hamiltonian model and circuit compilation, the :obj:`Noise` class in the current framework defines deviations of the real physical dynamics from the compiled one.
+It takes the compiled pulse-level description of the circuit (see also `Pulse representation`_) and adds noise elements to it, which allows defining noise that is correlated to the compiled pulses.
+In the following, we detail the three different noise models already available in the current framework.
 
-In the simulation, noise can be added to the processor at different levels:
+Noise in the hardware model
+---------------------------
 
-- The decoherence time T1 and T2 can be defined for the processor or for each qubit. When calculating the evolution, the corresponding collapse operators will be added automatically to the solver.
+The Hamiltonian model defined in the :obj:`.Model` class may contain intrinsic imperfections of the system and hence the compiled ideal pulse does not implement the ideal unitary gate.
+Therefore, building a realistic Hamiltonian model usually already introduces noise to the simulation.
+An example is the `superconducting qubits`_, where the physical qubit is represented by a multi-level system.
+Since the second excitation level is only weakly detuned from the qubit transition frequency, the population may leak out of the qubit subspace.
+Another example is an always-on ZZ type cross-talk induced by interaction with higher levels of the physical qubits :cite:`mundada2019suppression`, which is also implemented for the superconducting qubit model.
 
-- The noise of the physical parameters (e.g. detuned frequency) can be simulated by changing the parameters in the model, e.g. laser frequency in cavity QED. (This can only be time-independent since QuTiP open system solver only allows varying coefficients, not varying Hamiltonian operators.)
+Control noise
+-------------
 
-- The noise of the pulse intensity can be simulated by modifying the coefficients of the Hamiltonian operators or even adding new Hamiltonians.
+The control noise, as the name suggests, arises from imperfect control of the quantum system, such as distortion in the pulse amplitude or frequency drift.
+The simplest example is the random amplitude noise on the control coefficient :math:`c_j(t)` in `the control Hamiltonian`_.
 
-To add noise to a processor, one needs to first define a noise object :class:`.noise.Noise`. The simplest relaxation noise can be defined directly in the processor with relaxation time. Other pre-defined noise can be found as subclasses of  :class:`.noise.Noise`. We can add noise to the simulator with the method :meth:`.Processor.add_noise`.
-
-Below, we show two examples.
-
-The first example is a processor with one qubit under rotation around the z-axis and relaxation time :math:`T_2=5`. We measure the population of the :math:`\left| + \right\rangle` state and observe the Ramsey signal:
+The following example demonstrates a biased Gaussian noise on the pulse amplitude. For visualization purposes, we plot the noisy pulse intensity instead of the state fidelity. The three pulses can, for example, be a zyz-decomposition of an arbitrary single-qubit gate:
 
 .. plot::
 
-    import numpy as np
-    import matplotlib.pyplot as plt
-    from qutip import sigmaz, destroy, basis
-    from qutip_qip.device import Processor
-    from qutip_qip.operations import snot
-
-    a = destroy(2)
-    Hadamard = snot()
-    plus_state = (basis(2,1) + basis(2,0)).unit()
-    tlist = np.arange(0.00, 20.2, 0.2)
-
-    T2 = 5
-    processor = Processor(1, t2=T2)
-    processor.add_control(sigmaz(), label="sz")
-    processor.set_coeffs({"sz": np.ones(len(tlist))})
-    processor.set_tlist(tlist)
-    result = processor.run_state(
-        plus_state, e_ops=[a.dag()*a, Hadamard*a.dag()*a*Hadamard])
-
-    fig, ax = plt.subplots()
-    # detail about length of tlist needs to be fixed
-    ax.plot(tlist[:-1], result.expect[1][:-1], '.', label="simulation")
-    ax.plot(tlist[:-1], np.exp(-1./T2*tlist[:-1])*0.5 + 0.5, label="theory")
-    ax.set_xlabel("t")
-    ax.set_ylabel("Ramsey signal")
-    ax.legend()
-    ax.set_title("Relaxation T2=5")
-    ax.grid()
-    fig.tight_layout()
-    fig.show()
-
-The second example demonstrates a biased Gaussian noise on the pulse amplitude. For visualization purposes, we plot the noisy pulse intensity instead of the state fidelity. The three pulses can, for example, be a zyz-decomposition of an arbitrary single-qubit gate:
-
-.. plot::
-
-    import numpy as np
-    import matplotlib.pyplot as plt
     from qutip import sigmaz, sigmay
     from qutip_qip.device import Processor
     from qutip_qip.noise import RandomNoise
@@ -407,30 +360,94 @@ The second example demonstrates a biased Gaussian noise on the pulse amplitude.
         use_control_latex=False)
     axis2[0].step(qobjevo.tlist, noisy_coeff)
 
+Lindblad noise
+--------------
+
+The Lindblad noise originates from the coupling of the quantum system with the environment (e.g., a thermal bath) and leads to loss of information.
+It is simulated by collapse operators and results in non-unitary dynamics :cite:`breuer2002theory`:cite:`lidar2019lecture`.
+
+The most commonly used type of Lindblad noise is decoherence, characterized by the coherence time :math:`T_1` and :math:`T_2` (dephasing).
+For the sake of convenience, one only needs to provide the parameter ``t1``, ``t2`` to the processor and the corresponding operators will be generated automatically.
+Both can be either a number that specifies one coherence time for all qubits or a list of numbers, each corresponding to one qubit.
+
+For :math:`T_1`, the operator is defined as :math:`a/\sqrt{T_1}` with :math:`a` the destruction operator.
+For :math:`T_2`, the operator is defined as :math:`a^{\dagger}a\sqrt{2/T_2^*}`, where :math:`T_2^*` is the pure dephasing time given by :math:`1/T_2^*=1/T_2-1/(2T_1)`.
+In the case of qubits, i.e., a two-level system, the destruction operator :math:`a` is truncated to a two-level operator and is consistent with the Lindblad equation.
+Constant :math:`T_1` and :math:`T_2` can be provided directly when initializing the :obj:`.Processor`.
+Custom collapse operators, including time-dependent ones, can be defined through :obj:`.DecoherenceNoise`.
+For instance, the following code defines a collapse operator using :obj:`~qutip.sigmam()` and increases linearly as time:
+
+.. testcode::
+
+    from qutip_qip.device import LinearSpinChain
+    from qutip_qip.noise import DecoherenceNoise
+    tlist = np.linspace(0, 30., 100)
+    coeff = tlist * 0.01
+    noise = DecoherenceNoise(
+        sigmam(), targets=0,
+        coeff=coeff, tlist=tlist)
+    processor = LinearSpinChain(1)
+    processor.add_noise(noise)
+
+Similar to the control noise, the Lindblad noise can also depend on the control coefficient.
+
+Pulse representation
+====================
+
+.. _Pulse representation:
+
+As discussed before, in this simulation framework, we compile the circuit into pulse-level controls and simulate the time evolution of the physical qubits.
+In this subsection, we describe how the dynamics is represented internally in the workflow, which is useful for understanding the simulation process as well as defining custom pulse-dependent noise.
+
+A control pulse, together with the noise associated with it, is represented by a class instance of :obj:`.Pulse`.
+When an ideal control is compiled and returned to the processor, it is saved as an initialized :obj:`.Pulse`, equivalent to the following code:
+
+.. testcode:: [pulse]
+
+    from qutip_qip.pulse import Pulse
+    coeff = np.array([1.])
+    tlist = np.array([0., np.pi])
+    pulse = Pulse(
+        sigmax()/2, targets=0, tlist=tlist,
+        coeff=coeff, label="pi-pulse")
+
+This code defines a :math:`\pi`-pulse implemented using the term :math:`\sigma_x` in the Hamiltonian that flips the zeroth qubit specified by the argument ``targets``. The pulse needs to be applied for the duration :math:`\pi` specified by the variable ``tlist``. The parameters ``coeff`` and ``tlist`` together describe the control coefficient.
+Together with the provided Hamiltonian and target qubits, an instance of :obj:`.Pulse` determines the dynamics of one control term.
+
+.. note::
+
+   If the coefficients represent discrete pulse, the length of each array is 1 element shorter than ``tlist``. If it is supposed to be a continuous function, the length should be the same as ``tlist``.
+
+With a :obj:`.Pulse` initialized with the ideal control, one can define several types of noise, including the Lindblad or control noise as described in the `noise`_ section.
+An example of adding a noisy Hamiltonian as control noise through the method :obj:`.Pulse.add_control_noise` is given below:
+
+.. testcode:: [pulse]
+
+    pulse.add_control_noise(
+        sigmaz(), targets=[0], tlist=tlist,
+        coeff=coeff * 0.05)
+
+The above code snippet adds a Hamiltonian term :math:`\sigma_z`, which can, for instance, be interpreted as a frequency drift.
+Similarly, collapse operators depending on a specific control pulse can be added by the method :obj:`.Pulse.add_lindblad_noise`.
+
+In addition to a constant pulse, the control pulse and noise can also be provided as continuous functions. In this case, both ``tlist`` and ``coeff`` are given as NumPy arrays and a cubic spline is used to interpolate the continuous pulse coefficient.
+This allows using the compiled Cython version of the QuTiP solvers that have a much better performance than using a Python function for the coefficient.
+The option is provided as a keyword argument ``spline_kind="cubic"`` when initializing :obj:`.Pulse`.
+Similarly, the interpolation method can also be defined for :obj:`.Processor` using the same signature.
+
 
 Customize the simulator
 =======================
 
-The number of predefined physical models and compilers are limited.
-However, it is designed for easy customization and one can easily build customized model and compiling routines.
-For guide and examples, please refer to the tutorial notebooks
-at http://qutip.org/tutorials.html
-
-The workflow of the simulator
-=============================
-
-The following plot demonstrates the workflow of the simulator.
-
-.. image:: /figures//workflow.png
-
-The core of the simulator is :class:`.Processor`,
-which characterizes the quantum hardware of interest,
-containing the information such as the non-controllable drift Hamiltonian and
-the control Hamiltonian.
-Apart from the ideal system representing the qubits, one can also define
-hardware-dependent or pulse-dependent noise in :class:`.noise.Noise`.
-It describes how noisy terms such as imperfect control
-and decoherence can be added once the ideal control pulse is defined.
-When loading a quantum circuit, a :class:`.compiler.GateCompiler` compiles the circuit into a sequence of control pulse signals and schedule the pulse for parallel execution.
-For each control Hamiltonian, a :class:`.pulse.Pulse` instance is created that including the ideal evolution and associated noise.
-They will then be sent to the QuTiP solvers for the computation.
+As it is impractical to include every physical platform, we provide an interface that allows one to customize the simulators.
+In particular, the modular architecture allows one to conveniently overwrite existing modules for customization.
+
+To define a customized hardware model, the minimal requirements are a set of available control Hamiltonians :math:`H_j`, and a compiler, i.e., the mapping between native gates and control coefficients :math:`c_j`.
+One can either modify an existing subclass or write one from scratch by creating a subclass of the two parent classes :obj:`.Model` and :obj:`.GateCompiler`.
+Since different subclasses share the same interface, different models and compilers can also be combined to build new processors.
+
+Moreover, this customization is not limited to Hamiltonian models and compiler routines.
+In principle, measurement can be defined as a customized quantum gate and the measurement statistics can be extracted from the obtained density matrix.
+A new type of noise can also be implemented by defining a new :obj:`.Noise` subclass, which takes the compiled ideal :obj:`.Pulse` and adds noisy dynamics on top of it.
+
+For examples, please refer to the `tutorial notebooks <http://qutip.org/tutorials.html>`_.
diff --git a/doc/source/qip-simulator.rst b/doc/source/qip-simulator.rst
index d3de4939..08db9cd3 100644
--- a/doc/source/qip-simulator.rst
+++ b/doc/source/qip-simulator.rst
@@ -1,8 +1,8 @@
 .. _qip_simulator:
 
-*********************************
-Operator-level circuit simulation
-*********************************
+*****************************
+Gate-level circuit simulation
+*****************************
 
 Run a quantum circuit
 =====================
diff --git a/doc/source/qutip_qip.bib b/doc/source/qutip_qip.bib
index f2439a86..227204d6 100644
--- a/doc/source/qutip_qip.bib
+++ b/doc/source/qutip_qip.bib
@@ -1,7 +1,13 @@
-  @misc{li2021pulselevel,
-      title={Pulse-level noisy quantum circuits with QuTiP},
-      author={Boxi Li and Shahnawaz Ahmed and Sidhant Saraogi and Neill Lambert and Franco Nori and Alexander Pitchford and Nathan Shammah},
-      year={2021},
-      eprint={2105.09902},
-      archivePrefix={arXiv},
-      primaryClass={quant-ph}}
+@article{Li2022pulselevelnoisy,
+  doi = {10.22331/q-2022-01-24-630},
+  url = {https://doi.org/10.22331/q-2022-01-24-630},
+  title = {Pulse-level noisy quantum circuits with {Q}u{T}i{P}},
+  author = {Li, Boxi and Ahmed, Shahnawaz and Saraogi, Sidhant and Lambert, Neill and Nori, Franco and Pitchford, Alexander and Shammah, Nathan},
+  journal = {{Quantum}},
+  issn = {2521-327X},
+  publisher = {{Verein zur F{\"{o}}rderung des Open Access Publizierens in den Quantenwissenschaften}},
+  volume = {6},
+  pages = {630},
+  month = jan,
+  year = {2022}
+}
\ No newline at end of file
diff --git a/doc/source/references.bib b/doc/source/references.bib
new file mode 100644
index 00000000..f74e93e1
--- /dev/null
+++ b/doc/source/references.bib
@@ -0,0 +1,65 @@
+@article{Li2022pulselevelnoisy,
+  doi = {10.22331/q-2022-01-24-630},
+  url = {https://doi.org/10.22331/q-2022-01-24-630},
+  title = {Pulse-level noisy quantum circuits with {Q}u{T}i{P}},
+  author = {Li, Boxi and Ahmed, Shahnawaz and Saraogi, Sidhant and Lambert, Neill and Nori, Franco and Pitchford, Alexander and Shammah, Nathan},
+  journal = {{Quantum}},
+  issn = {2521-327X},
+  publisher = {{Verein zur F{\"{o}}rderung des Open Access Publizierens in den Quantenwissenschaften}},
+  volume = {6},
+  pages = {630},
+  month = jan,
+  year = {2022}
+}
+@article{magesan2020effective,
+  title={Effective Hamiltonian models of the cross-resonance gate},
+  author={Magesan, Easwar and Gambetta, Jay M},
+  journal={Physical Review A},
+  volume={101},
+  number={5},
+  pages={052308},
+  year={2020},
+  publisher={APS}
+}
+@article{blais2021circuit,
+  title={Circuit quantum electrodynamics},
+  author={Blais, Alexandre and Grimsmo, Arne L and Girvin, SM and Wallraff, Andreas},
+  journal={Reviews of Modern Physics},
+  volume={93},
+  number={2},
+  pages={025005},
+  year={2021},
+  publisher={APS}
+}
+@article{mundada2019suppression,
+  title={Suppression of qubit crosstalk in a tunable coupling superconducting circuit},
+  author={Mundada, Pranav and Zhang, Gengyan and Hazard, Thomas and Houck, Andrew},
+  journal={Physical Review Applied},
+  volume={12},
+  number={5},
+  pages={054023},
+  year={2019},
+  publisher={APS}
+}
+@book{breuer2002theory,
+  title={The theory of open quantum systems},
+  author={Breuer, Heinz-Peter and Petruccione, Francesco and others},
+  year={2002},
+  publisher={Oxford University Press on Demand}
+}
+@article{lidar2019lecture,
+  title={Lecture notes on the theory of open quantum systems},
+  author={Lidar, Daniel A},
+  journal={arXiv preprint arXiv:1902.00967},
+  year={2019}
+}
+@article{koch2007charge,
+  title={Charge-insensitive qubit design derived from the Cooper pair box},
+  author={Koch, Jens and Terri, M Yu and Gambetta, Jay and Houck, Andrew A and Schuster, David Isaac and Majer, Johannes and Blais, Alexandre and Devoret, Michel H and Girvin, Steven M and Schoelkopf, Robert J},
+  journal={Physical Review A},
+  volume={76},
+  number={4},
+  pages={042319},
+  year={2007},
+  publisher={APS}
+}
\ No newline at end of file
diff --git a/setup.cfg b/setup.cfg
index 40be6f8d..b5c5ac69 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -28,7 +28,7 @@ package_dir=
 packages = find:
 include_package_data = True
 install_requires =
-    numpy>=1.16.6
+    numpy>=1.16.6,<=1.21
     scipy>=1.0
     qutip>=4.6,<5
 setup_requires =
diff --git a/src/qutip_qip/algorithms/qft.py b/src/qutip_qip/algorithms/qft.py
index 493c1e70..48516c01 100644
--- a/src/qutip_qip/algorithms/qft.py
+++ b/src/qutip_qip/algorithms/qft.py
@@ -31,7 +31,7 @@ def qft(N=1):
     if N < 1:
         raise ValueError("Minimum value of N can be 1")
 
-    N2 = 2 ** N
+    N2 = 2**N
     phase = 2.0j * np.pi / N2
     arr = np.arange(N2)
     L, M = np.meshgrid(arr, arr)
diff --git a/src/qutip_qip/circuit.py b/src/qutip_qip/circuit.py
index 7c91431c..d06f7328 100644
--- a/src/qutip_qip/circuit.py
+++ b/src/qutip_qip/circuit.py
@@ -2403,8 +2403,6 @@ def initialize(self, state=None, cbits=None, measure_results=None):
         self.state = None
 
         if state is not None:
-            if state.shape[0] != 2 ** self.qc.N:
-                raise ValueError("dimension of state is incorrect")
             if self.mode == "density_matrix_simulator" and state.isket:
                 self.state = ket2dm(state)
             else:
diff --git a/src/qutip_qip/circuit_latex.py b/src/qutip_qip/circuit_latex.py
index 4c357972..d4f977cc 100644
--- a/src/qutip_qip/circuit_latex.py
+++ b/src/qutip_qip/circuit_latex.py
@@ -88,7 +88,6 @@ def _crop_pdf(filename):
         _force_remove(filename)
         os.rename(temporary, filename)
 
-
 else:
 
     def _crop_pdf(_):
@@ -248,7 +247,6 @@ def image_from_latex(code, file_type="png"):
                 os.chdir(previous_dir)
         return out
 
-
 else:
 
     def image_from_latex(*args, **kwargs):
diff --git a/src/qutip_qip/compiler/cavityqedcompiler.py b/src/qutip_qip/compiler/cavityqedcompiler.py
index a94742f2..7605bf10 100644
--- a/src/qutip_qip/compiler/cavityqedcompiler.py
+++ b/src/qutip_qip/compiler/cavityqedcompiler.py
@@ -36,7 +36,7 @@ class CavityQEDCompiler(GateCompiler):
 
     params: dict
         A Python dictionary contains the name and the value of the parameters.
-        See :meth:`.DispersiveCavityQED.set_up_params` for the definition.
+        See :obj:`.CavityQEDModel` for the definition.
 
     global_phase: float, optional
         Record of the global phase change and will be returned.
@@ -46,6 +46,33 @@ class CavityQEDCompiler(GateCompiler):
     gate_compiler: dict
         The Python dictionary in the form of {gate_name: decompose_function}.
         It saves the decomposition scheme for each gate.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from qutip_qip.circuit import QubitCircuit
+    >>> from qutip_qip.device import ModelProcessor, CavityQEDModel
+    >>> from qutip_qip.compiler import CavityQEDCompiler
+    >>>
+    >>> qc = QubitCircuit(2)
+    >>> qc.add_gate("ISWAP", targets=[0, 1])
+    >>>
+    >>> model = CavityQEDModel(2)
+    >>> processor = ModelProcessor(model=model)
+    >>> compiler = CavityQEDCompiler(2, params=model.params)
+    >>> processor.load_circuit(
+    ...     qc, compiler=compiler)  # doctest: +NORMALIZE_WHITESPACE
+    ({'sz0': array([   0.        , 2500.        , 2500.01315789]),
+    'sz1': array([   0.        , 2500.        , 2500.01315789]),
+    'g0': array([   0., 2500.]),
+    'g1': array([   0., 2500.])},
+    {'sz0': array([-0.5, -9.5]),
+    'sz1': array([-0.5, -9.5]),
+    'g0': array([0.01]),
+    'g1': array([0.01])})
+
+    Notice that the above example is equivalent to using directly
+    the :obj:`.DispersiveCavityQED`.
     """
 
     def __init__(
diff --git a/src/qutip_qip/compiler/circuitqedcompiler.py b/src/qutip_qip/compiler/circuitqedcompiler.py
index 2c64d83b..0c0992c7 100644
--- a/src/qutip_qip/compiler/circuitqedcompiler.py
+++ b/src/qutip_qip/compiler/circuitqedcompiler.py
@@ -26,6 +26,46 @@ class SCQubitsCompiler(GateCompiler):
         +-----------------+-----------------------+
         |``params``       | Hardware Parameters   |
         +-----------------+-----------------------+
+
+    Parameters
+    ----------
+    num_qubits: int
+        The number of qubits in the system.
+
+    params: dict
+        A Python dictionary contains the name and the value of the parameters.
+        See :meth:`.SCQubitsModel` for the definition.
+
+    Attributes
+    ----------
+    num_qubits: int
+        The number of the component systems.
+
+    params: dict
+        A Python dictionary contains the name and the value of the parameters,
+        such as laser frequency, detuning etc.
+
+    gate_compiler: dict
+        The Python dictionary in the form of {gate_name: decompose_function}.
+        It saves the decomposition scheme for each gate.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from qutip_qip.circuit import QubitCircuit
+    >>> from qutip_qip.device import ModelProcessor, SCQubitsModel
+    >>> from qutip_qip.compiler import SCQubitsCompiler
+    >>>
+    >>> qc = QubitCircuit(2)
+    >>> qc.add_gate("CNOT", targets=0, controls=1)
+    >>>
+    >>> model = SCQubitsModel(2)
+    >>> processor = ModelProcessor(model=model)
+    >>> compiler = SCQubitsCompiler(2, params=model.params)
+    >>> processor.load_circuit(qc, compiler=compiler);  # doctest: +SKIP
+
+    Notice that the above example is equivalent to using directly
+    the :obj:`.SCQubits`.
     """
 
     def __init__(self, num_qubits, params):
@@ -56,10 +96,10 @@ def _normalized_gauss_pulse(self):
         td = 2.4384880692912567
         sigma = 1 / 6 * td  # 3 sigma
         tlist = np.linspace(0, td, 1000)
-        max_pulse = 1 - np.exp(-((0 - td / 2) ** 2) / 2 / sigma ** 2)
+        max_pulse = 1 - np.exp(-((0 - td / 2) ** 2) / 2 / sigma**2)
         coeff = (
-            np.exp(-((tlist - td / 2) ** 2) / 2 / sigma ** 2)
-            - np.exp(-((0 - td / 2) ** 2) / 2 / sigma ** 2)
+            np.exp(-((tlist - td / 2) ** 2) / 2 / sigma**2)
+            - np.exp(-((0 - td / 2) ** 2) / 2 / sigma**2)
         ) / max_pulse
         return tlist, coeff
 
diff --git a/src/qutip_qip/compiler/gatecompiler.py b/src/qutip_qip/compiler/gatecompiler.py
index 9cb10f5d..3b4a5c58 100644
--- a/src/qutip_qip/compiler/gatecompiler.py
+++ b/src/qutip_qip/compiler/gatecompiler.py
@@ -349,6 +349,52 @@ def generate_pulse_shape(cls, shape, num_samples, maximum=1.0, area=1.0):
         ``Processor.pulse_mode`` to ``"continuous"``.
         Notice that finite number of sampling points will also make
         the total integral of the coefficients slightly deviate from ``area``.
+
+        Examples
+        --------
+        .. plot::
+            :context: reset
+
+            from qutip_qip.compiler import GateCompiler
+            import numpy as np
+            compiler = GateCompiler()
+            coeff, tlist= compiler.generate_pulse_shape(
+                "hann",  # Scipy Hann window
+                1000,  # 100 sampling point
+                maximum=3.,
+                # Notice that 2 pi is added to H by qutip solvers.
+                area= 1.,
+            )
+
+        We can plot the generated pulse shape:
+
+        .. plot::
+            :context: close-figs
+
+            import matplotlib.pyplot as plt
+            plt.plot(tlist, coeff)
+            plt.show()
+
+        The pulse is normalized to fit the area. Notice that due to
+        the finite number of sampling points, it is not exactly 1.
+
+        .. testsetup::
+
+            from qutip_qip.compiler import GateCompiler
+            import numpy as np
+            compiler = GateCompiler()
+            coeff, tlist= compiler.generate_pulse_shape(
+                "hann",  # Scipy Hann window
+                1000,  # 100 sampling point
+                maximum=3.,
+                # Notice that 2 pi is added to H by qutip solvers.
+                area= 1.,
+            )
+
+        .. doctest::
+
+            >>> round(np.trapz(coeff, tlist), 2)
+            1.0
         """
         coeff, tlist = _normalized_window(shape, num_samples)
         sign = np.sign(area)
diff --git a/src/qutip_qip/compiler/spinchaincompiler.py b/src/qutip_qip/compiler/spinchaincompiler.py
index 8828b569..f97129a2 100644
--- a/src/qutip_qip/compiler/spinchaincompiler.py
+++ b/src/qutip_qip/compiler/spinchaincompiler.py
@@ -35,7 +35,7 @@ class SpinChainCompiler(GateCompiler):
 
     params: dict
         A Python dictionary contains the name and the value of the parameters.
-        See :meth:`.SpinChain.set_up_params` for the definition.
+        See :obj:`.SpinChainModel` for the definition.
 
     setup: string
         "linear" or "circular" for two sub-classes.
@@ -70,6 +70,28 @@ class SpinChainCompiler(GateCompiler):
 
     global_phase: bool
         Record of the global phase change and will be returned.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from qutip_qip.circuit import QubitCircuit
+    >>> from qutip_qip.device import ModelProcessor, SpinChainModel
+    >>> from qutip_qip.compiler import SpinChainCompiler
+    >>>
+    >>> qc = QubitCircuit(2)
+    >>> qc.add_gate("RX", 0, arg_value=np.pi)
+    >>> qc.add_gate("RZ", 1, arg_value=np.pi)
+    >>>
+    >>> model = SpinChainModel(2, "linear", g=0.1)
+    >>> processor = ModelProcessor(model=model)
+    >>> compiler = SpinChainCompiler(2, params=model.params, setup="linear")
+    >>> processor.load_circuit(
+    ...     qc, compiler=compiler) # doctest: +NORMALIZE_WHITESPACE
+    ({'sx0': array([0., 1.]), 'sz1': array([0.  , 0.25])},
+    {'sx0': array([0.25]), 'sz1': array([1.])})
+
+    Notice that the above example is equivalent to using directly
+    the :obj:`.LinearSpinChain`.
     """
 
     def __init__(
diff --git a/src/qutip_qip/device/cavityqed.py b/src/qutip_qip/device/cavityqed.py
index 300bfac1..00084b17 100644
--- a/src/qutip_qip/device/cavityqed.py
+++ b/src/qutip_qip/device/cavityqed.py
@@ -28,7 +28,7 @@
 class DispersiveCavityQED(ModelProcessor):
     """
     The processor based on the physical implementation of
-    a dispersive cavity QED system.
+    a dispersive cavity QED system (:obj:`.CavityQEDModel`).
     The available Hamiltonian of the system is predefined.
     For a given pulse amplitude matrix, the processor can
     calculate the state evolution under the given control pulse,
@@ -51,6 +51,37 @@ class DispersiveCavityQED(ModelProcessor):
 
     **params:
         Hardware parameters. See :obj:`CavityQEDModel`.
+
+    Examples
+    --------
+
+    .. testcode::
+
+        import numpy as np
+        import qutip
+        from qutip_qip.circuit import QubitCircuit
+        from qutip_qip.device import DispersiveCavityQED
+
+        qc = QubitCircuit(2)
+        qc.add_gate("RX", 0, arg_value=np.pi)
+        qc.add_gate("RY", 1, arg_value=np.pi)
+        qc.add_gate("ISWAP", [1, 0])
+
+        processor = DispersiveCavityQED(2, g=0.1)
+        processor.load_circuit(qc)
+        result = processor.run_state(
+            qutip.basis([10, 2, 2], [0, 0, 0]),
+            options=qutip.Options(nsteps=5000))
+        final_qubit_state = result.states[-1].ptrace([1, 2])
+        print(round(qutip.fidelity(
+            final_qubit_state,
+            qc.run(qutip.basis([2, 2], [0, 0]))
+        ), 4))
+
+    .. testoutput::
+
+        0.9994
+
     """
 
     def __init__(
@@ -134,22 +165,59 @@ class CavityQEDModel(Model):
     """
     The physical model for a dispersive cavity-QED processor
     (:obj:`.DispersiveCavityQED`).
+    It is a qubit-resonator model that describes a system composed of
+    a single resonator and a few qubits connected to it.
+    The coupling is kept small so that the resonator is rarely
+    excited but acts only as a mediator for entanglement generation.
+    The single-qubit control Hamiltonians used are
+    :math:`\sigma_x` and :math:`\sigma_z`.
+    The dynamics between the resonator and the qubits is captured by
+    the Tavis-Cummings Hamiltonian,
+    :math:`\propto\sum_j a^\dagger \sigma_j^{-} + a \sigma_j^{+}`,
+    where :math:`a`, :math:`a^\dagger` are
+    the destruction and creation operators of the resonator,
+    while :math:`\sigma_j^{-}`, :math:`\sigma_j^{+}` are those of each qubit.
+    The control of the qubit-resonator coupling depends on
+    the physical implementation, but in the most general case
+    we have single and multi-qubit control in the form
+
+    .. math::
+
+        H=
+        \\sum_{j=0}^{N-1}
+        \\epsilon^{\\rm{max}}_{j}(t) \\sigma^x_{j} +
+        \\Delta^{\\rm{max}}_{j}(t) \\sigma^z_{j} +
+        J_{j}(t)
+        (a^\\dagger \\sigma^{-}_{j} + a \\sigma^{+}_{j}).
+
+    The effective qubit-qubit coupling is computed by the
+
+    .. math::
+
+        J_j = \\frac{g_j g_{j+1}}{2}(\\frac{1}{\\Delta_j} +
+        \\frac{1}{\\Delta_{j+1}}),
+
+    with :math:`\\Delta=w_q-w_0`
+    and the dressed qubit frequency :math:`w_q` defined as
+    :math:`w_q=\sqrt{\epsilon^2+\delta^2}`.
 
     Parameters
     ----------
     num_qubits : int
-        The number of qubits.
+        The number of qubits :math:`N`.
     num_levels : int, optional
         The truncation level of the Hilbert space for the resonator.
     **params :
         Keyword arguments for hardware parameters, in the unit of GHz.
         Qubit parameters can either be a float or a list of the length
-        ``num_qubits``.
+        :math:`N`.
 
         - deltamax: float or list, optional
-            The pulse strength of sigma-x control, default ``1.0``.
+            The pulse strength of sigma-x control,
+            :math:`\\Delta^{\\rm{max}}`, default ``1.0``.
         - epsmax: float or list, optional
-            The pulse strength of sigma-z control, default ``9.5``.
+            The pulse strength of sigma-z control,
+            :math:`\\epsilon^{\\rm{max}}`, default ``9.5``.
         - eps: float or list, optional
             The bare transition frequency for each of the qubits,
             default ``9.5``.
@@ -159,15 +227,13 @@ class CavityQEDModel(Model):
             The coupling strength between the resonator and the qubit,
             default ``1.0``.
         - w0 : float, optional
-            The bare frequency of the resonator. Should only be a float,
-            default ``0.01``.
+            The bare frequency of the resonator :math:`w_0`.
+            Should only be a float, default ``0.01``.
         - t1 : float or list, optional
             Characterize the amplitude damping for each qubit.
         - t2 : list of list, optional
             Characterize the total dephasing for each qubit.
 
-        The dressed qubit frequency is `wq` is computed by
-        :math:`w_q=\sqrt{\epsilon^2+\delta^2}`
     """
 
     def __init__(self, num_qubits, num_levels=10, **params):
diff --git a/src/qutip_qip/device/circuitqed.py b/src/qutip_qip/device/circuitqed.py
index c5968728..638b8d9c 100644
--- a/src/qutip_qip/device/circuitqed.py
+++ b/src/qutip_qip/device/circuitqed.py
@@ -16,7 +16,8 @@
 
 class SCQubits(ModelProcessor):
     """
-    A chain of superconducting qubits with fixed frequency.
+    A chain of superconducting qubits with fixed frequency
+    (:obj:`SCQubitsModel`).
     Single-qubit control is realized by rotation around the X and Y axis
     while two-qubit gates are implemented with Cross Resonance gates.
     A 3-level system is used to simulate the superconducting qubit system,
@@ -25,8 +26,6 @@ class SCQubits(ModelProcessor):
     system, as a demonstration and
     for simplicity, we only use a ZX Hamiltonian for
     the two-qubit interaction.
-    For details see https://arxiv.org/abs/2005.12667 and
-    https://journals.aps.org/pra/abstract/10.1103/PhysRevA.101.052308.
 
     Parameters
     ----------
@@ -39,6 +38,26 @@ class SCQubits(ModelProcessor):
         If ZZ cross-talk is included.
     **params:
         Hardware parameters. See :obj:`SCQubitsModel`.
+
+    Examples
+    --------
+
+    .. testcode::
+
+        import numpy as np
+        import qutip
+        from qutip_qip.circuit import QubitCircuit
+        from qutip_qip.device import SCQubits
+
+        qc = QubitCircuit(2)
+        qc.add_gate("RZ", 0, arg_value=np.pi)
+        qc.add_gate("RY", 1, arg_value=np.pi)
+        qc.add_gate("CNOT", targets=0, controls=1)
+
+        processor = SCQubits(2)
+        processor.load_circuit(qc)
+        init_state = qutip.basis([3, 3], [0, 0])
+        result = processor.run_state(init_state)
     """
 
     def __init__(self, num_qubits, dims=None, zz_crosstalk=False, **params):
@@ -63,6 +82,43 @@ class SCQubitsModel(Model):
     """
     The physical model for superconducting-qubit model processor
     (:obj:`.SCQubits`) with fixed frequency.
+    Each qubit is simulated by a multi-level Duffing model
+    :cite:`koch2007charge`,
+    in which the qubit subspace is provided by the ground state
+    and the first excited state.
+    By default, the creation and annihilation operators are
+    truncated at the third level, which can be adjusted.
+    The multi-level representation can capture the leakage of the population
+    out of the qubit subspace during single-qubit gates.
+    The single-qubit control is generated by two orthogonal quadratures
+    :math:`a_j^{\\dagger}+a_j` and :math:`i(a_j^{\\dagger}-a_j)`.
+    The interaction is possible only between adjacent qubits.
+    Although this interaction is mediated by a resonator, for simplicity,
+    we replace the complicated dynamics among two superconducting qubits
+    and the resonator with a two-qubit effective Hamiltonian
+    derived in Ref. :cite:`magesan2020effective`.
+
+    As an example, we choose the cross resonance interaction
+    in the form of :math:`\\sigma^z_{j} \\sigma^x_{j+1}`,
+    acting only on the two-qubit levels,
+    which is widely used, e.g., in fixed-frequency superconducting qubits.
+    We can write the Hamiltonian as
+
+    .. math::
+
+        H &= H_{\\rm{d}}
+         + \\sum_{j=0}^{N-1} \\Omega^x_{j} (a_j^{\\dagger} + a_j)
+        + \\Omega^y_{j} i(a_j^{\\dagger} - a_j) \\\\
+        &
+        + \\sum_{j=0}^{N-2} \\Omega^{\\rm{cr}1}_{j} \\sigma^z_j \\sigma^x_{j+1}
+        + \\Omega^{\\rm{cr}2}_{j} \\sigma^x_j \\sigma^z_{j+1}
+
+    where the drift Hamiltonian is defined as
+
+    .. math::
+
+        H_{\\rm{d}} = \\sum_{j=0}^{N-1}
+        \\frac{\\alpha_j}{2} a_j^{\\dagger}a_j^{\\dagger}a_j a_j
 
     Parameters
     ----------
@@ -90,7 +146,8 @@ class SCQubitsModel(Model):
             Anharmonicity for each superconducting qubit,
             default ``[-0.3]*num_qubits``.
         - omega_single : list, optional
-            Control strength for single-qubit gate,
+            The maximal control strength for single-qubit gate,
+            :math:`\\Omega^x` and :math:`\\Omega^y`,
             default ``[0.01]*num_qubits``.
         - omega_cr : list, optional
             Control strength for cross resonance gate,
@@ -103,7 +160,7 @@ class SCQubitsModel(Model):
 
     def __init__(self, num_qubits, dims=None, zz_crosstalk=False, **params):
         self.num_qubits = num_qubits
-        self.dims = dims
+        self.dims = dims if dims is not None else [3] * num_qubits
         self.params = {
             "wq": np.array(
                 ((5.15, 5.09) * int(np.ceil(self.num_qubits / 2)))[
@@ -130,7 +187,7 @@ def _set_up_drift(self):
             destroy_op = destroy(self.dims[m])
             coeff = 2 * np.pi * self.params["alpha"][m] / 2.0
             self._drift.append(
-                (coeff * destroy_op.dag() ** 2 * destroy_op ** 2, [m])
+                (coeff * destroy_op.dag() ** 2 * destroy_op**2, [m])
             )
 
     @property
diff --git a/src/qutip_qip/device/processor.py b/src/qutip_qip/device/processor.py
index f8eb10a9..6327da9d 100644
--- a/src/qutip_qip/device/processor.py
+++ b/src/qutip_qip/device/processor.py
@@ -27,9 +27,20 @@
 class Processor(object):
     """
     The noisy quantum device simulator using QuTiP dynamic solvers.
-    It compiles quantums circuit into a Hamiltonian model and then
+    It compiles quantum circuit into a Hamiltonian model and then
     simulate the time-evolution described by the master equation.
 
+    .. note::
+
+        This is an abstract class that includes the general API but
+        has no concrete physical model implemented.
+        In particular, it provides a series of low-level APIs that allow
+        direct modification of the Hamiltonian model and control pulses,
+        which can usually be achieved automatically using :obj:`.Model`
+        and build-in workflows.
+        They provides more flexibility but are not always the most
+        elegant approaches.
+
     Parameters
     ----------
     num_qubits : int, optional
@@ -221,7 +232,8 @@ def add_control(
         self, qobj, targets=None, cyclic_permutation=False, label=None
     ):
         """
-        Add a control Hamiltonian to the model.
+        Add a control Hamiltonian to the model. The new control Hamiltonian
+        is saved in the :obj:`.Processor.model` attributes.
 
         Parameters
         ----------
@@ -243,6 +255,21 @@ def add_control(
             If ``None``,
             it will be set to the current number of
             control Hamiltonians in the system.
+
+        Examples
+        --------
+        >>> import qutip
+        >>> from qutip_qip.device import Processor
+        >>> processor = Processor(1)
+        >>> processor.add_control(qutip.sigmax(), 0, label="sx")
+        >>> processor.get_control_labels()
+        ['sx']
+        >>> processor.get_control("sx") # doctest: +NORMALIZE_WHITESPACE
+        (Quantum object: dims = [[2], [2]], shape = (2, 2),
+        type = oper, isherm = True
+        Qobj data =
+        [[0. 1.]
+        [1. 0.]], [0])
         """
         targets = self._unify_targets(qobj, targets)
         if label is None:
@@ -268,6 +295,19 @@ def get_control(self, label):
         -------
         control_hamiltonian : tuple
             The control Hamiltonian in the form of ``(qobj, targets)``.
+
+        Examples
+        --------
+        >>> from qutip_qip.device import LinearSpinChain
+        >>> processor = LinearSpinChain(1)
+        >>> processor.get_control_labels()
+        ['sx0', 'sz0']
+        >>> processor.get_control('sz0') # doctest: +NORMALIZE_WHITESPACE
+        (Quantum object: dims = [[2], [2]], shape = (2, 2),
+        type = oper, isherm = True
+        Qobj data =
+        [[ 6.28318531  0.        ]
+        [ 0.         -6.28318531]], 0)
         """
         return self.model.get_control(label)
 
@@ -403,7 +443,6 @@ def set_coeffs(self, coeffs):
             )
 
     set_all_coeffs = set_coeffs
-    set_all_coeffs.__doc__ = "Equivalent to :obj:`Processor.set_coeffs`."
 
     def set_tlist(self, tlist):
         """
@@ -430,7 +469,6 @@ def set_tlist(self, tlist):
             self.pulses[pulse_dict[pulse_label]].tlist = value
 
     set_all_tlist = set_tlist
-    set_all_coeffs.__doc__ = "Equivalent to :obj:`Processor.set_tlist`."
 
     def get_full_tlist(self, tol=1.0e-10):
         """
diff --git a/src/qutip_qip/device/spinchain.py b/src/qutip_qip/device/spinchain.py
index 72d6170b..0a4e350e 100644
--- a/src/qutip_qip/device/spinchain.py
+++ b/src/qutip_qip/device/spinchain.py
@@ -102,6 +102,7 @@ def load_circuit(self, qc, setup, schedule_mode="ASAP", compiler=None):
 class LinearSpinChain(SpinChain):
     """
     Spin chain model with open-end topology.
+    For the control Hamiltonian please refer to :obj:`SpinChainModel`.
 
     Parameters
     ----------
@@ -115,6 +116,32 @@ class LinearSpinChain(SpinChain):
 
     **params:
         Hardware parameters. See :obj:`.SpinChainModel`.
+
+    Examples
+    --------
+
+    .. testcode::
+
+        import numpy as np
+        import qutip
+        from qutip_qip.circuit import QubitCircuit
+        from qutip_qip.device import LinearSpinChain
+
+        qc = QubitCircuit(2)
+        qc.add_gate("RX", 0, arg_value=np.pi)
+        qc.add_gate("RY", 1, arg_value=np.pi)
+        qc.add_gate("ISWAP", [1, 0])
+
+        processor = LinearSpinChain(2, g=0.1, t1=300)
+        processor.load_circuit(qc)
+        init_state = qutip.basis([2, 2], [0, 0])
+        result = processor.run_state(init_state)
+        print(round(qutip.fidelity(result.states[-1], qc.run(init_state)), 4))
+
+    .. testoutput::
+
+        0.994
+
     """
 
     def __init__(
@@ -159,6 +186,7 @@ class CircularSpinChain(SpinChain):
     """
     Spin chain model with circular topology. See :class:`.SpinChain`
     for details.
+    For the control Hamiltonian please refer to :obj:`SpinChainModel`.
 
     Parameters
     ----------
@@ -172,6 +200,31 @@ class CircularSpinChain(SpinChain):
 
     **params:
         Hardware parameters. See :obj:`.SpinChainModel`.
+
+    Examples
+    --------
+
+    .. testcode::
+
+        import numpy as np
+        import qutip
+        from qutip_qip.circuit import QubitCircuit
+        from qutip_qip.device import CircularSpinChain
+
+        qc = QubitCircuit(2)
+        qc.add_gate("RX", 0, arg_value=np.pi)
+        qc.add_gate("RY", 1, arg_value=np.pi)
+        qc.add_gate("ISWAP", [1, 0])
+
+        processor = CircularSpinChain(2, g=0.1, t1=300)
+        processor.load_circuit(qc)
+        init_state = qutip.basis([2, 2], [0, 0])
+        result = processor.run_state(init_state)
+        print(round(qutip.fidelity(result.states[-1], qc.run(init_state)), 4))
+
+    .. testoutput::
+
+        0.994
     """
 
     def __init__(
@@ -222,12 +275,28 @@ def topology_map(self, qc):
 class SpinChainModel(Model):
     """
     The physical model for the spin chian processor
-    (:obj:`CircularSpinChain` and :obj`LinearSpinChain`).
+    (:obj:`CircularSpinChain` and :obj:`LinearSpinChain`).
+    The interaction is only possible between adjacent qubits.
+    The single-qubit control Hamiltonians are :math:`\\sigma_j^x`$`,
+    :math:`\\sigma_j^z`, while the interaction is realized by
+    the exchange Hamiltonian
+    :math:`\\sigma^x_{j}\\sigma^x_{j+1}+\\sigma^y_{j}\\sigma^y_{j+1}`.
+    The overall Hamiltonian model is written as:
+
+    .. math::
+
+        H=
+        \\sum_{j=0}^{N-1}
+        \\Omega^x_{j}(t) \\sigma^x_{j} +
+        \\Omega^z_{j}(t) \\sigma^z_{j} + \\sum_{j=0}^{N-2}
+        g_{j}(t)
+        (\\sigma^x_{j}\\sigma^x_{j+1}+
+        \\sigma^y_{j}\\sigma^y_{j+1}).
 
     Parameters
     ----------
     num_qubits: int
-        The number of qubits.
+        The number of qubits, :math:`N`.
     setup : str
         "linear" for an open end and "circular" for a closed end chain.
     **params :
@@ -237,14 +306,16 @@ class SpinChainModel(Model):
         for each qubits.
 
         - sx : float or list, optional
-            The pulse strength of sigma-x control, default ``0.25``.
+            The pulse strength of sigma-x control, :math:`\\Omega^x`,
+            default ``0.25``.
         - sz : float or list, optional
-            The pulse strength of sigma-z control, default ``1.0``.
+            The pulse strength of sigma-z control, :math:`\\Omega^z`,
+            default ``1.0``.
         - sxsy : float or list, optional
-            The pulse strength for the exchange interaction,
+            The pulse strength for the exchange interaction, :math:`g`,
             default ``0.1``.
             It should be either a float or an array of the length
-            ``num_qubits-1`` for the linear setup or ``num_qubits`` for
+            :math:`N-1` for the linear setup or :math:`N` for
             the circular setup.
         - t1 : float or list, optional
             Characterize the amplitude damping for each qubit.
diff --git a/src/qutip_qip/operations/gates.py b/src/qutip_qip/operations/gates.py
index bfc45edb..509423f0 100644
--- a/src/qutip_qip/operations/gates.py
+++ b/src/qutip_qip/operations/gates.py
@@ -1356,7 +1356,7 @@ def globalphase(theta, N=1):
 
     """
     data = np.exp(1.0j * theta) * sp.eye(
-        2 ** N, 2 ** N, dtype=complex, format="csr"
+        2**N, 2**N, dtype=complex, format="csr"
     )
     return Qobj(data, dims=[[2] * N, [2] * N])
 
@@ -1706,7 +1706,7 @@ def qubit_clifford_group(N=None, target=0):
 
     X = sigmax()
     S = phasegate(np.pi / 2)
-    E = H * (S ** 3) * w ** 3
+    E = H * (S**3) * w**3
 
     for op in map(
         partial(reduce, mul),
diff --git a/src/qutip_qip/pulse.py b/src/qutip_qip/pulse.py
index 16e261c5..5ec95959 100644
--- a/src/qutip_qip/pulse.py
+++ b/src/qutip_qip/pulse.py
@@ -142,13 +142,13 @@ class Pulse:
     noise and the lindblad noise. The later two are lists of
     noisy evolution dynamics.
     Each dynamic element is characterized by four variables:
-    ``1`qobj``, ``targets``, ``tlist`` and ``coeff``.
+    ``qobj``, ``targets``, ``tlist`` and ``coeff``.
 
     See examples for different construction behavior.
 
     Parameters
     ----------
-    qobj : :class:'qutip.Qobj'
+    qobj : :class:`qutip.Qobj`
         The Hamiltonian of the ideal pulse.
     targets: list
         target qubits of the ideal pulse
@@ -285,7 +285,7 @@ def add_coherent_noise(self, qobj, targets, tlist=None, coeff=None):
 
         Parameters
         ----------
-        qobj: :class:'qutip.Qobj'
+        qobj: :class:`qutip.Qobj`
             The Hamiltonian of the pulse.
         targets: list
             target qubits of the pulse
@@ -316,7 +316,7 @@ def add_lindblad_noise(self, qobj, targets, tlist=None, coeff=None):
 
         Parameters
         ----------
-        qobj: :class:'qutip.Qobj'
+        qobj: :class:`qutip.Qobj`
             The collapse operator of the lindblad noise.
         targets: list
             target qubits of the collapse operator
@@ -505,7 +505,7 @@ def add_drift(self, qobj, targets):
 
         Parameters
         ----------
-        qobj: :class:'qutip.Qobj'
+        qobj: :class:`qutip.Qobj`
             The collapse operator of the lindblad noise.
         targets: list
             target qubits of the collapse operator
diff --git a/src/qutip_qip/qubits.py b/src/qutip_qip/qubits.py
index aeac776e..591dad9b 100644
--- a/src/qutip_qip/qubits.py
+++ b/src/qutip_qip/qubits.py
@@ -31,7 +31,7 @@ def qubit_states(N=1, states=[0]):
 
     return tensor(
         [
-            alpha * basis(2, 1) + sqrt(1 - alpha ** 2) * basis(2, 0)
+            alpha * basis(2, 1) + sqrt(1 - alpha**2) * basis(2, 0)
             for alpha in state_list
         ]
     )


Test patch:
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 8404390e..c2407d96 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -40,7 +40,7 @@ jobs:
     - name: Set up Python
       uses: actions/setup-python@v2
       with:
-        python-version: 3.6
+        python-version: 3.8
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
@@ -68,11 +68,11 @@ jobs:
     - name: Set up Python
       uses: actions/setup-python@v2
       with:
-        python-version: 3.6
+        python-version: 3.8
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
-        python -m pip install numpy scipy cython qutip sphinx matplotlib numpydoc sphinx_rtd_theme
+        python -mpip install -r doc/requirements.txt
         pip install .
     - name: Test code snippets in the documentation
       run: |
diff --git a/tests/test_circuit.py b/tests/test_circuit.py
index a4010431..f1fb442b 100644
--- a/tests/test_circuit.py
+++ b/tests/test_circuit.py
@@ -480,7 +480,7 @@ def test_N_level_system(self):
         """
         Test for circuit with N-level system.
         """
-        mat3 = rand_dm(3, density=1.)
+        mat3 = qp.rand_unitary_haar(3)
 
         def controlled_mat3(arg_value):
             """
@@ -495,7 +495,13 @@ def controlled_mat3(arg_value):
         qc.user_gates = {"CTRLMAT3": controlled_mat3}
         qc.add_gate("CTRLMAT3", targets=[1, 0], arg_value=1)
         props = qc.propagators()
-        np.testing.assert_allclose(mat3, ptrace(props[0], 0) - 1)
+        final_fid = qp.average_gate_fidelity(mat3, ptrace(props[0], 0) - 1)
+        assert pytest.approx(final_fid, 1.0e-6) == 1
+
+        init_state = basis([3, 2], [0, 1])
+        result = qc.run(init_state)
+        final_fid = qp.fidelity(result, props[0] * init_state)
+        assert pytest.approx(final_fid, 1.0e-6) == 1.
 
     @pytest.mark.repeat(10)
     def test_run_teleportation(self):