Spectral properties and breathing dynamics of a few-body Bose-Bose mixture in a 1D harmonic trap

# Spectral properties and breathing dynamics of a few-body Bose-Bose mixture in a 1D harmonic trap

Maxim Pyzh Zentrum für Optische Quantentechnologien, Universität Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany    Sven Krönke Zentrum für Optische Quantentechnologien, Universität Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany The Hamburg Centre for Ultrafast Imaging, Universität Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany    Christof Weitenberg Institut für Laserphysik, Universität Hamburg, 22761 Hamburg, Germany    Peter Schmelcher Zentrum für Optische Quantentechnologien, Universität Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany The Hamburg Centre for Ultrafast Imaging, Universität Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany
July 16, 2019
###### Abstract

We investigate a few-body mixture of two bosonic components, each consisting of two particles confined in a quasi one-dimensional harmonic trap. By means of exact diagonalization with a correlated basis approach we obtain the low-energy spectrum and eigenstates for the whole range of repulsive intra- and inter-component interaction strengths. We analyse the eigenvalues as a function of the inter-component coupling, covering hereby all the limiting regimes, and characterize the behaviour in-between these regimes by exploiting the symmetries of the Hamiltonian. Provided with this knowledge we study the breathing dynamics in the linear-response 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 inter-component coupling first leads to multi-mode oscillations in each relative coordinate, which turn into single-mode oscillations of the same frequency in the composite-fermionization regime.

## I Introduction

The physics of ultra-cold atoms has gained a great boost of interest since the first experimental realization of an atomic Bose-Einstein 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 lower-dimensional 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 ultra-cold gases the atoms are but weakly correlated and well described by a mean-field (MF) model, the well-known Gross-Pitaevskii equation (GPE), or in case of mixtures by coupled GPEs Ho and Shenoy (1996); Esry et al. (1997); Pu and Bigelow (1998a). Bose-Bose 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 inter-species interaction strengths, being experimentally tunable by e.g. Feshbach resonances (FR) Chin et al. (2010): the miscible (M), immiscible symmetry-broken (SB) or immiscible core-shell 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 low-lying excitations Egorov et al. (2013). In contrast to a single-species case the collective modes of mixtures exhibit new exciting phenomena: doublet splitting of the spectrum containing in-phase and out-of-phase oscillations, mode-softening for increasing inter-component 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 quasi-one-dimensional single-component 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 low-lying 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 quasi-1D set-up 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 Tonks-Girardeau limit Kinoshita et al. (2004); Paredes et al. (2004), where the bosons resemble a system of non-interacting 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 few-body systems. For instance, a profound investigation of the ground state phases of a few-body Bose-Bose 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 inter-component coupling. Only in the limit of infinite coupling the ground state becomes two-fold 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 time-independent problem of the simplest Bose-Bose mixture confined in a quasi-1D HO trap with two particles in each component, covering the whole parameter space of intra- and inter-species interactions, thereby complementing the analysis of some previous studies García-March 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 few-body Bose-Bose 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 inter-component 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 low-lying 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 Bose-Bose mixture containing two components, which are labelled by , confined in a highly anisotropic harmonic trap. We assume the low temperature regime, where the inter-particle interactions may be modelled via a contact potential, and strong transversal confinement allowing us to integrate out frozen degrees of freedom leading to a quasi-1D 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:

 H=∑σHσ+HAB, (1)

with single-component Hamiltonians and inter-component coupling .

 Hσ=2∑i=1(−12∂2∂x2σ,i+12x2σ,i)+gσδ(xσ,1−xσ,2), (2) HAB=gAB2∑i,j=1δ(xA,i−xB,j), (3)

where with the 3D s-wave scattering length and are effective (off-resonant) interaction strengths.

## Iii Methodology: Exact diagonalization in a correlated basis

To obtain information on the low-energy excitation spectrum we employ the well-established 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ía-March et al. (2014a), we pursue a different approach by using a correlated atom-pair basis. Thereby we can efficiently cover regimes of very strong intra- and inter-component 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 potential-optimized discrete variable representation (PO-DVR) Echave and Clary (1992). Here, one employs eigenstates of conveniently constructed one-dimensional 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 two-body interaction potential instead of an optimized basis for solving ultra-cold many-body problems Christensson et al. (2009); Rotureau (2013); Lindgren et al. (2014).

In order to construct a tailored basis, which already incorporates intra-component correlations, let us neglect for a moment the inter-component 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 equation111 with and under a delta potential constraint. Each eigenstate of this bosonic two-particle 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 Function222 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 delta-type constraint:

 g=−232Γ(1−μ2)Γ(−μ2). (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 .

 HRCM=−18∂2∂R2CM+2R2CM, (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 delta-type constraint.

 HRAB=−12∂2∂R2AB+12R2AB, (6) Hrσ=−∂2∂r2σ+14r2σ+gσδ(rσ), (7) H1=2∑i,j=1δ(RAB+(−1)irA2+(−1)jrB2). (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:

 ⟨RAB,rA,rB|k,l,m⟩=ΦABk(RAB)φAl(rA)φBm(rB), (9) E(0)k,l,m=k+μ(gA,l)+μ(gB,m)+32, (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 intra-component 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 inter-component 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 intra-component 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 four-dimensional 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 non-interacting intra-component scenario . For the Hamiltonian represents two uncoupled non-interacting bosonic species and we will label this regime as BEC-BEC following the nomenclature of García-March 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 inter-component 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 inter-component 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 two-fold 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 inter-component repulsion. The two-fold 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 non-integer 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 Bose-Bose 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ía-March et al. (2014a) and study the low-lying excitations for arbitrary interactions.

In Fig. 1(b) we show the impact of moderate but symmetric intra-component 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 real-valued 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 non-interacting regime (BEC-BEC 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 non-interacting case, it being also component-symmetric. 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 Wigner-von Neumann non-crossing rule von Neumann and Wigner (1929).

In Fig. 1(c) we asymmetrically increase the intra-component interactions as compared to the non-interacting 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 component-symmetric cases of Fig. 1(a) and (b) with the two-fold ground state degeneracy remaining untouched.

In Fig. 1(d) we choose very strong intra-component interaction strengths . When the -coupling is absent, we have two hard-core bosons in each component. The system can thus be mapped to a two-component mixture of non-interacting fermions Girardeau (1960); Yukalov and Girardeau (2005) and will be referred to as TG-TG limit. The PCFs become near degenerate with odd HO eigenstates, which again leads to integer-valued eigenenergies with equal spacings and the same degree of (near-)degeneracies as in the non-interacting case (Fig. 1(a)). The limit of strong inter-component coupling displays a completely different structure of the spectrum. The so-called full fermionization (FF) Girardeau and Minguzzi (2007) phase can be mapped to a non-interacting ensemble of four fermions with the ground state energy . However, in contrast to the single-component 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 delta-interaction. According to eq. (2) of Ref. Deuretzbacher et al. (2008), the snippet basis for distinguishable particles is defined as:

 ⟨x1,...,xN|Π⟩={√N!|ΨF0|ifxΠ(1)<...

where represents the ground state of non-interacting 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 hard-core limit.

To adjust this basis to our case one can follow a procedure described for a Bose-Fermi 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 Bose-Bose 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:

 ΨAABB ∝⟨→x|(e+(1,2))(e+(3,4)) (e)⟩, ΨBBAA ∝⟨→x|(e+(3,4))(e+(1,2)) (1,3)(2,4)⟩, ΨABAB ∝⟨→x|(e+(1,3))(e+(2,4)) (2,3)⟩, ΨBABA ∝⟨→x|(e+(2,4))(e+(1,3)) (1,4)⟩, ΨABBA ∝⟨→x|(e+(1,4))(e+(2,3)) (2,4)⟩, ΨBAAB ∝⟨→x|(e+(2,3))(e+(1,4)) (1,3)⟩,

where the notation describes the exchange of particles , is the identity permutation and for two permutations , , is defined as . The sub-index of the wave-function relates to the particle arrangements in the one-body density distribution.

The last case we discuss is the highly asymmetric case and in Fig. 1(e). For (BEC-TG) 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 non-degenerate ground state, the so-called phase separation (PS) phase García-March et al. (2014a), where the component occupies the centre of the harmonic trap, while the component, in order to reduce its intra-component interaction energy, forms a shell around the component.

#### Relative-frame densities

Let us now inspect the relative-frame probability densities instead of the usually studied one-body densities of the laboratory frame as e.g. in García-March 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 relative-frame 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:

 ρ(j)1(Yi)=∫∏p≠idYp |⟨→Y|Ej⟩|2, (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ía-March 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 BEC-BEC 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 TG-TG limit differs from the BEC-BEC 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 two-particle 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 counter-part to the TG-TG 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 one-body density distributions of the laboratory frame and one needs to additionally consider the two-body density function to verify this behaviour Zöllner et al. (2008). The density distribution of is more compressed compared to the BEC-BEC case due to the tighter confinement induced by the other component.

Finally, the PS phase corresponds to a core-shell 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 TG-TG 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 intra-component interaction energy by separating its particles. Secondly, the need to minimize the repulsive inter-component 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 BEC-BEC 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 single-component 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 few-body 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 component-symmetric quench of the HO trapping frequency, where our HO units will be given with respect to the post-quench system. The initial state for the time-propagation is the ground state of some pre-quench Hamiltonian with frequency . Then a sudden quench is performed to . The time evolution of this state is thus described as follows:

 |Ψ(t)⟩=e−iHt|E0⟩ω0≈n∑j=0cje−iEjt|Ej⟩, (13)

where and is the overlap between the initial state and the th eigenstate of the post-quench Hamiltonian . Since both the pre- and post-quench Hamiltonian are time-reversal symmetric, we assume their eigenstates and thereby also the overlap coefficients to be real-valued 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 component-symmetric case , states, which are antisymmetric w.r.t. the operation, have no overlap with the pre-quench ground state being symmetric under .

For the weakly coupled regime, the relative-frame coordinates turn out to be extremely helpful for characterizing the participating breathing modes. Therefore, we study in particular the reduced densities of the relative-frame coordinates. Employing the expansion in post-quench eigenstates from eq. (13), their time-evolution may be approximated within the linear-response regime as

 ρ1(Yi,t)≈c20ρ(0)1(Yi)+2n∑j=1c0cjρ(0,j)1(Yi)cos(Δjt), (14)

neglecting terms of the order for . So can be decomposed into the stationary background and time-dependent 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 relative-frame 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 relative-coordinate 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 component-asymmetric 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 component-symmetric case of Fig. 1(b).

### v.1 Component-asymmetric case

Because of the low amplitude quenching protocol, we will excite four breathing modes simultaneously in the component-asymmetric 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 up-shifting 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 single-particle 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:

1. ,

2. ,

3. ,

4. ,

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 relative-coordinate 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 relative-frame 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:

 ∑σ,ix2σ,i=4R2CM+R2AB+12r2A+12r2B. (15)

The expectation value of each observable with respect to the time-evolved state is directly related to the respective relative-frame density:

 ⟨Ψ(t)|Y2i|Ψ(t)⟩=∫dYi Y2i ρ1(Yi,t). (16)

Inserting the time-evolution of the relative frame density from eq. (14) one finds that the observables decompose into a stationary value and a time-dependent modulation as well. In particular, we are interested in the amplitudes of modulations, when the inter-component 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:

 ⟨Ej|Y2i|E0⟩=∫dYi Y2i ρ(0,j)1(Yi). (17)

In order to evaluate the overlaps with in terms of only the post-quench Hamiltonian eigenstates, we perform a Taylor approximation with respect to the weak quench strength , namely evaluated at , and arrive at

 cj≈−⟨Ej|ddω|E0⟩δω. (18)

Applying the (off-diagonal) Hellmann-Feynman theorem, one obtains:

 ⟨Ej|ddω|E0⟩=−ω⟨Ej|∑i,σx2i,σ|E0⟩Ej−E0. (19)

The overlaps are hence connected to the transition elements of each relative-frame 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 :

 Aj(Yi)=∣∣∣cjδω⟨Ej|Y2i|E0⟩∣∣∣, (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 relative-coordinate breathing-mode frequency in the single-component 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 relative-frame 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 non-zero only for one mode, namely the one with which oscillates in Fig. 3 (a)-(c). When we increase the inter-component 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 three-mode 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 relative-coordinate 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 dipole-mode like oscillations, see Fig. 3 (d).

### v.2 Component-symmetric case

Now we compare the above results with the component-symmetric case of Fig. 1 (b), where . Fig. 5 (a) depicts the possible breathing-mode 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 component-asymmetric 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 component-symmetry and in the uncoupled regime () they both oscillate solely with the frequency (see Fig. 5 (c) and (d)). Similarly to the component-asymmetric 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 component-asymmetric 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 inter-component coupling regime shows again the absolute dominance of the mode over the other lowest breathing modes besides .

## Vi Experimental realization

The few-body Bose-Bose mixture studied here should be observable with existing cold atom techniques. Quantum gas microscopes allow the detection of single particles in a well controlled many-body or few-body system Bakr et al. (2010); Sherson et al. (2010) and recent progress also allows for spin-resolved imaging in 1D systems using an expansion in the perpendicular direction Preiss et al. (2015); Boll et al. (2016). In these set-ups, the single-particle 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 time-resolved measurement of the breathing dynamics. Moreover, spin order was recently observed in very small fermionic bulk systems via spin-selective 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 spin-dependent transversal confinement as can be realized for heavier elements. Note that the longitudinal confinement needs to be spin-independent in order to ensure the same longitudinal trap frequencies and trap centres assumed in the calculations. The inter-component 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 inter-component background scattering lengths are negative Chin et al. (2010) leading to negative , but although not reported, inter-component FRs might exist. Alternatively, might be tuned via a CIR, which selectively changes at magnetic fields, where the intra-component 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 single-shots after a fixed hold time .

## Vii Discussion and Outlook

In this work we have explored a few-body problem of a Bose-Bose mixture with two atoms in each component confined in a quasi-1D 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 two-atom problem Busch et al. (1998) and include the information about the intra-species 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 low-lying energy spectrum as a function of the inter-species coupling for various fixed values of intra-species 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 relative-frame densities of some limiting cases to the known laboratory frame results García-March 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 relative-coordinate frame.

Finally, the obtained results were used to study the dynamics of the system under a slight component-symmetric quench of the trapping potential. We have derived expressions for the time evolution of the relative-frame 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 symmetry-induced selection rule: for the relative coordinates of both components and for the relative distance of the CMs of both components.

For not too strong inter-component coupling, each relative coordinate exhibits multi-mode 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 intra-component 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 relative-coordinate breathing frequency in the single-component 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 inter-component coupling and the symmetry of the intra-species interactions, which is in strong contrast to the single-component case Schmitz et al. (2013) as well as to the MF results, where two low-lying breathing modes can be obtained, namely an in-phase (out-of-phase) mode for a component-symmetric (component-asymmetric) quench Morise and Wadati (2000). Finally, we have argued that the experimental preparation of the considered few-body mixture and measurement of the predicted effects are in reach by means of state-of-the-art techniques.

This work serves as a useful analysis tool for future few-body 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 Bose-Fermi and Fermi-Fermi 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 Imaging-Structure, 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:

 ⟨a,b,c|H1|k,l,m⟩= 4∫drA∫drB ΦABa(rA+rB2)ΦABk(rA+rB2) φAb(rA)φAl(rA)φBc(rB)φBm(rB). (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 center-of-mass 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 so-called POTFIT algorithm Jäckle and Meyer (1996).

 χa,k(rA,rB) ≡ΦABa(rA+rB2)ΦABk(rA+rB2) ≈d∑i=0λ(a,k)iwA;(a,k)i(rA)wB;(a,k)i(rB).

Here333Note that in contrast to the usual convention we do not require the coefficients to be semi-positive without loss of generality., coincides with the th eigenvalue of the reduced one-body density matrix corresponding to the degree-of-freedom 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:

 ⟨a,b,c|H1|k,l,m⟩≈ 4d∑i=0λ(a,k)i∫drA w(a,k)i(rA)φAb(rA)φAl(rA) ∫drB w(a,k)i(rB)φBc(rB)φBm(rB). (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:

 ∫dr w(a,k)i(r) φσs(r) φσt(r), (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 intra-component 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 intra-component 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 time-consuming as one needs to calculate a bunch of 1D integrals for all the relevant quantum number configurations.

Now let’s analyse quantitatively the speed-up 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