Two-Site Bose-Hubbard Model

• The two-site Bose-Hubbard model is a minimal interacting bosonic system defined by two coupled modes, modeling double-well dynamics, Josephson oscillations, and cat states.
• Its Fock-space formulation maps the dynamics to a one-dimensional tight-binding chain with nonuniform couplings and site-dependent energies, allowing both exact and approximate treatments.
• Algebraic reformulations via SU(2) currents and Bethe ansatz integrability reveal rich quantum phases and dynamical regimes, including pair tunnelling, collapse-revival phenomena, and self-trapping.

The two-site Bose-Hubbard model is the minimal interacting bosonic lattice model with two modes, usually interpreted as a bosonic double well, a bosonic Josephson junction, or a two-mode condensate. In its standard form it combines intersite tunnelling with onsite interaction, while conserving the total particle number. Because the fixed-$N$ sector is finite dimensional and admits both exact and controlled approximate treatments, the model has become a canonical setting for analyzing Josephson oscillations, self-trapping, fragmentation, Schrödinger-cat formation, current algebras, Bethe-ansatz integrability, and semiclassical collapse-and-revival dynamics (Longhi, 2011, Filho, 2015, Veksler et al., 2014).

1. Hamiltonian structure and Fock-space formulation

A standard parameterization for $N$ bosons in two weakly coupled modes is

$$\hat{H} = -\hbar J\left(\hat{a}_1^\dagger \hat{a}_2+\hat{a}_2^\dagger \hat{a}_1\right) + \frac{\hbar U}{2} \left( \hat{a}_1^{\dagger 2}\hat{a}_1^2 + \hat{a}_2^{\dagger 2}\hat{a}_2^2 \right),$$

where $J>0$ is the hopping or tunnelling amplitude and $U$ is the onsite interaction strength, with $U>0$ corresponding to repulsive interactions. A biased form,

$$H=-\varepsilon \left(a_L^\dagger a_L - a_R^\dagger a_R\right) - J \left(a_L^\dagger a_R + a_L a_R^\dagger\right) + \frac{U}{2}\left( a_L^\dagger a_L^\dagger a_L a_L + a_R^\dagger a_R^\dagger a_R a_R \right),$$

is used to select one branch of the attractive problem, while the canonical Josephson Hamiltonian is often written as

$$\hat H=\frac{K}{8}(\hat N_1-\hat N_2)^2-\frac{\Delta\mu}{2}(\hat N_1-\hat N_2) -\frac{\mathcal E_J}{2}\left(\hat a_1^\dagger \hat a_2+\hat a_2^\dagger \hat a_1\right).$$

More general two-site Hamiltonians allow $K_{ij}\hat N_i\hat N_j$, on-well energies, chemical potentials, and asymmetric tunnelling amplitudes $\Omega_{ij}$ (Longhi, 2011, Melé-Messeguer et al., 2011, Filho, 2015, Filho, 2015, Filho, 2016).

The conserved total number $N$0 decomposes the Hilbert space into fixed-particle-number sectors. In the occupation basis $N$1, with $N$2 bosons in the left well and $N$3 in the right well, the state amplitudes $N$4 obey a tridiagonal system,

$N$5

with

$N$6

This representation is structurally decisive: the two-site problem becomes a one-dimensional tight-binding chain in Fock space, with nonuniform nearest-neighbour couplings and site-dependent diagonal energies. In the fixed-$N$7 sector the model therefore interpolates between a tunnelling-dominated delocalized regime and interaction-dominated occupation-number localization (Longhi, 2011).

For attractive interactions, one frequently introduces the dimensionless control parameter

$N$8

which organizes the crossover from a nearly coherent condensate to fragmented and cat-like regimes. In the large-$N$9 classical description the bifurcation occurs near $$\hat{H} = -\hbar J\left(\hat{a}_1^\dagger \hat{a}_2+\hat{a}_2^\dagger \hat{a}_1\right) + \frac{\hbar U}{2} \left( \hat{a}_1^{\dagger 2}\hat{a}_1^2 + \hat{a}_2^{\dagger 2}\hat{a}_2^2 \right),$$0, but for finite $$\hat{H} = -\hbar J\left(\hat{a}_1^\dagger \hat{a}_2+\hat{a}_2^\dagger \hat{a}_1\right) + \frac{\hbar U}{2} \left( \hat{a}_1^{\dagger 2}\hat{a}_1^2 + \hat{a}_2^{\dagger 2}\hat{a}_2^2 \right),$$1 and weak bias the crossover is smooth and shifted (Melé-Messeguer et al., 2011).

2. Algebraic structure, spin mappings, and current variables

The two-site model admits an $$\hat{H} = -\hbar J\left(\hat{a}_1^\dagger \hat{a}_2+\hat{a}_2^\dagger \hat{a}_1\right) + \frac{\hbar U}{2} \left( \hat{a}_1^{\dagger 2}\hat{a}_1^2 + \hat{a}_2^{\dagger 2}\hat{a}_2^2 \right),$$2-type reformulation in terms of collective currents. A particularly useful choice is

$$\hat{H} = -\hbar J\left(\hat{a}_1^\dagger \hat{a}_2+\hat{a}_2^\dagger \hat{a}_1\right) + \frac{\hbar U}{2} \left( \hat{a}_1^{\dagger 2}\hat{a}_1^2 + \hat{a}_2^{\dagger 2}\hat{a}_2^2 \right),$$3

respectively the imbalance current, the Josephson tunnelling current, and the coherent correlation tunnelling current. They satisfy the closed algebra

$$\hat{H} = -\hbar J\left(\hat{a}_1^\dagger \hat{a}_2+\hat{a}_2^\dagger \hat{a}_1\right) + \frac{\hbar U}{2} \left( \hat{a}_1^{\dagger 2}\hat{a}_1^2 + \hat{a}_2^{\dagger 2}\hat{a}_2^2 \right),$$4

With $$\hat{H} = -\hbar J\left(\hat{a}_1^\dagger \hat{a}_2+\hat{a}_2^\dagger \hat{a}_1\right) + \frac{\hbar U}{2} \left( \hat{a}_1^{\dagger 2}\hat{a}_1^2 + \hat{a}_2^{\dagger 2}\hat{a}_2^2 \right),$$5, $$\hat{H} = -\hbar J\left(\hat{a}_1^\dagger \hat{a}_2+\hat{a}_2^\dagger \hat{a}_1\right) + \frac{\hbar U}{2} \left( \hat{a}_1^{\dagger 2}\hat{a}_1^2 + \hat{a}_2^{\dagger 2}\hat{a}_2^2 \right),$$6, and $$\hat{H} = -\hbar J\left(\hat{a}_1^\dagger \hat{a}_2+\hat{a}_2^\dagger \hat{a}_1\right) + \frac{\hbar U}{2} \left( \hat{a}_1^{\dagger 2}\hat{a}_1^2 + \hat{a}_2^{\dagger 2}\hat{a}_2^2 \right),$$7, this is the standard angular-momentum algebra (Filho, 2015, Filho, 2015).

The corresponding Casimirs are the conserved total number and the quadratic invariant

$$\hat{H} = -\hbar J\left(\hat{a}_1^\dagger \hat{a}_2+\hat{a}_2^\dagger \hat{a}_1\right) + \frac{\hbar U}{2} \left( \hat{a}_1^{\dagger 2}\hat{a}_1^2 + \hat{a}_2^{\dagger 2}\hat{a}_2^2 \right),$$8

This gives the model a Bloch-sphere geometry: the dynamics in a fixed-$$\hat{H} = -\hbar J\left(\hat{a}_1^\dagger \hat{a}_2+\hat{a}_2^\dagger \hat{a}_1\right) + \frac{\hbar U}{2} \left( \hat{a}_1^{\dagger 2}\hat{a}_1^2 + \hat{a}_2^{\dagger 2}\hat{a}_2^2 \right),$$9 sector is constrained to an effective spin of length $J>0$0. In semiclassical treatments this becomes explicit because the rescaled spin operators satisfy commutators proportional to $J>0$1, so $J>0$2 and the limit $J>0$3 is the classical limit (Filho, 2015, Veksler et al., 2014).

A closely related identification appears in the hopping-only dimer. Restricting the Hamiltonian

$J>0$4

to the $J>0$5-particle sector $J>0$6, one obtains exactly the spin projection operator $J>0$7 for spin $J>0$8. In the symmetric and antisymmetric modes

$J>0$9

the hopping Hamiltonian becomes

$U$0

making the equally spaced spectrum in each fixed-$U$1 sector manifest (Berezowski et al., 6 May 2026).

3. Integrability, Bethe ansatz, and explicit eigenstate constructions

The two-site Bose-Hubbard model is one of the standard bosonic systems accessible to the algebraic Bethe ansatz. In one integrable parameterization,

$U$2

the construction starts from the rational $U$3-invariant $U$4-matrix, a bosonic Lax operator, and the monodromy matrix $U$5. The Bethe vector state is

$U$6

with Bethe roots satisfying

$U$7

A binomial-operator method rewrites $U$8 as $U$9, expands the Bethe vector in powers of $U>0$0, and yields a fully explicit combinatorial Fock-basis expression. The same construction gives scalar products, norms, and imbalance-current form factors as finite sums in symmetric polynomials of the Bethe roots and recursively generated coefficients $U>0$1 (Santos et al., 2015).

A generalized two-site Hamiltonian with interaction matrix $U>0$2, local terms $U>0$3, and asymmetric tunnelling $U>0$4,

$U>0$5

remains exactly solvable by a new parametrization of a bosonic Lax operator. The transfer matrix generates commuting conserved quantities, and the Hamiltonian can be written as a combination of them. In the no-interaction limit $U>0$6 with $U>0$7, the Bethe equations collapse to a sphere constraint,

$U>0$8

which becomes an $U>0$9-sphere $$H=-\varepsilon \left(a_L^\dagger a_L - a_R^\dagger a_R\right) - J \left(a_L^\dagger a_R + a_L a_R^\dagger\right) + \frac{U}{2}\left( a_L^\dagger a_L^\dagger a_L a_L + a_R^\dagger a_R^\dagger a_R a_R \right),$$0 if all roots are real (Filho, 2016).

For the attractive model, Quantum Inverse Scattering Method results can be converted into high-accuracy closed energy formulas. Using the auxiliary Hamiltonian

$$H=-\varepsilon \left(a_L^\dagger a_L - a_R^\dagger a_R\right) - J \left(a_L^\dagger a_R + a_L a_R^\dagger\right) + \frac{U}{2}\left( a_L^\dagger a_L^\dagger a_L a_L + a_R^\dagger a_R^\dagger a_R a_R \right),$$1

the exact energies are written in terms of Bethe roots $$H=-\varepsilon \left(a_L^\dagger a_L - a_R^\dagger a_R\right) - J \left(a_L^\dagger a_R + a_L a_R^\dagger\right) + \frac{U}{2}\left( a_L^\dagger a_L^\dagger a_L a_L + a_R^\dagger a_R^\dagger a_R a_R \right),$$2. The ground-state roots are approximately real, negative, and nearly equidistant,

$$H=-\varepsilon \left(a_L^\dagger a_L - a_R^\dagger a_R\right) - J \left(a_L^\dagger a_R + a_L a_R^\dagger\right) + \frac{U}{2}\left( a_L^\dagger a_L^\dagger a_L a_L + a_R^\dagger a_R^\dagger a_R a_R \right),$$3

which yields a simple approximation to the ground-state energy and, after mapping back to the physical Hamiltonian, accurate formulas for both the ground and first excited states over a broad parameter range. The first excited state changes character at $$H=-\varepsilon \left(a_L^\dagger a_L - a_R^\dagger a_R\right) - J \left(a_L^\dagger a_R + a_L a_R^\dagger\right) + \frac{U}{2}\left( a_L^\dagger a_L^\dagger a_L a_L + a_R^\dagger a_R^\dagger a_R a_R \right),$$4, reflecting a reorganization of the root pattern (Ermakov et al., 2017).

4. Dynamical regimes: Josephson oscillations, self-trapping, pair tunnelling, and revivals

In the tunnelling-dominated regime the basic observable is the population imbalance. For the $$H=-\varepsilon \left(a_L^\dagger a_L - a_R^\dagger a_R\right) - J \left(a_L^\dagger a_R + a_L a_R^\dagger\right) + \frac{U}{2}\left( a_L^\dagger a_L^\dagger a_L a_L + a_R^\dagger a_R^\dagger a_R a_R \right),$$5-boson Fock-space chain the normalized imbalance can be written as

$$H=-\varepsilon \left(a_L^\dagger a_L - a_R^\dagger a_R\right) - J \left(a_L^\dagger a_R + a_L a_R^\dagger\right) + \frac{U}{2}\left( a_L^\dagger a_L^\dagger a_L a_L + a_R^\dagger a_R^\dagger a_R a_R \right),$$6

For $$H=-\varepsilon \left(a_L^\dagger a_L - a_R^\dagger a_R\right) - J \left(a_L^\dagger a_R + a_L a_R^\dagger\right) + \frac{U}{2}\left( a_L^\dagger a_L^\dagger a_L a_L + a_R^\dagger a_R^\dagger a_R a_R \right),$$7, the onsite term $$H=-\varepsilon \left(a_L^\dagger a_L - a_R^\dagger a_R\right) - J \left(a_L^\dagger a_R + a_L a_R^\dagger\right) + \frac{U}{2}\left( a_L^\dagger a_L^\dagger a_L a_L + a_R^\dagger a_R^\dagger a_R a_R \right),$$8 vanishes and the effective Fock-space lattice becomes exactly self-imaging because the spectrum is equally spaced. If one initializes the system at one boundary of Fock space, the imbalance oscillates between $$H=-\varepsilon \left(a_L^\dagger a_L - a_R^\dagger a_R\right) - J \left(a_L^\dagger a_R + a_L a_R^\dagger\right) + \frac{U}{2}\left( a_L^\dagger a_L^\dagger a_L a_L + a_R^\dagger a_R^\dagger a_R a_R \right),$$9 and $$\hat H=\frac{K}{8}(\hat N_1-\hat N_2)^2-\frac{\Delta\mu}{2}(\hat N_1-\hat N_2) -\frac{\mathcal E_J}{2}\left(\hat a_1^\dagger \hat a_2+\hat a_2^\dagger \hat a_1\right).$$0 with period

$$\hat H=\frac{K}{8}(\hat N_1-\hat N_2)^2-\frac{\Delta\mu}{2}(\hat N_1-\hat N_2) -\frac{\mathcal E_J}{2}\left(\hat a_1^\dagger \hat a_2+\hat a_2^\dagger \hat a_1\right).$$1

As $$\hat H=\frac{K}{8}(\hat N_1-\hat N_2)^2-\frac{\Delta\mu}{2}(\hat N_1-\hat N_2) -\frac{\mathcal E_J}{2}\left(\hat a_1^\dagger \hat a_2+\hat a_2^\dagger \hat a_1\right).$$2 increases, the oscillations are distorted and eventually become self-trapped: the probability distribution remains localized near the initially occupied side, corresponding physically to interaction-induced persistence of imbalance (Longhi, 2011).

The $$\hat H=\frac{K}{8}(\hat N_1-\hat N_2)^2-\frac{\Delta\mu}{2}(\hat N_1-\hat N_2) -\frac{\mathcal E_J}{2}\left(\hat a_1^\dagger \hat a_2+\hat a_2^\dagger \hat a_1\right).$$3 sector isolates correlated pair transport. In the basis $$\hat H=\frac{K}{8}(\hat N_1-\hat N_2)^2-\frac{\Delta\mu}{2}(\hat N_1-\hat N_2) -\frac{\mathcal E_J}{2}\left(\hat a_1^\dagger \hat a_2+\hat a_2^\dagger \hat a_1\right).$$4, $$\hat H=\frac{K}{8}(\hat N_1-\hat N_2)^2-\frac{\Delta\mu}{2}(\hat N_1-\hat N_2) -\frac{\mathcal E_J}{2}\left(\hat a_1^\dagger \hat a_2+\hat a_2^\dagger \hat a_1\right).$$5, $$\hat H=\frac{K}{8}(\hat N_1-\hat N_2)^2-\frac{\Delta\mu}{2}(\hat N_1-\hat N_2) -\frac{\mathcal E_J}{2}\left(\hat a_1^\dagger \hat a_2+\hat a_2^\dagger \hat a_1\right).$$6, the amplitudes satisfy

$$\hat H=\frac{K}{8}(\hat N_1-\hat N_2)^2-\frac{\Delta\mu}{2}(\hat N_1-\hat N_2) -\frac{\mathcal E_J}{2}\left(\hat a_1^\dagger \hat a_2+\hat a_2^\dagger \hat a_1\right).$$7

As $$\hat H=\frac{K}{8}(\hat N_1-\hat N_2)^2-\frac{\Delta\mu}{2}(\hat N_1-\hat N_2) -\frac{\mathcal E_J}{2}\left(\hat a_1^\dagger \hat a_2+\hat a_2^\dagger \hat a_1\right).$$8 increases, the intermediate state $$\hat H=\frac{K}{8}(\hat N_1-\hat N_2)^2-\frac{\Delta\mu}{2}(\hat N_1-\hat N_2) -\frac{\mathcal E_J}{2}\left(\hat a_1^\dagger \hat a_2+\hat a_2^\dagger \hat a_1\right).$$9 is suppressed and the bosons tunnel together as a pair. This is the simplest setting in which the two-site model already exhibits second-order pair tunnelling rather than merely independent single-particle transfer (Longhi, 2011).

In the current-algebra formulation, exact second-order Heisenberg equations distinguish parameter regimes particularly clearly. In the extreme Rabi regime $K_{ij}\hat N_i\hat N_j$0, and practically $K_{ij}\hat N_i\hat N_j$1, the symmetric case $K_{ij}\hat N_i\hat N_j$2 gives uncoupled simple-harmonic motion for $K_{ij}\hat N_i\hat N_j$3 and $K_{ij}\hat N_i\hat N_j$4 with frequency $K_{ij}\hat N_i\hat N_j$5, while $K_{ij}\hat N_i\hat N_j$6 is conserved and the period is $K_{ij}\hat N_i\hat N_j$7. When $K_{ij}\hat N_i\hat N_j$8, $K_{ij}\hat N_i\hat N_j$9 remains harmonic with frequency $\Omega_{ij}$0, but $\Omega_{ij}$1 and $\Omega_{ij}$2 become coupled and interfere (Filho, 2015).

A semiclassical analysis of the weak-coupling regime $\Omega_{ij}$3 and large $\Omega_{ij}$4 shows that the occupation difference exhibits rapid oscillations together with collapses and revivals. For the initial condition with all particles on one site, the collapse time, revival time, and broader blurring time are

$\Omega_{ij}$5

The first-order $\Omega_{ij}$6 correction sets the revival period, while the $\Omega_{ij}$7 terms control the detailed width and shape of the revival peaks (Veksler et al., 2014).

Extensions with nonlinear atom-pair tunnelling add further dynamical structure. In the repulsive, symmetric case of the extended Hamiltonian

$\Omega_{ij}$8

the semiclassical fixed-point analysis identifies Josephson, phase-locking, and self-trapping states, with phase boundaries $\Omega_{ij}$9, $N$00, and $N$01 in the $N$02 plane (Rubeni et al., 2016).

5. Ground states, fragmentation, cat states, and entanglement

For the attractive biased two-site Hamiltonian, exact diagonalization reveals a sequence of qualitatively distinct ground-state regimes as $N$03 becomes more negative. For weak attraction, $N$04, the ground state remains close to a single mean-field condensate and the Fock distribution is roughly binomial. For $N$05, a fragmented or cat-like window appears roughly in the interval $N$06: the ground and first excited states become nearly degenerate, the Fock-space decomposition develops two peaks, and the ground and first excited states are approximately symmetric and antisymmetric combinations of the same bimodal structure. For stronger attraction the external bias selects one well, $N$07, $N$08, and the cat-like superposition disappears. A central point is that fragmentation begins before a fully developed cat state appears (Melé-Messeguer et al., 2011).

To capture this crossover, an improved variational ansatz uses a superposition of two exchanged mean-field states,

$N$09

Because $N$10 and $N$11 are independent, the ansatz interpolates continuously between a single condensate, a nearly symmetric cat-like state, and a bias-selected localized state. In the pre-bifurcation regime it captures the onset of fragmentation missed by a single mean-field state, and in the cat-like window it reproduces the exact Fock distribution, imbalance, fluctuations, and one-body density-matrix eigenvalues substantially better than earlier fixed-amplitude cat ansätze (Melé-Messeguer et al., 2011).

A complementary analytical treatment uses a permanent-based interpolation controlled by the single-particle overlap $N$12. In the repulsive regime this gives explicit formulas for the energy, visibility, Fisher information, number fluctuations, decay rates, and entanglement entropy in terms of hypergeometric polynomials. In the attractive regime it replaces naive symmetry-broken mean-field states by a symmetric superposition of two imbalanced coherent-like states, yielding a finite-$N$13 Schrödinger-cat description. Within that framework the semiclassical symmetry-breaking point

$N$14

is valid only in the limit of infinite number of bosons, while for large but finite $N$15 the relevant crossover scale is

$N$16

This distinction addresses a recurrent misconception: at finite $N$17, the exact ground state remains symmetric even when the semiclassical equations admit broken-symmetry solutions (Dell'Anna, 2011).

The same permanent-based approach yields the relation

$N$18

between Fisher information and on-site number variance for both repulsive and attractive interactions. The entanglement entropy is maximal close to $N$19 and exceeds its coherent value throughout the interval $N$20. In the NOON limit the variance reaches $N$21 and $N$22, while in the coherent limit $N$23 (Dell'Anna, 2011).

6. Optical realization, open-system variants, and generalized dimers

A particularly direct realization of the two-site Bose-Hubbard dynamics is obtained in classical optics. The fixed-$N$24 Fock-space equations are exactly reproduced by light propagation in an engineered finite waveguide array of $N$25 sites,

$N$26

under the identification $N$27 and $N$28. Each waveguide represents one occupation-number partition, the spacing $N$29 controls the coupling $N$30, and the index contrast $N$31 controls the propagation-constant shift $N$32. The mapping is exact at the level of Fock-space amplitudes, so the optical power distribution $N$33 directly visualizes $N$34 (Longhi, 2011).

This photonic construction makes several many-body phenomena spatially visible: Josephson-like oscillations, self-imaging at $N$35, the onset of self-trapping as $N$36 increases, and correlated pair tunnelling in the three-waveguide $N$37 realization. Its limitation is equally precise: the analogy is classical and linear, so it does not reproduce measurement backaction, coupling to an environment, or genuinely nonclassical photon statistics. It is exact only for the linear coupled-mode equations associated with a fixed-$N$38 Bose-Hubbard sector (Longhi, 2011).

Generalized two-site models preserve much of the same algebraic structure. Rewriting the Hamiltonian

$N$39

in current variables produces coefficients $N$40, $N$41, $N$42, $N$43, and $N$44 that separate interaction nonlinearity, interaction asymmetry, bias, and collective energy shifts. In the special case $N$45, $N$46 is conserved and $N$47 and $N$48 are independent harmonic oscillators, recovering pure Rabi/Josephson oscillations. When $N$49 but $N$50, $N$51 remains an independent oscillator while $N$52 and $N$53 stay coupled (Filho, 2015).

Open-system and nonlinear variants extend the model beyond closed Hamiltonian evolution. In a two-mode system with pair damping, thermo field dynamics and Hartree-Fock decoupling reduce the master equation to an $N$54 squeezing problem. The short-time dynamics generates two-mode squeezing, hence entanglement, while the purity-like quantity $N$55 decays faster as the damping increases (Chaitanya et al., 2013). In the hopping-only dimer, evolution under $N$56 rather than $N$57 generates two-component Schrödinger cat states; at time $N$58, a coherent state becomes a superposition of coherent states with opposite amplitudes (Berezowski et al., 6 May 2026).

Taken together, these constructions show why the two-site Bose-Hubbard model retains a special status. It is small enough to be analytically tractable, but large enough to exhibit interaction-driven localization, macroscopic coherence, fragmentation, near-degenerate cat-like structure, integrable root geometry, nontrivial current dynamics, and exact analog simulation in synthetic lattices.