Excitations in two-component Bose-gases

In this paper, we study a strongly correlated quantum system that has become amenable to experiment by the advent of ultracold bosonic atoms in optical lattices, a chain of two different bosonic constituents. Excitations in this system are first considered within the framework of bosonization and Luttinger liquid theory which are applicable if the Luttinger liquid parameters are determined numerically. The occurrence of a bosonic counterpart of fermionic spin-charge separation is signalled by a characteristic two-peak structure in the spectral functions found by dynamical DMRG in good agreement with analytical predictions. Experimentally, single-particle excitations as probed by spectral functions are currently not accessible in cold atoms. We therefore consider the modifications needed for current experiments, namely the investigation of the real-time evolution of density perturbations instead of single particle excitations, a slight inequivalence between the two intraspecies interactions in actual experiments, and the presence of a confining trap potential. Using time-dependent DMRG we show that only quantitative modifications occur. With an eye to the simulation of strongly correlated quantum systems far from equilibrium we detect a strong dependence of the time-evolution of entanglement entropy on the initial perturbation, signalling limitations to current reasonings on entanglement growth in many-body systems.


Introduction
One of the key proposals in the field of quantum computing, made by Feynman, is the idea to use one quantum system to simulate another one in order to circumvent the problem of the qualitatively different complexity of quantum systems and classical computers as the standard simulation tools of science. The advantage of Feynman's approach lies in the fact that the new simulating quantum system may have advantages in experimental control: preparation, manipulation of parameters and measurement.
Over the last decade, this proposal has been filled with life by the progress made in the preparation of dilute ultracold atom gases, highlighted by the now almost routinely achieved preparation of Bose-Einstein condensates. The use of Feshbach resonances or optical lattices has given us unprecedented control over interaction and dimensionality in strongly interacting quantum many-body systems of unique purity. This has been put to use in the creation of strongly correlated quantum systems that have been in the focus of interest in condensed-matter physics for a long time [1]. Examples are the observation of the superfluid to Mott insulator transition for Bose gases [2] and the fermionization of strongly interacting one dimensional bosons [3,4].
But one can do more: it is possible to create physical systems of their own interest that have no counterpart in conventional condensed matter physics. To take one example, due to the internal spin degree of freedom of electrons it is quite natural in solids to consider models of two components of fermions. The existence of two fermionic components is the driving force of important phenomena such as collective magnetism.
In this paper, we consider the bosonic equivalent of the two-component fermionic system. Strongly interacting two-component bosonic systems have no counterpart in condensed-matter physics, which is why they have found limited attention in the solidstate literature. Yet, they are quite easily implemented in the field of ultracold atom gases [5,6], and have generated quite some interest (for a review, see [7]). We focus on the case of one-dimensional two-component bosonic systems: on the one hand, they offer interesting physical phenomena such as a bosonic version [8,9,10] of fermionic spincharge separation which was discussed for cold atoms by [11,12,13,14]. On the other hand, very good control exists both analytically via bosonization [15] and numerically via (time-dependent) DMRG [16,17].
In this paper, we start out by considering the low-energy physics which is characterized, as all other critical one-dimensional quantum systems, by a few effective Luttinger liquid parameters, which we determine both by a mapping in a limiting regime and more generally by DMRG. In principle, the Luttinger liquid parameters determine completely the static and linear response properties of the system.
We move on to discuss the spectral functions which are the cleanest way to observe Luttinger liquid physics, in particular spin-charge separation and characterize the linearresponse behaviour of two-component bosonic systems. These are obtained in the framework of dynamical DMRG.
Experimentally, one will be confronted by certain limitations of ultracold atomic systems: currently, spectral functions are unavailable and one is restricted to monitor the time-evolution of excitations. These in fact show a separation of symmetric and antisymmetric density combinations ("charge" and "spin") of the individual species from which results in good agreement with those from the spectral functions can be derived. These observation still hold if we also take into account that in current implementations of two-component bosonic systems there are slight differences in the intraspecies interactions of the two components. Moreover, the presence of a confining harmonic trap potential does not qualitatively alter the results, as we can show by explicit simulation using time-dependent DMRG.
Last, but not least, we turn to the discussion of the time-evolution of the entropy of entanglement in the various out-of-equilibrium scenarios considered in this paper. This study was originally motivated by the fact that the efficiency of DMRG simulations is limited by entanglement growth which impacts exponentially on the numerical resources needed. Here, it turns out that there are surprisingly large variations in the generally accepted scenario of an essentially linear entanglement growth after quenches (as which all our scenarios can be interpreted). While we cannot give a general explanation for the phenomenon, we hope to provide a stimulus for further research to develop a more complete understanding of entanglement evolution.

Model
Cold atomic gases with two hyperfine species confined in optical lattices can be described by the two-component Bose-Hubbard model given by the following Hamiltonian [18]: where j is a site index and ν labels the two different flavours of bosons. J is the hopping strength, U ν the intraspecies and U 12 the interspecies onsite interaction. Here we assume that the two hyperfine species have the same mass and see the same lattice potential. Unless otherwise stated we use U 1 = U 2 ≡ U, and u = U/J, u 12 = U 12 /J. ε j,ν describes an external potential, i.e. given by a trap. In the following we denote the lattice spacing by a. We mainly consider incommensurable fillings with equal densities n 1 = n 2 ≡ n smaller than one. For vanishing interspecies interaction we recover the one component Bose-Hubbard model being in a superfluid phase. The superfluid phase remains stable for finite U 12 up to U 12 ∼ U. For U 12 > U the interspecies interaction becomes dominant and a demixing of the flavours occurs, i.e. a phase separation [8,19]. The physics corresponds then to the one of an itinerant ferromagnetic system with non-Luttinger liquid properties [20,21]. Experimentally, to test for Luttinger liquid properties, the relevant regime is given by U 12 ∼ U, with U 12 slightly below U to avoid the phase-separation regime. Strongly differing interaction parameters have not been realized experimentally so far; our typical interaction parameters are chosen accordingly.

Approximations
Continuum model In a weakly interacting superfluid phase in the low filling limit, the Bose-Hubbard model can be mapped to the continuous Lieb-Liniger model [22,23]. The mapping can be performed taking the limit a → 0 while leaving Ja 2 constant. The Hamiltonian for two bosonic species in the continuum is Here Ψ ( †) is the bosonic annihiliation (creation) operator, V the external potential, M the mass of the particles, and g and g 12 are the strengths of the intra-and interspecies interaction, respectively. The parameters of the continuum model and the lattice model are related by Ja 2 = 1 2M , the interaction strength to the δ-interaction strength by Ua = g and U 12 a = g 12 , and the density ρ to the filling factor n by ρa = n. In the hydrodynamic approximation [24] the sound velocities of this model are given by with v 0 = gρ/M.
The indices stand for charge and spin, corresponding to symmetric (charge) and antisymmetric (spin) combinations of the two bosonic fields. At the symmetric point, where the interspecies interaction is equal to the intraspecies interaction, the model is SU(2) symmetric and can be solved by the Bethe ansatz. For the special case g 12 = g the sound velocity has been determined using the Bethe ansatz [21,25]. The spin dispersion becomes quadratic at this point.
Bosonization The low energy physics of the Bose-Hubbard model can be described by the bosonization approach [15]. Two bosonic fields, θ ν and φ ν , are introduced in the continuum which are related to the phase and the amplitude of the original bosonic operator, respectively. To be more precise, in this representation the bosonic creation operator becomes b † ν (x) = √ ρ 0 p e i2p(πρ 0 x−φν (x)) e −iθν and the density operator Here ρ 0 is the average density, and 1 π ∂ x φ ν and θ ν are conjugate operators. The advantage of this representation is that the low-energy properties of the system are described by a quadratic Hamiltonian.
Introducing the "charge" and "spin" degrees of freedom by φ c = 1/ √ 2(φ 1 + φ 2 ) and φ s = 1/ √ 2(φ 1 − φ 2 ) it separates into two different part and is given by Here the parameters v ν are the velocities and K ν the Luttinger liquid parameters for the corresponding field. Due to this mapping the asymptotic physics of the Bose-Hubbard model is totally determined by the velocities and the Luttinger liquid parameters. In the case of intermediate interaction strengths, the relations between the microscopic parameters of the Bose-Hubbard model and the parameters of the bosonization approach are difficult to establish analytically. In the limit of small interaction strength and low filling they are given for the single species Bose-Hubbard by Janπ = v 0 K and Ua/2π = v 0 /K. In the limit of weak interspecies coupling the parameters for charge and spin degrees of freedom can be related to the parameters of the system without interspecies coupling by the perturbative expressions Outside the range of validity of these two approximations the parameters can be determined numerically which we do in the following for a broad range of parameters. The Luttinger parameter is determined by using The compressibility [15] is given by , where E 0 (L) is the energy of the ground state for a system of length L and the combined "charge" and "spin" densities are denoted n c,s = n 1 ± n 2 . The derivative can be computed numerically, yielding E(N 1 , N 2 ) is the ground state energy of the system with N 1 particle of the first flavour and N 2 of the second one. After a L → ∞ extrapolation, we can obtain the Luttinger parameter by using Eq. (9) and the numerically determined velocities. Note that the formulae given for K s and κ s in Equations (9) and (10) differ by (mutually cancelling) factors of 4 from the usual relationships in electronic models, where n s = (n 1 − n 2 )/2; the factor 1/2 is not meaningful in the model considered here and hence dropped. The final value of K s is of course unaffected. Fig. 1 shows the dependence of the Luttinger parameters on the interspecies interacting strength u 12 . The charge Luttinger parameter K c decreases with increasing u 12 , the spin Luttinger parameter increases. As expected, the bosonization result agrees with the numerical results for small u 12 and starts to deviate with increasing u 12 . The spin Luttinger parameter shows stronger deviations. Fig. 2 shows the Luttinger parameter for fixed interaction strength (u = 3, u 12 = 1.2, 2.7) and varying density n. Note that for the spin Luttinger parameter u 12 v 0 /(πK 0 ) is of the order of one, and thus the perturbative result (Eq. 8) is not reliable anymore. The DMRG results deviates significantly from the perturbative results; the deviations increase with density. A similar effect also occurs in the one-component Bose-Hubbard model where for densities above 0.5 lattice effects becomes relevant, see the Appendix for more details.

Spectral Functions
The key quantity to characterize single particle excitations in many-body systems is the dynamic single-particle spectral function, because it can be probed easily in solid-state setups. The single particle spectral function is defined as where |0 is the ground state with energy E 0 . For fermionic two-component Hubbard model this function signals the separation of spin and charge degrees of freedom by two distinct peaks at different frequencies.
In the following we determine this function for the two-species bosonic model using a variant of dynamical DMRG [26,27]. Our formulation of the algorithm is entirely based on matrix product states. This allows to avoid the targetting of multiple states in the DMRG algorithm, which is achieved only at substantial numerical cost. The new formulation of the algorithm saves more than an order of magnitude in time. As always, DMRG prefers open boundary conditions leading to spurious effects in the spectral functions; in this context, filtering procedures have been used [26]. In order to reduce boundary effects we use the quasi-momentum definition of the Fourier transform b q,ν = j sin qjπ L+1 b j,ν with the momentum |k| = qπ/(L + 1). These states create a single particle basis of eigenstates of the non-interacting system (u = u 12 = 0) with open boundary conditions. They have nodes at the edges of the systems and are therefore better suited to open boundary conditions than the usual definition of the Fourier transform using the eigenfunctions of the non-interacting system with periodic boundary conditions. Dynamical DMRG uses a correction vector method to calculate the spectral function: First define |lv(q, ω) = b † q,1 |0 . Then the correction vector is given by |cv(q, ω + iη) = 1 E 0 +ω+iη−Ĥ |lv(q, ω) and thus obey This complex equation for the correction vector can then be solved using the GMRES algorithm. The value of the spectral function is thus where the last line is only valid if |lv(q, ω) = |lv(q, ω) and |Ψ denotes the complex conjugate of |Ψ . Using the derivative, we can determine the spectral function with a smaller number of correction vectors and thus much more efficiently.
Note that in the following only the normalized spectral function is used, such that dω A norm (q, ω + iη) = 1 holds for every q. The full spectral function A norm (q, ω) is shown in Fig. 3. For the used system parameters L = 64 and η = 0.1 we needed up to 2000 states for each correction vector. One clearly sees two different branches with a linear dispersion relation ω ≈ v c,s q yielding two different velocities.
Thus this system exhibits spin-charge separation. Fig. 4 gives a more detailed view of the spectral function. For a number of given momenta q = kπ/(L + 1) the spectral function is plotted versus ωq −1 . Therefore the norm of the scaled spectral function obeys d (ωq −1 ) A norm (q, ω + iη) = q −1 and the norm decreases with increasing q. The position of the peak is roughly given by ω ≈ v c,s q, thus one expects two peaks at the two velocities v c and v s [28]. The spectral function shows a number of boundary effects similar to the one encountered for a single-component Bose-Hubbard model, see the Appendix for a detailed discussion of the one-component model. Here an addition peak at ω ≈ 0 occurs, visible for instance for q = 15π/65. Furthermore, for small values of q the charge peak is shifted to lower frequencies and above the charge peak an additional shoulder appears (i.e. for q = 10π/65). There are also some effects which do not occur for a single-component model: for large momenta the charge peak splits up and an additional peak below the main peak emerges. The plot also shows some effect of a nonlinear dispersion. Both spin and charge peak are shifted to higher frequencies for higher momenta, similar to the effects visible for a Bose-Hubbard model, see Fig. A2. In Fig. 3 momenta lower than q ≈ 0.25 are not plotted, because for the system sizes considered, finite size effects dominate there and obscure the physically relevant information.

Density perturbations and single-particle excitations
In ultracold atom experiments, spectral functions are hard to observe, but using suitably tuned and focused lasers it is easy to create local density perturbations. Theoretically the density perturbations can be created applying an external potential of the form ǫ j,cs ∼ exp(− (j − j 0 ) 2 /2σ j ) for times t < 0. For times t > 0 the potential is switched off.    Time-evolutions can be calculated by adaptive time-dependent DMRG most easily by using Trotter decompositions of short-ranged Hamiltonians leading to local timeevolutions [29,30,31]. For longer-ranged interactions or systems with large local state spaces, it is more efficient to consider global time-evolutions [32]. As we are dealing locally with products of two bosonic state spaces, we calculated the time-evolution of the density perturbations numerically in the latter framework. The algorithm was formulated using matrix product states and the global time evolution was done using a Krylov algorithm for exponentiation [33] with a fixed error bound per timestep [34]. The used error bound |Ψ(t + ∆t) − exp(−iĤ∆t)|Ψ(t) 2 is of the order of 10 −5 with a timestep of ∆t ≈ 0.1 − 0.2. Usually, 6 to 10 Krylov vector were used. For the Krylov vectors up to 3000 states were used in the case of the density perturbations. The time-evolutions of single particle excitations are much harder to perform, up to 7500 states were used there.
Snapshots of the density evolution of these excitations for one parameter set are shown in Fig. 5. One sees that the charge and spin perturbation which is created at time t = 0 splits up into two counter-propagating excitations. The velocity of the spin perturbation is much lower than the velocity of the charge perturbation. Even after separating into two perturbations the amplitude shows a decay. This decay is very slow for weak inter-species interaction (see Fig. 5). However, if the inter-species interaction is approaching the value of the intra-species interaction a strong broadening of the spin perturbation can be seen (see Fig. 6). Due to the very low velocity and the broadening the two peaks only separate at very long times. The strong broadening of the spin peak is expected since in the limit of equal interaction strength the system has a quadratic spin dispersion relation. Therefore the region of the linear regime of the dispersion for u 12 < u decreases if u 12 gets closer to u. The dips in front of the spin perturbation and behind the charge perturbations are due to the remaining interaction between the spin and the charge degrees of freedom and finite-size effects.
The results for a density perturbation connect directly to a setup where the time- evolution of densities is followed upon the creation of single-particle excitations in bosonic gases instead of a density perturbation. In a one-dimensional system these single particle excitations decay into a charge excitation and a spin excitation. In contrast in a higher-dimensional system single particle excitations have a finite life-time. Since in cold atomic gases one can realize systems of different dimensionality, this second setup, if realized, would give the possibility to directly confront the decay in a one-dimensional system with the finite life-time in a three dimensional system. Fig. 7 shows the time evolution of the density and bipartite entanglement entropy profiles for a creation and an annihilation of a particle at time t = 0. The initial excitation splits up into a right and a left moving part and one can clearly see the two different velocities for the charge and spin density. This eventually leads to the separation into a spin and a charge excitation. The entropy profiles exhibit a significant difference between the creation and the annihilation of a particle. After removing one particle the entropy stays almost constant (up to some small waves) with respect to time, the entropy strongly increases between the two counterprogating peaks after adding an addional particle. The reason for this behaviour is unknown; observations and conjectures will be discussed in more detail in Sec. 7.

Experimental constraints: Effects of unequal intraspecies interaction strengths and confining trap potentials
An experimental realisation of a two-component Bose-Hubbard model is given by the use of two different hyperfine states of 87 Rb, for instance |F = 2, m F = −1 and |F = 1, m F = 1 . The s-wave scattering lengths for these states are approximatly a 2 = 91.28a B , a 1 = 100.4a B , where a B is the Bohr radius [35]. Hence the intraspecies interactions are slightly distinct for the two bosonic species. By comparison, the interspecies scattering length a 12 is of the same order of magnitude and can be tuned by a Feshbach resonance [5,35]. Treating unequal intraspecies interactions in the bosonization approach results in a coupling of the spin and the charge part of the Hamiltonian. To decide if for the experimentally relevant parameters the coupling is already important we have calculated the time-evolution of a single particle excitation of a system with realistic parameters. The parameters have been determined using a lattice depth of V 0 = 4.3E R and a interspecies scattering length a 12 = 80a B , where E R is the recoil energy of the optical lattice. This yields the following Bose-Hubbard parameter, U 1 /J = 2.983, U 2 /J = 2.712, U 12 /J = 2.377. Fig. 8 shows the density profiles for different times. Even though a small coupling between the spin and the charge degree of freedom might be present, the single particle excitation splits up into a charge and a spin excitation. Therefore the separation in spin and charge survives a small experimental mismatch in the interspecies interaction strength.
A further complication of the ultracold atomic gases setup is the presence of an parabolic trapping potential. In the Bose-Hubbard model this trapping potential can be described by adding We calculated the time-evolution of a single particle excitation with a trapping potential, see Fig. 9. In contrast to the case without a trap the ground state density is not constant anymore. Therefore the velocity of the charge and the spin excitation now depend on the spatial position of the excitation. Still, two counter-propagating waves can be observed and the spin-charge separation is not qualitatively changed. The effect of a trapping potential has already been discussed in the context of spin-charge separation for two component fermionic systems [11,12,13], where similar robustness of the results was found.

Entropy of entanglement
The superposition of states that is characteristic of quantum mechanics implies for many-body systems the phenomenon of entanglement that is the key deviation of the quantum from the classical world and a key resource of quantum computing. In this Section, we focus on one measure of entanglement in bipartite systems, the entropy of entanglement, which is defined in the case of pure quantum states for an arbitrary bipartition of the system into a left part A and a right part B by cutting at bond i by forming the reduced density operators of A and B, The entropy of entanglement is then given by the von Neumann entropy of either reduced density operatorρ A orρ B : The identity of both definitions follows from the well-known observation that the eigenspectra of reduced density operators of a bipartition are identical. Entropy of entanglement has been studied in numerous contexts. In the present work, the focus is on the time-evolution of the entropy of entanglement in an out-ofequilibrium setting. On the one hand, this question is of fundamental interest for the understanding of coherent out-of-equilibrium quantum dynamics. On the other hand, it turns out that the time-evolution of the entropy of entanglement is closely related to the performance of time-dependent DMRG and related methods: the number of states (matrix dimensions) needed in such simulations are in a roughly exponential relationship with the entropy of entanglement. This means that e.g. a linear growth of entanglement entropy in time is reflected in an exponential growth in matrix dimensions, yielding the key limitation [36] for the time-scales accessible for such algorithms.
Building on the Lieb-Robinson theorem [37] it has been shown by Osborne [38] that for arbitrary short-ranged Hamiltonians in one dimension, entanglement growth is bounded linearly in time, S(t) ≤ S(0) + ct, reflecting a finite speed of propagation in such Hamiltonians, implying a potentially exponentially growth of matrix dimensions in time. In fact, such linear growth of entanglement has been observed [39,40,41] and been traced back to the fact that ouf of equilibrium quantum states show excitations propagating through the system, leading to a linearly expanding "light cone" where information is exchanged between bipartition parts A and B.
In this paper, we are considering three quite different types of time-evolution: (i) time-evolution after an insertion of a single particle, (ii) time-evolution after the extraction of a single particle, (iii) time-evolution after removing an external potential that created a density perturbation. All processes should generate excitations moving at the same maximal finite speed of propagation; one would therefore naively expect that all cases lead to a qualitatively similar linear growth of entanglement entropy in time, although we cannot expect identical results: for example, there is no particle-hole symmetry relating (i) and (ii). Fig. 7 shows the time-evolution of the entropy profiles after a single particle excitation by either insertion (left) or removal (right), i.e. cases (i) and (ii) where the entanglement entropy has been measured for all possible bipartitions. The drops at the ends are a natural consequence of the small dimension of eitherρ A orρ B there. The difference between (i) and (ii) is striking. In the case (i) entanglement entropy increases roughly linearly in time, as expected. However, in case (ii) entanglement entropy stays almost constant. This implies that simulating case (i) is much harder numerically than case (ii). For a comparison, we show case (iii), a small density perturbation evolving in time: the entanglement entropy is almost entirely unchanged under the time-evolution, see Fig. 10.
To summarise these observations, we show the maximal entanglement entropy on a chain versus time for cases (i)-(iii) in Fig. 11, with the ground state entanglement given as a reference. As already seen in Fig. 7, a enormous difference between the creation and the annihilation can be observed.
We can exclude that these observations are due to limitations of the numerical method that obviously neglects some of the entanglement entropy due to its inherent truncations. However, these results are converged in the used matrix dimensions; we also have calculated the fidelities between the initial states and the result of time-evolutions up to time t and then back to 0, finding fidelities that are essentially 1. The entanglement entropies shown can therefore be considered exact.
In the case of a weak density perturbation, we attribute the observed behaviour to the fact that the out-of-equilibrium wave function is essentially a weakly distorted ground state wave function that does not affect the entanglement structure of the ground state substantially. It is much less obvious to explain the difference between (i) and (ii) despite the absence of a particle-hole symmetry. If we consider Fig. 7, in the case of particle creation the increase of entanglement entropy is limited to the regions into which spin and charge excitations have already propagated that link ever larger regions of space. This is as expected. In the case of particle annihilation we might suppose that due to the relatively low density entanglement is removed because the chain is to some extent cut by the removal of a particle. In particular, if one thinks of the initial state as a superposition of many different Fock states, the application of the annihilation operator eliminates all states in which no particle occupies the site. Thereby the entanglement in the system is reduced. Indeed, a small dip is observed. The perturbation in the entanglement entropy propagates again linearly with the charge and spin excitations, but quantitatively it stays essentially unchanged. A possible explanation is suggested by perturbation theory. For short times the time-evolution operator can be approximated by an expression of the structure 1 + itJ (b † i b i±1 + h.c.) + it(U/2)n i (n i − 1), where we have suppressed the multi-component nature of the problem for illustrative purposes. While the time-evolution operator applied far from the perturbation conserves the state and thereby the entanglement entropy, we see that close to the perturbation changes occur. At the densities considered, the site where the annihilation operator is applied has  S t single particle excitation (b) single particle excitation (b + ) small density perturbation groundstate Figure 11. (Color online) Maximum of the bond entropy versus time for single particle excitations (both b † |0 and b|0 ) and for a small density perturbation. The small wave occuring for the b † |0 curve are an effect of the lattice. Note the enormous difference between the two single particle excitations. The system parameters are n 1,2 = 0.625, u = 3, u 12 = 2.1.
most likely occupation number 0 or 1 for the respective species in the contributing Fock states. After application of the operator, the Fock states either vanish or have occupation number 0 on the site, such that only one of the terms in the hopping operator couples to the state at the site of the perturbation. Rerunning the argument for the creation operator, no contribution vanishes (if we assume intermediate interaction strengths) and all hopping terms couple, which might indicate that changes in entanglement are much more pronounced in the latter case, as observed. This is obviously not rigorous at all, and at the moment we can only conclude that the simple picture of excitations transporting entanglement leading to linear entanglement growth needs substantial refinement.  boundary effects are visible: At ω ≈ 0 there exists an additional peak and the main peak is shifted to lower frequencies. Above the main peak there is also an additional shoulder. These are the same effects as we find for the two-component model. Let us note that all boundary effects become less pronounced for larger values of q.