Interaction of atomic systems with quantum vacuum beyond electric dipole approximation

The photonic environment can significantly influence emission properties and interactions among atomic systems. In such scenarios, frequently the electric dipole approximation is assumed that is justified as long as the spatial extent of the atomic system is negligible compared to the spatial variations of the field. While this holds true for many canonical systems, it ceases to be applicable for more contemporary nanophotonic structures. To go beyond the electric dipole approximation, we propose and develop in this article an analytical framework to describe the impact of the photonic environment on emission and interaction properties of atomic systems beyond the electric dipole approximation. Particularly, we retain explicitly magnetic dipolar and electric quadrupolar contributions to the light-matter interactions. We exploit a field quantization scheme based on electromagnetic Green’s tensors, suited for dispersive materials. We obtain expressions for spontaneous emission rate, Lamb shift, multipole-multipole shift and superradiance rate, all being modified with dispersive environment. The considered influence could be substantial for suitably tailored nanostructured photonic environments, as demonstrated exemplarily.

www.nature.com/scientificreports www.nature.com/scientificreports/ Furthermore, and potentially more interestingly, nanostructured environments may also open light-matter interaction channels beyond that one corresponding to coupling of electric field with electric dipoles. For example, light concentration in nanoscopic regions causes modulations of electromagnetic field at spatial distances comparable to the size of atomic systems, which may couple to electric-quadrupolar or higher-order moments. In addition, due to their large refractive index, dielectric nanomaterials offer the possibility of strong concentration of magnetic fields 17 . This prompts to consider not just electric-multipolar contributions but at the same time their magnetic counterparts. Until now, enhancement of the rate of a magnetic dipole emission by a nanostructure was considered 18 and reported experimentally in lanthanide ions [19][20][21][22] . Large enhancement of quadrupole [23][24][25] and even higher-order transitions was equally predicted. Transitions driven with several multipolar mechanisms have been considered 26 and observed 16,27 respectively in semiconductor quantum dots and transition-metal/lanthanide ions. Recently, it was pointed out that the simultaneous existence of several interaction channels enhanced in intensity with photonic nanostructures might open new avenues on the route to control spontaneous emission, related to their controlled interference 28 . To account for these effects, treatments beyond the point-dipole approximation were developed [29][30][31] , based on spatial integration of the transition elements expressed in terms of wavefunctions of the atomic system. Nonclassical corrections to these transition rates originating from effects such as nonlocality, electronic spill-out and Landau damping have been considered in ref. 32 . Although such approaches take into account all multipolar orders of light-matter interaction, they require information on the atomic system's wavefunctions, which are in general complex multidimentional structures. They are hardly accessible experimentally and their calculation may be a significant computational effort already in case of few-electron systems.
In this work, we develop an analytic theory based on the multipolar coupling Hamiltonian, in which the spatial extent of the atomic system is taken into account in terms of multipolar transition moments rather than wavefunctions. These moments can be evaluated with smaller computational effort than full wavefunctions, found in literature 33,34 or measured 21,35 . Contrary to the previous works mentioned above, we evaluate not just transition rates, but also energy shifts. In the case of an individual atomic system we evaluate Purcell enhancement of spontaneous emission and Lamb-shift modification in nanostructured dispersive environment. Additionally, we provide expressions to evaluate interaction strengths and collective emission rates of multiple atomic systems immersed in the same environment, while considering the magnetic-dipolar and electric-quadrupolar interaction channels. Our method exploits field quantization in dispersive media 36 , in which the optical properties of the environment are expressed in terms of the electromagnetic Green's tensor, defined by the spatial and spectral dependence of the electric permittivity. Our result is an extension to those results previously obtained in electric-dipole approximation in refs. 9,37 . Here, we include two next-order terms of the multipolar coupling Hamiltonian 38 , namely the magnetic dipole and the electric quadrupole. This unlocks qualitatively new routes to tailor light-matter coupling exploiting interference effects.
The work is organized as follows: in the next section we introduce the general framework used to describe light-matter coupling in the presence of dispersive and absorptive structured materials beyond electric-dipole approximation. Examples of application of the theory to specific geometries are given at the end of the Results section. Detailed numerical results and lengthy calculations are provided in the Supplementary Information available online.

Results
In this section, we introduce the general framework used to describe light-matter coupling in the presence of dispersive and absorptive structured materials beyond electric-dipole approximation. After giving a statement of the problem, we offer expressions for spontaneous emission rate and Lamb shift of a single atomic system in such environment at first. Next, we generalize these expressions to many-atom systems. Finally, examples of application of the theory to specific geometries are given. Lengthy calculations are documented in Supplementary Information.
Atomic system. We assume that the atomic system can be approximated by two active energy levels, separated by energy  0 ω , where  stands for the reduced Planck's constant. The corresponding ground and excited states are denoted by | ⟩ g and | ⟩ e , respectively. The system is fully described by a set of Pauli operators: the lowering operator g e σ = | ⟩⟨ | and the inversion operator | ⟩⟨ | | ⟩⟨ | e e g g z σ = − , following the usual commutation rules The free Hamiltonian of the system reads † 0 0 . This Hamiltonian can be generalized to the case of multiple emitters in a straightforward manner, as done in one of the following subsections.
Quantized electromagnetic field in dispersive media. In this work we follow the quantization scheme in dispersive and absorbing media, developed in refs. [39][40][41][42] . We restrict ourselves to nonmagnetic matter with relative permeability µ = 1. The constitutive equation relating the temporal Fourier components of the displacement field ω D r ( , ), the electric field E r ( , ) ω , and medium polarization P r ( , ) ω in an absorbing medium takes the form N describes a noise contribution to the polarization P r ( , ) ω arising from vacuum fluctuations in an absorbing medium. The related noise current density reads N N . Since we investigate vacuum-induced effects, we are interested in the case of a vanishing mean electric field. Then, the only field is related to noise current fluctuations and can be expressed as N 0 3 www.nature.com/scientificreports www.nature.com/scientificreports/ where µ 0 stands for vacuum permeability.
The dyadic tensor G r r ( , , ) ω ′ is the full electromagnetic Green's tensor characterizing the environment, determined from the Maxwell-Helmholtz equation with I representing the unit dyadic. The boundary condition for the Green's tensor at infinity reads G r r ( , , ) . The Green's tensor serves as the kernel that connects the electric field ω E r ( , ) with its source at position ′ r . The requirement for the canonical equal-time commutation relations of quantized fields to hold allows one to express the noise current in the form [39][40][41] Here, 0 ε represents the electric permittivity of vacuum and ε ε ε ω ω ω = +i r r r ( , ) Re ( , ) Im ( , ) is the relative permittivity of the dispersive and absorptive medium surrounding the atomic system. For simplicity, we assume isotropic media, so that the permittivity can be expressed as a scalar function. The bosonic operator fields take the form and e j is a unit vector in the j th direction. They obey the following commutation relations [39][40][41] f Finally, the electric field can be expressed in terms of bosonic operators as follows where c is the vacuum speed of light.
In the Schrödinger picture, the positive frequency part † . Similarly, the positive frequency part of the magnetic field is expressed given by , where the frequency components of the magnetic fields are connected to the corresponding electric field through coupling Hamiltonian. The electric dipole approximation is based on the assumption that the electric field can be approximated as constant across the spatial extent of the atomic system. In particular, this means that direct coupling of the emitter with the magnetic field is ignored and spatial modulations of the electric fields are neglected as well. Typically, the above assumption is very well met: the scale of spatial modulations of the electric field is set by its wavelength, usually in the optical or near-infrared range, while the modulations of the field's envelope are even slower and thus negligible. This assumption holds true in free space or traditional cavities, if the atomic system is represented by an atom or a molecule. However, the assumption might no longer be applicable if the emitter is positioned within a subwavelength electromagnetic hotspot, e.g. in photonic crystal cavities, in close vicinity of plasmonic nanostructures 16 , near picocavities or even in free space when geometrically large emitters like semiconductor quantum dots are considered 11,26,31 . For this reason, we include two higher-order terms of the multipolar coupling Hamiltonian, which include first-order spatial derivatives of the electric field 38 . The Hamiltonian is given in the rotating wave approximation, valid as long as the coupling strengths are small with respect to the transition frequency 0 ω where the fields are evaluated at the position of the atomic system r 0 , and for brevity we denote ≡ We assume that the electric field derivatives exist at this position, i.e. the atomic system should not be placed directly at an interface between two different media. Above,  = g e d d and = g e m m  are the electric and magnetic transition dipole moment elements, and  = g e Q Q is the electric transition quadrupole moment element, respectively. The dot denotes the standard scalar product of vectors, the double dot product of tensors is defined as In the case of a real electric quadrupole moment tensor, the quadrupolar contribution to the coupling Hamiltonian can be equivalently rewritten as . Please note that different degrees of freedom in the operators above are denoted as follows: • The degree of freedom related to the two-dimensional Hilbert space spanned by | ⟩ | ⟩ g e { , } is already included in the symbols of transition moments, and is relevant for Hermitian conjugation, e.g.
• the dipole moment vectors and the quadrupole moment tensor have elements corresponding to the x y z , , spatial directions, e.g. d i , Q ij , responsible for the orientation of the multipolar moment, • finally, the fields depend on position in space r, and each element of the Green's tensor is a function of the observation point r and the source point r′.
emission properties of single atomic system. To derive the spontaneous emission rate of an atomic system, we proceed as follows: First, Heisenberg equations are found for the atomic and the field operators. The equations for the field are then formally integrated, so that the field operators are expressed through the atomic ones. The result is then inserted into atomic equations, so that field variables are completely eliminated from the description. The resulting complicated integro-differential equations are simplified in the Markovian approximation, where the memory effects are neglected. As a result, we obtain effective dynamics of the atomic system alone. The procedure is a generalization of the one introduced in ref. 9 , where only the electric dipole interaction term was taken into account. Since the equations tend to be lengthy, we describe the consecutive steps in detail in Section 1 of the Supplementary Information. The effective evolution of the atomic system reads where the fields are always evaluated at r 0 , is the identity operator in the atomic Hilbert space, γ is the spontaneous emission rate, and δ stands for the analogue of Lamb shift, calculated beyond the electric dipole approximation. Explicit expressions for these quantities are given in the following. The "free" subscript corresponds to free fields, i.e. fields that are not influenced by atomic back-action. In the vacuum state, these fields account for vacuum fluctuations and their mean value vanishes. The influence of the photonic environment is taken into account through the modified Green's tensor. The emission rate γ includes contributions from the electric dipole, magnetic dipole, and electric quadrupole channels, as well as their interference. The rate is expressed as { , , }, G mn ″ denotes the imaginary part of an mn element of the Green's tensor and ε pkm is the Levi-Civita antisymmetric symbol. The derivatives in Eq. (10) should be evaluated at the position of the atomic system. Please note that the "generalized moment" is in fact a differential operator that acts on the Green's tensor, i.e. the "moment" combines atomic and field properties. We stress that D r becomes a purely atomic quantity only in the electric dipole approximation. Please note that in this case expression (9) is reduced to the well known form † . The Lamb shift in Eq. (7) is expressed through a principal-value integral: mn m r n r mn r r r r General comments. Before we move to the next section, we will make a few general comments.
The emission and frequency shifts originate from coupling to quantized electromagnetic vacuum and can be derived in the quantum model, albeit their enhancement with respect to free space can be expressed with classical quantities. The enhancement depends on the geometry of the surroundings of the emitter, which determines the properties of the field, here accounted for in terms of the Green's tensor. The tensor is governed by the classical www.nature.com/scientificreports www.nature.com/scientificreports/ electric permittivity. Such description does not take into account the electron tunneling effects which might appear between the atomic system and metallic nanoparticles for separations below 1 nm 45  can be decomposed into a sum of a homogeneous term G h and a scattered part G s 46 . The homogeneous term corresponds to the response in free space or a homogeneous medium, while the scattered one describes influence of scatterers in the environment, e.g. extended interfaces among different media, nanostructured particles, or photonic crystals. The contribution of the homogeneous part of the Green's tensor is already included in the homogeneous-medium spontaneous emission rate or the respective Lamb shift, while the contribution of the scattered part is of general interest. The scattered contribution is frequently finite at the origin, as ′ → r r r , 0 . From an analysis of the generalized multipole moment, we find that the electric dipole component depends on the imaginary part of Green's tensor, while the electric quadrupole and the magnetic dipole components are proportional to the sum and difference of the corresponding derivatives of the imaginary part of the Green's function: , respectively. A simple π 4 -rotation of the coordinate system shows that these can be independently tailored, since there are no general restrictions on the ratios of values of a function and its derivatives in different directions at a given point. This observation is a starting point to consider their full interference and engineer environments which might support it 28 .
Generalizing the expressions for the transition rate beyond the electric dipole approximation not only allows to consider corrections to the atomic systems' dynamics in cavities of extreme geometries. More importantly, it is a tool to describe, e.g. optical activity of chiral molecules for which an interplay between the electric and magnetic dipolar coupling of matter and light, i.e. interference of the two transition mechanisms, plays a crucial role. In centrosymmetric systems, parity is a good quantum number, allowing one to identify transitions either described by the electric dipole mechanism or a combination of electric quadrupole and magnetic dipole. Such systems could be considered sources of photons with well-defined parity, corresponding to a given transition mechanism. emission properties of multiple atomic systems. The same formalism can be applied to the case where multiple two-level atomic systems, indexed with α, share the same photonic environment. The systems do not need to be identical, but we assume the separations of their transition frequencies ω α to be small with respect to the scale of spectral modulations of the properties of the photonic environment. This assumption will be relevant for the Markovian approximation. We additionally assume that the systems do not directly interact. However, the shared environment can be a carrier of interatomic coupling, as we will see below.
In the case of multiple atomic systems, the Hamiltonian from Eq. (6) should be generalized to the form For a better understanding, it is useful to note that the same equations can be derived for a collection of atomic systems described by an effective Hamiltonian of the form are expressed through multipolar elements and differential operators, as given in Section 2 of the Supplementary Information. From Eq. S17 in Section 3 of the Supplementary Information it follows that if the transition frequency is sufficiently large, the expression for the coupling can be simplified to the form where the principal value integral is no longer present and the integration is now performed along the imaginary axis. There, the Green's tensor shows greater numerical stability due to its decaying rather than oscillating character.
The multipole-multipole interaction strength ξ αβ is a generalization of the dipole-dipole coupling which in free space scales as R 3 − , R being the interatomic distance. Contrary, a modified photonic environment, e.g., in a photonic crystal, near a nanoparticle or a nanowire, might allow not only for stronger interactions but also for extended interaction distances 9 . Due to the enhancement of off-diagonal elements of the Green's tensor, long-range coupling of multipoles or, in general, arbitrary nonparallel orientations may be enabled. Due to the strong field localization, corrections beyond the electric dipole approximation in such systems may be significant, and coupling of different multipoles is possible. Interference of different interaction components may lead to further enhancement or suppression of interaction strength, resulting in a corresponding modification of interaction distance. Please note that due to the large width of the peak in the density of states, assumed in Section 2 of the Supplementary Information, atomic systems with slightly different transition frequencies may in general be coupled. www.nature.com/scientificreports www.nature.com/scientificreports/ The dissipators in Eq. (14) include emission rates of an individual, th α atomic system γ αα (reducing to γ from the previous Section in case of a single system), and collective decay rates γ αβ . They arise because each atomic system is capable of modifying the photonic environment of the others and they are defined as www.nature.com/scientificreports www.nature.com/scientificreports/ examples. In this section we apply the formulas derived above first to the case of a homogeneous background dielectric and second to an exemplary selected nanostructured environment into which the atomic system is placed. In the first considered case, our goal is to retrieve familiar scaling of different multipolar components of spontaneous emission rates with different powers of refractive index 48,49 . In the latter case, we demonstrate that in a suitably engineered environment contributions beyond electric dipole may have a significant impact on the atomic system's dynamics.

Homogeneous background medium. In a homogeneous, isotropic, and infinitely extended medium the
Green's tensor takes the form is the distance between the source and observation point where the field is to be evaluated, and the wave number in the homogenous medium reads , ω n( ) being a position-independent refractive index.
In Section 4 of the Supplementary Information we show that away from the medium resonances the imaginary part of the Green's tensor is diagonal in the limit of → R 0 with R jk j k 0 3 for k j ≠ , and is exactly 0 for R 0 = . The diagonal elements however, are finite and read as Inserting the limit at R 0 → in Eq. (9), we retrieve the Weisskopf-Wigner result for the electric dipole contribution to spontaneous emission Higher-order terms may be evaluated based on derivatives found in Eqs. S27 & S28 of the Supplementary Information. As the result, we find the familiar result for transition rates via higher-order channels 48,51 , which scale with the third power of the refractive index 48,49 . Similarly, multipole-multipole interaction terms can be retrieved.
Pair of metallic nanospheres. To offer also an example of a structured photonic environment, we consider here a pair of silver nanospheres of 40 nm diameter, separated by a 6 nm gap inside of which an atomic system is positioned. The chosen coordinate frame is shown in the system sketch in Fig. 1. The Green's tensor was calculated using the MNPBEM toolbox for MATLAB 50 : a Maxwell equation solver based on the Boundary Element Method [52][53][54] . Full Maxwell equations were solved. Silver was modeled using data from ref. 55 . The tensor is calculated at the frequency 2 789 0 ω π = × − ps 1 on a 4 nm 3 nm × grid located in the symmetry plane = y 0 between the nanospheres, as marked by the rectangle in the system sketch in Fig. 1. The Green's tensor's elements and derivatives are presented in the Supplementary Information.
We now consider a two-level atomic system with transition frequency 0 ω . We choose a transition electric dipole moment to be oriented along the x axis, with a length of = d ea x 0 , where e stands for the elementary charge and a 0 is the Bohr radius. The magnetic transition dipole moment parallel to the z axis is µ = m i 2 z B , with µ B standing for Bohr magneton. The electric transition quadrupole moment in the xy plane is set to Q Q ea The chosen values correspond to all moments equal to 1 atomic unit, i.e. values characteristic for atoms and molecules. The transition rates depend on the position of the atomic system with respect to the nanospheres, and are shown in Fig. 1. We only consider the atomic system's positions in the rectangle from the system sketch in Fig. 1. In free space, the rates through respective channels are given by Eqs. (24)(25)(26) Hz. This means that the ED channel dominates the emission, and if it is present the visibility of the other channels is suppressed. In general, in the gap between the nanospheres the enhanced transition rates in all channels exceed the free-space values by several orders of magnitude: up to 4 in case of the ED and even 7 (8) for the MD (EQ), improving the impact of the higher order channels with respect to the dominant electric dipole. As a result, the channels beyond the electric dipole make a significant contribution www.nature.com/scientificreports www.nature.com/scientificreports/ to the overall transition rate. This contribution includes the interference channels: for the studied nanostructured geometry and orientations of the multipolar moments we find the interference between the MD and EQ to be relatively strong, while the interference with the ED channel is weaker. However, in general other geometries or orientations of the atomic systems could be considered where these channels, here weak, could be tuned 28 . Among possible applications this suggests potential for enhancement of optical activity, in particular circular dichroism. Among the two considered higher-order channels, the magnetic dipole transition channel dominates by two orders of magnitude over the electric quadrupole one. However, the latter is manifested through interference, which we find always destructive in the investigated region, and whose contribution to the transition rate in channels beyond the ED is of the order of 10%.
We now evaluate the correction of the frequency shift, i.e. the Lamb shift, of the atomic system coupled to the quantum vacuum near the two nanospheres. Please note that the full Green's tensor diverges at the origin. However, we are only interested to evaluate the correction of the shift with respect to free space, determined by the scattered part of the Green's tensor. The result obtained for selected positions within the gap is shown in Fig. 2.  www.nature.com/scientificreports www.nature.com/scientificreports/ The relative values of the shifts δ ω < − 10 3 0 in the investigated domain are dominated by the electric-dipole and the comparable magnetic-dipole channel. Please note the opposite sign of the interference term between the two.
Finally, we evaluate the interaction strengths ξ (Fig. 3) and collective emission rates γ 12 (Fig. 4) for a pair of atomic systems, one of which is positioned at the center of the domain in between the nanospheres, while the position of the other is swept across the simulation domain. The evaluated interaction strengths reach tens of THz throughout the simulation domain. Please note that atomic systems separated by very small distances might interact through mechanisms other than the optical one, e.g. through electron exchange. The total interaction strength and collective decay rate is a sum of contributions from all channels. Please note that in the case of interaction strength we find the contributions from all investigated multipoles to be significant. Contrary, the collective decay rates are dominated by the ED channel.

Discussion
We have studied dynamics of atomic systems coupled to a photonic environment in its vacuum state. The environment is described in terms of the electromagnetic Green's tensor and the interaction contributions beyond the paradigmatic electric dipole approximation have been included. The derived formalism allows us to evaluate dynamical parameters characterizing optical properties of atomic systems: both the individual and the collective contributions to energy shifts and decay rates. Inclusion of terms beyond the electric dipole approximation allows us to study the role of higher multipolar channels, including enhancement or suppression through interference of different interaction mechanisms. Examples of phenomena to be investigated are optical activity, multipole-multipole interactions between atomic systems and spontaneous emission suppression due to interference of different mutipolar channels in tailored photonic environments.

Methods
Our approach is based on quantization scheme introduced in refs. 39,42 . To derive equations presented in this work we apply the Markovian approximation, as described in detail in the Supplementary Information available online.
Green's tensor corresponding to the pair of metallic nanospheres in the Examples section was calculated using the MNPBEM toolbox for MATLAB 50 : a Maxwell equation solver based on the Boundary Element Method 52-54 . Full Maxwell equations were solved. Silver was modeled using data from ref. 55 . www.nature.com/scientificreports www.nature.com/scientificreports/