Combining density functional theory with macroscopic QED for quantum light-matter interactions in 2D materials

A quantitative and predictive theory of quantum light-matter interactions in ultra thin materials involves several fundamental challenges. Any realistic model must simultaneously account for the ultra-confined plasmonic modes and their quantization in the presence of losses, while describing the electronic states from first principles. Herein we develop such a framework by combining density functional theory (DFT) with macroscopic quantum electrodynamics, which we use to show Purcell enhancements reaching 107 for intersubband transitions in few-layer transition metal dichalcogenides sandwiched between graphene and a perfect conductor. The general validity of our methodology allows us to put several common approximation paradigms to quantitative test, namely the dipole-approximation, the use of 1D quantum well model wave functions, and the Fermi’s Golden rule. The analysis shows that the choice of wave functions is of particular importance. Our work lays the foundation for practical ab initio-based quantum treatments of light-matter interactions in realistic nanostructured materials.


Contents
In the case of the van der Waals heterostructure considered in this work, the interpretation of the f,M is as a sum over modes of the electromagnetic environment and final states in the upper subband, as defined by their in-plane momentum q and the dispersion relation of the subband. Noticeably, Eq. (1) describes only a single quanta excitation which can be in the excited state, or as a coherent superposition between many final states and electromagnetic modes. Using this expression for the wavefunction neglects high-order processes such as two-photon emission in addition to the counter-rotating terms. In the structure discussed in the paper, these contributions are weak and therefore neglected, though stronger interaction would force a wavefunction extension into additional terms.
Generally, the time dependent state of the system, Eq.(1), will evolve according to the time dependent Schrödinger equation. To determine its dynamics we write down the three term Hamiltonian of the coupled system: The emitter Hamiltonian can be written as: withb i,f andb † i,f representing the initial state annihilation and creation operators respectively, i, or final states, f , respectively. The field Hamiltonian can be expressed in two forms, symbolically as a sum over modes, M , with operatorsâ † M andâ M for creating and annihilating each mode respectively, or using the MQED notations so that each mode is represented as a dipole creation and annihilation (using the operatorsf † j (r, ω) andf j (r, ω) respectively) located at r, with spatial orientation j = x, y, z and frequency ω. The symbolic sum over modes should thus be understood in the following way, As discussed in the main text, the MQED vector potential in the Weyl gauge is given in terms of the dipole excitations of the combined light-matter system: 2 Finally, the Weyl gauge minimal coupling Hamiltonian becomes, where on one hand we use the electron mass, m, electron charge, e, the electric vector potential,Â(r), and momentum operatorp, and on the other hand the coupling constants, g i,f M , that couples the initial state with each optical excitation and allowed final state.
Inferring that the state has to obey the time dependent Schrodinger equation, the equations of motion for the state coefficients become: where we define ω if M = ω M +ω f −ω i . Using equations 8 and 9 leads to the closed expression for time evolution of the excited state coefficient: In the following we will show that for a stratified medium like the vdW-heterostructure considered in this work, equation 10 can be written in terms of a interaction kernel, K(q, ω), in the following way: wherehω if q = ε f q − ε i is the dispersive intersubband transition energy. As discussed in the main text, the kernel K(q, ω) is a function of frequency, ω, and in-plane momentum, q, and it holds all of the information about the coupling to the electromagnetic environment. The integral kernel is made up of three different contributions, arising from theÂ ·p and ∇ ·Â terms in the interaction Hamiltonian and their cross terms respectively. This follows from the expansion of the matrix elements in equation 10: where we have used the expression for the interaction Hamiltonian in Eq. 7 and further that after inserting the MQED vector potential, the modes of the electromagnetic field are labelled by ω and thus ω if M = ω − ω if , wherehω if is the energy difference between the initial and final state. To make the origin of the different terms clear, we have dubbed these KÂ ·p , Re[K cross ] and K ∇·Â respectively. Importantly, the two last terms only arise when the dielectric tensor of the material is anisotrpic. In an isotropic material, the condition ∇ · ( ↔ Â ) = 0 would mean that ∇ · (Â) = 0. We can therefore think of these extra terms as an anisotropic correction. For the TMDs considered in this work, the dielectric response is uniaxial and these terms therefore need to be taken into account. In the following, we will derive the contribution from each of these three terms to the interaction kernel.

TheÂ ·p contribution
We are first going to consider the dominantÂ ·p contribution to the interaction kernel, Inserting the MQED expansion of the vector potential, equation 6, into equation 14 one finds: where we have used the commutation relation f i (r, ω),f j (r , ω ) = δ(ω − ω )δ(r − r )δ ij 2 along with the Green's function identity: 2,3 Having arrived at Eq. 15, we next insert the explicit DGF and electronic wave functions. The subband states in the TMD are Bloch states and can be expressed as: where k and ρ are the in-plane momentum and location respectively, and A is the crystal area. φ j (r) inherits the in-plane periodicity of the lattice such that for any real space lattice . This leads to the following normalization condition for the meaning that we have to require, where, u.c. , u.c. and A u.c. denote the unit cell summation, the unit cell integral and the unit cell area respectively.
Combining the DGF from the methods section of the main text with the wave functions from equation 17, we find that: Assuming that the phase factor, e −i(q−k f −k f )·ρ , varies slowly inside of a single unit cell leads to conservation of in-plane momentum: This in turn means that equation 20 can be recast as: We next to turn to the anisotropic correction terms. In order to analyze these terms, we need the derivative of the MQED vector potential in Eq. 6, where we have used Einstein summation convention in the above. We note again here that in an isotropic media, it can be shown for the DGF discussed in the main text that In that case only the KÂ ·p term would thus survive. However, because of the uniaxial nature of the TMD's dielectric tensor we need to explicitly consider the extra terms.
By substituting the equation above into the expression of K ∇Â , we find where the symbol ∂ corresponds to the partial derivative with respect to the r coordinates. In order to continue we shall use two Green's function identities. First, we use ( , since the derivative is on a real parameter. Second, we use the identity Therefore, we find the compact formula, The next step is to insert the specific expression for the DGF and the wave functions. This procedure is very similar to the procedure for theÂ ·p term and in the end one ends up with the following result: By the same argument as for theÂ ·p term above, we can write the contribution to the integral kernel, K(q, ω), in equation 11 as: As expected, the kernel contribution in equation 29 is equal to zero in an isotropic medium where ⊥ = . We further note that the contribution is generally positive.
In order to find the final term of the kernel, K cross , we can use a similar procedure as was done above. Consequently, When substituting in the Green's function and the wave functions, we now find: Finally, this means that we can write the final contribution to the integral kernel as, Again, as expected, the kernel contribution in equation 32 is equal to zero in an isotropic medium where ⊥ = . We further note that this term can be negative, though the overall sum KÂ ·p + 2Re[K cross ] + K ∇·Â must be larger than 0.

Common approximations
The kernel function, K(q, ω) and equation 10 fully describe the interaction of the system with its environment and the corresponding dynamics. The work is therefore in principle done at this point. However, we now take few limits to compare with commonly used approximations.
Assuming that the emission takes place in the weak coupling regime, it is fair to assume that C i (t ) is slowly varying relative to the phase factor. We can therefore replace C i (t ) with C i (t) in equation 11 and extend the upper limit of the time integral to infinity (the Markov approximation). This leads to the following 1st order approximate expression for the rate: To investigate the effect of the nonlocality of the matrix elements and transition en- and the K ∇·Â terms. Further, in this local approximation, we do not consider the dispersion of the transition energy or the nonlocality of the matrix elements, that is ω if q = ω if q=0 . As a result, the coupling kernel becomes: where |f q=0 denotes the final state that would result from a vertical transition between the subbands.
The final approximation is what we denote as the sheet dipole approximation. In this approximation we again assume that the cosh(k z z) and sinh(k z z) factors can be replaced by their values in the middle of the slab. However, we now retain the dispersion of the transition energy and the nonlocality of the matrix elements. The sheet dipole approximation is therefore equivalent to calculations with the following kernel: When the dominant contributions to the decay rate come from excitations whose q values fulfill the condition q 1/d (d is the QW width), the coupling should be well described within the dipole approximation. Within the dipole approximation, the decay rates only depend on the wave functions of the subbands via their effective transition dipole moment.
Specifically, the rate depends proportionally on the absolute square of the effective dipole moments, |µ if | 2 . In general, the transition dipole moment is a nonlocal quantity and thus it has a dependence on the value of q. However, if all of the relevant excitations live at relatively small q-values, it should provide a decent approximation to say that the decay rates are directly proportional to the absolute square of the transition dipole moments at That means that if the we rescale the rates from the 1D wave function approximations by the ratio between the absolute square of their transition dipole moment and the absolute square of the transition dipole moment from the full 3D wave function approximations, the rates should agree, as long as q 1/d and thus |µ if q | 2 ≈ |µ if q=0 | 2 . In the following we thus investigate how rescaling the rates calculated within the 1D wave function approximations according to,  we observe that the rescaling of the Particle in a box states slightly overshoots the rates calculated using the 3D DFT wave functions, but that the qualitative trends agree well.
The reason for the overshoot is down to nonlocality of the transition dipole moment near where the plasmon intersects the transition energy dispersion. This nonlocality of the dipole moment near the plasmon is far more pronounced for the 2 layer slabs than for the thicker samples, because, as can be seen from figure 3 in the main text, the larger transition energies in the 2 layer slabs means that the plasmon intersects the transition energy dispersion at far larger q than what is the case for the thicker samples. This nonlocality doesn't manifest as much in the difference between the 1D and 3D DFT wave functions, because the difference in their effective dipoles in the out-of-plane direction is a constant factor independent of q.
This happens because the 1D DFT wave functions retain the out-of-plane character of the 3D DFT states, and therefore they differ from the full 3D states only in that they miss the in-plane overlap of the states.
Supplementary Figure 1: The scaling of the dipole approximated result to fit the full 3D DFT wavefunction transition. This figure shows the scaling of the first order decay rates by the square of the ratio between the vertical dipole transitions matrix elements, | f q=0 |p z |i | 2 , between the different wave function approximations.

Supplementary note 3: Reflection from a uniaxial material
In this section we will derive the Fresnel reflection coefficient off of the graphene sheet for a wave incident on the graphene sheet from a uniaxial medium. The incident wave with magnitude, E 0 , is defined as, where q is the wave vector in the x-direction and k z = x k 2 0 − x z q 2 is the out-of-plane wave number. With these definitions, the incident wave satisfies ∇ · · E = 0 and Defining the reflection-, r g , and transmission coefficient, t g , the expressions for the reflected and transmitted fields become: where the z-component of the wave vector for the transmitted field is k z,t = (k 2 0 − q 2 ) 1/2 .
Restricting the discussion to non-magnetic materials, µ r = 1, the magnetic field, H, and the magnetic flux density, B are trivially connected via H = 1 µ 0 B. This means that the magnetic field can be calculated directly from the electric field via Faraday's law, The next step is to use the electromagnetic boundary conditions. Specifically, the boundary conditions on the parallel parts of the fields. Noting that the unit normal of the interface isẑ and that the surface current can be related to the transmitted field via Ohm's law, where σ(q, ω) is the nonlocal conductivity of graphene as described in Ref., 4 the boundary conditions at the interface becomes: Inserting the electric fields into equation 45 leads to: and inserting the magnetic fields into equation 46 gives: Equations 47 and 48 represent two linear equations with two unknowns and they can thus be solved for the generalized reflection coefficient of the interface, r g . If we consider the quasi-static limit, which is well justified here since all of the relevant excitations are strongly confined, we can take k z ≈ i ⊥ q. Using this, one finds the following expression for the Fresnel reflection coefficient for a wave incident on a graphene sheet from a uniaxial medium: Supplementary note 4: Calculating the diamagneticÂ 2 term.
In this work, we neglect the A 2 term in the minimal coupling Hamiltonian. To provide at quantitative justification for this assumption, we here calculate it's value for the graphene-TMD-mirror cavity system and compare it to the other characteristic energies in the system.
Including the A 2 term would lead to the following addition to the Hamiltonian: To calculate the above term we need to use the MQED expansion of the vector potential from equation 2 in the main text. For the configuration considered in this work, the relevant terms in the above contribution comes when a field creation and annihilation operator are combined. In such cases, the terms lead to an interaction of the electronic state with itself via the electromagnetic field causing an energy shift of the electronic states. For the excited state this energy shift would take the form: where Tr Im ↔ G denotes the spatial trace of the DGF. Using equation 18 from the main text, one finds that: Inserting equation 53 into 52, while assuming azimuthal symmetry of the dispersion of the transition energy, leads to the final expression for the correction from the A 2 term: e, 0| e 2 2mÂ (r) 2 |e, 0 = αh 2 c π mA u.c dω dq q ω 2 Im r g e −2kzd − 1 r g e −2kzd + 1 d 3 r|ψ e (r)| 2 k z sinh 2 (k z z) + The fundamental equation that we start with is the integro-differential Eq. (5) in the main text that describes the full time dynamics of the excited state. We then describe this equation with new initial state amplitudesC i (t) = C i (t)e iδωt , to finḋ From here, the Markov approximation is used, which assumes thatC i is slowly varying in time compared to the rest of the time integrand, dω dqK(ω, q)e −i(ω−ω if (q)−δω)t . As a result, the time integral in the equation above can be taken to infinity. When comparing the real and imaginary parts of the resulted equation, we find the expression to the spontaneous emission rate which is equation 33 (Eq. (8) in the main text), and the energy renormalization as: In order to access the time dynamics of the excited state probability we have to numerically integrate equation (5) from the main text. However, equation (5) is a stiff differential equation, as it contains dynamic on two different time scales defined by the Rabi-oscillations and the overall decay. In order to ensure numerical stability of the numerical integration scheme it is therefore necessary to use a relatively small time spacing. In order to ensure stability, we employ a time spacing of dt = 10 −5 ps which is orders of magnitude smaller than the fastest time scale of the dynamics. Furthermore, we propagate everything 10 6 steps out to 10ps, which is long enough to ensure that all dynamics is well decayed for all configurations considered in this work.
In figures 4 and 5 we show the full time dependence of the excited state probabilities for a 2 and 5 layer WSe 2 slab at the graphene Fermi level resulting in the largest Purcell en-Another important parameter to converge is the number of grid points in the (q, ω) grid on which we compute the integration kernel, K(q, ω). Figure 6(a) shows the decay rates as a function of the number of points in a square (N, N ) grid of (q, ω) points, and figure 6(b) shows the deviation from between the different grids, defined relative to the rate at a grid resolution of N = 10000. The example given in the figure is for a 2 layer slab since this is the system that converges slowest with respect to the number of points. As can be seen in the figure, the decay rates are converged to within less than 0.2 % at N = 7000 and so we use this grid resolution for the calculations in paper. Further, we note that in principle the integration over both q and ω runs from zero to infinity. However, in practice we can cut the q-integration off at the q hz due to the conservation of in-plane momentum. The cut-off frequency is more arbitrary, but we find that a cut-off of ω cut = 2ω if q=0 gives well converged results.
Finally, for the wave function integrals, we compute the wave function integrals on a real space grid with a 0.1Å grid point spacing, which we find to be sufficient to obtain converged results. Because the electromagnetic excitations in this configuration carry significant momentum, we have to properly take the nonlocality of the matrix elements into account. To this end, we employ 10 q-points between q = 0 and q = q hz , and use the wave functions according to their q. With these parameters, we find that the values of the wave function integrals are converged.