{'fraction_non_alphanumeric': 0.03355177642980936, 'fraction_numerical': 0.008205155979202773, 'mean_word_length': 4.6242765763021625, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum optical master equation for solid-state quantum emitters

We provide an elementary description of the dynamics of defect centers in crystals in terms of a quantum optical master equation which includes spontaneous decay and a simplified vibronic interaction with lattice phonons. We present the general solution of the dynamical equation by means of the eigensystem of the Liouville operator and exemplify the usage of this damping basis to calculate the dynamics of the electronic and vibrational degrees of freedom and to provide an analysis of the spectra of scattered light. The dynamics and spectral features are discussed with respect to the applicability for color centers, especially for negatively charged nitrogen-vacancy centers in diamond.

We provide an elementary description of the dynamics of defect centers in crystals in terms of a quantum optical master equation which includes spontaneous decay and a simplified vibronic interaction with lattice phonons. We present the general solution of the dynamical equation by means of the eigensystem of the Liouville operator and exemplify the usage of this damping basis to calculate the dynamics of the electronic and vibrational degrees of freedom and to provide an analysis of the spectra of scattered light. The dynamics and spectral features are discussed with respect to the applicability for color centers, especially for negatively charged nitrogen-vacancy centers in diamond.

I. INTRODUCTION
Nowadays quantum optics avails itself not only of purely atomic sources of quantum light. The demand for easy-to-handle sources of non-classical light for modern experiments and quantum technological applications triggered a wide search for alternatives in several fields of physics. Among solid state systems, defect centers in diamond, primarily the negatively charged nitrogen-vacancy (NV − ) centers, have experienced an impressively successful development [1] during the last decade. These defect centers provide optically addressable electronic transitions between discrete quantum states energetically situated between the valence and conduction band. Antibunching of the emitted light as indicator for a single photon source was shown [2], Wheeler's delayed-choice experiment was realized [3] and coupling to optical microcavities was performed [4] to mention a few of the experimental breakthroughs. Moreover, applicability of the centers for quantum information technology using coupled electronic and nuclear spins has been shown [5][6][7], and the defect centers can serve as high precision magnetic sensors [8][9][10].
In contrast to atoms, solid state quantum emitter do not need to be laser cooled and trapped, as they are fixed in their solid state matrix. However, this embedding usually leads to strong coupling of the electronic degrees of freedom to lattice vibrations. While single atoms ideally only couple to the free radiation field, resulting in a damping due to spontaneous decay, their solid state counterparts experience among other effects, additional dephasing or non-radiative decay due to the coupling with crystal phonons. In the atomic case, the damping is typically described by a Markovian master equation. Solid state systems lack such a handy theoretical description based on a simplified physical picture. With the contribution we develop here, we hope to participate in filling this gap.
One of the first, most elementary and quantitative model for color centers was developed by Huang and Rhys [11] based on the Franck-Condon principle applied to color centers in crystals. It models the vibrational modes all at the same frequency in a continuum approximation, hence being essentially a single mode concept. Later developments generalized the approach to an arbitrary vibronic spectrum [12,13] and afterwards these ideas were applied to diamond [14], especially for the calculation of absorption and emission spectra. More recent models include pragmatic approaches [15], using somehow artificially but successfully a certain number of vibrational levels associated with the electronic ground state, with each transition to the excited state damped by an individual Markovian bath. A description based on a non-Markovian master equation is reported in [16] for general phononic spectral functions.
The numerous applications, especially of NV − centers, in quantum photonics call for a tractable theoretical description of the dipole-phonon interaction which renders the basics of the photon's spectral properties and is based on a clear physical picture. In this work we will present such a model by combining the Franck-Condon principle with a master equation formalism.
The basic idea is to describe the coupling to the vibrational modes by a single harmonic oscillator. Many defects in lattices are accompanied by the presence of strongly coupled localized modes, which have a exponentially decaying [13] amplitude away from the defect. These modes are represented by the single oscillator in our model. The delocalized phononic modes, in contrast, play the role of a temperature bath for the oscillator. We bring this idea into the form of a master equation and provide its full solution in terms of the damping basis, i.e. a spectral decomposition of the Liouville operator which generates the dynamics. In terms of this biorthogonal basis, we show how the absorption spectrum of a single defect center can be easily represented and analyzed. In such a way the basic features in the spectrum of light of solid state quantum emitters, in particular NV − centers, can be reproduced by a simple dynamical equation for the system's density operator. This paper has the following structure: In Sec. II we present the model and set up the master equation. FIG. 1. a) The Coulomb interaction of the ground state electronic orbital (shaded gray) of a defect center lets the surrounding atoms take on a certain lattice position. b) The excited electronic orbital pushes the adjacent crystal atoms into a displaced equilibrium configuration. c) Model of the Franck-Condon like lattice-defect center interaction in single mode and harmonic approximation including damping: In the excited state |e , the wavepacket of the single vibrational mode experiences a potential displaced by x0 compared to the potential in the ground state |g . The interaction causes a renormalized transition frequency between the states |e and |g , shifted by the relaxation frequency ωR. The single vibrational mode is coupled to a temperature bath modeling the phononic modes of the crystal and thermally relaxes with rate γ. The electronic two-level system decays with rate Γ due to spontaneous emission.
Sec. III is devoted to the damping basis as a solution of the dynamical equation. We derive the left-and right eigenelements of the Liouville operator together with its eigenvalues. In Sec. IV we use the damping basis to describe the time evolution and calculate the absorption spectrum of a single defect center. We discuss the features of the model in Sec. V, and put it into context with other descriptions. Finally, in Sec. VI we draw the conclusions. In the appendix we provide some details and techniques used for the calculations in the main text.

II. THE MODEL
In this section we develop a simplified theoretical description of a single color center coupled to a vibrational mode. This description employs the fact, that among the vibrational normal modes, often only a few localized modes strongly couple to the defect center, while delocalized phononic modes are considered here to play the role of a thermal bath causing the singled out localized mode to thermally relax.

A. Vibronic coupling
We consider two electronic states with transition frequency ω, a stable state |g and an energetically higher state |e . For an NV − center these states can be identified with the 3 A 2 and 3 E levels [17], respectively. The common understanding of the defect-vibration interaction assumes electronic-state dependent equilibrium positions of the atoms around the defect [13], see Fig. 1: Within the Franck-Condon assumption [18], atoms on the lattice sites around the defect experience a modified potential due to the reconfigured electronic orbitals, leading to an interaction term in the total Hamiltonian H when only a single vibrational mode with coordinate x is taken into account. Here the abbreviations σ g = |g g| and σ e = |e e| for the atomic projectors have been introduced for later convenience. For diamond, local vibration mode calculations [19] let it appear reasonable that only a few local modes strongly couple to the center, justifying the single mode assumption employed here. In the harmonic approximation and when expanding the potential around its new equilibrium position x 0 = F/M ν 2 up to first order in x, one finds We defined the Hamiltonian of a harmonic oscillator with frequency ν, momentum p and effective mass M associated with the local vibrational mode. The annihilation and creation operators b and b † are connected to the normal coordinate and momentum operators by with the ground state length scale ξ = /2M ν. In the second line of Eq. (3), we furthermore introduced the parts of the Hamiltonian representing the energy of the uncoupled two-level system and the interaction between the electronic and the vibrational degree of freedom with a coupling constant η = F ξ. Note that the coupling leads to an energy relaxation of ω R = F 2 /2M ν 2 , what can already be seen in the form (2) of the interaction. Up to this point, the model corresponds basically to the description of Huang and Rhys [1,11] of defect centers, describing the phononic degree of freedom by a single mode whose coordinate corresponds to the displacement of the surrounding lattice atoms from their equilibrium position. In this one-parametric model the dimensionless Huang-Rhys factor S = ω R /ν = (η/ν) 2 measures the vibronic interaction strength.

B. Spontaneous decay and phononic damping
To complete the description, we take electronic and vibrational relaxation processes into account using a master equation formalism. The time evolution of the system's density operator ̺, covering the electronic and vibrational degree of freedom, reads The first term on the right-hand side corresponds to the coherent evolution governed by the Hamiltonian H, Eq. (3). Spontaneous decay with the rate Γ is taken into account by the superoperator where we defined the lowering operator σ − = |g e| with the corresponding raising operator σ + = |e g| and used the abbreviation D[X]̺ = 2X̺X † − X † X̺ − ̺X † X to represent a Lindblad-form term. We further assume a finite lifetime of the vibrational excitation. Imperfections in the lattice structure or non-harmonic corrections may lead to a coupling between the vibrational modes and let the mode b interact with a bath of delocalized phononic modes with temperature T . In a Markovian approximation, such a coupling is described by and damps the localized mode irreversibly towards its stationary state with mean vibrational occupation numberm = [exp( ν/k B T )−1] −1 on a time scale given by the inverse decay rate 1/γ. Such a picture should be suitable for materials with few localized modes coupled strongly to the defect center compared to the delocalized modes forming the bath.

III. SOLUTION OF THE MASTER EQUATION: DAMPING BASE
An elegant form of providing a solution to the differential equation (9) is given in terms of the damping ba-sis [20], namely the spectral decomposition of the Liouville operator L. To this end, we seek for the left and right eigenelements of the Liouvillian, fulfilling The eigenelements are orthogonal by construction with respect to the scalar product Tr[̺ † λ̺ λ ′ ] = δ λ,λ ′ . If they form a complete set, the closure relation can be written as where the outer product is defined by its action (̺ λ ⊗ ̺ λ )X = Tr[̺ † λ X]̺ λ on an arbitrary operator X. The time evolution of an initial state ̺(t = 0) can then be written as where the coefficients are given by The remainder of this section is dedicated to the derivation of the explicit form of the damping basis of the problem at hand. The reader who is not interested in this technical part may skip the rest of this section without loss of coherence. The explicit formulas are summarized at the end of this section.
For later convenience we reformulate the master equation, Eq. (9), by introducing the Liouvillians describing a two-level system underlying spontaneous decay and a damped harmonic oscillator, respectively. Additionally we introduce and write the complete Liouvillian as The superoperators L tls and L osc exclusively act on the electronic or vibrational subspace, respectively. Later on we will hark back to the damping basis of L osc , which was derived in Ref. [20]. For the sake of completeness, we list the explicit expressions here: The normalized eigenelements read (where colons : ... : denote normal ordering, L l n (x) the associated Laguerre polynomials) with their corresponding eigenvalues

These eigenelements fulfill the completeness relation
In the following, we will express the eigenbasis of the coupled system in terms of these eigenelements of a damped quantum harmonic oscillator.

A. Right eigenelements
Inspired by the damping basis of a two-level system with spontaneous decay [21], the form of the Liouvillian, Eq. (9), suggests the ansatẑ for the right eigenelements. Here, µ are operators of the oscillator only and the eigenvalues λ depend on n, l and the superscript index of Eqs. (24)-(26).

Eigenelements connected to the atomic populations
We begin with Eq. (26), since it will also deliver the result for Eq. (24). Plugging the ansatz into Eq. (12) and using L in the form (19) yields where we defined the Liouvillian of a displaced harmonic oscillator with the displacement operator This quantity will appear throughout the rest of the treatment of this problem and is connected to a generalization of the Huang-Rhys parameter. One possibility to fulfill Eq. (27) is to chooseμ Then only a component in the σ g -subspace remains, obeying This obviously is the case from Eq. (24), and we consequently relabeled the operator in Eq. (30) correspondingly. The solution is given by Eq. (20) and the eigenvalues are the ones from Eq. (22). We now turn to the caseμ with the eigenvalues λ = λ nl − Γ.

Eigenelements connected to the atomic coherences
It remains to deal with the eigenelements associated with the electronic coherences, Eqs. (25). The action of L yields Here, we introduced To handle the eigenvalue equation, it is helpful to transform the expressions according tô and choose the displacement parameters α ± and β ± such that L ± µ → [b, µ], at the cost of an additional constant term. This turns out to be possible if one sets The condition on the curly bracket of Eq. (33) then transforms into the eigenvalue equation with the renormalized atomic transition frequency and linewidth are composed of the harmonic oscillator part and the renormalized electronic transition frequency and linewidth.

B. Left eigenelements
The derivation of the left eigenelements of L goes along the same lines as for their right counterparts. The action to the left [22] of L is represented by where the short notationĎ serves as an ansatz. Following an analogous treatment as in Sec. III A, we start with Eq. (45), as it will again deliver the solution for Eq. (47) too. One arrives ať to the eigenvalue λ = λ nl − Γ. Also the electronic coherence parts can be derived similarly as before, but this time using the transformatioň with α ± and β ± from Eqs. (35)-(38). One finally obtainš for the left eigenelements.

C. Summary
In the previous subsections we found that the eigenelements of the model Liouvillian L, Eq. (19), are assigned to the stable electronic subspace σ g , the electronic coherences σ ± and a subspace which can be associated with the population decay in the excited state at rate Γ, They fulfill the duality relation Tr[̺ both having the common eigenvalues for non-negative integers n and arbitrary integers l. The renormalized frequency ω and linewidth Γ of the two-level system are given by Eqs. (40) and (41), respectively. Moreover, α ± , β ± and ς are defined in Eqs. The dynamics of an arbitrary initial state in terms of the damping basis is given by Eq. (15). We focus on the case where the oscillator is initially in thermal equilibrium and the two-level system is prepared in an arbitrary state ρ 0 : ̺(t = 0) = ρ 0 µ th . For the time evolution of the reduced density operator ρ(t) = Tr osc [̺(t)] describing the dynamics of the atom, one finds The expansion coefficients c nl ̺(t = 0)] are determined by the initial state ̺(t = 0) and are calculated using the right eigenelements̺ where we defined the trace over the oscillator operators in̺ nl as The evaluation of these traces is conveniently carried out using phase space distributions, as it is sketched in App. B 1. The resulting functions depend only on β and m and are presented in formulas (B8)-(B9). We now turn to the expansion coefficients in Eq. (62). Because of Eqs. (63)-(64), for the two coefficients in the first line we only need the n = l = 0 case, giving unity and the initial excited state population ρ ee = e|ρ 0 |e . In the second line of Eq. (62), the coefficients for the coherences are required, where we have defined the initial values ρ ± = Tr[σ ∓ ρ 0 ] and Again we use phase space representations of the required eigenelements to evaluate the trace. The details are re-ported in App. B 2. It is clear from the form of the solution (62) that and A + nl B + nl = [A − n,−l B − n,−l ] * only appear as products. As a matter of course, for t = 0 the initial state must be reproduced by the solution (62) of the master equation. Hence, if the eigensystem is complete, the remaining sum over n and l has to fulfill this requirement. For finite times, the exponentials containing the eigenvalues generate the dynamics. The form of the eigenvalues is linear in n and l, such that the time dependency can be easily incorporated into Eq. (71). The sums can then be performed with the help of the power series representation and the generating function of the modified Bessel function I l (x) of order l [23]. After performing these steps one arrives at the solution with the time-dependent quantities The expression (73) indeed reproduces the initial electronic state for t = 0, as expected. For finite times, the population in the excited state decays with the natural rate Γ of the two-level system. It is the coherence between the ground and excited state which is modified by the coupled phonon mode. This influence on the temporal behavior is twofold: First, the operators σ ± , Eq. (75), rotate and decay with the renormalized frequencies from Eq. (60). Second, the evolution is modulated by the timedependency ofβ 2 (t), Eq. (74). The oscillator hence generates an effective dephasing of the two-level system. Although the dynamics of the oscillator appears to be less important for practical purposes, we nevertheless briefly delve into it in order to provide additional insight into the overall temporal behavior. We focus on the reduced density operator of the oscillator, µ(t) = Tr tls [̺(t)] with the same initial state ̺(t = 0) as above. The dynamics is described by Eq. (62) with the trace this time going over the electronic degrees of freedom, leaving only the terms in the first line. With the help of the completeness relation (23) and the overlaps of displaced harmonic oscillator eigenelements calculated in App. B 3, the expression for µ(t) can be recast into a form, where the time-dependent exponential factors can be incorporated into the displacement parameter Then, the time evolution of the oscillator's state can be expressed as The result consists of three contributions: (i) The thermal state corresponding to the |g -part of ̺(t = 0) showing no dynamics. (ii) A dynamically displaced thermal state component moving in the |e potential surface, whose population dies off with the electronic decay rate and (iii) the flow of the decayed population towards the undisplaced thermal state. In the asymptotic limit t → ∞ the oscillators state µ(t) ends in the thermal state µ th , as expected.
In Fig. 2 we illustrate the dynamics of both the electronic and the vibrational degree of freedom. The initial state of the two-level system is ρ 0 = |ψ 0 ψ 0 | with |ψ 0 = (|g + |e )/ √ 2 and a thermal state for the oscillator. Subplot a) shows the time evolution of the population in the excited state (dashed curve) which decays exponentially with rate Γ, as in the unperturbed twolevel system. The time evolution of the coherences (black solid curve) is modified by the oscillator: We show the absolute value of e|ρ(t)|g which exhibits a faster decay compared to the pure two-level case, superimposed by an oscillating modulation with frequency ν as described by Eq. (73). For larger temperature (gray solid curve) the coherences are damped away even faster. Fig. 2b) visualizes the temporal behavior of the oscillator's quantum state, Eq. (77), in phase space with the help of the Wigner function [24]. We show four snapshots of the Wigner function within one oscillation period τ = 2π/ν at the times marked in subfigure a) by empty circles. Initially, at t = 0, the phase space distribution corresponds to a thermal state, i.e. a Gaussian with variance proportional tom + 1/2. The portion of the distribution corresponding to the excited state starts to spirally circulate around the displaced equilibrium position at x/ √ 2ξ = Re β and p √ 2ξ/ = Im β and its population decays simultaneously at rate Γ, as described in the second term of Eq. (77). The black spiral visualizes the trajectory of the center of probability of this part while the blue (dashed) and red (solid) rings mark the variance of the Gaussian components moving in the harmonic potential of the ground and excited state, respectively. The contribution of the last term of Eq. (77) is hardly visible, but results, for example at t = 2τ /3, in a small asymmetry with respect to the axis connecting the centers of the two Gaussians. At t = τ , the components moving in the ground and excited state potentials strongly overlap again, leading to the partial revival of the electronic coherence in Fig. 2a).

B. Absorption and emission spectrum
The advantage of the eigenvalue decomposition of the Liouvillian becomes apparent most clearly when calcu-lating the spectral properties of the system. The method we use here was already applied for a systematic analysis of the fluorescence light of laser cooled atoms [25] and is now carried over to the case of solid state quantum emitters. For a weak probe laser of frequency ω L illuminating the defect center, the spectrum of absorption under stationary conditions is determined by [26,27] In this form the correlation function is evaluated without the probe. The time evolution of the operators is determined by the Liouvillian L, Eq. (19), and can be calculated using the quantum regression theorem [28] to . The stationary state of the undriven system is ̺ st = |g g|µ th , and hence, the second term of the commutator vanishes such that the spectrum takes on the form In the last step we expanded the spectrum in terms of the eigensystem of L, Eqs. (53)-(61), and performed the integration. The traces in Eq. (80) are readily calculated, we encountered them already in the previous subsection. They give with the explicit form of A − nl and B − nl specified in App. B and the product W nl = A − nl B − nl shown in Eq. (71). The decomposed spectrum is a superposition of curves associated with each eigenvalue λ (−) nl , weighted by the function W nl which also determines the shape: For real W nl the spectral contribution is a Lorentzian, for purely imaginary weight it is a Fano-like profile, and in the general case it is a mixture of both, positioned at ω L = Im λ (−) nl with a width given by Re λ (−) nl . From the l dependency it follows that the spectrum is built up by a group of peaks separated by the mechanical frequency ν and situated around the renormalized transition frequency ω. The peaks are distinguishable for Γ, γ ≪ ν and otherwise overlap. Each peak is a sum of several contributions of increasing width for rising n, whereby the n > 0 contributions vanish when m → 0. In this limit the weight factors take on the form for l ≤ 0 and otherwise zero. It goes over into a Poissonian distribution when β becomes real. For vanishing thermal occupation of the phonon mode, only phonons can be created during the absorption process and consequently, the spectral components for l > 0 disappear. Note that the absorption and excitation spectrum of the two-level system coincide in the situation discussed here.
We also remark that we had arrived to the same result by calculating the action of the propagator e Lt in Eq. (79) to the right using the solution (73). The integral representation can be advantageous in some cases, for example in the high temperature limit, when the expansion, Eq. (80), converges slowly. In Fig. 3 we show the absorption spectrum, Eq. (82), as a function of the probe frequency ω L for a temperature close to zero. The spectrum is composed by several resolved peaks. The zero-phonon line corresponding to a pure electronic transition without creation or annihilation of phononic excitations is shifted by the relaxation frequency ω R = ω − ω from the frequency of the bare transition |e ↔ |g . More peaks are visible at multiple distances of the oscillator frequency and extend towards higher frequencies, forming the vibronic sidebands corresponding to absorption processes accompanied by the creation of phonons. Phonons to be annihilated are only rarely present at low temperatures entailing strongly suppressed signals on the lower frequency side of the zerophonon line. The dashed curves show the most relevant terms of the sum in Eq. (82), that is the individual spectral components associated with a certain eigenvalue which compose together the total spectrum. Note that the single contributions can become negative, only the superposition generally yields a valid spectrum.
We conclude this section by discussing the spectral properties of a photon emitted from the defect center.
The emission spectrum reads [29] We assume here that the two-level system is initially excited and the oscillator is in the corresponding equilibrium µ 0 = D(β)µ th D † (β). This is a reasonable assumption if the lifetime of the excited state Γ −1 is much larger than the thermalization time scale γ −1 of the phonon mode. Again the quantum regression theorem and the completeness relation of the eigenelements are applied to bring the spectrum into the form The further evaluation uses and thus λ ′ = −Γ. The remaining trace in Eq. (85) can be calculated using the techniques of App. B and yields (−1) l B + nl e −iImβ 2 , viz. Eq. (B17). The spectrum contains the weight factors W ′ nl = W * n,−l , signifying the fact that the emission spectrum is a mirrored version of the absorption spectrum [30] with respect to the modified transition frequency ω.

V. DISCUSSION
We modeled the vibronic interaction of defect centers in crystals using a single-mode description including the radiative damping of the electronic degree of freedom as well as the damping of the vibrational mode by a heat bath with a temperature corresponding to a mean vibrational quantum numberm. The model can therefore be considered as an extension of the artificial but instructive single-mode Hamiltonian description, sometimes called the Huang-Rhys model [1] or Einstein model [30], which is the limiting case of our model when both linewidths, Γ and γ, tend to zero.
From the temporal behavior treated in Sec. IV A follows that the vibronic coupling results in an effective dephasing mechanism for the electronic subsystem: It leaves the dynamics of the population unaffected but expedites the decay of the electronic coherences. This decay of the coherences is superimposed by an oscillatory behavior whose origins can be traced back to the motion of the vibrational degree of freedom in Wigner phase space. In the spirit of master equations, the dynamics of the two-level system can be considered as being dephased by a structured environment with harmonic time evolution, when tracing out the vibrational coordinates. The different types of dynamics concerning the populations and coherences is formally resembled in the classes of eigenvalues, Eqs (59)-(61): The excited state simply decays with the natural rate Γ. In contrast, the electronic coherences oscillate with a renormalized frequency ω = ω − ω R and decay with the increased, temperaturedependent rate Γ = Γ + Γ R . The relaxation frequency ω R = |β| 2 ν goes over into ω R = Sν, proportional to the Huang-Rhys factor S, when the linewidth γ tends to zero. For finite linewidth, the damping rate γ additionally modifies the transition frequency of the two-level system, reflecting the fact that the electronic subsystem indirectly couples to the Markovian environment which damps the oscillator. Compared to the uncoupled case, the decay rate is increased by Γ R = |β| 2 γ(2m + 1).
Concerning the spectra, again a comparison with the Huang-Rhys model suggests itself. According to Ref. [31], this single-mode Hamiltonian model in the har-monic and Franck-Condon approximation, results in an absorption spectrum consisting of Dirac-δ-shaped spikes displaced from the renormalized frequency ω − ω R by multiples of the vibrational frequency ν, whose envelope (related to the intensities) is described by a Poissonian distribution I ℓ = e −S S ℓ /ℓ!, when ℓ numbers the vibronic sidebands. The generalized formulation applied here reflects these characteristics formally in Eq. (83), with the Huang-Rhys parameter replaced by the complex parameter β * 2 quantifying the vibronic interaction strength. For finite temperature, the Poissonian-like distribution breaks down and the intensity of the l-th vibronic sideband relative to the total spectral intensity is described by with θ l = Im β 2 + arg β 2l and x = ν/2k B T . This is a meaningful expression when the vibronic sidebands are resolved, that is for γ, Γ ≪ ν. The dependency of the absorption spectrum on the temperature is depicted in Fig. 4a). For low temperaturem ≈ 0 (solid line), the spectrum consists of the zero phonon line at ω L = ω − ω R with ω R ≈ 2ν for the parameters considered here. towards higher frequencies the vibronic sidebands are stretched out, starting to mutually overlap due to the swelling width of increasing orders. For higher temperature,m = 0.25, the individual lines are only hardly distinguishable and extend also towards lower frequencies on the red side of the zero-phonon line. The position of the zero phonon line is not temperature dependent, in contrast to its width given by Re λ (−) n0 , Eq. (60). Finally form = 1 only a broad spectral shape is present, centered approximately at the bare transition frequency ω of the two-level system. For ν ≫ γ, Γ this single peak approaches a Gaussian form around ω L = ω whose width is proportional to √m . Fig. 4b) shows the absorption spectrum at low temperature for three different vibronic coupling strengths, η = 0.9ν, 1.2ν and 1.5ν, corresponding to the dotted, dashed and solid curves, respectively. The peak strength follows the quasi-Poissonian distribution (83), with more peaks appearing for larger coupling: then higher order phonon processes become more likely. The position of the zero-phonon line moves towards lower frequencies for larger couplings in accordance with the renormalized frequency ω defined in Eq. (60) and the simplified physical picture of Fig. 1. Moreover, the curves in Fig. 4b) demonstrate that an increased coupling goes along with a broadening of the spectral components, as described by Γ, Eq. (60).
The master equation description of this work can also be related to the independent boson model [30], in which a linear coupling to essentially a continuum of phononic modes is taken into account. The phase diffusion in the temporal evolution of the phononic modes is at the heart of the dephasing mechanism there, being replaced in the present description by a Markovian bath. Consequently, the spectral function describing the frequency dependent coupling strength of the phononic continuum in that model is replaced here by a single parameter, namely the damping rate γ of the singled out phononic mode. Depending on the intended application, the more detailed independent boson model or the one-parametric simplified description presented here can be advantageous.
Finally we comment on the applicability of the model for NV − centers in diamond. As already mentioned, the presence of a few localized modes of similar frequencies suggests to highlight these modes -or ultimately a single mode -as it was done here. A temperature bath damping the mode has clear physical origins in the continuum of the delocalized phononic modes. The Markovian description, however, is a simplification which works reasonably well for rendering the basic properties of the NV − dynamics, but does not include the details of the temperature dependency of the zero-phonon line, for example. This would at least require to replace the flat spectral function implicitly applied here, if not more involved processes [32] have to be included, depending on the desired accuracy of the description. Of course, effects stemming from the dynamical Jahn-Teller effect are absent. Nevertheless, the model renders the basic properties of the spectra of NV − centers, especially in the room temperature regime and goes beyond simple pure dephasing models or artificial fitting models. Note that the room temperature spectra still correspond tom ≪ 1 due to the high Debye-temperature of diamond.

VI. CONCLUSIONS
We formulated a compact theoretical description for the dynamics of defect centers in crystals. This master equation description is based on a linear vibronic coupling in the Franck-Condon, single phonon mode and two-level approximation and includes radiative and phononic Markovian damping. The overall solution of the master equation is given in terms of the eigenelements and eigenvalues of the Liouville operator. This solution covers both, the electronic and phononic degree of freedom and can be further reduced to a master equation of the electronic two-level system only by tracing out the undesired degree of freedom. We demonstrated the calculation of the absorption and emission spectrum at different temperatures, which take on a insightful form and reproduce the basic features of real spectra, for example for NV − centers in diamond. The physical picture at the heart of the presented model was reviewed in detail by means of the dynamics of both involved degrees of freedom.
In the form presented here, the model can be easily extended by adding more vibrational modes in order to increase the accuracy of the theory, if desired. Moreover, such an extension would also be interesting from the point of view of a systematic approximation of open quantum dynamics [33]. It is also entirely conceivable to extend the treatment to include quadratic vibronic coupling or more electronic levels.
We believe that the model, based on an easy but clear physical picture, will be helpful for investigations of defect centers in quantum optical applications, for example as single photon sources, when the influence of the vibronic sidebands becomes important. With the analytic solution provided here it can deliver a basis for the description of coupled cavity-defect-center systems or similar hybrid setups.

VII. ACKNOWLEDGMENTS
We thank Giovanna Morigi, Christoph Becher and Ignacio Wilson-Rae for helpful comments and discussions. Furthermore, MB and RB acknowledge support from the German Research Foundation (DFG) within the project BI1694/1-1 and from the GradUS program of the Saarland University. JMT acknowledges support by CONA-CyT Grant No. 162781, by the BMBF project Q.com, and by CASED III.
The solution leads -after the substitution β → b, β * → b † in normal order -to the eigenelementsμ nl , Eq. (20) [20]. In our case, first the coordinate transformation is reversed, and eventuallŷ results as the right eigenelement of Eq. (A1), with eigenvalues given by λ nl , Eq. (22). Similarly the corresponding left eigenelements can be derived using the action to the left of L shift .

Appendix B: Overlaps between eigenelements
In Sec. IV we need explicit expressions for the formulas following from overlaps of eigenelements of the Liouvillian L, Eq. (19). Moreover, is helpful for the application of the damping basis. (B9) In this notation, the upper and lower line in the curly brackets correspond to l < 0 and l ≥ 0, respectively.