## Abstract

We investigate the process of light matter interaction in a spherical Mie nanolaser. We derive a rigorous theory based on a three dimensional vector set of Maxwell-Bloch equations and solve the resulting equations through a parallel Finite-Difference Time-Domain Maxwell-Bloch (FDTD-MB) code. Our results predicts a rich physical scenario, ranging from nontrivial vectorial energy matter interplay in the pre-lasing regime to mode competitions and dynamical frequency pulling phenomena. Application of these effects could favor the realization of largely-tunable, nonlinearly controlled nanolaser devices.

© 2008 Optical Society of America

## 1. Introduction

Since the pioneering investigations of Gustav Mie in the early days of the last century [1], the theoretical and experimental study of spherical resonators has deserved much attention from the scientific community. This interest stems from the analysis of fundamental processes such as scattering, energy propagation through disordered media and cavity quantum electrodynamics, and from the large number of applications in photonics, chemistry, meteorology, astronomy and sensing [2, 3, 4, 5, 6, 7, 8, 9]. Mie resonators display ultrahigh quality factors *Q* (≈10^{9}) [10], and are able to confine and store electromagnetic energy for long times within small dielectric volumes; Mie whispering galleries modes (WGM) can sustain enhanced stimulated emission processes so that ‘thresholdless’ lasing can be observed [11, 12, 13, 14]. Despite the wide-spread scientific production, the theoretical understanding of Mie lasers is still an open field of research. For what concerns laser emission, existing theoretical approaches either rely on a semi-classical treatment based on light interaction with a two level system [15, 16, 17] or on rate equations [18], with numerical analysis limited to one and two-dimensional media [19]. The former approach, as was recently observed [20], is not rigorous as the number of physical dimensions grows above one while the latter, accounting only for the contribution of atomic populations, misses important informations on quantum coherence [21], and cannot be effectively pursued to perform any realistic *ab initio* computation at ultra-fast time-scales (<ps). Furthermore, although low dimensional models permit a simplified analysis with respect to the fully vectorial set of Maxwell equations they left the general picture unknown. *Which is the outcome of a strongly nonlinear ultra-fast and multi-dimensional interaction on scales comparable with the wavelength?* Beside being central for the development of laser theory, this question has fundamental importance in quantum mechanics [22], soliton [23] and chaos theory [24]. However, to the best of our knowledge, neither *ab initio* investigations were reported nor quantum-mechanical models of energy-matter interaction for three-dimensional nanostructured lasing systems were developed. Owing to the large interest in the study of nano-scaled particle aggregates [25, 26, 27], the theoretical investigation of such dynamics is of considerable practical importance.

In this Letter we develop an *ab initio* rigorous theoretical model of light interaction in the presence of amplifying (or dissipative) materials, deriving a three dimensional vector set of Maxwell-Bloch (MB) equations within the real representation in terms of the *SU*(*n*) algebra [28]. We discretized the resulting equations on a Yee grid and numerically solved them within the Finite-Difference Time-Domain (FDTD) method [29]. The MB-FDTD approach is then applied to investigate the ultrafast dynamics of Mie nanoresonators. Specifically, we perform a series of numerical experiments by investigating the process of laser emission from a single nanosphere, covered by a layer of active material, for different pumping rates. This theoretical analysis yields the following scenario:

- In the transient before the steady state regime strong interactions between all components of the electromagnetic field and the atom take place. Such a
*vectorial*dynamics is the key ingredient for the reaching of a stable steady state. - When the lasing process settles up, modes possessing frequencies close to Mie resonances get excited and competition dynamics begins.
- As the pumping rate gets higher, a strongly multi-mode regime occurs. In this situation modes with higher quality factors tend to grow faster than the others and dynamical frequency pulling effects are observed. In particular, as the pumping rate increases the principal frequency of laser emission moves toward smaller wavelengths encompassing large shifts; a process radically different from the standard nonlinear frequency shift observed in lasers. The mechanism investigated here can be at the basis of the development of nonlinearly-controlled largely-tunable nanolaser devices.

## 2. Theory

We begin by writing Maxwell equations in an isotropic medium:

$$\frac{\partial \mathbf{E}}{\partial t}=\frac{1}{{\epsilon}_{0}}[\nabla \times \mathbf{H}-\frac{\partial \mathbf{P}}{\partial t}],$$

being **P** the material polarization. The latter can be decomposed as the sum of a linear contribution **P**
_{L}=*ε*
_{0}(*ε _{r}*-1)

**E**plus a nonlinear term

**P**

_{NL}=-

*eN*〈

_{a}**Q̂**〉 modeling the quantum-mechanical interaction between the electromagnetic field and the atoms being

*e*the electric charge,

*N*the density of polarizable atoms and 〈

_{a}**Q̂**〉 the expectation of the displacement operator

**Q̂**with respect to the quantum state |

*ψ*(

**r**,

*t*)〉. To derive a quantum mechanical model of the material, we consider a four-level system, with three degenerate levels, in which three electric mutually orthogonal dipole transitions can be excited by linearly polarized electromagnetic waves with energy equal to the atomic resonance

*ħω*

_{0}(Fig. 1-left). The total Hamiltonian operator

**Ĥ**of the system can be then expressed as

with the unperturbed Hamiltonian **Ĥ**
_{0} of components:

and the self-adjoint interaction Hamiltonian **Ĥ**
_{I}=*e*
**E**·**Q̂**. The displacement operator is

being *q*
_{0} the typical atomic length scale and **x**, **y**, **z** unit vectors along *x*, *y* and *z* axes, respectively. The dynamical evolution of the atomic system is expressed by the Liouville equation of motion of the density matrix operator *$\widehat{\rho}$*=|*ψ*〉〈*ψ*|:

with [*Â*, *B̂*]=*ÂB̂*-*B̂Â* defining the commutator between *Â* and *B̂*. By employing an homomorphism of the quantum mechanical Hilbert space and the Lie algebra *su*(4) [28], we introduce the coherence vector **S** of components *S _{j}*:

being Tr the trace operator and *ŝ _{j}* the

*j*-th generator of the

*SU*(4) algebra [30]. The time evolution of the coherence

**S**is then:

$${\Gamma}_{\mathrm{jh}}=\frac{i}{2\overline{h}}\mathrm{Tr}\left(\hat{\mathbf{H}}[{\hat{s}}_{j},{\hat{s}}_{h}]\right),$$

being *$\widehat{\gamma}$* a phenomenologically added [22] diagonal matrix of nonuniform relaxation rates with *γ _{ii}*=1/

*T*,

_{i}**S**

^{(0)}a vector containing the initial populations of the system. $\widehat{\Gamma}$ is skew and its nonzero components turn out to be:

$${\Gamma}_{\mathrm{1,7}}={\Gamma}_{\mathrm{2,6}}={\Gamma}_{\mathrm{3,5}}=\frac{{\Gamma}_{\mathrm{8,5}}}{\sqrt{3}}={\Gamma}_{\mathrm{14,9}}={\Gamma}_{\mathrm{10,13}}={\Omega}_{y},$$

$$\frac{{\Gamma}_{\mathrm{3,2}}}{2}={\Gamma}_{\mathrm{7,4}}={\Gamma}_{\mathrm{5,6}}={\Gamma}_{\mathrm{12,9}}={\Gamma}_{\mathrm{10,11}}={\Omega}_{x},$$

$${\Gamma}_{\mathrm{1,12}}={\Gamma}_{\mathrm{2,11}}={\Gamma}_{\mathrm{3,10}}={\Gamma}_{\mathrm{4,14}}={\Gamma}_{\mathrm{5,13}}=\frac{2\sqrt{2}{\Gamma}_{\mathrm{15,10}}}{\sqrt{3}}={\Omega}_{z},$$

with
${\Omega}_{i}=\frac{{E}_{i}\wp}{\overline{h}}$
and *℘*=*eq*
_{0}. Finally, the expectation value of the displacement operator **Q̂** expressed in terms of **S** is:

Equations (1)–(9) represent the three dimensional, full-wave vector set of Maxwell-Bloch FDTD equations modeling light-matter interaction in the presence of resonant quantum transitions.

## 3. Implementation

Equations (7)–(9) have been discretized on faces and sides of a cubic Yee grid while (1) on its center. Standard leapfrogging has been used as time marching algorithm for the FDTD, second order Crank-Nicholson for the solution of the MB set. A Uniaxial Perfectly Matched Layer (UPML) is employed to absorb outgoing waves. Following [19], we keep into account fluctuations by adding to the electric field a Markovian noise process with Gaussian statistic, the latter obtained from a Box-Müller random number generator algorithm. The code has been parallelized within the Message-Passing-Interface (MPI) standard.

We perform a series of numerical experiments by considering a spherical resonator (Fig. 1-right) of r=300 nm radius, while taking *ε _{r}*=8.41 as a realistic value for typical high-index materials (e.g.

*TiO*

_{2}), covered by a layer of active material of thickness

*d*=100 nm, with the same refractive index of the interior part and transition at

*ω*

_{0}=3.88·10

^{15}rad/s (λ

_{0}=483 nm) in proximity of Mie resonances, as detailed below. By assuming

*q*

_{0}≈0.1 nm, we obtain a coupling coefficient

*℘*=4.8·10

^{-28}C·m. The decaying constants were set to

*T*=10 fs (

_{i}*i*=1,2,4,5,6,7,9,10,11,12,13,14) and

*T*=1 ns (

_{j}*j*=3,8,15) according to the existing literature [19]. The pumping rate is settled by the density of polarizable atoms

*N*in the excited state, which is assumed to be constant in time and fixed by an external continuous excitation involving additional energy levels, as for standard laser systems. Furthermore, owing the small sphere radius, we assume a spatially uniform pumping of the active layer. The used discretization Δ

_{a}*x*=Δ

*y*=Δ

*z*=10 nm and Δ

*t*=0.006 fs guarantees both high accuracy (with 50 points per wavelength at λ=500 nm) as well as temporal stability.

## 4. Transient regime

We begin our analysis at a low pumping rate (*N _{a}*=10

^{26}m

^{-3}). To investigate the transient regime, we collect at each time step all the electric field components

*E*,

_{x}*E*,

_{y}*E*, sampled in proximity of the sphere center, and the output electromagnetic energy density

_{z}*ℰ*=

**E**·

**D**+

**H**·

**B**(Fig. 2). At

*t*≈600 fs the electromagnetic fields begin to coherently grow. Owing to the sphere-air high index contrast, the electromagnetic energy is well confined within the structure and the lasing mode experiences low losses (Fig. 2). This stage of evolution is characterized by a

*mutual*energy exchange between the atom and each electromagnetic component (Fig. 2a–c). The latter will eventually interact on the fs scale (fixed by

*T*=10fs) with each other owing to the presence of nonlinear non diagonal coupling terms in the density matrix (8). From a physical perspective, energy exchange originates from the presence of a triply degenerate atomic excited state (Fig. 1, left panel), which yields a nontrivial vectorial interplay between the electromagnetic field and the atom. Such a contribution, ignored by previous theories, is the essential ingredient for an isotropic laser three-dimensional dynamics. Without a four-level atomic system, in fact, no vectorial interaction is possible and no realistic steady state can be theoretically predicted.

_{i}## 5. Mode competition dynamics

Once the transient stage has been concluded (*t*≈2 ps), multimode lasing settles up. Even if nano-sized, the sphere has high quality factors that do not prevent lasing to occur (Fig. 3). To study the time evolution of the steady state regime we make use of a spectrogram (Fig. 3) and frequency demultiplexing (Fig. 4). The former shows a *mode-competition* dynamics between two modes in proximity of *ω*
_{0}. In Fig. 3b we plot the extinction factor *Q _{ext}* of the nanosphere as calculated from Mie theory [1], which shows the available modes (Fig. 3b solid line) in the gain bandwidth (Fig. 3b dotted line). The two modes (as shown in Figs. 3a–b) correspond to the nearest WGM modes with respect to the medium resonance

*ω*

_{0}. To characterize the mode-competition dynamics, we extract amplitudes

*a*

_{1}(

*t*),

*a*

_{2}(

*t*) and phases

*ϕ*

_{1}(

*t*),

*ϕ*

_{2}(

*t*) at carrier frequencies close to WGM modes by means of frequency demodulation (Fig. 4). In principle, a four dimensional phase-space [

*a*

_{1},

*a*

_{2},

*ϕ*

_{1},

*ϕ*

_{2}] supports a rich dynamics ranging from periodic motion to chaos. However, at this low pumping rate, no synchronization [24] occurs and each phase independently linearly grows, except in a small set of time instants where

*a*is close to zero (Fig. 4a–d). This is the “free-run” regime, as originally predicted by Lamb [17], where the phase variables can be averaged out from the coupled-mode equations and only the amplitudes remain as dynamical variables. As a result, we are left with a two-dimensional phase space which, according to the Poincaré-Bendixson theorem, can only sustain fixed points or periodic motion, as shown in Fig. 4a and in Fig. 4b, where we plot the map obtained from the Poincaré surface of section of (

_{i}*a*

_{1},

*a*

_{2}). The bi-modal steady state is characterized by a limit cycle where both amplitudes

*a*

_{1},

*a*

_{2}are trapped together.

## 6. Dynamical frequency pulling in multimodal regime

As the pumping rate *N _{a}* gets higher, a large number ofWGM is above the threshold for laser oscillation. To investigate the transition from two-mode to multi-mode regime, we calculate spectrograms obtained from MB-FDTD simulations for increasing pumping rates

*N*∈ [10

_{a}^{26},10

^{28}] m

^{-3}(Fig. 5). The dynamical evolution observed is characterized by two competing effects. As the frequency

*ω*increases with respect to

*ω*

_{0}: i) the gain experienced by the mode gets smaller; ii) the quality factor

*Q*of the WGM mode becomes higher. Therefore, although low wavelength modes experience small gain, they tend to grow faster with respect to modes at higher wavelengths. The result of this competition is a

*dynamical frequency pulling*towards smaller wavelength as the pumping rate grows (Fig. 5c). The predicted frequency shift is very large (≈ 100 nm, see Fig. 5c) and quantitatively it depends on the

*morphology*of the nanosphere resonance spectrum (see Fig. 3b), which fix the involved wavelenghts and the corresponding Q-factors; in principle it can be controlled by acting on the shape of the resonator. Dynamical frequency pulling is then expected to play a fundamental role in the nonlinear dynamics of laser systems, either based on a single sphere or on ordered and disordered resonator ensembles. The induced large shift can be exploited to develop largely-tunable nanosized laser devices operating at different wavelengths by acting on the input excitation; such applications would require a proper engineering of the cavity and will be the subject of future works.

## 7. Conclusion

In conclusion, we have developed a rigorous theory of electromagnetic energy-matter interaction in the presence of nonlinear resonant media. The derived equations are applied to investigate the process of laser emission of Mie nanospheres. By means of parallel simulations and frequency demultiplexing analysis, we discuss both the transient and the steady state regime of laser emission. During the transient, a nontrivial mutual energy exchange between the electromagnetic wave and the resonant medium occurs, thereby justifying a rigorous *ab-initio* investigation. On the other hand, the lasing state observed by increasing the pumping rate is characterized by a dynamical transition from two-mode to multi-mode regime, the latter accompanied by competition phenomena and dynamical frequency pulling. Such results have important implications in the theory of nano-lasers, in both ordered and disordered systems and are expected to stimulate new experiments and applications where single Mie resonators are employed, including nonlinear optics, colloidal physics, chemistry and quantum mechanics.

## Acknowledgements

We acknowledge support from the INFM-CINECA initiative for parallel computing and S. Trillo for useful discussions.

## References and links

**1. **G. Mie, “Beitrge zur Optik trber Medien, speziell kolloidaler Metallsungen,” Ann. Phys. **330**, 377–445 (1908). [CrossRef]

**2. **H. C. van de Hulst, *Light Scattering by Small Particles* (Dover, New York, 1981).

**3. **M. Xu and R. R. Alfano, “Random Walk of Polarized Light in Turbid Media,” Phys. Rev. Lett. **95**, 213901-(4) (2005). [CrossRef]

**4. **Y. Yamamoto, F. Tassone, and H. Cao, *Semiconductor Cavity Quantum Electrodynamics* (Springer, Berlin, 2000).

**5. **H. Yukawa, S. Arnold, and K. Miyano, “Microcavity effect of dyes adsorbed on a levitated droplet,” Phys. Rev. A **60**, 2491–2496 (1999). [CrossRef]

**6. **P. W. Barber and K. Chang, *Optical Effects Associated with Small Particles* (World Scientific, Singapore, 1988).

**7. **M. Kerker, *The Scattering of Light and Other Electromagnetic Radiation* (Academic, New York, 1969).

**8. **P. J. Wyatt, “Scattering of Electromagnetic Plane Waves from Inhomogeneous Spherically Symmetric Objects,” Phys. Rev. **127**, 1837–1843 (1962). [CrossRef]

**9. **A. Ashkin and J. M. Dziedzic, “Observation of Resonances in the Radiation Pressure on Dielectric Spheres,” Phys. Rev. Lett. **38**, 1351–1354 (1977). [CrossRef]

**10. **L. Yang and K. J. Vahala, “Gain functionalization of silica microresonators,” Opt. Lett. **28**, 592-(3) (2003). [CrossRef]

**11. **Y. Yamamoto and R. E. Slusher, “Optical processes in microcavities,” Phys. Today **46**, 66–73 (1993). [CrossRef]

**12. **U. Vietze, O. Krau, F. Laeri, G. Ihlein, F. Schth, B. Limburg, and M. Abraham, “Zeolite-Dye Microlasers,” Phys. Rev. Lett. **81**, 4628-(4) (1998). [CrossRef]

**13. **S. M. Spillane, T. J. Kippenberg, and K. J. Vahala, “Ultralow-threshold Raman laser using a spherical dielectric microcavity,” Nature **415**, 621-(3) (2002). [CrossRef]

**14. **V. Sandoghdar, F. Treussart, J. Hare, V. Lefvre-Seguin, and J. -M. Raimond, “Very low threshold whispering-gallery-mode microsphere laser,” Phys. Rev. A **54**, R1777-(4) (1996). [CrossRef]

**15. **E. T. Jaynes and F. W. Cummings, “Comparison of Quantum and Semiclassical Radiation Theory with Application to the Beam Maser,” Proc. IEEE. **51**, 89–109 (1963). [CrossRef]

**16. **G. W. Gardiner and P. Zoller, *Quantum noise* (Springer, Berlin, 2000).

**17. **W. E. Lamb, “Theory of an Optical Maser,” Phys. Rev. **134**, A1429–A1450 (1964). [CrossRef]

**18. **A. E. Siegman, *Lasers* (University Science Books, Sausalito, 1986).

**19. **G. Slavcheva, J. Arnold, and R. Ziolkowski, “FDTD simulation of the nonlinear gain dynamics in active optical waveguides and semiconductor microcavities,” IEEE J. Sel. Top. Quantum Electron. **10**, 1052–1062 (2004). [CrossRef]

**20. **G. Slavcheva, J. M. Arnold, I. Wallace, and R. W. Ziolkowski, “Coupled Maxwell-pseudospin equations for investigation of self-induced transparency effects in a degenerate three-level quantum system in two dimensions: Finite-difference time-domain study,” Phys. Rev. A **66**, 063418-(21) (2002). [CrossRef]

**21. **B. Bidégaray-Fesquet, F. Castella, and P. Degond, “From Bloch model to the rate equations,” Disc. Cont. Dyn. Sys. **11**, 1–26 (2004). [CrossRef]

**22. **L. Mandel and E. Wolf, *Optical Coherence and Quantum Optics* (Cambridge Press, New York, 1995).

**23. **S. L. McCall and E. L. Hahn, “Self-Induced Transparency,” Phys. Rev. **183**, 457–485 (1969). [CrossRef]

**24. **E. Ott, *Chaos in Dynamical Systems* (Cambridge University Press, Cambridge, 1997).

**25. **K. Li, M. I. Stockman, and D. J. Bergman, “Self-Similar Chain of Metal Nanospheres as an Efficient Nanolens,” Phys. Rev. Lett. **91**, 227402-(4) (2003). [CrossRef]

**26. **C. Conti, L. Angelani, and G. Ruocco, “Light diffusion and localization in three-dimensional nonlinear disordered media,” Phys. Rev. A **75**, 033812-(5) (2007). [CrossRef]

**27. **L. Rojas-Ochas, J. Mendez-Alcaraz, P. S. J.J. Saenz, and F. Scheffold, “Photonic Properties of Strongly Correlated Colloidal Liquids,” Phys.Rev.Lett. **93**, 073903-(4) (2004). [CrossRef]

**28. **F. T. Hioe and J. H. Eberly, “*N*-Level Coherence Vector and Higher Conservation Laws in Quantum Optics and Quantum Mechanics,” Phys. Rev. Lett. **47**, 838-(4) (1981). [CrossRef]

**29. **A. Taflove and S. C. Hagness, *Computational Electrodynamics: The Finite-Difference Time-Domain Method* (Artech House, Boston, 2005).

**30. **Y. Y. Lee and C. T. Chen-Tsay, “The fifteenfold way of the SU(4) symmetry scheme of strongly interacting particles,” Chinese Journal of Physics **3**, 45–68 (1965).