Dynamic transition of the generalized Jaynes–Cummings model: multi-particles and inter-particle interaction effects

Quantum dynamics lies at the center of the research of quantum science; compared with classical dynamics, its peculiarity offers novel prospective applications in quantum technologies Galindo and Martín-Delgado (2002); Degen et al. (2017); Pezzè et al. (2018); Amico et al. (2022). However, since experiments inevitably subject quantum systems to environments, it is crucial to investigate environmental effects on dynamic processes such as decoherence, dissipation, and entanglement Dobrovitski et al. (2003); Yao et al. (2007); Hanson et al. (2008); Orth et al. (2010); Strathearn et al. (2018); De Filippis et al. (2021), which often constrain the robustness of intended quantum operations.

A prototype model addressing the problem of decoherence due to a dissipative environment is the spin-boson model Leggett et al. (1987); Weiss (2012), whose Hamiltonian is usually given by $H_{\rm SB}=\Delta\sigma_{x}/2+\sigma_{z}\sum_{\nu}\lambda_{\nu}(a_{\nu}+a^{\dagger}_{\nu})+\sum_{\nu}\omega_{\nu}a_{\nu}^{\dagger}a_{\nu}$. This model derives from the familiar problem of a single particle tunneling between the two minima in a double well potential Weiss et al. (1987); Jelic and Marsiglio (2012). The eigenstates of $\sigma_{z}$, $\ket{\uparrow}$ and $\ket{\downarrow}$, correspond to the left and right minimum of the potential respectively. The off-diagonal term $\sim\Delta\sigma_{x}$ brings about the oscillation between the two eigenstates, equivalent to the particle tunneling between the two minima. The dissipative environment is modeled by the boson modes whose creation (annihilation) operators are $a_{\nu}^{\dagger}$ ($a_{\nu}$). The coupling between the spin and the bosons indicates that the environment constantly “monitors” the state (position) of the spin (particle). Thus it is anticipated if such a monitoring is strong enough, the environment collapses the spin state into either $\ket{\uparrow}$ or $\ket{\downarrow}$, and wipes out the quantum coherence.

Actually it has been shown in Ref. Leggett et al. (1987) that the environmental effects on the spin dynamics are encapsulated in the spectral function $J(\omega)$ of the bosons; at zero temperature, for $\Delta\to 0$, the spin dynamics is localized in the eigenstates of $\sigma_{z}$ for sub-Ohmic $J(\omega)$, and undergoes a damped oscillation for super-Ohmic $J(\omega)$, and transits from a damped oscillation to an incoherent relaxation and to localization for Ohmic $J(\omega)$ with increasing magnitude. Furthermore, for finite $\Delta$, numerical renormalization group calculations observed coherent dynamics even for sub-Ohmic $J(\omega)$ Bulla et al. (2003); Anders et al. (2007), and the numerically exact time-evolving matrix product operator method identified a new phase characterized by an aperiodic behavior for strong coupling Otterpohl et al. (2022). The spin-boson model not only sheds light on the fundamental question of how quantum phenomena transit to classical behavior Breuer and Petruccione (2002); Ingold (2002), but also provides a base to study topics ranging from electron transport to quantum information Zurek (2003); Kapral (2015); De La Torre et al. (2021); Cai et al. (2023); Hangleiter and Eisert (2023); Schuster and Yao (2023).

However, it is known that different specific coupling forms between systems and environments can lead to drastically different behaviors in system dynamics Breuer and Petruccione (2002); Zhang and Yu (2023). Other than the spin-boson model, another widely studied model is the Jaynes–Cummings model generalized to a continuum boson modes. The two models are related in the sense that the Hamiltonian of the latter is $\sigma_{y}H_{\rm SB}\sigma_{y}$ with only the system-environment coupling terms $\sim\sigma^{+}a_{\nu}$ and $\sigma^{-}a_{\nu}^{\dagger}$ retained. It can be shown analytically that if initially the particle is in its excited state and the bosons are in their vacuum state, the time evolution of the probability that the particle continues occupying the excited state exhibits a dynamic transition Pfeifer (1982); Cohen-Tannoudji et al. (1998): when the system-environment coupling is weak, the probability decays to zero monotonically; when the coupling is strong, a finite weight of the particle is localized in the excited state. The generalized Jaynes–Cummings model has been employed to explain stability of certain molecule states Pfeifer (1982); Bahrami and Bassi (2011); Taher Ghahramani and Shafiee (2013), and to explore non-Markovian behavior for a particle embedded in a boson bath John and Quang (1994); Florescu and John (2001); Burgess and Florescu (2022). However, our knowledge of the model system when applied to multiple particles with interaction included is rather limited.

In this work, we study the dynamic transition of the generalized Jaynes–Cummings model applied to multiple particles and a continuum boson bath; we access the effects of an all to all Ising type inter-particle interaction on the transition [see Eq. (1)]. Recently, this interaction is of particular interest partly due to its application in realizing the CNOT gate in quantum computation Mermin (2007). Our results are calculated by the non-perturbative time-dependent numerical renormalization group (TD-NRG) method Anders and Schiller (2006). Our numerical calculation agrees well with the benchmark provided by the analytic result for a single particle Pfeifer (1982); Cohen-Tannoudji et al. (1998). In the case of the number of particles $N\geq 2$, we find that even in the absence of the inter-particle interaction, the critical value of the system-environment coupling required for the transition decreases with larger $N$. Moreover, we find that the anti-ferromagnetic (ferromagnetic) Ising interaction suppresses (enlarges) the critical coupling value. Our results enrich the understanding of environmental effects on systems with intra-interactions.

We consider $N$ identical two-level particles coupled to a continuum of boson modes. The energy difference $E_{a}$ between the two internal levels of each particle $|e\rangle$ and $|g\rangle$ is close to the boson mode frequencies $\{\omega_{\nu}\}$. We take $\hbar=1$ throughout. We denote $\omega_{l}$ and $\omega_{c}$ as the lower bound and the upper bound cut-off of $\{\omega_{\nu}\}$. It is convenient for us to choose $\omega_{l}$ as the zero point for energyies. We consider an all to all Ising type interaction between the particles. This interaction can help to realize the CNOT gate in quantum computation Mermin (2007). The Hamiltonian of the combined system is given by

	$\displaystyle\tilde{H}=$	$\displaystyle\sum_{j=1}^{N}\Bigl{[}\frac{\Delta}{2}\sigma_{j}^{z}+\sum_{\nu}\lambda_{\nu}(\sigma_{j}^{-}a_{\nu}^{\dagger}+\sigma_{j}^{+}a_{\nu})\Bigr{]}+\sum_{\nu}\tilde{\omega}_{\nu}a_{\nu}^{\dagger}a_{\nu}$
		$\displaystyle+g\sum_{j<k}\sigma_{j}^{z}\sigma_{k}^{z},$		(1)

where $\Delta\equiv E_{a}-\omega_{l}$ and $\tilde{\omega}_{\nu}\equiv\omega_{\nu}-\omega_{l}$, and $\sigma_{j}^{z}=(|e_{j}\rangle\langle e_{j}|-|g_{j}\rangle\langle g_{j}|)/2$, $\sigma_{j}^{+}=|e_{j}\rangle\langle g_{j}|$ and $\sigma_{j}^{-}=|g_{j}\rangle\langle e_{j}|$, and $a_{\nu}^{\dagger}(a_{\nu})$ are creation (annihilation) operators of the boson modes with frequency $\omega_{\nu}$, $\lambda_{\nu}$ are the couplings between the particles and the boson modes, and $g$ is the interatomic interaction coupling. The Hamiltonian, Eq. (1), is a generalization of the Jaynes–Cummings model model to multiple particles and boson modes, plus the interatomic interaction introduced. We take the boson modes as a bath.

As in the case for the renowned spin-boson model, the effects of the boson bath on the particle dynamics take place via the spectral distribution Leggett et al. (1987); Weiss (2012)

$\displaystyle J(\omega)=\pi\sum_{\nu}\lambda_{\nu}^{2}\delta(\omega-\tilde{\omega}_{\nu});$

(2)

this can be clearly seen if one integrates out the boson bath in the path integral formalism. We consider the continuum limit such that $J(\omega)$ is a smooth function. To be specific, we assume the spectral distribution having the form

$\displaystyle J(\omega)=\begin{cases}2\pi\alpha(\omega_{c}-\omega_{l})^{1-s}\omega^{s},\,{\rm for}\,0<\omega<\omega_{c}-\omega_{l};\\ 0,{\rm otherwise}.\end{cases}$

(3)

The overall coupling of the atoms to the boson bath is characterized by the dimensionless strength $\alpha$; this form of $J(\omega)$ has been widely studied in the spin-boson model Leggett et al. (1987); Weiss (2012), in the context of which $s=1$ is called Ohmic with $s>1$ and $s<1$ called super-Ohmic and sub-Ohmic respectively. We focus our attention on the case $s=1$.

We are interested in the time evolution of $P_{e}(t)\equiv\sum_{j=1}^{N}\langle\sigma_{j}^{z}(t)+1\rangle/2N$ starting from the initial state that all the particles are in the excited state and the boson bath is not excited, i.e., $\ket{\psi(0)}=\prod_{j}\ket{e}_{j}\otimes\ket{0}$ state, where $\ket{0}$ is the vacuum state of the boson bath. We employ the non-perturbative time-dependent numerical renormalization group (TD-NRG) method to calculate $P_{e}(t)$ numerically Anders and Schiller (2006), and demonstrate that there exists a critical value $\alpha_{c}$ such that for long time, $P_{e}(t)$ decays to zero for $\alpha<\alpha_{c}$, and converges to a nonzero value for $\alpha>\alpha_{c}$. We show how $\alpha_{c}$ changes with the particle number $N$ and the interaction strength $g$.

Our calculation follows the algorithm of TD-NRG prescribed for boson baths Bulla et al. (2005). In the procedure of logarithmically discretizing the continuum boson spectrum, we take the discretization parameter $\Lambda=1.1$ and the $z$-trick parameter $N_{z}=4$. We use up to $N_{B}=200$ discretized boson modes and keep up to $N_{S}=1000$ lowest energy states in each TD-NRG iteration. For numerics, we use $\omega_{c}-\omega_{l}$ as the unit for energies and take $\omega_{c}-\omega_{l}$ to be unity, and set $\Delta=0.05$ ($\ll\omega_{c}-\omega_{l}$).

First we domenstrate that our numerical calculation agrees with the benchmark provided by the analytic solution for $N=1$. The initial state of the whole system is prepared as $\ket{\psi(0)}=\ket{e}\otimes\ket{0}$. Under the Hamiltonian (1), $\ket{\psi(0)}$ is coupled to the continuum of states $\{\ket{g}\otimes a_{\nu}^{\dagger}\ket{0}\}$.

Figure 1 shows the numerical results of $P_{e}(t)$, which by themselves indicate a transition in the dynamic behavior of $P_{e}(t)$ as $\alpha$ varies. The transition critical point turns out to occur at $\alpha_{c}=0.025$ for the numerics we take (see below). For $\alpha<\alpha_{c}$, Fig. 1 (a) and (b) show that $P_{e}(t)$ decays towards zero monotonically over time. When $\alpha$ is sufficiently small, one assumes the perturbation theory applicable and expects an exponential decay of $P_{e}(t)$ over time. As $\alpha$ increases further approaching $\alpha_{c}$, an early period in which $1-P_{e}(t)\sim t^{2}$ becomes more evident, and an exponential decay follows afterwards. When $\alpha$ is above $\alpha_{c}$, Fig. 1 (c) and (d) show that $P_{e}(t)$ behaves totally differently: $P_{e}(t)$ does not look to decay to zero any more for long time; phenomenologically it seems that a fraction of the particle’s weight is trapped in the excited state. Furthermore, an attenuating long time scale oscillation develops on top of the general trend of $P_{e}(t)$. This oscillation is more obvious with increasing $\alpha$.

The features of the above numerical results can be understood in an analytic way. For the present problem of a single particle coupled to the continuum bath, there exists an analytic expression of $P_{e}(t)$ Cohen-Tannoudji et al. (1998). Figure 1 shows that our numerical results agree well with the analytic calculation; the agreement benchmarks our numerical implementation of the TD-NRG algorithm.

As $P_{e}(t)=|\mathcal{U}_{e}(t)|^{2}$ with $\mathcal{U}_{e}(t)\equiv\bra{\psi(0)}e^{-i\tilde{H}t}\ket{\psi(0)}$, by the method of the Green’s function Cohen-Tannoudji et al. (1998), one can derive the Fourier transform $\mathcal{U}_{e}(\omega)$, defined by

$\displaystyle\ \mathcal{U}_{e}(t)=\int_{-\infty}^{\infty}d\omega\ \mathcal{U}_{e}(\omega)\ e^{-i\omega t},$

(4)

having the form

$\displaystyle\mathcal{U}_{e}(\omega)=\lim_{\eta\to 0^{+}}\frac{1}{\pi}\ \frac{J(\omega)+\eta}{[\omega-\Delta-2\alpha\Delta_{e}(\omega)]^{2}+[J(\omega)+\eta]^{2}}.$

(5)

Here $J(\omega)$ takes the role of the imaginary part of the self-energy. Since we take Eq. (3) with $s=1$ for $J(\omega)$, the real part of the self-energy

$\displaystyle\Delta_{e}(\omega)\equiv\mathcal{P}\int_{-\infty}^{\infty}\frac{d\omega^{\prime}}{2\alpha\pi}\frac{J(\omega^{\prime})}{\omega-\omega^{\prime}},$

(6)

is given by

$\displaystyle\Delta_{e}(\omega)=-1-\omega\ln\left|\frac{1-\omega}{\omega}\right|.$

(7)

Here $\mathcal{P}$ stands for principal value. Note due to the discontinuity of $J(\omega)$ at $\omega=1$, $\Delta_{e}(\omega)$ diverges correspondingly there.

Combining Eqs. (4) and (LABEL:uo), one can see that when $\alpha$ is sufficiently small, the contribution to Eq. (4) is mainly from the frequency domain where $\omega-\Delta-2\alpha\Delta_{e}(\omega)\approx 0$. Figure 2 shows the curve $\Delta_{e}(\omega)$ together with the straight line $(\omega-\Delta)/2\alpha$ of various $\alpha$ and correspondingly their intersections at frequency denoted by $\omega_{m}$, which is smaller than $\Delta$. Thus, for $\alpha\to 0$, the domain $|\omega-\Delta|\lesssim J(\Delta)$ dominates the contribution to Eq. (4), and one can approximate $\mathcal{U}_{e}(\omega)=\frac{1}{\pi}\ \frac{J(\Delta)}{[\omega-\Delta-2\alpha\Delta_{e}(\Delta)]^{2}+J^{2}(\Delta)}$ as a Lorentzian, and obtain $P_{e}(t)=e^{-2J(\Delta)t}$, an exponential decay as shown in Fig. 1 (a). Note although the straight line also intersects with $\Delta_{e}(\omega)$ at two other points close to $\omega=1$ for $\alpha\ll 1$, their contributions are negligible since $\Delta_{e}(\omega)\to-\ln|\omega-1|$ and $\partial_{\omega}\Delta_{e}$ is exponentially large there [also see Eq. (LABEL:uem)].

The deviation of the full $\mathcal{U}_{e}(t)$ from its Lorentzian approximation result shall be most noticable both for small and large $t$ Cohen-Tannoudji et al. (1998). Since for small $t$, if one expands $e^{-i\omega t}$ to the second order of $t$ in Eq. (4), one would conclude that $1-P_{e}(t)\sim t^{2}$. This discrepancy is due to the fact that the full $\mathcal{U}_{e}(\omega)$ decays faster than $|\omega|^{-2}$ for $|\omega|\to\infty$ as $J(\omega)$ also goes to zero in this limit. The insets of Fig. 1 (a) and (b) show that for small $t$, $P_{e}(t)$ does decrease quadratically, and this discrepency is, as expected, more obvious when $\alpha$ increases whereas the Lorentzian broadens. For long $t$, the difference between the full $\mathcal{U}_{e}(\omega)$ and its Lorentzian approximation at $\omega\to 0^{+}$ would matter. Note $\lim_{\omega\to 0^{+}}J(\omega)\to 0^{+}$. Our numerical results shown in Fig. 1 (a) and (b) seem not in that long $t$ regime yet.

A qualitative change occurs at the critical point $\alpha_{c}\equiv-\Delta/2\Delta_{e}(0)$ ($\alpha_{c}=0.025$ for $\Delta=0.05$); beyond this point, the straight line intersects $\Delta_{e}(\omega)$ at a negative frequency out of the range of the continuum. This intersection point corresponds to a discrete eigenstate with eigen-energy $\omega_{m}$ in the coupled system (see Appendix for details). Physically, the coupling has modified the initial discrete state of energy $\Delta$ so much that a dressed discrete state emerges below the continuum. This dressed discrete state is stable because it is not coupled to continuum states any more.

The emergence of the dressed discrete state brings about $\lim_{t\to\infty}P_{e}(t)\neq 0$. Since now for $\omega<0$, we have

$\displaystyle\mathcal{U}_{e}(\omega)=\frac{1}{|1-2\alpha\Delta_{m}^{{}^{\prime}}|}\delta(\omega-\omega_{m}),$

(8)

where $\Delta_{m}^{{}^{\prime}}\equiv\partial_{\omega}\Delta_{e}|_{\omega=\omega_{m}}$. This delta function transforms to an undamped term, i.e., $\frac{1}{|1-2\alpha\Delta_{m}^{{}^{\prime}}|}e^{-i\omega_{m}t}$, in $\mathcal{U}_{e}(t)$. Compared with Fig. 3 (a), Fig. 1 (c) shows that $P_{e}(t)$ does converge to $\frac{1}{|1-2\alpha\Delta_{m}^{{}^{\prime}}|^{2}}$ for long time.

The attenuating long time scale oscillation exhibited in $P_{e}(t)$ can be attributed to the interference between the contributions from Eq. (LABEL:uem) and the superposition of frequencies $\omega>0$ in Eq. (4). The oscillation period in Fig. 1 (c) for $\alpha=0.03$ and $0.04$ looks close to $1/|\omega_{m}|$ given in Fig. 3 (b). The nature of this oscillation can be readily elucidated in the large $\alpha$ limit, where the leading term in $\tilde{H}$ becomes the coupling Hamiltonian

$\displaystyle H_{a-b}=\sqrt{\sum_{k}\lambda_{k}^{2}}(\sigma^{-}A^{\dagger}+\sigma^{+}A),$

(9)

with $A\equiv\frac{1}{\sqrt{\sum_{k}\lambda_{k}^{2}}}\sum_{k}\lambda_{k}a_{k}$ and $[A,A^{\dagger}]=1$. In this limit, the initial state $\ket{\psi(0)}=\ket{e}\otimes\ket{0}$ is only coupled to another state $\ket{g}\otimes A^{\dagger}\ket{0}$ by $H_{a-b}$. Diagonalization of $H_{a-b}$ in the subspace spanned by these two states yields two eigen-energies $\pm\sqrt{\sum_{k}\lambda_{k}^{2}}$, which shall give rise to an oscillation in $P_{e}(t)$ of frequency $2\sqrt{\sum_{k}\lambda_{k}^{2}}$. By Eq. (2), for the form we assume for $J(\omega)$ with $s=1$, $\sum_{k}\lambda_{k}^{2}=\int d\omega^{\prime}J(\omega^{\prime})/\pi=\alpha$. These two eigen-energies correspond to the two intersections between the straight line and $\Delta_{e}(\omega)$ at $|\omega|\gg 1$ in the limit $\alpha\to\infty$ (see Fig. 2); since for $|\omega|\gg 1$, by Eq. (6), $\Delta_{e}(\omega)\approx\int d\omega^{\prime}J(\omega^{\prime})/2\alpha\pi\omega=1/2\omega$, and resultantly the straight line $(\omega-\Delta)/2\alpha$ intersects $\Delta_{e}(\omega)$ at two large magnitude frequencies which are approximately $\pm\sqrt{\alpha}$. The terms in Eq. (1) other than $H_{a-b}$ can be further treated perturbatively. Thus the oscillation must persist for a reasonably long time. Figure 1 (c) and (d) indicate that the attenuating long time scale oscillation sets in after passing the critical point and develops with increasing frequency all the way up to the large $\alpha$ limit.

Next we consider two particles, i.e., Eq. (1) with $N=2$, and the initial state is $|ee\rangle\otimes|0\rangle$. As the two particles couple to the same bath in the identical way, even in the absence of the inter-atomic interaction, i.e., $g=0$, difference in $P_{e}(t)$ from the case of one single particle is evident in the numerical results plotted in Fig. (4). For rather weak coupling strengths $\alpha$, Fig. 4 (a) shows that there seems to be two segments of exponential decay in $P_{e}(t)$. The decay rate of the first segment is about two times that of one single particle, and the decay rate of the second segment is further increased. We attribute the first segment to the decay process $|\psi(0)\rangle=|ee\rangle\otimes|0\rangle\to\{|\psi_{\nu}\rangle\equiv\frac{1}{\sqrt{2}}(|ge\rangle+|eg\rangle)\otimes a_{\nu}^{\dagger}|0\rangle\}$; the Fermi golden rule gives the decay rate $\Gamma$ to a resonant mode $\omega_{\nu}\approx\Delta$ proportional to $|\langle\psi_{\nu}|H_{a-b}|\psi(0)\rangle|^{2}=2|\langle g|\otimes\langle 0|a_{\nu}H_{a-b}|e\rangle\otimes|0\rangle|^{2}$. The further increased decay rate of the second segment may be due to the fact that for the following decay process $|\psi_{\nu}\rangle\to\{|\phi_{\nu\mu}\rangle\equiv|gg\rangle\otimes a_{\mu}^{\dagger}a_{\nu}^{\dagger}|0\rangle/\sqrt{1+\delta_{\nu\mu}}\}$, the rate is proportional to $|\langle\phi_{\nu\mu}|H_{a-b}|\psi_{\nu}\rangle|^{2}\approx(1+\delta_{\nu\mu})|\langle\psi_{\mu}|H_{a-b}|\psi(0)\rangle|^{2}$; if two excitations in the bath populate the same boson mode, the Bose enhancement would give rise to an additional factor $2$.

Figure 4 (b) and (c) indicate that there is a transition point we locate at $\alpha_{c}^{(2,0)}=0.013$; similar to the case of a single particle, now for $\alpha<\alpha_{c}^{(2,0)}=0.013$, $P_{e}(t)$ decays to zero for long time, and for $\alpha>\alpha_{c}^{(2,0)}=0.013$, $P_{e}(t)$ converges to a nonzero value for long time. Note that $\alpha_{c}^{(2,0)}$ is about half of $\alpha_{c}=0.025$ for a single particle. As $\alpha$ continues to increase, we also observe a long time scale oscillation in $P_{e}(t)$ as shown in Fig. 4 (d). Compared with that of a single particle, now the oscillation looks as a beat between different frequencies. This feature can be understood again in the large $\alpha$ limit: we diagonalize the dominant coupling Hamiltonian $H_{a-b}$ and obtain three eigenvalues $0,\pm\sqrt{3}\alpha$; to lowest order of $1/\alpha$, $P_{e}(t)=2/3+4\cos(\sqrt{3}\alpha t)/9-\cos(2\sqrt{3}\alpha t)/9$, consisting of two frequencies.

In the presence of the interaction ($g\neq 0$) for $N=2$, Fig. 5 shows the numerical results of $P_{e}(t)$ for $\alpha=0.012$ and varying $g$. Prominently, increasing the anti-ferromagnetic Ising interaction ($g>0$) can push $P_{e}(t)$ to undergo the dynamic transition. We identify that the critical point occurs at $g=0.01$ for fixed $\alpha=0.012$. In another word, for $g=0.01$, the critical value of $\alpha$ becomes $\alpha^{(2,g)}_{c}=0.012$, which is smaller than $\alpha^{(2,0)}_{c}=0.013$ for $g=0$. An approximation equivalent to choosing certain diagrams for self-energies also predicts that $\alpha^{(2,0)}_{c}<\alpha_{c}$ and $\alpha^{(2,g)}_{c}$ decreases (increases) with positive (negative) $g$ (see Appendix for details).

Our results of the critical value $\alpha^{(N,g)}_{c}$ for different values of $N$ and $g$ are summarized in Fig. 6. We find that the critical value $\alpha_{c}^{(N,0)}$ decreases with $N$; for fixed $N$, the critical value $\alpha_{c}^{(N,g)}$ generally increases for negative $g$ and decreases with positive $g$.

We studied the dynamics of the Jaynes–Cummings model generalized to multiple particles and a continuum boson bath. We also introduced an all to all Ising type inter-particle interaction. The dynamics starts with all the particles in their excited state and the boson bath in its vacuum state. The observable is the probability $P_{e}(t)$ that the particles remain in their excited state. We demonstrated that $P_{e}(t)$ exhibits a dynamic transition between decaying to zero and converging to a nonzero value when the system-environment coupling is tuned. We found that the critical coupling value decreases with the number of the particles $N$, and is suppressed (enlarged) by the anti-ferromagnetic (ferromagnetic) interaction ($g>0$). Our calculation is implemented via the non-perturbative time-dependent numerical renormalization group method, and agree with the benchmark for the case $N=1$. Our results reveal how the number of particles and their intra-interaction affect the dynamic transition.

This work is supported by the National Natural Science Foundation of China Grant No. 12474270.

To zero order of $\alpha$, the energies of $|\psi(0)\rangle$, $\{|\psi_{\nu}\rangle\}$ and $\{|\phi_{\nu\mu}\rangle\}$ are $E^{(0)}_{ee}=\Delta+g$, $E^{(0)}_{\nu}=-g+\omega_{\nu}$ and $E^{(0)}_{\nu\mu}=-\Delta+g+\omega_{\nu}+\omega_{\mu}$ respectively. The Hamiltonian $H_{a-b}$ couples $|\psi(0)\rangle$ with the band of states $\{|\psi_{\nu}\rangle\}$, and the band of states $\{|\psi_{\nu}\rangle\}$ with that of $\{|\phi_{\nu\mu}\rangle\}$. As $g$ varies, the band bottom of $\{|\psi_{\nu}\rangle\}$ moves with respect to that of $\{|\phi_{\nu\mu}\rangle\}$ and $E^{(0)}_{ee}$.

In the case of a single atom coupled to the continuum bath, the critical point is given by $\alpha_{c}\equiv-\Delta/2\Delta_{e}(0)$. We are going to show that the transition condition is the same as there is a stable dressed state whose energy $E$ touches the band bottom of the continuum bath, i.e., $E=-\Delta/2$. We expand the dressed state as

$\displaystyle|D_{1}\rangle=c_{e}|e\rangle\otimes|0\rangle+\sum_{\nu}c_{\nu}|g\rangle\otimes a_{\nu}^{\dagger}|0\rangle.$

(10)

By requiring $\tilde{H}|D_{1}\rangle=E|D_{1}\rangle$, we have

	$\displaystyle Ec_{e}=$	$\displaystyle\frac{\Delta}{2}c_{e}+\sum_{\nu}\lambda_{\nu}c_{\nu},$		(11)
	$\displaystyle Ec_{\mu}=$	$\displaystyle\left(-\frac{\Delta}{2}+\omega_{\mu}\right)c_{\mu}+\lambda_{\mu}c_{e}.$		(12)

From Eq. (12), we express $c_{\mu}$ in terms of $c_{e}$ and substitute $c_{\mu}$ in Eq. (11); we obtain

	$\displaystyle E-\frac{\Delta}{2}=$	$\displaystyle\sum_{\mu}\frac{\lambda_{\mu}^{2}}{E+\Delta/2-\omega_{\mu}}$
	$\displaystyle=$	$\displaystyle\int\frac{d\omega^{\prime}}{\pi}\frac{J(\omega^{\prime})}{E+\Delta/2-\omega^{\prime}}$
	$\displaystyle=$	$\displaystyle 2\alpha\Delta_{e}(E+\Delta/2),$		(13)

which determines $E$. In the absence of coupling, $\lambda_{\mu}=0$, $E=\Delta/2$. In the presence of nonzero $\lambda_{\mu}$, the continuum bath dresses the bare state $|e\rangle\otimes|0\rangle$. The dressed state becomes stable again once its energy drops out of the spectrum of the states $\{|g\rangle\otimes a_{\nu}^{\dagger}|0\rangle\}$, i.e., $E\leq-\Delta/2$. Thus the critical point is when $E=-\Delta/2$, and correspondingly

$\displaystyle-\Delta=2\alpha_{c}\Delta_{e}(0).$

(14)

One can carry out the same analysis for two particles. Likewise, we expand the dressed state as

	$\displaystyle|D_{2}\rangle=$	$\displaystyle c_{e}|ee\rangle\otimes|0\rangle+\sum_{\mu}c_{\mu}\left|\frac{ge+eg}{\sqrt{2}}\right\rangle\otimes a_{\mu}^{\dagger}|0\rangle$
		$\displaystyle+\sum_{\mu\nu}d_{\mu\nu}|gg\rangle\otimes a_{\mu}^{\dagger}a_{\nu}^{\dagger}|0\rangle.$		(15)

Note by commutation $d_{\mu\nu}=d_{\nu\mu}$. From $\tilde{H}|D_{2}\rangle=E|D_{2}\rangle$, we obtain

	$\displaystyle Ec_{e}=$	$\displaystyle(\Delta+g)c_{e}+\sqrt{2}\sum_{\rho}\lambda_{\rho}c_{\rho},$		(16)
	$\displaystyle Ec_{\mu}=$	$\displaystyle(\omega_{\mu}-g)c_{\mu}+\sqrt{2}\lambda_{\mu}c_{e}+\sqrt{2}\sum_{\rho}\lambda_{\rho}(d_{\mu\rho}+d_{\rho\mu}),$		(17)
	$\displaystyle(E-\omega_{\rho}$	$\displaystyle-\omega_{\sigma}+\Delta-g)(d_{\sigma\rho}+d_{\rho\sigma})=\sqrt{2}(\lambda_{\rho}c_{\sigma}+\lambda_{\sigma}c_{\rho}).$		(18)

Combining, Eqs. (17) and (18), we have

	$\displaystyle\left(E-\omega_{\mu}+g-2\sum_{\rho}\frac{\lambda_{\rho}^{2}}{E-\omega_{\rho}-\omega_{\mu}+\Delta-g}\right)c_{\mu}$
	$\displaystyle=\sqrt{2}\lambda_{\mu}c_{e}+2\lambda_{\mu}\sum_{\rho}\frac{\lambda_{\rho}c_{\rho}}{E-\omega_{\rho}-\omega_{\mu}+\Delta-g},$		(19)

which is essentially an “integral” equation. Since our numerical results show that at the transition $\alpha\ll 1$, we neglect the sum on the right hand side of Eq. (19), which is equivalent to an approximation of choosing a certain class of diagrams for the self energy. Now combining the approximation with Eq. (16), we have

	$\displaystyle E-\Delta-g=$	$\displaystyle 2\sum_{\mu}\frac{\lambda_{\mu}^{2}}{E-\omega_{\mu}+g-2\sum_{\rho}\frac{\lambda_{\rho}^{2}}{E-\omega_{\rho}-\omega_{\mu}+\Delta-g}}$
	$\displaystyle=$	$\displaystyle 2\int\frac{d\omega}{\pi}\frac{J(\omega)}{E-\omega+g-4\alpha\Delta_{e}(E-\omega+\Delta-g)}.$		(20)

The critical point is when $E=-\Delta+g$ (the band bottom set by the states $\{|gg\rangle\otimes a_{\mu}^{\dagger}a_{\nu}^{\dagger}|0\rangle\}$), which yields

$\displaystyle\Delta=\int\frac{d\omega}{\pi}\frac{J(\omega)}{\Delta+\omega-2g+4\alpha\Delta_{e}(-\omega)},$

(21)

determining the critical value of $\alpha$ for $N=2$. From Eq. (21), one can obtain the critical value $\alpha_{c}^{(2,0)}\approx 0.0247$, a little smaller than $\alpha_{c}=0.025$ for $N=1$, and can show that $\alpha_{c}^{(2,g)}$ decreases with positive $g$ and increases with negative $g$ (see Fig. 7). Qualitatively, $\alpha_{c}^{(2,g)}$ calculated from Eq. (21) differs from the numerical result presented in the main text.