Observing anyonization of bosons in a quantum gas

Experiment

The experiment started with an interaction-tunable 3D BEC of 1.3 × 10^{5 133}Cs atoms⁶⁰ prepared in the lowest magnetic hyperfine state \(| F,m_\rmF\rangle =| 3,3\rangle \equiv | \uparrow \rangle \), held in a crossed-beam dipole trap and levitated against gravity by a magnetic field gradient. The BEC is in the Thomas–Fermi regime with the 3D s-wave scattering length a_↑↑ tuned to a_↑↑ ≈ 220 a₀, corresponding to an offset magnetic field of B = 21.24(1) G. A 2D optical lattice, generated by two retro-reflected laser beams propagating in orthogonal directions, was gradually ramped up in 500 ms to a potential depth of 30E_r, with E_r = π²ħ²/(2ma²) the photon recoil energy, cutting the 3D system into an array of 1D tubes that are oriented along the vertical direction, as sketched in Fig. 1d. Here a = λ/2 is the lattice spacing with λ = 1,064.5 nm the wavelength of the lattice light. The longitudinal trapping frequency in the 1D tubes was 25.6(3) Hz. The magnetic field was then ramped up adiabatically to B = 35.1 G, tuning a_↑↑ to a_↑↑ ≈ 750 a₀, setting the Lieb–Liniger interaction parameter γ_↑↑ = mg_↑/(ħ²ρ) ≈ 14, where ρ  = N/L ≈ 1.33 μm⁻¹ is the average 1D density, and L is the average system length. The nominal value of the Fermi wavevector is given by k_F = πρ. Here g_↑ ≈ 2ħω_⊥a_↑↑ (ref. ⁵⁹), and ω_⊥ is the transversal trap frequency. For these values, the 1D systems are deeply in the fermionized Tonks–Girardeau regime^39,61.

The impurities were Cs atoms that were transferred to the Zeeman substrate \(| 3,2\rangle \equiv | \downarrow \rangle \) using a short radio-frequency pulse. Power and duration were set such that, on average, one impurity per tube was created. The pulse duration (15 μs) was much shorter than the Fermi time (t_F = 120 μs), ensuring that the spatial profile of the impurity closely matched that of the host gas. The number of impurities varied across the atomic density distribution. Because our detection was sensitive only to the impurity atoms, tubes with no impurities were irrelevant. For tubes with two impurities, the momentum distribution was not expected to be exactly anyonic, but the deviation was small because the impurities were still the minority component. The 3D scattering length between the impurity and host atoms a_↑↓ also varied with B (Fig. 1e). At B = 35.1 G, the host–impurity Lieb–Liniger parameter γ_↑↓ took the value γ_↑↓ ≈ 9. The impurity atoms in \(| \downarrow \rangle \) experienced a smaller levitating force and would be accelerated by F_↓ = mg/3. Such a comparatively strong force would lead to a non-adiabatic time evolution⁵⁹, populating the continuous spectrum of the gapless quantum liquid and pulling the system away from its ground state (see below). To avoid this, we adiabatically turned on optical levitation in 100 ms. Specifically, a 1,064-nm Gaussian beam with a 1/e² waist of σ_z ≈ 210 μm, positioned σ_z/2 above the atoms, generated a nearly linear optical potential gradient. A laser power of approximately 10 W indiscriminately levitated the host and impurity atoms when the magnetic force was off. A tunable force F_↓ on the impurity atoms while still fully levitating the host atoms can then be generated by adjusting the fraction of optical versus magnetic levitation.

Role of finite force and finite interaction

Here we studied the role of the finite force and finite interaction in our system. In Extended Data Fig. 1a, we show n_↓(k) at a fixed total momentum ħQ ≈ ħk_F for two different values of force F_↓. For a strong force F_↓ = mg/3, the distribution n_↓(k) was skewed and had a peak at around k = k_F. By contrast, for a relatively small force F_↓ = mg/18, the distribution was more symmetric and flat-top, as expected for a fermionic distribution. The simulations on the basis of sBHM were in good agreement with our experimental data. The residual asymmetry in the theoretical curve is attributed to the finite F_↓. The deviation between theory and experiment mainly resulted from the inhomogeneities of the experimental system in view of the distribution of k_F values for different tubes. From the observed broadening of the momentum distribution of the system, we estimated an upper bound on the Fermi wavevector variation across the tubes. The root mean square width corresponds to 0.2 ħk_F, where k_F corresponds to the mean value. Next, we turn to the effect of finite interaction strength on the momentum distribution. In Extended Data Fig. 1b, we showed n_↓(k) at a fixed Q ≈ 0.5k_F for three different values of interaction strength γ_↑↓. Close to the non-interacting point γ_↑↓ ≈ 0, the distribution resembled a bosonic distribution peaked around k = 0.5k_F and did not show any skewness. As we increased the interaction strength to a moderate value of γ_↑↓ ≈ 3, the height of the peak decreased, and n_↓(k) broadened to the left. Only for a sufficiently strong interaction did the distribution start to agree with the prediction from AHM. This confirmed that strong interactions are crucial for the emergence of anyonic correlations in our system. Note that the peak in the measured n_↓(k) was broader than that in the AHM predictions for a single tube. This is again attributed to the effect of inhomogeneities. Note that γ_↑↑ also varied when γ_↑↓ was changed, but it always stayed above 3.

Exchange symmetry engineering

We now elaborate on the way in which the emergence of a spin wave in the system led to the appearance of anyonic correlations on the original particles, as expressed by equation (1). Owing to the phenomenon of spin–charge separation, the exchange symmetry of the spatial part is dictated by the exchange symmetry of the spin part of the wavefunction. To obtain an exchange phase of θ in the spatial wavefunction, we needed to have an exchange phase of −θ on the spin wavefunction. To describe the system, we used the bosonic version of the approach described in ref. ⁴⁰, where the spinful bosonic system is replaced by a spinless bosonic charge sector and a spin chain, describing the spin of each atom. The unitary pairwise spin-exchange operators \(\hat\mathcalE_\ell ,\ell ^\prime \) exchange spin ℓ with spin \(\ell ^\prime \) in the spin chain. The set of \(\hat\mathcalE\) operators generates the symmetric group of permutations S_N. A fully anyonic wavefunction should be a simultaneous eigenstate of all \(\hat\mathcalE_\ell ,\ell ^\prime \), with the eigenvalue \(e^-i\theta \rmsgn(\ell -\ell ^\prime )\).

This state cannot exist for various reasons. Because \(\hat\mathcalE^2\) is the identity operator, the eigenvalues of \(\hat\mathcalE\) are ±1, corresponding to triplet (bosonic) and singlet (fermionic) wavefunctions. Furthermore, two exchange operators of the type \(\hat\mathcalE_\ell ,\ell ^\prime \) and \(\hat\mathcalE_\ell ^\prime ,\ell ^\prime\prime \) do not commute with each other, as can easily be verified. Simultaneous eigenstates of all pairwise exchange operators are therefore not easy to find, as a result of the fact that the group S_N for N larger than 2 is non-abelian. Nevertheless, certain observables in the form of correlation functions can be sensitive only to a subgroup of exchanges, as shown in the following.

We now try to find the common eigenstates of only a subgroup of S_N, with the required form of eigenvalues. In this sense, although this method cannot generate a fully anyonic wavefunction of the host–impurity system, it can at least give us direct access to specific observables of the anyonic gas. We look for a subgroup of S_N with elements that can have complex eigenvalues. The cyclic subgroups are abelian, and the eigenvalues of the different elements are given by the mth roots of unity if m is the size of the cycle. We concentrate on the cyclic group of maximal order, C_N, because this is the most relevant for us. The group generator \(\widehatC\) performs a cyclic rotation of the spin-chain configuraion of the system \(\widehatC| \sigma _1,…,\sigma _N\rangle =| \sigma _N,\sigma _1,…,\sigma _N-1\rangle \). The eigenvalues are given by e^−iθ, for θ = 2πn/N, with n = 0, …, N − 1, and the eigenstates are spin waves. Let us clarify the connection between the exchange phase and the eigenvalue of \(\widehatC\). One cyclic permutation corresponds to N − 1 backward binary exchanges. This can be seen by inspecting the effect of the operator on the state of the spin chain. To reproduce the behaviour of anyons with forward exchange phase −θ, the eigenvalue of \(\widehatC\) should correspond to e^iθ(N−1). This reduces to e^−iθ using the condition θN = 2πn, with \(n\in \mathbbZ\), which is necessary to keep the wavefunction single-valued. The allowed values of θ are therefore discretized but become dense in the thermodynamic limit. In our experiment, N ≃ 37, giving a discretization in steps of Δθ/π ≃ 0.03, which is below our uncertainty owing to inhomogeneities. In the case of a single impurity, the spin waves take the form

We wanted to identify the correlation functions that are well described by the \(\widehatC\) operator. The simplest example is the one-body correlation function of the impurity, for the single-impurity case. To see this connection, consider the action of the operator \(\widehatb_\downarrow ^\dagger (x)\widehatb_\downarrow (\,y)\) on the spin configuration of the 1D system. The destruction operator is only non-zero if the spin-down particle is found at position y, and the creation operator then places it at position x. As a result, the spin configuration of the system is shifted by exactly the amount \(\widehatN(x)-\widehatN(y)\), taking x > y. Here \(\widehatN(x)=\int _-\infty ^x\widehatn(y)dy\) counts the number of particles to the left of x. This corresponds to the application of the operator \(\widehatC^\widehatN(x)-\widehatN(y)\). We can therefore rewrite

where \(\widehatb\) is the destruction operator of spinless hardcore bosons in the charge sector, and \(\widehat\Pi _\downarrow (\widehatN(\,y))\) is the projector operator on spin down for the spin at position \(\widehatN(\,y)\) in the spin chain. If the spin state \(| \theta \rangle \) is prepared, we get

where in the last equivalence we used the Jordan–Wigner transformation \(\widehata=\widehatbe^i\theta \widehatN\). The factor 1/N results from the mean value of \(\widehat\Pi _\downarrow \) on the spin wave. It is easy to see how this argument can be generalized to the multi-impurity case, giving a family of anyonic correlation functions that can be exactly simulated with this method. Their explicit expression is given by

where the number m of creation (destruction) operators should match the number of spin σ particles in the spin wave. This demonstrates how, whenever the spin-wave state is realized, we can find correlation functions of the original spinful gas that map exactly onto the correlation functions of a system of N anyons, explaining why it is possible to access the momentum distribution of the anyons with measurements on the original spinful bosons. Making use of this equivalence in practice requires control of the spin state of the system, but it is completely independent of the state in the charge sector. It is therefore possible to directly measure the dynamics of the anyonic correlation functions, assuming that the spin wavefunction remains in a spin-wave state during evolution. In our system, we prepared a spin wave as the eigenstate of momentum with the lowest energy by slowly accelerating the impurity.

Emergence of anyons via spin–charge separation

We now turn to a lattice model to understand how the charge sector can be mapped onto an anyonic gas when the spinful hardcore bosons are prepared in a finite-momentum ground state. We consider the Hamiltonian \(\widehatH_\rmlat\) describing a gas of N spinful hardcore bosons:

Here \(\widehatb_\sigma \ell ^\dagger \) (\(\widehatb_\sigma \ell \)) are bosonic creation (annihilation) operators at site ℓ, σ = (↑, ↓) is the spin index and J is the hopping amplitude, and its value is specified below. We assume to be in the low-density limit N/L_S ≪ 1, where L_S is the number of lattice sites, and we impose periodic boundary conditions so that the conservation of momentum is assured. The operators \(\widehatb_\sigma \ell ^\dagger \) (\(\widehatb_\sigma \ell \)) are assumed to satisfy a no-double-occupancy constraint, \(\sum _\sigma \widehatb_\sigma \ell ^\dagger \widehatb_\sigma \ell \le 1\). Under this no-double-occupancy constraint, the spin and charge degrees of freedom separate, that is, the wavefunction \(| \varPsi \rangle \) can be written as \(| \varPsi \rangle =| \varphi \rangle \otimes | \chi \rangle \). Here \(| \varphi \rangle \) and \(| \chi \rangle \) denote the wavefunction for the charge and spin parts, respectively. The Hamiltonian \(\widehatH_\rmlat\) can be written in spin–charge separated form as⁶²

where \(\widehatf_\ell ^\dagger \) (\(\widehatf_\ell \)) is the spinless fermionic creation (annihilation) operator at site j, and \(\widehatC\) is the previously introduced spin permutation operator. Note that a bosonic description of the charge sector, with hardcore constraint, is also possible but has the disadvantage that the bosonic particles are still interacting so that diagonalization is not straightforward. The spin permutation operator \(\widehatC\) and the spinless fermionic operators can be diagonalized separately because they are independent of each other. The eigenstates of \(\widehatC\) are spin waves of the form

where \(| \sigma _1,…,\sigma _N\rangle \) is an arbitrary configuration of the spin chain, ν enumerates all possible disconnected spin blocks and N_ν corresponds to the total number of distinct elements of the form \(\widehatC^j| \sigma _1,…,\sigma _N\rangle \) in the νth block. The eigenvalues of \(\widehatC\) are given by e^−iθ, for θ = 2πn/N_ν, with n = 0, …, N_ν − 1. In the case of a single impurity N_↓ = 1, the eigenstates take the form of equation (3). By projecting \(\widehatH_\rmsc\) on the eigenspace of \(\widehatC\), we get an effective Hamiltonian for the charge sector

Here we see that the fermionic charge sector acquires an overall flux. This spin-generated flux is a collective effect, imposed by the spin waves onto the charge degrees of freedom. Note that the original Hamiltonian \(\widehatH_\rmlat\) does not break time-reversal symmetry. However, time-reversal symmetry is broken for the \(\widehatH_\rmeff\) governing the charge sector. This is a result of the projection onto a specific spin-wave subspace. Finally, we performed an anyonic transformation

The phase factor in the boundary term vanishes, (−1)^N−1e^iθe^{i(θ+π)(N−1)} = 1, and the Hamiltonian \(\widehatH_\rmeff\) can be mapped onto a system of hardcore anyons with a periodic boundary condition

As one can see, the anyonic model does not contain any concatenated flux. The transformation equation (11) is a generalized Jordan–Wigner transformation⁶³, and the anyons can be understood as composite particles in the charge sector^64,65. Each spin wave selects a specific value for the statistical phase. In the thermodynamic limit, this result also holds for any choice of boundary conditions. This justifies the use of fixed boundary conditions in the numerics.

Next, we turn to anyonic observables that can be measured experimentally. The real-space density of these anyons can be extracted by measuring the total density of the gas \(\langle \varphi | \widehata_\ell ^\dagger \widehata_\ell | \varphi \rangle =\langle \varphi | \,\widehatf_\ell ^\dagger \widehatf_\ell | \varphi \rangle \,=\) \(\sum _\sigma \langle \varPsi | \widehatb_\sigma \ell ^\dagger \widehatb_\sigma \ell | \varPsi \rangle \), where Ψ is the many-body wavefunction of the whole system. However, for hardcore anyons, the real-space density is independent of θ. The one-body correlator \(\langle \widehata_i^\dagger \widehata_j\rangle \), on the other hand, is very sensitive to θ. The Fourier transform of this gives the anyonic momentum distribution, which can be measured by measuring the momentum distribution of the impurity in our system through equation (1). Note that Hamiltonian (equation (10)) can be diagonalized exactly⁶². The momenta of the fermions correspond to the rapidities of the system.

Anyon Hubbard model

To benchmark the anyonic behaviour realized in the experiment, we next elaborated on the anyonic correlations of the paradigmatic AHM, which can be effectively simulated by using a bosonic model with density-dependent tunnelling. By using a fractional version of the Jordan–Wigner transformation, that is, the anyon–boson mapping

AHM from equation (2) can be expressed in terms of bosonic operators as

Here \(\widehatb_\ell \) are the bosonic annihilation operators at site ℓ.

Different from the bosonic one-body density correlation \(\langle \widehatb_\ell ^\dagger \widehatb_\ell ^\prime \rangle \), the correlator of anyons \(\langle \widehata_\ell ^\dagger \widehata_\ell ^\prime \rangle \) can be expressed as

For the data shown in Figs. 2 and 3, we have assumed N = 10 anyons in L_S = 40 lattice sites in the hardcore limit with open boundary condition. The effect of the boundary condition is negligible for large system sizes. The use of a reduced atom number speeds up considerably the numerics. We found a consistent and satisfactory agreement with the experimental data by considering a large system size at low filling (Supplementary Information).

Dynamical evolution with sBHM

In practice, a spin wave can be generated by slowly accelerating the impurity. To efficiently simulate such a dynamical process, we considered an sBHM on a 1D lattice:

Here \(\widehatb_\uparrow \ell \) and \(\widehatb_\downarrow \ell \) are the annihilation operators of the host particles and an impurity at site ℓ, respectively, with their hopping strength being denoted by J. We consider the hardcore limit of the intra-component interaction, that is, U_↑↑ → ∞ and U_↓↓ → ∞. The on-site interaction between the host particles and the impurity is denoted by U_↑↓. A constant force F_↓ is applied only to the impurity. We define the dimensionless force \(\mathcalF=\fracF_\downarrow m\hbar ^2\rho ^3\). At the low filling limit, any lattice model reduces to a continuum model with the effective mass given by

By setting the value of the effective mass to be equal to the particle’s mass, we fix the value of Ja². By defining the filling factor in a lattice n = N/L_S and a = L/L_S being the lattice constant, one obtains the following mapping between quantities:

In our simulation, the initial impurity distribution was defined by the ground state of the Hamiltonian (equation (16)) with F_↓ = 0 and a harmonic trapping potential V applied only for the impurity. At t = 0, we suddenly removed the traps and switched on the constant force F_↓. We simulated the quench dynamics by solving the time-dependent Schrödinger equation associated with the Hamiltonian (equation (16)) by using the time-dependent variational principle on the basis of matrix product states implemented using ITensors^45,66. The results are presented in Figs. 2 and 3. The parameters chosen were L_S = 40, N_↓ = 1, N_↑ = 20, U/J = 9.1 and F_↓a/J = 0.15 for numerical convenience. A more costly simulation by using a larger system size (for example, L_S = 120) at lower filling (for example, N_↑/L_S = 0.25) gives very similar results (Supplementary Information).

Swap model

Inspired by the central role of the spin wave in the emergence of anyonic behaviour of our system, we developed a toy model with a ground state that encodes the spin wave we are targeting:

with \(\widehatb_\uparrow \ell \) and \(\widehatb_\downarrow \ell \) being the annihilation operators of the host particles and the impurity at site ℓ, respectively, and their hopping strength is denoted by J. In the strongly interacting regime, the swapping strength J_ex is expected to be of the order of J²/U_↑↓. We encoded the spin-wave information by assigning the factor ^eiθ to the swapping terms. The ground state of the swap model is expected to effectively describe the lowest energy state of the spinful system for momentum ħQ = ħρθ (ref. ⁶⁷).

In a spin–charge separated representation, the one-body correlation function \(\langle \widehatb_\downarrow \ell ^\dagger \widehatb_\downarrow \ell ^\prime \rangle \) of the single impurity can be implemented by hopping of spinless particles and swapping \(\widehat\mathcalE\) on the spin chain^40,68. Taking \(\ell ^\prime \ge \ell \) as an example, we have

Here the Kronecker δ operators ensure that the \(\ell ^\prime \rmth\) site is occupied by the \(m^\prime \rmth\) spin, and after hopping, the ℓth site is occupied by the mth spin. In the case of a single impurity, the product of swap operators is related to the \(\widehatC\) operator, as shown previously, giving rise to a spin wave, which leads to the one-body correlator of the impurity shown in equation (1). Calculating the Fourier transform of the one-body correlator of the impurity and using the parameters L_S = 120, N_↓ = 1, N_↑ = 30 and J_ex/J = 0.01, we obtained the quasi-momentum distribution of the impurity shown in Figs. 2 and 3. Note that the small value of swapping strength J_ex is related to the strong host–impurity interaction, and the agreement with experimental data is found for a wide parameter regime (Supplementary Information).

Rapidity of anyons in one dimension

In Fig. 4d–f, we present the results of the simulation of the quench dynamics of anyonic gases after suddenly removing the harmonic trap in one dimension. The momentum distribution as a function of evolution time is expressed as

with the single-particle density matrix of hardcore anyons ρ_HCA(x, y; t). Following ref. ¹⁵, it can be efficiently computed as

where A_mn(x, y; t) are the matrix elements of \(\bfA(x,y;t)=(\bfP^-1)^T\det \bfP\), and the elements of matrix P(x, y; t) are \(P_mn(x,y;t)\,=\,\delta _mn\,-\) \((1-e^-i\theta \rmsgn(y-x))\rmsgn(y-x)\int _x^ydz\phi _m^* (z,t)\phi _n(z,t)\). Here, ϕ_n(x, 0) are the single-particle wavefunctions of the 1D harmonic oscillator, and ϕ_n(x, t) fulfill the time-dependent Schrödinger equation

with Heaviside step function Θ(t), which models a sudden quench ω(t) = ω₀Θ(−t). The solution was found to be \(\phi _n(x,t)\,=\) \(\phi _n(x/b(t),0)e^imx^2\mathopb\limits^./2b\hbar -iE_n\tau (t)/\hbar /\sqrtb(t)\), with the scaling factor \(b(t)=\sqrt1+\omega _^2t^2\), \(\tau (t)=\int _^tdt^\prime /b^2(t^\prime )\) and E_n = ħω₀(n + 1/2). In the experiment, the trapping frequency was set to ω₀ = 2π × 25.6(3) Hz, and the average Fermi time was \(t_\rmF=2m/\hbar k_\rmF^2\approx 0.12\) ms. Owing to the finite size of the optical levitation beam, the expansion time t_1D was limited to about 5 ms in the experiment.