Coherent Impurity Transport in a Binary Bose-Einstein condensate

We study the dynamics of a soliton-impurity system modeled in terms of a binary Bose-Einstein condensate. This is achieved by ‘switching oﬀ’ one of the two self-interaction scattering lengths, giving a two component system where the second component is trapped entirely by the presence of the ﬁrst component. It is shown that this system possesses rich dynamics, including the identiﬁcation of unusual ‘weak’ dimers that appear close to the zero inter-component scattering length. It is further found that this system supports quasi-stable trimers in regimes where the equivalent single-component gas does not, which is attributed to the presence of the impurity atoms which can dynamically tunnel between the solitons, and maintain the required phase diﬀerences that support the trimer state.


INTRODUCTION
Recent experimental progress has reached a state where the dynamics of a complex and thermally isolated quantum system can be studied for unprecedentedly long evolution times.In particular, advances in the field of ultra-cold atoms have allowed for such a high degree of controllability that, when combined with the absence of thermal phonons, studies of non-equilibrium coherent dynamics over timescales which are usually inaccessible in conventional condensed matter physics are possible [1,2].Not surprisingly, this has inspired a surge of theoretical interest and a growth of whole scientific communities which aim at the description of isolated, nonequilibrium, quantum systems [3][4][5][6][7].
Pioneering early experiments in this direction included the observation of the non-equilibrium dynamics of a one dimensional Bose gas (a paradigmatic integrable model) [8], which reopened foundational issues regarding thermalisation of observables in closed quantum systems [4,6,7].Perhaps the earliest experiment in this field was conducted by Greiner et al. [9], where a system was quenched across a superfluid to Mott-insulator transition and a coherent collapse and revival of the interference peaks in momentum space was observed in real time.This highly non-trivial non-equilibrium dynamics will be a central focus of this work and we aim at investigating its relationship to theoretical work which highlights the emergence of dynamical phase transitions (DPTs) in quenched dynamics.The idea of DPTs was first introduced by Heyl et al., who studied the vacuum persistence amplitude (survival probability) for certain quenches in the paradigmatic transverse Ising model [10].Through a well known mapping with the boundary partition function [11,12] they noticed that the rate function for certain quenches exhibits non-analyticities whenever the wave function becomes orthogonal to the initial state.According to Heyl et al. this behavior therefore indentifies a dynamical phase transition.Since the original inception, DPTs have been studied in a wide range of models  and while originally DPTs were believed to manifest when the system was quenched across an equilibrium phase transition, it is now known that they can occur even for quenches within the same phase [17,20,22].An exciting recent development is the observations of DPTs in experimental platforms such as ion trap architectures [39] and cold atom arrays [40].
Despite the range of models that have been investigated in relation to DPTs over the past years it is perhaps surprising that there have been little or no investigations of their manifestation in the original experiments which ignited the field, i.e. the breathing dynamics across the superfluid to Mott insulator transition [9] and dynamics in the Tonks-Girardeau gas [8].One central aim of this work is to fill that void.For this we first clarify the meaning of non-analyticities in the rate function proposed in [10] and show then that, in general, the orthogonality of the time evolved state to the initial state is not related to the temporal behaviour of local observables.Our first system of choice for this is an important continuum model, namely the Tonks-Girardeau gas [41] undergoing a pinning transition to an insulator by application of a commensurate lattice potential.This effect was first theoretically predicted by Büchler et al. [42] and later experimentally realized by Haller et al. [43].The dynamical quench problem was first studied by Lelas et al. in [44].In our calculations we provide the first evidence of periodically appearing nonanalyticities in the rate function for this process and explore the connection to the collapse/revival cycles in the dynamics of the momentum distribution.Both periodic cycles turn out to be connected only for deep quenches.
We then confirm this observation by presenting an exactly solvable discrete model which contains the same physical phenomenology i.e. hard-core bosons in a lattice at half filling with a staggered field.In this model analytic expressions can be found for the rate function and we compute the dynamics of the experimentally relevant parity operator and detail the connection with the rate function.
In the following we will first briefly review the basic ideas relating to DPTs and particular the connection with dynamical restoration of symmetry.We then first present our results for the continuum model and follow this with an in-depth discussion of the lattice model.After this we conclude with an overall discussion of some of the issues raised.

DYNAMICAL PHASE TRANSITIONS
The DPTs defined by Heyl et al. [10] are primarily centered around an object which is known as the Loschmidt amplitude and which has been exhaustingly studied under a number of guises in the past fifty years.This amplitude, following a Wick rotation z = it, can be thought of as a boundary partition function Z(z) = Ψ|e −zH |Ψ for z ∈ R [11,12].Exploiting this mapping, Heyl et al. noticed that, since the free energy density can be defined as f (z) = − lim L→∞ 1 L ln Z(z) for a system of size L, the Fisher zeros in this boundary partition function (corresponding to singularities in f (t)) coalesce into lines which can cross the real axis.This leads to the emergence of critical times t * n at which the so called rate function displays non-analyticities.According to the definition of DPTs, these singularities identify points at which the time evolved state is orthogonal to the initial one and in the following we will examine this definition for analyzing the dynamics in systems which contain a superfluid-Mott insulator transition.
It is interesting to note that in the presence of symmetry breaking one can also modify the concept of dynamical criticality as the dynamical restoration of symmetry rather than orthogonality [19].This can be seen by considering an initial condition which breaks an N s -fold symmetry of the Hamiltonian.Starting in |Ψ 0 and labeling the states obtained by repeated action of the symmetry operation as { Ψ j |} (j = 1, .., N s − 1), one can define the probability to remain in the ground-state manifold as This quantity turns out to have singularities not in the presence of temporal orthogonality but when the system crosses the boundary between two symmetry sectors.To demonstrate this let us consider for simplicity a twofold symmetry (like Z 2 ) and write according to eq. ( 2) where f 0 (t) and f 1 (t) correspond to the rate function in the two symmetry sectors.Let us now define the real valued rate function l(t) = 2 [f (t)].It is evident that at a certain time t * when the real parts of the rate function coincide, ie l 1 (t) = l 0 (t), the symmetry is dynamically restored, i.e. there is equal probability to be in both symmetry sectors.At all other times one has l 1 (t) > l 0 (t) or l 0 (t) < l 1 (t), which means that one of the two functions dominates P (t) because the L factor can be large in the exponentials.Therefore at the times t * cusp singularities appear in P (t) and a correspondence between DPTs and standard symmetry breaking in the steady state can be established [36].
It is therefore clear that great care must be taken when interpreting non-analyticities in the rate function as points of dynamical criticality.Strictly speaking such non-analytic points are times when the evolving state after the quench becomes orthogonal to the initial state.Since, in general, this has nothing to do with the restoration of a symmetry one would not expect the global orthogonality to be reflected in the dynamics of experimentally relevant observables.However, as pointed out in [10], there is a case when they can be interpreted to be the same: if the initial state is a Schrödinger cat state of the form i.e. a linear superposition of symmetry related ground-states of the initial Hamiltonian.Defining the generic rate functions via we get the Loschmidt amplitude Since in the thermodynamic limit this expression is dominated by the rate function with the smallest real part we have that i.e. the return probability calculated on a state like Eq.( 5) is equivalent to the probability to stay in the ground state manifold in the thermodynamic limit.Since the Fisher zeroes are singularities of the rate functions f jl (t), cusps in P (t) emerge when two rate functions have the same real part.The two objects therefore generally give different information about the state of the system and the question is whether this information can be extracted from local measurements or not.Indeed, P (t) can be shown to be connected to local symmetries of the Hamiltonian, since such symmetries are characterized by having local operators as generators.Furthermore, since the order parameter is an object which is in general not invariant under such local operations, the cusps in P (t) are naturally connected to zeroes of the order parameter since they indicate symmetry restoration.In turn singularities in G(t) and hence f (t) (or equivalently l(t)) indicate orthogonality.Since the ground states of Hamiltonians across a symmetry breaking phase boundary are orthogonal (in the thermodynamic limit), it is interesting to ask whether a connection between such singularities and the dynamics of local observables is present also in this case (see [45,46] for a related study of criticality in systems with long range interactions).This is what we will investigate below in the first of the two models where we focus on the emergence of non analyticities in the rate function f (t) in a highly experimentally relevant continuous model and explore their emergence with the dynamics of a measurable observable.

TEMPORAL ORTHOGONALITY IN THE TONKS-GIRARDEAU GAS
The first model we consider describes a one-dimensional system of N bosons confined in an external trapping potential.The Hamiltonian can be written as where g 1D is a parameter characterizing the sign and magnitude of the interaction and V b (x) is a box potential of length L with infinitely high walls.Let us assume an optical lattice potential of depth V 0 is applied in addition to the already existing trapping potential and is described by V (x) = V 0 cos 2 (k 0 x) where the wavevector is given by k 0 = M π/L and M is the number of wells in the lattice.When the strength of the lattice is much larger than the recoil energy, V 0 E r = ( k 0 ) 2 /(2m), the model above can be mapped onto the celebrated Bose-Hubbard model, which has a transition between a superfluid and insulating state [1].In the limit when V 0 E r , the Bose-Hubbard model is no longer applicable but interestingly it was shown in [42] that at low energies the model can be mapped on to the Sine-Gordon model and a phase transition between a superfluid and insulating state remains when the applied lattice is commensurate with the particle density.The transition was observed experimentally by Haller et al. in 2010 [43].
In this work we will consider the hard-core limit of the system, g 1D → ∞, where a pinning transition will occur for any infinitesimal lattice strength.In this limit the system is known as the Tonks-Girardeau gas, which allows for an exact solution due to the existence of the Fermi-Bose mapping theorem [41].The essence of the mapping is that the interaction term in Eq. ( 9) can be dealt with by imposing the following constraint The black solid lines are for times corresponding to the extrema of the momentum peak, while the red dotted lines are for times corresponding to the extrema of l(t).The grey solid line is the instantaneous momentum distribution of the insulating phase.
on the many-body wave-function: Ψ(x 1 , x 2 , . . ., x N ) = 0 if |x i − x j | = 0 for i = j and 1 ≤ i < j ≤ N .The system can then be mapped to free fermions subject to appropriate symmetry: ] is a Slater determinant of single particle states.This mapping theorem also holds time dependently and offers a convenient way to numerically calculate the real valued rate function l(t) = 2 [f (t)] from time evolving the single particle states in the quenched Hamiltonian H f as where the A mn (t) = ψ * m (x, 0)ψ n (x, t)dx are the matrix elements of the overlaps between the pre-and post-quench single particle states.This allows for a straightforward and numerically exact approach to the computation of the rate function.
The figure of merit we will consider is the time dependent momentum distribution n(p, t) which is routinely measured in cold atom setups.It is defined as the Fourier transform of the reduced single particle density matrix (RSPDM) where the time dependent RSPDM is ρ which is evaluated numerically using the technique developed in [47].
In the following we will study three types of quenches: switching the lattice on, switching the lattice off and changing the sign of the lattice potential.If the lattice potential is commensurate with the particle number, M = N , then switching on the lattice potential from an initial depth V i = 0 to a final depth V f > 0 allows one to observe temporal orthogonality occurring in a quench from a conducting to an insulating phase.The rate function for this quench is shown in Fig. 1(a) and non-analytic peaks can be seen to occur at times t/2 + α (where α is an integer) with a periodicity of T R = 4π/V f .In panel (b) the value of the momentum distribution at p = 0 is shown and for specific times the full momentum distribution is plotted in panel (c).The momentum distribution is initially sharply peaked at p = 0, which is characteristic for a Tonks-Girardeau gas trapped in an infinite well and which reflects the expected partial first order coherence due to the order present in the RSPDM .After the quench the sharp peak vanishes as the momentum distribution broadens, signaling the transition to the insulating phase.The magnitude of the zero momentum component therefore oscillates as the system moves between insulating and conducting phases, with the first minimum occurring at a time which is slightly earlier than the emergence of the non-analytic peak in l(t).For later times, this mismatch becomes more pronounced and the simulation clearly demonstrates that the timescale for non-analyticities in the rate function quantifying orthogonality and that for the collapse/revival cycles in the momentum distribution are close but not the same.However, the stronger the quench (V f > E R ), the more the two tend to coincide and we will explore this in more detail later when discussing the discrete model.FIG.2: Dynamics of a quench from insulator (V i = E R ) to superfluid (V f = 0) for several systems with different particle number, N .Note that the time axis is rescaled by the revival time in the box, N π/(2E r ), which has the implication that the non-analyticities will not be observed in the thermodynamic limit.
Let us now turn to the quench from insulator to superfluid, i.e. from V i > 0 to V f = 0.The behaviour of the rate function is shown in Fig. 2 for different system sizes on a time axis that is rescaled by N π/(2E r ).While one can observe a revival effect where at half the scaling time there is a type of transient criticality signalled by an apparent non-analyticity in l(t) at times t = α + 1/2 (α an interger), these non-analyticities do not signal the existence of DPTs, but rather are a result of the propagation of density waves from the box edges which then interfere at the box center.This is precisely the dynamical de-pinning effect that was studied by Cartarius et al recently in the same model [48].Therefore, this non-analyticity is the result of a finite size effect and does not exist in the thermodynamic limit.Instead the system undergoes a crossover from the insulating to the superfluid phase.This suggests that DPTs do not occur during dynamical de-pinning and we will explore this further in the discrete model in the next section.
Finally, we display in Fig. 3 the dynamics of the rate function and the momentum distribution for quenches within the insulating phase, for V i = V to V f = −V , which allows us to observe the post-quench dynamics on a timescale which is not governed by the lattice depth.Here the oscillations of both the rate function and momentum peak decay quickly, whilst it is clear that there is no simple relationship between nonanalyticities which emerge in the rate function and any features in the behavior of the momentum distribution.Let us attempt to understand in detail this phenomenology by studying a closely related exactly solvable model.

TIGHT BINDING MODEL
We consider a system of N hard core-bosons in a staggered onsite potential described by the Hamiltonian where b j are hard core bosons, J is the tunneling strength and V is the strength of the onsite potential.This model has the distinct advantage over the previous continuous model in that it is analytically solvable while retaining all the essential physics.The procedure is well known [49], using the Jordan-Wigner transformation b † j = e iπ l<j a † l a l a † j and using Fourier transformed variables, a j = 1 √ N k e −ikj a k the Hamiltonian can be written as where Ψ k = (a k , a k+π ) T and Ĥk = 2J cos(k)σ z + V σ y , where σ are the Pauli matrices.Notice that k = π(2n+1)/N , with n = 0, . . ., N/4 − 1.The Hamiltonian can be diagonalized in terms of the new variables Γ k = e iθ k σ y Ψ k , where tan(2θ k ) = V /(2 cos(k)), and the resulting spectrum is characterized by a dispersion k = (2J cos(k)) 2 + V 2 .For our purposes we will work at half filling where the spectrum is always gapped unless V = 0, in which case the gap at k = ±π/2 closes.Hence for V = 0 we have an insulating charge density wave phase while for V = 0 it is a "superfluid".In what follows we will consider three different types of quenched dynamics as we did in the previous section: quenches from the superfluid to insulator, quenches from the insulator to superfluid and then quenches within the superfluid phase.We note that the same model can also be solved in the presence of an external flux [50].
FIG. 4: l(t) vs. t/N for a quench from the insulator (V i = 0.3) to the superfluid (V f = 0) for different system sizes.
Fixing the tunneling strength J = 1 and considering a general quench from V i to V f , the Loschmidt amplitude can be computed using the Bogoliubov rotation connecting the old to the new quasiparticles and computing the time evolution one finally obtains Recalling that the Fisher zeroes are the roots of this complex valued function, one can solve them for G(z k ) = 0 and find the expression For quenches towards the insulating phase (V f > 0) it is evident that the Fisher zeros hit the real axis, hence corresponding to zeros of the Loschmidt amplitude (singularities of l(t)) A singularity at these momenta corresponds to a singularity in the rate function with a period For quenches towards the superfluid phase Fisher zeroes have always a finite imaginary part implying the absence of singularities in f (t) and therefore no DPTs are observed.However, keeping the system size finite and rescaling the time by it, one can observe a nice collapse and revival picture (see Fig.( 4)), as we previously discussed in the continuous model in Fig. 3(a).

ORTHOGONALITY AND OBSERVABLES
As discussed above, singularities in the rate function signal zeros at times when the time evolved state becomes orthogonal to the initial one.We now gain further understanding of why this occurrence is related to the time evolution of physical observables only for deep quenches.Notice that according to the calculation performed above, the overlap between the different ground states of the Hamiltonian Eq. ( 13) at different strengths of the staggered potential, V i and V f , is given by Therefore different ground states turn out to be orthogonal in the thermodynamic limit, with the overlap vanishing exponentially in the system size.This is suggestive, since if upon quenching V say from the superfluid V i = 0 to the insulator V f = 0 the system dynamics would result in consecutive collapses and revivals of the superfluid into the insulator, one could expect the system to attain orthogonality with the initial state at the farthest point from the superfluid, i.e. when the collapse into the insulator is complete.This intuition would be correct if the system would be able to dissipate the work done on it by the quench procedure.In the the present case of unitary dynamics, however, the fact that the system remains in a superposition of excited states of the post-quench Hamiltonian, makes the identification of the phenomena problematic.In other words it is only in the thermodynamic limit that ground states with different parameters are orthogonal.Hence only in that limit one could expect that, if the system was indeed able to collapse and turn back from one state to the other, one would get orthogonality when the superfluid fully collapses into a Mott insulator.An exception are deep quenches as we will now show.
In order to distinguish between the superfluid and the insulating phase, we choose the experimentally accessible parity operator which is an observable that characterizes charge density wave order In the fermionic representation this is given by The calculation of M gives Plotting this function in general in a situation where singularities in the rate function are present shows that while both quantities oscillate the time scales are typically very different (see Fig( 5)).There is however one instance in which the two quantities appear to have a correlated behavior (see Fig.( 6)), i.e. for quenches from V i ≤ 0 to a large V f > 0 (V f /J 1).In this case the parity operator oscillates between a positive and negative value periodically and each time a minimum is attained a cusp singularity is observed.This result is however simple to understand: the period of the oscillations of M (t) The two are clearly equal if V f V i , J in which case the dispersion is effectively flat and all k-modes oscillate with the same frequency.Therefore only in this case the orthogonality appears to be tied to oscillations of the order parameter.One might be tempted to argue that this is just the wrong operator to detect orthogonality.If however a different operator is used, such as for example the kinetic energy operator, it is easy to show that after a quench (2 cos(k)) cos(2θ k ) cos(2∆θ k ) + sin(θ k ) sin(∆θ k ) cos(2 k t) , which produces results similar to the ones presented above.

DISCUSSION
In this paper we have undertaken an extensive study of dynamical criticality in systems which contain a superfluid-Mott insulator transition in equilibrium.Furthermore, to our knowledge, this is the first numerical study of this type in a continuum model.We have found that although non-analyticities are present for certain quenches which signal temporal orthogonality this is only manifested in experimentally relevant observables for deep quenches.We studied numerically the dynamics of both the rate function and the momentum distribution following a quench in the Tonks-Girardeau gas across the pinning transition.In the discrete case we provided analytic calculations for the rate function and the dynamics of the parity operator.As known from state discrimination protocols in quantum information, it is an extremely difficult task to uncover global orthogonality from local measurements on pure states [51] and in the case of mixed states it is generically impossible [52], so we are lead to conjecture that in general it is not possible to detect orthogonality in the dynamics of the many-body state and hence non-analyticities in the rate functions by observing the dynamics of local observables alone.Nevertheless, we stress that one could still hope to detect such points through non-trivial order parameters [53] or perhaps even by extending ancilla based interferometry schemes which have been proposed [54][55][56][57] and experimentally implemented in local quenches in Fermi gases [58].In addition, studying the dynamics of the rate function and these experimentally relevant observables for quenches in critical models is interesting in its own right and we hope it will inspire further experiments in this direction.

FIG. 1 :
FIG. 1: Dynamics of a quench from the conducting (V i = 0) to the insulating phase (V f = E R ) for a system of N = 100 particles.(a) The rate function and (b) the height of the momentum peak as a function of time which is scaled with respect to T R .The solid vertical lines indicate the times that the non-analyticities appear in l(t), while the dashed vertical lines indicate the minima of l(t).For the times marked by symbols in (a-b) the momentum distribution is plotted in (c).The black solid lines are for times corresponding to the extrema of the momentum peak, while the red dotted lines are for times corresponding to the extrema of l(t).The grey solid line is the instantaneous momentum distribution of the insulating phase.

FIG. 3 :
FIG. 3: (a) The rate function and (b) the height of the momentum peak for N = 100 particles after a quench from the insulating phase (V i = 2E R ) to the insulating phase (V f = −2E R ).The solid vertical lines indicate the times at which the non-analyticities appear in l(t), while the dashed vertical line indicates its minima.

FIG. 5 :
FIG. 5: (a) l(t) and (b) M (t) vs. t for N = 100 and a quench fromV i = −1/6 to V f = 3.The correspondence between minima of M (t) and cusps of l(t) is not present in this case.