Huge congratulations to John Clarke, Michel H. Devoret, and John M. Martinis on the 2025 Nobel Prize in Physics “for the discovery of macroscopic quantum mechanical tunnelling and energy quantisation in an electric circuit.” Their superconducting Josephson-circuit experiments made quantum effects unmistakably visible at circuit scale, discrete, anharmonic energy levels and coherent tunnelling between macroscopically distinct states, laying key groundwork for modern superconducting qubits. In this short computational essay, we’ll walk through compact simulations that reproduce those signatures: a spectroscopy-style level map for the Cooper-pair box/transmon, and time-domain tunnelling dynamics with realistic decoherence to mirror the original observations.

Energy quantization in a Josephson circuit (Cooper-pair box / transmon toy)

A small superconducting island, linked by a Josephson junction and capacitively coupled to a gate, is described by conjugate variables: the superconducting phase across the junction and the integer number of Cooper pairs on the island. The key energy scales are the charging energy set by the island capacitance, the Josephson energy set by the junction’s critical current, and the offset charge controlled by the gate. In the charge basis the Hamiltonian becomes a real, symmetric, tridiagonal matrix with a diagonal charging term and nearest-neighbor couplings from the cosine Josephson term; diagonalization yields discrete, anharmonic energy levels. As the ratio of Josephson to charging energy increases, the device evolves from a charge-qubit regime with strong offset-charge dispersion and prominent avoided crossings near half-integer offset to a transmon regime with suppressed dispersion and a weakly anharmonic ladder whose anharmonicity is set by the charging scale. Microwave spectroscopy versus gate offset (or flux in related devices) reveals these discrete transitions and their anharmonicity, while linewidths report the underlying coherence times.
Consider the Hamiltonian as

H
=4
2
E
C


n
-
n
g



-
E
J
cos

ϕ

with

n
the Cooper-pair number with integer eigenvalues, and ϕ the superconducting phase across the junction. They are canonically conjugate:


ϕ
,

n
=
and

n
=-
∂
∂ϕ
. Additionally,
E
C
is the charging energy,
E
J
is the Josephson energy, and
n
g
is the gate offset charge. The Hamiltonian can be considered as a quantum rotor (phase) in a cosine potential, minimally coupled to a static “vector potential”
n
g
that shifts the momentum

n
.
The charge basis is denoted as
{|n〉}
nϵ
with

n
|n〉=n|n〉
. The phase operator acts as a ladder in this basis:
±

ϕ

|n〉=|n∓1〉
. Thus, this implies
cos

ϕ
=
1
2
∑
n
(|n〉〈n+1|+|n+1〉〈n|)
. Put all of these together, the Hamiltonian can be described in a truncated space, by setting
-
n
max
<=n<=
n
max
In[]:=
ClearAll[hamiltonian];​​hamiltonian[ratio_,ng_,nmax_:20]:=DiagonalMatrixTable4
2
(n-ng)
,{n,-nmax,nmax}-
ratio(*EJ/EC*)
2
SparseArray[{Band[{1,2}]1,Band[{2,1}]1},{2nmax+1,2nmax+1}]
Given above Hamiltonian (truncated version, find eigenvalues as a function of
n
g
):
In[]:=
​
E
J

E
C
E
1
E
2
E
3
E
4
E
5
E
6
E
7
E
8
E
9
E
10
E
1
E
2
The level curvature and avoided crossings reflect Josephson coupling, direct evidence of quantized levels in an electric circuit (a key piece of the laureates’ work).

Macroscopic quantum tunnelling (two-level double-well model)

Superconducting circuits can realize an effective double-well potential in phase or flux, hosting states localized in left and right wells that are macroscopically distinct. Restricting to the two lowest localized states gives a minimal model with a tunnel coupling, which mixes the wells, and a tunable bias, which tilts their relative energies. If the system is prepared in one well, the population coherently oscillates between wells with a frequency set by the combined effect of coupling and bias and with an amplitude that is maximal at symmetry and decreases as the bias grows. Pure dephasing damps the oscillation contrast, energy relaxation drives the system toward the lower-energy well’s stationary population, and finite temperature can repopulate the higher level; the observed coherence time reflects both relaxation and dephasing. Time-domain oscillations, spectroscopy of the bias-dependent splitting near symmetry, and escape-rate or switching-current measurements that crossover from thermal activation at high temperature to temperature-independent quantum tunnelling at low temperature together provide the experimental fingerprints of this macroscopic quantum behavior.
Let’s define a two-level system by the Hamiltonian
H=
ϵ
2
σ
z
+
Δ
2
σ
x
:
In[]:=
$Assumptions={ϵ,Δ}>0;​​h=
ϵ
2
PauliMatrix[3]+
Δ
2
PauliMatrix[2];​​h//MatrixForm
Out[]//MatrixForm=
ϵ
2
-
Δ
2
Δ
2
-
ϵ
2
Here ϵ represents the bias (tilt) between wells, and Δ is the tunnel splitting.
Find the eigenvalues of the Hamiltonian:
In[]:=
{
ℰ
1
,
ℰ
2
}=Eigenvalues[h]
Out[]=
-
1
2
2
Δ
+
2
ϵ
,
1
2
2
Δ
+
2
ϵ

Given above result,
Ω=
2
Δ
+
2
ϵ
is usually called the eigenfrequency.
In[]:=
{
e
1
,
e
2
}=Normalize/@Eigenvectors[h]
Out[]=

-ϵ+
2
Δ
+
2
ϵ

Δ
1+
2
Abs
-ϵ+
2
Δ
+
2
ϵ
Δ

,
1
1+
2
Abs
-ϵ+
2
Δ
+
2
ϵ
Δ

,-
ϵ+
2
Δ
+
2
ϵ

Δ
1+
2
Abs
ϵ+
2
Δ
+
2
ϵ
Δ

,
1
1+
2
Abs
ϵ+
2
Δ
+
2
ϵ
Δ


Set the initial state at the bottom of one of the wells (i.e., the eigenstate of
σ
z
):
In[]:=
ψ
0
={1,0};
Using the Schrodinger equation, find the final state at a given time:
In[]:=
ψ[t_]=
2
∑
i=1
Exp[-I
ℰ
i
t]
e
i
.
ψ
0
e
i
//FullSimplify
Out[]=
Cos
1
2
t
2
Δ
+
2
ϵ
-
ϵSin
1
2
t
2
Δ
+
2
ϵ

2
Δ
+
2
ϵ
,
ΔSin
1
2
t
2
Δ
+
2
ϵ

2
Δ
+
2
ϵ

Calculate the probability of tunneling:
In[]:=
{0,1}.ψ[t]{0,1}.ψ[t]//FullSimplify
Out[]=
2
Δ
2
Sin
1
2
t
2
Δ
+
2
ϵ

2
Δ
+
2
ϵ
The above behavior assumes coherent tunneling, i.e., no noise. In practice, environmental interactions introduce dephasing, relaxation, and other channels, so the dynamics deviate from the Schrödinger equation and the qubit becomes an open system. A common model for such open-system dynamics is the Lindblad master equation:
∂
t
ρ=-[H,ρ]+
∑
k
[
L
k
][ρ]
with
[
L
k
][ρ]=
L
k
ρ
†
L
k
-
1
2

†
L
k
L
k
,ρ
with
L
k
the corresponding Lindblad/jump operators.
Define the Lindblad terms:
In[]:=
ClearAll[]​​[L_,ρ_]:=L.ρ.ConjugateTranspose[L]-
1
2
(ConjugateTranspose[L].L.ρ+ρ.ConjugateTranspose[L].L)
For the case of energy relaxation (downward transitions), we have:
∂
t
ρ=-[H,ρ]+γ
σ
-
ρσ
+
-
1
2
{
σ
+
σ
-
,ρ}
with
γ=1/
T
1
For the case of pure dephasing (phase randomization without energy exchange), one gets:
∂
t
ρ=-[H,ρ]+
Γ
ϕ
2
(
σ
z
ρσ
z
-ρ)
with
Γ
ϕ
=1
T
ϕ
Putting them together, the total transverse coherence (T₂) combines both mechanisms
1
T
2
=
1
2
T
1
+
1
T
ϕ
Define a vector such that
ρ=
1
2
(+

r
.
σ
)
where 3-vector

r
is usually called the Bloch vector:
r={1,rx[t],ry[t],rz[t]};
Note the transformation from Bloch vector in the Cartesian coordinate to spherical and vice versa:
In[]:=
CoordinateTransformData"Cartesian""Spherical","Mapping",x.,y.,z.
Out[]=

2
x.
+
2
y.
+
2
z.
,ArcTanz.,
2
x.
+
2
y.
,ArcTan[x.,y.]
In[]:=
CoordinateTransformData"Spherical"->"Cartesian","Mapping",r.,θ,ϕ//FullSimplify
Out[]=
r.Cos[ϕ]Sin[θ],r.Sin[θ]Sin[ϕ],r.Cos[θ]
Find the density matrix:
ρ=
1
2
r.Table[PauliMatrix[j],{j,0,3}]//FullSimplify;​​ρ//MatrixForm
Out[]//MatrixForm=
1
2
(1+rz[t])
1
2
(rx[t]-ry[t])
1
2
(rx[t]+ry[t])
1
2
(1-rz[t])
Define the jump operator for relaxation:
In[]:=
σs={{0,0},{0,1}};
Put relaxation and pure dephasing together into dynamical equation:
In[]:=
dρt=-(h.ρ-ρ.h)+γ[σs,ρ]+
Γϕ
2
[PauliMatrix[3],ρ]//FullSimplify
Out[]=
-
1
2
Δrx[t],
1
4
(-((γ+2Γϕ+2ϵ)(rx[t]-ry[t]))+2Δrz[t]),
1
4
(-((γ+2Γϕ-2ϵ)(rx[t]+ry[t]))+2Δrz[t]),
1
2
Δrx[t]
Given the evolution of the density matrix, find the evolution of the Bloch vector:
In[]:=
drt=Table[Tr[dρt.PauliMatrix[j]]//FullSimplify,{j,3}]
Out[]=
-
1
2
(γ+2Γϕ)rx[t]-ϵry[t]+Δrz[t],ϵrx[t]-
1
2
(γ+2Γϕ)ry[t],-Δrx[t]
Set the dynamical equations and the initial conditions:
In[]:=
eqns=Thread[D[{rx[t],ry[t],rz[t]},t]==drt]​​ic={rx[0]==rx0,ry[0]==ry0,rz[0]==rz0}
Out[]=

′
rx
[t]-
1
2
(γ+2Γϕ)rx[t]-ϵry[t]+Δrz[t],
′
ry
[t]ϵrx[t]-
1
2
(γ+2Γϕ)ry[t],
′
rz
[t]-Δrx[t]
Out[]=
{rx[0]rx0,ry[0]ry0,rz[0]rz0}
Put all of them into a parametric DEs:
In[]:=
pnd=ParametricNDSolveValue[Join[eqns,ic],{rx,ry,rz},{t,0,T},{rx0,ry0,rz0,Δ,γ,Γϕ,ϵ,T}]
Out[]=
ParametricFunction
Expression: {rx,ry,rz}
Parameters: {rx0,ry0,rz0,Δ,γ,Γϕ,ϵ,T}

Visualize the parametric DEs, for different values of the corresponding parameters:
Out[]=
​
ϕ
θ
Δ
γ
Γϕ
ϵ
T
〈
σ
x
〉
〈
σ
y
〉
〈
σ
z
〉
Tunneling prob.

CITE THIS NOTEBOOK

Nobel Prize in Physics 2025: Macroscopic Quantum Effects and the Dawn of Quantum Computer​
by Mohammad Bahrami​
Wolfram Community, STAFF PICKS, October 8, 2025
​https://community.wolfram.com/groups/-/m/t/3558195