Fluxonium Qutrit Arrays for Quantum Simulation

• The paper demonstrates how fluxonium devices realize coherent qutrit operation by selecting three energy levels with strong anharmonicity, essential for simulating extended Bose–Hubbard models.
• Fluxonium qutrit arrays are superconducting circuits where external flux bias controls plasmonic and fluxonic excitations, enabling density-dependent single-particle and pair hopping.
• The design supports exploring many-body phenomena like superfluid, Mott insulator, and topologically ordered phases, while providing a testbed for lattice gauge theories.

Fluxonium qutrit arrays constitute a superconducting circuit platform where each site comprises a highly coherent fluxonium device engineered to realize a qutrit—an effective three-level system. The sites’ energy spectra and matrix elements are controlled by external flux bias, enabling operational regimes dominated by either plasmon-like or fluxon-like excitations. The interacting array realizes an extended Bose–Hubbard model with density-dependent single-particle hopping, correlated pair hopping, strong on-site interactions, and non-local couplings, offering a versatile testbed for quantum simulation of strongly correlated bosonic phases and lattice gauge models (Amelio et al., 29 Jan 2026, Sorokanich et al., 2024).

1. Single-Site Physics and Qutrit Encoding

A fluxonium device comprises a Josephson junction (energy $E_J$) shunted by a large linear inductance ("superinductor," energy $E_L$) and typical total capacitance $C_{\Sigma}$. The single-site Hamiltonian as a function of external flux $\Phi_\mathrm{ext}$ is:

$$H_{\mathrm{at}}(\Phi_\mathrm{ext}) = 4E_C\,\hat{n}^2 - E_J\cos(\hat{\phi} + 2\pi\Phi_\mathrm{ext}/\Phi_0) + \frac{E_L}{2}\hat{\phi}^2$$

where $E_C = e^2/(2C_{\Sigma})$, and $\Phi_0 = h/(2e)$. The eigenproblem $$H_{\mathrm{at}}(\Phi_\mathrm{ext})\,|\psi_a(\Phi_\mathrm{ext})\rangle = \omega_a(\Phi_\mathrm{ext})\,|\psi_a\rangle$$ is solved numerically.

Qutrit operation selects three local levels $|0\rangle$, $|1\rangle$, $|2\rangle$ such that $\omega_{21} \approx \omega_{10}$, with detuning $\Delta = \omega_{21} - \omega_{10} \approx 0$, ensuring strong anharmonicity with higher levels ($$\delta = \min |\omega_\mathrm{extra}-\omega_{10}| \gg g$$). Plasmonic excitations correspond to small oscillations around a single well minimum; fluxonic excitations entail 2$\pi$ phase slips between wells.

2. Qutrit Basis: Matrix Elements and Normalized Operators

Fock basis truncation to $\{|0\rangle,|1\rangle,|2\rangle\}$ defines local “photon number” $\hat{\rho}=0,1,2$. The relevant dipole matrix elements are $$n_{ab}(\Phi_\mathrm{ext}) = \langle a|\hat{n}|b\rangle$$ and $$\phi_{ab}(\Phi_\mathrm{ext}) = \langle a|\hat{\phi}|b\rangle$$. The normalized bosonic creation operator is

$$\hat{b}^\dagger = \sigma_{10} + \sqrt{2}\,\sigma_{21}, \quad \sigma_{ab}\equiv|a\rangle\langle b|$$

such that $\hat{b}^\dagger|{\rho}\rangle$ raises $\rho\to\rho+1$ (limited to $\rho\leq2$).

Two dimensionless parameters are central:

• $\alpha = |n_{21}/(\sqrt{2}n_{10})|$ (or analogously for $\phi$), which controls the density-dependence of the hopping,
• $\Delta = \omega_{21} - \omega_{10}$, which serves as an on-site interaction strength.

Depending on the plasmonic/fluxonic character, one finds $\alpha\sim1$ (plasmonic) or $\alpha\gg1,\alpha\ll1$ (fluxonic).

3. Operational Qutrit Regimes

Scanning the external flux and device parameters reveals four distinct regimes for qutrit transitions:

Regime	Excitation Type	Key Features
ΠΠ	plasmon–plasmon	$\alpha\approx1$, $P/J\sim0.3$; all transitions are intra-well oscillations
ΦΦ	fluxon–fluxon	$\alpha\gg1$, $P/J\gg1$; both transitions via phase slips, pair hopping dominates
ΠΦ	plasmon–fluxon	$\alpha\ll1$, $P/J\ll1$; 0⇄1 plasmon, 1⇄2 fluxon, single hopping suppressed
ΦΠ	fluxon–plasmon	$\alpha\gg1$, $P/J\sim O(1)$; 0⇄1 fluxon, 1⇄2 plasmon, moderate pair hopping

The qutrit subspace is robust since detuning to other levels is large compared to nearest-neighbor couplings.

4. Many-Body Model and Interactions

A chain or 2D array of fluxonium qutrits with nearest-neighbor capacitive $(g_C)$ or inductive $(g_L)$ couplings is captured by a generalized Bose–Hubbard Hamiltonian in the rotating-wave approximation:

$$H = - J \sum_{\langle ij\rangle} \alpha^{(\rho_i+\rho_j-1)} [\hat{b}_i^\dagger \hat{b}_j + \mathrm{h.c.}] + \frac{\Delta}{2} \sum_j (\hat{b}_j^\dagger)^2 \hat{b}_j^2 - \frac{P}{2} \sum_{\langle ij\rangle} [(\hat{b}_i^\dagger)^2 \hat{b}_j^2 + \mathrm{h.c.}] + \sum_{\langle ij\rangle} W(\rho_i, \rho_j)$$

with $J$ and $P$ determined by dipole matrix elements, $\alpha$ governing the occupation dependence of hopping, and a three-body hard core ($\rho_j\leq2$) truncation. Non-local interactions $W(\rho_i, \rho_j)$ originate primarily from persistent current matrix elements in inductive coupling.

Pair hopping, non-local interactions, and density-dependent hopping are tunable by external flux and circuit parameters, yielding rich physics beyond the canonical Bose-Hubbard model.

5. Array Mode Structure and Design Principles

The array’s linearized mode structure is solved exactly in terms of Chebyshev polynomial roots, with each array consisting of $N$ nonlinear phase nodes. Eigenmode frequencies are set by:

$\omega_\mu = \sqrt{8 E_J^a E_C^\mu}/\hbar$

where $E_C^\mu = e^2 \lambda_\mu/2$, and $\lambda_\mu$ is determined by the spectrum of the capacitance matrix, parameterized by array and grounding capacitances ($C_J, C_g^a, C_g^b$).

Eigenvectors have trigonometric spatial profiles (plane waves) after normalization. Approximations yield:

$$\omega_\mu \approx (\pi\mu/N)\sqrt{8E_J^a e^2/C_J}/\hbar + O(N^{-2}, C_g/C_J)$$

Design guidelines require the lowest array mode frequency $\omega_1$ well above thermal energy and drive frequencies ($\omega_1\gtrsim 10$ GHz for typical parameters), and low ground capacitance ($C_g^a/C_J\lesssim0.01$) for maximal anharmonicity. Parasitic couplings are minimized by suppressing $C_g^a$ and $C_g^b$, ensuring dispersive separation ($g_1\ll\Delta_\mathrm{qubit}-\omega_1$).

To mitigate mode degeneracy and crosstalk in arrays, spectral non-degeneracy is introduced by varying $N$ slightly across devices. Small-junction node shielding reduces stray capacitance and hybridization.

6. Phase Diagram and Dynamical Probes

At unit filling ($n=1$), the ground-state phase diagram encompasses superfluid (SF), pair superfluid (PSF), Mott insulator (MI), pair checkerboard (PCB), and clustered droplet (CL) phases. Order parameters include amplitudes $\psi_r=\langle\sigma_{rr}\rangle$, single-particle coherence $g^{(1)}$, and pair coherence $g_\mathrm{pair}$. Analytical mean-field phase boundaries are, for coordination $z$:

• MI–SF: $\Delta = z (1 + \sqrt{2}\,\alpha)^2 J$
• SF–PSF: $P = [2J/n (\sqrt{n-n^2/2}+n\alpha)^2 + \Delta/z]$

Dynamical experiments (“quench and watch”) are proposed: prepare local Fock state patterns, activate coupling, and monitor $\langle n_j(t)\rangle$ at each site. The space–time patterns distinguish SF (light-cone V-shapes), PSF (double V), heavy pair dispersion, cluster suppression, and checkerboard confinement.

7. Applications to Quantum Simulation and Topological Matter

The density-dependent hopping term, $\alpha^{(\rho_i+\rho_j-1)}$, emulates gauge–matter coupling, while pair hopping implements plaquette/ring-exchange operators critical for $Z_2$ or $U(1)$ lattice gauge theories. The three-body hard-core constraint and tunable $\Delta$ map directly onto models such as the bosonic Pfaffian (Moore–Read) state at filling $\nu=1$. The non-local term $W(\rho_i,\rho_j)$, particularly with synthetic gauge fluxes, enables simulation of anyon–Hubbard models and non-Abelian spin liquids, supporting explorations of topologically ordered regimes and lattice gauge dynamics beyond standard Bose–Hubbard physics.

This suggests fluxonium qutrit arrays offer a highly tunable, coherent, and theoretically tractable realization of complex quantum simulation platforms with extensibility toward gauge and topologically ordered phases.

• [Quantum Simulation with Fluxonium Qutrit Arrays, (Amelio et al., 29 Jan 2026)]
• [Exact and approximate fluxonium array modes, (Sorokanich et al., 2024)]