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

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 [1,2], where research topics such as collective modes [3][4][5], binary mixtures [6,7] and lower-dimensional geometries [8][9][10] 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 [11][12][13]. 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) [14]: the miscible (M ), immiscible symmetry-broken (SB ) or immiscible core-shell structure, also called phase separation (PS ) [15][16][17]. Comparing the experimentally obtained densities to numerical MF calculations [18,19] provides a sensitive probe for precision measurements of the scattering lengths or, if known, the magnetic fields used to tune them [16]. Another possibility to access the interaction regime and thus the scattering lengths is by exciting the system and extracting the frequencies of lowlying excitations [20]. In contrast to a single-species case the collective modes of mixtures exhibit new exciting phenomena: doublet splitting of the spectrum containing inphase and out-of-phase oscillations, mode-softening for * mpyzh@physnet.uni-hamburg.de † pschmelc@physnet.uni-hamburg.de 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 [21][22][23].
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 [24][25][26][27]. In the early theoretical investigations on quasi-one-dimensional single-component condensates [28] it was shown that different interaction regimes can be distinguished by the breathing mode frequency, which has been used in experiments [10,29,30]. Furthermore, the monopole mode provides indirect information on the ground state [31], its compressibility [32] and the low-lying energy spectrum such that an analogy has been drawn to absorption/emission spectroscopy in molecular physics [27].
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) [8] can be employed to realize the Tonks-Girardeau limit [33,34], where the bosons resemble a system of non-interacting fermions in many aspects. While this case can be solved analytically [35,36], 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 few-body Bose-Bose mixture [37,38] 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 [39], 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 [39][40][41]. 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 [42]. 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 [43][44][45][46][47][48], reporting a transition from a two mode beating of the center of mass Ω CM and relative motion Ω rel 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 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 Bose-Bose mixture containing two components, which are labelled by σ ∈ {A, B}, 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 N σ = 2 particles, which have the same mass m σ ≡ m 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 a ho = /(mω) one arrives at the sim-plified Hamiltonian: where α the 3D s-wave scattering length and α ∈ {A, B, AB} are effective (offresonant) 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. [40], we pursue a different approach by using a correlated atom-pair basis. Thereby we can efficiently cover regimes of very strong intra-and intercomponent interaction strengths g α 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) [49]. 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 [50,51]. 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 [52][53][54].
In order to construct a tailored basis, which already incorporates intra-component correlations, let us neglect for a moment the inter-component coupling H AB , 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 [55] and boils down to a coordinate transformation and solving a Weber differential equation 1 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 nor-malized as well as symmetrized Parabolic Cylinder Function 2 (PCF) ϕ n (r) ∝ D µ(g,n) (|r|) of the relative coordinate with µ(g, n) being a real valued quantum number depending on the interaction strength g and excitation level n ∈ N 0 , which is obtained by solving a transcendental equation derived from the delta-type constraint: Coming back to the binary mixture problem we apply a coordinate transformation to the relative frame Y ≡ (R CM , R AB , r A , r B ) T defined by: • total CM coordinate The Hamiltonian attains a new structure in this frame. Firstly, the total CM is separated H = H RCM +H rem and is simply a HO of mass M = 4.
with the spectrum E CM n = n + 1/2 and n ∈ N 0 . The remainder of the Hamiltonian can be decomposed as H rem = H 0 +g AB H 1 , where H 0 = H RAB + σ H rσ can be solved analytically and H 1 couples the eigenstates of H 0 . H RAB is HO of mass M = 1 and H rσ lead to the above mentioned Weber differential equations with delta-type constraint.
Now to diagonalize H rem we choose as basis the eigenvectors of H 0 , which we label as |k, l, m with k, l, m ∈ N 0 . The energy of a corresponding basis function and its spatial representation are given by: 2 Dµ(|r|) = 2 µ 2 e − r 2 4 U (− 1 2 µ, 1 2 , 1 2 r 2 ) with U (a, b, x) the Tricomi's hypergeometric function where Φ AB k are HO eigenstates and ϕ σ i (r) ∝ D µ(gσ ,i) (|r|). All ϕ σ i (r) 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 H 1 , 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 [56], 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 [57]. 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 H 0 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 H in the relative frame.

Symmetry analysis
First of all, H commutes with the individual parity operators P Yi of the relative frame coordinates, [H, P Yi ] = 0. The eigenvectors of P rσ are restricted to even parity because of the bosonic character of our components. Due to the decoupling of H RCM it is sufficient to consider only the ground state of the total CM motion, which is of even R CM parity, in the following. Then, the parity of the R AB degree of freedom completely determines the total parity of the eigenstates. Another symmetry arises, if one chooses equal intra-component interaction strengths g σ . Under these circumstances the Hamiltonian is invariant under r A ↔ r B exchange, which we define as the S r 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 g AB for various fixed values of g A and g B . The total CM is assumed to be in its ground state. Fig. 1(a)   The decoupled total CM is assumed to be in its ground state. The total parity is thus determined solely by the RAB parity and is marked by black lines (even states) and red lines (odd states). In Fig. a) b) d) solid curves correspond to symmetric (+1) and dashed to antisymmetric (−1) eigenstates under the Sr operation. The indicated (avoided) crossings are exemplary and simply outline specific features. In b) and c), we label the excited states which are relevant for the lowest monopole excitations as discussed in Sect. V by the corresponding excitation frequencies ΩAB etc. The labels for the ground-state phases in the limiting regimes BEC-BEC, TG-TG, etc. follow the nomenclature of [40], see also the main text. All quantities are given in HO units.

0.
For g AB = 0 the Hamiltonian represents two uncoupled non-interacting bosonic species and we will label this regime as BEC-BEC following the nomenclature of [40]. The eigenenergies are integers with equal spacings of w, which is 1 in our units. In this limit the PCFs are even HO eigenstates of mass m = 1/2. The eigenenergies are 0,k,l,m = k + 2l + 2m + 2. For g AB = 0, the ith eigenenergy corresponding to only even (odd) R AB -parity eigenstates is (i + 2)!/(2!i!) fold degenerate with i ∈ N 0 . 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 ith even or odd R AB -parity branch, respectively. Note that this grouping of the energy levels into branches will be used in the following for all values of g AB 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 R AB -parity incidentally cross, as they are of different symmetry and consequently not coupled by the H 1 perturbation.
For very strong g AB values, i.e. in the composite fermionization (CF ) limit [37,38], we observe a restoration of degeneracies, but in a different manner, namely the lowest states merge pairwise forming a two-fold de-generacy. 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 g AB → ∞ are non-integer eigenvalues and unequal energy spacings.
A very similar analysis concerning this specific choice of interactions (g σ = 0 and arbitrary g AB ) was performed in [39], 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 N = 10 particles. However, the analysis only covered a single line of the (g A , g B , g AB ) parameter space. Here, we extend the ground state analysis of [40] 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 g ≡ g σ = 2. Already in the uncoupled regime (g AB = 0) we observe fewer degeneracies compared to the g σ = 0 case. Nevertheless, we group the eigenstates into branches of even / odd R AB -parity also for finite g σ by continuously following the eigenenergies to the g σ → 0+ 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 E (0) 0,k,l,m = k+µ(g, l)+µ(g, m)+2. To roughly estimate the energetic ordering it is sufficient to know that the real-valued quantum number µ(g, n) fulfils for 0 < g < ∞ the following relations: • 2n < µ(g, n) < 2n + 1 • µ(g, n) + 1 < µ(g, n + 1) < µ(g, n) + 2 meaning that a single excitation of the relative motion r σ is energetically below a double excitation of the R AB degree of freedom. E.g. the first even R AB -parity branch in the uncoupled non-interacting regime (BEC-BEC in Fig.  1(a)) contains three degenerate states: |2, 0, 0 , |0, 1, 0 and |0, 0, 1 (eq. (9)). By choosing finite g σ values |2, 0, 0 acquires a higher energy than |0, 1, 0 and |0, 0, 1 leading to reduced degeneracies in the spectrum. Another striking feature is the appearance of additional crossings between states of the same R AB -parity due to the S r symmetry. States, which possess different quantum numbers concerning the S r transformation (+1 or −1), are allowed to cross as they 'randomly' do throughout the g AB variation. Of course, such crossings are also present in the previous non-interacting case, it being also componentsymmetric. An avoided crossing between a state of the first even R AB -parity branch and a state of the second even R AB -parity branch is worth mentioning, which is present for all values of g σ (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 noncrossing rule [58].
In Fig. 1(c) we asymmetrically increase the intracomponent interactions g σ as compared to the noninteracting case, namely to g A = 1 and g B = 2. For the uncoupled scenario (g AB = 0) all the degeneracies are lifted, because now PCFs of the A and the B components are different. The energy is E The energetic state ordering is far from obvious, which becomes apparent upon closer inspection of the µ(g, n) function. E.g. consider again the first even R AB -parity branch. Its lowest energy state is a single excitation of r B , followed by a single excitation of r A . The highest energy of this branch corresponds to a double R AB excitation. The ordering pattern for higher order branches is even more complicated. For intermediate values of g AB we observe that crossings from the previous scenario (with g σ = 2) between states of the same R AB -parity are replaced by avoided crossings because of the broken S r 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 g σ = 100. When the g AB -coupling is absent, we have two hard-core bosons in each component. The system can thus be mapped to a twocomponent mixture of non-interacting fermions [35,59] 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 E (0) 0,k,l,m ≈ k + 2l + 2m + 4 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 ) [36] phase can be mapped to a non-interacting ensemble of four fermions with the ground state energy N 2 /2 = 8. 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 N !/(N A !N B !) = 6-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 [60] of N distinguishable particles, which interact with each other by an infinite repulsive delta-interaction. According to eq. (2) of Ref. [60], the snippet basis for distinguishable particles is defined as: where Ψ F 0 represents the ground state of N noninteracting fermions and Π is a permutation of particle coordinates, which defines the sector where x 1 , ..., x N |Π has support. Every snippet state is one of N ! possible ground states in the hard-core limit.
The last case we discuss is the highly asymmetric case g A = 0 and g B = 100 in Fig. 1(e). For g AB = 0 (BEC-TG) one expects, based on the previous considerations, integer eigenvalues E (0) 0,k,l,m ≈ k + 2l + 2m + 3 and thus equal spacings as for the cases g σ = 0 and g σ = 100 depicted in Fig. 1(a) and (d), because the PCF of the A component is an even HO eigenstate and the PCF of the B 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 [40], where the A component occupies the centre of the harmonic trap, while the B component, in order to reduce its intra-component interaction energy, forms a shell around the A component.

Relative-frame densities
Let us now inspect the relative-frame probability densities ρ 1 (Y i ) instead of the usually studied one-body densities ρ 1 (x σ ) of the laboratory frame as e.g. in [40]. 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 H 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 g AB ≪ 1.
We define these quantities as follows: where |E j is the j-th eigenstate of H and we trace out all the degrees of freedom of the relative frame Y except for one. Let's compare our results concerning the ground state densities for some limiting cases to the ones obtained in [40]. In Fig. 2 we show the densities for all the degrees of   freedom except for R CM , 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 ρ (0) 1 (r σ ) 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 [55], where ρ (0) 1 (r) develops a minimum in the centre of the density distribution for finite g whose value tends to zero as g → ∞.
The CF phase is in some sense a complete counterpart to the TG-TG case. Now ρ (0) 1 (R AB ) 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 ρ (0) 1 (x σ ) of the laboratory frame and one needs to additionally consider the two-body density function ρ (0) 2 (x A , x B ) to verify this behaviour [37]. The density distribution of ρ    1 (r B ) 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 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 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 R AB 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 [24,43], 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 dx σ x 2 σ ρ 1 (x σ , t) where we have omitted the subtraction of the mean value dx σ x σ ρ 1 (x σ , t) squared, which vanishes due to the parity symmetry. From a theoretical point of view, it is fruitful to define a breathing observable as σ,i x 2 σ,i , 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 |Ψ(t = 0) = |E 0 ω0 of some pre-quench Hamiltonian with frequency ω 0 1. Then a sudden quench is performed to ω = 1. The time evolution of this state is thus described as follows: where |E j ≡ |E j ω and c j = E j |E 0 ω0 is the overlap between the initial state and the jth eigenstate |E j of the post-quench Hamiltonian H. Since both the preand post-quench Hamiltonian are time-reversal symmetric, we assume their eigenstates and thereby also the overlap coefficients c j to be real-valued without loss of generality. A small quench ensures that |c 0 | ≈ 1 and only the n lowest excited states are of relevance. Symmetry considerations further reduce the number of allowed contributions. E.g. states of odd R CM -parity or odd R AB -parity have zero overlap with |E 0 ω0 , because the initial state is of even R CM -and R AB -parity and the quench does not affect any of the symmetries discussed in the previous section. Similarly, in the component-symmetric case g A = g B , states, which are antisymmetric w.r.t. the S r operation, have no overlap with the pre-quench ground state being symmetric under S r . 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 relativeframe coordinates. Employing the expansion in postquench eigenstates from eq. (13), their time-evolution may be approximated within the linear-response regime as neglecting terms of the order c i c j for i, j > 0. So ρ 1 (Y i , t) can be decomposed into the stationary background ρ In the following, we regard the excitations of the first even R AB -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 g AB , 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 g A = g B (g A = g B ) 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 componentsymmetric case of Fig. 1(b).

A. Component-asymmetric case
Because of the low amplitude quenching protocol, we will excite four breathing modes simultaneously in the component-asymmetric case (g A = 1, g B = 2). Three of them stem from the first even R AB -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 ∆E = n with n ∈ N. 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 R AB -parity branch ensuring a considerable overlap with the initial state. The total CM trivially oscillates with the constant frequency Ω CM = 2 independent of any interactions g α it being a decoupled degree of freedom with the single-particle HO Hamiltonian (eq. (5)) [25,43]. 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 |n for a complete characterization of the involved states. States of higher order even R AB -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 g AB = 0, one can show analytically that each relative-coordinate density oscillates with a single frequency, each corresponding to exactly one eigenstate of the first even R AB -parity branch (see Fig.  3 (a)-(c)). E.g. for ρ 1 (R AB , t), the only transition density ρ (0,j) 1 (R AB ) which survives taking the partial trace is the one corresponding to |0 |2, 0, 0 , while the contributions from the remaining excited states vanish. This leads to the breathing motion in the R AB coordinate with a single frequency Ω AB . Analogously one can show that |0 |0, 1, 0 solely induces density modulation in ρ 1 (r A ) with the frequency Ω A , while ρ 1 (r B ) oscillates with Ω B exclusively due to |0 |0, 0, 1 . 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 ρ 1 (R CM , t), 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: The expectation value of each observable with respect to the time-evolved state |Ψ(t) is directly related to the respective relative-frame density: 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 jth mode is essentially composed of the overlap c j and of the transition element: In order to evaluate the overlaps c j = E j |E 0 ω0 with j = 0 in terms of only the post-quench Hamiltonian eigenstates, we perform a Taylor approximation with respect to the weak quench strength δω = ω − ω 0 , namely |E 0 ω0 ≈ |E 0 ω − δω d dω |E 0 ω evaluated at ω = 1, and arrive at Applying the (off-diagonal) Hellmann-Feynman theorem, one obtains: 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 A j , with which the j-th mode contributes to the oscillation of the observable Y 2 i : which may be interpreted as the susceptibility of the Y 2 i observable for the excitation of the state |E j . In Fig. 4 (a) we show the values of possible breathing mode frequencies, obtained from the spectrum of Fig. 1  (c). Ω CM = 2 does not depend on any interactions g α . In contrast to this, Ω AB is degenerate with Ω CM for g AB = 0, and when increasing g AB , decreases to a minimum first and then increases with the tendency to asymptotically reach Ω CM again. This behaviour strongly resembles the dependence of the relative-coordinate breathingmode frequency in the single-component case [43]. Ω A and Ω B have qualitatively akin curve shapes, varying much stronger with g AB . 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 R ABparity branch leads to a breathing motion of some specific relative-frame coordinate. Indeed, if we look at the amplitudes A j in Fig. 4 (b)-(d) in the decoupled regime (g AB = 0), we recognize that the amplitude for the coordinate Y i is non-zero only for one mode, namely the one with which ρ 1 (Y i , t) 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 g AB variation.
Another peculiarity worth noting arises in the strongly coupled regime: Ω σ oscillations become strongly suppressed for all the observables making Ω AB and Ω CM 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 Ω AB . We highlight that in the CF regime each density peak of the bimodal distribution ρ 1 (R AB , t) does not only breath periodically but also its maximum height position performs dipole-mode like oscillations, see Fig. 3 (d).

B. Component-symmetric case
Now we compare the above results with the component-symmetric case of Fig. 1 (b), where g A = g B = 2. Fig. 5 (a) depicts the possible breathing-mode frequencies. Here, two main differences arise: First, the Ω AB curve features two minima due to two avoided crossings with an eigenstate of the second even R ABparity branch. Second, the other two breathing mode frequencies of the relative coordinates are degenerate for g AB = 0, then separate with increasing g AB and approach one another asymptotically. Instead of using the labels Ω A and Ω B 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. S r , 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 r A and r B coordinates is identical due to the imposed component-symmetry and in the uncoupled regime (g AB = 0) they both oscillate solely with the Ω + frequency (see Fig. 5 (c) and (d)). Similarly to the component-asymmetric case, we see that when increasing g AB the Ω AB (Ω + ) mode contributes also to the observable r 2 σ (R 2 AB ). However, in the intermediate interaction regime 3.7 g AB 8.3 we observe a strong suppression of the Ω AB mode contribution for all observables, which is in stark contrast to the component-asymmetric case. Instead, a state of the second even R AB -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 g AB 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 Ω AB mode over the other lowest breathing modes besides Ω CM .

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 [62,63] and recent progress also allows for spin-resolved imaging in 1D systems using an expansion in the perpendicular direction [64,65]. In these set-ups, 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 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 [66].
Deterministic preparation of very small samples was demonstrated for fermions via trap spilling [67] and for bosons by cutting out a subsystem of a Mott insulator [68]. 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 7 Li, 39 K or 87 Rb. 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 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 [69]. In the case of 7 Li and 39 K, the inter-component back-ground scattering lengths are negative [14] leading to negative g AB , but although not reported, inter-component FRs might exist. Alternatively, g AB might be tuned via a CIR, which selectively changes g AB at magnetic fields, where the intra-component scattering lengths are very different.
In the following we give concrete numbers for a choice of 7 Li. The density distribution has structures on the scale of the HO unit a ho (Fig. 2). Choosing the trap parameters as in Ref. [67] as ω/(2π) = 1.5 kHz and ω ⊥ /(2π) = 15 kHz, yields a ho = 1 µm, i.e. larger than a typical optical resolution of 0.8 µm. At the same time, temperatures much lower than ω/k B = 72 nK are state of the art. The breathing dynamics will occur on a time scale of several 100 µs, 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 R 2 CM t , R 2 AB t , r 2 σ t by averaging the occurring relative coordinates over many single-shots after a fixed hold time t.

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 [55] 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 g A and g B and related the ground state relative-frame densities of some limiting cases to the known laboratory frame results [40]. 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 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 (g AB = 0) the density of each relative frame coordinate performs breathing oscillations with a single frequency corresponding to a specific excited state of the first even R AB -parity branch of the spectrum. The total CM coordinate performs breathing oscillations with the frequency Ω CM = 2 (HO units). For asymmetric choices of g σ values, three additional monopole modes participate in the dynamics, each of them corresponding to the motion of a particular relative coordinate: Ω A for the relative coordinate of the A component, Ω B for the relative coordinate of the B component and Ω AB for the relative distance of the CMs of both components. In contrast to this, the symmetric case g A = g B leads to only two additional modes because of a symmetry-induced selection rule: Ω + for the relative coordinates of both components and Ω AB 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: Ω AB and Ω CM . We have observed the same effect for the case of Phase Separation (results not shown). Interestingly, the dependence of Ω AB on g AB strongly resembles the behaviour of the relative-coordinate breathing frequency in the singlecomponent case [43]. 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 [43] as well as to the MF results, where two low-lying breathing modes can be obtained, namely an in-phase (outof-phase) mode for a component-symmetric (componentasymmetric) quench [23]. 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.