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$, $E_L$0 such that $E_L$1, with detuning $E_L$2, ensuring strong anharmonicity with higher levels ($E_L$3). Plasmonic excitations correspond to small oscillations around a single well minimum; fluxonic excitations entail 2$E_L$4 phase slips between wells.

2. Qutrit Basis: Matrix Elements and Normalized Operators

Fock basis truncation to $E_L$5 defines local “photon number” $E_L$6. The relevant dipole matrix elements are $E_L$7 and $E_L$8. The normalized bosonic creation operator is

$E_L$9

such that $C_{\Sigma}$0 raises $C_{\Sigma}$1 (limited to $C_{\Sigma}$2).

Two dimensionless parameters are central:

• $C_{\Sigma}$3 (or analogously for $C_{\Sigma}$4), which controls the density-dependence of the hopping,
• $C_{\Sigma}$5, which serves as an on-site interaction strength.

Depending on the plasmonic/fluxonic character, one finds $C_{\Sigma}$6 (plasmonic) or $C_{\Sigma}$7 (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	$C_{\Sigma}$8, $C_{\Sigma}$9; all transitions are intra-well oscillations
ΦΦ	fluxon–fluxon	$\Phi_\mathrm{ext}$0, $\Phi_\mathrm{ext}$1; both transitions via phase slips, pair hopping dominates
ΠΦ	plasmon–fluxon	$\Phi_\mathrm{ext}$2, $\Phi_\mathrm{ext}$3; 0⇄1 plasmon, 1⇄2 fluxon, single hopping suppressed
ΦΠ	fluxon–plasmon	$\Phi_\mathrm{ext}$4, $\Phi_\mathrm{ext}$5; 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 $\Phi_\mathrm{ext}$6 or inductive $\Phi_\mathrm{ext}$7 couplings is captured by a generalized Bose–Hubbard Hamiltonian in the rotating-wave approximation:

$\Phi_\mathrm{ext}$8

with $\Phi_\mathrm{ext}$9 and $$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$$0 determined by dipole matrix elements, $$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$$1 governing the occupation dependence of hopping, and a three-body hard core ($$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$$2) truncation. Non-local interactions $$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$$3 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 $$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$$4 nonlinear phase nodes. Eigenmode frequencies are set by:

$$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$$5

where $$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$$6, and $$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$$7 is determined by the spectrum of the capacitance matrix, parameterized by array and grounding capacitances ($$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$$8).

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

$$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$$9

Design guidelines require the lowest array mode frequency $E_C = e^2/(2C_{\Sigma})$0 well above thermal energy and drive frequencies ($E_C = e^2/(2C_{\Sigma})$1 GHz for typical parameters), and low ground capacitance ($E_C = e^2/(2C_{\Sigma})$2) for maximal anharmonicity. Parasitic couplings are minimized by suppressing $E_C = e^2/(2C_{\Sigma})$3 and $E_C = e^2/(2C_{\Sigma})$4, ensuring dispersive separation ($E_C = e^2/(2C_{\Sigma})$5).

To mitigate mode degeneracy and crosstalk in arrays, spectral non-degeneracy is introduced by varying $E_C = e^2/(2C_{\Sigma})$6 slightly across devices. Small-junction node shielding reduces stray capacitance and hybridization.

6. Phase Diagram and Dynamical Probes

At unit filling ($E_C = e^2/(2C_{\Sigma})$7), 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 $E_C = e^2/(2C_{\Sigma})$8, single-particle coherence $E_C = e^2/(2C_{\Sigma})$9, and pair coherence $\Phi_0 = h/(2e)$0. Analytical mean-field phase boundaries are, for coordination $\Phi_0 = h/(2e)$1:

• MI–SF: $\Phi_0 = h/(2e)$2
• SF–PSF: $\Phi_0 = h/(2e)$3

Dynamical experiments (“quench and watch”) are proposed: prepare local Fock state patterns, activate coupling, and monitor $\Phi_0 = h/(2e)$4 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, $\Phi_0 = h/(2e)$5, emulates gauge–matter coupling, while pair hopping implements plaquette/ring-exchange operators critical for $\Phi_0 = h/(2e)$6 or $\Phi_0 = h/(2e)$7 lattice gauge theories. The three-body hard-core constraint and tunable $\Phi_0 = h/(2e)$8 map directly onto models such as the bosonic Pfaffian (Moore–Read) state at filling $\Phi_0 = h/(2e)$9. The non-local term $$H_{\mathrm{at}}(\Phi_\mathrm{ext})\,|\psi_a(\Phi_\mathrm{ext})\rangle = \omega_a(\Phi_\mathrm{ext})\,|\psi_a\rangle$$0, 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)]