Strong multipolar transition enhancement with graphene nanoislands

During the past half century, a major approximation was natural in the field of light-matter interaction: the point-dipole model. It was assumed that the wavelength is much larger than the size of the emitting atom or molecule, so that the emitter can be described as a single or a collection of elementary dipoles. As it is legitimate for visible light, the approximation does no longer hold near plasmonic nanostructures, where the effective wavelength can drop below 10 nm. In that case deviations arise from the approximate model. First, the emitter spatial extent influences the far-field spectrum. Second, high-order transitions beyond the dipolar ones are not forbidden anymore. Going beyond the approximation requires intensive numerical efforts to compute the photonic response over the spatial extent of the emitter, since the complete Green's function is required. Here, we develop a general model that computes the multipolar transition rates of a quantum emitter in a photonic environment, by computing the Green's function through an eigenpermittivity modal expansion. We apply the method on graphene nanoislands, and we demonstrate a local breakdown of the selection rules, with quadrupolar transition rates becoming 100 times larger than dipolar ones.

the wave vector magnitude |k| by confining light in a nanophotonic structure. The wave vector can be written as k = η 0 ω/c, with the confinement factor η 0 the ratio between the vacuum and effective wavelength. In this case, higher-order transitions are enhanced by a factor η 0 to the power of the considered order (for example the octupolar transition is enhanced by a factor η 2 0 ). 5 Under these conditions, plasmonic nanoantennas are ideal candidates to enhance higher order transitions. For instance, in noble metals, forbidden quadrupolar transitions are enhanced for emitters close to tips, 6 interfaces, 7 nanowires, 8 nanogaps, 9 arrays 10 and gold dimers. 2,11,12 The strongly confined graphene plasmons (η 0 = 150 − 300, depending on the absorption losses 13,14 ) form an excellent platform for high-order transitions, which can occur efficiently, even similar to dipolar transitions. 5 In case of extremely high confinement (η 0 > 500) of plasmons in a two-dimensional material sheet, higher-order transitions rates can surpass lower order transitions, hence breaking the conventional selection rules. 5,15 Generally, accessing high-orders allows to probe a much larger range of the electronic energy level structure of an emitter, finding a way to a multiplex and broadband spectroscopy platform. 5,11 These higher-order transitions already play an important role in spectroscopy of many relevant chemical species, from individual atoms 16,17 to larger molecules with high symmetry, such as dihydrogen, carbon dioxide, methane, and benzene. 18,19 In photochemistry, enhancing the magnetic dipole transition in oxygen is interesting for photochemical reactions. 20 Finally, interference effects between multipolar orders can occur: the possibility of complete suppression of a certain transition through interference is required for many applications in the context of quantum computing, quantum storage, and quantum communication. 3,21 Despite its high potential and these developments, the field is currently limited by the difficulty of computing the electromagnetic environment of the emitter. Indeed, computing the spontaneous emission rates of a quantum emitter requires the knowledge of the electromagnetic field profile over the spatial extent of the wavefunctions of the emitter. 22 Usually, the problem is solved for absorption rates: in that case a plane wave excites a nanophotonic structure and the near-field is extracted. 1,6,9,10 This is a straightforward routine for conventional numerical methods such as the finite-element method 23,24 or the finite-difference time-domain method. 25 For spontaneous emission however, the knowledge of the vacuum field is essential. As a first approximation, one can resort to symmetric problems, 26 or consider only the relevant (properly quantized) modes of the structure for the process. In order to compute advanced photonic structures, a modal-based approach is very useful: a single simulation that determines the modes (e.g. of a cavity) is required to know the full spatial variation of the Green's function. 27 The eigenpermittivity mode expansion is particularly suited for the spontaneous emission of an emitter, for which the emission frequency is fixed. Eigenpermittivity modes have a permittivity eigenvalue that pertains only to a scattering element, which spans a finite portion of space. As a result, the normalization is trivial. Furthermore they are orthogonal and appear to form a complete set. 28,29 Once computed for a scatterer at a fixed wavelength, they straightforwardly give the optical response for any material constituting this scatterer. These modes have been derived in the seventies in the quasi-static approximation, and were used to derive bounds for scattering problems, 30 to study spasers, 31 disordered media, 32 and second harmonic generation. 33 The formalism (called GENOME for GEneralized NOrmal Mode Expansion) was recently extended beyond the quasistatic approximation by computing the electromagnetic fields and the associated Green's function of open and lossy electromagnetic systems, in particular for general nanoparticle configurations using a commercial software (COMSOL Multiphysics). 29 In this paper, we derive a general method to compute the transition rate of a quantum emitter influenced by its electromagnetic environment in the weak coupling regime (Sec-tion 2). The electromagnetic environment is embodied by the Green's function (computed with GENOME) and the quantum emitter is described by its wavefunction (a hydrogen-like emitter). Then, we apply this method to compute the electric dipolar, quadrupolar and octupolar transition rates of the emitter in the vicinity of graphene nanoislands for different geometries (triangle, square and crescent), showing strong enhancement of the transition rates (Section 3.1). Afterwards, we show that the graphene doping can be tuned to select particular transitions in Section 3.2, before demonstrating a local breakdown of the selection rules in Section 3.3.

Method
We consider the spontaneous emission of atomic hydrogen-like emitters into plasmons given by the minimal coupling Hamiltonian 34,35 with p i , r i and σ i the impulsion, position and spin of the i th electron, e the electronic charge, m e the electron mass, A and B the vector potential and magnetic field, H SO the spin-orbit coupling, and H e−e the electron-electron interaction. f † j (r, ω) and f j (r, ω) are the creation and annihilation operators, respectively.
For the interaction Hamiltonian, we neglect the ponderomotive potential (A 2 term) and the B term as the latter is negligible for non-magnetic structures. 36 We also neglect the p · A term since we use the Coulomb gauge (the scalar potential is identically 0) and therefore ∇·A = 0 except at an interface: With the atom-interface distance we consider, and the rapid decay of the atom wavefunctions, the contribution of this term will be negligible. Writing the vector potential with the Green's function of the system, and applying the Fermi Golden rule, one finds (for details, see 5 ) where ε 0 is the vacuum permittivity, c is the speed of light, ψ g and ψ e are the atomic wavefunctions of the ground and excited state of the emitter.Ḡ(r, r ′ , ω 0 ) is the Green's function of the Maxwell equations and satisfies ∇ × (∇ ×Ḡ) − ω 2 c 2 ε r (r, ω)Ḡ =Īδ(r − r ′ ), with ε r the relative permittivity and δ the delta function. 37 Equation 5 computes the transition rates of any emitter (via the atomic wavefunction) within any photonic environment (described by the Green's function) in the weak coupling regime. As mentioned, the spatial variation of the Green's function is known analytically for uniform media and for simple geometries. However, more complex structures need to be evaluated numerically, with high computational cost. 38 Equation 5 is also resource demanding since the integration is performed over 6 dimensions (r and r ′ ). In order to resolve these two issues, we resort to GENOME, 29 with the advantage that one modal computation allows the knowledge of the complete spatial Green's function.
In GENOME, the problem is written at a fixed frequency, and the mode-related eigenvalue is the permittivity. This formulation suits well the determination of spontaneous emission rates since the emission frequency is determined by the emitter. The modes E m (r) of the scatterer are solved with a commercial finite-element based software (COMSOL Multiphysics) 38 and the Green's function becomes where m is the mode number, k is the vacuum wave vector, ε m is the eigenpermittivity, ε i is the permittivity of the scatterer, ε b the permittivity of the background material and G 0 (|r − r ′ |) the Green's function of vacuum, which has an analytical form. 37 Inserting Equation 6 in Equation 5, we immediately see that the rate is a sum of two contributions Γ = Γ 0 + Γ s , with Γ 0 the decay rate in vacuum (based on the contribution ofḠ 0 (|r − r ′ |)), and Γ s depending on the modes and hence, the nanophotonic structure.
Focusing on this Γ s contribution, we can write where we defined γ m = . Note that the adjoint field (E † ) is the transposed vector, and there is no complex conjugate. 29 Since we can choose the wavefunctions to be real, the complex conjugate for the wavefunctions disappears, and we can integrate separately for r and r ′ . Both integrations give the same value, leading to: Im γ m ψ e (r)E m (r)∇ r ψ g (r)dr 2 (8) Finally, the transition rates are obtained with a 3 dimensional integration over the wavefunctions and the modes profiles, with a sum that can be truncated once the convergence is sufficient (40 modes in our case, see Supplementary Information). Note that when the integral is computed, the rate can be known for any material constituting the scatterer, enclosed in the parameter γ m . In that regards, graphene is the perfect candidate as it can be tuned to match a particular resonance (see Section 3.2).
In this work the graphene nanoislands are modeled with an effective thickness of t = 1 nm. The graphene permittivity (ε i ) is deduced from the surface optical conductivity (σ = σ intra + σ inter ) with ε i = 1+iσ/ωε 0 t. The optical conductivity is derived within the local random-phase approximation model 39,40 and is the sum of the two following contribution: with T the temperature, k B the Boltzmann constant and E F the doping level of graphene.
The scattering lifetime of electrons in graphene depends on the doping and is given by with the impurity-limited DC conductivity µ ≈ 10000 cm 2 /(V s) and v F = 10 6 m/s the graphene Fermi velocity. 41,42 The integration of Equation 8 is successfully compared to direct simulations of dipolar and quadrupolar transitions in the SI, showing great convergence with only 40 modes (1% relative error). In the next section, we implement Equation 8 to compute the rate of a H-like atom close to graphene nanoislands of varying geometry.

Results and discussion
We apply our method to compute the electric dipolar (E1), quadrupolar (E2) and octupolar (E3) transition rates of a H-like atom close to a graphene sheet with triangle, square and crescent geometry. We consider the transition series 6p, d, f → 4s, which are E1, E2 and E3 transitions, respectively. Note that we suppose that the angular magnetic number remains m = 0 during the transition, and we rotate the emitter wavefunctions to match the corresponding classical point-dipole orientation. The free-space wavelengths of the transitions are all λ = 2.63 µm, and in the whole paper the emitter is situated 5 nm above the graphene surface.
We then discuss the rate dependence on the graphene doping (Section 3.1) and we demonstrate the advantage of graphene tunability for multipolar transitions (Section 3.2). Finally, we optimize a configuration where the conventional selection rules break down, i.e. when the quadrupolar transition rate dominates the dipolar one (Section 3.3). Figure 1 shows the dipolar (Γ D ), quadrupolar (Γ Q ) and octupolar (Γ O ) transition rates of a H-like emitter in the vicinity of a graphene nanoisland for three geometries: square, triangle and crescent shape. The rates are normalized by the dipole emission rate in free space Γ D,0 = 4.484 × 10 5 s −1 . The latter is obtained by integrating Equation 5 in free space, and is in perfect agreement with the experimental values 43 (more information in SI). One can see, for example, that the octupolar rate is strongly enhanced with respect to vacuum: it is up to 300 times stronger than the dipolar rate in free space for the triangle geometry.

Transition rates
One observes that the strongest quadrupolar and octupolar rate enhancement appears along the edges and corners of the geometries. This is a consequence of the strong field gradients appearing along the graphene edge. 1,44 Second, for all geometries, the maximum quadrupolar rate is two orders of magnitude smaller than the maximum dipolar rate. That two-order magnitude difference compares fairly with the rates comparison obtained in Ref. [5] for a H-like emitter close to a non-structured two-dimensional material supporting plasmons confined with a factor η 0 ≈ 35 − 50 (corresponding to doping between 0.7 and 1 eV). The four-order magnitude difference between the dipolar and octupolar transition rates is also in agreement with the literature.
With the graphene nanoislands, we break the in-plane translational symmetry and the conventional dominance of the dipolar transition rate over the quadrupolar transition rate.
From the spatial maps, we observe that the maxima of the quadrupolar rate do not coincide with the maxima of the dipolar rate: by moving the emitter, one can find a position where the quadrupolar rate dominates the dipolar rate, breaking the conventional selection rules (see Section 3.3).
Note that the z-oriented emitter (out-of-plane direction) shows stronger rate enhancement, but the dipolar, quadrupolar and octupolar transitions are all three enhanced at the

Graphene tunability
In a spontaneous emission process the emission wavelength is fixed via the considered transition. Hence, as the frequency of the source is not a variable, a tuning knob is offered by the environment, e.g. the permittivity of the scatterer. The considered mode expansion is particularly well suited for this context as the permittivity is the eigenvalue of the problem. the following section, we show a particular doping of graphene where the quadrupolar rate dominates the dipolar transition rate, consequently breaking the conventional selection rules.

Local breakdown of conventional selection rules
At particular positions of the emitter, the quadrupolar transition rate overcomes the dipolar transition rate. This breakdown occurs at ultra-strong plasmon confinement (η 0 > 500) for planar two-dimensional materials, 5 which is experimentally achievable with graphene, but at the cost of considerable absorption losses. 13 The shape of the graphene nanoislands provides another degree of freedom to mold the field profile and break the selection rules.
We focus on the triangular graphene nanoisland of 50 nm side length, for which we computed the normalized dipolar and quadrupolar rates of an emitter 5 nm above its surface ( Figure 1). In Figure 3a, we plot the maximum of the ratio Γ Q /Γ D over all positions of the emitter, as a function of the graphene doping. This shows that the quadrupolar transition rate can be up to 100 times stronger than the dipolar transition rate at particular positions, breaking locally the conventional selection rules (the value is converged for 40 modes, as shown in the SI).
In Figures 3b and c, for the x-and y-oriented emitter, respectively, we observe enhancement where the field demonstrates strong gradients, i.e. at the corner of the triangle or along the edge. On the contrary, as observed in Figure 3a for the z-oriented emitter, the dipolar rate always dominates the quadrupolar one (the maximum rate enhancement of each order appear at the same position, as discussed in Section 3.1).
Note that the maximum is not a consequence of an inhibited dipolar transition: the quadrupolar rate is strongly enhanced. For example, for an x-oriented emitter (Figures 3b) at the left corner of the triangle, the dipolar transition remains enhanced (Γ D /Γ D,0 = 0.51×10 3 ), but its rate is weaker than the quadrupolar rate, which is 5.3×10 4 times the dipolar transition in free space.
Other areas further away from graphene seem to demonstrate a strong quadrupolar enhancement (for example position (-9;35) nm in Figure 3b). However, these are regions where the dipolar transition is poorly enhanced (Γ D /Γ D,0 = 16), as well as the quadrupolar transition (Γ Q /Γ D,0 = 161).

Conclusions and perspectives
We develop a numerical method that allows the computation of the transition rates of an emitter in the vicinity of any scatterer in the weak coupling regime, with the prerequisite knowledge of the wavefunctions of the excited and ground states of the emitter. With this method two important parameters are easily variable: the position of the emitter and the permittivity of the material, the latter being easily tunable for graphene. By optimizing both parameters we demonstrate a breakdown of the selection rules, with the quadrupolar transition rate 100 times stronger than the dipolar transition rate for a H-like emitter in the vicinity of a triangular graphene nanoisland. These results uncover interesting perspectives for applications in spectroscopy, photochemistry and quantum technologies.
In this paper, we focus on a particular photonic structure, i.e. a scatterer in free space. For future exploration, GENOME also allows the determination of the Green's function for more complex structures, such as a scatterer on a substrate, 29 non-uniform scatterers, 46 cluster of scatterers and finite periodic structures. 38 For experimental observation the investigation of these designs is important for the coupling of produced photons to the far-field. As an example, combining the near-field results (e.g. Figure 2) with the far-field out-coupling efficiency of the dominant mode, allows to select the graphene doping necessary to reach sufficient far-field emission. Such an analysis was carried out e.g. for two-photon emission processes near graphene nanoislands. 26 Other structures may be envisaged to enhance the coupling of a plane wave with a quadrupolar transition. 47 In parallel, our method allows for the computation of larger atoms and complex molecules by combining GENOME with time-dependent density functional theory techniques. Hence, controlling the emission rate of quantum dots 3,6 and Rydberg excitons 4 in complex electromagnetic environments is within reach. Furthermore, since the dipolar and quadrupolar rates compete, destructive interference effects can be observed and lead to suppression of particular transition channels, 3,21 leading to diverse quantum applications such as quantum computing, quantum storage, and quantum communication.

Supplementary Information
In this supplementary file, we discuss the method implementation. First, the results' convergence is shown for the different parameters in the integration. A verification of the dipolar transition rate in free-space is also proposed. Second, the integration is compared to full wave simulations and shows good agreement. Third, results are given for the quadrupolar transition rate dominance over the dipolar one for the square graphene nanoisland. Finally, supplementary figures complete Figure 1 of the main text.

S1 Convergence of the results
In this section, we show the convergence of the integration, with a cut-off set for 1% relative error. We implemented a spherical integration over the radial component and the azimuthal and polar angles. The wavefunctions are those of the hydrogen atom, and the field is obtained through an adjusted eigenfrequency study in COMSOL Multiphysics (see Ref. 38  The convergence is checked as a function of the number of subdivisions in the angular direction and in the radial direction, see Figure S1. A decent convergence is reached with 100 angular subdivisions and 50 radial subdivisions. The integration domain was truncated to 81a 0 (with a 0 the Bohr radius), as it reaches convergence for the dipolar transition. Finally, 40 plasmonic modes are necessary to reach a 1% relative error for the transition rate. Note that in Figure S1, the convergence is reached for a graphene doping of 0.7 eV: when the graphene doping decreases below 0.5 eV, more modes are necessary to reach 1% convergence (60 modes). Figure S2 shows that the ratio Γ Q /Γ D is also converged with 40 modes (so Figure  is Im Ḡ 0 (r, r, ω 0 ) = 1k/6π, we have drψ * e (r)∇ r ψ g (r) 2 (S1) The integration leads to Γ 0 = 4.484 × 10 5 s −1 in good agreement with the tabulated experimental values for the transition 6p → 4s (Γ D = 4.456 × 10 5 s −1 ), considering the 1% error on the integration.

S2 Comparison with direct COMSOL simulations
The comparison is conducted for the square graphene nanoisland. We compare the dipole and quadrupole transition rates evaluated with COMSOL direct simulations and with Equation 8 of the main text. In COMSOL, the structure is a square of 50 nm side length and 1 nm

S3 Square graphene nanoisland and selection rules
We focus on the square graphene nanoisland (50 nm side length), for which we computed the dipolar (Γ D ) and quadrupolar (Γ Q ) transition rates of an emitter 5 nm above its surface in Figure 1 in the main text. In Figure S4a rate is relatively weakly enhanced compared to its maximum value (Γ D /Γ D,0 = 2 × 10 6 ) : for example, at the maximum of the E F = 0.72 eV situation, Γ D /Γ D,0 = 1.2 × 10 3 .

S4 Supplementary data: yand z-oriented emitter
To complement Figure 1 of the main text, we plot here the dipolar, quadrupolar and octupolar transition rates of a y ( Figure S5) and a z ( Figure S6) oriented emitter 5 nm above square, triangle and crescent graphene nanoislands.  Figure S6: Dipolar, quadrupolar and octupolar transition rates of a z-oriented emitter close to graphene nanoislands of various geometries. The dipolar (left), quadrupolar (center), and octupolar (right) transition rates as a function of the emitter position 5 nm above the graphene nanoisland: 50 nm side length square (up), 50 nm side length triangle (middle), and 80 nm height crescent (bottom). The geometry is represented with a solid white line, and the rates are normalized by the dipolar emission rate in free space Γ 0 . For the triangular geometry, graphene is E F = 0.98 eV doped, for the square E F = 0.72 eV,and for the crescent E F = 0.88 eV.