# Does Medical Imaging learn different Convolution Filters?

Paul Gavrikov1 Janis Keuper1,2
1IMLA, Offenburg University, Germany
2Fraunhofer ITWM, Kaiserslautern, Germany
{paul.gavrikov<EMAIL_ADDRESS>

###### Abstract

Recent work has investigated the distributions of learned convolution filters
through a large-scale study containing hundreds of heterogeneous image models.
Surprisingly, on average, the distributions only show minor drifts in
comparisons of various studied dimensions including the learned task, image
domain, or dataset. However, among the studied image domains, medical imaging
models appeared to show significant outliers through “spikey” distributions,
and, therefore, learn clusters of highly specific filters different from other
domains. Following this observation, we study the collected medical imaging
models in more detail. We show that instead of fundamental differences, the
outliers are due to specific processing in some architectures. Quite the
contrary, for standardized architectures, we find that models trained on
medical data do not significantly differ in their filter distributions from
similar architectures trained on data from other domains. Our conclusions
reinforce previous hypotheses stating that pre-training of imaging models can
be done with any kind of diverse image data.

##  1 Introduction

Deep learning has accelerated the progress in many computer vision problems,
including various applications in medical imaging [1, 2, 3, 4]. Most deployed
architectures are based on convolutional neural networks (CNNs) which learn
large amounts of convolution filters to transform inputs into meaningful
feature spaces that can be further processed. Since the same CNN architectures
can learn many computer vision tasks and work on diverse image domains, the
question naturally arises whether there are shared similarities in the learned
convolution filters.

Previously, we introduced a large-scale dataset (CNN Filter DB) [5] that
consists of convolution filters with a kernel size of $3\times 3$ extracted
from hundreds of publicly available pre-trained CNNs trained for a multitude
of tasks. On that data, we performed a singular-value decomposition (SVD)
comparison of sets of filters aggregated by various dimensions of meta-data.
Therefore, sets of normalized filters were transformed to a common SVD-
obtained basis, and the shift between resulting distributions of the
coefficients ($c_{ij}$) was measured based on a variant of the Kullback-
Leibler Divergence (KL). Generally, individual models and model families (e.g.
ResNet) learn unique distributions of filters, but on average, distributions
segregated by the learned task, image domain, or dataset do not significantly
differ. Most observed shifts were due to different levels of “degenerated”
(i.e. sparse or highly repetitive) filters caused by the over-parameterization
of networks relative to the dataset. However, the most salient outlier in
shifts included medical imaging models, as this domain showed distinct
“spikey” coefficient distributions. In this paper, we set our focus to study
these models in detail and understand the possible causes of the shifts. The
models can be summarized in the following: CompNet [6] contains three
customized architectures trained for brain segmentation on the OASIS (MRI)
dataset [7]; LungMask [8] contains three UNet trained for lung segmentation on
the LTRC and some proprietary CT datasets; TorchXRayVision [9] contains seven
DensetNet-121 and one ResNet50 classifier trained to detect various
pathologies from various datasets. Additionally, the framework includes one
ResNet101 auto-encoder trained on a large aggregated dataset consisting of the
previously mentioned datasets; Lastly, one individual UNet [10] is included
that is also trained to perform brain segmentation on a public Kaggle MRI
dataset [11].

##  2 Analysis

Figure 1: KDEs of the coefficient distributions along every principal
component for all medical imaging models from CNN Filter DB.

In this section, we apply the same SVD-decomposition as defined in [5] to
visualize the coefficient distributions in kernel density estimation (KDE)
plots. We collect all filters for each medical model and collectively
visualize the resulting coefficients (Figure 1). First, we observe that the
salient “spikey” distributions are only seen in models of CompNet on all axes,
while all other models form less salient distributions. In the original
analysis, this model family was over-represented due to their relatively large
number of convolutions.

The CompNet architecture can be split into multiple encoder-decoder models:
source inputs are fed into an encoder that feeds into two parallel decoder
branches. One of those constructs a binary brain segmentation mask, and the
other computes a complementary output for the non-brain part of the image.
Finally, the outputs are aggregated and further processed by an auto-encoder
that learns a latent-space representation.

Figure 2: The first filters from a layer that learns only binary-like filters
found in CompNet.

Upon closer analysis, we observe that some layers learn filters with binary
weights (Figure 2) which explains the coefficient spikes in the KDEs, while
other layers learn a less conspicuous set of filters. We assume that the
binary filters are learned by the final auto-encoder as indicated by the
presence in later convolution stages, yet, we are unable to prove our
hypothesis due to the highly entangled implementation of the source model.
Still, it seems worthwhile to understand the cause of binary filters as it is
possible to learn them from a significantly lower dimensional and discrete
search space which should accelerate learning and require fewer data samples.
Interestingly, we find fewer such binary filter layers in the model trained on
the axial reconstruction plane.

Next, complementary to the findings in [5], we observe that models based on
the same architecture (e.g. DenseNet121) learn quite similar distributions
when compared to each other but differ substantially when compared to other
architectures: all CompNets learn similar coefficient distributions although
trained on different reconstruction planes; Albeit trained in different
regimes, all LungMask models also form similar distributions to the standalone
UNet. Still, some shifts in the distributions can be seen, in particular in
the DenseNet121 models from TorchXRayVision which were trained on different
datasets. The model trained on the Kaggle RSNA Pneumonia Challenge dataset
seems to overrepresent specific coefficients (and therefore filters clusters)
by a large amount, while the models trained on other datasets such as MIMIC-
CXR [12], NIH [13, 14], and CheXpert [15] tend to learn coefficient
distributions that are smoother and more similar to the distribution of the
model that was trained on all datasets. The most salient coefficient
distribution is seen in the ResNet50 which was trained on $512\times 512$ px
images instead of $256\times 256$ px. Spikes dominate the KDEs at 0 for all
axes, which indicates a highly sparse model. Indeed, by applying a structured
pruning on the kernels, we find that only approx. 1% of filters (N=12900) is
non-sparse. Due to a lack of access to the datasets, we cannot verify the
integrity of this model. However, if this model’s performance is on-par with
the other models, this opens the question of whether a) smaller and therefore
easier-to-train models would not be more suitable for this problem and perhaps
brain segmentation in general and b) higher resolution images may be more
suitable as they may contain more salient features for classification than
their down-sampled counter-parts.

Figure 3: Heatmaps showing the pair-wise distribution shift for UNet (left)
and DenseNet (right) trained on different datasets. Low values/dark colors
denote low shifts.

Finally, we compare the medical UNet segmentation and DenseNet121
classification models to similar models from CNN Filter DB trained on other
image domains (cars, seismic data, and natural images) to understand whether
medical models learn different convolution filters. We compute the KL in a
pair-wise manner between all model combinations and display the results in the
form of heatmaps (Figure 3). We observe that UNet learns strikingly different
representations on seismic data, but only minor differences when trained to
segment cars. DenseNet121 also learns slightly different representation when
trained on ImageNet [16]. Yet, it is worth noting that this is a significantly
larger dataset with 1.4 M samples and 1000 categories that occupy more of the
model capacity. Finally, the most noticeable shift is seen in the model
trained on the Kaggle dataset, although it originates from the same image
domain.

In all cases, it appears as if the shifts between filters of medical domains
and others are insignificant and may even vanish when compared to more models.

##  3 Conclusion and Future Work

The previously observed outliers in distributions [5] are directly caused by
specific branches in architectures, such as [6]. After all, it turns out that
medical imaging models do not learn fundamentally different filter
distributions than models of other image domains. In turn, this means that
pre-training with diverse image data originating from arbitrary domains should
accelerate the training of medical models (and others). The remaining shifts
to other domains may be explained by different levels of over-parameterization
and bear potential for future work which may study the filter quality through
metrics proposed in [5, 17]. For comparison, a larger sample size of models
should be used.

Additionally, we have observed that some layers learn highly specific clusters
of filters. Others, seem to develop extreme levels of sparsity. In those
cases, it would be reasonable to reduce the search space during training, and
in turn, accelerate training, and/or reduce the amount of necessary training
data. Given the scarcity of medical data and its acquisition cost, it seems
worthwhile to explore this direction in future work.

## Potential negative societal impact

Our paper is limited to the continued analysis of learned models, so we do not
believe our findings cause any negative societal impact. However, on contrary,
we hope that they can accelerate the learning of medical imaging models which
would result in reduced CO2 emissions during training and increase AI
democratization due to the reduced amount of necessary training data.

## Acknowledgments and Disclosure of Funding

Funded by the Ministry for Science, Research and Arts, Baden-Wuerttemberg,
Germany, Grant 32-7545.20/45/1 (Q-AMeLiA).

## References

  * [1] H.-C. Shin, H. R. Roth, M. Gao, L. Lu, Z. Xu, I. Nogues, J. Yao, D. Mollura, and R. M. Summers, “Deep convolutional neural networks for computer-aided detection: Cnn architectures, dataset characteristics and transfer learning,” _IEEE Transactions on Medical Imaging_ , vol. 35, no. 5, pp. 1285–1298, 2016.
  * [2] K. Suzuki, “Overview of deep learning in medical imaging,” _Radiol. Phys. Technol._ , vol. 10, no. 3, pp. 257–273, Sep. 2017.
  * [3] O. Ronneberger, P. Fischer, and T. Brox, “U-Net: Convolutional Networks for Biomedical Image Segmentation,” in _Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015_. Cham, Switzerland: Springer, Nov. 2015, pp. 234–241.
  * [4] P. Moeskops, M. A. Viergever, A. M. Mendrik, L. S. de Vries, M. J. N. L. Benders, and I. Išgum, “Automatic segmentation of mr brain images with a convolutional neural network,” _IEEE Transactions on Medical Imaging_ , vol. 35, no. 5, pp. 1252–1261, 2016.
  * [5] P. Gavrikov and J. Keuper, “CNN Filter DB: An Empirical Investigation of Trained Convolutional Filters,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2022, pp. 19 066–19 076.
  * [6] R. Dey and Y. Hong, “CompNet: Complementary Segmentation Network for Brain MRI Extraction,” in _Medical Image Computing and Computer Assisted Intervention – MICCAI 2018_ , A. F. Frangi, J. A. Schnabel, C. Davatzikos, C. Alberola-López, and G. Fichtinger, Eds. Cham: Springer International Publishing, 2018, pp. 628–636.
  * [7] D. S. Marcus, T. H. Wang, J. Parker, J. G. Csernansky, J. C. Morris, and R. L. Buckner, “Open Access Series of Imaging Studies (OASIS): cross-sectional MRI data in young, middle aged, nondemented, and demented older adults,” _Journal of Cognitive Neuroscience_ , vol. 19, no. 9, pp. 1498–1507, Sep. 2007.
  * [8] J. Hofmanninger, F. Prayer, J. Pan, S. Röhrich, H. Prosch, and G. Langs, “Automatic lung segmentation in routine imaging is primarily a data diversity problem, not a methodology problem,” _European Radiology Experimental_ , vol. 4, no. 1, pp. 1–13, 2020.
  * [9] J. P. Cohen, J. D. Viviano, P. Bertin, P. Morrison, P. Torabian, M. Guarrera, M. P. Lungren, A. Chaudhari, R. Brooks, M. Hashir, and H. Bertrand, “TorchXRayVision: A library of chest X-ray datasets and models,” in _Medical Imaging with Deep Learning_ , 2022. [Online]. Available: https://github.com/mlmed/torchxrayvision
  * [10] M. Buda, A. Saha, and M. A. Mazurowski, “Association of genomic subtypes of lower-grade gliomas with shape features automatically extracted by a deep learning algorithm,” _Computers in Biology and Medicine_ , vol. 109, 2019\.
  * [11] G. Shih, C. C. Wu, S. S. Halabi, M. D. Kohli, L. M. Prevedello, T. S. Cook, A. Sharma, J. K. Amorosa, V. Arteaga, M. Galperin-Aizenberg, R. R. Gill, M. C. B. Godoy, S. Hobbs, J. Jeudy, A. Laroia, P. N. Shah, D. Vummidi, K. Yaddanapudi, and A. Stein, “Augmenting the National Institutes of Health Chest Radiograph Dataset with Expert Annotations of Possible Pneumonia,” _Radiology: Artificial Intelligence_ , Jan. 2019. [Online]. Available: https://pubs.rsna.org/doi/full/10.1148/ryai.2019180041
  * [12] A. E. W. Johnson, T. J. Pollard, N. R. Greenbaum, M. P. Lungren, C.-y. Deng, Y. Peng, Z. Lu, R. G. Mark, S. J. Berkowitz, and S. Horng, “Mimic-cxr-jpg, a large publicly available database of labeled chest radiographs,” 2019. [Online]. Available: https://arxiv.org/abs/1901.07042
  * [13] X. Wang, Y. Peng, L. Lu, Z. Lu, M. Bagheri, and R. M. Summers, “Chestx-ray8: Hospital-scale chest x-ray database and benchmarks on weakly-supervised classification and localization of common thorax diseases,” in _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_. Los Alamitos, CA, USA: IEEE Computer Society, jul 2017, pp. 3462–3471. [Online]. Available: https://doi.ieeecomputersociety.org/10.1109/CVPR.2017.369
  * [14] A. Majkowska, S. Mittal, D. F. Steiner, J. J. Reicher, S. M. McKinney, G. E. Duggan, K. Eswaran, P.-H. Cameron Chen, Y. Liu, S. R. Kalidindi, A. Ding, G. S. Corrado, D. Tse, and S. Shetty, “Chest radiograph interpretation with deep learning models: Assessment with radiologist-adjudicated reference standards and population-adjusted evaluation,” _Radiology_ , vol. 294, no. 2, pp. 421–431, 2020, pMID: 31793848. [Online]. Available: https://doi.org/10.1148/radiol.2019191293
  * [15] J. Irvin, P. Rajpurkar, M. Ko, Y. Yu, S. Ciurea-Ilcus, C. Chute, H. Marklund, B. Haghgoo, R. Ball, K. Shpanskaya, J. Seekins, D. A. Mong, S. S. Halabi, J. K. Sandberg, R. Jones, D. B. Larson, C. P. Langlotz, B. N. Patel, M. P. Lungren, and A. Y. Ng, “Chexpert: A large chest radiograph dataset with uncertainty labels and expert comparison,” 2019. [Online]. Available: https://arxiv.org/abs/1901.07031
  * [16] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei, “Imagenet: A large-scale hierarchical image database,” in _2009 IEEE conference on computer vision and pattern recognition_. Ieee, 2009, pp. 248–255.
  * [17] P. Gavrikov and J. Keuper, “Adversarial Robustness Through the Lens of Convolutional Filters,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR) Workshops_ , June 2022, pp. 139–147.


# Floquet dynamics of Rabi model beyond the counterrotating hybridized
rotating wave method

Yingying Han<EMAIL_ADDRESS>Shenzhen Key Laboratory of Ultraintense
Laser and Advanced Material Technology, Center for Intense Laser Application
Technology, and College of Engineering Physics, Shenzhen Technology
University, Shenzhen 518118, China  Shuanghao Zhang  Shenzhen Key Laboratory
of Ultraintense Laser and Advanced Material Technology, Center for Intense
Laser Application Technology, and College of Engineering Physics, Shenzhen
Technology University, Shenzhen 518118, China  Meijuan Zhang  Shenzhen Key
Laboratory of Ultraintense Laser and Advanced Material Technology, Center for
Intense Laser Application Technology, and College of Engineering Physics,
Shenzhen Technology University, Shenzhen 518118, China  Institute of
Theoretical Physics, Shanxi University, Taiyuan 030006, China  Q. Guan
Department of Physics and Astronomy, Washington State University, Pullman,
Washington 99164-2814, USA  Wenxian Zhang  Key Laboratory of Artificial Micro-
and Nano-structures of Ministry of Education, and School of Physics and
Technology, Wuhan University, Wuhan, Hubei 430072, China  Weidong Li  Shenzhen
Key Laboratory of Ultraintense Laser and Advanced Material Technology, Center
for Intense Laser Application Technology, and College of Engineering Physics,
Shenzhen Technology University, Shenzhen 518118, China

###### Abstract

Monochromatically driven two-level systems (i.e., Rabi models) are ubiquitous
in various fields of physics. Though they have been exactly solved, the
physical pictures in these exact solutions are not clear. Recently,
approximate analytical solutions with neat physics have been obtained by using
the counterrotating hybridized rotating wave (CHRW) method, which has been
proven to be effective over a wider range of parameters than the previous
analytical solutions. However, the CHRW depends on a parameter $\xi$, which
has no solution in some regimes. Here we combine the double-unitary-
transformation approach with the generalized Van Vleck nearly degenerate
perturbation theory, and present approximate analytical results with clear
physics for almost all parameter regimes, which agree well with the numerical
solutions and the previous experimental results. Moreover, the dynamic
frequencies of the Rabi model are regular, and the frequency with the highest
Fourier amplitude changes from the Rabi frequency to $2n\omega$ with driving
frequency $\omega$ and integer n, as the driving intensity increases from weak
to deep-strong. In addition, we further explore the Floquet dynamics of the
dissipative open Rabi model. Remarkably, the dissipations are tunable in the
rotating frame, and the approximate analytical results obtained by our method
are in good agreement with the numerical results in the strong driving regime.
These results pave the way to quantum control using strong and deep-strong
driving with applications in quantum technologies.

##  I Introduction

Optical-control of materials based on Floquet engineering has been attracting
great interest ranging from the realization of novel phenomena [1, 2, 3, 4, 5,
6, 7, 8] to optical suppression of decoherence [9, 10, 11]. The prototype of
optically controlled quantum system is the Rabi model, describing a two-level
system driven by a monochromatic field with amplitude $A$ and frequency
$\omega$, which is common in different physical setups, ranging from quantum
optics to condensed matter and quantum information [12, 13]. For a weak and
near-resonance driving, one usually invokes the rotating-wave approximation
(RWA) and results show that the dynamics of the Rabi models are periodic
oscillations with Rabi frequency [12, 14]. In fact, the dynamics of the Rabi
model are much richer than those under the RWA, such as Rabi oscillations’
collapse and revival [15]. To observe these and other novel phenomena in the
lab, a strong (deep-strong) driving regime is required, that is, the driving
amplitude $A$ has to be comparable or larger than the system’s decoherence
rates (transition frequency of the system) [16, 17, 18, 19, 20, 21, 22, 23].
Beyond the weak driving regime, the RWA breaks down, and the frequency
components of the dynamics are no longer a single Rabi frequency and become
analytically difficult to predict.

An exact analytical solution to this Rabi model in terms of two known special
functions has been presented in Ref. [24], and an iterating approach for
strong-coupling periodically driven two-level systems is presented in Ref.
[25]. Though they all show analytical results that agree well with the
numerical results, it is difficult to analytically provide the frequency
components of the system dynamics, which are demanded in experiments for the
strongly driven system [26]. Therefore, analytical solutions with clear
physics are also needed. Counterrotating hybridized rotating wave (CHRW)
approach is beyond the traditional RWA and remains the RWA form with a
renormalized tunneling strength and a modified driving [27]. However, it is
based on a unitary transformation with a parameter $\xi$, which has no
solution in some situations. In this work, we present approximate analytical
results of the Rabi model with clear physics for almost all parameter regimes
by combining the double-unitary-transformation (DUT) approach with the
generalized Van Vleck (GVV) nearly degenerate perturbation theory. Our
analytical results agree well with the numerical solutions obtained by Floquet
theory and the previous experimental results [26]. Moreover, we further
investigate the effectiveness of the GVV method in the dissipative open Rabi
model and find that it is still applicable in the strong driving regime.

The paper is organized as follows. The formalism of Floquet dynamics is
presented in Sec. II. We explore the Floquet dynamics of the closed Rabi model
in Sec. III and the dissipative open Rabi model in Sec. IV. Conclusions are
given in the last section of the article.

##  II Floquet dynamics

In this work, we investigate the Rabi model in the Floquet picture: copy
infinite two-level systems and then shift the energy separation between each
copy by $\hbar\omega$ with Planck’s constant $\hbar$. Note that hereafter
$\hbar=1$. This correspondence allows one to view Rabi oscillations and
Floquet states from the simpler perspective of their time independent-problem
equivalent [28]. A periodic time-dependent Hamiltonian $H(t)=H(t+T)$ with
period $T$ can be transformed to an equivalent time-independent infinite-
dimensional Floquet matrix eigenvalue problem. The Hamiltonian $H(t)$ has
Fourier components of $\omega$ with $\omega=2\pi/T$,

$H(t)=\sum_{n}H^{[n]}\exp(-in\omega t),$ (1)

where the component $H^{[n]}$ is represented in the Floquet-state
$|\alpha,n\rangle=|\alpha\rangle\bigotimes|n\rangle$ with the system index
$\alpha$ (Note that in the generalized Floquet formalism $\alpha$ can be the
N-level system index, but in this work, the Rabi model is restricted to N=2)
and the Fourier index $n$ that runs from $-\infty$ to $\infty$ [29]. According
to the Floquet theory [30], the elements of the infinite-dimensional Floquet
matrix $H_{F}$ are defined by

$\langle\alpha,n|H_{F}|\beta,m\rangle=H_{\alpha\beta}^{[n-m]}+n\omega\delta_{\alpha\beta}\delta_{nm},$
(2)

with integers $n$ and $m$. The Floquet matrix is then diagonalized,

$H_{F}|\varepsilon_{\gamma l}\rangle=q_{\gamma l}|\varepsilon_{\gamma
l}\rangle,$ (3)

where $q_{\gamma l}$ is the quasienergy eigenvalue and $|\varepsilon_{\gamma
l}\rangle$ is the corresponding eigenvector. An initial state
$|\Psi(0)\rangle$ can be written as

$|\Psi(0)\rangle=\sum_{\gamma l}a_{\gamma l}|\varepsilon_{\gamma l}\rangle,$
(4)

with $a_{\gamma l}=\langle\varepsilon_{\gamma l}|\Psi(0)\rangle$. The wave
function at time $t$ is

$|\Psi(t)\rangle\\!=\\!\sum_{\gamma l}a_{\gamma l}e^{-iq_{\gamma
l}t}|\varepsilon_{\gamma l}\rangle\\!=\\!\sum_{\gamma l}e^{-iq_{\gamma
l}t|\varepsilon_{\gamma l}\rangle\langle\varepsilon_{\gamma
l}|\Psi(0)\rangle}.$ (5)

For a given initial state $|\Psi(0)\rangle=|\alpha,0\rangle$, the time-
averaged transition probability from $|\alpha\rangle$ to
$|\alpha^{\prime}\rangle$ can be calculated:

$\displaystyle P_{\alpha\rightarrow\alpha^{\prime}}(t)$ $\displaystyle=$
$\displaystyle|\sum_{n}\langle\alpha^{\prime},n|\Psi(t)\rangle|^{2}$ (6)
$\displaystyle=$ $\displaystyle|\sum_{n}\sum_{\gamma l}e^{-iq_{\gamma
l}t}\langle{\alpha^{\prime},}n|\varepsilon_{\gamma
l}\rangle\langle\varepsilon_{\gamma l}|\alpha,0\rangle|^{2}.$

According to the Floquet theory [31], the eigenvalues of the Floquet matrix
exhibit translational symmetry, i.e., $q_{\gamma n}=q_{\gamma 0}+n\omega$. As
a result, the probability to find the system in state
$|\alpha^{\prime}\rangle$ is expected to show oscillatory behavior with
frequencies $n\omega$ and $\pm\Delta\varepsilon+n\omega$. Here,
$\Delta\varepsilon=q_{\gamma n}-q_{\gamma^{\prime}n}$. In the next section, we
will focus on the analytical results of these frequencies.

##  III Closed Rabi model

The Rabi model describes a two-level system, denoted by ground state
$|0\rangle$ and excited state $|1\rangle$, driven by a harmonic driving with
amplitude $A$ and frequency $\omega$,

$H_{\text{Rabi}}(t)=\frac{\omega_{0}}{2}\sigma_{z}+\frac{A}{2}\cos(\omega
t)\sigma_{x},$ (7)

where $\sigma_{x,y,z}$ is the usual Pauli matrix and $\omega_{0}$ the
transition frequency of the system. Recently, strongly driven quantum systems
have attracted considerable attention [32, 21], and the traditional RWA is not
valid in this regime. Thus, we will analytically investigate the Floquet
dynamics of the Rabi model by using the CHRW [33] and GVV methods.

###  III.1  Counterrotating hybridized rotating wave

The essence of the CHRW method is a unitary transformation

$U_{1}(t)=\text{exp}\left[-i\frac{A\xi}{2\omega}\sin(\omega
t)\sigma_{x}\right],$ (8)

with parameter $\xi\in[0,1]$ to be determined later [27]. The Hamiltonian in
Eq. (7) after the transformation is

$\displaystyle H_{1}(t)$ $\displaystyle=$ $\displaystyle
U_{1}^{\dagger}(t)H_{\text{Rabi}}(t)U_{1}(t)-iU_{1}^{\dagger}(t)\frac{\partial
U_{1}(t)}{\partial t}$ (9) $\displaystyle=$
$\displaystyle\frac{\omega_{0}}{2}\left\\{\cos\left[\frac{A\xi}{\omega}\sin(\omega
t)\right]\sigma_{z}+\sin\left[\frac{A\xi}{\omega}\sin(\omega
t)\right]\sigma_{y}\right\\}$ $\displaystyle+$
$\displaystyle\frac{A}{2}(1-\xi)\cos(\omega t)\sigma_{x}.$

Note that after the transformation $U_{1}(t)$, the basis states of the
Hamiltonian in Eq. (9) are

$\displaystyle|s_{1}(t)\rangle$ $\displaystyle\\!=\\!$
$\displaystyle\cos\\!\left[\frac{A\xi}{2\omega}\sin(\omega
t)\right]\\!|1\rangle\\!+\\!i\sin\\!\left[\frac{A\xi}{2\omega}\sin(\omega
t)\right]\\!|0\rangle,$ $\displaystyle|s_{0}(t)\rangle$ $\displaystyle\\!=\\!$
$\displaystyle i\sin\\!\left[\frac{A\xi}{2\omega}\sin(\omega
t)\right]\\!|1\rangle\\!+\\!\cos\\!\left[\frac{A\xi}{2\omega}\sin(\omega
t)\right]\\!|0\rangle.$

Using $\text{exp}[iA\xi/\omega\sin(\omega
t)]=\sum_{n=-\infty}^{\infty}J_{n}(A\xi/\omega)\text{exp}(in\omega t)$ with
$n$th-order Bessel function of the first kind $J_{n}(\cdot)$, Eq. (9) can be
divided into three parts
$H_{1}(t)=H_{1}^{\prime}+H_{1}^{\prime\prime}(t)+H_{1}^{\prime\prime\prime}(t)$,
where

$\displaystyle H_{1}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{\omega_{0}}{2}J_{0}\left(\frac{A\xi}{\omega}\right)\sigma_{z},$
(11) $\displaystyle H_{1}^{\prime\prime}(t)$ $\displaystyle=$
$\displaystyle\frac{A}{2}(1-\xi)\cos(\omega
t)\sigma_{x}+\omega_{0}J_{1}\left(\frac{A\xi}{\omega}\right){\sin(\omega
t)}\sigma_{y},$

and
$H_{1}^{\prime\prime\prime}(t)=H_{1}(t)-H_{1}^{\prime}-H_{1}^{\prime\prime}(t)$
includes all higher-order harmonic terms $\sin(n\omega t)$ and $\cos(n\omega
t)$ with $n\geq 2$, which are ignored in the CHRW method. Therefore,
$H_{1}(t)\simeq H_{1}^{\prime}+H_{1}^{\prime\prime}(t)$ and can be rewritten
as

$H_{1}(t)=\frac{\omega_{0}}{2}J_{0}\left(\frac{A\xi}{\omega}\right)\sigma_{z}+\frac{\tilde{A}}{4}(e^{-i\omega
t}\sigma_{+}+e^{i\omega t}\sigma_{-}),$ (13)

where

$\frac{A}{2}(1-\xi)=\omega_{0}J_{1}\left(\frac{A\xi}{\omega}\right)\equiv\frac{\tilde{A}}{4},$
(14)

with $\xi\in[0,1]$. Clearly, Eq. (13) possesses a RWA-like form with a
renormalized transition frequency $J_{0}(A\xi/\omega)\omega_{0}$ and a
renormalized driving strength $\tilde{A}$.

For an initial state $|0\rangle$, the population of excited state $P_{1}(t)$
can be easily given as (see Appendix A)

$\displaystyle P_{1}(t)$ $\displaystyle=$ $\displaystyle
C_{0}+C_{1}\cos(\tilde{\Omega}t)+\sum_{n=1}^{\infty}[C_{2_{n}}\cos(2n\omega
t+\tilde{\Omega}t)$ (15) $\displaystyle+$ $\displaystyle
C_{3_{n}}\cos(2n\omega t-\tilde{\Omega}t)+C_{4_{n}}\cos(2n\omega t)],$

with effective Rabi frequency

$\tilde{\Omega}=\sqrt{\tilde{\Delta}^{2}+\tilde{A}^{2}/4},$ (16)

and renormalized detuning
$\tilde{\Delta}=J_{0}(A\xi/\omega)\omega_{0}-\omega$. The coefficients
$C_{0}$, $C_{1}$, $C_{2_{n}}$, $C_{3_{n}}$ and $C_{4_{n}}$ with
$n=1,2,3,\cdots,\infty$ in Eq. (15) are the Fourier amplitudes of the cosine
functions with frequencies $0$, $\tilde{\Omega}$, $2n\omega+\tilde{\Omega}$,
$2n\omega-\tilde{\Omega}$ and $2n\omega$ respectively, and they are shown in
Appendix A. Form Eq. (15), we find that the dynamic frequency components of
the Rabi model are infinite, and they exhibit translational symmetry with a
translation of $2n\omega$, which is different from the general $n\omega$
discussed in Sec. II. Although the even photon number (i.e., $2n\omega$), has
been discussed before [31, 34, 35, 36], here we present the analytical results
of all the dynamic frequencies and corresponding amplitudes, which, to our
knowledge, have not been presented directly before. The reason for the even
photon number is that the periodic-driven signal in the Rabi model (the cosine
signal) has additional symmetry [26].

Figure 1: Coefficients $C_{0}$, $C_{1}$, $C_{4_{1}}$, $C_{4_{2}}$,
$C_{2_{1}}$, $C_{2_{2}}$, $C_{3_{1}}$ and $C_{3_{2}}$ in Eqs. (15) as a
function of $A$ with $\omega/\omega_{0}=1$ (a) and $\omega/\omega_{0}=0.6$
(b), obtained from Eq. (45)-(49). The gray areas denote no results because
there are no reasonable $\xi$ based on Eq. (14), i.e., no solution for
$\xi\in[0,1]$ satisfies
$\frac{A}{2}(1-\xi)=\omega_{0}J_{1}\left(\frac{A\xi}{\omega}\right)$.

Figure 2:  The number of solutions to $\xi\in[0,1]$ that satisfy equation
$\frac{A}{2}(1-\xi)=\omega_{0}J_{1}\left(\frac{A\xi}{\omega}\right)$.

The magenta and red dashed lines are $\omega/\omega_{0}=0.6$ and
$\omega/\omega_{0}=1$, respectively. The yellow stars on the magenta and red
dashed lines represent the maximum driving strength (i.e.,
$A/\omega_{0}=4.18$) that were realized in experiments [26].

The effectiveness of CHRW has been demonstrated [27], and here we directly use
this method to analyze the major frequencies of system dynamics at different
driving intensities, which have not been studied previously. Figure 1 shows
the variation of Fourier amplitude [i.e., $C_{0}$, $C_{1}$, $C_{4_{1}}$,
$C_{4_{2}}$, $C_{2_{1}}$, $C_{2_{2}}$, $C_{3_{1}}$ and $C_{3_{2}}$ in Eq.
(15)] with the driving intensity $A$. For a resonant situation in Fig. 1(a),
when $A/\omega_{0}\rightarrow 0$, the Floquet dynamic is a constant of 0.5
plus a Rabi oscillation with an amplitude of 0.5 and a frequency of
$\tilde{\Omega}$. As $A$ increases, the Fourier amplitude of the Rabi
oscillation with frequency $\tilde{\Omega}$ decreases, while the Fourier
amplitudes of the high-frequency oscillation terms (i.e.,
$2n\omega\pm\tilde{\Omega}$ and $2n\omega$) increase and gradually exceed the
former, indicating that the RWA without high-frequency oscillation terms is no
longer valid. In addition, in most intervals of the deep-strong driving regime
(i.e., $A/\omega_{0}\gg 1$), the Fourier amplitudes of $2n\omega$ terms are
larger than those of other terms and become the main oscillation terms of the
Floquet dynamics. Moreover, we show the variation of the Fourier amplitudes
with $A$ for the off-resonance situation in Fig. 1 (b). When
$A/\omega\rightarrow 0$, the transition probability $P_{1}(t)\rightarrow 0$,
which is reasonable for non-resonant weak coupling. The amplitude of the Rabi
frequency $C_{1}$ increases with the increase of $A$, reaching its maximum at
about $A/\omega_{0}=1$, and then gradually decreases. Similarly, $2n\omega$
terms become the main oscillation terms of the Floquet dynamics in the most
intervals of deep-strong driving regime. Note that Eq. (15) has no solution
with some parameters in the gray areas as shown in Fig. 1, which we will
discuss in the next paragraph.

We present the valid conditions of the CHRW in this paragraph. From Eq. (14),
we see that the solution of $\xi$ is the key of the CHRW. Figure 2 shows the
number of solutions for the parameter $\xi\in[0,1]$ that satisfy equation
$A(1-\xi)/2=\omega_{0}J_{1}\left(A\xi/\omega\right)$. We find that with
parameters $\omega/\omega_{0}=0.27~{}(0.15),A/\omega_{0}=1.41~{}(1.35)$, there
are two (three) solutions, i.e., $\xi=0.253,0.6043$
($\xi=0.1691,0.2837,0.8256$), and so on. In these regions where $\xi$ has no
unique solution, the CHRW method is no longer applicable. Also in some
parameter ranges, no $\xi$ belongs to $[0,1]$ that satisfies
$A(1-\xi)/2=\omega_{0}J_{1}\left(A\xi/\omega\right)$, such as $3.84\leq
A/\omega_{0}\leq 7.01$ with $\omega/\omega_{0}=1$ (shown by the red dashed
line) and $2.30\leq A/\omega_{0}\leq 4.20$ with $\omega/\omega_{0}=0.6$ (shown
by the magenta dashed line). However, the maximum driving intensity
experimentally studied reached $A/\omega_{0}=4.18$, as the yellow stars shown
in Fig. 2. Therefore, analytical results with clear physics suitable for the
whole driving-strength range need to be solved urgently.

###  III.2 Double-unitary-transformation

In this section, we use the DUT approach to investigate the Floquet dynamics
of the Rabi model in Eq. (7). DUT is a powerful method for dealing with
transverse (off-diagonal) and periodically driven quantum systems (such as the
Rabi model), which makes it possible to obtain analytical results by using
perturbation theory. The first unitary transformation is
$U_{2}=\text{exp}(-i\pi\sigma_{y}/4)$, and the Hamiltonian in Eq. (7) after
this transformation is

$H_{2}(t)=\frac{A}{2}\cos(\omega t)\sigma_{z}-\frac{\omega_{0}}{2}\sigma_{x},$
(17)

with basis $|\pm\rangle=(|1\rangle\pm|0\rangle)/\sqrt{2}$. The second unitary
transformation is

$U_{3}(t)=\text{exp}\left[-i\frac{A}{2\omega}\sin(\omega t)\sigma_{z}\right],$
(18)

and the Hamiltonian in Eq. (17) after this transformation becomes

$\displaystyle H_{3}(t)$ $\displaystyle=$
$\displaystyle\frac{-\omega_{0}}{2}\sum_{n=-\infty}^{\infty}J_{n}e^{-in\omega
t}[(-1)^{n}\sigma_{+}+\sigma_{-}],$ (19)

where $J_{n}(A/\omega)$ abbreviate as $J_{n}$ and $\sigma_{\pm}=(\sigma_{x}\pm
i\sigma_{y})/2$. Note that the basis states of Eq. (19) are

$\displaystyle|s^{\prime}_{+}(t)\rangle=\frac{\text{exp}[\frac{iA\sin(\omega
t)}{2\omega}]|1\rangle+\text{exp}[\frac{-iA\sin(\omega
t)}{2\omega}]|0\rangle}{\sqrt{2}},$ (20)
$\displaystyle|s^{\prime}_{-}(t)\rangle=\frac{\text{exp}[\frac{iA\sin(\omega
t)}{2\omega}]|1\rangle-\text{exp}[\frac{-iA\sin(\omega
t)}{2\omega}]|0\rangle}{\sqrt{2}}.$ (21)

According to Eq. (1), the Floquet matrix block of Eq. (19) is

$\displaystyle H^{[n]}$ $\displaystyle=$
$\displaystyle\frac{-\omega_{0}}{2}J_{n}[(-1)^{n}\sigma_{+}+\sigma_{-}].$ (22)

Combined Eq. (22) with Eq. (2), the Floquet matrix of Eq. (19) is

$H_{F}=\left(\begin{array}[]{c c c|c c|c c|c c|c c
c}\ddots&&&&&&&&&&&\iddots\\\
\hline\cr&-2\omega&\frac{-\omega_{0}}{2}J_{0}&0&\frac{-\omega_{0}}{2}J_{1}&0&\frac{-\omega_{0}}{2}J_{2}&0&\frac{-\omega_{0}}{2}J_{3}&0&\frac{-\omega_{0}}{2}J_{4}&\\\
&\frac{-\omega_{0}}{2}J_{0}&-2\omega&\frac{\omega_{0}}{2}J_{1}&0&\frac{-\omega_{0}}{2}J_{2}&0&\frac{\omega_{0}}{2}J_{3}&0&\frac{-\omega_{0}}{2}J_{4}&0&\\\
\hline\cr&0&\frac{\omega_{0}}{2}J_{1}&-\omega&\frac{-\omega_{0}}{2}J_{0}&0&\frac{-\omega_{0}}{2}J_{1}&0&\frac{-\omega_{0}}{2}J_{2}&0&\frac{-\omega_{0}}{2}J_{3}&\\\
&\frac{-\omega_{0}}{2}J_{1}&0&\frac{-\omega_{0}}{2}J_{0}&-\omega&\frac{\omega_{0}}{2}J_{1}&0&\frac{-\omega_{0}}{2}J_{2}&0&\frac{\omega_{0}}{2}J_{3}&0&\\\
\hline\cr&0&\frac{-\omega_{0}}{2}J_{2}&0&\frac{\omega_{0}}{2}J_{1}&0&\frac{-\omega_{0}}{2}J_{0}&0&\frac{-\omega_{0}}{2}J_{1}&0&\frac{-\omega_{0}}{2}J_{2}&\\\
&\frac{-\omega_{0}}{2}J_{2}&0&\frac{-\omega_{0}}{2}J_{1}&0&\frac{-\omega_{0}}{2}J_{0}&0&\frac{\omega_{0}}{2}J_{1}&0&\frac{-\omega_{0}}{2}J_{2}&0&\\\
\hline\cr&0&\frac{\omega_{0}}{2}J_{3}&0&\frac{-\omega_{0}}{2}J_{2}&0&\frac{\omega_{0}}{2}J_{1}&\omega&\frac{-\omega_{0}}{2}J_{0}&0&\frac{-\omega_{0}}{2}J_{1}&\\\
&\frac{-\omega_{0}}{2}J_{3}&0&\frac{-\omega_{0}}{2}J_{2}&0&\frac{-\omega_{0}}{2}J_{1}&0&\frac{-\omega_{0}}{2}J_{0}&\omega&\frac{\omega_{0}}{2}J_{1}&0&\\\
\hline\cr&0&\frac{-\omega_{0}}{2}J_{4}&0&\frac{\omega_{0}}{2}J_{3}&0&\frac{-\omega_{0}}{2}J_{2}&0&\frac{\omega_{0}}{2}J_{1}&2\omega&\frac{-\omega_{0}}{2}J_{0}&\\\
&\frac{-\omega_{0}}{2}J_{4}&0&\frac{-\omega_{0}}{2}J_{3}&0&\frac{-\omega_{0}}{2}J_{2}&0&\frac{-\omega_{0}}{2}J_{1}&0&\frac{-\omega_{0}}{2}J_{0}&2\omega&\\\
\hline\cr\iddots&&&&&&&&&&&\ddots\end{array}\right)\begin{array}[]{c
c}\leftarrow&|s^{\prime}_{-},-2\rangle\\\
\leftarrow&|s^{\prime}_{+},-2\rangle\\\
\leftarrow&|s^{\prime}_{-},-1\rangle\\\
\leftarrow&|s^{\prime}_{+},-1\rangle\\\ \leftarrow&|s^{\prime}_{-},0\rangle\\\
\leftarrow&|s^{\prime}_{+},0\rangle\\\ \leftarrow&|s^{\prime}_{-},+1\rangle\\\
\leftarrow&|s^{\prime}_{+},+1\rangle\\\
\leftarrow&|s^{\prime}_{-},+2\rangle\\\
\leftarrow&|s^{\prime}_{+},+2\rangle\\\ \end{array}$ (23)

where $J_{k}=J_{k}(A/\omega)$. Then we perform a basis transformation with

$S=\left(\begin{array}[]{c c c|c c|c c|c c|c c c}\ddots&&&&&&&&&&&\iddots\\\
\hline\cr&1&1&0&0&0&0&0&0&0&0&\\\ &1&-1&0&0&0&0&0&0&0&0\\\
\hline\cr&0&0&1&1&0&0&0&0&0&0\\\ &0&0&1&-1&0&0&0&0&0&0\\\
\hline\cr&0&0&0&0&1&1&0&0&0&0&\\\ &0&0&0&0&1&-1&0&0&0&0&\\\
\hline\cr&0&0&0&0&0&0&1&1&0&0&\\\ &0&0&0&0&0&0&1&-1&0&0&\\\
\hline\cr&0&0&0&0&0&0&0&0&1&1&\\\ &0&0&0&0&0&0&0&0&1&-1&\\\
\hline\cr\iddots&&&&&&&&&&&\ddots\end{array}\right),$ (24)

and Eq. (23) becomes

$H^{\prime}_{F}=\left(\begin{array}[]{c c c|c c|c c|c c|c c
c}\ddots&&&&&&&&&&&\iddots\\\
\hline\cr&-2\omega-\frac{J_{0}\omega_{0}}{2}&0&0&\frac{\omega_{0}J_{1}}{2}&\frac{-\omega_{0}J_{2}}{2}&0&0&\frac{\omega_{0}J_{3}}{2}&\frac{-\omega_{0}J_{4}}{2}&0&\\\
&0&-2\omega+\frac{J_{0}\omega_{0}}{2}&\frac{-\omega_{0}J_{1}}{2}&0&0&\frac{\omega_{0}J_{2}}{2}&\frac{-\omega_{0}J_{3}}{2}&0&0&\frac{\omega_{0}J_{4}}{2}&\\\
\hline\cr&0&\frac{-\omega_{0}J_{1}}{2}&-\omega-\frac{J_{0}\omega_{0}}{2}&0&0&\frac{\omega_{0}J_{1}}{2}&\frac{-\omega_{0}J_{2}}{2}&0&0&\frac{\omega_{0}J_{3}}{2}&\\\
&\frac{\omega_{0}J_{1}}{2}&0&0&-\omega+\frac{J_{0}\omega_{0}}{2}&\frac{-\omega_{0}J_{1}}{2}&0&0&\frac{\omega_{0}J_{2}}{2}&\frac{-\omega_{0}J_{3}}{2}&0&\\\
\hline\cr&\frac{-\omega_{0}J_{2}}{2}&0&0&\frac{-\omega_{0}J_{1}}{2}&\frac{-\omega_{0}J_{0}}{2}&0&0&\frac{\omega_{0}J_{1}}{2}&\frac{-\omega_{0}J_{2}}{2}&0\\\
&0&\frac{\omega_{0}J_{2}}{2}&\frac{\omega_{0}J_{1}}{2}&0&0&\pagecolor{blue!20}\frac{\omega_{0}J_{0}}{2}&\pagecolor{blue!20}\frac{-\omega_{0}J_{1}}{2}&0&0&\frac{\omega_{0}J_{2}}{2}&\\\
\hline\cr&0&\frac{-\omega_{0}J_{3}}{2}&\frac{-\omega_{0}J_{2}}{2}&0&0&\pagecolor{blue!20}\frac{-\omega_{0}J_{1}}{2}&\pagecolor{blue!20}\omega-\frac{\omega_{0}J_{0}}{2}&0&0&\frac{\omega_{0}J_{1}}{2}&\\\
&\frac{\omega_{0}J_{3}}{2}&0&0&\frac{\omega_{0}J_{2}}{2}&\frac{\omega_{0}J_{1}}{2}&0&0&\omega+\frac{J_{0}\omega_{0}}{2}&\frac{-\omega_{0}J_{1}}{2}&0&\\\
\hline\cr&\frac{-\omega_{0}J_{4}}{2}&0&0&\frac{-\omega_{0}J_{3}}{2}&\frac{-\omega_{0}J_{2}}{2}&0&0&\frac{-\omega_{0}J_{1}}{2}&2\omega-\frac{J_{0}\omega_{0}}{2}&0&\\\
&0&\frac{\omega_{0}J_{4}}{2}&\frac{\omega_{0}J_{3}}{2}&0&0&\frac{\omega_{0}J_{2}}{2}&\frac{\omega_{0}J_{1}}{2}&0&0&2\omega+\frac{J_{0}\omega_{0}}{2}&\\\
\hline\cr\iddots&&&&&&&&&&&\ddots\end{array}\right).$ (25)
$\begin{array}[]{cccc c cccc cc c c c cc c c c c
c}&&&&&\uparrow&\uparrow&&&\uparrow&\uparrow&&\uparrow&\uparrow&\uparrow&\uparrow&&&\uparrow&\uparrow&\\\
&&&&&|s^{\prime}_{0},-2\rangle&|s^{\prime}_{1},-2\rangle&&&|s^{\prime}_{0},-1\rangle&|s^{\prime}_{1},-1\rangle&&|s^{\prime}_{0},0\rangle&|s^{\prime}_{1},0\rangle&|s^{\prime}_{0},1\rangle&|s^{\prime}_{1},1\rangle&&&|s^{\prime}_{0},2\rangle&|s^{\prime}_{1},2\rangle&\end{array}$
(26)

with Floquet basis

$\displaystyle|s^{\prime}_{1}(t)\rangle\\!$ $\displaystyle=$
$\displaystyle\\!\cos\left[\frac{A}{2\omega}\sin(\omega
t)\right]\\!|1\rangle\\!+\\!i\sin\left[\frac{A}{2\omega}\sin(\omega
t)\right]|0\rangle$ $\displaystyle|s^{\prime}_{0}(t)\rangle\\!$
$\displaystyle=$ $\displaystyle\\!i\sin\left[\frac{A}{2\omega}\sin(\omega
t)\right]\\!|1\rangle\\!+\\!\cos\left[\frac{A}{2\omega}\sin(\omega
t)\right]|0\rangle.$

From the matrix structure of $H^{\prime}_{F}$ in Eq. (25), one sees that
$|s^{\prime}_{1},0\rangle$ couples to $|s^{\prime}_{0},1\rangle$ via an off-
diagonal term of $-\omega_{0}J_{1}/2$, and this block matrix is shown in light
blue. By tuning the frequency of the driving field (i.e., $\omega$), the
Floquet state $|s^{\prime}_{1},0\rangle$ can become nearly degenerate with
$|s^{\prime}_{0},1\rangle$, namely, $\omega\approx J_{0}\omega_{0}$. Following
the standard GVV nearly degenerate perturbation theory and to the second order
[37, 38], Eq. (25) is reduced to a $2\times 2$ matrix by including all other
off-resonant coupling channels as perturbation terms (see Appendix B):

$\displaystyle H_{\rm GVV}=\left(\begin{array}[]{c c
}\frac{J_{0}\omega_{0}}{2}+\delta_{1^{\prime}}&\frac{-J_{1}\omega_{0}}{2}+\delta_{1^{\prime}0^{\prime}}\vspace{2mm}\\\
\frac{-J_{1}\omega_{0}}{2}+\delta_{0^{\prime}1^{\prime}}&\omega-\frac{J_{0}\omega_{0}}{2}+\delta_{0^{\prime}}\end{array}\right),$
(30)

where

$\displaystyle\delta_{1^{\prime}}$ $\displaystyle=$
$\displaystyle\sum_{\begin{subarray}{c}k=-\infty\\\ k\neq
0\end{subarray}}^{\infty}\left[\frac{(J_{2k+1}\omega_{0}/2)^{2}}{J_{0}\omega_{0}-(2k+1)\omega}+\frac{(J_{2k}\omega_{0}/2)^{2}}{-2k\omega}\right],$
(31) $\displaystyle\delta_{0^{\prime}}$ $\displaystyle\\!=\\!$
$\displaystyle\sum_{\begin{subarray}{c}k=-\infty\\\ k\neq
0\end{subarray}}^{\infty}\left[\frac{-(J_{2k-1}\omega_{0}/2)^{2}}{J_{0}\omega_{0}+(2k-1)\omega}\\!+\\!\frac{(J_{2k}\omega_{0}/2)^{2}}{-2k\omega}\right],$
(32) $\displaystyle\delta_{1^{\prime}0^{\prime}}$ $\displaystyle\\!=\\!$
$\displaystyle\sum_{\begin{subarray}{c}k=-\infty\\\ k\neq
0\end{subarray}}^{\infty}\left[\frac{-J_{2k-1}J_{2k}\omega_{0}^{2}/4}{J_{0}\omega_{0}\\!+\\!(2k-1)\omega}+\frac{J_{2k+1}J_{2k}\omega_{0}^{2}}{-8k\omega}\right],$
(33) $\displaystyle\delta_{0^{\prime}1^{\prime}}$ $\displaystyle\\!=\\!$
$\displaystyle\sum_{\begin{subarray}{c}k=-\infty\\\ k\neq
0\end{subarray}}^{\infty}\left[\frac{J_{2k+1}J_{2k}\omega_{0}^{2}/4}{J_{0}\omega_{0}-(2k+1)\omega}\\!+\\!\frac{J_{2k-1}J_{2k}\omega_{0}^{2}}{-8k\omega}\right],$
(34)

which result from the higher-order harmonic terms. Therefore, here our
approach includes the higher-order harmonic terms that are ignored in CHRW. In
Appendix B, we numerically prove that
$\delta_{0^{\prime}}=-\delta_{1^{\prime}}$,
$\delta_{1^{\prime}0^{\prime}}=\delta_{0^{\prime}1^{\prime}}$ within the
parameter range of our research and also demonstrate that the second-order
perturbation results are more suitable for the exact numerical results than
the first-order perturbation results.

Diagonalize Eq. (30) to obtain the eigenvalues, and the Rabi frequency is the
difference between the two eigenvalues:

$\Omega=\sqrt{4\delta_{1^{\prime}0^{\prime}}\delta_{0^{\prime}1^{\prime}}+(\delta_{1}^{\prime}-\delta_{0}^{\prime}-\omega)^{2}+B}$
(35)

with
$B=2(\delta_{1}^{\prime}\\!-\\!\delta_{0}^{\prime}-\omega)J_{0}\omega_{0}\\!-\\!2(\delta_{1^{\prime}0^{\prime}}+\delta_{0^{\prime}1^{\prime}})J_{1}\omega_{0}+(J_{0}^{2}+J_{1}^{2})\omega_{0}^{2}$.
Then all frequencies of the Floquet dynamics can be written as
$\pm\Omega+2n\omega$ and $2n\omega$. In Fig. 3, we compare the numerical,
analytical, and the previous experimental results of system dynamics
frequencies in the general off-resonance situation (i.e.,
$\omega/\omega_{0}=0.6$). The numerical results are obtained by Floquet
theory. The approximate analytical results are obtained by using CHRW and GVV,
and the experimental results extracted from Ref. [26]. One can see that the
GVV results are more suitable for the exact numerical results than CHRW,
because GVV results include higher-order harmonic terms through perturbation
theory, which are ignored in CHRW. When $A\in[2.30,4.20]$ or
$A\in[6.10,7.99]$, CHRW has no results as discussed in Sec. III.1. In a
contrast, the GVV results agree well with the experimental results, and the
slight discrepancy may be compensated by considering higher-order perturbation
terms. In addition, the GVV results are not only suitable for resonant and
near-resonant situations, but also for large detuning situations, as shown in
Fig. 4. In fact, the GVV results are suitable for almost all parameter
regimes, except $A/\omega_{0}\rightarrow 0$ and $\omega/\omega_{0}\rightarrow
0$ simultaneously. Therefore, the GVV results provide a better approximation
for the Rabi model with arbitrary strength and detuning.

Figure 3:  Frequency components of transition probability $P_{1}(t)$ as a
function of $A$ with $\omega/\omega_{0}=0.6$. The black lines plot
$\pm\Delta\epsilon+2n\omega$ and $2n\omega$, with quasienergy difference
$\Delta\epsilon$ determined numerically by Floquet matix (see Appendix C). The
red dashed lines plot $\pm\Omega+2n\omega$ and $2n\omega$ with $\Omega$ in Eq.
(35). The blue dash-dot lines plot $\pm\tilde{\Omega}+2n\omega$ and $2n\omega$
with $\tilde{\Omega}$ in Eq. (16). The green dots are the experimental results
extracted from Ref. [26].  Figure 4:  Frequency components of transition
probability $P_{1}(t)$ as a function of $A$ and $\omega$. Same as Fig. 3, the
black lines are the numerical results obtained by Floquet theory and the red
dashed lines are the approximate analytical results $\pm\Omega+2n\omega$ with
$\Omega$ in Eq. (35).

##  IV Dissipative open Rabi model

In this section, we further consider the dissipative open Rabi model with
population damping rate $\Gamma_{ij}$ from state $|i\rangle$ to $|j\rangle$
and the dephasing rate $\gamma_{ii}$ of state $|i\rangle$. With the assumption
of Markovian noise background, the density matrix for the two-level system
with basis $|i\rangle$ and $|j\rangle$ evolves according to the Lindblad
master equation [39, 40],

$\displaystyle\dot{\rho}$ $\displaystyle\\!=\\!$
$\displaystyle-i[H(t),\rho]\\!+\\!\frac{\Gamma_{ij}}{2}\mathcal{D}[|j\rangle\langle
i|]\rho\\!+\\!\frac{\Gamma_{ji}}{2}\mathcal{D}[|i\rangle\langle j|]\rho,$ (36)
$\displaystyle+$ $\displaystyle\gamma_{ii}\mathcal{D}[|i\rangle\langle
i|]\rho\\!+\\!\gamma_{jj}\mathcal{D}[|j\rangle\langle j|]\rho,$

with
$\mathcal{D}[\mathcal{O}]\rho=2\mathcal{O}\rho\mathcal{O}^{{\dagger}}-\mathcal{O}^{{\dagger}}\mathcal{O}\rho-\rho\mathcal{O}^{{\dagger}}\mathcal{O}$.
For the Rabi model, in the lab frame, we only consider the population damping
rate $\Gamma_{10}$ from the excited state $|1\rangle$ to the ground state
$|0\rangle$ and the dephasing rate $\gamma_{11}$ of state $|1\rangle$, i.e.,
$\Gamma_{01}=\gamma_{00}=0$.

The GVV method transforms the Rabi Hamiltonian with basis $|0\rangle$ and
$|1\rangle$ in Eq. (7) into an effective analytically solvable Hamiltonian in
Eq. (30) with basis $|s^{\prime}_{0,1}(t)\rangle$ through two unitary
transformations. The decay and excitation rate between
$|s^{\prime}_{0}(t)\rangle$ and $|s^{\prime}_{1}(t)\rangle$ can be easily
calculated by rewriting the lab frame dissipators (i.e., $\Gamma_{10}$ and
$\gamma_{11}$) in the new basis (see Appendix D),

$\displaystyle\\!\gamma_{s^{\prime}_{1}s^{\prime}_{1}}(t)$
$\displaystyle\\!=\\!$ $\displaystyle\sin^{2}\\!\left[\frac{A\sin(\omega
t)}{\omega}\right]\\!\frac{\Gamma_{10}}{8}\\!+\\!\cos^{4}\\!\left[\frac{A\sin(\omega
t)}{2\omega}\right]\\!\gamma_{11},$ (37)
$\displaystyle\\!\gamma_{s^{\prime}_{0}s^{\prime}_{0}}(t)$
$\displaystyle\\!=\\!$ $\displaystyle\sin^{2}\\!\left[\frac{A\sin(\omega
t)}{\omega}\right]\\!\frac{\Gamma_{10}}{8}\\!+\\!\sin^{4}\\!\left[\frac{A\sin(\omega
t)}{2\omega}\right]\\!\gamma_{11},$ (38)
$\displaystyle\\!\Gamma_{s^{\prime}_{1}s^{\prime}_{0}}(t)$
$\displaystyle\\!=\\!$ $\displaystyle\sin^{2}\\!\left[\frac{A\sin(\omega
t)}{\omega}\right]\\!\frac{\gamma_{11}}{2}\\!+\\!\cos^{4}\\!\left[\frac{A\sin(\omega
t)}{2\omega}\right]\\!\Gamma_{10},$ (39)
$\displaystyle\\!\Gamma_{s^{\prime}_{0}s^{\prime}_{1}}(t)$
$\displaystyle\\!=\\!$ $\displaystyle\sin^{2}\\!\left[\frac{A\sin(\omega
t)}{\omega}\right]\\!\frac{\gamma_{11}}{2}\\!+\\!\sin^{4}\\!\left[\frac{A\sin(\omega
t)}{2\omega}\right]\\!\Gamma_{10}.$ (40)

Equations (37)-(40) show that the decay parameters are time-dependent. One
finds that when $A$ is weak (i.e., $A/\omega\rightarrow 0$),
$\gamma_{s^{\prime}_{1}s^{\prime}_{1}}(t)\approx\gamma_{11}$,
$\gamma_{s^{\prime}_{0}s^{\prime}_{0}}(t)\approx 0$,
$\Gamma_{s^{\prime}_{1}s^{\prime}_{0}}(t)\approx\Gamma_{10}$ and
$\Gamma_{s^{\prime}_{0}s^{\prime}_{1}}(t)\approx 0$, because under this
condition $|s^{\prime}_{0}(t)\rangle\simeq|0\rangle$ and
$|s^{\prime}_{1}(t)\rangle\simeq|1\rangle$. When $A$ is deep-strong (i.e.,
$A/\omega\gg 1$) , the decay rates change greatly and periodically with time,
and sometimes the excitation rate between $|s^{\prime}_{0}\rangle$ and
$|s^{\prime}_{1}\rangle$ is larger than the decay rate between them, i.e.,
$\Gamma_{s^{\prime}_{0}s^{\prime}_{1}}>\Gamma_{s^{\prime}_{1}s^{\prime}_{0}}$.
Similarly, the dephasing rate of $|s^{\prime}_{0}\rangle$ can be larger than
that of $|s^{\prime}_{1}\rangle$, i.e.,
$\gamma_{s^{\prime}_{0}s^{\prime}_{0}}>\gamma_{s^{\prime}_{1}s^{\prime}_{1}}$.
These tunable decay parameters provide a new platform for quantum computing
[41].

Submitting the decay parameters in Eqs. (37)-(40) and the Hamiltonian in Eq.
(30) into the Lindblad master equation in Eq. (36), we obtain the dynamics of
the transition probability $P_{1}(t)$, as the red-dashed lines shown in Fig.
5. To verify the effectiveness of the GVV method, we also present the exact
numerical results obtained by solving the Floquet-Lindblad equation with
Hamiltonian in Eq. (7) and decay parameters $\Gamma_{10}/\omega=1$,
$\gamma_{11}/\omega=0.2$, $\Gamma_{01}/\omega=\gamma_{00}/\omega=0$. We find
that the GVV results agree well with the numerical results in the strong
driving regime, and the small deviations between them may result from higher-
order terms of perturbation theory and the loss of the fast oscillation terms
when rewriting the dissipation rates in Eqs. (37)-(40).

Figure 5:  Dynamics of transition probability $P_{1}$ with
$\omega/\omega_{0}=1$, $A/\omega_{0}=10$ (a), and $\omega/\omega_{0}=3$,
$A/\omega_{0}=10$ (b). The black lines denote the numerical results and the
red dashed lines are the results obtained by GVV. Here the initial state is
$\Psi(0)=|0\rangle$ and decay parameters are $\Gamma_{10}/\omega=1$,
$\gamma_{11}/\omega=0.2$ and $\Gamma_{01}/\omega=\gamma_{00}/\omega=0$.

##  V Discussion and conclusion

In conclusion, by combining DUT with GVV perturbation theories, we present an
approximate analytical result of the Rabi model for almost all parameter
regimes, except $A/\omega_{0}\rightarrow 0$ and $\omega/\omega_{0}\rightarrow
0$ simultaneously. The GVV results agree well with the numerical solutions
obtained by Floquet theory and previous experimental results, and are also
beyond the results obtained by using CHRW, which have no results within some
parameter ranges. Moreover, as the driving intensity increases from weak to
deep-strong, the main frequency component of the Floquet dynamics transitions
from Rabi frequency to $2n\omega$. In addition, we further consider the
dissipative open Rabi model with population damping and dephasing rates, and
the results obtained by GVV Hamiltonian with the regrouping decay rates agree
well with the numerical results in the strong driving regime. Our results
provide a desired theoretical method for studying the strongly driven closed
and open two-level systems. In particular, the development of ultrastrong
laser opens the doorway for light-matter interactions in the strong and deep-
strong coupling regimes.

The approximate analytical results of the Rabi model are helpful for the use
of strong driving for fast quantum gates, specifically qubit state
preparation. Starting with the qubit in its ground state, we apply a
monochromatic field with driving strength $A$ and frequency $\omega$. With the
help of our analytical results, the field strength or duration can be
quantitatively calculated to prepare the target state. In addition, the
tunable decay rates in the rotating frame may have applications in quantum
computing. In particular, the excitation rate between two states can be larger
than the decay rate between them, which is almost impossible in the lab frame.

###### Acknowledgements.

We thank Jianwen Jie for many helpful and intriguing discussions. This work is
supported by the National Natural Science Foundation of China under Grants
Nos. 12205199 and 12274331, and Natural Science Foundation of Top Talent of
SZTU(Grant No. GDRC202202 and GDRC202312).

##  Appendix A Derivation of Eq. (15)

The wave function in basis $|s_{1,0}(t)\rangle$ is
$|\Psi(t)\rangle=c_{1}(t)|s_{1}(t)\rangle+c_{0}(t)|s_{0}(t)\rangle$. The
initial condition after the unitary transformation is invariant,
$|s_{i}(0)\rangle=|i\rangle,i=1,0$. For an initial state in $|0\rangle$, we
can solve the Schr$\ddot{\text{o}}$dinger equation with Hamiltonian in Eq.
(13) easily as given in Ref. [12, 27],

$\displaystyle c_{1}(t)$ $\displaystyle=$
$\displaystyle-\frac{i\tilde{A}}{2\tilde{\Omega}}\sin\left(\frac{\tilde{\Omega}t}{2}\right)e^{i\frac{\omega
t}{2}}$ (41) $\displaystyle c_{0}(t)$ $\displaystyle=$
$\displaystyle\left[\cos\left(\frac{\tilde{\Omega}t}{2}\right)\\!+\\!\frac{i\tilde{\Delta}}{\tilde{\Omega}}\sin\left(\frac{\tilde{\Omega}t}{2}\right)\right]e^{-i\frac{\omega
t}{2}}.$ (42)

The population of excited state $|1\rangle$ can be expressed as

$P_{1}(t)\\!=\\!\left|\cos\left[\frac{A\xi\sin(\omega
t)}{2\omega}\right]c_{1}(t)\\!-\\!i\sin\left[\frac{A\xi\sin(\omega
t)}{2\omega}\right]c_{0}(t)\right|^{2}$ (43)

Substituting Eq. (41) and Eq. (42) into Eq. (43), and regrouping terms with
the same frequency, we have

$\displaystyle P_{1}(t)$ $\displaystyle=$ $\displaystyle
C_{0}+C_{1}\cos(\tilde{\Omega}t)+\sum_{n=1}^{\infty}[C_{2_{n}}\cos(2n\omega
t+\tilde{\Omega}t)$ (44) $\displaystyle+$ $\displaystyle
C_{3_{n}}\cos(2n\omega t-\tilde{\Omega}t)+C_{4_{n}}\cos(2n\omega t)],$

with

$\displaystyle C_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{2}-\frac{\tilde{\Delta}^{2}}{2\tilde{\Omega}^{2}}J_{0}\left(\frac{A\xi}{\omega}\right)+\frac{\tilde{\Delta}\tilde{A}}{4\tilde{\Omega}^{2}}J_{1}\left(\frac{A\xi}{\omega}\right),$
(45) $\displaystyle C_{1}$ $\displaystyle=$
$\displaystyle-\frac{\tilde{A}^{2}}{8\tilde{\Omega}^{2}}J_{0}\left(\frac{A\xi}{\omega}\right)-\frac{\tilde{\Delta}\tilde{A}}{4\tilde{\Omega}^{2}}J_{1}\left(\frac{A\xi}{\omega}\right),$
(46) $\displaystyle C_{2_{n}}$ $\displaystyle=$
$\displaystyle-\left[\frac{\tilde{A}^{2}}{8\tilde{\Omega}^{2}}+\frac{n\omega\tilde{A}}{2A\xi\tilde{\Omega}}\right]J_{2n}\left(\frac{A\xi}{\omega}\right)$
(47) $\displaystyle+$
$\displaystyle\frac{\tilde{\Delta}\tilde{A}}{8\tilde{\Omega}^{2}}\left[J_{2n-1}\left(\frac{A\xi}{\omega}\right)-J_{2n+1}\left(\frac{A\xi}{\omega}\right)\right],$
$\displaystyle C_{3_{n}}$ $\displaystyle=$
$\displaystyle-\left[\frac{\tilde{A}^{2}}{8\tilde{\Omega}^{2}}-\frac{n\omega\tilde{A}}{2A\xi\tilde{\Omega}}\right]J_{2n}\left(\frac{A\xi}{\omega}\right)$
(48) $\displaystyle+$
$\displaystyle\frac{\tilde{\Delta}\tilde{A}}{8\tilde{\Omega}^{2}}\left[J_{2n-1}\left(\frac{A\xi}{\omega}\right)-J_{2n+1}\left(\frac{A\xi}{\omega}\right)\right],$
$\displaystyle C_{4_{n}}$ $\displaystyle=$
$\displaystyle\frac{\tilde{\Delta}\tilde{A}}{4\tilde{\Omega}^{2}}\left[J_{2n+1}\left(\frac{A\xi}{\omega}\right)-J_{2n-1}\left(\frac{A\xi}{\omega}\right)\right]$
(49) $\displaystyle-$
$\displaystyle\frac{\tilde{\Delta}^{2}}{\tilde{\Omega}^{2}}J_{2n}\left(\frac{A\xi}{\omega}\right),$

where we used the formula $\cos[A\sin(\omega
t)/\omega]=\sum_{n=-\infty}^{\infty}J_{n}(A/\omega)\cos(n\omega t)$.

##  Appendix B Derivation of the $2\times 2$ effective Hamiltonian by the GVV
theory

In this section, we reduce the infinite-dimensional Floquet matrix in Eq. (25)
into a $2\times 2$ effective matrix by using GVV perturbation theory [42, 37].
Consider the Floquet states $|s^{\prime}_{1},0\rangle$ nearly degenerate with
$|s^{\prime}_{0},1\rangle$. According to the perturbation theory, we expand
the $2\times 2$ matrix $h$ and its eigenstates $\Phi$ in powers of
$\omega_{0}$. The zeroth and higher-order of $\Phi$ are given by

$\displaystyle\Phi_{1}^{\prime(0)}$ $\displaystyle\\!=\\!$
$\displaystyle|s_{1}^{\prime},0\rangle,~{}~{}\Phi_{0}^{(0)}=|s_{0}^{\prime},1\rangle,$
$\displaystyle\Phi_{1}^{\prime(1)}$ $\displaystyle\\!=\\!$
$\displaystyle\sum_{\begin{subarray}{c}k=-\infty\\\ k\neq
0\end{subarray}}^{\infty}\left[\frac{-J_{2k+1}\omega_{0}/2}{J_{0}\omega_{0}\\!-\\!(2k+1)\omega}|s_{0}^{\prime},2k\\!+\\!1\rangle\\!-\\!\frac{J_{2k}\omega_{0}}{4k\omega}|s_{1}^{\prime},2k\rangle\right],$
$\displaystyle\Phi_{0}^{(1)}$ $\displaystyle\\!=\\!$
$\displaystyle\sum_{\begin{subarray}{c}k=-\infty\\\ k\neq
0\end{subarray}}^{\infty}\left[\frac{-J_{2k-1}\omega_{0}/2}{J_{0}\omega_{0}\\!+\\!(2k-1)\omega}|s_{1}^{\prime},2k\rangle\\!+\\!\frac{J_{2k}\omega_{0}}{4k\omega}|s_{0}^{\prime},2k\\!+\\!1\rangle\right].$

The zeroth and higher-order of $h$ represented by $\Phi$ are

$h^{(0)}=\left(\begin{array}[]{ c c}\frac{J_{0}\omega_{0}}{2}&0\\\
0&\omega-\frac{J_{0}\omega_{0}}{2}\end{array}\right),$ (51)
$h^{(1)}=\langle\Phi^{(0)}|V^{\prime}|\Phi^{(0)}\rangle=\left(\begin{array}[]{
c c}0&-\frac{J_{1}\omega_{0}}{2}\\\
-\frac{J_{1}\omega_{0}}{2}&0\end{array}\right),$ (52) $\displaystyle h^{(2)}$
$\displaystyle=$
$\displaystyle\langle\Phi^{(0)}|V^{\prime}|\Phi^{(1)}\rangle-h^{(1)}\langle\Phi^{(0)}|\Phi^{(1)}\rangle$
(55) $\displaystyle=$ $\displaystyle\left(\begin{array}[]{ c
c}\delta_{1^{\prime}}&\delta_{1^{\prime}0^{\prime}}\\\
\delta_{0^{\prime}1^{\prime}}&\delta_{0^{\prime}}\end{array}\right),$

with $\delta_{1^{\prime}}$, $\delta_{0^{\prime}}$,
$\delta_{1^{\prime}0^{\prime}}$ and $\delta_{0^{\prime}1^{\prime}}$ given in
Eqs. (31)-(34). Equations (51) and (52) form a first-order perturbation
matrix:

$H_{\text{GRWA}}=\left(\begin{array}[]{ c
c}\frac{J_{0}\omega_{0}}{2}&-\frac{J_{1}\omega_{0}}{2}\vspace{2mm}\\\
-\frac{J_{1}\omega_{0}}{2}&\omega-\frac{J_{0}\omega_{0}}{2}\end{array}\right),$
(56)

which is similar to the rotating wave approximation and neglects all other
non-resonant coupling terms. Note that this RWA is different from the
conventional one subject to the transverse coupling, where RWA breaks down in
the strong field, and we refer it to generalized rotating wave
approximation(GRWA). Diagonalize Eq. (56) to obtain the eigenvalues, and the
Rabi frequency is the difference between the two eigenvalues:

$\Omega^{\prime}=\sqrt{(\omega-J_{0}\omega_{0})^{2}+J_{1}^{2}\omega_{0}^{2}}.$
(57)

Then all frequencies of the Floquet dynamics can be written as
$\pm\Omega^{\prime}+2n\omega$ and $2n\omega$.

Equations (51)-(55) form the second-order perturbation matrix shown in Eq.
(30). In Fig. 6, we plot $\delta_{1^{\prime}}$, $\delta_{0^{\prime}}$,
$\delta_{1^{\prime}0^{\prime}}$ and $\delta_{0^{\prime}1^{\prime}}$ as a
function of $A$ for resonant and off-resonant situations. We find that
$\delta_{0^{\prime}}=-\delta_{1^{\prime}}$,
$\delta_{1^{\prime}0^{\prime}}=\delta_{0^{\prime}1^{\prime}}$ with
$\omega/\omega_{0}=1$ and $\omega/\omega_{0}=0.6$. In Fig. 7, we compare the
numerical and analytic results of the frequencies obtained by GVV Hamiltonian
and GRWA Hamiltonian. Clearly, the GVV results fit better than the GRWA to the
exact numerical results, indicating the deviation of the GRWA and the validity
of the GVV.

Figure 6:  Comparison of $\delta_{1^{\prime}}$ and $\delta_{0^{\prime}}$,
$\delta_{1^{\prime}0^{\prime}}$ and $\delta_{0^{\prime}1^{\prime}}$ for
various $A$ with $\omega/\omega_{0}=1$(a),(b) and
$\omega/\omega_{0}=0.6$(c),(d). Figure 7:  Comparison of numerical results
obtained by Floquet theory and approximate analytical GRWA results [i.e.,
$\pm\Omega^{\prime}+2n\omega$ and $2n\omega$ with $\Omega$ in Eq. (57)] and
GVV results [i.e., $\pm\Omega+2n\omega$ and $2n\omega$ with $\Omega$ in Eq.
(35)] results of the frequency components for various $A$ with
$\omega/\omega_{0}=0.6$.

##  Appendix C Numerical results obtained by Floquet theory

According to the Floquet theory, the Rabi model in Eq. (7) can be rewritten as

$H=H^{[0]}+H^{[-1]}e^{i\omega t}+H^{[1]}e^{-i\omega t}$ (58)

with $H^{[0]}=\omega_{0}\sigma_{z}/2$ and $H^{[\pm 1]}=A\sigma_{x}/4$. Based
on Eq. (3), the Floquet matrix is

$H^{\prime\prime}_{F}=\left(\begin{array}[]{c c c|c c|c c|c c|c c
c}\ddots&&&&&&&&&&&\iddots\\\
\hline\cr&\frac{\omega_{0}}{2}-2\omega&0&0&\frac{A}{4}&0&0&0&0&0&0&\\\
&0&\frac{-\omega_{0}}{2}-2\omega&\frac{A}{4}&0&0&0&0&0&0&0&\\\
\hline\cr&0&\frac{A}{4}&\frac{\omega_{0}}{2}-\omega&0&0&\frac{A}{4}&0&0&0&0&\\\
&\frac{A}{4}&0&0&\frac{-\omega_{0}}{2}-\omega&\frac{A}{4}&0&0&0&0&0&\\\
\hline\cr&0&0&0&\frac{A}{4}&\frac{\omega_{0}}{2}&0&0&\frac{A}{4}&0&0&\\\
&0&0&\frac{A}{4}&0&0&\frac{-\omega_{0}}{2}&\frac{A}{4}&0&0&0&\\\
\hline\cr&0&0&0&0&0&\frac{A}{4}&\frac{\omega_{0}}{2}+\omega&0&0&\frac{A}{4}&\\\
&0&0&0&0&\frac{A}{4}&0&0&\frac{-\omega_{0}}{2}+\omega&\frac{A}{4}&0&\\\
\hline\cr&0&0&0&0&0&0&0&\frac{A}{4}&\frac{\omega_{0}}{2}+2\omega&0&\\\
&0&0&0&0&0&0&\frac{A}{4}&0&0&\frac{-\omega_{0}}{2}+2\omega&\\\
\hline\cr\iddots&&&&&&&&&&&\ddots\end{array}\right).\begin{array}[]{c c}\\\
\leftarrow&|1,-2\rangle\\\ \leftarrow&|0,-2\rangle\\\
\leftarrow&|1,-1\rangle\\\ \leftarrow&|0,-1\rangle\\\
\leftarrow&|1,0\rangle\\\ \leftarrow&|0,0\rangle\\\ \leftarrow&|1,+1\rangle\\\
\leftarrow&|0,+1\rangle\\\ \leftarrow&|1,+2\rangle\\\
\leftarrow&|0,+2\rangle\\\ \end{array}$ (59)

We obtain quasienergies by numerically diagonalize $H^{\prime\prime}_{F}$. In
our calculations, we truncate the matrix $H^{\prime\prime}_{F}$ with n ranging
from -30 to 30.

##  Appendix D Unitary transformation of decay parameters

In this section, we give the derivation of Eqs. (37)-(40). In the lab frame,
the master equation of the dissipative open Rabi model with population damping
rate $\Gamma_{10}$ and dephasing rate $\gamma_{11}$ is [43]

$\dot{\rho}=-i[H(t),\rho]+\frac{\Gamma_{10}}{2}\mathcal{D}[\sigma_{01}]\rho+\gamma_{11}\mathcal{D}[\sigma_{11}]\rho.$
(60)

According to Eq. (III.2), the rotation matrix is

$U_{4}(t)=\cos\left[\frac{A\sin(\omega
t)}{2\omega}\right]I+i\sin\left[\frac{A\sin(\omega
t)}{2\omega}\right]\sigma_{x},$ (61)

with $I$ a $2\times 2$ identity matrix. The decay and excitation rate between
$|s^{\prime}_{1}(t)\rangle$ and $|s^{\prime}_{0}(t)\rangle$ can be easily
calculated by rewriting the lab frame dissipators in the new basis [41],

$\displaystyle\frac{\Gamma_{10}}{2}\mathcal{D}[\sigma_{01}]\rho$
$\displaystyle=$
$\displaystyle\frac{\Gamma_{10}}{2}\mathcal{D}[U_{4}^{\dagger}(t)\sigma_{s^{\prime}_{0}s^{\prime}_{1}}U_{4}(t)]\rho$
$\displaystyle=$
$\displaystyle\frac{\Gamma_{10}}{2}\mathcal{D}\left[\zeta(\sigma_{s^{\prime}_{0}s^{\prime}_{0}}\\!-\\!\sigma_{s^{\prime}_{1}s^{\prime}_{1}})\\!+\\!\beta\sigma_{s^{\prime}_{1}s^{\prime}_{0}}\\!+\\!\eta\sigma_{s^{\prime}_{0}s^{\prime}_{1}}\right]\rho,$
$\displaystyle\gamma_{11}\mathcal{D}[\sigma_{11}]\rho$ $\displaystyle=$
$\displaystyle\gamma_{11}\mathcal{D}[U_{4}^{\dagger}(t)\sigma_{s^{\prime}_{1}s^{\prime}_{1}}U_{4}(t)]\rho$
$\displaystyle=$
$\displaystyle\gamma_{11}\mathcal{D}\left[\beta\sigma_{s^{\prime}_{0}s^{\prime}_{0}}\\!+\\!\eta\sigma_{s^{\prime}_{1}s^{\prime}_{1}}\\!+\\!\zeta(\sigma_{s^{\prime}_{1}s^{\prime}_{0}}\\!-\\!\sigma_{s^{\prime}_{0}s^{\prime}_{1}})\right]\rho$

with $\zeta=\frac{1}{2}\sin\left[\frac{A\xi\sin(\omega t)}{2\omega}\right]$,
$\beta=\sin^{2}\left[\frac{A\xi\sin(\omega t)}{2\omega}\right]$ and
$\eta=\cos^{2}\left[\frac{A\xi\sin(\omega t)}{2\omega}\right]$. Therefore, by
regrouping the above dissipators and dropping out the fast oscillating terms,
such as
$\sigma_{s^{\prime}_{0}s^{\prime}_{1}}\rho\sigma_{s^{\prime}_{0}s^{\prime}_{1}}$
and
$\sigma_{s^{\prime}_{1}s^{\prime}_{0}}\rho\sigma_{s^{\prime}_{1}s^{\prime}_{0}}$
etc., we obtain the effective decay rates in Eqs. (37)-(40).

## References

  * Else et al. [2016] D. V. Else, B. Bauer, and C. Nayak, Phys. Rev. Lett. 117, 090402 (2016).
  * Sacha and Zakrzewski [2018] K. Sacha and J. Zakrzewski, Rept. Prog. Phys. 81, 016401 (2018).
  * Sondhi [2020] F. H. R. S. R. Sondhi, Annu. Rev. Condens. Matter Phys. 11 (2020).
  * Shevchenko et al. [2010] S. Shevchenko, S. Ashhab, and F. Nori, Phys. Rep. 492, 1 (2010).
  * Han et al. [2019] Y. Han, X.-Q. Luo, T.-F. Li, W. Zhang, S.-P. Wang, J. S. Tsai, F. Nori, and J. Q. You, Phys. Rev. Applied 11, 014053 (2019).
  * Han et al. [2022a] Y. Han, W. Zhang, and W. Li, Opt. Express 30, 7987 (2022a).
  * Eckardt [2017] A. Eckardt, Rev. Mod. Phys. 89, 011004 (2017).
  * Luo et al. [2008] X. Luo, Q. Xie, and B. Wu, Phys. Rev. A 77, 053601 (2008).
  * Jing et al. [2015] J. Jing, L.-A. Wu, M. Byrd, J. Q. You, T. Yu, and Z.-M. Wang, Phys. Rev. Lett. 114, 190502 (2015).
  * Zhang et al. [2007] W. Zhang, V. V. Dobrovitski, L. F. Santos, L. Viola, and B. N. Harmon, Phys. Rev. B 75, 201302(R) (2007).
  * Zhang et al. [2008] W. Zhang, N. P. Konstantinidis, V. V. Dobrovitski, B. N. Harmon, L. F. Santos, and L. Viola, Phys. Rev. B 77, 125336 (2008).
  * Scully and Zubairy [1997] M. Scully and M. Zubairy, _Quantum Optics_ (Cambridge University Press, 1997).
  * Shore and Knight [1993] B. W. Shore and P. L. Knight, J. Mod. Opt. 40, 1195 (1993).
  * Grifoni and Hanggi [1998] M. Grifoni and P. Hanggi, Phys. Rep. 304, 229 (1998).
  * Eberly et al. [1980] J. H. Eberly, N. B. Narozhny, and J. J. Sanchez-Mondragon, Phys. Rev. Lett. 44, 1323 (1980).
  * Casanova et al. [2010] J. Casanova, G. Romero, I. Lizuain, J. J. García-Ripoll, and E. Solano, Phys. Rev. Lett. 105, 263603 (2010).
  * Agarwal et al. [2017] K. Agarwal, S. Ganeshan, and R. N. Bhatt, Phys. Rev. B 96, 014201 (2017).
  * Gagnon et al. [2017] D. Gagnon, F. Fillion-Gourdeau, J. Dumont, C. Lefebvre, and S. MacLean, Phys. Rev. Lett. 119, 053203 (2017).
  * Nori [2009] F. Nori, Science 325, 689 (2009).
  * Forn-Díaz et al. [2019] P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, and E. Solano, Rev. Mod. Phys. 91, 025005 (2019).
  * Ashhab et al. [2007] S. Ashhab, J. R. Johansson, A. M. Zagoskin, and F. Nori, Phys. Rev. A 75, 063414 (2007).
  * Barata and Wreszinski [2000] J. C. A. Barata and W. F. Wreszinski, Phys. Rev. Lett. 84, 2112 (2000).
  * Irish et al. [2005] E. K. Irish, J. Gea-Banacloche, I. Martin, and K. C. Schwab, Phys. Rev. B 72, 195410 (2005).
  * Xie and Hai [2010] Q. Xie and W. Hai, Phys. Rev. A 82, 032117 (2010).
  * Wu and Yang [2007] Y. Wu and X. Yang, Phys. Rev. Lett. 98, 013601 (2007).
  * Deng et al. [2015] C. Deng, J.-L. Orgiazzi, F. Shen, S. Ashhab, and A. Lupascu, Phys. Rev. Lett. 115, 133601 (2015).
  * Lü and Zheng [2012] Z. Lü and H. Zheng, Phys. Rev. A 86, 023831 (2012).
  * Merlin [2021] R. Merlin, Am. J Phys. 89, 26 (2021).
  * Chu and Telnov [2004] S. I. Chu and D. A. Telnov, Phys. Rep. 390, 1 (2004).
  * Kohler et al. [2005] S. Kohler, J. Lehmann, and P. Hanggi, Phys. Rep. 406, 379 (2005).
  * Shirley [1965] J. H. Shirley, Phys. Rev. 138, 979 (1965).
  * Martin et al. [2017] I. Martin, G. Refael, and B. Halperin, Phys. Rev. X 7, 041008 (2017).
  * Yan et al. [2015] Y. Yan, Z. Lü, and H. Zheng, Phys. Rev. A 91, 053834 (2015).
  * Reiss [1970a] H. R. Reiss, Phys. Rev. A 1, 803 (1970a).
  * Reiss [1970b] H. R. Reiss, Phys. Rev. Lett. 25, 1149 (1970b).
  * Ran et al. [2020] D. Ran, B. Zhang, Y. H. Chen, Z. C. Shi, and A. Gover, Opt. Lett. 45 (2020).
  * Han et al. [2020] Y. Han, X.-Q. Luo, T.-F. Li, and W. Zhang, Phys. Rev. A 101, 022108 (2020).
  * Son et al. [2009] S.-K. Son, S. Han, and S.-I. Chu, Phys. Rev. A 79, 032301 (2009).
  * Schnell et al. [2021] A. Schnell, S. Denisov, and A. Eckardt, Phys. Rev. B 104, 165414 (2021).
  * Schnell et al. [2020] A. Schnell, A. Eckardt, and S. Denisov, Phys. Rev. B 101, 100301(R) (2020).
  * Lu et al. [2017] Y. Lu, S. Chakram, N. Leung, N. Earnest, R. K. Naik, Z. Huang, P. Groszkowski, E. Kapit, J. Koch, and D. I. Schuster, Phys. Rev. Lett. 119, 150502 (2017).
  * Ho and Chu [1985] T.-S. Ho and S.-I. Chu, Phys. Rev. A 32, 377 (1985).
  * Han et al. [2022b] Y. Han, M. Qiao, X.-Q. Luo, T.-F. Li, W. Zhang, X.-H. Deng, J. Q. You, and D. Yu, Phys. Rev. A 105, 063724 (2022b).