Strong-coupling ansatz for the one-dimensional Fermi gas in a harmonic potential.

A major challenge in modern physics is to accurately describe strongly interacting quantum many-body systems. One-dimensional systems provide fundamental insights because they are often amenable to exact methods. However, no exact solution is known for the experimentally relevant case of external confinement. We propose a powerful ansatz for the one-dimensional Fermi gas in a harmonic potential near the limit of infinite short-range repulsion. For the case of a single impurity in a Fermi sea, we show that our ansatz is indistinguishable from numerically exact results in both the few- and many-body limits. We furthermore derive an effective Heisenberg spin-chain model corresponding to our ansatz, valid for any spin-mixture, within which we obtain the impurity eigenstates analytically. In particular, the classical Pascal's triangle emerges in the expression for the ground-state wave function. As well as providing an important benchmark for strongly correlated physics, our results are relevant for emerging quantum technologies, where a precise knowledge of one-dimensional quantum states is paramount.

One-dimensional (1D) systems occupy a unique place in strongly correlated many-body physics, as many are exactly solvable via methods such as the Bethe ansatz. Likewise, the exactly solvable harmonic oscillator plays a central role in quantum mechanics. However, when one combines these two fundamental models and considers interacting fermions in a 1D harmonic potential, there is no known solution in general (). Whereas the problem may be solved analytically for 2 particles () and numerically up to ~10 particles (–), the calculation rapidly becomes untenable beyond that. Recent theoretical works have proposed analytic forms of the lowest energy wave functions near the Tonks-Girardeau (TG) limit of infinitely strong contact interactions (–), but these do not match exact numerical studies (–) once the particle number exceeds three. Here, we present a novel, highly accurate ansatz for the wave function of the two-component 1D Fermi gas in a harmonic potential with strong repulsive interactions.

Harmonically confined 1D systems have received a considerable amount of interest, particularly since the experimental realization in ultracold atomic Bose (, ) and Fermi gases (–). Experimentalists have trapped fermionic 6Li atoms in a 1D waveguide, with a high degree of control over both the number of particles in two hyperfine states and the interspecies interaction strength. This allows one to study the evolution from few to many particles, as well as the possibility of magnetic transitions as the interactions are tuned through the TG limit. The approach we propose here provides a way to tackle the regime near the TG limit, which has been investigated in a recent experiment (). In this case, the ground-state manifold in the confined system consists of nearly degenerate states, where N↑ (N↓) is the number of spin-↑ (spin-↓) particles. For the “impurity” problem consisting of a single ↓ particle (N↓ = 1) in a sea of N↑ majority particles, we generate all these states in an essentially combinatorial manner. We show that the overlap between our ansatz wave function and exact states obtained by numerical calculations exceeds 0.9997 for N↑ ≤ 8. In particular, the overlap with the exact ground state for large repulsive interactions is found to extrapolate to a value ~0.9999 as N↑ → ∞. This remarkable accuracy shows that our ansatz effectively solves the strongly interacting single ↓ problem in a harmonic potential, from the few- to the many-body limit.

We have furthermore mapped the strongly interacting 1D problem onto an effective Heisenberg spin chain of finite length (, , ), and we derive an analytical expression for the Hamiltonian within which our ansatz is exact. In this case, our ansatz for the impurity wave functions in the spin basis corresponds to discrete Chebyshev polynomials. In accordance with the orthogonality catastrophe (), we find that the overlap with the noninteracting many-body ground state (that is, the quasiparticle residue) tends to zero in the thermodynamic limit N↑ → ∞. However, surprisingly, in this limit, the ground-state probability distribution of the impurity is a Gaussian only slightly broadened compared with the noninteracting ground state. As we argue, our effective spin model is expected to accurately describe any N↓, N↑, opening up the possibility of addressing the strongly interacting 1D Fermi gas with powerful numerical methods for lattice systems. We also discuss how our results can be extended to higher excited states and how they may be probed in cold-atom experiments. Our findings are of fundamental relevance to emerging quantum technologies where an accurate knowledge of 1D quantum states is needed, such as for engineering efficient state transfer () and for understanding relaxation and thermalization in out-of-equilibrium quantum systems ().

RESULTS

The model

We consider N↑ fermions in spin state ↑ and N↓ fermions in spin state ↓, both with mass m, confined in a 1D harmonic potential. The total number of particles is written as N↑ + N↓ = N + 1, which is convenient when we consider the impurity problem below. The Hamiltonian is thuswhere the coupling g quantifies the strength of the short-range interactions and ω is the harmonic oscillator frequency. Note that particles with the same spin do not interact because their wave function vanishes when |x − x | → 0 due to antisymmetry under particle exchange. Because commutes with the total spin operator, the eigenstates have well-defined spin projection S = (N↑ − N↓)/2 and total spin S. In the following, we use harmonic oscillator units where ω = m = ℏ = 1.

In the TG limit, the coupling strength g → ∞ and the system simplifies significantly owing to the form of the boundary conditions when two particles approach each other. Specifically, for a given wave function ψ(x), the infinite repulsion requires , with x ≡ x − x, the relative coordinate for any pair of fermions with opposite spin and x = (x0, x1,…, x). Identical fermions always obey this condition, as mentioned above. Because all particles experience the same boundary conditions, the ground-state manifold for a system with fixed S contains degenerate states, corresponding to the number of unique configurations of ↑ and ↓ particles.

The simplest eigenstate of the Hamiltonian (Eq. 1) is the fully ferromagnetic state, corresponding to the maximum total spin S = (N + 1)/2. In this case, the spin part of the wave function is always symmetric regardless of S, and thus, the wave function in real space must be antisymmetric. Specifically, the wave function takes the form of a Slater determinant of single-particle harmonic oscillator wave functions and can be written ()with normalizationThe ferromagnetic state corresponds to that of N + 1 identical fermions, and thus, its energy is , where we subtract the center-of-mass zero point energy. Furthermore, it is an eigenstate for all g because the wave function antisymmetry guarantees that it vanishes when x → 0 so that it does not experience the particle-particle interaction.

For a given S (corresponding to fixed N↑ and N↓ ), the remaining eigenstates with the same energy E0 in the TG limit may be characterized by other values of S. For instance, for N↑ = N↓ = 1, there are two states, characterized by either S = 0 or S = 1. However, for general particle number, spin alone is not sufficient to determine the states with S < (N +1)/2 because the degeneracy of the ground-state manifold is , whereas the number of different S for a given S is 1 + min(N↑, N↓). Thus, to construct a unique orthogonal basis of eigenstates in the TG limit, we must consider how the states in the ground-state manifold evolve as g → ∞, as illustrated in Fig. 1. We mostly focus on the impurity problem where we have one ↓ particle at position x0 and N ↑ particles at positions x with 1 ≤ i ≤ N . In this case, we have N eigenstates with spin S = S = (N − 1)/2, in addition to the ferromagnetic state.

Fig. 1

Energy levels in the Tonks-Girardeau limit.

We display the exact energies (ψ | |ψ ) (red) and the result of our ansatz 〈〉 (blue dashed), given by E0 − in this limit, for the case of one ↓ particle and N = 3 ↑ particles. For g < 0, there is also a two-body bound state at negative energies (not shown).

Ground-state manifold in the TG limit

To construct the wave functions for the impurity problem with S = Sz = (N − 1)/2, it is useful to define a complete (but not orthogonal) set of basis functions involving ϕ0 = ψ0(x) and the N states:where s ≡ sign (x). For simplicity, we omit the dependence of ϕ on the coordinates. Each sign function simply replaces a zero-crossing in the Slater determinant (Eq. 2) with a cusp at the position where the impurity meets a majority (↑) particle (x = 0). As an example, for N = 2, we have basis functions:The basis functions are clearly degenerate with the ferromagnetic state (Eq. 2) when g → ∞ because the interaction energy vanishes, whereas the energy of motion in the harmonic potential is the same for all ϕ. This can be shown by noting that for any ordering of the particles (say x0 < x1 < … < x ) we have ϕ ∝ ψ0 (x). Thus, all eigenstates of the ground-state manifold in the TG limit must be linear combinations of the basis functions. Note that, alternatively, we could have chosen a basis set whose functions are non-zero only for a particular ordering of particles as in ().

The central question we address here concerns the nature of the eigenstates in the vicinity of the TG limit; that is, we wish to know the wave functions and energies perturbatively in the small parameter 1/g. This allows one to uniquely define the eigenstates at g → ∞ as being those that are adiabatically connected to the states at finite g. Before proceeding with degenerate perturbation theory, it is instructive to consider the structure of the exact eigenstates ψ (up to corrections of order 1/g) for N = 1, (), and N = 2, (),The subscripts on the wave functions order these in terms of decreasing energy for small but positive 1/g. Note that the eigenstates split into two orthogonal sets, which are even or odd with respect to parity because the Hamiltonian commutes with the parity operator. Referring to Fig. 1 and focusing on the repulsive case g > 0, we see that the ferromagnetic state ψ0 has the maximum energy within the manifold (), whereas the ground state ψ has the lowest total spin, that is, S = Sz, in accordance with the Lieb-Mattis theorem (). Physically, the cusps in the wave function for g → ∞ can easily be shifted from zero, which decreases the kinetic energy and thus leads to a lower energy compared with the ferromagnetic state. Indeed, we see two patterns emerging: ψ1 only contains states with one cusp, and only the ground state ψ contains the state with the maximal number of cusps. These observations suggest that the system may lower its energy by successively acquiring more cusps in the wave function.

Inspired by the above considerations, we now propose the following strong-coupling ansatz for the impurity eigenstates of the ground-state manifold in the vicinity of the TG limit:

For any N, the exact wave function ψ essentially corresponds to , a superposition of the basis functions ϕ restricted to k ≤ l.

In other words, the wave functions are obtained by a Gram-Schmidt orthogonalization scheme on the set of basis functions {ϕ}: is obtained by adding one cusp to ψ0, is obtained by adding one more cusp and then orthogonalizing it to , and so on. We emphasize that the ansatz allows one to obtain the entire ground-state manifold using linear algebra manipulations only. Thus, one can go far beyond the limit N ≲ 9 of current state-of-the-art calculations (). In fact, we will see that the ansatz allows us to obtain analytic expressions for all wave functions in the ground-state manifold for any N. We will show that the ansatz is remarkably accurate compared with exact numerical results, and that it allows one to calculate several observables analytically, even in the many-body limit.

The procedure for constructing our ansatz wave functions as outlined above can be performed straightforwardly even for large N, by noting that the inner products of the basis functions (Eq. 3), Φ ≡ 〈ϕ|ϕ〉, may be calculated combinatorially (see Methods).

Perturbation theory around the TG limit

To demonstrate the accuracy of our ansatz, we now turn to the explicit solution of the Schrödinger equation in the vicinity of the TG limit. Because here there are N + 1 degenerate states, we apply degenerate perturbation theory and obtain the ground-state manifold by means of finite-matrix diagonalization, in a similar manner to (, ). The energy can be written as E ≃ E0 − , where is the 1D contact density (–). From the Hellmann-Feynman theorem, we then obtainwhich defines the perturbation due to a non-zero 1/g. The state |Ψ〉 is a linear combination of the basis states {ϕ} of the ground-state manifold: |Ψ〉 = ∑α|ϕ〉. To obtain the eigenstates, we require |Ψ〉 to be a stationary state, that is, , resulting in the matrix equation , with the eigenvalue (contact density) of the state |Ψ〉. The matrix elements of are(see Methods).

For N = 1 (N = 2), the evaluation of is straightforward and yields , whereas all other elements vanish. Thus, we findwhere ≡ 〈ψ||ψ〉 is the contact coefficient corresponding to the state ψ. All the eigenstates for N ≤ 2 are uniquely determined by the two symmetries of parity and spin, so that the ratios of and the general structure of the wave functions in Eq. 4 hold for any confining potential that preserves parity and spin. However, these symmetries alone are not sufficient to determine the eigenstates for N > 2, and therefore N = 3 will provide a nontrivial test of our ansatz. In this case, the coefficients of and the eigenstates may still be evaluated analytically, but their form is sufficiently complicated that we relegate these to Methods. Converting long analytical expressions into numerical values for brevity, we obtainwhereas the contact coefficients areThese contact coefficients determine the energy splitting shown in Fig. 1 and they agree with those obtained in (). For N ≥ 4, we resort to a numerical evaluation of the matrix elements of , which may be calculated efficiently using a novel method outlined in Methods.

Our ansatz, however, is far simpler than the brute-force approach above. Whereas the evaluation of the multidimensional integral in Eq. 5 quickly becomes untenable as N increases above ~10, the implementation of our ansatz is a basic exercise in linear algebra: Applying our Gram-Schmidt orthogonalization scheme, we find the statesComparing Eq. 10 with Eq. 8, we see that our ansatz is extremely accurate, with only a minute deviation from the exact result for and (ψ0 and ψ2 are determined exactly from parity and spin). We note that our proposed wave functions are identical to those obtained numerically in (), illustrating that the results of our ansatz are essentially indistinguishable from numerical calculations.

We now explicitly demonstrate that the very high accuracy of our ansatz also holds for higher particle number N and that it even seems to hold in the many-body limit. A natural measure of its accuracy is the wave function overlap |〈ψ|〉| between the exact eigenstates ψ and our proposed ones . Writing the wave functions as , the overlap is simply . For the two nonexact states with N = 3 discussed above, we then find this quantity to be 0.999993, where we remind the reader that this is the numerical value of an analytic result (see Methods). Strikingly, we find that the overlap exceeds 0.9997 for all states up to N = 8, with the error being largest for the states “intermediate” between the ferromagnetic state ψ0 and the ground state ψ In Fig. 2, we illustrate how the wave function overlaps for the and always exceed 0.99994. In addition, the overlap in the ground state appears to extrapolate to a value ~0.9999 as N → ∞. Our ansatz is therefore essentially indistinguishable from “numerically exact” methods, even in the many-body limit. This shows that our ansatz effectively solves the strongly interacting 1D impurity problem for general N.

Fig. 2

Accuracy of the ansatz.

Overlaps between our ansatz and the exact wave functions ψ for majority particle numbers N ≤ 8. For the ferromagnetic state, this always equals 1 (black line). The red and blue dots depict the overlap for and (ground state), respectively. These are both 1 for N = 2 because these states are uniquely determined by spin and parity, whereas they are both 0.999993 for N = 3. The extrapolations (dashed lines) are least-squares fits of the data points to cubic polynomials. Inset: The wave function overlap of Girardeau’s proposed state () with the exact ground state.

Of particular interest is the state ψ, the ground state for repulsive interactions. Girardeau proposed () that this state is simply given by the state with the maximum number of cusps inserted, that is, ψ = ϕ. As shown in the inset of Fig. 2, the overlap of Girardeau’s proposed state with the exact ground state is 76% for N = 8, and it most likely tends to zero as N → ∞. Thus, our ansatz is a significant improvement compared to previous proposals for the ground-state wave function.

We now turn to the contact coefficients of the N + 1 energy levels in the ground-state manifold, that is, the splitting of the spectrum at finite coupling. In Fig. 3, we show how the energy takes the following approximate form:Comparing with Eqs. 6, 7, and 9, we see that this expression is exact for N = 1 and 2, whereas it holds to within 3.0% for N ≤ 8. We show in the next section that the spectrum given by Eq. 11 is intimately linked with an effective Heisenberg spin Hamiltonian within which our ansatz is exact.

Fig. 3

Contact coefficients of the ground-state manifold.

The contact coefficients (blue dots) of the exact eigenstates ψ in the Tonks-Girardeau regime, controlling the splitting of the energy levels at finite but large coupling. The gray lines represent the approximate relationship ; see Eq. 11.

Effective Heisenberg spin chain

We now discuss how the 1D problem can be mapped onto a Heisenberg spin model (, , ). This enables us to determine the states analytically and also allows us to generalize our ansatz for the impurity problem to any N↓.

In the limit g → ∞, the system consists of impenetrable particles because the wave function must vanish when two particles approach each other. Thus, if the particles are placed in a particular order, they should retain that ordering as long as the repulsion is infinite. This allows us to consider the system in the TG limit as a discrete lattice of finite length N + 1, where the particle furthest to the left is at site i = 0, the next particle is at site i = 1, and so on. A small but finite value of 1/g then allows neighboring particles to exchange position, introducing a nearest-neighbor spin interaction in the lattice picture. We can thus write the Hamiltonian in the lattice aswhere S is the spin operator at site i and J is the nearest neighbor exchange constant, which can in general depend on i (, ). Subtracting the constant in each term of the sum ensures that the ferromagnetic state has energy E0. The Hamiltonian is valid to linear order in 1/g and the general form holds for any external potential.

The couplings J in the Heisenberg model (Eq. 12) can be determined by considering the single ↓ impurity problem in a new basis of position states |↓〉 with 0 ≤ i ≤ N. The lattice position i corresponds to the position of the impurity relative to the N majority particles. The position states are orthonormal, 〈↓|↓〉 = δ, and can be related back to ϕ (see Methods). The perturbation may then be evaluated in the position basis of the impurity by inserting a complete set of eigenstates, yieldingThe matrix elements in Eq. 13 provide an explicit construction of the Heisenberg Hamiltonian, Eq. 12.

We now determine the Heisenberg Hamiltonian within which our strong-coupling ansatz for the eigenstates is exact. Proceeding via “reverse engineering,” we form the effective Hamiltonian by replacing ψ with our ansatz wave functions in Eq. 13. By inspection of the effective Hamiltonian for all N ≤ 100, we find that we must use the approximation from Eq. 11 in Eq. 13 to obtain a Hamiltonian restricted to nearest-neighbor interactions, and we obtain the couplingsThe exchange constants take the form of an inverted parabola and are thus reminiscent of the real space harmonic oscillator potential (see Fig. 4). The form of the coefficients means that the impurity at small positive 1/g may minimize its energy by occupying primarily the center of the spin chain while alternating the sign of the wave function on the different sites. Contrast this with the ferromagnetic state, which is a completely symmetric function of the impurity position . This is equivalent to the state obtained by applying the total spin lowering operator to the spin polarized state with S = (N + 1)/2. Note that this symmetric spin function corresponds to an antisymmetric wave function in real space.

Fig. 4

Exchange constants and ground state of the Heisenberg model.

Illustration of the nearest-neighbor exchange constants (Eq. 14) of the spin Hamiltonian (Eq. 12) for N = 100. We also show the ground-state wave function (Eq. 18) within our ansatz (green dots).

The Heisenberg model obtained from our ansatz is exact for N = 1 and N = 2, whereas it is approximate for larger N. In particular, for N = 3, our ansatz yieldswhich should be compared with the result obtained by using the exact eigenstates and energies in Eq. 13, yieldingThe error in the coefficients is thus less than 1%. Note that the values in Eq. 16 agree with those obtained in (). For larger N ≤ 8, we find that the error in the coefficients at the central sites remains ≲0.3%, whereas the error at the edges of the spin chain remains ≲5%. This shows that our ansatz is most accurate when the impurity is near the center of the harmonic potential, which is always the case for the ground-state wave function as we demonstrate below.

Our effective “harmonic” Heisenberg model allows us to determine the general solution for the single ↓ impurity within our ansatz analytically. We obtainfor the eigenstates in the ground-state manifold, where is a normalization constant. This result may be verified by direct application of the Hamiltonian (Eq. 12), and follows from the basis functions ϕ being discrete polynomials of the variable (i − N/2) of maximum order l in the spin chain. The Gram-Schmidt procedure of our ansatz then yields the orthonormal discrete polynomials with maximal order l in the variable (i − N/2). The functions in Eq. 17 are well known in the field of approximation theory as discrete Chebyshev polynomials—see, for example, (). The analytical form for the ansatz wave functions provides a simple solution to the Gram-Schmidt procedure for general N. In particular, the ground-state wave function is simply a (sign-alternating) Pascal’s triangle:Note that, in real space, this wave function does not change sign under the exchange of the impurity with a majority particle.

From the analytical expression (Eq. 18), we can determine the probability that the impurity is at position i relative to the majority particles in the ground state. We obtain . This prediction is dramatically different from the constant probability distribution P(i) = 1/(N + 1) predicted by Girardeau’s proposed ground state, which in the spin-chain model takes the form Indeed, we see that as N → ∞. Thus, ψ is inaccurate for the ground state in the harmonic potential.

Finally, we emphasize that the mapping to the effective Heisenberg model allows us to find solutions for any N↑ and N↓: one simply needs to calculate the eigenstates of the Hamiltonian (Eq. 12) with coefficients given by Eq. 14. For instance, in the case of N↑ = N↓ = 2, the ground-state manifold is spanned by the six states with i ≠ j. Within this basis, we find overlaps ≳0.99998 between exact and approximate eigenstates, and our results are in excellent agreement with the wave functions obtained numerically in (). Furthermore, within our ansatz, the contact coefficients of the six states take the form where is the contact coefficient from the (N↑, N↓) = (3, 1) problem—see Eq. 9. As expected, because the Hamiltonian commutes with the spin operator, the spectrum contains that of the single-impurity problem and, in accordance with the Lieb-Mattis theorem (), the ground state has S = 0.

Approaching the many-body limit

The fact that the wave function overlaps appear to extrapolate to a numerical value very close to 1 (see Fig. 2) indicates that our ansatz is also highly accurate in the many-body limit. Thus, we now investigate the limit N → ∞ for the impurity ground state (Eq. 18) at large repulsion. We focus on properties that depend on the impurity probability distribution in the bulk of the system.

The first such quantity is the contact coefficient, as shown in Fig. 5. We compare it with the expression corresponding to McGuire’s exact solution to the single impurity problem in free space () mapped onto the harmonically confined system using the local density approximation (). We see that our prediction for the contact appears to extrapolate to the Bethe ansatz result in the many-body limit, thus implying that the local density mapping is valid for the single-impurity ground-state energy. Indeed, this is consistent with the fact that the impurity ground-state wave function is confined to the central region of the trap (see Fig. 4) where the density of majority particles is highest.

Fig. 5

Contact of the ground state in the few- and many-body limit.

Contact coefficient of the ground state at small positive 1/g as a function of N. The dots are the analytical results for N ≤ 3 and numerical results for 4 ≤ N ≤ 8. We do not show a comparison between the ground-state contact and the perturbation evaluated within our approximate states, , because the relative error between these is less than 0.05% for N ≤ 8. The dashed line is McGuire’s free-space solution mapped to the harmonic potential using the local density approximation—see the discussion in the main text. Inset: The ground-state contact coefficient in units of N3/2 and plotted as a function of 1/N to illustrate the possible convergence to McGuire’s prediction (marked by a triangle). The dashed line is a cubic fit to our data. We also compare with the expectation value of Girardeau’s proposed ground state (green squares).

We next calculate the probability density of the impurity in real space, P(x0) = ∫dx1⋯dxN| ψ(x)|2. This is very complicated to evaluate for general N, but in the thermodynamic limit N → ∞, the probability distribution of the approximate ground-state wave function (Eq. 18) may be converted into P(x0). The distribution of majority particles is unaffected by the presence of the impurity in the thermodynamic limit, and according to the local density approximation, it is , where μ0 is the chemical potential at the center of the harmonic potential. This, in turn, yields . The lattice index i in the Heisenberg model corresponds to the number of majority particles to the left of the impurity; thus, it may be related to the position in real space via . Because , we can then write:where we have taken the central part of the harmonic potential with . Substituting this into Stirling’s approximation to the ground-state probability distribution, P(i) ≈ 2(πN)−1/2 exp[−(2i − N)2/N], finally yields the probability density of the impurity particle in the thermodynamic limit:Remarkably, upon tuning the system from the noninteracting ground state at g = 0+ with probability density , the impurity wave function has only spread out slightly in the TG limit. The broader distribution may be viewed as an increase in the harmonic oscillator length by a factor or as a decreased effective mass. Note that our predicted probability density is completely different from that of the Girardeau state, n(x0)/N, which equals that of the ferromagnetic state. In particular, our predicted distribution retains a width for any N, whereas for the ferromagnetic state the width is . The narrow width of the ansatz ground state compared with the length of the spin chain, in turn, implies that its overlap with the exact ground-state wave function could approach 1 in the limit of large N. For instance, the overlap with a Gaussian of the same width indeed converges to 1 in the thermodynamic limit.

The small change in the ground-state impurity probability distribution from the noninteracting to the TG limit appears to suggest that the wave function of the system is only weakly perturbed by infinite interactions. On the other hand, it is well known that the system encounters the orthogonality catastrophe in the thermodynamic limit, where the state of the system has no overlap with the noninteracting state (). To reconcile these points, it is necessary to consider the impact of the interactions on the majority fermions, which reshuffle the Fermi sea. This is embedded in the residue Z = |〈ψ|ψNI〉|2, that is, the squared overlap of the ground-state wave function with the noninteracting ground state at g = 0:We compute the residue using a numerical method similar to that outlined in Methods, and the result is shown in Fig. 6. By fitting, we find that the residue decreases with N as ~1/. Intriguingly, the same scaling of the residue with particle number was predicted for a massive impurity immersed in a 1D Fermi gas in uniform space ().

Fig. 6

Emergence of the orthogonality catastrophe.

Residue of the wave function ψ as a function of N. For N = 1 and N = 2, we find the analytic results 2/π and 81/(16π2), respectively. The dashed line is 0.89/ . We do not show a comparison between the residue of the ground state in our approximation scheme and that of the exact ground state because the relative error is less than 0.07% for N ≤ 7.

Higher energy manifolds and breathing modes

We have demonstrated that our ansatz is extremely accurate for the (N + 1)–dimensional ground-state manifold of the impurity problem to order 1/g. We now show how it can be extended to states in higher energy manifolds. It is known that the 2D version of the Hamiltonian (Eq. 1) is part of a “spectrum generating” SO(2,1) algebra connected with scale transformations x → x/λ (). In 1D, this symmetry is broken for a finite interaction strength because the scaling of the interaction, gδ(x) → λgδ(x), is different than in 2D. The key point, however, is that the SO(2,1) symmetry is recovered in the TG limit g → ∞. We can then use the technology developed for the SO(2,1) symmetry, suitably adapted to the 1D case. In particular, one can define an operator B so that if |n〉 is an eigenstate with energy E, then |n + 1〉 = B†|n〉 is an eigenstate with energy E = E + 2 (see Methods). The spectrum in the TG limit therefore consists of towers of states separated by twice the harmonic potential frequency, where B|0〉 = 0 for the lowest state in each tower.

Away from the TG limit, each level in these towers is shifted in energy, and Eq. 11 gives the energy shift of the ground-state manifold to a very good approximation. Each state in the ground-state manifold represents the lowest state in a separate tower of modes, and we now use the SO(2,1) algebra to calculate the energy shift of the excited states in each tower. Within first-order perturbation theory, the energy shift δE of the nth excited mode |n〉 away from its value in the TG limit is given byNote that this result is exact to order 1/g. The expectation values in Eq. 22 can be calculated using operator algebra only, once the so-called scaling dimension of is known. Using Eq. 5, we findwhere (x) = ϕ(x/λ)/λ(. It follows that the scaling dimension of is = 3 in 1D. One can now use Eq. 22 to express the energy shift of the state |n〉 as a function of the energy shift of all lower states in the tower (see Methods). The simplest case is the energy shift δE1 of the first excited breathing state |1〉 = B†|0〉, given bywhere δE0 is the energy shift of the N + 1 particle ground state away from the value E0 = N(N + 2)/2 in the TG limit. Equation 24 predicts that the energy shift of the first excited mode is larger than the shift of the state in the ground-state manifold. Physically, this means that the excited state energy approaches its noninteracting value faster than the ground state as one moves away from the TG limit. Used in combination with Eq. 11, Eq. 24 generalizes our ansatz for the spectrum in the single impurity problem to higher energy manifolds. However, we emphasize that our results for the excited manifolds only depend on the energy levels in the ground-state manifold and are not limited to the impurity problem.

We can compare the prediction of Eq. 24 with the exact solution to the two-body problem. In 1D, the exact two-body energies E are determined by the equation (). Close to the TG limit g → ∞, we have E = 3/2 + δE0 for a state in the ground-state manifold and E = 7/2 + δE1 for the first excited state in the tower, with δEi/E ≪ 1. Expanding the Γ functions yields δE1/δE0 = 3/2, which is identical to the result obtained from Eq. 24 when using E0 = 3/2. This demonstrates explicitly that Eq. 24, valid for any N, recovers the exact two-body theory close to the TG limit.

For large N, it immediately follows from Eq. 24 that the correction to the energy shift of the first excited manifold goes as 1/N2. Moreover, this holds for any n ≪ N (see Methods). Thus, in the thermodynamic limit, we find that the dynamic SO(2,1) symmetry extends to finite interactions, up to order 1/g.

The lowest breathing mode frequency has been measured in several atomic gas experiments (–). The results in this section can therefore be tested experimentally providing a sensitive probe of interactions in the few-body system.

We note that the approach described in this section is exact to lowest order in 1/g, and it is completely general. It would for instance be interesting to apply it to a system of 1D bosons close to the TG limit, where the frequency of the lowest breathing mode was recently calculated using a mapping to an effective fermionic Hamiltonian ().

Radio-frequency spectroscopy

Our results may be probed directly in cold atomic gases using radiofrequency (RF) spectroscopy, as already applied in the recent experiment (). Consider a homogeneous RF probe with frequency ωrf, which flips the impurity atom from the hyperfine state |a〉 to the state |b〉. Within linear response, the RF signal is proportional towhere Pi(Pf) is the probability of occupation of the initial |i〉 (final |f〉) many-body state, and ψσ is the field operator for the hyperfine state |σ〉.

Assume that the system is initially in a definite state |i〉 and that all final states are empty. There are two kinds of RF spectroscopy. In direct RF spectroscopy, a = ↓ and the impurity atom interacts with the ↑ atoms in the initial state, which belongs to the interacting many-body ground-state manifold, whereas the final hyperfine state |b〉 of the impurity atom does not interact with the majority atoms. There will then be a peak at ωrf = −E0 in the RF spectrum in the TG limit, and the reduction of the height of the peak from its noninteracting value gives the quasiparticle residue of the initial state. There will also be peaks at ωrf = −E0 + 2n with n = 1, 2,… because the initial interacting wave function has components in excited noninteracting states with the same parity. The shift of the peak position away from ωrf = −E0 gives the energy shift of the many-body ground state when 1/g > 0. In inverse RF spectroscopy, the initial state |a〉 of the impurity atom does not interact with the majority atoms, whereas the final state does with b = ↓ (). There will then be a peak at ωrf = E0 in the RF spectrum in the TG limit and the shift in position when 1/g > 0 again gives the many-body energy shift directly. The reduction of the height of the peak from its noninteracting value gives the quasiparticle residue. There will also be RF peaks at higher frequencies corresponding to flipping into the excited interacting states.

DISCUSSION

In this work, we investigated in detail the properties of a single impurity immersed in a Fermi sea of N majority particles near the TG limit. By comparing the properties with exact numerical results, we have demonstrated the impressive accuracy of our strong-coupling ansatz for arbitrary N. We have furthermore identified the effective Heisenberg Hamiltonian within which our ansatz is exact, and this has allowed us to analytically evaluate the entire ground-state manifold, yielding the discrete Chebyshev polynomials. In particular, the ground-state wave function from our ansatz at strong repulsion is a sign-alternating Pascal’s triangle in the spin chain. Because its overlap with the exact ground-state wave function extrapolates to a value ~0.9999 for N → ∞, we believe that Eq. 18 is essentially indistinguishable from the result of numerically exact approaches.

In addition to the static properties considered here, our ansatz provides a framework for investigating impurity dynamics in a harmonic potential because we have determined the entire spectrum of the ground-state manifold and associated excited states related via a scale transformation. The impurity dynamics in 1D gases have recently been investigated theoretically in (, ) and experimentally in (, ).

Our results also extend far beyond the single-impurity problem because the effective Heisenberg Hamiltonian (Eq. 12) accurately describes any number of ↑, ↓ particles in the strongly coupled regime, as explicitly demonstrated for N↑ = N↓ = 2. For larger N↑ and N↓, where the number of states grows dramatically, our Hamiltonian can be tackled with numerical tools developed for lattice systems, such as the density matrix renormalization group () or matrix product states (). Extending our approach beyond the two-component Fermi gas to other quantum mixtures would also enable us to address open problems in the context of quantum magnetism, such as the nature of correlations and dynamical quantum phase transitions. In particular, it would be interesting to investigate whether our ansatz may be extended to harmonically confined N-component fermions, a scenario that has relevance to SU(N) magnetism () and that has recently been experimentally realized ().

Finally, the exceptional accuracy of our simple ansatz and the suggestive form of the impurity spectrum, E − E0 ∝ l(l + 1)/g, lead us to speculate that our results are the manifestation of a hidden approximate symmetry. This raises the tantalizing possibility that the 1D Fermi gas in a harmonic potential becomes near-integrable in the strongly interacting limit.

METHODS

Manipulations of the basis functions

To calculate the overlaps of the basis functions ϕ, it is useful to introduce an alternative formulation of the problem. First, we note that, for a given ordering of particles, all basis functions ϕ (and consequently any superposition of these) are proportional to ψ0. We may then define the complete and orthonormal set of basis states:where 1 ≤ j ≤ N. {x} corresponds to any set of i spin-↑ particles, and the step function Θ is 1 if precisely i of the N majority particles are to the left of the impurity, and zero otherwise. Clearly, the states described by Eq. 25 do not overlap. To see that they are properly normalized, considerHere, we used the fact that the integral over the (normalized) ferromagnetic state ψ0 does not depend on the ordering of particles. Thus, restricting the integral to a particular ordering, the result is 1/(N + 1)!. Now, if i particles are to the left of the impurity, there are ways of choosing these, with i! × (N − i)! ways of ordering those to the left and right of the impurity. Gathering the terms, the basis states {|↓〉} are thus seen to form an orthonormal basis.

Recall now the definition of the basis states ϕ:with s ≡ sign(xj − x0). Assuming there are exactly i majority particles to the left of the impurity, we may evaluate the sum of sign functions:This simply counts how many combinations exist with k [l − k] majority particles involved in the sign functions out of the i [N − i] majority particles located to the left [right] of the impurity, respectively. Thus, we arrive at the projection of the two sets of basis states:The prefactor arises from the same arguments that led to Eq. 25.

The inner products Φ = 〈ϕ|ϕ〉 then simply follow from Eq. 29,We also note that the matrix Φ is bisymmetric, that is, Φ = Φ = Φ.

Numerical evaluation of .

The matrix elements of can be evaluated aswhere we have used the fact that . This quantity can be non-zero because of the presence of cusps in the basis functions. Note that for the ferromagnetic state and therefore . This implies that ϕ0 is an eigenstate with = 0, as expected.

We now show that the (N + 1)–dimensional integrals appearing in the matrix elements of , Eq. 30, may conveniently be expressed in terms of a Taylor expansion of a suitable function. We begin by noting that the bracketed term of the wave functiononly depends on the relative coordinates. It is then convenient to introduce the center-of-mass coordinate and the relative coordinates y0 = (N + 1)(x0 − xcm) and y = x (1 ≤ i ≤ N). The constraint may be enforced through the use of a δ-function, which in turn may be converted into an extra integral, yieldingThis allows us to decouple the center-of-mass coordinate. The integrand of Eq. 30 contains a factor of (following from the form of ψ0), which may be written as . Because this is the only term in which the center-of-mass coordinate appears in Eq. 30, xcm may be integrated out to give .

Next, consider the part of the integral involving both k and y0:where in the last step we have integrated by parts. Collecting all factors, we get an expression where the integrals over the y have been decoupled:where . We have made use of the fact that the integral is independent of i so that we can differentiate with respect to just y and multiply by N. We also used the relation . Denoting the quantity in square parentheses by h(k) and introducing its Taylor expansion around , we see that the desired matrix element is converted into a quickly convergent series, .

For the example of N = 3, the function appearing in square parentheses for the matrix element 11 iswhereas for the matrix element , it isIn both cases, we see how the integrals over the relative coordinates between the impurity and the majority particles separate. It is also easy to see that h11(k) = h33(k), as expected.

Analytic solution for N = 3

We start by considering the explicit solution of the Schrödinger equation perturbatively for 1/|g| ≪ 1. First, we note that the matrix elements of can be written aswhere we have the integralthat is, we integrate over i of the majority particles to the left of the impurity, and N − i − 1 to the right. Note that the integral only couples basis states with the same parity. One can also show that

For N = 3, it is possible to evaluate analytically, thus givingCombining this with the matrix of inner products Φ then allows us to solve the eigenvalue equation and determine and ψ analytically. The resulting expressions are rather cumbersome and yield the numerical values shown in the main text.

SO(2,1) symmetry and excited states

We define the scaling operator S(λ) asIt can be written as (, )The raising and lowering operators can then be defined as ()for the N + 1 particle state. Here, is twice the harmonic potential, and with the total momentum and the center-of-mass coordinate. One can show that, in the TG limit, from which it follows that if |n〉 is an eigenstate with energy E, then |n + 1〉 = B†|n〉 is an eigenstate with energy E = E + 2. The operator B† excites breathing modes, and the spectrum in the TG limit consists of towers of states separated by twice the harmonic potential frequency, where B|0〉 = 0 for the lowest state in each tower.

To calculate the scaling dimension, we take the derivative of Eq. 23 with respect to λ, setting λ = 1, from which it follows that the scaling dimension is . The calculation of the energy shift, Eq. 22, is now rather long and cumbersome but straightforward because all necessary commutators are known (). We obtainwhere we have defined . Here, E is the internal energy of the jth state, that is, the energy minus the zero point energy of the center-of-mass. The last term in Eq. 36 is zero for n = 1. Equation 36 gives the shift δE of the energy of the state |n〉 away from its value in the TG limit in terms of the energy shifts of the lower modes.