QuTiP 2— Rabi Oscillations in Python with QuTiP

Simulating coherent quantum control in a two-level system

Infographic by NotebookLM

In the previous article, QuTiP for Quantum Technology: Getting Started, we built a first working example of quantum simulation in Python. The aim was not to cover every feature of QuTiP but to make the basic workflow concrete: define a quantum system, choose an initial state, evolve it, and visualise the result. In this second article we will use the same workflow to simulate one of the simplest and most important examples of coherent quantum dynamics — the Rabi oscillation.

A Rabi oscillation describes a two-level quantum system being driven back and forth between two states. In physical terms, those two states could represent selected energy levels in an atom, the spin-up and spin-down states of a spin, an ion qubit driven by a laser, or a superconducting qubit controlled by a microwave pulse. This makes the Rabi oscillation a natural second step after the introductory article: it is simple enough to simulate in a few lines of Python, yet important enough to connect directly with real quantum technologies.

To make it easy to follow along, I have prepared a companion Google Colab notebook that contains all of the code in this article. You can open it at: https://colab.research.google.com/drive/13ab5HspJ4_5WXAjeWN4p3ujR8Y7_6IHm and run each cell as you read, adjusting parameters such as the Rabi frequency to see how the dynamics change.

1. What we will build

We will simulate a two-level quantum system whose two states are the ground state |0⟩ and the excited state |1⟩. The system will start in |0⟩, and we will then apply a Hamiltonian that drives transitions between the two levels. As time evolves, the probability of finding the system in each state changes periodically: the population flows from |0⟩ to |1⟩ and then back again.

By the end of the article we will have:

• built a two-level quantum system in QuTiP;
• defined a Hamiltonian that drives coherent oscillations;
• used mesolve to calculate the time evolution;
• visualised the state probabilities with Plotly;
• interpreted the same motion using σz and the Bloch sphere;
• changed the Rabi frequency to see how the control speed changes.

2. Rabi oscillations as a model of coherent control

A two-level system is the simplest useful model in quantum technology. It is not always the whole physical system, but it often captures the part we care about: in an atom we may focus on two selected energy levels, in a spin system on the spin-up and spin-down states, and in a qubit on the computational basis states we label |0⟩ and |1⟩.

A Rabi oscillation describes what happens when we coherently drive such a system. If the system begins in |0⟩, the drive can move population into |1⟩; but if the drive continues, the system does not simply stay there — it moves back again, and the probability flows between the two states in a regular oscillation. This is a central idea behind quantum control. A pulse of the right duration can move the system from |0⟩ to |1⟩, while a shorter pulse can create a superposition. In practical systems the same idea underlies state preparation, qubit rotations, spectroscopy, and atomic-clock interrogation.

For this article we will begin with the cleanest possible version: a closed, resonantly driven two-level system. Later articles can add the effects that make real systems more complicated.

3. The model: a driven two-level system

We use the Hamiltonian:

H = (Ω / 2) σx

Here Ω is the Rabi angular frequency, σx is the Pauli-X operator, and the factor of one half is conventional when writing qubit rotations. The Pauli-X operator couples the two basis states; it is the operator that turns |0⟩ towards |1⟩ and |1⟩ towards |0⟩. Starting from |0⟩, this Hamiltonian rotates the state around the x-axis of the Bloch sphere, and in the probability picture that rotation appears as population moving back and forth between the two basis states.

4. Setting up Python in Google Colab

If you are using Google Colab, start by installing QuTiP and Plotly.

!pip -q install qutip plotly

Then import the libraries we need.

import numpy as np
import plotly.graph_objects as go

from qutip import basis, sigmax, sigmay, sigmaz, mesolve

We will use numpy for time arrays and numerical values, plotly for interactive visualisation, basis to create |0⟩ and |1⟩, sigmax, sigmay and sigmaz for the Pauli operators, and mesolve for the time evolution. The overall workflow is the same as before:

define states → define operators → define Hamiltonian → solve dynamics → visualise

5. Defining the basis states and projectors

A two-level system has a Hilbert space of dimension two. In QuTiP we can define the two basis states as follows.

ket_0 = basis(2, 0)
ket_1 = basis(2, 1)

Here basis(2, 0) represents |0⟩ and basis(2, 1) represents |1⟩. We choose the initial state to be the ground state, so that the system starts in |0⟩.

psi0 = ket_0

Next we define the projection operators.

P0 = ket_0 * ket_0.dag()
P1 = ket_1 * ket_1.dag()

Mathematically these are P0 = |0⟩⟨0| and P1 = |1⟩⟨1|. The expectation value of P0 gives the probability of measuring the system in |0⟩, and the expectation value of P1 gives the probability of measuring it in |1⟩. This is a useful QuTiP pattern to keep in mind:

state → projector → expectation value → measurement probability

6. Defining the Rabi Hamiltonian

Now we define the Hamiltonian.

Omega = 2 * np.pi  # Rabi angular frequency

H = 0.5 * Omega * sigmax()

With this choice the Rabi angular frequency is Ω = 2π, so the population oscillation has a period of one unit of time. Before running the simulation we can already predict the main behaviour:

• t = 0.0 — mostly |0⟩
• t = 0.5 — mostly |1⟩
• t = 1.0 — back to |0⟩

This kind of prediction is worth making, because it gives us a sanity check against which to compare the numerical result.

7. Solving the time evolution with mesolve

Now we create a time grid that spans three full Rabi cycles.

tlist = np.linspace(0, 3, 600)

Next we solve the time evolution.

result = mesolve(
    H,
    psi0,
    tlist,
    c_ops=[],
    e_ops=[P0, P1, sigmax(), sigmay(), sigmaz()]
)

This is the central QuTiP line in the article. The arguments mean the following:

mesolve(
    H,                  # Hamiltonian
    psi0,               # initial state
    tlist,              # times at which to evaluate the system
    c_ops=[],           # no collapse operators
    e_ops=[...]         # observables to calculate
)

The empty c_ops=[] list means that we are not including relaxation, dephasing, or any other dissipative process; this is therefore a closed-system calculation. The e_ops list tells QuTiP which expectation values to return, and here we request five quantities in a single solver call:

• P0 — probability of |0⟩
• P1 — probability of |1⟩
• σx — x-component of the Bloch vector
• σy — y-component of the Bloch vector
• σz — z-component of the Bloch vector

We then extract the results, which the result object stores in the same order as e_ops.

prob_0 = np.real(result.expect[0])
prob_1 = np.real(result.expect[1])

sx = np.real(result.expect[2])
sy = np.real(result.expect[3])
sz = np.real(result.expect[4])

8. Plotting the Rabi oscillation

Now we plot the two probabilities.

fig = go.Figure()

fig.add_trace(
    go.Scatter(
        x=tlist,
        y=prob_0,
        mode="lines",
        name="Probability of |0⟩"
    )
)

fig.add_trace(
    go.Scatter(
        x=tlist,
        y=prob_1,
        mode="lines",
        name="Probability of |1⟩"
    )
)

fig.update_layout(
    title="Rabi oscillations in a two-level quantum system",
    xaxis_title="Time",
    yaxis_title="Probability",
    yaxis=dict(range=[0, 1.05]),
    template="plotly_white"
)
fig.show()

The result should show two oscillating curves that mirror each other: when the probability of |0⟩ is high, the probability of |1⟩ is low, and vice versa. This is the simplest visual picture of coherent population transfer. The system is not randomly jumping between states; it is evolving smoothly according to the Hamiltonian.

9. Reading the result

Because this is a closed two-level system, the total probability should remain equal to one throughout the evolution. It is worth checking this directly rather than trusting the plot alone.

max_probability_error = np.max(np.abs(prob_0 + prob_1 - 1))

print("Maximum probability error:", max_probability_error)

The value should be very small, close to numerical precision. This is a good habit when writing simulations: do not only look at the plot, but also check a conservation rule whenever one is available. In this example the rule is simply that P(|0⟩) + P(|1⟩) = 1.

Get Ken Moriwaki’s stories in your inbox

Join Medium for free to get updates from this writer.

Remember me for faster sign in

We also calculated σz. In QuTiP’s convention |0⟩ has a σz expectation value of +1 and |1⟩ has a σz expectation value of −1, so the σz expectation value is exactly the population difference:

⟨σz⟩ = P(|0⟩) − P(|1⟩)

Let us plot it.

fig = go.Figure()

fig.add_trace(
    go.Scatter(
        x=tlist,
        y=sz,
        mode="lines",
        name="⟨σz⟩"
    )
)

fig.update_layout(
    title="Population difference during Rabi oscillation",
    xaxis_title="Time",
    yaxis_title="⟨σz⟩",
    yaxis=dict(range=[-1.05, 1.05]),
    template="plotly_white"
)

fig.show()

The probability plot shows the occupation of each state, while the σz plot shows the population imbalance. They are two views of the same dynamics.

10. Visualising the Bloch-vector trajectory

A two-level quantum state can also be represented on the Bloch sphere. For a pure qubit state, the expectation values of the three Pauli operators give the Bloch-vector coordinates x = ⟨σx⟩, y = ⟨σy⟩ and z = ⟨σz⟩. We already calculated these values in the original mesolve call, so there is no need to solve the dynamics again.

fig = go.Figure()

fig.add_trace(
    go.Scatter3d(
        x=sx,
        y=sy,
        z=sz,
        mode="lines",
        name="State trajectory"
    )
)

fig.add_trace(
    go.Scatter3d(
        x=[sx[0]],
        y=[sy[0]],
        z=[sz[0]],
        mode="markers",
        name="Initial state"
    )
)

fig.update_layout(
    title=dict(text="Rabi oscillation as a Bloch-vector trajectory", y=0.97),
    scene=dict(
        xaxis_title="⟨σx⟩",
        yaxis_title="⟨σy⟩",
        zaxis_title="⟨σz⟩",
        xaxis=dict(range=[-1, 1]),
        yaxis=dict(range=[-1, 1]),
        zaxis=dict(range=[-1, 1])
    ),
    template="plotly_white",
    legend=dict(
        orientation="h",
        yanchor="bottom",
        y=1.0,
        xanchor="center",
        x=0.5
    ),
    margin=dict(l=0, r=0, t=70, b=0)
)

fig.show()

This view connects the probability oscillation to the geometry of a qubit. Because the Hamiltonian is proportional to σx, the state rotates around the x-axis; and since the initial state is |0⟩, it starts near the north pole of the Bloch sphere and traces a circle in the y–z plane during the oscillation. The probability plot and the Bloch-vector plot are not separate simulations but two representations of the same quantum evolution.

11. Changing the Rabi frequency

The parameter Ω controls how quickly the system oscillates. To explore this, we can wrap the simulation in a small function.

def simulate_rabi(Omega, tlist, initial_state, excited_projector):
    H = 0.5 * Omega * sigmax()
    result = mesolve(
        H,
        initial_state,
        tlist,
        c_ops=[],
        e_ops=[excited_projector]
    )
    return np.real(result.expect[0])

Now we compare different Rabi frequencies.

tlist = np.linspace(0, 3, 600)

Omega_values = [
    2 * np.pi * 0.5,
    2 * np.pi * 1.0,
    2 * np.pi * 2.0
]

fig = go.Figure()

for Omega_value in Omega_values:
    prob_excited = simulate_rabi(
        Omega=Omega_value,
        tlist=tlist,
        initial_state=psi0,
        excited_projector=P1
    )
    fig.add_trace(
        go.Scatter(
            x=tlist,
            y=prob_excited,
            mode="lines",
            name=f"Ω / 2π = {Omega_value / (2 * np.pi):.1f}"
        )
    )

fig.update_layout(
    title="Effect of changing the Rabi frequency",
    xaxis_title="Time",
    yaxis_title="Probability of |1⟩",
    yaxis=dict(range=[0, 1.05]),
    template="plotly_white"
)

fig.show()

As Ω increases, the oscillation becomes faster. This matches the physical interpretation: a stronger resonant drive produces faster coherent population transfer. It is here that the simulation begins to look like a control problem, because the duration of the pulse determines what operation we perform:

• a π pulse moves population from |0⟩ to |1⟩;
• a π/2 pulse creates an equal superposition;
• a 2π pulse returns the observable populations to their starting values.

Strictly speaking, for a spin-1/2 state vector a 2π rotation introduces a global phase of −1. This phase does not change the probabilities or the Bloch vector, so it is invisible in the plots above; for the observables we measure here, the system has returned to the same physical state.

12. What this example leaves out

This simulation is intentionally idealised. There is no relaxation, no dephasing, no detuning, and no pulse envelope; the drive is perfectly resonant and perfectly coherent. That is not a weakness of the example but its purpose. Before studying imperfect quantum systems we need to understand the ideal motion, and a Rabi oscillation in a closed two-level system gives us a clean reference case against which more realistic simulations can be compared.

Later, we can add a number of effects that bring the model closer to real systems:

• detuning — the drive is not exactly resonant;
• relaxation — the excited state decays;
• dephasing — phase coherence is lost;
• pulse shapes — the drive turns on and off smoothly;
• open-system dynamics — the system interacts with its environment.

QuTiP is well suited to this because the same basic modelling structure can be extended. In this article the collapse-operator list was empty; in a later article we can add collapse operators to model relaxation and dephasing.

13. Complete code

Here is the full code in one block.

#!pip -q install qutip plotly

import numpy as np
import plotly.graph_objects as go

from qutip import basis, sigmax, sigmay, sigmaz, mesolve

# Basis states
ket_0 = basis(2, 0)
ket_1 = basis(2, 1)

# Initial state
psi0 = ket_0

# Projectors
P0 = ket_0 * ket_0.dag()
P1 = ket_1 * ket_1.dag()

# Rabi frequency
Omega = 2 * np.pi

# Hamiltonian
H = 0.5 * Omega * sigmax()

# Time grid
tlist = np.linspace(0, 3, 600)

# Solve dynamics once and request all observables we need
result = mesolve(
    H,
    psi0,
    tlist,
    c_ops=[],
    e_ops=[P0, P1, sigmax(), sigmay(), sigmaz()]
)

# Extract expectation values
prob_0 = np.real(result.expect[0])
prob_1 = np.real(result.expect[1])
sx = np.real(result.expect[2])
sy = np.real(result.expect[3])
sz = np.real(result.expect[4])

# Plot probabilities
fig = go.Figure()

fig.add_trace(
    go.Scatter(
        x=tlist,
        y=prob_0,
        mode="lines",
        name="Probability of |0⟩"
    )
)

fig.add_trace(
    go.Scatter(
        x=tlist,
        y=prob_1,
        mode="lines",
        name="Probability of |1⟩"
    )
)

fig.update_layout(
    title="Rabi oscillations in a two-level quantum system",
    xaxis_title="Time",
    yaxis_title="Probability",
    yaxis=dict(range=[0, 1.05]),
    template="plotly_white"
)

fig.show()

# Check probability conservation
max_probability_error = np.max(np.abs(prob_0 + prob_1 - 1))
print("Maximum probability error:", max_probability_error)

# Plot σz expectation value
fig = go.Figure()

fig.add_trace(
    go.Scatter(
        x=tlist,
        y=sz,
        mode="lines",
        name="⟨σz⟩"
    )
)

fig.update_layout(
    title="Population difference during Rabi oscillation",
    xaxis_title="Time",
    yaxis_title="⟨σz⟩",
    yaxis=dict(range=[-1.05, 1.05]),
    template="plotly_white"
)

fig.show()

# Bloch-vector trajectory
fig = go.Figure()
fig.add_trace(
    go.Scatter3d(
        x=sx,
        y=sy,
        z=sz,
        mode="lines",
        name="State trajectory"
    )
)

fig.add_trace(
    go.Scatter3d(
        x=[sx[0]],
        y=[sy[0]],
        z=[sz[0]],
        mode="markers",
        name="Initial state"
    )
)

fig.update_layout(
    title=dict(text="Rabi oscillation as a Bloch-vector trajectory", y=0.97),
    scene=dict(
        xaxis_title="⟨σx⟩",
        yaxis_title="⟨σy⟩",
        zaxis_title="⟨σz⟩",
        xaxis=dict(range=[-1, 1]),
        yaxis=dict(range=[-1, 1]),
        zaxis=dict(range=[-1, 1])
    ),
    template="plotly_white",
    legend=dict(
        orientation="h",
        yanchor="bottom",
        y=1.0,
        xanchor="center",
        x=0.5
    ),
    margin=dict(l=0, r=0, t=70, b=0)
)

fig.show()

# Compare different Rabi frequencies
def simulate_rabi(Omega, tlist, initial_state, excited_projector):
    H = 0.5 * Omega * sigmax()
    result = mesolve(
        H,
        initial_state,
        tlist,
        c_ops=[],
        e_ops=[excited_projector]
    )
    return np.real(result.expect[0])

Omega_values = [
    2 * np.pi * 0.5,
    2 * np.pi * 1.0,
    2 * np.pi * 2.0
]
fig = go.Figure()

for Omega_value in Omega_values:
    prob_excited = simulate_rabi(
        Omega=Omega_value,
        tlist=tlist,
        initial_state=psi0,
        excited_projector=P1
    )
    fig.add_trace(
        go.Scatter(
            x=tlist,
            y=prob_excited,
            mode="lines",
            name=f"Ω / 2π = {Omega_value / (2 * np.pi):.1f}"
        )
    )

fig.update_layout(
    title="Effect of changing the Rabi frequency",
    xaxis_title="Time",
    yaxis_title="Probability of |1⟩",
    yaxis=dict(range=[0, 1.05]),
    template="plotly_white"
)

fig.show()

14. Summary

In this article we simulated Rabi oscillations using QuTiP. We started with a two-level quantum system, placed it in the |0⟩ state, and applied a Hamiltonian proportional to σx, which caused the system to move coherently between |0⟩ and |1⟩. The main workflow was to define the basis states, define the projectors, define the Hamiltonian, run mesolve, extract the expectation values, and finally plot and interpret the dynamics.

The key idea is that a Rabi oscillation is a minimal model of coherent quantum control: it shows how a quantum system can be driven predictably between states. This prepares us for the next step — realistic quantum dynamics. Once the ideal oscillation is clear, we can ask what happens when the system is noisy, dissipative, or coupled to its environment.

Further reading

• Nation, P.D. (2026). QuTiP — Quantum Toolbox in Python. [online] Qutip.org. Available at: https://qutip.org
• ‌QuTiP developers and contributors (2026). QuTiP: Quantum Toolbox in Python — QuTiP 5.3 Documentation. [online] Available at: https://qutip.readthedocs.io
• Lambert, N., Giguère, E., Menczel, P., Li, B., Hopf, P., Suárez, G., Gali, M., Lishman, J., Rushiraj Gadhvi, Agarwal, R., Galicia, A., Shammah, N., Nation, P., Johansson, J.R., Ahmed, S., Cross, S., Pitchford, A. and Nori, F. (2025). QuTiP 5: The Quantum Toolbox in Python. Physics Reports, 1153, pp.1–62. doi:https://doi.org/10.1016/j.physrep.2025.10.001.
• Johansson, J.R., Nation, P.D. and Nori, F. (2012). QuTiP: An open-source Python framework for the dynamics of open quantum systems. Computer Physics Communications, 183(8), pp.1760–1772. doi:https://doi.org/10.1016/j.cpc.2012.02.021.
• Johansson, J.R., Nation, P.D. and Nori, F. (2013). QuTiP 2: A Python framework for the dynamics of open quantum systems. Computer Physics Communications, 184(4), pp.1234–1240. doi:https://doi.org/10.1016/j.cpc.2012.11.019.