Dynamics of a degenerate Cs-Yb mixture with attractive interspecies interactions

We probe the collective dynamics of a quantum degenerate Bose-Bose mixture of Cs and $^{174}$Yb with attractive interspecies interactions. Specifically, we excite vertical center of mass oscillations of the Cs condensate, and observe significant damping for the Cs dipole mode, due to the rapid transfer of energy to the larger Yb component, and the ensuing acoustic dissipation. Numerical simulations based on coupled Gross-Pitaevskii equations provide excellent agreement, and additionally reveal the possibility of late-time revivals (beating) which are found to be highly sensitive to the Cs and Yb atom number combinations. By further tuning the interaction strength of Cs using a broad Feshbach resonance, we explore the stability of the degenerate mixture, and observe collapse of the Cs condensate mediated by the attractive Cs-Yb interaction when $a_{\mathrm{Cs}}<50 \, a_0$, well above the single-species collapse threshold, in good agreement with simulations.

Superfluid mixtures formed of ultracold atoms offer a versatile playground for studies of binary fluid dynamics in both Bose-Bose and Bose-Fermi mixtures. Interplay between superfluids can lead to counterflow instabilities [1][2][3][4], and quantum turbulence [1,5]. Previous work has explored collective excitations in the form of scissors mode oscillations in a repulsive Bose-Bose mixture [6], and dipolar oscillations [7], breathing dynamics [8], and the exchange of angular momentum [9] in Bose-Fermi mixtures. Experimental control in ultracold atomic mixtures allows the balance of the inter-and intraspecies interactions to be tuned, enabling further manipulation of the coupled dynamics and stability of the system. This can result in novel bound states such as bright-dark solitons [10], and exotic vortex lattices [11,12], and enable miscibility studies [13]. Prior research on dynamical instabilities associated with attractive interspecies interactions has focused on Bose-Fermi mixtures, where the bosonic species mediates collapse and subsequent density-dependent atom loss in the fermionic species [14][15][16]. However, recent studies have used Bose-Bose mixtures to probe the region where the net mean-field interaction is close to zero, resulting in the observation of quantum droplets where quantum fluctuations provide a repulsive force that stabilizes against collapse [17][18][19].
In this Letter, we probe the collective dynamics and stability of a quantum degenerate Bose-Bose mixture of Cs and 174 Yb with an attractive interspecies interaction characterised by the scattering length a CsYb = −75 a 0 [20]. The Cs intraspecies scattering length a Cs is a tunable experimental parameter via a broad, lowfield magnetic Feshbach resonance [21,22], whilst the * These two authors contributed equally. † s.l.cornish@durham.ac.uk 174 Yb intraspecies scattering length is a Yb = 105 a 0 . We use vertical Cs center-of-mass (CoM) oscillations to probe the collective dynamics of the system, observing an upwards mean-field shift of the CoM dipole mode frequency. We also find that the attractive interspecies interaction results in complicated coupled dynamics, as Cs pulls along Yb in the region of Cs-Yb overlap. We probe the stability of the mixture by tuning the Cs interaction strength across the collapse threshold, demonstrating Yb-mediated collapse of the Cs Bose-Einstein condensate (BEC), well above the single-species collapse threshold. We prepare dual-degenerate mixtures of 133 Cs in the |F = 3, m F = +3 hyperfine state and 174 Yb using the apparatus and procedures described in Refs. [23][24][25][26][27]. Our dual-species evaporation relies on sympathetic cooling of Cs with Yb, resulting in an imbalanced mixture with typical condensate atom numbers of (N Yb , N Cs ) = (50, 5) × 10 3 with up to 20% shot-to-shot variation [27]. The condensates are confined in a bichromatic optical dipole trap (BODT), using laser beams with wavelengths of 532 nm and 1070 nm aligned as shown in Fig. 1(a). We define a coordinate system (x, y, z) with x along the copropagating BODT beams and z the vertical direction. The axesx andỹ for the Cs trap are rotated 20 • from the Yb trap axes. The final trap frequencies are (ω x , ω y , ω z ) Yb = 2π × (10, 120, 80) Hz for Yb and (ωx, ωỹ, ω z ) Cs = 2π × (40, 70, 260) Hz for Cs. The vertical trap frequencies are very sensitive to the overlap of the BODT beams and can vary by up to 10% with day-to-day realignment of the trap [27].
We first explore the coupled dynamics of the mixture by exciting Cs CoM oscillations in the vertical direction. At the end of the dual-species evaporation sequence we adiabatically ramp on a magnetic field gradient of 9.0 G/cm in the z-direction. This additional magnetic potential causes a 0.96 (7) µm vertical shift of the Cs trap minimum, without directly acting on Yb which has no magnetic moment. The gradient is rapidly switched off, exciting the vertical CoM dipole mode. We wait a variable hold time before turning off the trapping potential, then perform absorption imaging of the Cs condensate after 40 ms of levitated time of flight (ToF). We hold a Cs constant at 275 a 0 for both in trap oscillations and the subsequent ToF. To account for the high sensitivity of the vertical trap frequencies to alignment, we reference the coupled Cs-Yb dynamics to the Cs CoM oscillations in the absence of Yb by taking both measurements in quick succession. We obtain the Cs-only reference by selectively removing Yb with a 5 ms pulse of light resonant with the 1 S 0 → 1 P 1 transition at 399 nm prior to exciting the Cs CoM oscillations.
The measured vertical CoM oscillations of the Cs BEC in the presence of the Yb BEC are shown by the green circles in Fig. 1(b). These oscillations are strikingly different from the corresponding Cs-only oscillations (dashed green line), which occur at the Cs natural trap frequency ω Cs,z = 2π × 251(1) Hz. In the presence of Yb, the Cs CoM oscillations experience a frequency shift upwards of 4(1)%, with significant damping. This shift is in good agreement with a simple mean-field prediction that assumes a Thomas-Fermi profile for the Yb BEC and neglects the back-action of the Cs BEC on the Yb atoms due to the large number imbalance [28]. Under such assumptions, the interspecies interaction results in an additional attractive quadratic potential for Cs that leads to a 5% upwards frequency shift.
We perform numerical simulations based on coupled 3D Gross-Pitaevskii equations (GPE) [28] to gain further insight into the rich oscillatory dynamics. A numerical visualization of the initial geometry with (N Yb , N Cs ) = (40, 6.5) × 10 3 atoms, ω Cs,z = 2π × 250 Hz, and an imposed +0.95 µm vertical shift of the Cs trap position is shown in Fig. 1(c). The numerical simulations reveal that both components start the coupled dynamics from non-equilibrium positions with respect to their trap centers. The attractive Cs-Yb interaction induces a density 'bulge' in the Yb BEC in the region of Cs-Yb overlap, visible in the z-axis density profile shown in Fig. 1(c). The off-center bulge leads to a vertical Yb CoM shift upwards of ∆z Yb,CoM = 0.1 µm. Similarly, the Cs CoM shift reduces to ∆z Cs,CoM = 0.88 µm.
The simulated vertical CoM oscillations of the Cs BEC are shown by the red line in Fig. 1(b). The simulations indicate that the CoM shifts and ensuing coupled dynamics are highly sensitive to the relative atom numbers [28]. Therefore, to account for shot-to-shot variations in both Cs and Yb atom numbers, the simulation result is based on the average of 206 atom number combinations, randomly selected such that their mean and standard deviations are (N Yb , N Cs ) = (40, 6.5) × 10 3 and 20%, respectively, in broad agreement with the values measured experimentally. As shown in Fig. 1(b), this leads to excellent agreement with the experimental observations. A video of the coupled dynamics is available in the Supplemental Materials (SM) [28]. Closer analysis of the  (2) Hz, and β = 0.10(2) ms −1 . The dashed green line shows an equivalent fit to a measurement performed in the absence of Yb, yielding ω = 2π × 251(1) Hz [28]. The numerical plot is based on an average over 206 particle number sets, with the red band indicating the standard deviation. Atom numbers are randomly selected such that their mean is (N Yb , NCs) = (40, 6.5) × 10 3 , with a 20% standard deviation.
simulations reveals the dynamics are complicated by the different vertical trap frequencies experienced by the two components (ω Cs,z ≈ 3 ω Yb,z ), the large atom number imbalance (N Yb /N Cs ≈ 6), and the broader spatial extent of the Yb BEC in the weakly-confined direction. Accordingly, in what follows we focus our analysis on the central bulge in the Yb BEC that directly interacts with the Cs BEC. A striking feature of Fig. 1(b) is the fast damping of the Cs CoM oscillations in the presence the Yb BEC, with the amplitude falling to 1/e its initial value in 10(2) ms. We explore the origin of this damping in Fig. 2. Figure 2(a) shows the numerical CoM motions out to 180 ms, whilst Fig. 2(b) shows the time evolution of the relevant energies. Looking first at the energies in Fig. 2(b), we see that there is a net loss of ∼ 85% of the initial excess Cs potential energy over the first ∼ 20 ms. This energy is transferred to the Yb component through a coupling of the CoM motions to low-lying excitations. Although the initial excitations are produced along the z axis of relative motion, an additional secondary mode coupling between the z and y directions gradually develops owing to the comparable Yb trap frequencies in these directions. This is evident in Fig. 2(b) as an out-of-phase oscillation in  , and the main corresponding single-particle energy changes along the z and (for Yb only) the y axis. Note the more gradual initial increase along y, consistent with a secondary excitation mechanism, and the initial out-of-phase oscillations between z and y directions for Yb. (c) Time evolution of the Cs velocity profile vCs = ∇ arg[ψCs(0, 0, z)]/mCs. (d) Corresponding Yb velocity profile. Velocity plots include a grey mask to filter out regions corresponding to atomic densities lower than 2 µm −3 . E Yb,y and E Yb,z , correlated with oscillations of the condensate widths. The end result is that after about 20 ms, the initial excess Cs single-particle energy concentrated along z has been converted to Yb single-particle energy across both the z and y directions in comparable fractions. Smaller, secondary couplings along the x direction (the weak axis of the trap) occur at longer timescales [28].
Earlier related studies of repulsively-coupled nonlinear superfluids [8,[29][30][31] have identified two distinct excitation mechanisms. Firstly, in a number-imbalanced system like ours, the motion of the smaller superfluid through the larger component can lead to single-component excitations in the larger superfluid. For the attractive interaction in our system, this is analogous to stirring a single-component superfluid with a focused reddetuned laser beam [32,33]. Secondly, hydrodynamic considerations in the homogeneous limit reveal the existence of a dynamical instability whenever the relative motion of two superfluids exceeds some characteristic speed. In the limit of small interspecies interactions, this speed approaches the sum of the sound velocities of the two (uncoupled) superfluids. This counterflow instability leads to the simultaneous generation of pairs of excitations in both superfluids [29]. We believe both these mechanisms are at play in our experiment [28].
The temporal evolution of the local velocity fields shown in Fig. 2(c) for Cs and Fig. 2(d) for Yb is useful for visualizing the two excitation regimes. The velocity fields along the z axis for each component α = {Cs, Yb} are obtained from v α,z = ∇ arg[ψ α (0, 0, z)]/m α , where m is the mass and ψ is the condensate wave function. During the initial half oscillation (0 < t < T Cs /2), Cs exhibits coherent harmonic motion, maintaining a constant phase profile across its vertical width. In contrast, during this same initial period, incoherent features gradually emerge in the Yb component in the form of counterpropagating density waves as the relative velocity between the two components increases, consistent with the single-excitation picture. Here the Cs BEC is both smaller than Yb, and executes faster oscillations. Therefore, it acts as a perturbing attractive barrier, which travels through the Yb BEC, and drags along the central Yb density bulge. As shown in Fig. 2(d), the bulge begins to track the downward Cs motion at t ∼ 1 ms (T Cs /4), which corresponds to the point of maximum downward Cs velocity. This in turn causes Cs to slow down as it continuously transfers energy to Yb.
Counterflow dynamics occur when Cs reverses its motion at t ∼ 2 ms (T Cs /2) whilst the central Yb bulge continues its downward motion, thus creating a direct counterflow between the two components, whose magnitude increases over the next quarter-cycle of the Cs oscillation (until t = 3T Cs /4). This leads to the onset of counterflow instabilities, which are particularly pronounced at the edges of the mixture's central overlap region (|z| 2 µm), where the relative speed exceeds the sum of the local critical velocities of the two components. Moreover, the moving Cs component induces a secondary, at times opposing in direction, flow within the central Yb bulge, thus giving rise to the co-existence of flows in the same and opposing directions (the extension of co-moving and counter-moving modes for two superfluids which are not in the same trap [13]). The ensuing Yb motion, which extends beyond the region of CsYb overlap, feeds back into Cs by exciting phonons at the edges of the Cs component where its density is lowest [28], leading to a complicated energy transfer and phonon excitation mechanism in both components. We have also verified that smaller initial amplitude oscillations give rise to reduced damping and fewer excitations.
Although the dominant dynamics occur during the experimentally probed timescale of 22 ms, closer inspection of the long-time coupled evolution for a fixed total atom number reveals smaller scale energy transfers between the Cs and Yb components. This results in partial revivals in the Cs CoM oscillations, as shown in Fig. 2(a). Such beating arises from the existence of two closelylocated dominant frequencies for Cs (both within a few % of the natural Cs frequency), the origin of which can be attributed to the co-existing secondary dynamical micromotions (combinations of appropriate co-moving and counter-moving modes) in the central z region. The combination of rapid damping and frequency sensitivity to particle number combinations make it hard to characterize such secondary effects in our current experiments. More details of these complicated dynamics are given in the SM [28].
In a separate study, we explore the dynamical instabilities of the mixture further by tuning the relative meanfield contributions through changes to the Cs scattering length, a Cs . The attractive interspecies interaction between Cs and 174 Yb gives rise to mean-field collapse of the Cs BEC mediated by the presence of the Yb BEC, as shown in Fig. 3. In both single-component and dualspecies BECs, the onset of collapse is related to a mechanical instability arising from the population of imaginary Bogoliubov modes [34]. In a single-component BEC, collapse [35][36][37][38][39], along with its associated density-dependent atom loss, occurs when the intraspecies interaction becomes sufficiently attractive to overcome the zero-point kinetic energy associated with the harmonic trapping potential. For a given atom number N , the scattering length a crit associated with the onset of collapse is given by |a crit | = C /(mω)/N , whereω is the geometric mean of the trap frequencies, and C is a numerical constant weakly dependent on the trap geometry [37,40,41]. In contrast, the dual-species collapse threshold primarily depends on the balance of the inter-and intraspecies interactions, and the critical point for collapse is determined by the parameter δg = g 12 + √ g 11 + g 22 which describes the balance of mean-field interactions in the system. Here, g ij = 2π 2 a ij (m i + m j )/(m i m j ) is the interaction coupling constant. When δg < 0, the modes of the density branch of the Bogoliubov excitation spectrum become imaginary, leading to the onset of mechanical instabilities, and subsequent collapse [14][15][16].
In our experiment, we approach the collapse instability by varying a Cs via a magnetic Feshbach resonance. We begin with a dual-species BEC formed with a Cs = 147 a 0 , with atom numbers (N Yb , N Cs ) = (50, 4) × 10 3 . We ramp a Cs to the desired value over 10 ms, hold for 30 ms, and then ramp back to a Cs = 147 a 0 over 10 ms prior to ToF expansion and imaging. As shown in Fig. 3, we observe the onset of collapse as a marked decrease in Cs atom number due to density dependent 3-body losses (K 3,Cs ∼ 1 × 10 −27 cm 6 s −1 [26,42] for the range of magnetic fields explored). lengths, both with (top row) and without (bottom row) Yb present. The presence of Yb causes the collapse instability to occur at a Cs > 0, where we would expect a stable single-component Cs BEC. In addition to atom loss, we also observe a change in aspect ratio for the Cs cloud. We quantify the collapse in Fig. 3 We model the experimental protocol numerically using two-component coupled GPEs with a density dependent 3-body loss term for Cs, and ω Cs,z = 2π × 260 Hz [28]. The simulation results are shown in Fig. 3(b) for Cs+Yb (green squares) and Cs only (red squares). For Cs scattering lengths lower than the range shown, the simulations do not converge at the highest resolutions we could feasibly use, indicative of mean-field collapse [28]. We find good agreement with the experimental results using a thermal Cs three-body loss coefficient of K 3,Cs = 1.8 × 10 −27 cm 6 s −1 . However, we note that the Ybmediated collapse is broader than predicted by the GPE simulations. This may be indicative of beyond-meanfield effects, such as the formation of quantum droplets recently observed in both K spin mixtures [17,18] and Rb-K mixtures [19], and is a direction for future studies.
In conclusion, we have probed the collective dynamics of a quantum degenerate Bose-Bose mixture of Cs and 174 Yb with attractive interspecies interactions. Observations of the CoM oscillations of Cs in the presence of Yb have revealed a measurable increase in the dipole mode frequency and significant damping as energy is transferred to Yb. We find good agreement with numerical simulations which highlight the complexity of the coupled dynamics. We have also investigated the stability of the degenerate mixture by crossing the transition for dualspecies collapse and find good agreement with coupled-GPE simulations that include a Cs three-body loss term.
This opens up new prospects for the study of the BEC to quantum droplet phase transition. In this context, a key advantage of Cs-Yb mixtures is that the Cs scattering length may be precisely tuned without affecting Yb. Additionally, the very different atomic structure of Cs and Yb enables species-specific optical traps with low photon scattering rates. These advantages may prove invaluable for precise studies of binary fluid dynamics in both Bose-Bose and Bose-Fermi Cs-Yb mixtures.

Supplemental Materials
In these Supplemental Materials we present details of our numerical modelling, offering further evidence to support the findings made in the main text.
We start by describing our theoretical model which is based on a three-dimensional mean-field description by the coupled Gross-Pitaevskii equations (GPEs), and where ψ Yb(Cs) is the mean-field wavefunction of Yb (Cs) trapped in the harmonic potential, V Yb(Cs) (r) and M µ = 2m Yb m Cs /(m Yb + m Cs ) is the reduced mass. The Cs three-body loss coefficient for the condensate is K (C) 3,Cs while the corresponding thermal coefficient K 3,Cs is a factor of 3! higher [43]. The wavefunctions are normalized to the total atom numbers, The trapping potentials are considered as and where the 'tilde' notationx andỹ axes in V Cs indicates that the longitudinal (x) and one of the transverse (y) axes of Cs are rotated by 20 degrees with respect to the x and y axes of Yb while z axis remains coincided, in agreement with the experimental setup, and ∆z is the trap displacement in Cs.
In the following, we investigate the coupled oscillations and collapse via GPE modeling. The GPEs are simulated by the Fourier pseudo-spectral methods with the Runge-Kutta method. In the simulations, the trapping frequencies are as described in the experimental section and we take a Yb = 105 a 0 and a YbCs = −75 a 0 , where a 0 is the Bohr radius. In the coupled oscillation section, we consider fixed a Cs = 275 a 0 with time-dependent ∆z and ignore three-body loss, whereas in the collapse section, a Cs is considered as a time-dependent variable with ∆z = 0, and we include the Cs condensate three-body loss term.

I. COUPLED OSCILLATIONS
The coupled oscillations of Cs and Yb BECs triggered by the sudden removal of the trap displacement in Cs, are modeled by the GPEs without the inclusion of the threebody loss term, based on an initial Cs trap displacement ∆z = 0.95 µm, and ω Cs,z = 2π × 250 Hz. The initial trap displacement in Cs is implemented in the stationary initial condition, generated by solving Eq. (2) with the imaginary-time propagation method. Subsequently, ∆z is set to zero for all t ≥ 0 to mimic the experimental procedure. The required set of particle numbers for the initial conditions is reached by a self-consistent choice of the chemical potentials µ Cs(Yb) [44] during the imaginarytime evolution.

A. Time-of-Flight Analysis
In the experiments, the center of mass (CoM) motion is measured by time-of-flight (ToF) imaging. Since a full ToF simulation is numerically forbidden due to the huge computational size constraints imposed by the V Cs(Yb) anisotropy, our numerical simulations focus on the characterization of the in situ CoMs, with an effective mapping to the ToF measurement along the kick axis. The in situ CoM along the z direction is defined by The presence of the attractive Cs-Yb interaction implies that the actual initial shift exhibited by the CoM of Cs is reduced, while simultaneously the Yb density is locally distorted by Cs, acquiring an initial CoM shift. All the results presented here are based on a definite particle number combination. For our 'primary' atom number combination of (N Yb , N Cs ) = (40, 6.5) × 10 3 we find such initial shifts to be z CoM Cs (t = 0) ∼ 0.88 µm, in comparison to the imposed trap shift (∆z = 0.95 µm), and z CoM Yb (t = 0) ∼ 0.10 µm. This 'primary' atom number combination is in broad agreement with the values measured experimentally, and leads to excellent agreement with the experimental observations [see Fig. 1 Direct comparison to the experimental ToF data is facilitated by mapping the numerical in-situ results to the predicted ToF measurement by based on an expansion time τ ToF = 40 ms, and on the assumption that the condensate carries its moving velocity during the ToF expansion. The time derivative is computed by the central finite difference method, namely, dz CoM Cs (t)/dt ≈ [z CoM Cs (t + 0.5∆t) − z CoM Cs (t − 0.5∆t)]/∆t, with ∆t = 0.5/ω Cs,z . This mapping provides an excellent agreement of numerical predictions with experimental observations as shown in Fig. SM.1.
Specifically, Fig. SM.1(a) shows the Cs CoM oscillation in the absence of the Yb component both based on experimental measurements (green points/line) and corresponding numerical predictions (red line), which yields (atom-number-independent) oscillations at the expected trap frequency ω Cs,z . The corresponding comparison for the Cs CoM oscillations in the presence of Yb is shown in Fig. SM.1 (b)-(e). Specifically, the same set of experimental data (green points/line) are compared against numerical simulations (red lines) using 4 different atom number combinations, chosen here to reflect number variations of up to ≈ 1.25σ in either Cs (N Cs = (6.5 ± 1.5) × 10 3 ) or Yb (N Yb = (40 ± 10) × 10 3 ). These reveal very good overall agreement over the range of available experimental data, but with a sensitivity of the late-time [t > (15 − 20) ms] oscillations to the precise atom numbers (and in particular to the Yb atom number, which is the largest component). [Recall that Fig. 1(b) of the main paper is based on weighted combinations of different atom numbers, statistically distributed around the 'primary' mean value.] The effect of the tilted axes in Cs (α = 20 • ) is also investigated (dashed lines in Fig. SM.1 (b)-(e)), and it is found to have a marginal influence on the CoM motion (at least over the probed timescales), clearly ruling out such axes mismatch as a primary reason for the observed dynamics.

B. Excitation Dynamics Analysis
Here we analyze the dominant collective mode frequencies, provide information regarding the numericallyobserved excitations and explain the physical origin leading to the beating at longer times.

Excitation Profile
As discussed in the main text (see also Fig. 2(c) of the main paper) the dominant dynamics, energy transfer and dissipation manifest themselves over the first ∼ 22 ms. In this subsection we try to shed more light on the complex key dynamical features during this evolution phase. In order to do this, we present in Fig. SM.2 the study of a number of relevant observables, both over the first 22 ms dominant timescale (left panels), and also zoomed in to examine the initial system behaviour over the first 5 ms (right panels), which enable a more clear visualization over the first complete oscillation cycle, which seeds the subsequent dynamics.
Given that Cs has a smaller spatial extent than Yb along both z and x axes and thus only interacts with the central part of Yb, we can expect minor modifications between the central part of the Yb component (which is co-located with Cs) and the corresponding behaviour of the entire Yb cloud, due to excitations propagating towards the Yb trap edges along the z direction. We can characterize this by looking at the CoM of Yb when measured either as a whole (as done in the main text), or when instead focusing only on the relevant region of direct interactions with Cs. Specifically, we contrast the earlier findings to the 'local' Yb CoM oscillations extracted within the central region |x| < 10 µm, |y| < 6 µm and |z| < 3 µm. These features, plotted in Fig. SM.2 (a), reveal some evidence for secondary (non-central) dynamics in the Yb CoM, as a result of micro-motions in Yb within and beyond its region of co-existence with Cs. Such micro-motions, which we will show can even feature oscillations in opposing directions for the central and total Yb component, are associated with the emerging excitations.
To determine the difference between the dynamical and initial (static) profiles for each component (which is a measure of the extent of excitations present in such systems [45]), we measure their overlap for each component through the 'Fidelity' defined here as Here ψ ref Cs(Yb) (r(t)) = ψ Cs(Yb) r − ∆z CoM Cs(Yb) (t)ẑ, t = 0 is constructed from the initial wavefunction of each component [ψ Cs(Yb) (r, t = 0)], by shifting it along the z direction by an amount dependent on its CoM difference ∆z CoM Cs(Yb) (t) = z CoM Cs(Yb) (t) − z CoM Cs(Yb) (0), and computed numerically via interpolation.
Uncoupled oscillations for the two components give a Fidelity ≈ 1 (with some undamped small amplitude oscillations), as demonstrated by the thin/light pink/blue lines in Fig. SM.2(b). Contrary to this, the Fidelity begins to drop rapidly for both components after ≈ 1 ms, with an evident few % drop maintained at all subsequent evolution times, clearly demonstrating the impact of incoherent excitation processes during the first coupled oscillation cycle.
We attribute this feature to the time-evolving relative speed between the two components. At the outset, the two components are moving in the same direction towards decreasing values of z, but starting with different offsets and with highly distinct trap frequencies (specifically along the direction of excitation, the frequencies differ by a factor of ≈ 3), signalling very different oscillation dynamics. As such, the Cs CoM acquires negative values already around ≈ 1 ms (as opposed to ≈ 2.7 ms for Yb), with Cs reversing its direction of motion already at ≈ 1. the local velocity fields along the z direction for Cs and Yb, respectively. These plots include a grey mask to filter out the information for densities lower than 2 µm −3 , as this enables us to better highlight the most relevant features.
Initially, and for the first half oscillation cycle (up to ≈ 1.9 ms), Cs exhibits coherent motion, with a uniform (but time-dependent) velocity profile. Contrary to this, the Yb component demonstrates gradual emergence of incoherent features during this early period, in the form of counter-propagating density waves arising from the relative motion between the two components. During this timescale, we can thus visualize the initial excitation process as a single-excitation in Yb, caused by the coherent Cs component moving through the central Yb region, which acts much like a red-detuned laser being accelerated through Yb.
At ≈ 1 ms, the Cs CoM reaches the trap centre and the coherent motion in Cs reaches the maximum velocity for   its harmonic oscillation. At 1.9 ms, while Yb is still mov-ing towards lower z values, the Cs CoM reverses its direc-tion, thus creating a direct counterflow between the two components. The magnitude of the counterflow is gradually increased and maintained at high values over the next ≈1 ms when the Cs oscillates to its center for a second time at 2.9 ms. This is clearly visible in Fig. SM.2(e), which depicts the spatial-temporal evolution of the (conveniently scaled) relative velocities. Here, the local velocities along the z axis for each component α (= Cs, Yb), are obtained from v α = ∇arg[ψ α ]/m α . Their relative velocity is shown here scaled to the sum of the instantaneous local speeds of sound, which can be visualized as a potential upper bound for the critical velocity of the coupled system. This choice is based on the predicted critical velocity of counterflow for two weakly-coupled superfluids of (c Cs + c Yb ) [29,31], where c α = g α n α /m α with g α = 4π 2 a αα /m α is the critical sound speed for each individual component, and the density n α is measured locally at each moment in time.
The establishment of direct counterflow between the Cs and Yb components leads to the emergence of incoherent features in the propagation of Cs, and coinciding secondary excitations of Yb. Such larger relative speeds are generally reached towards the edges of the Cs component, so typically for |z| 3µm, suggesting such secondary excitation processes originate at the edges of the Cs component. These features therefore appear consistent with the generation of coupled excitations in the two components, as a result of the counterflow instability near the edges of the component overlap where the densities, and thus speeds of sound, are appropriately reduced.
Interestingly, at ≈ 3.7 ms, when the Cs CoM next reverses its direction, we see a stark difference manifesting itself in the behaviour of the central 'local' Yb CoM from the corresponding CoM behaviour of the full Yb system [ Fig. SM.2(a)(ii)]. This indicates the existence of opposing flows in Yb along the z direction, and can be thought of as the extension of co-moving and countermoving modes to our complicated system of different trap frequencies. The existence of opposing flows in Yb, in turn gradually couples to the Cs component, leading to a rather complicated excitation pattern after the first few oscillations, as evident from the left panels of Fig. SM.2(c)-(d). Such features can also be seen in the subsequent Fig. SM.3 revealing co-existence of dynamical components for both Cs and Yb across the same dominant frequencies.
Our findings can be reinforced by a detailed analysis of the changes ∆E in energy exhibited by the various single-particle energy contributions along the different directions [ Fig. SM.2(f)-(i)]. During the period ≈ (2 − 4) ms, when Cs and Yb are primarily moving in opposing directions, we see the gradual establishment of an additional important mechanism of internal energy transfer within the Yb component, specifically from its directly excited (by the moving Cs component) z direction, to the other transversal y direction of the elongated cloud (ω Yb,y = 1.5ω Yb,z ). This feature is evident in Fig. SM.2(h)(ii). Such irregularly-oscillating energy ex-change across the two tight Yb directions (visible already in Fig. 2(c) of the main paper) saturates after ≈ 20 ms, signalling the end of the primary phase of rapid evolution.
During this period, excitations along the other directions are present [see Fig. SM.2(i)], but suppressed in comparison, with only modest energy transfer on a much longer timescale (due to the smaller natural frequencies). This explains why we did not consider energy transfer along such channels as important in the main text.
The above energy calculations across the three directions x, y, and z (collectively labelled by ν) are based on the single-particle Hamiltonian contributions wherê includes both kinetic energy (KE) and potential energy (V) contributions, with intra-and inter-component interaction energies respectively evaluated by and The abundant excitations in the system begin to gradually thermalize, causing further damping, minimizing the subsequent beat amplitudes, the origin of which is explained in more detail in Sec. I B 2 below.

Excitation Spectrum and Origin of Beating
The frequency shift of dynamical collective modes of coupled superfluids due to the cross interaction between components has been studied in the context of both Bose-Bose and Bose-Fermi superfluid mixtures [8,31,[46][47][48]. In the Cs-Yb mixture under consideration here, the Cs particle number is much smaller than the Yb number. We therefore may neglect the back action of the Cs on the Yb along the z axis. As a consequence, the Yb component oscillates with the harmonic trap frequency ω Yb,z . Moreover, within the Thomas-Fermi approximation, the Yb density profile, given by n Yb (z) = [µ Yb − V Yb,z (z)] /g Yb , can be used to estimate the effective potential acting on Cs due to interspecies interactions. The effective Cs potential along the z axis is then V eff Cs,z (z) ≈ V Cs,z + (4π 2 a CsYb /M µ )n Yb (z). This potential has an effective harmonic frequency of This simple formula predicts ω eff Cs,z ≈ 1.05ω Cs,z , which captures the effective frequency for the Cs in the presence of a static Yb cloud in the following discussion.
The previous subsection studied the dominant dynamical regime of the first ≈ 22 ms, where most of the Cs excess energy is transferred to Yb through the generation of excitations due to their relative motion. Detailed numerical consideration of the heavily damped coupled oscillations indicates both the existence of multiple dynamical modes, and the emergence of beating on longer timescales. Here we shed more light on both these features, and the dominant oscillation frequencies of the coupled mixture, by analyzing the mixture's excitation spectrum, through the Fourier spectra of z CoM Cs(Yb) . These are shown by the red (Cs) and blue (Yb) solid lines in Fig. SM.3.
We observe a multitude of excitation frequencies in both components, consistent with the complicated excitation generation discussed above. As expected, the dominant dynamics for each component arise at frequencies close to their natural frequency, slightly upshifted due to the mutual attractive mean-field potential, as discussed above. Such peaks, and the general underlying structure for each component, reflect well the excitation spectrum that would be expected when considering the modification imposed by the attractive mean-field potential, but ignoring the dynamical back action of the other component. To better characterize such an effect of the attractive mean-field potential, we have also performed corresponding dynamical numerical simulations for Cs (and Yb) when coupled to the static other component, i.e. respectively Yb (and Cs). This has been done by adding to our simulations a static mean-field contribution from the second component, set to its configuration at t = 0 in our simulations, whose role is to change the  Beyond the dominant contributions for each component discussed above, the dynamical excitation frequencies reveal a multitude of other relevant modes. Firstly, we find each component to be additionally responding (in their central overlap region) to the mutual attractive mean field by co-moving at the dominant frequency of the other superfluid. In other words, parts of the (centrallylocated) Yb are dragged along by Cs at ≈ ω Cs,z , while correspondingly there is a notable component of Cs being dragged along by Yb at ≈ ω Yb,z . We also find a clear (broad) peak for both components around their frequency difference ≈ (ω Cs,z − ω Yb,z ) ≈ 0.7ω Cs,z [which is also close to the average frequency (ω Cs,z +ω Yb,z )/2]. The simultaneous presence of dominant oscillation modes for the two components both close to the natural frequencies of each component, and close to the difference (and average value) of such frequencies, indicates the strong presence of co-moving and counter-moving oscillating components. Here the dynamics are complicated by the different masses, atom numbers, interaction strengths and trapping potentials for the two components.
Beyond such modes, our spectra clearly demonstrate the emergence of splitting of the dominant upshifted frequencies of each component into two (or multiple) closeby frequency peaks. Such frequencies are generally found to be upshifted from the natural frequency, and to lie on either side of the static-mean-field contribution, thus highlighting the dynamical role of the other component. Specifically for Cs, we find two dominant frequencies, of broadly comparable importance. For the 'primary' atom number combination considered here, these are ω (1) Cs,z ∼ 1.063ω Cs,z and ω (2) Cs,z ∼ 1.005ω Cs,z . The difference in these two frequencies establishes a timescale 2π/ ω (1) Cs,z − ω  Fig. 2(c) of the main text. The period of this beating cycle is highly sensitive to the particular atom numbers considered. Additional features seen within each beating cycle can be attributed to the intricate interplay between the various relevant modes discussed above, which will in general have different relative phases and damping rates, and lead to dynamics on a faster, few ms, timescale.

C. Density Snapshots
To better visualize the complicated dynamics, we present in Fig. SM.4 representative 3D density isosurfaces and 2D sliced density snapshots during the first ≈ 22 ms of coupled evolution, showing also the corresponding 1D Shown are (i) 3D profiles corresponding to density isosurfaces for the 5% peak density of Yb at t = 0 with corresponding slices shown along each axis and 2D sliced density snapshots across the (ii) x − y (z = 0) plane, (iii) x − z (y = 0) plane and (iv) y − z (x = 0) plane.
slices. The full evolution, which also includes later times can be found in the accompanying Supplementary Video.
From these, we can see the different stages of the evolution, including the Yb component being dragged by Cs over the first ≈ 2 ms, the subsequent double-peak distribution of Yb along the z axis, and the gradual excitations along different directions.

II. COLLAPSE PROPERTIES AND DYNAMICS
We model the experiment using a three-dimensional mean-field description, expressed mathematically as twocomponent coupled GPEs with an added three-body loss term [49] for Cs as detailed in Eqs. 1-5. Here we have ∆z = 0, and ω z = 2π × 260 Hz. We take a Yb = 105 a 0 and a CsYb = −75 a 0 . Our numerical simulations of these coupled GPEs, described below, use Fourier pseudospectral methods with adaptive Runge-Kutta timestepping and are implemented using the XMDS2 software [50].

A. Ground states and static collapse analysis
One can attempt to find ground state solutions [∂ t ψ Yb = ∂ t ψ Cs = 0] to the coupled GPEs with a fixed Cs scattering length and zero three-body loss K Cs , collapse occurs [51] and mean-field ground state solutions cease to exist. We use imaginary time propagation to find ground state solutions numerically. For Cs only, with N Cs = 4 × 10 3 , we find a (crit) Cs ≈ −2.15 a 0 . We found that an approximate variational calculation [52,53] assuming a Gaussian ground state profile (and ignoring rotation of trap axes) yields a fairly similar result (a (crit) Cs ≈ −2.48 a 0 ) without the computational effort of solving the GPE. With Yb also present, with atom numbers N Cs = 4 × 10 3 and N Yb = 5 × 10 4 , we find a (crit) Cs ≈ 51.25 a 0 using imaginary time propagation. In this case, a Gaussian variational calculation gives incorrect results (a (crit) Cs ≈ −1.75 a 0 ) since the Yb condensate is far from having a Gaussian profile.

B. Dynamical collapse analysis
To model the dynamic collapse experiments reported in the main text, we perform real-time simulations of the coupled GPEs that are intended to match the experimental protocol. We begin from the Cs+Yb ground state solution, found as described above, with a Cs (t = 0) = 147 a 0 , N Cs = 4 × 10 3 , and N Yb = 5 × 10 4 . When simulating the Cs-only dynamic collapse experiments we continue to find the Cs+Yb ground state, but we remove Yb atoms immediately at time t = 0 by setting ψ Yb = 0, in analogy to the experimental protocol.
We then perform a real time simulation of the coupled GPEs in which the Cs three-body loss coefficient is set to a constant, non-zero, value and the Cs scattering length is varied in time to match the experimental protocol. Specifically, we implement in sequence a hold at a Cs (t) = 147 a 0 for 5 ms, a 10 ms linear ramp to a Cs (t) = a (jump) Cs , a hold at this value for 30 ms, a 10 ms ramp back to a Cs (t) = 147 a 0 , and a final 1 ms hold at this value. At this time in the experimental protocol a timeof-flight measurement begins; in the simulations we simply read out the final atom numbers. We repeat this process for varying values of a (jump) Cs , obtaining the results reported in the main text. All parameters are set to experimentally measured values rather than fitted to the data, except for the Cs three-body loss coefficient. The latter is less precisely constrained by experimental measurements. In practice we repeated the simulations for a few trial values and found that K 3,Cs = 1.8 × 10 −27 cm 6 s −1 gives a reasonable quantitative match to the results of the collapse experiment, while being consistent with other experimental measurements [26,42]. We have not searched for a precise best-fit value of K (T) 3,Cs (e.g., by least-squares fitting) because of the long computational times needed for the simulations.
By default our data points are from simulations in a box of dimension (L x , L y , L z ) = (96, 24, 24) µm with grid resolution of 3/2 points per µm in each direction. We checked convergence of the final N Cs value to ≈ 1% by repeating selected simulations at higher resolution (3 points per µm) and in larger boxes in the most-constrained direction (L y = 48 µm). Our data points with the lowest a (jump) Cs value in each simulation were obtained from higher-resolution simulations, for which the final N Cs value agreed to within ≈ 1% with the larger-box simulations. For values of a (jump) Cs −0.7 a 0 (Cs-only), or 51 a 0 (Cs+Yb), the three-body loss coefficient is not sufficient to prevent a mean-field collapse, even in our higher-resolution simulations. When this occurs the final Cs number is not converged and changes with the grid resolution. Quantitative modeling of this collapse regime would require a beyond-mean-field description, which is beyond the scope of this work.