NUMERICAL INVESTIGATION OF COLLECTIVE EFFECTS IN QUANTUM OPTICAL MEDIA VIA INTEGRAL OPERATOR ELECTRIC FIELDS By Elliot Lu A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Physics—Doctor of Philosophy Electrical and Computer Engineering—Dual Major 2023 ABSTRACT The study of quantum optics is principally concerned with investigating light-matter interactions. Within the discipline, computational simulation is a burgeoning field that can lend new insights into optical phenomena previously uncovered by theory or experiment. Collective emission effects such as superradiance serve as one prominent example. In contrast to ordinary emissions, superradiance involves dipolar coupling within optical ensembles and produces a coherent burst of radiation whose intensity scales with the square of the number of emitters. Whereas theoretical results involving superradiance are often shoehorned into small, ideal systems, numerical simulations permit the examination of much larger realistic systems, and can further aid in verifying experimental results. Studies of other phenomena, such as polarization enhancement, inhomogenenous broadening, and subradiance, benefit similarly. To design new systems that exploit quantum optical effects, we devise in this thesis a new numerical approach that can faithfully simulate dynamics of optical active media. Such material are characterized by their ability to modify and re-emit radiation. Nanoscale semiconductor particles known as quantum dots serve as a prime example. Their larger dipole moments–compared to atoms– enable them to experience strong interactions with radiation fields, and permit the observation of a variety of optical phenomena, including superradiance. Despite this merit, numerical simulation of large ensembles of quantum dots–and for long time periods–is challenging. In contrast to previous counterparts, our computational model, which involves the solution to the Maxwell-Bloch equations via integral operator electric fields, is massively scalable in both time and space. This is facilitated by the Adaptive Integral Method (AIM), which effects FFT-based convolutions to evaluate the field. This allows us to perform large scale simulations that reproduce optical effects such as superradiance. To demonstrate the fidelity of our approach, we evaluate the rate of photon emission from our ensemble and show that it reproduces the quadratic scaling of superradiance. In simulations of medium-sized (𝑁 = 50 − 300) ensembles of quantum dots in a Gaussian cloud, we confirm this quadratic scaling by subtracting independent emissions from total emissions. We also observe anisotropy of emission–another hallmark of superradiance–in the field radiated by the Gaussian cloud. Subradiance is revealed in steady state plots of the population excitation, which display diminished emissions. This effect is amplified by inhomogeneous broadening, which induces greater disorder and thus interference within the ensemble, but diminished by the presence of collective Lamb shifts. Additionally, we compare the results of this calculation to those using another formalism, the Master equation. By applying zero-averaging random initial conditions to the polarization, we achieve strong numerical agreement between the two approaches. We observe both superradiant scaling, and destructive interference among dots separated by half-wavelengths. We remark, however, that the Maxwell-Bloch model is superior to the Master equation in resolving time delays and capturing propagation and memory effects. Hence, simulations involving ensembles of emitters separated far apart in space should opt for the Maxwell-Bloch approach to accurately account for delay effects. ACKNOWLEDGEMENTS I would like to thank all members of the guidance committee, C. Piermarocchi, B. Shanker, M. Maghrebi, S. Tessmer, and E. Rothwell, for their time and diligence in reviewing this dissertation. iv TABLE OF CONTENTS LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Active media and secondary radiation . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Collective emission effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 CHAPTER 2 DYNAMICS OF QUANTUM DOTS VIA MAXWELL-BLOCH EQUATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1 Semiclassical Maxwell-Bloch model . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 Computational approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3 Acceleration techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.4 Large scale simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.5 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 CHAPTER 3 COLLECTIVE EMISSION EFFECTS IN ENSEMBLES OF QUANTUM DOTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.1 Dicke model of superradiance . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2 Quantization of electric field . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.3 Maxwell-Bloch equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.4 Photon emission rate, Superradiance . . . . . . . . . . . . . . . . . . . . . . . 31 3.5 Numerical simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.6 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 CHAPTER 4 MASTER EQUATION DESCRIPTION OF EMISSION DYNAMICS . . 45 4.1 Master equation formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.2 Comparison of formulations . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.3 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 CHAPTER 5 CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 APPENDIX A INTEGRAL OPERATOR ELECTRIC FIELDS . . . . . . . . . . . . . 63 APPENDIX B MAXWELL-BLOCH EQUATIONS . . . . . . . . . . . . . . . . . . . 68 APPENDIX C DERIVATION OF LINDBLAD MASTER EQUATION . . . . . . . . . 72 v LIST OF TABLES Table 2.1 Simulation parameters for Section 3.5. Here 𝑐, 𝑒, and 𝑎 0 denote the vacuum speed of light, elementary charge, and Bohr radius respectively. The decoherence times here, while shorter than those typical of optical resonance experiments, afford a shorter computational time but preserve dynamical emission phenomena . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 Table 2.2 AIM parameters for Section 3.5. . . . . . . . . . . . . . . . . . . . . . . . . . . 20 Table 3.1 Simulation parameters for Section 3.5. Here 𝑐 0 and 𝑒 are the vacuum speed of light and elementary charge, respectively. The decoherence (spontaneous emission) time 𝑇1 = 1/Γ0 is obtained as the inverse of the spontaneous emission parameter Γ0 (2.13), hence is proportional to 𝑑 −2 . . . . . . . . . . . . . . . . . . 34 vi LIST OF FIGURES Figure 1.1 Illustration of superradiance: post-excitation, dipoles establish correlation via their radiation fields and emit a superradiant pulse. For this to occur, the separation between emitters must be sufficiently small compared to the wavelength of excitation (not to scale). . . . . . . . . . . . . . . . . . . . . . . 3 Figure 1.2 Illustration of post-excitation emission rates in a hypothetical scenario. Uncoupled emissions display a simple exponential decay. In contrast, coupled emissions occur in a short window and decay transiently. Additionally, whereas uncoupled emissions scale linearly with number of emitters 𝑁, the superradiant contribution to coupled emissions scales as 𝑁 2 (not shown). . . . 4 Figure 2.1 Bloch sphere representation of a two-level quantum system. Equation (2.2) describes the evolution of |𝜓⟩ on the surface of this unit sphere. The “poles" |0⟩ and |1⟩ represent the ground and excited levels, respectively [1]. . . . . . . . 6 Figure 2.2 Precession of a Bloch vector s about a torque 𝛀 lying in the 𝑢𝑤 plane. The presence of field excitations tilts 𝛀 away from the 𝑤-axis, causing the 𝑤-component of s to depart from its initial value and inducing a transition between states. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 Figure 2.3 A set of Lagrange polynomials for 𝑘 = 0, ..., 𝑝 and 𝑝 = 3, serving as time basis functions for the polarization source P̃(r, 𝑡). . . . . . . . . . . . . . . . . 12 Figure 2.4 Illustration of the nearfield criterion for a third order expansion 𝑀 = 3 corresponding to (3 + 1) 2 = 16 gridpoints in this 2D illustration, and 𝛾 = 2. The dashed line indicates the complete nearfield of the box associated with r0 —i.e. all boxes that have an expansion point within 𝛾Δ𝑠 (infinity norm) of the expansion around r0 . Consequently, all of the sℓ (r) within the central dark blue square have a pairwise interaction with the sℓ′ (r) inside the dashed box. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 Figure 2.5 𝑙2 error of the Rabi frequency magnitude | 𝜒| with respect to grid spacing Δ𝑠 for expansion orders 𝑀 = 2 through 6 and 𝛾 = 1, using source and observer boxes of volume 𝜆3 separated by Δ𝑟 = 2𝜆( 𝑥ˆ + 𝑦ˆ + 𝑧ˆ), each containing 64 randomly generated dots. For an expansion order 𝑀 one expects the overall error to scale as 𝑂 (Δ𝑠 𝑀−1 ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 Figure 2.6 FFT runtime (excluding setup time) using a third-order expansion. (Top) 1024 timesteps with Δ𝑠 = 𝜆/400 and prescribed polarizations in the fixed frame. (Bottom) 1000 timesteps with Δ𝑠 = 𝜆/10 and Liouville-dynamics polarization in the rotating frame. Both cases have a quasi-quadratic scaling in the direct calculation, whereas the FFT-accelerated calculation performs slightly worse than linear. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 vii Figure 2.7 ẑ-distribution of polarization | 𝜌˜ 01 | for a 10 000-dot cylindrical simulation, replicating the parameters in [2]. The AIM calculation recovers the oscillatory long-range pattern that we obtained using a direct calculation [2]. . . . . . . . . 19 Figure 2.8 Coloration of | 𝜌˜ 01 | as an indicator of | P̃| at 𝑡1 = 6.0 ps (top), 𝑡 2 = 7.0 ps (middle), 𝑡3 = 8.0 ps (bottom) relative to the peak of a 1 ps-wide pulse, for a system of 100 000 dots. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 Figure 2.9 Scatterplots of | 𝜌˜ 01 | corresponding to the bottom two color plots of fig. 2.8. There exists a single preferred polarization, represented by the linear region of greatest density, arising from dots whose transition dipole moments align or anti-align with the laser field. Radiative coupling produces polarizations that exceed this value. The changes in the red line reflect pulse propagation. . . 22 Figure 3.1 Comparison of photon emission rates for a coupled and uncoupled pair of emitters. Superradiance is described by the enhanced emission rate for the case of coupled emitters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 Figure 3.2 Population excitation 𝜌11 for a Gaussian distribution of 𝑁 dots with 𝑑/𝑑0 = 1.0 and various 𝑁. (a) Per-dot time evolution for 𝑁 = 200. We show here 𝜌11 𝑖 , with 𝑖 corresponding to ten randomly chosen dots. Data for the three dots closest to the center are shown in color. (b) Dot-averaged values ⟨𝜌11 ⟩ for different 𝑁. (c) Immediately post-excitation, a faster decay is observed for larger 𝑁, as evidenced by the crisscrossing of curves. . . . . . . . . . . . . . . 37 Figure 3.3 Total emission rates 𝛾(𝑡) (dashed) and off-diagonal emission rates 𝛾 𝐼 (𝑡) (solid) calculated from (3.43), for the simulations of Figure 3.2. Only the latter, which is associated with superradiant emissions, displays quadratic scaling with number of emitters 𝑁. The superradiant burst lies between 𝑡 ≈ 5 and 15 ps, independently of 𝑁. Overall, the curves bear a strong resemblance to Figure 1.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 Figure 3.4 Time evolution of space-averaged population excitation 𝜌11 for 𝑁 = 200 at various dipole strengths. For larger 𝑑 we observe: (a) a faster decay post- excitation and (b) a reduction in the number of subradiant, slowly-decaying states, as suggested by smaller populations at 1 ns. . . . . . . . . . . . . . . . . 38 Figure 3.5 Log-plot of averaged-over-trials population ⟨𝜌11 ⟩ at 1 ns as a function of dipole strength and 𝑁. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 viii Figure 3.8 Plot of logarithmic radiated field intensities for 𝑁 = 100, 𝑑/𝑑0 = 1.0 (a, b), 𝑁 = 200, 𝑑/𝑑0 = 1.0 (c), 𝑑/𝑑0 = 2.0 (d) in time and space. Evident are not only spatial but temporal oscillations, becoming enhanced for larger values of the dipole moment. It is also evident that the intensity amplitude increases with the number of emitters. Groups of dots in the cloud (𝑦 ∼ 0) undergo emission synchronization leading to periodic oscillations that become more irregular as the density increase. Finally, (b) displays emission enhancement in the laser propagation direction (cf. Fig 3). . . . . . . . . . . . . . . . . . . . 39 Figure 3.6 Colormaps of logarithmic field intensity (field norm squared) for a configuration with 𝑁 = 100 and 𝑑/𝑑0 = 1.00 after 20 ps, on the 𝑥 − 𝑦 (top) and 𝑥 − 𝑧 (bottom) planes. The spatial oscillations occur with a period about half the wavelength of 253 nm. Also note the enhancement along the laser propagation direction in the positive 𝑧 axis. . . . . . . . . . . . . . . . . . 40 Figure 3.7 Plot of logarithmic field intensity along the 𝑦 axis for 𝑑/𝑑0 = 1.0 (Blue) and 𝑑/𝑑0 = 2.0 (Magenta), corresponding to Fig. 3 (a). The Fourier data is normalized against the 𝑑/𝑑0 = 2.0 series. . . . . . . . . . . . . . . . . . . . . 41 Figure 3.9 Temporal (top) and spatial (bottom) Fourier plots of field intensity 𝐼 (𝑦, 𝑡) = |E(𝑦, 𝑡)| 2 along the 𝑦 axis (as in Fig. 4) for 𝑁 = 200 and various dipole strengths, normalized by maximum intensity: 𝐼˜1 ( 𝑓 ) = 𝐼˜(𝑘 = 0, 𝑓 )/ 𝐼˜(0, 0) (top), 𝐼˜2 (𝑘) = 𝐼˜(𝑘, 𝑓 = 0)/ 𝐼˜(0, 0) (bottom). Spectral broadening with increasing dipole strength is evident in both cases. . . . . . . . . . . . . . . . . 42 Figure 3.10 Averaged population excitation curves for 𝑁 = 200 including and excluding the collective Lamb shift. In the latter (Magenta dashed), the occurrence of purely subradiant states is suggested by the complete suppression of emission following the period of superradiant decay. For reference, a plot for uncoupled emitters is also shown, exhibiting only spontaneous decay. . . . . . . . . . . . . 43 Figure 3.11 Off-diagonal emission rates 𝛾 𝐼 (𝑡) calculated from (3.43) for 𝑁 = 200, including and excluding the collective Lamb shift. We observe that superradiant emissions (𝑡 ≈ 5 and 15 ps) are diminished by the shift. . . . . . . 43 Figure 3.12 Time evolution of space-averaged population excitation ⟨𝜌11 ⟩ for 𝑁 = 200, for different inhomogeneous broadening values 𝛿 (in meV). (Inset) Excitation values as a function of 𝛿, at 𝑡 = 500 ps (blue circles) and 𝑡 = 1000 ps (red squares). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 Figure 4.1 Comparison of the off-diagonal term of the emission rate, 𝛾 𝐼 (𝑡) using the semiclassical Maxwell Bloch (solid) and the Master equation (dashed). Eight √ dots are initially prepared in the maximum polarization state (|0⟩ + quantum |1⟩)/ 2. Curves for different dot separations 𝑠 compared to 𝜆 are shown. . . . . 49 ix Figure 4.2 Comparison of 𝛾 𝐼 (𝑡) between the semiclassical Maxwell Bloch and the Master equation. Eight quantum dots equally spaced by 0.1𝜆 are considered. Each dot is initially in the maximum population state with the addition of a random polarization (the average has been taken over 500 random initial conditions). Additionally, the dots are divided into 𝑁𝑔 groups each receiving a different random polarization phase. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 Figure 4.3 The second term of the emission rate, 𝛾 𝐼 (𝑡), via the semiclassical Maxwell Bloch equations is shown, for the Gaussian distributed dot configurations of Sect. III and all dots starting in the excited state (no pulse, 𝑁𝑔 = 1). Solid curves are calculated directly from (4.12), while dashed curves are calculated by subtracting independent emissions from the total emission as in (3.43). The average has been taken over 500 random initial conditions. . . . . . . . . . 51 Figure 4.4 Comparison of short time behavior of the population excitation ⟨𝜌01 ⟩ using the Master equation and Maxwell-Bloch equations with and without delays. 𝑁 = 2 (top) and 𝑁 = 8 (bottom) quantum dots are separated by 1.0𝜆, with the first dot initialized in the excited state with a random initial polarization, and the rest lying initially in the ground state. Shown is the population inversion ⟨𝜌01 ⟩ for the last dot in the chain, i.e. the one that would receive a retarded signal the latest. The delay in excitation corresponds to the time 𝑅/𝑐 for the signal to propagate to the last dot (for a total separation 𝑅 = 1.0𝜆, 𝑅/𝑐 ≈ 0.157 𝑇1 ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 x CHAPTER 1 INTRODUCTION 1.1 Active media and secondary radiation In contrast to conventional matter which scatters light directly, the term active media refers to matter which modify propagating waves by absorbing and re-emitting radiation. Such media consist of nanoscale elements, typically atoms and semiconductor quantum dots, characterized by strong optical resonances. These materials see applications as diverse as gain media for lasers, optical amplifiers, and components for quantum technologies. In this respect, quantum dots possess salient advantages over their atomic counterparts. Their nanoscale size tightly confines electrons and induces energy quantization. By further tuning their size, shape, and composition, dots can be engineered to emit light at a wide range of desired frequencies. For instance, smaller sizes yield larger differences in energy between the valence and conduction bands, which increases the frequency of secondary radiation. Additionally, dots experience stronger dipolar transitions than atoms; this permits more pronounced coupling via secondary radiated (re-emitted) fields and enables the observation of nonlinear effects, such as Rabi oscillations [3, 4, 5], at lower laser intensities. In light of these merits, we focus our analysis of active media in this thesis on quantum dots. The consequences of secondary radiation in quantum dots have been studied recently in both experimental [6] and theoretical/computational [7, 8] works. In the latter, Maxwell-Bloch equations [9] are typically used to describe interactions among an ensemble of dots, each modeled as a many- level quantum system. Such a description may be semiclassical, wherein the dots absorb and re-emit classical Maxwellian fields, or the fields may be fully quantized. Regardless, the innate coupling of quantum mechanics to electrodynamics results in a set of equations which are highly formidable to solve. Theoretical calculations are thus limited to systems with very small numbers (< 10) of emitters, and neglect effects due to the quantization of the electromagnetic fields. Numerical solutions, on the other hand, have witnessed an evolution from continuum models [10, 11] to models employing spatial homogeneity [12] and mesh-based PDE solvers [13, 14]. These methods share 1 common shortcomings: they are neither able to resolve short-ranged effects, nor are practically scalable to larger systems. In Chapter 2, we describe a computation approach to the semiclassical Maxwell-Bloch equations that remedies these shortcomings, enabling the simulation of a variety of radiation-induced phenomena with ensembles of more than 104 dots. 1.2 Collective emission effects A major consequence of secondary radiation is the existence of collective emission, wherein an ensemble of active media cooperates to emit radiation coherently. In contrast to light from independent emitters, superradiance is defined by a greater intensity, as well as anisotropy [9, 10] of radiation. The dual phenomena, in which destructive coupling diminishes the intensity of radiation, is termed subradiance. Following Dicke’s seminal work [15] in the 1950s, superradiance has been examined extensively in the theoretical literature [16, 17, 18, 19, 20]. In theory, the occurrence of superradiance is deeply connected to the nonlinear coupling of quantum states with their radiated fields. It involves the coherent “alignment" of dipoles, producing a short radiation burst whose intensity scales as the square of the number of emitters. Thus, the phenomenon can be regarded as a process in which a disordered system builds up order via the correlation amongst dipoles. It can also be viewed as a process in which an initially quantum mechanical state takes on a classical character, i.e. in which a set of correlated dipoles emits radiation akin to that of a classical antenna [9]. In this regard, the representation of the radiation as a classical electric field is a fitting one, and sufficient to capture the intricacies of collective emission for many applications. 2 Figure 1.1 Illustration of superradiance: post-excitation, dipoles establish correlation via their radiation fields and emit a superradiant pulse. For this to occur, the separation between emitters must be sufficiently small compared to the wavelength of excitation (not to scale). Experimental observation of superradiance became widespread after the invention of dye-laser systems, which produce short, narrow bursts that enable excitation of a collection of media to specific energy levels [9]. Early studies in the 1970s revealed superradiance in atoms [21, 22] and gas molecules [23]; more recently, it has been observed in lead halide perovskite lattices [24], atomic lattices [25], free-electron lasers [26], and isolated quantom dots [27]. In a typical study, a pulsed laser that is resonant with the transition frequency of the material is used to initially excite the system. The resulting secondary radiation is emitted with characteristic intensity curves (Fig. 1.2). Whereas independent emissions decay exponentially, superradiant emissions occur in a transient time window. Thereafter, one typically observes a diminished emission rate. This characterizes subradiance, in which dipoles anti-align to cause destructive interference and the suppression of emission. In Chapter 3, we shall numerically explore both collective effects in large ensembles of quantum dots. 3 5 Normalized emission rate 4 3 2 1 0 Time Figure 1.2 Illustration of post-excitation emission rates in a hypothetical scenario. Uncoupled emissions display a simple exponential decay. In contrast, coupled emissions occur in a short window and decay transiently. Additionally, whereas uncoupled emissions scale linearly with number of emitters 𝑁, the superradiant contribution to coupled emissions scales as 𝑁 2 (not shown). Despite the breadth of theoretical and experimental literature, numerical simulation of collective emission in large systems is a relatively unexplored area. To induce transitions, simulations typically mimic experiments and incorporate a resonant pulsed laser. Modeling the subsequent re-radiation is a much more arduous task. In general, this can be approached from the semiclassical Maxwell- Bloch or quantum Master equation models. Each formulation poses distinct challenges while offering certain advantages. The former involves the solution of a system of 𝑁 coupled, nonlinear, delayed differential equations, each describing the evolution of a simple 2 × 2 density matrix representing the state of each quantum dot (𝑁 is the total number of dots in the system). As such, it is readily scalable to large systems of dots. In contrast, the Master equation is a linear differential equation for the enormous 2𝑁 ×2𝑁 density matrix of the entire ensemble, which impedes large-scale simulations. Moreover, the Born-Markov approximation implies that delays are absent from the Master equation, a feature that simplifies calculations but makes this formulation ill-suited for modeling propagation and memory effects in cases when the electromagnetic field is far from equilibrium. We elaborate on these remarks in Chapters 3 and 4, which are devoted to simulations involving the Maxwell-Bloch equations and Master equation, respectively. 4 CHAPTER 2 DYNAMICS OF QUANTUM DOTS VIA MAXWELL-BLOCH EQUATIONS A theory involving the coupling of quantum mechanics to electrodynamics has several realizations. Fully quantized theories of the electromagnetic field are able to capture quantum statistics, and faithfully describe such phenomenon as sub-Poissonian photon emissions. However, in analyses of dense ensembles of discrete emitters such as quantum dots, which involve measurements of macroscopic quantities such as field intensity or even spontaneous emission [28], a classical treatment of the field is often sufficient. The interplay between a classical field and an ensemble of dots represented by many-level quantum states, defines the semiclassical model and is the crux of this chapter. The semiclassical model permits us to draw upon techniques from classical electrodynamics to analyze the electric field. It represents the field in terms of propagating waves emanating from a source distribution. Mathematically, we write the field as the action of a wave equation Green’s function on the source. This integral operator approach–so called because of the integral kernel that appears in the Maxwell-Bloch equations–has several merits. It readily admits a numerical discretization in terms of spatial and temporal basis functions. The former encodes the spatial distribution of dots, while capturing their mesoscopic local structure; the latter resolves time derivatives and captures retardation. Our method also delineates near- and farfield regions of interaction, and facilitates the tracking of data in both regimes. Compared to predecessors such as continuum models [10, 11] and mesh-based PDE solvers [13, 14], nearfield interactions can be resolved at a practical computational cost. As such, our computational algorithm [29] enables the simulation of coupling effects due to large numbers (𝑁 > 100) of dots separated by distances far smaller than the excitation wavelength. Despite the superior scalability of the integral operator approach, a direct implementation is still hampered by the 𝑂 (𝑁 2 ) cost of computing matrix-vector products in convolutions. In Section 2.3, we describe an acceleration technique for reducing the computation cost to 𝑂 (𝑁 log 𝑁). In this Adaptive Integral Method (AIM), source data is projected onto a set of auxiliary basis functions 5 representing grid points, through which interactions are evaluated and then projected back onto observers. We demonstrate both the expected 𝑂 (𝑁 log 𝑁) performance scaling, and convergence of the algorithm based on error scaling with grid density and expansion order. By applying AIM to the integral-operator Maxwell-Bloch equations, we run large scale (𝑁 > 10000) simulations displaying effects due to secondary radiation, including oscillatory patterns, polarization modulations, and ultimately super- and subradiant emissions. 2.1 Semiclassical Maxwell-Bloch model In what follows, we outline a self-consistent model of quantum dot ensembles, and later provide insight into the physics. In the context of this thesis, we restrict ourselves to two-level systems (a generalization to three-level systems will briefly be discussed in the Conclusion) with associated transition frequency 𝜔0 . In a typical process, stimulation by a resonant light pulse induces electric dipole transitions from the ground to excited state, and generates a coherent polarization, which, in turn, induces re-emission of secondary radiation. The two-level system admits a geometric interpretation in the Bloch vector s = (𝑢, 𝑣, 𝑤)𝑇 , with the configuration space of Bloch vectors known as the Bloch sphere (Fig. 2.1). ŵ = |1⟩ |𝜓⟩ 𝜃 v̂ 𝜑 û −ŵ = |0⟩ Figure 2.1 Bloch sphere representation of a two-level quantum system. Equation (2.2) describes the evolution of |𝜓⟩ on the surface of this unit sphere. The “poles" |0⟩ and |1⟩ represent the ground and excited levels, respectively [1]. 6 An equivalent formulation in terms of a 2 × 2 density matrix 𝜌ˆ is given by the one-to-one correspondence: © 𝜌00 𝜌01 ª 1 © 1 − 𝑤 𝑢 + 𝑖𝑣 ª 𝜌 = ­­ ®= ­ ® 2­ ® ® (2.1) 𝜌10 𝜌11 𝑢 − 𝑖𝑣 1 + 𝑤 « ¬ « ¬ In the 𝜌 basis, the Liouville equation describes the time evolution of each dot: d𝜌 −𝑖 = [H (𝑡), 𝜌] − D [𝜌]. (2.2) d𝑡 ℏ where the dot Hamiltonian © 0 −ℏ𝜒(𝑡) ª H (𝑡) ≡ ­­ ® (2.3) ∗ ® −ℏ𝜒 (𝑡) ℏ𝜔0 « ¬ governs the interaction with an electric field E(r, 𝑡). It consists of diagonal terms representing the internal energies ℏ𝜔0 of each state, and off-diagonal terms containing the interaction with fields via 𝜒(𝑡) = d · E(r, 𝑡)/ℏ. Here d is the dipole moment of the transition, which determines the strength and direction of the induced polarization. Additionally, one typically uses a decoherence matrix © (𝜌00 − 1)/𝑇1 𝜌01 /𝑇2 ª D [𝜌] ≡ ­­ ® ® (2.4) 𝜌10 /𝑇2 𝜌11 /𝑇1 « ¬ to describe the effects of spontaneous emission on each dot. Here 𝑇1 and 𝑇2 are time constants, related to the spontaneous decay rate Γ0 via 𝑇1 = 1/Γ0 , 𝑇2 = 2𝑇1 = 2/Γ0 . Equation (2.2) has a physical interpretation in terms of Bloch vectors. Using eq. (2.1) and assuming real fields, rewriting eq. (2.2) in terms of Bloch components yields:   𝑢¤ = −𝜔0 𝑣 − 𝑢/𝑇2         𝑣¤ = 𝜔0 𝑢 − 2𝜒𝑤 − 𝑣/𝑇2 (2.5)       𝑤¤ = 2𝜒𝑣 − (𝑤 + 1)/𝑇1    Disregarding the damping terms, this can be written as a precession equation: s¤ = 𝛀 × s (2.6) 7 where the effective torque 𝛀 = (2𝜒, 0, 𝜔0 )𝑇 acts as the axis of rotation of the Bloch vector s. In the absence of electric fields (𝜒 = 0), 𝛀 is parallel to the 𝑤-axis, and the system rotates in a plane of constant 𝑤; thus there are no transitions between states. The presence of fields tilts 𝛀 away from the 𝑤-axis, stimulating transitions (see Fig. 2.2). These oscillate between the two states at a rate equal to the Rabi frequency 𝜒, hence the term Rabi oscillations. ŵ 𝛀 s v̂ û Figure 2.2 Precession of a Bloch vector s about a torque 𝛀 lying in the 𝑢𝑤 plane. The presence of field excitations tilts 𝛀 away from the 𝑤-axis, causing the 𝑤-component of s to depart from its initial value and inducing a transition between states. To stimulate transitions, E(r, 𝑡) typically includes an incident laser field E 𝐿 (r, 𝑡) oscillating at frequency 𝜔 𝐿 . In all relevant scenarios, this laser is nearly resonant with the transition frequency: 𝜔0 ∼ 𝜔 𝐿 , which lies in the optical frequency band (∼ 2278 ps−1 ), for the type of dots considered in this thesis). Resolving such fast oscillations requires the choice of a very small timestep; hence integrating eq. (2.2) directly on the timescale of typical Rabi dynamics (∼ 1 ps) becomes computationally laborious. We therefore introduce 𝜌˜ = 𝑈 𝜌𝑈 † where 𝑈 = diag(1, 𝑒𝑖𝜔 𝐿 𝑡 ) and 𝑈 † denotes the conjugate transpose of 𝑈, transforming to a frame rotating with the frequency of the incident laser. In this rotating frame, (2.2) becomes: d 𝜌˜ −𝑖   = H̃ (𝑡), 𝜌˜ − D [ 𝜌]. ˜ (2.7) d𝑡 ℏ where the rotating frame Hamiltonian: © 0 −ℏ𝜒(𝑡)𝑒 −𝑖𝜔 𝐿 𝑡 ª ˜ 𝐻 (𝑡) = ­ ­ ® (2.8) ∗ ® −ℏ𝜒 (𝑡)𝑒 𝐿 ℏ(𝜔0 − 𝜔 𝐿 ) 𝑖𝜔 𝑡 « ¬ 8 contains only terms proportional to 𝑒𝑖(𝜔0 ±𝜔 𝐿 )𝑡 if E 𝐿 (𝑡) ∼ Ẽ 𝐿 (𝑡) cos(𝜔 𝐿 𝑡). The rotating wave approximation (RWA) neglects the high frequencies 𝜔0 + 𝜔 𝐿 , assuming that terms containing them will integrate to zero in solving (2.7) over appreciable timescales. In the presence of interactions among dots, their secondary radiated fields E𝑟𝑎𝑑 are added to the laser field in calculating the total electric field experienced by each dot: E(r, 𝑡) ≡ E 𝐿 (r, 𝑡) + E𝑟𝑎𝑑 (r, 𝑡) (2.9) Here we assume the radiation field E𝑟𝑎𝑑 (r, 𝑡) = 𝔉{P(r, 𝑡)} arises from a polarization density Í  𝑖 P(r, 𝑡) = 𝑖 𝛿(r − r𝑖 ) d Re 2𝜌01 (𝑡)) , where 𝑖 labels each dot located at r𝑖 . The radiated field can then be written in integral operator form as (see Appendix A):   𝔉{P(r, 𝑡)} ≡ −𝜇0 𝜕𝑡2 − 𝑐2 ∇∇ 𝑔(r, 𝑡) ★𝑠𝑡 P(r, 𝑡) ∫ " # 𝜕𝑡2 P(r′, 𝑡 𝑅 ) 𝜕𝑡 P(r′, 𝑡 𝑅 ) P(r′, 𝑡 𝑅 )  −1 = (I − r̂ ⊗ r̂) · 2 + (I − 3r̂ ⊗ r̂) · 2 + 3 𝑑 3 r′ 4𝜋𝜖 𝑐 𝑅 𝑐𝑅 𝑅 (2.10) Here, 𝑔(r, 𝑡) = 𝛿(𝑡 𝑅 )/𝑅 is the retarded wave equation Green’s function, R = r − r′, r̂ = R/𝑅, ⊗ is a tensor product, and 𝑡 𝑅 = 𝑡 − 𝑅/𝑐. As the polarization density P(r, 𝑡) arises from the off-diagonal elements of 𝜌, ˆ this formulation effectively couples the evolution of all quantum dots via their radiated fields 𝔉{P(r, 𝑡)}. Furthermore, it depicts the propagation of electric signals through space with finite velocity, such that each dot receives the field radiated by another dot at the retarded time 𝑡𝑅. In the rotating frame, the source distribution must also be transformed. Writing P(r, 𝑡) = P̃(r, 𝑡)𝑒𝑖𝜔 𝐿 𝑡 and substituting into (2.10) gives the rotating frame equivalent of the radiated field: ∫ " −1 𝜕𝑡2 P̃(r′, 𝑡 𝑅 ) + 2𝑖𝜔 𝐿 𝜕𝑡 P̃(r′, 𝑡 𝑅 ) − 𝜔2𝐿 P̃(r′, 𝑡 𝑅 ) 𝔉{P̃(r, 𝑡)} ≡ (I − r̂ ⊗ r̂) · 4𝜋𝜖 𝑐2 𝑅 # (2.11) 𝜕𝑡 P̃(r′, 𝑡 𝑅 ) + 𝑖𝜔 𝐿 P̃(r′, 𝑡 𝑅 ) P̃(r′, 𝑡 𝑅 )  −𝑖𝜔 𝐿 𝑅/𝑐 3 ′ + (I − 3r̂ ⊗ r̂) · + 𝑒 𝑑 r 𝑐𝑅 2 𝑅3 In most scenarios, the incident laser E 𝐿 is sharply peaked and decays transiently in a time 2𝑡 0 . Thereafter, the polarization density P̃(𝑡) varies appreciably slower than the laser. In that case, we 9 may disregard the derivative terms and write for the Rabi energy of a dot at r and for 𝑡 ≫ 2𝑡0 : 𝜒(r, 𝑡) ≡ d · 𝔉{P̃(r, 𝑡)}/ℏ −P̃ˆ (r′, 𝑡 𝑅 ) 𝑖 P̃ˆ (r′, 𝑡 𝑅 ) P̃ˆ (r′, 𝑡 𝑅 ) ∫ " !# 3Γ0 ≈− (I − r̂ ⊗ r̂) · + (I − 3r̂ ⊗ r̂) · + 𝑒 −𝑖𝜔 𝐿 𝑅/𝑐 𝑑 3 r′ 4 𝑘𝑅 (𝑘 𝑅) 2 (𝑘 𝑅) 3 (2.12) where 𝑘 = 𝜔 𝐿 /𝑐, P̃ = P̃ˆ 𝑑 and: 𝑑 2 𝜔3𝐿 Γ0 ≡ (2.13) 3𝜋𝜖ℏ𝑐3 is the spontaneous decay rate (we have further assumed that all dipole strengths are equal). 2.1.1 Radiation reaction field Expressions (2.10-2.12) appear singular in the limit 𝑅 → 0. However, it is possible to approximate this limit such that the resulting quantity is finite. This leads to the radiation reaction field, a self-field which an emitter experiences due to its own emitted radiation. We may derive such a field E 𝑅𝑅 from (2.10) via Taylor expansion (see Appendix A for details). The result in the rotating frame is [30, 31]: ∫ 1 2  3 ′ 2 ′ 2 ′ 3 ′  ′ E 𝑅𝑅 (r, 𝑡) ≡ − 𝜕𝑡 P̃(r , 𝑡) + 3𝑖𝜔 𝐿 𝜕𝑡 P̃(r , 𝑡) − 3𝜔 𝜕𝑡 P̃(r , 𝑡) − 𝑖𝜔 P̃(r , 𝑡) 𝑑r (2.14) 4𝜋𝜖 3𝑐3 𝐿 𝐿 As before, the last term in the above expression will dominate after the laser decays in 2𝑡 0 . Substitution this term into (2.7) via (2.8) and (2.9) yields the following term to the RHS: ©− (1 − 2𝜌00 ) 2 − 1 2𝜌01 (1 − 2𝜌00 ) ª  𝛽 𝜌¤̃ 𝑅𝑅 ≈ ­ ® (2.15) 2 ­ 2𝜌 (1 − 2𝜌 ) (1 − 2𝜌00 ) − 1 2 ® 10 00 « ¬ Here 𝛽 = (1/4𝜋𝜖)(2𝑑 2 𝜔3𝐿 /3ℏ𝑐3 ) is a parameter that measures the strength of radiation reaction [32]. As suggested by the simple relation 𝛽 = Γ0 /2, the radiation reaction field functions similarly to (2.4) and serves to dampen the system, albeit in a nonlinear manner and with an overall magnitude half that of the spontaneous decay. The self-field E 𝑅𝑅 further induces a “residual" nearfield around the vicinity of each dot, i.e. we may write for the total nearfield: E𝑁 𝐹 (r, 𝑡) = E 𝑅𝑅 (r, 𝑡) + E𝑁 𝐹−𝑅𝑅 (r, 𝑡) (2.16) 10 where the support of E𝑁 𝐹−𝑅𝑅 extends beyond the location of the point source (see Appendix A). Thus, the inclusion of the radiation reaction field constitutes a model of superradiance where this field strongly couples dots separated by distances much smaller than the excitation wavelength. This idea has been explored in several works such as [12, 32]. However our simulations have not definitively shown that including E 𝑅𝑅 at the location of each dot necessarily reproduces superradiance. Nonetheless, it is illustrative of an interesting physical effect that can be incorporated into our integral operator model. 2.2 Computational approach We now consider the numerical solution of (2.7). To start, let us expand P̃(r, 𝑡) in terms of 𝑁 𝑠 space and 𝑁𝑡 time basis functions such that: 𝑠 −1 𝑁 𝑁∑︁ ∑︁𝑡 −1 P̃(r, 𝑡) ≈ Ãℓ(𝑚) sℓ (r)𝑇 (𝑡 − 𝑚 Δ𝑡), (2.17) ℓ=0 𝑚=0 Assuming that the dots are point sources, we take s𝑙 (r) = 𝛿(r − r𝑙 ). Here 𝛿(r) is a 3-D delta function, d𝑙 and r𝑙 denote respectively the dipole moment and position of the ℓ th dot, while Ãℓ(𝑚) = d𝑙 𝜌˜ ℓ,01 (𝑚 Δ𝑡)) represents its polarization at the 𝑚 th timestep. The time basis functions 𝑇 (𝑡) interpolate the function of interest at each timestep, and are required to have finite support and obey the causality clause 𝑇 (𝑡) = 0 if 𝑡 < −Δ𝑡. Additionally they must be at least twice-differentiable to recover the time derivatives in (2.10). To this end we elect to use shifted Lagrange polynomials (𝜏 = 𝑡/Δ𝑡) [33]: 𝑘 𝑝−𝑘 𝜏 − ( 𝑝 − 𝑘) − 𝑖 Ö 𝜏 − ( 𝑝 − 𝑘) + 𝑖   Ö × 𝑘 −1 ≤ 𝜏 < 𝑘    , 𝑝   𝐿 𝑘 (𝑡) = 𝑖=1 −𝑖 𝑖=1 +𝑖 (2.18)    0, otherwise    𝑝 of order 𝑝 = 3 or higher (see Fig. 2.3), in terms of which 𝑇 (𝑡 − 𝑚Δ𝑡) = 𝐿 𝑘 (𝑡) for 𝑚 = 𝑝 − 𝑘. 11 1 0.5 0 0 0.5 1 1.5 2 2.5 3 𝑡/Δ𝑡 Figure 2.3 A set of Lagrange polynomials for 𝑘 = 0, ..., 𝑝 and 𝑝 = 3, serving as time basis functions for the polarization source P̃(r, 𝑡). Substituting (2.17) into (2.9) via (2.10) and projecting the result onto the 𝛿(𝑡 − 𝑚 Δ𝑡)sℓ (r) basis yields a set of discrete convolution equations: 𝑚 ∑︁ ′ ′ Ẽ (𝑚) = ẼL(𝑚) + F̃ (𝑚−𝑚 ) · à (𝑚 ) (2.19) 𝑚 ′ =0 where Ẽℓ(𝑚) ≡ sℓ (r), Ẽ(r, 𝑚 Δ𝑡) ; 0 ⩽ ℓ < 𝑁𝑠 (2.20a) (𝑚) ẼL,ℓ ≡ sℓ (r), ẼL (r, 𝑚 Δ𝑡) ; 0 ⩽ ℓ < 𝑁𝑠 (2.20b) and F̃ (𝑘) gives a sparse matrix of dimension 𝑁 𝑠 × 𝑁 𝑠 such that F̃ℓℓ(𝑘) ˜ ′ ≡ sℓ (r), 𝔉{sℓ ′ (r)𝑇 (𝑘 Δ𝑡)} . (2.20c) Note that due to the finite support of the 3-D retarded potential, F̃ (𝑘) has a sparse, Toeplitz lower-triangular structure. This facilitates a cost complexity of O (𝑁𝑡 ), as only a fixed number of multiplications need to be performed at each timestep 𝑚Δ𝑡 [34]. Equation (2.19) defines a marching-on-time scheme for evaluating the total electric field. The determination of the polarization à (𝑚+1) thus proceeds from integrating the equation of motion (2.7) from 𝑡𝑖 = 𝑚Δ𝑡 to 𝑡 𝑓 = (𝑚 + 1)Δ𝑡 for every quantum dot. To solve this system, we use a numerical predictor-corrector approach derived in [35]. Defining 𝑡 𝑚 ≡ 𝑚Δ𝑡 and approximating 𝜌(𝑡) ˜ as a weighted sum of exponentials, the predictor-corrector scheme proceeds with an extrapolation 12 predictor step: 𝑊−1 ∑︁ 𝜌˜ ℓ (𝑡 𝑚+1 ) ← P𝑤(0) 𝜌˜ ℓ (𝑡 𝑚−𝑤 ) + P𝑤(1) 𝜕𝑡 𝜌˜ ℓ (𝑡 𝑚−𝑤 ), (2.21) 𝑤=0 and iterated corrector steps: 𝑊−1 ∑︁ 𝜌˜ ℓ (𝑡 𝑚+1 ) ← C𝑤(0) 𝜌˜ ℓ (𝑡 𝑚−𝑤 ) + C𝑤(1) 𝜕𝑡 𝜌˜ ℓ (𝑡 𝑚−𝑤 ) (2.22) 𝑤=−1 Such an integrator has significantly better convergence properties than Runge-Kutta integrators for equations of the type seen in (2.7) and naturally accommodates basis functions within 𝑐 Δ𝑡 of each other. A self-consistent solution to eq. (2.7) has the following prescription for any timestep: (i) determine Ãℓ(𝑚) = 2 Re( 𝜌˜ ℓ,01 (𝑚 Δ𝑡)) from the known history of the system, (ii) compute Ẽℓ(𝑚) using eq. (2.19), (iii) find 𝜕𝑡 𝜌˜ ℓ,01 (𝑚 Δ𝑡) using eq. (2.7), and (iv) correct 𝜌˜ ℓ,01 (𝑚 Δ𝑡) and iterate steps (ii) through (iv) until converged [34]. 2.3 Acceleration techniques It is apparent that the discrete convolution in (2.19) involves a matrix-vector product with spatial complexity 𝑂 (𝑁 2 ). This step bottlenecks the naive algorithm, which has limited scalability in terms of number of emitters. Techniques for effecting a sub-quadratic complexity are varied, but general fall into two regimes. In kernel-type solvers, nearby sources are aggregated and propagated far away via an appropriate series representation of the Green’s function. In contrast, source-type solvers transform the source distribution independently of the kernel. As they are adaptable to a much wider range of kernels, such as the rotating frame kernels encountered in Section 2.2, source-type solvers will be employed to accelerate the simulation of dot dynamics. The Time Domain Adaptive Integral Method (TD-AIM, or simply AIM) is a source-type solver that represents the discrete convolution of (2.19) in terms of near- and farfield contributions. Nearfield interactions are computed directly, where transformation of the sources would incur undue approximation errors. Meanwhile, farfield sources are transformed by projecting onto auxiliary basis functions representing points on a regular Cartesian grid. Farfield interactions are then computed via auxiliary sources, which offer two advantages: (i) they compress the interaction 13 matrix by representing sources within the same spatial region in terms of the same auxiliary set, and (ii) they impose a Toeplitz structure on the resulting interaction matrix that lends itself to efficient diagonalization through application of a multidimensional FFT. The projection to and from the auxiliary basis is accomplished via auxiliary matrices Λ and Λ† . In terms of this basis, they may be presented as: ∑︁ † 𝜓ℓ (r) ≈ Λℓu 𝛿(r − u). (2.23) u∈𝐶ℓ Here 𝜓ℓ (r) ∈ {sℓ (r) · x̂, sℓ (r) · ŷ, sℓ (r) · ẑ} and 𝐶ℓ denotes the collection of grid points within the expansion region of sℓ (r) (Figure 2.4 illustrates the relevant geometry). For an expansion of order 𝑀, the sum (2.23) contains (𝑀 + 1) 3 terms corresponding to the (𝑀 + 1) 3 grid points nearest to sℓ (r). Additionally, the border parameter 𝛾 separates near- and farfield interactions based upon 𝑅ℓℓ′ , which gives the minimum distance in units of the grid spacing between expansion regions enclosing dot ℓ and ℓ′ : grid 𝑅ℓℓ′ = min{∥u − u′ ∥ ∞ | u ∈ 𝐶ℓ , u′ ∈ 𝐶ℓ′ }. (2.24) † The matrix elements Λℓu are then determined via the multipole moment matching scheme: ∫ " # ∑︁ † (𝑥 − 𝑥 0 ) 𝑚 𝑥 (𝑦 − 𝑦 0 ) 𝑚 𝑦 (𝑧 − 𝑧0 ) 𝑚 𝑧 𝜓ℓ (r) − Λℓu 𝛿(r − u) d3 r = 0. (2.25) u∈𝑐 ℓ (r0 ≡ 𝑥0 x̂ + 𝑦 0 ŷ + 𝑧0 ẑ denotes the origin about which we compute the multipoles) and solving for those elements via least squares. In short, the algorithm effects the convolution (2.19) via: Ẽ (𝑚) ≈ ẼL(𝑚) + Ẽdirect (𝑚) (𝑚) + ẼFFT (2.26) where  (𝑚−𝑚 ′ ) (𝑚)  𝑚 ′ =0 F̃direct · à (𝑚) − ẼFFT  Í𝑚 𝑅ℓℓ′ ⩽ 𝛾   (𝑚)  Ẽdirect =  0 otherwise,     (2.27)  (𝑚)  Fℓℓ′ 𝑅ℓℓ′ ⩽ 𝛾   (𝑚)  F̃direct,ℓℓ ′ =  0 otherwise,    14 𝛾 Δ𝑠 r𝑏 𝛾 Δ𝑠 𝛾 Δ𝑠 r0 𝛾 Δ𝑠 (𝛾 + 1) Δ𝑠 𝛾 Δ𝑠 r𝑎 r𝑐 Figure 2.4 Illustration of the nearfield criterion for a third order expansion 𝑀 = 3 corresponding to (3 + 1) 2 = 16 gridpoints in this 2D illustration, and 𝛾 = 2. The dashed line indicates the complete nearfield of the box associated with r0 —i.e. all boxes that have an expansion point within 𝛾Δ𝑠 (infinity norm) of the expansion around r0 . Consequently, all of the sℓ (r) within the central dark blue square have a pairwise interaction with the sℓ′ (r) inside the dashed box. (𝑚) Farfield interactions are evaluated by computing the FFT-facilitated convolution ẼFFT between auxiliary sources, while nearfield interactions are evaluated by computing interactions between sources directly (as in (2.19)), then subtracting the contribution to the nearfield interaction due to (𝑚) ẼFFT . This procedure mitigates approximation errors between nearby auxiliary sources by replacing their contribution with contributions from the original sources. 2.3.1 Error analysis A simulation was run to analyze errors incurred in evaluating the convolution (2.19). We set up two domains with sufficient separation such that the interactions between these occur only via AIM. Each domain contains 64 randomly distributed dots, we prescribe the temporal variation of the polarization of each dot, and we measure the total radiated field at each dot. Finally, we fix the 15 temporal interpolation basis order at 3 and the polarization of each dots varies as (𝑡 −𝑡0 ) 2 − 𝑃(𝑡) = 𝑒 2𝜎 2 . (2.28) The simulation runs for 1024 timesteps of size Δ𝑡 = 0.1 ps, the width of the Gaussian 𝜎 = 1024 Δ𝑡/12 and its center 𝑡0 = 1024 Δ𝑡/2. This method admits a straightforward analytic solution via (2.10). ℓ2 norm differences between analytic and numerical solutions are then calculated as a function of grid size Δ𝑠 for different expansion orders. As Figure 2.5 suggests, we observe excellent convergence; indeed, the error scales as 𝑂 (Δ𝑠 𝑀−1 ), consistent with the results of previous convergence analysis [34]. 2.3.2 Performance analysis Next, we present a set of experiments that demonstrate the 𝑂 (𝑁 log(𝑁)) complexity scaling of AIM. For this, we perform simulations in both the fixed frame with prescribed polarizations, and the rotating frame with full Liouville equation dynamics. To ensure proper examination of computational complexity, we start with a box of side length 6 Δ𝑠 (chosen to minimize the number of nearfield pairs), and filled with dots at random locations. We obtain each successive value of 𝑁 by doubling the sidelength and in effect, increasing the number of dots by a factor of eight. We use a third order expansion 𝑀 = 3 with AIM spacings Δ𝑠 = 𝜆/400 and Δ𝑠 = 𝜆/10 for the fixed and rotating wave cases, respectively. Timesteps mirror those used in Section 2.3.1. Figure 2.6 gives runtimes for both cases, demonstrating that the two FFT-accelerated simulations outpace their direct counterparts near 𝑁 = 1000 and 𝑁 = 2000, respectively. 16 10−1 10−2 ℓ2 error (ps−1 ) 10−3 10−4 10−5 10−6 0.83𝑥 1.0 0.28𝑥 2.0 0.08𝑥 2.9 0.13𝑥 4.0 10−7 0.19𝑥 5.0 0.05 0.1 0.15 0.2 Δ𝑠/𝜆 Figure 2.5 𝑙2 error of the Rabi frequency magnitude | 𝜒| with respect to grid spacing Δ𝑠 for expansion orders 𝑀 = 2 through 6 and 𝛾 = 1, using source and observer boxes of volume 𝜆3 separated by Δ𝑟 = 2𝜆( 𝑥ˆ + 𝑦ˆ + 𝑧ˆ), each containing 64 randomly generated dots. For an expansion order 𝑀 one expects the overall error to scale as 𝑂 (Δ𝑠 𝑀−1 ). 17 105 104 103 Runtime (s) 102 101 100 10−1 3.671 44 × 10−5 𝑁 𝑠2.0701 (direct) 10−2 5.917 79 × 10−2 𝑁 𝑠1.00663 (nearfield + FFT) 10−3 101 102 103 104 Number of basis functions 106 105 104 Runtime (s) 103 102 101 100 5.224 31 × 10−4 𝑁 𝑠1.9852 (direct) 10−1 3.683 39 × 10−1 𝑁 𝑠1.1014 (nearfield + FFT) 10−2 101 102 103 104 Number of basis functions Figure 2.6 FFT runtime (excluding setup time) using a third-order expansion. (Top) 1024 timesteps with Δ𝑠 = 𝜆/400 and prescribed polarizations in the fixed frame. (Bottom) 1000 timesteps with Δ𝑠 = 𝜆/10 and Liouville-dynamics polarization in the rotating frame. Both cases have a quasi-quadratic scaling in the direct calculation, whereas the FFT-accelerated calculation performs slightly worse than linear. 2.4 Large scale simulations Having analyzed both its error scaling and performance complexity, we apply AIM to the simulation of large ensembles (𝑁 > 10000) of dots, where AIM performance outpaces that of the basic algorithm. Fig (2.7) shows a simulation employing AIM that reproduces features displayed in 18 a previous work. The figure displays the polarization of each dot in a cylinder as a function of their z-coordinate (the axis of the cylinder), under the effect of a resonant 𝜋 pulse. Each of the dots has an identical (fixed) dipole moment (see [2] for the details of the simulation parameters). Evidently, the secondary radiation produces random shifts in the polarization due to short-range effects in the local neighborhood of each dot. Also present is oscillation of the polarization due to long-range collective effects. This oscillation reflects the role of boundary conditions in the confinement of the macroscopic electric field in the system. 0.5 0.4 0.3 |𝜌01 | 0.2 0.1 0 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 𝑍 (𝜇𝑚) Figure 2.7 ẑ-distribution of polarization | 𝜌˜ 01 | for a 10 000-dot cylindrical simulation, replicating the parameters in [2]. The AIM calculation recovers the oscillatory long-range pattern that we obtained using a direct calculation [2]. In Figures 2.8-2.9 we examine the response of a system of 100 000 dots—randomly distributed throughout a cuboid—to an applied laser pulse of the form: (k·r− 𝜔 𝐿 (𝑡 −𝑡0 ) ) 2 − EL (r, 𝑡) = 𝐸 0 𝑒 2𝜎 2 cos(k · r − 𝜔 𝐿 𝑡) x̂. (2.29) The laser amplitude 𝐸 0 is chosen to produce a 𝜋-pulse on each dot in the absence of interactions. The transition dipole moment of each dot has a fixed magnitude but random orientation: d = 𝑑 ê, where ê is sampled from the uniform distribution on the unit sphere. Tables 2.1-2.2 list additional simulation parameters. 19 Quantity Symbol Value Speed of light 𝑐 300 𝜇m ps −1 Dipole moment 𝑑/𝑒 10 𝑎 0 Transition frequency 𝜔0 1500 meV/ℏ Laser frequency 𝜔𝐿 1500 meV/ℏ Laser wavevector k 𝜔 𝐿 /𝑐 ẑ Pulse width 𝜎/𝜔 𝐿 1 ps Pulse delay 𝑡0 5 ps Decoherence times 𝑇1 , 𝑇2 10 ps, 20 ps Table 2.1 Simulation parameters for Section 3.5. Here 𝑐, 𝑒, and 𝑎 0 denote the vacuum speed of light, elementary charge, and Bohr radius respectively. The decoherence times here, while shorter than those typical of optical resonance experiments, afford a shorter computational time but preserve dynamical emission phenomena Quantity Symbol Value Simulation timestep Δ𝑡 0.02 ps AIM spacing Δ𝑠 0.040𝜆 = 33.06 nm AIM expansion order 𝑀 5 Nearfield border parameter 𝛾 1 Transverse domain length — 16 Δ𝑠 = 529 nm Longitudinal domain length — 1500 Δ𝑠 = 49.59 𝜇m Table 2.2 AIM parameters for Section 3.5. Fig 2.8 displays a color map of | 𝜌˜ 01 | as an indicator of the polarization | P̃| of each quantum dot at different timesteps after the pulse peak. The figure shows only dots located in a central segment of about 4 µm of the entire cuboid. The random orientation of the dipole moments creates a variation in the amplitude of the polarization with dots whose dipole moments align or anti-align with the laser field having greatest amplitude. In addition, despite each dot resonantly coupling to the pulse, inhomogeneity arises due to the inter-dot coupling. These simulation can resolve inhomogeneities at the microscopic level, taking into account the orientation of the transition dipole moment of each quantum dot, as well as the effect of local secondary fields. 20 Figure 2.8 Coloration of | 𝜌˜ 01 | as an indicator of | P̃| at 𝑡1 = 6.0 ps (top), 𝑡2 = 7.0 ps (middle), 𝑡3 = 8.0 ps (bottom) relative to the peak of a 1 ps-wide pulse, for a system of 100 000 dots. To visualize long-range effects, Fig. (2.9) shows | 𝜌˜ 01 | as a function of the z coordinate of each dot, corresponding to the color plots of fig. 2.8. Here we show the entire cuboid having sides of 20 µm. In contrast to the results of fig. 2.7, no oscillations due to confinement are observed, since the length of the system far exceeds the radiation wavelength. Moreover, we observe a dispersion of the polarization due to the random orientation of the transition dipoles. Since the strength of the coupling scales with d · E ∝ cos(𝜃), the distribution peaks at the value of | 𝜌˜ 01 | when 𝜃 = 0 or 𝜃 = 𝜋, with a tail corresponding to all the intermediate values. Only a few dots, for which the secondary fields constructively interfere, have a polarization larger than the peak value. Finally, note how the value of the peak polarization slightly increases from left to right due to pulse propagation. 2.5 Concluding remarks In this chapter, we described a flexible, self-consistent model for simulating the dynamics of quantum dots. This model, wherein an integral operator electric field couples to density matrices representing the dot dynamics, permits the simulation of dense ensembles, characterized by strongly 21 0.5 0.48 0.46 0.44 0.42 |𝜌01 | 0.4 0.38 0.36 0.34 0.32 0.3 −30 −25 −20 −15 −10 −5 0 5 10 15 20 25 30 𝑍 (𝜇𝑚) 0.5 0.48 0.46 0.44 0.42 |𝜌01 | 0.4 0.38 0.36 0.34 0.32 0.3 −30 −25 −20 −15 −10 −5 0 5 10 15 20 25 30 𝑍 (𝜇𝑚) Figure 2.9 Scatterplots of | 𝜌˜ 01 | corresponding to the bottom two color plots of fig. 2.8. There exists a single preferred polarization, represented by the linear region of greatest density, arising from dots whose transition dipole moments align or anti-align with the laser field. Radiative coupling produces polarizations that exceed this value. The changes in the red line reflect pulse propagation. 22 nonlinear optical effects. To examine these effects in large ensembles, we applied the Adaptive Integral Method (AIM) to facilitate computing matrix-vector products involved in the discretization of the electric field. Simulations revealed both short range oscillatory behavior resulting from confinement of the field, and long range effects, including polarization dispersion. 23 CHAPTER 3 COLLECTIVE EMISSION EFFECTS IN ENSEMBLES OF QUANTUM DOTS A major application of nonlinear quantum optics is the analysis of collective emission effects such as superradiance and subradiance. Theoretical descriptions of these phenomena often rely on effective Hamiltonians, such as the Dicke model, where interaction with one or few cavity modes is assumed and emitters are homogeneous. For extended systems, Maxwell-Bloch equations can be used where the electric field couples to a continuous local-averaged polarization field. However, to understand these collective phenomena, it is essential to consider the role of the emitters’ local configuration and their spatial and energy distribution. In fact, the disorder in the local distribution of the emitters can strongly affect the superradiant and subradiant dynamics. The coupling between emitters resulting from the exchange of virtual photons, known as van-der-Waals coupling [9], has been recognized for a long time to be an obstacle to the experimental observation of superradiant behavior [36]. Far from being a limitation, subradiance has been recently proposed as a mechanism for photon storage in quantum memories [37, 38]. As superradiance and subradiance are highly sensitive to underlying parameters, yet their analysis provides crucial insights for other research areas, it becomes important to devise methods to replicate these effects numerically. In this chapter, we detail the theory of superradiance and subradiance, starting with Dicke’s original model, then elaborating into a quantized field Maxwell-Bloch formalism. From this formalism, we derive an expression for the photon emission rate as a summation of pairwise emission terms, and show that this implies superradiance. Using the Maxwell-Bloch model of the previous chapter, we run simulations that reproduce superradiance, which we verify by calculating the photon emission rate and showing that it scales quadratically with the number of emitters. Simultaneously, we observe subradiance in the steady-state solution of these simulations. We then examine patterns in the radiated field that result from collective coupling, and conclude with studies on the collective Lamb shift and inhomogeneous broadening, and their effects on emission dynamics. 24 3.1 Dicke model of superradiance In Dicke’s original model, an ensemble of 𝑁 identical emitters are represented by 𝑁 two-level states |𝜓𝑖 ⟩, with associated transition energy ℏ𝜔0 and transition operators 𝜎𝑖± and diagonal operator 𝜎𝑖𝑧 . In terms of transition operators, the electric dipole of the 𝑖th emitter can be represented as: d̂𝑖 = (𝜎𝑖+ + 𝜎𝑖− ) di (3.1) Rather than viewing each two-level state separately, we may consider the tensor product of all Ë𝑁 states as constituting a single ensemble state |𝜓⟩ with initial value |𝜓(𝑡 = 0)⟩ = 𝑖=1 |𝜓𝑖 (𝑡 = 0)⟩. In doing so, we make two crucial assumptions: (i) the volume in which the emitters are embedded has dimensions much smaller than the transition wavelength 𝜆 0 = 2𝜋𝑐/𝜔0 , and (ii) the coupling of the emitters to the radiation field, which governs the evolution of the ensemble, is symmetric with respect to the pairwise exchange of emitters. The first assumption ensures emitters are confined to sufficiently small a volume that they experience non-negligible effects due to their radiation fields. This strong coupling assumption drives the coherent alignment of dipoles that ultimately yields superradiance. The second assumption imposes an equivalence relation on the Hilbert space of all combinations of tensor products of states |𝜓𝑖 ⟩. The system is thus confined to evolve within a Hilbert subspace invariant under permutations of emitters, i.e. the emitters remain indiscernible with respect to their emissions at all times. Under this assumption, it is permissible to represent each two-level state as a spin-1/2 state, and the entire 𝑁-state ensemble as a symmetric combination of 𝑁 such states, i.e. a spin angular momentum eigenstate with maximum angular momentum 𝑆 = 𝑁/2. Defining the spin angular momentum operators: ∑︁ 𝑆ˆ± = 𝜎𝑖± (3.2) 𝑖 ∑︁ 𝑆ˆ𝑧 = 𝜎𝑖𝑧 (3.3) 𝑖 𝑆 = ( 𝑆ˆ+ 𝑆ˆ− + 𝑆ˆ− 𝑆ˆ+ )/2 + 𝑆ˆ2𝑧 ˆ2 (3.4) 25 these eigenstates may be labeled |𝑆, 𝑀⟩ with associated eigenvalues given by: 𝑆ˆ2 |𝑆, 𝑀⟩ = 𝑆(𝑆 + 1) |𝑆, 𝑀⟩ 𝑆ˆ𝑧 |𝑆, 𝑀⟩ = 𝑀 |𝑆, 𝑀⟩ where 𝑀 ∈ {−𝑆, −𝑆 + 1, ..., 𝑆 − 1, 𝑆}. Suppose that we initialize all emitters in the excited state |↑⟩, i.e. |𝜓𝑖 (𝑡 = 0)⟩ = |↑⟩ for all 𝑖. The ensemble thus begins in the state |𝑆, 𝑀 = 𝑆⟩ and descends down the “chain" of 2𝑆 + 1 discrete energy levels, emitting a photon of energy ℏ𝜔0 at each transition. The rate of emission may be calculated by treating each emitter as an independent dipole, so that the ensemble has total dipole moment: ∑︁ ∑︁ d̂ = d̂𝑖 = (𝜎𝑖+ + 𝜎𝑖− ) di (3.5) 𝑖 𝑖 The emission rate for a single emitter is given by 𝑊𝑖 = Γ0 ⟨𝜎𝑖+ 𝜎𝑖− ⟩, where Γ0 is the spontaneous emission rate (as in 2.13) and ⟨·⟩ denotes an expectation value. Using expressions (3.2-3.4), we find a corresponding expression for the emission rate from an ensemble of 𝑁 emitters: 𝑊𝑁 = Γ0 ⟨𝑆ˆ+ 𝑆ˆ− ⟩ = Γ0 (𝑆 + 𝑀)(𝑆 − 𝑀 + 1) (3.6) wherein it is evident that the emission rate starts at the low value 𝑁 Γ0 for 𝑀 = 𝑆, increases to a maximum Γ0 𝑆(𝑆 + 1) at 𝑀 = 0, then decreases to zero at the ground state 𝑀 = −𝑆. It is also informative to calculate correlations between different emitters due to their dipole moments. Assuming symmetry under pairwise exchanges, we may write ⟨𝜎𝑖+ 𝜎 𝑗− ⟩ = ⟨𝜎 + 𝜎 − ⟩ for all 𝑖, 𝑗. Thus: ∑︁ ⟨𝑆ˆ+ 𝑆ˆ− ⟩ = ⟨𝜎𝑖+ 𝜎𝑖− ⟩ + 𝑁 (𝑁 − 1)⟨𝜎 + 𝜎 − ⟩ (3.7) 𝑖 It follows that: 𝑆2 − 𝑀 2 ⟨𝜎 + 𝜎 − ⟩ = (3.8) 𝑁 (𝑁 − 1) which shows that dipole-dipole correlations grow as the system decays to the half-excited state 𝑀 = 0, then decrease as further energy is lost down to the ground state. Comparing to (3.6), we see that the enhancement of emission rates parallels the buildup of correlations amongst dipoles. 26 Notice, however, that the symmetry assumption is crucial to the derivation of this result, without which expression (3.7) would be invalid. For the simple case of 𝑁 = 2, explicit expressions for 𝑊 (𝑡) can be found, assuming either independent or correlated emissions. The result, as shown by Haroche [9], can be written: 𝑊2 (𝑡) = 2Γ0 𝑒 −2Γ0 𝑡 (1 + 2Γ0 𝑡) (3.9) 𝑊2(𝑖𝑛𝑑) (𝑡) = 2Γ0 𝑒 −Γ0 𝑡 (3.10) Curves for both types of emission are shown in Figure 3.1, where it is apparent that the coupled pair of emitters experiences a higher initial emission rate compared to the uncoupled pair. 1 0.8 Emission rate 0.6 0.4 0.2 0 0 0.5 1 1.5 2 2.5 3 3.5 4 𝑡 Figure 3.1 Comparison of photon emission rates for a coupled and uncoupled pair of emitters. Superradiance is described by the enhanced emission rate for the case of coupled emitters. 3.1.1 Subradiance Subradiance, in which destructive interference among emitters induces suppression of emission, may be expressed as a generalization of the Dicke model. Consider the state of all possible single excitation states |𝑘⟩ where only the 𝑘th emitter is excited: 1 ∑︁ + 1 ∑︁ |𝑆, 𝑀 = −𝑆 + 1⟩ = √ 𝜎𝑘 |↓⟩ = √ |𝑘⟩ (3.11) 𝑁 𝑘 𝑁 𝑘 27 which is clearly a spin angular momentum eigenstate, i.e. a Dicke state. However, the subradiant states: ! 𝑁+ 𝑁 1 ∑︁ ∑︁ |𝑆, 𝑀 = −𝑆 + 1⟩ − = √ |𝑘⟩ − |𝑘⟩ (3.12) 𝑁 𝑘=1 𝑘=𝑁0 +1 with 𝑁 − 𝑁+ destructively interfering emitters, do not lie in the Dicke subspace. More generally, consider the state of all possible excitation states with 𝑆 + 𝑀 emitters excited: ! 𝑆+𝑀 1 ∑︁ + 1 ∑︁ |𝑆, 𝑀⟩ = √ 𝜎𝑘 |↓⟩ = √ |k⟩ (3.13) N 𝑘 N k 𝑁! where k = (𝑘 1 , 𝑘 2 , ..., 𝑘 𝑆+𝑀 ) labels each of the N = (𝑆+𝑀)!(𝑆−𝑀)! possible states. The corresponding subradiant states are: N+ N ! 1 ∑︁ ∑︁ |𝑆, 𝑀⟩ − = √ |k⟩ − |k⟩ (3.14) N k k and are similarly no longer an element of the eigenspace. Indeed, the states |𝑆, 𝑀⟩ − form a non- orthogonal basis for the space of subradiant states. Moreover, as the Dicke states are characterized by their symmetry under pairwise exchange of particles, the emergence of subradiant states constitutes a form of symmetry breaking. Compared to the emission rate of a Dicke state with 𝑁 emitters (𝛾 = 𝑁 Γ0 ), the emission rate of a subradiant state is diminished: 𝛾 = (2𝑁+ − 𝑁)Γ0 and is precisely zero when 𝑁+ = 𝑁/2, i.e. when there are an equal number of constructively and destructively interfering emitters. In systems exhibiting superradiance, subradiance dominates the steady state, as the components of the superradiant eigenstate decay on transient timescales. This behavior, which is mediated by Lamb shifts (see Subsection 3.4.1), shall be replicated in the numerical simulations of Section 3.5. 3.2 Quantization of electric field The Dicke model offers a theoretical explanation of superradiance posed in the familiar terms of addition of angular momentum. Yet as a simplified model, it neither captures the intricate interaction between a system of emitters and their radiated fields, nor faithfully represents mixed state ensembles of emitters. To account for these deficiencies, we return to the Maxwell-Bloch formalism, and analyze the coupling of density matrices to a quantized electric field. As our 28 computational model is based on classical fields, our ultimate goal is to derive an expression for the photon emission rate, analogous to (3.7) of the Dicke model. This would establish the consistency of the two models, and allow us to quantify the observation of superradiance in terms of the emission rate. We prescribe a quantized electric field as the superposition of modes with wavevector k and polarization n̂, defined in terms of photon creation and annihilation operators 𝑎 †k,n̂ , 𝑎 k,n̂ : ∑︁ E+ (r) = −𝑖 E 𝑘,n̂ 𝑎 k,n̂ 𝑒𝑖k·r (3.15) k,n̂ ∑︁ E− (r) = 𝑖 E 𝑘,n̂ 𝑎 †k,n̂ 𝑒 −𝑖k·r (3.16) k,n̂ √︃ ℏ𝑐𝑘 where E 𝑘,n̂ = 𝐸 𝑘 n̂ and 𝐸 𝑘 = 2𝜖𝑉 is the per-photon field amplitude, and 𝑉 is a quantization volume. The total Hamiltonian may then be separated into three terms H = 𝐻𝑑𝑜𝑡 + 𝐻𝑟𝑎𝑑 + 𝐻𝑖𝑛𝑡 representing the dots, radiated field, and interactions respectively: ∑︁ 𝐻𝑑𝑜𝑡 = ℏ𝜔0 𝜎𝑖𝑧 (3.17) 𝑖 ∑︁ 𝐻𝑟𝑎𝑑 = ℏ 𝜔 𝑘 (𝑎 †k,n̂ 𝑎 k,n̂ + 1/2) (3.18) k,n̂ ∑︁ 𝐻𝑖𝑛𝑡 = − d𝑖 · (E+ (r𝑖 ) + E− (r𝑖 )) ⊗ (𝜎𝑖+ + 𝜎𝑖− ) 𝑖 (3.19) RWA ∑︁ + ≈− d𝑖 · (E (r𝑖 ) ⊗ 𝜎𝑖+ − + E (r𝑖 ) ⊗ 𝜎𝑖− ) 𝑖 where the last step invokes the Rotating Wave Approximation, retaining only energy-conserving terms (those in which a photon is created from decay to the ground state, and vice versa). The similarity with the Hamiltonian (2.3) of Section 2.1 is apparent, with 𝐻𝑑𝑜𝑡 comprising the diagonal elements and 𝐻𝑖𝑛𝑡 the off-diagonal coherences (The field Hamiltonian 𝐻𝑟𝑎𝑑 is absent, as the photon modes are decoupled from the density matrices in the semiclassical model). 29 3.3 Maxwell-Bloch equations Rather than tracking the evolution of the density operators, we may transfer the time dependence onto the observables, yielding the Heisenberg picture. It is in this representation that the Maxwell- Bloch equations can be derived, consistent with the semiclassical model. As the Heisenberg picture essentially “forgets" about the global configuration, it permits the definition of observables E± (r, 𝑡) representing local fields. To determine the evolution of the field, note that the time dependence of E± (r, 𝑡) is transferred in the Heisenberg picture onto the operators 𝑎(𝑡) and 𝑎 † (𝑡). Inserting these operators into Heisenberg’s equation of motion for observable 𝐴: d𝐴 𝑖 = − [ 𝐴, H ] (3.20) d𝑡 ℏ with H = 𝐻𝑑𝑜𝑡 + 𝐻𝑟𝑎𝑑 + 𝐻𝑖𝑛𝑡 as defined by eqs. (3.17-3.19) gives: d𝑎 k,n̂ 𝑖 = − [𝑎 k,n̂ (𝑡), 𝐻𝑟𝑎𝑑 + 𝐻𝑖𝑛𝑡 ] d𝑡 ℏ" # 𝑖 ∑︁ † ∑︁ ∑︁ † ′ = − 𝑎 k,n̂ , ℏ 𝜔 𝑘 ′ 𝑎 k′ ,n̂′ 𝑎 k′ ,n̂′ − 𝑖 d𝑖 · E 𝑘 ′ ,n̂′ 𝑎 k′ ,n̂′ 𝑒 −𝑖k ·r𝑖 𝜎𝑖− (3.21) ℏ k′ ,n̂′ 𝑖 k′ ,n̂′ 1 ∑︁ = −𝑖𝜔 𝑘 𝑎 k,n̂ (𝑡) − d𝑖 · E 𝑘,n̂ 𝑒 −𝑖k·r𝑖 𝜎𝑖− (𝑡) ℏ 𝑖 and similarly for 𝑎 † (𝑡). Additionally, let us define the local polarization and population inversion operators P(r, 𝑡) and 𝑊 (r, 𝑡): ∑︁ P± (r, 𝑡) = d 𝛿(r − r𝑖 ) 𝜎𝑖± (𝑡) (3.22) 𝑖 ∑︁ 𝑊 (r, 𝑡) = 𝛿(r − r𝑖 ) 𝜎𝑖𝑧 (𝑡) (3.23) 𝑖 from which the motion of P± (r, 𝑡) is given by: 𝜕 ± 𝑖 P (r, 𝑡) = − ([P± , 𝐻𝑑𝑜𝑡 ] + [P± , 𝐻𝑖𝑛𝑡 ]) 𝜕𝑡 ℏ 𝑖 = − (−ℏ𝜔0 P± (r, 𝑡)) ∓ 2d(d · E∓ (r, 𝑡))𝑊 (r, 𝑡) (3.24) ℏ 2𝑖 = 𝑖𝜔0 P± (r, 𝑡) ± d(d · E∓ (r, 𝑡))𝑊 (r, 𝑡) ℏ 30 Deriving (3.21-3.24) again with respect to time and equating terms yields the Maxwell-Bloch equations in the Heisenberg picture (see Appendix B for a derivation): 1 𝜕2   ± 1 ∓ 1 2 2 2 − ∇ 2 E (r, 𝑡) = ∇ × ∇ × P (r, 𝑡) = − (∇ − ∇∇)P∓ (r, 𝑡) (3.25) 𝑐 𝜕𝑡 𝜖 𝜖 This resembles the classical equations (from which (2.10) follows, see Appendix A): 1 𝜕2 1 1 𝜕2     2 − ∇ E(r, 𝑡) = − − ∇∇ P(r, 𝑡) (3.26) 𝑐2 𝜕𝑡 2 𝜖 𝑐2 𝜕𝑡 2 with the difference 1/𝑐2 𝜕𝑡2 → ∇2 on the RHS, a discrepancy which has been explored, for instance in [39]. This discrepancy vanishes when the fields are cast in integral operator form (2.10), which is equivalent to either set of equations (see Appendices B) . Additionally, the mathematical objects representing the fields E+ and E− (and polarization sources P+ and P− ) in (3.25) become non-commuting operators acting on photon states. 3.4 Photon emission rate, Superradiance We now derive an expression for the emission rate, analogous to expression (3.6) of the Dicke model, via the Heisenberg equation of motion. To start, reconsider (3.21) and transform to the rotating basis 𝑎 𝑘 (𝑡) = 𝑎˜ 𝑘 (𝑡)𝑒 −𝑖𝜔 𝑘 𝑡 and 𝜎𝑖− (𝑡) = 𝜎 ˜ 𝑖 (𝑡)𝑒 −𝑖𝜔 𝐿 𝑡 . The first term is eliminated: d𝑎˜ k,n̂ ∑︁ 𝑖 −𝑖(k·r𝑖 −(𝜔 𝑘 −𝜔 𝐿 )𝑡) − =− 𝜒𝑘, n̂ 𝑒 𝜎˜ 𝑖 (𝑡) (3.27) d𝑡 𝑖 where 𝜒𝑘,𝑖 = d · E /ℏ is the Rabi energy of the 𝑖th emitter due to a radiated field of mode (𝑘, n̂). n̂ 𝑖 𝑘,n̂ Integrating, and assuming that no photons exist initially (𝑎˜ k,n̂ (𝑡 = 0) = 0) yields: ∑︁ ∫ 𝑡 ′ 𝑎˜ k,n̂ (𝑡) = − 𝜒𝑘,𝑖 n̂ 𝑒 −𝑖k·r𝑖 𝑑𝑡 ′ 𝑒𝑖(𝜔 𝑘 −𝜔 𝐿 )𝑡 𝜎 ˜ 𝑖− (𝑡 ′) (3.28) 𝑖 0 The rate of photon emission may then be equated to the gain in number of photons, defined by the photon number operator: ∑︁ 𝑁 𝑝 (𝑡) = 𝑎 †k,n̂ 𝑎 k,n̂ (3.29) k,n̂ 31 Substituting in (3.28) gives (R = r𝑖 − r 𝑗 ) ∑︁ ∑︁ ∫ 𝑡 ∫ 𝑡 −𝑖k·R ′ −𝑖(𝜔 𝑘 −𝜔 𝐿 )𝑡 ′ + ′ ′′ 𝑁 𝑝 (𝑡) = 𝑖 ∗ 𝑗 ( 𝜒𝑘,n̂ ) 𝜒𝑘,n̂ 𝑒 𝑑𝑡 𝑒 𝜎˜ 𝑖 (𝑡 ) 𝑑𝑡 ′′ 𝑒𝑖(𝜔 𝑘 −𝜔 𝐿 )𝑡 𝜎 ˜ 𝑗− (𝑡 ′′) k,n̂ 𝑖, 𝑗 0 0 ∫ ∫ 𝑡 ∫ 𝑡 𝑐 ∑︁ −𝑖k·R ′ ′ ′′ ) = d𝑖 · 𝑑k 𝑘 (I − k̂ ⊗ k̂)𝑒 · d𝑗 𝑑𝑡 𝑑𝑡 ′′ 𝑒 −𝑖(𝜔 𝑘 −𝜔 𝐿 ) (𝑡 −𝑡 × 16𝜋 3 𝜖ℏ 𝑖, 𝑗 0 0 𝜎˜ 𝑖+ (𝑡 ′) 𝜎 ˜ 𝑗− (𝑡 ′′) (3.30) Í 𝑉 ∫ where we have replaced k → (2𝜋) 3 𝑑k and used the completeness relation: ∑︁ n̂n̂ = I − k̂ ⊗ k̂ (3.31) n̂ ∇∇ Further replacing I − k̂ ⊗ k̂ → I + 𝑘2 and integrating over angles gives: ∞ 𝜔3𝑘 ∫   ∫ 𝑡 ∫ 𝑡 1 ∑︁ ∇∇ sin 𝑘 𝑅 ′ ′ ′′ 𝑁 𝑝 (𝑡) = 3 d𝑖 · 𝑑𝜔 𝑘 I+ 2 · d𝑗 𝑑𝑡 𝑑𝑡 ′′ 𝑒 −𝑖(𝜔 𝑘 −𝜔 𝐿 ) (𝑡 −𝑡 ) × 8𝜋 𝜖ℏ 𝑖, 𝑗 0 𝑐3 𝑘 𝑅 0 0 ˜ 𝑖+ (𝑡 ′) 𝜎 𝜎 ˜ 𝑗− (𝑡 ′′) (3.32) To eliminate one time integral, the Wigner-Weisskopf approximation is invoked, in which the term ˜ 𝑖+ (𝑡 ′) 𝜎 𝜎 ˜ 𝑗− (𝑡 ′′) is assumed to vary on a timescale much smaller than 𝜔−1 𝐿 , thus is effectively constant. Thus: ∫ 𝑡 ∫ ∞ ′′ −𝑖(𝜔 𝑘 −𝜔 𝐿 ) (𝑡 ′ −𝑡 ′′ ) ′ ′′ ) 𝑑𝑡 𝑒 ˜ 𝑖+ (𝑡 ′) 𝜎 ˜ 𝑗− (𝑡 ′′) 𝜎 ≈ ˜ 𝑖+ (𝑡 ′) 𝜎 ˜ 𝑗− (𝑡 ′) 𝜎 𝑑𝑡 ′′ 𝑒 −𝑖(𝜔 𝑘 −𝜔 𝐿 ) (𝑡 −𝑡 0 0 (3.33) ≈ 𝜋𝛿(𝜔 𝑘 − 𝜔𝐿 ) 𝜎 ˜ 𝑖+ (𝑡 ′) 𝜎˜ 𝑗− (𝑡 ′) where the Cauchy Principle Value (PV) of the integral over 𝑡 ′′ yields a Lamb shift, and will be considered in the following subsection. Ignoring the PV, the result is: 𝜔3𝐿 ∑︁   ∫ 𝑡 ∇∇ sin 𝑘 𝑅 𝑁 𝑝 (𝑡) = 2d𝑖 · I + 2 · d𝑗 𝑑𝑡 ′ 𝜎 ˜ 𝑖+ (𝑡 ′) 𝜎˜ 𝑗− (𝑡 ′) (3.34) 4𝜋𝜖ℏ𝑐3 𝑖, 𝑗 𝑘 𝑅 0 The photon emission rate is thus: 𝜔3𝐿 ∑︁   d𝑁 𝑝 (𝑡) ∇∇ sin 𝑘 𝑅 𝛾(𝑡) = = 3 2d𝑖 · I + 2 · d𝑗 𝜎 ˜ 𝑖+ (𝑡) 𝜎˜ 𝑗− (𝑡) d𝑡 4𝜋𝜖ℏ𝑐 𝑖, 𝑗 𝑘 𝑅 ∑︁ (3.35) = 2 𝜖ˆ𝑖 · 𝚪𝑖 𝑗 · 𝜖ˆ 𝑗 ⟨𝜎𝑖+ (𝑡)𝜎 𝑗− (𝑡)⟩ 𝑖𝑗 32 where d𝑖 = 𝑑 𝜖ˆ𝑖 , and the expectation value is taken with respect to the initial state of the system |Ψ0 ⟩: ⟨𝜎𝑖+ (𝑡)𝜎 𝑗− (𝑡)⟩ = ⟨Ψ0 |𝜎𝑖+ (𝑡)𝜎 𝑗− (𝑡)|Ψ0 ⟩ (3.36) Additionally:    3Γ0 sin 𝑘 𝑅 cos 𝑘 𝑅 sin 𝑘 𝑅 𝚪𝑖 𝑗 = (1 − r̂ ⊗ r̂) + (1 − 3r̂ ⊗ r̂) − (3.37) 4 𝑘𝑅 (𝑘 𝑅) 2 (𝑘 𝑅) 3 𝑑 2 𝜔3𝐿 𝜇0 𝑑 2 𝜔3𝐿 Γ0 = = (3.38) 3𝜋𝜖ℏ𝑐3 3𝜋ℏ𝑐 where we have recovered expression (2.13) for the spontaneous decay rate Γ0 , and 𝚪𝑖 𝑗 is a correlation tensor that describes the pairwise rate of emission induced by the interaction between the 𝑖th and 𝑗th emitters. This is related to the total emission rate of the dot ensemble via (3.35). To see that (3.35) implies superradiance, we separately write the diagonal and off-diagonal elements of the RHS: ∑︁ ∑︁ 𝛾(𝑡) = 2 Γ𝑖𝑖 ⟨𝜎𝑖+ 𝜎𝑖− ⟩ +2 Γ𝑖 𝑗 ⟨𝜎𝑖+ 𝜎 𝑗− ⟩ (3.39) 𝑖 𝑖≠ 𝑗 where Γ𝑖 𝑗 = 𝜖ˆ𝑖 · 𝚪𝑖 𝑗 · 𝜖ˆ 𝑗 and it can easily be shown that lim 𝑅→0 Γ𝑖𝑖 = Γ0 /2. Thus: ∑︁ ∑︁ 𝛾(𝑡) = Γ0 ⟨𝜎𝑖+ 𝜎𝑖− ⟩ + 2 Γ𝑖 𝑗 ⟨𝜎𝑖+ 𝜎 𝑗− ⟩ 𝑖 𝑖≠ 𝑗 (3.40) = 𝛾0 (𝑡) + 𝛾 𝐼 (𝑡) which is reminiscent of expression (3.7) of the Dicke model. Eq. 3.40 asserts that the total emission consists of independent emissions 𝛾0 , plus emission from coupling due to secondary radiated fields 𝛾 𝐼 . Moreover, the magnitude of this off-diagonal term scales approximately as 𝑁 (𝑁 − 1); hence we uncover the characteristic quadratic scaling of superradiance. 3.4.1 Collective Lamb shift Previously, we had discarded the Cauchy Principle Value of the Wigner-Weisskopf approximated integral (3.33), as we were solely interested in calculating the photon emission rate. Evaluating the PV by contour integration yields the following expression (see also Appendix C): 33 Quantity Symbol Value Refractive index 𝑛 1.54 Speed of light 𝑐 = 𝑐 0 /𝑛 194.6704 𝜇m ps −1 (Reference) Dipole moment 𝑑0 /𝑒 2.536 × 10−3 𝜇m Transition frequency 𝜔0 3175 meV/ℏ Laser frequency 𝜔𝐿 3175 meV/ℏ Laser wavevector k 𝜔 𝐿 /𝑐 ẑ Pulse width 𝜎/𝜔 𝐿 0.3446 ps Pulse delay 𝑡0 5 ps Decoherence time 𝑇1 8.31 ps Table 3.1 Simulation parameters for Section 3.5. Here 𝑐 0 and 𝑒 are the vacuum speed of light and elementary charge, respectively. The decoherence (spontaneous emission) time 𝑇1 = 1/Γ0 is obtained as the inverse of the spontaneous emission parameter Γ0 (2.13), hence is proportional to 𝑑 −2 .    3Γ0 cos 𝑘 𝑅 sin 𝑘 𝑅 cos 𝑘 𝑅 𝛀𝑖 𝑗 (r) = −(1 − r̂ ⊗ r̂) + (1 − 3r̂ ⊗ r̂) + (3.41) 4 𝑘𝑅 (𝑘 𝑅) 2 (𝑘 𝑅) 3 It can be shown that (3.41) induces a collective energy shift between pairs of emitters 𝑖 and 𝑗, corresponding to an additional term for the emitter Hamiltonian (3.17), which is a Lamb shift. This will become more apparent in Chapter 4, where we recast this phenomenon in the formalism of the Master equation. In the context of collective emissions, the Lamb shift has the effect of degrading subradiance [40]; that is, the purely subradiant states (3.12 or 3.14) are split into a combination of subradiant and superradiant states. The latter decay on transient timescales, leaving the slowly decaying subradiant states which are observed in the steady state. In Subsection 3.5.3, we exemplify purely subradiant states that occur in the absence of Lamb shifts. 3.5 Numerical simulation Here we describe the results of simulations that, using the Maxwell-Bloch model expounded thus far, demonstrate the consequences of collective emission. As in Section 2.4, we apply an incident field with the shifted Gaussian waveform: (k·r− 𝜔 𝐿 (𝑡 −𝑡0 ) ) 2 − EL (r, 𝑡) = 𝐸 0 𝑒 2𝜎 2 cos(k · r − 𝜔 𝐿 𝑡) x̂. (3.42) to excite an ensemble of dots lying initially in the ground state ( 𝜌˜ 00 , 𝜌˜ 01 )| 𝑡=0 = (1, 0). The dots are Gaussian distributed in space with a standard deviation along each dimension of 0.5𝜆. The 34 laser amplitude 𝐸 0 is adjusted to produce a 𝜋-pulse on each dot in the absence of interactions. The systems of 𝑁 quantum dots are assumed to be embedded in a NaCl medium with refractive index 𝑛. Each dot has an identical dipole moment d = 𝑑 x̂, which is varied based on a reference dipole moment 𝑑0 . For 𝑑 = 𝑑0 , decay times 𝑇1 and 𝑇2 = 2𝑇1 are chosen, and modified for other values of 𝑑 to satisfy the 𝑑 −2 dependence [41]. Values for simulation parameters, which mirror those in an experimental study [42], are shown in Table 3.1. 3.5.1 Population excitation Fig. 3.2a depicts the time behavior of a set of ten dots chosen from a simulation with 𝑁 = 200 dots, portraying a rich phenomenology of oscillations following the initial excitation. The figure shows the excited state population of each dot as a function of time, i.e. 𝜌11 𝑖 (𝑡). These oscillations result from local energy shifts induced by randomly distributed neighboring dots and are dominated by the 1/𝑅 3 contribution from Eq. (2.10). After excitation, we observe a superradiant behavior, which occurs on the time scale of 𝑇1 . This can be seen in Fig. 3.2b where we show the population dynamics averaged over all dots, ⟨𝜌11 ⟩ for different 𝑁. The figure and inset show a faster re-emission after pulse excitation in configurations with a greater density of emitters. To better quantify this behavior, we evaluate the off-diagonal contribution 𝛾 𝐼 to the emission rate via Eq. (3.40). However, as the response to the incident pulse obfuscates a direct calculation, it is more enlightening to obtain 𝛾 𝐼 by subtracting independent emissions 𝛾0 from the total emission rate 𝛾, which we equate to the rate of decay in total number of excited states 𝑁 (𝑡) = 𝑖 ⟨𝜎𝑖+ (𝑡)𝜎𝑖− (𝑡)⟩ = 𝑖 𝜌11 Í Í 𝑖 (𝑡). That is: ! d𝑁 (𝑡) ∑︁ d𝜌𝑖 (𝑡) 11 𝑖 𝛾 𝐼 (𝑡) = − − 𝛾0 (𝑡) = − + Γ0 𝜌11 (𝑡) (3.43) d𝑡 𝑖 d𝑡 Figure 3.3 shows the results of calculating both the total and off-diagonal emission rates. Apparent is the quadratic scaling of superradiant emissions with 𝑁. In contrast, the total emission does not scale quadratically, but becomes increasingly dominated by the superradiant contribution. As the superradiant states decay transiently, the system settles into a subradiant regime where re- emission slows down. This transition is visible in both per-dot (Fig. 3.2a) and averaged (Fig. 3.2b) plots. Higher emitter density also leads to more subradiance, resulting in a larger population 35 decaying at long times. Increasing the dipole strength 𝑑, while increasing the strength of interactions, induces shorter decay times 𝑇1 . The overall effect, as captured by Figure 3.4, is an enhancement of superradiance, at the cost of population of slowly decaying states. Figure 3.5 summarizes the effects of 𝑁 and 𝑑 on the average population at 1 ns. The subradiant slowly decaying states are enhanced by larger 𝑁 and lower 𝑑, corresponding to dense systems of weakly interacting dots. 3.5.2 Radiated field patterns Other signatures of superradiance, such as anisotropy, emerge when examining the radiated field of Eq. (2.10). Field intensities in the region external to the dot cloud show patterns with a period of half the wavelength (Fig. 3.7-3.6). The intensity along the 𝑦-axis, shown in Figure 3.7, indicates that the patterns persist for stronger dipole values. As the excitation pulse propagates along 𝑧, we also observe a clear enhancement of the emission along the positive 𝑧 axis, a well-known signature of collective radiative emission [10]. It is worth noting that the three terms of Eq. (2.10) are not commensurate. While the near-field interaction terms produce a random pattern near the origin, the far-field 1/𝑅 term produces a regular pattern characterized by phase shifts in some directions. Time-space plots reveal the synchronization of groups of dots into temporal oscillations in Fig. 3.8, where we plot the intensity along the 𝑦-axis as a function of time. These oscillations become more pronounced and irregular with increased dipole strength and dot density. This effect is captured in Fig. 5, which illustrates temporal and spatial Fourier transforms of the field intensity along the 𝑦-axis for different dipole strengths. These plots display considerable spectral broadening with increasing dipole strength. In time, this broadening is suggested by the emergence of additional peaks corresponding to new oscillation periods in the dot ensemble. Peaks corresponding to characteristic lengths also appear in space, which are nonetheless strongly dependent on the random spatial configuration chosen. 36 100 10−2 𝜌11 10−4 10−6 10−8 (a) 0 100 200 300 400 500 600 700 800 900 1,000 100 10−1 10−2 10−3 ⟨𝜌11 ⟩ 10−4 10−5 10−6 (b) 10−7 0 100 200 300 400 500 600 700 800 900 1,000 100 10−1 ⟨𝜌11 ⟩ 10−2 (c) 10−3 0 10 20 30 40 50 60 70 Time ( 𝑝𝑠) Figure 3.2 Population excitation 𝜌11 for a Gaussian distribution of 𝑁 dots with 𝑑/𝑑0 = 1.0 and various 𝑁. (a) Per-dot time evolution for 𝑁 = 200. We show here 𝜌11 𝑖 , with 𝑖 corresponding to ten randomly chosen dots. Data for the three dots closest to the center are shown in color. (b) Dot-averaged values ⟨𝜌11 ⟩ for different 𝑁. (c) Immediately post-excitation, a faster decay is observed for larger 𝑁, as evidenced by the crisscrossing of curves. 37 500 Emission rate 𝛾 𝐼 (𝑡) (units of Γ0 ) 400 300 200 100 0 0 5 10 15 20 25 30 35 40 Time ( 𝑝𝑠) Figure 3.3 Total emission rates 𝛾(𝑡) (dashed) and off-diagonal emission rates 𝛾 𝐼 (𝑡) (solid) calculated from (3.43), for the simulations of Figure 3.2. Only the latter, which is associated with superradiant emissions, displays quadratic scaling with number of emitters 𝑁. The superradiant burst lies between 𝑡 ≈ 5 and 15 ps, independently of 𝑁. Overall, the curves bear a strong resemblance to Figure 1.2. 100 10−1 10−2 10−3 𝜌11 10−4 10−5 10−6 10−7 0 100 200 300 400 500 600 700 800 900 1,000 Time ( 𝑝𝑠) Figure 3.4 Time evolution of space-averaged population excitation 𝜌11 for 𝑁 = 200 at various dipole strengths. For larger 𝑑 we observe: (a) a faster decay post-excitation and (b) a reduction in the number of subradiant, slowly-decaying states, as suggested by smaller populations at 1 ns. 38 2 −4 1.5 −5 𝑑/𝑑0 1 −6 0.5 50 100 150 200 250 300 350 400 log⟨𝜌11−7 ⟩ 𝑁 Figure 3.5 Log-plot of averaged-over-trials population ⟨𝜌11 ⟩ at 1 ns as a function of dipole strength and 𝑁. log |E| 2 −10 −5 0 5 10 y () 2 2 𝑦 (𝜇m) 0 𝑦 (𝜇m) 0 −2 (a) −2 (c) 0 50 100 150 200 0 50 100 150 200 2 2 𝑧 (𝜇m) 0 𝑦 (𝜇m) 0 −2 (b) −2 (d) 0 100 50 150 200 0 50 100 150 200 Time (ps) Time (ps) Figure 3.8 Plot of logarithmic radiated field intensities for 𝑁 = 100, 𝑑/𝑑0 = 1.0 (a, b), 𝑁 = 200, 𝑑/𝑑0 = 1.0 (c), 𝑑/𝑑0 = 2.0 (d) in time and space. Evident are not only spatial but temporal oscillations, becoming enhanced for larger values of the dipole moment. It is also evident that the intensity amplitude increases with the number of emitters. Groups of dots in the cloud (𝑦 ∼ 0) undergo emission synchronization leading to periodic oscillations that become more irregular as the density increase. Finally, (b) displays emission enhancement in the laser propagation direction (cf. Fig 3). 39 log |E| 2 y (𝜇m) −20 −10 0 10 20 1 0.5 y (𝜇m) 0 −0.5 −1 −1 −0.5 0 0.5 1 1 0.5 z (𝜇m) 0 −0.5 −1 −1 −0.5 0 0.5 1 x (𝜇m) Figure 3.6 Colormaps of logarithmic field intensity (field norm squared) for a configuration with 𝑁 = 100 and 𝑑/𝑑0 = 1.00 after 20 ps, on the 𝑥 − 𝑦 (top) and 𝑥 − 𝑧 (bottom) planes. The spatial oscillations occur with a period about half the wavelength of 253 nm. Also note the enhancement along the laser propagation direction in the positive 𝑧 axis. 40 5 0 −5 −10 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 𝑦 (𝜇m) Figure 3.7 Plot of logarithmic field intensity along the 𝑦 axis for 𝑑/𝑑0 = 1.0 (Blue) and 𝑑/𝑑0 = 2.0 (Magenta), corresponding to Fig. 3 (a). The Fourier data is normalized against the 𝑑/𝑑0 = 2.0 series. 3.5.3 Collective Lamb shift Next we consider the effect of disregarding the pairwise Lamb shift, represented by (3.41). Comparing (3.41) to (2.12) we notice that 𝛀𝑖 𝑗 arises from the real part of the Rabi energy via: ˜ P̃(r 𝑗 )} = 𝜖ˆ𝑖 · (𝛀𝑖 𝑗 + 𝑖𝚪𝑖 𝑗 ) · 𝜖ˆ 𝑗 𝜒 𝑗 (r𝑖 ) = d · 𝔉{ (3.44) where 𝜖ˆ𝑖 is the unit dipole vector for the 𝑖th dot. Hence writing our radiated field as E𝑟𝑎𝑑 (r, 𝑡) = ˜ P̃(r, 𝑡)} , where 𝔉{ ˜ P̃(r, 𝑡)} is defined as in (2.11), effectively excludes the shift. Results  𝑖 Im 𝔉{ of this simulation are shown in Figures 3.10 and 3.11; parameters follow those in Subsection 3.5.1. Evident is the appearance of purely subradiant states in the absence of Lamb shifts, characterized by a complete suppression of emission, thus a higher population of excited states at long times. We further observe from Figure 3.11 that superradiance is degraded by the shift as well. 3.5.4 Inhomogeneous broadening Lastly, we studied the effect of inhomogeneous broadening by considering dots with energy ℏ𝜔𝑖0 that follow a Gaussian distribution of width 𝛿 centered at 𝜔 𝐿 . Increasing 𝛿 affects the dynamics in two ways: (1) the excitation induced by the 𝜋 pulse is less efficient, so the population inversion decreases, and (2) the population of subradiant modes increases due to increased disorder. These two competing effects can be seen in Figure 3.12. The average excitation of the dots right 41 100 10−1 | 𝐼˜1 | 10−2 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 𝑓 (ps−1 ) 100 10−1 | 𝐼˜2 | 10−2 10−3 10−4 −200 −150 −100 −50 0 50 100 150 200 𝑘 (𝜇m−1 ) Figure 3.9 Temporal (top) and spatial (bottom) Fourier plots of field intensity 𝐼 (𝑦, 𝑡) = |E(𝑦, 𝑡)| 2 along the 𝑦 axis (as in Fig. 4) for 𝑁 = 200 and various dipole strengths, normalized by maximum intensity: 𝐼˜1 ( 𝑓 ) = 𝐼˜(𝑘 = 0, 𝑓 )/ 𝐼˜(0, 0) (top), 𝐼˜2 (𝑘) = 𝐼˜(𝑘, 𝑓 = 0)/ 𝐼˜(0, 0) (bottom). Spectral broadening with increasing dipole strength is evident in both cases. 42 100 10−2 10−4 𝜌11 10−6 10−8 10−10 0 100 200 300 400 500 600 700 800 900 1,000 Time ( 𝑝𝑠) Figure 3.10 Averaged population excitation curves for 𝑁 = 200 including and excluding the collective Lamb shift. In the latter (Magenta dashed), the occurrence of purely subradiant states is suggested by the complete suppression of emission following the period of superradiant decay. For reference, a plot for uncoupled emitters is also shown, exhibiting only spontaneous decay. 150 Emission rate 𝛾 𝐼 (𝑡) (units of Γ0 ) 100 50 0 −50 0 5 10 15 20 25 30 35 40 Time ( 𝑝𝑠) Figure 3.11 Off-diagonal emission rates 𝛾 𝐼 (𝑡) calculated from (3.43) for 𝑁 = 200, including and excluding the collective Lamb shift. We observe that superradiant emissions (𝑡 ≈ 5 and 15 ps) are diminished by the shift. 43 after excitation decreases with 𝛿, but at longer times, its dependence on the inhomogeneity is characterized by a peak around ℏ𝛿 = 0.1 meV. The resonant pulse is chosen to peak at 𝑡0 = 5 ps. 100 10−1 10−2 10−3 10−4 ⟨𝜌11 ⟩ 10−5 10−6 10−7 10−8 10−9 0 100 200 300 400 500 600 700 800 900 1,000 Time ( 𝑝𝑠) Figure 3.12 Time evolution of space-averaged population excitation ⟨𝜌11 ⟩ for 𝑁 = 200, for different inhomogeneous broadening values 𝛿 (in meV). (Inset) Excitation values as a function of 𝛿, at 𝑡 = 500 ps (blue circles) and 𝑡 = 1000 ps (red squares). 3.6 Concluding remarks Using the Maxwell-Bloch model, we discovered signatures of collective emissions in both the population excitation and radiated fields of medium-sized ensembles of quantum dots. A crucial step of our study involved calculation of the photon emission rate, which we used to measure the scaling of superradiant emission with number of dots. Our expression for the emission rate (3.40) was derived in the Heisenberg picture, in which the fields are assumed to be time dependent, while the density operators representing the ensemble state are held constant. Of course, it is equally valid to transfer the time dependence to the density matrices, leading to the Schrödinger picture. While the two pictures offer an equivalent description of superradiance in terms of secondary radiated fields, their regimes of validity are dependent on both the spatial and temporal characteristics of the system. Details are considered in the next chapter, where we compare superradiant systems governed by the Maxwell-Bloch equations (3.25) to that under the Master equation–originating from the Schrödinger picture. 44 CHAPTER 4 MASTER EQUATION DESCRIPTION OF EMISSION DYNAMICS In the last chapter, we presented the theory and simulation of collective emission in the language of the Maxwell-Bloch equations. As it is grounded on Heisenberg’s equation of motion for operators representing photons and local polarization, the Maxwell-Bloch approach is innately well-suited for capturing memory and propagation effects due to electromagnetic fields arising from localized emitters. In contrast, a global description of all emitters and their aggregate quantum state requires an alternative approach based on the Schrödinger picture. This motivates the Master Equation, which prescribes the evolution of the global 2𝑁 × 2𝑁 density matrix representing the quantum state of all emitters in a large Hilbert space. In this chapter, we first derive the Master Equation from the Schrödinger equation of motion. This entails the Born-Markov approximation [43, 44], in which correlations between quantum dots and their radiated fields are assumed to occur on a short timescale. Afterward, we compare the results of simulations using both approaches. The purpose of this is threefold. First, we demonstrate that our numerical model is consistent with established algorithms for solving the Master equation. Second, we numerically demonstrate that the Born-Markov approximation is equivalent to disregarding effects due to the finite propagation speed of the electromagnetic field. Lastly, we elucidate a method for inducing decays in perfectly inverted states, which corresponds to an unstable equilibrium for the Maxwell-Bloch model. We conclude with a discussion of merits and drawbacks of both approaches. 4.1 Master equation formulation Ë𝑁 Consider a global density matrix 𝜚 for 𝑁 dots with initial state 𝜚(𝑡0 ) = 𝑖=1 𝜌𝑖 (𝑡0 ), which also tracks the degrees of freedom of the field: 𝜚(𝑡) = 𝜌(𝑡) ⊗ 𝜌𝑟𝑎𝑑 (𝑡) (4.1) In the Schrödinger picture, the Liouville equation of motion for 𝜚 is given by: d𝜚 −𝑖 = [H (𝑡), 𝜚] (4.2) d𝑡 ℏ 45 where the total Hamiltonian H = 𝐻𝑑𝑜𝑡 + 𝐻𝑟𝑎𝑑 + 𝐻𝑖𝑛𝑡 is defined by (3.17-3.19). Via the unitary operator 𝑈 = exp(−𝑖(𝐻𝑑𝑜𝑡 + 𝐻𝑟𝑎𝑑 )𝑡/ℏ), we transform to the interaction picture: d 𝜚˜ −𝑖 = [𝐻𝑖𝑛𝑡 (𝑡), 𝜚] ˜ (4.3) d𝑡 ℏ Tracing out over field variables, we obtain the reduced density matrix: ˜ = Tr𝑟𝑎𝑑 𝜚(𝑡) 𝜌(𝑡) ˜ which enables us to examine the evolution of the emitters independently of the degrees of freedom of the field. Integrating eq. (4.2), and iterating the procedure to second-order yields: ∫ 𝑡 d −1 ˜ = 2 Tr𝑟𝑎𝑑 𝜌(𝑡) 𝑑𝑡 ′ [𝐻𝑖𝑛𝑡 (𝑡), [𝐻𝑖𝑛𝑡 (𝑡 − 𝑡 ′), 𝜚(𝑡 ˜ − 𝑡 ′)]] (4.4) d𝑡 ℏ 0 As the right-hand-side still contains the total density operator 𝜚, ˜ the Born approximation is invoked, wherein the coupling between emitters and field is assumed sufficiently weak that the field density operator is constant, and 𝜚(𝑡) ˜ ⊗ 𝜌˜𝑟𝑎𝑑 . Integration still requires knowledge of the state of the ˜ = 𝜌(𝑡) system at all previous times, so a further simplification is made. Assuming that the field correlation lasts much shorter than the time in which the system varies appreciably, the Markov approximation replaces 𝑡 − 𝑡 ′ → 𝑡, yielding the Born-Markov (or Redfield) master equation: ∫ ∞ d −1 ˜ = 2 Tr𝑟𝑎𝑑 𝜌(𝑡) 𝑑𝑡 ′ [𝐻𝑖𝑛𝑡 (𝑡), [𝐻𝑖𝑛𝑡 (𝑡 − 𝑡 ′), 𝜌(𝑡) ˜ ⊗ 𝜌˜𝑟𝑎𝑑 ]] (4.5) d𝑡 ℏ 0 where we have assumed that the total system is initialized as the product of a field in its vacuum state with all emitters excited: Ö ˜ 𝜚(0) = 𝜌(0) ˜ ⊗ 𝜌˜𝑟𝑎𝑑 = |+⟩𝑖 ⟨+|𝑖 ⊗ |0⟩ ⟨0| (4.6) 𝑖 We may now substitute in the expression for 𝐻𝑖𝑛𝑡 in terms of field modes and transform back to the Schrödinger picture via 𝜌(𝑡) = 𝑈 𝜌(𝑡)𝑈 ˜ † , which ultimately yields the Lindblad master equation. Steps of this derivation, which follow those in [44] and involve projection onto the eigenbasis of 𝐻𝑑𝑜𝑡 , are detailed in Appendix C. 46 The result is: 𝑖 ∑︁ Γ𝑖 𝑗 𝜌¤ = − [𝐻𝑑𝑜𝑡 + 𝐻𝑠ℎ𝑖 𝑓 𝑡 , 𝜌] − ({𝜎𝑖+ 𝜎 𝑗− , 𝜌} − 2𝜎 𝑗− 𝜌) ℏ 𝑖𝑗 2 (4.7) 𝑖 ∑︁ Γ𝑖 𝑗 + − = − [𝐻𝑑𝑜𝑡 + 𝐻𝑠ℎ𝑖 𝑓 𝑡 , 𝜌] − (𝜎𝑖 𝜎 𝑗 𝜌 − 2𝜎 𝑗− 𝜌𝜎𝑖+ + 𝜌𝜎𝑖+ 𝜎 𝑗− ) ℏ 𝑖𝑗 2 where:    3Γ0 sin 𝑘 𝑅 cos 𝑘 𝑅 sin 𝑘 𝑅 Γ𝑖 𝑗 (r) = 𝜖ˆ𝑖 · (1 − r̂ ⊗ r̂) + (1 − 3r̂ ⊗ r̂) − · 𝜖ˆ 𝑗 (4.8) 4 𝑘𝑅 (𝑘 𝑅) 2 (𝑘 𝑅) 3    3Γ0 cos 𝑘 𝑅 sin 𝑘 𝑅 cos 𝑘 𝑅 Ω𝑖 𝑗 (r) = 𝜖ˆ𝑖 · −(1 − r̂ ⊗ r̂) + (1 − 3r̂ ⊗ r̂) + · 𝜖ˆ 𝑗 (4.9) 4 𝑘𝑅 (𝑘 𝑅) 2 (𝑘 𝑅) 3 ∑︁ 𝐻𝑠ℎ𝑖 𝑓 𝑡 = Ω𝑖 𝑗 𝜎𝑖+ 𝜎 𝑗− (4.10) 𝑖𝑗 and the spontaneous decay rate Γ0 is defined as in (2.13) or (3.38). Expression (4.8) describes the pairwise rate of emission induced by the interaction between the 𝑖th and 𝑗th emitters, and reproduces (3.37). Similarly, (4.9) describes a pairwise Lamb shift, and reproduces (3.41); furthermore, it is now apparent that the Lamb shift contributes an additional term (4.10) to the Hamiltonian 𝐻𝑑𝑜𝑡 . Equivalently, Ω and Γ can be derived from the real and imaginary parts of the Rabi energy of a classical radiated field (cf. 2.10) due to a polarization density P(r, 𝑡). Explicitly we may write: ˜ P̃(r 𝑗 )}/ℏ = −(Ω𝑖 𝑗 + 𝑖𝛾𝑖 𝑗 ) 𝜒 𝑗 (r𝑖 ) = d · 𝔉{ (4.11) where 𝜒 𝑗 (r𝑖 ) labels the Rabi energy of the 𝑖th dot due to the polarization of the 𝑗th dot. In this way, we establish the equivalence of the Master equation model and semiclassical model of Chapter 2. In Section 3.4, it was shown that the pairwise emission terms Γ𝑖 𝑗 give rise to superradiant emissions that scale quadratically as the number of emitters. As Γ𝑖 𝑗 plays an identical role in the Master equation, that result is equally valid here: ∑︁ ∑︁ 𝛾(𝑡) = Γ0 ⟨𝜎𝑖+ 𝜎𝑖− ⟩ + 2 Γ𝑖 𝑗 ⟨𝜎𝑖+ 𝜎 𝑗− ⟩(𝑡) 𝑖 𝑖≠ 𝑗 (4.12) = 𝛾0 (𝑡) + 𝛾 𝐼 (𝑡) As previously shown, the term 𝛾 𝐼 (𝑡), which describes emission from off-diagonal interactions among dots, corresponds to superradiant emissions, whereas 𝛾0 measures independent emissions. 47 4.2 Comparison of formulations We now describe results of simulations of quantum dot ensembles that employ both Maxwell- Bloch and Master equation approaches. As a differential equation for a 2𝑁 × 2𝑁 density matrix, the Master equation incurs an exponential computational cost and prohibits the simulation of large (𝑁 > 100) numbers of emitters, and for extended periods of time. Nonetheless, this global density matrix contains much more information on the quantum state of the system than what is typically measured in experiments. As such, we focus on transient dynamics of smaller systems of emitters. Using expression (4.12), we compute the off-diagonal contribution to the photon emission rate 𝛾 𝐼 (𝑡) under both approaches. In the Master equation approach, the time-dependent expectation values ⟨𝜎𝑖+ 𝜎 𝑗− ⟩ are calculated as Tr[𝜌(𝑡) 𝜎𝑖+ 𝜎 𝑗− ] after solving for 𝜌(𝑡) using (4.7). In the Maxwell- Bloch approach, ⟨𝜎𝑖+ 𝜎 𝑗− ⟩ are calculated as in (3.36): ⟨𝜎𝑖+ 𝜎 𝑗− ⟩ = ⟨Ψ0 |𝜎𝑖+ (𝑡)𝜎 𝑗− (𝑡)|Ψ0 ⟩ (4.13) where |Ψ0 ⟩ is the state of the system at 𝑡 = 0 and the 𝜎𝑖+ (𝑡) are operators in the Heisenberg representation. This can be approximated in the semiclassical model as ⟨Ψ0 |𝜎𝑖+ (𝑡)𝜎𝑖− (𝑡)|Ψ0 ⟩ ∼ 𝑖 (𝑡) and ⟨0|𝜎 + (𝑡)𝜎 − (𝑡)|0⟩ ∼ 𝜌𝑖 (𝑡) 𝜌 (𝑡), where 𝜌𝑖 (𝑡) is the two-level density matrix for the 𝑖th 𝜌11 𝑗 𝑖 𝑗 10 01 dot. Fig.4.1 show a comparison of 𝛾 𝐼 (𝑡) calculated for an initial state of the system in which each dot √ is in the (|0⟩ + |1⟩)/ 2 state, corresponding to the maximum initial polarization. Eight dots are in a chain along the 𝑦-axis, equally separated by distances given by 𝑠/𝜆, in a configuration studied in Ref. [45] with the Master equation method. There is a tight agreement between the approaches, and slight deviations are visible only in the 𝑑/𝜆 = 0.1 case. Also evident is the destructive interference that occurs when the dots are separated by a half-wavelength. As expected, for dots separated by 𝜆 the off-diagonal term becomes negligible in both cases. A special treatment has to be made if every dot is initially in the |1⟩ state. In this case, the initial polarization of the system is zero and remains zero at all times. This initial condition is equivalent to setting the Bloch vector of each dot in the equilibrium “up" initial position, which is unstable. The semiclassical approximation would then give no contribution to 𝛾 𝐼 (𝑡) in (4.12). This limitation 48 6 Emission rate 𝛾 𝐼 (𝑡) (units of Γ0 ) 4 2 0 0 2 4 6 8 10 12 14 16 Time ( 𝑝𝑠) Figure 4.1 Comparison of the off-diagonal term of the emission rate, 𝛾 𝐼 (𝑡) using the semiclassical Maxwell Bloch (solid) and the Master equation√(dashed). Eight quantum dots are initially prepared in the maximum polarization state (|0⟩ + |1⟩)/ 2. Curves for different dot separations 𝑠 compared to 𝜆 are shown. has been addressed by Haake et al. [46] who showed that semiclassical Maxwell-Bloch equations can describe superradiant pulses if Gaussian-distributed zero-averaging random initial conditions for the polarization are used. The random initial conditions provide the necessary tipping angle leading to spontaneous emission, and can be seen as the effect of a polarization measurement on the fully inverted state, which gives noise since the population and polarization operators do not commute. Fig. 4.2 shows a comparison of 𝛾 𝐼 (𝑡) for an initial state of the system in which each dot is in the |1⟩ state, corresponding to the maximum initial population inversion. We added a complex √︁ zero-averaging polarization following a Gaussian distribution with 𝜎 = 1/ 𝑁/𝑁𝑔 according to the method of Ref. [46]. Here 𝑁𝑔 is the number of groups into which the dots are divided, each receiving a different random polarization. In general the profile of the emission, consisting of a sharp rise post-excitation followed by decay to a steady low-emission state, resembles that of the Master equation. The methodology of Ref. [46] becomes exact in the limit 𝑁/𝑁𝑔 ≫ 1 and 𝑁𝑔 ≫ 1, therefore the discrepancies observed are due to the fact that we are considering a small system with 𝑁 = 8. Small values of 𝑁𝑔 lead to spurious non-zero values for 𝛾 𝐼 (𝑡) at 𝑡 = 0. However, we observe 49 Emission rate 𝛾 𝐼 (𝑡) (units of Γ0 ) 3 2 1 0 0 5 10 15 20 25 30 35 40 Time ( 𝑝𝑠) Figure 4.2 Comparison of 𝛾 𝐼 (𝑡) between the semiclassical Maxwell Bloch and the Master equation. Eight quantum dots equally spaced by 0.1𝜆 are considered. Each dot is initially in the maximum population state with the addition of a random polarization (the average has been taken over 500 random initial conditions). Additionally, the dots are divided into 𝑁𝑔 groups each receiving a different random polarization phase. that the averaged initial emission rate correctly tends to zero as the number of groups is increased. Using the same zero-averaging of initial conditions, emission rates were finally calculated for the Gaussian distributed dot configurations of Section III. In Fig. 4.3, superradiant emission is evidenced by both the characteristic rise and decay profile, as well as the 𝑁 2 scaling, of the emission rate. We also observed high variability of the emission for different random configurations in the Gaussian cloud. Of course, a comparison to the Master equation here is infeasible due to the exponential computation cost of that approach. 4.2.1 Time retardation effects Whereas the Maxwell-Bloch equations are a set of delayed differential equations that express their solution in terms of the polarization at previous times, the solution to the Master equation is essentially local in time. This feature, which leads to causality violation on small scales, is ultimately due to the Markov approximation. A simple simulation demonstrates that the Markov approximation amounts to disregarding delays in the Maxwell-Bloch equations. Reconsider a chain of dots, but where only the first dot is initially excited with the addition of random polarization. In the absence of delays, we would expect 50 12 Emission rate 𝛾 𝐼 (𝑡) (units of Γ0 ) 10 8 6 4 2 0 0 5 10 15 20 25 30 35 40 Time ( 𝑝𝑠) Figure 4.3 The second term of the emission rate, 𝛾 𝐼 (𝑡), via the semiclassical Maxwell Bloch equations is shown, for the Gaussian distributed dot configurations of Sect. III and all dots starting in the excited state (no pulse, 𝑁𝑔 = 1). Solid curves are calculated directly from (4.12), while dashed curves are calculated by subtracting independent emissions from the total emission as in (3.43). The average has been taken over 500 random initial conditions. instantaneous excitation of the other dots by re-radiation from the first dot, a feature which violates causality. Figure 4.4 shows the results of this simulation. The Markov approximated solution is in strong agreement with the Maxwell-Bloch solution without delay: we observe excitation at 𝑡 = 0 as this dot immediately receives a signal from the first dot. In contrast, the delayed Maxwell-Bloch solution is briefly unexcited, as the signal propagates down the chain, exciting each dot in an inverted “domino" fashion before arriving at the final dot. However, as this solution is convergent to that of the Master equation, the effect of the delay is negligible for simulation times much longer than the propagation time of the ensemble. Likewise, the Markov approximation becomes valid in that regime. Simulations of superradiance, for instance, which occur on the time scale 𝑇1 of the spontaneous decay, and involve tightly packed ensembles of dots, generally fall under this regime. 51 10−8 10−9 10−10 ⟨𝜌01 ⟩ 10−11 10−12 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 ·10−3 10−9 10−10 ⟨𝜌01 ⟩ 10−11 10−12 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time (units of 𝑇1 ) ·10−2 Figure 4.4 Comparison of short time behavior of the population excitation ⟨𝜌01 ⟩ using the Master equation and Maxwell-Bloch equations with and without delays. 𝑁 = 2 (top) and 𝑁 = 8 (bottom) quantum dots are separated by 1.0𝜆, with the first dot initialized in the excited state with a random initial polarization, and the rest lying initially in the ground state. Shown is the population inversion ⟨𝜌01 ⟩ for the last dot in the chain, i.e. the one that would receive a retarded signal the latest. The delay in excitation corresponds to the time 𝑅/𝑐 for the signal to propagate to the last dot (for a total separation 𝑅 = 1.0𝜆, 𝑅/𝑐 ≈ 0.157 𝑇1 ). 52 4.3 Concluding remarks In simulations of quantum dot ensembles, we have shown that the Maxwell-Bloch and Master equation approaches yield similar results and solid quantitative agreement. An inherent limitation of the Maxwell-Bloch model is the existence of an unstable equilibrium solution, corresponding to a state with perfect population inversion. This was overcome by applying zero-averaging of initial random polarizations. Indeed, when the initial polarization is non-zero, or when an exciting pulse induces a non-zero polarization, the Maxwell-Bloch solution never converges to this unstable equilibrium. In contrast, the Master equation is limited by the Markov approximation, which was shown to be equivalent to assuming the absence of delays in the Maxwell-Bloch model. However, this effect is negligible for simulations lasting much longer than the time for an electromagnetic signal to propagate through the medium. Additionally, we invoked the Born approximation, and assumed that the field is in “thermal" equilibrium and is weakly coupled to the ensemble. This is a reasonable assumption after strong excitations (such as a pulsed laser) have passed and the emitters are allowed to equilibrate with their radiated fields. 53 CHAPTER 5 CONCLUSION 5.1 Summary We began by outlining a theoretical framework for analyzing the interaction of semiconductor quantum dots with their radiated fields. We expressed such fields as classically propagating waves, and the dots as point polarization sources which generate these waves. The fields are coupled to density matrices representing each dot as a quantum two-level system, which evolve according to the Liouville equation. Numerically integrating this equation involved discretization in terms of spatial and temporal basis functions, and application of a Predictor-Corrector algorithm. Our computational model afforded us several advantages: (a) it permits fields to be resolved within near vicinities of each other, (b) it is scalable as 𝑂 (𝑁 log 𝑁) to large ensembles of dots with the assistance of the Adaptive Integral Method (AIM), which effects FFT-based convolutions to evaluate the field. This enabled large scale (𝑁 > 104 ) simulations that revealed both short range oscillatory behavior resulting from confinement of the field, and long range effects such as polarization dispersion. Having demonstrated the fidelity of our computational model, we focused in Chapter 3 on simulations of superradiance and subradiance. To enrich the analysis, we described the theory of quantized electric fields as a generalization of Dicke’s model. From this, we derived the Maxwell- Bloch equations and in turn, an expression for the photon emission rate. We showed that this consists of two terms, representing independent and coupled emissions respectively, and concluded that the latter corresponds to superradiance emissions, which scale quadratically with number of emitters. In simulations of medium-sized ensembles of quantum dots in a Gaussian cloud, we confirmed this quadratic scaling by subtracting independent emissions from total emissions. We also observed anisotropy of emission–another hallmark of superradiance–in the field radiated by the Gaussian cloud. Subradiance was discovered in steady state plots of the population excitation, which displayed diminished emissions. This effect was amplified by inhomogeneous broadening, which induced greater disorder and thus interference within the ensemble, but diminished by the presence of collective Lamb shifts. 54 In Chapter 4, we sought to validate our theoretical model by running in tandem simulations using an alternative formulation, the Master equation. Here we recovered the expression for the pairwise and total emission rates–first derived from the Maxwell-Bloch equations in Chapter 3. By applying zero-averaging random initial conditions to the polarization, we achieved strong numerical agreement between the two approaches. We observed both superradiant scaling, and destructive interference among dots separated by half-wavelengths. Lastly, we compared merits and drawbacks of both approaches. In particular, we showed that the Markov approximation, implicit in the Master equation, is equivalent to neglecting delays due to finite propagation speeds. Hence, simulations involving ensembles of emitters separated far apart in space should opt for the Maxwell-Bloch approach to accurately account for delay effects. 5.2 Future work 5.2.1 Multi-level systems In outlining our semiclassical model, we explicitly restricted ourselves to two-level systems, as this is sufficient to accurately describe a range of phenomena–such as superradiance–while simplifying calculations. For other scenarios, however, it may be convenient–or even necessary–to consider states with more than two energy levels. For instance, a three-level state can be represented by the density matrix: © 𝜌00 𝜌01 𝜌02 ª ­ ® ­ ∗ 𝜌ˆ ≡ ­ 𝜌01 (5.1) ® 𝜌11 𝜌12 ® ­ ® ­ ® ∗ 𝜌02 ∗ 𝜌12 𝜌22 « ¬ The Hamiltonian (2.3) and dissipation matrix (2.4) in the Maxwell-Bloch model generalizes naturally to: © 0 −ℏ𝜒01 (𝑡) −ℏ𝜒02 (𝑡) ª ­ ® Ĥ (𝑡) ≡ ℏ ­−ℏ𝜒01 ­ ∗ (𝑡) 𝜔1 ® −ℏ𝜒12 (𝑡) ® (5.2) ­ ® ­ ® −ℏ𝜒02∗ (𝑡) −ℏ𝜒 ∗ (𝑡) 𝜔 12 2 « ¬ 55 and: 𝜌 𝛾 𝜌 (𝛾 +𝛾 ) © 𝜌11 𝛾1 + 𝜌22 𝛾2 − 012 1 − 02 22 3 ª ­ ∗ 𝛾 ® −𝜌01 1 −𝜌 (𝛾 +𝛾 +𝛾 ) D̂ [ 𝜌] ˆ ≡­ 12 1 2 3 (5.3) ­ ® 2 −𝜌11 𝛾1 + 𝜌22 𝛾3 2 ® ­ ® ­ −𝜌∗ (𝛾2 +𝛾3 ) ∗ −𝜌12 (𝛾1 +𝛾2 +𝛾3 ) ® 02 2 2 −𝜌22 𝛾2 − 𝜌22 𝛾3 « ¬ where 𝜒𝑖 𝑗 labels the Rabi energy associated with transitions between the 𝑖th and 𝑗th levels, and 𝛾𝑖 is the spontaneous decay rate of the 𝑖th level. The generalization to three-levels can be applied to the study of quantum information. For instance, it can be used to model higher harmonic generation and electromagnetic induced transparency (EIT), which can be applied to the design of quantum devices. 5.2.2 Collective effects in waveguides and cavities Our prior investigations considered an ensemble of quantum dots embedded in vacuum or some homogeneous dielectric medium. In such free space, the re-emitted radiation can be represented as an infinite summation of eigenmodes with distinct wavevectors, propagating out to infinity. This enabled us to observe a variety of optical effects in a controlled, idealistic setting, while disregarding other effects that unnecessarily burden the analysis. A more realistic scenario has the ensemble contained within some waveguide or other closed conductive object. There, boundary conditions impose restrictions on the possible eigenmodes, and the form of the fields is highly dependent on the geometry of the system. The analysis therefore becomes far more complex, often requiring numerical Green’s functions [47] or other numerical techniques to simulate such fields. Nonetheless, the representation of the emitters remains intact. Within the Maxwell-Bloch formalism, the Hamiltonian for each emitter can be written in the form (2.3) or (5.2). Alternatively, it can be written in the context of the Master equation for 𝑁 emitters: ∑︁𝑁 ∑︁ Ĥ (𝑡) ≡ 𝜔𝑖 𝜎𝑖+ 𝜎𝑖− − 𝜎𝑖+ 𝜎 𝑗− 𝜒𝑖 𝑗 (𝑡) (5.4) 𝑖 𝑖≠ 𝑗 where 𝜎𝑖± are the usual raising and lowering operators for the 𝑖th emitter. In both cases, ℏ𝜒𝑖 𝑗 = d𝑖 · G(r𝑖 , r 𝑗 ) · d 𝑗 describes the Rabi energy of the 𝑖th emitter due to radiation from the 𝑗 emitter. As mentioned, the form of the Green’s function G is constrained by the underlying boundary conditions, and admits an analytical solution only for the simplest geometries. For instance, for a 56 perfectly cylindrical waveguide along 𝑧 with a single guided mode with longitudinal wavenumber 𝑘 𝑧 , the Green’s function takes the form [48]:  ′  𝐸 1 (𝜌)𝐸 2 (𝜌′) 𝑒𝑖𝑘 𝑧 (𝑧−𝑧 ) , 𝑧 > 𝑧′     ′ 𝐺 (r, r ) = 𝑖𝑔0 × (5.5)  ′ −𝑖𝑘 (𝑧−𝑧 ′) ′  𝐸 1 (𝜌)𝐸 2 (𝜌 ) 𝑒   𝑧 , 𝑧<𝑧  where 𝜌 = (𝑥, 𝑦) labels transverse coordinates, and 𝑔0 is a constant [49]. Such a mode propagates if |𝑘 𝑧 | is smaller than 𝜔/𝑐 where 𝑐 = 𝑐 0 /𝑛 is the speed of light in the material; otherwise it is attenuated transverse to the direction of the waveguide. In that case, emissions are completely suppressed for an infinite array, or are strongly subradiant for finite arrays [49]. This simple example shows that confinement of propagating waves can strongly affect the presence and degree of coupled emissions. Recent experiments have also observed confinement of atoms in waveguides due to their dipole radiation [50]. Such developments motivate the numerical simulation of optical effects in waveguides. 57 BIBLIOGRAPHY [1] Connor Glosser. The Quest for Active Media Models: A Self-Consistent Framework for Simulating Wave Propagation in Nonlinear Systems. PhD thesis, Michigan State University, 2018. [2] Connor Glosser, B Shanker, and Carlo Piermarocchi. Collective rabi dynamics of electromagnetically coupled quantum-dot ensembles. Physical Review A, 96(3):033816, 2017. [3] T. H. Stievater, Xiaoqin Li, D. G. Steel, D. Gammon, D. S. Katzer, D. Park, C. Piermarocchi, and L. J. Sham. Rabi oscillations of excitons in single quantum dots. Phys. Rev. Lett., 87:133603, Sep 2001. [4] H. Kamada, H. Gotoh, J. Temmyo, T. Takagahara, and H. Ando. Exciton rabi oscillation in a single quantum dot. Phys. Rev. Lett., 87:246401, Nov 2001. [5] H. Htoon, T. Takagahara, D. Kulik, O. Baklenov, A. L. Holmes, and C. K. Shih. Interplay of rabi oscillations and quantum interference in semiconductor quantum dots. Phys. Rev. Lett., 88:087401, Feb 2002. [6] Kenta Asakura, Yasuyoshi Mitsumori, Hideo Kosaka, Keiichi Edamatsu, Kouichi Akahane, Naokatsu Yamamoto, Masahide Sasaki, and Naoki Ohtani. Excitonic rabi oscillations in self-assembled quantum dots in the presence of a local field effect. Phys. Rev. B, 87:241301, Jun 2013. [7] G. Ya. Slepyan, S. A. Maksimenko, A. Hoffmann, and D. Bimberg. Quantum optics of a quantum dot: Local-field effects. Phys. Rev. A, 66:063804, Dec 2002. [8] G. Ya. Slepyan, A. Magyarov, S. A. Maksimenko, A. Hoffmann, and D. Bimberg. Rabi oscillations in a semiconductor quantum dot: Influence of local fields. Phys. Rev. B, 70:045320, Jul 2004. [9] M. Gross and S. Haroche. Superradiance: An essay on the theory of collective spontaneous emission. Physics Reports, 93(5):301 – 396, 1982. [10] Nicholas E. Rehler and Joseph H. Eberly. Superradiance. Phys. Rev. A, 3:1735–1751, May 1971. [11] J. C. MacGillivray and M. S. Feld. Theory of superradiance in an extended, optically thick medium. Phys. Rev. A, 14:1169–1189, Sep 1976. [12] C. R. Stroud, J. H. Eberly, W. L. Lama, and L. Mandel. Superradiant effects in systems of two-level atoms. Phys. Rev. A, 5:1094–1104, Mar 1972. 58 [13] C. Vanneste and P. Sebbah. Selective excitation of localized modes in active random media. Phys. Rev. Lett., 87:183903, Oct 2001. [14] A Fratalocchi, C Conti, and G Ruocco. Mode competitions and dynamical frequency pulling in Mie nanolasers: 3d ab-initio Maxwell-Bloch computations. Optics express, 16(12):8342– 8349, 2008. [15] R. H. Dicke. Coherence in spontaneous radiation processes. Phys. Rev., 93:99–110, Jan 1954. [16] G. S. Agarwal. Master equations in phase-space formulation of quantum optics. Phys. Rev., 178:2025–2035, Feb 1969. [17] R. H. Lehmberg. Radiation from an 𝑛-atom system. i. general formalism. Phys. Rev. A, 2:883–888, Sep 1970. [18] F. T. Arecchi and Eric Courtens. Cooperative phenomena in resonant electromagnetic propagation. Phys. Rev. A, 2:1730–1737, Nov 1970. [19] G. L. Lamb. Analytical descriptions of ultrashort optical pulse propagation in a resonant medium. Rev. Mod. Phys., 43:99–124, Apr 1971. [20] C Leonardi, J S Peng, and A Vaglica. Beats in dicke superradiant emission. Journal of Physics B: Atomic and Molecular Physics, 15(21):4017, nov 1982. [21] M. Gross, C. Fabre, P. Pillet, and S. Haroche. Observation of near-infrared dicke superradiance on cascading transitions in atomic sodium. Phys. Rev. Lett., 36:1035–1038, Apr 1976. [22] Ph. Cahuzac, H. Sontag, and P.E. Toschek. Visible superfluorescence from atomic europium. Optics Communications, 31(1):37–41, 1979. [23] N. Skribanowitz, I. P. Herman, J. C. MacGillivray, and M. S. Feld. Observation of dicke superradiance in optically pumped hf gas. Phys. Rev. Lett., 30:309–312, Feb 1973. [24] Gabriele Rainò, Michael A Becker, Maryna I Bodnarchuk, Rainer F Mahrt, Maksym V Kovalenko, and Thilo Stöferle. Superfluorescence from lead halide perovskite quantum dot superlattices. Nature, 563(7733):671–675, 2018. [25] Shoichi Okaba, Deshui Yu, Luca Vincetti, Fetah Benabid, and Hidetoshi Katori. Superradiance from lattice-confined atoms inside hollow core fibre. Communications Physics, 2:1–10, 2019. [26] L. Giannessi, M. Artioli, M. Bellaveglia, F. Briquez, E. Chiadroni, A. Cianchi, M. E. Couprie, G. Dattoli, E. Di Palma, G. Di Pirro, M. Ferrario, D. Filippetto, F. Frassetto, G. Gatti, M. Labat, G. Marcus, A. Mostacci, A. Petralia, V. Petrillo, L. Poletto, M. Quattromini, J. V. Rau, J. Rosenzweig, E. Sabia, M. Serluca, I. Spassovsky, and V. Surrenti. High-order-harmonic 59 generation and superradiance in a seeded free-electron laser. Phys. Rev. Lett., 108:164801, Apr 2012. [27] Michael Scheibner, Thomas Schmidt, Lukas Worschech, Alfred Forchel, Gerd Bacher, Thorsten Passow, and Detlef Hommel. Superradiance of quantum dots. Nature Physics, 3(2):106–110, Feb 2007. [28] Hideaki Taniyama, Hisashi Sumikura, and Masaya Notomi. Simulation technique of quantum optical emission process from multiple two-level atoms based on classical numerical method. Opt. Express, 27(9):12070–12079, Apr 2019. [29] Elliot Lu, Connor Glosser, Carlo Piermarocchi, and B. Shanker. QuEST:doi.org/10.5281/zenodo.6499556, May 2022. [30] Elliot Lu, Carlo Piermarocchi, and B Shanker. Modeling Radiation Reaction Induced Superradiance in Quantum Dot Systems. In 2020 IEEE International Symposium on Antennas and Propagation and North American Radio Science Meeting, pages 1653–1654. IEEE, 2020. [31] Elliot Lu, Connor Glosser, Carlo Piermarocchi, and B. Shanker. Numerical simulations of laser pulse propagation in quantum active media: Using a semiclassical model. IEEE Antennas and Propagation Magazine, 64(5):8–15, 2022. [32] C.R. Stroud and E.T. Jaynes. Long-term solutions in semiclassical radiation theory. Phys. Rev. A, 1:106–121, Jan 1970. [33] Balasubramaniam Shanker, Mingyu Lu, Jun Yuan, and Eric Michielssen. Time domain integral equation analysis of scattering from composite bodies via exact evaluation of radiation fields. Antennas and Propagation, IEEE Transactions on, 57:1506 – 1520, 06 2009. [34] C. Glosser, E. Lu, T.J. Bertus, C. Piermarocchi, and B. Shanker. Acceleration techniques for semiclassical maxwell–bloch systems: An application to discrete quantum dot ensembles. Computer Physics Communications, 258:107500, 2021. [35] Andreas Glaser and Vladimir Rokhlin. A new class of highly accurate solvers for ordinary differential equations. Journal of Scientific Computing, 38(3):368–399, 2009. [36] R Friedberg, SR Hartmann, and JT Manassah. Limited superradiant damping of small samples. Physics Letters A, 40(5):365–366, 1972. [37] Ana Asenjo-Garcia, M Moreno-Cardoner, Andreas Albrecht, HJ Kimble, and Darrick E Chang. Exponential improvement in photon storage fidelities using subradiance and “selective radiance” in atomic arrays. Physical Review X, 7(3):031024, 2017. [38] A. Piñeiro Orioli, J. K. Thompson, and A. M. Rey. Emergent dark states from superradiant dynamics in multilevel atoms in a cavity. Phys. Rev. X, 12:011054, Mar 2022. 60 [39] Anatoly A. Svidzinsky, Xiwen Zhang, and Marlan O. Scully. Quantum versus semiclassical description of light interaction with atomic ensembles: Revision of the maxwell-bloch equations and single-photon superradiance. Phys. Rev. A, 92:013801, Jul 2015. [40] Han Cai. Superradiance and Topological Quantum Optics in Atomic Medium and Cavity QED. PhD thesis, Texas A&M University, 2017. [41] Elliot Lu, B. Shanker, and Carlo Piermarocchi. Transient dynamics of subradiance and superradiance in open optical ensembles. Phys. Rev. A, 107:043703, Apr 2023. [42] Kensuke Miyajima, Y Kagotani, S Saito, M Ashida, and T Itoh. Superfluorescent pulsed emission from biexcitons in an ensemble of semiconductor quantum dots. Journal of Physics: Condensed Matter, 21(19):195802, 2009. [43] Girish S Agarwal. Quantum optics. Cambridge University Press, 2012. [44] H.-P. Breuer and F. Petruccione. The Theory of Open Quantum Systems. Oxford University Press, 2002. [45] Stuart J. Masson and Ana Asenjo-Garcia. Universality of Dicke superradiance in arrays of quantum emitters. Nature Communications, 13(1):2285, April 2022. [46] Fritz Haake, Harald King, Guntram Schröder, Joe Haus, and Roy Glauber. Fluctuations in superfluorescence. Physical Review A, 20(5):2047, 1979. [47] Chen-To Tai. Dyadic Green Functions in Electromagnetic Theory 2nd ed. IEEE Press, 1993. [48] Sina Saravi, Alexander N. Poddubny, Thomas Pertsch, Frank Setzpfandt, and Andrey A. Sukhorukov. Atom-mediated spontaneous parametric down-conversion using bandgap modes in nonlinear periodic waveguides. In Advanced Photonics 2018 (BGPP, IPR, NP, NOMA, Sensors, Networks, SPPCom, SOF), page NpTh3I.5. Optica Publishing Group, 2018. [49] Alexandra S. Sheremet, Mihail I. Petrov, Ivan V. Iorsh, Alexander V. Poshakinskiy, and Alexander N. Poddubny. Waveguide quantum electrodynamics: Collective radiance and photon-photon correlations. Rev. Mod. Phys., 95:015002, Mar 2023. [50] A. Goban, C.-L. Hung, S.-P. Yu, J. D. Hood, J. A. Muniz, J. H. Lee, M. J. Martin, A. C. McClung, K. S. Choi, D. E. Chang, O. Painter, and H. J. Kimble. Atom-light interactions in photonic crystals. Nature Communications, 5(1):3808, May 2014. [51] L.D. Landau and E. M. Lifschitz. The Classical Theory of Fields 3rd ed. Pergamon Press Ltd., 1971. [52] Carlos M. Bustamante, Esteban D. Gadea, Tchavdar N. Todorov, and Damián A. Scherlis. Tailoring cooperative emission in molecules: Superradiance and subradiance from first- 61 principles simulations. The Journal of Physical Chemistry Letters, 13(50):11601–11609, 2022. PMID: 36480910. 62 APPENDIX A INTEGRAL OPERATOR ELECTRIC FIELDS A.1 Fundamentals The following forms for Maxwell’s equations will be assumed:   ∇ · E = 𝜖𝜌0         ∇ · B = 0    (A.1) ∇ × E = − 𝜕𝜕𝑡B         1 𝜕E   ∇ × B = 𝜇0 J + 𝑐2 𝜕𝑡    Additionally, Fourier transforms are assumed to be unitary: ∫ ∫ 1 −𝑖𝜔𝑡 1 𝐹 (𝜔) = √ 𝑓 (𝑡) 𝑒 𝑑𝑡 ⇐⇒ 𝑓 (𝑡) = √ 𝐹 (𝜔) 𝑒𝑖𝜔𝑡 𝑑𝜔 (A.2) 2𝜋 2𝜋 A.2 Electric field due to polarization sources Starting from Maxwell’s equations (A.1), one may express a time-varying electric field in terms of scalar and vector potentials 𝜙 and A as: 𝜕A E(r, 𝑡) = −∇𝜙 − (A.3) 𝜕𝑡 Differentiating with respect to time, and applying the Lorenz gauge condition: 1 𝜕𝜙 ∇·A+ =0 (A.4) 𝑐2 𝜕𝑡 to eliminate 𝜙 yields: 𝜕2A 𝜕2   𝜕E 2 2 = −∇(−𝑐 ∇ · A) − 2 = 𝑐 ∇∇ − 2 A (A.5) 𝜕𝑡 𝜕𝑡 𝜕𝑡 where ∇∇ is a dyadic operator. The retarded potential due to a current source J(r, 𝑡) is given by: 𝛿(𝑡 𝑅 − 𝑡 ′) 𝛿(𝑡 𝑅 − 𝑡 ′) 𝜕P ′ ′ ∫ ∫ 𝜇0 𝜇0 A(r, 𝑡) = ′ J(r′, 𝑡 ′) 𝑑r′ 𝑑𝑡 ′ = 𝑑r 𝑑𝑡 (A.6) 4𝜋 |r − r | 4𝜋 |r − r′ | 𝜕𝑡 ′ 63 where 𝛿 is a delta function, 𝑡 𝑅 = 𝑡 − |r − r′ |/𝑐, P is a polarization density, related to the current density by J = 𝜕𝑡 P. This can be written in terms of the Green’s function for the wave operator   1 𝜕2 2 : 𝑐2 𝜕𝑡 2 − ∇ 𝛿(𝑡 − 𝑟/𝑐) 𝑔(r, 𝑡) = (A.7) 4𝜋𝑟 as: A(r, 𝑡) = 𝜇0 𝑔(r, 𝑡) ★𝑠𝑡 𝜕𝑡 P(r, 𝑡) (A.8) where ★𝑠𝑡 denotes both spatial and temporal convolution. Substituting this expression back into (A.5) and canceling one time derivative we obtain: E(r, 𝑡) ≡ −𝜇0 (𝜕𝑡2 − 𝑐2 ∇∇)𝑔(r, 𝑡) ★𝑠𝑡 P(r, 𝑡) (A.9) 𝛿(𝑡 𝑅 − 𝑡 ′)  2 ∫ 𝜕 = −𝜇0 2 − 𝑐 2 ∇∇ ′ P(r′, 𝑡 ′) 𝑑r′ 𝑑𝑡 ′ 𝜕𝑡 4𝜋|r − r | Alternatively, one may start by deriving eq. (3.26) from Maxwell’s equations as follows: 𝜕2P 1 𝜕2E   2 ∇·P 𝜕 ∇ × ∇ × E = ∇(∇ · E) − ∇ E = ∇ − − ∇2 E = − ∇ × B = −𝜇0 2 − 2 2 (A.10) 𝜖0 𝜕𝑡 𝜕𝑡 𝑐 𝜕𝑡 where again J = 𝜕𝑡 P and we assume only a bound charge density 𝜌 = −∇ · P exists. Rearranging yields (3.26), reproduced here: 1 𝜕2    2  2 𝜕 2 − ∇ E(r, 𝑡) = −𝜇0 − 𝑐 ∇∇ P(r, 𝑡) (A.11) 𝑐2 𝜕𝑡 2 𝜕𝑡 2 As the wave operator appears on the LHS and has associated Green’s function 𝑔(r, 𝑡), the result (A.9) then follows after observing that the operators involved commute due to their linearity. This approach has several merits. First, it passes directly from source to field, thus eschews gauge fixing; second, it draws upon a set of equations (A.11) which have a counterpart in the quantum Heisenberg picture, namely (3.25). The derivation of that set of equations is the subject of Appendix B. While we have obtained an expression for the field radiated by a polarization source P, direct evaluation of the convolution (A.9) is intractable, as it involves taking spatial derivatives of a delta function. To obtain an explicit analytic expression, we make a foray into Fourier space. For the problem of dipole radiation, we assume time-harmonic fields, i.e. E(r, 𝑡) = E (r)𝑒𝑖𝜔𝑡 and similarly 64 for the potential and sources. Substitution into (A.5) yields a Fourier space representation of the field:   2 ∇∇ 𝑖𝜔E(r, 𝜔) = 𝜔 I+ 2 A (A.12) 𝑘 The time-harmonic equivalent of (A.6) is: ∫ 𝑖𝑘 |r−r′ | ∫ 𝑖𝑘 |r−r′ | 𝜇0 𝑒 ′ ′ 𝑖𝜔𝜇0 𝑒 A(r, 𝜔) = ′ J(r ) 𝑑r = ′ P(r′) 𝑑r′ (A.13) 4𝜋 |r − r | 4𝜋 |r − r | Inserting back into (A.5) we obtain: ′ 𝑒𝑖𝑘 |r−r |  ∫ ∇∇ 2 E(r, 𝜔) = 𝜔 𝜇0 I + 2 P(r′, 𝜔) 𝑑r′ ≡ 𝜔2 𝜇0 𝐺 (r, 𝜔) ★ P(r, 𝜔) (A.14) 𝑘 4𝜋|r − r′ | where ★ denotes spatial convolution and:   ∇∇ 𝑒𝑖𝑘𝑟 𝐺 (r, 𝜔) = I + 2 (A.15) 𝑘 4𝜋𝑟 is the dyadic Green’s function for the Helmholtz equation. By explicitly expanding the Helmholtz Green’s function (A.15) in spherical coordinates we obtain (r̂ = r/𝑟):      𝑖𝑘𝑟 ∇∇ 𝑒𝑖𝑘r 1 1 𝑒 I+ 2 = (𝐼 − r̂ ⊗ r̂) − (𝐼 − 3r̂ ⊗ r̂) + 2 (A.16) 𝑘 4𝜋r 𝑖𝑘𝑟 (𝑘𝑟) 4𝜋𝑟 which facilitates evaluating the convolution in (A.14). The result is: ∫    𝑖𝑘 𝑅 𝜇0 1 1 𝑒 E(r, 𝜔) = − (𝐼 − r̂ ⊗ r̂) − (𝐼 − 3r̂ ⊗ r̂) + 2 · P(r′, 𝜔) 𝑑r′ (A.17) 4𝜋 𝑖𝑘 𝑅 (𝑘 𝑅) 𝑅 where R = r − r′ and r̂ = R/𝑅. Transforming back to the time domain by replacing each factor of 𝑖𝜔 by a time derivative, and 𝑒𝑖𝑘 𝑅 /𝑅 by 𝛿(𝑡 − 𝑅/𝑐)/𝑅 = 𝛿(𝑡 𝑅 )/𝑅 gives: 𝜕𝑡2 𝛿(𝑡 𝑅 − 𝑡 ′) 𝑐𝜕𝑡 𝛿(𝑡 𝑅 − 𝑡 ′) 𝑐2 𝛿(𝑡 𝑅 − 𝑡 ′) ∫    𝜇0 E(r, 𝑡) = − (𝐼 − r̂ ⊗ r̂) − (𝐼 − 3r̂ ⊗ r̂) + · 4𝜋 𝑅 𝑅2 𝑅3 P(r′, 𝑡 ′) 𝑑r′ 𝑑𝑡 ′ (A.18) Linearity permits us to transfer the derivatives acting on the 𝛿-functions to the polarization. The final result is (cf. (2.10)): ∫ " # 𝜕𝑡2 P(r′, 𝑡 𝑅 ) 𝜕𝑡 P(r′, 𝑡 𝑅 ) P(r′, 𝑡 𝑅 )  −1 E(r, 𝑡) = (I − r̂ ⊗ r̂) · + (I − 3r̂ ⊗ r̂) · + 𝑑 3 r′ (A.19) 4𝜋𝜖 𝑐2 𝑅 𝑐𝑅2 𝑅3 65 A.3 Radiation reaction field The singularity of (A.19) at 𝑅 = 0 can be resolved by Taylor expanding (A.14) and retaining specific terms. This ultimately yields the radiation reaction field. To start, let us separate out the two terms in parentheses from (A.14): ′ ′ 𝑒𝑖𝑘 |r−r | 𝜔2 𝜇0 𝑒𝑖𝑘 |r−r | ∫ ∫ ′ ′ E(r, 𝜔) = 𝜔 𝜇02 P(r , 𝜔) 𝑑r + ∇∇ P(r′, 𝜔) 𝑑r′ (A.20) 4𝜋|r − r′ | 𝑘2 4𝜋|r − r′ | and Taylor expand each in 𝑘𝑟 in the limit of small 𝑘𝑟. Expanding the first term to first order we obtain: 𝑖𝜔2 𝜇0 ∫ E 𝑅𝑅1 (r, 𝜔) ≈ 𝑘 P(r′, 𝜔) 𝑑r′ (A.21) 4𝜋 For the second term, the two spatial derivatives eliminate two powers of 𝑟, so we expand to third order: 𝜔2 𝜇0 ∇∇ 𝑖𝑘 3 |r − r′ | 2 ∫   E 𝑅𝑅2 (r, 𝜔) ≈ − P(r′, 𝜔) 𝑑r′ 4𝜋 𝑘2 3! (A.22) 𝑖𝜔2 𝜇0 𝑘 ∫ =− P(r′, 𝜔) 𝑑r′ 4𝜋 3 The total nearfield is thus: 𝑖𝜔2 𝜇0 2𝑘 ∫ E 𝑅𝑅 (r, 𝜔) = P(r′, 𝜔) 𝑑r′ (A.23) 4𝜋 3 Transforming back to time by replacing 𝑖𝜔 with time derivatives, we arrive at the radiation reaction field (cf. [51], note the minus sign convention): ∫ 1 2 E 𝑅𝑅 (r, 𝑡) ≡ 𝜕𝑡3 P(r′, 𝑡) 𝑑r′ (A.24) 4𝜋𝜖 3𝑐3 The rotating frame equivalent is: ∫ 1 2 𝜕𝑡3 P̃(r′, 𝑡) + 3𝑖𝜔 𝐿 𝜕𝑡2 P̃(r′, 𝑡) − 3𝜔2𝐿 𝜕𝑡 P̃(r′, 𝑡) − 𝑖𝜔3𝐿 P̃(r′, 𝑡) 𝑑r′ (A.25)   E 𝑅𝑅 (r, 𝑡) = − 4𝜋𝜖 3𝑐3 As these expressions were obtained by assuming the limit 𝑘𝑟 → 0, which is no longer valid as one moves away from the vicinity of an emitter, the radiation reaction field is felt exclusively by the emitter and those in its nearfield. The force due to the radiation reaction field is the well-known 66 Abraham-Lorentz force. This parallels a similar derivation in the context of Heisenberg’s equations [52]. In the derivation of (A.22) we have ignored the contributions of the zeroth through second order terms, which are obvious non-negligible in the limit of small 𝑘𝑟. Evaluating their contribution leads to an expression for the leftover nearfield: 𝜔2 𝜇0 I ∇∇ 1 + 𝑖𝑘𝑟 − (𝑘𝑟) 2 /2    E𝑁 𝐹−𝑅𝑅 (r, 𝜔) = + 2 ★ P(r, 𝜔) (A.26) 4𝜋 𝑟 𝑘 𝑟 Via the expansion: ∇∇ 1 + 𝑖𝑘𝑟 − 𝑘 2𝑟 2 /2       1 1 1 1 = 2 − 3 (I − 3r̂ ⊗ r̂) − (I − r̂ ⊗ r̂) (A.27) 𝑘2 𝑟 𝑘 𝑟 2 𝑟 one obtains, after substituting back into (A.26): P(r′, 𝜔) P(r′, 𝜔) ∫   1 E𝑁 𝐹−𝑅𝑅 (r, 𝜔) = − 2 (𝑖𝜔) (I + r̂ ⊗ r̂) + (I − 3r̂ ⊗ r̂) 𝑑r′ (A.28) 4𝜋𝜖 2𝑐2 𝑅 𝑅3 Replacing 𝑖𝜔 by time derivatives, restoring the temporal convolution, and applying the resulting 𝛿-functions to P(r′, 𝑡) gives: 𝜕𝑡2 P(r′, 𝑡) P(r′, 𝑡) ∫   1 E𝑁 𝐹−𝑅𝑅 (r, 𝑡) = − (I + r̂ ⊗ r̂) 2 + (I − 3r̂ ⊗ r̂) 3 𝑑r′ (A.29) 4𝜋𝜖 2𝑐 𝑅 𝑅 The total nearfield is then approximated as: E𝑁 𝐹 (r, 𝑡) = E 𝑅𝑅 (r, 𝑡) + E𝑁 𝐹−𝑅𝑅 (r, 𝑡) (A.30) Note that for both terms the retarded time 𝑡 𝑅 has been replaced by the present time 𝑡, which is a suitable approximation when the source and observer are sufficiently close. 67 APPENDIX B MAXWELL-BLOCH EQUATIONS B.1 Derivation We start by deriving the Maxwell-Bloch equations: 1 𝜕2   2 1 1 2 2 − ∇ E(r, 𝑡) = ∇ × ∇ × P(r, 𝑡) = − (∇2 − ∇∇)P(r, 𝑡) (B.1) 𝑐 𝜕𝑡 𝜖 𝜖 Recall from Section (3.2) the definitions of the quantized fields in the Heisenberg picture: ∑︁ E+ (r, 𝑡) = −𝑖 E 𝑘,n̂ 𝑎 k,n̂ (𝑡)𝑒𝑖k·r (B.2) k,n̂ ∑︁ E− (r, 𝑡) = 𝑖 E 𝑘,n̂ 𝑎 †k,n̂ (𝑡)𝑒 −𝑖k·r (B.3) k,n̂ √︃ ℏ𝑐𝑘 (here E 𝑘,n̂ = 𝐸 𝑘 n̂ = 2𝜖𝑉 n̂ is the vector amplitude for a mode with wavevector k and polarization n̂). For simplicity let us assume that the interaction is between the field and a single emitter located at r0 (which easily generalizes to the case of 𝑁 emitters). Hence the Hamiltonians (without assuming RWA) are: 𝐻𝑑𝑜𝑡 = ℏ𝜔 𝜎 𝑧 (B.4) ∑︁ 𝐻𝑟𝑎𝑑 = ℏ 𝜔 𝑘 (𝑎 †k,n̂ 𝑎 𝑘,n̂ + 1/2) (B.5) k,n̂ 𝐻𝑖𝑛𝑡 = −d · (E+ (r0 ) + E− (r0 )) ⊗ (𝜎 + + 𝜎 − ) (B.6) Heisenberg’s equation of motion then implies that (cf. 3.21): d𝑎 𝑘,n̂ 1 = −𝑖𝜔 𝑘 𝑎 𝑘,n̂ (𝑡) − d · E 𝑘,n̂ 𝑒 −𝑖k·r0 (𝜎 + (𝑡) + 𝜎 − (𝑡)) (B.7) d𝑡 ℏ Writing 𝑎 𝑘 (𝑡) = 𝑎˜ 𝑘 (𝑡)𝑒 −𝑖𝜔 𝑘 𝑡 eliminates the first term: d𝑎˜ k,n̂ 1 = − d · E 𝑘,n̂ 𝑒 −𝑖(k·r0 −𝑖𝜔 𝑘 𝑡) (𝜎 + (𝑡) + 𝜎 − (𝑡)) (B.8) d𝑡 ℏ 68 in terms of which (B.2) can be written: ∑︁ E+ (r, 𝑡) = −𝑖 E 𝑘,n̂ 𝑎˜ 𝑘,n̂ 𝑒𝑖(k·r−𝜔 𝑘 𝑡) (B.9) k,n̂ and similarly for E− (r, 𝑡).   1 𝜕2 Applying the wave operator 𝑐2 𝜕𝑡 2 − ∇2 to E+ (r, 𝑡) yields   2 1 𝜕2     2 + ∑︁ 1 d 𝑎˜ k,n̂ d 𝑎˜ k,n̂ 2 2 2 2 − ∇ E (r, 𝑡) = −𝑖 E 𝑘,n̂ 2 2 − 2𝑖𝜔 𝑘 − 𝜔 𝑘 𝑎˜ k,n̂ + 𝑘 𝑎˜ k,n̂ 𝑒𝑖(k·r−𝜔 𝑘 𝑡) 𝑐 𝜕𝑡 𝑐 d𝑡 d𝑡 k,n̂ d2 𝑎˜ k,n̂   𝑖 ∑︁ d 𝑎˜ k,n̂ 𝑖(k·r−𝜔 𝑘 𝑡) =− 2 E 𝑘,n̂ − 2𝑖𝜔 𝑘 𝑒 𝑐 k,n̂ d𝑡 2 d𝑡 (B.10) Substituting in (B.8) yields: 1 𝜕2     + 𝑖 ∑︁ d + 2 − ∇ E (r, 𝑡) = 2 (d · E 𝑘,n̂ )E 𝑘,n̂ 𝜎 + 𝜎 − 𝑖𝜔 𝑘 (𝜎 + 𝜎 ) 𝑒𝑖k·(r−r0 ) (B.11) − + − 𝑐2 𝜕𝑡 2 ℏ𝑐 k,n̂ d𝑡 Repeating the procedure for E− (r, 𝑡), and applying the wave operator to the total field E(r, 𝑡) = E+ (r, 𝑡) + E− (r, 𝑡) we have: 1 𝜕2   2 2 ∑︁ − ∇ E(r, 𝑡) = (d · E 𝑘,n̂ )E 𝑘,n̂ 𝑐2 𝜕𝑡 2 ℏ𝑐2 k,n̂   d + − + − × − 𝜎 + 𝜎 sin(k · (r − r0 )) + 𝜔 𝑘 (𝜎 + 𝜎 ) cos(k · (r − r0 )) d𝑡 2 ∑︁ = 2 (d · E 𝑘,n̂ )E 𝑘,n̂ 𝜔 𝑘 (𝜎 + + 𝜎 − ) cos(k · (r − r0 )) ℏ𝑐 k,n̂ (B.12) where the last step is justified since (d · E 𝑘,n̂ )E 𝑘,n̂ sin(k · (r − r0 )) is an odd function of k. Using √︃ ℏ𝑐𝑘 the definitions 𝜔 𝑘 = 𝑐𝑘, E 𝑘,n̂ = 𝐸 𝑘 n̂ = 2𝜖𝑉 n̂, this can be rewritten: 1 𝜕2   1 ∑︁ 2 2 − ∇ 2 E(r, 𝑡) = (dn̂) n̂ 𝑘 2 (𝜎 + + 𝜎 − ) cos(k · (r − r0 )) (B.13) 𝑐 𝜕𝑡 𝜖𝑉 k,n̂ Í 𝑉 ∫ Í Replacing k → (2𝜋) 3 𝑑k and using the completeness relation n̂ n̂n̂ = I − k̂ ⊗ k̂ (cf. (3.31)) we arrive at: 1 𝜕2   ∫  1 𝑉 2 − ∇ E(r, 𝑡) = d 𝑑k(𝑘 − k ⊗ k) cos(k · (r − r0 )) (𝜎 + + 𝜎 − ) 2 (B.14) 𝑐2 𝜕𝑡 2 𝜖 (2𝜋) 3 69 The term in brackets can be written: ∫ ∫ 2 2 𝑑k(𝑘 −k⊗k) cos(k · (r − r0 )) = (−∇ +∇∇) 𝑑k cos(k · (r − r0 )) = (2𝜋) 3 (−∇2 +∇∇)𝛿(r−r0 ) (B.15) Thus: 1 𝜕2   2 1 2 + − + −  − ∇ E(r, 𝑡) = −d ∇ 𝛿(r − r 0 )(𝜎 + 𝜎 ) + d ∇∇𝛿(r − r 0 )(𝜎 + 𝜎 ) 𝑐2 𝜕𝑡 2 𝜖 1 (B.16) = − (∇2 − ∇∇)(P+ (r, 𝑡) + P− (r, 𝑡)) 𝜖 1 2 = − (∇ − ∇∇)P(r, 𝑡) 𝜖 using the definition (3.22) for the polarization operators. B.2 Equivalence to classical equations The solution of (B.16) is: 1 2 1 2  ∫ 𝛿(𝑡 − 𝑡 ′) P(r′, 𝑡 ′) 𝑑r′ 𝑑𝑡 ′ (B.17) 𝑅 E(r, 𝑡) ≡ − (∇ − ∇∇)𝑔(r, 𝑡) ★𝑠𝑡 P(r, 𝑡) = − ∇ − ∇∇ ′ 𝜖 𝜖 4𝜋|r − r | which differs from the classical solution (A.9) in that the temporal Laplace operator acting on the polarization density is replaced by a spatial one, i.e. 1/𝑐2 𝜕𝑡2 P → ∇2 P. However, passage to the Fourier domain shows that the solution is unaffected by this discrepancy. Indeed, we may rewrite the Maxwell-Bloch solution in terms of the vector potential operator A as (cf. A.5): 𝜕E   = 𝑐2 ∇∇ − ∇2 A (B.18) 𝜕𝑡 which, assuming a time harmonic solution E(r, 𝑡) = E (r)𝑒𝑖𝜔𝑡 , has the Fourier space representation:   𝑖𝜔E(r, 𝜔) = 𝑐2 ∇∇ − ∇2 A (B.19) wherein E(r, 𝜔) can be written in terms of the dyadic Green’s function 𝐺 (r, 𝜔) as: E(r, 𝜔) ≡ 𝜔2 𝜇0 𝐺 (r, 𝜔) ★ P(r, 𝜔) (B.20) where: 1  2  𝑒𝑖𝑘𝑟 𝐺 (r, 𝜔) = 2 −∇ + ∇∇ (B.21) 𝑘 4𝜋𝑟 70 𝑒 𝑖𝑘𝑟 As 𝑓 (𝑟) = 4𝜋𝑟 is the solution to: (∇2 + 𝑘 2 ) 𝑓 (𝑟) = 𝛿(𝑟) (B.22) we may write (for 𝑟 ≠ 0): 1  2  𝑒𝑖𝑘𝑟   ∇∇ 𝑒𝑖𝑘𝑟 𝐺 (r, 𝜔) = 2 𝑘 + ∇∇ = I+ 2 (B.23) 𝑘 4𝜋𝑟 𝑘 4𝜋𝑟 which recovers (A.15), and completes the proof of the equivalence of the two expressions (A.9) and (B.17) representing a time-harmonic field E in terms of the polarization density P. 71 APPENDIX C DERIVATION OF LINDBLAD MASTER EQUATION In Section 4.1, we saw that invoking the Born-Markov approximation yields an integro-differential master equation for the density matrix in the interaction picture: ∫ ∞ d −1 ˜ = 2 Tr𝑟𝑎𝑑 𝜌(𝑡) 𝑑𝑡 ′ [𝐻𝑖𝑛𝑡 (𝑡), [𝐻𝑖𝑛𝑡 (𝑡 − 𝑡 ′), 𝜌(𝑡) ˜ ⊗ 𝜌˜𝑟𝑎𝑑 ]] (C.1) d𝑡 ℏ 0 where 𝐻𝑖𝑛𝑡 (𝑡) is the interaction Hamiltonian. Now observe that the interaction Hamiltonian in the Schrödinger picture: ∑︁ 𝐻𝑖𝑛𝑡 = − d𝑖 · (E+ (r𝑖 ) + E− (r𝑖 )) ⊗ (𝜎𝑖+ + 𝜎𝑖− ) 𝑖 ∑︁ ∑︁ (C.2) =− d̂𝑖 · E𝑖 = 𝐴𝛼 ⊗ 𝐵 𝛼 𝑖 𝛼 Í ˆ Í can be written as a product of operators 𝐴𝛼 = − 𝑖 ( 𝑑 𝛼 )𝑖 and 𝐵𝛼 = 𝑖 (𝐸 𝛼 )𝑖 acting on the dots and field, respectively (here 𝛼 indexes each of the three Cartesian components of d̂𝑖 or E𝑖 ). Let us now decompose the 𝐴𝛼 in terms of the eigenstates of 𝐻𝑑𝑜𝑡 , which are precisely the Dicke states |𝑆, 𝑀⟩: 𝐻𝑑𝑜𝑡 |𝑆, 𝑀⟩ = ℏ𝜔𝜎𝑧 |𝑆, 𝑀⟩ = 𝑀ℏ𝜔 |𝑆, 𝑀⟩ (C.3) It therefore suffices to project 𝐴𝛼 onto subspaces defined by an energy difference ℏ𝜔 (or equivalently, a difference in number of excited states 𝑀 ′ − 𝑀 = 1): ∑︁ 𝐴𝛼 (𝜔) = ⟨𝑆, 𝑀 | 𝐴𝛼 |𝑆, 𝑀 ′⟩ (C.4) 𝑀 ′ −𝑀=1 By the completeness of the eigenstates of 𝐻𝑑𝑜𝑡 we may write: ∑︁ ∑︁ 𝐴𝛼 = 𝐴𝛼 (𝜔) = 𝐴𝛼 (𝜔) † (C.5) 𝜔 𝜔 Thus: ∑︁ ∑︁ 𝐻𝑖𝑛𝑡 = 𝐴𝛼 (𝜔) ⊗ 𝐵𝛼 = 𝐴𝛼 (𝜔) † ⊗ 𝐵†𝛼 (C.6) 𝛼,𝜔 𝛼,𝜔 72 To transform 𝐻𝑖𝑛𝑡 to the interaction picture, recall that 𝑈 (𝑡) = exp(−𝑖(𝐻𝑑𝑜𝑡 + 𝐻𝑟𝑎𝑑 )𝑡/ℏ), so 𝑈 † (𝑡) 𝐴𝛼 (𝜔)𝑈 (𝑡) = 𝑒 −𝑖𝜔𝑡 𝐴𝛼 (𝜔), and: ∑︁ ∑︁ 𝐻𝑖𝑛𝑡 (𝑡) = 𝑒 −𝑖𝜔𝑡 𝐴𝛼 (𝜔) ⊗ 𝐵𝛼 (𝑡) = 𝑒𝑖𝜔𝑡 𝐴𝛼 (𝜔) † ⊗ 𝐵†𝛼 (𝑡) (C.7) 𝛼,𝜔 𝛼,𝜔 where 𝐵𝛼 (𝑡) = 𝑒𝑖𝐻𝑟 𝑎𝑑 𝑡 𝐵𝛼 𝑒 −𝑖𝐻𝑟 𝑎𝑑 𝑡 . We can finally substitute this expression into the Born-Markov master equation (C.1). After expanding out the commutators one obtains: ∫ ∞ d 1 ˜ = 2 Tr𝑟𝑎𝑑 𝜌(𝑡) 𝑑𝑡 ′ [𝐻𝑖𝑛𝑡 (𝑡 − 𝑡 ′) 𝜌(𝑡) ˜ 𝜌˜𝑟𝑎𝑑 𝐻𝑖𝑛𝑡 (𝑡) − 𝐻𝑖𝑛𝑡 (𝑡)𝐻𝑖𝑛𝑡 (𝑡 − 𝑡 ′) 𝜌(𝑡) ˜ 𝜌˜𝑟𝑎𝑑 ] + etc. d𝑡 ℏ 0 ∑︁ ∑︁ ′ = 𝑒 −𝑖(𝜔−𝜔 )𝑡 𝛾𝛼𝛽 (𝜔) [ 𝐴 𝛽 (𝜔) 𝜌(𝑡) ˜ 𝐴𝛼 (𝜔′) † − 𝐴𝛼 (𝜔′) † 𝐴 𝛽 (𝜔) 𝜌(𝑡)] ˜ + etc. 𝜔,𝜔′ 𝛼,𝛽 (C.8) where: ∫ ∞ ∫ ∞ 1 ′ 𝑖𝜔𝑡 ′ ′ 1 ′ 𝛾𝛼𝛽 = 2 Tr𝑟𝑎𝑑 𝑑𝑡 𝑒 𝐸 𝛼 (𝑡)𝐸 𝛽 (𝑡 − 𝑡 ) 𝜌˜𝑟𝑎𝑑 = 2 𝑑𝑡 ′ 𝑒𝑖𝜔𝑡 ⟨𝐸 𝛼 (𝑡)𝐸 𝛽 (𝑡 − 𝑡 ′)⟩ (C.9) ℏ 0 ℏ 0 defines elements of a spectral correlation tensor. To reduce the double sum over frequencies, the secular approximation is invoked, in which all terms with 𝜔 ≠ 𝜔′ are neglected due to their fast oscillations averaging out their contribution (this step is equivalent to invoking the Rotating Wave Approximation to the interaction Hamiltonian (C.2)). The result is: d ∑︁ ∑︁ ˜ = 𝜌(𝑡) ˜ 𝐴𝛼 (𝜔) † − 𝐴𝛼 (𝜔) † 𝐴 𝛽 (𝜔) 𝜌(𝑡)] 𝛾𝛼𝛽 (𝜔) [ 𝐴 𝛽 (𝜔) 𝜌(𝑡) ˜ + etc. (C.10) d𝑡 𝜔 𝛼,𝛽 Í Í 𝑉 ∫ After taking the continuum limit 𝜔 = 𝑘 → (2𝜋) 3 𝑑k, we must write 𝛾𝛼𝛽 in terms of the fields defined by (3.15 - 3.16), then evaluate the time and k-space integrals. This mirrors steps in the derivation of the photon emission rate (see Section 3.4). However, in the evaluation of the Wigner-Weisskopf approximated integral (3.33), the contribution from the Cauchy Principle Value can in general be retained, giving rise to a Lamb shift. Assuming all dipole strengths are equal (𝑑𝑖 = 𝑑 𝑗 = 𝑑 for all 𝑖, 𝑗), the result is:   1 ∑︁ 𝚪𝑖 𝑗 𝛾= 2 + 𝑖𝛀𝑖 𝑗 (C.11) 𝑑 𝑖𝑗 2 73 where:    3Γ0 sin 𝑘 𝑅 cos 𝑘 𝑅 sin 𝑘 𝑅 𝚪𝑖 𝑗 (r) = (1 − r̂ ⊗ r̂) + (1 − 3r̂ ⊗ r̂) − (C.12) 4 𝑘𝑅 (𝑘 𝑅) 2 (𝑘 𝑅) 3    3Γ0 cos 𝑘 𝑅 sin 𝑘 𝑅 cos 𝑘 𝑅 𝛀𝑖 𝑗 (r) = −(1 − r̂ ⊗ r̂) + (1 − 3r̂ ⊗ r̂) + (C.13) 4 𝑘𝑅 (𝑘 𝑅) 2 (𝑘 𝑅) 3 2𝑑 2 𝜔3 Γ0 = (C.14) 3𝜋𝜖 0 ℏ𝑐3 Substituting this back into (C.10), transforming back to the Schrödinger picture via 𝜌(𝑡) = 𝑈 𝜌(𝑡)𝑈 ˜ †, and using: ∑︁ ∑︁ A=− d̂𝑖 = −𝑑 𝜖ˆ𝑖 (𝜎𝑖+ + 𝜎𝑖− ) (C.15) 𝑖 𝑖 for each component 𝐴𝛼 yields the Lindblad master equation: d 𝑖 ∑︁ Γ𝑖 𝑗 𝜌(𝑡) = − [𝐻𝑑𝑜𝑡 + 𝐻𝑠ℎ𝑖 𝑓 𝑡 , 𝜌] − (𝜎𝑖+ 𝜎 𝑗− 𝜌 − 2𝜎 𝑗− 𝜌𝜎𝑖+ + 𝜌𝜎𝑖+ 𝜎 𝑗− ) (C.16) d𝑡 ℏ 𝑖𝑗 2 where Γ𝑖 𝑗 = 𝜖ˆ𝑖 · 𝚪𝑖 𝑗 · 𝜖ˆ 𝑗 , Ω𝑖 𝑗 = 𝜖ˆ𝑖 · 𝛀𝑖 𝑗 · 𝜖ˆ 𝑗 and: ∑︁ 𝐻𝑠ℎ𝑖 𝑓 𝑡 = Ω𝑖 𝑗 𝜎𝑖+ 𝜎 𝑗− (C.17) 𝑖𝑗 describes a Lamb (energy) shift. (It is apparent that the Cauchy Principle value term Ω is responsible for this shift.) 74