Quantum magnetism in strongly interacting one-dimensional spinor Bose systems

Strongly interacting one-dimensional quantum systems often behave in a manner that is distinctly different from their higher-dimensional counterparts. When a particle attempts to move in a one-dimensional environment it will unavoidably have to interact and ‘push’ other particles in order to execute a pattern of motion, irrespective of whether the particles are fermions or bosons. A present frontier in both theory and experiment are mixed systems of different species and/or particles with multiple internal degrees of freedom. Here we consider trapped two-component bosons with short-range inter-species interactions much larger than their intra-species interactions and show that they have novel energetic and magnetic properties. In the strongly interacting regime, these systems have energies that are fractions of the basic harmonic oscillator trap quantum and have spatially separated ground states with manifestly ferromagnetic wave functions. Furthermore, we predict excited states that have perfect antiferromagnetic ordering. This holds for both balanced and imbalanced systems, and we show that it is a generic feature as one crosses from few- to many-body systems.


Results
Energetics and wave functions. Our  m , and energy, ω  , which we will use throughout (here  is Planck's constant divided by π 2 ). We assume short-range interactions between A and B particles that we model by a Dirac delta-function parameterized by an interaction strength, g, i.e.
where x and y denote the coordinates of A and B particles, respectively. The intraspecies interaction strengths are assumed to be much smaller than g and we will therefore neglect such terms. To access the quantum mechanical properties of our system we must solve the N -body Schrödinger equation. This will be done using novel analytical tools and using exact diagonalization. In the latter case we have adapted an effective interaction approach that has recently been succesfully applied to fermions in harmonic traps 42,41 (see the Methods section for further details). The analytical and numerical methods allow us to address up to ten particles, which is larger than most previous studies not based on stochastic or Monte Carlo techniques. The simplest non-trivial case is the three-body system which has two A and one B particle. The energy spectrum is shown in Fig. 1 as a function of g. The most interesting feature to notice is the ground state behaviour as / → + g 1 0 . Here, an odd and an even state become degenerate at an energy of ω .  2 5 . This should be contrasted to the behaviour of single-component bosons or two-component fermions which will always have energies that are an integer times ω  when / → g 1 0. Furthermore, we notice how the two states that merge at / = g 1 0 become two excited state branches on the attractive side of the resonance but the even parity state remains the lower one. This is opposite to the behaviour of fermions 43 where the hierarchy of states is inverted at / = g 1 0. The ground state for large and negative g is very different as it contains deeply bound molecules, which we will not consider further. The fractional energies in the spectrum can be explained by a schematic three-body model and in stochastic variational calculations 29 . This provides a hint that larger systems could also display fractional energy states in the strongly interacting limit and begs the question as to what spatial configurations such states correspond to.
We will now show that the fractional energy states are generic for strongly interacting two-component bosons in 1D and, importantly, for the ground state they realize perfect ferromagnetic behaviour irrespective of whether the system is balanced ( = N N A B ) or not. The term perfect ferromagnetic behaviour implies that we have a full spatial separation of the two components in the exact ground state wave function of the system, i.e. the probability to find only A on one side and only B on the other side of the system is not just dominant, it is exactly unity. The ground state has only a single 'domain wall' at which an A and a B particle are neighbours. As a consequence, imagine that you detect an A particle on the left (right) side of the system, then you can immediately conclude that all the B particles reside to the right (left) of this A particle. Balanced systems. We first consider a four-body system that has two A and two B particles. The energy spectrum for > g 0 is shown in Fig. 2a). A striking feature is the two-fold degenerate ground state for / → g 1 0 that has a non-integer energy similar to the three-body problem. In this strongly repulsive limit, the system realizes a perfect spatially ferromagnetic quantum state as we will now demonstrate analytically.
First we note that the center-of-mass motion of the four-body system can be separated and thus ignored. This leaves three Jacobi coordinates to describe the system. The details of these reductions can be found in the Methods section below. In Fig. 2b) we show the space of the Jacobi coordinates and highlight all the planes at which an AB (solid planes), an AA or a BB (checkerboard planes) pair of particles overlap. The main observation is that as / → g 1 0, the wave function must vanish on all the solid planes in Fig. 2b) and we arrive at the disconnected regions shown with different colours. The particles become effectively impenetrable and we may characterize the wave function by specifying the amplitudes of all possible spatial configurations of the four particles. These regions correspond to specific orderings on a line of the four particles. In particular, the large (red) region dominating the figure corresponds to spatial configurations AABB or BBAA. The green region occupies half the spatial volume of the red and corresponds to ABBA or BAAB, while the yellow region has one-fourth the volume of the red region and corresponds to ABAB or BABA configurations. A wave function that vanishes on all AB interfaces may now be constructed in each of these regions. However, it is immediately clear that it will have lower energy when it can spread over a larger volume. We thus conclude that the doubly degenerate ground state at / → g 1 0 has the structure ± AABB BBAA (taking into account the parity invariance of the Hamiltonian).
As discussed in the Methods section, one may solve a simple wave equation in the red region and obtain the ground state energy to arbitrary precision. The triangles at the / = g 1 0 line in Fig. 2a) show the energies obtained in this manner. We reproduce both the ground state and a set of excited states. All of these have fractional energies and all of them are perfectly ferromagnetically ordered. The remaning states of the spectrum can be obtained by solving in the other regions of Fig. 2b). Note that states with amplitude exclusively in the yellow regions are perfectly spatial antiferromagnetic, ± ABAB BABA, and have energies (n + 1/2) ω  with integer n. They are the only parts of the spectrum which can be constructed by starting from identical fermions using Girardeau's mapping techniques 44 . The arguments presented here are neither restricted to = N 4 nor to a harmonic trapping potential and hold for any N Figure 1. Three-body spectral flow. The energy spectrum of two A and one B particle as a function of interaction strength, g, obtained by numerical calculations. In the limit / → g 1 0, the ground state becomes doubly degenerate and has half-integer energy. The contribution from center-of-mass motion has been removed. For visibility, we have dimmed states from the attractive side that diverge to large negative energies close to / → g 1 0.
Scientific are adjacent and connected. This implies that the volume for ABBA is half as large.
In order to further study the correlation between the A and B particle subsystems, one can use the pair-correlation function. The pair-correlation function measures the probability of finding a particle from one subsystem (say of type A) at position x 2 given that we know the position, x 1 , of a particle in the other subsystem (say of type B). The pair-correlation function for the balanced case is shown in Fig. 3a) for the case where g is zero and therefore no separation and hence no specific ordering is present, but as g becomes strongly repulsive evidence of separation is seen as illustrated in Fig. 3b,c). In particular, Fig. 3c) shows that if we have a given particle from one subsystem situated on the negative x-axis ( < x 0 1 ) then a particle from the other subsystem is most likely to be found on the positive x-axis ( > x 0 2 ) and vice versa. These numerical calculations are also supported by our analytical method for the four-body system at / = g 1 0. Fig. 3d) shows the analytical result that is virtually identical to the numerical results of Fig. 3c.
Larger systems may in principle be handled in similar fashion by solving wave equations with proper boundary conditions and obtaining the fractional energies in the limit / → g 1 0. However, the increase in dimension of the problem makes this very difficult in practice. In order to further demonstrate that balanced systems have perfect ferromagnetic ground states irrespective of particle number, we have numerically computed the ground state densities for systems with as shown in Fig. 4. Evidence of the separation of A and B can be seen in the total density already in Fig. 4a as g increases (note the perfect agreement with the analytical result in the limit / → g 1 0). We expect the two degenerate ground states to have structure … … ± … … A AB B B BA A. In order to prove this perfect ferromagnetic behaviour, we consider the odd and even superposition of the two degenerate states which we expect will yield states with exclusively A or B particles on either side of the system As the ground state for / → g 1 0 is spatially separated, one may speculate that it can be understood physically as two ideal Bose gases or 'condensates' sitting on either side of the system even in this strongly interacting limit. In Fig. 4d) we plot the densities in a rescaled fashion where we multiply by / case indicates that the system does behave as two ideal Bose gases as the particle number grows. In the limit , we would expect the overlap of the two gases to vanish as the energy cost of overlap goes to infinity. We therefore expect that the occupied mode in this large system limit is the first excited state of the harmonic trap which vanishes at the center. The dashed line in Fig. 4d) shows this state rescaled to . This analytical guess displays the same features as the numerical densities and we conclude that already for ten particles the many-body properties are emerging. Imbalanced systems. The extremely imbalanced limit, where and N B varies, provides a realization of a strongly interacting Bose polaron in 1D, i.e. an impurity that interacts strongly with an ideal Bose gas. In Fig. 5 we plot the densities of systems with . We see that the impurity sits at the edge of the system (Fig. 5a), while the majority component tends to occupy the center (Fig. 5d). We confirm the numerical results by employing an analytical model, which shows excellent agreement. The details can be found in the Methods section. To confirm that the wave function of the strongly interacting ground state has intrinsic phase separation, i.e. has the form … ± … AB B B BA, we plot the densities for a sum of the nearly degenerate ground states in Fig. 5b,c). As in the balanced case above, we find a perfectly separated ground state behaviour. For instance, if we locate the single A particle on one side of the trap, we would thus immediately know that all the B particles reside on the other side, and vice versa. This behaviour is opposite to the case where the B particles are identical fermions where the impurity resides mainly in the center of the system 41 . We have confirmed that this structure is , and it is therefore a generic feature that the two species are perfectly spatially separated (ferromagnetic) in the ground state for strong interactions.
A remarkable feature of the densities in Fig. 5b,c) is the movement of the centroids of the peaks with particle number. We clearly see the majority moving into the center and the impurity being pushed toward the edge. This demonstrates how an ideal condensate is being built in the center. For large N B the energy per particle goes to ω /  1 2 , which implies a single-mode condensate that is becoming macroscopically occupied (see Methods section for details). The relative deviation between numerical and analytical energies is below three percent for . We thus have an analytic model for the crossover between the few-and many-body limit for the bosonic polaron in one dimension. This includes the external trap that is a reality of most experiments.

Discussion
We have shown that a mixture of two ideal Bose systems in one dimension has unusual properties when the inter-species interactions is strong. The systems have energies that are non-integer multiples of ω  . In Fig. 6 we show the ground state energies for / = . g 1 0 01 in systems with ten particles or less relative to a ground state with only B particles. Driving an A to B transition via radiofrequency spectroscopy would be a possible way to confirm the predicted energies in Fig. 6. This technique has been demonstrated for fermions in recent few-body experiments in 1D 17 . The results in Fig. 6 show that the energy per particle tends to saturate for large systems and that this happens faster the more imbalanced the system is. We also see that the balanced case has an almost linear energy dependence.
The ferro-and antiferromagnetic states can also be detected by measuring momentum distributions. In Fig. 7a,d) we show the distributions for / = .
, respectively. The purple solid line is the even parity ground state while the solid green line is the excited state with antiferromagnetic ordering. The striking difference of the two distributions implies that they should be easily identifiable in experiments. For comparison, the solid black line with a multi-peak structure shows the distribution for a system of identical fermions. Comparing the solid green and black curves in Fig. 7a,d) clearly demonstrates that, in spite of the fact that these two states have equal energy, the correlations are very different. In Fig. 7a) the dashed black line corresponds to the Tonks-Girardeau hard-core boson state 44 , which is also seen to be very different from the states discussed here. For the density has been rescaled to a total density with . The dashed line corresponds to the density expected in the many-body limit imbalanced systems, we find that measuring the impurity momentum distribution, Fig. 7b), yields information about the parity of the state. On the other hand, the majority distributions in Fig. 7c) are identical in the two opposite parity ground states. A characteristic feature seen in Fig. 7b) is the development of oscillatory structure as the number of majority particles increases and pushes the impurity further out  in the trap, see also Fig. 5b). A third technique for experimentally addressing the systems we study is controlled tunneling as the trap is gradually lowered 16 (see Methods below).
The separation of components in the ground state for strong interactions is intrinsic to both balanced and imbalanced mixtures, as is the presence of other spatial configurations in specific excited states. Furthermore, the magnetic behaviour discussed above is not connected to the harmonic confinement and should be seen in an arbitrary confining geometry. While we have studied the balanced and the extremely imbalanced limits here, we have checked numerically that the spatial separation of components is an intrinsic feature of the system for systems with ten or less particles. We therefore infer that this will hold also for larger systems, regardless of the population ratio. A simple physical picture can be given in terms of domain walls, i.e. points at which the two components interface. The system tends to minimize the number of domain walls and this principle can be used to understand the ferromagnetic ground state and predict the ordering in energy of other configurations.
In the paradigmatic two-component (spin 1/2) Fermi system, the ground state is never purely ferro-or antiferromagnetic for strong interactions 43 , and Bose mixtures therefore provide a unique set of quantum ground states for exploring and exploiting magnetic behaviour. The description of these systems clearly goes beyond the famous Bose-Fermi mappings 44,45 ] and we provide not only numerical but also new analytical tools to fill this gap. Importantly, we demonstrate that the crossover from few-to many-body physics can be studied already at the level of ten particles.

Methods
Numerical method. We solve numerically the many-body Schrödinger equation by exact diagonalization with the full Hamiltonian projected onto a finite basis constructed from harmonic oscillator single-particle states. Each many-body basis state is written as a product of symmetrized states of N A and N B particles. The model space truncation is defined by an upper limit of the total energy.
Instead of the bare zero-range interaction in (1), we consider an effective two-body interaction in order to speed up the convergence of the eigenstates with respect to the size of the many-body basis. The effective potential is constructed in a truncated two-body space, and is designed such that its solutions correspond to exact two-body solutions given by the Busch formula 46 . As explained in detail in Refs. 41,42, this is achieved using a unitary transformation that involves the lowest eigensolutions given by the Busch formula. By construction, this unitary transformation approach will reproduce exact bare Hamiltonian results for the many-body system (both energy spectrum and wave functions) in the limit of infinite model space.
The excellent convergence property of this effective-interaction approach was demonstrated in Ref. 41 and is key to the quality of our numerical results and to our conclusions. In the construction of the effective interactions we benefit from having access to the exact two-body solutions for short-range interactions in harmonic traps. However, we stress that using numerical two-body solutions this approach can be generalized to study many-body systems in higher dimensions with finite-range interactions and in any trapping potential.
Density and pair-correlation profiles. In the second quantization formalism the density profile is cal- . In our analytical results we use the following expressions to obtain the density and pair-correlation from the N -body wave function, ψ ( , …, , , …, ) . To get the density of A particles in the system we need to calculate Likewise, to obtain the density of a B particle we integrate over all variables except y 1 instead of x 1 . The AB pair-correlation is obtained by calculating Analytics for balanced systems. Here we outline the calculational procedures required to obtain the exact solutions for the four-body system in the / → g 1 0 limit. The method can in principle be extended to larger systems, but it becomes increasingly difficult. In the next subsection we provide an alternative method that works well for larger systems in the imbalanced case. with g the AB interacting coupling constant and we assume that the AA and BB interactions vanish. We now perform an orthogonal coordinate transformation is the generalized Laguerre polynomial. When = g 0, the angular functions, θ φ ( , ) Y lm , are the usual spherical harmonic functions with l the total and m the projection angular momentum quantum number. However, in the limit / → g 1 0 we have to enforce non-trivial boundary conditions whenever A and B particles overlap. Let us focus on the region > z 0 by restricting to θ π ≤ ≤ / 0 2 (solutions for < z 0 may be obtained by symmetry arguments or by considering instead π θ π / ≤ ≤ 2 ). The arguments of the Dirac delta-functions in Eq. (6) vanish when tan sin 0 and θ φ ± =  1 tan cos 0. The regions defined by these conditions are illustrated in Fig. 2b). The solid red planes show exactly where the arguments of the interaction Dirac delta-functions have to vanish.
We now make the simple transformation tan cos . In these new , u v variables, the boundaries are simply ± = u 1 0 and ± = v 1 0, i.e. the function must vanish on the boundary of a square. Finally, one must transform the angular part of the Laplacian into the new variables which yields By the procedure outlined above we have transformed the problem of solving a harmonic oscillator problem in a non-trivial geometry, i.e. the red region in Fig. 2b), into solving a very simple boundary value problem We write the eigenvalue in this way so it matches the usual 3D angular eigenvalue ( + ) l l 1 . The equation for ( , ) f u v may be straighforwardly solved by using a two-dimensional Fourier expansion of the wave function. This will produce some spurious solutions as we must also impose bosonic symmetry among the two A and two B particles separately. This translates to the requirement that the solution be symmetric when reflected across the two diagonals. Notice that for each eigenvalue of this problem, λ, we obtain a whole class of solutions with energies ω λ / = + + /  E n 2 3 2 as we may add radial excitations. The low-lying solutions of Eq. (9) are given in Table. 1. State number 1, 4, 6, 11, 13, and 15 have the required bosonic symmetries for the balanced system. A number of doubly degenerate states in the spectrum may be used to construct eigenfunctions for a four-body system with = N 3 A and = N 1 B (or vice versa). In this case the wave function must vanish on the diagonal of the ( , ) u v square domain which is achievable by taking proper linear combinations. The states marked 'fermions' in Table. 1 are antisymmetric across the two diagonals in the ( , ) u v square and provide allowed states for all four-body two-component Fermi systems, i.e. 2+ 2, 3+ 1 or four identical fermions. The eigenenergies of the fermionic states have the exact values 7.5, 10.5 and 11.5. Our results differ by ⋅ − 4 10 5 which attests to the accuracy of our method. All states in Table. 1 have been obtained using a modest 400 Fourier basis states. Notice that even though the lowest perfectly antiferromagnetic state for bosons is at the same energy as the fermionic state number 5 in Table. 1, they are not related since states with the configurations ABAB or BABA solve a different boundary value problem (corresponding to the yellow regions in Fig. 2b).
The energies obtained using this (semi)-analytical approach for the = = N N 2

A B
system are given in Fig. 2a) as triangles at / = g 1 0. The two lowest triangles correspond to the angular ground state (lowest λ value) with = n 0 and = n 1. The two upper triangles are the first and second excited angular solutions both with = n 0. All four solutions have the spatial structure ± AABB BBAA. The blue dots in Fig. 4a) show the ground state density obtained by the transformation method discussed here. The rest of the spectrum at / → g 1 0 can be obtained by solving the boundary value problem in the green ( ± ABBA BAAB) and yellow areas ( ± ABAB BABA). In the latter case a fermionized (totally antisymmetric) wave function is a solution. Our main interest here is to understand the ground state so we leave the remaining states and regions for future investigations.
Analytics for imbalanced systems. The analytics provided here is applied for the Bose polaron, where ω is a 1D harmonic oscillator. The x 1 coordinate denotes the single A particle, the 'impurity' , while y i denotes the coordinates of the majority B particles. We introduce an adiabatic decomposition of the total wave function of the form which depends parametrically on x 1 . This expansion can be related to the Born-Oppenheimer approximation in which case one may consider x 1 the 'slow' variable (typically the nuclear coordinate in molecular physics). In the limit of interest / → g 1 0, we impose the condition that the total wave function vanishes for where.
The subscript y on the brackets denote integration over all , …, y y N Note that the energy E 0 provides a variational upper bound to the exact energy. The energies computed via this method for the polaron are shown in Fig. 6 and agree with the numerical results to within a few percent for the largest particle numbers in the figure. We expect the agreement to become better for even larger particle numbers. Furthermore, since we obtain the full wave function in an analytical form we may also compute the densities of both impurity and majority components. In Fig. 5a) we show the impurity densities for = N 3 while Fig. 5d) shows the corresponding majority density. We see a striking agreement between the numerical results and the analytically tractable model presented here. The model presented here can be extended to excited states and also to systems Tunneling experiments. For 1D few-body systems it is now possible to experimentally access the ratio of tunnelling probabilities for particles with different spins 16 . For strongly interacting fermionic systems this ratio can be obtained in a rather simple way without knowing the parameters of the experiment (e.g. height of the barrier), as it is determined solely by the probability for the particle to be at the edge of the trap 43 . For strongly interacting bosonic systems this is unfortunately no longer the case and the detailed parameters of the experiment are needed. To illustrate this consider a system with one impurity in a sea of N majority bosons. For  N 1 we have shown that for the ground state the impurity is pushed to the edge of the trap whereas the N bosons can be accurately described using only the lowest energy level of the harmonic trap. Let us assume that the trap is lowered on one side. If the impurity sits on this side we will detect an impurity after some time τ which may be very short, since the impurity is very close to the barrier. If, on the contrary, the impurity sits on the other side of the trap we will detect a majority particle after time τ 1 . This time may be exponentially enhanced since tunnelling is from State λ + .  Table 1. Low-lying spectrum of Eq. (9). The last column denotes the specific system for which the particular solution is an allowed eigenstate using the notation + N N A B for bosons. See the text for details.
the ground state of the trap and the barrier is consequently large. The only system where the ratio of probabilities does not depend on the geometry of the experiment is the = N N A B system where symmetry dictates that one will detect A particles as frequently as B particles. As the perfectly antiferromagnetic state is only present for , its tunneling signature is similar to many other states in the spectrum and detection of the state via tunneling would be quite difficult. For the excited states the ratio of probabilities is again very dependent on the parameters of the experiment. For example for the impurity and N bosons system if the experiment is constructed such that only one particle can tunnel then there generally will be many states where only majority particles will be detected. This happens for the states where the impurity sits closer to the center of the trap.