Kinetic equations for two-photon light in random media

The propagation of light in random media is usually considered within the setting of classical optics [1]. However, there has been considerable recent interest in phenomena where quantum effects play a key role [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. Of particular importance is understanding the impact of multiple scattering on entangled two-photon states, with an eye towards characterizing the transfer of entanglement from the field to matter [13, 16, 5, 17, 18, 19, 20]. Progress in this direction can be expected to lead to advances in spectroscopy [21], imaging [22, 23, 24, 25, 26, 27, 28, 29, 30, 31] and communications [32, 33, 34] .

The propagation of two-photon light is generally investigated either in free space or, in some cases, with account of diffraction [35, 36] or scattering [37]. In this paper, we consider the propagation of two-photon light in a random medium. A step in this direction was taken in [38], where a model in which the field is quantized and the medium is treated classically was investigated. The main drawback of that work is that it is does not allow for the transfer of entanglement between the field and the atoms or between the atoms themselves. Instead, we treat the problem from first principles, employing a model in which the field and the matter are both quantized. We show that for a medium consisting of two-level atoms, the field and atomic probability amplitudes for a two-photon state obey a system of nonlocal partial differential equations with random coefficients. Using this result, we find that at long times and large distances, the corresponding average probability densities in a random medium can be determined from the solutions to a system of kinetic equations. These equations follow from the multiscale asymptotics of the average Wigner transform of the amplitudes in a suitable high-frequency limit [39, 40, 41, 42]. This formulation of the problem builds on earlier research by the authors on collective emission of single photons in a random medium of two-level atoms [43]. In that work, we employed a formulation of quantum electrodynamics in which the field is quantized in real space, thus allowing for the field and atomic degrees of freedom to be treated on the same footing.

This paper is organized as follows. In section 2 we formulate our model for the propagation of two-photon light and derive the equations governing the dynamics of the field and atomic degrees of freedom. Section 3 is concerned with the application of the real-space formalism to the problem of stimulated emission by a single atom. In section 4 we study the dynamics of a two-photon state in a homogeneous medium. In section 5 we discuss the problem of stimulated emission in a random medium and obtain the governing kinetic equations. Section 6 takes up the general problem of two-photon transport in random media and presents the relevant kinetic equations. Our conclusions are formulated in section 7. The appendices contain the details of calculations.

We consider the following model for the interaction between a quantized massless scalar field and a system of two-level atoms. The atoms are taken to be identical, stationary and sufficiently well separated that interatomic interactions can be neglected. The overall system is described by the Hamiltonian

$$H=H_{F}+H_{A}+H_{I}\ ,$$

(1)

where $H_{F}$ is the Hamiltonian of the field, $H_{A}$ is the Hamiltonian of the atoms and $H_{I}$ is the interaction Hamiltonian. In order to treat the atoms and the field on the same footing, it is useful to introduce a real-space representation of $H$ [43]. The Hamiltonian of the field is of the form

$\displaystyle H_{F}=\hbar c\int d^{3}x(-\Delta)^{1/2}\phi^{\dagger}({\bf x})\phi({\bf x})\ ,$

(2)

where $\phi$ is a Bose scalar field that obeys the commutation relations

	$\displaystyle[\phi({\bf x}),\phi^{\dagger}({\bf x}^{\prime})]$	$\displaystyle=\delta({\bf x}-{\bf x}^{\prime})\ ,$		(3)
	$\displaystyle[\phi({\bf x}),\phi({\bf x}^{\prime})]$	$\displaystyle=0\ .$		(4)

where we have neglected the zero-point energy. The nonlocal operator $(-\Delta)^{1/2}$ is defined by the Fourier integral

	$\displaystyle(-\Delta)^{1/2}f({\bf x})$	$\displaystyle=\int\frac{d^{3}k}{(2\pi)^{3}}e^{i{\bf k}\cdot{\bf x}}|{\bf k}|\tilde{f}({\bf k})\ ,$		(5)
	$\displaystyle\tilde{f}({\bf k})$	$\displaystyle=\int d^{3}xe^{-i{\bf k}\cdot{\bf x}}f({\bf x})\ .$		(6)

We note that (2) is equivalent to the usual oscillator representation of $H_{F}$.

To facilitate the treatment of random media, it will prove convenient to introduce a continuum model of the atomic degrees of freedom [43]. The Hamiltonian of the atoms is given by

$\displaystyle H_{A}=\hbar\Omega\int d^{3}x\rho({\bf x})\sigma^{\dagger}({\bf x})\sigma({\bf x})\ ,$

(7)

where $\Omega$ is the atomic resonance frequency, $\rho({\bf x})$ is the number density of the atoms, and $\sigma$ is a Fermi field that obeys the anticommutation relations

	$\displaystyle\{\sigma({\bf x}),\sigma^{\dagger}({\bf x}^{\prime})\}$	$\displaystyle=\frac{1}{\rho({\bf x})}\delta({\bf x}-{\bf x}^{\prime})\ ,$		(8)
	$\displaystyle\{\sigma({\bf x}),\sigma({\bf x}^{\prime})\}$	$\displaystyle=0\ .$		(9)

The Hamiltonian describing the interaction between the field and the atoms is taken to be

$\displaystyle H_{I}=\hbar{g}\int d^{3}x\rho({\bf x})\left(\phi^{\dagger}({\bf x})\sigma({\bf x})+\sigma^{\dagger}({\bf x})\phi({\bf x})\right)\ ,$

(10)

where ${g}$ is the strength of the atom-field coupling. Here we have made the Markovian approximation, in which the coupling constant is independent of the frequency of the photons or positions of the atoms, and have imposed the rotating wave approximation (RWA).

We suppose that the system is in a two-photon state of the form

	$\displaystyle|\Psi\rangle=\int d^{3}x_{1}d^{3}x_{2}$	$\displaystyle\left(\psi_{2}({\bf x}_{1},{\bf x}_{2},t)\phi^{\dagger}({\bf x}_{1})\phi^{\dagger}({\bf x}_{2})+\psi_{1}({\bf x}_{1},{\bf x}_{2},t)\rho({\bf x}_{1})\sigma^{\dagger}({\bf x}_{1})\phi^{\dagger}({\bf x}_{2})\right.$		(11)
	$\displaystyle+$	$\displaystyle\left.a({\bf x}_{1},{\bf x}_{2},t)\rho({\bf x}_{1})\rho({\bf x}_{2})\sigma^{\dagger}({\bf x}_{1})\sigma^{\dagger}({\bf x}_{2})\right)|0\rangle\ ,$

where $|0\rangle$ is the combined vacuum state of the field and the ground state of the atoms. Here the atomic amplitude $a({\bf x}_{1},{\bf x}_{2},t)$ is the probability amplitude for exciting two atoms at the points ${\bf x}_{1}$ and ${\bf x}_{2}$ at time $t$, the one-photon amplitude $\psi_{1}({\bf x}_{1},{\bf x}_{2},t)$ is the probability amplitude for exciting an atom at ${\bf x}_{1}$ and creating a photon at ${\bf x}_{2}$, and the two-photon amplitude $\psi_{2}({\bf x}_{1},{\bf x}_{2},t)$ is the probability amplitude for creating photons at ${\bf x}_{1}$ and ${\bf x}_{2}$. The functions $\psi_{2}$ and $a$ are symmetric and antisymmetric, respectively:

$$\psi_{2}({\bf x}_{2},{\bf x}_{1},t)=\psi_{2}({\bf x}_{1},{\bf x}_{2},t)\ ,\quad a({\bf x}_{2},{\bf x}_{1},t)=-a({\bf x}_{1},{\bf x}_{2},t)\ ,$$

(12)

consistent with the bosonic and fermionic character of the corresponding fields.

The state $|\Psi\rangle$ is the most general two-photon state within the RWA. In addition, $|\Psi\rangle$ is normalized so that $\langle\Psi|\Psi\rangle=1$. It follows from (3) and (8) that the amplitudes obey the normalization condition

$\displaystyle\int d^{3}x_{1}d^{3}x_{2}\left(2|\psi_{2}({\bf x}_{1},{\bf x}_{2},t)|^{2}+\rho({\bf x}_{1})|\psi_{1}({\bf x}_{1},{\bf x}_{2},t)|^{2}+2\rho({\bf x}_{1})\rho({\bf x}_{2})|a({\bf x}_{1},{\bf x}_{2},t)|^{2}\right)=1\ .$

(13)

If the amplitudes $a({\bf x}_{1},{\bf x}_{2},t)$, $\psi_{1}({\bf x}_{1},{\bf x}_{2},t)$ and $\psi_{2}({\bf x}_{1},{\bf x}_{2},t)$ are factorizable as functions of ${\bf x}_{1}$ and ${\bf x}_{2}$, then there are no quantum correlations and the state $|\Psi\rangle$ is not entangled. Otherwise $|\Psi\rangle$ is entangled. If $\psi_{2}$ alone is not factorizable, we say that $|\Psi\rangle$ is an entangled two-photon state.

The dynamics of $|\Psi\rangle$ is governed by the Schrodinger equation

$\displaystyle i\hbar\partial_{t}|\Psi\rangle=H|\Psi\rangle\ .$

(14)

Projecting onto the states $\phi^{\dagger}({\bf x})|0\rangle$ and $\sigma^{\dagger}({\bf x})|0\rangle$ and making use of (3) and (8), we arrive at the following system of equations obeyed by $a$, $\psi_{1}$ and $\psi_{2}$:

		$\displaystyle i\partial_{t}\psi_{2}({\bf x}_{1},{\bf x}_{2},t)=c(-\Delta_{{\bf x}_{1}})^{1/2}\psi_{2}({\bf x}_{1},{\bf x}_{2},t)+c(-\Delta_{{\bf x}_{2}})^{1/2}\psi_{2}({\bf x}_{1},{\bf x}_{2},t)$
		$\displaystyle+\frac{g}{2}(\rho({\bf x}_{1})\psi_{1}({\bf x}_{1},{\bf x}_{2},t)+\rho({\bf x}_{2})\psi_{1}({\bf x}_{2},{\bf x}_{1},t))\ ,$		(15)
		$\displaystyle\rho({\bf x}_{1})i\partial_{t}\psi_{1}({\bf x}_{1},{\bf x}_{2},t)=2{g}\rho({\bf x}_{1})\psi_{2}({\bf x}_{1},{\bf x}_{2},t)+\rho({\bf x}_{1})\left[c(-\Delta_{{\bf x}_{2}})^{1/2}+\Omega\right]\psi_{1}({\bf x}_{1},{\bf x}_{2},t)$
		$\displaystyle-2{g}\rho({\bf x}_{1})\rho({\bf x}_{2})a({\bf x}_{1},{\bf x}_{2},t)\ ,$		(16)
		$\displaystyle\rho({\bf x}_{1})\rho({\bf x}_{2})i\partial_{t}a({\bf x}_{1},{\bf x}_{2},t)=\frac{g}{2}\rho({\bf x}_{1})\rho({\bf x}_{2})(\psi_{1}({\bf x}_{2},{\bf x}_{1},t)-\psi_{1}({\bf x}_{1},{\bf x}_{2},t))+2\Omega\rho({\bf x}_{1})\rho({\bf x}_{2})\,a({\bf x}_{1},{\bf x}_{2},t)\ .$		(17)

The overall factors of $\rho({\bf x})$ in the above will be cancelled as necessary. The details of the calculation are presented in Appendix A.

The above model views the atomic degrees of freedom as fermionic. In Appendix D we consider the bosonic case.

In this section we consider the problem of stimulated emission by a single atom. We assume that the atom is located at the origin and put $\rho({\bf x})=\delta({\bf x})$. Thus the system (2)–(17) becomes

	$\displaystyle i\partial_{t}\psi_{2}({\bf x}_{1},{\bf x}_{2},t)=c(-\Delta_{{\bf x}_{1}})^{1/2}\psi_{2}({\bf x}_{1},{\bf x}_{2},t)+c(-\Delta_{{\bf x}_{2}})^{1/2}\psi_{2}({\bf x}_{1},{\bf x}_{2},t)$
	$\displaystyle+\frac{g}{2}(\delta({\bf x}_{1})\psi_{1}({\bf x}_{1},{\bf x}_{2},t)+\delta({\bf x}_{2})\psi_{1}({\bf x}_{2},{\bf x}_{1},t))\ ,$		(18)
	$\displaystyle\delta({\bf x}_{1})i\partial_{t}\psi_{1}({\bf x}_{1},{\bf x}_{2},t)=2{g}\delta({\bf x}_{1})\psi_{2}({\bf x}_{1},{\bf x}_{2},t)+\delta({\bf x}_{1})\left[c(-\Delta_{{\bf x}_{2}})^{1/2}+\Omega\right]\psi_{1}({\bf x}_{1},{\bf x}_{2},t)\ ,$		(19)

where the term $\delta({\bf x}_{1})\delta({\bf x}_{2})a({\bf x}_{1},{\bf x}_{2},t)$ does not contribute due to the antisymmetry of $a$. We assume that initially there is only one photon present in the system, which means that the initial conditions for the amplitudes $\psi_{1}$ and $\psi_{2}$ are given by

	$\displaystyle\psi_{1}(0,{\bf x},0)$	$\displaystyle=e^{i{\bf k}_{0}\cdot{\bf x}}\ ,$
	$\displaystyle\psi_{2}({\bf x}_{1},{\bf x}_{2},0)$	$\displaystyle=0\ ,$

where ${\bf k}_{0}$ is the wavevector of the photon. Note that the amplitude of $\psi_{1}$ is set to unity for convenience. To proceed, we take the Laplace transform with respect to $t$ and the Fourier transforms with respect to ${\bf x}_{1}$ and ${\bf x}_{2}$ of (18) and (19). We thus obtain

	$\displaystyle iz\psi_{2}({\bf k}_{1},{\bf k}_{2},z)$	$\displaystyle=\left[c|{\bf k}_{1}|+c|{\bf k}_{2}|\right]\psi_{2}({\bf k}_{1},{\bf k}_{2})+\frac{g}{2}\left[\psi_{1}({\bf k}_{1},z)+\psi_{1}({\bf k}_{2},z)\right]\ ,$		(20)
	$\displaystyle i\left[z\psi_{1}({\bf k},z)-\delta({\bf k}-{\bf k}_{0})\right]$	$\displaystyle=2{g}\int d^{3}k^{\prime}\psi_{2}({\bf k}^{\prime},{\bf k},z)+\left[c|{\bf k}|+\Omega\right]\psi_{1}({\bf k},z)\ .$		(21)

Here we have defined the Laplace transform by

$\displaystyle a(z)=\int_{0}^{\infty}dte^{-zt}a(t)\ ,$

(22)

where we denote a function and its Laplace transform by the same symbol. Solving (20) and (21) leads to an integral equation for $\psi_{1}({\bf k},z)$ of the form

$\displaystyle\psi_{1}({\bf k},z)=\frac{\delta({\bf k}-{\bf k}_{0})}{z+i(c|{\bf k}|+\Omega)-i\Sigma({\bf k},z)}+\frac{i{g}^{2}}{z+i(c|{\bf k}|+\Omega)-i\Sigma({\bf k},z)}\int d^{3}k^{\prime}\frac{\psi_{1}({\bf k}^{\prime},z)}{(c|{\bf k}^{\prime}|+c|{\bf k}|)-iz}\ ,$

(23)

where

$\displaystyle\Sigma({\bf k},z)=g^{2}\int d^{3}k^{\prime}\frac{1}{(c|{\bf k}^{\prime}|+c|{\bf k}|)-iz-i\epsilon}\ ,$

(24)

with $\epsilon\to 0^{+}$. In order to evaluate the integral in (23), we make the pole approximation where we replace $\Sigma({\bf k},z)$ with $\Sigma({\bf k},-i(\Omega+c|{\bf k}|))$. We note that this quantity is independent of ${\bf k}$ and $z$, and so we will denote it by $\Sigma$. For consistency, we also replace $z$ by $-i(\Omega+c|{\bf k}|)$ under the integral in (23). Note that this approximation arises in the Wigner-Weisskopf theory of spontaneous emission. In addition, we split $\Sigma$ into its real and imaginary parts:

	$\displaystyle\mathrm{Re}\,\Sigma$	$\displaystyle=\delta\omega\ ,$		(25)
	$\displaystyle\mathrm{Im}\,\Sigma$	$\displaystyle=\Gamma/2\ .$		(26)

We can calculate $\Gamma$ and $\delta\omega$ by making use of the identity

$\displaystyle\frac{1}{c|{\bf k}|-\Omega-i\epsilon}=P\frac{1}{c|{\bf k}|-\Omega}+i\pi\delta(c|{\bf k}|-\Omega)\ .$

(27)

We find that

	$\displaystyle\Gamma$	$\displaystyle=2g^{2}\pi\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}}\,\delta(c|{\bf k}|-\Omega)$		(28)
		$\displaystyle=\frac{g^{2}\Omega^{2}}{\pi c^{3}}\ ,$		(29)

and

$$\delta\omega=\frac{g^{2}}{2\pi^{2}}\int_{0}^{2\pi/\Lambda}\frac{k^{2}dk}{ck-\Omega}\ ,$$

(30)

where we have introduced a high-frequency cutoff to regularize the divergent integral. Using the above results, (23) becomes

$\displaystyle\psi_{1}({\bf k},z)=\frac{\delta({\bf k}-{\bf k}_{0})}{z+i(c|{\bf k}|+\Omega)-i\Sigma}+\frac{ig^{2}}{z+i(c|{\bf k}|+\Omega)-i\Sigma}\int d^{3}k^{\prime}\frac{\psi_{1}({\bf k}^{\prime},z)}{c|{\bf k}^{\prime}|-\Omega}\ .$

(31)

Since $\psi_{1}$ appears on both the left- and right-hand sides of (31), upon iteration we obtain an infinite series for $\psi_{1}$, which to order $O(g^{2})$ is given by

	$\displaystyle\psi_{1}({\bf k},z)$	$\displaystyle=\frac{\delta({\bf k}-{\bf k}_{0})}{z+i(c|{\bf k}|+\Omega)-i\Sigma}+\frac{ig^{2}}{z+i(\Omega-c|{\bf k}^{\prime}|)-i\Sigma}\int d^{3}k^{\prime}\frac{\delta({\bf k}^{\prime}-{\bf k}_{0})}{(\Omega-c|{\bf k}^{\prime}|)(z+i(c|{\bf k}^{\prime}|+\Omega)-i\Sigma)}$		(32)
		$\displaystyle=\frac{\delta({\bf k}-{\bf k}_{0})}{z+i(c|{\bf k}|+\Omega)-i\Sigma}+\frac{ig^{2}}{z+i(c|{\bf k}|+\Omega)-i\Sigma}\frac{1}{(c|{\bf k}_{0}|-\Omega)(z+i(c|{\bf k}_{0}|+\Omega)-i\Sigma)}\ .$		(33)

Inverting the Laplace transform, we find that $\psi_{1}({\bf k},t)$ is given by

$\displaystyle\psi_{1}({\bf k},t)$

$\displaystyle=e^{-i(\Omega-\delta\omega)t-\Gamma t/2}\left[\delta({\bf k}-{\bf k}_{0})e^{-ic|{\bf k}|t}+\frac{i{g}^{2}}{((c|{\bf k}_{0}|-\Omega)(c|{\bf k}_{0}|-c|{\bf k}|)}\left(e^{-ic|{\bf k}|t}-e^{-ic|{\bf k}_{0}|t}\right)\right]\ .$

(34)

Note that $\psi_{1}$ decays exponentially at long times ($\Gamma t\gg 1$) and that $\psi_{2}$ can be obtained from (20).

In this section we consider the case of a homogeneous medium and set $\rho({\bf x})=\rho_{0}$, where $\rho_{0}$ is constant. It is useful to define the function $\tilde{\psi_{1}}({\bf x}_{1},{\bf x}_{2},t)=\psi_{1}({\bf x}_{2},{\bf x}_{1},t)$ and write (17) as a $4\times 4$ symmetric system. If we further define the vector

$\displaystyle{\bf\Psi}({\bf x}_{1},{\bf x}_{2},t)=\begin{bmatrix}\sqrt{2}\psi_{2}({\bf x}_{1},{\bf x}_{2},t)\\ \sqrt{\rho_{0}/2}\psi_{1}({\bf x}_{1},{\bf x}_{2},t)\\ \sqrt{\rho_{0}/2}\tilde{\psi_{1}}({\bf x}_{1},{\bf x}_{2},t)\\ \sqrt{2}\rho_{0}a({\bf x}_{1},{\bf x}_{2},t)\end{bmatrix}\ ,$

(35)

then (17) can be written as

$\displaystyle i\partial_{t}{\bf\Psi}=A{\bf\Psi}\ ,$

(36)

where

$\displaystyle A=\begin{bmatrix}c(-\Delta_{{\bf x}_{1}})^{1/2}+c(-\Delta_{{\bf x}_{2}})^{1/2}&{g}\sqrt{\rho_{0}}&{g}\sqrt{\rho_{0}}&0\\ {g}\sqrt{\rho_{0}}&c(-\Delta_{{\bf x}_{2}})^{1/2}+\Omega&0&-{g}\sqrt{\rho_{0}}\\ {g}\sqrt{\rho_{0}}&0&c(-\Delta_{{\bf x}_{1}})^{1/2}+\Omega&{g}\sqrt{\rho_{0}}\\ 0&-{g}\sqrt{\rho_{0}}&{g}\sqrt{\rho_{0}}&2\Omega\end{bmatrix}.$

(37)

This definition of ${\bf\Psi}$ has the advantage that the matrix $A$ is symmetric. The Fourier transform of (36) with respect to the variables ${\bf x}_{1}$ and ${\bf x}_{2}$ is given by

$\displaystyle i\partial_{t}\hat{{\bf\Psi}}=\hat{A}\hat{{\bf\Psi}}\,$

(38)

where

$\displaystyle\hat{A}({\bf k}_{1},{\bf k}_{2})=\begin{bmatrix}c|{\bf k}_{1}|+c|{\bf k}_{2}|&{g}\sqrt{\rho_{0}}&{g}\sqrt{\rho_{0}}&0\\ {g}\sqrt{\rho_{0}}&c|{\bf k}_{2}|+\Omega&0&-{g}\sqrt{\rho_{0}}\\ {g}\sqrt{\rho_{0}}&0&c|{\bf k}_{1}|+\Omega&{g}\sqrt{\rho_{0}}\\ 0&-{g}\sqrt{\rho_{0}}&{g}\sqrt{\rho_{0}}&2\Omega\end{bmatrix}.$

(39)

The matrix $\hat{A}$ has eigenvalues

	$\displaystyle\lambda_{1}$	$\displaystyle=b_{1}-\frac{\sqrt{b_{2}-2\sqrt{b_{3}}}}{2}\,,$
	$\displaystyle\lambda_{2}$	$\displaystyle=b_{1}+\frac{\sqrt{b_{2}-2\sqrt{b_{3}}}}{2}\,,$
	$\displaystyle\lambda_{3}$	$\displaystyle=b_{1}-\frac{\sqrt{b_{2}+2\sqrt{b_{3}}}}{2}\,,$
	$\displaystyle\lambda_{4}$	$\displaystyle=b_{1}+\frac{\sqrt{b_{2}+2\sqrt{b_{3}}}}{2}\,,$

where

	$\displaystyle b_{1}$	$\displaystyle=\Omega+\frac{c\,|{\bf k}_{1}|}{2}+\frac{c\,|{\bf k}_{2}|}{2}\,,$		(40)
	$\displaystyle b_{2}$	$\displaystyle=2\,\Omega^{2}+8\,g^{2}\rho_{0}+c^{2}\,{|{\bf k}_{1}|}^{2}+c^{2}\,{|{\bf k}_{2}|}^{2}-2\,\Omega\,c\,|{\bf k}_{1}|-2\,\Omega\,c\,|{\bf k}_{2}|\,,$		(41)
	$\displaystyle b_{3}$	$\displaystyle=\Omega^{4}-2\,\Omega^{3}\,c\,|{\bf k}_{1}|-2\,\Omega^{3}\,c\,|{\bf k}_{2}|+\Omega^{2}\,c^{2}\,{|{\bf k}_{1}|}^{2}+4\,\Omega^{2}\,c^{2}\,|{\bf k}_{1}|\,|{\bf k}_{2}|+\Omega^{2}\,c^{2}\,{|{\bf k}_{2}|}^{2}$
		$\displaystyle+8\,\Omega^{2}\,g^{2}\rho_{0}-2\,\Omega\,c^{3}\,{|{\bf k}_{1}|}^{2}\,|{\bf k}_{2}|-2\,\Omega\,c^{3}\,|{\bf k}_{1}|\,{|{\bf k}_{2}|}^{2}-8\,\Omega\,c\,g^{2}\rho_{0}\,|{\bf k}_{1}|$
		$\displaystyle-8\,\Omega\,c\,g^{2}\rho_{0}\,|{\bf k}_{2}|+c^{4}\,{|{\bf k}_{1}|}^{2}\,{|{\bf k}_{2}|}^{2}+4\,c^{2}\,g^{2}\rho_{0}\,{|{\bf k}_{1}|}^{2}+4\,c^{2}\,g^{2}\rho_{0}\,{|{\bf k}_{2}|}^{2}.$		(42)

The components of the associated eigenvectors ${\bf v}_{i}({\bf k}_{1},{\bf k}_{2})$ are given by

	$\displaystyle v_{i1}$	$\displaystyle=-\frac{2\,\Omega^{3}-2\,\Omega\,g^{2}\rho_{0}+2\,\Omega^{2}\,c\,|{\bf k}_{1}|+2\,\Omega^{2}\,c\,|{\bf k}_{2}|-c\,g^{2}\rho_{0}\,|{\bf k}_{1}|-c\,g^{2}\rho_{0}\,|{\bf k}_{2}|+2\,\Omega\,c^{2}\,|{\bf k}_{1}|\,|{\bf k}_{2}|}{g^{2}\rho_{0}\,\left(c\,|{\bf k}_{1}|-c\,|{\bf k}_{2}|\right)}$
		$\displaystyle-\frac{\left(4\,\Omega+c\,|{\bf k}_{1}|+c\,|{\bf k}_{2}|\right)\lambda_{i}^{2}}{g^{2}\rho_{0}\,\left(c\,|{\bf k}_{1}|-c\,|{\bf k}_{2}|\right)}+\frac{\lambda_{i}\left(5\,\Omega^{2}-2\,g^{2}\rho_{0}+3\,\Omega\,c\,|{\bf k}_{1}|+3\,\Omega\,c\,|{\bf k}_{2}|+c^{2}\,|{\bf k}_{1}|\,|{\bf k}_{2}|\right)}{g^{2}\rho_{0}\,\left(c\,|{\bf k}_{1}|-c\,|{\bf k}_{2}|\right)}$
		$\displaystyle+\frac{\lambda_{i}^{3}}{g^{2}\rho_{0}\,\left(c\,|{\bf k}_{1}|-c\,|{\bf k}_{2}|\right)}\,,$		(43)
	$\displaystyle v_{i2}$	$\displaystyle=\frac{\lambda_{i}^{2}}{g\sqrt{\rho_{0}}\,\left(c\,|{\bf k}_{1}|-c\,|{\bf k}_{2}|\right)}+\frac{2\,\left(\Omega^{2}+c\,|{\bf k}_{1}|\,\Omega-g^{2}\rho_{0}\right)}{g\sqrt{\rho_{0}}\,\left(c\,|{\bf k}_{1}|-c\,|{\bf k}_{2}|\right)}-\frac{\left(3\,\Omega+c\,|{\bf k}_{1}|\right)\,\left(\lambda_{i}\right)}{g\sqrt{\rho_{0}}\,\left(c\,|{\bf k}_{1}|-c\,|{\bf k}_{2}|\right)}\,,$		(44)
	$\displaystyle v_{i3}$	$\displaystyle=\frac{\lambda_{i}^{2}}{g\sqrt{\rho_{0}}\,\left(c\,|{\bf k}_{1}|-c\,|{\bf k}_{2}|\right)}+\frac{2\,\left(\Omega^{2}+c\,|{\bf k}_{2}|\,\Omega-g^{2}\rho_{0}\right)}{g\sqrt{\rho_{0}}\,\left(c\,|{\bf k}_{1}|-c\,|{\bf k}_{2}|\right)}-\frac{\left(3\,\Omega+c\,|{\bf k}_{2}|\right)\,\left(\lambda_{i}\right)}{g\sqrt{\rho_{0}}\,\left(c\,|{\bf k}_{1}|-c\,|{\bf k}_{2}|\right)}\,,$		(45)
	$\displaystyle v_{i4}$	$\displaystyle=1.$		(46)

It follows that the solution to (38) is

$\displaystyle{\bf\Psi}({\bf x}_{1},{\bf x}_{2},t)=\sum_{i=1}^{4}\int\frac{d^{3}k_{1}}{(2\pi)^{3}}\frac{d^{3}k_{2}}{(2\pi)^{3}}e^{i{\bf x}_{1}\cdot{\bf k}_{1}+i{\bf x}_{2}\cdot{\bf k}_{2}}C_{i}({\bf k}_{1},{\bf k}_{2})e^{-i{g}\sqrt{\rho_{0}}\lambda_{i}({\bf k}_{1},{\bf k}_{2})t}{\bf v}_{i}({\bf k}_{1},{\bf k}_{2})\,,$

(47)

where the values $C_{i}$ are solutions to the linear system

$\displaystyle\hat{{\bf\Psi}}({\bf k}_{1},{\bf k}_{2},0)=\sum_{i=1}^{4}{\bf v}_{i}({\bf k}_{1},{\bf k}_{2})C_{i}({\bf k}_{1},{\bf k}_{2}).$

(48)

In order to study the emission of photons, we assume that initially there are two localized volumes of excited atoms of linear size $l_{s}$ centered at the points ${\bf r}_{1}$ and ${\bf r}_{2}$. The initial amplitudes are taken to be

	$\displaystyle\psi_{2}({\bf x}_{1},{\bf x}_{2},0)$	$\displaystyle=0\,,$		(49)
	$\displaystyle\psi_{1}({\bf x}_{1},{\bf x}_{2},0)$	$\displaystyle=0\,,$		(50)
	$\displaystyle a({\bf x}_{1},{\bf x}_{2},0)$	$\displaystyle=\left(\frac{1}{\pi l_{s}^{2}}\right)^{3/2}\left(e^{-|{\bf x}_{1}-{\bf r}_{1}|^{2}/2l_{s}^{2}}e^{-|{\bf x}_{2}-{\bf r}_{2}|^{2}/2l_{s}^{2}}-e^{-|{\bf x}_{2}-{\bf r}_{1}|^{2}/2l_{s}^{2}}e^{-|{\bf x}_{1}-{\bf r}_{2}|^{2}/2l_{s}^{2}}\right).$		(51)

Figure 1 illustrates the time dependence of $\psi_{1},\psi_{2}$ and $a$, where we have set the dimensionless quantities ${\Omega}/{\sqrt{\rho_{0}}{g}}={c}/{l_{s}\sqrt{\rho_{0}}{g}}=1$. Note that the atomic amplitude is antisymmetric, consistent with (8) and that it corresponds to an entangled state. We observe that the amplitudes oscillate and decay in time. Moreover, the atomic and one-photon amplitudes are out of phase with one another, and the two-photon amplitude is an order of magnitude smaller.