\section{The Polychromatic Jaynes-Cummings Model in the Fock basis} \label{subsec:fockbasisproblem}

\begin{figure*}[ht]
	\subfloat[Degenerate basis \label{fig:basisdiagram}]{
			\begin{picture}(280,250)
				\includegraphics[width=280pt]{degeneratebasisdiagram}
			\end{picture}
	}
	\subfloat[Energy levels \label{fig:energylevels}]{
			\begin{picture}(210,250)
				\put(100,5){\includegraphics[width=110pt]{energylevels}}
				\put(57,12){$\ket{10-2\uparrow}$}
				%
				\put(4,41){$\ket{10-1\downarrow}$}
				\put(45,41){$\ket{1-1-1\uparrow}$}
				\put(56,53){$\ket{01-2\uparrow}$}
				%
				\put(15,76){$\ket{1-10\downarrow}$}
				\put(56,76){$\ket{1-20\uparrow}$}
				\put(15,88){$\ket{01-1\downarrow}$}
				\put(56,88){$\ket{00-1\uparrow}$}
				%
				\put(68,110){$\ket{000\downarrow}$}
				\put(27,110){$\ket{0-10\uparrow}$}
				\put(50,122){$\ket{-11-1\uparrow}$}
				%
				\put(23,146){$\ket{-100\uparrow}$}
				\put(60,146){$\ket{-110\downarrow}$}
				\put(15,157){$\ket{0-11\downarrow}$}
				\put(56,157){$\ket{0-21\uparrow}$}
				%
				\put(60,180){$\ket{-101\downarrow}$}
				\put(12,180){$\ket{-1-11\uparrow}$}
				\put(60,192){$\ket{-210\uparrow}$}
				%
				\put(60,220){$\ket{-201\uparrow}$}
			\end{picture}
	}
\caption{\label{fig:degeneratebasisdiagram} Feynman diagram (a) and dressed state energies (b) of a two-level atom interacting with a three frequency electromagnetic field. The Feynman diagram, which illustrates states with different spin and photon number (vertexes) and their interactions (wavy lines), demonstrates the complexity of the multifrequency Jaynes-Cummings model in the Fock-spin basis. The energy level diagram highlights the erroneous calculation of the dressed energies in this basis (dashed lines) in contrast with the correct result (solid lines). 
The field frequencies are $j \omega_f$ and $(j\pm1) \omega_f$, and the energy levels are shown when the spin is resonant with $\omega_0=j\omega_f$.
In (a) the open (filled) circles are spin up (down) states with $N+\Delta N$, each labelled by their energy difference from the state in the centre of the diagram. The orientation of the wavy lines indicates the mode of interaction, and the shaded bands contain states of constant photon number in mode 1, 2 or 3. States such as $A$ and $B$ are degenerate, which cause problems in perturbative or numerical calculations of the eigenstates of \eqnref{eq:polychromaticJCM}. This can be seen in (b) by the incorrect appearance of singlets and triplets, and the disagreement of the level shifts. The states in (b) are labelled with shorthand $\ket{\Delta n_1, \Delta n_2, \Delta_n , m}$ for $\ket{n_1+\Delta n_1, n_2 + \Delta n_2, n_3 + \Delta n_3, m}$ such that $\ket{000 m} = \ket{n_1,n_2,n_3 m}$.
}
\end{figure*}


We consider a two-level system, such as an atom, interacting with a multi-frequency field with fundamental frequency $\omega_f$ \footnote{The frequency $\omega_f$ need only be the greatest common denominator of the frequencies under consideration, but in principle could be the fundamental frequency $2 \pi c/L$ when one imposes periodic boundary conditions to quantise the electromagnetic field. Thus, requiring that the field frequencies are rationally related in this model is not an overly restrictive assumption.},
described by the polychromatic Jaynes-Cummings model with Hamiltonian
\begin{equation} \label{eq:polychromaticJCM}
	H = \tfrac12 \omega_0 \sigma_z + \sum_i \left[k_i \omega_f a_{k_i}^{\dagger} a_{k_i} + g_{k_i} \left( a_{k_i} \sigma_+ + a_{k_i}^{\dagger} \sigma_-\right) \right].
\end{equation}
The rotating-wave approximation has been made and this form is relevant for the systems considered in \secref{sec:calculations} and \secref{sec:effectiveHs}; the Rabi model without this approximation is treated in \secref{sec:quantumclassicalboundary}.
In this model $\pm \tfrac12 \omega_0$ are the spin eigenenergies, the field modes are indexed by $i$ and the $i^{\mathrm{th}}$ mode has frequency $k_i \omega_f$ where $k_i$ are integers and need not be commensurate. $g_{k_i}$ are the coupling spin-field constants. More generally the mode frequencies do not need to be commensurate.

A natural basis to describe the coupled atoms and photons comprises of Fock-spin states $\ket{n_{k_1},n_{k_2},... ,m}$, where $n_{k_i}$ is the number of photons in mode $k_i$ and $m \hbar$ the atom's spin projection. These are eigenstates of $H_0 = \tfrac12 \omega_0 \sigma_z + \sum_{\{k \}} k \omega_f a_{k}^{\dagger} a_{k}$ and an extension of the basis $\{ \ket{n,m} \vert n \in \mathds Z^+ , m=\pm\tfrac12 \}$ used extensively for the single-frequency Jaynes-Cummings model. However, calculations performed in this basis are often divergent. The problems associated with this basis for the multi-frequency model are caused by (\emph{i}) the nature of multi-frequency interaction $V= \sum_{\{k \}} g_k \left( a_k \sigma_+ + a_k^{\dagger} \sigma_- \right)$ and (\emph{ii}) degeneracy of these states. In section \ref{sec:basis} we introduce a new non-degenerate basis upon which the dynamics of the multi-frequency model \eqnref{eq:polychromaticJCM} can be solved using a wide range of conventional methods. In the remainder of this section we elaborate on the key differences between the single and multi-frequency Jaynes-Cummings model and the problem of working with Fock-spin states.

In striking contrast to the single-frequency model is that a vast number of different states can now interact with each other. For a single-frequency field only pairs of states $\ket{n,-\tfrac12}$ and $\ket{n-1,\tfrac12}$ are permitted to interact - under the rotating wave approximation - by the terms $a \sigma_+$ and $a^{\dagger} \sigma_-$. For multiple frequencies, atoms can climb the Jaynes-Cummings ladder by interacting on alternate modes.
The extended range of interactions, together with energy degeneracies of the Fock states $\ket{n_1,n_2,...}$ cause the standard numerical or perturbative approaches to breakdown due to divergences in their solutions.

It is well known for example, that an energy degeneracy between a pair of interacting states $\ket{\psi_1}$ and $\ket{\psi_2}$ causes the perturbative solutions to the Schr\"odinger equation to diverge due to the energy denominator $E_1-E_2=0$. Such degeneracies are common for multi-frequency fields of \eqnref{eq:polychromaticJCM}. Consider for example a field with three commensurate frequencies (i.e. $k_1$, $k_2=k_1+1$ and $k_3=k_1+2$), and the interaction between states $\ket{n_1,n_2,n_3,m}$ and $\ket{n_1+1,n_2-2,n_3+1,m}$. These states are degenerate since $k_1-2k_2+k_3 = 0$. They also interact via a fourth order interaction $\propto a_{k_1}^{\dagger} \sigma_- a_{k_2} \sigma_+ a_{k_3}^{\dagger} \sigma_- a_{k_2} \sigma_+$. The energy shift caused by this interaction incorrectly appears to diverge when calculated perturbatively. So called `non perturbative' analytic methods such as the resolvent formalism employed in \secref{sec:effectiveHs} also suffer from such degeneracies for the same reason.

Simple numerical approaches also suffer from these degeneracies.
To diagonalise \eqnref{eq:polychromaticJCM} numerically we must truncate the basis over a finite range of values $n_k$ for each mode. An efficient truncation method starts with state $\ket{n_1,n_2,n_3,-\tfrac12}$ for a three mode field for example, and includes any states which can be connected by applying $a_k \sigma_+$ and $a_k^{\dagger} \sigma_-$ up to some maximum number $N$ times. For example, including up to $N=1$ interactions, the basis consists of $\ket{n_1,n_2,n_3,-\tfrac12}$ together with $\ket{n_1-1, n_2, n_3,\tfrac12}, \ \ket{n_1, n_2-1
, n_3,\tfrac12}$ and $\ket{n_1, n_2, n_3-1,\tfrac12}$. For $N=2$ our basis includes six further states connected to $\ket{n_{k_1}, n_{k_2}, n_{k_3},-\tfrac12}$ via $a_{k_j}^{\dagger} \sigma_- a_{k_i} \sigma_+$ for all combinations $i\neq j$. The truncated basis is shown diagrammatically in \figref{fig:basisdiagram} for $N=3$ along with their energies. 
Even in this limited basis degeneracies are abundant; There are seven pairs of states which have the same energy and spin.

The problem encountered numerically is that the basis must be truncated and subsequently the degenerate states are not all surrounded by a similar set of interactions. Look for example at the state spin up state $\ket{n_1,n_2-1,n_3,\tfrac12}$ with energy $0$ which lies near the centre of the basis in \figref{fig:basisdiagram} and is marked $A$. Its degenerate partner $\ket{n_1-1,n_2+1,n_3-1,\tfrac12}$ , marked $B$ lies on the left side of this figure. The energy of $\ket{n_1,n_2-1,n_3,\tfrac12}$ is predominantly shifted by the interactions with the three neighbouring states, while $\ket{n_1-1,n_2+1,n_3-1,\tfrac12}$ only interacts directly with two neighbouring state in this truncated basis. Thus, the dressed energies of these states, which ought to be degenerate, are not. This can be seen by the dressed energies shown in \figref{fig:energylevels}. This problem cannot be solved simply by expanding the basis as further degenerate states will also near the edge of the extended basis - extending the basis eventually leads to bands of dressed energies with erroneous intra-band avoided crossings.

The solution to the problems caused by degenerate multi-frequency Fock states is to describe the field using sets of non-degenerate states $ \{ \ket{N} \}$ which we derive below. Each state is labelled by $N$ and has energy $N \omega_f$. Each set spans a non-degenerate subspace of Fock space which is sufficient to describe a multimode field in a product of coherent-states. We derive the multi-frequency number operator \eqnref{eq:diagfieldenergy}, and show that the action of the field creation and annihilation operators are $a_k \ket{N} = \alpha_k (\gamma_{N-k}/\gamma_N) \ket{N-k}$ and $a_k^{\dagger} \ket{N} = a_k^* (\gamma_N/\gamma_{N+k}) \ket{N+k}$ (\eqnref{eq:offdiagmatrixelements}).
These relations can be understood from the conservation of energy since $a_k$ annihilates a photon with energy $k \omega_k$.
We define $\gamma_N$ in \eqnref{eq:defgammaN} below, and show that in the limit of large photon numbers $\gamma_N/\gamma_{N\pm k} \approx 1$ using \eqnref{eq:gammaapprox}. 
With these relations one can study the quantum electrodynamics of multi-frequency atom photon interactions in coherent fields using the non-degenerate basis $\{ N \}$. 
The matrix elements of the Jaynes-Cummings Hamiltonian for example are readily evaluated in this basis using \eqnref{eq:approxoffdiagmatrixelements} and \eqnref{eq:diagonalmatrixelements}.
We explore these dynamics for strong fields in \secref{sec:effectiveHs}, and weak fields in \secref{sec:quantumclassicalboundary}.

\section{A non-degenerate dressed basis for a spin-half particle in a multi-frequency field} \label{sec:basis}

In this section we define the non-degenerate basis and derive the action of the field operators in this basis.
We assume each frequency mode of the field is in a coherent-state
\begin{equation}
	\ket{\alpha_k} = e^{- \abs{\alpha_k}^2 / 2} \sum_{n_k=0}^{\infty} \frac{\alpha_k^{n_k}}{\sqrt{n_k!}} \ket{n_k}
\end{equation}
Thus, the field is in the tensor product of coherent-states defined by the set of complex numbers $\alpha_k$, which we write as
\begin{equation}
	\ket{ \{ \alpha_k \}} = {\prod_k}^{\otimes} \ket{\alpha_k}.
\end{equation}
We next express such field states in a non-degenerate basis $\{\ket{N} \}$.

The Hilbert space of the field $\mathscr H_F$ is spanned by the set of Fock states $\{ \ket{n_k} \}$. We partition this space into distinct subspaces $\mathscr E_N$ of energy $N \omega_f$ (see \appref{ap:partition} for proof). Formally, $\mathscr E_N = \left \{ \ket{\psi} \vert H_F\ket{\psi} = N \hbar \omega_f \ket{\psi} \right \}$, where $N$ can be any positive integer.
The Hamiltonian $H_F = \sum_k k \hbar \omega_f a_k^{\dagger} a_k$ is the energy of the field minus the vacuum energy.

Let $P_N:\mathscr H_F \rightarrow \mathscr E_N$ be the projector onto the $N^{\text{th}}$ subspace. Given the state $\ket{\{\alpha_k \}}$ we define the state $\ket{N}$ as
\begin{equation}
	\ket{N} = \frac{P_N \ket{\{ \alpha_k \}}}{\sqrt{\bra{\{\alpha_k \}} P_N \ket{\{ \alpha_k\} } }},
\end{equation}
which has norm $\braket{ N } = 1$ since $P_N^2=P_N$. 

Since $\{ \mathscr E_N \}$ is a partition of $\mathscr H_F$, we can write the identity operator $\mathds 1 = \sum_{N=0}^{\infty} P_N$. Applying this to coherent-state $\ket{\{ \alpha_k \}}$ expands the field's state in terms of states $\ket{N}$,
\begin{equation} \label{eq:Nexpansion}
	\ket{ \{ \alpha_k\} } = \sum_{N=0}^{\infty} \gamma_N \ket{N}
\end{equation}
where 
\begin{equation} \label{eq:defgammaN}
\gamma_N = \sqrt{\bra{\{\alpha_k \}} P_N \ket{\{ \alpha_k\}}} \in \mathds R.
\end{equation} 
Importantly, this expansion of $\ket{\{ \alpha_k\}}$ is on a non-degenerate basis, unlike its expansion in terms of Fock state products. We note that the basis $\{ \ket{N} \}$ is specific to the set $\{ \alpha_k \}$ of complex numbers which describe the coherent-state of the field. Furthermore, the basis does not span $\mathscr H_F$ but is sufficient to describe this specific state of the field. This subspace $\{\ket{N} \} \subset \mathscr H_F$ is closed under the action of $a_j$, a property which is inherited from the coherent-state of the field. Although it is not strictly closed under $a_j^{\dagger}$, it approximately closes to the degree of $\braket{\alpha,1}{\alpha}$, the over-lap between a coherent-state and a photon added coherent-state \cite{agarwal1991}, which falls off rapidly for $\alpha$ larger than unity.

We now turn to the matrix elements of \eqnref{eq:polychromaticJCM} (and \eqnref{eq:shrodrabi}) in this basis. The field energy in the diagonal terms is 
\begin{equation} \label{eq:diagfieldenergy}
	\bra{N} \sum_k k \hbar \omega_f a_k^{\dagger} a_k \ket{N}=N \hbar \omega_f,
\end{equation}
since $\ket{N} \in \mathscr E_N$. For the off-diagonal terms we need to evaluate $\bra{N} a_j \ket{M}$ and its Hermitian conjugate $\bra{M} a_j^{\dagger} \ket{N}$. To do this we first prove that $a_j \ket{M} \in \mathscr E_{M-j}$.
Applying $H_F$ to $a_j \ket{M}$m,
	\begin{IEEEeqnarray*}{rCl}
		H_F a_j \ket{N} 
			&=& \hbar \omega_f \sum_{l=1}^{\infty} l a_l^{\dagger} a_l a_j \ket{N} \\
			&=& \hbar \omega_f \left[j a_j^{\dagger} a_j a_j+ a_j \sum_{l \neq j} l a_l^{\dagger} a_l \right] \ket{N} \\
			&=& \hbar \omega_f \left[j \left (a_j a_j^{\dagger} - 1\right ) a_j + a_j \sum_{l \neq j} l a_l^{\dagger} a_l \right] \ket{N} \\
			&=& \hbar \omega_f a_j \left[ -j + \sum_{l=1}^{\infty} l a_l^{\dagger} a_l \right] \ket{N} \\
			&=& (N-j) \hbar \omega_f a_j \ket{N}.
	\end{IEEEeqnarray*}
This shows that the state $a_j \ket{N}$ is an eigenstate of $H_F$ with energy $(N-j)\hbar \omega_f$. Thus $a_j \ket{N} \in \mathscr E_{N-j}$.

From this we note that $\bra{N} a_j \ket{M} = 0 \forall M \neq N+j$. We next evaluate $\bra{N} a_j \ket{\{ \alpha_k \}}$ in two ways. On the one hand,
\begin{IEEEeqnarray*}{rCl}
	\bra{N} a_j \ket{\{ \alpha_k \}} 
		&=& \bra{N} \alpha_j \ket{\{ \alpha_k \}} \\
		&=& \alpha_j \bra{N} \sum_M \gamma_M \ket{M} \\
		&=& \alpha_j \gamma_N. \IEEEyesnumber \label{eq:Najalphak1}
\end{IEEEeqnarray*}
On the other hand,
\begin{IEEEeqnarray*}{rCl}
	\bra{N} a_j \ket{\{ \alpha_k \} }
		&=& \bra{N} a_j \sum_M \gamma_M \ket{N} \\
		&=& \gamma_{N+j} \bra{N} a_j \ket{N+j} \IEEEyesnumber \label{eq:Najalphak2}
\end{IEEEeqnarray*}
since only $a_j \ket{M} = a_j \ket{N+j}$ is in the same subspace $\mathscr E_N$ as $\ket{N}$. Equating the two expression \eqnref{eq:Najalphak1} and \eqnref{eq:Najalphak2} then rearranging, we find
\begin{IEEEeqnarray*}{rCl} \IEEEyesnumber \label{eq:offdiagmatrixelements}
	\bra{N} a_j \ket{M} &=& \frac{\gamma_N}{\gamma_{N+j}} \alpha_j \delta_{M,N+j}. \IEEEyessubnumber
\end{IEEEeqnarray*}
Taking the Hermitian conjugate and switching labels $N$ and $M$,
\begin{IEEEeqnarray*}{rCl}
	\bra{N} a_j^{\dagger} \ket{M} &=& \frac{\gamma_{N-j}}{\gamma_N} \alpha_j^* \delta_{M,N-j}.\IEEEyessubnumber
\end{IEEEeqnarray*}
To evaluate these expression we must first evaluate $\gamma_N^2$.

The probability distribution of states $\ket{N}$ is given by $\bra{\{\alpha_k \}} P_N \ket{ \{ \alpha_k \} } = \gamma_N^2$. This distribution has mean value $\bar N = \sum_k k \abs{\alpha_k}^2$ since $\langle \hat N\rangle = \sum_k \langle  k a_k^{\dagger} a_k \rangle = \sum_k k \abs{\alpha}_k^2$, where $\langle \hat O \rangle$ is the expectation value of an operator $\hat O$. The variance of $\gamma_N^2$ is $\sigma_N^2 = \sum_k k^2 \abs{\alpha_k}^2$ using $\langle N^2 \rangle - \bar N^2$ and the commutator relation $[a_k,a_k^{\dagger}]=1$.
We show in \appref{ap:gammaapprox} that the coefficients $\gamma_N^2$ tend towards the Gaussian
\begin{equation} \label{eq:gammaapprox}
	\gamma_N^2 \approx \frac{k_0}{\sqrt{2 \pi} \sigma_N} \exp \left[ - \frac{(N-\bar N)^2}{2 \sigma_N^2} \right],
\end{equation}
in the limit that the excitation of each mode $\abs{\alpha_k}^2$ is large, where $k_0$ is the greatest common denominator of $\{k\}$.

From the distribution \eqnref{eq:gammaapprox} we see that the states with non-negligible population are those distributed within $\sim \sigma_N$ of $\bar N$. The value of $\sigma_N$ is bounded such that $\sigma_N \leq \sqrt {\bar N} \sum_k k$ since $\lambda_k \leq 1/k$. Thus the standard deviation $\sigma_N \ll \bar N$ provided that $\bar N \gg 1$ and $\sum_k k \ll \sqrt{\bar N}$. In this case, we see from \eqnref{eq:gammaapprox} that the ratio $\gamma_N/\gamma_{N+j}\approx 1$. Under this approximation the off diagonal matrix elements ($N \neq M$) of \eqnref{eq:polychromaticJCM} are
\begin{IEEEeqnarray*}{rCl} \IEEEyesnumber \label{eq:approxoffdiagmatrixelements}
	\bra{N,1/2} H \ket{M,-1/2} &=& \frac{g}{\sqrt 2}\alpha_j e^{i \mathbf k_j \cdot \mathbf r}\delta_{M,N+j},\IEEEyessubnumber \\
	\bra{N,-1/2} H \ket{M,1/2} &=& \frac{g}{\sqrt 2} \alpha_j^* e^{i \mathbf k_j \cdot \mathbf r}  \delta_{M,N-j}.\IEEEyessubnumber
\end{IEEEeqnarray*}
The diagonal terms are given exactly by
\begin{equation} \label{eq:diagonalmatrixelements}
	\bra{N,\pm 1/2} H \ket{N,\pm 1/2} = \pm \frac{1}{2}\hbar \omega_0 + N \hbar \omega_f.
\end{equation}

\section{Eigenenergies calculated in the non-degenerate basis} \label{sec:calculations}

\begin{figure}[!t]
	\subfloat[Nondegenerate basis \label{fig:ndbasis}]{
	\begin{picture}(150,300)
		\put(40,10){\includegraphics[width=50pt]{nondegeneratebasisdiagram}}
		\put(0,13){$\ket{N-3,\downarrow}$}
		\put(95,13){$\ket{N-3-k_2,\uparrow}$}
		\put(0,58){$\ket{N-2,\downarrow}$}
		\put(95,58){$\ket{N-2-k_2,\uparrow}$}
		\put(0,103){$\ket{N-1,\downarrow}$}
		\put(95,103){$\ket{N-1-k_2,\uparrow}$}
		\put(0,148){$\ket{N,\downarrow}$}
		\put(95,148){$\ket{N-k_2,\uparrow}$}
		\put(0,193){$\ket{N+1,\downarrow}$}
		\put(95,193){$\ket{N+1-k_2,\uparrow}$}
		\put(0,237){$\ket{N+2,\downarrow}$}
		\put(95,237){$\ket{N+2-k_2,\uparrow}$}
		\put(0,283){$\ket{N+3,\downarrow}$}
		\put(95,283){$\ket{N+3-k_2,\uparrow}$}
	\end{picture}
	}
	\subfloat[Energy levels \label{fig:ndlevels}]{
	\begin{picture}(69,295)
		\put(0,0){\includegraphics[width=69pt]{ndbelevels}}
	\end{picture}
	}
\caption{\label{fig:ndbasisandlevels}
Feynman diagram (a) of a two-level atom in a three frequency field represented in the non-degenerate basis, and (b) the ladder of eigenenergies calculated using this basis. The Feynman diagram here demonstrates the simplicity of working in this basis, compared to \figref{fig:basisdiagram} which represents the same problem in the Fock-spin basis. Furthermore, the eigenenergies calculated in a truncated basis correctly form a set of doublet states, unlike those in \figref{fig:energylevels}. In (a) atom-photon states are represented by circles at each vertex as labelled, and the diagram is truncated after states that interact with $\ket{N, \downarrow}$ or $\ket{N-j,\uparrow}$ via a minimum of three photons.
The wavy lines show the interactions between them through the field with frequencies $\omega_1, \ \omega_2$ and $\omega_3$ equal to $(j-1) \omega_f, \ j \omega_f$ and $(j+1)\omega_f$. The state energy increases going up the page with resonant states shown at the same height. In this example $\omega_2$ is chosen to be resonant with the spin-half.
Thus, wavy horizontal lines represent interactions $\propto a_2 \sigma_+$ or $a_2^{\dagger} \sigma_-$. Diagonal wavy lines in one direction represent $a_1 \sigma_+$ and $a_1^{\dagger} \sigma_-$, and the other $a_3 \sigma_+$ and $a_3^{\dagger} \sigma_-$. At this resonance the energies shown in (b) of pairs of states $\ket{N, \downarrow}$ and $\ket{N-{k_2}, \uparrow}$ etc. are equal and form singlets. These are split into doublets in the dressed state basis as is typical in the Jaynes-Cummings model. Solid lines show state energies calculated in this truncated non-degenerate basis larger basis and are in very close agreement with more precise calculations in a large truncated basis which includes states up to 5 interactions from $\ket{N, \downarrow}$.
}
\end{figure}

Immediately one finds a multi-frequency model's Hamiltonian matrix can be diagonalised without the problems caused by level degeneracies.
The set of basis states can be efficiently truncated in the same way that produced the basis in \figref{fig:basisdiagram}. 
We consider again the three mode field with frequencies $(j-1)\omega_f, \ j \omega_f$ and $(j+1)\omega_f$ for some integer $j$. These modes are the same as in \secref{subsec:fockbasisproblem} when $j=1$.
Starting from a state $\ket{N,-\tfrac12}$ we include all states connected to this by up to $N$ interactions.

Figure \ref{fig:ndbasis} shows the diagram for $N=3$, which is significantly simpler than the degenerate diagram in \figref{fig:basisdiagram}. In fact the two diagrams become topologically equivalent if one identifies the degenerate states with one another.
Figure \ref{fig:ndlevels} shows some of the eigenenergies calculated in the non-degenerate basis truncated at $N=5$. In this case the energies do not suffer the problems seen in the degenerate basis. On the left hand side of this diagram lie the energy levels of pairs of states $\ket{n,-\tfrac12}$ and $\ket{n-j,\tfrac12}$, $\ket{n+1,-\tfrac12}$ and $\ket{n-j+1,\tfrac12}$, ... Each pair are resonant since we have taken $\omega_0-j \omega_f=0$. These dressed state energies (far right of \figref{fig:ndlevels}) are shifted apart by the interaction between them forming the doublets. The level shifts are calculated most accurately for states near the centre of the basis. When $N=3$ states near the edge of the basis (the highest or lowest energies in this example) can be seen to deviate from more accurate calculations with $N=5$.
\begin{figure}
	\includegraphics[height = 5.8cm]{evals3fplot.eps}
\caption[]{\label{fig:evals3f} 
Energy levels of a spin-half in a three frequency field versus detuning between the spin-half and the mean field frequency. Each field is resonant at either $\Delta_0 = -1, \ \Delta_0 = 0$ or $\Delta_0 = 1$.
The dashed lines show energies $N\omega_f \pm \tfrac12 \Delta_0$ of the states $\ket{N, \pm \tfrac12}$. The interactions between these states causes avoided crossings in the dressed state energies (solid lines) at resonance $\Delta_0=0$ and $\pm1$. The Rabi frequencies are equal between the fields; $\Omega = 0.2 \omega_f, \ 0.5 \omega_f$ and $0.8 \omega_f$ for the black, grey and light grey lines.
The avoided crossings are shifted towards $\Delta_0$ due multi-frequency processes involving the off resonant fields. At $\Delta_0=\pm2$ further avoided crosses emerge due to three photon interactions, clearly visible for stronger Rabi frequencies. The light grey line shows how the the central three avoided crossings merge to form a single resonance as the Rabi frequencies approaches the field's frequency spacing $\omega_f$.
}
\end{figure}

Figure \ref{fig:evals3f} shows a generalisation of this energy level diagram, plotted as a function of $\Delta_0 = \omega_0-j \omega_f$ along the horizontal axis. The nondegenerate basis states (dashed lines) $\ket{N,\tfrac12}$ increase linearly with $\Delta_0$ while $\ket{N,-\tfrac12}$ decrease linearly. These basis states cross whenever $\Delta_0/\omega_f$ takes an integer value. We label each crossing by this integer value, $q$. For $q=0$ the states $\ket{N,-\tfrac12}$ and $\ket{N+j,\tfrac12}$ are resonant as they are in \figref{fig:ndbasisandlevels}. Generally the $q^{\mathrm{th}}$ resonance is between $\ket{N-\tfrac12}$ and $\ket{N-(j+q),\tfrac12}$.


The dressed levels (solid lines) are quite clearly forced apart for $q=-1,0,1$ predominantly due to the first order interactions with modes $j-1, \ j$ and $j+1$ respectively.
The $q=\pm1$ resonances are also shifted inwards slightly due to the level shift caused by second order interactions involving the off resonant modes. These shifts are the multi-frequency equivalent to the Bloch-Siegert shift. The contributing shifts to the $q=0$ resonance cancel when the amplitudes of modes $j\pm1$ are balanced.
Less prominent avoided crossings occur for $q=\pm2$. These become noticeable as the field amplitudes increase (grey solid line) and are caused by resonant three photon interactions. Using our non-degenerate basis we are able to investigate these multi-frequency resonances and level shifts analytically in the section below.

\section{Multi-photon processes of a spin-half particle in strong multi-chromatic field} \label{sec:effectiveHs}

We now derive effective Hamiltonians which describe the system near each of the crossings of the levels of $H_0$. For this we use the resolvent formalism, allowing us to account for the level shifts and interaction processes which couple these levels to an arbitrary precision.

Briefly, we introduce the formalism of the resolvent $G(z)=(z-H)^{-1}$ which is the advanced Green's function for the Shro\"odinger equation. Its Fourier transform is the time evolution operator and we can interpret $G(z)$ as the propagator in the complex frequency ($z$) space.
We are able to describe the dynamics between two energy states of $H_0$ (close to where they cross) using the projection of the resolvent $PG(z)P=(z-H_{\mathrm{eff}})^{-1}$ onto subspace $\mathscr E$ these states span. $H_{\mathrm {eff}}$ is the effective Hamiltonian between the two states and is given by
\begin{IEEEeqnarray}{rCl} \label{eq:Heffgeneral}
	H_{\mathrm{eff}} &=& P H P + PR(z)P 
\end{IEEEeqnarray}
where $R(z)$ has the series expansion \cite{Cohen92}
\begin{equation} \label{eq:Rgeneral}
	R(z) = VQG_0QV + VQG_0QVQG_0QV +  \cdots .
\end{equation}
$P$ is the projection from $\mathscr H$ to $\mathscr E$ and $Q$ the projector onto its compliment $\bar{\mathscr E}$ and $G_0=(z-H_0)^{-1}$ is the free propagator. The series expansion of $R(z)$ converges provided the energies within $\mathscr E$ are closer to each other than the energy of any state in $\bar{\mathscr E}$ with which they interact. Thus, the non-degenerate basis introduced in \secref{sec:basis} is also necessary here to calculate multi-frequency effects using the resolvent formalism.
The multi-frequency effects are encompassed by $R(z)$ where each term describes an interaction in $\mathscr E$ mediated via virtually excited states in $\bar {\mathscr E}$. These terms are easily evaluated analytically when the two states of interest are close to resonance with each other, but far from resonance with any other state, in which case $z$ is well approximated by their mean energy.

We interpret the effective Hamiltonian in powers of $V$.
The terms of $H_{\mathrm{eff}}$ zeroth order in $V$ are merely the energies of the unperturbed levels of $H_0$. The first order terms in $V$ describe the direct interactions between states of $\mathscr E$. The second order terms, for the interaction discussed in this work, describe interactions which return to the same state they started, via an intermediate state in $\bar{\mathscr E}$. These self interactions shift the levels energy levels of $PH_0P$. Third order terms are between different initial and final states in $\mathscr E$ via two intermediate states in $\bar{\mathscr E}$.


To define the effective Hamiltonian at a given resonance we must find the relevant multi-photon processes which modify $PHP$. These processes, described generally by \eqnref{eq:Rgeneral}, are those that begin and end with states in $\mathscr E$ but have all intermediate states in $\bar{\mathscr E}$. Transitions between initial, intermediate and final states can only be those driven by the field modes present. 
Provided the Rabi frequencies are smaller than $\omega_f$, the leading terms are those lowest order in $V$. For three frequencies the processes up to third order are shown diagramatically in \figref{fig:q0diagrams} for the $q=0$ resonance.
The term for each multi-photon process is evaluated from matrix elements of individual transitions involved using \eqnref{eq:approxoffdiagmatrixelements} and the free propagators ($G_0$) of the intermediate states. The later are known from the intermediate state energies evaluated using \eqnref{eq:diagonalmatrixelements}.

We calculate effective Hamiltonians up to third order in $V$, first for the three frequency example considered of \figref{fig:evals3f}, then more generally for a large number of evenly spaced off-resonant frequencies. For convenience we abbreviate the resonant states $\ket{N-(j+q),\tfrac12}$ and $\ket{N, -\tfrac12}$ to $\ket{a}$ and $\ket{b}$ respectively. All energies are expressed relative to the mean energy $E_0=(E_a+E_b)/2=\left [N-(j+q)/2 \right ] \omega_f$ at the resonance under consideration.

\begin{figure}[!t]
	\begin{picture}(250,215)
	
		\put(75,202){$V_{ab}=$}
		\put(105,200){\includegraphics[width=40pt]{vab}}
		\put(104,190){$\ket{b}$}
		\put(137,190){$\ket{a}$}
	
		\put(0,128){$R^{(2)}_{aa}\quad=$}
		\put(56,138){$\ket{a}$}
		\put(25,93){\includegraphics[width=40pt]{raa1}}
		\put(80,128){$+$}
		\put(75,125){\includegraphics[width=40pt]{raa2}}
		\put(106,115){$\ket{a}$}
	
		\put(130,128){$R^{(2)}_{bb}=$}
		\put(160,125){\includegraphics[width=40pt]{rbb1}}
		\put(158,115){$\ket{b}$}
		\put(185,128){$+$}
		\put(203,138){$\ket{b}$}
		\put(205,93){\includegraphics[width=40pt]{rbb2}}
	
		\put(30,43){$R_{ab}^{(3)}\quad=$}
		\put(75,40){\includegraphics[width=40pt]{rab1}}
		\put(73,30){$\ket{b}$}
		\put(106,30){$\ket{a}$}
		\put(124,42){$+$}
		\put(138,52){$\ket{b}$}
		\put(171,52){$\ket{a}$}
		\put(140,7){\includegraphics[width=40pt]{rab2}}
	\end{picture}
\caption{\label{fig:q0diagrams} 
Diagrams for one, two and three photon processes for a spin-half in a three frequency field. These diagrams represent the contributions to the effective Hamiltonian describing the resonance between the spin and the middle frequency field mode. The effective Hamiltonian, expanded up to third order interactions, is written $H_{\mathrm{eff}} = PH_0P + PVP + PVQG_0VP + PVQG_0VQG_0VP$. Figure (a) represents $P VP$ between $ \ket{a}=\ket{N, \downarrow}$ and $\ket{b}=\ket{N-j,\uparrow}$. Figure (b) represents the different two photon processes of $R^{2} = PVQG_0VP$ starting and returning to either $\ket{a}$ and $\ket{b}$. Figure (c) represents the three photon processes $R^{3} = PVQG_0VQG_0VP$ between $\ket{a}$ and $\ket{b}$. Note that we have not shown $R^{3}_{ba}$ which is the Hermitian conjugate of $R^{3}_{ab}$ and has the same diagrammatic representation.
}
\end{figure}

We begin by looking at the crossing $q=0$, where the initial and final states are $\ket a =\ket{N-j,\tfrac12}$ and $\ket b=\ket{N,-\tfrac12}$.
We evaluate \eqnref{eq:Heffgeneral} for the processes shown in \figref{fig:q0diagrams}.
Up to first order
\begin{equation} \label{eq:q0firstorderHeff}
	H_{\mathrm{eff}} = P H P = \frac12 \left(
		\begin{array}{cc}
			\Delta_0 & \Omega_j \\
			\Omega_j^* & -\Delta_0
		\end{array}
	\right).
\end{equation}
where $\Omega_j = g_j (\gamma_{N-j}/\gamma_N) \alpha_j\approx g_j \alpha_j$. We have labelled $\Delta$ with the subscript $q=0$ to show which resonance this detuning is from.
At this level of approximation the system behaves like it is coupled only by a single-frequency $j \omega_f$.

The second order term of $R(z)$, $PVQG_0(z)QVP$, shifts the energies of $\ket{a}$ and $\ket{b}$ due virtual excitation of the intermediate states $\ket{N\pm1,-\tfrac12}$ and $\ket{N-j\pm, \tfrac12}$ respectively. 
Thus, at second order the diagonal elements of \eqnref{eq:q0firstorderHeff} are shifted by
\begin{IEEEeqnarray*}{rCl} \label{eq:qzerolevelshifts}
	R^{(2)}_{aa}(z) &=&
		\frac{\abs{\Omega_{j-1}}^2/4}{z+\omega_f+\tfrac12 \Delta_j}
		+\frac{\abs{\Omega_{j+1}}^2/4}{z-\omega_f+\tfrac12 \Delta_j} \IEEEyesnumber \IEEEyessubnumber\\
	R^{(2)}_{bb}(z) &=& 
		\frac{\abs{\Omega_{j-1}}^2/4}{z-\omega_f-\tfrac12 \Delta_j}
		+\frac{\abs{\Omega_{j+1}}^2/4}{z+\omega_f-\tfrac12 \Delta_j}, \IEEEyessubnumber
\end{IEEEeqnarray*}
where $z$ can be approximated by zero.
The resonance is therefore shifted to $\Delta_0=R_{bb}^{(2)}(0)-R_{aa}^{(2)}(0)$.
However, when the amplitudes of the modes $j \pm 1$ are equal the two terms which contribute to each shift cancel to zero when $\Delta_j=0$, so in this case the resonance is unchanged.

The third order contribution modifies the interaction between $\ket{a}$ and $\ket{b}$, coupling them additionally via two intermediate states $\ket{N-j+1,\tfrac12}$ and $\ket{N+1,-\tfrac12}$ or $\ket{N-j-1,\tfrac12}$ and $\ket{N-1,-\tfrac12}$. 
The matrix elements for these terms are
\begin{IEEEeqnarray*}{rCr}
	R^{(3)}_{ab}(z) &=& \frac{\Omega_{j-1}\Omega_j^* \Omega_{j+1}}{8}
		\left[
			\frac{1}{(z-\omega_f-\tfrac12\Delta_j)(z-\omega_f+\tfrac12 \Delta_j)}
		\right.\\
	&&
		\left. +\frac{1}{(z+\omega_f-\tfrac12\Delta_j)(z+\omega_f+\tfrac12 \Delta_j)}
		\right] \IEEEeqnarraynumspace \IEEEyesnumber
\end{IEEEeqnarray*}
and its Hermitian conjugate. At resonance, this third order interaction increases the splitting between $E_a$ and $E_b$ further than the value $\Omega_j$ expected for a single-frequency field.

The same approach can be adopted for the resonances $q=\pm1$ with similar results. To lowest order the matrix elements are $\pm \tfrac12 \Delta_q$, where $\Delta_q=\omega_0-(j+q)\omega_f$. The first order contribution to the interaction matrix elements are $\tfrac12 \Omega_{j+q}$ and its Hermitian conjugate.
The result \eqnref{eq:qzerolevelshifts} for the leading order level shift at $q=0$ can be generalised to
\begin{IEEEeqnarray*}{rCl}
	R^{(2)}_{aa} &=& 
		\sum_{p\neq q, \abs{p}\leq 1} \frac{\abs{\Omega_{j+p}}^2}{4} \frac{1}{z-\tfrac12 \Delta_{j+q}-(q-p)\omega_f},\\
	R^{(2)}_{bb} &=&
		\sum_{p\neq q, \abs{p}\leq 1} \frac{\abs{\Omega_{j+p}}^2}{4} \frac{1}{z+\tfrac12 \Delta_{j+q}+(q-p)\omega_f}.
\end{IEEEeqnarray*}
which not only applies to the cases $q=\pm1$ but for any $q$. Unlike the $q=0$ case, these level shifts do not cancel on resonance. For $q=1$ these terms shift the resonance to negative value of $\Delta_{j+1}$, while for $q=-1$ the shift is to a positive value of $\Delta_{j-1}$.

\begin{figure*}
	\subfloat[Energy Levels \label{fig:q2resonance}]{
		\includegraphics[height=157pt]{q2avoidedcrossingplot}
		}
	\subfloat[Excitation Spectrum \label{fig:q2spectrum}]{
		\includegraphics[height=157pt]{q2spectrumplot}
		}
	\caption{{\label{fig:q2resonance}} 
Dressed state energies (a) and excitation probabilities (b) for a three photon resonance of a spin-half in a three frequency field.
This figure demonstrates the accuracy of the effective Hamiltonian discussed in \secref{sec:effectiveHs}, where dashed lines plot the analytic results of the effective theory, and solid lines the exact numerical results. 
In (a) the three photon avoided crossing becomes larger and shifts left as the Rabi frequencies increase. Here the field frequencies are $j \omega_f$ and $(j\pm1) \omega_f$. The results are plotted against the detuning $\omega_0-j\omega_f$ such that the unshifted resonance occurs at $\Delta = 2$ (i.e $q=2$). All fields have the same Rabi frequencies or equivalently their coherent amplitudes $\alpha_1, \alpha_2$ and $\alpha_3$ are the same. Curves are labelled by these Rabi frequencies in units of $\omega_f$. In (b) this manifests as a broadening and shift left of the peak excitation probability. The excitation probability was calculated for a square pulse of all three frequencies with duration $\tilde \Omega_{\mathrm{eff}} t = \pi$ where $\tilde \Omega_{\mathrm{eff}}=2 R_{ab}^{(3)}$ is the effective three-photon Rabi frequency given in \eqnref{eq:q2omegaeff}.
}
\end{figure*}

For level crossings where $\abs{q}>1$ there is no direct coupling between the states $\ket{a}$ and $\ket{b}$ within the resonant subspace $\mathscr E$. For $q=\pm2$ or $q=\pm3$ the resonant states are only coupled by third and higher order processes. In general we find that at the $q^{\mathrm{th}}$ resonance, between states $\ket{a}=\ket{N-(j+q),\tfrac12}$ and $\ket{b}=\ket{N,-\tfrac12}$, the lowest order interaction is of order $q$ in $V$ when $q$ is odd, and of order $q+1$ when $q$ is even. For odd $q$ the effective interaction matrix element between the resonant states is
\begin{IEEEeqnarray*}{rCl}
	R_{ab}^{(q)} &= 
		\frac{\left({\Omega_{j+1}^{\frac12(q+1)}}\right)^* \Omega_{j-1}^{\frac12 (q-1)}}{2^q} 
		&\prod_{n=1}^{\tfrac12(q-1)} \left(\frac{1}{z+\tfrac12 \Delta_{j+q}+2 n \omega_f} \right.\\
	&& \left.
		\times \frac{1}{z-\tfrac12 \Delta_{j+1} - 2 n \omega_f} \right). \IEEEyesnumber
\end{IEEEeqnarray*}
On resonance ($\Delta_q=0$), and approximating $z$ with zero, this effective interaction matrix element simplifies to
\begin{equation}
	R_{ab}^{(q)} \approx 
		\left(\frac{\Omega_{j+1}^*}{2}\right)^{\frac12(q+1)}  \left(\frac{-\Omega_{j-1}}{2} \right)^{\frac12 (q-1)}
		\frac{1}{(2 \omega_f)^{q-1} \left(\frac{q-1}{2}!\right)^2}
\end{equation}
This shows the multi-photon resonant interaction strength scales generally as
\begin{equation}
	R_{ab}^{(q)} \sim 8 \Omega \left(\frac{\Omega}{2 \omega_f}\right)^q \left(\frac{q}{q!}\right)^2
\end{equation}
with a characteristic Rabi frequency $\Omega\sim \Omega_j \forall j$.

For even $q$ there are $q+1$ different processes which contribute to the interaction matrix element, each of order $q+1$ in $V$. While the interaction processes are more numerous for even $q$ resonances, the coupling is typically weaker than that for $q-1$ unless the number of processes exceeds the relevant ratio $\Omega / \Delta$. Due to the multiple processes involved, the lowest order multi-photon matrix element for even $q$ is given by the more complicated expression
\begin{widetext}
\begin{IEEEeqnarray*}{rCl}
	R_{ab}^{(q+1)} =
		&\sum_{s=1}^{q/2} &\left \{
			 \frac{\Omega_j^*}{2} \frac{\Omega_{j+1}/2}{z-(q-2s+1)\omega_f-\tfrac12 \Delta_{j+q}} 
			\left[ \prod_{n=1}^{s-1} \frac{\Omega_{j-1}^*/2}{z+2 n \omega_f+\tfrac12 \Delta_{j+q}}
				\frac{\Omega_{j+q}/2}{z-q \omega_f + (2n+1)\omega_f-\tfrac12 \Delta_{j+q}} \right] \right. \\
	&& \left. \times
			\frac{\Omega_{j+1}/2}{z+(2s-1)\omega_f+\tfrac12 \Delta_{j+q}} 
			\left[
				\prod_{n=1}^{q/2-s-1} \frac{\Omega_{j-1}^*/2}{z+q \omega_f - (2n-1)\omega_f+\tfrac12 \Delta_{j+q}}
				\frac{\Omega_{j+1}/2}{z-2 n\omega_f-\tfrac12 \Delta_{j+q}}
			\right] \right\} \\
	&+ \sum_{s=0}^{q/2} & \left\{
		\frac{\Omega_j}{2}
		\left[
			\prod_{n=1}^{s} \frac{\Omega_{j+1}/2}{z+2 n \omega_f+\tfrac12 \Delta_{j+q}}
			\frac{\Omega_{j-1}^*/2}{z-q \omega_f + (2n-1)\omega_f-\tfrac12 \Delta_{j+q}}
		\right] \right. \\
	&&\left. \times
		\left[ \prod_{n=s}^{q/2 -1} \frac{\Omega_{j+1}/2}{z+(2n+1)\omega_f+\tfrac12 \Delta_{j+q}}
			\frac{\Omega_{j-1}^*/2}{z-q \omega_f +2 n \omega_f - \tfrac12 \Delta_{j+q}}
		\right]
		\right\}.
\end{IEEEeqnarray*}
\end{widetext}
This expression is greatly simplified when applied to an individual resonance. For example the third order interaction for $q=2$ has
\begin{IEEEeqnarray*}{rCl}
	R_{ab}^{(3)} &=& \frac{\Omega_j}{2} \frac{\Omega_{j-1}^*}{2} \frac{\Omega_{j+1}}{2} 
		\left[ \frac{-1}{(\Delta_2+2\omega_f)(\Delta_2+\omega_f)}\right. \IEEEyesnumber \label{eq:q2omegaeff} \\
	&& \left.
		+  \frac{-1}{(\Delta_2+\omega_f)(\Delta_2+2 \omega_f)} +  \frac{-1}{(\Delta_2+\omega_f)(\Delta_2+ \omega_f)} \right],
\end{IEEEeqnarray*}
with $z\approx0$.
On resonance this simplifies to $\Omega_j \Omega_{j-1}^* \Omega_{j+1}/4$.

The effective Hamiltonian for each two dimensional subspace is easily solved to find the dressed states and their energies, which in general are
\begin{equation} \label{eq:Heffdressedenergies}
	E_{\pm} = \frac12 \sqrt{\left[\Delta_q-(R_{bb}-R_{aa})\right]^2+(2 R_{ab})^2}.
\end{equation}
Thus the effective Rabi splitting at the multiphoton resonances is $2 R_{ab}$. 
An atom, initially in its ground state therefore undergo Rabi oscillations
\begin{equation} \label{eq:rabiosc}
	\left \langle \tfrac12 \vert  \psi (t) \right \rangle = 
		\frac{\Omega_{\mathrm{eff}}}{\widetilde \Omega_{\mathrm{eff}}} \sin \left(\widetilde \Omega_{\mathrm{eff}} t/2 \right),
\end{equation}
where $\Omega_{\mathrm{eff}}=2R_{ab}$ and $\widetilde \Omega_{\mathrm{eff}} = \sqrt{\left[\Delta_q-(R_{bb}-R_{aa})\right]^2+\Omega_{\mathrm{eff}}^2}$.
Figure \ref{fig:q2resonance} shows the dressed energies calculated in this way agree well with the numerical approach given in \secref{sec:calculations}.
Figure \ref{fig:q2spectrum} shows the Rabi oscillations of this two-level subsystem calculated from the square modulus of \eqnref{eq:rabiosc}, and for comparison, the same result calculated numerically. As $\Omega_k$ decrease the analytic approximation becomes more accurate. For larger $\Omega_k$ there is a leakage of the population to states outside the subspace considered, which is accounted for in the numerical results. In other work, we find simple analytic expression which accurately describe the time evolution in this case \cite{yuen2018analytic}.

\section{Weak fields and the boundary between classical and quantum electromagnetic fields} \label{sec:quantumclassicalboundary}

\begin{figure*}
	\begin{tabular}{lr}
		\subfloat[Energy Levels \label{fig:alphaplot1}]{
			\includegraphics[width=240pt]{alphaplot1}
			\put(-70,60){$\alpha = 2, \ \omega_0=2$}
		} 
		& 
		\subfloat[Energy Levels \label{fig:g0plot1}]{
		\includegraphics[width=240pt]{g0plot1}
		\put(-90,60){$\omega_0=2, \ g_0/\bar \omega = 10^{-2}$}
		} 
		\\
		\subfloat[Energy Levels \label{fig:alphaplot2}]{
			\includegraphics[width=240pt]{alphaplot2}
			\put(-70,60){$\alpha = 1, \ \omega_0=11$}
		} 
		& 
		\subfloat[Energy Levels \label{fig:g0plot2}]{
			\includegraphics[width=240pt]{g0plot2}
			\put(-90,60){$\omega_0=2, \ g_0/\bar \omega = 10^{0.5}$}
		} 
		\\
				\subfloat[Energy Levels \label{fig:alphaplot3}]{
			\includegraphics[width=240pt]{alphaplot3}
			\put(-70,85){$\alpha = \frac{1}{\sqrt{2}}, \ \omega_0=3$}
		} 
		& 
		\subfloat[Energy Levels \label{fig:g0plot3}]{
			\includegraphics[width=240pt]{g0plot3}
			\put(-90,85){$\omega_0=9, \ g_0/\bar \omega = 10^{0.5}$}
		} 
	\end{tabular}
	\caption{\label{fig:popinversioncollapse} 
Multi-frequency collapse and revivals of a spin-half system driven by two frequencies. Each plot shows the spin-half population inversion calculated in the non-degenerate basis  (orange lines) and Fock basis (blue lines) for comparison. The multi-frequency dynamics are clearly more complex than for a single-frequency coherent driving field. Collapse and revivals are predicted by the non-degenerate formalism because it accounts for the state-dependent action of the field operator, unlike semi-classical approaches such as Floquet theory. Results are shown for a variety of different coherent amplitudes $\alpha_1=\alpha_2=\alpha$, coupling strengths $g_0/\bar \omega$ (normalised to the mean field freqeuncy), and mean field frequencies. The different plots show that as one moves away from the parameter regime $\alpha \sim 1, \ g_0/\bar \omega \sim 1$ or towards higher $\bar \omega_0$ that the accuracy of results calculated in the non-degenerate basis increases. Calculations in this basis are significantly more efficient than in the Fock basis as the system size does not grow exponentially with the number of field modes. A full comparison between the non-degenerate basis and Fock basis (exact results) is given in \figref{fig:l2comparisons}.
Energy is in units of $\omega_f$ which is also the spacing between the two frequency modes, and time is shown in units of $1/g_0$.
}
\end{figure*}

\begin{figure*}
	\subfloat[\label{fig:alphavsf0plot}]{
		\includegraphics[width=250pt]{alphavsf0plot.eps}
		}
	\subfloat[ \label{fig:g0vsf0plt}]{
		\includegraphics[width=250pt]{g0vsf0plt.eps}
		}
	\caption{{\label{fig:l2comparisons}} 
This figure shows there are two regimes of multi-frequency quantum dynamics for a spin-half driven by a two frequency field.
These regimes occur in the parameter space of the mean photon number in each mode ($\abs{\alpha}^2$), mean field frequency $\bar \omega_0$, and interaction strength ($g_0/\bar \omega$) scaled by this frequency.
These regimes are highlighted by the discrepancy (indicated by colour) between the population inversion calculated in the non-degenerate and Fock bases. In region II, the non-degenerate basis produces accurate results as the state of the field remains separable into different modes. In region I the discrepancy occurs because the field modes become entangled, requiring a full description in the Fock basis. The discrepancy was quantified using the square difference between the population inversion results up to time $g_0 t = 1$, as given by \eqnref{eq:shrodrabi}. The colour indicates this discrepancy on log scale with base 10. The boundary between the regions was drawn at a 1\% discrepancy. In (a) $g_0/\bar \omega $ is constant at 0.1, and in (b) $\abs{\alpha}^2 = 5$. Examples of the population inversion at six different points in the parameter space shown here are given in \figref{fig:popinversioncollapse}.
}
\end{figure*}

The electromagnetic fields considered here are in a product of coherent-states, which are the quantum states closest to coherent classical fields. We know however, that they differ from classical fields, and this differenced is most pronounced for weak fields where the mean photon number in each mode is small $\abs{\alpha_k} \sim 1$. For a single-frequency interacting with a spin-half system, a collapse of the Rabi oscillations is observed due dependence on the interaction strength on the square root of the photon number. Some time later revivals occur, but never complete since the frequency with which different pairs of atom-photon states oscillate is irrationally related \cite{Eberly1980}.
Thus, collapse and revivals are one hallmark of a quantised coherent field.
We investigate the limits in which our non-degenerate theory accurately encapsulates the collapse of Rabi oscillations. Remarkably, we find that the our theory accurately predicts the collapse in the deep strong coupling regime in addition to the limit where the coherent fields are highly excited.

When multiple frequencies are present additional mechanisms effect the process of collapse. Additional frequencies change how the different components of oscillating spin-field system dephase. When additional frequencies are far from resonance compares to the interaction strength they drive virtual excitations which shift and the resonantly interacting levels providing further dephasing. 
When the interaction strength becomes comparable to the detuning between the field modes, the spin system is driven by many modes without any single mode dominating its dynamics. 
In turn the different modes of the field become entangled via the spin.
As the interaction strength is increase beyond the mode spacing the multiple modes start to drive the spin collectively, each acting in a similar way.

To investigate the multi-mode dynamics and validate our theoretical formalism we calculate the time evolution for the multi-mode Rabi model numerically, both in the Fock basis and our non-degenerate basis, and compare the time evolution of the spin inversion $\rho_{ee}-\rho_{gg}$. The multimode Rabi model considered is
\begin{IEEEeqnarray}{C} \label{eq:shrodrabi}
	i \partial_t \ket{\psi} = \left[\tfrac12 \omega_0 \sigma_z + \sum_i a_i^{\dagger} a_i \omega_i + g_i (a_i^{\dagger}+a_i) \sigma_x \right] \ket{\psi} \IEEEeqnarraynumspace
\end{IEEEeqnarray}
with coupling strength $g_i = \sqrt{\omega_i/\omega_0} g_0$, and $\omega_0$ is the lowest field frequency which we set resonant spin. Thus the coupling strength grows with $\sqrt{\omega}$ as is the case for the electromagnetic field, and is parameterised by $g_0$. We calculate $\ket{\psi(t)}$ for fields with two three frequencies with mode space $\Delta=1$ in dimensionless units. The initial state is a product of coherent-states, each with equal amplitude $\alpha$, with spin up. We compare $\rho_{ee}-\rho_{gg}$ calculated in the Fock and non-degenerate bases for a range of $g_0, \omega_0 $ and $\alpha$. We quantify the comparison between results in each basis using the $L_2$ measure
\begin{IEEEeqnarray}{rCl} \label{eq:l2measure}
	L_2[f(t),g(t)] = g_0 \int_0^{1/g_0}  \abs{f(t)-g(t)}^2 \mathrm d t.
\end{IEEEeqnarray}
This integral is over the interval $t=0$ to $t=1/g$ - the characteristic time of the collapse for single-frequency Rabi oscillations \cite{puri1986}.

Figure \ref{fig:l2comparisons} compares the spin population inversion calculated in the Fock and non-degenerate basis over a range of field amplitudes $\alpha$, coupling strengths $g_0$, and field frequencies $\bar \omega_0 = \omega_0 + \Delta/2$. Figure \ref{fig:popinversioncollapse} shows typical population inversions for different $\alpha$ and $\omega_0$ with $g_0=\bar \omega/10$ in figures (a) to (c)), and different $g_0=\bar \omega$ and $\omega_0$ with $\abs{\alpha}=5$ in figures (d) to (f).
In all cases there there is close agreement between Fock and non-degenerate basis calculations for times $g_0 t \ll 1$. When $\bar \omega \sim 1$, $\alpha \sim 1$ and $g0/\bar \omega_0 \sim 10^{-1}$ (ultra strong coupling), we find the results are quite different already by $g_0 t \sim 1$. Typical examples are shown in \figref{fig:alphaplot1} and \ref{fig:g0plot1}. We refer to this region of parameter space as region I, and the region outside this as region II. The level of agreement is quantified by their $L_2$ measure which is plotted in \figref{fig:l2comparisons} and the boundary between regions I and II is indicated by the contour $L_2 = 10^{-3}$. For parameters in region II, the agreement becomes remarkably good even for longer times. Figures \ref{fig:alphaplot2}, \ref{fig:g0plot2}, \ref{fig:alphaplot3} and \ref{fig:g0plot3} show the population inversion typical of region II.

The strongest discrepancy between the Fock and non-degenerate basis calculations is for small $\omega_0$ where it is comparable to the frequency difference $\Delta$, small $g_0/\Delta$, and mean photon numbers $\abs{\alpha} \sim 1$. This arises because the non-degenerate basis does not cover Fock space since it represents many non-degenerate Fock states by a single state. The action of the field operators on this state is derived from the action on each underlying Fock state together with the set of amplitudes with which each contributes to the non-degenerate state. The action is correct provided that the field remains in a coherent-state. Gradually, the Fock states which contribute to a non-degenerate state will dephase. The rate of dephasing will be most significant when both field modes interact strongly with the spin but the two modes remain distinct such that their dynamics are different. 


When the dephasing between components of the non-degenerate states is slow the non-degenerate basis yields accurate results.
This is the case when $g_0$ becomes much larger than the mode spacing such that the difference in detuning between the spin and field modes becomes irrelevent. The dephasing is also slows as $\bar \omega_0$ increase for fixed $\Delta$. A perturbative expansion of the interaction shows that there are fewer low order processes which resonantly couple multi-mode Fock states when their frequencies approach.

In the limit of large photon number $\bar n \gg 1$ it is common to ignore the $n$-dependence of the interaction between different Fock states since its variance $\mathrm {var}(n)=\sqrt{\bar n}\ll \bar n$. While this is widely accepted, its universal application is limited to short times much shorted than $1/g_0$ - the timescale of Rabi oscillation collapse and is caused by the $n$-dependence of interactions. For multiple frequencies we give a similar approximation where the ratios $\gamma_N/\gamma_{N\pm k}$ are ignored, the variance of $N$ for this approximation being $\sigma_N$ and not $\sqrt{N}$. Figure \ref{fig:alphavsf0plot} suggests that on a timescale $1/g_0$ the non-degenerate basis is remains suitable for accurate estimation of the spin population in the limit of large $\bar N$, but that this limit is approached slowly, compared with the increase in accuracy of the non-degenerate basis calculations for increasing frequency or coupling strength.

The comparison between the spin-dynamics calculated in the Fock and non-degenerate bases highlight two distinct regimes of multi-frequency light-matter interactions with coherent fields. Regime I: For weak fields and small $g_0/\Delta$ the dephasing between degenerate Fock states becomes important and leads to entanglement between the field modes via the spin. Regime II: For strong fields, or for weak fields and large $g_0/\Delta$, each mode remains in a coherent-state, and although the field and spin become entangled, the entanglement between modes is negligible. In the later regime, the dynamics of the field can be accurately represented in the non-degenerate basis. The semi-classical limit applies only for short times compared to $1/g_0$. The quantum dynamics then emerge for longer times, in one of these two regimes. In regime I the entanglement between the modes require a full quantum picture. Regime II lies between the quantum and classical limits, whereby its dynamics depend on the level of excitation of the field, but the field remains factorizable.

\section{Conclusion} \label{sec:conclusion}

We have developed a formalism which tackles the complexity of multi-frequency atom-photon interactions when each field mode is well described by a coherent-state.
We demonstrate its utility for strong strong fields and remarkably, for a wide range of parameters for weak fields where the quantum mechanical nature of the dynamics are important.

This formalism utilises a subspace of Fock space, spanned by a non-degenerate basis. This basis addresses the issue of level degeneracy encountered when working in the Fock basis. We have derived the effect of the photon creation and annihilation operators on these basis states including the asymmetry between them associated with quantum fluctuations. 
For highly excited coherent-states, where the photon number in each mode is large compared to unity, we show these fluctuations are small and can be ignored.
We have demonstrated that by working in this basis, accurate calculations of the multi-frequency dressed states can be performed. Applying the resolvent formalism in this basis we have derived effective Hamiltonians which includes analytic expression for the level shifts and multi-photon interactions for an atom in a multi-frequency field.

For weak fields, where quantum effects are significant, we have compared the long time dynamics predicted by our formalism to exact calculations in the Fock basis. By investigating the parameter range for which our formalism is applicable we have identified two different quantum regimes for the collapse and revival of Rabi oscillations of a spin-half system in a multi-frequency field. One of these covers a large region in parameter space where some quantum aspects of the field are essential for the collapse, but the fields do not become entangled. Our formalism is uniquely capable of addressing problems in this regime.

\begin{acknowledgments}
This work was supported through the UK National Quantum Technologies Programme (NQIT hub, EP/M013243/1) and the EU H2020 Collaborative project QuProCS (Grant Agreement No. 641277). 
The authors would like to thank Axel Kuhn, Candadi Sukumar, Kathrin Luksch and Rian Hughes for helpful discussions during the preparation of this work, and Ezra Kassa for comments on the manuscript.
\end{acknowledgments}



