Spectral properties and breathing dynamics of a fewbody BoseBose mixture in a 1D harmonic trap
Abstract
We investigate a fewbody mixture of two bosonic components, each consisting of two particles confined in a quasi onedimensional harmonic trap. By means of exact diagonalization with a correlated basis approach we obtain the lowenergy spectrum and eigenstates for the whole range of repulsive intra and intercomponent interaction strengths. We analyse the eigenvalues as a function of the intercomponent coupling, covering hereby all the limiting regimes, and characterize the behaviour inbetween these regimes by exploiting the symmetries of the Hamiltonian. Provided with this knowledge we study the breathing dynamics in the linearresponse regime by slightly quenching the trap frequency symmetrically for both components. Depending on the choice of interactions strengths, we identify 1 to 3 monopole modes besides the breathing mode of the center of mass coordinate. For the uncoupled mixture each monopole mode corresponds to the breathing oscillation of a specific relative coordinate. Increasing the intercomponent coupling first leads to multimode oscillations in each relative coordinate, which turn into singlemode oscillations of the same frequency in the compositefermionization regime.
I Introduction
The physics of ultracold atoms has gained a great boost of interest since the first experimental realization of an atomic BoseEinstein condensate Anderson et al. (1995); Davis et al. (1995), where research topics such as collective modes Jin et al. (1996); Mewes et al. (1996); Stringari (1996), binary mixtures Myatt et al. (1997); Hall et al. (1998) and lowerdimensional geometries Olshanii (1998); Görlitz et al. (2001); Moritz et al. (2003) were in the focus right from the start. In most of the experiments on ultracold gases the atoms are but weakly correlated and well described by a meanfield (MF) model, the wellknown GrossPitaevskii equation (GPE), or in case of mixtures by coupled GPEs Ho and Shenoy (1996); Esry et al. (1997); Pu and Bigelow (1998a). BoseBose mixtures exhibit richer physics compared to their single component counterpart. For instance, different ground state profiles can be identified depending on the ratios between the intra and interspecies interaction strengths, being experimentally tunable by e.g. Feshbach resonances (FR) Chin et al. (2010): the miscible (M), immiscible symmetrybroken (SB) or immiscible coreshell structure, also called phase separation (PS) Papp et al. (2008); Tojo et al. (2010); McCarron et al. (2011). Comparing the experimentally obtained densities to numerical MF calculations Pattinson et al. (2013); Lee et al. (2016) provides a sensitive probe for precision measurements of the scattering lengths or, if known, the magnetic fields used to tune them Tojo et al. (2010). Another possibility to access the interaction regime and thus the scattering lengths is by exciting the system and extracting the frequencies of lowlying excitations Egorov et al. (2013). In contrast to a singlespecies case the collective modes of mixtures exhibit new exciting phenomena: doublet splitting of the spectrum containing inphase and outofphase oscillations, modesoftening for increasing intercomponent coupling, onset of instability of the lowest dipole mode leading to the SB phase as well as minima in the breathing mode frequencies w.r.t. interaction strength Pu and Bigelow (1998b); Gordon and Savage (1998); Morise and Wadati (2000).
The breathing or monopole mode, characterized by expansion and contraction of the atomic density, has in particular proven to be a useful tool for the diagnostics of static and dynamical properties of physical systems. It is sensitive to the system’s dimensionality, spin statistics as well as form and strength of interactions Bauch et al. (2009, 2010); Abraham et al. (2012); Abraham and Bonitz (2014). In the early theoretical investigations on quasionedimensional singlecomponent condensates Menotti and Stringari (2002) it was shown that different interaction regimes can be distinguished by the breathing mode frequency, which has been used in experiments Moritz et al. (2003); Haller et al. (2009); Fang et al. (2014). Furthermore, the monopole mode provides indirect information on the ground state McDonald et al. (2013), its compressibility Altmeyer et al. (2007) and the lowlying energy spectrum such that an analogy has been drawn to absorption/emission spectroscopy in molecular physics Abraham and Bonitz (2014).
From a theoretical side, those of the above experiments which are concerned with quasi1D setup are in particular interesting, since correlations are generically stronger, rendering MF theories often inapplicable. Here, confinement induced resonances (CIR) Olshanii (1998) can be employed to realize the TonksGirardeau limit Kinoshita et al. (2004); Paredes et al. (2004), where the bosons resemble a system of noninteracting fermions in many aspects. While this case can be solved analytically Girardeau (1960); Girardeau and Minguzzi (2007), strong but finite interactions are tractable only to numerical approaches, which limits the analysis to fewbody systems. For instance, a profound investigation of the ground state phases of a fewbody BoseBose mixture Zöllner et al. (2008); Hao and Chen (2009) showed striking differences to the MF calculations: for coinciding trap centres, a new phase with bimodal symmetric density structure, called composite fermionization (CF), is observed while SB is absent for any finite intercomponent coupling. Only in the limit of infinite coupling the ground state becomes twofold degenerate enabling to choose between CF and SB representations Dehkharghani et al. (2015), while the MF theory predicts the existence of SB already for finite couplings. This observation accentuates the necessity to include correlation effects.
In this work we solve the timeindependent problem of the simplest BoseBose mixture confined in a quasi1D HO trap with two particles in each component, covering the whole parameter space of intra and interspecies interactions, thereby complementing the analysis of some previous studies GarcíaMarch et al. (2014a, b); Dehkharghani et al. (2015). To accomplish this, an exact diagonalization method based on a correlated basis is introduced. We unravel how the distinguishability of the components renders the spectrum richer and complexer compared to a single component case Zöllner et al. (2007). Furthermore, these results are used to investigate the breathing dynamics of the composite system. While the breathing spectrum of a single component was recently investigated comprehensively in Schmitz et al. (2013); Tschischik et al. (2013); Chen et al. (2015); Choi et al. (2015); Gudyma et al. (2015); Atas et al. (2016), reporting a transition from a two mode beating of the center of mass and relative motion frequencies for few atoms to a single mode breathing for many particles, the breathing mode properties of fewbody BoseBose mixtures are not characterized so far. For this reason, we analyse the number of breathing frequencies and the kind of motion to which they correspond in dependence on the intra and intercomponent interaction for the binary mixture at hand.
This work is structured as follows. In Sec. II we introduce the Hamiltonian of the system. In Sec. III we perform a coordinate transformation to construct a fast converging correlated basis. Using exact diagonalization with respect to this basis we study in Sec. IV the lowlying energy spectrum for various interaction regimes. Sec. V is dedicated to the breathing dynamics within the linear response regime. An experimental realization is discussed in Sec. VI and we conclude the paper with a summary and an outlook in Sec. VII.
Ii Model
We consider a BoseBose mixture containing two components, which are labelled by , confined in a highly anisotropic harmonic trap. We assume the low temperature regime, where the interparticle interactions may be modelled via a contact potential, and strong transversal confinement allowing us to integrate out frozen degrees of freedom leading to a quasi1D model. Our focus lies on a mixture of particles, which have the same mass and trapping frequencies , in the transversal, longitudinal direction, respectively. This can be realized by choosing different hyperfine states of the same atomic species. By further rescaling the energy and length in units of and one arrives at the simplified Hamiltonian:
(1) 
with singlecomponent Hamiltonians and intercomponent coupling .
(2)  
(3) 
where with the 3D swave scattering length and are effective (offresonant) interaction strengths.
Iii Methodology: Exact diagonalization in a correlated basis
To obtain information on the lowenergy excitation spectrum we employ the wellestablished method of exact diagonalization. However, concerning the choice of basis, instead of taking bosonic number states w.r.t. HO eigenstates as in e.g. GarcíaMarch et al. (2014a), we pursue a different approach by using a correlated atompair basis. Thereby we can efficiently cover regimes of very strong intra and intercomponent interaction strengths and achieve convergence for relatively small basis sizes of about 700.
Actually, the idea of choosing optimized basis sets to speed up the convergence with respect to the size of basis functions can be also seen in the context of the potentialoptimized discrete variable representation (PODVR) Echave and Clary (1992). Here, one employs eigenstates of conveniently constructed onedimensional reference Hamiltonians in order to incorporate more information on the actual Hamiltonian into the basis compared to the standard DVR technique Light and Carrington Jr (2000); Beck et al. (2000). Another approach, stemming from nuclear physics, uses an effective twobody interaction potential instead of an optimized basis for solving ultracold manybody problems Christensson et al. (2009); Rotureau (2013); Lindgren et al. (2014).
In order to construct a tailored basis, which already incorporates intracomponent correlations, let us neglect for a moment the intercomponent coupling , which leaves us with two independent bosonic components, each consisting of two particles. The problem of two particles in a harmonic trap can be solved analytically Busch et al. (1998) and boils down to a coordinate transformation and solving a Weber differential equation^{1}^{1}1 with and under a delta potential constraint. Each eigenstate of this bosonic twoparticle problem turns out to be a tensor product of a HO eigenfunction of the center of mass (CM) coordinate and a normalized as well as symmetrized Parabolic Cylinder Function^{2}^{2}2 with the Tricomi’s hypergeometric function (PCF) of the relative coordinate with being a real valued quantum number depending on the interaction strength and excitation level , which is obtained by solving a transcendental equation derived from the deltatype constraint:
(4) 
Coming back to the binary mixture problem we apply a coordinate transformation to the relative frame defined by:

total CM coordinate
, 
relative CM coordinate
, 
relative coordinate for each component
.
The Hamiltonian attains a new structure in this frame. Firstly, the total CM is separated and is simply a HO of mass .
(5) 
with the spectrum and . The remainder of the Hamiltonian can be decomposed as , where can be solved analytically and couples the eigenstates of . is HO of mass and lead to the above mentioned Weber differential equations with deltatype constraint.
(6)  
(7)  
(8) 
Now to diagonalize we choose as basis the eigenvectors of , which we label as with . The energy of a corresponding basis function and its spatial representation are given by:
(9)  
(10) 
where are HO eigenstates and . All are of even parity because of the bosonic nature of particles of each component.
The main challenge now is the calculation of the matrix elements of , which are complicated 2D integrals at first sight and need to be tackled numerically. In the Appendix we provide a circumvention of this problem via the Schmidt decomposition Nielsen and Chuang (2000), allowing us to replace one 2D integral by multiple 1D integrals, which results in faster computation times. In quantum chemistry, the algorithm for achieving such a representation is known as POTFIT Jäckle and Meyer (1996). We will point out several symmetries to avoid the calculation of redundant terms and outline (dis)advantages of the whole method.
To summarize, the coordinate transformation to the chosen relative frame (i) decouples the CM motion and (ii) naturally guides us to employ the analytically known eigenstates of as the basis states in order to incorporate intracomponent correlations into our basis.
Iv Stationary Properties
By means of the correlated basis introduced above and an efficient strategy for calculating the Hamiltonian matrix to be diagonalized (see Appendix A), we can easily obtain the static properties of our system for a huge variety of different intra and intercomponent interaction strengths. Before going into the details, we need to address the symmetries of in the relative frame.
Symmetry analysis
First of all, commutes with the individual parity operators of the relative frame coordinates, . The eigenvectors of are restricted to even parity because of the bosonic character of our components. Due to the decoupling of it is sufficient to consider only the ground state of the total CM motion, which is of even parity, in the following. Then, the parity of the degree of freedom completely determines the total parity of the eigenstates. Another symmetry arises, if one chooses equal intracomponent interaction strengths . Under these circumstances the Hamiltonian is invariant under exchange, which we define as the transformation. It should be noted, that translating all these transformations to the laboratory frame leads to certain proper and improper rotations of the fourdimensional coordinate space.
Energy spectra
In Fig. 1 we show the total energy spectrum as a function of for various fixed values of and . The total CM is assumed to be in its ground state. Fig. 1(a) depicts the noninteracting intracomponent scenario . For the Hamiltonian represents two uncoupled noninteracting bosonic species and we will label this regime as BECBEC following the nomenclature of GarcíaMarch et al. (2014a). The eigenenergies are integers with equal spacings of , which is 1 in our units. In this limit the PCFs are even HO eigenstates of mass . The eigenenergies are thus .
For , the th eigenenergy corresponding to only even (odd) parity eigenstates is fold degenerate with . Already a small intercomponent coupling lifts all these degeneracies such that branches of eigenenergies arise. In the following, we label these resulting branches as the th even or odd parity branch, respectively. Note that this grouping of the energy levels into branches will be used in the following for all values of and in particular for the analysis of the breathing dynamics in section V. As we further increase the intercomponent coupling strength, we observe that states corresponding to branches of opposite parity incidentally cross, as they are of different symmetry and consequently not coupled by the perturbation.
For very strong values, i.e. in the composite fermionization (CF) limit Zöllner et al. (2008); Hao and Chen (2009), we observe a restoration of degeneracies, but in a different manner, namely the lowest states merge pairwise forming a twofold degeneracy. In this regime the two components spatially separate for the ground state, where one component locates on the left side of the trap, while the other is pushed to the right side due to the strong intercomponent repulsion. The twofold degeneracy of the ground state reflects actually the two possible configurations: A left B right and A right B left. This behaviour can be observed in the relative frame densities, discussed later in this section. Another striking peculiarity for are noninteger eigenvalues and unequal energy spacings.
A very similar analysis concerning this specific choice of interactions ( and arbitrary ) was performed in Dehkharghani et al. (2015), where an effective interaction approach was employed to greatly improve the convergence properties of exact diagonalization in order to access properties of a BoseBose mixture up to particles. However, the analysis only covered a single line of the parameter space. Here, we extend the ground state analysis of GarcíaMarch et al. (2014a) and study the lowlying excitations for arbitrary interactions.
In Fig. 1(b) we show the impact of moderate but symmetric intracomponent interactions of strength . Already in the uncoupled regime () we observe fewer degeneracies compared to the case. Nevertheless, we group the eigenstates into branches of even / odd parity also for finite by continuously following the eigenenergies to the limit. The reason for the reduced degeneracies is that PCFs are not HO eigenstates any more, while PCFs of both components are still the same. The energy is . To roughly estimate the energetic ordering it is sufficient to know that the realvalued quantum number fulfils for the following relations:
meaning that a single excitation of the relative motion is energetically below a double excitation of the degree of freedom. E.g. the first even parity branch in the uncoupled noninteracting regime (BECBEC in Fig. 1(a)) contains three degenerate states: , and (eq. (9)). By choosing finite values acquires a higher energy than and leading to reduced degeneracies in the spectrum. Another striking feature is the appearance of additional crossings between states of the same parity due to the symmetry. States, which possess different quantum numbers concerning the transformation ( or ), are allowed to cross as they ’randomly’ do throughout the variation. Of course, such crossings are also present in the previous noninteracting case, it being also componentsymmetric. An avoided crossing between a state of the first even parity branch and a state of the second even parity branch is worth mentioning, which is present for all values of (see the exemplary arrow in Fig. 1 (b) or (d)). States of the same symmetry obviously do not cross according to the Wignervon Neumann noncrossing rule von Neumann and Wigner (1929).
In Fig. 1(c) we asymmetrically increase the intracomponent interactions as compared to the noninteracting case, namely to and . For the uncoupled scenario () all the degeneracies are lifted, because now PCFs of the and the components are different. The energy is . The energetic state ordering is far from obvious, which becomes apparent upon closer inspection of the function. E.g. consider again the first even parity branch. Its lowest energy state is a single excitation of , followed by a single excitation of . The highest energy of this branch corresponds to a double excitation. The ordering pattern for higher order branches is even more complicated. For intermediate values of we observe that crossings from the previous scenario (with ) between states of the same parity are replaced by avoided crossings because of the broken symmetry. The strong coupling regime displays less degeneracies as compared to the componentsymmetric cases of Fig. 1(a) and (b) with the twofold ground state degeneracy remaining untouched.
In Fig. 1(d) we choose very strong intracomponent interaction strengths . When the coupling is absent, we have two hardcore bosons in each component. The system can thus be mapped to a twocomponent mixture of noninteracting fermions Girardeau (1960); Yukalov and Girardeau (2005) and will be referred to as TGTG limit. The PCFs become near degenerate with odd HO eigenstates, which again leads to integervalued eigenenergies with equal spacings and the same degree of (near)degeneracies as in the noninteracting case (Fig. 1(a)). The limit of strong intercomponent coupling displays a completely different structure of the spectrum. The socalled full fermionization (FF) Girardeau and Minguzzi (2007) phase can be mapped to a noninteracting ensemble of four fermions with the ground state energy . However, in contrast to the singlecomponent case of four bosons, we need to take into account that the components are distinguishable. The degeneracy of the ground state is expected to be fold and corresponds to the different possibilities of ordering the laboratory frame coordinates while keeping in mind the indistinguishability of particles of each component. A profound study of these ground states can be done by employing a snippet basis Deuretzbacher et al. (2008) of distinguishable particles, which interact with each other by an infinite repulsive deltainteraction. According to eq. (2) of Ref. Deuretzbacher et al. (2008), the snippet basis for distinguishable particles is defined as:
(11) 
where represents the ground state of noninteracting fermions and is a permutation of particle coordinates, which defines the sector where has support. Every snippet state is one of possible ground states in the hardcore limit.
To adjust this basis to our case one can follow a procedure described for a BoseFermi mixture in Fang et al. (2011). The only modification is to replace the antisymmetric exchange symmetry of the fermionic species by a symmetric one. With respect to the spatial projection we will combine the sectors of the distinguishable case into ground states of our BoseBose mixture. Each ground state consists of different sectors, reflecting the indistinguishability within each component, meaning that by exchanging the spatial order of identical particles we switch between these sectors. By exchanging the spatial order of distinguishable particles we switch between the different ground states. Thus, the permutations can be decomposed in transpositions which exchange identical particles and transpositions which exchange distinguishable particles, leading to the following ground state configurations:
where the notation describes the exchange of particles , is the identity permutation and for two permutations , , is defined as . The subindex of the wavefunction relates to the particle arrangements in the onebody density distribution.
The last case we discuss is the highly asymmetric case and in Fig. 1(e). For (BECTG) one expects, based on the previous considerations, integer eigenvalues and thus equal spacings as for the cases and depicted in Fig. 1(a) and (d), because the PCF of the component is an even HO eigenstate and the PCF of the component is degenerate with an odd HO eigenstate. Very peculiar is the strong coupling case, where we observe a nondegenerate ground state, the socalled phase separation (PS) phase GarcíaMarch et al. (2014a), where the component occupies the centre of the harmonic trap, while the component, in order to reduce its intracomponent interaction energy, forms a shell around the component.
Relativeframe densities
Let us now inspect the relativeframe probability densities instead of the usually studied onebody densities of the laboratory frame as e.g. in GarcíaMarch et al. (2014a). We will see that these quantities can be used to identify regions of most probable relative distances and provide a more detailed picture of particle arrangements than their laboratory frame counterparts. Moreover, in the quench dynamics study, the subject of the next section, an occupation of a certain eigenstate of will lead to the breathing oscillation of only one relativeframe density, making it possible to connect different breathing modes to specific relative motions within the system, at least for the weakly coupled case .
We define these quantities as follows:
(12) 
where is the th eigenstate of and we trace out all the degrees of freedom of the relative frame except for one.
Let’s compare our results concerning the ground state densities for some limiting cases to the ones obtained in GarcíaMarch et al. (2014a). In Fig. 2 we show the densities for all the degrees of freedom except for , which trivially obeys a Gaussian distribution. In the BECBEC case all the densities are characterized by a Gaussian density profile, since the Hamiltonian consists of completely decoupled HOs for each degree of freedom. The TGTG limit differs from the BECBEC case in the distributions featuring two maxima and a minimum in between, a result of strong repulsion within each component. This behaviour reflects actually the already known results of the analytical twoparticle solution Busch et al. (1998), where develops a minimum in the centre of the density distribution for finite whose value tends to zero as .
The CF phase is in some sense a complete counterpart to the TGTG case. Now features two maxima and a minimum in between, a result of A and B strongly repelling each other. This feature is blurred in the onebody density distributions of the laboratory frame and one needs to additionally consider the twobody density function to verify this behaviour Zöllner et al. (2008). The density distribution of is more compressed compared to the BECBEC case due to the tighter confinement induced by the other component.
Finally, the PS phase corresponds to a coreshell structure, where and show a more pronounced peak, while obeys a bimodal distribution with two density peaks being much further apart than both in the CF and in the TGTG case. This can be understood in the following way: firstly, the fact that A locates in the trap centre and not the other way around is because B needs to minimize its repulsive intracomponent interaction energy by separating its particles. Secondly, the need to minimize the repulsive intercomponent energy pushes the B particles even further along the harmonic trap at the cost of increased potential energy until these two energies balance themselves out. The two A particles are compressed to closer distances as compared to the BECBEC case because of a tighter trap induced by B, while at the same time A modifies the HO potential to a double well for B. This results in stronger localization of particles, which leads to a more pronounced peak in the distribution.
V Breathing Dynamics
The spectral properties discussed above can be probed by slightly quenching a system parameter such that the lowest lying collective modes are excited. Here, we focus on a slight quench of the trapping frequency in order to excite the breathing or monopole modes being characterized by a periodic expansion and compression of the atomic density. While in the singlecomponent case two lowest lying breathing modes of in general distinct frequencies exist, being associated with a motion of the CM and the relative coordinates Bauch et al. (2009); Schmitz et al. (2013), respectively, the number of breathing modes, their frequencies and the associated “normal coordinates” are so far unknown for the more complex case of a binary fewbody mixture and shall be the subject of this section.
Experimentally, breathing oscillations can be studied by measuring the width of species density distribution where we have omitted the subtraction of the mean value squared, which vanishes due to the parity symmetry. From a theoretical point of view, it is fruitful to define a breathing observable as , whose expectation value is essentially the sum of the widths of the A and the B component.
To study the breathing dynamics we will perform a slight and componentsymmetric quench of the HO trapping frequency, where our HO units will be given with respect to the postquench system. The initial state for the timepropagation is the ground state of some prequench Hamiltonian with frequency . Then a sudden quench is performed to . The time evolution of this state is thus described as follows:
(13) 
where and is the overlap between the initial state and the th eigenstate of the postquench Hamiltonian . Since both the pre and postquench Hamiltonian are timereversal symmetric, we assume their eigenstates and thereby also the overlap coefficients to be realvalued without loss of generality. A small quench ensures that and only the lowest excited states are of relevance. Symmetry considerations further reduce the number of allowed contributions. E.g. states of odd parity or odd parity have zero overlap with , because the initial state is of even  and parity and the quench does not affect any of the symmetries discussed in the previous section. Similarly, in the componentsymmetric case , states, which are antisymmetric w.r.t. the operation, have no overlap with the prequench ground state being symmetric under .
For the weakly coupled regime, the relativeframe coordinates turn out to be extremely helpful for characterizing the participating breathing modes. Therefore, we study in particular the reduced densities of the relativeframe coordinates. Employing the expansion in postquench eigenstates from eq. (13), their timeevolution may be approximated within the linearresponse regime as
(14) 
neglecting terms of the order for . So can be decomposed into the stationary background and timedependent modulations of the form , further called transition densities, with oscillation frequency .
In the following, we regard the excitations of the first even parity branch as the lowest monopole modes and show that each monopole mode is directly connected to the breathing modulation of a single relativeframe density, if the two components are but weakly coupled. This behaviour changes for increasing , where each coordinate begins to exhibit an oscillation with more than one frequency. By inspecting the modulations of the variances of each relativecoordinate and taking the excitation amplitudes into account, we show that four (three) breathing modes are excited for () in the weakly coupled regime, while only two breathing modes are of relevance in the strongly coupled regime. First, we inspect the componentasymmetric case of Fig. 1(c) in detail to illustrate some peculiarities of involved breathing modes, since it contains the most relevant features. Thereafter, we unravel differences to the componentsymmetric case of Fig. 1(b).
v.1 Componentasymmetric case
Because of the low amplitude quenching protocol, we will excite four breathing modes simultaneously in the componentasymmetric case (, ). Three of them stem from the first even parity branch of Fig. 1(c). Remember, however, that the total CM was assumed to be in the ground state to keep the spectrum discernible. One obtains the full spectral picture by including all CM excitations, meaning duplicating and upshifting depicted energy curves by with . This reveals a forth mode, namely a double total CM excitation. It features the same parity symmetries and is energetically of the same order as the states from the first even parity branch ensuring a considerable overlap with the initial state. The total CM trivially oscillates with the constant frequency independent of any interactions it being a decoupled degree of freedom with the singleparticle HO Hamiltonian (eq. (5)) Bauch et al. (2010); Schmitz et al. (2013). The other three modes, which are excited, are known analytically, when there is no coupling between the components, and we label the corresponding mode frequencies as:

,

,

,

,
where we have prepended the CM eigenstate for a complete characterization of the involved states. States of higher order even parity branches as well as higher excitations of the CM coordinate are negligible due to small overlaps with the initial state.
In the uncoupled regime , one can show analytically that each relativecoordinate density oscillates with a single frequency, each corresponding to exactly one eigenstate of the first even parity branch (see Fig. 3 (a)(c)). E.g. for , the only transition density which survives taking the partial trace is the one corresponding to , while the contributions from the remaining excited states vanish. This leads to the breathing motion in the coordinate with a single frequency . Analogously one can show that solely induces density modulation in with the frequency , while oscillates with exclusively due to . Thereby, the relativeframe coordinates render “normal coordinates” in the uncoupled regime, which is also a valid picture for extremely weak couplings.
By introducing a larger coupling between the components one observes that each relative frame density, except for , begins to oscillate with up to three frequencies simultaneously. So all the modes begin to contribute to the density modulation of each relative coordinate. However, there are some peculiarities we observe, for the visualization of which the densities are not well suited any more. Instead, we will transform the breathing observable to the relative frame and consider the expectation values of individual terms it decomposes into:
(15) 
The expectation value of each observable with respect to the timeevolved state is directly related to the respective relativeframe density:
(16) 
Inserting the timeevolution of the relative frame density from eq. (14) one finds that the observables decompose into a stationary value and a timedependent modulation as well. In particular, we are interested in the amplitudes of modulations, when the intercomponent coupling is varied, since they determine how many frequencies are of essential relevance for the considered motion. The amplitude of the th mode is essentially composed of the overlap and of the transition element:
(17) 
In order to evaluate the overlaps with in terms of only the postquench Hamiltonian eigenstates, we perform a Taylor approximation with respect to the weak quench strength , namely evaluated at , and arrive at
(18) 
Applying the (offdiagonal) HellmannFeynman theorem, one obtains:
(19) 
The overlaps are hence connected to the transition elements of each relativeframe breathing observable (eq. (15)), weighted with the inverse of the mode frequency, which leads to a damping of contributions from higher order branches. This relation enables us to calculate the amplitude , with which the th mode contributes to the oscillation of the observable :
(20) 
which may be interpreted as the susceptibility of the observable for the excitation of the state .
In Fig. 4 (a) we show the values of possible breathing mode frequencies, obtained from the spectrum of Fig. 1 (c). does not depend on any interactions . In contrast to this, is degenerate with for , and when increasing , decreases to a minimum first and then increases with the tendency to asymptotically reach again. This behaviour strongly resembles the dependence of the relativecoordinate breathingmode frequency in the singlecomponent case Schmitz et al. (2013). and have qualitatively akin curve shapes, varying much stronger with . In particular, we note that these frequencies reach values below the frequency of the CM dipole mode being equal to unity.
The breathing mode frequencies discussed above are labelled according to the peculiarity of the uncoupled regime, where each eigenstate from the first even parity branch leads to a breathing motion of some specific relativeframe coordinate. Indeed, if we look at the amplitudes in Fig. 4 (b)(d) in the decoupled regime (), we recognize that the amplitude for the coordinate is nonzero only for one mode, namely the one with which oscillates in Fig. 3 (a)(c). When we increase the intercomponent coupling, the eigenstates cease to be simple product states in the relative coordinate frame resulting in contamination of each density modulation with the frequencies from the other modes as well, which leads to a threemode oscillation. Nevertheless, we label the frequencies corresponding to the uncoupled case and follow the states continuously throughout the variation.
Another peculiarity worth noting arises in the strongly coupled regime: oscillations become strongly suppressed for all the observables making and the main contributors to the density modulations. In Fig. 3 (d)(f) we show that all the relativecoordinate densities oscillate with the same frequency . We highlight that in the CF regime each density peak of the bimodal distribution does not only breath periodically but also its maximum height position performs dipolemode like oscillations, see Fig. 3 (d).
v.2 Componentsymmetric case
Now we compare the above results with the componentsymmetric case of Fig. 1 (b), where . Fig. 5 (a) depicts the possible breathingmode frequencies. Here, two main differences arise: First, the curve features two minima due to two avoided crossings with an eigenstate of the second even parity branch. Second, the other two breathing mode frequencies of the relative coordinates are degenerate for , then separate with increasing and approach one another asymptotically. Instead of using the labels and as in the componentasymmetric case, we label these two modes with and , since the corresponding excited eigenstates are symmetric and antisymmetric w.r.t. , respectively (see Fig. 1 (b)). The latter frequency, however, does not give any contribution to the breathing dynamics, it being symmetry excluded.
The motion of the and coordinates is identical due to the imposed componentsymmetry and in the uncoupled regime () they both oscillate solely with the frequency (see Fig. 5 (c) and (d)). Similarly to the componentasymmetric case, we see that when increasing the () mode contributes also to the observable ().
However, in the intermediate interaction regime we observe a strong suppression of the mode contribution for all observables, which is in stark contrast to the componentasymmetric case. Instead, a state of the second even parity branch gains relevance. These two eigenstates actually participate in the avoided crossing exemplary indicated by an arrow in the spectrum (see Fig. 1 (b)) and thereby exchange their character. By further increasing we observe another exchange of roles, which is attributed to the presence of a second avoided crossing (see Fig. 5 (a)), such that the strong intercomponent coupling regime shows again the absolute dominance of the mode over the other lowest breathing modes besides .
Vi Experimental realization
The fewbody BoseBose mixture studied here should be observable with existing cold atom techniques. Quantum gas microscopes allow the detection of single particles in a well controlled manybody or fewbody system Bakr et al. (2010); Sherson et al. (2010) and recent progress also allows for spinresolved imaging in 1D systems using an expansion in the perpendicular direction Preiss et al. (2015); Boll et al. (2016). In these setups, the singleparticle sensitivity relies on pinning the atoms in a deep lattice during imaging and experiments have so far focused on lattice systems. However, bulk systems might be imaged with high spatial resolution by freezing the atomic positions in a lattice before imaging. For fast freezing, this would allow a timeresolved measurement of the breathing dynamics. Moreover, spin order was recently observed in very small fermionic bulk systems via spinselective spilling to one side of the system Murmann et al. (2015).
Deterministic preparation of very small samples was demonstrated for fermions via trap spilling Serwane et al. (2011) and for bosons by cutting out a subsystem of a Mott insulator Islam et al. (2015). The tight transverse confinement for a 1D system can be obtained from a 2D optical lattice, while the axial confinement would come from an additional optical potential, which can be separately controlled to initialize the breathing mode dynamics.
Choosing two hyperfine states of the same atomic species ensures the same mass of the two bosonic species. Possible choices include , or . While the former have usable FRs to tune the interaction strengths, the latter allows selective tuning via CIR in a spindependent transversal confinement as can be realized for heavier elements. Note that the longitudinal confinement needs to be spinindependent in order to ensure the same longitudinal trap frequencies and trap centres assumed in the calculations. The intercomponent interaction strength can be tuned via a transverse spatial separation as obtained e.g. from a magnetic field gradient Schachenmayer et al. (2015). In the case of and , the intercomponent background scattering lengths are negative Chin et al. (2010) leading to negative , but although not reported, intercomponent FRs might exist. Alternatively, might be tuned via a CIR, which selectively changes at magnetic fields, where the intracomponent scattering lengths are very different.
In the following we give concrete numbers for a choice of . The density distribution has structures on the scale of the HO unit (Fig. 2). Choosing the trap parameters as in Ref. Serwane et al. (2011) as kHz and kHz, yields , i.e. larger than a typical optical resolution of . At the same time, temperatures much lower than nK are state of the art. The breathing dynamics will occur on a time scale of several , which is easily experimentally accessible. Choosing a smaller would make the imaging easier, but would impose stricter requirements on the temperature.
For the observation of the breathing mode dynamics, one would record the positions of all four particles in each experimental image and obtain the widths , , by averaging the occurring relative coordinates over many singleshots after a fixed hold time .
Vii Discussion and Outlook
In this work we have explored a fewbody problem of a BoseBose mixture with two atoms in each component confined in a quasi1D HO trapping potential by exact diagonalization. By applying a coordinate transformation to a suitable frame we have constructed a rapidly converging basis consisting of HO and PCF eigenfunctions. The latter stem from the analytical solution of the relative part of the twoatom problem Busch et al. (1998) and include the information about the intraspecies correlations, which renders our basis superior to the common approach of using HO eigenstates as basis states.
We have then explored the behaviour of the lowlying energy spectrum as a function of the interspecies coupling for various fixed values of intraspecies interaction strengths. Hereby we have covered the strongly coupled limiting cases of Composite Fermionization, Full Fermionization and Phase Separation, studied also intermediate symmetric and asymmetric values of and and related the ground state relativeframe densities of some limiting cases to the known laboratory frame results GarcíaMarch et al. (2014a). We have discussed the evolution of degeneracies and explained appearing (avoided) crossings in terms of the symmetries of the Hamiltonian, which become directly manifest in the chosen relativecoordinate frame.
Finally, the obtained results were used to study the dynamics of the system under a slight componentsymmetric quench of the trapping potential. We have derived expressions for the time evolution of the relativeframe densities within the linear response regime and observed that in the uncoupled regime () the density of each relative frame coordinate performs breathing oscillations with a single frequency corresponding to a specific excited state of the first even parity branch of the spectrum. The total CM coordinate performs breathing oscillations with the frequency (HO units). For asymmetric choices of values, three additional monopole modes participate in the dynamics, each of them corresponding to the motion of a particular relative coordinate: for the relative coordinate of the A component, for the relative coordinate of the B component and for the relative distance of the CMs of both components. In contrast to this, the symmetric case leads to only two additional modes because of a symmetryinduced selection rule: for the relative coordinates of both components and for the relative distance of the CMs of both components.
For not too strong intercomponent coupling, each relative coordinate exhibits multimode oscillations and we have explored their relevance for the density modulations by analysing the behaviour of suitably chosen observables as one gradually increases the coupling between the components for symmetric and asymmetric choices of intracomponent interactions strengths. Thereby, we have found that for strong couplings, where Composite Fermionization takes place, the () modes become highly suppressed, leaving only two monopole modes in this regime: and . We have observed the same effect for the case of Phase Separation (results not shown). Interestingly, the dependence of on strongly resembles the behaviour of the relativecoordinate breathing frequency in the singlecomponent case Schmitz et al. (2013). All in all, we have obtained 2 to 4 monopole modes for the quench dynamics depending on the strength of the intercomponent coupling and the symmetry of the intraspecies interactions, which is in strong contrast to the singlecomponent case Schmitz et al. (2013) as well as to the MF results, where two lowlying breathing modes can be obtained, namely an inphase (outofphase) mode for a componentsymmetric (componentasymmetric) quench Morise and Wadati (2000). Finally, we have argued that the experimental preparation of the considered fewbody mixture and measurement of the predicted effects are in reach by means of stateoftheart techniques.
This work serves as a useful analysis tool for future fewbody experiments. Measurements of the monopole modes can be mapped to the effective interactions within the system such that precise measurements of the scattering lengths or external magnetic fields can be performed. The numerical method used here can be applied to BoseFermi and FermiFermi mixtures with two particles in each component simplifying the numerics, because the PCFs have to be replaced by odd HO eigenstates, if a bosonic component is switched to a fermionic one, which significantly accelerates the calculation of integrals. Further, it would be interesting to see how the frequencies and the amplitudes of the monopole modes vary for an increasing number of particles. Exploring the spectrum for negative values of interaction parameters is also a promising direction of future research.
Acknowledgements.
The authors acknowledge fruitful discussions with H.D. Meyer, J. Schurer, K. Keiler and J. Chen. C.W. and P.S. gratefully acknowledge funding by the Deutsche Forschungsgemeinschaft in the framework of the SFB 925 ”Light induced dynamics and control of correlated quantum systems”. S.K. and P.S. gratefully acknowledge support for this work by the excellence cluster ”The Hamburg Centre for Ultrafast ImagingStructure, Dynamics and Control of Matter at the Atomic Scale” of the Deutsche Forschungsgemeinschaft.Appendix A
In the following, we discuss how to efficiently calculate matrix elements of the coupling operator from eq. (8) with respect to the basis (9). Because of the already mentioned even parity of one can make simple substitutions of the form to show that each delta in the sum of gives the same contribution, such that after performing an integral over one obtains:
(21) 
At this point it is important to notice that the integral vanishes for odd because of the parity symmetries of and , which can be seen by transforming to relative and centerofmass coordinates , . In the following, we assume that the quantum numbers for the coordinate are restricted to and for the coordinates to . Now our computation strategy consists of three steps:
First, we circumvent evaluating the 2D integral from eq. (21) by viewing the product of the two HO eigenstates as a pure, in general not normalized state depending on the two coordinates and and applying the Schmidt decomposition Nielsen and Chuang (2000) or, equivalently, the socalled POTFIT algorithm Jäckle and Meyer (1996).
Here^{3}^{3}3Note that in contrast to the usual convention we do not require the coefficients to be semipositive without loss of generality., coincides with the th eigenvalue of the reduced onebody density matrix corresponding to the degreeoffreedom and denotes the corresponding eigenvector, which can be shown to feature a definite parity symmetry. We remark that (i) we may choose because of without loss of generality and that (ii) the decomposition of eq. (A) becomes exact for . Having ordered the coefficients in decreasing sequence w.r.t. to their modulus, we choose such that only terms with are taken into account, which results in an accurate approximation to . We perform this decomposition for all the relevant HO quantum numbers , meaning with even. This procedure is independent of any interactions and needs to be executed only once.
By inserting eq. (A) into eq. (21) we obtain the following expression:
(23) 
As one can see, the 2D integral is replaced by a sum of products of 1D integrals. In order to greatly overcome the numerical effort for computing a 2D integral should be preferably a small number (see below).
The second step consists in the calculation of 1D integrals and here we provide an efficient strategy to circumvent redundant computations. Consider the integral:
(24) 
with . The PCFs are of even parity and thus only even actually contribute allowing to reduce the number of expansion terms in eq. (23) to . Both PCFs and depend on the same interaction strength meaning that the above integral does not distinguish between the two subsystems such that the index can be dropped for the moment. Now we fix the PCFs by specifying the strength of intracomponent interaction and further we fix the HO quantum numbers , which determine the functions , as well as PCF quantum numbers . We loop over all and save the integral values, labelled as . This procedure is performed for a set of multiple values of we are interested in and for all the relevant quantum number configurations with even and .
In the last step we calculate the matrix elements from eq. (23). For HO quantum numbers we extract all the expansion coefficients , obtained in the first step, and for the chosen interactions parameters and PCF quantum numbers we pick the appropriate integral values corresponding to and . The advantage of this procedure is that not only symmetric choices of intracomponent interaction strengths are accessible, but also an arbitrary asymmetric combination . Additionally, the proposed scheme can be easily parallelized. However, adding new values to the set is in general very timeconsuming as one needs to calculate a bunch of 1D integrals for all the relevant quantum number configurations.
Now let’s analyse quantitatively the speedup obtained by our algorithm in contrast to the straightforward evaluation of 2D integrals. Since the energy spacing of the PCF modes is approximately twice the energy spacing of the HO modes corresponding to the motion, we assume with an odd to keep the number of even and odd parity basis states the same. The number of 1D integrals one needs to compute for each in order to construct the matrix is approximately , where is an average number of terms in (23), as the criterion requires more terms for larger values of . The number of 2D integrals amounts to . For checking the convergence (see below), we have chosen and , i.e. basis states. The number of expansion terms varies in the interval resulting in . Thus, we need to either evaluate 1D integrals or 2D integrals. Not only is the number of 1D integrals smaller, the computation of one 2D integral takes also significantly longer than of one 1D integral, especially for higher quantum numbers. Moreover, in order to build the Hamiltonian matrix for all , the 2D integrals (21) would have to be evaluated for the distinct combinations , where denotes the cardinality of , while the 1D integrals (24) must be calculated only for all , i.e. distinct values, which renders this approach much more efficient.
For the spectra shown in section IV we have chosen and , i.e. basis states. Since we know, that basis states of different symmetry do not couple, we can split the matrix into subspaces of even and odd