Emergence of junction dynamics in a strongly interacting Bose mixture

We study the dynamics of a one-dimensional system composed of a bosonic background and one impurity in single- and double-well trapping geometries. In the limit of strong interactions, this system can be modeled by a spin chain where the exchange coefficients are determined by the geometry of the trap. We observe non-trivial dynamics when the repulsion between the impurity and the background is dominant. In this regime, the system exhibits oscillations that resemble the dynamics of a Josephson junction. Furthermore, the double-well geometry allows for an enhancement in the tunneling as compared to the single-well case.


Introduction
The experimental investigation of ultracold atomic systems has made possible the realization of several celebrated models in quantum mechanics and condensed matter. These experiments are characterized by the rigorous control over the parameters of the system and by precise measurement techniques. The manipulation of optical traps, for instance, allows for the construction of different confining geometries [1], from one-dimensional tubes [2,3] to lattice systems [4,5]. Traps consisting of two wells, in particular, have been extensively employed in recent experiments [6,7,8,9,10]. They are of special interest in the study the Josephson effect [11] in cold atoms, where Bose-Einstein condensates placed in such potentials [12,13] are considered in analogy to superconductors [14,15]. The high degree of precision demonstrated in these experiments also extends to the number of particles under consideration and to the strength of the interactions between them. Two-component fermionic systems with only a few atoms [16,17,18] and paradigmatic models such as the infinitely repulsive Tonks-Girardeau gas [19,20] can also be realized and studied in the lab. Combinig these features, other experiments have recently explored multicomponent strongly correlated gases [21], which have shown several exotic properties [22].
The limit of strong interactions has been a particularly favored starting point in the theoretical study of one-dimensional systems with internal degrees of freedom, due to the possibility of mapping the Hamiltonian to an effective spin chain [23,24,25,26,27]. One of the key features of this mapping is that the exchange coefficients of the spin chain are solely determined by the trapping potential, and powerful numerical methods to calculate these coefficientes are now available [28,29]. Approaching the problem of few atoms in a trap in the strongly interacting regime has also provided knowledge of the fundamental properties of quantum magnetism [30,31,32,33,34,35]. On the other hand, the many-body case in the limit of total population imbalance − the "impurity" problem − also presents many interesting features, such as quantum flutter [36] and Bloch oscillations in the absence of a lattice [37,38,39,40,41] (the latter having been recently observed in experiments [42]).
In this work we study a strongly interacting system composed of an impurity and a bosonic background in single-and double-well potentials. Different methods have been employed to study the many-body problem of bosons in a double-well, even outside the mean-field regime [43,44,45,46,47]. Addressing the subject from a few-body perspective [48,49,50,51,52], however, might lead to new insight on the properties of these systems. Here we show that, in the regime where the repulsion between the impurity and the background is dominant, the system can exhibit non-trivial dynamical effects: the impurity undergoes Josephson-like oscillations when initialized at the edge of the system, and can have its tunneling enhanced when a barrier is present. These effects provide new perspectives in the study of spin state transfer and quantum transport in one-dimensional systems, and should be observable using current experimental techniques.

System description and Hamiltonian
We consider the problem of an impurity confined in the presence of a background of strongly interacting bosons. We assume the impurity is a boson in an internal state defined by |↓ , while the remaining identical bosons are described by |↑ . Consequently, all atoms have the same mass m. Two-component Bose gases can be realized experimentally using, for instance, 87 Rb atoms in different hyperfine states such as |F = 2, m F = −1 and |F = 1, m F = 1 [53,54]. The number of identical bosons is given by N ↑ , while the total system size is N = N ↑ + 1. The Hamiltonian for this problem is written as where the first sum involves the single-particle Hamiltonian H 0 (see below) which is the same for both components, while the remaining terms account for the contact interactions. The coordinates are denoted by x ↓ for the impurity and x ↑i for the remaining bosons. The interaction strength is defined by g for the impuritybackground interactions, and as κg for the background-background interactions. Those parameters can be experimentally manipulated using Feshbach [55] or confinement induced resonances [56]. The single-particle Hamiltonian in Eq. 1 is given by Here, V (x) is a double-well potential (see Fig. 1) expressed as where ω is the trapping frequency. The parameterb denotes the displacement of the two minima of the wells with respect to the origin, and also defines the size of the barrier at this point as V (0) = 1 2 mω 2b2 . We then refer tob as the "barrier parameter". By makingb = 0 we naturally recover the harmonic single-well potential. Although this form of potential has analytical solutions in terms of parabolic cylinder functions [57], we obtain the single-particle wave functions and energies through numerical diagonalization (see Appendix A.1). We will focus on the behavior of the spatial distributions and the impurity dynamics in the repulsive case (g, κ > 0), for different choices of the intraspecies interaction parameter κ and the barrier parameter b. While cases of attractive interactions (g, κ < 0) can in principle be explored, the properties of the system in this regime reproduce only highly excited states related to the so-called Super Tonks-Girardeau gas [58,59]. Simulating the dynamics of systems with attractive interactions would likely require taking into account the formation of bound pairs, an effect which is beyond the scope of the formalism employed here. Throughout this work, we will consider all quantities in harmonic oscillator units; therefore, length, energy and time are given in units of l = /mω, ω and ω −1 , respectively. While the intraspecies interaction parameter κ is dimensionless, the parameter g is considered in units of 2 /ml. For simplicity, we also make the barrier parameter dimensionless by rescaling it as b =b/l. In our calculations, we set g = 20 and assume that = ω = m = 1. Figure 1: (color online) Sketch of the 1D bosonic system with an impurity in the doublewell potential. The parameterb sets the position of the minimum of each well and also the size of the barrier between them at x = 0. In the limit of strong interactions, the system can be mapped to a spin chain where the coefficients α are determined by the shape of the trap.
In the limit of strong interactions (g 1), Hamiltonian 1 can be written, up to linear order in 1/g, as the XXZ spin chain (see Appendix A.2 for details) where E 0 is the energy of the system in the limit of infinite repulsion and σ i denotes a Pauli matrix acting on site i. The coefficients α i , often called geometric coefficients, are calculated using the wave function for a system of N spinless fermions, which is constructed as the Slater determinant of the lowest occupied orbitals in the trap (see Appendix A.3). The spatial part of the wave function for a bosonic system is then obtained by means of the Fermi-Bose mapping [60]. A comparative study of the spatial distributions for a strongly interacting few-body bosonic system in the double-well has been presented in [61].

Spin densities
To obtain the probability densities for each component in the system, we must combine the spatial distributions of the atoms in the trap with the probability of magnetization of the corresponding site for an eigenstate of the spin chain described by Eq. 2. The spin densities are therefore given by , where ρ i (x) describes the individual atomic densities (see Appendix A.4), while m i ↑(↓) denotes the probability of each site in a spin wave function |χ having spin up or down. The quantities ρ ↓ (x) and ρ ↑ (x) thus describe the spatial distributions of the impurity and the background bosons, respectively. In Fig. 2 a) we show the results for the spin densities of a 3+1 system for b = 0 (single-well) and b = 2 (double-well) at different values of the intraspecies parameter κ. The cases where b = 0 correspond to the results expected for the harmonic trap, which have been broadly covered for bosonic and fermionic systems in previous works [62,63,32,64,25]. For κ < 1 the repulsion between the impurity and the background is larger than the background repulsion. This causes the impurity to be pushed to the edges of the system, which is an effect also found in the case of a weakly interacting background [65]. In this regime, the system exhibits Ising-type ferromagnetic correlations, and is characterized by a nearly degenerate ground state [35]. At κ = 1 all interactions are equal and the spin densities show the Heisenberg-type ferromagnetic profiles expected for isospin bosons [31]. In this case both distributions display the same characteristic Tonks-Girardeau spatial densities, but scaled to the number of particles in each species. When κ > 1, the repulsion between the background bosons dominates, and we observe predominantly antiferromagnetic correlations, where the impurity is placed near the center of the trap. In the double-well potential, all densities are depleted in the center of the system, but this is not the only relevant effect. For the case of κ < 1 the impurity has now a larger probability of being near the center of the trap (as compared to the single-well case), since the background density is strongly reduced in this region. It is even possible to expect a configuration (specially for a larger number of background bosons) where the impurity is completely localized near the barrier between the wells (see Fig.2 b)). This effect is directly related to the imbalance in the numerical values of the geometrical coefficients at the edges and near the center of the system as b is increased. For the cases of κ = 1 , again the impurity and the background densities have the same shape, aside from normalization. At κ > 1, we observe a similar configuration, with a small bias of the impurity toward the center of the trap.

Dynamics
We now turn to the dynamics of the impurity after being initialized at the left edge of the system. The corresponding spin state is therefore given by |↓↑↑↑ . Since this is not an eigenstate of the spin chain, we can expect the spin state to evolve in time governed by Eq. 2, and we denote it by |χ(t) . A thorough study of spin state transfer in traps of different shapes has been done by Volosniev et. al. in [24]. It has been shown [66,67,68] that transfer is optimized by considering κ = 2 (which turns Eq. 2 into an XX Hamiltonian) with a set of exchange coefficients where α j ∝ j(N − j). Here we focus on the tunneling times for the impurity between the wells (or between the left and right sides of the system in the case of a single well) when the background repulsion is smaller than the repulsion between the impurity and the background (κ < 1). We point out that, since we do not consider any other external perturbations, like trap or interaction quenches, we can assume that the spatial distributions remain in the ground state. This also allows us to consider only the manifold of Eq. 2 with lowest energy. To quantify the dynamics of the impurity, we calculate its average position as where ρ ↓ (x, t) is the time dependent spin density calculated with the spin state |χ(t) .
When considering the regime of κ < 1, we observe that the projections of the initial wave function on the two-lowest eigenstates are dominant when compared to the case of higher excited states. This allows us to attempt a two-level description for the time evolution of the spin wave function; we thus write |ψ(t) = c g e −iω R t |g + c e |e , where |g and |e denote the two egeinstates of the spin chain with lower energy, and c g and c e are the projections of the initial wave function over these states. The frequency ω J = E e − E g is given by the gap between the energies of the first excited state E e and the ground state E g . In Fig. 3, we present the results for x ↓ (t) in the single-well (b = 0) and doublewell potentials (b = 2), with two choices of κ < 1, also showing a comparison with the two-level description in each of these cases. At b = 0, we notice that the motion of the impurity between the edges of the system is very well captured by the two-level approximation. This behavior clearly resembles the oscillations in population expected for a bosonic Josephson junction described as a many-body system in a double-well. In the present context, however, the barrier is composed by the repulsive background gas. We expect these results to hold even in the case of more than one impurity, provided that the system is imbalanced (that is, the background gas must have a larger number of particles). In this situation, an initial state decribed by the minority species completely localized at either side of the trap should have its time evolution governed mainly by the two lowest energy states. In the single-well case with weaker intraspecies interaction (κ = 0.2) the tunneling of the impurity is suppressed. Here, the the behavior of the background approaches that of an ideal Bose gas, where the atoms tend to "bunch up" in the center of the trap. Now, comparing the single-and double-well cases, we see that, for κ = 0.5, the presence of the barrier slows down the tunneling of the impurity. Furthermore, we observe oscillations on a smaller scale, due to a larger overlap between the initial state and the excited states of the spin chain Hamiltonian. At κ = 0.2, however, we get an enhanced tunneling of the impurity when considering a double-well as opposed to the single-well case. This effect has been also found with a different choice of double-well potential [24]. One might interpret it as a splitting of the background gas by the barrier in such a way that the impurity is able to tunnel through faster than it would in the absence of the barrier. However, if we consider a single-particle problem where an atom is initialized in the left well, it is clear that increasing the barrier size would only lead to exponential suppression in the tunneling frequencies. We therefore conclude that the accelerated tunneling observed in the regimes we consider is only possible due to the presence of the bosonic background, and thus constitutes a many-body effect.
To get an understanding of this behavior over a larger parameter space, we plot in Fig. 4 the energy gap between the ground state and the first excited state for several values of κ and b. The non-monotonic behavior of the gap as a function of b indicates that, for small κ, there is some choice of barrier size that increases that energy gap, and therefore enables a higher tunneling frequency between the wells. As κ increases, however, we see that this behavior disappears and the presence of the barrier only reduces the gap, thus making the tunneling slower.

Increasing N ↑
As a final example, we consider a case where we increase the number of background bosons to N ↑ = 5 (to observe similar effects as in the case of N ↑ = 3, we choose to maintain an even total number of atoms). The initial state is once again defined by the impurity placed at the left edge of the system, that is, |↓↑↑↑↑↑ . We keep the intraspecies repulsion parameter fixed at κ = 0.25. In Fig. 5 a) we once again show the results for the average position of the impurity as a function of time. For the single-well, the tunneling times are so long that the impurity is effectively frozen at the left edge of the system. In this case, an analogy can be drawn to the self-trapped regime in a few-body system as presented in [48]. For b = 3, however, we again notice that a faster motion of the impurity from the left to the right well is induced. The difference in the results with and without the barrier is even clearer than in the case of N ↑ = 3. This can also be seen in the energy gap between the two lowest states, as presented in Fig. 5  b): a very pronounced curve shows the increase in this quantity for small κ and b > 1. We point out that the final time (t = 10 4 in units of ω −1 ) considered in the present case is five times larger than in the case of N ↑ = 3. Time scales for harmonically trapped systems are set by the inverse frequency, which in present experiments with few-body cold atoms is of around 100µs [18]. The total times obtained in experimental setups can therefore be decreased by considering tighter traps.

Conclusions
We have studied the static and dynamic properties of an impurity in the presence of a background of bosons in single-well and double-well geometries. The ground state spin densities are described by a combination of the spatial distributions in the limit of infinite repulsion and the eigenstates of a spin chain. We have shown that the dynamics of an impurity initialized at the left edge of the system displays oscillations similar to the ones observed in Josephson junctions. Additionally, for weaker background interactions, the tunneling of the impurity can in fact be enhanced by the introduction of a barrier. This many-body effect is only possible in the presence of a background. We interpret it as the increase of the gap between the two lowest energy states, which governs the low-frequency dynamics of the system. Our results open new perspectives on the study of quantum transport in one-dimensional systems, hinting at the possibility of realizing a bosonic Josephson junction in the complete absence of an artificial barrier. Moreover, the inclusion and manipulation of double-well potentials and even lattices may allow for the optimized transfer of spin states.

Acknowledgments
The authors thank Xiaoling Cui and Artem Volosniev for reading and commenting on the manuscript, and Tomasz Sowiński for fruitful discussions. The following agencies -Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), the Danish Council for Independent Research DFF Natural Sciences and the DFF Sapere Aude program -are gratefully acknowledged for financial support.

Appendix A. Exact solutions, geometrical coefficients and spatial distributions in the limit of infinite repulsion
In this Appendix we present the solutions obtained with numerical diagonalization for the single-particle in a double-well, mapping between a strongly interacting onedimensional system and the spin chain Hamiltonian described by Eq. 2, expressions for the geometrical coefficients and densities in the impenetrable limit, and the analytical form of the eigenstates and eigenvalues of the spin chain in the 3+1 case. The singleparticle energies and numerical values for the exchange coefficients are calculated using the open-source code CONAN [29]. Some of these quantities depend on the spinless fermion wave function Φ(x 1 , ..., x N ), which is constructed as the Slater determinant of the N lowest-lying orbitals of the trapping potential. The energy E 0 from Eq. 2 is then simply the sum of the energies of these individual states.
Appendix A.1. Single-particle solutions in the double-well The eigenvalues and eigenstates of a particle in a double-well were obtained through numerical diagonalization of the hamiltonian H 0 using the 50 lowest states of the harmonic oscillator trap (b = 0) as basis. In Fig. A1 we present these solutions for different values of the barrier parameter b. In panel a), we show how each pair of states becomes degenerate as the barrier size is increased. This is reflected in the eigenstates shown in panel b): at larger values of b, the ground state and the first excited state have the same probability distribution, differing only in parity.

Appendix A.2. Mapping the strongly interacting system to a spin chain Hamiltonian
In this section we show details of the mapping between Hamiltonian 1 and the spin chain described by Eq.2. Although different approaches have been used to describe this mapping [23,32,33,26], we focus on the one presented in [24]. We start by considering a more general bosonic Hamiltonian of the form where the total number of particles is given by N = N ↑ + N ↓ . In the limit of infinite interactions (g → ∞), the eigenstates of Hamiltonian 1 can be described by where the sum runs over the L(N ↑ , N ↓ ) = N ↑ +N ↓ N ↑ permutations of the coordinates, and P k is the permutation operator. In this expression, Φ 0 is simply the wave function in the impenetrable limit, with coordinates fixed as {x ↑i } and {x ↓j }, with i = 1, ..., N ↑ and j = 1, ..., N ↓ . To investigate the behavior of the energy at very strong (but finite) interactions, we use the Hellmann-Feynman theorem, which gives where the first term on the right-hand side accounts for interactions between different bosons, while the remaining terms arise from interactions between identical bosons. The conditions for the derivatives at the contact point between two particles are given by for identical bosons and ∂Ψ ∂x ↑i − ∂Ψ ∂x ↓j for a distinguishable pair. In the expressions above we have where the denominator introduces a normalization factor. Integrating with respect to g we obtain the following energy functional where E 0 is the energy in the limit of infinite repulsion, and we neglect terms of higher order in (1/g). By introducing the wave function described by Eq.A.2 in the expression above, we obtain where a ↑↓ ik represents the coefficients in Eq. A.2 multiplying terms with nighboring ↑ and ↓ particles at position i and i + 1, while the remainig terms have the same role, for ↓↑, ↑↑ and ↓↓ pairs. The purpose of such terms is to account for the energy contribution of exchanging two neighboring particles with particular spin projections. The coefficients α i are now independet of spin, and can be written as where Φ 0 (x 1 , · · · , x i , · · · , x N ) is again the wave function present in Eq. A.2 (where we have ommited the spin indices). Since this wave function is defined in the region determined by a particular order of the coordinates, it is enough to calculate the integrals in one particular sector, such as x 1 < x 2 · · · < x N − 1. Now, let us consider a spin chain Hamiltonian defined as ) is the operator that exchanges neighboring spins with different projections and Π i,i+1 ↑↑ = Π i,i+1 ↓↓ = 1 2 (1 + σ i z σ i+1 z ) have the same action, but for identical spins. A generic spin state can now be written as where once again the sum runs over the permutations of the N ↑ and N ↓ spins. If we calculate the expected value of Hamiltonian A.11 as χ|H|χ , we obtain where the coefficients A ↑↓ ik , A ↑↑ ik and A ↓↓ ik have the same meaning as in Eq. A.9, and L(N ↑ ,N ↓ ) k=1 a 2 k introduces a normalization factor. It becomes clear that the energy functionals given by Eqs. A.8 and A.13 are identical if J i = α i /g. Furthermore, by rewriting Eq. A.11 in terms of the Pauli matrices, we obtain which is the spin chain Hamiltonian described in Eq. 2 of the main text. We conclude then that the eigenvalue problems defined with Eqs. A.8 and A.13 are identical, which validates, for a strongly interacting system, the mapping between Hamiltonians 1 and 2.

Appendix A.3. Geometrical coefficients
Different methods have been developed for calculating the coefficients in Eq. A.10, in particular exploring the determinant form of Φ 0 (x 1 , ..., x N ). Here, we use the opensource code CONAN [29] to calculate them, considering the double-well potential with different barrier parameters. An equivalent approach, based on Chebyshev polynomials, has been published in [28]. In Fig. A2, we show results for the geometrical coefficients in the cases of N = 4, 5 and 6. Due to the parity invariance of the trapping potentials considered here, we have α N = α 1 , α N − 1 = α 2 , ..., which means we must calculate, at most, N/2 coefficients in each case. In the single-well cases (b = 0), the coefficients are simply the ones expected for the given number of particles in a harmonic trap. Particularly, for N = 4 we have α 1 ≈ 1.78, α 2 ≈ 2.34. As b is increased, the major change is observed for even N in the central coefficient, which vanishes for b 3. In this limit, we have the two sides of the system described by almost completely separate harmonic traps with N/2 in each well. The values of the coefficients for harmonic traps in this situation have analytical expressions, given by α = π/2 for N = 2 and α = 3 3 /(2 3 √ 2π) for N = 3. The situation is not the same for odd N . In this case, there is no central coefficient to vanish as the barrier is increased. We would expect such a system to exhibit different spatial densities and dynamics.
Appendix A.4. Spatial correlations in the impenetrable limit In this section we describe the spatial densities for a given number of atoms N in the impenetrable limit. These densities reproduce the results expected for a Tonks-Girardeau gas or a gas of spinless fermions, and are characterized by a chain of localized atoms. The individual distributions are calculated with the following expression [32] where the integral is performed in the region Γ = x 1 < x 2 < ... < x N . The quantity ρ i (x) then gives the spatial distribution of the atom with index i. A formula for calculating Figure A3: Spatial distributions in the impenetrable limit for the case of a) N = 4 and b) N = 6.