Strongly interacting confined quantum systems in one dimension

Abstract

In one dimension, the study of magnetism dates back to the dawn of quantum mechanics when Bethe solved the famous Heisenberg model that describes quantum behaviour in magnetic systems. In the last decade, one-dimensional (1D) systems have become a forefront area of research driven by the realization of the Tonks-Girardeau gas using cold atomic gases. Here we prove that 1D fermionic and bosonic systems with strong short-range interactions are solvable in arbitrary confining geometries by introducing a new energy-functional technique and obtaining the full spectrum of energies and eigenstates. As a first application, we calculate spatial correlations and show how both ferro- and antiferromagnetic states are present already for small system sizes that are prepared and studied in current experiments. Our work demonstrates the enormous potential for quantum manipulation of magnetic correlations at the microscopic scale.

ARTICLE

Received 1 Sep 2014 | Accepted 18 Sep 2014 | Published 4 Nov 2014

DOI: 10.1038/ncomms6300

Strongly interacting conned quantum systems in one dimension

A.G. Volosniev1, D.V. Fedorov1, A.S. Jensen1, M. Valiente2 & N.T. Zinner1

In one dimension, the study of magnetism dates back to the dawn of quantum mechanics when Bethe solved the famous Heisenberg model that describes quantum behaviour in magnetic systems. In the last decade, one-dimensional (1D) systems have become a forefront area of research driven by the realization of the TonksGirardeau gas using cold atomic gases. Here we prove that 1D fermionic and bosonic systems with strong short-range interactions are solvable in arbitrary conning geometries by introducing a new energy-functional technique and obtaining the full spectrum of energies and eigenstates. As a rst application, we calculate spatial correlations and show how both ferro- and antiferromagnetic states are present already for small system sizes that are prepared and studied in current experiments. Our work demonstrates the enormous potential for quantum manipulation of magnetic correlations at the microscopic scale.

1 Department of Physics and Astronomy, Aarhus University, Ny Munkegade 120, Aarhus C DK-8000, Denmark. 2 SUPA, Institute of Photonics and Quantum Sciences, Heriot-Watt University, Edinburgh EH14 4AS, UK. Correspondence and requests for materials should be addressed to N.T.Z. (email: mailto:[email protected]

Web End [email protected] ).

NATURE COMMUNICATIONS | 5:5300 | DOI: 10.1038/ncomms6300 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

& 2014 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms6300

Strongly interacting quantum particles are ubiquitous in nature and play a vital role in superuidity, superconductivity and magnetism, and low-dimensional magnetic

systems have great potential to deliver the key insights into fundamental properties of the materials used in modern technology1. The study of magnetism in one-dimension (1D) goes back to the dawn of quantum mechanics when Bethe solved the Heisenberg model by introducing the famous Bethe ansatz method2. During the 1960s the technique was used to solve several paradigmatic models including the LiebLiniger model of repulsive bosons3, impurity problems with fermions by McGuire4,5 and Yangs solutions of homogeneous two-component Fermi systems6 that led to Lieb and Wus solution of the 1D Hubbard model in 1968 (ref. 7). Common to the theoretical models is that periodic or open boundaries are essential ingredients, and the Bethe ansatz cannot be applied for general connement, which is nevertheless the reality in many state-of-the-art set-ups that have for instance been used to realize the TonksGirardeau8,9 gas using cold bosonic atoms1012.

We present a new functional method that is specically designed to overcome these problems and include the external conning geometry exactly in the strongly interacting regime to linear order in the inverse interaction strength. Furthermore, it can be used to interpolate from few to mesoscopic particle numbers and address how spatial correlations emerge and evolve1316. Our basic example below is a four-body system of two spin up and two spin down fermions in an external trap that is solved exactly for the rst time here. Direct access to the wave function allows us to see ferro- and antiferromagnetic correlations in the eigenstates and give exact probabilities for these congurations. This is an extremely important feature of this new approach as multicomponent systems become highly degenerate for strong interactions, and thus the energy itself yields little information about the system. Furthermore, we solve exactly the impurity or polaron problem of one spin down interacting strongly with a number of spin up particles in a harmonic trap, a set-up that has been realized experimentally16. This allows us to compare different conning potentials and show that correlations are strongly inuenced by the geometry. This has ramications on density functional approaches for strong interactions, and our scheme can provide invaluable benchmarks of procedures where Bethe ansatz solutions obtained with periodic boundary conditions are supplemented by the local density approximation. For nite particle numbers, periodic boundary conditions is a strong assumption whose justication can now be addressed using exact solutions. Ultimately, our method could realize the goal of answering the key question: how many particles does it take to make a strongly interacting many-body system?

ResultsDerivation of the energy functional. The general system we consider has N particles of mass m with coordinates x1,y, xN and

is described by the Hamiltonian

H Xip2i=2m Vxi

p

=mo

, and energy, E o. Here o is the oscillator

angular frequency and h/2p is the reduced Plancks constant.

For the hard wall, we have l L and E 2/2mL2. Henceforth, we

will measure lengths, energies and g in these units. Our focus will be on fermions with two internal spin states, up and down. Systems with more than two internal states can be addressed in a similar manner (see Methods).

A general eigenstate wave function has the form C(x1,y, xN).

For simplicity, we omit the coordinates from now on. The zero-range interaction implies that C obeys the boundary conditions

1

2g

" #

@C

@xj

@C

@xi

@C

@xi

@C

@xj

Cxi xj; 2

where the subscripts indicate the limits xi xj-0, that is,

they are derivatives from each side of the point xi xj. In the limit

where 1/g-0, the boundary conditions and the Pauli principle imply that C must vanish whenever xi xj for any i and j. Such

functions can be constructed from the eigenstates of the single-particle Hamiltonian (the rst term in equation (1)) by taking the antisymmetrized product of N states. This state we denote CA. Its

energy, EA, is a sum of the occupied single-particle energies. However, the boundary conditions allow us to write a more general state on the basis of CA9,2022,

C

XN !k1akyxPk1; . . . ; xPkNCAx1; . . . ; xN; 3

where we sum over the N! permutations, Pk, of the N coordinates, and y(x1,y,xN) 1 when x1ox2oyoxN and zero otherwise.

Fortunately, symmetries reduce the number of independent ak coefcients. For the present case of two-component fermions, the Pauli principle dictates that there are only M N!/(Nm!Nk!)

degrees of freedom. This is the number of degenerate states at 1/g 0, which shows that the functions in equation (3) constitute

a basis. The basic idea is now very simple. To linear order in 1/g we can write E EA K/g, where K K(a1,y,aM) is a functional

of the ak coefcients and is independent of g by the Hellmann Feynman theorem. One can now prove that K has the simple form (see Methods)

K

Xi4jdxi xj; 1

where pi is the momentum operator of particle i and V(xi) is an external conning potential. We assume a short-range two-body interaction that we model by a Dirac delta function of strength g. In the following, we are interested in the strongly interacting limit where g-N (or 1/g-0). For simplicity, our focus will be the repulsive case (g40), although our results can be extended linearly to the attractive side (go0) of 1/g 0. The deeply bound

states for g - N (refs 17,18) are irrelevant for our arguments

and will not be addressed here. The external potential produces an energy scale, E, and a length scale, l, in which we will express all other quantities. In the examples below, we will consider a double-well potential, V(x) mo2(|x| b)2/2, with barrier

parameter, b, and a hard wall potential of length L that vanishes for 0oxoL and has innite strength for xo0 and x4L. A double-well potential has recently been realized experimentally15,19. For b 0, the double well reduces to the

harmonic oscillator potential from which we adopt our units of length, l

Pk;p ak;pak ap2

Pk a2k; 4

where k and p run from 1 to M, and ak,p are matrix elements that depend only on the single-particle potential, V(x). The eigenfunctions and eigenenergies to order 1/g can now be obtained by variation of K with respect to ak and diagonalizing the resulting matrix. K is equivalent to Tans contact parameter23 in 1D24 and we compute it exactly for 1/g-0. Furthermore, the derivation can be easily extended to multicomponent bosons, fermions or mixtures, and it provides an effective Hamiltonian that can be used to study perturbations in the strongly interacting limit (see Methods for details).

Remarkably, in the strongly interacting regime, the effective Hamiltonian can always be written as a spin model. In the important case of N spin 1/2 fermions or two-component bosons

g

2 NATURE COMMUNICATIONS | 5:5300 | DOI: 10.1038/ncomms6300 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications & 2014 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms6300 ARTICLE

governed by the Hamiltonian in equation (1), it is a Heisenberg model of the form

PN 1i1 JiSi Si1 (ref. 25,26), which is a favourite starting point for research into quantum magnetism27.

Here Si is the spin operator of particle i. It has been shown20 that for a half-lled Hubbard model, Ji J is constant. This is obtained

by using the Bethe ansatz and by the same method one can also prove that Ji J for particles in hard wall (box) connement28.

Our approach generalizes these important results and not only do we nd that Ji depends on the external connement but also provide a procedure for computing these nearest-neighbour interaction coefcients exactly. We may therefore use the external connement to tailor the Ji coefcients into desirable spin models and thus manipulate static and dynamic quantum magnetic correlations26. Note also that this is not only true for the ground state manifold but also for higher manifolds as illustrated in Fig. 1a. Each manifold will have its own unique set of Ji coefcients that we can compute exactly. In the language of the Hubbard model, one can think of higher manifolds as belonging to higher bands.

Four-body systems. A central example is the hitherto unsolved four-body problem since it illustrates the method and it allows us to address magnetic correlation physics. We take Nm Nk 2

with an M 4!/(2!)2 sixfold degeneracy at 1/g 0 as shown in

Fig. 1a for the case of a harmonic trap (b=0). In general, the spectrum around 1/g-0 has the form of a ladder of manifolds each of which contains an M-fold fan of states as illustrated in Fig. 1a. For the ground state manifold we also show the adiabatic connection of states from weak to strong coupling where the third and fourth excited states are initially degenerate at g 0. The

parity invariance of the double well and hard wall potentials means that the three types of spatially correlated states shown in Fig. 1b,c, ferromagnetic, antiferromagnetic and mixed, completely

specify all solutions at 1/g 0. Figure 1b shows the conguration

probabilities for C1 and the state C4 for the hard wall, which turns out to have exact opposite ferro- and antiferromagnetic probabilities compared with the ground state. In Fig. 1d we show the double-well probabilities as functions of b for C1 and C3, again picked as examples because of their signicantly different correlations. In both cases we nd a ground state that is dominantly spatially antiferromagnetic, and perhaps more remarkably we nd excited states that are dominantly spatially ferromagnetic. Preparing different states at g 0 and then tuning

to 1/g 0 (refs 14,16) would thus produce completely different

correlation patterns. Note that if one considers two-component bosons instead of fermions the results are very different (see Methods).

One can understand intuitively what is going on by looking at equation (4). The antiferromagnetic conguration is favoured since (ai aj)2 is large for ai and aj differing in both sign and

magnitude. The functional approach presented here thus provides a very precise mathematical insight into the preference for domain walls of opposite spin in the strongly interacting regime and provides a spatial explanation of antiferromagnetism in repulsively interacting 1D systems. Moreover, our results also demonstrate the potential for manipulating correlations by state preparation and trap shape modulation. A step in this direction was recently reported using anisotropic optical lattices15. Most often one discusses ferromagnetism induced by symmetry breaking. We clearly have the presence of a degenerate manifold of states to induce such breaking and a small spin gradient is enough to drive the system into a purely ferromagnetic state. However, our direct access to the exact wave function demonstrates the presence of intrinsic magnetic correlations even without breaking the spin symmetry. More generally, our method can be used to study the correlations that drive quantum phase transitions in larger systems using exact wave functions.

d

0.8

1, (AF)

1, (F) 3, (AF)

1, (M)

3, (F)

3, (M)

Excited state manifold

E1

E

0.6

Ground state manifold

P 0.4

A

E0

5

3, 4

2

1

0.2

1/g 0

Repulsive Attractive

0 0.5 1 1.5 2

b

Hard wall trap AF

0.333

Mixed

Double well

trap

Antiferro

magnetic

F

Ferro

magnetic

0.045 Mixed

P(1)

P(4)

0.622

0.045 0.333 0.622

Figure 1 | Four-body system with two spin up and two spin down particles. (a) Schematic spectrum of eigenstates showing the ground state and an excited state manifold. The slopes at innite g will generally be different around E0 and E1. For the ground state manifold, we indicate the adiabatic connection between the strongly and weakly interacting regime for a harmonic trap (b 0). The red lines are positive, while the green lines are negative

parity states. The structure around 1/g 0 is the same for both double-well and hard wall traps, although the slopes are different. (b) Hard wall trap

congurations and their probabilities for C1 and C4 that have opposite antiferromagnetic (AF) and ferromagnetic (F) contributions. The probabilities include both the congurations shown and their parity inverted partners that are equal. (c) Same as b for the double-well trap. (d) The probability to nd the three congurations in c as function of the barrier parameter, b, for the ground state, C1, and an excited state, C3.

NATURE COMMUNICATIONS | 5:5300 | DOI: 10.1038/ncomms6300 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

& 2014 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms6300

Impurity problems. As another demonstration of the nature of the strong coupling regime, we consider the case of a single spin down (impurity) interacting with a variable number of spin up fermions. An impressive recent experiment has considered this system for Nr6 (ref. 16). As the energy at 1/g-0 is degenerate, further insight into the strongly interacting regime has to come from correlations in the N-body wave functions. Here we consider the probability for the impurity to tunnel out of the trap as shown in Fig. 2a. In a simple model, we assume that because of the strong repulsion only the particle on the far right can tunnel out as the barrier is lowered (see Methods for details). Since the wave function in equation (3) contains a superposition of states with the impurity in different positions, the probability is simply given by the amplitude for it to be on the right, aN. In

Fig. 2b we plot the impurity tunnelling probability, Pk |aN|2, for

N 2,y, 10 for the ground state of both a harmonic trap and a

hard wall potential (see Methods for details). The results for N 2 and 3 are trap-independent, while for NZ4 we see a clear

geometrical dependence. In particular, we nd that the scaling with N is completely different; whereas by a t we nd approximately PkpN 3 for the hard wall, the harmonic trap is not a power law but rather closer to a 1/N! behaviour as seen in the inset of Fig. 2b. The exact results thus allow us to conclude that geometry has a strong effect on correlations in the wave function. The exact wave functions also show that the impurity has a peak in its probability density at the centre of the trap. This is already true for N 3 and shows that for strongly interacting

systems, studying the few-body limit gives insight into the behaviour of larger systems. Also note that the combination of McGuires solution (using periodic boundary conditions)4 and the local density in the trap29 can only capture the energy in the strongly interacting regime but does not reproduce the energy slope to order 1/g. The 1/N line in Fig. 1b applies to the

non-interacting state CA (which is often referred to as the fermionized state). More importantly, 1/N is also the probability obtained if the spin up particles had instead been strongly interacting bosons (see Methods for further details). This demonstrates a very strong deviation from the common perception of the similarity of strongly interacting fermions and bosons in 1D.

DiscussionThe strongly interacting regime is very difcult to access both numerically and analytically because of many (nearly) degenerate energy eigenstates. The present approach nds the exact solution to linear order in 1/g in a manner that automatically yields eigenstates that are adiabatically connected to the eigenstates for smaller values of g. Combining our analytical approach with numerical techniques that perform well in the weakly and intermediate strength regime will allow us to access the quantitative and qualitative behaviour of 1D systems in arbitrary conning geometries. Furthermore, our approach provides a necessary starting point for including higher orders in 1/g (refs 30,31) and it represents an essential benchmark for numerical calculations. References 17,18,25 contain recent numerical calculations very close to the 1/g 0 limit. The results presented

here are in agreement with those calculations. A leading numerical method used for 1D systems is the density matrix renormalization group technique. As this is intrinsically a variational approach, it is not clear how well it will perform in strongly interacting limit for multicomponent systems where the degeneracy of the spectrum is large. In addition, our method can provide a fairly simple benchmark of density matrix renormalization group by using the exact slopes of the energy as 1/g-0.

It is extremely important to note that with no extra effort our method obeys the LiebMattis theorem32, which states that the energies may be ordered according to total spin, S, with the ground state having the minimal S. It is tempting to conclude that the exact solution can be obtained by constructing eigenstates with well-dened total spin at 1/g 0, and many previous

attempts to solve this problem have been based on spin algebra and spin mappings in one way or another21,22. However, this construction is not unique and our method demonstrates that the condition is insufcient to determine the eigenstates to order 1/g for N43 for multicomponent systems where M43. Insisting on states that are eigenstates of the total spin provides just one extra constraint and this is not enough to determine all the coefcients ak. A more direct way to see that the spin algebra approach is incomplete is to note that the construction of the spin functions ignores the connement V(x)22. While it does yield eigenstates in the strict limit 1/g 0, the spin states obtained for N43 are

generally not adiabatically connected to states at large but nite g. The method presented here overcomes this naturally by extremizing the slope of the energy as the criteria that determines the eigenstates.

Several previous papers have introduced BoseFermi and FermiFermi duality mappings for interacting 1D systems3336. This is a very nice mathematical tool for transforming between strongly and weakly interacting systems; however, it must be stressed that these techniques cannot be used to solve the problem considered here. Many of these techniques start from the antisymmetrized noninteracting state and then multiplies by a factor that ensures that under exchange of two particles the sign comes out correctly (plus for two bosons and minus for two fermions). As we have clearly shown, the coefcients (ak) are

generally not integer. For multicomponent systems and for Bose Fermi mixtures, duality mappings must be supplemented by knowledge of the solution on either the fermion or boson side of

0.50.450.40.350.3

100

102

P

104

106

0.250.20.150.1

0.0502 3 4 5 6

N

N

1082 10

P

7 8 9 10

Figure 2 | Tunnelling of a strongly interacting impurity in a Fermi sea. (a) Illustration of a tunnelling experiment where the trap is opened on one side and the out-going particle and its spin is detected14 for the case of Nm 3 and Nk 1. (b) Probability to nd the spin down impurity on the far

right for Nm 1,?, 9 in a b 0 harmonic trap (solid) and a hard wall

trap (dashed) for the ground state. The dotted line is 1/N and is the probability in both traps for the non-interacting state, CA. It is also the result if the Nm majority particles are strongly interacting identical bosons. The inset shows the same data on a double-log plot. The dash-dotted line is 1/N! for comparison.

4 NATURE COMMUNICATIONS | 5:5300 | DOI: 10.1038/ncomms6300 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications & 2014 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms6300 ARTICLE

the duality transformation. Without our solution one would merely be mapping into another unsolved problem. On the other hand, combining our technique and duality transformations we expand the class of solvable systems.

Methods

Details of the energy-functional derivation. As noted in the main text, in the strongly interacting limit, 1/g-0, the boundary conditions in equation (2) imply that the total wave function C must vanish whenever xi xj. For two identical

fermions, this is trivial since the Pauli principle dictates that both sides of equation (2) vanish. For non-identical particles, the wave function must still vanish when they overlap; however, the derivatives from each side can generally be different since the Pauli principle provides no restrictions. The basic idea of our method is now very simple. First construct an antisymmetric function, CA,

with energy EA using single-particle states as described in the main text. The most general N-body wave function is shown in equation (3). The number of independent coefcients in equation (3) can be deduced from Pauli symmetry to be M N!/S, where the symmetry factor is calculated according to the number of

groups of identical particles in the system; S N1!N2!yNn! if there are n groups of

identical particles with N1,N2,y,Nn in each group. Now we can construct an energy-functional, K(a1,y,aM), such that E EA K/g. Subsequently, we vary K

with respect to ak and diagonalize the resulting linear system to obtain the exact eigenstates and slopes of the energy to linear order in 1/g. Intuitively, the functional gives the slopes of the energy so that the ground state on the repulsive side (g40)

around 1/g 0 will maximize K, the rst excited state will be the next extreme

point, and so on.

The proof of equation (4) is an exercise in application of techniques from standard quantum mechanics. Using either perturbation theory or the HellmannFeynman theorem, we have

K

@E@g 1 g2

@E

@g lim

g!1

g2

Pi4j

R QN

k1 dxk j C j 2 dxi xj

hC j Ci

; 5

where the Dirac bracket /C|CS denotes the normalization integral. The dependence on g can be eliminated by using equation (2) in equation (5) that yields

K

1 4

C

@@x

x x !0

@ @x

@@x

x x !0

@ @x

Pi4j

R QN

i1 dxidxi xj

2

hC j Ci

;

6

where it is important that one rst evaluates the derivatives and then integrates out the delta function. Note that if i and j are the indices of identical fermions, then the Pauli principle requires antisymmetry under permutation and we get a vanishing contribution to K. We can now split this integral into N! sectors with different particle orderings, that is,

K

1 4

R

G

QNi1 dxidxi xj@@x @@x

x x !0

C

@@x

x x !0

Pioj

Pk

@ @x

2

hC j Ci

;

7

where we sum over permutation of the coordinates, Pk, and the integration regions, Gk, are such that xP 1

oxP 2

o:::oxP N

. This is very much in the spirit of the

Bethe ansatz of course. In the Bethe ansatz, the assumption is that the two-body potential scatters without diffraction37. In our case, the very notion of scattering is compromised by the presence of the external trap and not even asymptotically can we talk about free particles in a general conning 1D geometry. Here we use instead the local properties of the two-body interaction.

From every boundary where two particles coincide, we will obtain a factor (ak ap)2 times a derivative of CA. If they are identical fermions, then ak ap and

the term vanishes; however, to keep it general we do not make such an assumption here. K can thus be written as a sum of quadratic differences of the ak coefcients.

Likewise, we may use the normalization C C

j

h i

Pk a2k (corresponding to unit normalization on each of the M sectors in the expansion equation (3)). Therefore,

K can be written as

K

Pk;p

8 >

>

>

>

>

<

>

>

>

>

>

:

ak ap2ak;p

Pa2k; 8 where k and p run over the number of independent coefcients M N!/S and we

have used the antisymmetry of CA to eliminate the factor 1/4. This is equation (4) of the main text. The quantity ak,p is dened as

ak;p

ZG

dx1 . . . dxNdxi xj

10

where we have xed x1 and x2 to be spin projection up, while x3 and x4 have spin projection down. The topology of each conguration is indicated by the arrows. Note that we have only written the independent pieces of the wave function; the remaining terms are dictated by the Pauli principle. In the example we consider the ground state manifold (lowest energy at 1/g 0), meaning that CA(x1, x2, x3, x4) is

the antisymmetric function formed by occupying the four lowest states in the conning potential. The functional for the energy around 1/g 0 using this basis

becomes

K

a1;2

@CA @xi

a1 a22 a2;3a2 a32 a2;4a2 a42 a3;5a3 a52 a4;5a4 a52 a5;6a5 a62

P6ki a 2 k

: 9

In addition, we rst have to take the derivative before integrating over the delta function. Here Pk is a permutation of the coordinates that has the property that xi

and xj are next to each other so that they can interact, while p denotes a permutation, Pp, of the same kind but with xi and xj in reverse orders. This shows why we do not need to put the index p explicitly on the right-hand side sinceit is uniquely specied for given k, i and j. Note that these integrals will generally also depend on the ordering of all the other particles besides xi and xj, which is specied by Pk.

The decisive observation is that the ground state in the vicinity of 1/g-0 will be the state that maximizes the slope K. In fact, all sets of ak that extremize K dene a wave function that is an eigenstate around 1/g-0, and these will be orthogonal. This is proved as follows. First dene a basis of states given by setting ak 1 and

ap 0 for pak; this denes a set of M so-called bump functions that all have

energy EA through CA. We now apply degenerate perturbation theory to rst order, which yields a secular matrix (to be discussed below) whose eigenvalues are the slopes Ki and eigenstates are the correct eigenfunctions for 1/g-0. The result now follows from the linear variation method, which states that the extremizing combinations are orthogonal eigenstates. We have just shown that we can use either degenerate perturbation or variation to nd the exact wave functions for the ground state for 1/g-0 and we obtain the slopes, K, automatically. It is straightforward to argue for the adiabatic connection of the ground state at g 0

and the ground state around 1/g-0 for the lowest EA value possible where CA constructed by one particle in each of the N lowest single-particle states. This follows from the LiebMattis theorem and the fact that the largest total spin state is uniquely dened. For higher states, one must be more careful in connecting the states, and symmetry classications at both weak and strong interactions is a very useful tool38. However, we stress that symmetries (permutation group, parity invariance and so on) cannot be used to determine the ak coefcients themselves for general external connement and arbitrary number of particles.

The determination of the amplitudes, ak, now proceeds by linearizationof the functional, qK/qak 0. This produces an eigenvalue equation of the form

Au Ku, where u is a vector with ak as entries, while A is a symmetric matrix

containing combinations of the ak,p coefcients. By diagonalizing A we obtain the orthogonal and complete set of eigenstates. This completes the proof of the solvability of the strongly interacting problem in an arbitrary conning potential in 1D to linear order in 1/g. The derivation above allows us to write down the effective strong interaction Hamiltonian Heff EAIM (1/g)A, where IM is an

M-dimensional identity matrix. Heff can be used to study additional perturbations

on the system such as external electromagnetic elds in the strongly interacting regime.

The exact solutions generically have different ak coefcients (in fact a subset of coefcients can even vanish in a given eigenstate). This explains why it is very difcult to achieve convergent results in the strongly interacting limit using numerical techniques that are not optimized to take this into account. By taking different CA with different energies EA we can now build the entire spectrum that will consist of a ladder of states each with an M-fold degeneracy at 1/g 0 and

determine their slopes, K, around 1/g 0. This is illustrated for ground state and

excited state manifolds in Fig. 1 in the main text.

Fermions with I42 internal states or colours (such as an SU(I) model) are solved by exactly the same method but with a different M, depending on the number of such colours. A minor adjustment for strongly interacting bosons is that when two identical bosons are interchanged in permutations Pk and Pp of equation (4), we must take ak ap to compensate the antisymmetry of CA and

consequently add a term 4akpa2k to the numerator in equation (8) to account for the interactions in equation (1). Mixtures of fermions and bosons run along the same lines. Our only assumptions are that the particles have equal masses, the same interaction strength g between all components and a conning potential, V(x), which is the same for all particles.

Details for four-body systems. The two spin up and two spin down systems discussed in the text have the general wave function

C

a1CA for x1ox2ox3ox4 ""##

a2CA for x1ox3ox2ox4 "#"#

a3CA for x3ox1ox2ox4 #""#

a4CA for x1ox3ox4ox2 "##"

a5CA for x3ox1ox4ox2 #"#"

a6CA for x3ox4ox1ox2 ##""

2

;

11

where we notice the absence of a term with a3,4 since those congurations have no matching boundaries. By using parity invariance, one sees that a1,2 a5,6 a and

NATURE COMMUNICATIONS | 5:5300 | DOI: 10.1038/ncomms6300 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

& 2014 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms6300

a2,4 a3,5 a4,5 b, so that we have only two independent coefcients. By

variation of K with respect to ak, we obtain the A matrix that has the rather simple form

A

a a 0 0 0 0 a a 2b b b 0 0

0 b 2b 0 b 0

0 b 0 2b b 0

0 0 b b a 2b a

0 0 0 0 a a

The different coefcients of the matrix are

ak;k 1

ZG Y

N

l1

a2,3

2

6666664

3

7777775

dxldxk1 x1 j

@CA @xk1 j 2

; 18

for k 1,y, N/2 for N even and k 1,y, (N 1)/2 for N odd. Here the integral is

over the region G x2ox3oyox1oxk1oyoxN. In addition, we note that

these constants are not equal for general potentials. However, once more the hard wall connement is a truly special case. There one can prove that the ak,k1 are equal and one recovers the Bethe ansatz results28.

The results above can be applied to the case where x2,y, xN are identical bosons instead of fermions. For the Hamiltonian in equation (1) with a single coupling, g, for all pair-wise interactions, the ground state is the one found by Girardeau many years ago9 since the Hamiltonian does not distinguish between the identical bosons and the impurity. Alternatively, by counting interacting pairs in each conguration one can show that in a hard wall connement the two-component boson case is obtained by adding a diagonal matrix to A of the form diag(2(N 2), 2(N 3), 2(N 3),y, 2(N 3), 2(N 2)). Note that the inter- and

intraspecies interactions are the same for these two-component bosons, and the coupling constant goes to innity (1/g-0). This follows from the structure of our Hamiltonian that does not distinguish between the two species. An interesting example of non-identical inter- and intraspecies interactions can be found inref. 39.

For general potentials with nonconstant ak,k1, the matrix to add will still be

diagonal but now the entries are sums over a set of ak,k1 with a contribution from

each adjacent pairs of identical bosons for a given conguration. Let us denote the two species as A and B and consider the example with three A-type bosons and one B-type (acting as the impurity). In this case the four diagonal terms (corresponding to congurations BAAA, ABAA, AABA and AAAB) are 2(a2,3 a3,4) for BAAA and

also for AAAB, and 2a1,2 for ABAA and also for AABA (there is symmetry around the centre of the diagonal because of the parity symmetry of the potential). Thus, in the limit |g|-N, we have the ground state wave function C |CA|. This has

probability 1/N for the impurity to sit on the right-hand side of the system. This is a clear demonstration of the differences between fermions and bosons when there are multiple internal states. Our formalism can also be applied to a mixture with N1 single component fermions and N2 single component bosons. However, if the BoseFermi and BoseBose zero-range interactions have different strength parameters, then the situation changes drastically and the ground state depends on the ratio of these strength parameters as they diverge to innity39.

Tunnelling. The tunnelling problem considered in the main text is considered within a simple model that can be justied in several different ways. Here we will use arguments based on the many-body tunnelling theory of Bardeen40. It has recently been shown41 that the Bardeen theory can reproduce experimental results on strongly interacting two-body systems14, and we assume that this remains true for more particles as studied in the main text.

We assume that the exact eigenstate we have obtained is the initially stationary state that starts to decay as the barrier is lowered. This initial state can be written as a sum of eigenstates of the new Hamiltonian with a lowered barrier, and we arrive at a dynamical problem. While this may sound very complicated, the new Hamiltonian is still assumed to be strongly interacting and particles are still not allowed to exchange positions. The particle that was initially closest to the barrier before it is lowered will therefore provide the dominant contribution to the ux going out of the trap. The initial state therefore contains the information about the tunnelling rates. This can be made quantitative by computing very accurate matrix elements for the rates. Since our main purpose is to point out the difference of different initial states in tunnelling experiments, we postpone further tunnelling calculations for future studies.

References

1. Deshpande, V. V., Bockkrath, M., Glazman, L. I. & Yacoby, A. Electron liquids and solids in one dimension. Nature 464, 209216 (2010).

2. Bethe, H. A. Zur Theorie der Metalle. I. Eigenwerte und Eigenfunktionen der linearen Atomkette. Z. Phys. 71, 205226 (1931).

3. Lieb, E. H. & Liniger, W. W. Exact analysis of an interacting Bose gas. The general solution and the ground state. Phys. Rev. 130, 16051616 (1963).4. McGuire, J. B. Interacting Fermions in one dimension. I. Repulsive potential.J. Math. Phys. 6, 432439 (1965).5. McGuire, J. B. Interacting Fermions in one dimension. II. Attractive potential.J. Math. Phys. 7, 123132 (1966).6. Yang, C. N. Some exact results for the many-body problem in one dimension with repulsive delta-function interaction. Phys. Rev. Lett. 19, 13121315 (1967).

7. Lieb, E. H. & Wu, F. Y. Absence of Mott transition in an exact solution of the short-range, one-band model in one dimension. Phys. Rev. Lett. 20, 14451448 (1968).

8. Tonks, L. W. The complete equation of state of one, two and three-dimensional gases of hard elastic spheres. Phys. Rev. 50, 955963 (1936).

9. Girardeau, M. D. Relationship between systems of impenetrable bosons and fermions in one dimension. J. Math. Phys. 1, 516523 (1960).

; 12

and acts on the vector u [a1, a2, a3, a4, a5, a6]T, where the superscript T denotes

the transpose. Note how parity symmetry is explicitly featured in A since it is symmetric when reading it from the top left corner to the bottom right corner or vice versa. The explicit forms of a and b are

a

Z

x ox ox

dx1dx2dx3dx4dx2 x3 j

@CA

@x3 j 2 13

@x4 j 2; 14

where the difference between the two quantities is the ordering of coordinates in the integration region. Integrating out the delta function and renaming the variables in a (x4-x3) show that we have x1ox2ox3 and x1ox3ox2, respectively.

An easy way to think about the difference is that b has the pair that interacts on the side of the system (left and right sides are equivalent because of parity invariance), while for a the interacting pair is in the middle. These two are generally not the same. However, for the particular case of the hard wall connement, it turns out that a b. In this case our results are consistent with the Bethe ansatz approach28.

For the double-well potential, aab. This is perhaps a subtle point even to experts in the eld and may explain a number of previous failed attempts to obtain the exact solution. By diagonalization of A we obtain the slopes and wave functions that allows us to determine the structures and probabilites discussed in themain text.

For the case where the two pairs of identical particles are instead bosons, there are additional interaction terms when the two identical bosons are next to each other. The modication to the A matrix is very simple and consists of adding a diagonal matrix diag(4b, 0, 2a, 2a, 0, 4b) to A. The ground state is then of the simple form u [1, 1,1,1, 1,1]T, which is easily seen to correspond to C |CF|

where CF is the wave function for spinless fermions, that is, the Girardeau wave function9. Thus, we get uniform probability of 1/3 for each of the three congurations shown in Fig. 1. In addition, we notice the differences between strongly interacting bosons and fermions.

Details for polaron systems. The fermionic polaron system where a single spin down (often called an impurity for obvious reasons) interacts with a number of identical spin up fermions is handled in a similar manner to the four-body system, although the A matrix has an even simpler structure. Let the single spin down particle have coordinate x1, and the N 1 spin up particles have coordinates x2,y,

xN. There are N!/(N 1)! N independent ak coefecients. The wave function is

now given by

C

b

Z

x ox ox

dx1dx2dx3dx4dx2 x4 j

@CA

8 >

>

>

>

>

>

<

>

>

>

>

>

>

:

a1CA for x1ox2ox3o . . . oxN 1oxN #"" . . . ""

a2CA for x2ox1ox3o . . . oxN 1oxN "#" . . . ""

a3CA for x2ox3ox1o . . . oxN 1oxN ""# . . . ""

::aNCA for x2ox3o . . . oxN 1oxNox1 """ . . . "#

15

It is now straightforward to obtain the slope, which is given by

K

PN 1 k1

ak;k 1ak ak12

PN k1 a2k

: 16

As before, we can use parity invariance of the connement to conclude that a1,2 aN 1,N, a2,3 aN 2,N 1 and so forth. This means that the number of

independent ai,j is N/2 for N even and (N 1)/2 for N odd. The corresponding

A matrix is tridiagonal with entries

A

a1;2

2

6666664

a1;2 0 0 . . . 0 0 0 a1;2 a1;2 a2;3 a2;3 0 . . . 0 0 0

: : : : . . . : : :: : : : . . . : : :0 0 0 0 . . . a2;3 a2;3 a1;2 a1;2

0 0 0 0 . . . 0 a1;2 a1;2

:

17

3

7777775

6 NATURE COMMUNICATIONS | 5:5300 | DOI: 10.1038/ncomms6300 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications & 2014 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms6300 ARTICLE

10. Paredes, B. et al. Tonks-Girardeau gas of ultracold atoms in an optical lattice. Nature 429, 277281 (2004).

11. Kinoshita, T., Wenger, T. & Weiss, D. S. Observation of a one-dimensional Tonks-Girardeau gas. Science 305, 11251128 (2004).

12. Haller, E. et al. Realization of an excited, strongly correlated quantum gas phase. Science 325, 12241227 (2009).

13. Serwane, F. et al. Deterministic preparation of a tunable few-fermion system. Science 332, 336338 (2011).

14. Zrn, G. et al. Fermionization of two distinguishable fermions. Phys. Rev. Lett. 108, 075303 (2012).

15. Greif, D., Uehlinger, T., Jotzu, G., Tarruell, L. & Esslinger, T. Short-range quantum magnetism of ultracold fermions in an optical lattice. Science 340, 13071310 (2013).

16. Wenz, A. N. et al. From few to many: observing the formation of a Fermi sea one atom at a time. Science 342, 457 (2013).

17. Gharashi, S. E. & Blume, D. Correlations of the upper branch of 1D harmonically trapped two-component Fermi gases. Phys. Rev. Lett. 111, 045302 (2013).

18. Lindgren, E. J., Rotureau, J., Forssn, C., Volosniev, A. G. & Zinner, N. T. Fermionization of two-component few-fermion systems in a one-dimensional harmonic trap. New J. Phys. 16, 063003 (2014).

19. Jochim, S. Talk at ITAMP (Cambridge, MA, USA, 27 March 2013).20. Ogata, M. & Shiba, H. Bethe-ansatz wave function, momentum distribution, and spin correlation in the one-dimensional strongly correlated Hubbard model. Phys. Rev. B 41, 2326 (1990).

21. Deuretzbacher, F. et al. Exact solution of strongly interacting quasi-one-dimensional spinor bose gases. Phys. Rev. Lett. 100, 160405 (2008).

22. Guan, L., Chen, S., Wang, Y. & Ma, Z.-Q. Exact solution for innitely strongly interacting Fermi gases in tight waveguides. Phys. Rev. Lett. 102, 160402 (2009).

23. Tan, S. Energetics of a strongly correlated Fermi gas. Ann. Phys. 323, 29522970 (2008).

24. Barth, M. & Zwerger, W. Tan relations in one dimension. Ann. Phys. 326, 25442565 (2011).

25. Deuretzbacher, F. et al. Quantum magnetism without lattices in strongly-interacting one-dimensional spinor gases. Phys. Rev. A 90, 013611 (2014).

26. Volosniev, A. G. et al. Engineering the dynamics of effective spin-chain models for strongly interacting atomic gases. Preprint at http://arxiv.prg/abs/1408.3414

Web End =http://arxiv.prg/abs/1408.3414 (2014).

27. Auerbach, A. Interacting Electrons and Quantum Magnetism (Springer-Verlag, 1998).

28. Oelkers, N., Batchelor, M. T., Bortz, M. & Guan, X.-W. Bethe ansatz study of one-dimensional Bose and Fermi gases with periodic and hard wall boundary conditions. J. Phys. A Math. Gen. 39, 10731098 (2006).

29. Astrakharchik, G. E. & Brouzos, I. Trapped one-dimensional ideal Fermi gas with a single impurity. Phys. Rev. A 88, 021602(R) (2013).

30. Busch, T., Englert, B.-G., Rzazewski, K. & Wilkens, M. Two cold atoms in a harmonic trap. Found. Phys. 28, 549559 (1998).

31. Sen, D. Perturbation theory for singular potentials in quantum mechanics. Int.J. Mod. Phys. 14, 1789 (1999).32. Lieb, E. H. & Mattis, D. Theory of ferromagnetism and the ordering of electronic energy levels. Phys. Rev. 125, 164 (1962).

33. Cheon, T. & Shigehara, T. Fermion-Boson duality of one-dimensional quantum particles with generalized contact interactions. Phys. Rev. Lett. 82, 2536 (1999).

34. Girardeau, M. D. & Olshanii, M. Theory of spinor Fermi and Bose gases in tight atom waveguides. Phys. Rev. A 70, 023608 (2004).

35. Girardeau, M. D. Two super-Tonks-Girardeau states of a trapped one-dimensional spinor Fermi gas. Phys. Rev. A 82, 011607(R) (2010).

36. Girardeau, M. D. Tonks-Girardeau and super-Tonks-Girardeau states of a trapped one-dimensional spinor Bose gas. Phys. Rev. A 83, 011601(R) (2011).

37. Sutherland, B. Beautiful Models (World Scientic, 2004).38. Harshman, N. Spectroscopy for a few atoms harmonically trapped in one dimension. Phys. Rev. A 89, 033633 (2014).

39. Zinner, N. T. et al. Fractional energy states of strongly-interacting bosons in one dimension. Europhys. Lett. 107, 60003 (2014).

40. Bardeen, J. Tunnelling from a many-particle point of view. Phys. Rev. Lett. 6, 57 (1961).

41. Rontani, M. Tunneling theory of two interacting atoms in a trap. Phys. Rev. Lett. 108, 115302 (2012).

Acknowledgements

We thank J. Arlt, D. Blume, G. M. Bruun, F. Mila, K. Riisager, A. Svane andS. Gammelmark for discussion and feedback on the manuscript. This work was funded by the Danish Council for Independent Research DFF Natural Sciences, the DFF Sapere Aude program and EPSRC grant No. EP/J001392/1.

Author contributions

A.G.V. and N.T.Z. devised the project and developed the formalism in collaboration with M.V., A.S.J. and D.V.F. The numerical calculations were carried out by A.G.V. under the supervision of D.V.F. and by N.T.Z. The initial draft of the manuscript was written by A.G.V., N.T.Z. and M.V. and all authors contributed to the revisions that led to the nal version.

Additional information

Competing nancial interests: The authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/

Web End =http://npg.nature.com/ http://npg.nature.com/reprintsandpermissions/

Web End =reprintsandpermissions/

How to cite this article: Volosniev A. G. et al. Strongly interacting conned quantum systems in one dimension. Nat. Commun. 5:5300 doi: 10.1038/ncomms6300 (2014).

NATURE COMMUNICATIONS | 5:5300 | DOI: 10.1038/ncomms6300 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 7

& 2014 Macmillan Publishers Limited. All rights reserved.