Signatures of Non-Markovianity in Cavity-QED with Color Centers in 2D Materials

Light-matter interactions of defects in two dimensional materials are expected to be profoundly impacted by strong coupling to phonons. In this work, we combine ab initio calculations of a defect in hBN, with a fully quantum mechanical and numerically exact description of a cavity-defect system to elucidate this impact. We show that even at weak light-matter coupling, the dynamical evolution of the cavity-defect system has clear signatures of non-markovian phonon effects, and that linear absorption spectra show the emergence of hybridised light-matter-phonon states in regimes of strong light-matter coupling. We emphasise that our methodology is general, and can be applied to a wide variety of material/defect systems.

Single-photon emission has been observed from a broad range of two dimensional materials (2DM) [1][2][3][4][5][6][7].Of particular interest to quantum technologies is the emission from defect complexes in hexagonal Boron Nitride (hBN) [8][9][10][11][12][13], owing to the materials wide band-gap [14] and stability of the emitters [10,15].Not only does the precise nature of these defect complexes remain under scrutiny [10,11,16], but also their quantum optical properties remain largely unexplored, particularly in regimes of strong light-matter coupling [17,18].Crucially, the optical properties of condensed matter systems cannot be understood through the standard theory of quantum optics, since strong system-environment interactions lead to a breakdown of the underpinning peturbative methods [19,20].Of particular relevance to the optical properties of defect complexes in 2DM are strong electronphonon interactions.Specifically, defect emission spectra typically show sharp and well resolved phonon sidebands [21][22][23][24], a hallmark of vibronic state formation and long lived correlations between electronic and vibrational degrees of freedom.It is therefore crucial that any theory describing the dynamical or optical behaviour of a 2DM defect complex, accurately accounts for electron-phonon interactions.
To this end, we study the cavity quantum electrodynamics (cQED) of defects in 2DM [17,18], focusing on two colour centres in hBN proposed as singlephoton emitters, the C B V N [10,25] and C 2 C N defect complexes [26,27].We accurately account for electronphonon interactions by combining ab initio calculations of the defects, with a fully quantum mechanical and numerically exact description of the dynamics of the defect states interacting with a single mode optical cavity.The former provides an atomistic description of the phonon modes of the material system.The latter employs the time evolved matrix product operator (TEMPO) algorithm [28], where non-Markovian influences of the phonon environment are encoded within a tensor network.When combined with tensor compression methods, TEMPO provides an efficient approach for calculating the reduced dynamics of the system [29][30][31], which we extend to extract the linear response spectrum of the defect-cavity system.This hybrid approach is inspired by [32] and allows us to maintain the accuracy and insight of first-principles calculations, while enabling coherent dynamics to be captured in previously inaccessible regimes of light-matter coupling, information that is typically inaccessible to first-principle methods.
A central finding of this work is that electron-phonon interactions in the studied defect complexes lead to highly non-Markovian dynamics, even in regimes of moderate light-matter coupling strengths, which are accessible to state-of-the-art experiments.We attribute this behaviour to coherent coupling between the electronic states of the defect and high-quality (Q) factor phonon modes, leading to long-lived correlations between the system and its phonon environment.This provides a mechanism to optically manipulate mechanical modes of a condensed matter system that would otherwise not directly couple to light.Furthermore, from linear response spectra, we see evidence of hybridisation between cavitydefect polaritons and high-Q phonon modes when entering regimes of strong light-matter coupling.These states bare resemblance to recently predicted exciton-photonphonon hybridisation in hBN [33], however in this instance they occur in a continuum of modes coupled to a defect state, and therefore will inherit non-linear features from the localised electronic states.For simplicity we focus on the C B V N defect complex in the manuscript, leaving the discussion of C 2 C N to the supplementary information (SI).The electronic structure of the C B V N defect consists of two separate manifolds [34]: the first consists of the (2) 1 A 1 and (1) 1 B 1 excited states, and a (1) 1 A 1 singlet ground state; the second contains the triplet states (1) 3 B 1 and (2) 3 B 1 [16], separated in energy by ω e = 2 eV [34].By assuming the inter-system crossing due to spin-orbit coupling occurs on a much longer timescale (10 − 100 ns) than cavityenhanced optical transitions (< 1 ns), we restrict our attention to the triplet state manifold, which reduces to a two level system with ground and exited states, |g and |e respectively [35].The electronic transition interacts with a single mode optical cavity as shown schematically in Fig. 1a, with a cavity resonance Ω c and described by a Jaynes-Cummings type interaction: where σ = |g e| is the transition dipole operator, a (a † ) is the annihilation (creation) operator of the cavity mode, and g is the light-matter coupling strength.Here we limit ourselves to the single excitation subspace, spanned by the basis {|g, 0 , |e, 0 , |g, 1 }.The composite cavityemitter system interacts with both vibrational and electromagnetic environments, such that the global Hamiltonian is written as The emitter-cavity system couples to an external electromagnetic environment either through direct emission from the defect into free space, or leakage from the cavity mode.In the rotating wave approximation the interaction to the electromagnetic field takes the form Â † j c j,l +h.c.), where Â1 = σ and Â2 = a, and c j,l is the annihilation operator for the l th mode of the j th electromagnetic environment.In the limit of weak coupling to these fields, we can make a flat spectrum approximation, such that f 1,1 = √ 2πΓ and f 2,1 = √ 2πκ [36,37].Here Γ is the spontaneous emission rate from the TLE into free space, and κ is the cavity leakage rate.Using these approximations and tracing over the external EM fields, we can describe the optical contribution to the system evolution through the superoperator V t = exp(t L S ), where L S ρ = −i[H S , ρ]+ΓL σ ρ+κL a ρ, and we have introduced the Lindblad dissipators L o ρ = oρo † − {o † o, ρ}/2.For a full derivation and discussion of these approximations please refer to the SI.
We take a linear electron-phonon interaction [38] of the form , where b k is the annihilation operator for a phonon mode with wave vector k, with the corresponding free field evolution , and g k is the electron-phonon coupling strength of the k th phonon mode.The electronphonon coupling can be fully characterised in terms of the spectral density is the partial Huang-Rhys factor of mode k.The total Huang-Rhys parameter can then be recovered by integrating the spectral density through S Tot = ∞ 0 dωJ Ph (ω).To determine the partial Huang-Rhys factors, and therefore the spectral density, we treat the electronphonon coupling from first principles, using density functional theory (DFT) to calculate the normal modes of the ground-and excited states of the hBN lattice with the considered defect complex.Importantly, the normal modes of the ground-and excited states of the system can be very different and can even be of different types.To properly account for this, we employ the method outlined by Duschinsky [39] to calculate the partial Huang-Rhys factors.For completeness we also checked the generating function approach [40] and found almost identical results, meaning that the modes in the ground and excited states for the system studied in the present work are in fact very similar.The DFT calculations of the phonons were performed for periodically repeated defect complexes in hBN monolayers using the Vienna Ab Initio Simulation Package (VASP) [41].In order to avoid interactions between the periodic copies of the defects, we used a 9x9x1 supercell for the calculations, which we relax on a 3x3x1 K-point grid, to a force convergence of 10 −3 eV Å−1 using a plane wave cut-off of 700 eV.All calculations were performed with the HSE06 exchange correlation functional, which is essential for an accurate description of the electron-phonon coupling [42].The normal modes and the dynamical matrix was calculated at the Γ-point.Note that we relax the defect in the inplane constrained C 2 V configuration.
The resulting spectral density is shown in Fig. 1(b), where we see multiple sharp peaks present.These peaks correspond to high-Q phonon modes present in hBN; notably from the atomistic simulations, we can assign the dominant contributions to the peak ∼ 125 meV from de-localized defect breathing modes, in which the atoms surrounding the defect oscillate along the dipole direction of the defect.
To understand the influence of these modes will have on the cavity-defect system, we can calculate the bath correlation function shown in Fig. 1(c), where T is the temperature, k B is the Boltzmann constant, and we have assumed the phonon bath to be initially in thermal equilibrium.This correlation function encodes the timescales over which memory effects last between the system and environment [43], for a full derivation we refer the reader to the SI.The sharp peaks in the ab intio spectral density lead to long lived oscillation in the correlation function, which can be attributed to phonon modes with large Q-factors.
To highlight the importance of this structure we can replace the spectral density with a broadened function with the same total Huang-Rhys parameter, as shown by the dashed curve in Fig. 1(b).The corresponding correlation function, shown in Fig. 1(c), has only small oscillations which are rapidly damped.
To study the time evolution of the composite cavityemitter system, we employ the time-evolving matrix product operator (TEMPO) algorithm initially developed by Strathearn et al [28].TEMPO is a powerful method for the study of open quantum systems in strong coupling regimes [29,30], and has been applied to study quantum thermodynamics in the strong coupling regime [44], as well as non-additive phenomena in quantum optics [31], and optimal control [45].The starting point for the TEMPO formalism, alongside other numerically exact path integral methods [46][47][48], is the Trotter decomposition [49], where for a sufficiently small time increment δt, the propagator for the open sys-tem can be factorised such that ).The superoperator V δt is as defined above, and captures the evolution of the system and dissipation through the external electromagnetic fields.The interaction between the system and phonon environment is captured through W δt = exp(δtL B ), and is given by . This partitioning allows one to construct a discretetime influence functional of Feynman-Vernon type [46,47], which captures the influence of the environment to all orders in the interaction strength.The reduced state after k-time-steps can be expressed as: where we have introduced the compound indices α k = (s k , r k ) and β k = (t k , u k ), yielding the density matrix elements ρ α k = r k | ρ |s k , and the system superoperator is given by , the trace over the environmental degrees of freedom can be done analytically [46,47].
Crucially, the influence tensor scales exponentially in the number of time-steps taken [46,47]; while applying a finite time memory approximation can reduce the computational cost of propagating the reduced state of the system out to long times [46,47,50], it limits one to scenarios where key dynamics occur on short timescales.A key insight of Strathearn et al [28] was that the influence tensor may be represented in matrix product operator (MPO) form.This allows one to encode the exponentially growing tensor as a tensor network, and apply tensor compression [51] to reduce the rank of the elements in the network, thereby circumventing exponential scaling.The reduced state of the electronic system is then calculated by contracting the network down after each time-step.The convergence of the TEMPO algorithm is sensitive to taking a sufficiently small timestep δt, and the degree to which the tensors are compressed.Further details of TEMPO and its convergence properties are discussed in the SI.
We now consider the dynamics of the cavity-defect system initialised in the state ρ(0) = |e, 0 e, 0|.Fig. 2 shows the time evolution of the excited state population and the cavity mode occupation for different cavity coupling parameters.We find that across all parameter regimes TEMPO predicts complex oscillations in the cavity occupation.We attribute these oscillations to high-Q phonon modes in the environment, which lead to longlived oscillations in the bath correlation function shown in Fig. 1c, indeed these oscillations follow closely those observed in the bath correlation function.We confirm this by repeating the TEMPO calculations with the effective width spectral density shown in Fig. 1b.Since this spectral density has the same total Huang-Rhys parameter, naïvely we might expect similar dynamics to emerge.However, the resulting dynamics in Fig. 2 show no oscillations at weak light-matter coupling, suggesting it is indeed the structure in the spectral density that has the dominant contribution to non-Markovian behaviour.Interestingly, phonon induced oscillations are present in the cavity dynamics across all parameter regimes studied, but only emerge in the emitter dynamics in the strong light-matter coupling regime.This is a consequence of phonon processes inducing fluctuations to the defects' excited state energy; the cavity is sensitive to these fluctuations through the dipole coupling with the emitter, leading to the observed oscillations in the occupation.However, these effects do become increasingly visible in the emitter dynamics with increasing coupling strength, with the onset of vacuum Rabi oscillations.
We now turn our attention to the linear spectrum of the cavity-defect system.The linear response of an electronic system coupled to a quantized mode can be formulated in two ways, by an external field coupling to the electronic degrees of freedom, or by an external current pumping the cavity mode [52,53].In both cases the driving can be included as a perturbation at the Hamiltonian level, such that H (t) = H + E(t) • µ.Here the system transition operator can take on two values µ = {σ † + σ, a † + a}.The first is the dipole operator of the defect, describing the scattering of light directly off the two level transition, and forms the basis of standard linear response theory [54].The second is the quadrature operator of the cavity field, which induces a static polarisation of the cavity mode and excites real photons into the field [53].E(t) then corresponds to the external field or the external current respectively.
Treating the driving as weak, we can extract the linear response spectrum using density matrix perturbation theory [54].Taking the semi-impulsive limit such that E(t) ∝ e iω D t δ(t), where ω D is the driving frequency, we obtain the absorption spectrum: where we have introduced the first-order response function S (1)  All other parameters are the same as Fig. 2 cavity parameters.To understand the role of phonons in these spectra, we artificially include a scaling parameter to spectral density J Ph (ω) → α HRF J Ph (ω), where α HRF ∈ [0, 1], equivalent to scaling the total Huang-Rhys parameter.Fig. 3a shows absorption spectrum at the onset of the strong coupling regime, where g = κ.
In the absence of phonons (α HRF = 0) we see a clear Rabi splitting, with the separation of the peaks determined by the the light-matter coupling g.As α HRF increases, the separation of these peaks is reduced.We can understand this reduction by appealing to the polaron formalism commonly used to study the behaviour semiconductor quantum dots [20,55]: here, the light-matter coupling strength is reduced by the Frank-Condon factor F = exp(−α HRF S Tot /2), which accounts for the geometric difference of the phonon modes associated to the emitter ground and excited states.Fig. 3d compares the peak separation as a function of the scaling parameter α HRF for the three parameter regimes.For the parameters associated with Fig. 3a, the peak separation follows the Frank-Condon factor F (solid curve).
In regimes of stronger light-matter coupling show in Fig. 3b and c, a more complex picture emerges: signifi-cant structure becomes apparent in the spectra as α HRF increases.Of particular interest is the departure from the well understood polaronic physics seen in Fig. 3a, which is highlighted in Fig. 3d.Here we see that at stronger coupling regimes (g = 50 meV), the renormalisation no longer follows the Frank-Condon factor, and in regime of strong light-matter coupling (g = 100 meV), we in-fact see the splitting increase with α HRF .We attribute this behaviour to a hybridisation of the light-matter polariton and high-Q phonon modes.The resultant state is a tri-partite quasi-particle with characteristics of light, matter, and vibrations.This can be seen most clearly in the upper-polariton of Fig. 3c, where an additional splitting emerges when α HRF → 1.This interpretation is supported by the C 2 C N calculations shown in the SI; this defect complex has a spectral density with little low energy structure, such that at g = 50 meV no hybridisation occurs, and the renormalisation of the Rabi frequency follows closely the Frank-Condon factor F .It is only at higher light-matter coupling strengths (g = 100 meV) for C 2 C N , when the polariton splitting approaches resonance with a high-Q phonon mode that we observe a departure from Frank Condon physics, heralding hybridised polariton-polaron states.
Conclusion-In this letter we have combined atomistic simulations of a defect complex in hBN with TEMPO, a numerically exact and fully quantum mechanical simulation method.Our hybrid approach allows us to study realistic emitters beyond phenomenological and approximate treatments [22], providing a complete description of electron-phonon interactions in condensed matter single photon emitters in optical or plasmonic cavities, with direct relevance to on-going experiments.Furthermore, by considering the cavity quantum electrodynamics of a defect across the weak and strong light-matter coupling regimes, we have shown that strong coupling to high-Q vibrational modes play a significant role in determining the dynamics of the cavity-defect system, even in the weak light-matter coupling regime.At strong lightmatter coupling, the absorption spectrum show clear signatures of hybridisation between the light-matter polaritons and phonon modes inherent to hBN.
The method we present here is general, and not restricted to the material system or specific defect complexes studied here, with potential application to organic polaritons [56].It is worth noting that we have restricted ourselves to low temperatures in the above discussion; at elevated temperatures phonon processes beyond linear electron-phonon coupling become important through mechanisms such as the Jahn-Teller effect [57].

I. MODEL AND DYNAMICAL DESCRIPTION
In this section we will derive the dynamical description of an emitter interacting with electromagnetic and phonon environments used in the manuscript.We start by considering the Hamiltonian of the emitter and environments, given by: The time evolution of the emitter and the single mode the optical cavity is generated by the Jaynes-Cummings Hamiltonian H S = ω e σ † σ + g(σ † a + σa † ) + Ω c a † a, with the parameters defined in the manuscript.In the dipole and rotating wave approximations, the interaction with the electromagnetic environment can be written as [36]: where Â1 = σ and Â2 = a, and c j,l is the annihilation operator for the l th mode of the j th electromagnetic environment and f j,1 the corresponding coupling strength.The free evolution of the fields are then determined by H EM B = j,l ω j,l c † j,l c j,l .The electron-phonon interaction is assumed to be linear [38] and takes the form, where the parameters are defined in the manuscript.In the subsequent sections of this supplement, we will derive the dynamical description of system subject to these interactions, starting with the optical field, before introducing the Time Evolved Matrix Product Operator formalism (TEMPO) [28].

A. Coupling to the electromagnetic field
To derive the impact of the electromagnetic fields, we first consider its influence in isolation from the electron-phonon interactions.By assuming the coupling between the cavity-emitter system and the external electromagnetic is weak, the resulting dynamics may be described by second order Born Markov master equation, which in the interaction picture takes the form [36,43]: where S j,1 = A † j and S j,2 = A j are the system coupling operators given in Eq. 2 corresponding to the j th -optical environment.Similarly C (j) αβ = B j,α (τ )B j,β is the correlation function of the j th -environment, with B j,1 = l f j,l c j,l and B j,2 = B † j,1 .Taking the optical fields to be in their vacuum state, the only non-zero correlation function is C (j) 21 (τ ) = l |f j,l | 2 e iω j,l τ .We can now make two further simplifications: the first is the flat spectrum approximation arXiv:2207.10630v2[quant-ph] 22 Jul 2022 ∀ l, which is valid in the weak coupling regime, and in situations in which the bath spectral density does not vary appreciably over the system energy scales.The second is to take the continuum limit, extending the lower limit of integration to −∞ [36].In combination these approximations simplify the bath correlation functions to δ-functions, such that C (j) 21 (τ ) ≈ |f j | 2 δ(τ )/2π.We thus, end up with a master equation of the form: where we have defined the lindblad dissipator Moving back into the Schrodinger picture, and including system Hamiltonian, we obtain

B. Deriving the influence functional and augmented density operator
To derive the TEMPO algorithm, we follow Refs.[28].We start by considering the Liouville superoperator that generates the evolution of the system L = L 0 + L E , where •] accounts for the dynamics of the phonon bath.The dynamics of this system are generated by the superoperator U(t) = exp(Lt), such that the reduced state of the system is given by ρ S (t) = tr E (U(t)ρ 0 ⊗ τ B ), where ρ 0 is the initial state of the system, and is the Gibbs state of the phonon bath, with Z the corresponding partition function.
We start by discretising the time interval, such that t = kδt for integer k and δt 1.The propagator then becomes U(t) = U N δt , with U δt = exp(Lδt), allowing us to apply the symmetric Trotter splitting such that: )δt .For an initially uncorrelated system and environment, we can decompose the initial state of the system and environment into the system basis, such that χ(t = 0) = r0s0 ρ r0s0 |r 0 s 0 | ⊗ τ B , where ρ s0r0 = r 0 | ρ 0 |s 0 .Applying the Trotterised propagator to the initial density operator yields: where we have introduced the compound indices α i = (r i , s i ) and β i = (t i , u i ).By assuming the interaction term is diagonal in the system basis, the superoperators can be expanded such that Applying this expression iteratively and tracing over the environmental degrees of freedom, we have after k timesteps: where |s k and we have defined the influence functional as ).For an environment in a Gibbs state, the trace in the expression for the influence functional can be done analytically [46,47], and written in product form as, where we have defined the influence tensors: Here λ s are the diagonal elements of the coupling operator, in this case σ † σ, and we have introduced the discretised memory kernel with bath autocorrelation function, It is clear from Eq. 6, that to calculate the state of the system at some time step k, requires the construction of the (2k + 1)-index object: this object is referred to as the augmented density tensor (ADT), and, for arbitrarily small timesteps δt, is an exact representation of the system-environment interaction for the open quantum system of interest.The size of the ADT scales exponentially with number of timesteps taken.To circumvent this scaling, we follow Strathearn et al [28] and represent the the ADT in matrix product operator form, utilising standard tensor compression methods [51] to reduce the rank of the elements in the tensor network.We first start by considering the product of influence tensors in Eq. 11, . We can write this tensor in matrix product operator (MPO) form by defining the higher-rank tensor: where we have promoted the influence tensors to rank 4 tensors using the Kronecker-δ: As discussed in [28] and [29], each time step corresponds to an application of the MPO B, which is represented as a tensor network.Each node in the network consists of the rank-4 influence tensors defined above.We then compress the network after each timestep using a singular value decomposition (SVD) [51] at each node, and make a low rank approximation, whereby singular values less than a threshold value ε C are discarded.This reduces the internal bond dimension of the MPS, providing a controllable way to reduce the size of the tensors.Alongside the size of the timestep dt, the singular values cut-off will be a principle convergence parameter, which will be discussed in Section IV.

II. LINEAR RESPONSE THEORY WITH TEMPO
In this section we extend TEMPO to calculate linear response functions, and subsequently the absorption spectrum.Following the example of Flick et al. [53], we consider a system weakly perturbed by by a driving field, H = H +H D (t), where H D (t) = E(t) • µ.The above Hamiltonian can represent two distinct physical scenarios: probing the emitter degrees of freedom using a classical driving field E(t), or probing the cavity fields with E(t) representing an external current.In the two different cases, the system transition operator take on two values, µ = σ † +σ, or a † +a respectively.The first is the dipole operator of the defect, describing the scattering of light directly off the two level transition, and forms the basis of standard linear response theory [54].The second is the quadrature operator of the cavity field, which induces a static polarisation of the cavity mode and excites real photons into the field [53].Conceptually we can separate these two situations into different modes of driving, where the external field perturbs the cavity mode or the two level transition directly.
Treating the driving as weak, we follow the density matrix perturbation theory formalism outlined by Mukamel [54], and treat H D (t) perturbatively.Moving into the interaction picture with respect to the global Hamiltonian H, the density operator may be written as, where χ represents a state in the full Hilbert space of the open quantum system, and The equation of motion for the state is then governed by the interaction picture von-Neumann equation χI (t) = −i [H D (t), χ I (t)].We can formally write the solution to this equation in terms of the series expansion , where in the Schrödinger picture, the n th -order contribution takes the form We choose the initial state of the open quantum system χ(−∞), to be the equilibrium state such that it does not evolve in time.Using this series expansion, we can consider the consider the effect the peturbative driving has on the polarisation of the field P (t) = tr(µχ(t)) = ∞ n=0 P (n) (t).By substituting in the expression for the driving Hamiltonian H D (t), and making a change of variables for the integrals, we obtain where we have introduced n th -order response function, For this work, we will focus on the linear response of the system, that is, the first order response function: For a defect coupled to phonons and an optical field, the equilibrium state is simply , that is a product state of the system in its ground state, the phonons in Gibbs state τ B , and the electromagnetic field in the vacuum |{0} = l |0 l .The trace in Eq. 18 with the global equilibrium state χ(−∞) means that even in the case of the bosonic quadrature operator µ = a + a † we never leave the single excitation Hilbert space in the case of linear response.We can therefore consider the linear response of the coupled cavity-TLE system within the single excitation subspace without loss of generality.
To obtain the linear spectrum we make a further simplification by working in the semi-impulsive limit, that is, we assume E(t) ≈ E 0 e iω D t δ(t), where E 0 and ω D are respectively the amplitude and frequency of the driving field.In this case, the linear polarisation takes the simple form P (1) (t) = E 0 S (1) (t), and we can link this to the field scattered by the system simply as E Out (t) ∝ iP (1) (t).The general expression for the (heterodyne) linear absorption spectrum can be found as [54]: Therefore, to obtain the linear spectrum, we simply need to calculate tr (µ(t)µ(0)χ(−∞)).This can be done using TEMPO by propagating a system in the unphysical initial state (0) = µ(0) |g, 0 g, 0|, and tracing with relevant driving operator µ at time t.Formally, one would have to propagate the system to t = ∞.However, in practice the important thing to ensure is that the system propagated to equilibrium.

III. AB INITIO SPECTRAL DENSITIES AND C2CNCALCULATIONS
In addition to the C B V N defect, we also investigated the effect of the cavity on the C 2 C N defect complex, which is considered another strong candidate as the source of single photon emission in hBN.In this section, we outline how we extract the spectral densities from the output of the atomistic simulations, before presenting the results for C 2 C N .Note that we shall not give a detailed description of the atomistic simulation, focusing only on the aspects important to the manuscript.

A. Spectral function of the electron phonon coupling
As discussed in the manuscript, the electron phonon coupling of the defect complexes studied are treated from first-principles using density function theory (DFT).Specifically, the normal modes of the lattice are calculated independently for the defect complex in its ground and excited state configuration.Partial Huang Rhys parameters, S k , are then calculated using the Duschinsky method [39], which accounts for the possible geometric differences of phonon modes in the ground and excited state configurations.
We can define the spectral density of the electron phonon coupling in terms of these partial Huang Rhys parameters, J Ph (ν) = k S k δ(ν − ν k ).To account for the finite lifetime of these phonon modes, and to obtain a smooth function for the spectral density, we approximate the δ-function with a Gaussian, such that J Ph (ν) = k S k f (ν − ν k ), where we have introduced the broadening function f , with broadening parameter ς. ς was chosen to 2.5 meV for C B V N and 10 meV for C 2 C N , similar to the values used in Refs.[25] and [16] respectively.Fig. 1 shows the partial Huang-Rhys factors in black, and the broadened spectral function of the electron-phonon coupling in orange, for the C B V N (a) and C 2 C N (b) defect complexes respectively.In contrast to the C B V N spectral density, the phonon environment of C 2 C N has very little structure at low energies, with a dominant contribution at ω ≈ 187 meV.This mode is an out-of-phase phonon mode and discussed in Ref. [16].N defect complex as a function of time for the same cavity parameters as considered for C B V N in the main manuscript.We clearly observe some structure in the population dynamics, but it is clear that the lack of a dominant phonon peak at small detuning from the ZPL means that the population dynamics of the C 2 C N defect are significantly less affected by the non-Markovian effects than was the case for C B V N .
Following the procedure in the manuscript, Figs. 3 (a-c) show the absorption spectra when an artificial scaling of the Huang Rhys factor (HRF), α HRF , is introduced to the phonon spectral density for C 2 C N .Again, we see that the impact of phonons for g = 50 meV is reduced in comparison to those found in the example of C B V N .This is emphasised in Fig. 3 (d), where we plot the renormalised Rabi spltting as a function of the HRF scaling parameter.Here we see that for g = 5 meV and g = 50 meV the scaling matches the trend predicted by the Frank Condon factor F , defined in the manuscript.We do, however, see clearly the emergence of the hybridised polaron-polariton state in Fig. 3 (c), which is reflected in the departure from Frank Condon physics for g = 100 meV and κ = 50 meV, where the polariton splitting approaches the high-Q mode of C 2 C N at ∼ 100 meV.

IV. COMPUTATIONAL DETAILS OF THE TEMPO CALCULATIONS
The two parameters that need to be converged in TEMPO calculations are the time step and the singular value decomposition cut-off used in the truncation of the bond dimension of the MPS.For all configurations considered in this work, we find that a SVD cut-off of λ c = 10 −6 is sufficient to ensure converged calculations.For each set of parameters, we separately converge the time step and we end up using a time step of δt = 5 eV −1 for the simulations with g = 5, 15 meV, and δt = 3 eV −1 and 2 eV −1 for g = 50 meV and 100 meV respectively.
For the spectrum calculations, it is important to ensure that the system is propagated long enough to capture the full dynamics, which we ensure in our calculations.However, because the frequency space resolution after the FFT is inversely proportional to the number of elements in the time array, it makes sense to pad with zeros to improve the resolution of the spectrum.We emphasize that since we ensure that all dynamics is properly accounted for by propagating the system to equilibrium, this padding has no effect on the physics of the system.We pad all time arrays such that the full array includes 2 15 time steps.Note that it is crucial that the total number of array elements is a power of 2 to ensure that the FFT needed to calculate the spectrum runs efficiently.FIG. 3. Effect of the electron-phonon coupling on the interaction with the cavity mode for the C2CN defect.(a-c) the linear absorption spectra when probing the cavity mode, for various values of the scaling parameter α.Here structure in the absorption spectra becomes increasingly important for large coupling, and broad cavity widths, but the effect is less pronounced than for CBVN.(d) shows the change in peak position as a function of α (points) for the cavity parameters studied in (a-c)

FIG. 1 .
FIG. 1. (a)A schematic figure of the CBVNdefect in hBN interacting with a single quantised field mode.This mode could be plasmonic or photonic in nature.(b) A plot of the ab initio (solid) and effective width (dashed) spectral densities.Using the atomistic simulation, we ascribe the dominant peak at ∼ 125 meV to a phonon breathing mode.(c) shows the bath correlation function for the ab intio (solid) and effective width (dashed) spectral densities.The sharp peaks of the ab initio spectral density leads to long lived oscillations in the correlation function, which are washed out when an effective width approximation is used.

3 FIG. 2 .
FIG. 2. Time dependent emitter population (top) and cavity occupation (bottom) for different coupling strengths and cavity widths.In contrast to the effective width spectral density (dashed), the ab initio spectral density (solid orange) predicts significant oscillations in the cavity occupation over all parameter regimes.Parameters used are T = 4 K, Γ = 4 meV, and Ωc = ωe − λ where λ = ∞ 0 dωω −2 J Ph (ω) is the reorganisation energy.The step size used to obtain convergence was dt = 3.2 fs, with SVD cut-off C = 10 −6 .

FIG. 3 .
FIG.3.(a-c) the linear absorption spectra when probing the cavity mode, for various values of the scaling parameter αHRF.Here structure in the absorption spectra becomes increasingly important for large coupling, and broad cavity widths.(d) shows the change in peak position as a function of αHRF (points) for the cavity parameters studied in (a-c).All other parameters are the same as Fig.2 where we have introduced the superoperators V δt [ρ] = exp(L 0 δt/2)[ρ], and W δt [ρ] = U δt ρU † δt with the unitary environment propagator defined as U δt = exp −i(H PH I + H PH B

3 FIG. 2 .
FIG. 2. Time dependent emitter population (top) and cavity occupation (bottom) for different coupling strengths and cavity widths for the C2CN defect.Parameters used are T = 4 K, Γ = 4 meV, and Ωc = ωe − λ where λ = ∞ 0 dωω −2 J Ph (ω) is the reorganisation energy.The step size used to obtain convergence was dt = 3.2 fs, with SVD cut-off C = 10 −6 , the same as for CBVN in the main manuscript. 2