Pump-probe spectroscopy of Bose polarons: Dynamical formation and coherence

We propose and investigate a pump-probe spectroscopy scheme to unveil the time-resolved dynamics of fermionic or bosonic impurities immersed in a harmonically trapped Bose-Einstein condensate. In this scheme a pump pulse initially transfers the impurities from a noninteracting to a resonantly interacting spin state and, after a ﬁnite time in which the system evolves freely, the probe pulse reverses this transition. This directly allows us to monitor the nonequilibrium dynamics of the impurities as the dynamical formation of coherent attractive or repulsive Bose polarons and signatures of their induced interactions are imprinted in the probe spectra. We show that for interspecies repulsions exceeding the intraspecies ones a temporal orthogonality catastrophe occurs, followed by enhanced energy redistribution processes, independently of the impurity’s ﬂavor. This phenomenon takes place over the characteristic trap timescales. For much longer timescales a steady state is reached characterized by substantial losses of coherence of the impurities. This steady state is related to eigenstate thermalization and it is demonstrated to be independent of the system’s characteristics.


I. INTRODUCTION
Time-resolved spectroscopy is an established technique for the characterization of the dynamical response of a wide range of physical systems [1].The general idea underlying a pump-probe spectroscopy (PPS) scheme is that a pump pulse prepares a nonstationary state of the system under consideration, which is then interrogated by a time-delayed probe pulse.This allows for simultaneous spectral and temporal resolution of the induced dynamical processes, exposing the energy redistribution of the selectively triggered excitations [2,3], in sharp contrast to time-independent spectroscopic techniques like injection spectroscopy [4][5][6][7].Applications of the PPS protocol range from two-and three-level atomic systems [8][9][10][11][12][13][14] to the ultrafast dynamics of photoexcited quantum materials [15][16][17][18][19][20][21][22].Such a time-domain analysis has been proven to be a powerful tool for resolving the ultrafast molecular dynamics allowing, for instance, for a coherent control of bound excited-state dimers over long timescales [23,24].PPS has also been utilized for studying the paircorrelation dynamics of ultracold Bose gases [25], offering a potential connection between ultrafast and ultracold physics [24,26].
Operating in the ultracold regime, in this work we propose a PPS scheme as a toolkit for investigating in a time-resolved manner the impurity problem and the related formation of Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. and interactions between quasiparticles .Understanding the physics of quasiparticles is important beyond cold atom settings in semiconducting [65] and superconducting devices [66].Additionally, interactions among quasiparticles in liquid helium mixtures [67,68] and cuprates [69,70] are considered to be responsible for conventional and high-T c superconductivity [71][72][73][74][75][76][77].Here we consider a Bose-Einstein condensate (BEC) with one or two impurities of either bosonic or fermionic nature immersed into it and track the emergent Bose polaron formation [42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][61] with a PPS radiofrequency protocol analogous to the one used in the experiment of Ref. [78].This allows us to probe and control the coherence properties of the quasiparticles.Our results pave the way for transferring the knowledge regarding the ultrafast dynamics of condensed matter systems [79][80][81][82] to the ultracold atomic realm.
In our investigation, an intense pump pulse transfers the initially free bosonic or fermionic impurities to an attractively or repulsively interacting state with the environment.After a variable dark time, during which the system evolves freely, a probe pulse of weaker intensity is applied, which deexcites the impurities.As the formation of well-defined attractive and repulsive Bose polarons in this many-body (MB) system is captured in the probe spectrum, this process allows us to monitor the dynamics.In systems where the interaction strength between the impurity and the background is not larger than the interaction strength within the background gas, polaronic excitations can have long lifetimes.However, beyond that limit substantial losses of coherence occur with a temporal orthogonality catastrophe (TOC) [59][60][61]83] being imprinted in the probe spectrum.The TOC emerges due to the relaxation of the quasiparticles into energetically lower-lying, phaseseparated states.This process is independent of the number of the impurities or their statistics.Remarkably, for timescales longer than the characteristic confinement one, the probe spectrum unveils evidence toward eigenstate thermalization [84][85][86][87], where the impurities reside in an incoherent state characterized by a large effective temperature.This relaxation dynamics [88,89] is found to be independent of the size of the bath, the number and nature of the impurities, and their interaction strengths and mass.
Our work is structured as follows.Section II introduces the setup under consideration and briefly comments on the employed variational approach to tackle the nonequilibrium dynamics of Bose polarons.In Sec.III we discuss the utilized PPS scheme and demonstrate the resulting Bose polaron spectrum for short and long evolution times with a particular focus on the impurity-impurity-induced interactions, coherence properties, and thermalization processes.In Sec.IV we elaborate on the emergent energy redistribution processes, while in order to gain further insights into the spectroscopically observed relaxation dynamics we invoke in Sec.V the eigenstate thermalization hypothesis (ETH).We summarize our results and provide an outlook including future perspectives in Sec.VI.Appendix A presents in detail the used radiofrequency spectroscopy scheme and Appendix B explicates briefly the predictions of a Ramsey protocol for strong impurity-medium interactions.The dimensional reduction of our MB Hamiltonian from three to one dimension (3D to 1D) is showcased in Appendix C. Finally, Appendix D deals with the variational method employed herein so as to simulate the PPS protocol and Appendix E delineates the convergence of the presented results.

II. MODEL SETUP
Our model is a highly particle imbalanced mixture.It consists of N I = 1, 2 bosonic or fermionic impurities (I) having a spin-1/2 degree of freedom [90] being immersed in a bosonic bath of N B = 100 structureless bosons (B).The mixture is assumed to be mass balanced, m B = m I ≡ m (unless stated otherwise), while both species are harmonically confined in the same one-dimensional potential.Details of the dimensional reduction of our system are discussed in Appendix C. The MB Hamiltonian reads denote the noninteracting Hamiltonian of the BEC and the impurities, respectively, while a ∈ {↑, ↓}.Additionally, ˆ B (x) [ ˆ a (x)] is the field operator of the BEC (spin-a impurities).We further consider that the dominant interaction is an s-wave one since we operate in the ultracold regime.As such both intra-(g BB , g II ) and interspecies (g BI ) interactions are adequately described by a contact potential [91], see also Appendix C. Furthermore, , with a, a ∈ {↑, ↓}, correspond to the contact intraspecies interaction terms of the bosonic bath and the impurities, respectively.Note that only the spin-↑ component of the impurities interacts with the BEC while the spin-↓ one is noninteracting.The relevant interspecies interaction term reads ĤBI Finally, and β = ν β − ν 0 referring to the bare Rabi frequency and the detuning of the radiofrequency pulse when the bosonic bath is absent, see Appendix A for further details.Here β ∈ {pump, probe, dark}.Moreover, the total spin operators are given by Ŝi = dx ab ˆ a (x)σ i ab ˆ b (x), with σ i ab denoting the Pauli matrix i ∈ {x, y, z}.
It is worth mentioning at this point that the onedimensional description adopted holds under the conditions k B T hω h2 m [ρ (1)  B (x = 0)] 2 ≈ 3 4/3 16 ( 1 [92,93].In these expressions, a BB is the three-dimensional s-wave scattering length between the particles of the medium, and α = √ h/(mω) and α ⊥ = √ h/(mω ⊥ ) correspond to the axial and transversal length scales.ρ (1)  B (x = 0) is the initial one-body density of the environment at x = 0, k B is the Boltzmann constant, and T refers to the temperature of the bosonic bath.To provide a concrete example, assuming ω ≈ 2π × 100 Hz and considering a 87 Rb gas with g BB = 0.5 ( h3 ω)/(m) ≈ 3.55 × 10 −38 Jm our 1D setting can be realized for transverse frequencies ω ⊥ ≈ 2π × 5.1 kHz.Accordingly, the 1D treatment is valid since N B a BB α ⊥ /α 2 ≈ 0.07 1 and temperature effects are negligible for k B T 316 hω ≈ 1.5 μK.
To access the time-resolved spectral response of bosonic and fermionic impurities immersed in the BEC bath the multilayer multiconfiguration time-dependent Hartree method for atomic mixtures is utilized [94][95][96].The latter is a nonperturbative approach that uses a variationally optimized timedependent basis which spans the optimal subspace of the Hilbert space at each time instant and allows for tackling all interatomic correlations [59].In particular, the MB wave function is expressed as a truncated Schmidt decomposition using D species functions for each component [Eq.(D1) in Appendix D].Next, each of these species functions is expanded in a basis of d B and d I single-particle functions for the BEC background and the impurities, respectively [Eq.(D2)].These single-particle functions utilize a time-independent primitive basis that is a tensor product of basis states regarding the spatial and the spin degrees of freedom [Eq.(D3)].Then, by following a variational principle, we arrive at a set of coupled nonlinear integrodifferential equations of motion [94][95][96].A detailed description of our MB variational approach and the ingredients of our numerical simulations are provided in Appendices D and E, respectively.

III. PUMP-PROBE SPECTROSCOPY SCHEME
We prepare the multicomponent system in its ground state with fixed g BB and g II = 0.The impurities are in their spin-↓ state and thus ĤBI = 0. To trigger the dynamics, an intense, pump

= 10ω
ω, rectangular pump pulse drives the noninteracting spin-↓ impurities to their interacting with the bath spin-↑ state for −t e < t < 0 (where t e denotes the exposure time) [Fig.1(a)].The condition pump R0 ω ensures that the duration of the pump pulse is much smaller than the time interval in which the polarons form.Accordingly, the polaron formation can only occur after the termination of the pump pulse and therefore it can be captured by the subsequent probe pulse.To ensure the resonance condition of the pump pulse, namely pump = + , and to optimize t e = π/ pump R , the fraction of impurity atoms that have been successfully transferred to the spin-↑ state, N ↑ (t = 0) /N I , is monitored for variable pump .The resulting pump spectrum features a coherent atomic resonance [34,36,39,57,78] at pump = + .The latter, for N I = 1 and g BI = ±0.5, 1.5 h3 ω/m, is clearly visible in Fig. 1(b).Notice also that secondary peaks possessing an intensity of the order of 12% with respect to the dominant ones also emerge due to the rectangular shape of the pump pulse (see also Appendix A).
After the initial pump sequence the remaining population of the spin-↓ state is annihilated by employing an optical blast that projects the impurities to the |↑ state (as described in Appendix A) and subsequently the spin-↑ atoms are left to evolve for fixed g BI and dark R0 = 0 but variable dark time t d .The polaronic states can form within 0 t t d while at t = t d a probe pulse is applied.This pulse is characterized by probe R0 = ω pump R0 so as to enhance the spectral resolution of the signal obtained by the fraction of impurity atoms transferred to the spin-↓ state, N ↓ (t d ) /N I for variable probe .For the same reason the duration of the probe pulse is fixed to t e = π/ probe R , where probe R is the resonant Rabi frequency of the probe pulse at probe = + , t d = 0, and N I = 1.
Concluding within the PPS scheme, polaronic states can be identified in the probe spectrum as well-defined peaks with amplitude N↓ (t d ) /N I < 1.For our purposes (accounting for the finite fidelity resulting after the probe pulse) we employ the criterion N↓ (t d ) /N I < 0.96 in order to identify the polaronic resonances.Interestingly, a peak with N↓ (t d ) /N I ≈ 1 does not correspond to a polaron as it implies that the accessed MB state is equivalent to a noninteracting state.Accordingly, the peaks exactly at t d = 0 and N↓ (t d ) /N I ≈ 1, that will appear later, do not indicate the formation of polarons.No-tice that polaronic peaks with N↓ (t d ) /N I < 1 can occur for strong impurity-BEC interactions g BI > g BB , even for t d = 0, demonstrating fast energy transfer to the polaronic states for t d < ( probe R0 ) −1 .

A. Short-time dynamics of Bose polarons
The short-time (t ∼ ω −1 ) dynamics of few, N I = 1, 2, fermionic or bosonic impurities with g II = 0 immersed in a BEC bath of N B = 100 atoms is captured by the probe spectra for distinct attractive [Figs.1(c)-1(e)] and repulsive [Figs.1(f)-1(k)] interspecies interactions g BI .Focusing on the attractive side, a well-defined polaron at t d = 0 with central peak location at + = −8.7ω[Fig.1(c)], +,B = −8.9ω[Fig.1(d)], and +,F = −8.6ω[Fig.1(e)] is identified in the cases of N I = 1 and N I = 2 bosonic and fermionic impurities respectively.These polarons show a nonsizeable shift for all the different evolution times t d as long as N I = 1.However, a clear shift can be inferred for N I = 2 [Fig.1(d)].This shift, being of about 10%, is a consequence of the energy redistribution between the bosonic impurities and the BEC as demonstrated in Ref. [57] and it is further related to the fact that for g BI < 0 attractive induced interactions are significantly enhanced [57].Additionally, due to the pronounced induced interactions a collisional broadening [97] of the spectral line is clearly observed for t d = 1, 11ω −1 .Indeed, since the two-body state of the impurities evolves rapidly during the probe sequence the spectral resolution of the measurement decreases, giving rise to a wide background in the PPS spectrum for −10 < probe /ω < 2. The imprint of induced interactions in the spatiotemporal evolution of the one-body density is the dephasing of the breathing oscillations [hardly visible in the inset of Fig. 1(d)] within the time interval 10 < ωt d < 15, which is absent for the single impurity [see the inset in Fig. 1(c)].In contrast to the above dynamics, the time-resolved evolution of fermionic impurities closely resembles the single impurity one with the two fermions undergoing at short times a coherent breathing motion, as is evident in the inset of Fig. 1(e).This latter result can be easily understood by the fact that attractive induced-interactions between fermionic impurities are known to be suppressed, providing in turn a nonsizeable shift of the respective atomic peak resonance captured by the probe spectra [36,39].
Switching to repulsive interactions, the dynamical evolution of the system changes dramatically.Independently of flavor and concentration the motion of the impurities, as detected by the one-body density evolution for g BI = 0.5 h3 ω/m, is apparently qualitatively similar [insets in Figs.1(f From the very early stages of the nonequilibrium dynamics the density filamentizes with recurrences of an almost central density peak occurring at the collision points, i.e., around t d ≈ 16ω −1 and t d ≈ 19ω −1 for the bosonic and fermionic impurities, respectively [Fig.1(g) and Fig. 1(h)].However, in all three cases a clean quasiparticle peak is monitored in the respective probe spectra indicating the existence of welldefined polarons for t ∼ ω −1 .The dominant peak location appears to be shifted for t d = 0 when compared to + , while an overall broadening of the spectrum is also inferred.The observed shift depends strongly on the impurity's nature and it is measured to be of about 10% for bosonic but dropping down to almost 5% for fermionic impurities [Figs.1(g) and 1(h)].Interestingly, the difference in the spectrum between the N I = 1 and the N I = 2 bosonic case [compare Figs.1(f) and 1(g)] is negligible, suggesting that attractive induced interactions cannot be directly unveiled by the observed shift.A result that complements earlier predictions indicating that attractive induced interactions are suppressed in the repulsive regime [60,98].Indeed, for g BI > 0 the density of the BEC is less distorted compared to g BI < 0 [60,61].Accordingly, the impurities are less attracted to these distortions and consequently to each other.In turn, by inspecting the oscillatory tail of the probe spectra for g BI > 0 interference phenomena associated with the filamentation process can be identified.Indeed, already from the single impurity [Fig.1(f)] the amplitude, A( probe ), of the secondary peak appearing in the spectra, e.g., at t d = 8ω −1 , A( probe ≈ 4.6ω) = 0.352 is comparable with the dominant one A( probe ≈ 7.5ω) = 0.755.Notice that the intensity ratio of the secondary to the dominant peak is larger that 12% and thus cannot be attributed to the rectangular shape of the probe pulse.The latter, directly reflects the coherence between the filaments formed in the one-body density (see Discussion).However, as the number of impurities increases significant losses of coherence take place.Indeed, the secondary peak at the two bosonic impurities [Fig.1(j)] and to a slightly shifted but significantly broadened peak at +,F = 25.6ω for the fermionic ones [Fig.1(k)].The latter broadening is attributed to the fermion statistics.Indeed, fermionic impurities occupy higher momenta and as such couple stronger to the BEC excitations.However, at t d = 2ω deformation of the central peak is present and a highly oscillatory tail is seen in all cases.To appreciate the aforementioned degree of coherence already indicated by the probe spectra we next invoke the spatial first-order coherence function, Here t ) is the MB wave function, and ρ (1)  σ (x; t ) is the one-body density, see also Appendix D. Importantly, |g (1)  σ (x, x ; t )| ∈ [0, 1] indicates the spatially resolved deviation of a MB wave function from a corresponding product state.Specifically, if |g (1)  σ (x, x ; t )| = 1, then the system is termed fully coherent;, otherwise, coherence losses occur, signifying the buildup of correlations [60,99].Indeed, the instantaneous |g (1)  ↑ (x, x ; t d = 2ω)| for g BI = 1.5 h3 ω/m clearly dictates that the quasiparticle remains adequately coherent since, e.g., |g (1)  ↑ (x = −5 2(a)].Finally, notice that for t d > 6ω any quasiparticle notion is lost as detected by the probe spectra.This outcome, being consistent with recent works [39,57], is also supported by the diffusive behavior of the corresponding one-body density evolution [insets in Figs.1(i cles probe = 0 emanates in the PPS spectrum, referring to a phase separation between the impurity and the BEC.It is also worth mentioning that by employing a corresponding Ramsey scheme (see the discussion in Appendix B) it is not possible to conclude the emergence of the TOC within the same time interval since the structure factor is still finite and drops close to zero for substantially longer evolution times.

B. Long-time Bose polaron dynamics
Next let us study the evolution of the system at longer times, 100 < ωt d < 300.Note that for a typical axial confinement ω ≈ 2π × 100 Hz, the interval 100 < t d < 300 corresponds to 0.16 < t d < 0.48 s.As time evolves one expects that significant losses of coherence signaling the buildup of correlations will take place in the MB evolution of the system [Fig.2(b)].A powerful asset of exposing the latter is the temporal average which depends solely on the eigenstate properties of the interacting system [100].This allows us to infer the relaxation tendency of the impurities in the framework of the ETH [84,86], see also the discussion in Sec.V. Evidently, ḡ(1) ↑ (x, x ) reveals that for g BI = 0.5 h3 ω/m the impurity is largely coherent [Fig.2(c)] while at g BI = 1.5 h3 ω/m any coherence property is lost [Fig.2(d)].This outcome is further supported by the time-resolved probe spectra illustrated for longer times in Figs.3(a)-3(i).A strong suppression of the interaction shift with respect to + is found to persist at long times, which, together with the weakly fluctuating amplitude A( probe ≈ + ) observed in the course of time, verifies the longevity of coherent single and two polarons irrespectively of their flavor and for both attractive and moderate repulsive g BI = ±0.5 h3 ω/m [Figs.3(a)-3(f)].Alterations come into play for g BI > g BB , where, as per our discussion above, losses of coherence, as captured by g (1)  ↑ (x, x ; t d ), become significant and the polaron picture breaks down.Here a two-humped distribution appears in our probe spectra independently of the number of the impurities and their nature.The most pronounced feature of Figs.3(g)-3(i) is the peak located at probe = 0.The latter findings suggest that a relaxed state is reached characterized by incoherent impurities being unpredicted so far.
Indeed, by fitting ḡ(1) ↑ (x, x ) to the corresponding prediction of the N I -particle Gibbs ensemble we obtain large effective temperatures.These refer to k B T eff = 8.45 hω for N I = 1 and k B T eff = 8.58 hω (k B T eff = 5.89 hω) in the case of two bosons (fermions) showcasing their tendency to approach an incoherent thermalized state, see also our detailed discussion in Secs.V and IV.Notice that the initial state of fermions involves higher momenta than bosons, while the critical velocity of the BEC is the same [101][102][103].Therefore, fermions couple stronger to the BEC excitations losing a larger portion of their energy, implying a smaller T eff .Further evidences supporting the observed thermalization [88,89] are provided by the temporal evolution of the variance, with S denoting the relevant spatial region in which the impurities reside and ḡ(1) ↑ ∈ (0, 1).Remarkably a tendency toward thermalization is seen [Fig. 2

(e)], with ḡ(1)
↑ saturating at long times irrespectively of the size of the BEC cloud and whether one or two, noninteracting or interacting, impurities are present and what their nature is.

IV. IMPURITY-MEDIUM INTERACTION ENERGY
To further support the thermalization tendency of the multicomponent system for strong impurity-BEC interactions at longer times of the nonequilibrium dynamics, we next inspect the behavior of the interspecies interaction energy.The latter quantity is defined as ĤBI (t ) ≡ (t )| ĤBI | (t ) , where the operator of the interspecies interactions is ĤBI = . Also, ˆ σ (x) and ˆ † σ (x) denote the σ -species field operator that annihilates and creates, respectively, a σ -species particle at position x.
The time evolution of ĤBI (t ) /N I is illustrated in Fig. 4 upon considering a pumping that drives the atoms to the spin-↑ state with g BI = 1.5 h3 ω/m.Specifically, we consider different settings consisting of N I = 1, N I = 2 spin-polarized bosons or fermions as well as a few-body system containing N B = 10 bosons and N I = 2 noninteracting impurities.In all cases we observe that ĤBI (t ) /N I decreases up to t = 100ω −1 while for later times, and in particular for t > 200ω −1 , it shows a saturation trend to a certain value depending on both N I and N B .Recall that this saturation effect allowed for the derivation of Eq. ( 11) within the ETH scheme.Notice also that the saturation value of ĤBI (t ) /N I is smaller for the few-body system (N B = 10, N I = 2) while for the N B = 100 setups ĤBI (t ) /N I acquires its smaller value for two fermionic impurities and takes almost the same value for N I = 1 and N I = 2 noninteracting bosonic impurities.Finally, let us note that the overall decreasing behavior of ĤBI (t ) /N I suggests a transfer of energy from the impurities to the bosonic gas as it has been also demonstrated in Refs.[59,61].This energy transfer process, identified by the decreasing rate FIG. 4. Evolution of the impurity-BEC interaction energy per impurity particle applying a pump pulse to drive the impurities to the spin-↑ state with g BI = 1.5 h3 ω/m for a single (N I = 1) and two (N I = 2) bosonic or fermionic impurities and also for a few-body bath consisting of N B = 10 particles (see legend).In all cases g BB = 0.5 h3 ω/m and g II = 0. of ĤBI (t ) /N I , seems to be enhanced for N B = 10 while for the N B = 100 case it is more pronounced for the fermionic impurities.We remark that a saturation trend at long timescales being in turn suggestive of the thermalization tendency of the system occurs also for other observables.These include, for instance, the dynamical structure factor (Appendix B) as well as entropic measures [104,105] such as the von-Neumann entropy [61,99,106] quantifying the degree of entanglement (results not shown for brevity).

V. CHARACTERIZATION OF THE RELAXATION DYNAMICS
We next explicate our method for characterizing the relaxation dynamics occurring in our setup during the hold time t d .To achieve this, we employ the ETH [84,86].Within this framework it is assumed that after a quench, the finite subsystems of a larger extended system relax to a steady state reminiscent of thermal equilibrium.Here, by fitting the time-averaged one-body density of the impurities to a thermal equilibrium one, we show that this thermalization process explains the relaxed state of the impurities emanating for long times after the orthogonality catastrophe of the polarons.
The relaxation of an isolated (closed) system is understood in terms of the principle of local equivalence [107].Within this framework as the thermodynamical limit is approached, i.e., the system size tends to infinity, the reduced density matrices of the involved few-particle subsystems at long times can be calculated in terms of the density matrix of a (generalized) Gibbs ensemble at thermal equilibrium.Indeed, if the only conserved quantity of the Hamiltonian is the total energy, then it can be shown that these reduced density matrices can be calculated in terms of the equilibrium density matrix within the Gibbs ensemble ρGibbs In Eq. ( 5) Z is the partition function stemming from the normalization of the density matrix, i.e., Tr[ ρGibbs ] = 1.Ĥ refers to the MB Hamiltonian and T eff , k B correspond to the effective temperature and the Boltzmann constant, respectively.Of course, our setup exhibits also other conserved quantities than the total energy.Below we resort to the approximation of no further symmetries as it is the only case that explicit results showing the relaxation dynamics of the system are available within ETH [107].As we shall show later on, the aforementioned choice leads to an excellent agreement between our numerical findings and the relevant estimates provided by applying Eq. (5).Within this approximation the effective temperature, T eff , is fixed by the conserved value of the energy per particle in the thermodynamic limit (TL), Here BEC → TL is defined as the limit where N B → ∞, g BB → 0, N B g BB = const, and g BI /g BB = const.Notice, however, that Eqs. ( 5) and ( 6), are impractical for calculations since the eigenvalues and eigenstates of the full interacting Hamiltonian, Ĥ , are required for the evaluation of the Gibbs ensemble of the (N B + N I ) MB ensemble which are difficult if not impossible to obtain.For this reason, we simplify the above-mentioned set of equations so as to obtain explicit results which can be subsequently compared with those obtained by the time evolution of the (N B + N I ) MB system within the multilayer multiconfiguration time-dependent Hartree (ML-MCTDHX) approach.
Since we intend to employ the thermodynamic limit where the MF Gross-Pitaevskii treatment of the BEC is exact in the weak interaction limit [92], it is reasonable to assume that the corresponding density operator of the Gibbs ensemble acquires the product form ρGibbs = ρ(N B ) B;Gibbs ⊗ ρ(N I ) ↑;Gibbs .Recall that during the dark time all of the impurities are in their spin-↑ state.In this case, the form of ρ↑;Gibbs is similar to Eq. ( 5), namely where H eff I is an effective Hamiltonian that acts only on the impurity.Equation ( 7) greatly simplifies the description of our system, since the density matrix of the impurity depends only on the eigenvectors and eigenvalues of an N I -body effective Hamiltonian.To proceed further we specify the form of Ĥeff ↑ .Within a zeroth-order approximation we assume that the BEC acts solely as a potential barrier for the impurities, and as consequence their effective Hamiltonian reads Notice that this approximation for the effective potential is a simplification of the impurity problem.First, it neglects, among others, the renormalization of the impurity's mass, m I → m eff I , due to the presence of the BEC [59].Second, the presence of induced interactions between the impurities cannot be captured [60].
The time dependence of the Hamiltonian of Eq. ( 8) implies a nonstationary state for the impurities.However, it is well known that following an interaction quench the density of the BEC is only slightly perturbed by the motion of the impurities [59][60][61].The latter justifies the substitution of the effective Hamiltonian by its time-averaged value Ĥeff B (x i ; t ) for the density of the bath.By incorporating the above-mentioned approximations we obtain explicit forms for the one-body density of the impurity within the Gibbs ensemble, namely ρ (1)  I;Gibbs (x, x ; Here n i (T eff ) is the distribution function of the N I particles and φ i refers to the eigenstates of Ĥeff ↑ .Due to the small number of impurities considered herein (N I = 1, 2) both the fermionic and the bosonic impurities do not follow the appropriate, for N I → ∞, Fermi-Dirac or Bose-Einstein distributions.Instead, it can be shown that the relevant distribution in the case of a single particle or two bosons is the Boltzmann distribution, with i being the eigenvalues of Ĥeff For two fermions the corresponding distribution reads where Note also that the one-body density of the impurity, Eq. ( 9), depends only on a single parameter, namely the effective temperature, T eff .
As shown earlier, the one-body density matrix of the impurities, ρ (1)  ↑ (x, x ; t ), saturates to its timeaveraged value, i.e., ρ (1)  ↑ (x, x ; t → ∞) ≈ ρ(1) (1)  ↑ (x, x ; t ), for long hold times, as it is evident in the relaxation dynamics of ḡ (1) [see also Fig. 2(e)].In order to facilitate the comparison of our results to the ETH prediction [Eq.( 9)] we fit the averaged one-body density matrix, ρ(1) ↑ (x, x ), obtained within ML-MCTDHX to the corresponding Gibbs ensemble, ρ (1)  ↑;Gibbs (x, x ; T eff ), and extract the value of T eff .Our results for the best-fitted parameters are shown in Fig. 5.We remark that the fitting is performed on the level of ρ(1) ↑ (x, x ), while only the diagonal ρ(1) ↑ (x) ≡ ρ(1) ↑ (x, x) is presented in Fig. 5 in order to enhance the visibility of the obtained results.By comparing the time-averaged one-body density and the fitted Gibbs ensemble prediction a very good agreement is observed for both one [Fig.5(a)] and two impurities of either bosonic [Fig.5(b)] or fermionic [Fig.5(c)] nature.This result holds equally also in the case of mass-imbalanced mixtures composed, for instance, of heavy bosonic impurities, m I = 133/78m B [see, e.g., Fig. 5(d)].The above findings indicate that despite the employed approximations the ETH scheme is able to capture the main features exhibited by the relaxed state of the impurities within our correlated MB system.Regarding the effective temperature we find rather large values of T eff for one and two bosonic impurities that is of the order of T eff ≈ 8 hω/k B .While for two fermions and two heavier bosonic impurities the temperature is slightly smaller, possessing values of the order of T eff ≈ 6 hω/k B .These large values of the effective temperature are indicative of the incoherent character of the impurities after the probe pulse [see also Fig. 2(d)].To advance further the correspondence between the ETH model and the correlated MB results we estimate the effective temperature of the relaxed state by expressing Eq. ( 6) only in terms of the impurity's degrees of freedom.Notice that the energy of the impurity is not conserved during the MB evolution of our system due to the presence of energy exchange processes between the impurity and the bath.However, as evidenced in Fig. 4, the energy of the impurities saturates for large times (see also Sec.IV).Indeed, by taking advantage of this observation we can cast Eq. ( 6) in the form where Ē↑ is the time-averaged impurity energy.In the case of a single impurity, Eq. ( 12) gives an estimation for the effective temperature of T eff = 8.56 hω/k B which is in good agreement with the effective temperature obtained by fitting T eff = 8.45 hω/k B .Note also here that the effective Hamiltonian of Eq. ( 8) is known to overestimate the zero-point energy of the impurity since it neglects its dressing by the excitations of the BEC [59].This in turn explains the higher T eff obtained via Eq.( 12).In contrast, in the case of two impurities the related estimates for T eff are much higher than the ones obtained by the fitting of ρ(1) ↑ (x, x ).Specifically, Eq. ( 12) yields T eff = 10.28 hω/k B and T eff = 6.89 hω/k B for the two bosonic and the two fermionic impurities, respec-tively.The observed discrepancy is attributed to the presence of induced interactions between the impurities that are more prevalent in the case of bosonic impurities than fermionic ones [60].However, their effect is neglected within the effective Hamiltonian of Eq. ( 8).

VI. CONCLUSIONS
We have developed a PPS scheme to study the timeresolved dynamics of fermionic and bosonic impurities immersed in a harmonically confined BEC.Coherence properties and induced interactions are encoded in the probe spectra for both attractive and repulsive interactions.Moreover, longlived attractive and repulsive polarons exist up to g BI ≈ g BB .For g BI > g BB , with the dynamics being dominated by energy redistribution processes, a rather rapid temporal orthogonality catastrophe occurs.To explicitly showcase that energy redistribution processes take place we have discussed the behavior of the corresponding interspecies interaction energy which decreases for short evolution times and thus captures the energy transfer from the impurities to the environment.Furthermore, it shows a saturation tendency for large evolution times, a behavior that is indicative of a relaxation process of the impurities.Indeed, at longer times (t d > 100ω −1 ), where any coherence information is lost, a thermalized state is reached.To further characterize the aforementioned relaxation dynamics we have resorted to an effective ETH model which has enabled us to identify that for strong interspecies couplings and long evolution times the impurities acquire an effective temperature.This effective temperature is found to be smaller for the fermionic impurities than the bosonic ones.Importantly, we have found that the thermalization process is independent of the size of the bath and the impurity concentration, the interacting nature of the impurities as well as their flavor and mass.
It would be intriguing to utilize PPS at finite temperature [108,109] and also in higher dimensions to explore the timeresolved formation of quasiparticles.As further perspectives, PPS could be exploited to unravel recondensation dynamics [110,111] in excited bands of optical lattices and the dynamics of vibrational states of ultra-long-range Rydberg molecules [112,113] to infer their lifetime.
Let us elaborate on the model that allows for the simulation of the MB dynamics under the influence of radiofrequency fields [5,36].This model has been employed in the main text for the characterization of the coherence properties of polaronic quasiparticles in the context of PPS.
In our case few atomic impurities are immersed in a BEC environment close to an interspecies magnetic Feshbach resonance [114].The case of bosonic impurities possessing equal mass to the BEC atoms can be realized by employing different hyperfine states of a particular isotope, e.g., 85 Rb or 87 Rb.For fermionic impurities the equal-mass scenario occurs approximately, e.g., for 173 Yb impurities immersed in a 174 Yb BEC with mass ratio of m B /m I ≈ 1.006.Different masses for the impurities and the BEC atoms can be realized by invoking different atomic species, e.g., considering 87 Rb and 133 Cs [115].Typically broad Feshbach resonances occur at magnetic fields of the order of 800G [114], referring to a regime where the atoms experience a sizable Zeeman shift.This sizable Zeeman shift allows us to address selectively the distinct m F transitions provided that the intensity of the radiofrequency pulse results in a Rabi frequency R much smaller than the Zeeman splitting of the involved hyperfine levels.The latter is typically of the order of a few tenths of MHz.This large splitting of the different m F levels implies that magnetic phenomena such as spin-exchange interactions can be safely neglected for these values of the magnetic field, see also Appendix C.
In this work we consider two hyperfine levels of the impurity atoms, denoted as |↑ and |↓ .These states can be identified and resonantly coupled for a frequency ν 0 , corresponding to the Zeeman splitting between the two levels, when a BEC environment is absent.Due to the harmonic confinement of the atoms each of the hyperfine levels is further divided into states of different atomic motion.The average spacing between these sublevels is of the order of the harmonic trap frequency, hω lying within the range of a few tenths of h × Hz to a few h × kHz in typical ultracold atom experiments [78,117].In the vicinity of a Feshbach resonance the energy of these sublevels strongly depends on the interspecies interaction strength, g BI , between the impurity atoms in the resonantly interacting hyperfine state and their BEC environment.Accordingly, the energy of each motional state shifts by + (g BI ) from the corresponding noninteracting one.As is also made obvious within the main text [see Fig. 1(b)] this shift is of the order of ω to ∼10ω.Therefore, due to the separation of the different involved energy scales it suffices to treat the impurities as two-level atoms.Furthermore, even for β R0 + ∼ kHz where the regime of strong intense pulses is accessed, β R0 ν 0 ∼ 10 MHz, allowing us to invoke the rotating wave approximation.Notice here that β ∈ {pump, probe, dark}.Within this approximation and in the interaction picture of the |↑ ↔ |↓ transition, the Hamiltonian for the internal state of the impurities reads 2 Ŝx , which is exactly the form employed in Eq. ( 1) of the main text.β R0 and β = ν β − ν β 0 refer respectively to the (bare) Rabi frequency and the detuning with respect to the resonance of the |↑ ↔ |↓ transition at g BI = 0. We remark that the |↑ and |↓ states in the Schrödinger and interaction pictures are equivalent, so our conclusions are invariant under this frame transformation [118].
To populate the polaronic states we employ a pump pulse of rectangular shape as depicted in Fig. 6(a).The system is initialized in the noninteracting ground state where the impurity atoms are spin polarized in their |↓ state.The pump pulse is characterized by frequency ν pump , and a detuning pump is employed.This pulse is further characterized by an exposure time t e and a bare Rabi-frequency pump R0 .Different realizations utilize different detunings pump and exposure times t e but the same pump R0 .In the duration of the pulse the system undergoes Rabi oscillations which for strong-enough pulses pump R0 ω are well characterized by a Rabi frequency R ( pump ) = ( pump R+ ) 2 + ( pump − + ) 2 , where pump R+ and + are the corresponding resonance values.By fitting the spectroscopic signal, which is the fraction of atoms transferred to the |↑ hyperfine state, to the theoretical line shape for rectangular pulses reading N↑ (t ) these resonance values of pump R+ and + can be obtained.Note here that the line shape [Eq.(A1)] exhibits an infinite sequence of peaks at the locations, pump n , n = 0, ±1, . . ., given by the solutions of R+ t e ) and A ±1 ≈ 0.01179( pump R+ t e ) 2 .In order to achieve a high spectroscopic signal, N↑ (t d ) /N I , we set the exposure time to t e = π/ pump R+ (to the obtained fitting accuracy) ensuring that A 0 ≈ 1.This choice implies that the peaks at pump ±1 are clearly imprinted in the obtained spectrum possessing an amplitude A ±1 ≈ 0.116438.Indeed, these side peaks can be clearly identified in Fig. 1(b).
To infer about the coherence properties of the polaronic states we employ PPS, see Fig. 6(a).Initially, we prepare the system in the same noninteracting ground state as in the previously examined protocol and apply a rectangular π pulse, with pump R0 = 10ω and t e = π/ pump R+ on a polaronic resonance where we have identified the resonant pump R+ and pump + as explained above.This sequence transfers the atoms from the ground state to the polaronic state in a very efficient manner, see Fig. 6(b).Then the impurity atoms are projected to the spin-↑ state by employing an optical burst transition on the lowest hyperfine state | ↓ to an available P electronic level at t = 0 which essentially ejects all the spin-↓ atoms from the trap.This procedure has been simulated by the application of the operator ĤP = −i h dx ˆ † ↓ (x) ˆ ↓ (x) over a short time interval t b .We can numerically verify that for large > 100ω and small t b ω −1 (corresponding to the experimentally relevant values) the action of ĤP to the state after the pump pulse is equivalent to the projection of the impurity to the spin-↑ configuration.For this reason and for computational simplicity we employ the state || P| (t=0 − ) || as an initial state for the subsequent time evolution t > 0. Note that this sequence for pump R ω is approximately equivalent to an interaction quench, as the pump-pulse has almost no spectral selectivity due to the pronounced power broadening of the radiofrequency transition, as pump R ∼ pump + and the fact that the out-of-equilibrium dynamics is effectively frozen due to the small timescale t e = π/ pump R ω −1 .Indeed, these properties of the pump-pulse have been verified numerically for the selected parameters pump R0 and pump as the MB state after this pulse is found to possess a fidelity in excess of 90% to the initial one.
After the pump sequence is completed we let the system evolve in the absence of radiofrequency fields, dark R0 = 0, for a dark time, t d .Finally, we apply a second probe π pulse with a smaller probe R0 = 1ω to the first one and varying probe + to transfer the atoms from the polaronic to the initial ground state.The employed spectroscopic signal is the fraction of atoms that have been deexcited by the probe pulse (recall that within our scheme all of the particles are at t d in the spin-↑ state) divided by the total number of impurities, N ↓ (t d ) N I .Note that a smaller probe R0 , when compared to pump R0 , is employed in order to reduce the power broadening during the probe sequence and subsequently increase the resolution in terms of detuning.However, this value cannot be arbitrarily lowered since for decreasing probe intensities the frequency resolution is increased at the expense of lower temporal resolution.For such low intensities the motional state of the spin-↑ impurities is significantly altered during the application of the probe pulse.As a heuristic argument the relation δνδt ≈ 1 that connects the temporal (δt) and spectral (δν) resolution is commonly employed [119].The value of probe R0 = 1ω is selected within this work as it consists an adequate trade-off between the spectral and the temporal resolution.Finally, due to the rectangular shape of the probe pulse the exhibited line shape of NI (t d ) N I is given by Eq. (A1) as long as the impurity is coherent, i.e., |g (1) (x, x ; t d )| ≈ 1. Accordingly, in our analysis we attribute all fringes appearing in the spectra to the line shape of a single resonance if the ratio of the amplitude of two neighboring peaks satisfies A n+1 A n < 0.12.

APPENDIX B: COMPARISON WITH RAMSEY SPECTROSCOPY
Next we demonstrate the advantage of utilizing the PPS scheme in comparison to Ramsey spectroscopy.In particular, we explicitly showcase the differences between the predictions of the PPS and the Ramsey schemes for intriguing phenomena exhibited by our system including the TOC and the thermalization process.To achieve this comparison we have simulated the Ramsey response of our system following the scheme described in Refs.[59,60].The main facet of this Ramsey protocol is that by applying an intense radiofrequency pulse to the initially noninteracting with the bath spin-↓ impurities we transfer them into a superposition state |↑ +|↓ √ 2 , where the state |↑ interacts with the bosonic medium.Thus, the time-evolved MB wave function, e.g., of a single impurity is given by BI is the spatial part of the initial MB wave function (see also Appendix D), and E 0 refers to the corresponding eigenenergy.In this protocol the structure factor, |S(t )| = | 0 BI |e iE 0 t/ he −i ĤR t/ h| 0 BI |, of the system is monitored by inspecting the magnitude of the impuritys' spin | Ŝ(t ) |.
According to our discussion in Sec.III it becomes apparent that time-dependent phenomena such as the TOC and the consecutive thermalization of the impurities can be clearly tracked in an experimentally relevant fashion via employing the temporarily resolved PPS scheme.Indeed, Figs. 1(i)-1(k) in the main text reveal that TOC takes place already for t d = 2ω −1 , with the presence of quasifree impurities at probe ≈ 0 being also readily imprinted in the probe spectrum.On the contrary, the only information that Ramsey spectroscopy conveys is the value of the structure factor, |S(t )|.Importantly, it does not deliver any further insights about the physical origin of its decreasing tendency and thus the underlying physical processes, see in particular Figs.7(a) and 7(b).Indeed, for t < 10ω −1 |S(t )| oscillates having a minimum value of 0.4 when g BI = 1.5 h3 ω/m and thus does not provide any clear signature for the emergence of the TOC identified using PPS.Along the same lines for g BI = 1.5 h3 ω/m > g BB , namely after the TOC manifests itself, a thermalization tendency is clearly imprinted in the probe spectrum [see Figs.3(g)-3(i)] with a predominant peak appearing at probe ≈ 0. In other words, the dynamics after the decay of the strongly (g BI > g BB ) repulsive Bose polarons leads to a quasistationary state of the MB system with respect to the energy redistribution among the different dynamical modes, providing this way strong evidences toward a thermalized state.This mechanism cannot be even suggested by invoking the contrast com- puted within the above-discussed Ramsey scheme.Indeed, the Ramsey scheme only indicates the tendency toward thermalization due to the decreasing behavior of the structure factor which, however, fluctuates within the depicted time interval.

APPENDIX C: DIMENSIONAL REDUCTION OF THE MANY-BODY HAMILTONIAN FROM THREE TO ONE DIMENSIONS
We consider an ensemble of confined ultracold atoms in three different hyperfine states, denoted as B, ↑, and ↓.State B is occupied by bosonic bath particles and the ↑ and ↓ states constitute a pseudo-spin-1/2 subsystem.We assume that state B belongs to a different manifold of hyperfine states with respect to the total angular-momentum quantum number, F , than the other two pseudospin states.The system is optically trapped and therefore all of the above hyperfine states experience the same confinement.
The ab initio Hamiltonian of this multicomponent system reads Ĥ = Ĥ0 + ĤSD + ĤI .The spin-independent part Ĥ0 is given by where m is the mass of the chemical element and V 0 (r) refers to the confining potential.By imposing a homogeneous magnetic field, along the zdirection, the energy of the magnetic sublevels characterized by different m F shift due to the Zeeman effect.Accordingly, the state-dependent part of the Hamiltonian is expressed as where σ z αβ corresponds to the spin-z Pauli matrix and E σ is the energy of the atomic state of species σ ∈ {B, ↑, ↓}.The typical energy difference between hyperfine levels possessing different F is of the order of several h × GHz.In particular, for 87 Rb the hyperfine splitting between the two lowest hyperfine manifolds F = 1 and F = 2 is E F =2 − E F =1 ≈ h × 6.83 GHz [116].Additionally, the amplitude of the Zeeman energy shifts is also of the order of h× MHz/G.For instance, in 87 Rb this amplitude is of the order of ∼0.7 h× MHz/G [116].Furthermore, in the same species quadratic Zeeman shifts that lead to a nonequidistant distribution of magnetic sublevels possessing an amplitude of several h × MHz can be observed already for magnetic fields of the order of ∼10 G [116].Typical ultracold atom experiments involve interaction energies ranging from hundreds of h × Hz to a few h × kHz generating this way interaction energy shifts and spin-exchange processes characterized by energies of the same order of magnitude.Therefore, except for the case where the magnetic field applied is of the order of few Gauss, the spin-exchanging collisions are strongly suppressed [120].
Operating in the ultracold limit dominated by s-wave scattering [92] the interaction Hamiltonian can be expressed as [121] The scattering lengths a σ σ , with σ ∈ {B, ↑, ↓}, can be tuned via a Fano-Feshbach resonance between two distinct hyperfine levels [114].
In order to effectively reduce the dimensionality of the above system from 3D to 1D a strong confinement along the two perpendicular spatial directions is usually employed [93].Then the confining potential reads where ω ⊥ ω holds for the transverse and longitudinal trapping frequencies.Note that the potential of Eq. (C4) can be realized either by a single optical dipole trap [120] or by applying a deep two-dimensional optical lattice potential [42,122].To access the 1D regime, the frequency of the transverse confinement ω ⊥ has to be selected such that the excited states of the harmonic trap along the transverse directions (y, z) are not populated.The condition for a 1D BEC is well known [93] and reads N B a BB α ⊥ /α 2  1, where α ⊥ = √ h/mω ⊥ and α = √ h/mω.In the few atom case, referring to the impurity species, a sufficient condition for accessing the 1D limit is ω ⊥ Nω [123].Indeed, under this assumption it is known that even in the strong interaction limit the system behaves as a Tonks-Girardeau gas of hard-core bosons sharing some characteristics with a gas of free fermions of the same particle number [124,125].Properties of such fermionized 1D bosons have been probed experimentally in Refs.[126,127].
Accordingly, the corresponding 3D field operators can be expressed in terms of 1D ones as follows: By employing Eq. (C5) we can then evaluate straightforwardly the reduced 1D effective Hamiltonian for Ĥ0 + ĤSD which takes the form Here all terms contributing to the total energy shift for constant N I and N B are dropped while the Ŝz operator reads The dimensional reduction of ĤI is, however, more complicated.In particular, the phenomenon of the confinement induced resonance [91,128] occurs when α ⊥ = √ h/(mω ⊥ ) is comparable to a σ σ .This implies that the actual 1D coupling constant deviates from g MF σ σ = 2 h2 a σ σ ma 2 ⊥ [93], which is obtained by evaluating the integrals appearing in Eq. (C3) along the transverse (y, z) directions.Detailed theoretical and experimental investigations [91,128] reveal that the 1D coupling strength g σ σ possesses a simple analytical form, i.e., ) −1 , and the effective 1D interaction Hamiltonian simplifies to Note here that due to the double counting for intraspecies interaction terms in Eq. (C8), the parameter g BB appearing in the latter is two times larger than the corresponding one that is involved in Eq. ( 1).According to the above discussion, for the experimental implementation of the setup described in the main text the interaction parameters, g BB , g BI , and g II , used herein are related to the corresponding 3D scattering lengths as follows: Here gσσ = g σ σ / h3 ω/m refers to the dimensionless interaction strength and η = α/α ⊥ = √ ω ⊥ /ω is the aspect ratio.
Furthermore, BECs involving particle numbers of the order of N B ∼ 100 are already accsessible by current state-of-theart experimental settings, e.g., in optical lattice experiments [42,122].Finally, it is worth commenting that three-body recombination processes are highly suppressed for alkali ultracold atomic vapors, as the one considered herein.For instance, for a 87 Rb BEC and in the presence of three-body recombination a lifetime of 14.8 s has been reported [129].Note also that the rate of three-body recombination scales with the cube of the density and, as a consequence, this effect is negligible for the mesoscopic system under consideration which involves low densities.

APPENDIX D: THE MANY-BODY VARIATIONAL METHODOLOGY: ML-MCTDHX
To track the stationary properties and, most importantly, the MB quantum dynamics of the multicomponent system addressed in the main text, we resort to the ML-MCTDHX [94][95][96].It constitutes an ab initio variational method for solving the time-dependent MB Schrödinger equation of atomic mixtures possessing either bosonic [57,99,130] or fermionic [39,106,118,131] spinor components.A major advantage of this approach is the expansion of the total MB wave function with respect to a time-dependent and variationally optimized basis (see below).This allows us to capture all the relevant inter-and intraspecies correlations of a multicomponent system in an efficient manner at each time instant by utilizing a reduced number of basis states when compared to expansions relying on a time-independent basis.
The system considered in the main text consists of a bosonic bath (B) with N B = 100 atoms and either one (N I = 1) or two (N I = 2) impurity (I) atoms.Most importantly, the impurities being either bosons or fermions possess an internal pseudospin-1/2 degree of freedom [59,90].To account for interspecies correlations, the MB wave function | (t ) is expressed according to a truncated Schmidt decomposition [99,130,132] where the time-dependent Schmidt weights λ k (t ) are also known as the natural species populations of the kth species function and provide information about the degree of entanglement between the individual subsystems.For instance, if two different λ k (t ) are nonzero, then | (t ) is a linear superposition of two states and therefore the system is entangled [132,133] or interspecies correlated.On the other hand, in the case of λ 1 (t ) = 1, λ k>1 (t ) = 0, the wave function is a direct product of two states and the system is nonentangled.Next, in order to include intraspecies correlations into our MB wave function ansatz, each species function | σ k (t ) is further expanded on a time-dependent number-state basis set | n(t ) σ .Namely where A σ k; n (t ) denote the underlying time-dependent expansion coefficients.Moreover, each number state | n(t ) σ corresponds to a permanent for bosons or a determinant for fermions building on d σ time-dependent variationally optimized single-particle functions (SPFs), i.e., |φ σ l (t ) , with l = 1, 2, . . ., d σ , being characterized by occupation numbers n = (n 1 , . . ., n d σ ).Additionally, the SPFs are expanded with respect to a time-independent primitive basis.For the majority species, this primitive basis corresponds to an M dimensional discrete variable representation denoted in the following by {|q }.However, for the impurities the primitive basis refers to the tensor product {|q, s } of the discrete variable representation basis regarding the spatial degrees of freedom and the two-dimensional pseudospin-1/2 basis {|↑ , |↓ }.Consequently, each SPF of the impurities acquires the following Here B I jq↑ (t ) [B I jq↓ (t )] are the time-dependent expansion coefficients of the pseudospin-↑ and ↓ respectively, see also Refs.[59,118].
Having exemplified the MB wave-function ansatz and in order to address the time evolution of the (N B + N I )-body wave function | (t ) obeying the Hamiltonian of Eq. ( 1) provided in the main text we then numerically solve the so-called ML-MCTDHX equations of motion [94].These equations are determined by following the Dirac-Frenkel [134,135] variational principle for the generalized ansatz of Eqs.(D1), (D2), and (D3).In this way, we obtain a set of D 2 linear differential equations of motion for the λ k (t ) coefficients coupled to D( (N B +d B −1)! N B !(d B −1)! + (N I +d I −1)! N I !(dI −1)! ) nonlinear integrodifferential equations for the species functions and d B + d I nonlinear integrodifferential equations for the SPFs.

APPENDIX E: CONVERGENCE OF THE MANY-BODY SIMULATIONS
The Hilbert space truncation within the ML-MCTDHX method is determined by the considered orbital configuration space, i.e., C = (D; d B ; d I ).In this notation, D = D B = D I and d B and d I denote the number of species functions and SPFs respectively for each species [Eqs.(D1) and (D2)].Moreover, within our numerical calculations we employ a primitive basis based on a sine discrete variable representation for the spatial part of the SPFs with M = 600 grid points.This sine discrete variable representation intrinsically introduces hard-wall boundary conditions at both edges of the numerical grid which in our case are located at x ± = ±50 √ h/mω.We assured that the location of the hard-wall boundaries does not impact our findings since no significant density portion occurs beyond x ± = ±20 √ h/mω.The eigenstates of the multicomponent system are obtained by utilizing the so-called improved relaxation method [94][95][96] within ML-MCTDHX.To address the corresponding nonequilibrium dynamics, we numerically solve the ML-MCTDHX equations of motion using the MB wave function [Eq.(D1)] under the influence of the Hamiltonian (1) of the main text.
To testify the convergence of the MB results we ensured that all observables of interest are to a certain level of accuracy insensitive for a varying orbital configuration space, C = (D; d B ; d I ).Note that for the MB simulations discussed in the main text we relied on the orbital configuration C = (10; 3; 8).To infer the convergence of our results we exemplarily showcase below the behavior of the spatially integrated one-body coherence function, g (1) (x, x ; t ), for different number of species and SPFs in the course of time.In particular we calculate its normalized absolute deviation between the C = (10; 3; 8) and other orbital configurations C = (D; d B ; d I ), namely G C,C (t ) = dxdx g (1)  C (x, x ; t ) − g (1)  C (x, x ; t ) dxdx g (1)  C (x, x ; t ) . (E1) The dynamics of G C,C (t ) is illustrated in Fig. 8 for the multicomponent bosonic system consisting of N B = 100 atoms and N I = 2 noninteracting impurities upon considering the pump spectroscopic sequence introduced in Sec.A from g BI = 0 either to g BI = 0.5 h3 ω/m [Fig.8 at long evolution times t > 150ω −1 .It is also worth mentioning at this point that for all other observables and interspecies interaction strengths discussed in the main text a similar degree of convergence takes place (results not shown here for brevity).

5 FIG. 1 .
FIG. 1.(a) Schematic illustration of the PPS pulse sequences used.(b) Spectral response of the pump pulse N↑ (t = 0) /N I versus its detuning pump for g BB = 0.5 h3 ω/m, N B = 100, N I = 1 and varying g BI .Vertical dashed lines indicate the resonant detunings + .[(c)-(k)]Time-resolved probe spectra at different g BI , bosonic (B) or fermionic (F) impurity numbers N I = 1, 2 with g II = 0, and for various short dark times, t d (see legend).In all cases insets illustrate the spatiotemporal evolution of the impurity's one-body density and dashed lines mark the instants for which the probe spectrum is provided.

FIG. 3 .
FIG. 3. [(a)-(i)]Probe spectra for different g BI , N I , and impurity flavors at various dark times t d deep in the evolution (see legend).The remaining system parameters are the same as in Fig.1.

FIG. 6 .
FIG. 6.(a) Schematic illustration of the employed PPS pulse sequences.The involved spin configuration at each state of the dynamics is also provided.(b) Expected time evolution of the population of spin-↑ atoms during the PPS sequence.The inset depicts the evolution of N↑ (t ) /N I during the application of the optical burst pulse (blue line) and the approximation of employing the projection operator P at t = 0 (dark red line).

1 FIG. 7 .
FIG. 7. (a) Time evolution of the structure factor of one impurity at different impurity-BEC interaction strengths (see legend).(b) Dynamics of the structure factor depicted in (a) of one impurity for g BI = 1.5 h3 ω/m over a longer timescale.The bosonic bath contains N B = 100 bosons.The system is initialized in its ground state with g BI = 0.