Many-body quantum dynamics and induced correlations of Bose polarons

We consider a highly particle number imbalanced Bose–Bose mixture composed of N_I = 2 bosonic impurities (I) possessing an additional pseudospin-1/2 degree of freedom [76], which are immersed in a bosonic gas of N_B = 100 structureless bosons (B). Moreover, the mixture is assumed to be mass-balanced, namely m_B = m_I ≡ m and each species is confined in the same one-dimensional external harmonic oscillator potential of frequency ω_B = ω_I = ω. Such a system can be experimentally realized by considering e.g. a ⁸⁷Rb BEC where the majority species resides in the hyperfine state $| F=2,{m}_{F}=1\rangle$ and the pseudospin degree of freedom of the impurities refers for instance to the internal states $| \uparrow \rangle \equiv | F=1,{m}_{F}=1\rangle$ and $| \downarrow \rangle \equiv | F=1,{m}_{F}=-1\rangle$ [77, 78]. Alternatively, it can be realized to a good approximation by a mixture of isotopes of ⁸⁷Rb for the bosonic gas and two hyperfine states of ⁸⁵Rb for the impurities. The underlying MB Hamiltonian of this system reads

$$\begin{eqnarray}&&\hat{H}={\hat{H}}_{B}^{0}+\sum _{a=\uparrow ,\downarrow }{\hat{H}}_{a}^{0}+\sum _{a=\uparrow ,\downarrow }{\hat{H}}_{{aa}}^{{\rm{int}}}+{\hat{H}}_{\uparrow \downarrow }^{{\rm{int}}}+{\hat{H}}_{{BB}}^{{\rm{int}}}+{\hat{H}}_{{BI}}^{{\rm{int}}}.\end{eqnarray} \tag{ 1 }$$

The non-interacting Hamiltonian of the bosonic gas is ${\hat{H}}_{B}^{0}=\int {\rm{d}}{x}\,{\hat{{\rm{\Psi }}}}_{B}^{\dagger }(x)\left(-\tfrac{{{\hslash }}^{2}}{2m}\tfrac{{{\rm{d}}}^{2}}{{{\rm{d}}{x}}^{2}}+\tfrac{1}{2}m{\omega }^{2}{x}^{2}\right){\hat{{\rm{\Psi }}}}_{B}(x)$, while for the impurities it reads ${\hat{H}}_{a}^{0}=\int {\rm{d}}{x}\,{\hat{{\rm{\Psi }}}}_{a}^{\dagger }(x)\left(-\tfrac{{{\hslash }}^{2}}{2m}\tfrac{{{\rm{d}}}^{2}}{{{\rm{d}}{x}}^{2}}+\tfrac{1}{2}m{\omega }^{2}{x}^{2}\right){\hat{{\rm{\Psi }}}}_{a}(x)$ with $a=\left\{\uparrow ,\downarrow \right\}$ being the indices of the spin components. Here ${\hat{{\rm{\Psi }}}}_{\sigma }(x)$ refers to the bosonic field-operator of either the bosonic gas (σ = B) or the impurity ($\sigma =a=\left\{\uparrow ,\downarrow \right\}$) atoms. Furthermore, we operate in the ultracold regime where s-wave scattering is the dominant interaction process. Therefore both the intra- and the intercomponent interactions can be adequately modeled by contact ones. The contact intraspecies interaction of the BEC component is modeled by ${\hat{H}}_{{BB}}^{{\rm{int}}}={g}_{{BB}}\int {\rm{d}}{x}\,{\hat{{\rm{\Psi }}}}_{B}^{\dagger }(x){\hat{{\rm{\Psi }}}}_{B}^{\dagger }(x){\hat{{\rm{\Psi }}}}_{B}(x){\hat{{\rm{\Psi }}}}_{B}(x)$ and between the impurities via ${\hat{H}}_{{aa}^{\prime} }^{{\rm{int}}}$ = ${g}_{{aa}^{\prime} }\int {\rm{d}}{x}{\hat{{\rm{\Psi }}}}_{a}^{\dagger }(x){\hat{{\rm{\Psi }}}}_{a^{\prime} }^{\dagger }(x){\hat{{\rm{\Psi }}}}_{a^{\prime} }(x){\hat{{\rm{\Psi }}}}_{a}(x)$ where either $a=a^{\prime} =\uparrow ,\downarrow$ or $a=\uparrow$, $a^{\prime} =\downarrow$. Note also that we assume ${g}_{\uparrow \uparrow }\,=\,{g}_{\downarrow \downarrow }\,=\,{g}_{\uparrow \downarrow }\equiv {g}_{{II}}$. Most importantly, we consider that only the pseudospin-$\uparrow$ component of the impurities interacts with the bosonic gas while the pseudospin-$\downarrow$ is non-interacting. The resulting intercomponent interaction is ${\hat{H}}_{{BI}}^{{\rm{int}}}={g}_{{BI}}\int {\rm{d}}{x}\,{\hat{{\rm{\Psi }}}}_{B}^{\dagger }(x){\hat{{\rm{\Psi }}}}_{\uparrow }^{\dagger }(x){\hat{{\rm{\Psi }}}}_{\uparrow }(x){\hat{{\rm{\Psi }}}}_{B}(x)$, where ${g}_{{BI}}\equiv {g}_{B\uparrow }$ and ${g}_{B\downarrow }\,=\,0$.

In all of the above-mentioned cases, the effective one-dimensional coupling strength [79] is given by ${g}_{\sigma \sigma ^{\prime} }=\tfrac{2{{\hslash }}^{2}{a}_{\sigma \sigma ^{\prime} }^{s}}{\mu {a}_{\perp }^{2}}{\left(1-\left|\zeta (1/2)\right|{a}_{\sigma \sigma ^{\prime} }^{s}/\sqrt{2}{a}_{\perp }\right)}^{-1}$, where $\sigma ,\sigma ^{\prime} =B,\uparrow ,\downarrow$ and $\mu =\tfrac{m}{2}$ is the reduced mass. The transversal length scale is ${a}_{\perp }=\sqrt{{\hslash }/\mu {\omega }_{\perp }}$ with ω_⊥ being the transversal confinement frequency and ${a}_{\sigma \sigma ^{\prime} }^{s}$ denotes the three-dimensional s-wave scattering length within (σ = σ') or between ($\sigma \ne \sigma ^{\prime}$) the components. In a corresponding experiment, g_σσ' can be tuned either via ${a}_{\sigma \sigma ^{\prime} }^{s}$ with the aid of Feshbach resonances [15, 16] or by adjusting ω_⊥ using confinement-induced resonances [79]. In the following, the MB Hamiltonian of equation (1)is rescaled with respect to ℏω. As a consequence, length, time, and interaction strengths are given in units of $\sqrt{\tfrac{{\hslash }}{m\omega }}$, ω⁻¹ and $\sqrt{\tfrac{{{\hslash }}^{3}\omega }{m}}$ respectively.

To study the quench dynamics, the above-described multicomponent system is initially prepared in its ground state configuration for fixed g_BB = 0.5 and g_BI = 0 and either g_II = 0 or g_II = 0.2. In this way, the case of two non-interacting and that of weakly interacting impurities are investigated. This initial (ground) state emulates a system prepared in the $| 1,-1\rangle =| \downarrow {\rangle }_{1}\otimes | \downarrow {\rangle }_{2}$ configuration for the spin degree of freedom i.e. where the impurity-BEC interaction is zero. Note that the spinor part of the wavefunction is expressed in the basis of the total spin i.e. $| S,{S}_{z}\rangle$ [80]. Accordingly, the spatial part $| {{\rm{\Psi }}}_{{BI}}^{0}\rangle$ of the ground state of the system obeys the following eigenvalue equation $\left(\hat{H}-{\hat{H}}_{{BI}}\right)| {{\rm{\Psi }}}_{{BI}}^{0}\rangle | 1,-1\rangle ={E}_{0}| {{\rm{\Psi }}}_{{BI}}^{0}\rangle | 1,-1\rangle$, with E₀ being the corresponding eigenenergy and ${\hat{H}}_{{BI}}| {{\rm{\Psi }}}_{{BI}}^{0}\rangle | 1,-1\rangle =0$. To trigger the dynamics we carry out an interspecies interaction quench from g_BI = 0 to a finite positive or negative value of g_BI at t = 0 and monitor the subsequent time-evolution. In a corresponding experiment, this quench protocol can be implemented by using a radiofrequency π/2 pulse with an exposure time much smaller than ω⁻¹ [19]. The pulse acts upon the spin degree of freedom of the impurity, which maps the pseudospin-$\downarrow$ impurities to the superposition state $| {\psi }_{S}{\rangle }_{i}\equiv \tfrac{| \uparrow {\rangle }_{i}+| \downarrow {\rangle }_{i}}{\sqrt{2}}$ with i = 1, 2 [18]. The corresponding MB wavefunction of the system, $| {\rm{\Psi }}(t)\rangle ={{\rm{e}}}^{-{\rm{i}}\hat{H}t/{\hslash }}\left[| {{\rm{\Psi }}}_{{BI}}^{0}\rangle (| {\psi }_{S}{\rangle }_{1}\otimes | {\psi }_{S}{\rangle }_{2})\right]$, is then given by

$$\begin{eqnarray}&&\begin{array}{l}| {\rm{\Psi }}(t)\rangle =\displaystyle \frac{1}{\sqrt{2}}{{\rm{e}}}^{-{\rm{i}}\hat{H}t/{\hslash }}\left[| {{\rm{\Psi }}}_{{BI}}^{0}\rangle | 1,0\rangle \right]+\displaystyle \frac{1}{2}\left({{\rm{e}}}^{-{{iE}}_{0}t/{\hslash }}| {{\rm{\Psi }}}_{{BI}}^{0}\rangle | 1,-1\rangle +{{\rm{e}}}^{-{\rm{i}}\hat{H}t/{\hslash }}| {{\rm{\Psi }}}_{{BI}}^{0}\rangle | 1,1\rangle \right).\end{array}\end{eqnarray} \tag{ 2 }$$

The setup and processes addressed in our work can be experimentally realized utilizing radiofrequency spectroscopy [9, 18, 22, 23, 43] and Ramsey interferometry [18].

To calculate the stationary properties and to track the MB non-equilibrium quantum dynamics of the multicomponent bosonic system discussed above we employ the ML-MCTDHX method [70–72]. This is an ab initio variational method for solving the time-dependent MB Schrödinger equation of atomic mixtures and it is based on the expansion of the total MB wavefunction with respect to a time-dependent and variationally optimized basis tailored to capture both the intra- and the interspecies correlations of a multicomponent system [35, 65, 81, 82].

To include the interspecies correlations, the MB wavefunction ($| {\rm{\Psi }}(t)\rangle$) is first expanded in terms of D distinct species functions, $| {{\rm{\Psi }}}_{i}^{\sigma }(t)\rangle$, for each component σ = B, I, and then expressed according to a truncated Schmidt decomposition [83] of rank D, namely

$$\begin{eqnarray}&&| {\rm{\Psi }}(t)\rangle =\sum _{k=1}^{D}\sqrt{{\lambda }_{k}(t)}| {{\rm{\Psi }}}_{k}^{{B}}(t)\rangle | {{\rm{\Psi }}}_{k}^{{I}}(t)\rangle .\end{eqnarray} \tag{ 3 }$$

Here the time-dependent expansion coefficients λ_k(t) are the Schmidt weights and will be referred to in the following as the natural populations of the kth species function. Evidently, the system is entangled [84] or interspecies correlated when at least two different λ_k(t) possess a non-zero value. If this is not the case, i.e. for λ₁(t) = 1, λ_k>1(t) = 0, the wavefunction is a direct product of two states.

Therefore, in order to account for intraspecies correlations, each of the above-mentioned species functions is expressed as a linear superposition of time-dependent number-states, $| \vec{n}(t){\rangle }^{\sigma }$, with time-dependent coefficients ${A}_{i;\vec{n}}^{\sigma }(t)$ as

$$\begin{eqnarray}&&| {{\rm{\Psi }}}_{i}^{\sigma }(t)\rangle =\sum _{\vec{n}}{A}_{i;\vec{n}}^{\sigma }(t)| \vec{n}(t){\rangle }^{\sigma }.\end{eqnarray} \tag{ 4 }$$

Each number state $| \vec{n}(t){\rangle }^{\sigma }$ is a permanent building upon d^σ time-dependent variationally optimized single-particle functions (SPFs) $\left|{\phi }_{l}^{\sigma }(t)\right\rangle$, l = 1, 2, ..., d^σ with occupation numbers $\vec{n}=({n}_{1},\,\ldots ,\,{n}_{{d}^{\sigma }})$. Consecutively, the SPFs are expanded on a time-independent primitive basis. The latter refers to an ${ \mathcal M }$ dimensional discrete variable representation (DVR) for the majority species and it is denoted by $\{\left|k\right\rangle \}$. For the impurities this corresponds to the tensor product $\{\left|k,s\right\rangle \}$ of the DVR basis for the spatial degrees of freedom and the two-dimensional pseudospin-1/2 basis $\{| \uparrow \rangle ,| \downarrow \rangle \}$. Accordingly, each SPF of the impurities is a spinor wavefunction of the form

$$\begin{eqnarray}&&| {\phi }_{j}^{{I}}(t)\rangle =\sum _{k=1}^{{ \mathcal M }}\left({B}_{{jk}\uparrow }^{{I}}(t)| k\rangle | \uparrow \rangle +{B}_{{jk}\downarrow }^{{I}}(t)| k\rangle | \downarrow \rangle \right),\end{eqnarray} \tag{ 5 }$$

with ${B}_{{jk}\uparrow }^{{I}}(t)$ $[{B}_{{jk}\downarrow }^{{I}}(t)]$ being the time-dependent expansion coefficients of the pseudospin-$\uparrow$ $[\downarrow ]$ (see also [35, 82] for a more detailed discussion).

The time-evolution of the (N_B + N_I)-body wavefunction $\left|{\rm{\Psi }}(t)\right\rangle$ governed by the Hamiltonian of equation (1) is obtained via solving the so-called ML-MCTDHX equations of motion [70]. The latter are determined by utilizing e.g. the Dirac–Frenkel [85, 86] variational principle for the generalized ansatz introduced in equations (3)–(5). This procedure results in a set of D² linear differential equations of motion for the λ_k(t) coefficients which are coupled to $D\left(\tfrac{({N}_{B}+{d}^{B}-1)!}{{N}_{B}!({d}^{B}-1)!}+\tfrac{({N}_{I}+{d}^{I}-1)!}{{N}_{I}!({d}^{I}-1)!}\right)$ nonlinear integrodifferential equations for the species functions and d^B + d^I nonlinear integrodifferential equations for the SPFs.

A main aspect of the ansatz outlined above is the expansion of the system's MB wavefunction with respect to a time-dependent and variationally optimized basis. The latter allows to efficiently take into account the intra- and intercomponent correlations of the system using a computationally feasible basis size. In the present case the Bose gas consists of a large number of weakly interacting particles and therefore its intracomponent correlations are suppressed. As a consequence they can be adequately captured by employing a small number of orbitals, d^B < 4. Additionally, the number of impurities, N_I < 3, is small giving rise to a small number of integrodifferential equations allowing us to employ many orbitals, d_I, and thus account for strong impurity–impurity and impurity-BEC correlations. Therefore, the number of the resulting equations of motion that need to be solved is numerically tractable. Since our method is variational, its validity is determined upon examining its convergence. For details on the precision of our simulations see appendix.

To study the quench-induced dynamics of each species at the single-particle level we calculate the one-body reduced density matrix for each species [87, 88]

$$\begin{eqnarray}&&{\rho }_{\sigma }^{(1)}(x,x^{\prime} ;t)=\langle {\rm{\Psi }}(t)| {\hat{{\rm{\Psi }}}}_{\sigma }^{\dagger }(x){\hat{{\rm{\Psi }}}}_{\sigma }(x^{\prime} )| {\rm{\Psi }}(t)\rangle .\end{eqnarray} \tag{ 6 }$$

Here, ${\hat{{\rm{\Psi }}}}_{\sigma }(x)$ is the σ-species bosonic field operator acting at position x and satisfying the standard bosonic commutation relations [89]. For simplicity, we will use in the following the one-body densities for each species i.e. ${\rho }_{\sigma }^{(1)}(x;t)\equiv {\rho }_{\sigma }^{(1)}(x,x^{\prime} =x;t)$, which is a quantity that is experimentally accessible via averaging over a sample of single-shot images [65, 90, 91]. We remark that the eigenfunctions and eigenvalues of ${\rho }_{\sigma }^{(1)}(x,x^{\prime} ;t)$ are termed natural orbitals ${\varphi }_{i}^{\sigma }(x;t)$ and natural populations ${n}_{i}^{\sigma }(t)$ [65, 70] respectively. In this sense, each bosonic subsystem is called intraspecies correlated if more than a single natural population possess a non-zero contribution. Otherwise, i.e. for ${n}_{1}^{\sigma }(t)=1$ and ${n}_{i\gt 1}^{\sigma }(t)=0$, the corresponding subsystem is said to be fully coherent and the MB wavefunction (equations (3), (5)) reduces to a mean-field product ansatz [92, 93].

To unveil the role of impurity–impurity correlations following the interspecies interaction quench we calculate the time-evolution of the corresponding diagonal of the two-body reduced density matrix

$$\begin{eqnarray}&&{\rho }_{{aa}^{\prime} }^{(2)}({x}_{1},{x}_{2};t)=\langle {\rm{\Psi }}(t)| {\hat{{\rm{\Psi }}}}_{a}^{\dagger }({x}_{1}){\hat{{\rm{\Psi }}}}_{a^{\prime} }^{\dagger }({x}_{2}){\hat{{\rm{\Psi }}}}_{a^{\prime} }({x}_{2}){\hat{{\rm{\Psi }}}}_{a}({x}_{1})| {\rm{\Psi }}(t)\rangle ,\end{eqnarray} \tag{ 7 }$$

where $a,a^{\prime} =\uparrow ,\downarrow$. The two-body reduced density matrix refers to the probability of finding simultaneously one pseudospin-a boson at x₁ and a pseudospin-a' boson at x₂ [65, 66]. Moreover, it provides insights into the spatially resolved dynamics of the two impurities with respect to one another. Indeed, the impurities are dressed by the excitations of the bosonic gas forming quasiparticles which in turn can move independently or interact, and possibly form a bound state [8, 10, 42, 94].

To capture the emerging effective interactions between the two bosonic impurities we monitor their relative distance [9, 42] given by

$$\begin{eqnarray}&&\left\langle {r}_{{aa}}(t)\right\rangle =\displaystyle \frac{\int {{\rm{d}}{x}}_{1}{{\rm{d}}{x}}_{2}| {x}_{1}-{x}_{2}| {\rho }_{{aa}}^{(2)}({x}_{1},{x}_{2};t)}{\left\langle {\rm{\Psi }}(t)| {\hat{N}}_{a}\left({\hat{N}}_{a}-1\right)| {\rm{\Psi }}(t)\right\rangle }.\end{eqnarray} \tag{ 8 }$$

Here, ${\hat{N}}_{a}$ with $a=\uparrow ,\downarrow$ is the number operator that measures the number of bosons in the spin-a state. Experimentally, $\left\langle {r}_{{aa}}(t)\right\rangle$ can be probed via in situ spin-resolved single-shot measurements on the spin-a state [91]. More precisely, each image gives an estimate of $\left\langle {r}_{{aa}}(t)\right\rangle$ between the bosonic impurities if their position uncertainty is assured to be adequately small [91]. Subsequently, $\left\langle {r}_{{aa}}(t)\right\rangle$ is obtained by averaging over several such images.

To examine the quench-induced dynamics of the two spinor bosonic impurities we first determine the time-evolution of the total spin polarization (contrast) $| \left\langle \hat{{\bf{S}}}(t)\right\rangle | =\sqrt{{\left\langle {\hat{S}}_{x}(t)\right\rangle }^{2}+{\left\langle {\hat{S}}_{y}(t)\right\rangle }^{2}}$ which enables us to infer the dressing of the impurities during the dynamics [18]. Note that $\left\langle {\hat{S}}_{z}(t)\right\rangle =\left\langle {\hat{S}}_{z}(t=0)\right\rangle =0$ since $\left[{\hat{S}}_{z},\hat{H}\right]=0$ and the spin operator in the kth direction (k = x, y, z) is given by ${\hat{S}}_{k}=(1{/N}_{I})\int {\rm{d}}{x}{\sum }_{{ab}}{\hat{{\rm{\Psi }}}}_{a}^{\dagger }(x){\sigma }_{{ab}}^{k}{\hat{{\rm{\Psi }}}}_{b}(x)$, with ${\sigma }_{{ab}}^{k}$ denoting the Pauli matrices. The contrast for a single impurity has been extensively studied [25, 39, 43, 97] and it is related to the so-called Ramsey response [18] and therefore the structure factor. The time-dependent overlap between the interacting and the non-interacting states is given by

$$\begin{eqnarray}&&| \left\langle \hat{{\bf{S}}}(t)\right\rangle {| }^{2}=| \left\langle {\tilde{{\rm{\Psi }}}}_{{BI}}^{0}| {{\rm{e}}}^{{\rm{i}}{\tilde{E}}_{0}t/{\hslash }}{{\rm{e}}}^{-{\rm{i}}\hat{\tilde{H}}t/{\hslash }}| {\tilde{{\rm{\Psi }}}}_{{BI}}^{0}\right\rangle {| }^{2}\equiv {\left|{S}_{1}(t)\right|}^{2},\end{eqnarray} \tag{ 10 }$$

where $| {\tilde{{\rm{\Psi }}}}_{{BI}}^{0}\rangle$ is the spatial part of the MB ground state wavefunction of a single impurity with energy ${\tilde{E}}_{0}$ when g_BI = 0. $\hat{\tilde{H}}=\hat{P}\hat{H}\hat{P}$ with $\hat{P}$ being the projector operator to the spin-$\uparrow$ configuration, and $\hat{H}$ denotes the postquench Hamiltonian (equation (1)). Note also that the contrast is chosen here to take values in the interval [0, 1]. From equation (10) zero contrast implies that the overlap between the interacting and the non-interacting states vanishes signifying an orthogonality catastrophe phenomenon [52, 97]. On the other hand, if $| \left\langle \hat{{\bf{S}}}(t)\right\rangle {| }^{2}=1$ then the non-interacting and the interacting states coincide and no quasiparticle is formed. Therefore only in the case that $0\lt | \left\langle \hat{{\bf{S}}}(t)\right\rangle {| }^{2}\lt 1$ we can infer the dressing of the impurity and the formation of a quasiparticle.

When increasing the number of impurity atoms to N_I > 1, $| \left\langle \hat{{\bf{S}}}(t)\right\rangle {| }^{2}$ is more complex since additional spin states contribute to the MB wavefunction (see equation (2)). To understand the interpretation of $| \left\langle \hat{{\bf{S}}}(t)\right\rangle {| }^{2}$ during the dynamics we therefore first discuss it for the case of two impurities. The contrast of two pseudospin-1/2 bosonic impurities reads

$$\begin{eqnarray}&&| \left\langle \hat{{\bf{S}}}(t)\right\rangle {| }^{2}=\displaystyle \frac{1}{4}{\left|A(| 1,0\rangle ;| 1,-1\rangle )+{A}^{* }(| 1,0\rangle ;| 1,1\rangle )\right|}^{2},\end{eqnarray} \tag{ 11 }$$

where the spatial overlap between two different spin configurations namely $| S,{S}_{z}\rangle$ and $| S^{\prime} ,{S}_{z}^{{\prime} }\rangle$ is defined as [80]

$$\begin{eqnarray}\begin{array}{rcl}A(| S,{S}_{z}\rangle ;| S^{\prime} ,{S}_{z}^{{\prime} }\rangle ) & \equiv & \left[\langle S,{S}_{z}| \langle {{\rm{\Psi }}}_{{BI}}^{0}| \right]{{\rm{e}}}^{{\rm{i}}\hat{H}t/{\hslash }}| S,{S}_{z}\rangle \langle S^{\prime} ,{S}_{z}^{{\prime} }| {{\rm{e}}}^{-{\rm{i}}\hat{H}t/{\hslash }}\left[| {{\rm{\Psi }}}_{{BI}}^{0}\rangle | S^{\prime} ,{S}_{z}^{{\prime} }\rangle \right]\\ & = & \ \displaystyle \int {{\rm{d}}{x}}^{{N}_{B}}{{\rm{d}}{x}}^{{N}_{I}}{{\rm{\Psi }}}_{S,{S}_{z}}^{* }({\vec{x}}^{B},{\vec{x}}^{I};t){{\rm{\Psi }}}_{S^{\prime} ,{S}_{z}^{{\prime} }}({\vec{x}}^{B},{\vec{x}}^{I};t),\end{array}\end{eqnarray} \tag{ 12 }$$

with ${{\rm{\Psi }}}_{S,{S}_{z}}({\overrightarrow{x}}^{B},{\overrightarrow{x}}^{I};t)={\textstyle \tfrac{ \langle {\overrightarrow{x}}^{B},{\overrightarrow{x}}^{I}| \langle S,{S}_{z}| {\rm{\Psi }}(t) \rangle }{{\left|\left| \langle S,{S}_{z}| {\rm{\Psi }}(0) \rangle \right|\right|}^{2}}}$ referring to the spatial wavefunction corresponding to the spin configuration $| S,{S}_{z}\rangle$ and $| {{\rm{\Psi }}}_{{BI}}^{0}\rangle$ being the spatial part of the initial MB state for two impurities. Also, ${\overrightarrow{x}}^{B}=({x}_{1}^{B},\ldots ,{x}_{{N}_{B}}^{B})$ and ${\overrightarrow{x}}^{I}=({x}_{1}^{I},\ldots ,{x}_{{N}_{B}}^{I})$ refer to the coordinates of each bath and impurity particle, respectively. In particular in our case we consider two pseudospin-1/2 bosons where $| 1,1\rangle \equiv | \uparrow {\rangle }_{1}\otimes | \uparrow {\rangle }_{2}$, $| 1,-1\rangle \equiv | \downarrow {\rangle }_{1}\otimes | \downarrow {\rangle }_{2}$, $| 1,0\rangle \equiv \tfrac{| \uparrow {\rangle }_{1}\otimes | \downarrow {\rangle }_{2}+| \downarrow {\rangle }_{1}\otimes | \uparrow {\rangle }_{2}}{\sqrt{2}}$.The relevant overlaps read $A(| 1,0\rangle ;| 1,-1\rangle )$ = ${{\rm{e}}}^{-{{\rm{i}}E}_{0}t/\hslash }\int {{\rm{d}}x}^{{N}_{B}}{d}^{{N}_{B}}{\overrightarrow{x}}^{B}{d}^{{N}_{I}}{\overrightarrow{x}}^{I}{{\rm{\Psi }}}_{1,0}^{\ast }({\overrightarrow{x}}^{B},{\overrightarrow{x}}^{I};t){{\rm{\Psi }}}_{BI}^{0}({\overrightarrow{x}}^{B},{\overrightarrow{x}}^{I};0)$ and $A(| 1,0\rangle ;| 1,1\rangle )$ = $\int {d}^{{N}_{B}}{\overrightarrow{x}}^{B}{d}^{{N}_{I}}{\overrightarrow{x}}^{I}{{\rm{\Psi }}}_{1,0}^{\ast }({\overrightarrow{x}}^{B},{\overrightarrow{x}}^{I};t){{\rm{\Psi }}}_{1,1}({\overrightarrow{x}}^{B},{\overrightarrow{x}}^{I};t)$. Recall that a quasiparticle is a free particle that is dressed by the excitations of a bosonic bath via their mutual interactions. As a consequence, ${{\rm{\Psi }}}_{{BI}}^{0}({\vec{x}}^{B},{\vec{x}}^{I})$ refers to the wavefunction where no polaron quasiparticle exists since it is the ground state wavefunction of the system with ${g}_{B\uparrow }\,=\,0$. Moreover, ${{\rm{\Psi }}}_{\mathrm{1,0}}({\vec{x}}^{B};{\vec{x}}^{I})$ and ${{\rm{\Psi }}}_{\mathrm{1,1}}({\vec{x}}^{B};{\vec{x}}^{I})$ denote the wavefunctions where a single and two impurities respectively interact with the bosonic gas and therefore describe the formation of a single and two polarons, respectively. Accordingly, $A(| 1,0\rangle ;| 1,-1\rangle )$ provides the overlap between the state of a single and no impurities interacting with the bath, while $A(| 1,0\rangle ;| 1,1\rangle )$ is the overlap between a single and two impurities interacting with the bath.

As a result, $| \left\langle \hat{{\bf{S}}}(t)\right\rangle {| }^{2}=1$ means that $A(| 1,0\rangle ;| 1,-1\rangle )=A(| 1,0\rangle ,| 1,1\rangle )={{\rm{e}}}^{{\rm{i}}\varphi }$ where φ is a phase factor. The fact that $\left|A(| 1,0\rangle ;| 1,-1\rangle )\right|=1$ implies that the spatial state of a single impurity interacting with the bath is the same as the non-interacting one, except for a possible phase factor, and therefore a quasiparticle is not formed. Moreover since also $\left|A(| 1,0\rangle ,| 1,1\rangle )\right|=1$ it holds that the state of a single pseudospin-$\uparrow$ interacting impurity coincides with the state of two pseudospin-$\uparrow$ impurities interacting with the bath and as a consequence with a bare particle due to $\left|A(| 1,0\rangle ;| 1,-1\rangle )\right|=1$. Thus, $| \left\langle \hat{{\bf{S}}}(t)\right\rangle {| }^{2}=1$ implies that there is no quasiparticle formation. On the contrary for $| \left\langle \hat{{\bf{S}}}(t)\right\rangle {| }^{2}=0$ either $A(| 1,0\rangle ;| 1,-1\rangle )=A(| 1,0\rangle ,| 1,1\rangle )=0$ or $A(| 1,0\rangle ;| 1,-1\rangle )=-{A}^{* }(| 1,0\rangle ,| 1,1\rangle )$ should be satisfied. In the former case we can deduce the occurrence of an orthogonality catastrophe phenomenon as in the single impurity case while the latter scenario is given by the destructive interference of the $A(| 1,0\rangle ;| 1,-1\rangle )$ and $A(| 1,0\rangle ,| 1,1\rangle )$ terms However, for $0\lt | \left\langle \hat{{\bf{S}}}(t)\right\rangle {| }^{2}\lt 1$ the corresponding overlaps acquire finite values and a quasiparticle can be formed.

Notice also that in the special case of ${g}_{\uparrow \downarrow }\,=\,0$ and ${g}_{\downarrow \downarrow }\,=\,0$ (but g_↑↑ arbitrary) it can be shown that $A(| 1,0 \rangle ;| 1,-1 \rangle )= \langle {\tilde{{\rm{\Psi }}}}_{BI}^{0}| {{\rm{e}}}^{{\rm{i}}{\hat{P}}_{0}\hat{H}{\hat{P}}_{0}t/\hslash }{{\rm{e}}}^{{\rm{i}}{\tilde{E}}_{0}t/\hslash }| {\tilde{{\rm{\Psi }}}}_{BI}^{0} \rangle \equiv {S}_{1}(t)$, where ${\hat{P}}_{0}$ refers to the projection operator to the spin state $| 1,0 \rangle$. The latter is exactly the contrast or the structure factor of a single impurity (equation (10)). Indeed $| {{\rm{\Psi }}}_{{BI}}^{0}\rangle =| {\tilde{{\rm{\Psi }}}}_{{BI}}^{0}\rangle \otimes | {\psi }_{I}^{0}\rangle$ for ${g}_{\downarrow \downarrow }\,=\,0$ holds where $| {\psi }_{I}^{0}\rangle$ is the single-particle ground state of the impurity while $| {\tilde{{\rm{\Psi }}}}_{{BI}}^{0}\rangle$ and $| {{\rm{\Psi }}}_{{BI}}^{0}\rangle$ refer to the spatial part of the MB ground state wavefunction of a single (energy ${\tilde{E}}_{0}$) and two impurities (energy E₀), respectively. Additionally $\hat{H}$ is the postquench Hamiltonian given by equation (1). Consequently, the contrast in this special case acquires the simplified form

$$\begin{eqnarray}&&| \left\langle \hat{{\bf{S}}}(t)\right\rangle {| }^{2}=\displaystyle \frac{1}{4}{\left|{S}_{1}(t)+{A}^{* }(| 1,0\rangle ;| 1,1\rangle )\right|}^{2}.\end{eqnarray} \tag{ 13 }$$

Evidently, here $| \left\langle \hat{{\bf{S}}}(t)\right\rangle |$ depends explicitly on the structure factor S₁(t) of a single impurity allowing for a direct interpretation of the dynamical dressing of the two impurities with respect to the single impurity case discussed in [35]. In the following, ${g}_{\uparrow \uparrow }\,=\,{g}_{\downarrow \downarrow }\,=\,{g}_{\uparrow \downarrow }\equiv {g}_{{II}}$ and as a consequence ${g}_{\uparrow \downarrow }\,=\,0$, ${g}_{\downarrow \downarrow }\,=\,0$ is encountered for g_II = 0 while the general case of equation (12) applies for the case of g_II = 0.2 analyzed below.

The dynamics of the two particle contrast $| \left\langle \hat{{\bf{S}}}(t)\right\rangle |$ is presented in figures 3(a)–(c) for both attractive and repulsive postquench interspecies interactions ${g}_{B\uparrow }$. In particular, $| \left\langle \hat{{\bf{S}}}(t)\right\rangle |$ is shown for either two non-interacting (figure 3(a)) or interacting (figure 3(b)) impurities and N_B = 100 as well as for a few-body bosonic gas with N_B = 10 and g_II = 0 (figure 3(c)). In all cases, six different dynamical regions with respect to ${g}_{B\uparrow }$ can be identified marked as R_I, R_II, R_III, R_IV, ${R}_{{II}}^{{\prime} }$ and ${R}_{{III}}^{{\prime} }$. Focusing on the system with N_B = 100 and g_II = 0 these regions correspond to $-0.2\leqslant {g}_{B\uparrow }^{{R}_{I}}\lt 0.2$, $0.2\leqslant {g}_{B\uparrow }^{{R}_{{II}}}\lt 0.4$, $0.4\leqslant {g}_{B\uparrow }^{{R}_{{III}}}\lt 1$, $1\leqslant {g}_{B\uparrow }^{{R}_{{IV}}}\lt 5$, $-0.5\leqslant {g}_{B\uparrow }^{{R}_{{II}}^{{\prime} }}\lt -0.2$ and $-1\leqslant {g}_{B\uparrow }^{{R}_{{III}}^{{\prime} }}\lt -0.5$ respectively (figure 3(a)). Specifically, within the very weakly interacting region R_I the contrast is essentially unperturbed remaining unity in the course of the time-evolution and therefore there is no quasiparticle formation. For postquench interactions lying within R_II or ${R}_{{II}}^{{\prime} }$ the contrast performs small and constant amplitude oscillations, weakly deviating from $| \left\langle \hat{{\bf{S}}}(t=0)\right\rangle | =1$ (figure 3(f)). This behavior indicates the generation of two long-lived coherent quasiparticles (see also section 4.3). Entering the intermediate repulsive interaction region R_III, $| \left\langle \hat{{\bf{S}}}(t)\right\rangle |$ exhibits large amplitude ($0\lt | \left\langle \hat{{\bf{S}}}(t)\right\rangle | \lt 1$) multifrequency temporal oscillations (figure 3(f)). The latter signifies the dynamical formation of two Bose polarons which are coupled with higher-order excitations of the bosonic bath when compared to regions R_II and ${R}_{{II}}^{{\prime} }$ as we shall expose in section 4.4.1. For intermediate attractive interactions (region ${R}_{{III}}^{{\prime} }$) $| \left\langle \hat{{\bf{S}}}(t)\right\rangle |$ undergoes large amplitude oscillations taking values in the interval $0\lt | \left\langle \hat{{\bf{S}}}(t)\right\rangle | \lt 1$ (figure 3(f)). This response of $| \left\langle \hat{{\bf{S}}}(t)\right\rangle |$ again signals quasiparticle formation. However, in addition to this dynamical dressing the destructive ($| \left\langle \hat{{\bf{S}}}(t)\right\rangle | =0$) and the constructive ($| \left\langle \hat{{\bf{S}}}(t)\right\rangle | \approx 1$) interference between the states of a single and two Bose polarons can be seen (see also equation (11) and its interpretation in section 4.1).

Zoom In Zoom Out Reset image size

For strong repulsive interactions lying within R_IV the contrast shows a fastly decaying amplitude at short evolution times (0 < t < 2) and subsequently fluctuates around zero (figure 3(f)). This latter behavior of $| \left\langle \hat{{\bf{S}}}(t)\right\rangle {| }^{}\to 0$ is a manifestation of an orthogonality catastrophe phenomenon of the spontaneously generated short-lived (0 < t < 2) Bose polarons. It is a consequence of the spatial phase separation between the impurity and the bosonic bath (see also figure 5(h) and the discussion in section 4.4.1), where the impurity prefers to reside at the edges of the BEC background, see also figure 2(c). Note that this behavior is also supported by the effective potential of the impurities, see equation (9). Most importantly this process results in an energy transfer from the impurity to the BEC, which prohibits the revival of the dynamical state of the impurity to its initial one, implying $| \left\langle \hat{{\bf{S}}}(t)\right\rangle {| }^{}\ll 1$. Such a mechanism has been also identified to occur for the case of a single impurity, see [35].

The emergence of the different dynamical regions in the evolution of the contrast holds equally when the size of the bath decreases to N_B = 10 (figure 3(c)). For such a few-body scenario region R_II, where coherently long-lived quasiparticles are formed, becomes slightly wider, i.e. $0.2\leqslant {g}_{B\uparrow }^{{R}_{{II}}}\lt 0.6$, compared to the N_B = 100 case. The most notable difference between the few and the many particle bath takes place in the intermediate interaction region R_III. The latter, occurs now at $0.6\leqslant {g}_{B\uparrow }^{{R}_{{III}}}\lt 1.8$, with $| \left\langle \hat{{\bf{S}}}(t)\right\rangle |$ performing large amplitude multifrequency oscillations implying in turn the formation of highly excited polaronic states. Note that the amplitude of the oscillations of $| \left\langle \hat{{\bf{S}}}(t)\right\rangle |$ here is larger than in the N_B = 100 case (figure 3(a)). Additionally, we observe that $| \left\langle \hat{{\bf{S}}}(t)\right\rangle |$ decreases smoothly as ${g}_{B\uparrow }$ increases, which is in sharp contrast to the N_B = 100 case. Recall that such a smooth behavior occurring in the few-body scenario has already been identified in our discussion of the ground state properties and in particular when inspecting the relative distance between the impurities. Also, the oscillations of $| \left\langle \hat{{\bf{S}}}(t)\right\rangle |$($0\lt | \left\langle \hat{{\bf{S}}}(t)\right\rangle | \lt 1$) for intermediate attractive interactions (region ${R}_{{III}}^{{\prime} }$) being a consequence of the destructive ($| \left\langle \hat{{\bf{S}}}(t)\right\rangle | =0$) and constructive ($| \left\langle \hat{{\bf{S}}}(t)\right\rangle | \approx 1$) interference between the states of a single and two Bose polarons are much more prevalent and regular for N_B = 10 as compared to the N_B = 100 case. Concluding, we can infer that the overall phenomenology of the dynamical formation of quasiparticles as imprinted in the contrast is similar for N_B = 10 and N_B = 100.

To test the effect of the number of impurities on the interaction intervals of quasiparticle formation we also consider the case of N_I = 3 non-interacting, g_II = 0, bosons immersed in a few-body bath of N_B = 10 atoms. The dynamics of the corresponding contrast for this system following a quench from ${g}_{B\uparrow }\,=\,0$ to a finite either attractive or repulsive ${g}_{B\uparrow }$ is illustrated in figure 3(d). As it can be seen, $| \left\langle \hat{{\bf{S}}}(t)\right\rangle |$ shows a similar behavior to the case of two impurities (figure 3(c)) but the regions of finite contrast become narrower. Particularly, the intermediate repulsive interaction region here occurs for $0.5\leqslant {g}_{B\uparrow }^{{R}_{{III}}}\lt 1.5$ instead of $0.6\leqslant {g}_{B\uparrow }^{{R}_{{III}}}\lt 1.8$ for N_I = 2. Additionally, $| \left\langle \hat{{\bf{S}}}(t)\right\rangle |$ acquires lower values within the regions R_III and ${R}_{{III}}^{{\prime} }$ for more impurities. Moreover, for N_I = 3 within ${R}_{{III}}^{{\prime} }$ we observe a pronounced dephasing of the contrast which is absent for the N_I = 2 case, see figures 3(e₁), (e₂). As a consequence, we can deduce that the basic characteristics of the regions of dynamical polaron formation do not significantly change for a larger number of impurities in the regime N_I ≪ N_B.

Finally, we discuss $| \left\langle \hat{{\bf{S}}}(t)\right\rangle |$ for weakly interacting impurities. Comparing the temporal evolution of $| \left\langle \hat{{\bf{S}}}(t)\right\rangle |$ for g_II = 0.2 (figure 3(b)) to the one for g_II = 0 (figure 3(a)) we observe that the extent of the above-described dynamical regions (R_I, R_II, R_III, R_IV, ${R}_{{II}}^{{\prime} }$ and ${R}_{{III}}^{{\prime} }$) can be tuned via g_II. For instance, region R_II occurs at $0.2\leqslant {g}_{B\uparrow }^{{R}_{{II}}}\lt 0.4$ for g_II = 0.2 instead of $0.2\leqslant {g}_{B\uparrow }^{{R}_{{II}}}\lt 0.5$ when g_II = 0, while region R_III takes place at $0.4\leqslant {g}_{B\uparrow }^{{R}_{{III}}}\lt 1.3$ if g_II = 0.2 and within $0.5\leqslant {g}_{B\uparrow }^{{R}_{{III}}}\lt 1$ in the non-interacting scenario. Also region R_IV where the orthogonality catastrophe takes place is shifted to slightly larger interactions for g_II = 0.2 compared to the g_II = 0 case. Interestingly we observe that the contrast within R_III and ${R}_{{III}}^{{\prime} }$ exhibits a decaying tendency for long evolution times t > 50 in the presence of weak impurity–impurity interactions, a behavior which is absent when g_II = 0.

To quantify the excitation spectrum of the impurity we calculate the spectrum of the contrast, namely

$$\begin{eqnarray}&&A({\omega }_{f})=\displaystyle \frac{1}{\pi }\left|{\int }_{0}^{\infty }{\rm{d}}{t}\,{{\rm{e}}}^{{\rm{i}}{\omega }_{f}t}\,| \ \left\langle \hat{{\bf{S}}}(t)\right\rangle \ | {{\rm{e}}}^{{\rm{i}}{\tan }^{-1}\tfrac{\left\langle {\hat{S}}_{x}(t)\right\rangle }{\left\langle {\hat{S}}_{y}(t)\right\rangle }}\right|.\end{eqnarray} \tag{ 14 }$$

Recall that at low impurity densities and weak interspecies interactions it has been shown that $| \left\langle \hat{{\bf{S}}}(t)\right\rangle |$ is proportional to the so-called spectral function of quasiparticles [18, 97, 98]. Figure 4 presents A(ω_f) in the case of a single and two either non-interacting (g_II = 0) or weakly interacting (g_II = 0.2) impurities when N_B = 100 for different interspecies couplings of either sign. Evidently, for weak ${g}_{B\uparrow }$ belonging either to region R_II with ${g}_{B\uparrow }\,=\,0.25$ (figure 4(a)) or ${R}_{{II}}^{{\prime} }$ with ${g}_{B\uparrow }\,=\,-0.25$ (figure 4(d)) we observe a single peak in A(ω_f) located at ω_f ≈ 4.27 and ω_f ≈ −4.39 respectively. This single peak occurs independently of the number of impurities and their intraspecies interactions. Therefore, this peak at small ${g}_{B\uparrow }\,=\,\pm 0.25$ corresponds to the long-time evolution of a well-defined repulsive or attractive Bose polaron respectively. Within region R_III e.g. at ${g}_{B\uparrow }\,=\,0.5$ two dominant peaks occur in A(ω_f) (figure 4(b)) at frequencies ω_f ≈ 8.42 and ω_f ≈ 8.79 for both the N_I = 1 and N_I = 2 cases. Accordingly, these two peaks suggest the formation of a quasiparticle dressed, for higher frequencies, by higher-order excitations of the BEC background.

Zoom In Zoom Out Reset image size

Entering the strongly interspecies repulsive region R_IV a multitude of frequencies are imprinted in the impurity's excitation spectrum e.g. at ${g}_{B\uparrow }\,=\,1.5$, see figure 4(c). The number of the emerging frequencies is larger for the two compared to the single impurity but does not significantly depend on g_II for N_I = 2. For instance, when N_I = 1 mainly three predominant peaks centered at ω_f ≈ 23.75, ω_f ≈ 25.13, and ω_f ≈ 26.26 appear in A(ω_f) whilst for N_I = 2 and g_II = 0 five dominantly contributing frequencies located at ω_f ≈ 22.31, ω_f ≈ 23.81, ω_f ≈ 25.2, ω_f ≈ 26.39 and ω_f ≈ 27.52 occur. These frequency peaks correspond to even higher excited states of the quasiparticle than the ones within the region R_III. We note that for values of ${g}_{B\uparrow }$ deeper in R_IV a variety of low amplitude but large valued frequency peaks occur in A(ω_f). This fact indicates that the impurities tend to populate a multitude of states indicating the manifestation of the polaron orthogonality catastrophe as discussed in [35, 75] (results not shown here).

Turning to intermediate attractive interactions lying within ${R}_{{III}}^{{\prime} }$ such as ${g}_{B\uparrow }\,=\,-0.5$ a single frequency peak can be seen in A(ω_f) whose frequency is shifted towards more negative values for N_I = 2 compared to N_I = 1 and also for increasing g_II (figure 4(e)). Specifically, when N_I = 1 the aforementioned peak occurs at ω_f ≈ −8.79 while for N_I = 2 and g_II = 0 [g_II = 0.2] it lies at ω_f ≈ −8.92 [ω_f ≈ −8.86]. This peak indicates the generation of an attractive Bose polaron. A further increase of the attraction, e.g. ${g}_{B\uparrow }\,=\,-1$, leads to the appearance of three quasiparticle peaks in A(ω_f) when N_I = 2 and either g_II = 0 or g_II = 0.2, centered at ω_f ≈ −18.1, ω_f ≈ −18.35 and ω_f ≈ −18.98, but only one for N_I = 1 with ω_f ≈ −17.91, as shown in figure 4(f). This change of A(ω_f) for increasing N_I within the regions ${R}_{{II}}^{{\prime} }$ and ${R}_{{III}}^{{\prime} }$ demonstrates the prominent role of induced interactions for attractive interspecies ones. More specifically for N_I = 2, A(ω_f) possesses additional quasiparticle peaks as compared to the N_I = 1 case. Indeed, according to equation (11) we can predict at least two peaks at positions ${\omega }_{f}={{\rm{\Delta }}}_{+}^{{N}_{I}=1}=-17.96$ and ${\omega }_{f}=2{{\rm{\Delta }}}_{+}^{{N}_{I}=2}-{{\rm{\Delta }}}_{+}^{{N}_{I}=1}=-18.98$ explaining two of the above identified peaks. The third dominant peak at ω_f = 18.35 appearing in the spectrum is attributed to the occupation of an excited state with S_z = 1 (see also equation (2)) according to equation (11). Recall that the $| 1,1\rangle$ spin state in the time-evolved wavefunction (equation (2)) corresponds to the two polaron case while $| 1,0\rangle$ contains only one polaron and the $| 1,-1\rangle$ describes impurities that do not interact with the bath and thus no polarons. The aforementioned population of the additional polaronic states for N_I = 2 is a clear evidence of impurity–impurity induced interactions.

The overall behavior of the excitation spectrum $A({\omega }_{f};{g}_{B\uparrow })$ for N_I = 2 and g_II = 0 is shown in figure 4(g) with varying ${g}_{B\uparrow }$. Evidently, the position of the dominant quasiparticle peak in terms of ω_f increases almost linearly for larger ${g}_{B\uparrow }$. This behavior essentially reflects the linear increase of the energy of the initial state $| {\rm{\Psi }}(0)\rangle$ (equation (2)) directly after the quench. Moreover, comparing the position of the dominant quasiparticle peak with ${{\rm{\Delta }}}_{+}^{{N}_{I}=1}$ reveals that for ${g}_{B\uparrow }\,\gt \,0.5$, while the latter saturates, the former increases and additional peaks appear in the spectrum $A({\omega }_{f};{g}_{B\uparrow })$. These peaks correspond to excited states of the system and already for ${g}_{B\uparrow }\,\gt \,1$ the ground states corresponding to ${{\rm{\Delta }}}_{+}^{{N}_{I}}$ cease to be populated during the dynamics. In a similar fashion, such additional quasiparticle peaks occur also for attractive interactions, see figure 4(g) for ${g}_{B\uparrow }\,\lt \,-0.5$. In this case the additional quasiparticle peaks stem from the induced interactions resulting in the presence of a peak at ${\omega }_{f}=2{{\rm{\Delta }}}_{+}^{{N}_{I}=2}-{{\rm{\Delta }}}_{+}^{{N}_{I}=1}\ne {{\rm{\Delta }}}_{+}^{{N}_{I}=1}$ and other ones which correspond to the occupation of higher-lying excited polaronic states with S_z = 1(equation (2)). Note that such an almost linear behavior of the polaronic spectrum is reminiscent of the corresponding three-dimensional scenario but away from the Feshbach resonance regime. The latter corresponds in one-dimension to an interspecies Tonks–Girardeau interaction regime which is not addressed in the present work. We remark that in one-dimension there is no molecular bound state occurring for repulsive interactions.

Summarizing, we can infer that the quasiparticle excitation spectrum depends strongly on the value of the postquench interspecies interaction strength and also on the number of impurities outside the weakly attractive and repulsive coupling regimes [98]. However, this behavior is also slightly altered when going from two non-interacting to two weakly interacting impurities. For a relevant discussion on the lifetime of the above-described spectral features we refer the interested reader to [99]. It is also important to mention that in the weakly interacting impurity-BEC regime where the contrast is finite in the course of the evolution the spectral function A(ω_f) corresponds to the injection spectrum in the framework of the reverse rf spectroscopy [2, 26].

Below we further analyze the dynamical response of the multicomponent system, and especially of the impurities, following an interspecies interaction quench from ${g}_{B\uparrow }\,=\,0$ to ${g}_{B\uparrow }\,\gt \,0$ within the above identified dynamical regions of the contrast. In particular, we explore the dynamics of the system on both the single- and the two-body level and further develop an effective potential picture to provide a more concrete interpretation of the emergent phenomena. We mainly focus on the nonequilibrium dynamics of two non-interacting impurities (g_II = 0) and subsequently discuss whether possible alterations might occur for weakly interacting (g_II = 0.2) impurities. Also, in the following, only the temporal-evolution of the pseudospin-$\uparrow$ part of the impurities is discussed since the pseudospin-$\downarrow$ component does not interact with the bosonic medium.

4.4.1. Density evolution and effective potential

To visualize the spatially resolved dynamics of the system on the single-particle level we first inspect the time-evolution of the σ-species single-particle density ${\rho }_{\sigma }^{(1)}(x;t)$ (equation (6)) illustrated in figure 5. For weak postquench interspecies repulsions lying within the region R_II e.g. ${g}_{B\uparrow }\,=\,0.25$, such that ${g}_{B\uparrow }\,\lt \,{g}_{{BB}}$, the impurities (see figure 5(b)) exhibit a breathing motion of frequency ${\omega }_{{\rm{br}}}^{I}\approx 1.44$ inside the bosonic medium [73, 74]. Moreover, at initial evolution times (t < 60) the amplitude of the breathing is almost constant whilst later on (t > 60) it shows a slightly decaying tendency, see for instance the smaller height of the density peak at t = 70 compared to t = 20 in figure 5(b). This decaying amplitude can be attributed to the build up of impurity–impurity correlations in the course of the evolution [42] due to the presence of induced interactions discussed later on, see also figure 7(a). The breathing motion of the impurities is directly captured by the periodic contraction and expansion in the shape of the instantaneous density profiles of ${\rho }_{\uparrow }^{(1)}(x;t)$ depicted in figure 6(b). On the other hand, the bosonic gas remains essentially unperturbed (figure 5(a)) throughout the dynamics, showing only tiny distortions from its original Thomas–Fermi cloud due to its interaction with the impurity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

An intuitive understanding of the observed dynamics of the impurities is provided with the aid of an effective potential picture. Indeed, the impurity-BEC interactions can be taken into account, to a very good approximation, by employing a modified external potential for the impurities. The latter corresponds to the time-averaged effective potential created by the harmonic oscillator and the density of the bosonic gas [35, 51, 73, 75] namely

$$\begin{eqnarray}&&{\bar{V}}_{I}^{{\rm{eff}}}(x)=\displaystyle \frac{1}{2}m{\omega }^{2}{x}^{2}+\displaystyle \frac{{g}_{{BI}}}{T}{\int }_{0}^{T}{\rm{d}}{t}{\rho }_{B}^{(1)}(x;t).\end{eqnarray} \tag{ 15 }$$

The averaging process aims to eliminate the emergent very weak distortions on the instantaneous density of the BEC ${\rho }_{B}^{(1)}(x;t)$, and it is performed herein over T = 100. These distortions being a consequence of the motion of the impurities within the BEC are imprinted as a slow and very weak amplitude breathing motion of ${\rho }_{B}^{(1)}(x;t)$ with ${\omega }_{{\rm{br}}}^{B}\approx 1.82$, hardly visible in figure 5(a). They are canceled out in our case for T > 20. Note that ${\omega }_{{\rm{br}}}^{B}\lt 2$ is attributed to the repulsive character of the BEC background which negatively shifts its breathing frequency from the corresponding non-interacting value [100]. At ${g}_{B\uparrow }\,=\,0.25$ this ${\bar{V}}_{I}^{{\rm{eff}}}(x)$ takes the form of a modified harmonic oscillator potential illustrated in figure 6(a) together with the densities of its first few single-particle eigenstates. Furthermore, assuming the Thomas–Fermi approximation for ${\rho }_{B}^{(1)}(x,t)$ the effective trapping frequency of the impurities corresponds to ${\omega }_{{\rm{eff}}}=\omega \sqrt{1-\tfrac{{g}_{B\uparrow }}{{g}_{{BB}}}}$. Therefore their expected effective breathing frequency would be ${\omega }_{{\rm{br}}}^{{\rm{eff}},I}=2{\omega }_{{\rm{eff}}}\approx 1.41$ which is indeed in a very good agreement with the numerically obtained ${\omega }_{{\rm{br}}}^{I}$. The discrepancy between the prediction of the effective potential and the MB approach is attributed to the approximate character of the effective potential which does not account for possible correlation induced shifts to the breathing frequency. Moreover, in the present case the impurities which undergo a breathing motion within ${\bar{V}}_{I}^{{\rm{eff}}}(x)$ reside predominantly in its energetically lowest-lying state E₁, see figure 6(a). It is also important to mention that this effective potential approximation is adequate only for weak interspecies interactions where the impurity-BEC entanglement is small [35, 75]. Note also that the inclusion of the Thomas–Fermi approximation in the effective potential of equation (15) can not adequately describe the impurities dynamics when they reach the edges of the bosonic cloud, see [36] for more details. However in this case ${\rho }_{\uparrow }^{(1)}(x;t)$ lies within ${\rho }_{B}^{(1)}(x;t)$ throughout the evolution indicating the miscible character of the dynamics for ${g}_{B\uparrow }\,\lt \,{g}_{{BB}}$ [35, 65]. Furthermore, for these weak postquench interspecies repulsions a similar to the above-described dynamics takes place also for two weakly (g_II = 0.2) repulsively interacting impurities as shown in figure 5(c). The impurities undergo a breathing motion within the bosonic medium in the course of the time-evolution exhibiting a slightly larger oscillation frequency than for the g_II = 0 case but with the same amplitude (hardly visible by comparing figures 5(b) and (c)).

For larger postquench interaction strengths ${g}_{B\uparrow }\,=\,0.5$ (region R_III), i.e. close to the intraspecies interaction of the bosonic bath g_BB, the impurities show a more complex dynamics compared to the weak interspecies repulsive case (figure 5(e)). Also, the BEC medium performs a larger amplitude breathing motion (figure 5(d)) compared to the ${g}_{B\uparrow }\,=\,0.25$ scenario but again with a frequency ${\omega }_{{\rm{br}}}^{{B}}\approx 1.82$. Focusing on the impurities motion, we observe that at short evolution times (0 < t < 5) after the quench ${\rho }_{\uparrow }^{(1)}(x;t)$ expands and then splits into two counterpropagating density branches with finite momenta that travel towards the edges of the bosonic cloud, see figure 5(e) and the profiles shown in figure 6(d). The appearance of these counterpropagating density branches is a consequence of the interaction quench which imports energy into the system. Reaching the edges of ${\rho }_{B}^{(1)}(x;t)$ the density humps of ${\rho }_{\uparrow }^{(1)}(x;t)$ are reflected back towards the trap center (x = 0) where they collide around t ≈ 15 forming a single density peak (figure 6(d)). The aforementioned impurity motion repeats itself in a periodic manner for all evolution times (figure 5(e)). Here, the underlying time-averaged effective potential (equation (15)) corresponds to a highly deformed harmonic oscillator possessing an almost square-well like profile as illustrated in figure 6(c). Moreover, a direct comparison of the densities of the lower-lying single-particle eigenstates of ${\bar{V}}_{I}^{{\rm{eff}}}(x)$ (figure 6(c)) with the density profile snapshots of ${\rho }_{\uparrow }^{(1)}(x;t)$ of the MB dynamics (figure 6(d)) reveals that the impurities predominantly reside in a superposition of the two lower-lying excited states (E₁ and E₂) of ${\bar{V}}_{I}^{{\rm{eff}}}(x)$. Additionally in the case of two weakly repulsively interacting impurities, shown in figure 5(f), the impurities' motion remains qualitatively the same. However, due to the inclusion of intraspecies repulsion the impurities possess a slightly larger overall oscillation frequency and the collisional patterns at the trap center appear to be modified as compared to the g_II = 0 case.

Turning to strong postquench repulsions, i.e. ${g}_{B\uparrow }\,=\,1.5\gg {g}_{{BB}}$ which belongs to R_IV, the dynamical response of the impurities is greatly altered and the bosonic gas exhibits an enhanced breathing dynamics as compared to the weak and intermediate interspecies repulsions discussed above. Initially ${\rho }_{\uparrow }^{(1)}(x;t=0)$ consists of a density hump located at the trap center which, following the interaction quench, breaks into two density fragments, as illustrated in figure 5(h), each of them exhibiting a multihump structure (see also figure 6(f)). Note that the density hump at the trap center remains the dominant contribution of ${\rho }_{\uparrow }^{(1)}(x;t)$ until it eventually fades out for t > 5, see figure 5(h). This multihump structure building upon ${\rho }_{\uparrow }^{(1)}(x;t)$ is clearly captured in the instantaneous density profiles depicted in figure 6(f). Remarkably, the emergent impurity density fragments that are symmetrically placed around the trap center (x = 0) perform a damped oscillatory motion in time around the edges of the Thomas–Fermi radius of the bosonic gas, see in particular figures 5(g), (h).

The emergent dynamics of the impurities can also be interpreted to lowest order approximation (i.e. excluding correlation effects) by invoking the corresponding effective potential which for these strong interspecies repulsions has the form of the double-well potential shown in figure 6(e). Comparing the shape of the densities of the eigenstates of ${\bar{V}}_{I}^{{\rm{eff}}}(x)$ (figure 6(e)) with the density profiles ${\rho }_{\uparrow }^{(1)}(x;t)$ (figure 6(f)) obtained within the MB dynamics simulations it becomes evident that the impurities reside in a superposition of higher-lying states of the effective potential. Furthermore the double-well structure of ${\bar{V}}_{I}^{{\rm{eff}}}(x)$ suggests that each of the observed density fragments of the impurities is essentially trapped in each of the corresponding two sites of ${\bar{V}}_{I}^{{\rm{eff}}}(x)$. Of course, as already mentioned above, for these strong interactions ${\bar{V}}_{I}^{{\rm{eff}}}(x)$ provides only a crude description of the impurity dynamics since it does not account for both intra- and interspecies correlations that occur during the MB dynamics. However ${\bar{V}}_{I}^{{\rm{eff}}}(x)$ enables the following intuitive picture for the impurity dynamics. Namely, the damped oscillations of ${\rho }_{\uparrow }^{(1)}(x;t)$ designate that the pseudospin-$\uparrow$ impurities at initial times are in a superposition state of a multitude of highly excited states (see e.g. figure 6(f) at t = 8) while for later times they reside in a superposition of lower excited states (see e.g. figure 6(f) at t = 15). We should also remark that a similar overall dynamical behavior on the single-particle level has been reported in the case of a single spinor impurity and has been also related to an enhanced energy transfer from the impurity to the bosonic bath [35, 68, 69, 75]. Such an energy transfer process takes place also in the present case (results not shown here). Another important feature of the observed dynamical response of the impurities is the fact that they are not significantly affected by the presence of weak intraspecies interactions. This can be seen by inspecting figure 5(i) which shows the time-evolution of ${\rho }_{\uparrow }^{(1)}(x;t)$ for g_II = 0.2. Here, the most noticeable difference when compared to the g_II = 0 scenario is that the splitting of ${\rho }_{\uparrow }^{(1)}(x;t)$ into two branches occurs at shorter time scales (compare figures 5(h), (i)) due to the additional intraspecies repulsion.

4.4.2. Dynamics of the two-body reduced density matrix

To investigate the development of impurity–impurity correlations during the quench dynamics we next resort to the time-evolution of the pseudospin-$\uparrow$ impurity intraspecies two-body reduced matrix ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ (equation (7)). Recall that ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ provides the probability of finding at time t a pseudospin-$\uparrow$ boson at location x₁ and a second one at x₂ [65, 66]. Most importantly, it allows us to monitor the two-body spatially resolved dynamics of the impurities and infer whether they move independently or correlate with each other [8, 10, 42].

Figure 7 shows ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ at specific time-instants of the evolution of two non-interacting (figures 7(a₁)–(b₆) and (d₁)–(d₆)) as well as weakly interacting (figures 7(c₁)–(c₆)) impurities for different postquench interspecies repulsions. To reveal the role of induced impurity–impurity correlations via the bath we mainly focus on two initially non-interacting impurities where ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t=0)={\rho }_{\uparrow }^{(1)}({x}_{1},t=0){\rho }_{\uparrow }^{(1)}({x}_{2},t=0)/2$ since g_II = 0 and initially ${g}_{B\uparrow }\,=\,0$. As already discussed in section 4.4.1 for weak interspecies postquench repulsions, namely ${g}_{B\uparrow }\,=\,0.25$ (region R_II), the impurities perform a breathing motion on the single-particle level (figure 5(b)) exhibiting a decaying amplitude for large evolution times. Accordingly, inspecting ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ (figures 7(a₁)–(a₆)) we observe that the impurities are likely to reside together close to the trap center since ${\rho }_{\uparrow \uparrow }^{(2)}(-2\lt {x}_{1}\lt 2,-2\lt {x}_{2}\lt 2;t)$ is mainly populated throughout the evolution. In particular, at initial times ${\rho }_{\uparrow \uparrow }^{(2)}(-2\lt {x}_{1}\lt 2,-2\lt {x}_{2}\lt 2;t)$ shows a Gaussian-like distribution which contracts (figure 7(a₂)) and expands (figures 7(a₃), (a₄)) during the dynamics as a consequence of the aforementioned breathing motion. Deeper in the evolution ${\rho }_{\uparrow }^{(1)}(x;t)$ decays and ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ is deformed along its diagonal (figures 7(a₄), (a₆)) or its anti-diagonal (figure 7(a₅)) indicating that the impurities tend to be slightly apart or at the same location respectively. This is indicative of the admittedly weak induced interactions as the breathing mode along the anti-diagonal of ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ (relative coordinate breathing mode) does not possess exactly the same frequency as the breathing along the diagonal (center-of-mass breathing mode).

Zoom In Zoom Out Reset image size

For larger interspecies repulsions e.g. for ${g}_{B\uparrow }\,=\,0.5$ (region R_III) the two-body dynamics of the impurities is significantly altered, see figures 7(b₁)–(b₆). At the initial stages of the dynamics the impurities reside together in the vicinity of the trap center as ${\rho }_{\uparrow \uparrow }^{(2)}(-3\lt {x}_{1}\lt 3,-3\lt {x}_{2}\lt 3;t)$ is predominantly populated. However for later times two different correlation patterns appear in ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ in a periodic manner. Recall that for these interactions ${\rho }_{\uparrow }^{(1)}(x;t)$ splits into two counterpropagating density branches traveling towards the edges of the bosonic bath and then are reflected back to the trap center where they collide (figure 5(e)). Consequently, when the two density fragments appear in ${\rho }_{\uparrow }^{(1)}(x;t)$ the impurities reside in two different two-body configurations (figures 7(b₂), (b₄) and (b₆)). Namely the bosonic impurities either lie together at a certain density branch (see the diagonal elements of ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t))$ or they remain spatially separated with one of them residing in the left and the other in the right density branch (see the anti-diagonal elements of ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$). Moreover, during their collision at x = 0 the impurities are very close to each other as it is evident by the enhanced two-body probability in the neighborhood of x₁ = x₂ = 0 (figures 7(b₃), (b₅)). The dynamics of two weakly repulsive (g_II = 0.2) impurities shows similar two-body correlation patterns to the non-interacting ones, as it can be seen by comparing figures 7 (b₁)–(b₆) to (c₁)–(c₆). This behavior complements the similarities already found at the single-particle level (see section 4.4.1). The major difference on the two-body level between the g_II = 0.2 and g_II = 0 scenario is that in the former case ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ is more elongated along its anti-diagonal when the impurities collide at x = 0 (figures 7(c₁), (c₃)). Therefore weakly interacting impurities tend to be further apart compared to the g_II = 0 case, a result that reflects their direct repulsion. Other differences observed at the same time-instant in ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ between the interacting and the non-interacting cases are due to the repulsive s-wave interaction that directly competes with the attractive induced interactions emanating in the system. For instance, shortly after a collision point e.g. at t = 55, shown in figures 7(b₅) and (c₅), we observe that due to the repulsive s-wave interactions the attractive contribution between the impurities, see the diagonal of ${\rho }_{\uparrow \uparrow }^{(2)}(-2\lt {x}_{1}\lt 2,-2\lt {x}_{2}\lt 2;t)$ in figure 7(b₅) disappears (figure 7(c₅)).

Turning to very strong repulsions, e.g. for ${g}_{B\uparrow }\,=\,1.5$ lying in region R_IV, the correlation patterns of the two non-interacting impurities (figures 7(d₁)–(d₆)) show completely different characteristics compared to the ${g}_{B\uparrow }\,\leqslant \,{g}_{{BB}}$ regime. Note here that in the dynamics of ${\rho }_{\uparrow }^{(1)}(x;t)$ the initially formed density hump breaks into two density fragments (figure 5(h)) possessing a multihump shape (see also figure 6(f)). Subsequently, the fragments lying symmetrically with respect to x = 0 perform a damped oscillatory motion in time residing around the edges of the Thomas–Fermi radius of the bosonic gas. The corresponding two-body reduced density matrix shows a pronounced probability peak around x₁ = x₂ = 0 (figure 7(d₁)) indicating that at the initial stages of the dynamics the impurities are mainly placed together in this location. As time evolves, the impurities predominantly move as a pair towards the edge of the Thomas–Fermi background, see in particular the diagonal of ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ in figures 7(d₂), (d₃), and simultaneously they start to exhibit a delocalized behavior as can be deduced by the small values of the off-diagonal elements of ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2}\ne {x}_{1};t)$. Entering deeper in the evolution the aforementioned delocalization of the impurities becomes more enhanced since ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ disperses as illustrated in figures 7(d₄)–(d₆). This dispersive behavior of ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ is inherently related to the multihump structure of ${\rho }_{\uparrow }^{(1)}(x;t)$ and suggests from a two-body perspective the involvement of several excited states during the impurity dynamics. It is also worth mentioning that at specific time instants the diagonal of ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ is predominantly populated (figures 7(d₂), (d₃), (d₅)) which is indicative of the presence of induced interactions.

4.4.3. Two-body dynamics within the effective potential picture

To further expose the necessity of taking into account the intra- and the interspecies correlations of the system in order to accurately describe the MB dynamics of the impurities we next solve the time-dependent Schrödinger equation that governs the system's dynamics relying on the previously introduced effective potential picture (equation (15)) via exact diagonalization⁶ . Thus our main aim here is to test the validity of ${\bar{V}}_{I}^{{\rm{eff}}}(x)$ at least to qualitatively capture the basic features of the emergent non-equilibrium dynamics of the two impurities. We emphasize again that ${\bar{V}}_{I}^{{\rm{eff}}}$ does not include any interspecies correlation effects that arise in the course of the temporal-evolution of the impurities. Within this approximation the effective Hamiltonian that captures the quench-induced dynamics of the impurities reads

$$\begin{eqnarray}&&{H}^{{\rm{eff}}}=\displaystyle \int {\rm{d}}{x}\,{\hat{{\rm{\Psi }}}}_{\uparrow }^{\dagger }(x)\left(-\displaystyle \frac{{{\hslash }}^{2}}{2m}\displaystyle \frac{{{\rm{d}}}^{2}}{{{\rm{d}}{x}}^{2}}+{\bar{V}}_{I}^{{\rm{eff}}}\right){\hat{{\rm{\Psi }}}}_{\uparrow }(x)+{g}_{\uparrow \uparrow }\displaystyle \int {\rm{d}}{x}{\hat{{\rm{\Psi }}}}_{\uparrow }^{\dagger }(x){\hat{{\rm{\Psi }}}}_{\uparrow }^{\dagger }(x){\hat{{\rm{\Psi }}}}_{\uparrow }(x){\hat{{\rm{\Psi }}}}_{\uparrow }(x),\end{eqnarray} \tag{ 16 }$$

where ${\hat{{\rm{\Psi }}}}_{\uparrow }(x)$ is the bosonic field-operator of the pseudospin-$\uparrow$ impurity and ${g}_{\uparrow \uparrow }$ denotes the intraspecies interactions between the two pseudospin-$\uparrow$ impurity atoms. Recall that the intercomponent contact interaction of strength ${g}_{B\uparrow }$ and the intraspecies interaction between the bath atoms are inherently embedded into ${\bar{V}}_{I}^{{\rm{eff}}}$ (equation (15)). In particular, within ${\bar{V}}_{I}^{{\rm{eff}}}$ we account for the correlated Thomas–Fermi profile of the BEC since ${\rho }_{B}^{(1)}(x;t)$ is determined from the MB approach. Below, we exemplarily study the dynamics of two non-interacting impurities and therefore we set ${g}_{\uparrow \uparrow }\,=\,0$ in equation (16). Moreover, in order to trigger the non-equilibrium dynamics we consider an interspecies interaction quench from ${g}_{B\uparrow }\,=\,0$ (t = 0) to a finite repulsive value of ${g}_{B\uparrow }$. Such a sudden change is essentially taken into account via a deformation of ${\bar{V}}_{I}^{{\rm{eff}}}$ (equation (15)).

The corresponding instantaneous two-body reduced density matrix of the impurities within H^eff is depicted in figure 8 for distinct values of ${g}_{B\uparrow }$. Focusing on weak postquench interactions, e.g. ${g}_{B\uparrow }\,=\,0.25$, we observe that at the initial times the two-body dynamics of the impurities is adequately described within H^eff (compare figures 7(a₁)–(a₃) to figures 8(a₁)–(a₃)). Indeed, in this time-interval only some minor deviations between the heights of the peaks of ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ obtained within the MB and the H^eff approach are observed. However, for longer times H^eff (figures 8(a₄)–(a₆)) fails to capture the correct shape of ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ and more precisely its deformations occurring along its diagonal or anti-diagonal (see figures 8(a₄)–(a₆)) which stem from the build up of higher-order correlations during the dynamics.

Zoom In Zoom Out Reset image size

Increasing the repulsion such that ${g}_{B\uparrow }\,=\,0.5$, deviations between the effective potential approximation and the correlated approach become more severe. For instance, at the initial times the sharp two-body probability peak of ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ in the vicinity of x₁ = x₂ = 0 arising in the MB dynamics (figure 7(b₁)) becomes smoother within H^eff (figure 8(b₁)) although the overall shape of ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ remains qualitatively similar. Moreover, the observed elongations along the diagonal of ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ exhibited due to the presence of correlations are not captured in the effective picture, e.g. compare figures 7(b₃), (b₅) with figures 8(b₃), (b₅). Remarkably, the two-body superposition identified in ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ of two different two-body configurations occurring at specific time-instants is also predicted at least qualitatively via H^eff, see figures 8 (b₂), (b₄) and (b₆). We remark that the differences in the patterns of ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ between H^eff and the correlated approach are even more pronounced when g_II = 0.2 (results not shown).

Strikingly for strongly repulsive interactions, ${g}_{B\uparrow }\,=\,1.5$, H^eff completely fails to capture the two-body dynamics of the impurities. This fact can be directly inferred by comparing ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ within the two approaches, see figures 7(c₁)–(c₆) and figures 8(c₁)–(c₆). Even at the initial stages of the dynamics the effective potential cannot adequately reproduce the correct shape of ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$, compare figure 8(c₁) with figure 7(d₁). Note, for instance, the absence of the central two-body probability peak in the region −2 < x₁, x₂ < 2 within H^eff which demonstrates the correlated character of the dynamics. More precisely, ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ obtained via H^eff shows predominantly the development of two different two-body configurations. The first pattern suggests that the impurities either reside together at the same edge of the BEC background or each one is located at a distinct edge of the Thomas–Fermi profile, see e.g. figures 8(c₁), (c₅). However, at different time-instants ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ indicates that the impurities lie in the vicinity of the trap center as illustrated e.g. in figures 8(c₂), (c₄) and (c₆), an event that never occurs for t > 5 in the MB dynamics (see figure 5(h)). It is also worth mentioning that the observed dispersive character of ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ in the MB dynamics (see e.g. figures 7(d₄)–(d₆)) is a pure correlation effect and a consequence of the participation of a multitude of excited states in the impurity dynamics which is never captured within H^eff.

Next we discuss the dynamical behavior of both the BEC medium and the bosonic impurities on both the one- and the two-body level after an interspecies interaction quench from ${g}_{B\uparrow }\,=\,0$ to the attractive regime of ${g}_{B\uparrow }\,\lt \,0$. To explain basic characteristics of the dynamics of the impurities an effective potential picture is also employed. As in the previous section we first examine the emergent time-evolution of two non-interacting impurities (g_II = 0) and then compare our findings to that of two weakly interacting (g_II = 0.2) ones.

4.5.1. Single-particle dynamics and effective potential

To investigate the spatially resolved dynamics of the multicomponent system after an interaction quench from ${g}_{B\uparrow }\,=\,0$ to ${g}_{B\uparrow }\,\lt \,0$, we first analyze the spatio-temporal evolution of the σ-species single-particle density ${\rho }_{\sigma }^{(1)}(x;t)$. The dynamical response of ${\rho }_{\sigma }^{(1)}(x;t)$ triggered by the quench is presented in figure 9 for postquench interspecies attractions ${g}_{B\uparrow }\,=\,-0.5$ (figures 9(a)–(c)) and ${g}_{B\uparrow }\,=\,-1$ (figures 9(d)–(f)).

Zoom In Zoom Out Reset image size

Inspecting the dynamics of two non-interacting impurities at ${g}_{B\uparrow }\,=\,-0.5$ (region ${R}_{{III}}^{{\prime} }$), shown in figures 9(a), (b), we deduce that ${\rho }_{\uparrow }^{(1)}(x;t)$ undergoes a breathing motion inside ${\rho }_{B}^{(1)}(x;t)$ characterized by a predominant frequency ${\omega }_{{br}}^{I}\approx 2.76$ and a secondary one ${\omega }_{{\rm{br}}}^{{\prime} I}\approx 2.88$ thus producing a beating pattern. These two distinct frequencies stem from the center-of-mass and relative coordinate breathing modes of the impurities, whose existence originates from the presence of attractive induced interactions in the system. We remark that the breathing frequency of the center-of-mass can be estimated in terms of the corresponding effective potential of the impurities, see also equation (17). In particular for ${g}_{B\uparrow }\,=\,-0.5$, ${\omega }_{{\rm{br}}}^{I}=2\sqrt{2.06}\approx 2.87$ (see also the comment in⁷ )which is in very good agreement with ${\omega }_{{\rm{br}}}^{{\prime} I}$. The relevant contraction of ${\rho }_{\uparrow }^{(1)}(x;t)$ can be inferred by its increasing amplitude that takes place from the very early stages of the non-equilibrium dynamics (figure 10(b)). The beating pattern can be readily identified e.g. by comparing the maximum height of ${\rho }_{\uparrow }^{(1)}(x;t)$ during its contraction at initial and later stages of the dynamics, see e.g. ${\rho }_{\uparrow }^{(1)}(x;t)$ at t = 10 and t = 40 in figure 9(b). Moreover, as a consequence of the motion of the impurity and the relatively weak interspecies attraction, i.e. ${g}_{B\uparrow }\,=\,-0.5$, the Thomas–Fermi cloud of the bosonic gas becomes slightly distorted. In particular, a low amplitude density hump is imprinted on ${\rho }_{B}^{(1)}(x;t)$ exactly at the position of ${\rho }_{\uparrow }^{(1)}(x;t)$ as shown by the white colored region in figure 9(a) in the vicinity of x = 0 [75]. An almost similar effect to the above-mentioned breathing dynamics is present also for the case of two weakly interacting impurities (figure 9(c)). Here, the secondary frequency manifests itself at later evolution times resulting in turn in a slower beating of ${\rho }_{\uparrow }^{(1)}(x;t)$ compared to the g_II = 0 scenario (hardly visible in figure 9(c)). This delayed occurrence is attributed to the presence of intraspecies repulsion which competes with the attractive induced interactions.

Zoom In Zoom Out Reset image size

For a larger negatively valued interspecies coupling, e.g. for g_BI = −1 within region ${R}_{{III}}^{{\prime} }$, ${\rho }_{\uparrow }^{(1)}(x;t)$ becomes more spatially localized and again performs a decaying amplitude breathing motion, the so-called beating identified above, but with a larger major frequency, ${\omega }_{{\rm{br}}}^{I}\approx 3.2$, compared to the ${g}_{B\uparrow }\,=\,-0.5$ case (figure 9(e)). Notice that the observed beating motion of the impurities persists while being more dramatic for this stronger attraction (compare figures 9(b) and (e)). This enhanced attenuation of the breathing amplitude together with the strong localization of the impurities is a direct effect of the dominant presence of interspecies attractions between the impurity and the bath, see also [75]. Also, due to the stronger ${g}_{B\uparrow }$ and the increased spatial localization of ${\rho }_{\uparrow }^{(1)}(x;t)$, the density hump building upon ${\rho }_{B}^{(1)}(x;t)$ at the instantaneous position of the impurities is much more pronounced than that found for ${g}_{B\uparrow }\,=\,-0.5$ (figure 9(d)). Note that the density hump appearing in ${\rho }_{B}^{(1)}(x;t)$ is essentially an imprint of the impurities presence and motion within the bosonic medium. Indeed, ${\rho }_{\uparrow }^{(1)}(x;t)$ exhibits a sech-like form tending to be more localized for a larger interspecies attractions ${g}_{B\uparrow }$, see e.g. ${\rho }_{\uparrow }^{(1)}(x;t)$ at a fixed time-instant for ${g}_{B\uparrow }\,=\,-0.5$ and ${g}_{B\uparrow }\,=\,-1$ in figures 9(b) and (e) respectively, a behavior that also holds for the consequent density hump in ${\rho }_{B}^{(1)}(x;t)$ (figures 9(a), (d)). We should remark that for large negative ${g}_{B\uparrow }$ the system becomes strongly correlated and the BEC is highly excited. The latter is manifested by the development of an overall weak amplitude breathing motion of the bosonic gas, see figure 9(d). Furthermore, the inclusion of weak intraspecies repulsions between the impurities does not significantly alter their dynamics (figure 9(f)). Indeed, a faint increase of their expansion magnitude takes place and the corresponding amplitude of the beating decays faster (compare figures 9(d) and (f)).

The above-mentioned dynamics can also be qualitatively explained in terms of a corresponding effective potential approximation [35, 73, 75]. Yet again, the effective potential experienced by the impurities consists of the external harmonic oscillator V(x) and the single-particle density of the BEC background. Importantly, since ${\rho }_{B}^{(1)}(x;t)$ is greatly distorted from its original Thomas–Fermi profile due to the motion of the impurities, we invoke a time-averaged effective potential. Consequently, the effective potential of the impurity reads

$$\begin{eqnarray}&&{\bar{V}}_{I}^{{\rm{eff}}}(x)=V(x)-\displaystyle \frac{\left|{g}_{{BI}}\right|}{T}{\int }_{0}^{T}{\rm{d}}{t}{\rho }_{B}^{(1)}(x;t),\end{eqnarray} \tag{ 17 }$$

where T = 100 denotes the corresponding total propagation time. We remark that for the considered negative values of ${g}_{B\uparrow }$ the shape of ${\bar{V}}_{I}^{{\rm{eff}}}(x)$ does not significantly change after averaging over T = 60. A schematic illustration of ${\bar{V}}_{I}^{{\rm{eff}}}(x)$ and the densities of its first few single-particle eigenstates at ${g}_{B\uparrow }\,=\,-1$ is presented in figure 10(a), see also remark (see footnote 7). The observed localization tendency of ${\rho }_{\uparrow }^{(1)}(x;t)$ around the aforementioned potential minimum is essentially determined by the strongly attractive behavior of ${\bar{V}}_{I}^{{\rm{eff}}}(x)$. Remarkably, the distinct dynamical features of the impurities for an increasing interspecies attraction can be partly understood with the aid of ${\bar{V}}_{I}^{{\rm{eff}}}(x)$. Indeed, for increasing $\left|{g}_{{BI}}\right|$ the effective frequency of ${\bar{V}}_{I}^{{\rm{eff}}}(x)$ is larger and ${\bar{V}}_{I}^{{\rm{eff}}}(x)$ becomes more attractive. The former property of ${\bar{V}}_{I}^{{\rm{eff}}}(x)$ accounts for the increasing breathing frequency of the impurity wavepacket for larger $\left|{g}_{{BI}}\right|$. Additionally, the increasing attractiveness of ${\bar{V}}_{I}^{{\rm{eff}}}(x)$ is responsible for the reduced width of ${\rho }_{\uparrow }^{(1)}(x;t)$ for a larger $\left|{g}_{{BI}}\right|$ and thus its increasing localization tendency.

4.5.2. Two-body correlation dynamics and comparison to the effective potential approximation

Having described the time-evolution of the impurities on the single-particle level, we next analyze the dynamical response of the pseudospin-$\uparrow$ component by invoking the corresponding two-body reduced density matrix ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ (see also equation (7)).

The time-evolution of ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ is depicted in figures 11(a₁)–(a₅) for two non-interacting (g_II = 0) impurities following an interspecies interaction quench from ${g}_{B\uparrow }\,=\,0$ to ${g}_{B\uparrow }\,=\,-0.5$ (region ${R}_{{III}}^{{\prime} }$). Before the quench the impurities lie together in the vicinity of the trap center since ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1}=0,{x}_{2}=0;t=0)$ shows a high probability peak (figure 11(a₁)). However as time evolves the two bosons start to occupy a relatively smaller spatial region as can be deduced by the shrinking of the central two-body probability peak across the diagonal at t = 10 in figure 11(a₂). Then they move either opposite to each other (see the elongated anti-diagonal in figures 11(a₃), (a₅)) or tend to bunch together at the same location (see the pronounced diagonal of ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2}={x}_{1};t=60)$ in figure 11(a₄)). This latter behavior of the impurities is the two-body analog of their wavepacket periodic expansion and contraction (relative coordinate breathing motion) discussed previously on the single-particle level (figure 9(b)).

Zoom In Zoom Out Reset image size

The dynamics of two weakly repulsively interacting (g_II = 0.2) impurities (figures 11 (b₁)–(b₅)) shows similar characteristics to the above-described non-interacting scenario. Indeed, initially (figure 11(b₁)) and at short times (figure 11(b₂)) the impurities reside close to the trap center while later on they repel (see e.g. figure 11(b₃)) or attract (figure 11(b₄)) each other as a result of their breathing dynamics (see also figure 9(c)). The major difference between the weakly interacting and the non-interacting impurities is that their distance which is given by the anti-diagonal distribution of their two-body reduced density matrix is slightly different, see figures 11(e₁)–(e₅). For instance at t = 40 the non-interacting impurities are further apart from each other as compared to the case of interacting impurities, while this situation is reversed at t = 90. The aforementioned difference owes its existence to the distinct relative coordinate breathing frequencies. This can be directly inferred from the fact that ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ possesses a larger spatial distribution when g_II = 0.2 and it is attributed to their underlying mutual repulsion. For instance, even initially ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t=0)$ for g_II = 0.2 (figure 11(b₁)) is slightly deformed towards its anti-diagonal compared to the g_II = 0 case (figure 11(a₁)). This behavior persists also during the evolution independently of the expansion or the contraction of the impurity cloud, as can be seen by comparing figures 11(b₄) to (a₄) and figures 11(b₅) to (a₅).

To reveal the importance of both intra- and interspecies correlations for the impurity dynamics we then utilize the effective potential, ${\bar{V}}_{I}^{{\rm{eff}}}(x)$, introduced in equation (17) and solve numerically the time-dependent Schrödinger equation of the impurities via exact diagonalization. We remark once more that ${\bar{V}}_{I}^{{\rm{eff}}}$ neglects the interspecies correlations of the multicomponent system but includes the density profile of the BEC determined by the MB approach. In particular, we construct the effective Hamiltonian H^eff of equation (16) but using the ${\bar{V}}_{I}^{{\rm{eff}}}(x)$ of equation (17). For brevity we focus on the case of ${g}_{\uparrow \uparrow }\,=\,0.2$ and analyze the dynamics after an interspecies interaction quench from ${g}_{B\uparrow }\,=\,0$ (t = 0) to ${g}_{B\uparrow }\,=\,-0.5$. As explained in section 4.4.3 within the effective potential picture this quench scenario accounts for the deformation of ${\bar{V}}_{I}^{{\rm{eff}}}$. Snapshots of ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ when g_II = 0.2 and ${g}_{B\uparrow }\,=\,-0.5$ obtained within H^eff are illustrated in figures 11(c₁)–(c₅). As it can be seen by comparing ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ for the MB approach (figures 11(b₁)–(b₅)) and H^eff (figures 11(c₁)–(c₅)) significant deviations occur between the two methods. Indeed, during the time-evolution the correlation patterns visible in ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ calculated via H^eff exhibit similar overall characteristics to the ones taking place in the correlated approach but at completely different time-scales. In fact, ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ shows elongated shapes along its diagonal (figure 11(c₃)) or anti-diagonal (figure 11(c₄)) implying that the impurities tend to be relatively close or apart from one another respectively. The latter is again a manifestation of the breathing motion of the impurities at the two-body level. However H^eff fails in general to adequately capture the correct spatial shape of ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$, since e.g. it predicts a repulsion of the impurities (figure 11(c₄)) when in the presence of correlations they attract each other (figure 11(b₄)) and vice versa (compare figures 11(c₃) and (b₃)). This difference is caused by the failure of the effective potential to account for induced interactions emanating within the MB setting.

Finally, turning to strong postquench attractions within ${R}_{{III}}^{{\prime} }$, e.g. for ${g}_{B\uparrow }\,=\,-1$ presented in figures 11(d₁)–(d₅), we observe that the two-body dynamics of the impurities is drastically altered with respect to the weakly attractive case ${g}_{B\uparrow }\,=\,-0.5$ described above. Initially, at t = 0, the two bosons bunch together in the vicinity of the trap center since ${\rho }_{\uparrow \uparrow }^{(2)}(-1\lt {x}_{1}\lt 1,-1\lt {x}_{2}\lt 1;t=0)$ is predominantly populated (figure 11(d₁)). Subsequently the two-body distribution of ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ spatially shrinks exhibiting a highly intense peaked structure around −0.2 < x₁, x₂ < 0.2 as shown in figures 11(d₂), (d₃). For longer evolution times ${\rho }_{\uparrow \uparrow }^{(2)}({x}_{1},{x}_{2};t)$ deforms possessing an elongated shape across its diagonal (see figures 11(d₄), (d₅)) which indicates that the impurities experience a mutual attraction. This latter behavior suggests the appearance of attractive induced interactions between the impurities mediated by the bosonic gas.