# Role of generalized parity in the symmetry of fluorescence spectrum from
two-level systems under periodic frequency modulation

Yiying Yan<EMAIL_ADDRESS>Department of Physics, School of Science,
Zhejiang University of Science and Technology, Hangzhou 310023, China  Zhiguo
Lü<EMAIL_ADDRESS>Key Laboratory of Artificial Structures and Quantum
Control (Ministry of Education), Department of Physics and Astronomy, Shanghai
Jiao Tong University, Shanghai 200240, China  Collaborative Innovation Center
of Advanced Microstructures, Nanjing 210093, China  JunYan Luo  Department of
Physics, School of Science, Zhejiang University of Science and Technology,
Hangzhou 310023, China  Hang Zheng<EMAIL_ADDRESS>Key Laboratory of
Artificial Structures and Quantum Control (Ministry of Education), Department
of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai 200240,
China  Collaborative Innovation Center of Advanced Microstructures, Nanjing
210093, China

###### Abstract

We study the origin of the symmetry of the fluorescence spectrum from the two-
level system subjected to a low-frequency periodic modulation and a near-
resonant high-frequency monochromatic excitation by using the analytical and
numerical methods based on the Floquet theory. We find that the fundamental
origin of symmetry of the spectrum can be attributed to the presence of the
generalized parity of the Floquet states, which depends on the driving
parameters. The absence of the generalized parity can lead to the asymmetry of
the spectrum. Based on the generalized parity, the conditions for the symmetry
and asymmetry of the spectrum can be derived, which succeeds in predicting
symmetry and asymmetry of the spectrum for the harmonic, biharmonic, and
multiharmonic modulations. Moreover, we find that the secular approximation
widely used in the analytical calculation may lead to artifact symmetry of the
spectrum that vanishes when such approximation is avoided. The present study
provides a significant perspective on the origin of the symmetry of the
spectrum.

##  I Introduction

Resonance fluorescence, arising from a quantum emitter driven by an external
field and coupled to a radiative reservoir Mollow (1969); Scully and Zubairy
(1997); Cohen-Tannoudji  _et al._ (1998), is not only an important concept in
quantum optics but also has potential application in quantum information
technology, for instance, it plays an important role in realizing the single-
photon source He  _et al._ (2013); Santana  _et al._ (2017); Kiršanskė  _et
al._ (2017). Particularly, the resonance fluorescence of two-level systems has
attracted much interest and been studied in various aspects such as spectrum
Ficek and Freedhoff (1993); Agarwal  _et al._ (1991); Ficek and Freedhoff
(1996); Ficek and Rudolph (1999); Peiris  _et al._ (2014); Konthasinghe  _et
al._ (2014); He  _et al._ (2015); Toyli  _et al._ (2016), squeezing Carmichael
(1985); Grünwald and Vogel (2012, 2013), photon statistics Kimble  _et al._
(1977); D’Souza  _et al._ (1990); Nazir (2008); Pastukhov  _et al._ (2014),
photon antibunching Itano  _et al._ (1988); Ficek  _et al._ (1984); Damanet
_et al._ (2018), and so on. The line shape of the spectrum is found to depend
strongly on the external field that interacts with the quantum emitters as
well as the reservoirs to which the quantum emitters are coupled. As is well-
known, for a sufficiently strong monochromatic field, the spectrum has a
symmetric three-peak structure, known as the Mollow triplet Mollow (1969).
More recently, the bi- and multi-chromatically driven quantum systems are of
interest Kryuchkyan  _et al._ (2017); Antón  _et al._ (2017); Yan  _et al._
(2018); Saiko  _et al._ (2018). In such systems, the spectrum turns out to
have a complicated multipeak structure Ficek and Freedhoff (1993); Agarwal
_et al._ (1991); Ficek and Freedhoff (1996); Ficek and Rudolph (1999); Peiris
_et al._ (2014); Konthasinghe  _et al._ (2014); He  _et al._ (2015), which can
be either symmetric or asymmetric. In principle, the physical origin of the
triplet and multipeak structures can be understood in terms of the transitions
between the quantum dressed states Cohen-Tannoudji  _et al._ (1998) or in
terms of the transitions between the semiclassical Floquet states Breuer and
Petruccione (1997); Yan  _et al._ (2016a). The studies on the resonance
fluorescence have enriched the physics concerning the light-matter
interaction.

The origin of the symmetry of the spectrum has been found in the case of the
monochromatic field. Specifically, it is the detailed balance condition that
guarantees the symmetry of the Mollow triplet Cohen-Tannoudji  _et al._
(1998). As is well-known, the breakdown of such a condition leads to the
asymmetry of the spectrum, for instance, in the presence of a pure dephasing
reservoir Roy and Hughes (2012); McCutcheon and Nazir (2013) or the counter-
rotating terms of the external field under certain conditions Browne and
Keitel (2000); Yan  _et al._ (2013, 2016a). The dephasing-induced asymmetric
Mollow triplet has been experimentally observed in the quantum dots (the pure
dephasing arises because of the interaction between the quantum dot and its
solid-state environment) Ulrich  _et al._ (2011); Ulhaq  _et al._ (2013). For
the bi- and multi-chromatic fields, the origin of the symmetry of the spectrum
is rarely discussed, owing to the fact that the physically transparent
spectrum is hardly analytically derived, and has not been comprehensively
understood.

Recent studies show that the fluorescence spectrum from a driven two-level
system with a modulated transition frequency is symmetrically multipeaked for
the vanishing detuning while asymmetrically multipeaked for the finite
detuning Yan  _et al._ (2016b); Kryuchkyan  _et al._ (2017); Antón  _et al._
(2017); Yan  _et al._ (2018). Such an exotic bichromatically driven two-level
system with coexistence of the longitudinal and transversal coupling between
the system and the applied fields has been experimentally studied in the
superconducting qubits Li  _et al._ (2013); Pan  _et al._ (2017), single
molecule Brunel  _et al._ (1998), and nitrogen-vacancy spin qubits Rohr  _et
al._ (2014). The quantum systems under frequency modulation are also of
interest in theoretical studies Kibis  _et al._ (2009); Macovei and Keitel
(2014); Zhao  _et al._ (2015); Silveri  _et al._ (2013); Macovei  _et al._
(2015), the intriguing phenomena of which were reviewed recently Silveri  _et
al._ (2017). It is worthwhile to note that the bichromatically driven two-
level system with frequency modulation differs from those considered in Refs.
Agarwal  _et al._ (1991); Ficek and Freedhoff (1993), where the two-level
systems are transversely driven by a bichromatic field. In such a case, the
symmetry of the fluorescence spectrum is found to depend on the average
detuning if the strengths of the two components of the bichromatic field are
the same; the pronounced asymmetry of the spectrum is revealed when the
average detuning is finite and/or the strengths of the two components of the
field are unequal Agarwal  _et al._ (1991); Ficek and Freedhoff (1993). For a
bichromatically amplitude-modulated field, the spectrum is also found to be
symmetric and asymmetric for the vanishing and finite detuning, respectively
Wilkens and Rza¸ewski (1989). So far the fundamental origin of such a
detuning-dependent symmetry remains obscure.

In this work, we use both analytical and numerical methods based on the
Floquet theory to study the fundamental origin of the symmetry of the
fluorescence spectrum from the two-level system under a low-frequency periodic
modulation and a near-resonant monochromatic excitation. We address the
symmetry and asymmetry of the spectrum by considering the generalized parity
of Floquet states rather than the behaviors of the bare-state or dressed-state
populations as considered in Refs. Das and Macovei (2013); Macovei  _et al._
(2015); Antón  _et al._ (2017). The generalized parity is found to guarantee
the symmetry of the spectrum while the breaking of such a parity can yield
pronouncedly asymmetric spectrum even in the vanishing detuning case. Based on
the generalized parity, the conditions for the symmetric and asymmetric
spectra are derived, which are not given in the previous works and cannot be
derived from the behaviors of the bare or dressed state population. The
generalized-parity-induced symmetry of the spectrum is verified and
illustrated in the context of the biharmonic modulation by the comparison
between the analytical and numerical results. The analytical results are found
to be in agreement with the numerically exact results in the regimes where the
perturbation theory and secular approximation can be justified. In addition,
we find that the spectrum with the secular approximation may have artifact
symmetry under certain conditions, i.e., the spectrum with secular
approximation is symmetric while the numerically exact calculation shows
asymmetric spectra because of the broken parity. The present finding simply
interprets the detuning-dependent symmetry in the harmonic modulation case and
can also be extended to analyze the symmetry and asymmetry of the spectrum in
the multiharmonic modulation cases. Our results suggest that it is feasible to
control the symmetry and asymmetry of the spectrum via engineering the
generalized parity of the Floquet states.

The rest of the paper is organized as follows. In Sec. II, we first discuss
the generalized-parity-induced symmetry of the fluorescence spectrum without
the secular approximation and further elucidate the symmetry of the spectrum
with a physically transparent formal spectrum with the secular approximation.
In Sec. III, we analytically and numerically calculate the fluorescence
spectrum in the context of the biharmonic modulation to verify the symmetry
and asymmetry of the spectrum predicted based on the generalized parity. In
the last section, the conclusions are given.

##  II Fluorescence spectrum and generalized parity

We consider that the transition frequency of the two-level system is modulated
periodically via a low-frequency external field $f(t)$ and the two-level
system is also excited by a near-resonant monochromatic field, which is
described by the following Hamiltonian ($\hbar=1$)

$H(t)=\frac{1}{2}[\omega_{0}+f(t)]\sigma_{z}+\frac{\Omega_{x}}{2}(\sigma_{+}e^{-i\omega_{x}t}+\sigma_{-}e^{i\omega_{x}t}),$
(1)

where $\sigma_{z(x,y)}$ is the usual Pauli matrix, $\omega_{0}+f(t)$ is the
modulated transition frequency, $\sigma_{\pm}=(\sigma_{x}\pm i\sigma_{y})/2$
are the raising and lowering operators, and $\Omega_{x}$ ($\omega_{x}$) is the
strength (frequency) of the monochromatic driving. Here we choose
$f(t)=f(t+T)$ with $T$ being the fundamental period of the modulation and much
greater than $2\pi/\omega_{x}$. This is a generalized model as compared with
the previous one considered in Refs. Yan  _et al._ (2016b); Kryuchkyan  _et
al._ (2017); Antón  _et al._ (2017).

To study the emission processes, we need to take account of the spontaneous
decay. Thus, the time evolution of the driven two-level system under study is
modeled by the Lindblad master equation. In the frame rotating at the
frequency $\omega_{x}$, the Lindblad master equation takes the form

$\frac{d}{dt}\tilde{\rho}(t)={\cal L}(t)\tilde{\rho}(t),$ (2)

where $\tilde{\rho}(t)$ is the reduced density matrix in the rotating frame
and the superoperator ${\cal L}(t)$ is given by ${\cal
L}(t)\tilde{\rho}(t)=-i[\tilde{H}(t),\tilde{\rho}(t)]-\kappa/2[\\{\sigma_{+}\sigma_{-},\tilde{\rho}(t)\\}-2\sigma_{-}\tilde{\rho}(t)\sigma_{+}]$
with $\kappa$ being the radiative decay rate. $\tilde{H}(t)$ is the effective
Hamiltonian and reads

$\tilde{H}(t)=\frac{\Omega_{x}}{2}\sigma_{x}+\frac{1}{2}[\delta+f(t)]\sigma_{z},$
(3)

with $\delta=\omega_{0}-\omega_{x}$ being the detuning between the bare
transition frequency and monochromatic excitation frequency. This master
equation is actually a set of first-order differential equations with periodic
coefficients. It can be directly solved by the so-called Floquet-Liouville
(FL) approach with a desire accuracy Ho  _et al._ (1986); Yan  _et al._
(2016b). Although such a Floquet-theory-based numerical method is simple and
efficient, it is not physically transparent to analyze the role of generalized
parity of Floquet states in the symmetry of the fluorescence spectrum. We use
an alternative method which is developed in our previous works Yan  _et al._
(2016a, 2018) to solve the master equation and calculate the fluorescence
spectrum. We first calculate the Floquet states for $\tilde{H}(t)$ and use
them as the bases to reformulate Eq. (2) and derive its analytical formal
solutions with the aid of the secular approximation in the Floquet picture.

###  II.1 The symmetry of fluorescence spectrum without secular approximation

The steady-state fluorescence spectrum is given by the Fourier transform of
the time-averaged first-order correlation function Mollow (1969); Ho  _et al._
(1986)

$S(\Delta)\propto{\rm
Re}\frac{1}{T}\int_{0}^{\infty}\int_{0}^{T}\lim_{t^{\prime}\rightarrow\infty}\left\langle\tilde{\sigma}_{+}(t^{\prime}+\tau)\tilde{\sigma}_{-}(t^{\prime})\right\rangle
e^{-i\Delta\tau}dt^{\prime}d\tau,$ (4)

where $\Delta=\omega-\omega_{x}$ and
$\left\langle\tilde{\sigma}_{+}(t^{\prime}+\tau)\tilde{\sigma}_{-}(t^{\prime})\right\rangle$
is the first-order correlation function and the tilde indicates that it is
evaluated in the rotating frame. In general, it is difficult to derive an
exact analytical spectrum. Nevertheless, we find that it is possible to show
that the spectrum is exactly symmetric about $\Delta=0$ when
$\delta+f(t)=-[\delta+f(t+T/2)]$ by realizing the fact that the driven two-
level system possesses a generalized parity symmetry, i.e.,

$\sigma_{x}\tilde{H}(t+T/2)\sigma_{x}=\tilde{H}(t).$ (5)

Here, the generalized parity transformation consists of an exchange between
the up and down states of two-level system
($\sigma_{z}\rightarrow-\sigma_{z}$) and a time shift of half period of the
modulation ($t\rightarrow t+T/2$).

We state briefly how the generalized parity guarantees the symmetry of the
spectrum. Owing to Eq. (5), we can construct a generalized parity
transformation in the Liouville space, the details of which can be found in
Appendix A. When $\delta+f(t)=-[\delta+f(t+T/2)]$, the superoperator ${\cal
L}(t)$ is similarly found to be invariant under the generalized parity
transformation. Based on this property, it can be derived from the master
equation (2) without the secular approximation that in the steady-state limit,
the time-averaged first-order correlation function is a real-valued function
in the rotating frame. As a result, the fluorescence spectrum is symmetric
about $\Delta=0$. This finding shows that the symmetry of the spectrum occurs
when $\delta+f(t)=-[\delta+f(t+T/2)]$ and results from the generalized parity.
We will numerically verify the generalized-parity-induced symmetry in Sec.
III.

###  II.2 The symmetry of fluorescence spectrum with secular approximation

To further elucidate the role of the generalized parity in determining the
symmetry of the spectrum, we calculate the spectrum in the Floquet picture
which allows us to derive a physically transparent formal spectrum with the
aid of the secular approximation.

According to the Floquet theory Shirley (1965); Sambe (1973), the time-
dependent Schrödinger equation governed by $\tilde{H}(t)$ possesses a set of
formal solutions
$|\tilde{\psi}_{\alpha}(t)\rangle=|\tilde{u}_{\alpha}(t)\rangle
e^{-i\tilde{\varepsilon}_{\alpha}t}$, where
$|\tilde{u}_{\alpha}(t)\rangle=|\tilde{u}_{\alpha}(t+T)\rangle$ is Floquet
state and $\tilde{\varepsilon}_{\alpha}$ is the corresponding real-valued
quasienergy. The index $\alpha$ labels independent Floquet states.
Substituting the formal solution into the Schrödinger equation, one readily
finds that

$[\tilde{H}(t)-i\partial_{t}]|\tilde{u}_{\alpha}(t)\rangle=\tilde{\varepsilon}_{\alpha}|\tilde{u}_{\alpha}(t)\rangle.$
(6)

On solving this equation, one obtains the Floquet states and quasienergies of
the driven two-level system.

We use $|\tilde{u}_{\alpha}(t)\rangle$ ($\alpha=\pm$) as the basis to
reformulate the master equation (2) and invoke the secular approximation Yan
_et al._ (2016a, 2018), yielding

$\displaystyle\frac{d}{dt}\tilde{\rho}_{++}(t)$ $\displaystyle=$
$\displaystyle-\Gamma_{{\rm rel}}\tilde{\rho}_{++}(t)+\Gamma_{{\rm s}},$ (7)
$\displaystyle\frac{d}{dt}\tilde{\rho}_{+-}(t)$ $\displaystyle=$
$\displaystyle-(i\Delta_{+-}+\Gamma_{{\rm deph}})\tilde{\rho}_{+-}(t),$ (8)

where
$\tilde{\rho}_{\alpha\beta}(t)=\langle\tilde{u}_{\alpha}(t)|\tilde{\rho}(t)|\tilde{u}_{\beta}(t)\rangle$
is the element of density operator,
$\Delta_{+-}=\tilde{\varepsilon}_{+}-\tilde{\varepsilon}_{-}$ is the
difference of two quasienergies, and $\Gamma_{{\rm
s}}=\kappa\sum_{l}|x_{-+,l}^{(+)}|^{2}$, where $x^{(+)}_{\alpha\beta,l}$ is a
time-averaged transition matrix element defined as follows:

$x^{(\pm)}_{\alpha\beta,l}=\frac{1}{T}\int^{T}_{0}\langle\tilde{u}_{\alpha}(t)|\sigma_{\pm}|\tilde{u}_{\beta}(t)\rangle
e^{-i2\pi lt/T}dt.$ (9)

The relaxation rate $\Gamma_{{\rm rel}}$ and dephasing rate $\Gamma_{{\rm
deph}}$ are given by

$\displaystyle\Gamma_{{\rm rel}}$ $\displaystyle=$
$\displaystyle\kappa\sum_{l}(|x_{+-,l}^{(+)}|^{2}+|x_{-+,l}^{(+)}|^{2}),$ (10)
$\displaystyle\Gamma_{{\rm deph}}$ $\displaystyle=$
$\displaystyle\frac{\kappa}{2}\sum_{l}(|x_{+-,l}^{(+)}|^{2}+|x_{-+,l}^{(+)}|^{2}+4|x_{++,l}^{(+)}|^{2}).$
(11)

The analytical formal solutions in the Floquet picture can be easily found as
follows:

$\displaystyle\tilde{\rho}_{++}(t)$ $\displaystyle=$
$\displaystyle\tilde{\rho}_{++}(0)e^{-\Gamma_{{\rm
rel}}t}+\tilde{\rho}_{++}^{{\rm ss}}(1-e^{-\Gamma_{{\rm rel}}t}),$ (12)
$\displaystyle\tilde{\rho}_{+-}(t)$ $\displaystyle=$
$\displaystyle\tilde{\rho}_{+-}(0)e^{-(\Gamma_{{\rm deph}}+i\Delta_{+-})t},$
(13)

where

$\tilde{\rho}_{++}^{{\rm ss}}=\frac{\Gamma_{{\rm s}}}{\Gamma_{{\rm
rel}}}=\frac{\sum_{l}|x_{-+,l}^{(+)}|^{2}}{\sum_{l}(|x_{+-,l}^{(+)}|^{2}+|x_{-+,l}^{(+)}|^{2})}$
(14)

is the steady-state population of the Floquet state. These solutions together
with the quantum regression theory enable us to derive a physically
transparent spectrum function Yan  _et al._ (2016a, 2018)

$\displaystyle S(\Delta)$ $\displaystyle\propto$
$\displaystyle\sum_{l}\bigg{\\{}\pi|x_{++,l}^{(+)}|^{2}(\tilde{\rho}_{++}^{{\rm
ss}}-\tilde{\rho}_{--}^{{\rm ss}})^{2}\delta(\Delta-l\omega_{z})$
$\displaystyle+4|x_{++,l}^{(+)}|^{2}\tilde{\rho}_{++}^{{\rm
ss}}\tilde{\rho}_{--}^{{\rm ss}}\frac{\Gamma_{{\rm rel}}}{\Gamma_{{\rm
rel}}^{2}+(\Delta-l\omega_{z})^{2}}$
$\displaystyle+|x_{+-,l}^{(+)}|^{2}\tilde{\rho}_{++}^{{\rm
ss}}\frac{\Gamma_{{\rm deph}}}{\Gamma_{{\rm
deph}}^{2}+(\Delta-l\omega_{z}-\Delta_{+-})^{2}}$
$\displaystyle+|x_{-+,l}^{(+)}|^{2}\tilde{\rho}_{--}^{{\rm
ss}}\frac{\Gamma_{{\rm deph}}}{\Gamma_{{\rm
deph}}^{2}+(\Delta-l\omega_{z}+\Delta_{+-})^{2}}\bigg{\\}},$

It is evident that the accuracy of Eq. (LABEL:eq:sfunsa) is limited by the
secular approximation when the transition matrix elements
$x^{(+)}_{\alpha\beta,l}$ and quasienergies are exactly calculated. As is
well-known, the secular approximation can be justified under the strong
driving condition, i.e., $\Delta_{+-}\gg\kappa$. In general, we can calculate
the quasienergies and transition matrix elements based on both analytical and
numerical diagonalization (ND) of the Floquet Hamiltonian
$\tilde{H}(t)-i\partial_{t}$ in the Sambe space Shirley (1965); Sambe (1973) ,
yielding the analytical and semianalytical spectra, respectively.

Next, we discuss the parity phenomenon of the Floquet states resulting from
Eq. (5). We consider the behavior of the Floquet states under the generalized
parity transformation ${\cal P}_{G}$, which is defined as

${\cal
P}_{G}|\tilde{u}_{\alpha}(t)\rangle:=\sigma_{x}|\tilde{u}_{\alpha}(t+T/2)\rangle.$
(16)

By differentiating $\sigma_{x}|\tilde{u}_{\alpha}(t+T/2)\rangle$ with respect
to $t$, we readily obtain

$\left[\sigma_{x}\tilde{H}\left(t+T/2\right)\sigma_{x}-i\partial_{t}\right]\sigma_{x}\left|\tilde{u}_{\alpha}\left(t+T/2\right)\right\rangle=\tilde{\varepsilon}_{\alpha}\sigma_{x}\left|\tilde{u}_{\alpha}\left(t+T/2\right)\right\rangle.$
(17)

When $\delta+f(t)=-[\delta+f(t+T/2)]$,
$\sigma_{x}|\tilde{u}_{\alpha}(t+T/2)\rangle$ satisfies the same differential
equation as $|\tilde{u}_{\alpha}(t)\rangle$ because of Eq. (5). Recalling the
uniqueness of solutions of the differential equations, in such cases we must
have

$\sigma_{x}\left|\tilde{u}_{\alpha}\left(t+T/2\right)\right\rangle=\lambda_{\alpha}|\tilde{u}_{\alpha}(t)\rangle,$
(18)

where $\lambda_{\alpha}$ is a constant. Furthermore, we have
$\lambda_{\alpha}=\pm 1$ because of ${\cal
P}_{G}^{2}|\tilde{u}_{\alpha}(t)\rangle=\lambda_{\alpha}^{2}|\tilde{u}_{\alpha}(t)\rangle=|\tilde{u}_{\alpha}(t)\rangle.$
Specifically, when $\delta+f(t)=-[\delta+f(t+T/2)]$, the Floquet states may be
even or odd functions under the generalized parity transformation, which is
referred to as the generalized parity of the Floquet states. The generalized
parity has been previously investigated in other phenomena such as the
coherent destruction of tunneling Grossmann  _et al._ (1991) and the laser-
induced electronic transport Lehmann  _et al._ (2003).

Clearly, if $\delta+f(t)\neq-\left[\delta+f\left(t+T/2\right)\right]$, Eq.
(18) cannot hold as
$\sigma_{x}\tilde{H}\left(t+T/2\right)\sigma_{x}\neq\tilde{H}(t)$, i.e., the
effective Hamiltonian is no longer invariant under the generalized parity
transformation. Consequently, the Floquet states also do not have the
generalized parity.

Figure 1: The incoherent components of the fluorescence spectrum for $p=3$,
$\Omega_{x}=10\kappa$, $\delta=0$, $\Omega_{z}=\omega_{z}=40\kappa$, $r=1$,
and various phase. “Ana.” and “Num.” denote the analytical and the FL
numerical results, respectively.

We show that the symmetry of the spectrum may be a consequence of the
generalized parity of the Floquet states. By using Eq. (18) and
$x_{\alpha\beta,l}^{(+)}=\left[x_{\beta\alpha,-l}^{(-)}\right]^{\ast}$, it is
straightforward to show the following identity for arbitrary integer $l$ from
the definition (9) of the transition matrix element:

$x_{\alpha\beta,l}^{(+)}=(-1)^{l}\lambda_{\alpha}\lambda_{\beta}\left[x_{\beta\alpha,-l}^{(+)}\right]^{\ast},$
(19)

provided $\delta+f(t)=-[\delta+f(t+T/2)]$. It follows that

$|x_{\alpha\beta,l}^{(+)}|=|x_{\beta\alpha,-l}^{(+)}|$ (20)

also holds for any integer $l$. We emphasize that the relation (20) can be
deduced from relation (19), however, the relation (19) cannot be derived from
relation (20). With the relation (20), it is straightforward to show that the
spectrum (LABEL:eq:sfunsa) is symmetric about $\Delta=0$ Yan  _et al._ (2018).
Specifically, since $|x_{++,l}^{(+)}|=|x_{++,-l}^{(+)}|$, the emission lines
at $\Delta=\pm l\omega_{z}$ (the positions are symmetric about $\Delta=0$)
have the equal weights. Moreover, since $|x_{+-,l}^{(+)}|=|x_{-+,-l}^{(+)}|$,
we also have $\tilde{\rho}_{++}^{{\rm ss}}=\tilde{\rho}_{--}^{{\rm ss}}$
according to Eq. (14), leading to $|x_{+-,l}^{(+)}|^{2}\tilde{\rho}_{++}^{{\rm
ss}}=|x_{-+,-l}^{(+)}|^{2}\tilde{\rho}_{--}^{{\rm ss}}$. That is to say, the
emission lines at $\Delta=\pm(l\omega_{z}+\Delta_{+-})$ (the positions are
symmetric about $\Delta=0$) have the same weights. It turns out that the
symmetry of the spectrum fundamentally originates from the generalized parity
of the Floquet states when $\delta+f(t)=-[\delta+f(t+T/2)]$. Conversely, one
may expect that the symmetry of the spectrum may break when such a parity is
absent. However, it is a formidable task to analytically prove that the
spectrum is asymmetric in the absence of the generalized parity.

Let us discuss what happens to the formal spectrum if
$\delta+f(t)\neq-[\delta+f(t+T/2)]$. Under such a condition, the generalized
parity is absent, and thus we cannot have the relation (19). In principle, the
absence of the generalized parity will result in two possible situations. One
is that the spectrum becomes asymmetric about $\Delta=0$ because the relation
$|x^{(+)}_{\alpha\beta,l}|\neq|x^{(+)}_{\beta\alpha,-l}|$ can be derived at
least for a certain $l$. The other is that the spectrum is symmetric because
the equality $|x^{(+)}_{\alpha\beta,l}|=|x^{(+)}_{\beta\alpha,-l}|$ still
holds for any $l$, originating from other kinds of identities between the
transition matrix elements rather than the generalized-parity-induced identity
(19). Apparently the first situation is more trivial than the second one. Most
importantly, the present analysis suggests that the formal spectrum may be
symmetric even without the generalized parity. Consequently, we cannot
conclude from the formal spectrum (LABEL:eq:sfunsa) that the symmetry of the
spectrum breaks as long as the generalized parity is absent.

To end this section, we give some remarks on the above findings based on the
formal spectrum. First, we find that the symmetry of the spectrum may result
from the generalized parity and requires $\delta+f(t)=-[\delta+f(t+T/2)]$.
This is consistent with the analysis above without the secular approximation.
Moreover, the generalized parity is found to be an important underlying cause
of the relation (20), which was numerically found in harmonic modulation case
Yan  _et al._ (2018). It turns out here that the relation (20) can be
established due to the generalized parity in the bi- and multi-harmonic cases.
Second, without the generalized parity, namely, when
$\delta+f(t)\neq-[\delta+f(t+T/2)]$, the formal spectrum can be either
trivially asymmetric or nontrivially symmetric. The symmetry requires the
relation (20) in the absence of the generalized parity, namely, Eq. (19).
Third, the formal spectrum is derived with the secular approximation and thus
the present analysis needs further verification. In what follows we consider a
concrete biharmonic modulation to verify whether the generalized parity
guarantees the symmetry of the spectrum when the secular approximation is not
invoked and we also check whether the relation (20) can be established without
the generalized parity and whether such relations lead to the symmetry of the
spectrum without the secular approximation.

Figure 2: The incoherent components of the fluorescence spectrum for $p=3$,
$\delta=5\kappa$, $\Omega_{x}=10\kappa$, $\Omega_{z}=\omega_{z}=40\kappa$,
$r=1$, and various phase.

##  III Verification of symmetry and asymmetry of the spectrum

To calculate fluorescence spectrum, without loss of generality, we mainly
consider the biharmonic modulation in this work, namely, the modulation
consists of two harmonics

$f(t)=\Omega_{z}[\cos(\omega_{z}t)+r\cos(p\omega_{z}t+\phi)],$ (21)

where $\Omega_{z}$ and $\omega_{z}=2\pi/T$ are the amplitude and fundamental
frequency of the modulation, respectively, $p$ is a positive integer, $r$ is
the ratio of the amplitude of the second harmonic to that of the first one,
and $\phi$ is a relative phase. Since $\frac{1}{T}\int^{T}_{0}f(t)dt=0$, the
condition for the presence of the generalized parity
$\delta+f(t)=-[\delta+f(t+T/2)]$ is equivalent to $\delta=0$ and
$f(t)=-f(t+T/2)$. The condition for the absence of the generalized parity
$\delta+f(t)\neq-[\delta+f(t+T/2)]$ is simply divided into three cases:

$\left\\{\begin{array}[]{c}\delta\neq 0\,{\rm and}\,f(t)=-f(t+T/2);\\\
\delta=0\,{\rm and}\,f(t)\neq-f(t+T/2);\\\ \delta\neq 0\,{\rm
and}\,f(t)\neq-f(t+T/2).\end{array}\right.$ (22)

It is noted that for the biharmonic modulation (21), both $f(t)=-f(t+T/2)$ and
$f(t)\neq-f(t+T/2)$ can be realized by setting $p$ odd and even numbers,
respectively. To verify above analysis, we calculate the numerically exact
fluorescence spectrum from master equation (2) with the FL formalism Ho  _et
al._ (1986); Yan  _et al._ (2016b), which is compared with the analytical and
semianalytical results from Eq. (LABEL:eq:sfunsa). The analytical and
semianalytical results are obtained by using the transition matrix elements
and quasienergies calculated with the Van Vleck perturbation theory and the ND
of the Floquet Hamiltonian, respectively. The detailed analytical calculation
is presented in Appendix B. In addition, we just focus on the incoherent
components of the fluorescence spectrum, which is of interest in the
experiments. In principle, similar analysis is applicable to the coherent
components. In this work, we mainly consider the parameters regime
$\omega_{z}\sim\Omega_{z}\gg\Omega_{x}\gg\kappa$, in which case both the Van
Vleck perturbation theory (up to second order in $\Omega_{x}$) and secular
approximation can be justified. Importantly, this regime is experimentally
accessible in the artificial atoms, e.g., the transmon qubit Li  _et al._
(2013). We should emphasize that if the perturbation theory is inapplicable,
we can obtain the transition matrix elements and quasienergies by the ND of
the Floquet Hamiltonian.

We first verify whether the generalized parity guarantees the symmetry of the
spectrum. In Fig. 1, we display the incoherent component of fluorescence
spectra obtained by the FL numerical method (solid line) and analytical result
(dashed line) for $p=3$, $\delta=0$, and various values of $\phi$. Apparently
the spectra are symmetric as expected. The analytical results are in agreement
with the FL results. These results also show that the spectrum depends weakly
on the relative phase $\phi$. In addition, it is straightforward to verify
that for other driving parameters, the spectrum is symmetric as well when $p$
is an odd number and $\delta=0$. In Appendix C, we show that when $\delta=0$
and $p$ is odd, the transition matrix elements indeed satisfy Eq. (19), which
guarantees the symmetry of the spectrum. The present results suggest that the
symmetry of the spectrum appears as long as $\delta=0$ and $f(t)=-f(t+T/2)$
and fundamentally originates from the generalized parity of the Floquet states
in such a situation.

Figure 3: The incoherent components of the fluorescence spectrum for $p=2$,
$\delta=0$, $\Omega_{x}=10\kappa$, $\Omega_{z}=\omega_{z}=40\kappa$, $r=1$,
and various phases.

We move to examine whether the symmetry of the spectrum breaks when the
generalized parity is absent, namely, under the conditions
$\delta+f(t)\neq-[\delta+f(t+T/2)]$. We calculate the spectra with the
parameters being the same as in Fig. 1 except for the detuning
$\delta=5\kappa$, corresponding to the case of $\delta\neq 0$ and
$f(t)=-f(T+T/2)$. In Fig. 2, the analytical and FL numerical spectra agree
with each other and are found to be asymmetric for the finite detuning,
indicating that in spite of $f(t)=-f(t+T/2)$, the asymmetry of spectrum
appears when $\delta\neq 0$.

Let us consider the case of $\delta=0$ and $f(t)\neq-f(t+T/2)$ by setting $p$
being even. We calculate the spectrum for $p=2$ and the other parameters being
the same as in Fig. 1. Figure 3 displays that the analytical and numerical
spectra are pronouncedly asymmetric even though $\delta=0$ except for
$\phi=\pi/2$ in which case the analytical spectrum is found to be strictly
symmetric (see discussion below) while the numerical spectrum is slightly
asymmetric [in particular, the intensities of emission lines at
$\Delta=\pm\omega_{z}$ are unequal as shown in Fig. 6(a)]. These results
confirm that the formal spectrum (LABEL:eq:sfunsa) may be symmetric without
the generalized parity of the Floquet states. However, the numerically exact
spectrum is asymmetric in the absence of the generalized parity. This shows
that the generalized parity plays an important role in determining the
symmetry of the exact spectrum. We will further analyze such discrepancy
between the analytical and numerical results later. In addition, we find that
in contrast with $p=3$, the spectrum is found to depend strongly on relative
phase $\phi$ when $p=2$.

Finally we calculate the spectra for $\delta\neq 0$ and $f(t)\neq-f(t+T/2)$.
Figure 4 shows the spectra obtained for the detuning $\delta=5\kappa$ and the
other parameters being the same as in Fig. 3. The spectra are still
asymmetric. In general, it is straightforward to verify the asymmetry of the
spectrum under the condition that $\delta+f(t)\neq-[\delta+f(t+T/2)]$. All in
all, it turns out that the symmetry of the spectrum breaks in the absence of
the generalized parity. Conversely, we can say that the symmetry of the
spectrum can be fully attributed to the presence of the generalized parity. In
contrast to the previous studies, we ascribe the asymmetry to the breaking of
the generalized parity rather than the unequal populations of dressed states
Antón  _et al._ (2017) or the breakdown of relation (20) Yan  _et al._ (2018).

Figure 4: The incoherent components of fluorescence spectrum for $p=2$,
$\delta=5\kappa$, $\Omega_{x}=10\kappa$, $\Omega_{z}=\omega_{z}=40\kappa$,
$r=1$, and various phase.

Let us explore how the analytical spectrum becomes symmetric in the absence of
the generalized parity of the Floquet states. To this end, we show that the
relation (20) can originate from the identities different from Eq. (19). Based
on the results from the Van Vleck perturbation theory, we analytically derive
the identities for the transition matrix elements in the case of vanishing
detuning and even $p$. The derivation are given in Appendix C. When $p$ is
even, $\delta=0$, and $\phi=\left(1/2+n\right)\pi$ ($n=0,\pm 1,\pm 2,\ldots$),
we find that the following relations hold for arbitrary integer $l$:

$\displaystyle x^{(+)}_{++,-l}$ $\displaystyle=$
$\displaystyle(-1)^{l}x^{(+)}_{++,l},$ (23) $\displaystyle x^{(+)}_{-+,-l}$
$\displaystyle=$ $\displaystyle-(-1)^{l}e^{-i2\theta_{0}}x^{(+)}_{+-,l},$ (24)

where $\theta_{0}$ is a phase defined in Eq. (86). Although the relations (23)
and (24) are derived based on the perturbation theory, it is straightforward
to show that they hold in the nonperturbative regimes. In Fig. 5, we calculate
$x^{(+)}_{++,l}$ $(l=\pm 1,\pm 2)$ with the variation of $\Omega_{x}$ by using
the analytical and ND methods. We see that the deviation between the
analytical and numerical results becomes larger and larger as $\Omega_{x}$
increases, which is due to the breakdown of the perturbation calculation.
Nevertheless, $x^{(+)}_{++,l}$ obtained by the ND method still satisfies Eq.
(23). This suggests that the relations (23) and (24) are not limited to the
perturbative regimes. More importantly, it follows from the identities (23)
and (24) that $|x_{\alpha\beta,l}^{(+)}|=|x_{\beta\alpha,-l}^{(+)}|$, which
leads to the symmetry of the formal spectrum (LABEL:eq:sfunsa). That is to
say, without the generalized parity of the Floquet states, the relation (20)
can also be established from other kinds of the identities for the transition
matrix elements instead of the generalized-parity-induced identity (19) under
certain conditions.

The discrepancy in the symmetry predicted by the analytical and numerical
methods shown in Fig. 3(b) indicates that the relations (23) and (24) cannot
guarantee the symmetry of the spectrum without the secular approximation. To
further verify this, in Fig. 6, we use semianalytical and FL numerical methods
to calculate the weights of the emission lines at $\Delta=\pm\omega_{z}$ as
the increasing of $\Omega_{x}$ for $p=2$, $\delta=0$, and two values of
$\phi$. It is evident that the weights calculated from the semianalytical
method (solid and dashed lines) are the same while the weights from the
numerical method (dot-dashed and dotted lines) are unequal, indicating that
the semianalytical spectrum is symmetric but the numerical spectrum is not
symmetric. The present results illustrate that that provided the relation (20)
is established in the absence of the generalized parity, the secular
approximation can induce artifact symmetry that vanishes if such approximation
is not invoked.

Figure 5: Transition matrix elements $x^{(+)}_{++,l}$ versus driving strength
$\Omega_{x}$, calculated from the analytical method and the numerical method
based on the ND of the Floquet Hamiltonian for $p=2$, $\delta=0$,
$\Omega_{z}=\omega_{z}=40\kappa$, $\phi=\pi/2$, and $r=1$.

Apart from the biharmonic modulation, we find that the conditions for the
symmetry and asymmetry of the spectrum, which are derived based on the
generalized parity, are applicable to the simple harmonic and multiharmonic
modulation cases. For the simple harmonic modulation
$f(t)=\Omega_{z}\cos(\omega_{z}t)$, $f(t)=-f(t+T/2)$ is met. Therefore, the
symmetry and asymmetry of the spectrum is uniquely controlled by the detuning
$\delta$, which simply interprets the detuning-dependent symmetry of the
spectrum. Specifically, the spectrum is expected to be symmetric when
$\delta=0$ and asymmetric when $\delta\neq 0$. This is consistent with the
findings of previous studies Yan  _et al._ (2016b); Antón  _et al._ (2017);
Yan  _et al._ (2018). For the multiharmonic modulation
$f(t)=\sum_{p=1}^{N}\Omega_{z,p}\cos(p\omega_{z}t+\phi_{p})]$, where
$\Omega_{z,p}$ and $\phi_{p}$ are the amplitude and phase of the $p$th
harmonic, respectively, either $f(t)=-f(t+T/2)$ or $f(t)\neq-f(t+T/2)$ can be
met, similarly to the biharmonic case. We have calculated the spectrum with
the FL and semianalytical methods for the cases of $N=3$, $N=4$, and $N=5$.
The results (not shown here) further confirm that the symmetry and asymmetry
of spectrum fundamentally originate from the presence and absence of the
generalized parity of the Floquet states, respectively.

Figure 6: Weights of emission lines at $\Delta=\pm\omega_{z}$ versus driving
strength $\Omega_{x}$, calculated from the semianalytical method and the FL
method, for $p=2$, $\delta=0$, $\Omega_{z}=\omega_{z}=40\kappa$, $r=1$, and
two values of $\phi$. “Semiana.” denotes the semianalytical result.

##  IV Conclusions

In summary, we have studied the fundamental origin of the symmetry of the
resonance fluorescence from the two-level system subjected to a periodic
frequency modulation and a near-resonant high-frequency monochromatic
excitation by using both analytical and numerical methods based on the Floquet
theory. In such a driven two-level system, we have found that the generalized
parity of Floquet states plays a fundamental role in the symmetry of the
spectrum. Specifically, the generalized parity guarantees the symmetry of the
spectrum. On the other hand, when the generalized parity is broken, the
spectrum becomes asymmetric. This has been illustrated in the context of the
biharmonic modulation, the parameters of which can be tuned to induce or break
the generalized parity. For the biharmonic modulation, we find that when
$\delta=0$ and $f(t)=-f(t+T/2)$, the generalized parity exists and the
spectrum is symmetric. When $\delta+f(t)\neq-[\delta+f(t+T/2)]$, the
generalized parity is broken and the spectrum is found to be asymmetric.
Interestingly, we can obtain pronouncedly asymmetric spectrum by requiring the
modulation $f(t)\neq-f(t+T/2)$ even though $\delta=0$. Moreover, these
conditions for the symmetry and asymmetry of the spectrum are found to be
applicable to the simple harmonic and multiharmonic modulation cases. In
addition, we illustrated that the secular approximation may induce artifact
symmetry that vanishes if the secular approximation is avoided under certain
conditions. The present study gives a deep insight into the origin of the
symmetry of the spectrum and reveals a simple relation between the symmetry of
the spectrum and the generalized parity of the Floquet states.

###### Acknowledgements.

This work was supported by the National Natural Science Foundation of China
(Grants No. 11647082, No. 11774311, No. 11774226, and No. 11874260).

##  Appendix A Derivation of symmetry of the spectrum without the secular
approximation

The master equation can be rewritten in a matrix form

$\frac{d}{dt}\vec{\tilde{\rho}}(t)={\cal L}(t)\vec{\tilde{\rho}}(t).$ (25)

Here the vector is defined as

$\vec{\tilde{\rho}}(t)=(\langle\tilde{\sigma}_{+}(t)\rangle,\langle\tilde{\sigma}_{-}(t)\rangle,\langle\tilde{\pi}_{+}(t)\rangle,\langle\tilde{\pi}_{-}(t)\rangle)^{{\rm
T}},$ (26)

where $\pi_{\pm}=(1\pm\sigma_{z})/2$ and
$\langle\tilde{\hat{o}}(t)\rangle\equiv{\rm Tr}[\hat{o}\tilde{\rho}(t)]$. The
superoperator ${\cal L}(t)$ in the Liouville space spanned by the matrix bases
$\\{\sigma_{\pm},\pi_{\pm}\\}$ is given by

${\cal
L}(t)=\left(\begin{array}[]{cccc}i[\delta+f(t)]-\frac{\kappa}{2}&0&-\frac{i\Omega_{x}}{2}&\frac{i\Omega_{x}}{2}\\\
0&-i[\delta+f(t)]-\frac{\kappa}{2}&\frac{i\Omega_{x}}{2}&\frac{-i\Omega_{x}}{2}\\\
\frac{-i\Omega_{x}}{2}&\frac{i\Omega_{x}}{2}&-\kappa&0\\\
\frac{i\Omega_{x}}{2}&\frac{-i\Omega_{x}}{2}&\kappa&0\end{array}\right).$ (27)

If $\delta+f(t)=-[\delta+f(t+T/2)]$, in which case the Hamiltonian is
invariant under the generalized parity transformation, one readily finds that

${\cal T}{\cal L}(t+T/2){\cal T}={\cal L}(t),$ (28)

where the transformation matrix is given by

${\cal T}=\left(\begin{array}[]{cccc}0&1&0&0\\\ 1&0&0&0\\\ 0&0&-1&0\\\
0&0&0&-1\end{array}\right),$ (29)

and ${\cal T}^{2}=I$ with $I$ being the identity matrix. Similarly to the
Hamiltonian, the matrix ${\cal L}(t)$ is invariant under the transformation
defined in Eq. (28), which can be regarded as the generalized parity
transformation in the Liouville space, similarly to that defined in Eq. (16)
of the main text.

Let us derive the specific property of the steady state in the long-time limit
[as $\det{\cal L}(t)=0$, there exists a nontrivial steady state]. It follows
from Eq. (25) that

$\frac{d}{dt}\vec{\tilde{\rho}}(t+T/2)={\cal
L}(t+T/2)\vec{\tilde{\rho}}(t+T/2),$ (30)

which leads to

$\displaystyle\frac{d}{dt}{\cal T}\vec{\tilde{\rho}}(t+T/2)$ $\displaystyle=$
$\displaystyle{\cal T}{\cal L}(t+T/2){\cal T}{\cal
T}\vec{\tilde{\rho}}(t+T/2)={\cal L}(t){\cal T}\vec{\tilde{\rho}}(t+T/2),$
(31)

which means that ${\cal T}\vec{\tilde{\rho}}(t+T/2)=c\vec{\tilde{\rho}}(t)$,
owing to the uniqueness of solutions of the differential equation. On using
the fact that $\vec{\tilde{\rho}}(t)=\vec{\tilde{\rho}}(t+T)$ as
$t\rightarrow\infty$ because of ${\cal L}(t)={\cal L}(t+T)$, we find that $c$
may be either $+1$ or $-1$. It is easy to prove by contradiction that $c=-1$.
Suppose that $c=1$, yielding
$\langle\tilde{\pi}_{+}(t+T/2)\rangle=-\langle\tilde{\pi}_{+}(t)\rangle$.
However, if one considers $\delta+f(t)=0$ in which case ${\cal L}(t)$ is time
independent while Eq. (28) still holds, the steady state becomes time
independent and one gets
$\langle\tilde{\pi}_{+}(t)\rangle=\langle\tilde{\pi}_{+}(t+T/2)\rangle$. By
contradiction, one finds that $c=-1$. Consequently, in the steady-state limit,
we have

${\cal
T}\vec{\tilde{\rho}}(t+T/2)=-\vec{\tilde{\rho}}(t)\quad(t\rightarrow\infty).$
(32)

Next, let us derive the property of the principal matrix solution
$\Pi(t,t^{\prime})$ of the master equation, which solves the differential
equation

$\frac{d}{dt}\Pi(t,t^{\prime})={\cal L}(t)\Pi(t,t^{\prime}),$ (33)

with the initial condition $\Pi(t^{\prime},t^{\prime})=I$. It is
straightforward to show that

$\displaystyle\frac{d}{dt}{\cal T}\Pi(t+T/2,t^{\prime}+T/2){\cal T}$
$\displaystyle=$ $\displaystyle{\cal T}{\cal L}(t+T/2){\cal T}{\cal
T}\Pi(t+T/2,t^{\prime}+T/2){\cal T}={\cal L}(t){\cal
T}\Pi(t+T/2,t^{\prime}+T/2){\cal T},$ (34)

namely, ${\cal T}\Pi(t+T/2,t^{\prime}+T/2){\cal T}$ satisfies the same
differential equation and the same initial condition as $\Pi(t,t^{\prime})$.
As a result, we simply have

${\cal T}\Pi(t+T/2,t^{\prime}+T/2){\cal T}=\Pi(t,t^{\prime}).$ (35)

According to the quantum regression theory Mollow (1969), the two-time
correlation functions

$\vec{\tilde{g}}(t,t^{\prime})=(\langle\tilde{\sigma}_{+}(t)\tilde{\sigma}_{-}(t^{\prime})\rangle,\langle\tilde{\sigma}_{-}(t)\tilde{\sigma}_{-}(t^{\prime})\rangle,\langle\tilde{\pi}_{+}(t)\tilde{\sigma}_{-}(t^{\prime})\rangle,\langle\tilde{\pi}_{-}(t)\tilde{\sigma}_{-}(t^{\prime})\rangle)^{{\rm
T}}$ (36)

satisfy the same equation as $\vec{\tilde{\rho}}(t)$, however, with a
different initial condition

$\vec{\tilde{g}}(t^{\prime},t^{\prime})=(\langle\tilde{\pi}_{+}(t^{\prime})\rangle,0,0,\langle\tilde{\sigma}_{-}(t^{\prime})\rangle)^{{\rm
T}}.$ (37)

Similarly, another set of two-time correlation functions

$\vec{\tilde{G}}(t,t^{\prime})=(\langle\tilde{\sigma}_{+}(t^{\prime})\tilde{\sigma}_{+}(t)\rangle,\langle\tilde{\sigma}_{+}(t^{\prime})\tilde{\sigma}_{-}(t)\rangle,\langle\tilde{\sigma}_{+}(t^{\prime})\tilde{\pi}_{+}(t)\rangle,\langle\tilde{\sigma}_{+}(t^{\prime})\tilde{\pi}_{-}(t)\rangle)^{{\rm
T}}$ (38)

also satisfy the same differential equation as $\vec{\tilde{g}}(t,t^{\prime})$
but with the initial condition

$\vec{\tilde{G}}(t^{\prime},t^{\prime})=(0,\langle\tilde{\pi}_{+}(t^{\prime})\rangle,0,\langle\tilde{\sigma}_{+}(t^{\prime})\rangle)^{{\rm
T}}.$ (39)

Using Eq. (32), we have

${\cal T}\vec{\tilde{g}}(t^{\prime},t^{\prime})=\left(\begin{array}[]{c}0\\\
\langle\tilde{\pi}_{+}(t^{\prime})\rangle\\\ 0\\\
-\langle\tilde{\sigma}_{-}(t^{\prime})\rangle\end{array}\right)=\left(\begin{array}[]{c}0\\\
\langle\tilde{\pi}_{+}(t^{\prime}+T/2)\rangle\\\ 0\\\
\langle\tilde{\sigma}_{+}(t^{\prime}+T/2)\rangle\end{array}\right)=\vec{\tilde{G}}\left(t^{\prime}+\frac{T}{2},t^{\prime}+\frac{T}{2}\right)\quad(t^{\prime}\rightarrow\infty).$
(40)

In the steady-state limit, the correlation functions are found to have the
following relation

$\displaystyle\vec{\tilde{g}}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\Pi(t,t^{\prime})\vec{\tilde{g}}(t^{\prime},t^{\prime})$ (41)
$\displaystyle=$ $\displaystyle{\cal
T}\Pi\left(t+\frac{T}{2},t^{\prime}+\frac{T}{2}\right){\cal
T}\vec{\tilde{g}}(t^{\prime},t^{\prime})$ $\displaystyle=$ $\displaystyle{\cal
T}\Pi\left(t+\frac{T}{2},t^{\prime}+\frac{T}{2}\right)\vec{\tilde{G}}\left(t^{\prime}+\frac{T}{2},t^{\prime}+\frac{T}{2}\right)$
$\displaystyle=$ $\displaystyle{\cal
T}\vec{\tilde{G}}\left(t+\frac{T}{2},t^{\prime}+\frac{T}{2}\right)\quad(t^{\prime}\rightarrow\infty).$

It follows that as $t^{\prime}\rightarrow\infty$,

$\displaystyle\langle\tilde{\sigma}_{+}(t)\tilde{\sigma}_{-}(t^{\prime})\rangle$
$\displaystyle=$
$\displaystyle\langle\tilde{\sigma}_{+}(t^{\prime}+T/2)\tilde{\sigma}_{-}(t+T/2)\rangle$
(42) $\displaystyle=$
$\displaystyle\langle\tilde{\sigma}_{+}(t+T/2)\tilde{\sigma}_{-}(t^{\prime}+T/2)\rangle^{\ast}.$

In the steady-state limit, the first-order correlation function depends
explicitly on time $t^{\prime}$, however, the $t^{\prime}$ dependence can be
eliminated by setting $t=\tau+t^{\prime}$ and integrating over $t^{\prime}$
(because the contributions of $t^{\prime}$-dependent terms are negligible to a
long-time observation), yielding the $\tau$-dependent first-order correlation
function

$\displaystyle\bar{\tilde{g}}_{1}(\tau)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{T}\int_{0}^{T}\lim_{t^{\prime}\rightarrow\infty}\langle\tilde{\sigma}_{+}(\tau+t^{\prime})\tilde{\sigma}_{-}(t^{\prime})\rangle
dt^{\prime}$ (43) $\displaystyle=$
$\displaystyle\frac{1}{T}\int_{0}^{T}\lim_{t^{\prime}\rightarrow\infty}\langle\tilde{\sigma}_{+}(\tau+t^{\prime}+T/2)\tilde{\sigma}_{-}(t^{\prime}+T/2)\rangle^{\ast}dt^{\prime}$
$\displaystyle=$
$\displaystyle\frac{1}{T}\int_{T/2}^{T+T/2}\lim_{t^{\prime}\rightarrow\infty}\langle\tilde{\sigma}_{+}(\tau+t^{\prime})\tilde{\sigma}_{-}(t^{\prime})\rangle^{\ast}dt^{\prime}$
$\displaystyle=$
$\displaystyle\frac{1}{T}\int_{0}^{T}\lim_{t^{\prime}\rightarrow\infty}\langle\tilde{\sigma}_{+}(\tau+t^{\prime})\tilde{\sigma}_{-}(t^{\prime})\rangle^{\ast}dt^{\prime}$
$\displaystyle=$ $\displaystyle\bar{\tilde{g}}_{1}^{\ast}(\tau),$

where we used relation (42) and the fact that
$\langle\tilde{\sigma}_{+}(\tau+t^{\prime}+T)\tilde{\sigma}_{-}(t^{\prime}+T)\rangle^{\ast}=\langle\tilde{\sigma}_{+}(\tau+t^{\prime})\tilde{\sigma}_{-}(t^{\prime})\rangle^{\ast}$
as $t^{\prime}\rightarrow\infty$. This means that the generalized parity
guarantees that the correlation function is a real-valued function of $\tau$
in the rotating frame and thus results in the symmetry of the spectrum when
$\delta+f(t)=-[\delta+f(t+T/2)]$. This is consistent with the prediction from
the spectrum (LABEL:eq:sfunsa).

In general, it is a formidable task to show that the spectrum is asymmetric
when $\delta+f(t)\neq-[\delta+f(t+T/2)]$ with or without the secular
approximation. Nevertheless, from the above derivation, one readily notes that
the generalized parity plays an important role in determining the symmetry of
the spectrum. Consequently, if such parity breaks, it is not difficult to
imagine that the symmetry of the spectrum also breaks trivially if there is no
other symmetry-inducing mechanism.

##  Appendix B Analytical calculation of quasienergies and transition matrix
elements in the biharmonic modulation case

We use the Van Vleck perturbation theory Cohen-Tannoudji  _et al._ (1998);
Hausinger and Grifoni (2010) to analytically calculate the quasienergies and
transition matrix elements $x_{\alpha\beta,l}^{(+)}$ for the biharmonic
modulation, which leads to the analytical fluorescence spectrum. Since we are
interested in the regime of $\Omega_{z},\,\omega_{z}\gg\Omega_{x}$, which is
accessible in the experiment Li  _et al._ (2013), we use $\Omega_{x}$ as the
perturbation parameter. We first transform Eq. (6) with the unitary
transformation

$e^{S(t)}[\tilde{H}(t)-i\partial_{t}]e^{-S(t)}e^{S(t)}|\tilde{u}_{\alpha}(t)\rangle=\tilde{\varepsilon}_{\alpha}e^{S(t)}|\tilde{u}_{\alpha}(t)\rangle,$
(44)

where

$S(t)=i\frac{\Omega_{z}}{2\omega_{z}}\left\\{\sin(\omega_{z}t)+\frac{r}{p}[\sin(p\omega_{z}t+\phi)-\sin\phi]\right\\}\sigma_{z}.$
(45)

We can define the transformed Floquet states and transformed Hamiltonian as
follows:

$|u_{\alpha}^{\prime}(t)\rangle=e^{S(t)}|\tilde{u}_{\alpha}(t)\rangle,$ (46)
$\displaystyle H^{\prime}(t)$ $\displaystyle=$ $\displaystyle
e^{S(t)}[\tilde{H}(t)-i\partial_{t}]e^{-S(t)}$ (47) $\displaystyle=$
$\displaystyle\frac{1}{2}\delta\sigma_{z}+\frac{1}{2}\sum_{l}(f_{l}\sigma_{+}+f_{-l}^{\ast}\sigma_{-})e^{il\omega_{z}t},$

where

$f_{l}=\Omega_{x}F_{l},$ (48)

and

$\displaystyle F_{l}$ $\displaystyle=$
$\displaystyle\frac{1}{T}\int_{0}^{T}e^{i\frac{\Omega_{z}}{\omega_{z}}\left\\{\sin(\omega_{z}t)+\frac{r}{p}[\sin(p\omega_{z}t+\phi)-\sin\phi]\right\\}-il\omega_{z}t}dt$
(49) $\displaystyle=$ $\displaystyle
e^{-i\Theta}\sum_{k}J_{k}\left(\frac{r\Omega_{z}}{p\omega_{z}}\right)J_{l-kp}\left(\frac{\Omega_{z}}{\omega_{z}}\right)e^{ik\phi},$

with $\Theta=\frac{r\Omega_{z}}{p\omega_{z}}\sin\phi$ and $J_{k}(z)$ being the
Bessel function of the first kind. To proceed, we introduce an extended
Hilbert space in which the time-dependent Floquet Hamiltonian
$H^{\prime}(t)-i\partial_{t}$ becomes time independent Sambe (1973). One
readily introduces the Fourier basis $|l\rangle\equiv\exp(il\omega_{z}t)$ and
inner product $\langle
l|n\rangle\equiv\frac{1}{T}\int_{0}^{T}\exp[i(n-l)\omega_{z}t]dt=\delta_{l,n}$,
where $\delta_{l,n}$ is the Kronecker delta function. Denoting
$|\uparrow\rangle$ and $|\downarrow\rangle$ as the eigenstates for
$\sigma_{z}$ with the eigenvalues $+1$ and $-1$, respectively, one gets the
composite bases
$|\uparrow(\downarrow),l\rangle=|\uparrow(\downarrow)\rangle\otimes|l\rangle$.
In the extended Hilbert space spanned by such bases, we can obtain the
explicit form of the Floquet Hamiltonian, which is written as

$\displaystyle H_{{\cal F}}^{\prime}$ $\displaystyle=$ $\displaystyle
H^{\prime}(t)-i\partial_{t}$ (50) $\displaystyle=$
$\displaystyle\frac{1}{2}\delta\sigma_{z}+\sum_{n}n\omega_{z}|n\rangle\langle
n|+\frac{1}{2}\sum_{n,l}(f_{l}\sigma_{+}+f_{-l}^{\ast}\sigma_{-})$
$\displaystyle\otimes|n+l\rangle\langle n|.$

The Floquet Hamiltonian has an infinite size and is difficult to be
diagonalized exactly in analytical calculation. To carry out perturbation
calculation, we transform the Floquet Hamiltonian with a further unitary
transformation with the Hermitian generator $K$, leading to

$\displaystyle H_{{\cal F}}^{\prime\prime}$ $\displaystyle=$ $\displaystyle
e^{iK}H_{{\cal F}}^{\prime}e^{-iK}$ (51) $\displaystyle=$ $\displaystyle
H_{\cal F}^{\prime}+[iK,H_{{\cal F}}^{\prime}]+\frac{1}{2!}[iK,[iK,H_{{\cal
F}}^{\prime}]]+\ldots,$

where the explicit form of $K$ is to be determined by requiring $H_{{\cal
F}}^{\prime\prime}$ to be block diagonal. The generator is expanded as

$K=K^{(1)}+K^{(2)}+K^{(3)}+\ldots,$ (52)

where the superscripts indicate the orders in the perturbation. We use
$H_{0}=\frac{1}{2}\delta\sigma_{z}+\sum_{n}n\omega_{z}|n\rangle\langle n|$ and
$V=\frac{1}{2}\sum_{n,l}(f_{l}\sigma_{+}+f_{-l}^{\ast}\sigma_{-})\otimes|n+l\rangle\langle
n|$ as the dominate and perturbation components, respectively. Up to the
second order in $\Omega_{x}$, we have

$\displaystyle H_{{\cal F}}^{\prime\prime}$ $\displaystyle\simeq$
$\displaystyle H_{0}+V+[iK^{(1)},H_{0}]+[iK^{(1)},V]+[iK^{(2)},H_{0}]$ (53)
$\displaystyle+\frac{1}{2}[iK^{(1)},[iK^{(1)},H_{0}]].$

Next, we discuss under which condition the transformed Hamiltonian may
reasonably be block diagonal. For the dominate component $H_{0}$, we simply
have
$H_{0}|\uparrow(\downarrow),n\rangle=[+(-)\delta/2+n\omega_{z}]|\uparrow(\downarrow),n\rangle\equiv\tilde{\varepsilon}^{(0)}_{+(-),n}|\uparrow(\downarrow),n\rangle$.
Provided that
$\tilde{\varepsilon}^{(0)}_{+,n}-\tilde{\varepsilon}^{(0)}_{-,n+m}=\delta-m\omega_{z}\approx
0$, we have a subspace spanned by two almost degenerate unperturbed states
$|\uparrow,n\rangle$ and $|\downarrow,n+m\rangle$, where $n$ is an arbitrary
integer and $m$ is the integer nearest to $\delta/\omega_{z}$. The projection
onto such a subspace is realized by the operator:

$\Pi_{n}=|\uparrow,n\rangle\langle\uparrow,n|+|\downarrow,n+m\rangle\langle\downarrow,n+m|.$
(54)

The eigenvalues of the dominate component $H_{0}$ in the $n$th subspace are
well-separated from those in the $(n+l)$th subspace as long as
$|l\omega_{z}|\gg|\delta-m\omega_{z}|$ for any $l\neq 0$. Moreover, if we
assume that

$|\langle\uparrow,n|V|\downarrow,n+l+m\rangle|\ll|\tilde{\varepsilon}^{(0)}_{+,n}-\tilde{\varepsilon}^{(0)}_{-,n+l+m}|,$
(55)

which is simply $|f_{-l-m}/2|\ll|l\omega_{z}|$, the transitions between the
states in the different subspaces can be neglected up to a certain order in
the perturbation Cohen-Tannoudji  _et al._ (1998), yielding the following
condition

$\Pi_{n}H_{{\cal F}}^{\prime\prime}\Pi_{l}=0,$ (56)

for $n\neq l$. Therefore, $H_{{\cal F}}^{\prime\prime}$ is block diagonal. The
second condition that $K$ cannot have matrix elements inside each subspace of
two almost degenerate states is also assumed, i.e.,

$\Pi_{n}K\Pi_{n}=0.$ (57)

The generator can now be fully determined via Eqs. (56) and (57). The
nonvanishing elements of $K^{(1)}$ and $K^{(2)}$are given by

$\langle\uparrow,n|iK^{(1)}|\downarrow,l\rangle=\frac{1}{2}\frac{f_{n-l}}{\delta+(n-l)\omega_{z}},$
(58)
$\langle\downarrow,l|iK^{(1)}|\uparrow,n\rangle=-\frac{1}{2}\frac{f_{n-l}^{\ast}}{\delta+(n-l)\omega_{z}},$
(59)

for $n-l\neq-m$, and

$\displaystyle\langle\uparrow,n|iK^{(2)}|\uparrow,l\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{4(n-l)\omega_{z}}\left\\{\sum_{k\neq
n+m,l+m}\frac{f_{n-k}f_{l-k}^{\ast}}{2}\left[\frac{1}{\delta+(n-k)\omega_{z}}+\frac{1}{\delta+(l-k)\omega_{z}}\right]\right.$
(60)
$\displaystyle\left.+\frac{f_{l-n-m}^{\ast}f_{-m}}{\delta+(l-n-m)\omega_{z}}+\frac{f_{n-l-m}f_{-m}^{\ast}}{\delta+(n-l-m)\omega_{z}}\right\\},$
$\displaystyle\langle\downarrow,n|iK^{(2)}|\downarrow,l\rangle$
$\displaystyle=$ $\displaystyle-\frac{1}{4(n-l)\omega_{z}}\left\\{\sum_{k\neq
l-m,n-m}\frac{f_{k-n}^{\ast}f_{k-l}}{2}\left[\frac{1}{\delta+(k-n)\omega_{z}}+\frac{1}{\delta+(k-l)\omega_{z}}\right]\right.$
(61)
$\displaystyle+\left.\frac{f_{l-n-m}^{\ast}f_{-m}}{\delta+(l-n-m)\omega_{z}}+\frac{f_{n-l-m}f_{-m}^{\ast}}{\delta+(n-l-m)\omega_{z}}\right\\},$

for $n\neq l$. The rest elements of $K^{(1)}$ and $K^{(2)}$ are vanishing.

The transformed Hamiltonian have the $2\times 2$ submatrix $H_{{\cal
F}}^{\prime\prime(n)}$ in the diagonal, which reads Cohen-Tannoudji  _et al._
(1998)

$\displaystyle H_{{\cal F}}^{\prime\prime(n)}$ $\displaystyle=$ $\displaystyle
H_{0}\Pi_{n}+\Pi_{n}V\Pi_{n}+\frac{1}{2}\Pi_{n}[iK^{(1)},V]\Pi_{n}$ (64)
$\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}\frac{\delta}{2}+n\omega_{z}+\sum_{j\neq-m}\frac{|f_{j}|^{2}}{4(\delta+j\omega_{z})}&\frac{f_{-m}}{2}\\\
\frac{f_{-m}^{\ast}}{2}&-\frac{\delta}{2}+(n+m)\omega_{z}-\sum_{j\neq-m}\frac{|f_{j}|^{2}}{4(\delta+j\omega_{z})}\end{array}\right).$

One can diagonalize the submatrix $H_{{\cal F}}^{\prime\prime(n)}$
analytically. Its eigenvalues (quasienergies) are

$\tilde{\varepsilon}_{\pm,n}=\frac{1}{2}\left(m\omega_{z}\pm\Omega_{m}\right)+n\omega_{z},$
(65)

where

$\Omega_{m}=\sqrt{\left[\delta-m\omega_{z}+\sum_{j\neq-m}\frac{|f_{j}|^{2}}{2(\delta+j\omega_{z})}\right]^{2}+|f_{-m}|^{2}}.$
(66)

The eigenvectors are given by

$\displaystyle|\Psi_{+,n}^{\prime\prime}\rangle$ $\displaystyle=$
$\displaystyle u|\uparrow,n\rangle+v|\downarrow,n+m\rangle,$ (67)
$\displaystyle|\Psi_{-,n}^{\prime\prime}\rangle$ $\displaystyle=$
$\displaystyle v|\uparrow,n\rangle-u^{\ast}|\downarrow,n+m\rangle,$ (68)

with

$\displaystyle u$ $\displaystyle=$
$\displaystyle\frac{f_{-m}}{|f_{-m}|}\sqrt{\frac{1}{2}\left[1+\frac{1}{\Omega_{m}}\left(\delta-m\omega_{z}+\sum_{j\neq-m}\frac{|f_{j}|^{2}}{2(\delta+j\omega_{z})}\right)\right]},$
(69) $\displaystyle v$ $\displaystyle=$
$\displaystyle\sqrt{\frac{1}{2}\left[1-\frac{1}{\Omega_{m}}\left(\delta-m\omega_{z}+\sum_{j\neq-m}\frac{|f_{j}|^{2}}{2(\delta+j\omega_{z})}\right)\right]}.$
(70)

The eigenvectors for $H_{{\cal F}}^{\prime}$ can be derived as follows:

$|\Psi_{\pm,n}^{\prime}\rangle=e^{-iK}|\Psi_{\pm,n}^{\prime\prime}\rangle\simeq\left(1-iK^{(1)}-iK^{(2)}+\frac{1}{2!}iK^{(1)}iK^{(1)}\right)|\Psi_{\pm,n}^{\prime\prime}\rangle.$
(71)

It is straightforward to derive the explicit form of the eigenvectors, which
reads

$\displaystyle|\Psi_{+,n}^{\prime}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{{\cal{\cal N}}}\left\\{uB|\uparrow,n\rangle-\sum_{j\neq
0}P_{j}|\uparrow,n+j\rangle+vB|\downarrow,n+m\rangle+\sum_{j\neq
0}Q_{j}|\downarrow,n+m+j\rangle\right\\},$ (72)
$\displaystyle|\Psi_{-,n}^{\prime}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{{\cal{\cal N}}}\left\\{vB|\uparrow,n\rangle+\sum_{j\neq
0}Q_{-j}^{\ast}|\uparrow,n+j\rangle-u^{\ast}B|\downarrow,n+m\rangle+\sum_{j\neq
0}P_{-j}^{\ast}|\downarrow,n+m+j\rangle\right\\},$ (73)

where

$B=1-\frac{1}{8}\sum_{l\neq-m}\frac{|f_{l}|^{2}}{(\delta+l\omega_{z})^{2}},$
(74) $\displaystyle P_{j}$ $\displaystyle=$
$\displaystyle\frac{f_{j-m}}{2[\delta+(j-m)\omega_{z}]}\left(v+\frac{uf_{-m}^{\ast}}{2j\omega_{z}}\right)+\frac{u}{4j\omega_{z}}\sum_{k\neq-m}\frac{f_{k+j}f_{k}^{\ast}}{\delta+k\omega_{z}},$
(75) $\displaystyle Q_{j}$ $\displaystyle=$
$\displaystyle\frac{f_{-j-m}^{\ast}}{2[\delta-(j+m)\omega_{z}]}\left(u+\frac{vf_{-m}}{2j\omega_{z}}\right)+\frac{v}{4j\omega_{z}}\sum_{k\neq-m}\frac{f_{k-j}^{\ast}f_{k}}{\delta+k\omega_{z}},$
(76)

and ${\cal N}=\sqrt{B^{2}+\sum_{j\neq 0}(|P_{j}|^{2}+|Q_{j}|^{2})}$ is the
normalization factor. The Floquet states $|u_{\alpha,n}^{\prime}(t)\rangle$
with the quasienergy $\tilde{\varepsilon}_{\alpha,n}$ can be derived from
$|\Psi_{\alpha,n}^{\prime}\rangle$ by replacing $|n\rangle$ with
$e^{in\omega_{z}t}$.

With above results at hand, we can analytically calculate the transition
matrix element

$\displaystyle x_{\alpha\beta,l}^{(+)}$ $\displaystyle=$
$\displaystyle\frac{1}{T}\int_{0}^{T}\langle\tilde{u}_{\alpha}(t)|\sigma_{\pm}|\tilde{u}_{\beta}(t)\rangle
e^{-il\omega_{z}t}dt=\frac{1}{T}\int_{0}^{T}\langle
u_{\alpha}^{\prime}(t)|e^{S(t)}\sigma_{+}e^{-S(t)}|u_{\beta}^{\prime}(t)\rangle
e^{-il\omega_{z}t}dt$ (77) $\displaystyle=$
$\displaystyle\sum_{n}\frac{1}{T}\int_{0}^{T}F_{n}\langle
u_{\alpha}^{\prime}(t)|\sigma_{+}|u_{\beta}^{\prime}(t)\rangle
e^{i(n-l)\omega_{z}t}dt=\sum_{n}F_{n+l}\langle\Psi_{\alpha,0}^{\prime}|\sigma_{+}|\Psi_{\beta,n}^{\prime}\rangle,$

and

$\displaystyle\langle\Psi_{+,0}^{\prime}|\sigma_{+}|\Psi_{+,n}^{\prime}\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{{\cal
N}^{2}}\left\\{u^{\ast}vB^{2}\delta_{n,-m}-\sum_{j\neq
0,n+m}P_{j}^{\ast}Q_{j-n-m}+(u^{\ast}Q_{-n-m}-vP_{n+m}^{\ast})B(1-\delta_{n,-m})\right\\},$
(78)
$\displaystyle\langle\Psi_{+,0}^{\prime}|\sigma_{+}|\Psi_{-,n}^{\prime}\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{{\cal
N}^{2}}\left\\{-(u^{\ast})^{2}B^{2}\delta_{n,-m}-\sum_{j\neq
0,n+m}P_{j}^{\ast}P_{n+m-j}^{\ast}+2u^{\ast}P_{n+m}^{\ast}B(1-\delta_{n,-m})\right\\},$
(79)
$\displaystyle\langle\Psi_{-,0}^{\prime}|\sigma_{+}|\Psi_{+,n}^{\prime}\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{{\cal
N}^{2}}\left\\{v^{2}B^{2}\delta_{n,-m}+\sum_{j\neq
0,n+m}Q_{-j}Q_{j-n-m}+2vQ_{-n-m}B(1-\delta_{n,-m})\right\\},$ (80)
$\displaystyle\langle\Psi_{-,0}^{\prime}|\sigma_{+}|\Psi_{-,n}^{\prime}\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{{\cal
N}^{2}}\left\\{-u^{\ast}vB^{2}\delta_{n,-m}+\sum_{j\neq
0,n+m}P_{j}^{\ast}Q_{j-n-m}+(vP_{n+m}^{\ast}-u^{\ast}Q_{-n-m})B(1-\delta_{n,-m})\right\\},$
(81)

where $(1-\delta_{n,-m})$ indicates that the term vanishes for $n=-m$.
Clearly, the validity of the perturbation theory is limited to the condition
(55). For $\delta\approx 0$, roughly speaking, the above results can be
justified when $r\sim 1$ and $\omega_{z}\sim\Omega_{z}\gg\Omega_{x}$.

##  Appendix C Equalities for transition matrix elements in the vanishing
detuning case

For the biharmonic modulation, we show the equalities that the transition
matrix elements satisfy under the vanishing detuning condition ($\delta=0$)
using the above analytical results, which helps us to understand the symmetry
of the spectrum in the main text. It follows from Eq. (49) that

$\displaystyle F_{-l}$ $\displaystyle=$ $\displaystyle
e^{-i\Theta}\sum_{k}J_{k}\left(\frac{r\Omega_{z}}{p\omega_{z}}\right)J_{-l-kp}\left(\frac{\Omega_{z}}{\omega_{z}}\right)e^{ik\phi}$
(82) $\displaystyle=$
$\displaystyle(-1)^{l}e^{-i\Theta}\sum_{k}J_{k}\left(\frac{r\Omega_{z}}{p\omega_{z}}\right)(-1)^{k(p+1)}$
$\displaystyle\times
J_{l-kp}\left(\frac{\Omega_{z}}{\omega_{z}}\right)e^{-ik\phi},$

where we used the relation $J_{-n}(z)=(-1)^{n}J_{n}(z)$. It is evident that
when $p$ is an odd number, $p+1$ is even and thus $(-1)^{k(p+1)}=1$, leading
to

$F_{-l}=(-1)^{l}e^{-i2\Theta}F_{l}^{\ast}.$ (83)

When $p$ is an even number, $(-1)^{k(p+1)}=(-1)^{k}$ may be either $+1$ or
$-1$. Nevertheless, we can obtain a simple relation between $F_{l}$ and
$F_{-l}$ by setting

$(-1)^{k}e^{-ik\phi}=e^{ik\phi},$ (84)

which yields that $\phi=\left(1/2+n\right)\pi$ $(n=0,\pm 1,\pm 2,\ldots)$.
With an even $p$ and such values of phase, we have

$F_{l}=(-1)^{l}F_{-l}.$ (85)

We should emphasize that Eqs. (83) and (85) hold under different conditions.
The former is available when $p$ is odd and regardless of $\phi$ while the
latter is established when $p$ is even and $\phi=(1/2+n)\pi$.

Provided that $\delta=0$, we get $m=\delta/\omega_{z}=0$. We define the phase
of $F_{0}$ via

$F_{0}=e^{-i\theta_{0}}|F_{0}|.$ (86)

Together with Eqs. (69) and (70), we simply have

$v=ue^{i\theta_{0}}$ (87)

with the aid of Eq. (83) or (85). Such an equality between $u$ and $v$ is
valid only for $\delta=0$ and in the valid regime of Eq. (83) or (85).

###  C.1 Odd $p$

We consider that $p$ is an odd number. It follows from Eq. (49) that
$\theta_{0}=\Theta$. Using $\delta=0$ and Eqs. (83) and (87), one readily gets
from Eqs. (75) and (76) that

$\displaystyle Q_{j}$ $\displaystyle=$
$\displaystyle-\frac{f_{-j}^{\ast}}{2j\omega_{z}}\left(u+\frac{vf_{0}}{2j\omega_{z}}\right)+\frac{v}{4j\omega_{z}}\sum_{k\neq
0}\frac{f_{k-j}^{\ast}f_{k}}{k\omega_{z}}$ (88) $\displaystyle=$
$\displaystyle\frac{(-1)^{j+1}e^{i2\Theta}f_{j}}{2j\omega_{z}}\left(u+\frac{vf_{0}^{\ast}e^{-i2\Theta}}{2j\omega_{z}}\right)+\frac{v}{4j\omega_{z}}\sum_{k\neq
0}\frac{f_{-k-j}^{\ast}f_{-k}}{-k\omega_{z}}$ $\displaystyle=$
$\displaystyle\frac{(-1)^{j+1}e^{i\Theta}f_{j}}{2j\omega_{z}}\left(v+\frac{uf_{0}^{\ast}}{2j\omega_{z}}\right)+\frac{e^{i\Theta}u}{4j\omega_{z}}\sum_{k\neq
0}\frac{(-1)^{j+1}f_{k+j}f_{k}^{\ast}}{k\omega_{z}}$ $\displaystyle=$
$\displaystyle(-1)^{j+1}e^{i\Theta}P_{j}.$

From this relation and Eqs. (77)-(80), it is straightforward to show that

$\displaystyle\left[x_{-+,-l}^{(+)}\right]^{\ast}$ $\displaystyle=$
$\displaystyle\sum_{n}\frac{F_{n-l}^{\ast}}{{\cal
N}^{2}}\left\\{v^{2}B^{2}\delta_{n,0}+\sum_{n\neq
0,n}Q_{-j}^{\ast}Q_{j-n}^{\ast}+2vBQ_{-n}^{\ast}(1-\delta_{n,0})\right\\}$
(89) $\displaystyle=$ $\displaystyle\sum_{n}\frac{F_{-n-l}^{\ast}}{{\cal
N}^{2}}\left\\{v^{2}B^{2}\delta_{n,0}+\sum_{j\neq
0,-n}Q_{-j}^{\ast}Q_{j+n}^{\ast}+2vBQ_{n}^{\ast}(1-\delta_{n,0})\right\\}$
$\displaystyle=$
$\displaystyle\sum_{n}\frac{(-1)^{n+l}F_{n+l}e^{i2\Theta}}{{\cal
N}^{2}}\left\\{v^{2}B^{2}\delta_{n,0}+\sum_{j\neq
0,n}Q_{j}^{\ast}Q_{n-j}^{\ast}+2vBQ_{n}^{\ast}(1-\delta_{n,0})\right\\}$
$\displaystyle=$
$\displaystyle\sum_{n}\frac{(-1)^{n+l}F_{n+l}e^{i2\Theta}}{{\cal
N}^{2}}\left\\{v^{2}B^{2}\delta_{n,0}+\sum_{j\neq
0,n}(-1)^{n}e^{-i2\Theta}P_{j}^{\ast}P_{n-j}^{\ast}+2vB(-1)^{n+1}e^{-i\Theta}P_{n}^{\ast}(1-\delta_{n,0})\right\\}$
$\displaystyle=$ $\displaystyle(-1)^{l}\sum_{n}\frac{F_{n+l}}{{\cal
N}^{2}}\left\\{(u^{\ast})^{2}B^{2}\delta_{n,0}+\sum_{j\neq
0,n}P_{j}^{\ast}P_{n-j}^{\ast}-2u^{\ast}BP_{n}^{\ast}(1-\delta_{n,0})\right\\}$
$\displaystyle=$ $\displaystyle-(-1)^{l}x_{+-,l}^{(+)}.$

Similarly, we find that
$\left[x^{(+)}_{++,-l}\right]^{\ast}=(-1)^{l}x^{(+)}_{++,l}$. Not
surprisingly, due to the generalized parity of the Floquet states, the
transition matrix elements satisfy Eq. (19) as long as
$\delta+f(t)=-[\delta+f(t+T/2)]$. For the biharmonic modulation, such
equalities are established when $p$ is odd and $\delta=0$.

###  C.2 Even $p$

We move to consider that $p$ is an even number. In such a case, the
generalized parity of the Floquet states is broken even if $\delta=0$. Thus,
we cannot expect that the transition matrix elements satisfy Eq. (19).
However, we have another type of equality. With Eqs. (85) and (87), one gets

$\displaystyle Q_{j}$ $\displaystyle=$
$\displaystyle\frac{f_{-j}^{\ast}}{-2j\omega_{z}}\left(u+\frac{vf_{0}}{2j\omega_{z}}\right)+\frac{v}{4j\omega_{z}}\sum_{k\neq
0}\frac{f_{k-j}^{\ast}f_{k}}{k\omega_{z}}$ (90) $\displaystyle=$
$\displaystyle\frac{(-1)^{j+1}f_{j}^{\ast}}{2j\omega_{z}}\left(u+\frac{vf_{0}}{2j\omega_{z}}\right)+\frac{v}{4j\omega_{z}}\sum_{k\neq
0}\frac{(-1)^{j+1}f_{j-k}^{\ast}f_{-k}}{-k\omega_{z}}$ $\displaystyle=$
$\displaystyle\frac{(-1)^{j+1}e^{-i\theta_{0}}f_{j}^{\ast}}{2j\omega_{z}}\left(v+\frac{u^{\ast}f_{0}}{2j\omega_{z}}\right)+\frac{e^{-i\theta_{0}}u^{\ast}}{4j\omega_{z}}\sum_{k\neq
0}\frac{(-1)^{j+1}f_{j+k}^{\ast}f_{k}}{k\omega_{z}}$ $\displaystyle=$
$\displaystyle(-1)^{j+1}e^{-i\theta_{0}}P_{j}^{\ast}.$

It is straightforward to derive Eqs. (23) and (24) in the main text via Eqs.
(77)-(80) and (90). We stress that the conditions for establishing such
relations require that $p$ is even, $\phi=(1/2+n)\pi$, and $\delta=0$.

## References

  * Mollow (1969) B. R. Mollow, Phys. Rev. 188, 1969 (1969).
  * Scully and Zubairy (1997) M. O. Scully and M. S. Zubairy,  _Quantum optics_ (Cambridge University Press, 1997).
  * Cohen-Tannoudji  _et al._ (1998) C. Cohen-Tannoudji, J. Dupont-Roc, G. Grynberg, and P. Thickstun,  _Atom-photon interactions: basic processes and applications_ (Wiley, 1998).
  * He  _et al._ (2013) Y.-M. He, Y. He, Y. Wei, D. Wu, M. Atature, C. Schneider, S. Hofling, M. Kamp, C. Lu, and J. Pan, Nat. Nanotech. 8, 213 (2013).
  * Santana  _et al._ (2017) T. S. Santana, Y. Ma, R. N. E. Malein, F. Bastiman, E. Clarke, and B. D. Gerardot, Phys. Rev. B 95, 201410(R) (2017).
  * Kiršanskė  _et al._ (2017) G. Kiršanskė, H. Thyrrestrup, R. S. Daveau, C. L. Dreeßen, T. Pregnolato, L. Midolo, P. Tighineanu, A. Javadi, S. Stobbe, R. Schott, A. Ludwig, A. D. Wieck, S. I. Park, J. D. Song, A. V. Kuhlmann, I. Söllner, M. C. Löbl, R. J. Warburton, and P. Lodahl, Phys. Rev. B 96, 165306 (2017).
  * Ficek and Freedhoff (1993) Z. Ficek and H. S. Freedhoff, Phys. Rev. A 48, 3092 (1993).
  * Agarwal  _et al._ (1991) G. S. Agarwal, Y. Zhu, D. J. Gauthier, and T. W. Mossberg, J. Opt. Soc. Am. B 8, 1163 (1991).
  * Ficek and Freedhoff (1996) Z. Ficek and H. S. Freedhoff, Phys. Rev. A 53, 4275 (1996).
  * Ficek and Rudolph (1999) Z. Ficek and T. Rudolph, Phys. Rev. A 60, R4245 (1999).
  * Peiris  _et al._ (2014) M. Peiris, K. Konthasinghe, Y. Yu, Z. C. Niu, and A. Muller, Phys. Rev. B 89, 155305 (2014).
  * Konthasinghe  _et al._ (2014) K. Konthasinghe, M. Peiris, and A. Muller, Phys. Rev. A 90, 023810 (2014).
  * He  _et al._ (2015) Y. He, Y.-M. He, J. Liu, Y.-J. Wei, H. Y. Ramírez, M. Atatüre, C. Schneider, M. Kamp, S. Höfling, C.-Y. Lu, and J.-W. Pan, Phys. Rev. Lett. 114, 097402 (2015).
  * Toyli  _et al._ (2016) D. M. Toyli, A. W. Eddins, S. Boutin, S. Puri, D. Hover, V. Bolkhovsky, W. D. Oliver, A. Blais, and I. Siddiqi, Phys. Rev. X 6, 031004 (2016).
  * Carmichael (1985) H. J. Carmichael, Phys. Rev. Lett. 55, 2790 (1985).
  * Grünwald and Vogel (2012) P. Grünwald and W. Vogel, Phys. Rev. Lett. 109, 013601 (2012).
  * Grünwald and Vogel (2013) P. Grünwald and W. Vogel, Phys. Rev. A 88, 023837 (2013).
  * Kimble  _et al._ (1977) H. J. Kimble, M. Dagenais, and L. Mandel, Phys. Rev. Lett. 39, 691 (1977).
  * D’Souza  _et al._ (1990) R. D’Souza, A. S. Jayarao, and S. V. Lawande, Phys. Rev. A 41, 4083 (1990).
  * Nazir (2008) A. Nazir, Phys. Rev. B 78, 153309 (2008).
  * Pastukhov  _et al._ (2014) V. M. Pastukhov, Y. V. Vladimirova, and V. N. Zadkov, Phys. Rev. A 90, 063831 (2014).
  * Itano  _et al._ (1988) W. M. Itano, J. C. Bergquist, and D. J. Wineland, Phys. Rev. A 38, 559 (1988).
  * Ficek  _et al._ (1984) Z. Ficek, R. Tanaś, and S. Kielich, Phys. Rev. A 29, 2004 (1984).
  * Damanet  _et al._ (2018) F. Damanet, J. Kübler, J. Martin, and D. Braun, Phys. Rev. A 97, 023832 (2018).
  * Kryuchkyan  _et al._ (2017) G. Y. Kryuchkyan, V. Shahnazaryan, O. V. Kibis, and I. A. Shelykh, Phys. Rev. A 95, 013834 (2017).
  * Antón  _et al._ (2017) M. A. Antón, S. Maede-Razavi, F. Carreño, I. Thanopulos, and E. Paspalakis, Phys. Rev. A 96, 063812 (2017).
  * Yan  _et al._ (2018) Y. Yan, Z. Lü, J. Y. Luo, and H. Zheng, Phys. Rev. A 97, 033817 (2018).
  * Saiko  _et al._ (2018) A. P. Saiko, S. A. Markevich, and R. Fedaruk, Phys. Rev. A 98, 043814 (2018).
  * Breuer and Petruccione (1997) H.-P. Breuer and F. Petruccione, Phys. Rev. A 55, 3101 (1997).
  * Yan  _et al._ (2016a) Y. Yan, Z. Lü, and H. Zheng, Ann. Phys. 371, 159 (2016a).
  * Roy and Hughes (2012) C. Roy and S. Hughes, Phys. Rev. B 85, 115309 (2012).
  * McCutcheon and Nazir (2013) D. P. S. McCutcheon and A. Nazir, Phys. Rev. Lett. 110, 217401 (2013).
  * Browne and Keitel (2000) D. E. Browne and C. H. Keitel, J. Mod. Opt. 47, 1307 (2000).
  * Yan  _et al._ (2013) Y. Yan, Z. Lü, and H. Zheng, Phys. Rev. A 88, 053821 (2013).
  * Ulrich  _et al._ (2011) S. M. Ulrich, S. Ates, S. Reitzenstein, A. Löffler, A. Forchel, and P. Michler, Phys. Rev. Lett. 106, 247402 (2011).
  * Ulhaq  _et al._ (2013) A. Ulhaq, S. Weiler, C. Roy, S. M. Ulrich, M. Jetter, S. Hughes, and P. Michler, Opt. Express 21, 4382 (2013).
  * Yan  _et al._ (2016b) Y. Yan, Z. Lü, H. Zheng, and Y. Zhao, Phys. Rev. A 93, 033812 (2016b).
  * Li  _et al._ (2013) J. Li, M. P. Silveri, K. S. Kumar, J. M. Pirkkalainen, A. Vepsäläinen, W. C. Chien, J. Tuorila, M. A. Sillanpää, P. J. Hakonen, E. V. Thuneberg, _et al._ , Nat. Commun. 4, 1420 (2013).
  * Pan  _et al._ (2017) J. Pan, H. Z. Jooya, G. Sun, Y. Fan, P. Wu, D. A. Telnov, S.-I. Chu, and S. Han, Phys. Rev. B 96, 174518 (2017).
  * Brunel  _et al._ (1998) C. Brunel, B. Lounis, P. Tamarat, and M. Orrit, Phys. Rev. Lett. 81, 2679 (1998).
  * Rohr  _et al._ (2014) S. Rohr, E. Dupont-Ferrier, B. Pigeau, P. Verlot, V. Jacques, and O. Arcizet, Phys. Rev. Lett. 112, 010502 (2014).
  * Kibis  _et al._ (2009) O. V. Kibis, G. Y. Slepyan, S. A. Maksimenko, and A. Hoffmann, Phys. Rev. Lett. 102, 023601 (2009).
  * Macovei and Keitel (2014) M. Macovei and C. H. Keitel, Phys. Rev. A 90, 043838 (2014).
  * Zhao  _et al._ (2015) Y.-J. Zhao, Y.-L. Liu, Y.-x. Liu, and F. Nori, Phys. Rev. A 91, 053820 (2015).
  * Silveri  _et al._ (2013) M. Silveri, J. Tuorila, M. Kemppainen, and E. Thuneberg, Phys. Rev. B 87, 134505 (2013).
  * Macovei  _et al._ (2015) M. Macovei, M. Mishra, and C. H. Keitel, Phys. Rev. A 92, 013846 (2015).
  * Silveri  _et al._ (2017) M. P. Silveri, J. A. Tuorila, E. V. Thuneberg, and G. S. Paraoanu, Rep. Prog. Phys. 80, 056002 (2017).
  * Wilkens and Rza¸ewski (1989) M. Wilkens and K. Rza¸ewski, Phys. Rev. A 40, 3164 (1989).
  * Das and Macovei (2013) S. Das and M. A. Macovei, Phys. Rev. B 88, 125306 (2013).
  * Ho  _et al._ (1986) T.-S. Ho, K. Wang, and S.-I. Chu, Phys. Rev. A 33, 1798 (1986).
  * Shirley (1965) J. H. Shirley, Phys. Rev. 138, B979 (1965).
  * Sambe (1973) H. Sambe, Phys. Rev. A 7, 2203 (1973).
  * Grossmann  _et al._ (1991) F. Grossmann, T. Dittrich, P. Jung, and P. Hänggi, Phys. Rev. Lett. 67, 516 (1991).
  * Lehmann  _et al._ (2003) J. Lehmann, S. Kohler, P. Hänggi, and A. Nitzan, J. Chem. Phys. 118, 3283 (2003).
  * Hausinger and Grifoni (2010) J. Hausinger and M. Grifoni, Phys. Rev. A 81, 022117 (2010).