text:
\section{Introduction}

Solid-state photonic qubits based on the fermion
systems coupled  to a quantized electromagnetic (EM) field in a plasmonic or dielectric nanocavity are promising for a variety of quantum information and quantum sensing applications \cite{haroche,lodahl2015,degen2017}.  Their benefits include compatibility with semiconductor technology, scalability, and potential for operation at temperatures much higher than the alternative platforms based on superconducting qubits or trapped ions. Indeed, strong coupling to single quantum emitters in dielectric nanocavities was demonstrated in various systems, for example color centers \cite{lukin2016} or quantum dots (QDs) Refs.~\cite{deppe,reithmaier}. In plasmonic cavities, strong coupling to single molecules  Refs.~\cite{chikkaraddy2016,benz2016, park2016} and colloidal QDs \cite{pelton2018, gross2018, park2019} has been achieved at room temperature; see, e.g., Refs.~\cite{lodahl2015,thorma2015,bitton2019,kono2019,may2020} for recent reviews.

While the quantum dynamics of entangled nanophotonic systems is interesting by itself, many applications would benefit from to control and modify the qubit states by time-dependent variation of certain parameters, while taking into account various processes of decoherence and dissipation. There is of course a large body of work related to cavity quantum electrodynamics (QED) with time-variable parameters. For example, the dynamics of nanophotonics systems with periodic modulation of some parameter, such as the cavity size or the position of a quantum emitter in a cavity, has been studied extensively in the burgeoning fields of cavity optomechanics \cite{aspelmeyer2014,meystre2013,pirkkalainen2015} and quantum acoustics \cite{chu2017,hong2017, arriola2019}. In this case the most interesting new element added to the nanophotonic system is the parametric resonance or the dressing of the electron-photon coupling by mechanical oscillations. Near the parametric resonance, the system can be mapped to an exactly solvable time-independent Hamiltonian  within the rotating-wave approximation  \cite{tokman2020}.

There is a class of time-dependent Hamiltonians for which the nonstationary Schr\"{o}dinger equation can be solved exactly in the analytic form, notably multistate Landau-Zener Hamiltonians and driven Tavis-Cummings Hamiltonians \cite{sinitsyn2018, chernyak2018}; see also \cite{li2018} where this technique was applied to the quantum annealing problem. Here we are interested in the nanophotonic applications, so we have to consider open multimode photonic systems with an arbitrary time dependence of the parameters. Therefore, we restrict ourselves to the adiabatic dynamics, for which the analytic solution can be found for a broad variety of systems with time-dependent cavity or fermion emitter parameters, and with dissipation included at the level of the Lindblad formalism. We find that the condition of adiabaticity is not that restrictive; in particular it still allows one to consider the parameter variation at a rate comparable to or faster than the generalized Rabi frequency in strongly coupled systems, which may be required for qubit manipulation.

We will also stick to the rotating wave
approximation (RWA)  \cite{Scully1997}.
The use of RWA restricts the coupling strength to the values much lower than
the characteristic energies in the system, such as the optical transition or
photon energy. The emerging studies of the so-called ultra-strong
coupling regime \cite{kono2019} have to go beyond the RWA. Nevertheless, for the
vast majority of experiments, including nonperturbative strong coupling
dynamics and entanglement, the RWA is adequate and provides some crucial
simplifications that allow one to obtain analytic solutions.

In particular, within Schr\"{o}dinger's description, the equations of motion
for the components of an infinitely dimensional state vector $\left\vert
\Psi \right\rangle $ that describes a coupled fermion-boson system can be
split into the blocks of low dimensions if the RWA is applied. This is true
even if the dynamics of the fermion subsystem is nonperturbative, e.g. the
effects of saturation are important. Note that there is no such
simplification in the Heisenberg representation, except within the perturbation theory; see e.g. \cite{Scully1997}. This is  because boson operators are defined on a basis of infinite dimension and truncation of their dynamics into blocks of small dimensions is generally not possible (see also \cite{tokman2020}).  The Schr\"{o}dinger's approach also leads to fewer equations
for the state vector components than the approach based on the von Neumann
master equation for the elements of the density matrix. This is especially true for a system with many degrees of freedom, e.g., many electron states coupled to multiple boson field modes.

Obviously, the Schr\"{o}dinger equation in its standard form cannot be
applied to describe open systems coupled to a dissipative reservoir. In this
case the stochastic versions of the equation of evolution for the state vector
have been developed, e.g. the method of quantum jumps \cite{Scully1997,Plenio1998}. This method is optimal for numerical analysis in
the Monte-Carlo type schemes. Here we formulate the stochastic equation which is more
conducive to the analytic treatment. In \cite{tokman2020} we showed that the stochastic equation of evolution for the state vector can be derived directly from the Heisenberg-Langevin formalism.

The paper is structured as follows. Section II formulates the model and the Hamiltonian for two-level electron system and a quantized EM field in a nanocavity with time-variable parameters. It treats a single-mode cavity in detail as a particular case and describes simple manipulations with a single cavity mode coupled to a single fermionic qubit. Section III considers the dynamics of two time-modulated cavity modes coupled to a single quantum emitter and Sec.~IV treats the case of a variable frequency of the optical transition in a fermion qubit.  Section V solves the quantum dynamics for an open time-dependent system with the coupling to dissipative reservoirs taken into account. An interesting phenomenon of modulation-induced transparency is analyzed. Numerical estimations for various nanophotonic systems reported in the literature are presented. Conclusions are in Section VI. Appendix A describes the quantization procedure for a plasmon cavity field with strongly subwavelength localization.  Appendix B summarizes the main properties of the stochastic equation of evolution and compares with the Lindblad density-matrix formalism.

\section{Cavity QED with time-variable parameters}

\subsection{Standard cavity QED Hamiltonian for a quantized field coupled to a two-level emitter}

For reference, we start from summarizing basic textbook  facts about a quantized electron system resonantly coupled to the quantum multimode EM field of a
nanocavity without any time dependence, and then consider the time-dependent models in the next sections.

\begin{figure}[htb]
\centering

\begin{subfigure}[b]{0.6\textwidth}
\includegraphics[width=1\linewidth]{fig1a.pdf}
\caption{ }
\label{fig1}
\end{subfigure}

\begin{subfigure}[b]{0.6\textwidth}
\includegraphics[width=1\linewidth]{fig1b.pdf}
\caption{}
\label{fig1}
\end{subfigure}

\caption { (a) A sketch of a quantum emitter (e.g. a quantum dot or a single molecule) in a nanocavity with time-dependent parameters created by
a metallic nanotip of the scanning probe and a metallic substrate. The profiles of the electric potential $ \Phi\left( \mathbf{r}, t\right)$ for the symmetric and antisymmetric mode (see Appendix A) are sketched. Other parameters are the transition energy $W(t)$ for a quantum emitter, the optical field frequency $\omega(t)$, the cavity height $d(t)$, and the relaxation constants of the cavity field, $\mu$, and a quantum emitter, $\gamma$.  (b) A quantum emitter coupled to the cavity surface plasmon field supported by graphene. The mode frequency $\omega(t)$ can be varied by applying variable voltage $V(t)$ which modifies the charge density in graphene. }

\end{figure}

Consider the simplest version of the fermion subsystem: two electron states $
\left\vert 0\right\rangle $ and $\left\vert 1\right\rangle $ with energies $
0 $ and $W$, respectively. We will call it an
``atom'' for brevity, although it can be electron states of a molecule, a
quantum dot, a defect in a semiconductor, or any other electron system. Introduce creation and
annihilation operators of the excited state $\left\vert 1\right\rangle $,
$
\hat{\sigma}=\left\vert 0\right\rangle \left\langle 1\right\vert$, $\hat{\sigma
}^{\dagger }=\left\vert 1\right\rangle \left\langle 0\right\vert$,
which satisfy standard commutation relations for fermions:
\begin{equation*}
\hat{\sigma}^{\dagger }\left\vert 0\right\rangle =\left\vert 1\right\rangle ,%
\hat{\sigma}\left\vert 1\right\rangle =\left\vert 0\right\rangle ,\hat{\sigma%
}\hat{\sigma}=\hat{\sigma}^{\dagger }\hat{\sigma}^{\dagger }=0;\left[ \hat{%
\sigma},\hat{\sigma}^{\dagger }\right] _{+}=\hat{\sigma}\hat{\sigma}%
^{\dagger }+\hat{\sigma}^{\dagger }\hat{\sigma}=1.
\end{equation*}%
The Hamiltonian of an atom is
\begin{equation}
\hat{H}_{a}=W\hat{\sigma}^{\dagger }\hat{\sigma}.  \label{ha}
\end{equation}%
We will also need the dipole moment operator, $
\mathbf{\hat{d}}=\mathbf{d}\left( \hat{\sigma}^{\dagger }+\hat{\sigma}%
\right)$,
where $\mathbf{d=}\left\langle 1\right\vert \mathbf{\hat{d}}\left\vert
0\right\rangle $ is a real vector. For a finite motion we can always choose
the coordinate representation of stationary states in terms of real
functions.

We assume that an atom is placed in a nanocavity and is resonantly coupled to the electric field of quantized cavity modes. Figure 1 sketches two out of many possible geometries of a time-variable nanocavity, e.g. formed by the nanotip of the scanning probe and the metallic substrate (Fig.~1a), similar to the recent experiments with strong coupling to single quantum emitters \cite{park2016, pelton2018, gross2018, park2019}. Of course many other cavity geometries are possible, such as the one in Fig.~1b where the quantum emitter is coupled to the cavity surface plasmon field supported by graphene.  Here the optical transition energy $W(t)$, the photon mode frequency $\omega(t)$, and field amplitudes described by an electric potential $ \Phi\left( \mathbf{r}, t\right)$  are all subject to external modulation by e.g.~varying the tip distance to the substrate, the position of a quantum emitter in a cavity, or a variable voltage applied to graphene or to a QD in a semiconductor nanostructure, but we will start from the Hamiltonian without any time dependence for future comparison.

We use a standard representation for the electric field operator in a cavity:%
\begin{equation}
\mathbf{\hat{E}}=\sum_{i}\left[ \mathbf{E}_{i}\left( \mathbf{r}\right) \hat{c%
}_{i}+\mathbf{E}_{i}^{\ast }\left( \mathbf{r}\right) \hat{c}_{i}^{\dagger }%
\right] ,  \label{ef}
\end{equation}%
where $\hat{c}_{i}^{\dagger },\hat{c}_{i}$ are creation and annihilation
operators for photons at frequency $\omega _{i}$; the functions $%
\mathbf{E}_{i}\left( \mathbf{r}\right) $ describe the spatial structure of
the EM modes in a cavity. The relation between the modal frequency $\omega
_{i}$ and the function $\mathbf{E}_{i}\left( \mathbf{r}\right) $ can
be found by solving the boundary-value problem of the classical
electrodynamics \cite{Scully1997}. The normalization conditions \cite{Tokman2016}
\begin{equation}
\int_{V}\frac{\partial \left[ \omega_i^{2}\varepsilon \left( \omega_i,\mathbf{r%
}\right) \right] }{\omega_i \partial \omega_i }\mathbf{E}_{i}^{\ast }\left(
\mathbf{r}\right) \mathbf{E}_{i}\left( \mathbf{r}\right) d^{3}r=4\pi \hbar
\omega _{i}  \label{nc-const}
\end{equation}%
ensure correct bosonic commutators $\left[ \hat{c}_{i},\hat{c}_{i}^{\dagger }%
\right] =\delta _{ij}$ and the field Hamiltonian in the form
\begin{equation}
\hat{H}_{em}=\hbar \sum_{i}\omega _{i}\left( \hat{c}_{i}^{\dagger }\hat{c}%
_{i}+\frac{1}{2}\right) .  \label{hem}
\end{equation}%
Here $V$ is a quantization volume  and $\varepsilon \left( \omega ,\mathbf{r}%
\right) $ is the dielectric function of a dispersive medium that fills the
cavity.

 Equation (\ref{nc-const}) is true for any fields satisfying Maxwell's equations as long as intracavity losses can be neglected and the flux of the Poynting vector through the total cavity surface is zero; see, e.g., Refs.~\cite{Tokman2016,Tokman2013, Tokman2015, Tokman2018}. Of course the photon losses are always important when calculating the decoherence rates and fluctuations. What matters for Eq.~(\ref{nc-const}) is that the effect of losses on the {\it spatial structure} of the cavity modes is insignificant. The latter is true as long as it makes sense to talk about cavity modes at all, which means in practice that the cavity Q-factor is at least around 10 or greater.

In many experiments involving strong coupling to a single quantum emitter the plasmonic cavities of nanometer size and even below 1 nm are used. The quantization procedure for a strongly subwavelength plasmon field has its peculiarities. We describe it in detail in Appendix A.

Adding
the interaction Hamiltonian with a EM cavity mode in the electric dipole
approximation, $-\mathbf{\hat{d}}\cdot \mathbf{\hat{E}}$, the Hamiltonian of
an atom coupled to a single mode EM field is
\begin{equation}
\hat{H}=\hat{H}_{em}+\hat{H}_{a}-\mathbf{d}\left( \hat{\sigma}^{\dagger }+%
\hat{\sigma}\right) \cdot \left[ \mathbf{E}\left( \mathbf{r}\right) \hat{c}+%
\mathbf{E}^{\ast }\left( \mathbf{r}\right) \hat{c}^{\dagger }\right] _{%
\mathbf{r}=\mathbf{r}_{a}},  \label{h}
\end{equation}%
where $\mathbf{r}=\mathbf{r}_{a}$ denotes the position of an atom inside the
cavity. This can be rewritten as
\begin{eqnarray}
\hat{H} &=&\hat{H}_{em}+\hat{H}_{a} - \left( \chi \hat{\sigma}
^{\dagger }\hat{c}+\chi ^{\ast }\hat{\sigma}\hat{c}^{\dagger }+\chi \hat{
\sigma}\hat{c}+\chi ^{\ast }\hat{\sigma}^{\dagger }\hat{c}^{\dagger }\right)
 \label{toh}
\end{eqnarray}
where
$ \chi =\left( \mathbf{d}\cdot \mathbf{E}\right) _{\mathbf{r}=
\mathbf{r}_{a}}$.

The best conditions for entanglement are realized in the vicinity of an atom-field resonance, where one can apply the rotating wave approximation (RWA). The RWA Hamiltonian is obtained by dropping the last two terms in Eq.~(\ref{toh}).
Note that we can always take the functions $\mathbf{E}\left( \mathbf{r}%
\right) $ to be real at
the position of an atom. This single-mode model is of course the Jaynes-Cummings Hamiltonian \cite{JC1963}.

\subsection{Quantized electromagnetic field in a time-variable cavity}

In a standard approach to quantization of the EM field based on Eqs.~(\ref{ef})-(\ref{hem}), a set of mode
frequencies $\omega _{i}$ and the relation between the frequency $%
\omega _{i}$ and the spatial structure of the field mode $\mathbf{E}
_{i}\left( \mathbf{r}\right) $ are determined by solving the boundary-value
problem for the classical EM field. Let's assume that the solution of this
boundary-value problem depends on a certain parameter $p$, for example the
cavity height $d(t)$ in Fig.~1 or the position of the emitter with respect to the field
distribution. In this case $\omega _i(p)$ and $%
\mathbf{E}_{i}\left( \mathbf{r,}p\right) $ are functions of $p$. Of course
the solution depends on many parameters of the cavity, but we consider the
situation when this particular parameter is adiabatically changing with
time. As usual, ``adiabatically'' means that the change can be
arbitrary (e.g. periodic or not) but it should be slow as compared to typical frequencies of all subsystems when the parameters are constant, such as
the modal frequencies and the transition frequency of a quantum emitter. It is important that the rate of change of parameters can be arbitrary as compared to characteristic frequency scales which determine the interaction between subsystems, such as the Rabi frequency, as long as these scales are smaller than the modal or transition frequencies \cite{kruskal1962,tokman2011}.

For an adiabatically varying parameter Eqs.~(\ref{ef})-(\ref{hem}) depend on the
instantaneous value of the parameter,
\begin{equation}
\mathbf{\hat{E}}=\sum_{i}\left[ \mathbf{E}_{i}\left( \mathbf{r,}p\right)
\hat{c}_{i}+\mathbf{E}_{i}^{\ast }\left( \mathbf{r,}p\right) \hat{c}%
_{i}^{\dagger }\right] ,  \label{elef}
\end{equation}
\begin{equation}
\hat{H}\ =\sum_{i}\hat{H}_{i}\ \ \ \hat{H}_{i}=\hbar \omega _i(p) \left( \hat{c}_{i}^{\dagger }\hat{c}_{i}+\frac{1}{2}\right)
,   \label{Hem}
\end{equation}
\begin{equation}
\int_{V\left( p\right) }\frac{\partial \left[ \omega _{i}^{2}\varepsilon
\left( \omega _{i},\mathbf{r}\right) \right] }{\omega _{i}\partial \omega
_{i}}\mathbf{E}_{i}^{\ast }\left( \mathbf{r,}p\right) \mathbf{E}_{i}\left(
\mathbf{r,}p\right) d^{3}r=4\pi \hbar \omega _i (p).
\label{NC}
\end{equation}%
The solution of the Schr\"{o}dinger
equation $i\hbar \frac{\partial }{\partial t}\left\vert \Psi
_{i}\right\rangle =\hat{H}_{i}\left\vert \Psi _{i}\right\rangle $ for a
given field mode is
\begin{equation}
\left\vert \Psi _{i}\right\rangle =\sum_{n=0}^{\infty
}   C_{n}\left\vert n\right\rangle  \label{Sol of SE}
\end{equation}%
where $C_n = C_n^0 e^{-i(n + \frac{1}{2}) \int_{0}^{t}\omega _{i}\left( \tau \right) d\tau
}$, $\omega _{i}\left( t\right) \equiv $ $\omega _{i%
}\left( p\left( t\right) \right) $, and $\left\vert n\right\rangle $ are Fock
states. For a bosonic field described by a standard quantized harmonic oscillator, if we choose the coordinate representation expressed via
Hermite polynomials, the parameters of the polynomials will be
time-dependent. One can easily see that the above solution conserves the
adiabatic invariant $\frac{\left\langle \Psi _{i}\right\vert \hat{H}%
_{i}\left\vert \Psi _{i}\right\rangle }{\omega _{i}\left( t\right) }$
, just like in a classical slowly time-varying harmonic oscillator problem \cite{LL1}.

\subsection{Quantum emitter coupled to a quantized EM field
with a time-variable amplitude}

Let a two-level electron system (an atom)  be located at the point $\mathbf{r}=\mathbf{0}$
inside the cavity. The Hamiltonian of the system including the coupling of
an atom to a particular cavity mode and its adiabatic modulation can be
described within the RWA,
\begin{equation}
\hat{H}=\hbar \omega(t) \left( \hat{c}^{\dagger }%
\hat{c}+\frac{1}{2}\right) +W\hat{\sigma}^{\dagger }\hat{\sigma}-\left[
\mathbf{d}\cdot \mathbf{E}\left( t\right) \hat{\sigma}^{\dagger }\hat{c}+%
\mathbf{d}^{\ast }\cdot \mathbf{E}^{\ast }\left( t\right) \hat{\sigma}\hat{c}%
^{\dagger }\right] ,  \label{H of sys with RWA}
\end{equation}%
where $\mathbf{E}\left( \mathbf{0},t\right) =\mathbf{E}\left( t\right) $. The time dependence of the field
amplitude follows from the parameter modulation.

The wave function of the coupled photon-electron state can be written as
\begin{equation}
\label{wavefun}
\Psi =\sum_{n=0}^{\infty } \left( C_{n0}\left\vert n\right\rangle \left\vert
0\right\rangle +C_{n1}\left\vert n\right\rangle \left\vert 1\right\rangle \right),
\end{equation}
 where we will maintain the same order, $\left\vert
photon\right\rangle \left\vert fermion\right\rangle $ of the state products
everywhere. Substituting it in the Schr\"{o}dinger equation and taking into
account the time variation of the parameter, we obtain the equation for the ground energy state,
\begin{equation}
\dot{C}_{00}+i\omega _{00}\left( t\right) C_{00}=0,  \label{Eq for les}
\end{equation}
and a pair of equations for ``resonant'' states,
\begin{equation}
\dot{C}_{n0}+i\omega _{n0}\left( t\right) C_{n0}-i\Omega _{R}^{\ast }\left(
t\right) C_{\left( n-1\right) 1}=0,  \label{Cn0 dot}
\end{equation}%
\begin{equation}
\dot{C}_{\left( n-1\right) 1}+i\omega _{\left( n-1\right) 1}\left( t\right)
C_{\left( n-1\right) 1}-i\Omega _{R}\left( t\right) C_{n0}=0,
\label{Cn-11 dot}
\end{equation}%
where
\begin{equation*}
\omega _{n0}\left( t\right) =\omega \left( t\right) \left( n+\frac{1}{2}%
\right) ,\ \ \omega _{n1}\left( t\right) =\omega _{n0}\left( t\right) +\frac{%
W}{\hbar },\ \ \Omega _{R}\left( t\right) =\frac{\mathbf{d\cdot E}\left(
t\right) }{\hbar }\sqrt{n}.
\end{equation*}%

Equations (\ref{Cn0 dot}), (\ref{Cn-11 dot}) can be written in a more convenient form after making a substitution
\begin{equation}
\left(
\begin{array}{c}
C_{n0} \\
C_{\left( n-1\right) 1}%
\end{array}%
\right) =\left(
\begin{array}{c}
G_{n0}e^{-i\int_{0}^{t}\omega _{n0}\left( \tau \right) d\tau } \\
G_{\left( n-1\right) 1}e^{-i\int_{0}^{t}\omega _{\left( n-1\right) 1%
}\left( \tau \right) d\tau }%
\end{array}%
\right) ,  \label{Cn0 and Cn-11}
\end{equation}%
which gives
\begin{equation}
\dot{G}_{n0}-i\Omega _{R}^{\ast }\left( t\right) e^{i\int_{0}^{t}\delta
\left( \tau \right) d\tau }G_{\left( n-1\right) 1}=0,
\label{Gn0 dot}
\end{equation}
\begin{equation}
\dot{G}_{\left( n-1\right) 1}-i\Omega _{R}\left( t\right)
e^{-i\int_{0}^{t}\delta \left( \tau \right) d\tau }G_{n0}=0,
\label{Gn-11 dot}
\end{equation}%
where $\delta \left( t\right) =\omega \left( t\right) -\frac{W}{\hbar }$.

When there is no modulation, i.e. $\delta$, $\omega$, and $\Omega_R$ are constant, Eqs.~(\ref{Gn0 dot}), (\ref{Gn-11 dot}) have a simple solution $G_{(n-1)1},G_{n0} \propto e^{-i\nu t}$, where the eigenvalues are
\begin{equation}
    \nu_{1,2} = \frac{\delta}{2} \pm \sqrt{\frac{\delta^2}{4} + |\Omega_R|^2 },
    \label{nu12}
\end{equation}
and the eigenvectors satisfy
\begin{equation}
    K_{1,2} = \left[ \frac{G_{n0}}{G_{(n-1)1}} \right]_{1,2} = \frac{\nu_{1,2} e^{i \delta t}}{\Omega_R},
\end{equation}
where $K_1 K_2^{\ast} = -1$. The eigenvalues $\nu_{1,2}$ as a function of detuning $\delta$ are shown in Fig.~2. It is easy to verify that in the region $\delta \ll -|\Omega_R|$   the eigenvalue $\nu_1$ corresponds to the dominant state $\left\vert n-1\right\rangle \left\vert 1\right\rangle$, whereas in the region $\delta \gg |\Omega_R|$   this eigenvalue corresponds to dominant state $\left\vert n\right\rangle \left\vert 0\right\rangle$. For the eigenvalue $\nu_2$ it is exactly the opposite.

\begin{figure}[htb]
\centering
\begin{subfigure}[b]{0.8\textwidth}
\includegraphics[width=1\linewidth]{fig2.pdf}
\caption{ }
\label{fig2}
\end{subfigure}

\caption { Frequency eigenvalues $\nu_{1,2}$ from Eq.~(\ref{nu12}) as a function of detuning $\delta$ from the  resonance, $\delta = \omega - \frac{W}{\hbar}$. All frequencies are in units of the average Rabi frequency $\overline{\Omega _{R}}$.}
\end{figure}

When a cavity parameter is modulated, for example, a cavity height $d(t)$ in Fig.~1, both frequencies and field amplitudes $\mathbf{E}_{i}\left(\mathbf{r,}p\right)$ get modulated; see Eq.~(\ref{NC}). Therefore, the Rabi frequency $\Omega_R$ gets modulated. For a periodic modulation, the function $ \Omega _{R}\left( t\right)
e^{-i\int_{0}^{t}\delta \left( \tau \right) d\tau }$ is periodic and can be expanded in the Fourier series,
\begin{equation}
    \Omega _{R}\left( t\right)
e^{-i\int_{0}^{t}\delta \left( \tau \right) d\tau } = \sum_{n = -\infty}^{\infty} R_n e^{-i n \Omega t},
\label{fourier}
\end{equation}
where $\Omega$ is the modulation frequency. The explicit expressions for the Fourier amplitudes $R_n$ can be obtained  for any specific model of a cavity; see, e.g., Appendix A for the plasmonic cavity, which shows specific examples of the cavity mode frequencies, field amplitudes, and their modulation.

When the modulation frequency and amplitude of the eigenmode frequencies are small enough, one can neglect the modulation of the Rabi frequency in Eqs.~(\ref{Gn0 dot}), (\ref{Gn-11 dot}). This corresponds to the WKB approximation and one can see it by
 taking the time derivative of Eq.~(\ref{Gn-11 dot}):
\begin{equation}
\frac{d^{2}G_{\left( n-1\right) 1}}{dt^{2}} + \left( i\delta
\left( t\right) - \frac{1}{\Omega _{R}\left(
t\right) }  \frac{d\Omega _{R}}{dt} \right) \frac{dG_{\left( n-1\right) 1}}{dt} +\left\vert \Omega
_{R}\left( t\right) \right\vert ^{2}G_{\left( n-1\right) 1}=0.
\label{Gn-11 double dot}
\end{equation}
Now we can estimate the order of magnitude of different terms in Eq.~(\ref{Gn-11 double dot}).
Assume that the cavity mode frequency is modulated as $\omega_i(t) = \bar{\omega} + \delta(t) $. Since Eq.~(\ref{NC}) defines a certain dependence $\Omega_R(\omega_i)$, one can estimate $\left|\frac{\Omega_R^{-1} \dot{\Omega}_R}{\delta}\right| \sim \left| \frac{\dot{\delta}}{\Omega_R \delta} \frac{d \Omega_R}{d\omega_i}\right|_{\omega_i = \bar{\omega}}$ and $\left\vert \Omega _{R}-\overline{\Omega _{R}}\right\vert \sim \left| \frac{d\Omega_R}{d\omega_i} \delta \right|$. For estimations we take $\frac{d\Omega_R}{d\omega_i} \sim \frac{\overline{\Omega _{R}}}{\bar{\omega}} $, $\delta \sim \Delta \omega$, and $\dot{\delta} \sim \Omega \Delta \omega$ where $\overline{\Omega _{R}} = \Omega_R(\bar{\omega})$ and $\Delta \omega$ is the frequency change over the time $\Omega^{-1}$. This gives $\left|\frac{\Omega_R^{-1} \dot{\Omega}_R}{\delta}\right| \sim \frac{\Omega}{ \bar{\omega}}$ and $\left\vert \Omega _{R}-\overline{\Omega _{R}}\right\vert \sim \overline{%
\Omega _{R}}\frac{\Delta \omega }{\bar{\omega}}$. If $\Delta \omega ,\Omega
\ll \bar{\omega}$, Eq.~(\ref{Gn-11 double dot}) becomes
\begin{equation}
\frac{d^{2}G_{\left( n-1\right) 1}}{dt^{2}}+i\delta \left( t\right) \frac{%
dG_{\left( n-1\right) 1}}{dt}+\left\vert \overline{\Omega _{R}}\right\vert
^{2}G_{\left( n-1\right) 1}=0.  \label{Gn-11 double dot for sc}
\end{equation}%
Equation~(\ref{Gn-11 double dot for sc}) corresponds to the set of Eqs.~(\ref{Cn0 dot}),~(\ref{Cn-11 dot}) with $\Omega _{R} = $ const $ = \overline{\Omega _{R}}$.

If we consider for definiteness a sinusoidal
modulation of the frequency of a given mode,
$
\omega \left( t\right) =\bar{\omega}-\Delta \omega \cdot \sin \left( \Omega
t\right)$, and take into account that $\Delta \omega, \Omega, \overline{\Omega_R} \ll \overline{\omega}$, the Fourier amplitudes in Eq.~(\ref{fourier}) can be expressed through the Bessel functions,
\begin{equation}
    R_0 = \overline{\Omega _{R}} J_0\left(\frac{\Delta \omega}{\Omega}\right), \; R_n = (-i)^{|n|}  \overline{\Omega _{R}} J_{|n|}\left(\frac{\Delta \omega}{\Omega}\right).
    \label{amplitudes}
\end{equation}

The decoherence processes can be added within the stochastic equation of evolution for the state vector, which is derived in Appendix B. However, we postpone doing this until we consider a more complex case of two quantized modes interacting with a quantum emitter.

\subsection{Simple manipulations with a qubit coupled to a single-mode field}

A single emitter coupled to a single-mode field in a time-variable cavity permits simple manipulations: a slow or fast sweep through the resonance $\omega(t) = W/\hbar$, bringing an electron-photon system in and out of entanglement by changing the values of coefficients in Eq.~(\ref{wavefun}), transduction of the excitation between an atom and the EM field, e.g., between  $\left\vert 0\right\rangle \left\vert
1 \right\rangle$ and $\left\vert 1\right\rangle \left\vert
0\right\rangle$ states, etc.

Note that the rate of modulation or parameter variation has to be slow only as compared to the optical frequency. It does not have to be slow as compared to the average Rabi frequency $\overline{\Omega _{R}}$. Therefore, in the strong coupling regime a desired switching can be completed faster than the Rabi oscillations and decoherence rates.

Let's look at some of these control operations in more detail. The sweep through resonance can be calculated exactly for each specific time dependence $\delta(t)$, but the limiting cases are well understood from the vast amount of literature on the linear coupling of the optical modes, Landau-Zener-type problems, etc \cite{kruskal1962,tokman2011,kochar1983,hallin1995,yokomizo2014}.

For a slow sweep, $\left| \frac{d\delta}{dt} \right| \ll \overline{|\Omega_R|^2}$, the system will follow each eigenvalue branch plotted in Fig.~2 without jumping between them: for example, if the system starts from $\nu_1$ at $\delta \ll -\overline{|\Omega_R|}$, it will stay on $\nu_1$ as it moves through resonance to $\delta \gg \overline{|\Omega_R|}$. This means that the quantum state of the system will be switched from $\left\vert n-1\right\rangle \left\vert 1\right\rangle$ to $\left\vert n\right\rangle \left\vert 0\right\rangle$.

In the opposite limit of a fast sweep, $\left| \frac{d\delta}{dt} \right| \gg \overline{|\Omega_R|^2}$, as the system moves through resonance from $\delta \ll -\overline{|\Omega_R|}$ to $\delta \gg \overline{|\Omega_R|}$ it jumps from one eigenvalue branch to another. As a result, the quantum state stays unchanged.

In the intermediate region $\left| \frac{d\delta}{dt} \right| \sim \overline{|\Omega_R|^2}$, by varying the sweep rate or the Rabi frequency  $\overline{|\Omega_R|}$ one can get any desired combination of the quantum states at the output. In particular, for linear variation of the detuning, $\delta(t) = \beta t$ where $\beta$ is a constant, one can obtain an exact analytic solution of Eq.~(\ref{Gn-11 double dot for sc}) to predict the evolution of the system:
\begin{equation}
G_{(n-1)1}(t) =\displaystyle  e^{-\frac{i \beta t^2}{4}} \left[ c_1 D_{i \frac{\overline{|\Omega_R|}^2}{\beta}} \left( \sqrt{\beta} e^{-\frac{i \pi}{4}} t \right) + c_2 D_{-i \frac{\overline{|\Omega_R|}^2}{\beta}-1} \left( i \sqrt{\beta} e^{-\frac{i \pi}{4}} t \right)  \right],
\label{linear}
\end{equation}
where $D_{\nu}$ are the parabolic cylinder functions \cite{bateman} and $c_{1,2}$ are arbitrary constants determined by initial conditions. This solution can be used, for example, to calculate the efficiency of the $\left\vert n-1\right\rangle \left\vert 1\right\rangle$ quantum state tunneling, i.e., the probability of the transition from the top to bottom branch in Fig.~2 as the detuning $\delta(t)$ varies from $-\infty$ to $+\infty$:
$$
\left| C_{(n-1)1} \right|_{\delta \rightarrow \infty}^2 \approx \displaystyle e^{-\frac{\overline{2 \pi |\Omega_R|}^2}{\beta}} \left| C_{(n-1)1} \right|_{\delta \rightarrow - \infty}^2.
$$
As expected, the probability is approaching 1 when $\left| \frac{d\delta}{dt} \right| = \beta \gg \overline{|\Omega_R|}^2$ and becomes exponentially small in the opposite limit.

\section{Dynamics of two modulated cavity modes coupled to a
quantum emitter}

In order to perform more complex operations on the photonic qubits and get more functionality, we need to add one more quantized degree of freedom to the system. Here we consider {\it two} cavity modes in a time-variable cavity,
\begin{equation}
\mathbf{\hat{E}}=\mathbf{E}_{a}\left( \mathbf{r,}t\right) \hat{a}+\mathbf{E}%
_{a}^{\ast }\left( \mathbf{r,}t\right) \hat{a}^{\dagger }+\mathbf{E}%
_{b}\left( \mathbf{r,}t\right) \hat{b}+\mathbf{E}_{b}^{\ast }\left( \mathbf{%
r,}t\right) \hat{b}^{\dagger }.  \label{two mode ef}
\end{equation}%

We assume that the modulation of both frequencies has a small amplitude and average frequencies of both modes $\bar{\omega}_{a,b}$
are close to the transition frequency. In this case the RWA Hamiltonian for
an atom + field system is
\begin{equation}
\hat{H}=\hbar \omega _{a}\left( t\right) \left( \hat{a}^{\dagger }%
\hat{a}+\frac{1}{2}\right) +\hbar \omega _{b}\left( t\right) \left(
\hat{b}^{\dagger }\hat{b}+\frac{1}{2}\right) +W\hat{\sigma}^{\dagger }\hat{%
\sigma}-\left[ \hat{\sigma}^{\dagger }\left( \chi _{a}\hat{a}+\chi _{b}\hat{b%
}\right) +\hat{\sigma}\left( \chi _{a}^{\ast }\hat{a}^{\dagger }+\chi
_{b}^{\ast }\hat{b}^{\dagger }\right) \right] ,  \label{toh in RWA}
\end{equation}%
where $\chi _{a,b}\left( t\right) =\mathbf{d\cdot E}_{a,b}\left( t\right) $.

The Schr\"{o}dinger equation can be solved analytically within the RWA
\cite{tokman2020}. As a simple example, we include only the transitions between
the states with lowest energies, namely $\left\vert 0_{a}\right\rangle
\left\vert 0_{b}\right\rangle \left\vert 0\right\rangle ,\left\vert
0_{a}\right\rangle \left\vert 0_{b}\right\rangle \left\vert 1\right\rangle
,\left\vert 1_{a}\right\rangle \left\vert 0_{b}\right\rangle \left\vert
0\right\rangle ,\left\vert 0_{a}\right\rangle \left\vert 1_{b}\right\rangle
\left\vert 0\right\rangle $, i.e. we seek the solution in the form%
\begin{equation}
\Psi =C_{000}\left\vert 0_{a}\right\rangle \left\vert 0_{b}\right\rangle
\left\vert 0\right\rangle +C_{001}\left\vert 0_{a}\right\rangle \left\vert
0_{b}\right\rangle \left\vert 1\right\rangle +C_{100}\left\vert
1_{a}\right\rangle \left\vert 0_{b}\right\rangle \left\vert 0\right\rangle
+C_{010}\left\vert 0_{a}\right\rangle \left\vert 1_{b}\right\rangle
\left\vert 0\right\rangle.   \label{sol of SE}
\end{equation}

For arbitrary coefficients $C$ the state (\ref{sol of SE}) is a tripartite entangled state which can be reduced to standard GHZ states by  local operations  \cite{dur2000, cunha2020}, e.g. by rotations on the Bloch sphere of each qubit. In most cases discussed in the literature the GHZ states are made of identical subsystems, e.g., photons \cite{shalm2012, agusti2020}. In our case the subsystems are of different nature: a fermionic electron system and bosonic EM field modes. This makes their rotations more complicated, but on the other hand, enables other interesting applications. For example, one can determine the statistics of atomic excitations by measuring the statistics of photons, or change the entangled state of coupled photon modes by changing the atomic state with a classical control field.

Similarly to \cite{tokman2020}, the equations for the coefficients are
\begin{equation}
\dot{C}_{000}+i\frac{\omega _{a}\left( t\right) +\omega _{b}\left( t\right)
}{2}C_{000}=0;  \label{c000 dot}
\end{equation}%
\begin{equation}
\dot{C}_{001}+i\left( \frac{1}{2}\omega _{a}\left( t\right) +\frac{1}{2}%
\omega _{b}\left( t\right) +\frac{W}{\hbar }\right) C_{001}-i\Omega
_{Ra}(t) C_{100}-i\Omega _{Rb}(t) C_{010}=0,  \label{c001 dot}
\end{equation}%
\begin{equation}
\dot{C}_{100}+i\left( \frac{3}{2}\omega _{a}\left( t\right) +\frac{1}{2}%
\omega _{b}\left( t\right) \right) C_{100}-i\Omega _{Ra}^{\ast }(t)C_{001}=0,
\label{c100 dot}
\end{equation}%
\begin{equation}
\dot{C}_{010}+i\left( \frac{1}{2}\omega _{a}\left( t\right) +\frac{3}{2}%
\omega _{b}\left( t\right) \right) C_{010}-i\Omega _{Rb}^{\ast }(t)C_{001}=0,
\label{c010 dot}
\end{equation}%
where $\Omega _{Ra,b}=\frac{\chi _{a,b}}{\hbar }$. Making the substitution%
\begin{equation}
\left(
\begin{array}{c}
C_{001} \\
C_{100} \\
C_{010}%
\end{array}%
\right) =\left(
\begin{array}{c}
G_{0}\exp \left[ -i\int_{0}^{t}\left( \frac{1}{2}\omega _{a}\left( \tau
\right) +\frac{1}{2}\omega _{b}\left( \tau \right) +\frac{W}{\hbar }\right)
d\tau \right] \\
G_{a}\exp \left[ -i\int_{0}^{t}\left( \frac{3}{2}\omega _{a}\left( \tau
\right) +\frac{1}{2}\omega _{b}\left( \tau \right) \right) d\tau \right] \\
G_{b}\exp \left[ -i\int_{0}^{t}\left( \frac{1}{2}\omega _{a}\left( \tau
\right) +\frac{3}{2}\omega _{b}\left( \tau \right) \right) d\tau \right]%
\end{array}%
\right) ,  \label{substitution}
\end{equation}%
we obtain%
\begin{equation}
\dot{G}_{0}-i\Omega _{Ra}(t)G_{a}\exp \left[ -i\int_{0}^{t}\left( \omega
_{a}\left( \tau \right) -\frac{W}{\hbar }\right) d\tau \right] -i\Omega
_{Rb}(t) G_{b}\exp \left[ -i\int_{0}^{t}\left( \omega _{b}\left( \tau \right) -%
\frac{W}{\hbar }\right) d\tau \right] =0,  \label{g0 dot}
\end{equation}%
\begin{equation}
\dot{G}_{a}-i\Omega _{Ra}^{\ast }(t) G_{0}\exp \left[ i\int_{0}^{t}\left( \omega
_{a}\left( \tau \right) -\frac{W}{\hbar }\right) d\tau \right] =0,
\label{ga dot}
\end{equation}%
\begin{equation}
\dot{G}_{b}-i\Omega _{Rb}^{\ast }(t) G_{0}\exp \left[ i\int_{0}^{t}\left( \omega
_{b}\left( \tau \right) -\frac{W}{\hbar }\right) d\tau \right] =0,
\label{gb dot}
\end{equation}%

In Fig.~\ref{Fig:eigen_2photon}, we show the eigenstates of the system described by Eqs.~\eqref{c001 dot}, \eqref{c100 dot} and \eqref{c010 dot} as a function of frequency detuning defined as $\omega_a-W/\hbar$.  Here we assumed that $\Omega_{Ra}=\Omega_{Rb}\equiv \Omega_R$ and kept the difference  $\omega_b - \omega_a = 5\Omega_{R}$ constant, which can be achieved either by varying $W/\hbar$ while keeping constant $\omega_{a,b}$ or by varying $\omega_a$ and $\omega_b$ at the same rate while keeping $W/\hbar$ constant. The anticrossings are clearly seen in the plot of eigenfrequencies, when either $\omega_a$ or $\omega_b$ is resonant with the optical transition of an atom. As compared to Fig.~2, Fig.~\ref{Fig:eigen_2photon}(b) shows more possibilities for switching between the three product states as the detuning is swept through the two resonances at the rate slower than the Rabi frequencies and the generation of both bipartite and tripartite entangled states in the vicinity of resonances if the sweeping rate is comparable to the Rabi frequencies.

\begin{figure}[htb]
	\centering
	\begin{subfigure}{0.45\textwidth}
		\centering
		\includegraphics[width=\linewidth]{fig3a.pdf}
	\end{subfigure}
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\linewidth]{fig3b.pdf}
	\end{subfigure}
	\caption{ The eigenstates of the system described by Eqs.~\eqref{c001 dot}, \eqref{c100 dot} and \eqref{c010 dot} as a function of frequency detuning defined as $\omega_a-W/\hbar$, whereas the difference of modal frequencies $\omega_b - \omega_a = 5\Omega_{R}$ is kept constant. The eigenfrequencies are shown in (a), and the amplitudes of the eigenstates are shown in (b), in which the amplitudes of $C_{001}$, $C_{100}$ and $C_{010}$ are represented by the solid, dashed and dotted lines, respectively. The eigenfrequencies are shifted by $\left.\left( \frac{1}{2}\omega_a + \frac{1}{2}\omega_b + \frac{W}{\hbar} \right)\right|_{\omega_a=W/\hbar}$. }
	\label{Fig:eigen_2photon}
\end{figure}

Since the functions $\omega_{a,b}(t)$ and $\Omega_{a,b}$ are periodic with period $2 \pi/\Omega$, we can use the expansion (\ref{fourier}) in Eqs.~(\ref{g0 dot})-(\ref{gb dot}).

If we keep only the resonant terms, assuming for example the following
resonances, $\bar{\omega}_{a}=\frac{W}{\hbar }$ and $\bar{\omega}%
_{b}+m\Omega =\frac{W}{\hbar }$, where $m$ is the number of a particular Fourier harmonic, the equations get simplified,
\begin{equation}
\frac{d}{dt}\left(
\begin{array}{c}
G_{0} \\
G_{a} \\
G_{b}%
\end{array}%
\right) +\left(
\begin{array}{ccc}
0 & -i R_{a0} &
-i R_{bm} \\
-iR_{a0}^{\ast } & 0 & 0 \\
-i R_{bm}^{\ast} & 0 & 0
\end{array}%
\right) \left(
\begin{array}{c}
G_{0} \\
G_{a} \\
G_{b}%
\end{array}%
\right) =0.  \label{simp eq}
\end{equation}%

Other (nonresonant) harmonics can be neglected only if $\Omega_{Ra,b} \ll \Omega$, see \cite{tokman2020}. When the modulation amplitude is zero, $R_{a0} = \Omega_{Ra}$ and $R_{bm} = 0$. In this case one of the eigenvalues $\Gamma_0$ corresponds to the decoupled state  $\left\vert 0_{a}\right\rangle \left\vert 1_{b}\right\rangle \left\vert
0\right\rangle$. Two other eigenvalues $\Gamma_{1,2}$ describe the solution with Rabi oscillations between states $\left\vert 0_{a}\right\rangle \left\vert 0_{b}\right\rangle \left\vert
1\right\rangle$ and $\left\vert 1_{a}\right\rangle \left\vert 0_{b}\right\rangle \left\vert
0\right\rangle$.This is an obvious limit since frequency $\omega_a$ is in resonance with the transition frequency, whereas  $\omega_b$ is out of resonance.

Assuming a sinusoidal
modulation of the partial frequencies of both cavity modes as an example,
\begin{equation}
\omega _{a,b}\left( t\right) =\bar{\omega}_{a,b}-\Delta \omega _{a,b}\cdot
\sin \left( \Omega t\right) ,
\label{two mode sin modul of freq}
\end{equation}%
 and using the well-known expansion in series of the
harmonics of the modulation frequency $\Omega $, with coefficients expressed
in terms of Bessel functions,
\begin{equation}
e^{-i\frac{\Delta \omega }{\Omega }\cos \left( \Omega t\right) }=J_{0}\left(
\frac{\Delta \omega }{\Omega }\right) +2\sum_{n=1}^{\infty }\left( -i\right)
^{n}J_{n}\left( \frac{\Delta \omega }{\Omega }\right) \cos \left( n\Omega
t\right),  \label{coe}
\end{equation}
we can express Fourier amplitudes in Eq.~(\ref{simp eq}) through Bessel functions:
\begin{equation}
R_{a0} = \overline{\Omega _{Ra}} J_0\left(\frac{\Delta \omega_a}{\Omega}\right), \; R_{bm} = (-i)^{|m|}  \overline{\Omega _{Rb}} J_{|m|}\left(\frac{\Delta \omega_b}{\Omega}\right).
    \label{amplitudes2}
\end{equation}

Note that the modulation amplitudes in Eq.~(\ref{two mode sin modul of freq}) \textit{can be of
the order of the modulation frequency}, $\frac{\Delta \omega _{a,b}}{\Omega }%
\sim 1$, despite the requirement $\Delta \omega _{a,b}\ll $ $\bar{\omega}%
_{a,b}$.

As usual, to solve Eq.~(\ref{simp eq}) one has to find the
eigenvalues $\Gamma _{0,1,2}$ and eigenvectors of the matrix of
coefficients. The characteristic equation for the eigenvalues is $\Gamma
\left( \Gamma ^{2}+\Omega _{R\Sigma }^{2}\right) =0$, where the cumulative Rabi frequency is
\begin{equation}
\label{cum-rabi}
\Omega
_{R\Sigma }=\sqrt{|R_{a0}|^2 + |R_{bm}|^2 }.
\end{equation}
 The result is
\begin{equation}
\Gamma _{0}=0,\ \Gamma _{1,2}=\pm i\Omega _{R\Sigma }.  \label{eigenvalues}
\end{equation}

Figure 4 shows one example of the cumulative Rabi frequency $\Omega
_{R\Sigma }$ as a function of $\Delta \omega  = \Delta \omega _{a} = \Delta \omega _{b}$ for $m = 1$ and $\Omega _{Ra} = \Omega _{Rb}$. As expected, $\Omega_{R\Sigma }$ decays with detuning from resonances but the decay is nonmonotonic and depends on the order of harmonic resonances.

\begin{figure}[htb]
	\includegraphics[width=0.5\textwidth]{fig4.pdf}
\caption { Cumulative Rabi frequency $\Omega
_{R\Sigma }$ as a function of $\Delta \omega  = \Delta \omega _{a} = \Delta \omega _{b}$ for $m = 1$ and $\Omega _{R} = \Omega _{Ra} = \Omega _{Rb}$.}
\end{figure}

The eigenvalue $\Gamma _{0}$ (i.e. the solution $\propto e^{-\Gamma _{0}t}$
) corresponds to the eigenvector $\left(
\begin{array}{c}
0 \\
1 \\
- \frac{R_{a0}}{R_{bm}}
\end{array}%
\right) $, whereas eigenvalues $\Gamma _{1,2}$ (i.e. the solution behaving
as $\propto e^{-\Gamma _{1,2}t}$) correspond to the eigenvectors $\left(
\begin{array}{c}
\frac{\pm \Omega _{R\Sigma }}{R_{a0}^{\ast }} \\
1 \\
\frac{R_{bm}^{\ast}}{R_{a0}^{\ast}}.
\end{array}%
\right) $. Here the eigenvectors are not normalized to $1$. The resulting
solution is
\begin{equation}
\left(
\begin{array}{c}
G_{0} \\
G_{a} \\
G_{b}%
\end{array}%
\right) =A\left(
\begin{array}{c}
0 \\
1 \\
- \frac{R_{a0}}{R_{bm}} %
\end{array}%
\right) +B e^{-i\Omega _{R\Sigma }t}  \left(
\begin{array}{c}
\frac{ \Omega _{R\Sigma }}{R_{a0}^{\ast }} \\
1 \\
\frac{R_{bm}^{\ast}}{R_{a0}^{\ast}}.
\end{array}%
\right)
+C e^{i\Omega _{R\Sigma }t} \left(
\begin{array}{c}
\frac{- \Omega _{R\Sigma }}{R_{a0}^{\ast }} \\
1 \\
\frac{R_{bm}^{\ast}}{R_{a0}^{\ast}}.
\end{array}%
\right) ,  \label{sol of simp eq}
\end{equation}%
where the constants $A$, $B$ and $C$ are determined by the initial
conditions.

For an arbitrary initial state vector
\begin{equation}
\Psi =C_{000}(0)\left\vert 0_{a}\right\rangle \left\vert 0_{b}\right\rangle
\left\vert 0\right\rangle +C_{001}(0)\left\vert 0_{a}\right\rangle \left\vert
0_{b}\right\rangle \left\vert 1\right\rangle +C_{100}(0)\left\vert
1_{a}\right\rangle \left\vert 0_{b}\right\rangle \left\vert 0\right\rangle
+C_{010}(0)\left\vert 0_{a}\right\rangle \left\vert 1_{b}\right\rangle
\left\vert 0\right\rangle,
   \label{WF-init}
\end{equation}
satisfying the normalization condition
$$
|C_{000}(0)|^2 + |C_{001}(0)|^2 + |C_{100}(0)|^2 + |C_{010}(0)|^2 = 1,
$$
the constants in Eq.~(\ref{sol of simp eq}) are
\begin{align}
    & A = \frac{C_{100}(0) \frac{|R_{bm}|^2}{|R_{a0}|^2} - C_{010}(0) \frac{R_{bm}}{R_{a0}} }{1 + \frac{|R_{bm}|^2}{|R_{a0}|^2}}, \nonumber \\
    & B = \frac{1}{2} \left( \frac{C_{100}(0) + C_{010}(0) \frac{R_{bm}}{R_{a0}} }{1 + \frac{|R_{bm}|^2}{|R_{a0}|^2}} + C_{001}(0) \frac{R_{a0}^{\ast}}{\Omega_{R\Sigma}} \right), \nonumber \\
    & C = \frac{1}{2} \left( \frac{C_{100}(0) + C_{010}(0) \frac{R_{bm}}{R_{a0}} }{1 + \frac{|R_{bm}|^2}{|R_{a0}|^2}} - C_{001}(0) \frac{R_{a0}^{\ast}}{\Omega_{R\Sigma}} \right).
 \label{abc}
    \end{align}

Let's consider some examples of the initial conditions to
illustrate this solution.

\subsection{ An atom is excited; both modes are in the vacuum state:}

The initial state vector is $\Psi(0)
= \left\vert 0_{a}\right\rangle \left\vert 0_{b}\right\rangle
\left\vert 1\right\rangle$. In this case Eq.~(\ref{abc}) gives $A=0,B=-C=\frac{R_{a0}^{\ast} }{2\Omega
_{R\Sigma }}$. The full expression for the state vector at any moment of time becomes
\begin{eqnarray}
\Psi &=& e^{-i\int_{0}^{t}\omega _{001}\left( \tau \right) d\tau
}\cos \left( \Omega _{R\Sigma }t\right) \left\vert 0_{a}\right\rangle
\left\vert 0_{b}\right\rangle \left\vert 1\right\rangle -i\frac{R_{a0}^{\ast }}{\Omega
_{R\Sigma }} e^{-i\int_{0}^{t}\omega _{100}\left( \tau \right) d\tau
}\sin \left( \Omega _{R\Sigma }t\right) \left\vert 1_{a}\right\rangle
\left\vert 0_{b}\right\rangle \left\vert 0\right\rangle  \notag \\
&&- i \frac{R_{bm}^{\ast } }{\Omega _{R\Sigma }}e^{-i\int_{0}^{t}\omega _{010%
}\left( \tau \right) d\tau }\sin \left( \Omega _{R\Sigma }t\right)
\left\vert 0_{a}\right\rangle \left\vert 1_{b}\right\rangle \left\vert
0\right\rangle,  \label{sv for excited atom}
\end{eqnarray}
where
\begin{equation*}
\omega _{001}\left( t\right) =\frac{1}{2}\omega _{a}\left( t\right) +\frac{1%
}{2}\omega _{b}\left( t\right) +\frac{W}{\hbar },\ \ \omega _{100}\left(
t\right) =\frac{3}{2}\omega _{a}\left( t\right) +\frac{1}{2}\omega
_{b}\left( t\right) ,\ \ \ \ \omega _{010}\left( t\right) =\frac{1}{2}\omega
_{a}\left( t\right) +\frac{3}{2}\omega _{b}\left( t\right) .
\end{equation*}%

As we see, an initial atomic excitation decays into a pair of
electromagnetic modes. Their frequencies are modulated due to the modulation
of the cavity geometry and are split by the cumulative Rabi frequency.
In the absence of dissipation the excitation energy oscillates back and forth between an atom and the field modes at the cumulative Rabi frequency.

\subsection{Both cavity modes are excited; the atom is in the ground state:}

The initial state vector is $\Psi(0) =C_{100}(0)\left\vert
1_{a}\right\rangle \left\vert 0_{b}\right\rangle \left\vert 0\right\rangle
+C_{010}(0) \left\vert 0_{a}\right\rangle \left\vert 1_{b}\right\rangle
\left\vert 0\right\rangle$.
In this case the state vector is
\begin{eqnarray}
\Psi &=& - 2i B \frac{\Omega_{R\Sigma}}{R_{a0}^{\ast}} e^{-i\int_{0}^{t}\omega _{001}\left( \tau \right) d\tau
}\sin \left( \Omega _{R\Sigma }t\right) \left\vert 0_{a}\right\rangle
\left\vert 0_{b}\right\rangle \left\vert 1\right\rangle + \left( A + 2B\cos\left( \Omega _{R\Sigma }t\right) \right) e^{-i\int_{0}^{t}\omega _{100}\left( \tau \right) d\tau
} \left\vert 1_{a}\right\rangle
\left\vert 0_{b}\right\rangle \left\vert 0\right\rangle  \notag \\
&& + \left( 2 B \frac{R_{bm}^{\ast}}{R_{a0}^{\ast}}\cos \left( \Omega _{R\Sigma }t\right) -A\frac{R_{a0}}{R_{bm}} \right)  e^{-i\int_{0}^{t}\omega _{010
}\left( \tau \right) d\tau }
\left\vert 0_{a}\right\rangle \left\vert 1_{b}\right\rangle \left\vert
0\right\rangle,  \label{sv-excited-modes}
\end{eqnarray}
where
\begin{align}
    & A = \frac{C_{100}(0) \frac{|R_{bm}|^2}{|R_{a0}|^2} - C_{010}(0) \frac{R_{bm}}{R_{a0}} }{1 + \frac{|R_{bm}|^2}{|R_{a0}|^2}}, \;
    B = C = \frac{1}{2}  \frac{C_{100}(0) + C_{010}(0) \frac{R_{bm}}{R_{a0}} }{1 + \frac{|R_{bm}|^2}{|R_{a0}|^2}}.
 \label{abc2}
    \end{align}

\begin{figure}[htb]
\centering

\begin{subfigure}[a]{0.4\textwidth}
\includegraphics[width=1\linewidth]{fig5a.pdf}
\caption{ }
\end{subfigure}

\begin{subfigure}[b]{0.4\textwidth}
\includegraphics[width=1\linewidth]{fig5b.pdf}
\caption{}
\end{subfigure}

\begin{subfigure}[c]{0.4\textwidth}
\includegraphics[width=1\linewidth]{fig5c.pdf}
\caption{}
\end{subfigure}

\caption {  (a) The average normalized energy of an atom, (b) the number of quanta in mode $a$, and (c) the number of quanta in mode $b$ as a function of normalized time. The initial conditions are $C_{000}(0) = 0$, $C_{001}(0) = 0$, $C_{100}(0) = 1/2$, and $C_{001}(0) = \sqrt{3}/2$; i.e., the two modes are initially excited with different amplitudes whereas an atom is in the ground state. Other parameters are $ \Delta \omega _{a} = \Delta \omega _{b} = \Omega$, $m = 1$, and $\Omega _{Ra} = \Omega _{Rb} =\Omega _{R}$. }
\label{fig5}
\end{figure}

An atom, originally in its ground state, will get excited through resonant coupling to the EM field, as is obvious from Eq.~(\ref{sv-excited-modes}). The resulting dynamics of the averaged normalized energy of an atom $\left\langle \Psi \right\vert \hat{H}_{a}\left\vert \Psi \right\rangle /W$ and the numbers of quanta in mode $a$, $|C_{100}(t)|^2$ and mode $b$, $|C_{100}(t)|^2$  is shown in Fig.~\ref{fig5} for one generic set of initial conditions. Due to the presence of three coupled degrees of freedom, the evolution is more complicated than single-sinusoidal Rabi oscillations. Moreover, there is one particular choice of initial conditions,   $C_{010}(0) = - C_{100}(0) \frac{R_{a0}}{R_{bm}}$, which corresponds to $B = C = 0$ and $A = C_{100}(0)$, where the normalization condition gives $|C_{100}(0)|^2 = \left( 1 + \frac{|R_{a0}|^2}{|R_{bm}|^2} \right)^{-1}$. This gives the following state vector,
\begin{equation}
\left(
\begin{array}{c}
C_{001} \\
C_{100} \\
C_{010}%
\end{array}%
\right) = C_{100}(0) \left(
\begin{array}{c}
0 \\
e^{-i\int_{0}^{t}\omega _{100}\left( \tau \right) d\tau } \\
-\frac{R_{a0}}{R_{bm}}
e^{-i\int_{0}^{t}\omega _{010}\left( \tau \right) d\tau }%
\end{array}%
\right)   .  \label{coe for one excited cav mode}
\end{equation}%

It corresponds to
the solution in which an atom stays in the ground state and \textit{is not excited by the electromagnetic
field despite being in resonance}. It happens because of destructive
interference between two frequency-modulated electromagnetic modes. In this case the three quantities shown in Fig.~\ref{fig5} become constant in time, with the average atomic energy being zero at all times. This effect is discussed in more detail in Sec.~V where the dissipation is taken into account.

\section{Dynamics of two cavity modes coupled to a time-variable atom}

Consider now the situation in which the cavity is not changing with time
whereas the transition energy of an atom depends on the parameter $p$ which
is adiabatically modulated. For example, it could be an optical transition
in a semiconductor nanostructure under an applied time-variable bias. The
Hamiltonian of such an atom can be written as $\hat{H}=W\left( p\right) \hat{%
\sigma}^{\dagger }\hat{\sigma}$. The dynamics of an isolated atom conserves
the adiabatic invariant $\frac{\left\langle \Psi \right\vert \hat{H}%
\left\vert \Psi \right\rangle }{\omega \left( t\right) }$ , where $\omega
\left( t\right) =\frac{W\left( p\left( t\right) \right) }{\hbar }$.

The dipole moment of the transition is also modulated, $\left\langle
1\right\vert \mathbf{\hat{d}}\left\vert 0\right\rangle =\mathbf{d}\left(
p\left( t\right) \right) $, because atom wave functions in the coordinate
representation depend on the parameter $p$. We again consider small enough amplitude of modulation of the transition energy. In this case, using the arguments
similar to those in Sec.~II B, we can show that the dependence $\mathbf{d}(t) $ can be neglected; it is the dependence $W(t) $
which is important for the evolution of a coupled atom-field system. The
RWA Hamiltonian which describes such a system is
\begin{equation}
\hat{H}=\hbar \omega _{a}\left( \hat{a}^{\dagger }\hat{a}+\frac{1}{2}%
\right) +\hbar \omega _{b}\left( \hat{b}^{\dagger }\hat{b}+\frac{1}{2%
}\right) +W\left( t\right) \hat{\sigma}^{\dagger }\hat{\sigma}-\left[ \hat{%
\sigma}^{\dagger }\left( \chi _{a}\hat{a}+\chi _{b}\hat{b}\right) +\hat{%
\sigma}\left( \chi _{a}^{\ast }\hat{a}^{\dagger }+\chi _{b}^{\ast }\hat{b}%
^{\dagger }\right) \right] .  \label{RWA H}
\end{equation}
Consider again a sinusoidal modulation of the transition energy,

\begin{equation}
W\left( t\right) =\bar{W}-\hbar \Delta \omega \sin \left( \Omega t\right) ,
\label{Wt}
\end{equation}%
where $\frac{\bar{W}}{\hbar }\gg \Delta \omega $.

The Schr\"{o}dinger equation with this Hamiltonian allows analytic
solutions. For simplicity, we again consider the basis states with lowest
energies: $\left\vert 0_{a}\right\rangle \left\vert 0_{b}\right\rangle
\left\vert 0\right\rangle$, $\left\vert 0_{a}\right\rangle \left\vert
0_{b}\right\rangle \left\vert 1\right\rangle$, $\left\vert 1_{a}\right\rangle
\left\vert 0_{b}\right\rangle \left\vert 0\right\rangle$, and $\left\vert
0_{a}\right\rangle \left\vert 1_{b}\right\rangle \left\vert 0\right\rangle$.
The corresponding wave function is
\begin{equation}
\Psi =C_{000}\left\vert 0_{a}\right\rangle \left\vert 0_{b}\right\rangle
\left\vert 0\right\rangle +C_{001}\left\vert 0_{a}\right\rangle \left\vert
0_{b}\right\rangle \left\vert 1\right\rangle +C_{100}\left\vert
1_{a}\right\rangle \left\vert 0_{b}\right\rangle \left\vert 0\right\rangle
+C_{010}\left\vert 0_{a}\right\rangle \left\vert 1_{b}\right\rangle
\left\vert 0\right\rangle,
   \label{WF}
\end{equation}%
where the coefficients obey the equations%
\begin{equation}
\dot{C}_{000}+i\frac{\omega _{a}+\omega _{b}}{2}C_{000}=0;  \label{C000 dot}
\end{equation}%
\begin{equation}
\dot{C}_{001}+i\left( \frac{1}{2}\omega _{a}+\frac{1}{2}\omega _{b}+\frac{%
W\left( t\right) }{\hbar }\right) C_{001}-i\Omega _{Ra}C_{100}-i\Omega
_{Rb}C_{010}=0,  \label{C001 dot}
\end{equation}%
\begin{equation}
\dot{C}_{100}+i\left( \frac{3}{2}\omega _{a}+\frac{1}{2}\omega _{b}\right)
C_{100}-i\Omega _{Ra}^{\ast }C_{001}=0,  \label{C100 dot}
\end{equation}%
\begin{equation}
\dot{C}_{010}+i\left( \frac{1}{2}\omega _{a}+\frac{3}{2}\omega _{b}\right)
C_{010}-i\Omega _{Rb}^{\ast }C_{001}=0.
  \label{C010 dot}
\end{equation}%
After the substitution%
\begin{equation}
\left(
\begin{array}{c}
C_{001} \\
C_{100} \\
C_{010}%
\end{array}%
\right) =\left(
\begin{array}{c}
G_{0}\exp \left[ -i\left( \frac{1}{2}\omega _{a}+\frac{1}{2}\omega
_{b}\right) t+\int_{0}^{t}\frac{W\left( \tau \right) }{\hbar }d\tau \right]
\\
G_{a}\exp \left[ -i\left( \frac{3}{2}\omega _{a}+\frac{1}{2}\omega
_{b}\right) t\right]  \\
G_{b}\exp \left[ -i\left( \frac{1}{2}\omega _{a}+\frac{3}{2}\omega
_{b}\right) t\right]
\end{array}%
\right) ,  \label{Substitution}
\end{equation}%
we obtain%
\begin{equation}
\dot{G}_{0}-i\Omega _{Ra}G_{a}\exp \left[ -i\left( \omega _{a}t-\int_{0}^{t}%
\frac{W\left( \tau \right) }{\hbar }d\tau \right) \right] -i\Omega
_{Rb}G_{b}\exp \left[ -i\left( \omega _{b}t-\int_{0}^{t}\frac{W\left( \tau
\right) }{\hbar }d\tau \right) \right] =0,  \label{G0 dot}
\end{equation}%
\begin{equation}
\dot{G}_{a}-i\Omega _{Ra}^{\ast }G_{0}\exp \left[ -i\left( \omega
_{a}t-\int_{0}^{t}\frac{W\left( \tau \right) }{\hbar }d\tau \right) \right]
=0,  \label{Ga dot}
\end{equation}%
\begin{equation}
\dot{G}_{b}-i\Omega _{Rb}^{\ast }G_{0}\exp \left[ -i\left( \omega
_{b}t-\int_{0}^{t}\frac{W\left( \tau \right) }{\hbar }d\tau \right) \right]
=0,  \label{Gb dot}
\end{equation}%
Similarly to the previous section, we expand the exponents in Eqs.~(\ref{G0
dot})-(\ref{Gb dot}) over the harmonics of the modulation frequency $\Omega
$ using Eq.~(\ref{coe}) and keep only the resonant terms, assuming for
definiteness that $\omega _{a}=\frac{\bar{W}}{\hbar }$ and $\omega
_{b}+m\Omega =\frac{\bar{W}}{\hbar }$. We again obtain Eq.~(\ref{simp eq}), where $R_{a0} = \Omega_{Ra} J_0\left(\frac{\Delta \omega}{\Omega}\right)$, $R_{bm} = (-i)^{|m|}  \Omega _{Rb} J_{|m|}\left(\frac{\Delta \omega}{\Omega}\right)$.
Therefore, the
modulation of the atomic transition and the cavity parameters leads to a
similar dynamics.

\section{Dynamics of open time-dependent cavity QED systems}

\subsection{The stochastic evolution of the state vector}

Consider again the dynamics of two adiabatically varying cavity modes
coupled to an atom, but this time we include the processes of relaxation and decoherence in an open
system, which is (weakly) coupled to a dissipative reservoir. We will use the approach based on the
stochastic evolution of the state vector; see Appendix B and
\cite{tokman2020}.
This is basically the Schr\"{o}dinger equation modified by
adding a linear relaxation operator and the noise source term with
appropriate correlation properties. The latter are related to the parameters
of the relaxation operator, which is a manifestation of the
fluctuation-dissipation theorem \cite{Landau1965}. In Appendix B we outlined the main properties of the stochastic equation of evolution
and showed how physically reasonable constraints on the observables
determine the properties of the noise sources. We also demonstrated the
relationship between our approach and the Lindblad method of solving the
master equation.

Within our approach the system is described by a state vector which has a
fluctuating component: $\left\vert \Psi \right\rangle =\overline{\left\vert
\Psi \right\rangle }+\widetilde{\left\vert \Psi \right\rangle }$, where the
straight bar means averaging over the statistics of noise and the wavy bar
denotes the fluctuating component. This state vector is of course very
different from the state vector obtained by solving a standard Schr\"{o}%
dinger equation for a closed system. In fact, coupling to a dissipative
reservoir leads to the formation of a mixed state, which can be described by
a density matrix $\hat{\rho}=\overline{\left\vert \Psi \right\rangle }\cdot
\overline{\left\langle \Psi \right\vert }+\overline{\widetilde{\left\vert
\Psi \right\rangle }\widetilde{\left\langle \Psi \right\vert }}$. However, the density matrix equations are more cumbersome for the analytic solution as compared to the formalism used in this paper.

One can view the stochastic equation approach as a convenient formalism for
calculating physical observables which allows one to obtain analytic solutions for the evolution of a coupled system in the presence of dissipation and decoherence. When the Markov approximation is applied, the results are equivalent to those obtained within the Lindblad master equation formalism. Within the Markov approximation, the relaxation operator in the stochastic equation for the state vector is obtained simply by summing up partial Lindbladians for all subsystems, whatever they are (in our case these are a fermion emitter and two EM cavity modes). Then the noise source term is determined unambiguously by conservation of the norm of the state vector and the requirement that the system should approach thermal equilibrium when the external perturbation is turned off. This immediately gives Eqs.~(\ref{modified c000 dot})-(\ref{modified c010 dot}) below.

Following the derivation in Appendix B, equations~(\ref{c000 dot})-(\ref{c010 dot}) are modified due
to the terms with relaxation constants $\gamma _{000}$,$\gamma _{001}$,$\gamma _{010}$, and $\gamma _{100}$ which are originated from the Lindladians,  and
the noise sources,
\begin{equation}
\left( \frac{\partial }{\partial t}+\gamma _{000}\right) C_{000}+i\frac{%
\omega _{a}\left( t\right) +\omega _{b}\left( t\right) }{2}C_{000}=-\frac{i}{%
\hbar }\mathfrak{R}_{000};  \label{modified c000 dot}
\end{equation}%
\begin{equation}
\left( \frac{\partial }{\partial t}+\gamma _{001}\right) C_{001}+i\left(
\frac{1}{2}\omega _{a}\left( t\right) +\frac{1}{2}\omega _{b}\left( t\right)
+\frac{W}{\hbar }\right) C_{001}-i\Omega _{Ra}C_{100}-i\Omega _{Rb}C_{010}=-%
\frac{i}{\hbar }\mathfrak{R}_{001},  \label{modified c001 dot}
\end{equation}%
\begin{equation}
\left( \frac{\partial }{\partial t}+\gamma _{100}\right) C_{100}+i\left(
\frac{3}{2}\omega _{a}\left( t\right) +\frac{1}{2}\omega _{b}\left( t\right)
\right) C_{100}-i\Omega _{Ra}^{\ast }C_{001}=-\frac{i}{\hbar }\mathfrak{R}_{100},
\label{modified c100 dot}
\end{equation}%
\begin{equation}
\left( \frac{\partial }{\partial t}+\gamma _{010}\right) C_{010}+i\left(
\frac{1}{2}\omega _{a}\left( t\right) +\frac{3}{2}\omega _{b}\left( t\right)
\right) C_{010}-i\Omega _{Rb}^{\ast }C_{001}=-\frac{i}{\hbar }\mathfrak{R}_{010}.
\label{modified c010 dot}
\end{equation}%

We assume that  noise terms in Eq.~(\ref{modified c000 dot})-(\ref{modified c010 dot}) become equal to zero after averaging over the noise statistics. The averages of the
quadratic combinations of noise source terms are nonzero and we assume here that they are delta-correlated in time (the Markov approximation),
\begin{equation}
\overline{\mathfrak{R}_{\beta }^{\ast }\left( t+\xi \right) \mathfrak{R}_{\alpha }\left( t\right) }%
=\overline{\mathfrak{R}_{\beta }^{\ast }\left( t\right) \mathfrak{R}_{\alpha }\left( t+\xi \right)
}=\hbar ^{2}\delta \left( \xi \right) D_{\alpha \beta }.  \label{correlator}
\end{equation}%
Here the indices $\alpha$ and $\beta$ span a set of the lowest-energy states $\left\vert 0_{a}\right\rangle \left\vert 0_{b}\right\rangle
\left\vert 0\right\rangle$, $\left\vert 0_{a}\right\rangle \left\vert
0_{b}\right\rangle \left\vert 1\right\rangle$, $\left\vert 1_{a}\right\rangle
\left\vert 0_{b}\right\rangle \left\vert 0\right\rangle$, and $\left\vert
0_{a}\right\rangle \left\vert 1_{b}\right\rangle \left\vert 0\right\rangle$.  Including the
noise sources is crucial for consistency of the formalism: it ensures the
conservation of the norm of the state vector and leads to a physically
meaningful equilibrium state.

Consider the case of zero temperatures for all reservoirs, which means in
practice that these temperatures in energy units are much lower than the
atomic transition energy and the cavity mode frequencies. In this case the relaxation constants are greatly simplified as compared to the general expressions given in Appendix B,
\begin{equation}
\gamma_{000} = 0, \, \gamma_{001} = \frac{\gamma }{2}, \,
\gamma _{100}=\frac{\mu _{a}}{2}, \,
\gamma _{010}= \frac{\mu
_{b}}{2},  \label{gamma at 0 temp}
\end{equation}
where $\gamma $ is the inelastic relaxation rate for an isolated atom, $\mu _{a,b}$
are relaxation rates of the EM modes determined by the cavity Q-factor; these ``partial'' relaxation constants are determined by couplings to their respective dissipative reservoirs. Appendix B outlines how to include elastic decoherence processes.

In this limit we can drop the noise terms in the right-hand side of all
equations for the components of the state vector, except the term $\mathfrak{R}_{000}$
in the equation for $C_{000}$; see Appendix B. This noise term ensures
conservation of the norm,
\begin{equation*}
 \overline{
\left\vert C_{000}\right\vert ^{2}}+\overline{\left\vert
C_{001}\right\vert ^{2}}+ \overline{\left\vert
C_{010}\right\vert ^{2}}+ \overline{\left\vert
C_{100}\right\vert ^{2}} =0,
\end{equation*}%
if its correlator is given by
\begin{equation*}
\overline{\mathfrak{R}_{000}\left( t+\xi \right) \mathfrak{R}_{000}^{\ast }\left( t\right) }%
= 2\hbar ^{2}\delta \left( \xi \right) \left( \gamma _{100}\overline{\left\vert C_{100}\right\vert
^{2}} +\gamma _{001}\overline{
\left\vert C_{001}\right\vert ^{2}} +\gamma _{010}\overline{
\left\vert C_{010}\right\vert ^{2}} \right).
\end{equation*}

As an example, consider a high-quality cavity and neglect the cavity losses as compared to the atomic decay.
In this case, and for a low temperature of an atomic reservoir, Eqs.~(\ref%
{modified c001 dot})-(\ref{modified c010 dot}) take the form%
\begin{equation}
\left( \frac{\partial }{\partial t}+\frac{\gamma }{2}\right) C_{001}+i\left(
\frac{1}{2}\omega _{a}\left( t\right) +\frac{1}{2}\omega _{b}\left( t\right)
+\frac{W}{\hbar }\right) C_{001}-i\Omega _{Ra}C_{100}-i\Omega _{Rb}C_{010}=0,
\label{modified c001 dot at low temp}
\end{equation}%
\begin{equation}
\frac{\partial }{\partial t}C_{100}+i\left( \frac{3}{2}\omega _{a}\left(
t\right) +\frac{1}{2}\omega _{b}\left( t\right) \right) C_{100}-i\Omega
_{Ra}^{\ast }C_{001}=0,  \label{modified c100 dot at low temp}
\end{equation}%
\begin{equation}
\frac{\partial }{\partial t}C_{010}+i\left( \frac{1}{2}\omega _{a}\left(
t\right) +\frac{3}{2}\omega _{b}\left( t\right) \right) C_{010}-i\Omega
_{Rb}^{\ast }C_{001}=0.  \label{modified c010 dot at low temp}
\end{equation}%
Using the substitution of variables in Eq.~(\ref{substitution}) and
repeating the same derivation as in Sec.~III, we arrive at
\begin{equation}
\frac{d}{dt}\left(
\begin{array}{c}
G_{0} \\
G_{a} \\
G_{b}%
\end{array}%
\right) +\left(
\begin{array}{ccc}
\frac{\gamma}{2} & -i R_{a0} &
-i R_{bm} \\
-iR_{a0}^{\ast } & 0 & 0 \\
-i R_{bm}^{\ast} & 0 & 0
\end{array}%
\right) \left(
\begin{array}{c}
G_{0} \\
G_{a} \\
G_{b}%
\end{array}%
\right) =0.
\label{modified simp eq}
\end{equation}%
Its solution is determined by the eigenvalues and eigenvectors of the matrix
in Eq.~(\ref{modified simp eq}). The eigenvalues are given by%
\begin{equation*}
\Gamma \left[ \left( \Gamma -\frac{\gamma }{2}\right) \Gamma +\Omega
_{R\Sigma }^{2}\right] =0,
\end{equation*}%
which yields%
\begin{equation}
\Gamma _{0}=0,\ \Gamma _{1,2}=\frac{\gamma }{4}\pm i\sqrt{\Omega _{R\Sigma
}^{2}-\frac{\gamma ^{2}}{16}}.  \label{modified eigenvalues}
\end{equation}%
The eigenvector corresponding to the eigenvalue $\Gamma _{0}=0$ is the same
as in the absence of dissipation (see Sec.~III B), whereas the expressions for
the eigenvectors corresponding to eigenvalues $\Gamma _{1,2}$ can be
obtained from \textquotedblleft dissipationless\textquotedblright\
expressions by replacing $\pm $ $\Omega _{R\Sigma }\longrightarrow \pm \sqrt{%
\Omega _{R\Sigma }^{2}-\frac{\gamma ^{2}}{16}}-i\frac{\gamma }{4}$. As a
result, we obtain the following expression for the state vector,%
\begin{eqnarray}
\left(
\begin{array}{c}
C_{001} \\
C_{100} \\
C_{010}%
\end{array}%
\right) &=&A\left(
\begin{array}{c}
0 \\
e^{-i\int_{0}^{t}\omega _{100}\left( \tau \right) d\tau } \\
- \frac{R_{a0}}{R_{bm}}
e^{-i\int_{0}^{t}\omega _{010}\left( \tau \right) d\tau }%
\end{array}%
\right)  \notag \\
+ &&B e^{\left( -i\sqrt{\Omega _{R\Sigma }^{2}-\frac{\gamma ^{2}}{16}}-\frac{%
\gamma }{4}\right) t}\left(
\begin{array}{c}
\frac{\sqrt{\Omega _{R\Sigma }^{2}-\frac{\gamma ^{2}}{16}}-i\frac{\gamma }{4}%
}{R_{a0}^{\ast} }%
e^{-i\int_{0}^{t}\omega _{001}\left( \tau \right) d\tau } \\
e^{-i\int_{0}^{t}\omega _{100}\left( \tau \right) d\tau } \\
 \frac{R_{bm}^{\ast}}{R_{a0}^{\ast}}  e^{-i\int_{0}^{t}\omega _{010}\left( \tau \right) d\tau }%
\end{array}%
\right)  \notag \\
&& + C e^{\left( i\sqrt{\Omega _{R\Sigma }^{2}-\frac{\gamma ^{2}}{16}}-\frac{%
\gamma }{4}\right) t}\left(
\begin{array}{c}
\frac{-\sqrt{\Omega _{R\Sigma }^{2}-\frac{\gamma ^{2}}{16}}-i\frac{\gamma }{4%
}}{R_{a0}^{\ast}}  e^{-i\int_{0}^{t}\omega _{001}\left( \tau \right) d\tau } \\
e^{-i\int_{0}^{t}\omega _{100}\left( \tau \right) d\tau } \\
 \frac{R_{bm}^{\ast}}{R_{a0}^{\ast}}  e^{-i\int_{0}^{t}\omega _{010}\left( \tau \right) d\tau }%
\end{array}%
\right) .  \label{sol for modified simp eq}
\end{eqnarray}%
Where the constants $A$, $B$ and $C$ are given by initial conditions. In the limit $\Omega_{R\Sigma} \gg \gamma$ their dependence on the initial values $C_{100}(0)$, $C_{010}(0)$, and $C_{001}(0)$ is given by Eqs.~(\ref{abc}) from the previous section,  whereas their dependence on $C_{000}(0)$ is determined by the normalization condition.

\subsection{Modulation-induced transparency}

Note
again the existence of the solution with $B = C = 0$ in which an atom initially in the ground state \textit{is decoupled
from the electromagnetic field} and stays in the ground state because of destructive interference between the EM modes.
There is however an interesting difference as compared to the dissipationless case discussed in Sec.~IIE. For arbitrary initial conditions, when $A,B,C$ are not equal to zero, part of the field energy will be resonantly transferred to the atom and dissipate through the atomic decay. However, the terms with $B$ and $C$ factors in Eq.~(\ref{sol for modified simp eq}) decay exponentially with time, and  the solution to Eq.~(\ref{sol for modified simp eq}) at $ t \gg 1/\gamma$ will acquire the same form as  in the case of $B = C = 0$:
\begin{align}
\Psi & =A \left( e^{-i\int_{0}^{t}\omega _{100}\left( \tau \right)
d\tau }\left\vert 1_{a}\right\rangle \left\vert 0_{b}\right\rangle
\left\vert 0\right\rangle  - \frac{R_{a0}}{R_{bm}}   e^{-i\int_{0}^{t}\omega _{010%
}\left( \tau \right) d\tau }\left\vert 0_{a}\right\rangle \left\vert
1_{b}\right\rangle \left\vert 0\right\rangle \right)  +C_{000}\left\vert
0_{a}\right\rangle \left\vert 0_{b}\right\rangle \left\vert 0\right\rangle
\label{modified wf}.
\end{align}

The value of $C_{000}$ at $ t \gg 1/\gamma$ is determined by the noise term $-\frac{i}{\hbar} \mathfrak{R}_{000}$ in the right-hand side of Eq.~(\ref{modified c000 dot}) and satisfies $\overline{C_{000}} = 0$, $\overline{|C_{000}|^2} = 1- |A|^2 \left[ 1 + \frac{|R_{a0}|^2}{|R_{bm}|^2} \right]$ (see Appendix B).

The value of $|A|^2$ is given by
\begin{equation}
    |A|^2 = \displaystyle \frac{1-|C_{000}(0)|^2 - |C_{001}(0)|^2}{1+|Z|^2} \left[ \frac{\left| \frac{|R_{bm}|^2}{|R_{a0}|^2} - Z \frac{R_{bm}}{R_{a0}} \right|^2 }{\left( 1+ \frac{|R_{bm}|^2}{|R_{a0}|^2} \right)^2 } \right],
\end{equation}
where $Z = \displaystyle \frac{C_{010}(0)}{C_{100}(0)}$.
The value of $|A|^2$ reaches a maximum when Arg$[Z] = \pi - {\rm Arg}\left[\frac{R_{bm}}{R_{a0}} \right]$ and $|Z| = \left| \frac{R_{a0}}{R_{bm}} \right|$, which corresponds to $C_{010}(0) = - C_{100}(0) \frac{R_{a0}}{R_{bm}}$ and
\begin{equation}
    \label{a2}
|A|^2 = \displaystyle \frac{1-|C_{000}(0)|^2 - |C_{001}(0)|^2}{1+ \frac{|R_{a0}|^2}{|R_{bm}|^2}}.
\end{equation}
This equation has a simple interpretation. According to Eq.~(\ref{modified wf}), the average steady-state number of quanta in both modes is
\begin{equation}
|C_{100}|^2 + |C_{010}|^2 = \displaystyle  |A|^2 \left(1+ \frac{|R_{a0}|^2}{|R_{bm}|^2} \right).
\label{avnumber}
\end{equation}
Comparing Eq.~(\ref{avnumber}) and Eq.~(\ref{a2}), one can see that despite the presence of dissipation, when the value of $|A|^2$ reaches a maximum given by Eq.~(\ref{a2}) the average steady-state number of field quanta given by Eq.~(\ref{avnumber}) is equal to its initial value:  $ |C_{100}(0)|^2 + |C_{010}(0)|^2 = 1-|C_{000}(0)|^2 - |C_{001}(0)|^2 $.

\begin{figure}[htb]
	\includegraphics[width=0.5\textwidth]{fig6.pdf}
\caption { The contour plot of the normalized average number of quanta $\overline{N_q}$ on the complex $Z$ plane for $m = 1$, $\Omega_{Ra} = \Omega_{Rb}$ and $\Delta \omega_a = \Delta \omega_b = \Omega$.}
\label{asquared}
\end{figure}

The contour plot of the  average steady-state number of quanta normalized by its initial value,
$$ \overline{N_q} = \displaystyle  \frac{|A|^2 \left(1+ \frac{|R_{a0}|^2}{|R_{bm}|^2} \right)}{ |C_{100}(0)|^2 + |C_{010}(0)|^2  }$$
on the complex $Z$ plane is shown in Fig.~\ref{asquared} for $m = 1$, $\Omega_{Ra} = \Omega_{Rb}$ and $\Delta \omega_a = \Delta \omega_b = \Omega$. For this particular choice of parameters, the maximum of  the number of quanta is reached at Arg$[Z] = -\pi/2$, i.e., it is located on the imaginary axis as shown in the figure. For $m = 2$ the maximum will be on the real axis. At its maximum, the average number of quanta is equal to its initial value, i.e. it remains constant.

For any initial conditions other than those corresponding to the maximum of $\overline{N_q}$, a part of the EM field energy will dissipate through interaction with an atom, and eventually only the part which corresponds to the combination of modes completely decoupled from an atom due to destructive interference survives. This will result in smaller values of $\overline{N_q}$.  Of course, eventually the finite cavity losses will kick in and the
field will dissipate to the level of quantum and thermal fluctuations.

Finally, for the initial state $\Psi
\left( 0\right) =\left\vert 0_{a}\right\rangle \left\vert 0_{b}\right\rangle
\left\vert 1\right\rangle $ (only the atom is excited) we have $A=0$, i.e.
the system goes into the ground state as expected.

To summarize, the modulated system of an atom resonantly coupled to two EM
cavity modes demonstrates an interesting effect of modulation-induced
transparency. In the absence of modulation, the presence of an atom
experiencing an incoherent decay leads to the dissipation of the EM field
 even if the empty cavity is ideal,
i.e. has zero losses. However, low-frequency modulation of the cavity or of
the transition frequency of an atom creates the EM field distribution which
is completely decoupled from an atom due to destructive interference between
the cavity modes, even at resonance between the atomic transition and the
cavity mode frequencies. Therefore, the atom will remain in the ground state and the field will experience no
dissipation in the absence of cavity losses. For a classical field, such a destructive interference effect which switches off the
field dissipation in resonant medium by introducing low-frequency modulation
was considered, in particular, in Ref.~\cite{rad2006} for acoustically modulated
two-level atoms. Similar effects in the interaction of classical fields with atoms are discussed in the introduction of Ref.~\cite{rad2020}.

\subsection{Prospects for strong coupling and quantum entanglement in various nanophotonic systems}

Expressions in this section and more general expressions for the relaxation rates in Appendix B (see, e.g., Eqs.~(\ref{gamma0}),(\ref{gamma1})) allow one to calculate the effective decoherence rates from the known  ``partial'' relaxation rates
for individual subsystems: EM cavity modes and any kind of a fermionic qubit. One can compare the decoherence rates with characteristic Rabi frequencies which enter the solution for the evolution equations such as Eqs.~(\ref{modified c001 dot at low temp})-(\ref{modified c010 dot at low temp}) in order to determine if the strong coupling regime and quantum entanglement in the electron-photon system can be achieved. For any specific application, one should also compare the effective relaxation times with relevant operation times (gate transition time, read/write time etc.) In the discussion below, we rely on the parameters obtained from Refs.~\cite{haroche}-\cite{aspelmeyer2014} which we already cited in the Introduction. Many of them are recent reviews and one can find further references there. We don't attempt to overview here the vast and rapidly growing amount of literature on the subject.

In electron-based quantum emitters the largest oscillator strengths in the visible/near-infrared range have been observed for excitons in organic molecules, followed by perovskites and more conventional inorganic semiconductor quantum dots. The typical variation of the dipole matrix element of the optical transition which enters the Rabi frequency is from tens of nm to a few Angstrom. The dipole moment grows with increasing wavelength. The relaxation times are strongly temperature and material quality-dependent, varying from tens or hundreds of ps for single quantum dots at 4 K to the $\mu$s range for defects in semiconductors and diamond at mK temperatures. At room temperature the typical decoherence rates for the optical transition are in the $\sim 10$ meV range.

The photon decay times are longest for dielectric micro- and nano-cavities: photonic crystal cavities, nanopillars, distributed Bragg reflector mirrors, microdisk whispering gallery mode cavities, etc. Their quality factors are typically between $10^3 - 10^7$, corresponding to photon lifetimes from sub-ns to $\mu$s range. However, the field localization in the dielectric cavities is diffraction-limited, which limits the attainable Rabi frequency values to hundreds of $\mu$eV. The effective decoherence rate in dielectric cavity QED systems is typically limited by the relaxation in the fermion quantum emitter subsystem,

In plasmonic cavities, field localization on a nm and even sub-nm scale has been achieved, but the photon losses are in the ps or even fs range and therefore, they dominate the overall decoherence rate. Still, when it comes to strong coupling at room temperature to a single quantum emitter such as a single molecule or a quantum dot, the approach utilizing plasmonic cavities has seen more success so far. In these systems the Rabi splitting of the order of 100-200 meV has been observed. In plasmonic systems it may be beneficial to consider longer-wavelength emitters with the optical transition at the mid-infrared and even terahertz wavelengths. Indeed, with increasing wavelength the plasmon losses go down, the matrix element of a dipole-allowed transition increases, whereas the plasmon localization stays largely the same.

Another factor that has to be taken into account when choosing a nanophotonic system for a specific application is the rate with which the modulation of the cavity or emitter parameters has to be performed. For example, if the modulation at the rate comparable to the Rabi frequency or operation with $\pi$- or $\pi/2$ pulses is required, the plasmonic-based systems run into a problem: they would require $\sim 10-100$ fs pulses for modulation, which obviously can be achieved only with fs lasers. All electronic operations typically have a cutoff at tens of GHz. Applications of nanophotonics to quantum computing are especially challenging, because  computations require at least 99.99\% fidelity, i.e. at least $10^4$ ``flops'' before decoherence kicks in.

\section{Conclusions}

In conclusion, we developed the analytic theory describing the dynamics and control  of strongly coupled nanophotonic systems with time-variable parameters. The coupling of the fermion and photon subsystems to their dissipative reservoirs are  included within the stochastic equation of evolution approach, which is equivalent to the Lindblad approximation in the master equation formalism. Our analytic solution is valid in the approximation that the rate of parameter modulation and the amplitude of the frequency modulation are much smaller than the optical transition frequencies. At the same time, they can be arbitrary with respect to the generalized Rabi oscillations frequency which determines the coherent dynamics. Therefore, we can describe an arbitrary modulation of the parameters, both slower and faster than the Rabi frequency, for complete control of the quantum state. For example, one can turn on and off the entanglement between the fermionic and photonic degrees of freedom, swap between the quantum states, or decouple  the fermionic qubit from the cavity field via modulation-induced transparency.

\begin{acknowledgments}
This work has been supported in part by the Air Force Office for Scientific Research
Grant No.~FA9550-17-1-0341, National Science Foundation Award No.~1936276, and Texas A\&M University through STRP, X-grant and T3-grant programs.
 M.E.~and M.T.~acknowledge the support from RFBR Grant No. 20-02-00100.

\end{acknowledgments}

redpajama_set_name: RedPajamaArXiv


text:
Cryptocurrency: checking the future of Bitcoin, Ethererum & Co
After the cryptocurrency market was able to stabilize slightly after the losses in the first half of 2022, Bitcoin fell to its lowest level in 3 months in mid-September. The second largest cryptocurrency, Ether, is also affected. Ethereum only switched to the more energy-efficient proof-of-stake mechanism in early September — an update that many investors had been hoping for. COMPUTER BILD took a close look at the top coins and explained which cryptocurrency has a future and which coin is mostly driven by hype.
Bitcoin: does cryptocurrency have a future?
The oldest, largest and most valuable cryptocurrency is the undisputed number one. Technically, however, we note the pioneering role of the Bitcoin protocol. With its proof-of-work mechanism to validate transactions on the blockchain, the currency has gained a reputation as a climate polluter in recent years – mining and sending Bitcoin is resource-intensive. There are now many other cryptocurrencies that rely on an efficient mechanism. In real-world application, Ethereum, Cardano, Solana, and other coins have overtaken Bitcoin.
However, the proof-of-work mechanism also has advocates, whose voices have been heard again in recent months. While the proof-of-stake mechanism ultimately centralizes cryptocurrencies in the hands of a few, proof-of-work enables true decentralization. Ethereum and other cryptocurrencies are issued by foundations or corporations, which typically hold a significant amount of coins in circulation, giving them de facto major influence over decisions affecting the future of cryptocurrency. The Bitcoin developer is anonymous as of this writing and does not appear publicly.
Whether Bitcoin has a future as a cryptocurrency or will be replaced by other coins depends on the role cryptocurrencies will play in a financial system of the future. If decentralization is really happening, Bitcoin is the only major cryptocurrency that truly meets this requirement. However, an integration of blockchain technology into the existing financial system is also possible – as Ethereum & Co. are striving for this with smart contracts. A cryptocurrency like Bitcoin would have no future here due to its slow architecture.
Trade stocks and ETFs for free on Finanzen.net Zero
Stablecoins: future of digital payments?
Bitcoin has largely said goodbye to the original idea of a decentralized means of payment. Instead, this role is increasingly fulfilled by stablecoins, which provide the stability and efficiency needed for financial transactions. But stablecoin is not just stablecoin. True stablecoins are backed by fiat currencies and therefore always offer a relatively stable exchange rate around one US dollar. Algorithmic stablecoins trick investors into directing the market to settle around $1. The collapse of stablecoin Terra has shown just how fragile this system is.
DeFi applications, i.e. cryptocurrencies used in decentralized financial products, will continue to change the financial landscape regardless of the fate of individual cryptocurrencies. stablecoins like attached or the USD Coins have a good chance of continuing to play an important role in the future of cryptocurrencies. Whether this trend will affect the future cryptocurrency already in 2022 and 2023 remains to be seen. Stablecoins are hardly interesting for investors – at least without the other currencies with some volatility.
Ripple: competition for the SWIFT system
Despite the ongoing digitization, classic SEPA transfers still often take a whole day. Ripple decided to change this situation. With the Ripple protocol, the developers want to establish a global decentralized digital transfer system in the future, in which every transaction is carried out securely and quickly. Although Ripple is counted among cryptocurrencies, XRP is not based on blockchain technology. Ripple processes transactions through independent validation servers operated by banks, institutions or individuals. The underlying protocol is released under an open source license, so banks can easily adapt it to their needs. However, Ripple is still under pressure right now, as the US Securities and Exchange Commission (SEC) wants to treat Ripple Coin as a security and subject it to the same regulatory conditions. The future of this cryptocurrency essentially depends on the outcome of the lawsuit. A verdict is expected for 2022 – but has been repeatedly postponed.
Simply copy top performing traders (including crypto trading) to eToro
*78% of private CFD accounts lose money
Smart Contracts: where are Ethereum, Cardano and Solana in 2030?
Smart contracts represent payment agreements on a blockchain and work like a contract: two parties agree to exchange a service for money at a specific time. If a party fails to fulfill its obligation, the laws or courts regulate the consequences. With smart contracts, contracts are automatically validated via the blockchain without third parties. Since the majority of contracts govern the simple exchange of goods/services and money – for example, leases, loan agreements or insurance policies – there are potentially significant efficiencies to be gained through to the use of smart contracts.
Close to Ethereum also represent other major cryptocurrencies such as Solana (SOL) Where Gimbal (ADA) focus on smart contracts. It has not yet been decided which cryptocurrency will prevail. More recently, Ethereum completed the long-awaited update to the proof-of-stake mechanism, which makes transactions faster and cheaper.
Trade shares at a lower cost (0 to 4 euros) with Smartbroker
Will there still be a cryptocurrency in 2030?
The countries of the world are only slowly catching up with technical developments and starting to take care of the crypto market. You currently approach the phenomenon very differently:
In Australia In mid-September, a parliamentarian tabled a bill to regulate crypto exchanges, stablecoins, and China's central bank cryptocurrency. Adoption is not expected before the end of 2022 at the earliest. The draft will first be submitted to Parliament for deliberation in early November.
China completely banned cryptocurrency trading for locals in 2021. The background is the development of its own cryptocurrency – the digital yuan. By centralizing digital money flows, the Chinese government would not only have full visibility of all transactions, but also control. How the global crypto market will deal with a digital yuan is unclear. In any case, the Bitcoin & Co. ban caused strong setbacks.
Out of Russia In late 2021/early 2022, there were also considerations for banning cryptocurrencies. A law has been passed to regulate it. Russian crypto companies now pay a 13% tax.
Also from the United States comes a piece of news that may be game-changing for the future of cryptocurrency in 2030. The US Federal Reserve has raised interest rates to 2.5%. Two years ago, it was still 0.13%. Further increases are to be expected. The decision of the US Federal Reserve has an impact on financial markets around the world. Bonds will become significantly more attractive in the medium term – to the detriment of cryptocurrencies.
On the other hand, a different approach is pursued El Salvador. President Nakib Bukele is a supporter of Bitcoin and introduced cryptocurrency as an official and legal currency. Moreover, the Head of State speculates publicly with public funds. According to the President's statements on Twitter, El Salvador has invested around $225 million in Bitcoin. With the current price drop, the country has lost over US$60 million.
On offer: trade stocks and bitcoins for 1 euro at Trade Republic
Conclusion: Which cryptocurrency has a future?
There is no clear answer to the question "Which cryptocurrency has a future?" Blockchain technology has long taken root in all industries and will continue to play an important role in the financial world and in everyday life. Cryptocurrencies will probably still exist in 2030. However, this does not necessarily mean that prices will continue to rise. So far, cryptocurrencies have been used in the economy, which can ensure transactions are as fast and cheap as possible with a stable exchange rate – the most valuable cryptocurrency is not one of them. .
After some countries wanted to adopt a restrictive approach towards Bitcoin, investors are now mainly focused on the United States. How they shape their monetary policy will not only significantly influence the future of cryptocurrency in 2022, but also the future of cryptocurrency.
commercial republic
Finanzen.net Zero
BSDEX
Deposit/Withdrawal Options
Possible Crypto Savings Plans
Coin transfer via wallet possible
intuitive app interface
Bank transfer, credit card, Apple Pay
1 euro third-party fees per transaction, spread
intuitive app, transparent pricing model
high fees for small volume
relatively complex due to the social trading function
PayPal, credit card, bank transfer, other services
1% of order volume + spread, $5 withdrawal fee, $10/month inactivity fee after 12 months without login
large selection, social commerce
relatively high fees
none, but 1 euro surcharge for small amounts for transactions below 500 euro, plus spreads
intuitive app, crypto and stocks in one wallet
high fees for small amounts
clear interface
no custody fees, no order fees, minimum spread of 0.3%
no order commission or custody fees
50 euros minimum order
intuitive and clearly structured app
0.2% of trading volume for maker orders (min. EUR 0.01); 0.35% for taker orders (min. EUR 0.01)
no minimum order volume, clear application
only 7 tradable crypto assets
*Our independent experts regularly deal with suppliers of products and services. We will make the resulting articles available to you free of charge. COMPUTER BILD receives a small commission if you click on a link or enter into a contract with a linked provider. Note: The content of computerbild.de does not constitute specific investment recommendations and contains only general information. The authors, publishers and sources cited are not responsible for any losses resulting from the purchase or sale of the securities or financial products mentioned in the articles.
Experiences with RXK Capital (rxkcapital.org)? Fraud or does the payment arrive?
NFT Technologies action: sensational development ()

redpajama_set_name: RedPajamaCommonCrawl


Don't miss the last in the Thursday Night Concerts series tomorrow night at Discovery Green this fall .
This week, dance the night away to R &B divas Nikki Hill and Annika Chambers!
The concert begins at 7. Picnics, blankets, and lawn chairs are welcome, but glass containers and alcohol are prohibited. Don't worry—food, wine, and beer are available for purchase from The Lake House.


Absolutely, the New is Here!
A Rental Renovation – and newest team members!
Our 1 Year "Move to Mississippi" Anniversary!
This has been THE best day! We've worked hard, played hard, with so much gratitude because 1 year ago today is when Rip and I drove into Madison, Mississippi and today we can't THANK God enough!


The team at Dickson Plumbing and Gas Fitting would love to renovate your home and transform your kitchen, bathroom or laundry into the latest modern living space that would make you feel proud. We offer free consultation and advice before you begin the design process to avoid costly mistakes and enrich the overall experience.
Updating your bathroom can be very rewarding and it can be a straight forward operation if you plan ahead. Getting the right advice and guidance as to buying the correct plumbing fixtures and doing the set out will make the transformation go smoothly. The Plumbers at Dickson Plumbing regularly work on Bathroom Renovations, so we have the experience to make it easy for you.
Replacing a Kitchen is done in two stages, before the the old cupboards are removed we come in for a short time to disconnect the services,once the new kitchen has been install we do a Re-connect .
In steps for your bathroom or laundry we would begin by removing all your current PC items (toilet/tapware/basins, etc), reposition pipework, if needed and prepare the new canvas for the upgrade to commence. Once the tiling has been completed we then install all new upgraded PC items to complete the transformation of your new space.
Dickson Plumbing and Gas Fitting are fully qualified and experienced tradesmen with a vast range of skill sets. This combination of skills ensures that all Taps and fixtures are fitted off without any issues. Preparation is the key and it important to get it right also in a timely manner and with attention to detail.
Renovating you Laundry can be one of the smaller jobs to renovate.
We take pride in our work and find renovations a pleasurable experience. Like you, we love seeing your vision unfold before our eyes.
Finding an experienced plumber with home renovation knowledge can sometimes be difficult. Dickson Plumbing and Gas offer licensed and insured services for your project.


text:
Covid-19 Vaccine involvement
· Edited: Jan 4
At Alpine Air Express, we are proud to announce our involvement in the distribution of the Pfizer-BioNTech's COVID-19 vaccine, which is the first COVID-19 vaccine to receive FDA authorization for emergency use. On Monday, December 14, 2020, we launched our first shipments of the vaccine.
Our priority is to move the COVID-19 vaccines from our base hubs to our outstation areas as quickly and efficiently as possible. These outstation areas include the intermountain regions of Montana, North and South Dakota, Idaho, Wyoming, Colorado, Utah and Minnesota. Once the cargo arrives at its respective location, the vaccines will be distributed to hospitals, clinics and other medical facilities intended for front line healthcare workers and those working in close proximity to the virus.
"We will get through it together" has been a popular slogan in 2020; and, we here at Alpine are excited to do our part. As front line workers, we are teamed up with the brave doctors, nurses, other healthcare workers and first responders. Alpine's team is proud to work together helping each other during these unprecedented times.


text:
Władimir Pietrowicz Worobjow (ros. Владимир Петрович Воробьёв, ur. 28 czerwca 1876 w Odessie, zm. 31 października 1937 w Charkowie) – radziecki anatom. W 1924 wraz z biochemikiem prof. B.I. Zbarskim dokonał mumifikacji ciała Lenina, pomimo sprzeciwu jego rodziny. Inicjatywa ta wyszła od Biura Politycznego RKP(b).

Przypisy

Bibliografia
 Владимир Петрович Воробьев. [dostęp: 2010-11-17].
 Воробьёв Владимир Петрович (1876-1937). [dostęp: 2010-11-17].

Rosyjscy naukowcy
Urodzeni w 1876
Zmarli w 1937

redpajama_set_name: RedPajamaWikipedia


text:
Home Sports News Devils Porridge Museum holds first Commonwealth Sport's Day
Devils Porridge Museum holds first Commonwealth Sport's Day
On Saturday June 22nd, the Devil's Porridge Museum hosted its first Commonwealth Sport's Day. There were egg and spoon races, sack races, three legged races and even tug of war. A party atmosphere was created by local young musician Jay who played the piano while Chameleon Face painters provided the beautiful face art. To top it all, the Vice Chair of the Museum, Graham Croft, cooked delicious burgers and hot dogs on the barbecue.
"It was a really fun day" commented Judith Hewitt, Museum Manager. "We recently came across a photo of girls in World War One having a tug of war somewhere in Eastriggs. This got us thinking that we could do something similar. We also have a document in our archive which tells of a Sports Day in 1919 which was held on Eastriggs playing ground and included workers from HM Factory Gretna getting involved in different sports. We currently have an exhibition in the museum called 'The Contribution of the Commonwealth' which explores why Eastriggs is known as the Commonwealth village and it seemed like an ideal time to have a mini Commonwealth games of our own!"
The Devil's Porridge Museum has a full programme of summer activities for young people from July 1st to August 17th. All activities are free.
Devils Porridge first Commonwealth Sport's Day
Previous articleChance to Get South of Scotland on The Map With OVO Energy Tour of Britain National Land Art competition
Next article200 D&G health and social care staff sing about 'working together'


One major challenge USACE and other government agencies face is finding the right balance among the often conflicting concerns our society has related to our water resources. We want these resources to help fuel economic growth (navigation, hydropower). We want them to provide social benefits (recreation). And finally, we want to be sure that they are available for future generations (environmental protection and restoration).
USACE seeks to achieve this balance through an integrated approach to water resources management that focuses on regional solutions involving an array of stakeholders (i.e. other government agencies, environmental groups, businesses and private organizations). More information about this approach is available in the Civil Works Strategic Plan (pdf, 2.46 MB).
For more information read our Civil Works Program Brochure (pdf, 3.90 MB) or browse through the information for each business line.
2016 Fast Facts are available!


Mene modeled many of the popular dog breeds of his day and this alert dog is most likely a Miniature Pinscher or Min Pin as they are known today.
This sculpture of a Miniature Pinscher is one of his more obscure works. He has managed to capture the essence of this feisty breed in the bronze. The dog is in a playful pose and this is a charming piece which would be a classic addition to any collection of the artist's work. The dog is standing on an Oriental rug and the pattern can be seen on the base of the sculpture. This is a model that is seldom encountered. It was cast by Mene in his own foundry and is finished in a rich brown patina.
Animals in Bronze by Christopher Payne, pg. 234.
Les Animaliers by Jane Horswell, pg. 136.


Ahmedabad is a 5000-year-old city and has a lot to offer to a travelling chef.
Old 'Amdavad' is a culturally (culinary) rich and an architecturally opulent city. The old city bustles in a mish-mash of narrow streets and 'pols' enclosed in a 12-gate fort. Each Pol typically has a corner shop selling Chawana and farsan, and is a world in itself because these super narrow streets house some of the most beautiful and humongous havelis that you will ever see.
In one such Haveli I met Abhay Mangaldas, who has taken upon himself to restore these architectural gems and lure the people of "new" Ahmedabad into the old city by opening restaurants and cafes in the restored Havelis. One food conversation led to another, and a lunch seemed inevitable. At the lunch table Abhay made me taste a salad that stuck in my head forever.
A raw turmeric kachumber that was fresh, aromatic, mildly astringent and divinely combined with 'fafda', a pulled lentil crisp that gave it a perfect balance of flavours and a perfect contrast of textures. Surprisingly, raw turmeric and peanuts - even though underground veggies - are acceptably used in Jain food. There's a lot more to this rhizome that one finds out on the slightest scratch, literally.
Native to Tamil Nadu, this pre-Aryan spice has been known to mankind for over 4,000 years. The term 'Haridra' (haldi) is believed to have a Munda air to it; Munda is an ancient aborigine dialect. Even in the early Vedic times the only four spices recorded are mustard (Baja),a sour Citrus (Jambira), Turmeric (Haridra) and long pepper (Pippali).
Out of these, turmeric was regarded as the most auspicious because it was the most useful for the entire body. Turmeric stands tall as a cure in Ayurveda, Chinese, Unani and Siddha medicine thanks to Curcumin, a compound that is now also widely accepted in the West as a cancer buster.
Yet all I say is forget the past and experience the present, grate some fresh turmeric and raw papaya together, toss it up with chopped green coriander, lemon juice and black pepper and top a papad with this salad for a taste you will savour forever. A chef's promise!
1. In a pan put water, coconut milk, turmeric soaked sago and cook. Now, add mango pulp. 2. Take mango and coconut sago mixture in a glass. Keep it in the fridge to cool.
3. When cooled take it out from the fridge. Then in the same glass put mango pulp.
4. Cut a mango wedge and keep it on the glass.