# Supplemental Material: Decoherence and momentum relaxation in Fermi-polaron
Rabi dynamics: a kinetic equation approach

Tomasz Wasak<EMAIL_ADDRESS>Institute of Physics, Faculty of Physics,
Astronomy and Informatics, Nicolaus Copernicus University in Toruń,
Grudzia̧dzka 5, 87-100 Toruń, Poland  Max-Planck-Institut für Physik komplexer
Systeme, Nöthnitzer Str. 38, 01187 Dresden, Germany  Matteo Sighinolfi  INO-
CNR BEC Center and Dipartimento di Fisica, Università di Trento, 38123 Trento,
Italy  Johannes Lang  Max-Planck-Institut für Physik komplexer Systeme,
Nöthnitzer Str. 38, 01187 Dresden, Germany  Institut für Theoretische Physik,
Universität zu Köln, Zülpicher Straße 77, 50937 Cologne, Germany  Francesco
Piazza<EMAIL_ADDRESS>Theoretical Physics III, Center for Electronic
Correlations and Magnetism, Institute of Physics, University of Augsburg,
86135 Augsburg, Germany  Max-Planck-Institut für Physik komplexer Systeme,
Nöthnitzer Str. 38, 01187 Dresden, Germany  Alessio Recati
<EMAIL_ADDRESS>INO-CNR BEC Center and Dipartimento di Fisica,
Università di Trento, 38123 Trento, Italy  Trento Institute for Fundamental
Physics and Applications, INFN, 38123, Trento, Italy

In this Supplemental Material (SM), we present derivations of the kinetic
equations (3) of the main text. To this aim, we write the theory in the
Keldysh formalism and then we approximately solve the Dyson equations to
obtain the desired kinetic equations. We also include a derivation of the
collision integral $I_{\alpha}$ and of the polaron decoherence rate
$\Gamma_{\alpha}^{\mathrm{dec}}$. Finally, we also discuss the repulsive
polaron case. For more details on the Keldysh formalism and on its usage, we
refer the reader to [1, 2, 3] and references therein. A more detailed
derivation of the kinetic equations can be found in Ref. [4].

##  I Keldysh formulation of the theory

In our work, we use the Keldysh formalism as presented in [1] to describe our
theory. In this way, we can compute the expectation value of any operator
using the path integral formalism.

The total Hamiltonian $\hat{H}$ of the system is given by Eq. (1) of the main
text. The Keldysh partition function ${\cal Z}$ of $\hat{H}$ is

${\cal Z}=\int\\!{\cal D}\psi_{1}{\cal D}\psi_{2}{\cal
D}\psi_{3}\leavevmode\nobreak\
e^{iS[\psi_{1+},\psi_{1-},\bar{\psi}_{1+},\bar{\psi}_{1-},\psi_{2+},\psi_{2-},\bar{\psi}_{2+},\bar{\psi}_{2-},\psi_{3+},\psi_{3-},\bar{\psi}_{3+},\bar{\psi}_{3-}]},$
(1)

where the integration measures are in the form ${\cal
D}\psi_{i}=\prod_{\sigma=\pm}{\cal D}\bar{\psi}_{i\sigma}{\cal
D}\psi_{i\sigma}$. The fields $\psi_{i\sigma}$ are Grassmann anticommuting
fields for the fermions in level $i$ and the $+(-)$ index labels fields on the
forward (backward) branch of the Keldysh time contour [1]. The total action
can be decomposed, analogously to $\hat{H}$, as
$S=\sum_{i}S_{i}+S_{\Omega}+S_{\rm int}$, i.e., free particle actions $S_{i}$,
the Rabi action $S_{\Omega}$ and the interaction part $S_{\rm int}$. After the
Keldysh rotation of the fermionic fields defined as [1, 5]

$\displaystyle\psi_{i1}$
$\displaystyle=\frac{1}{\sqrt{2}}(\psi_{i+}+\psi_{i-}),\quad\psi_{i2}=\frac{1}{\sqrt{2}}(\psi_{i+}-\psi_{i-})$
(2) $\displaystyle\bar{\psi}_{i1}$
$\displaystyle=\frac{1}{\sqrt{2}}(\bar{\psi}_{i+}-\bar{\psi}_{i-}),\quad\bar{\psi}_{i2}=\frac{1}{\sqrt{2}}(\bar{\psi}_{i+}+\bar{\psi}_{i-}),$
(3)

the free action component is (summation over repeated indices $a,b$ is
understood from now on; each index takes components $1$ or $2$ in the Keldysh
space)

$S_{i}=\iint\\!dxdx^{\prime}\leavevmode\nobreak\
\bar{\psi}_{ia}(x)\hat{G}^{-1}_{0i,ab}(x,x^{\prime})\psi_{ib}(x^{\prime}),$
(4)

with $dx=d{\bf r}dt$ and $\hat{G}^{-1}_{0i}$ is the $2\times 2$ bare inverse
matrix of Green functions (GF); summation over repeated indices is implied.

For our fermionic system, the causal structure of this matrix is given by,

$\hat{G}^{-1}_{0i}=\begin{pmatrix}G^{-1,R}_{0i}&G^{-1,K}_{0i}\\\
0&G^{-1,A}_{0i}\end{pmatrix}.$ (5)

In the last expression, the dependence on $(x,x^{\prime})$ is understood and
R, A, K label respectively the retarded, advanced and Keldysh components of
the GF. With the same convention, the Rabi action is

$S_{\Omega}=\frac{\Omega}{2}\int\\!dx\leavevmode\nobreak\
\big{[}\bar{\psi}_{3a}(x)\psi_{2a}(x)+\bar{\psi}_{2a}(x)\psi_{3a}(x)\big{]}.$
(6)

Note that, differently from the main text, here we reabsorb the detuning
$\Delta$ of level 2 inside the free particle action $S_{2}$. The non-Keldysh-
rotated interaction action is

$S_{\rm int}=-U\sum_{\sigma=\pm}\int\\!dx\leavevmode\nobreak\
\bar{\psi}_{3\sigma}(x)\bar{\psi}_{1\sigma}(x)\psi_{1\sigma}(x)\psi_{3\sigma}(x).$
(7)

$S_{\rm int}$ is quadratic in the bath fields $\psi_{1}$ and in order to
obtain an action that is linear in $\psi_{1}$ we introduce via a Hubbard-
Stratonovich transformation (HST) an auxiliary molecular field
$\Delta_{\sigma}(x)=U\psi_{1\sigma}(x)\psi_{3\sigma}(x)$ with bosonic
statistics. The HST is explicitly written as

$e^{iS_{\rm int}}=\int
D[\bar{\Delta},\Delta]e^{iS_{\Delta}+iS_{\Delta\psi_{1}\psi_{3}}},$ (8)

with the actions $S_{\Delta}$ and $S_{\Delta\psi_{1}\psi_{3}}$ defined on the
Keldysh contour ${\cal C}$

$\displaystyle S_{\Delta}=\int_{\cal C}\\!dx\leavevmode\nobreak\
\bar{\Delta}(x)\frac{1}{U}\Delta(x)$ (9) $\displaystyle
S_{\Delta\psi_{1}\psi_{3}}=-\int_{\cal C}\\!dx\leavevmode\nobreak\
[\bar{\psi}_{3}\bar{\psi}_{1}\Delta+\bar{\Delta}\psi_{1}\psi_{3}].$ (10)

Now $S$ is linear in the bath field $\psi_{1}$ as desired, therefore they can
be integrated out using standard Gaussian integrals. The final effective
action, after the proper Keldysh rotation, reads

$S_{\Delta 3}=-\iint\\!dx\leavevmode\nobreak\ dx^{\prime}\leavevmode\nobreak\
\bar{\zeta}_{a}(x)\hat{G}_{01,ab}(x,x^{\prime})\zeta_{b}(x^{\prime}),$ (11)

where the fields $\zeta$ are defined as

$\zeta_{1}=\frac{\bar{\psi}_{3,1}\Delta_{2}+\bar{\psi}_{3,2}\Delta_{1}}{\sqrt{2}}\quad\zeta_{2}=\frac{\bar{\psi}_{3,1}\Delta_{1}+\bar{\psi}_{3,2}\Delta_{2}}{\sqrt{2}}\quad\bar{\zeta}_{1}=\frac{\bar{\Delta}_{1}\psi_{3,2}+\bar{\Delta}_{2}\psi_{3,1}}{\sqrt{2}}\quad\bar{\zeta}_{2}=\frac{\bar{\Delta}_{1}\psi_{3,1}+\bar{\Delta}_{2}\psi_{3,2}}{\sqrt{2}},$
(12)

with the Keldysh rotation of $\Delta$ defined as in Eq.(2), while
$\bar{\Delta}$ is transformed by

$\bar{\Delta}_{1}=\frac{1}{\sqrt{2}}(\bar{\Delta}_{+}+\bar{\Delta}_{-})\quad\bar{\Delta}_{2}=\frac{1}{\sqrt{2}}(\bar{\Delta}_{+}-\bar{\Delta}_{-}).$
(13)

The different rotations for $\bar{\Psi}$ and $\bar{\Delta}$ reflect different
statistics of the fields. In Eq. (11), we now have an interaction between
$\psi_{3}$ and $\Delta$ mediated by the fermions in state 1. This describes a
scattering between an impurity in level 3 and a molecule, as shown in Fig.1.

Figure 1: Vertex of $S_{\Delta 3}$ in position space as represented in Eq.
(11). Each line carries a Keldysh index, with $\alpha,\,\beta=cl,\,q$ and
$c,\,d=1\,2$.

##  II Derivation of Dyson equations for Keldysh Green functions

Dressed GFs for molecules and impurities can be derived from Dyson equations
in the Keldysh formalism, with the key difference that the $2\times 2$
structure of the GF matrix in this formalism gives 3 different equations to
solve for the R, A, K components. The retarded GF is needed to describe the
excitation spectrum, while the Keldysh GF is related to the distribution
function $F$ via $G^{K}=G^{R}\circ F-F\circ G^{A}$, where the symbol $\circ$
stands for the convolution operator [1].

The Dyson equations for molecules and for atoms in levels 2 and 3 read

$\displaystyle(\hat{G}_{0\Delta}^{-1}-\hat{\Sigma}_{\Delta})\cdot\hat{G}_{\Delta}$
$\displaystyle=\mathbb{1}$ (14)
$\displaystyle(\hat{G}_{0,ij}^{-1}-\hat{\Sigma}_{ij})\cdot\hat{G}_{jk}$
$\displaystyle=\mathbb{1}\delta_{ik},\quad$ (15)

where $i,j,k=2,3$ and the dot product denotes the matrix product combined with
a convolution. In the last equation, we adopted for the bare and dressed GF of
atoms the same spin notation as in Ref. [1]. The bare $G_{0}^{-1\,K}$ are
regularization factors and can therefore be neglected for both atoms and
molecules. The bare retarded and advanced Green’s functions are

$\displaystyle G_{0\Delta}^{-1\,R(A)}(x_{1},x_{2})$
$\displaystyle=\delta(x_{1}-x_{2})U$ (17) $\displaystyle
G_{0,jk}^{-1\,R(A)}(x_{1},x_{2})$
$\displaystyle=\delta(x_{1}-x_{2})\left[\left(i\partial_{t_{2}}+\frac{1}{2m}\nabla^{2}_{\mathbf{r}_{2}}\right)\delta_{jk}-\mathbf{H}\cdot\hat{s}_{jk}+i0^{+}\delta_{jk}\right]\,,$
(18)

where $\hat{s}$ is the vector of Pauli matrices and
$\mathbf{H}=(\Omega/2,0,\Delta/2)$, consistent with the definition of
$H_{\Omega}$. Note that $\Omega$ is in general time dependent. The time-
dependence has to be chosen in order to match the experimental conditions
described in Ref. [6] and the related Supplemental Material. This means that
first the detuning $\Delta$ is tuned in resonance with the desired polaron
branch and then the Rabi frequency is switched on at the initial time of the
dynamics. To solve the Dyson equation we compute the 1-loop self-energies
$\hat{\Sigma}_{jk}$ and $\hat{\Sigma}_{\Delta}$. These are obtained by
contracting the vertex of $S_{\Delta 3}$ with a Green’s function. Following
the summation over internal Keldysh indices one finds

$\displaystyle\Sigma_{\Delta}^{R}(p)$
$\displaystyle=\frac{i}{2V}\sum_{q}\left(G_{01}^{R}(p-q)G_{0,33}^{K}(q)+G_{01}^{K}(p-q)G_{0,33}^{R}(q)\right)$
(19) $\displaystyle\Sigma_{\Delta}^{K}(p)$
$\displaystyle=\frac{i}{2V}\sum_{q}\left\\{G_{01}^{K}(p-q)G_{0,33}^{K}(q)+[G_{01}^{R}(p-q)-G_{01}^{A}(p-q)][G_{0,33}^{R}(q)-G_{0,33}^{A}(q)]\right\\}$
(20) $\displaystyle\Sigma_{33}^{R}(x,p)$
$\displaystyle=-\frac{i}{2V}\sum_{q}\left(G_{01}^{K}(q)G_{\Delta}^{R}(x,p+q)+G_{01}^{A}(q)G_{\Delta}^{K}(x,p+q)\right)$
(21) $\displaystyle\Sigma_{33}^{K}(x,p)$
$\displaystyle=-\frac{i}{2V}\sum_{q}\left(G_{01}^{K}(q)G_{\Delta}^{K}(x,p+q)+G_{01}^{A}(q)G_{\Delta}^{R}(x,p+q)+G_{01}^{R}(q)G_{\Delta}^{A}(x,p+q)\right),$
(22)

with $x=(x_{1}+x_{2})/2$ and $p=(\omega,{\bf p})$ is the relative momentum,
i.e., the conjugate coordinate of $x_{1}-x_{2}$ in the sense of the Fourier
transform. In our approximation scheme, we compute $\Sigma_{\Delta}$ using the
bare GF, Eqs. (17)-(18): since the bare GF are time and space-translation
invariant, the molecular self energy depends only on the relative momentum. We
then use the molecular self energies obtained in this way to derive — as
explained below — dressed molecular propagators $G_{\Delta}$ that we use to
compute $\Sigma_{33}$, Eqs. (21)-(22), that also depends on $x$. In the
impurity limit, i.e., when the density of the impurities is vanishingly small
compared to the density of the bath, this scheme is equivalent to the non-
self-consistent $T$-matrix approximation for $\Sigma_{33}$. A representation
of the self-energies is shown in Fig. 2.

Figure 2:  Schematic representation of the molecular and impurity self
energies, $\Sigma_{\Delta}$ and $\Sigma_{33}$, defined in Eqs. (19)-(22). Wavy
lines denote $G_{01}$, full lines $G_{33}$ and double lines $G_{\Delta}$,
while the labels $R$, $A$, $K$ denote the components of the propagators inside
the loop.

##  III Retarded Green functions and polaron energy

The first quantity we need to know, because it is needed for
$\hat{\Sigma}_{33}$ in our approximation scheme, is $G_{\Delta}^{R}$. From the
Dyson equation for the molecules, Eq. (14), we obtain

$G_{\Delta}^{R}(p)=\frac{1}{U^{-1}-\Sigma_{\Delta}^{R}(p)}$ (23)

The explicit form of $G_{\Delta}^{R}$ is obtained through the relations
$G_{i}^{K}(p)=F_{i}(p)(G_{i}^{R}(p)-G_{i}^{A}(p))$. For the bath at
equilibrium, $F_{1}(p)=1-2n_{1}^{\rm eq}({\bf k})$, with $n_{1}^{\rm eq}$ the
Fermi-Dirac distribution at temperature $T$, while in the impurity limit
$N_{3}\ll N_{1}$, $F_{3}\simeq 1$. The form of $G_{\Delta}$ obtained in this
way is

$G_{\Delta}^{R}(\omega,\mathbf{k})=\left(\frac{m}{4\pi
a_{13}}-\frac{1}{V}\sum_{\mathbf{q}}\frac{1-n_{1}^{\rm
eq}(\mathbf{k}-\mathbf{q})}{\omega-\varepsilon_{3}(\mathbf{q})-\varepsilon_{1}(\mathbf{k}-\mathbf{q})+i0^{+}}+\frac{m}{\mathbf{q}^{2}}\right)^{-1},$
(24)

where the contact interaction strength $U$ is related to the $s$-wave
scattering length $a_{13}$ by the Lippmann-Schwinger equation, so that
$U^{-1}=m/4\pi a_{13}+\sum_{\mathbf{q}}m/\mathbf{q}^{2}$ [7].

$1/k_{F}a$ |  $E_{0}$ $[\varepsilon_{F}]$ | $m_{*}/m$ | $Z_{a}(0)$ | $\partial_{\omega}\text{Im}\Sigma_{33}^{R}(0)$
---|---|---|---|---
0 | -0.625 | 1.16 | 0.775 | -69.0 $\times 10^{-4}$
0.25 | -0.858 | 1.29 | 0.673 | -33.5 $\times 10^{-4}$
Table 1: Attractive polaron parameters at $T=0.135T_{F}$.

$1/k_{F}a$ |  $E_{0}$ $[\varepsilon_{F}]$ | $m_{*}/m$ | $Z_{r}(0)$ | $\partial_{\omega}\text{Im}\Sigma_{33}^{R}(0)$
---|---|---|---|---
1.27 | 0.423 | 1.13 | 0.823 | -2.7 $\times 10^{-4}$
2.63 | 0.188 | 1.07 | 0.963 | -0.3 $\times 10^{-4}$
Table 2: Repulsive polaron parameters at $T=0.135T_{F}$

$1/k_{F}a$ | $1-\sum_{\alpha}Z_{\alpha}$
---|---
0.0 | 0.225
0.25 | 0.140
1.27 | 0.037
2.63 | 0.018
Table 3: Sum of quasiparticle residues at $T=0.135$.

The same approach applied to $G_{33}^{R}$ gives

$G_{33}^{R}(x,p)=\left(\omega-\varepsilon_{3}(\mathbf{k})-\Sigma_{33}^{R}(x,p)-\frac{\Omega^{2}/4}{\omega-\varepsilon_{2}(\mathbf{k})+\Delta+i0}\right)^{-1},$
(25)

where the last term describes the Rabi coupling to state 2. The polaron energy
is the pole of the propagator, and the Rabi coupling splits each branch of the
polaron spectrum in two sub-branches. Indeed, the $\alpha$-polaron dispersion
relation is

$\varepsilon_{\alpha\,1,2}(\mathbf{k})=\varepsilon_{3}(\mathbf{k})+\frac{\text{Re}\Sigma_{33}^{R}(\varepsilon_{\alpha\,1,2}(\mathbf{k}),\mathbf{k})-\Delta}{2}\pm\frac{1}{2}\sqrt{(\text{Re}\Sigma_{33}^{R}(\varepsilon_{\alpha\,1,2}(\mathbf{k}),\mathbf{k})+\Delta)^{2}+\Omega^{2}}.$
(26)

However, we will not consider the dressing effect of the Rabi coupling on the
branches in our treatment, leaving it for future development. Therefore, the
polaron energy dispersion relation has the usual form

$\varepsilon_{\alpha}(\mathbf{k})=\varepsilon_{3}(\mathbf{k})+\text{Re}\Sigma_{33}^{R}(\varepsilon_{\alpha}(\mathbf{k}),\mathbf{k}).$
(27)

Finally, a comment on $\text{Re}\Sigma_{33}^{R}$ is also useful: we derive it
in the impurity limit, where the non-self consistent $T$-matrix approximation
is recovered as known for equilibrium [8]. Within this approximation polaron
dispersion relation, quasiparticle residue $Z_{\alpha}(\mathbf{k}=0)$ and
effective mass $m_{*}$ are derived. Results obtained in the non-self
consistent $T$-matrix approximation are known to be in good agreement with
experiments and quantum Monte Carlo calculations, see e.g. Ref. [9] and
references therein. The results used in the numerical simulation are listed in
Tabs. 3 and 3. As we will show below, we have to go beyond the impurity limit
to describe the dynamics of populations.

##  IV Kinetic equations, collision integrals and decoherence rates

The last component of the Dyson equation is the Keldysh one, and it is
necessary to derive kinetic equations. We start our derivation from the
molecular Dyson equation, where the Keldysh component is

$G_{\Delta}^{-1\,R}\circ G_{\Delta}^{K}\circ
G_{\Delta}^{-1\,A}=\Sigma_{\Delta}^{K}-G_{0,\Delta}^{-1\,K}.$ (28)

We then use the relation $G_{\Delta}^{K}=G_{\Delta}^{R}\circ
F_{\Delta}-F_{\Delta}\circ G_{\Delta}^{A}$ for the distribution function
$F_{\Delta}$ and we insert it into Eq. (28) to obtain an equation for the
distribution function of molecules that reads

$F_{\Delta}\circ G_{\Delta}^{-1\,A}-G_{\Delta}^{-1\,R}\circ
F_{\Delta}=\Sigma_{\Delta}^{K}-G_{0,\Delta}^{-1\,K}.$ (29)

We then apply the Wigner transformation [1, 10] to simplify the equation for
$F_{\Delta}$. We consider only the zeroth order of the expansion in gradients,
treating the dynamics of the molecules as subordinated to the dynamics of the
polarons. In this approximation, Eq. (29) reads

$2iF_{\Delta}\text{Im}\Sigma_{\Delta}^{R}=\Sigma_{\Delta}^{K},$ (30)

and the following equations hold:

$\displaystyle
F_{\Delta}(G_{\Delta}^{R}-G_{\Delta}^{A})=|G_{\Delta}^{R}|^{2}\Sigma_{\Delta}^{K}$
(31) $\displaystyle F_{\Delta}{\cal
A}_{\Delta}=|G_{\Delta}^{R}|^{2}i\Sigma_{\Delta}^{K}$ (32) $\displaystyle{\cal
A}_{\Delta}=-2|G_{\Delta}^{R}|^{2}\text{Im}\Sigma_{\Delta}^{R},$ (33)

where ${\cal A}_{\Delta}=-2\text{Im}G_{\Delta}^{R}$.

Now we consider the impurities. From Eq. (15) we have 4 kinetic equations for
$F_{22}$, $F_{33}$, $F_{23}$ and $F_{32}$ but only 3 of them are relevant,
because $F_{32}=F_{23}^{*}$. It is convenient to define
$\tilde{G}_{0,jk}^{-1\,R(A)}$ as

$\tilde{G}_{0,jk}^{-1\,R(A)}(x_{1},x_{2})=\delta(x_{1}-x_{2})\delta_{jk}\left(i\partial_{t_{2}}+\frac{1}{2m}\nabla^{2}_{\mathbf{r}_{2}}\right),$
(34)

so that
$G_{0,jk}^{-1\,R(A)}=\tilde{G}_{0,jk}^{-1\,R(A)}-\mathbf{H}\cdot\hat{s}_{jk}$

From the matrix structure of the Dyson equation for the impurities we have
that

$\displaystyle
F_{22}\circ\tilde{G}_{0,22}^{-1\,A}-\tilde{G}_{0,22}^{-1\,R}\circ
F_{22}+\frac{\Omega}{2}\circ F_{23}^{*}-F_{23}\circ\frac{\Omega}{2}=0,$ (35)
$\displaystyle
F_{33}\circ\tilde{G}_{0,33}^{-1\,A}-\tilde{G}_{0,33}^{-1\,R}\circ
F_{33}+\frac{\Omega}{2}\circ
F_{23}-F_{23}^{*}\circ\frac{\Omega}{2}+\Sigma_{33}^{R}\circ
F_{33}-F_{33}\circ\Sigma_{33}^{A}=\Sigma_{33}^{K},$ (36) $\displaystyle
F_{23}\circ\tilde{G}_{0,33}^{-1\,A}-\tilde{G}_{0,22}^{-1\,R}\circ
F_{23}+\frac{\Omega}{2}\circ F_{33}-F_{22}\circ\frac{\Omega}{2}+\Delta
F_{23}=F_{23}\circ\Sigma_{33}^{A}.$ (37)

We now perform a Wigner transformation at linear order in gradients on the
previous equations. Note that the free fermionic energies $\varepsilon_{2}$
and $\varepsilon_{3}$ are the same and this leads to (the dependence
$F_{ij}(x,k)$ is not shown)

$\displaystyle\left\\{\omega-\varepsilon_{2}(\mathbf{k}),\,F_{22}\right\\}-i\frac{\Omega}{2}(F_{23}-F_{23}^{*})-\frac{1}{2}\left\\{\frac{\Omega}{2},F_{23}+F_{23}^{*}\right\\}=0$
(38) $\displaystyle\
\left\\{\omega-\varepsilon_{3}(\mathbf{k})-\text{Re}\Sigma_{33}^{R}(x,k),\,F_{33}\right\\}+i\frac{\Omega}{2}(F_{23}-F_{23}^{*})-\frac{1}{2}\left\\{\frac{\Omega}{2},F_{23}+F_{23}^{*}\right\\}=i\Sigma_{33}^{K}(x,k)+2\text{Im}\Sigma_{33}^{R}(x,k)F_{33}$
(39)
$\displaystyle\left\\{\omega-\varepsilon_{3}(\mathbf{k})-\frac{1}{2}\text{Re}\Sigma_{33}^{R}(x,k),\,F_{23}\right\\}-i(\text{Re}\Sigma_{33}^{R}(x,k)-\Delta)F_{23}-i\frac{\Omega}{2}(F_{22}-F_{33})-\frac{1}{2}\left\\{\frac{\Omega}{2},F_{22}+F_{33}\right\\}+$
$\displaystyle\hskip
156.49014pt+\frac{i}{2}\\{\text{Im}\Sigma_{33}^{R}(x,k),F_{23}\\}=\text{Im}\Sigma_{33}^{R}(x,k)F_{23},$
(40)

where the Poisson brackets are defined as
$\\{A,B\\}=\partial_{x}A\partial_{k}B-\partial_{k}A\partial_{x}B$ and
$\partial_{x}A\partial_{k}B=\nabla_{\mathbf{r}}A\nabla_{\mathbf{k}}B-\partial_{t}A\partial_{\omega}B$.

Note that the right-hand side of the Eq. (39) defines the collision integral
$\tilde{I}_{\rm coll}=i\Sigma_{33}^{K}+2\text{Im}\Sigma_{33}^{R}F_{33}$, while
in the right-hand side of Eq. (40) we only have the imaginary part of the
retarded self-energy, so this term is a decay rate. We also point out that we
are not in the impurity limit: indeed a small but finite concentration of
impurities is necessary to compare with experimental results. Note that in
Eqs. (38)-(40) the explicit time dependence of the Rabi frequency is
$\Omega(t)=\theta(t)\Omega$, but it will be not considered form now on 111Even
if this time dependence is included, it does not change the final kinetic
equations because it will be removed by the procedure adopted..

The equations above describe the distribution functions $F_{ij}$ for any
$\omega$ and we do not have the equations for populations yet. In systems with
a Rabi coupling, the connection between $F_{ij}$ and the populations of the
levels is not trivial to find, as opposed to the common situation in
literature [1, 12, 10].

A way to find this connection is to work with the Kadanoff-Baym equation [13]
for the Keldysh component of the GF. This equation is, in matrix form

$\displaystyle D_{0}G^{\rm K}-G^{\rm K}D_{0}+\frac{i}{2}(\Gamma G^{\rm
K}+G^{\rm K}\Gamma)+\frac{i}{2}\left[\left\\{D_{0},G^{\rm
K}\right\\}-\left\\{G^{\rm
K},D_{0}\right\\}\right]-\frac{1}{4}\left[\left\\{\Gamma,G^{\rm
K}\right\\}+\left\\{G^{\rm K},\Gamma\right\\}\right]=$ $\displaystyle\hskip
85.35826pt=\Sigma^{\rm K}G^{\rm A}-G^{\rm R}\Sigma^{\rm
K}+\frac{i}{2}\left[\left\\{\Sigma^{\rm K},G^{\rm A}-\left\\{G^{\rm
R},\Sigma^{\rm K}\right\\}\right\\}\right],$ (41)

where

$\displaystyle D_{0}$
$\displaystyle=\begin{pmatrix}&\omega-\varepsilon_{2,{\bf p}}&-\Omega/2\\\
&-\omega/2&\omega-\varepsilon_{3,{\bf
p}}-\text{Re}\Sigma_{33}^{R}+\Delta\end{pmatrix}$ (42) $\displaystyle\Gamma$
$\displaystyle=\begin{pmatrix}&0&0\\\
&0&-2\text{Im}\Sigma_{33}^{R}\end{pmatrix}=\begin{pmatrix}&0&0\\\
&0&\Gamma_{33}\end{pmatrix}.$ (43)

When written in components, Kadanoff-Baym equation leads to equation that are
similar to the ones for the distribution functions

$\displaystyle\partial_{t}G^{\rm K}_{22}+i\frac{\Omega}{2}(G^{\rm
K}_{23}-G^{\rm K}_{32})=0$ (44)
$\displaystyle(1-\partial_{\omega}\text{Re}\Sigma_{33}^{R})\partial_{t}G^{\rm
K}_{33}-i\frac{\Omega}{2}(G^{\rm K}_{23}-G^{\rm K}_{23})=i\Sigma_{33}^{\rm
K}(G_{33}^{\rm R}-G_{33}^{\rm A})$ (45)
$\displaystyle\left(1-\frac{1}{2}\partial_{\omega}\text{Re}\Sigma_{33}^{R}\right)\partial_{t}G^{\rm
K}_{23}-i(\text{Re}\Sigma_{33}^{R}-\Delta)G^{\rm
K}_{23}+i\frac{\Omega}{2}(G^{\rm K}_{22}-G^{\rm K}_{33})+$
$\displaystyle\hskip
71.13188pt+\frac{i}{2}\partial_{\omega}\text{Im}\Sigma_{33}^{R}\partial_{t}G^{\rm
K}_{23}=-\frac{\Gamma_{33}}{2}G^{\rm K}_{23}+iG_{23}^{\rm R}\Sigma_{33}^{K}.$
(46)

At linear order in gradients of the WT, the matrix $\hat{G}^{\rm K}$ can be
parametrized as

$\hat{G}^{\rm K}=-\frac{i}{2}(\hat{\cal A}\hat{F}-\hat{F}\hat{\cal A}),$ (47)

where $\hat{\cal A}$ is the matrix of spectral functions. In the limit of
small $\Omega$, i.e. when the Rabi coupling is a probe, the energies of levels
$\left|2\right\rangle$ and $\left|3\right\rangle$ are not modified and the
components of $\hat{G}^{\rm K}$ are

$\displaystyle G_{22}^{\rm K}$ $\displaystyle=-i{\cal A}_{22}F_{22}$ (48)
$\displaystyle G_{23}^{\rm K}$ $\displaystyle=-\frac{i}{2}({\cal A}_{22}+{\cal
A}_{33})F_{23}$ (49) $\displaystyle G_{32}^{\rm K}$
$\displaystyle=-\frac{i}{2}({\cal A}_{22}+{\cal A}_{33})F_{32}$ (50)
$\displaystyle G_{33}^{\rm K}$ $\displaystyle=-i{\cal A}_{33}F_{33},$ (51)

and in the same limit the spectral functions are

$\displaystyle{\cal A}_{22}$
$\displaystyle=2\pi\delta(\omega-\varepsilon_{2}({\bf p}))$ (52)
$\displaystyle{\cal A}_{33}$
$\displaystyle=\sum_{\alpha=a,r}Z_{\alpha}2\pi\delta(\omega-\varepsilon_{\alpha}({\bf
p})),$ (53)

where $\varepsilon_{2}$ is the bare energy of level $\left|2\right\rangle$ and
$\varepsilon_{\alpha}$ is given by the polaron energy not modified by the Rabi
coupling, Eq. (27). With the above results for spectral functions and
components of $\hat{G}^{\rm K}$, it is possible to perform and integration
over $\omega$ and obtain the so-called on-shell equation for the distribution
functions. The advantage of using the Kadanoff-Baym equation instead of the
Dyson’s equation is that for the former the on-shell projection is naturally
defined, because all the distribution functions in the equations are
multiplied by a spectral function. In the on-shell equations, also the terms
with $\partial_{\omega}\text{Im}\Sigma_{33}^{R}$ is neglected because they are
small, see Tables 3 and 3. The final result is

$\displaystyle\partial_{t}F_{22}^{(2)}-i\frac{\Omega}{2}\left[\frac{1}{2}(F_{23}^{(2)}-F_{32}^{(2)})+\frac{1}{2}\sum_{\alpha}Z_{\alpha}(F_{23}^{(\alpha)}-F_{32}^{(\alpha)})\right]$
(54)
$\displaystyle\sum_{\alpha}\partial_{t}F_{33}^{(\alpha)}+i\frac{\Omega}{2}\left[\frac{1}{2}(F_{23}^{(2)}-F_{32}^{(2)})+\frac{1}{2}\sum_{\alpha}Z_{\alpha}(F_{23}^{(\alpha)}-F_{32}^{(\alpha)})\right]=\sum_{\alpha}(i\Sigma_{33}^{{\rm
K}\,(\alpha)}+2\text{Im}\Sigma_{33}^{{\rm R}\,(\alpha)}F_{33}^{(\alpha)})$
(55)
$\displaystyle\left(1-\frac{1}{2}\partial_{\omega}\text{Re}\Sigma_{33}^{{\rm
R}\,(2)}\right)\frac{1}{2}\partial_{t}F_{23}^{(2)}+\sum_{\alpha}\tilde{Z}_{\alpha}^{-1}Z_{\alpha}\frac{1}{2}\partial_{t}F_{23}^{(\alpha)}-i\left[(\text{Re}\Sigma_{33}^{{\rm
R}\,(2)}-\Delta)\frac{1}{2}F_{23}^{(2)}+\sum_{\alpha}\times\right.$
$\displaystyle\left.\times(\text{Re}\Sigma_{33}^{{\rm
R}\,(\alpha)}-\Delta)\frac{1}{2}Z_{\alpha}F_{23}^{(\alpha)}\right]+i\frac{\Omega}{2}(\sum_{\alpha}Z_{\alpha}F_{33}^{(\alpha)}-F_{22}^{(2)})=\frac{1}{2}\left(\text{Im}\Sigma_{33}^{(R)\,(2)}F_{23}^{(2)}+\sum_{\alpha}Z_{\alpha}\text{Im}\Sigma_{33}^{(R)\,(\alpha)}F_{23}^{(\alpha)}\right),$

where in the last equation the term $iG_{23}^{\rm R}\Sigma_{33}^{\rm K}$
vanishes on shell. The notation

$\displaystyle F_{ij}^{(n)}$ $\displaystyle\equiv
F_{ij}(\omega=\varepsilon_{n})$ (57) $\displaystyle Z_{\alpha}$
$\displaystyle=(1-\partial_{\omega}\text{Re}\Sigma_{33}^{R}|_{\omega=\varepsilon_{\alpha}})^{-1}$
(58) $\displaystyle\tilde{Z}_{\alpha}$
$\displaystyle=\left(1-\frac{1}{2}\partial_{\omega}\text{Re}\Sigma_{33}^{R}|_{\omega=\varepsilon_{\alpha}}\right)^{-1},$
(59)

has been used. Now, from the above equations for the on-shell distribution
functions it is possible to obtain the ones for the populations. Indeed, in
general one has that

$\displaystyle i\int_{\omega}G_{22}^{\rm K}$ $\displaystyle\equiv
1-2n_{22}^{(2)}=F_{22}^{(2)}$ (60) $\displaystyle i\int_{\omega}G_{33}^{\rm
K}$
$\displaystyle\equiv\sum_{\alpha}1-2n_{33}^{(\alpha)}=\sum_{\alpha}Z_{\alpha}F_{33}^{(\alpha)}$
(61) $\displaystyle i\int_{\omega}G_{23}^{\rm K}$ $\displaystyle\equiv
f_{23}^{(2)}+\sum_{\alpha}f_{23}^{(\alpha)}=f_{23}=\frac{1}{2}(F_{23}^{(2)}+\sum_{\alpha}Z_{\alpha}F_{23}^{(\alpha)})$
(62)

With the above definition, the equations for the populations in
$\left|2\right\rangle$ and $\left|3\right\rangle$ are easily obtained

$\displaystyle\partial_{t}n_{22}^{(2)}+i\frac{\Omega}{2}(f_{23}-f_{23}^{*})=0$
(63)
$\displaystyle\sum_{\alpha}\frac{1}{Z_{\alpha}}\partial_{t}n_{33}^{(\alpha)}-i\frac{\Omega}{2}(f_{23}-f_{23}^{*})=\sum_{\alpha}I_{\rm
coll}^{\alpha}.$ (64)

From the above equations, the conservation of particles is obtained, i.e.,

$\partial_{t}(N_{2}+N_{3})=\int_{{\bf
p}}\partial_{t}n_{22}^{(2)}+\sum_{\alpha}\frac{1}{Z_{\alpha}}\partial_{t}n_{33}^{(\alpha)}=0.$
(65)

The conservation of particles inside the kinetic equations is a consequence of
the assumption made on ${\cal A}_{33}$, Eq. (53), where the non-coherent part
of the spectral function is totally neglected. Finally, in order to close the
kinetic equations, only a single polaron branch has to be considered. In this
case, both $n_{ii}$ and $f_{23}$ are projected on the same energy up to an
effective mass correction and the kinetic equations become

$\displaystyle\partial_{t}n_{2}+i\frac{\Omega}{2}(f_{23}-f_{23}^{*})=0$ (66)
$\displaystyle\partial_{t}n_{\alpha}-iZ_{\alpha}\frac{\Omega}{2}(f_{23}-f_{23}^{*})=I_{\rm
coll}^{\alpha}$ (67)
$\displaystyle\partial_{t}f_{23}+i\tilde{Z}_{\alpha}\delta_{\alpha}f_{23}+i\frac{\Omega}{2}(n_{\alpha}-n_{2})=-\frac{\Gamma_{\alpha}^{\rm
dec}}{2}f_{23},$ (68)

with $n_{\alpha}=n_{33}^{(\alpha)}$, $\delta_{\alpha}({\bf
p})=\varepsilon_{\alpha}({\bf p})-\varepsilon_{2}({\bf p})-\Delta$ and
$\Gamma_{\alpha}^{\mathrm{dec}}$ is the decoherence rate. The absence of the
tilde accent over collisional integral and decoherence rate indicates that
they are written in terms of populations and not of distribution functions.

To better understand the kinetic equations, it is necessary to derive the
explicit form of the collisional integral and of the decoherence rate. The
strategy is the following: one starts from the definitions, uses the kinetic
equation for the dynamics of the molecules, projects on the polaron energy and
then write everything in terms of the populations.

Starting from the definition of the collision integral $\tilde{I}_{\rm coll}$
given below Eq. (40), where Eqs. (31-33) and the definitions of self-energies,
Eqs. (19-22) are used, we obtain the expression

$\displaystyle\tilde{I}_{\rm coll}$
$\displaystyle=\left(\frac{1}{2V}\right)^{2}\sum_{q,q^{\prime}}|G_{\Delta}^{R}(q)|^{2}{\cal
A}_{\phi}(q-p){\cal A}_{\phi}(q-q^{\prime}){\cal A}_{33}(q^{\prime})\times$
$\displaystyle\times[(F_{33}(p)+F_{\phi}(q-p))(1+F_{\phi}(q-q^{\prime})F_{33}(q^{\prime}))-(1+F_{33}(p)F_{\phi}(q-p))(F_{33}(q^{\prime})+F_{\phi}(q-q^{\prime}))],$

After the integration over energies, that can be performed using the spectral
functions in Eq. (LABEL:I_coll_1), we obtain

$\displaystyle\tilde{I}_{\rm coll}$
$\displaystyle=\frac{2\pi}{(2V)^{2}}\sum_{\mathbf{q},\mathbf{q}^{\prime}}\sum_{\beta={\rm
a,r}}|G_{\Delta}^{R}(\omega+\varepsilon_{\phi}(\mathbf{q}-{\bf
p}),\mathbf{q})|^{2}Z_{\beta}(\mathbf{q}^{\prime})\delta(\omega+\varepsilon_{\phi}(\mathbf{q}-{\bf
p})-\varepsilon_{\beta}(\mathbf{q}^{\prime})-\varepsilon_{\phi}(\mathbf{q}-\mathbf{q}^{\prime}))\times$
$\displaystyle\times[(F_{33}(p)+F_{\phi}(\mathbf{q}-{\bf
p}))(1+F_{\phi}(\mathbf{q}-\mathbf{q}^{\prime})F_{33}^{\beta}(\mathbf{q}^{\prime}))-(1+F_{33}(p)F_{\phi}(\mathbf{q}-{\bf
p}))(F_{33}^{\beta}(\mathbf{q}^{\prime})+F_{\phi}(\mathbf{q}-\mathbf{q}^{\prime}))].$
(70)

Note that the full $G_{\Delta}^{R}$ is equivalent to the $T$-matrix, so the
notation $G_{\Delta}^{R}=T_{\rm sc}$ is used from now on. After the on-shell
projection, the term inside $\sum_{\alpha}$ in the right-hand side of Eq. (55)
becomes

$\displaystyle\tilde{I}_{\rm coll}^{\alpha}$
$\displaystyle=\frac{2\pi}{(2V)^{2}}\sum_{\mathbf{q},\mathbf{q}^{\prime}}\sum_{\beta={\rm
a,r}}|T_{\rm sc}(\varepsilon_{\alpha}({\bf
p})+\varepsilon_{\phi}(\mathbf{q}-{\bf
p}),\mathbf{q})|^{2}Z_{\beta}(\mathbf{q}^{\prime})\delta(\varepsilon_{\alpha}({\bf
p})+\varepsilon_{\phi}(\mathbf{q}-{\bf
p})-\varepsilon_{\beta}(\mathbf{q}^{\prime})-\varepsilon_{\phi}(\mathbf{q}-\mathbf{q}^{\prime}))\times$
$\displaystyle\times[(F_{33}^{\alpha}({\bf p})+F_{\phi}(\mathbf{q}-{\bf
p}))(1+F_{\phi}(\mathbf{q}-\mathbf{q}^{\prime})F_{33}^{\beta}(\mathbf{q}^{\prime}))-(1+F_{33}^{\alpha}({\bf
p})F_{\phi}(\mathbf{q}-{\bf
p}))(F_{33}^{\beta}(\mathbf{q}^{\prime})+F_{\phi}(\mathbf{q}-\mathbf{q}^{\prime}))].$
(71)

The physical interpretation of the above is clear: the $T$-matrix plays the
role of a cross section, the Dirac delta enforces the energy conservation in
the collision and the terms inside the curly brackets are the in and out
processes of the collision and $\sum_{\beta}$ is responsible for both the
relaxation and the conversion processes between the polaron branches. Now it
is possible to rewrite everything in terms of populations, obtaining the final
form of collision integral inside Eq. (67)

$I^{\alpha}_{\rm coll}=\sum_{\beta}I_{\alpha\beta},$ (72)

where $\sum_{\beta}I_{\alpha\beta}$ describes a scattering process between
polarons in $\alpha$ and $\beta$. $I_{\alpha\beta}$ can be written as

$I_{\alpha\beta}({\bf p})=\frac{1}{V}\sum_{\mathbf{q}^{\prime}}\left\\{W_{{\bf
p}\mathbf{q}^{\prime}}^{\alpha\beta}[1-n_{\alpha}({\bf
p})]n_{\beta}(\mathbf{q}^{\prime})-W_{\mathbf{q}^{\prime}{\bf
p}}^{\beta\alpha}[1-n_{\beta}({\bf
p}^{\prime})]n_{\alpha}(\mathbf{q})\right\\}.$ (73)

In this form, the collisional integral is expressed in terms of the transition
rate $W_{\mathbf{q}^{\prime}{\bf p}}^{\beta\alpha}$ (from the $\alpha$ polaron
with momentum ${\bf p}$ to the $\beta$ polaron with momentum
$\mathbf{q}^{\prime}$) and of its complementary process. The transition rates
follow the Fermi golden rule and are given by

$\displaystyle W_{{\bf p}\mathbf{q}^{\prime}}^{\alpha\beta}=$
$\displaystyle\frac{2\pi}{V}\\!\\!\sum_{\mathbf{q}}|T_{\mathrm{sc}}(\varepsilon_{\beta}(\mathbf{q}^{\prime})\\!+\\!\varepsilon_{\phi}(\mathbf{q}-\mathbf{q}^{\prime}),\mathbf{q})|^{2}Z_{\alpha}(\mathbf{q})Z_{\beta}(\mathbf{q}^{\prime})\times\delta(\varepsilon_{\alpha}({\bf
p})+\varepsilon_{\phi}(\mathbf{q}-{\bf
p})-\varepsilon_{\phi}(\mathbf{q}-\mathbf{q}^{\prime})-\varepsilon_{\beta}(\mathbf{q}^{\prime}))$
$\displaystyle\times n_{\phi}^{\rm
eq}(\mathbf{q}-\mathbf{q}^{\prime})[1-n_{\phi}^{\rm eq}(\mathbf{q}-{\bf
p})],.$ (74)

where $n_{\phi}^{\rm eq}$ is the thermal distribution of the majority atoms
with chemical potential $\mu$ and energy
$\varepsilon_{\phi}(\mathbf{k})=\mathbf{k}^{2}/2m$,
$n_{\phi}^{\mathrm{eq}}(\mathbf{k})=1/[e^{(\varepsilon_{\phi}(\mathbf{k})-\mu)/T}+1]$.
A very similar procedure to what has been done for the collisional integral
can be applied also to the decoherence rate,
$\tilde{\Gamma}_{33}^{\text{dec}}=-2\text{Im}\Sigma_{33}^{R}$, obtaining

$\displaystyle\tilde{\Gamma}_{33}^{\rm dec}$
$\displaystyle=\frac{2\pi}{V^{2}}\sum_{\mathbf{q},\mathbf{q}^{\prime}}\sum_{\beta={\rm
a,r}}|T_{\rm
sc}(\varepsilon_{\beta}(\mathbf{q}^{\prime})+\varepsilon_{\phi}(\mathbf{q}-\mathbf{q}^{\prime}),\mathbf{q})|^{2}Z_{\beta}(\mathbf{q}^{\prime})\delta(\omega+\varepsilon_{\phi}(\mathbf{q}-{\bf
p})-\varepsilon_{\phi}(\mathbf{q}-\mathbf{q}^{\prime})-\varepsilon_{\beta}(\mathbf{q}^{\prime}))$
$\displaystyle\times[(1-n_{\phi}^{\rm
eq}(\mathbf{q}-\mathbf{q}^{\prime}))(1-n_{\beta}(\mathbf{q}^{\prime}))n_{\phi}^{\rm
eq}(\mathbf{q}-{\bf p})+(1-n_{\phi}^{\rm eq}(\mathbf{q}-{\bf p}))n_{\phi}^{\rm
eq}(\mathbf{q}-\mathbf{q}^{\prime})n_{\beta}(\mathbf{q}^{\prime})].$ (75)

After the projection on polaron energies and the substitution of distribution
functions with population it is straightforward to obtain the final form of
the decoherence rate

$\Gamma_{\alpha}^{\rm dec}=\frac{\tilde{Z}_{\alpha}({\bf p})}{Z_{\alpha}({\bf
p})}\frac{1}{V}\sum_{\mathbf{q}^{\prime}}\sum_{\beta={\rm a,r}}W_{{\bf
p}\mathbf{q}^{\prime}}^{\alpha\beta}n_{\beta}(\mathbf{q}^{\prime})+W_{\mathbf{q}^{\prime}{\bf
p}}^{\beta\alpha}(1-n_{\beta}(\mathbf{q}^{\prime})),$ (76)

where the transition rates defined in Eq. (74) have been used. In this
expression, the first term is related to the in processes, while the second is
related to the out ones. Note that while both of them are originating from
collisions of the impurity atom with the bath, collisional integral and
decoherence rate are different in nature. The first one is indeed related to
the imbalance of populations at different momenta and eventually at different
branch, the second one affects coherences. In the long time, the collisional
integral drives the population of minority atoms in thermal equilibrium with
the bath.

In the impurity limit, where $n_{\alpha}(\mathbf{k})\approx 0$, the
decoherence rate takes the form

$\Gamma^{\mathrm{dec}}_{\alpha}(\mathbf{k})\\!=\\!\frac{\tilde{Z}_{\alpha}(\mathbf{k})}{Z_{\alpha}(\mathbf{k})}\\!\frac{1}{V}\\!\\!\sum_{\mathbf{k}^{\prime},\beta}W_{\mathbf{k}^{\prime}\mathbf{k}}^{\beta\alpha},$
(77)

while $I_{\alpha\beta}$ can now be written as

$I_{\alpha\beta}(\mathbf{k})=\frac{1}{V}\sum_{\mathbf{k}^{\prime}}\bigg{[}W_{\mathbf{k}\mathbf{k}^{\prime}}^{\alpha\beta}n_{\beta}(\mathbf{k}^{\prime})-W_{\mathbf{k}^{\prime}\mathbf{k}}^{\beta\alpha}n_{\alpha}(\mathbf{k})\bigg{]}.$
(78)

As briefly mentioned above, it is important to notice that, due to the
omission of the back-flow term, the kinetic equations Eqs. (66-68) conserve
the sum $N_{3}(t)+\tilde{N}_{2}(t)$, with
$\tilde{N}_{2}(t)=\sum_{\mathbf{k}}Z_{\alpha}n_{2}(\mathbf{k},t)$ and
$N_{3}(t)=\sum_{\beta,\mathbf{k}}n_{\beta}(t,\mathbf{k})$, where the polaron
$\alpha$ is on resonance with the non-interacting state. Since at the initial
time, a fraction $1-Z_{\alpha}$ is missing from $\tilde{N}_{2}$ with respect
to the total $N_{2}$, in order to implement the conservation law we assume for
consistency that this missing part is not evolving in time and directly add it
to $\tilde{N}_{2}(t)$ to obtain $N_{2}(t)$. Therefore, the experimental
observable $N_{2}/(N_{2}+N_{3})$ can be rewritten as
$1-{\sum_{\alpha,\mathbf{k}}}n_{\alpha}(t,\mathbf{k})/N_{2}(t=0)$. We note
that the same form can be obtained by using the total number conservation to
rewrite $N_{2}(t)/(N_{2}(t)+N_{3}(t))=1-N_{3}(t)/(N_{2}(t)+N_{3}(t))$. By
adding to kinetic equations back-flow terms and gradient corrections to the
collision integral encoding memory effects, proper particle-number
conservation can be achieved [14], but at the cost of largely increasing the
complexity of the equations as well as the computational demand of their
solution. This is beyond the scope of the present work and will be the subject
of future investigations.

##  V Repulsive polaron

In this section, we provide a comparison between the observed Rabi oscillation
for the repulsive polaron, see Ref. [6], and the theoretical predictions
derived using the kinetic equations, Eqs. (66-68). To do so, we have to
compute $\Gamma_{\mathrm{r}}^{\mathrm{dec}}$ and $I_{\mathrm{r}}$ and insert
them into kinetic equations. For the decoherence rate, we cannot simply
replace a$\rightarrow$r in Eq. (7) of the main text, because
$\Gamma_{\mathrm{r}}^{\mathrm{dec}}$ and $I_{\mathrm{r}}$ also contain terms
$\propto W^{\mathrm{a}\mathrm{r}}$, describing collisions that result in de-
excitation to the attractive branch. However, the observed decay rate [6]
suggests that such a process is absent at weak interactions, where it is
instead replaced by a much slower three-body recombination to the molecular
state. Since we will present results for $k_{F}a>0$ and far from unitarity, we
will not consider the de-excitation process.

Our predictions for the Rabi oscillations of repulsive polarons, i.e., the
case when the laser drive is resonant with the transition between state
$|2\rangle$ and the repulsive polaron branch, i.e.,
$\delta_{\mathrm{r}}(0)=0$, are reported in Fig. 3, where we plot the same
quantity of the main text, i.e.,

$\frac{N_{2}(t)}{N_{2}(t)+N_{3}(t)}=1-\frac{\sum_{\alpha,\mathbf{k}}n_{\alpha}(t,\mathbf{k})}{N_{2}(t=0)}.$
(79)

We consider the experimentally relevant cases $1/(k_{F}a)=2.63$ and 1.27,
corresponding to the polaron parameters
$(E_{\mathrm{r}}/\varepsilon_{F},\,m_{*}/m,\,Z_{\mathrm{r}}(0))$ =
$(0.442,\,1.13,\,0.823)$ and $(0.188,\,1.07,\,0.963)$, respectively, see Tab.
3.

Figure 3: Rabi oscillations of the repulsive polarons for $1/(k_{F}a)=$ 2.63
and 1.27 (upper and lower panels, respectively). The blue solid line is
calculated from Eq. (LABEL:eq:keq), whereas the red dots are from the
experiment [6]. The parameters are $\Omega=0.7\varepsilon_{F}$,
$T=0.135\varepsilon_{F}$, $\delta_{\mathrm{r}}(\mathbf{k}=0)=0$. The shaded
blue region (mostly unnoticeable) corresponds to a $20\%$ uncertainty in the
temperature of the majority component.

The dots correspond to the experimental results from Ref. [6]. To simplify the
problem, as before, we assumed that no population is coherently transferred to
the attractive branch, i.e., $n_{\mathrm{a}}(t,{\bf k})=0$. Neglecting the
decay to the attractive branch, the decoherence rate of the repulsive polarons
$\Gamma_{\mathrm{r}}^{\mathrm{dec}}$ is not noticeable on the timescale
considered. In Fig. 3, we observe indeed coherent Rabi oscillations given by
the renormalized frequency $\sqrt{Z_{\mathrm{r}}}\Omega$ without any
significant decay of the signal, in good agreement with the experimental data.
For these parameter regimes, the uncertainty in the temperature (shaded
region) is almost unnoticeable.

For the parameters considered here, the collision integral plays only a minor
part, and neglecting it is a very good approximation for the observables shown
in Fig. 3, which are dominated by the coherent Rabi transfer between the
states and the decoherence rate $\Gamma_{\alpha}^{\mathrm{dec}}$. The
smallness of the collision integral also implies slow thermalization of the
impurity.

Interestingly, as pointed out in the experiment [6], it appears that the
theoretically predicted strong decay of the repulsive polaron to the
attractive branch does not play a role in the Rabi dynamics (see also the
discussion in [15]).

Figure 4: The map of the violation of the detailed balance quantified by the
parameter $B(k_{1},k_{2})$ from Eq. (81). The dashed lines mark the contours
where $B=0.9$, 1, 1.1; notice that $B=1$ also on the diagonal $k_{1}=k_{2}$.
The interaction strengh is $1/(k_{F}a)=-0.25$ (left column), 0 (middle
column), +0.25( right column). Each row corresponds to a minimum (the first
and third rows) or maximum (the second and fourth rows) of the Rabi
oscillations shown in Fig. 2 of the main text. The parameters:
$\Omega=0.7\varepsilon_{F}$, $T=0.135\varepsilon_{F}$,
$\delta_{\mathrm{a}}(0)=0$.

##  VI Thermalization and bound on the violation of detailed balance

In this Section we provide additional results for the time-dependent violation
of the detailed balance in the approach of the polaron to a thermal state. We
also derive a bound on the average detailed balance violation given by the
Kullback-Leibler divergence as mentioned in the main text. Finally, we also
discuss the shape of the distributions of the polarons.

###  VI.1 Detailed balance violation

The detailed balance in the main text is probed by the following quantity:

$B\equiv\frac{1-n_{\mathrm{a}}(\mathbf{k})}{n_{\mathrm{a}}(\mathbf{k})}\frac{1-n_{1}(\mathbf{k}_{1})}{n_{1}(\mathbf{k}_{1})}\frac{n_{\mathrm{a}}(\mathbf{k}^{\prime})}{1-n_{\mathrm{a}}(\mathbf{k}^{\prime})}\frac{n_{1}(\mathbf{k}_{1}^{\prime})}{1-n_{1}(\mathbf{k}_{1}^{\prime})}.$
(80)

If the environment $n_{1}$ is in equilibrium, this object simplifies to the
following form, which depends only on two momenta:

$B(\mathbf{k},\mathbf{k}^{\prime})=\frac{f(\mathbf{k})}{f(\mathbf{k}^{\prime})},$
(81)

where, defining $\beta=1/T$,

$f(\mathbf{k})=\frac{1-n_{\mathrm{a}}(\mathbf{k})}{n_{\mathrm{a}}(\mathbf{k})}e^{-\beta\varepsilon_{\mathrm{a}}(\mathbf{k})}.$
(82)

In this notation, we dropped the time dependence of $n_{\mathrm{a}}$ for
brevity. If the polarons are in a thermal state,
$B(\mathbf{k},\mathbf{k}^{\prime})\equiv 1$ for all momenta, and the deviation
from 1 signals the violation of detailed balance. If, on the other hand,
$B(\mathbf{k},\mathbf{k}^{\prime})>1$, the rate of the _in_ scattering
processes, which populate the mode $\mathbf{k}$ as a result of collisions of
polarons with momentum $\mathbf{k}^{\prime}$, is larger than the _out_
processes, in which polarons with momentum $\mathbf{k}$ are scattered out to
momentum $\mathbf{k}^{\prime}$.

In Fig. 4, we plot $B(\mathbf{k}_{1},\mathbf{k}_{2})$ for times corresponding
to extrema of Fig. 2 in the main text for $1/(k_{F}a)=$ -0.25 (left column), 0
(middle column), 0.25 (right column). The dashed lines mark the values 0.9,
1.0 (notice also $B(k_{1},k_{1})=1$ on the diagonal), 1.1. Here, we observe
that, for the times corresponding to the minima of Rabi oscillations, when the
polaron population is maximal, the region in momentum space (see the first and
third rows), where the detailed balance is $B\approx 1$, is large. This is
also manifested as minima of the main panel in Fig. 3 in the main text. On the
contrary, a large region with significant violation of the detailed balance is
seen for times when the population is very small (see the second and fourth
rows), i.e., at the maxima of the oscillations.

From Fig. 4, in the second and fourth rows, which correspond to minima of the
polaron densities in the Rabi oscillations, a simple structure forms. Namely,
for all $k_{2}<k_{1}$ we have $B(k_{1},k_{2})<1$, which means that the _out_
scattering from $k_{1}$ to $k_{2}<k_{1}$ overcomes the _in_ collisions, and we
observe an effective cooling of the system. This effect is also visible in the
inset of Fig. 3 of the main text, where the inferred temperature $T_{K}$
decreases at times corresponding to maxima of the Rabi oscillations (see Fig.
2 in the main text).

If the system is close to thermal equilibrium, $B\approx 1$ and the deviation
$B-1$ quantifies the degree of detail balance violation. In the main text, we
average this deviation over the system, to obtain a single measure of the
violation. Thus, we have

$(\Delta B)^{2}\equiv\langle[B(\mathbf{k},\mathbf{k}^{\prime})-1]^{2}\rangle,$
(83)

where the average is taken over the normalized distribution $\propto
n_{\mathrm{a}}(\mathbf{k})n_{\mathrm{a}}(\mathbf{k}^{\prime})$. Since Eq. (81)
is a product of two factors, we obtain:

$(\Delta B)^{2}=1+\langle f(\mathbf{k})^{2}\rangle\langle
f(\mathbf{k})^{-2}\rangle-2\langle f(\mathbf{k})\rangle\langle
f(\mathbf{k})^{-1}\rangle,$ (84)

where $\langle\cdot\rangle\equiv C\sum_{\mathbf{k}}(\cdot
n_{\mathrm{a}}(\mathbf{k}))$, and $C$ is the normalization constant.

Now, we assume that $n_{\mathrm{a}}(\mathbf{k})=n_{T}(\mathbf{k})+\delta
n_{\mathbf{k}}$, where
$n_{T}=1/(e^{\beta(\varepsilon_{\mathrm{a}}(\mathbf{k})-\mu_{a})}+1)$ is a
thermal state, and the chemical potential $\mu_{a}$ is adjusted in order to
match the density of the polarons.

Expanding now Eq. (84) to second order in $\delta n_{\mathbf{k}}$, we obtain

$\displaystyle\langle f(\mathbf{k})^{2}\rangle$ $\displaystyle=$
$\displaystyle\langle
1+2f_{1}(\mathbf{k})+(f_{1}^{2}(\mathbf{k})+2f_{2}(\mathbf{k}))\rangle,$ (85)
$\displaystyle\langle f(\mathbf{k})^{-2}\rangle$ $\displaystyle=$
$\displaystyle\langle
1-2f_{1}(\mathbf{k})+(3f_{1}^{2}(\mathbf{k})-2f_{2}(\mathbf{k}))\rangle,$ (86)
$\displaystyle\langle f(\mathbf{k})\rangle$ $\displaystyle=$
$\displaystyle\langle 1+f_{1}(\mathbf{k})+f_{2}(\mathbf{k})\rangle,$ (87)
$\displaystyle\langle f(\mathbf{k})^{-1}\rangle$ $\displaystyle=$
$\displaystyle\langle
1-f_{1}(\mathbf{k})+(f_{1}^{2}(\mathbf{k})-f_{2}(\mathbf{k}))\rangle,$ (88)

where $f_{1}(\mathbf{k})=-(\delta
n_{\mathbf{k}}e^{-\beta\varepsilon_{\mathrm{a}}(\mathbf{k})})/n_{T}^{2}(\mathbf{k})$
and $f_{2}(\mathbf{k})=(\delta
n_{\mathbf{k}}^{2}e^{-\beta\varepsilon_{\mathrm{a}}(\mathbf{k})})/n_{T}^{3}(\mathbf{k})$

Up to the second order in $\delta n_{\mathbf{k}}$ we obtain

$(\Delta B)^{2}=2\big{[}\langle f_{1}^{2}\rangle-\langle
f_{1}\rangle^{2}\big{]}\geqslant 0.$ (89)

Since we can drop the second term in order to get the bound $\Delta
B^{2}\leqslant 2\langle f_{1}^{2}\rangle$, we have

$\langle
f_{1}^{2}\rangle=C\sum_{\mathbf{k}}n_{\mathrm{a}}(\mathbf{k})\bigg{(}-\frac{\delta
n_{\mathbf{k}}e^{-\beta\varepsilon_{\mathrm{a}}(\mathbf{k})}}{n_{T}^{2}(\mathbf{k})}\bigg{)}^{2}\approx
C\sum_{\mathbf{k}}\frac{\delta
n_{\mathbf{k}}^{2}}{[1-n_{T}(\mathbf{k})]^{2}n_{T}(\mathbf{k})},$ (90)

where we replaced $n_{\mathrm{a}}$ by $n_{T}$, since we are concerned with the
second order in $\delta n_{\mathbf{k}}$.

Now, we use
$1/[n_{T}(1-n_{T})^{2}]=(1+e^{-\beta\varepsilon_{\mathrm{a}}})^{2}/n_{T}$, and
we obtain:

$\langle f_{1}^{2}\rangle=C\sum_{\mathbf{k}}\frac{\delta
n_{\mathbf{k}}^{2}}{n_{T}(\mathbf{k})}\bigg{(}1+e^{-\beta\varepsilon_{\mathrm{a}}(\mathbf{k})}\bigg{)}^{2}\leqslant\bigg{(}1+e^{-\beta\varepsilon_{\mathrm{a}}(0)}\bigg{)}^{2}C\sum_{\mathbf{k}}\frac{\delta
n_{\mathbf{k}}^{2}}{n_{T}(\mathbf{k})}.$ (91)

Finally, since up to the second order, the Kullback-Leibler divergnce has the
form $D(\tilde{n}_{\mathrm{a}}|\tilde{n}_{T})\approx C\sum_{\mathbf{k}}\delta
n_{\mathbf{k}}^{2}/(2n_{T}(\mathbf{k}))$, here tilde indicates the
normalization, we obtain

$(\Delta B)^{2}\leqslant
4\bigg{(}1+e^{-\beta\varepsilon_{\mathrm{a}}(0)}\bigg{)}^{2}D(\tilde{n}_{\mathrm{a}}|\tilde{n}_{T}),$
(92)

which has the form mentioned in the main text, i.e.,
$\langle(B-1)^{2}\rangle\leqslant\zeta
D(\tilde{n}_{\mathrm{a}}|\tilde{n}_{T})$ with
$\zeta=4(1+e^{-\beta\varepsilon_{\mathrm{a}}(0)})^{2}>0$.

Figure 5: The numerically determined occupations of
$n_{\mathrm{a}}(t,\mathbf{k})$ (solid) and the inferred equilibrium
distribution $n_{\mathrm{a}}^{\mathrm{eq}}(T_{K}(t)),\mathbf{k})$ (dashed) as
a function of momentum $k=|\mathbf{k}|$ (in the units of $k_{F}$) for various
times corresponding to the points indicated in the insets. The color code of
the curves in the main panel is the same as the times corresponding to the
colors of the points in the inset. Main panel: from top to bottom:
$1/(k_{F}a)=$ -0.25 (panel a), 0 (b), 0.25 (c). The inset: shows the inferred
temperature $T_{K}$ as a function of time $t$ (in units of $t_{F}$) as
described in the main text. The parameters: $\Omega=0.7\varepsilon_{F}$,
$T=0.135\varepsilon_{F}$.

###  VI.2 Occupation of attractive polarons

In Fig. 5, we show the occupations of the momenta of the attractive polaron
branch, i.e., $n_{\mathrm{a}}(t,\mathbf{k})$, for $1/(k_{F}a)=$ -0.25 (panel
a), 0 (b), 0.25 (c) at some characteristic times $t$. We have chosen the times
$t$ corresponding to the minima (red and pink points in the insets),
$T_{K}/T=1$ (orange and blue) and maxima (green and purple) of the inferred
temperature $T_{K}$ which varies with time as shown in the corresponding
insets. The solid curves show the results of the numerical simulations of
$n_{\mathrm{a}}(t,\mathbf{k})$ as a function of $k=|\mathbf{k}|$ while the
dashed curves are the optimal thermal distribution functions
$n_{\mathrm{a}}^{\mathrm{eq}}(T_{K},\mathbf{k})$ corresponding to the
optimized temperatures $T_{K}$.

The presented results show that the overall shape of the distributions
$n_{\mathrm{a}}(t,\mathbf{k})$ is captured by the found equilibrium
distribution, which explains the fact that the inferred $T_{K}$ is of the
correct order as $T$. All the nearest equilibrium distributions are, however,
overestimated at low momenta [apart from the dip in the distribution for the
green curve in (b) corresponding to the first maximum of $T_{K}/T$]. This low
momentum region, however, is weakly represented in the Kullback-Leibler
divergence due to the jacobian $k^{2}$ in the integration. The worst agreement
is found for low densities of the polarons, which correspond to the peaks of
$T_{K}$ (green and purple points in the insets), since the distribution is
very far from equilibrium. Here, we note that the non-Boltzmann-like shape of
the distributions confirms the need for the non-equilibrium theory of the
Rabi-coupled dynamics of polarons.

Finally, we comment on the timescales of the approach to the steady state.
Since it is challenging to extract its form analytically due to non-trivial
dependence of the $T$-matrix on momentum $\mathbf{k}$, we resorted to its
simple estimation from the dynamics of the violation of detailed balance.
Assuming an exponential decay, we fitted an exponential function $\propto
e^{-t/\tau_{\mathrm{th}}}$ to the peaks that are shown in Fig. 3 in the main
text. We have found $\tau_{\mathrm{th}}=$ 22 $t_{F}$ for $1/(k_{F}a)=$ -0.25,
29 $t_{F}$ for $1/(k_{F}a)$ = 0 and 14 $t_{F}$ for $1/(k_{F}a)$ = +0.25. The
markedly shorter time for $1/(k_{F}a)$ = +0.25 captures the behaviour of the
red dotted curve from Fig. 3 (main text, main panel) which approaches zero on
a faster timescale. Interestingly, these simple results show non-monotonic
behaviour of thermalization across the unitary limit $1/(k_{F}a)=0$, which is
also visible in Fig. 3 of the main text.

## References

  * Kamenev [2011] A. Kamenev,  _Field theory of non-equilibrium systems_ (Cambridge University Press, 2011) pp. 208–209.
  * Wasak  _et al._ [2021a] T. Wasak, R. Schmidt, and F. Piazza, Physical Review Research 3, 013086 (2021a).
  * Wasak  _et al._ [2021b] T. Wasak, F. Pientka, and F. Piazza, arXiv preprint arXiv:2103.14040 (2021b).
  * Sighinolfi [2022] M. Sighinolfi,  _Open Quantum Systems and Ultracold atoms_ (2022).
  * Larkin and Ovchinnikov [1975] A. Larkin and Y. Ovchinnikov, Sov. Phys. JETP 41, 960 (1975).
  * Scazza  _et al._ [2017] F. Scazza, G. Valtolina, P. Massignan, A. Recati, A. Amico, A. Burchianti, C. Fort, M. Inguscio, M. Zaccanti, and G. Roati, Phys Rev. Lett. 118, 083602 (2017).
  * Chevy [2006] F. Chevy, Phys. Rev. A 74, 063628 (2006).
  * Combescot  _et al._ [2007] R. Combescot, A. Recati, C. Lobo, and F. Chevy, Physical review letters 98, 180402 (2007).
  * Massignan  _et al._ [2014] P. Massignan, M. Zaccanti, and G. M. Bruun, Reports on Progress in Physics 77, 034401 (2014).
  * Sieberer  _et al._ [2016] L. M. Sieberer, M. Buchhold, and S. Diehl, Reports on Progress in Physics 79, 096001 (2016).
  * Note [1] Even if this time dependence is included, it does not change the final kinetic equations because it will be removed by the procedure adopted.
  * Ruckenstein and Lévy [1989] A. E. Ruckenstein and L. P. Lévy, Phys. Rev. B 39, 183 (1989).
  * Kadanoff and Baym [2018] L. P. Kadanoff and G. Baym,  _Quantum statistical mechanics: Green’s function methods in equilibrium and nonequilibrium problems_ (CRC Press, 2018) pp. 102–120.
  * Knoll  _et al._ [2001] J. Knoll, Y. Ivanov, and D. Voskresensky, Annals of Physics 293, 126 (2001).
  * Adlong  _et al._ [2020] H. S. Adlong, W. E. Liu, F. Scazza, M. Zaccanti, N. D. Oppong, S. Fölling, M. M. Parish, and J. Levinsen, Phys. Rev. Lett. 125, 133401 (2020).