Which solver do you call, what do you hand it, and how do you know the answer is right?
Installation first, with the exact error messages you are likely to meet, and a map of where each program’s own documentation lives. Then QuTiP in two topics (objects, then solvers) and scqubits in two (a qubit, then coupled systems and sweeps). The last of them joins the two programs, which is how they are used together. The code is meant to be copied and changed.
sesolve. Losses: mesolve, which is exact. Losses and a space too large for a density matrix: mcsolve, which averages random histories and has a sampling error.Pick an experiment (a Rabi oscillation, a T1 decay, a Ramsey fringe, an echo), set its numbers, and run it. The plot is drawn by the same solver the earlier topics verified against QuTiP; below it is the exact Python that would give you the same curve, with your numbers filled in, and a button to copy it. Turn on “random histories” to see what mcsolve does: each history is one qubit that either jumped or did not, and their average is the smooth curve.
Honest limits: this is one qubit. It is the same equation QuTiP solves, but it will not replace QuTiP for a larger system, and the generated code is a starting point you should read before you trust it.
Every solver is asked the same question: here is a Hamiltonian, here is where the system starts, here are the times I care about; tell me what I would measure. The three you need are sesolve for a closed system, mesolve for a system with losses, and mcsolve, which solves the lossy problem a different way. What you get back is a Result object. If you passed e_ops, a list of operators, it holds their expectation values at each time in result.expect; if you did not, it holds the state itself at every time in result.states. Asking for expectation values is almost always the right choice, because the states are large and the expectation values are the numbers an experiment reports.
Losses. mesolve takes the collapse operators of topic 28 as a list. It evolves the density matrix, which has d² entries for a d-level system, and its answer is exact for the Lindblad equation. mcsolve takes the same list but evolves a single state vector through random histories, in each of which the system either jumps or does not, and averages them. It needs d numbers instead of d², which matters when d is large, but its answer carries a sampling error that shrinks only as one over the square root of the number of histories. With 400 histories the error is about 0.02 to 0.03 in a probability, which is what the example below shows. Pass seeds to make a run repeatable.
A pulse that changes in time. A constant Hamiltonian is not an experiment: pulses start and stop. QuTiP’s list format handles it. You write the Hamiltonian as a list whose entries are either an operator, which is constant, or a pair [operator, function of time]. The solver multiplies each operator by its function at every instant. A Gaussian pulse whose area under the curve equals π is a π pulse, whatever its shape, in the rotating frame: the example shows it flips the qubit completely.
How to know the answer is right. The solvers are good, so errors come from what you gave them. Run a case you can do by hand first (the closed-system line below has the answer sin²(Ωt/2)). Check that a quantity that must be conserved is conserved (the trace of a density matrix, the total probability). Tighten the tolerance with options={"atol": 1e-10, "rtol": 1e-8} and see that the answer does not move. And check the units: times in ns need rates in rad per ns (topic 25). A solver that says it has finished has not said it is right.
| Solver | Starting state | Losses | Returns |
|---|---|---|---|
sesolve(H, psi0, times, e_ops=[...]) | a ket | none | expectation values, or states |
mesolve(H, rho0, times, c_ops, e_ops=[...]) | a ket or a density matrix | list of collapse operators | expectation values, or states |
mcsolve(H, psi0, times, c_ops, e_ops=[...], ntraj=N, seeds=s) | a ket | list of collapse operators | averaged expectation values, with the individual runs kept |
Other solvers in the same family are brmesolve (the Bloch–Redfield equation, for noise with a spectrum), fmmesolve and the Floquet tools (for periodic drives), steadystate (the state a lossy system settles into) and the Krylov and stochastic solvers. The QuTiP 5 paper (Lambert et al., Physics Reports 1153, 2026) describes them, and the user guide chapters “Time Evolution and Quantum System Dynamics” and “Environments of Open Quantum Systems” are where to read next.
The worked example. A 10 MHz Rabi drive with T1 = 120 ns and T2 = 90 ns. Closed, the excited-state population at 0, 50, 100, 150 and 200 ns is 0, 1, 0, 1, 0. With the losses it is 0, 0.789, 0.3038, 0.6021 and 0.4187, a flipping whose peaks sink towards a half as the qubit decoheres. 400 random histories (seed 7) give 0, 0.8, 0.32, 0.61 and 0.4, within 0.02 of the exact answer. A Gaussian pulse of width 15 ns and area π leaves the qubit excited with probability 1.0000.
Try it in Python. Run on QuTiP 5.2.0, 25 October 2026:
import numpy as np
import qutip as qt
TAU = 2 * np.pi
Omega = TAU * 0.01 # a 10 MHz Rabi frequency, in rad per ns
H = 0.5 * Omega * qt.sigmax() # resonant drive, rotating frame
T1, T2 = 120.0, 90.0 # ns
gphi = 1 / T2 - 1 / (2 * T1)
c_ops = [np.sqrt(1 / T1) * qt.destroy(2), np.sqrt(gphi / 2) * qt.sigmaz()]
times = np.linspace(0, 200, 5) # 0, 50, 100, 150, 200 ns
e_ops = [qt.num(2)] # what to record: the excited-state population
closed = qt.sesolve(H, qt.basis(2, 0), times, e_ops=e_ops) # no losses: a pure state
open_ = qt.mesolve(H, qt.fock_dm(2, 0), times, c_ops, e_ops=e_ops) # with losses: a density matrix
print("sesolve:", np.round(closed.expect[0], 4))
print("mesolve:", np.round(open_.expect[0], 4))
# the same losses, played out as single random histories and averaged (Monte Carlo)
mc = qt.mcsolve(H, qt.basis(2, 0), times, c_ops, e_ops=e_ops, ntraj=400, seeds=7, options={"progress_bar": False})
print("mcsolve:", np.round(mc.expect[0], 2), " largest gap to mesolve:", round(float(np.max(np.abs(mc.expect[0] - open_.expect[0]))), 2))
# a pulse that is not constant: a Gaussian whose area is pi. QuTiP's list format: [operator, function of time]
sigma = 15.0 # ns
area_pi = np.pi / (sigma * np.sqrt(TAU)) # peak rate that makes the area under the pulse equal pi
pulse = lambda t: area_pi * np.exp(-((t - 100.0) ** 2) / (2 * sigma ** 2))
H_t = [0.5 * qt.sigmax(), pulse]
res = qt.sesolve(H_t, qt.basis(2, 0), np.linspace(0, 200, 401), e_ops=e_ops)
print("after the pi pulse, P(excited) =", round(float(res.expect[0][-1]), 4))sesolve: [0. 1. 0. 1. 0.]
mesolve: [0. 0.789 0.3038 0.6021 0.4187]
mcsolve: [0. 0.8 0.32 0.61 0.4 ] largest gap to mesolve: 0.02
after the pi pulse, P(excited) = 1.0An echo cancels only what is static: a detuning that is the same during the first and the second half of the sequence. The Lindblad model of this part has noise with no memory, so nothing in it is static, and its echo curve is the plain decay e−t/T2: the envelope of the Ramsey fringes, without the fringes. The explorer for this topic shows it, and the Python it writes reproduces it (the verifier checks the two against each other). An echo that outlives Ramsey in a real device is evidence of slow noise, and simulating that needs a noise model with memory, such as a detuning that is itself a random function of time, not a larger T2.
If your own echo decays faster than e−t/T2, check three things: that the π pulse acts on the density matrix as σxρσx, that it comes exactly half way, and that T2 was entered no larger than 2T1. The same question appears in QuTiP’s Discussions as Hahn Echo in a qubit decoherence simulation giving a T2 echo too short.
And for a pulse that changes in time (Methods for writing time-dependent operators, the other long thread there): use the list format of the example above, [operator, function of time]. A function that also takes args still runs but raises a FutureWarning (topic 33 lists the message); write it as f(t).
Every pulse-level and open-system result in a hardware paper
The curves in a calibration or coherence figure are, in most simulation sections, the output of one of these three calls with a Hamiltonian written as in the earlier topics and collapse operators taken from measured T1 and T2. Pulse shaping (a Gaussian or DRAG envelope, a flat-top for a two-qubit gate) is a list-format Hamiltonian. Qiskit Dynamics does a similar job inside the Qiskit toolkit: its documentation (version 0.6.0 when checked) lists a Solver class, a LindbladModel and a tutorial on the Lindblad dynamics of a qubit chain, with JAX as an optional backend. Qiskit 2.0 removed Qiskit’s own Pulse module, so check its pages before you rely on pulse-schedule features. Choosing between them is mostly a question of which toolkit the rest of your work lives in.
Check yourself
You run mcsolve with 400 histories and mesolve on the same system. How should the two answers compare?
Averaged over enough histories, mcsolve converges to the Lindblad answer, but each run carries noise that falls as 1/√N: about 0.02 to 0.03 in a probability at 400 histories in the topic’s example.
The tempting wrong answers.