To describe the spontaneous emission of an atomic gas inside a structured environment, we consider the probability that an excited atom emits a photon into a free mode during a time interval τ. From Eq. 8 the probability for this process to occur, regardless of the final vibrational state, is given by the integral P s → s ′ γ = ∫ 0 τ d t ∑ B e i ω s s ′ + ω A B − ω γ t ⟨ Φ B | κ s s ′ γ X | Φ A ⟩ 2 , (12) thus posing frequency missmatch conditions on the modes the atom interacts most strongly with. In general these conditions depend on a vibrational shift ω _AB , but, when the electronic transition frequency is much larger than the relevant CM transition frequencies ω _ss′ ≫ ω _AB , this shift can be neglected [30]. By removing these shifts from the equation the states |Φ _B ⟩ can be averaged out using the completeness of the vibrational states, leading to P s → s ′ γ ≃ S ω a b − ω γ ⟨ Φ A | κ s s ′ γ † X κ s s ′ γ X | Φ A ⟩ , (13) with S ( ω a b − ω γ ) a sharp spectral function that satisfies S ( ω ) ∼ 2 π τ δ ( ω ) as τ → ∞ [31].

Equation 13 should be read as a probability distribution that weights the decay process. The atom interacts with many modes of the environment such that the spontaneous decay rates is obtained from the sum Γ s s ′ ; A = ∑ γ ∂ ∂ τ P s → s ′ γ (14) consistent with the Born and Markov approximations. By removing the correlations that build-up between field and atom we have made the Born approximation, and by extending the time integral without accounting for selfconsistent exchanges we have performed the Markov approximation. In this sense, the spread ω _ss + ω _AB introduced by vibrational states has to be much smaller than the free mode density in order to be neglected. This is usually achieved in free space, and is a good approximation to half cavities where modes acquire a linewidth.

It is now possible to identify the decay rate of an atom inside an structured environment. The transition rates between electronic states incorporate information about the spatial region explored by the atom and the vectorial nature of the radiated field. By using the angular spectrum of the field defined in Eq. 6 the spontaneous rates are found to be [25]. Γ a b ; A ed = ∑ σ = , ± , 0 ∑ γ S ω a b − ω γ d a b ∗ − σ T γ ; A ed σ σ d a b − σ , (15a) Γ a b ; A md = ∑ σ = , ± , 0 ∑ γ S ω a b − ω γ m a b ∗ − σ T γ ; A md σ σ m a b − σ , (15b) Γ a b ; A eq = ∑ σ , σ ′ = , ± , 0 ∑ γ S ω a b − ω γ q a b ∗ − σ − σ ′ T γ ; A eq σ σ ′ σ σ ′ q a b − σ − σ ′ , (15c) where we have used square brackets [x] _σ to denote the σ component of vectors and tensors and adopted the circular polarization basis {e _σ } = {e _± = e _x ± i e _y , e ₀ = e _z }. The contribution of each mode becomes T γ ; A ed σ σ = 2 π ℏ 2 ∫ d Ω k ∫ d Ω k ′ g A k , k ′ ; ω γ f γ ∗ θ k , φ k σ f γ θ k ′ , φ k ′ σ , (16a) T γ ; A md σ σ = 2 π c 2 ℏ 2 ω γ 2 ∫ d Ω k ∫ d Ω k ′ g A k , k ′ ; ω γ k × f γ ∗ θ k , φ k σ k ′ × f γ θ k ′ , φ k ′ σ , (16b) T γ ; A eq σ σ ′ σ σ ′ = 2 π c 2 ℏ 2 ω γ 2 ∫ d Ω k ∫ d Ω k ′ g A k , k ′ ; ω γ k σ f γ ∗ θ k , φ k σ ′ k ′ σ f γ θ k ′ , φ k ′ σ ′ . (16c)

The effect of the vibrational states for all multiple moments is given entirely by a form factor g A k , k ′ ; ω γ = ∫ d 3 X Φ A ∗ X e i k ̂ − k ̂ ′ ⋅ ω γ X / c Φ A X . (17)

The form factor shows that the center-of-mass motion performs an average over the field distributions. The average depends on the atomic center of mass initial state. For an atom prepared in the ground state of an harmonic trap with angular frequency components Λ _x,y,z centered at the position X ₀—which does not need to coincide with the origin used to describe the EM field—the form factor is g A k , k ′ ; ω γ = e i k ̂ − k ̂ ′ ⋅ ω γ X 0 / c ∏ i = x , y , z e − 2 π η i k ̂ i − k ̂ i ′ 2 2 (18)

The average carries information on the trap through the Lamb-Dicke parameters η i = ℏ ω γ 2 / 2 M c 2 Λ i = l i / λ γ relating the ground-state size of the trap l i = 2 π ℏ / 2 M Λ i to the wavelength of the electromagnetic mode λ _γ . In the Lamb-Dicke limit η _i < 1, the atoms are confined below the photon wavelength and will not be heated by light scattering [32]. Notice that even if ω _γ ≫Λ _i , the condition is achieved for ω γ ≪ 2 M c 2 / ℏ taking into account the value of the atom rest energy Mc ². The Lamb-Dicke regime has been achieved for trapped ions [33] and neutral atoms [34].

Through Eqs 15a, 15b, 15c, 16a, 16b, 16c, Eq. 17 we have shown that decay rates depend on the EM mode density evaluated at the atomic resonance frequency and also on the average value determined by the vibrational states. The electromagnetic modes γ that participate in Eqs 15a, 15b, 15c can be constrained by imposing physical boundaries, as done in experiments with optical cavities. Then, the enhancement [29] or inhibition [35, 36] of the spontaneous emission rate depends on the electromagnetic structure of the environment and the location and trapping strength of the atom. Notice that this control is expected not only for electric dipole transitions, but for any multipole transition as we have just illustrated for electric quadrupole and magnetic dipole cases.

It is worth mentioning that for atoms in free-space the EM field is homogeneous and the dipolar and quadrupolar spontaneous emission decay rates are Γ s s ′ ed = 4 | d s s ′ | 2 ω s s ′ 3 3 ℏ c 3 , (19a) Γ s s ′ md = 4 | m s s ′ | 2 ω s s ′ 3 3 ℏ c 3 , (19b) Γ s s ′ eq = | q s s ′ | 2 ω s s ′ 5 15 ℏ c 5 . (19c)

This can be shown by evaluating the T γ ; A ( mult ) integrals with angular spectrum f γ free θ k , ϕ k = 1 | sin θ k | δ θ k − θ γ δ ϕ k − ϕ γ for each polarization, and using the completeness relation of polarizations ϵ _k,λ and wavevectors k. For a systematic approach beyond the quadrupolar interaction, however, it is more convenient to use the spherical modes.

The structure of the radiated field incorporates the atomic symmetries that arise from the central field model. In the mean field scheme, individual radiative electronic transitions involve a single electron that changes its orbital and yields a non-null electromagnetic multipole for the atom as a whole. An EM multipole transition of an atom can be either described as the emission or absorption of a single photon in an appropriate spherical vectorial mode, or as divided along multiple photon channels with wavevectors k and specific angular distribution probabilities. Spherical vectorial modes then lead to a more efficient transfer than their wavevector counterparts.

This efficient transfer is already suggested by the form of the spherical vectorial modes. In free space a monochromatic spherical wave of frequency ω and amplitude A has an angular spectrum f γ = A Y ̃ j m P θ k , φ k . (20)

The spectrum describes the coupling of orbital and polarization angular momenta of photons to yield a total angular momentum {jm} as described by the functions Y ̃ j m ( P ) where P = E, M refers to either transverse magnetic or transverse electric waves. These functions are written explicitly in the Supplementary Appendix SAI in terms of standard spherical harmonic functions Y _ℓm (θ, φ)—that account for the orbital angular momentum—and spherical polarizations e _± = e _x ± i e _y , e ₀ = e _z . Spherical vectorial modes are then characterized by the parameters γ: {ω/c = k, j, m, P} that can be compared to the electronic degrees-of-freedom involved during a transition. In particular to the atomic angular momentum that includes both orbital and spin contributions of the electronic configurations and the angular momentum of the nucleus. Note that the electronic and electromagnetic fields should be described in the same reference frame. In particular, the same quantization z-axis, which is either arbitrarily selected or predetermined by anisotropic environments. The latter can be, e.g., external electric and magnetic fields selected to manipulate the internal or external degrees of the atom, or they could refer to the geometry of cavities designed to control the classical or quantum features of the electromagnetic field.

The radiative transitions that connect two-atomic states are ruled by strict conservation laws for energy, linear momentum, and angular momentum of the atom-radiation system as a whole. These conservation laws are naturally satisfied by spherical modes where electric and magnetic multipole transitions involve a single Y j m ( E ) or Y j m ( M ) photon [37]. As such, they can be used to connect the atomic internal state to the most probable decay processes and the specific EM modes involved. For example, the spatial pattern of magnetic spherical waves Y j m ( M ) is obtained from Eq. 6 as ∫ d Ω k e − i k ⋅ r Y ̃ j m M θ k , φ k = 4 π i − j j j k r Y j m M θ , φ , (21) with j _j (kr) the spherical Bessel functions with k = ω/c and the position r = {r, θ, ϕ}. This form was used to obtain the radiative decay of Eqs 19a, 19b, 19c, but, by portraying the wave as a single vectorial mode, there is no need to account for the contribution of each wavevector and polarization. Similar descriptions are found for the spherical electric modes. Relevant, well-known mathematical properties of spherical vectorial modes are summarized in Supplementary Appendix SAI including a connection to vectorial plane waves. We show explicitly the simple structure of the modes in k-space that is used to define a scalar product from which the EM field can be quantized. These properties emphasize that the mode structure results from the direct coupling of orbital and polarization angular momenta of the field.

Figures 2–4 are used to illustrate the rich spatial patterns of vectorial modes. In an effort to show the symmetries involved, isointensity surfaces are drawn for each polarization e _±,0, and for the total intensity of the spherical waves. As anticipated in the Introduction, the polarization and configuration structure in this subwavelength region is complex. High gradients of the intensity and vortices along a dislocation line for certain polarizations can be found as shown by the phase structure in the XY plane.

We begin with Figure 2 where each polarization component of the electric dipole spherical wave is described by a combination of spherical functions Y _ℓm with ℓ = 0 and 2, as we now discuss. Figure 2 shows the Y 11 ( E ) spherical wave. The σ ₊ polarization component is a combination of both Y ₂₀ and Y ₀₀, the latter being independent of θ and φ and not null at r = 0. Neither Y ₂₀ nor Y ₀₀ components exhibit orbital vortices. For comparison, the σ ₋ component is proportional Y ₂₂ and displays an optical vortex with topological charge 2 along the Z-axis. The π component is proportional to Y ₂₁ with a unit topological charge vortex. Consider next the quadrupole electric spherical wave Y 22 ( E ) shown in Figure 3. An analogous description follows but through a combination of Y _ℓm functions with ℓ = 1 and 3. This implies an odd behavior for quadrupole waves in contrast to the even behavior of dipole electric waves with respect to the parity transformation. Finally, Figure 4 illustrates a Y 11 ( M ) spherical wave. In general, each circular polarization σ _±(π) component of Y j m ( M ) is proportional to the Y _jm∓1(Y _jm ) spherical harmonic. The parity of magnetic spherical modes is even (odd) for even (odd) values of j.

The results are to be compared with standard paraxial optics. For paraxial optics, light polarization is approximately a global concept. There is a main direction of propagation and the polarization vectors are approximately perpendicular to it. This facilitates the identification of processes where only either σ ₊, σ ₋ or π transitions occur; π-transitions require a main direction of propagation perpendicular to a quantization-axis defined by the environment as mentioned in the beginning of this Subsection. For each type of transition the atomic internal states experience a change of the internal magnetic number Δm = 1, −1, 0 respectively. If a realization of a few atomic levels model is desired, a search of an atomic- EM field configuration is performed to maximize the relevance of predetermined internal atomic sublevels that participate in the nonlinear optical process.

Having revisited the connection between electronic processes and nonlinear optical processes, we can now exploit the connection to Maxwell equations. For localized sources which admit an efficient multipole expansion, the radiated field operator E ̂ is known to satisfy the equation [31, 37, 40] − ∇ ⋅ ∇ + 1 c 2 ∂ 2 ∂ t 2 E ̂ = 4 π c ∂ ∂ t J ̂ (32) where the sources J ̂ = J ̂ ( ed ) + J ̂ ( eq ) (33) divide into electric dipole J ̂ ( ed ) and quadrupole J ̂ ( eq ) operators that satisfy J ( ed ) = 1 c ∂ P ̂ ∂ t , J ( eq ) = − 1 c ∂ ∂ t ∇ ⋅ Q ̂ . (34)

It is possible to extend the formalism to include magnetic dipole transitions through a magnetization term J ̂ ( m d ) = ∇ × M ̂ . (35)

In our case, Eq. 32 stems from the interaction Hamiltonian beyond the dipole approximation written in Eq. 7 and the decomposition of the electric fields [4]. The form of the interaction led naturally to the vectorial spherical waves associated to each atomic radiative transition and then to source terms of this form. The particular atomic evolution follows from a master equation evolution (as defined in Supplementary Appendix SAII), but the form is universal.

Note that radiation can also be generated by so-called “free” sources. An example is radiation scattered by structureless ideal charges such as electrons, i.e., the Compton effect, which has frequently been described using a plane wave expansion of the quantized EM field [31]. Recent theoretical studies concern the Compton effect with photons associated to structured beams [41], and the implementation of the inverse Compton effect for twisted electron beams is an active area of research [42].

For the three-wave mixing process considered here, the induced moments oscillate at frequencies ω ₁, ω ₂ and ω ₃. As such, the evolution of the field amplitudes can be decomposed into a Fourier series that leads to − ∇ ⋅ ∇ − ω 1 2 c 2 E γ 1 = − 4 π ω 1 2 c 2 ∇ ⋅ Q − ω 1 , (36a) − ∇ ⋅ ∇ − ω i 2 c 2 E γ i = 4 π ω i 2 c 2 P s − ω i , i = 2,3 (36b) whenever the temporal phase matching condition is imposed ω ₁ = ω ₂ + ω ₃. From them, conditions yielding a high efficiency of parametric processes based on atomic coherence can be obtained [43]. Note that the full quantum treatment implicit in Eq. 32 allows to study correlations of the quadrature equations of the EM modes, including squeezing conditions [44].

Most implementations of nonlinear optics processes consider input paraxial beams driving the atomic gas. The set of paraxial modes { γ i p x } neither matches the atomic symmetry nor provides any guarantee of an adequate description of the EM field in the radiation zone [45]. Within the paraxial regime, the theoretical description of the classical and quantum properties of light is incorporated using the slowly varying amplitude approximation [1]. It states that the relative change in the amplitude per wavelength is small and there is a main direction of propagation z of the light fields so that the dominant spatial variation of the amplitude can be approximated by ( ω / c ) ∂ z E ( γ i ) . This is not valid when the electromagnetic beams are focused, but, as we now show, the vectorial waves provide a form to describe the evolution.

We can now use this formalism to describe the phase matching conditions that rule the underlying nonlinear processes. For plane waves phase matching conditions establish the relations between the wave vectors of the involved EM modes that guarantee the highest efficiency of an optical nonlinear process. In the quantum realm, these equations are interpreted as the conservation of linear momentum of the participating photons. For a homogeneous atomic gas, the three-wave mixing phase matching conditions are k 1 = k 2 + k 3 . (40)

For beams exhibiting a common dislocation line—and correspondingly a non-trivial local orbital momentum along that line—the phase matching conditions correlate the topological charge of the vortices so that the angular momentum of the photons is conserved [47, 48]. For an isotropic atomic sample and for modes exhibiting optical vortices of topological charge m _i along a common axes, for our process m 1 = m 2 + m 3 . (41)

In this Section we show how they naturally emerge from the dyadic treatment of the scattered field. To that end it is just necessary the integrate the Maxwell Equation Eqs 38a, 38b over the spatial variables of the localized source ∫ d 3 X ′ j ℓ ω r / c Y ̃ ℓ m P ∗ θ ′ , φ ′ ⋅ J X ′ , t ′ (42)

If the 3D sample is in an isotropic trap, the angular integration acquires an analytic expression (some formulae useful for the calculation of the spatial derivatives of the spherical waves can be found in the Supplementary Appendix SAI). The integrals are, in general, a linear combination of the Wigner 3-j symbols ∫ d Ω Y ℓ 1 m 1 θ , φ Y ℓ 2 m 2 θ , φ Y ℓ 3 m 3 θ , φ = 2 ℓ 1 + 1 2 ℓ 2 + 1 2 ℓ 3 + 1 4 π ℓ 1 ℓ 2 ℓ 3 0 0 0 ℓ 1 ℓ 2 ℓ 3 m 1 m 2 m 3 (43) available in most numerical platforms. The phase matching conditions result from the identification of non-null Wigner 3-j symbols, m i ϵ − ℓ i , … , ℓ i , i = 1,2,3 , m 1 + m 2 + m 3 = 0 , | ℓ 1 − ℓ 2 | ≤ ℓ 3 ≤ ℓ 1 + ℓ 2 .

These phase matching conditions can be interpreted as a conservation of the total angular momentum and not just its z-component for the photons involved in the three-wave mixing process.