Simulation technique of quantum optical emission process from multiple two-level atoms based on classical numerical method

In this paper, we report a numerical method for analyzing optical radiation from a two-level atom. The proposed method can consistently consider the optical emission and absorption process of an atom, and also the interaction between atoms through their interaction with a radiation field. The numerical model is based on a damping oscillator description of a dipole current, which is a classical model of atomic transition and is implemented with a finite-difference time-domain method. Using the method, we successfully simulate the spontaneous emission phenomena in a vacuum, where the interaction between an atom and a radiated field plays an important role. We also simulate the radiation from an atom embedded in a photonic crystal (PhC) cavity. As a result, an atom-cavity field interaction is sucessfuly incorporated in the simulation, and the enhancement of the optical emission rate of an excited atom is explained. The method considers the effect of the interaction between atoms through the radiated field. We simulate the optical emission process of the multiple atoms and show that an enhancement of the emission rate can occur owing to the an atom-atom interaction (superradiance)(R. H. Dicke, Phys. Rev. {\bf 93}, 99[1954]). We also show that the emission rate is suppressed by the effect of the destructive dipole-dipole interaction under an out-of-phase excitation condition (subradiance).


Introduction
The spontaneous emission rate is a fundamental physical quantity that governs light-matter interaction. It is well known that the radiation process of an atom depends on its environment [1]. This means that the radiation process can be controlled artificially by changing the environment. Many research studies have been carried out aiming at the control of the radiation process in combination with a micro optical resonator, such as PhC cavities [2] and whispering-gallery cavities [3].
The finite-difference time-domain (FDTD) method solving Maxwell's equation has become a powerful tool in computational electrodynamics. The method has been applied to studies of the optical process in actual photonic nanostructures. Recently, an FDTD method has been applied to phenomena including the spontaneous emission process of an excited two-level atom located in an optical cavity, which is called cavity quantum electrodynamics (cQED) [4,5]. However, their numerical model adopts dipole current decay in time at a given a priori rate. In their calculation, the decay rate was chosen to match the spontaneous emission rate of an atom in a vacuum. Due to their formulation, which is based on a dipole model assuming attenuation beforehand, their application is limited. It cannot be applied to pure spontaneous emission phenomenon in a vacuum.
Recently, we reported the numerical technique based on FDTD method to analyze optical process [6]. In the report, we successfully simulated an emission process in an optical cavity.
However, the numerical model used in the calculation assumed a non-radiative decay term. In this paper, we propose a numerical method that can analyze a spontaneous emission process even in a vacuum. A two-level atom is modeled based on an equation of motion for dipole oscillation, which radiates an electromagnetic field and interacts with the radiated electromagnetic field, which means a consideration of the absorption process of the electromagnetic field. A dipole exchanges electromagnetic field with another dipole by absorbing the radiated field, which leads to an dipole-dipole interaction. The method does not require a priori knowledge of the decay rate of the dipole current caused by the spontaneous emission in a vacuum. A spontaneous emission is purely quantum mechanical phenomenon. However, a classical oscillating current dipole radiates an electromagnetic field and is decelerated by a radiated field, which, in a sense, constitutes a classical description of quantum phenomena. This model makes no assumptions about the spontaneous emission rate of an dipole in a vacuum. Therefore, the spontaneous emission process can be consistently analyzed using this model. Using this model, we calculate the dipole current decay in a vacuum, where the calculated value of the spontaneous emission rate agrees well with that of the theory. Rabi oscillation and the enhancement of the emission rate of a dipole embedded in a PhC cavity are also obtained. The enhancement of the optical emission rate caused by coupling with the cavity's optical mode is successfully simulated. The advantage of our model is that it incorporates the absorption of the field as well as the emission automatically. This means the dipole-dipole interaction is included as a fundamental process of the formulation [7]. The radiation process with a dipole-dipole interaction can also be calculated using this method.
A numerical method that describes the time evolution of an electromagnetic field with a dipole is explained in Sec. 2. Using the method, we present the calculated result of spontaneous emission in a vacuum in Sec. 3. Calculated results for the radiation from two excited dipoles and the effect of the dipole-dipole interaction between multiple dipoles are also described. The aim of this study is to analyze the radiation process of a photon from an excited dipole. In this paper, we assume an atom with two energy levels and model its atomic transition as an oscillating dipole (Fig. 1). An oscillating dipole radiates an electromagnetic field, which can be expressed by the current term in Maxwell's equation. The time evolution of the electromagnetic field and dipole is expressed by the following Maxwell equation,

Numerical model of spontaneous emission
where the generation of a magnetic field by an oscillating dipole current is included. The equation of motion of the dipole current is expressed as where ω p is the polarization frequency, and ∆ε = ε s, p − ε ∞, p is the change in the relative permittivity, where ε s, p is the static relative permittivity and ε ∞, p is the relative permittivity at infinite frequency. As can be understood from the equation, ∆ε expresses the strength of coupling between dipole and electric field. The equation describes the time evolution of a dipole accelerated or decelerated by an electric field. In the equation, E(t) expresses the total electric field radiated by the dipole current and also an external field. We solve equations (1) -(3) using the FDTD method described in [6]. Momentum conservation implies that the particle recoils when it radiates an electromagnetic field [8]. This effect is called the radiation reaction, and it is consistently included in Eq. (3). In other words, a two-level system in an excited state makes a transition to a fundamental state by emitting a photon with a certain probability. A transition from a fundamental state to an excited state also occurs when a photon is absorbed. This could occur even in a vacuum through the re-absorption of an emitted photon. Both processes are included in Eq. (3). When we calculate the time evolution of an initially excited dipole located in a vacuum it emits an electromagnetic field. The radiated field accelerates or decelerates the motion of the dipole. Then the amplitude of the dipole current oscillation decays over time. This can be considered as a classical model of the spontaneous emission process. As can be understood from Eq. (3), we assume no artificial damping term of the dipole in our formulation. An amplitude of the radiated field is proportional to the dipole current as ∼ dj(t)/dt. It means the time differentiation of the dipole current term is implicitly included in E(t) of Eq. (3). If we express this term explicitly, we obtain the following equation where the second term expresses friction and E ex (t) on the right hand side expresses the external field. In the equation, 2δ is a proportionality constant. Eq. (4) has the same form as that of Lorentz-dispersive media, with the optical susceptibility given as In the standard FDTD method for Lorentz-dispersive media, the friction term has been thought to correspond to a non-radiative recombination process [9].

Spontaneous emission in a vacuum
Spontaneous emission is a fundamental and purely quantum mechanical phenomenon. Its exact theoretical description can be achieved using a quantum mechanical expression. However, classical theory can also provide approximate results for spontaneous emission. The classical description is based on a classical antenna model. It describes a quantum state transition caused by an oscillating dipole current. The model reasonably describes electromagnetic field radiation by a dipole. The only difference from quantum mechanical theory is that the classical theory gives the half values of the spontaneous emission rate [10,11]. The fact that a quantum mechanical analysis gives twice the value is attributed to the effect of vacuum fluctuation. In classical analysis, the vacuum fluctuation is ignored. Using the expression of optical susceptibility, and Eq. (5), we obtain the following relation, In quantum mechanical theory, the spontaneous emission rate is given as where d is the transition dipole moment. From this relation, the spontaneous emission rate is expressed using FDTD parameters as As a result, we obtain a classical expression of a spontaneous emission rate,  To check the validity of our proposed method, we simulate the radiation process from a dipole located in a vacuum. As an initial condition of the calculation, we set the dipole in an excited state. There is no electromagnetic field at a beginning of the calculation. The dipole radiates an electromagnetic field according to the Eq. (2). Without right side of Eq. (3), a dipole oscillation doesn't decay. The right side of Eq. (3) causes the dipole oscillation to be accelerated or decelerated by the electromagnetic field. A deceleration force make a dipole oscillation decay, which can be understood as recoil force of emission process. If the dipole feels no recoil force, its amplitude does not decay during the simulation. However, the simple exponential decay of the dipole current amplitude occurs, which means the recoil force of emission is included in the simulation. The decay time is estimated by the exponential fitting the time evolution of the dipole current amplitude. The calculated decay time of the dipole current is shown in Fig. 2. Theoretical values obtained with Eq. (10) are also shown in the figure. As seen in the figure, the results agree well. Our proposed numerical model can well describe the spontaneous emission process in a vacuum.

Spontaneous emission in a PhC cavity
An advantage of our model that we wish to emphasize is that it can naturally incorporate a dipoledipole interaction through a radiation field. In the proposed numerical model, dipole oscillation can absorb as well as radiate a field. If there are multiple dipoles, a field radiated by a dipole can be absorbed by another dipole, and vice versa. This means an interaction between dipoles, which in turn means that a dipole-dipole interaction is naturally included in the formulation.

Spontaneous emission of two excited dipoles
Planar PhC cavities are attracting growing interest due to their potential for confining light in extremely small volumes and their high Q factor. The density of state for a photon is greatly modified in a PhC cavity. The probability of transition from an excited state to the fundamental state is dependent on its environment, namely the photon density of states. This is called the Purcell effect and it suggests the possibility of controlling the amount of optical nonlinear effect and also the probability of radiation from excited dipoles artificially.
As an example of a PhC cavity, a PhC cavity with three hole missing (L3) is assumed in the calculation. The thickness of the PhC slab is t = 210 nm, the lattice constant is a = 420 nm, and the radius of the air hole is r = 115.5 nm, and the refractive index of the slab material is assumed to be 3.475. For the fundamental mode, the amplitude of the electric field distribution is shown in Fig. 3 along with the missing holes of the cavity. The electric field has an even symmetry and the magnetic field has an odd symmetry in this direction. The electric field has its maximum value at the center of the cavity and the magnetic field has a node there. We locate one or two dipoles in the cavity at the center of cavity structure and the middle of slab. The two dipoles are located one numerical grid (30 nm) apart. The eigenfrequencies of the dipoles are assumed to be the same. When there are two dipoles, we assume that both are in excited states with the same phase as the initial condition. No electromagnetic field is assumed to exist in the cavity at the beginning of the calculation. Figure 4 shows the dipole current intensity as a function of time for ∆ε = 5.0 × 10 −5 and ω p = 1.19865 × 10 15 Hz. At this condition, dipole's eigen frequency corresponds to resonance condition. The results for one and two dipoles are shown in the figure. As shown in the previous section, the radiated field affects the dipole radiation process. Because of Purcell's effect the cavity mode strongly affects the dipole radiation process in a cavity. The time dependence of the dipole current does not exhibit simple exponential behavior. This is because the formation of a cavity mode by a radiated field needs a certain amount of time. In other words, a radiated field does not formulate a cavity mode instantly and gradually affects the radiation process in a cavity.
With two dipoles, the dipole current decays much faster than that of a single dipole, although simple exponential fitting is not possible. Of course, the decay rate of a dipole strongly depends on the resonant condition of the dipole frequency in the cavity. To check a resonant condition effect on the radiation process from the dipole, the decay rates of two dipoles are calculated for several eigenfrequency conditions. The results for both one and two dipoles are shown in Fig. 5, which reveals clear resonance in the decay rate. This indicates a strong enhancement of the radiative decay when the dipole frequency coincides with the cavity's resonance. For two dipoles, the decay rate in a resonant condition is about twice that of a single dipole. We consider that the enhancement is induced by the interaction between two dipoles. We considered that a  phenomenon similar to superradiance [7] occurs in this case, although the simulation is purely classical. Dicke's superradiance process consists of a two-step process. First, the phases are substantially aligned by the formation of maximally-entangled state through the spontaneous emission, and then the light-emitting body having in-phase acts like a phased array. Here our model deals with the latter of them.  In the previous calculation, two dipoles are located at the center of cavity with one numerical grid (30 nm) apart, where the electric field of the fundamental mode has its maximum value and their initial phases coincide. The radiated field and dipole oscillation always positively interfere in this condition. This interference depends on the relative positions and also on the cavity mode. To check the dipole location dependence on the radiation process, we performed numerical calculations for different dipole conditions as shown in Fig. 3. The results are shown in table 1. (a) are calculations performed for the dipole locations at the center of the PhC cavity, where the fundamental mode of an electric field has it maximum value, (b) correspond to its node, (c) are at its second maximum points, and (d) are around the middle of the electric field amplitude. In each condition, two dipoles are located symmetrically opposite in terms of the cavity center and in the same initial phase. The result indicates that the decay rate of a dipole becomes large at the position where cavity mode has a large amplitude. In contrast, the enhancement of the decay rate for dipoles at the node position is very small.

Effect of multiple dipole interaction in the emission process in a PhC cavity
The cooperative spontaneous emission from a collection of dipoles was studied by Dicke in 1954 [7]. When the dipoles are coherently radiating in phase, the dipoles radiate N times faster than in an incoherent emission, where N is the number of dipoles. This phenomenon is known as superradiance. Superradiance in Dicke's theory is quantum mechanical, but classical theory gives an approximate description of the phenomenon. To check if the proposed method can simulate this phenomenon, we calculate the time evolution of initially excited dipoles in phase.   The dipoles are located in the center of the PhC cavity. The time evolution of the dipole current is calculated and the time dependence of the amplitude decay is shown in Fig. 6. As expected, the calculated results show an enhancement in the radiation rate from the dipole. The amount of enhancement is larger for a larger number of dipoles. Dicke's work indicated destructive interference as well as constructive interference. The spontaneous emission of a dipole can be seriously suppressed by destructive interference. This phenomenon is called subradiance.
To simulate the subradiance that occurs as a result of destructive interference, we simulate the radiation process of two dipoles located in the center of a cavity. As an initial condition, the outof-phase excitation of two dipoles is assumed and compared with two dipoles in-phase excitation and with single dipole excitation. The result is shown in Fig. 7. The figure shows that the decay rate is greatly suppressed with out-of-phase excitation and agrees well with the prediction for subradiance. As a result, we can say that the model can approximately and reasonably describe the dipole-dipole interaction of radiation processes.

Conclusion
In summary, we have presented a numerical technique that allows us to analyze fundamental optical processes, such as the spontaneous emission process of an atom in a vacuum and atoms located in a cavity. The multi-dipole interaction through a radiated electromagnetic field is included in the method. The influence of a dipole-dipole interaction on the emission process in a cavity is simulated and both an enhancement and a suppression of the radiation rate are exhibited. Although the model is based on classical theory, it helps us to better understand the phenomenon. We believe that the numerical method proposed here would be useful for a better understanding of fundamental physical phenomena such as spontaneous emission processes. We thank A. Shinya and H. Gotoh for their encouragement throughout this work.