Photon Emission Rate Engineering using Graphene Nanodisc Cavities

In this work, we present a systematic study of the plasmon modes in a system of vertically stacked pair of graphene discs. Quasistatic approximation is used to model the eigenmodes of the system. Eigen-response theory is employed to explain the spatial dependence of the coupling between the plasmon modes and a quantum emitter. These results show a good match between the semi-analytical calculation and full-wave simulations. Secondly, we have shown that it is possible to engineer the decay rates of a quantum emitter placed inside and near this cavity, using Fermi level tuning, via gate voltages and variation of emitter location and polarization. We highlighted that by coupling to the bright plasmon mode, the radiative efficiency of the emitter can be enhanced compared to the single graphene disc case, whereas the dark plasmon mode suppresses the radiative efficiency.


Introduction
Technological advances in the field of nanofabrication have provided a powerful tool to tailor light-matter interaction. Metallic nanoparticles support surface plasmon resonances where collective oscillations of electron and photons can result in a localization of electromagnetic fields into subwavelength scales [1]. Apart from localizing the incident plane waves, these plasmon modes strongly modify spontaneous emission properties of quantum emitters [2], such as quantum dots, placed close to them. In particular, radiative decay rates of such fluorescent particles can be tuned, depending on whether the emitter couples to a radiative or non-radiative plasmon mode [3]. This approach of tunable fluorescence quenching and enhancement finds its uses in applications such as molecular imaging [4]. A number of geometries have been explored for such decay rate engineering, for instance metallic planar surface [5], photonic crystals [6] and various collections of metal nanoparticles [3]. In particular, collections of nanoparticles, such as dimers provide an additional parameter, namely, separation and orientation of the individual particles with respect to each other to tune the local electromagnetic density of states. In certain dimer systems with inversion symmetry, symmetric or antisymmetric dipole modes can be excited, even using a plane wave, in accordance with the plasmon hybridization model. The antisymmetric or dark mode has a cancellation of the two induced electric dipole moments hence suppresses far field radiation. The converse is true for the symmetric or bright mode where the electric dipole moments add up constructively. Fig. 1. Geometry for studying decay rate engineering: A) In this geometry, the emitter can excite only one of the two modes, depending on its location and polarization B) In this geometry, the emitter can excite both the dark as well as the bright modes. Hence it is suited to studying comparatively, the effect of these modes on the the decay rate of the quantum emitter.
The numerical values of the disc separation D and the radius R, which were used for the BEM simulation are shown here.
Metals, as plasmonic materials, however suffer from a number of drawbacks, such as high losses [7] and limited tunability of electronic carrier concentration. In recent years, graphene [8] has emerged as a very efficient plasmonic material in the far infrared and terahertz range [9]. Because of its unique bandstructure, electrons in graphene behave as Dirac fermions. A consequence of this is that backscattering of electrons from impurities is forbidden [10], which results in graphene plasmons being much less lossy in the far infrared compared to metals in the visible range. In addition to chemical doping, the carrier density or equivalently, the Fermi level in graphene can be tuned via electrostatic gating [11]. Thus graphene is an excellent candidate for tuning light-matter interaction in this wavelength range.
In this work, we discuss how decay rate can be engineered via plasmon modes in a dimer of vertically stacked graphene nanodiscs. Plasmons in a single Graphene disc have been shown to provide very high total decay rate enhancements [12]. A dimer system of nanodiscs, while still having these advantages, provides a route to engineer radiative decay rates via excitation of dark and bright modes. Such bright dipolar modes have recently been experimentally observed in the case of graphene micro-disc dimers [13].
We will firstly solve for the eigenmodes of the systems using a quasistatic approximation. Secondly, a general recipe using the eigenresponse theory will be provided which can be used to model the spatial and polarization dependence of the local density of states. Some of these results for the lowest dipolar mode will be compared with full-wave boundary element simulations.
Finally we use the example of the lowest dipolar mode to show that the fluorescence quantum yield can be tuned by modifying the Fermi level of the graphene nanodiscs and the possibility of obtaining vacuum Rabi splitting in the cavity.

Methods
In this section, we briefly mention the various analytical and numerical techniques that were employed to arrive at the results in this paper.
The first term in Eq. (1) represents intraband contribution and the remaining terms are contributions of the interband transitions to the total graphene conductivity. Here τ is the electron relaxation time. The relaxation time typically has contributions from 1) scattering from impurities in infinite graphene, 2) coupling to acoustic and optical phonons (hω OPh = 0.2 eV) in graphene and phonon modes of polar substrates and 3) edge scattering in the case of finite nanostructures, such as the one discussed in the present paper( [15] and references therein). In literature, relaxation times as high as 1000 fs have been reported [16,17]. For frequencies larger than the optical phonon frequency of graphene, τ ∼ 50 fs [9]. For the results of this paper, we use the conservative lower end of this range. For all the results concerning excitation of the disc modes via local emitters, we use τ = 50 fs and T = 300 K. For the plane-wave excitation result, a larger τ of 100 fs is used since for smaller values, the extinction peak for the dark mode is too broad to be separated out from the background, dominated by the bright mode. However, both these values we used for τ are on the very conservative side of the range of experimentally measured values.

Simulation of graphene plasmon modes
All the simulations of this paper were performed using a Boundary Element Method (BEM) code, SCUFF-EM suite [18], a free, open-source software implementation of the boundaryelement method that implements specialized algorithms for efficient computation of scattered and absorbed power in scattering problems [19]. For simulations we consider a very small thickness "effective graphene" [20]. This now becomes a 3D structure, whose conductivity is given by dividing the 2D conductivity by the thickness of the effective graphene. This allows us to define a permittivity for effective graphene using Maxwells Equations: where ∆ is the thickness of the effective graphene. Convergence tests were performed with ∆ as a parameter and a value of ∆ = 0.25 nm was chosen as the appropriate thickness for the specific range of frequencies and lengthscales of our problem.

Calculation of decay rates
In the second half of the paper, we will discuss the decay rates of emitters placed close to or inside the nanodisc cavity. An ideal dipole emitter can get rid of its energy through two pathways: 1) radiatively into free space propagating modes or 2) non-radiatively into material absorption. The decay rate into the plasmon mode is mostly dominated by absorption. We calculate these decay rates as follows. The total decay rate is calculated using the scattered electric field at the location of the emitter [21]: where Γ 0 is the decay rate of the emitter, if it were in free space. The radiative decay rate is given by Γ rad /Γ 0 = P rad /P 0 = 1 + P sca /P 0 , where P rad is the power radiated to the far field, P sca is the total scattered power and P 0 is the power radiated by the emitter when placed in free-space. The non-radiative decay rate, which is the dominant contribution from the decay into the plasmon mode is given by Γ abs /Γ 0 = P abs /P 0 ≈ Γ plasmon /Γ 0 , where P abs is the power absorbed in the graphene nanostructure.
Vacuum Rabi splitting calculation: The most common way for describing atom-cavity interaction quantum mechanically is through the Jaynes-Cummings Model (JCM) Hamiltonian. The JCM Hamiltonian, in the rotating wave approximation (RWA) is given by: where ω 0 is resonant frequency of the quantum emitter and ω is the frequency of the plasmon mode. σ + (σ − ) are the atomic raising (lowering) operators and a † (a) are the creation (annihilation) operators for a cavity photon.
For the present problem we use an open quantum system approach, in order to incorporate absorption and radiative decay. Thus the evolution of the density matrix is given by [23]: where κ is the rate of decay of the plasmon mode. κ contains both the radiation as well and absorption mechanisms for broadening the plasmon resonance [24]. However, it is usually dominated by absorption. Γ ′ is the decay rate of the quantum emitter into free space modes, modified by geometrical effects. However, in this paper we assume the Γ ′ to be equal to Γ 0 , the spontaneous emission rate in free space, in accordance with Wigner-Weisskopf theory. The system density operator evolves according to Eq. (5). In the single excitation manifold, only the states {|g |1 , |e |0 , |g |0 } need to be retained. Here |g and |e are the ground and excited states of the atom and |0 and |1 denote the number of photons in the cavity mode. It can then be shown [25] that for Rabi oscillations to exist, on resonance (∆ = 0), one needs the condition |g/(κ − Γ ′ )| > 1/2.
In the JCM, the coupling strength g is determined by the details of the cavity field mode and the atomic dipole matrix element. For our purpose, we determine g classically, using the limit of a low finesse cavity. In this limit, g satisfies the following equation: where ∆ is the detuning between the resonant frequency of the plasmon mode and that of the quantum emitter. For the present work, the typical spontaneous emission rate of the emitter is much smaller than the cavity line-width. Thus, Eq. (6) suggests that on resonance, Γ tot /κ = (g/κ) 2 . Hence the g factor can be determined. This expression also points out that Rabi oscillations should exist when Γ tot /κ > 1/4.

Calculation of the eigen-modes in the quasistatic limit
Firstly, we discuss the mathematical formulation of the eigenvalue equation in terms of an electrostatic potential. This is essentially a solution of the Laplace equation for the disc geometry.
This section is divided into two subsections. In the first subsection, we repeat the derivation for the case of single disc, which had been worked out by Fetter [26] in 1986. In the second subsection, we formulate the eigenvalue problem for the case of a stack of two discs.
In the second section, we discuss a numerical framework to solve the eigenvalue equations we obtained in the previous section. In the third section, we summarize the results of the calculation, by providing details of the eigen-mode plots and a comparison to full wave boundary element simulations.

Mathematical framework for a stacked dimer of discs
The general strategy is to solve the Poisson equation for the electrostatic potential Φ(r), with the surface charge boundary condition due to the graphene discs. The surface charge density itself can be related to the electrostatic potential via the continuity equation and the surface conductivity of graphene. This leads to an eigenvalue equation with Φ on both sides. Subsequently, numerical techniques are used to solve this eigenvalue problem to get the resonant frequencies as well as the potential profile. This potential can then be used to calculate various other quantities of interest such as surface polarization and surface current density.
The calculation for the single disc case, was carried out in the Fetter's paper [26] on magnetoplasmons in disk geometries of 2DEG. For the sake of comparison the notation from a recent paper [27] on edge plasmons in a single graphene disc, has been used here.
It should also be noted that for a single disc in the quasistatic regime, closed form solution is possible [28]. However, we use a numerical approach here which can be easily extended to stacks consisting of arbitrary number of discs, where closed form solution becomes very cumbersome.
The geometry consists of two identical discs each of radius R, stacked vertically with a separation D in between (see Fig. 1). The location of the discs in our chosen coordinate system is z = ±D/2. The approach that will be presented here can easily incorporate the case where the two discs are non-identical. However, for the sake of clarity for our specific case, we will only consider identical discs for now.
First note that because of circular symmetry, the potential can be expressed as Φ(r) = Φ(r, z)e iLφ , in cylindrical coordinates.
We follow a two step procedure to get to the eigenvalue equation: • Express the surface potential Φ(r, z = ±D/2) in terms of the surface charge density σ b , using the Laplace equation and the normal electric field boundary condition • Express surface charge density σ b in terms of the surface potential Φ(r, z = ±D/2), using the continuity equation and the current-field relation Expressing Φ in terms of σ b : The Poisson equation in this case is given by: where σ b is the surface charge density (and not the surface conductivity, which is represented by σ ). One way of solving such problems is to write the general form of the solution in the regions on either side of the the boundaries and then match the boundary conditions. We will use this approach.
Thus we write the solution for the Laplace equation for z = 0 and then use the boundary condition for the normal electric field. To be specific, these equations are given below: where ε m is the relative permittivity of the medium in between the discs. Note that the e iLφ dependence was suppressed in the boundary condition equation. (Note that there is another boundary condition which is the continuity of the potential across z = ±D/2.) Now let us express Φ in terms of its Hankel transform component: where the Hankel transform is only taken in the radial coordinate of the cylindrical system. Now we substitute Eq. (11) into Eq. (8). After some manipulation and using Bessel's differential equation, we obtain the following simplified form: For z = ±D/2, Eq. (8) holds for the potential Φ(r) in real coordinates. Equivalently, for z = ±D/2, Eq. (12) holds for the Hankel transformed potential. We can write down the form of the solution in the three different regions as follows: There are four unknowns We also have four equations, two for the continuity of the potential across the discs and the other two for the normal electric field boundary condition.
It is quite straightforward to solve for the general case of different relative permittivities. In the following, we choose ε u = ε m = ε d = ε for simplicity. Solving the above linear system of equations, we get the solution for the Hankel-transformed potential on the discs: Now we go to real space, by taking the inverse Hankel transform on each side of the above equation. For brevity, we denote the Hankel transform operator asĤ(p; r ′ ) = ∞ 0 dr r ′ J L (pr ′ ) and its inverse operator asĤ −1 (r; p) = ∞ 0 d p pJ L (pr). With this notation, we can express the real space solution as: Expressing σ b in terms of Φ: There are two equations that we need to express σ b in terms of Φ. One is the continuity equation for surface current density and the other is the relation between surface current density and the electric field. These equations are given below: Now using the relation E || = −∇ || (Φ(r, z = 0)e ıLφ ), we arrive at the relation: As the last step, we move to normalized coordinates x → r/R and express the final equation in operator form: where η = σ (ω)/ıωε 0 εR. The reader is reminded that hereÎ uu L (x; x ′ ) andÎ dd L (x; x ′ ) are associated with the on-site term for the upper and lower disc respectively, whereas the offdiagonal terms represent the interaction between the two discs. In terms of normalized coordinates, Φ(r, z = ±D/2) → φ u,d (x)e ıLφ . Eq. (20) is an eigenvalue problem in the parameter η = σ (ω)/(iωε 0 εR), which can be related to the resonant frequencies of the modes. Note that in the normalized coordinates, the exponential term in the off-diagonal kernel depends on the ratio D/R instead of just D.
The solution to Eq. (20), will give us the resonant frequencies and the mode-profiles of the plasmons in stacked dimer of graphene discs.
It should be noted that this kind of approach is easily extensible to more than two discs or discs with different radii or surface conductivities.
We solve the eigenvalue problems for the single disc and the stacked dimer of discs case using the standard method of polynomial expansion. The results for the eigen-frequencies are shown in Fig. 2.

Comparison with full wave simulation
We compared the resonant frequencies obtained using the quasi-static solution to those obtained using a full-wave boundary element simulation (BEM). For this comparison we only choose the L = 1, n = 1 mode since that is the mode that we will be concerned with in the rest of the paper, when talking about photon emission rate engineering. It should also be noted that for the simulations, we use a realistic absorption in the graphene conductivity. The comparison is presented in Fig.2. Figure 2 suggests that there is a good overall match between the resonant frequencies found from the BEM and the quasistatic result. The resonant frequencies in the simulation were obtained from the LDOS spectrum.
Two features in Fig. 2 are worth highlighting. Firstly, the resonant frequencies of both the modes increase with E F . This is due to the fact that increasing E F results in an increase in the carrier density, which in turn causes an increase in the restoring force. This explanation is similar to how the plasma frequency in noble metals increases with carrier concentration. Secondly, the frequency splitting ω B − ω D increases with E F . This is due to the fact that at higher E F , the plasmon modes of individual discs are more leaky. This results in the interaction between the two discs being even stronger, resulting in a larger splitting.

Eigen-response theoretic framework and calculation of overlap
To understand dependence of the spatial and polarization dependence of the LDOS, we resort to an eigen-response theory [29]. In the present case of a system with a symmetric and an antisymmetric mode, the polarization density can be expressed as: where α A,L and α S,L are the eigen-polarizabilities of the antisymmetric and symmetric modes and P A,L and P S,L are the eigenmodes. The excitability of the modes is related to the overlap terms, which are the inner products of the mode profile and the excitation.
Relation between overlap and LDOS: The spatial dependence of the LDOS is contained in the overlap term, since the eigen-polarizability is usually only frequency dependent. In general the projected LDOS in Eq. (3) can be written as: where |i is a quantity proportional to the surface polarization, for each mode. In the following section we present the calculation of the surface polarization, which will help us calculate these overlap terms. Calculation of surface polarization: From our quasistatic approach, we determined the surface potential on both the discs. Using the surface potential, it is straightforward to obtain the surface polarization. The surface polarization P s can be related to the potential on the discs in the following way. In the absence of magnetization, the surface current J s can be related to P s as follows: J s can also related to the electrostatic potential using the relation: Thus we have the relation between P s and Φ: The potential Φ obtained by solving the eigenvalue equation as mentioned in an earlier note, can be plugged in Eq. (25) to calculate P s and subsequently the overlap terms. We will consider the source dipole generating a field which excite various infinitesimal dipoles on the disc surface. For this purpose, we will need the Green's tensor which is defined as [21]:Ĝ In order to keep the mathematical framework completely general, we consider the geometry shown in Fig. 1. For this geometry, r 0 = X dx + Y dŷ + Z dẑ and the location of the infinitesimal dipoles r = xx + yŷ + zẑ. The overlap terms for the case of a single disc are shown in Fig. 3.

Decay rate engineering
Single graphene nanodiscs, have been shown to provide very high Purcell factors [12]. Such high enhancement factors are possible due to the very small plasmon mode volume, which is a general characteristic of graphene films and nanostructures in the far infrared and terahertz range [22]. Using a vertical nanodisc dimer cavity, should provide an additional degree of freedom for engineering the decay rates. For instance, other than applying gate voltage or changing the radius of the discs, there is now an additional parameter, which is the separation of the discs, that can be used to tune the resonances [30]. The dark dipolar plasmon mode only weakly couples to plane waves hence its excitation using plane wave is not very efficient. However, if we use local emitters, such as quantum dots, it is indeed possible to excite the dark mode very efficiently [30]. The coupling of quantum emitters to such plasmon modes is reflected in the modification of the decay rates of the former. Depending on the nature of the plasmon mode, radiative decay of an emitter can be suppressed (quenching) or enhanced. In the case of graphene nanodisc dimer cavities, one can easily modify the radiative or non-radiative processes, by using gate voltage or disc separation to tune the dark or the bright mode to be resonant to the quantum dot. In this context, location and polarization of the emitter is another variable [31], which will be discussed in this work.
We consider two situations (Fig. 1) for studying the coupling between the emitter and the plasmon modes. Firstly, we study the case when the emitter placed inside the cavity. If the Fig. 3. Semi-analytically calculated overlap terms for a Single Disc: The emitter is located at a vertical distance of 15nm above the disc and moves along Y d = 0. The colors correspond to different polarizations of the emitter:x(-),ŷ(-) andẑ (-) emitters are located at the inversion plane then it allows us to excite either the dark or bright plasmon mode, depending on the emitter polarization and position. We will use the eigenresponse theory to calculate the overlap terms as a function of the position and polarization of the emitter. We will also present a comparison with LDOS for the lowest dipolar mode calculated via BEM simulation and show that the shape of the LDOS spectrum can be well explained by the calculated overlap terms. Secondly, we place the emitter outside the cavity, resulting in it being able to excite both the modes, for the same location and orientation of the emitter's dipole moment vector. This enables us to directly compare the bright and dark modes, in terms of their efficiency in suppressing or enhancing the decay pathways and quantum efficiency of the emitter, under similar excitation conditions. To avoid repetition, in this section we will only present simulation results for the radiative decay rates.

Case A: Emitters inside the stacked disc cavity
In principle, one can calculate the decay rates for all emitter positions inside the cavity. A randomly oriented emitter should then couple to both the dark as well as a bright mode. However, if the quantum dots are placed inside the cavity, on the inversion symmetry plane parallel to the plane of the disks, then depending of the alignment of the dipole matrix element, the quantum dot will only be able to couple to either the dark or bright mode, but not both. Hence in order to understand the physics of the problem, this is a convenient choice of emitter location.
Position and polarization dependence of the LDOS: Having worked out the matrix elements for the case of the single disc, we now move on to consider the cases of the dimer of discs. Because of mirror symmetry we can predict that there are even and odd modes. In literature, these are often called the bright and dark modes respectively. If we label the upper disc by U and the lower disc by L, then the following relations hold: • Bright mode: dp U = dp L , z U = −z L = d/2 • Dark mode: dp U = −dp L , z U = −z L = d/2 Now to calculate the total overlap for the two modes, we add the contributions from the two discs, taking into account the appropriate sign changes as mentioned in the relations given above: Bright Mode: Since the sign of the infinitesimal dipole moment does not change, we only need to consider the sign change in the z coordinate of the two discs. This results in the total sum of the overlap term for the z = ±d/2 giving a zero for the z−polarization of the emitter. The other two terms for the x and y polarizations survive and are basically twice of the corresponding value for the single disc case. Dark Mode: In this case, the sign of the infinitesimal dipole moment does change, and so does the sign change in the z coordinate of the two discs. This results in the total sum of the overlap term for the z = ±d/2 giving a zero for the x− and y−polarizations of the emitter. The only nonzero term is the one for the z−polarizations and is just twice of the corresponding value for the single disc case. The overlap terms are calculated by evaluating integrals of the form dp * mode · E exc , as discussed earlier. The calculated overlap terms are presented in Fig. 4 for the bright mode and Fig.  5 for the dark mode.
To verify our approach, we perform BEM simulations for a dimer of graphene discs, each 100 nm in diameter and separated vertically by 30 nm. The frequency range for the simulation is chosen to highlight the contribution of the lowest dipolar mode. The semi-analytical calculations using the eigen-response theory gives the dips in the LDOS spectrum ( Fig. 4 for the Similarly for the bright mode, there is a contribution from the z-polarized emitter. This effect is due contributions from neighbouring resonances. in fact this effect can be easily taken in to account in the eigen-response theory, if we include losses in the eigen-polarizabilities. These calculations will be published elsewhere. In the next section, we point out how radiative decay rates can be engineered using these dark and bright modes. In order to compare quantitatively, the effect of the dark and bright modes on the quantum emitter decay rate, we now study a geometry in which the emitter can excite both modes. Any location of the emitter other than the inversion plane is a valid choice. However, for simplicity, we chose to consider the emitter located outside the cavity as shown in Fig. 1. This might also be a convenient scheme, in as far as experimental realization is concerned. We can study all three orientations of the emitter dipole moment as done in the previous section. However, since in this section our main aim is to highlight the engineering of radiative decay rates, we will consider only one one polarization of the emitter, namely the x-polarization. To gain some qualitative insight into the excitability of the modes, we analytically calculate the overlap terms. Based on this calculation (not shown), we find that for an emitter close to the center X d = 0 nm, both the modes are excitable, with highest probability. Radiative decay rate engineering: We then performed simulations to calculate the radiative, non-radiative and total decay rates, when the quantum emitter, located at (X d , 0, Z d ) couples to either of the modes at their respective resonant frequencies. An example spectrum, for X d = 0 nm is shown in Fig. 7(a). In this case, the total decay rate enhancement, which is close to the non-radiative part, is almost the same for the two modes. This is consistent with our analytical calculation of the overlap terms which show that at X d = 0 nm both the modes are equally excitable. However, there is a difference between the radiative decay rate enhancement. This situation results because of partial cancellation of the induced moments on two discs, for the dark mode.
Further, we demonstrate the tunability of the radiative efficiency as a function of carrier concentration. Our calculations in Fig. 7(b) show that the contrast in radiative efficiency increases with increasing E F . The dark mode becomes less, and the bright mode, more radiative, as the carrier concentration increases. As mentioned before, this is because at higher carrier concentrations, the modes of the two discs can interact more strongly.
We also study the dependence of the radiative efficiency, as a function of position X d of the emitter, in Fig. 8(a). Clearly, the bright mode has a higher radiative efficiency compared to the . Z d is fixed such that the emitter is located 15 nm above the closest disc. The bright mode is found to enhance radiative decay rate, but the dark mode does not. b) Vacuum Rabi splitting versus Fermi level: Normalized vacuum Rabi splitting for the dark and bright modes when an x-polarized emitter is located at (0, 0, Z d ) outside the cavity. Z d is fixed such that the emitter is located 15 nm above the closest disc. The dashed line represents g/κ = 1/2 below which splitting will not be observed. (κ ≈10 meV).
dark mode, for various locations of the emitter. We find that the radiative efficiency drops from a maximum at X d = 0 nm to a local minimum, as the emitter approaches a certain horizontal distance ( X d = 40 nm in the specific case of Fig. 8(a)) and then rises again. This behaviour can be qualitatively explained by looking at the overlap term (not shown), which predicts that the mode excitability being maximum at X d = 0 nm, drops to zero at a certain X d to subsequently rise again. Similar pattern in observed for the dark mode.
In the analysis provided here, one must note that the dominant decay pathway for the emitter in this geometry is still non-radiative. Hence the overall quantum efficiency when the emitter couples to either the dark or the bright mode, is rather small for both cases.
Strong coupling regime: Strong coupling regime in a system composed of a single graphene disc and a quantum emitter was predicted in [12]. Here we briefly mention that the same can be obtained in the dimer system as well. However, so far we have not carried out a complete study in this direction.
In our classical electrodynamic simulations carried out for the graphene nanodisc cavity, we obtain Γ tot ≈ Γ plasmon > κ on resonance. This condition ensures the existence of coherent coupling between the plasmon and the emitter and the possibility of observing vacuum Rabi splitting [23].
Using the approach mentioned in the Methods section, we calculated g factors for different Fermi-level energies, ranging from 0.2 − 0.8 eV, for different positions of the emitter both inside and outside the cavity. Our calculations for the normalized Rabi splitting are shown in Fig. 8(b), for an example case of the emitter placed outside the cavity (same geometr y as considered earlier in this section). We find that for a range of values of the doping level, we do obtain g/κ > 1/2 and we predict Rabi splitting values g of up to 10 meV at room temperature. Here for calculating the absolute value of g, we have used Γ 0 ≈ 5 × 10 7 s −1 [12]. We note that it should possible to increase the vacuum Rabi splitting by using higher quality graphene samples, so that τ is larger. Another important trend suggested by Fig. 8(b) is that the value of g/κ decreases with increasing E F . This can be qualitatively understood as resulting from the increasing leakiness of the plasmon mode as the carrier concentration is increased. This results in a larger mode volume V , which is expected to cause a decrease in the decay rate Γ tot . Since κ ∼ 1/τ is almost independent of Fermi Level E F , therefore g/κ ≈ Γ tot /κ decreases with increasing E F .

Conclusion
We have performed a complete study of the plasmon modes of a vertically stacked dimer of graphene discs. The eigenmodes and resonant frequencies of the modes were calculated semianalytically using a quasistatic approximation. We showed a convincing match between full wave BEM simulations and the quasistatic approach. Secondly we explained the position and polarization dependence of the LDOS using the framework of eigen-response theory. In this case also, results were consistent with simulations.
Subsequently, we focused on the dark and bright plasmon modes formed out of dipolar modes of each constituent disc of vertically stacked graphene disc dimer cavity. Due to the different symmetry of these two modes, completely different behaviour is observed in the far field response as well as decay rates. We have shown that it is possible to engineer the decay rates of a quantum emitter placed inside and near this cavity, using Fermi level tuning, via gate voltages and variation of emitter location and polarization. We highlighted that by coupling to the bright plasmon mode, the radiative efficiency of the emitter can be enhanced compared to the single graphene disc case, whereas the dark plasmon mode suppresses the radiative efficiency. Such a system can offer new degrees of freedom for controlling radiative and non-radiative emission properties of quantum emitters.