Analytical Study of the Non-Hermitian Semiclassical Rabi Model

The Rabi model, describing the interaction between a two-level system and a classical driving field, is a cornerstone of quantum optics and atomic physics  [1, 2, 3, 4, 5]. In recent years, there has been growing interest in extending the Rabi model to non-Hermitian systems, where the Hamiltonian is no longer self-adjoint  [6, 7, 8, 9, 10, 11, 12]. The non-Hermitian systems, usually accompanied with complex eigenvalues, can exhibit rich and novel phenomena such as non-Hermitian phase transition  [13, 14, 15], non-Hermitian quantum Hall effect  [16, 17], non-reciprocal dynamics [18, 19], non-Hermitian quantum decoherence  [20, 21], and so on. The non-Hermitian Hamiltonian provides a framework for investigating open systems and facilitates a deep understanding of the phenomena associated with non-Hermitian physics  [22].

In particular, one of the most intriguing features of non-Hermitian systems is the concept of parity-time ($\mathcal{PT}$) symmetry, where the combined operation of parity (spatial inversion) and time-reversal symmetry allows for real eigenvalues even in the absence of Hermiticity  [23, 24, 25, 26, 27]. When the system parameters cross certain critical values, the $\mathcal{PT}$ symmetry is broken, leading to the emergence of complex eigenvalues and the formation of exceptional points (EPs), where the eigenvalues and eigenstates coalesce  [28]. The $\mathcal{PT}$-symmetric non-Hermitian systems display numerous novel phenomena, such as $\mathcal{PT}$-symmetric phase transition  [29, 25], sensitivity enhancement  [30, 31, 32], topological insulators  [33, 34, 35, 36], and it thus offers some new physics which is absent in the Hermitian systems  [37, 38, 39].

Although the time-independent non-Hermitian systems has been well developed, the question remains challenging when extended to time-dependent situations  [40]. Several recent works have focused on time evolution and have uncovered a big amount of intriguing dynamical behaviors  [41, 7, 8, 42, 43, 9, 44, 45, 11]. Some efforts have been made to investigate the nature of the non-Hermitian Rabi model (NHRM)  [7, 8, 9, 46]. The NHRM can also be regarded as a periodically driven non-Hermitian two-level system  [24]. The time-independent non-Hermitian two-level system  [24] has been extensively studied in the literature and can be mapped to many important non-Hermitian systems. Hoverer, an exact analytical solution to the NHRM remains elusive, and the physical picture remains unclear. In this study, we aim to find a solvable analytical effective Hamiltonian and gain deeper understanding of the $\mathcal{PT}$-symmetric NHRM.

The paper is structured as follows: In Sec. II, by using a similarity transformation, we transform the NHRM into the rotating-wave approximation (RWA) form with renormalized parameters, making the NHRM effectively solvable. In Sec. III, the analytical $\mathcal{PT}$ spectrum is obtained analytically from the effective solvable model. A comparison with the numerically exact results is performed. We introduce the Floquet parity operator and explain how the $\mathcal{PT}$-broken phase emerges in Sec. IV. Section V is devoted to an analytical analysis of dynamics. The Fourier spectrum analysis is also presented in this section. The analytical Bloch-Siegert shift is presented and analyzed in Sec. VI. The last section contains some concluding remarks.

The Hamiltonian of the NHRM with the purely imaginary coupling strength  [7] can be expressed as

where $\Delta$ represents the atomic frequency, $A$ denotes the coupling strength between the atom and the external field, and $\omega$ stands for the field frequency. The operators $\hat{\sigma}_{k}$ represent the components of the Pauli matrix, where $k=x,y,z$. Notably, this Hamiltonian exhibits $\mathcal{PT}$ symmetry, which can be demonstrated as

where the parity operator $\mathcal{P}=\hat{\sigma}_{z}$ acts as $\hat{\sigma}_{x}\rightarrow-\hat{\sigma}_{x}$, and the time-reversal operator $\mathcal{T}$ transforms $i\rightarrow-i$ and $t\rightarrow-t$. In the $\mathcal{PT}$-unbroken phase, all quasi-energies are real, whereas in the $\mathcal{PT}$-broken phase, the quasi-energies become complex.

In recent experiments, optical and radio-frequency fields provide an effective tool for manipulating atom loss, enabling the implementation of the NHRM Hamiltonian through time-periodic coupling between two spins in a dissipative system of ultracold atoms [46]. The two lowest ${}^{2}S_{1/2}$ hyperfine levels in a non-interacting Fermi gas of ⁶Li atoms are prepared. A resonant optical beam excites the atoms from the lower level to the $2P$ excited state, generating atom loss in the lower level. This state-dependent, purely lossy system allows the lossy Hamiltonian to be mapped onto the $\mathcal{PT}$-symmetric NHRM (1).

The structure of the phase diagram of the NHRM has been theoretically studied using perturbation theory  [3], as discussed in Refs.  [6, 7]. For the Hermitian semi-classical Rabi model, an analytically exact solution remains unavailable. However, some analytical studies based on perturbation theory  [47] and unitary transformations  [48, 49, 50] have been proposed, in addition to numerical solutions based on Floquet theory. For the NHRM, developing an accurate analytical scheme to bridge the gap between perturbation theory and the numerical exact solution would be valuable.

To facilitate the analytical study, we transform the above Hamiltonian in Eq. (1) into an effective RWA. Similar to the unitary transformation in the Hermitian semi-classical RM  [48, 51], we introduce a similarity transformation with a generator in the present NHRM,

where $\alpha$ is a coefficient to be determined later. We can then obtain the following transformed Hamiltonian

$\hat{H}^{\prime}(t)$ can be divided into three parts $\hat{H}^{\prime}(t)=H_{0}^{\prime}+H_{1}^{\prime}+H_{2}^{\prime}\$

where $J_{n}(x)$ is the $n$th-order Bessel function of the first kind, and the following identities have been used

Note that $H_{0}$ is the time-independent component, while $H_{1}$ oscillates with a single-harmonic frequency. The term $H_{2}^{\prime}$ includes all higher-order harmonic terms and Bessel functions with orders $n\geq 2$. By neglecting $H_{2}^{\prime}$, we can derive the effective Hamiltonian that retains only a single-harmonic frequency, which is $H^{\prime}\simeq H_{0}^{\prime}+H_{1}^{\prime}$.

At this point, the coefficient $\alpha$ can be determined by solving the equation,

where $I_{n}(x)$ is $n$th-order modified Bessel function, and the identity $J_{n}(ix)=i^{n}I_{n}(x)$ has been used. The transformed Hamiltonian now takes an RWA-type form with a non-Hermitian coupling strength,

Applying a unitary transformation $\hat{R}(t)=\exp\left(i\omega t\hat{\sigma}_{z}/2\right)$, we obtain a time-independent Hamiltonian as

where

Here, the renormalized atom characteristic frequency $\tilde{\Delta}$ can be interpreted as an effective detuning, while the effective coupling strength $\tilde{A}$ is given by by Eq.  (9). Thus, by performing a single similarity transformation and neglecting higher-order harmonic terms, we derive an effective non-Hermitian two-level atomic Hamiltonian, which forms the foundation of our subsequent work.

An essential characteristic of $\mathcal{PT}$-symmetric models is their spectrum. Despite the non-Hermitian nature of the Hamiltonian, all eigenvalues become real under certain parameter conditions. We designate the system as being in a $\mathcal{PT}$-unbroken phase when these parameters guarantee real eigenvalues. Conversely, when the parameters result in the appearance of imaginary parts in the eigenvalues, we characterize the system as being in a $\mathcal{PT}$-broken phase.

In the time-periodic model, the eigenvalues of the instantaneous Hamiltonian do not have the same physical meaning as those in time-independent scenarios. Therefore, we must study the quasi-eigenvalues. Fortunately, we derived a solvable effective time-independent Hamiltonian in Eq. (10) in the previous section, allowing us to readily obtain its eigenvalues.

The exceptional points are given by $\tilde{\Delta}^{2}=\tilde{A}^{2}/4$, as derived from Eq. (11), which directly yield the $\mathcal{PT}$ phase diagram directly in Fig. 1. When $\tilde{\Delta}^{2}>\tilde{A}^{2}/4$, the model is in a $\mathcal{PT}$-unbroken phase, characterized by real eigenvalues. Conversely, when $\tilde{\Delta}^{2}<\tilde{A}^{2}/4$, the eigenvalues become complex and conjugate in pairs. The black dashed line in Fig. 1 separates the two phases, which are just the EPs.

To solve the time-periodic NHRM numerically, we apply Floquet theory to obtain the quasi-eigenvalues. For the NHRM Hamiltonian (1), the Floquet Hamiltonian is given by $\hat{\mathcal{F}}=\hat{H}(t)-i\frac{\partial}{\partial t}$. As detailed in Appendix A, we perform exact diagonalization (ED) on the matrix form of the Floquet Hamiltonian. Using the numerically obtained quasi-eigenvalues, we determine the boundary between the $\mathcal{PT}$-broken and $\mathcal{PT}$-unbroken phases, as shown in Fig. 1, marked by the blue-shaded region.

Interestingly, the present analytical phase boundary (black dashed lines) agrees excellently with with the exact numerical results over a wide atomic frequency range, up to $\Delta/\omega\approx 3$ (the three-photon resonance), which constitutes the primary region for the $\mathcal{PT}$-broken phase. When $\Delta/\omega>3$, additional $\mathcal{PT}$-broken phases intervene the $\mathcal{PT}$-unbroken region. These secondary $\mathcal{PT}$-broken phases, however, are not captured by the current analytical scheme. Is the present analytical scheme unable to provide any information about the secondary $\mathcal{PT}$-broken phases? The answer is ”no”.

We present the eigenvalues from both the analytical scheme ( Eq. (11), open circles) and the Floquet theory (thick curves) in Fig. 2 for $\Delta/\omega=2.5$ and $3.5$. These two values are marked by two vertical red dashed lines in Fig. 1, both cross the primary $\mathcal{PT}$-broken phase boundary, but $\Delta/\omega=3.5$ line also crosses the secondary $\mathcal{PT}$-broken phase boundary. According to the Floquet theorem, for a time-periodic Hamiltonian, the quasi-eigenvalues shifted by $n\omega$ ($n$ is an integer) do not change the physical state, similar to a generalized Brillouin zone in the Bloch band theory. For later reference, we also display some shifted quasi-eigenvalues for both the analytical scheme (Eq.  (11)) and Floquet theory.

Figure 2 (a) shows that the analytical results closely agree with the numerical quasi-energies for the coupling strengths up to $A\simeq 6\omega$ for $\Delta/\omega=2.5$. The analytical EPs, indicated by dashed lines, are nearly identical to the numerical ones.

As shown in Fig. 2 (b) for$\ \Delta/\omega=3.5$, the analytical quasi-eigenvalues agree with the numerical results, except in the small, ellipse-shaped region around quasi-eigenvalues $\varepsilon/\omega=(2n+1)/2$. In this region, the analytical scheme can only provide real quasi-eigenvalues, and level crossings occur, while the numerical results yield complex quasi-eigenvalues. Further analysis will be presented in the next section to explain these observations, based on the uncovered symmetry of the NHRM within the framework of Floquet theory.

To explain the $\mathcal{PT}$-symmetry breaking and the EPs discussed in the previous section, we introduce the Floquet parity operator

where $\hat{G}=\sum_{n=-\infty}^{\infty}n|n\rangle\langle n|$ is constructed using the state vector $|n\rangle$ from an orthonormal, complete Fourier basis set $\exp\left(-in\omega t\right)$ in the Hilbert space of $2\pi/\omega$-periodic function. It can be shown that $\hat{\Pi}\mathcal{F}\hat{\Pi}^{\dagger}=\mathcal{F}$, where $\mathcal{F}=\hat{H}-i\partial/\partial t$ is the Floquet Hamiltonian. This implies that the Hilbert space of the Floquet Hamiltonian can be decomposed into two independent subspaces, corresponding to even and odd Floquet parity, with eigenvalues $\Pi=\pm 1$, respectively.

For odd (even) parity ($\Pi=\mp 1$), the subspace basis should be $\{\cdots|\mp,-1\rangle,|\pm,0\rangle,|\mp,1\rangle,\cdots\}$ , so that the Floquet Hamiltonian can be expressed as

We can also numerically diagonalize this matrix for both even and odd parity to calculate the corresponding quasi-eigenvalues.

The quasi-eigenvalues, denoted by the thick curves in Fig. 2, correspond to the odd parity subspace. Quasi-eigenvalues shifted by $2n\omega$, marked in black, also belong to this subspace. In contrast, the quasi-eigenvalues of even parity are shifted by $(2n+1)\omega$, as marked in red. In principle, quasi-eigenvalues with different parities may cross. However, within the same parity subspace, they cannot intersect without an additional symmetry. Fortunately the present NHRM just possesses an alternative $\mathcal{PT}$ symmetry, which enables quasi-eigenvalues within the same parity subspace to intersect, leading to EPs. This phenomenon underlies the $\mathcal{PT}$ symmetry breaking and the emergence of EPs.

The picture above illustrating the $\mathcal{PT}$ symmetry breaking fully describes the features of the quasi-eigenvalue spectrum shown in Fig. 2. Two branches of the thick curve (representing real values), both with the same odd parity, approach each other as $A/\omega$ increases and intersect at the vertical dashed line, forming the primary EPs. These EPs can be accurately captured by $\tilde{\Delta}^{2}=\tilde{A}^{2}/4$ in Eq. (11).

As mentioned above, the eigenvalue $\varepsilon_{-}$, shifted by $2n\omega$, retains the same parity. If the shifted $\varepsilon_{-}$ cross with $\varepsilon_{+}$, i.e., $\varepsilon_{+}=\varepsilon_{-}+2n\omega$, then we have

This equation gives the coupling strength $A/\omega$ for a given $\Delta/\omega$, at which the real eigenvalues cross, thus corresponding to the $\mathcal{PT}$-broken region. The primary EP corresponds to $n=0$ in this equation. We can also determine the crossing points where real eigenvalues with the same parity intersect, for any $n>0$.

The lines representing the crossing points for $n=1$ and $n=2$ are shown as red solid lines in Fig. 1, both lying within the $\mathcal{PT}$-broken region. Using perturbation theory, Lee et al. predicted a line of maximal $\mathcal{PT}$ symmetry breaking (where the imaginary eigenvalue is maximized) in Ref. [7], shown as a blue dash-dotted line in Fig.1 for comparison. The two lines agree well for small $A/\omega$ but begin to deviate as $A/\omega$ increases. Noticeably, their line for $n=2$ does not lie within the $\mathcal{PT}$-broken region, while ours does. Therefore, our analytical scheme can more accurately predict the onset of $\mathcal{PT}$ symmetry breaking, beyond the primary $\mathcal{PT}$-broken phase, particularly when $\Delta$ exceeds $5\omega$.

If $A$ is very small, and higher-order terms in $A$ are neglected, we can derive a closed-form solution to Eq. (13),

Based on the similarity transformation in Sec. II, this result is exact in the limit as $A/\omega\to 0$. Since $\Delta=(2n+1)\omega$, the crossing point occurs at $A/\omega=0$. This indicates that the $\mathcal{PT}$-broken line intersect the $A/\omega=0$ axis at these resonance points.

As shown in Fig. 1 and in Ref. [7], the second and third $\mathcal{PT}$-broken phases, for $n=1$ and $2$, are expected to be extremely narrow at small values of $A/\omega$ near 0, but they do not extend to the $A/\omega=0$ axis according to numerical calculations. This discrepancy arises from the numerical calculations, which require truncation, while the exact solution can only be obtained using an infinite basis. At the single-photon resonance $\Delta=\omega$ for $n=0$, Eq. (14) gives $A/\omega=0$, consistent with the primary $\mathcal{PT}$-broken region, which touches $A/\omega=0$ axis at $\Delta/\omega=1$, as shown in Fig. 1.

We now focus on the more complicated Floquet quasi-eigenvalues spectrum in Fig. 2 (b) for $\Delta=3.5\omega$. Specifically, within the small, ellipse-shaped region around quasi-eigenvalues $\varepsilon=(2n+1)\omega/2$, the present analytical results from Eq. (11), denoted by the open circles, show real level crossings. The crossing point at $A/\omega=2.39$ can be calculated from Eq.  (13) for $n=1$. This is also the point where the first red solid line intersects the second vertical red dashed line in Fig.1. However, the two real quasi-eigenvalues with the same parity cannot cross. As a result, an imaginary part emerges to resolve this issue, signaling the new PT symmetry breaking. The analytical results in Eq. (11) cannot fully describe, but can predict this $\mathcal{PT}$-broken region. In contrast, real spectrum lines with different parities can cross at $\varepsilon=n\omega$.

When $\Delta/\omega=(2n+1)$, i.e. $(2n+1)$ -photon resonance, the two quasi-eigenvalues with the same parity coincide at $A/\omega=0$, as proven in Eq.  (14) for any $n$. If $\Delta$ increases slightly, the two quasi-eigenvalues will intersect. However, they cannot cross because they share the same parity. As a result, a complex quasi-eigenvalue must emerge to prevent the crossing of real quasi-eigenvalues, leading to the formation of EP points.

If $3<\Delta/\omega<5$, as demonstrated in Fig. 2 (b), the two adjacent quasi-eigenvalues obtained from the numerical diagonalization of the matrix (12), which share the same parity (both black or both red lines), should avoid crossing each other at a finite value of $A/\omega$. Consequently, a $\mathcal{PT}$-broken phase emerges, indicated by small, ellipse-shaped blue curves representing the imaginary part of the quasi-eigenvalues. This explains the origin of the second $\mathcal{PT}$-broken region in Fig. 1. However, two adjacent real quasi-eigenvalues from the analytical scheme cross in the $\mathcal{PT}$-broken region. In other words, the analytical scheme cannot fully describe this $\mathcal{PT}$-broken region, but can predict its location. If we further increase $\Delta/\omega$ to the range $5<\Delta/\omega<7$, both adjacent and next-adjacent numerical quasi-eigenvalues avoid the real quasi-eigenvalue crossing. Consequently, two $\mathcal{PT}$-broken phases emerge, as illustrated in Fig. 1 by the second and third $\mathcal{PT}$-broken regions. At large $A/\omega$, the widths of the different $\mathcal{PT}$-broken regions can become large enough to merge into a single, contiguous $\mathcal{PT}$-broken region.

So far, we have explained the origin of the $\mathcal{PT}$-broken phases and accounted for the phase diagram using the conserved Floquet parity operator. Additionally, the Floquet parity operator can be applied to the Hermitian semi-classical Rabi model. Within the same parity subspace, real quasi-energies cannot cross but instead exhibit an avoided-crossing feature  [3]. As shown in Appendix A, quasi-energies with different parities can cross, whereas quasi-energies with the same parity can only avoid each other. The present non-Hermitian model exhibits richer and more complex properties than its Hermitian counterpart, due to the additional $\mathcal{PT}$ symmetry.

Our current analytical scheme, based on a single similarity transformation, accurately describes the phase diagram for $\Delta/\omega<3$ and can also predict the new $\mathcal{PT}$-broken region at $\Delta/\omega>3$. By combining the primary $\mathcal{PT}$-broken region for $n=0$, determined by Eq. (11), with the $\mathcal{PT}$-broken line from Eq. (13) for $n>0$, as shown in Fig. 1, our analytical scheme can produce a highly accurate, nearly exact phase diagram within $A/\omega<2$ for arbitrary $\Delta/\omega$, except for the small region of the second $\mathcal{PT}$-broken phase. In the next section, we will examine whether it accurately captures the dynamics.

In this section, We explore the time-evolution of the two-level system. The time dependent wavefunction can be represented as $|\psi(t)\rangle=(c_{+}(t)\ c_{-}(t))^{\text{T}}$, where $c_{+}(t)$ and $c_{-}(t)$ denote the probability amplitudes of the upper and lower states, respectively, in a two-level atom. Using the Schrödinger equation, we can derive the following two coupled first-order differential equations:

The exact dynamics can be determined by numerically solving two coupled dynamical equations, as an analytical solution is highly challenging. An effective approximation involves applying the RWA to decouple the probability amplitudes. This approach simplifies the calculation and provides the population evolution of the excited state as

where $\Omega_{\mathrm{RWA}}=\sqrt{(\Delta-\omega)^{2}-(A/2)^{2}}$. The population grows exponentially when $\Omega_{\mathrm{RWA}}$ is imaginary and can exceed unity due to the non-conservation of probability in non-Hermitian systems.

The above quasi-eigenvalues in Eq.  (11) allow us to analytically study the dynamical properties of the NHRM. Using the similarity transformation (2) and the unitary transformation $\hat{R}(t)$, we can define the time evolution operator as

by which we can derive analytically the time-dependent wavefunction $|\psi(t)\rangle$ starting from an initial state $|\psi_{i}\rangle$. The time evolution operator in Eq.  (18) applies a sequence of transformations to the initial state, moving it to a new gauge. The state evolves under the time-independent Hamiltonian $\tilde{H}$ until time $t$. Finally, we arrive at the evolved state to the original frame using the transformation $e^{-\hat{S}(t)}\hat{R}^{\dagger}(t)$, yielding the final state $|\psi(t)\rangle=\hat{U}(t)|\psi_{i}\rangle$.

The atom is initially prepared in the ground state, so $|\psi_{i}\rangle=(0\ 1)^{\text{T}}$. In our analysis, we focus on the atomic population of the excited state, denoted as $P_{e}(t)=|c_{+}(t)|^{2}$, representing the probability of finding the atom in its excited state at time $t$. This initial state has a definite parity, so its evolution involves only eigenstates in the same parity subspace. The analytical result can be written as:

where

Here, $\Omega_{R}=\left(\varepsilon_{+}-\varepsilon_{-}\right)/2$ is the Rabi frequency from Eq. (11). $2\Omega_{R}$ measures the gap between the two branches of the quasi-eigenvalues for the same $n$ in the eigenvalue spectrum, as illustrated in Fig. 2 (b) by the black double arrow connecting the open circles. Adjacent quasi-eigenvalues with the same parity are shifted by $2\omega$. The quasi-eigenvalue differences dominating the dynamics can include $2\omega$, $2\omega\pm\Omega_{R}$, as well as other combinations and multiples of $2\omega$ and $2\Omega_{R}$.

In Fig. 3, we show the time evolution of the atom at three representative points: A ($\Delta/\omega=2.5$, $A/\omega=1$), B ($\Delta/\omega=3.5$, $A/\omega=3$), and C ($\Delta/\omega=2.5$, $A/\omega=4$). These points lie in different regions of the phase diagram in Fig. 1, and we analyze them using the present analytical approach (19), numerical simulations, and the RWA (17). As illustrated in Fig. 3 (c), at point C in the $\mathcal{PT}$-broken phase, dominated by single-photon resonance, the population in the excited state grows roughly exponentially, eventually surpassing unity. This feature is fully described by Eq.  (19). This intriguing emergent phenomenon has significant implications for non-Hermitian dynamics. Unlike Hermitian systems, our model equivalently treats open systems, leading to a breakdown in probability conservation. Interestingly, the present analytical results for population evolution perfectly agree with the numerical results. The RWA can only accurately describe population evolution in the short-time regime. While the RWA can predict the dominant exponential behavior (cf. Eq. (17)), it underestimates the growth rate. Additionally, we find that the RWA yields a smoothed evolution without the oscillatory feature due to the neglect of counter-rotating terms.

In the $\mathcal{PT}$-unbroken phase, represented in Fig. 3 (a) and (b), stable oscillations with multiple frequencies and varying amplitudes are observed, in contrast to the $\mathcal{PT}$-broken phase, where exponential dynamics dominate. The present analytical results, given by Eq. (19), perfectly describe these oscillations at point A for $\Delta/\omega<3$. At point B, the analytical results match the exact dynamics only at the very early stage and deviate significantly throughout the rest of the process. Nevertheless, they qualitatively capture the multi-oscillation feature, including oscillations with multiple frequencies.

To analyze the dynamics and stable oscillations in the $\mathcal{PT}$-unbroken phase in greater detail, we perform a Fourier spectral analysis of the population time evolution. The Fourier transform, a fundamental technique for frequency decomposition, is defined as follows:

where $\mathcal{F}(\nu)$ represents the spectrum as a function of frequency $\nu$, and $F(t)$ denotes the observed time evolution as a function of time $t$.

The frequency spectrum of the atomic population time evolution is shown in Fig. 4. At the point A, as shown in Fig. 4 (a), the numerical frequency spectrum exhibits three peaks, precisely corresponding to $2\Omega_{R}$, $2\omega$, and $2\omega+2\Omega_{R}$, all these three frequencies measure the gaps between the two quasi-eigenvalues within the same parity in the analytical scheme. These frequencies also correspond to the lowest three frequencies in Eq. (19), further demonstrating that our analytical method accurately describes the time evolution in the $\mathcal{PT}$-unbroken phase for $\Delta/\omega<3$.

As shown in Fig. 4 (b) for point B, the exact frequency spectrum exhibits five peaks, while our analytical results yield four main peaks, missing only the last one. The analytical frequency peaks are marked with arrows, indicating the differences between two quasi-eigenvalues of the same parity. In both Fig. 4 (a) and (b), the main peaks corresponding to $2\Omega_{R}$, $2\omega$, and $2\omega+2\Omega_{R}$ also appear in the exact dynamics. In panel (b), however, the additional peaks around $2\omega-2\Omega_{R}$ and $4\omega$ in the numerically exact dynamics are only weakly captured in the analytical results. At point B, the atomic frequency $\Delta/\omega=3.5$ lies beyond the three-photon resonance region, so additional higher-order frequency terms are required. Therefore, to obtain a more precise description, additional similarity transformations should be added in Eq.  (2). For qualitative studies, the current analytical scheme, which uses a single similarity transformation, can capture the main features of the dynamics.

We have also extensively examine the dynamics for large atomic frequencies $\Delta/\omega\geq 4$, in the $\mathcal{PT}$-unbroken phase, as demonstrated in Appendix B. Our analytical results agree well with the exact ones as long as the coupling strength is not too strong. The RWA can only provide the primary oscillation frequency, known as the Rabi frequency, and thus cannot fully capture the complex features of the time evolution, as shown in Figs. 3,  4, and 7. The RWA Rabi frequency also deviates from the exact values at strong coupling, whereas the analytical Rabi frequency agrees well with the exact values in all cases.

In this section, we will discuss the Bloch-Siegert (BS) shift in the NHRM, which is a key issue in its Hermitian counterpart  [3, 49, 52].

In the Hermitian system, the BS shift represents the offset of the resonance frequency $\Delta_{\mathrm{res}}$ from the original driven frequency, i.e., $\delta\omega_{\mathrm{BS}}=\Delta_{\mathrm{res}}-\omega$. The resonance frequency is determined by the maximum time-averaged transition probability, corresponding to the largest oscillation amplitude. However, in the non-Hermitian system, when the frequency reaches resonance, the probability diverges rapidly. In this case, the quasienergies exhibit the largest imaginary component. The condition of the resonance is

where $\partial\alpha/\partial\omega$ can be obtained from Eq. (9) as

Solving Eq.  (20) yields $\Delta_{\mathrm{res}}$, allowing us to compute the BS shift in the present analytical scheme.

A simple analytical expression for the BS shift, expanded in powers of the coupling strength $A/\omega$, can also be derived. In this process, we expand $I_{0}\left(\frac{A}{\omega}\alpha\right)$ and $I_{1}\left(\frac{A}{\omega}\alpha\right)$ up to the fourth order in $A$,

Using Eq. (9), we expand the coefficient $\alpha$ up to second order in $A$,

The minus sign in the second term arises from the imaginary coupling $iA$, which changes the sign in the Hermitian counterpart  [48]. Therefore, the modulated Rabi frequency, up to fourth order in $A$, is given by

We can now obtain the resonance frequency,

Note that the $A^{2}$ term is consistent with the one found in perturbation theory [7]. Additionally, by applying the concise similarity transformation method, we can derive the higher-order solution for the BS shift.

The analytical BS shift, based on Eq.  (20), is shown by the red line in Fig. 5 . Interestingly, our results match the numerical data very well across a wide range of coupling strengths, up to $A/\omega=5$. Note that the perturbation result for the BS in Ref.  [7], marked in blue, begins to deviate from the exact values at $A/\omega=1.5$. Surprisingly, the presence of the $\mathcal{PT}$-broken phase in this non-Hermitian model does not significantly affect the BS shift across a broad range of coupling strengths, even when passing through multi-photon resonance. This may be due to the fact that the BS shift can be expressed as a function of $(iA)^{2}$.

This work studies the NHRM analytically. First, we derive a transformed NHRM in the RWA form using a single similarity transformation. Next, we obtain an effective Hamiltonian for the non-Hermitian two-level system by applying a rotating operation, which explicitly yields two eigenvalues and an EP. Thus, the quasi-eigenvalues can thus be obtained from the two eigenvalues by shifting them by a multiple of the field frequency. The exact quasi-eigenvalues can be determined numerically using Floquet theory. Notably, our analytical results for the quasi-eigenvalues agree excellently with the numerical results in the single-photon resonance regime. The primary $\mathcal{PT}$-broken phase boundary is also detected in a wide parameter regime, extending even beyond the multi-photon resonance region.

We also uncover a Floquet parity symmetry that commutes with the Floquet Hamiltonian of the NHRM. The conserved parity decomposes the entire Floquet space into two subspaces with different parity. The quasi-eigenvalues are associated with one of two parities. Quasi-eigenvalues shifted by even multiples of the field frequency, $2n\omega$, belong to the same parity. When the atomic frequency exceeds an odd multiple of the field frequency, i.e., $\Delta>(2n+1)\omega$, two adjacent quasi-eigenvalues within the same parity will intersect. However, they cannot truly cross without an additional symmetry. Therefore, a complex quasi-eigenvalue emerges to resolve this issue, leading to the secondary $\mathcal{PT}$-broken phase. Interestingly, in this region, our analytical real quasi-eigenvalues cross, predicting the $\mathcal{PT}$-broken location, though it cannot describe the $\mathcal{PT}$-broken.

The dynamics are also investigated within the analytical scheme. In the single-photon resonance regime, i.e., $\Delta/\omega<3$, the time evolution of the atomic population, starting from the atomic upper state, as predicted by the analytical scheme, closely matches the numerical results. Even in the higher-photon resonance regime, Fourier transform analysis shows that the analytical scheme can captures multi-oscillations, with the dominant frequencies aligning with those observed in the exact dynamics. Additionally, the BS shift, as derived analytically in the present scheme, agrees well with the exact result over a wide coupling regime. This further highlights that the simple analytical scheme effectively captures the main features of this non-Hermitian qubit-cavity coupling model.

To accurately describe the secondary $\mathcal{PT}$-broken phase and dynamics in higher-order resonance regimes, additional similarity transformations are required. The research along this direction is in progressing. The proposed analytical method for the NHRM may offers a valuable alternative for exploring open systems and atom-field interactions. Its versatility extends to more complex models, including those with multiple cavity frequencies, higher-order interactions, or additional qubits.

In this Appendix, we provide a brief description of the Floquet picture for the NHRM. A complete set of orthonormal basis states $|\alpha,m(t)\rangle\rangle$ of space $\mathscr{F}$ can be constructed by combining a complete set of orthogonal and normalized basis states of $\hat{H}$, $|\alpha\rangle=|\pm\rangle$, with the complete set of time-periodic functions $e^{im\omega t}$ labeled by the integer $m$,

Using the definition of the scalar product in the extended Floquet Hilbert space, the matrix elements can be obtained $\mathcal{F}_{m^{\prime}m}^{\alpha^{\prime}\alpha}=\langle\langle\alpha^{\prime},m^{\prime}|\hat{\mathcal{F}}|\alpha,m\rangle\rangle$ of the Floquet Hamiltonian $\hat{\mathcal{F}}=\hat{H}(t)-i\frac{\partial}{\partial t}$ with respect to the basis $|\alpha,m\rangle\rangle$,

where

In the matrix form, the Floquet Hamiltonian can be expressed as

The quasi-eigenvalues are the eigenvalues of this infinite matrix. The exact numerical diagonalization of the matrix (24) is performed using a large truncated space.

For comparison, we calculate the Floquet quasi-energies for the Hermitian counterpart using the same model parameters as in Fig. 2. The results are presented in Fig. 6. In the time-periodic Hermitian model, quasi-energies with the same parity exhibit only the avoided-crossing feature, in contrast to its non-Hermitian counterpart. Quasi-energies with different parities can cross. This is consistent with the Floquet parity symmetry.

We can examine the dynamics using our analytic approach in $\mathcal{PT}$-unbroken phase over a broader parameter regime. As shown in Fig. 7, we compare the dynamics for the different atomic frequencies: (a) $\Delta=4\omega$, (b) $\Delta=5\omega$, and (c) $\Delta=6\omega$, all at the same coupling strength $A=2\omega$. We observe that the population of the excited state maintains stable oscillations, and our analytical method matches remarkably well with the numerical results.

As shown in Fig. 7 (c), even when the atomic frequency lies above the higher-order $\mathcal{PT}$-broken regions, such as the five-photon resonance region, our analytical approach still accurately captures the dynamic behavior. This suggests that our analytic method is available for a wide range of coupling strengths and arbitrary atomic frequencies. Interestingly, as the frequency of atoms increases, the analytical results align more closely with the numerical ones.