DERIVATION OF THE STATIONARY FLUORESCENCE SPECTRUM FORMULA FOR MOLECULAR SYSTEMS FROM THE PERSPECTIVE OF QUANTUM ELECTRODYNAMICS

Lithuanian Journal of Physics, Vol. 61, No. 2, pp. 110–123 (2021)
© Lietuvos mokslų akademija, 2021

Numerical simulations of stationary fluorescence spectra of molecular systems usually rely on the relation between the photon emission rate and the system’s dipole–dipole correlation function. However, research papers usually take this relation for granted, and standard textbook expositions of the theory of fluorescence spectra also tend to leave out this important relation. In order to help researchers with less theoretical training gain a deeper understanding of the emission process, we perform a step-by-step derivation of the expression for the fluorescence spectrum, focusing on rigorous mathematical treatment and the underlying physical content. Right from the start, we employ quantum description of the electromagnetic field, which provides a clear picture of emission that goes beyond the phenomenological treatment in terms of the Einstein A coefficient. Having obtained the final expression, we discuss the relation of the latter to the present level of theory by studying a simple two-level system. From the technical perspective, the present work also aims at familiarizing the reader with the density matrix formalism and with the application of the double-sided Feynman diagrams.

Keywords: quantum molecular electrodynamics, dipole–dipole correlation function, Einstein coefficients, Feynman diagrams

1. Introduction

Fluorescence (or emission) spectroscopy is the key tool for obtaining information about the bright excited states in molecular systems [1]. The concepts behind the experiment and its implementation are not difficult, and the simplest theoretical description follows from the usual textbooks on quantum mechanics. On the other hand, researchers are often confronted by the need to understand fluorescence spectra of condensed matter systems. These systems are complex entities, usually described within the theoretical framework of open quantum systems. It provides a way of calculating the evolution of the transition dipole moment operator and, consequently, its autocorrelation function. This function is usually called the dipole–dipole correlation function, and it fully describes the stationary photon emission process. Thus, in order to understand, and, more importantly, to simulate the fluorescence spectra of molecular systems, one requires a deep understanding of the relation between the fluorescence spectra and the dipole–dipole correlation function. Unfortunately, this is rarely detailed in textbooks or monographs.

Indeed, while the relation of the absorption spectra of condensed matter systems to the corresponding dipole–dipole correlation function is quite easily understood and is often presented in literature [2–4], that is not the case for fluorescence spectroscopy. The key difficulty lies in the fact that classical description of the field is sufficient when analysing absorption, but quantum description must be used for fluorescence. Another issue is that the system under consideration has to be described by the density operator, rather than the wave function. That is why most quantum mechanics textbooks describe fluorescence only in terms of the Einstein coefficients [5–7]. The same is true for more specialized textbooks pertaining to spectroscopy [8, 9], while books specifically on fluorescence are more focused on concepts and experimental issues, neglecting the theoretical description [1, 10]. The seminal textbook by Mukamel does contain the description of the fluorescence in terms of a quantum field, but it focuses on the time-resolved fluorescence, complicating the description [11]. Research papers, on the other hand, usually either just state the end result without much justification [12–16] or simply cite Mukamel’s book [17, 18]. Therefore, researchers new to the field face significant difficulties when they try to understand the theoretical concepts and reasoning behind the results stated in literature.

In this paper, we present a detailed derivation of the relation between the stationary (or relaxed) fluorescence spectrum of a condensed matter system and the relevant dipole–dipole correlation function. This is done by employing the quantum description of the electromagnetic field and the density operator formalism. Even though the derivation is somewhat lengthy, the mathematical procedures are not difficult, as long as they are done carefully. Moreover, it is instructive to demonstrate how one may deal with rather bulky expressions, which may sometimes be done partly by following the physical reasoning. Having derived the formula for the fluorescence spectrum, we apply it to a simple two-level system and show its relation to the Einstein A coefficient. We believe that our work will be beneficial for researchers seeking to bridge the gap between textbook knowledge and research papers.

2. Preliminaries

In its most general form, the total Hamiltonian of a condensed matter system interacting with the electromagnetic field may be written in the following form,

{\hat{H}}_{T} = {\hat{H}}_{s} + {\hat{H}}_{SF} + {\hat{H}}_{F}, (1)

where the three terms describe, correspondingly, the system, its interaction with the electromagnetic field, and the electromagnetic field itself. We will discuss each term in more detail below.

The molecular system under consideration described by the Hamiltonian Ĥ_Smay either be isolated or it can be ‘open’ – that is, interacting with its environment, such as a solvent. In the latter case, the Ĥ_S operator is usually partitioned even further, separating the Hamiltonians of the system (usually the electronic degrees of freedom of the relevant molecules) and its environment (usually the vibrational degrees of freedom of both the relevant system and its environment), together with an additional Hamiltonian describing the system–environment interaction. The exact form of Ĥ_S will not be relevant for the present derivation as the results obtained below will be rather general, but the important point is that the final expression is valid for open systems.

The system–field interaction Hamiltonian is taken to be of the form

{\hat{H}}_{S F} = - \sum_{α = 1}^{N_{m o l}} {\hat{μ}}_{α} \cdot \hat{E} (r_{α}), (2)

where $\hat{μ}$ _α is the transition dipole moment of the αth molecule, N_mol is the number of molecules in the system, and r_α is the position of the αth molecule. The operator Ê(r_α) is the electric field strength operator that describes the field confined within the volume V. It is given by [11, 19]

\hat{E} (r_{α}) = \sum_{k, s} (i N_{k} p_{k s} {\hat{a}}_{_{k s}} e^{i k \cdot r_{α}} - i N_{k} p_{k s}^{*} {\hat{a}}_{k s}^{†} e^{i k \cdot r_{α}}), (3)

where summation over all possible wave vectors k and polarizations s is performed. Operators â_k_s and ${\hat{a}}_{k s}^{†}$ are bosonic annihilation and creation operators, and $N_{k} = \sqrt{ℏ c | k | / 2 ε_{0} V}$ is the normalization factor with V being the quantization volume. The method of complex amplitudes is being used here, so the polarization vector p_k_s is a complex quantity (similarly to the case when, for example, elliptically polarized waves are studied [19]). Assuming that V corresponds to a cube of side 𝓁, components of k are restricted to k_j = 2πn_j/𝓁, with j = x, y, z denoting the Cartesian components of the vector, and n_j being integers. Since p_k_s ⊥ k, summation over s means summation over two polarizations. It should be noted that the operator $\hat{μ}$ _α acts only on the system states, while the operator Ê(r_α) acts only on the field states. Therefore, the tensor product of these operators is to be understood in Eq. (2), i.e.

\begin{array}{l} {\hat{μ}}_{α} \cdot \hat{E} (r_{α}) \equiv {\hat{μ}}_{α}^{x} \otimes {\hat{E}}^{x} (r_{α}) \\ + {\hat{μ}}_{α}^{y} \otimes {\hat{E}}^{y} (r_{α}) + {\hat{μ}}_{α}^{z} \otimes {\hat{E}}^{z} (r_{α}), (4) \end{array}

where the symbol “⊗” explicitly denotes the tensor product. The Hamiltonian (2) corresponds to the dipole approximation [11, 19], whereby it is assumed that each molecule may be adequately described as a single electric dipole, neglecting higher order effects related to the magnetic dipole moment and electric quadrupole moment. Therefore, the electric field is assumed not to change appreciably over the dimensions of the molecules. This approximation is valid as long as the wavelength of the radiation is much greater than the size of the molecule. For example, the dimensions of the porphyrin ring of a chlorophyll molecule are ≲ 1 nm, clearly much less than the wavelength of visible light. Finally, the field Hamiltonian is [11, 19]

{\hat{H}}_{F} = \sum_{k, s} ℏ ω_{k} ({\hat{a}}_{k s}^{†} {\hat{a}}_{_{k s}} + \frac{1}{2}), (5)

where the mode frequency ω_k = c|k|. Note that by employing the Hamiltonians defined by Eqs. (2) and (5) we are using the second quantization notation, with the field being quantized in addition to the molecular system. As a consequence, the electric field strength (3) is described by an operator rather than a function, and it has no time dependence since we are working in the Schrödinger picture at this point. The (harmonic) time dependence will emerge naturally once we switch to the interaction picture and transform the creation and annihilation operators accordingly (see Eqs. (45, 46)).

3. Derivation

We seek to obtain the formula for the stationary (or relaxed) fluorescence spectrum of molecular systems, and we begin by discussing the assumptions involved. The relevant molecular system may be of arbitrary size, ranging from single isolated molecules to large molecular aggregates in a solvent or protein environment. The system is assumed to be excited by an infinitely short laser pulse, and the subsequent ‘equilibration’ – relaxation to the lowest electronic levels of the excited electronic state manifold – is assumed to be instantaneous as well. Therefore, we will not take into account any photons emitted during the equilibration period, but rather restrict the analysis to the photons emitted thereafter. Since radiative lifetime of the molecular excited states is often on the order of nanoseconds and the equilibration time of the excited states is a few to tens of picoseconds [3], our assumption of instantaneous excitation and following equilibration corresponds to this separation of timescales. Finally, we ignore the multiphoton relaxation and any non-radiative relaxation channels, such as thermal decay or internal conversion [3].

The fluorescence spectrum f(ω, t) of a molecular system is effectively a distribution function of photon frequencies which shows the number of photons per unit frequency range (around ω) that the system emits per unit time at time t after excitation. We aim at describing the relaxed fluorescence spectrum, therefore we will be interested in the steady-state distribution f(ω). As will be shown below, f(ω) is easily related to the single-mode distribution f_k_s(t) that describes emission of photons with frequencies corresponding to a specific wave vector, |k| = ω_k/c, and a specific polarization s. Thus now we will seek an expression for f_k_s(t).

From a purely quantum perspective, the rate of emission of a photon with a wave vector k and polarization s may be calculated as the rate of change of the occupation number corresponding to the (k, s) mode of the electromagnetic wave [11, 19]. The occupation number operator is given by

\hat{𝒩} = {\hat{a}}_{k s}^{†} {\hat{a}}_{k s} . (6)

We must now find the operator f̂_k_s that corresponds to the emission rate observable f_k_s. By identifying the photon emission rate with the emission spectrum [11], we adopt one of the several possible quantum definitions of the spectrum [20]. Such an intuitive definition turns out to be adequate for the present analysis, as the spectra calculated using the final formula agree with the measured spectra reasonably well [15, 16].

It follows from our definition of the emission rate that

f_{k s} (t) = \frac{d}{d t} 𝒩_{k s} (t), (7)

where 𝒩_k_s is the mode occupation number observable which operator ${\hat{𝒩}}_{k s}$ corresponds to. Meanwhile, it is known from standard courses in quantum mechanics [5, 21, 22] that if an observable A has a corresponding (time-independent) operator Â, then an observable B whose quantum average equals to the average of $\frac{d}{d t} A (t)$ has a corresponding operator

\hat{B} = \frac{i}{ℏ} [\hat{H}, \hat{A}], (8)

where Ĥ is the Hamiltonian of this exemplary system. Therefore, it follows that

{\hat{f}}_{k s} = \frac{i}{ℏ} [{\hat{H}}_{T}, {\hat{𝒩}}_{k s}] . (9)

We recall that Ĥ_T consists of three terms [see Eq. (1)], but the occupation number operator acts only in the Hilbert space of the field states, and it therefore commutes with Ĥ_S. Moreover, looking at Eqs. (5) and (6), we notice that ${\hat{𝒩}}_{k s}$ commutes with Ĥ_F as well. Therefore, [Ĥ_T, ${\hat{𝒩}}_{k s}$ ] = [Ĥ_SF, ${\hat{𝒩}}_{k s}$ ]. Further, the fact that the creation and annihilation operators corresponding to different modes mutually commute allows us to leave only a single (k, s) mode when substituting Eq. (3) into Eq. (9). Thus,

{\hat{f}}_{k s} = \frac{i}{ℏ} [{\hat{H}}_{k s}, {\hat{N}}_{k s}], (10)

where we define a single-mode system–field interaction Hamiltonian,

{\hat{H}}_{k s} = - \sum_{α = 1}^{N_{mol}} {\hat{μ}}_{α} \cdot (i N_{k} p_{k s} {\hat{a}}_{k s} e^{ik \cdot r_{α}} - i N_{k} p_{k s}^{*} {\hat{a}}_{k s}^{†} e^{-ik \cdot r_{α}}) . (11)

The aim is now to find the quantum average of the emission rate operator. This may be conveniently approached with the help of the density matrix formalism, so we start by introducing the total density operator for $\hat{ρ}$ _T(t) that describes both the system and the field. The key advantage of working with density operators rather than wave functions is that the final result will be suitable for open quantum systems as well. Let t₀ denote the excitation time; then, as the system state for $\hat{ρ}$ _S(t₀) has just been prepared, we may assume that initially there are no correlations between the system and the field. Hence, the total density operator is initially a tensor product of the system and the field parts [23],

{\hat{ρ}}_{T} (t_{0}) = {\hat{ρ}}_{S} (t_{0}) \otimes {\hat{ρ}}_{F} (t_{0}) . (12)

Once the system has been prepared, it starts interacting with the electromagnetic field, and the relation no longer holds at later times, with the evolution of the total density matrix being controlled by the total Hamiltonian (1). The required quantum average is then given by

\begin{array}{l} f_{k s} (t) = {Tr}_{T} ({\hat{f}}_{k s} {\hat{ρ}}_{T} (t)) \\ = {Tr}_{T} (\frac{i}{ℏ} [{\hat{H}}_{k s}, {\hat{𝒩}}_{k s}] {\hat{ρ}}_{T} (t)), (13) \end{array}

where the symbol Tr stands for the trace, and the index T (total) emphasizes that the trace is taken with respect to both the system and the field states. For the interpretation of subsequent equations, it is beneficial to use the trace identity (A.4) to recast (13) into

f_{k s} (t) = {Tr}_{T} (- \frac{i}{ℏ} {\hat{𝒩}}_{k s} [{\hat{H}}_{k s}, {\hat{ρ}}_{T} (t)]) . (14)

We now have to find the expression for the density operator. The density operator satisfies the Liouville–von Neumann equation [23],

\frac{d}{d t} {\hat{ρ}}_{T} (t) = - \frac{i}{ℏ} [{\hat{H}}_{o}, {\hat{ρ}}_{T} (t)] - \frac{i}{ℏ} [{\hat{H}}_{SF}, {\hat{ρ}}_{T} (t)], (15)

where we have partitioned the total Hamiltonian as

{\hat{H}}_{T} = {\hat{H}}_{0} + {\hat{H}}_{SF}, {\hat{H}}_{o} = {\hat{H}}_{S} + {\hat{H}}_{F} . (16)

This partitioning will allow us to treat the system–field interaction as a small perturbation and develop a perturbation expansion with respect to the electric field strength [11],

{\hat{ρ}}_{T} (t) = {\hat{ρ}}_{T}^{(0)} (t) + {\hat{ρ}}_{T}^{(1)} (t) + .... (17)

The initial condition is chosen such that at t = t₀ (when the fluorescence starts) the system has reached the excited state equilibrium. Further, we assume that the resulting density operator ${\hat{ρ}}_{S}^{eq} \equiv {\hat{ρ}}_{S} (t_{0})$ is known, while the field is in the vacuum state $| 0_{F} ⟩ ⟨ 0_{F} |$ ,

{\hat{ρ}}_{T}^{e q} \equiv {\hat{ρ}}_{T} (t_{0}) = {\hat{ρ}}_{S}^{e q} \otimes | 0_{F} ⟩ ⟨ 0_{F} | . (18)

Fluorescence thus corresponds to the spontaneous light emission process. Formally, the vacuum state is a tensor product of the zero-occupation states of all the modes. Specifically, let |0_k′s′⟩ denote the zero-occupation state of the (k′, s′) mode; then,

| 0_{F} ⟩ = \prod_{k^{'}, s^{'}} | 0_{k^{'}, s^{'}} ⟩, (19)

where the product symbol is to be understood as a tensor product. Hence, â_k_s|0_F⟩ = 0 for all k and s.

The zeroth order term in (17) corresponds to the unperturbed system with Ĥ_SF = 0, hence it corresponds to the equilibrium state,

{\hat{ρ}}_{T}^{(0)} (t) = {\hat{ρ}}_{T}^{(e q)}, t > t_{0} . (20)

For the first order term, Eq. (15) takes the form

\frac{d}{d t} {\hat{ρ}}_{T}^{(1)} (t) = - \frac{i}{ℏ} [{\hat{H}}_{o}, {\hat{ρ}}_{T}^{(1)} (t)] - \frac{i}{ℏ} [{\hat{H}}_{SF}, {\hat{ρ}}_{T}^{(0)} (t)], (21)

As we show in Appendix B, the formal solution is

{\hat{ρ}}_{T}^{(1)} (t) = - \frac{i}{ℏ} \int_{0}^{t - t_{0}} d t_{1} e^{- (i / ℏ) {\hat{H}}_{0} t_{1}} [{\hat{H}}_{SF}, {\hat{ρ}}_{T}^{eq}] e^{(i / ℏ) {\hat{H}}_{0} t_{1}} . (22)

Substituting ${\hat{ρ}}_{T}^{(0)} (t) + {\hat{ρ}}_{T}^{(1)} (t)$ for $\hat{ρ}$ _T(t) in Eq. (14) and using Eqs. (20) and (22), we obtain

\begin{array}{l} f_{k s} (t) = {Tr}_{T} ({\hat{f}}_{k s} {\hat{ρ}}_{T}^{eq}) + {(- \frac{i}{ℏ})}^{2} \\ \times \int_{0}^{t - t_{0}} d t_{1} {Tr}_{T} ({\hat{𝒩}}_{ks} [{\hat{H}}_{k s}, e^{- (i / ℏ) {\hat{H}}_{0} t_{1}} [{\hat{H}}_{SF}, {\hat{ρ}}_{T}^{eq}] e^{(i / ℏ) {\hat{H}}_{0} t_{1}}]) . (23) \end{array}

The first term on the right-hand side of Eq. (23) represents the time-independent zeroth order contribution to the fluorescence. However, the interaction with the field is totally neglected at the zeroth order, therefore, the said term is equal to zero (as can be verified directly). Since we are interested in the steady state emission rate, we will now set t₀ → –∞ (the system started to fluoresce in the infinite past) and use the cyclic invariance of the trace to rearrange the exponential factors. Equation (23) becomes

f_{k s} (t) = {(- \frac{i}{ℏ})}^{2} \times \int_{0}^{\infty} d t_{1} {Tr}_{T} ({\hat{N}}_{ks} [{\hat{H}}_{k s}^{I}, (t_{1}), [{\hat{H}}_{SF}, {\hat{ρ}}_{T}^{eq}]]) . (24)

where we have defined the single-mode system– field interaction Hamiltonian in the interaction picture as

{\hat{H}}_{k s}^{I} (t_{1}) = e^{(i/ ℏ) {\hat{H}}_{0} t_{1}} {\hat{H}}_{k s} e^{- (i/ ℏ) {\hat{H}}_{o} t_{1}} . (25)

As we can see from Eq. (24), the photon emission rate becomes time-independent, as it should be in the case of the stationary fluorescence. The presence of a double commutator and two Hamiltonians indicate that two interactions with the field have to occur for the system to emit a photon. The integration variable t₁ in Eq. (24) represents the delay between these two interactions. This delay may be arbitrary, therefore, we have to perform integration over all possible delays spanning from 0 to infinity in order to obtain the stationary value. This explains the physics behind the integral in Eq. (24).

Now let us analyse the trace in Eq. (24). The trace identity (A.6) allows us to change the order of the operators under the trace as follows,

\begin{array}{l} {Tr}_{T} ({\hat{N}}_{k s} [{\hat{H}}_{k s}^{I} (t_{1}), [{\hat{H}}_{SF}, {\hat{ρ}}_{T}^{eq}]]) \\ = {({Tr}_{T} ({\hat{ρ}}_{T}^{eq} [{\hat{H}}_{SF}, [{\hat{H}}_{k s}^{I} (t_{1}), {\hat{𝒩}}_{k s}]]))}^{*} . (26) \end{array}

We now have a commutator of Ĥ_SF, which includes all modes of the field, with a single-mode operator [ ${\hat{H}}_{k s}^{I}$ (t₁), 𝒩̂_k_s]. Since operators corresponding to different modes mutually commute, we may use the single-mode Hamiltonian Ĥ_k_s instead of Ĥ_SF here and, consequently, in the initial Eq. (24). Also noting that, according to Eq. (25), ${\hat{H}}_{k s}^{I}$ (0) = Ĥ_k_s, we rewrite the trace in a slightly more consistent way,

T \equiv {Tt}_{T} ({\hat{N}}_{k s} [{\hat{H}}_{k s}^{I} (t_{1}), [{\hat{H}}_{k s}^{I} (0), {\hat{ρ}}_{T}^{eq}]]) . (27)

Upon expansion of the commutators, we will obtain four terms in total,

T = T_{1} + T_{2} + T_{3} + T_{4}, (28)

T_{1} = {Tt}_{T} ({\hat{𝒩}}_{k s} {\hat{H}}_{k s}^{I} (t_{1}) {\hat{H}}_{k s}^{I} (0) {\hat{ρ}}_{T}^{eq}), (29)

T_{2} = {Tt}_{T} (- {\hat{𝒩}}_{k s} {\hat{H}}_{k s}^{I} (t_{1}) {\hat{ρ}}_{T}^{eq} {\hat{H}}_{k s}^{I} (0)), (30)

T_{3} = {Tt}_{T} (- {\hat{𝒩}}_{k s} {\hat{H}}_{k s}^{I} (0) {\hat{ρ}}_{T}^{eq} {\hat{H}}_{k s}^{I} (t_{1})), (31)

T_{4} = {Tt}_{T} ({\hat{𝒩}}_{k s} {\hat{ρ}}_{T}^{eq} {\hat{H}}_{k s}^{I} (0) {\hat{H}}_{k s}^{I} (t_{1})) . (32)

These four terms may be presented graphically using the double-sided Feynman diagrams [11] shown in Fig. 1. The thick vertical lines represent the elements of the total density matrix, while the dots and the wavy arrows represent interactions with the field; time runs from below up. A Hamiltonian operator appearing to the left (right) of the density matrix in the formula corresponds to an interaction with the field of the left (right) vertical line in the diagram. An incoming wavy arrow represents photon absorption (annihilation), and an outgoing wavy arrow represents photon emission (creation). Hence, the directions of the wavy arrows are dictated by the state of the system and the field: if the field is in the vacuum state and the system is excited, then photon emission occurs; conversely, if the field is in a singly excited state and the system is in the ground state, photon absorption takes place. The diagrams in Fig. 1 depict all possible relaxation pathways (restricted to only two interactions with the field, as required in a single-photon process), and we encourage the reader to use the stated rules and check the correspondence between the diagrams and the equations. An important point here is that knowing the rules for constructing the diagrams we could have first drawn them and then immediately written down Eqs. (29–32). In fact, this is the way the two-sided Feynman diagrams are actually used in practice (just like the Feynman diagrams for the processes in the Standard Model [24]). For example, we can construct (29) from Fig. 1(a), by starting from ${\hat{ρ}}_{T}^{eq}$ , acting from the left with ${\hat{H}}_{k s}^{I}$ (0), and again acting from the left with ${\hat{H}}_{k s}^{I}$ (t₁). The final interaction is from the left with ${\hat{𝒩}}_{k s}$ , and we take the trace over all degrees of freedom. Other terms can be constructed similarly, with every interaction from the right introducing a minus sign.

We are only concerned with a single (k, s) mode in this discussion, and the initial density matrix of this mode is |0_k_s⟩⟨0_k_s|. The equilibrium density matrix of the system, on the other hand, is generally a certain superposition state,

{\hat{ρ}}_{S}^{e q} = \sum_{a, b} ({\hat{ρ}}_{S}^{eq}) a b | a ⟩ ⟨ b |, (33)

where |a⟩ is an eigenvector of Ĥ_S. (In contrast to the field state vectors, such as |0_k_s⟩, the system state vectors will not have any indices in our notation.) The total density matrix in the equilibrium is thus

\begin{array}{l} {\hat{ρ}}_{S}^{eq} = \sum_{a, b} {({\hat{ρ}}_{S}^{eq})}_{a b} | a ⟩ ⟨ b | \otimes | 0_{k s} ⟩ ⟨ 0_{k s} | \\ = \sum_{a, b} {({\hat{ρ}}_{S}^{eq})}_{a b} | 0_{k s} ⟩ | a ⟩ ⟨ b | ⟨ 0_{k s} |, (34) \end{array}

and each of the diagrams in Fig. 1 represents a single element of this sum. By also showing the field state in addition to the system state, we extend the usual Feynman diagram nomenclature.

The diagrams enable us to discard the terms that will not contribute to the fluorescence. For example, in Fig. 1(a) corresponding to T₁, we can see that only the ket field vector undergoes changes (since the Hamiltonians act on ${\hat{ρ}}_{T}^{eq}$ from the left in Eq. (29)), and the field remains in the ground state after the interactions. Therefore, the sequence of interactions shown in Fig. 1(a) will not contribute to fluorescence, so T₁ should equal to zero. A similar situation is seen in Fig. 1(d), which implies that T₄= 0. Only the terms T₂ and T₃ will contribute, where one interaction changes the ket and another one changes the bra, as shown in Figs. 1(b, c).

Fig. 1. The double-sided Feynman diagrams corresponding to Eqs. (29–32).

Looking closely at the terms in Eqs. (29–32) we notice that, since all the operators under the traces are Hermitian,

T_{1} = T_{4}^{*}, T_{2} = T_{3}^{*} . (35)

This follows from the trace identity (A.4). The relations given in Eq. (35) are also encoded in the diagrams in Fig. 1, because Figs. 1(a, d) as well as Figs. 1(b, c) are mirror images of each other.

We can now support our claim that T₁ = T₄ = 0 by a direct calculation. Substituting Eq. (18) into the expression for T₄ in Eq. (32), we obtain

T_{4} = {Tr}_{T} ({\hat{𝒩}}_{k s} | 0_{F} ⟩ ⟨ 0_{F} | {\hat{ρ}}_{S}^{eq} {\hat{H}}_{k s}^{I} (0) {\hat{H}}_{k s}^{I} (t_{1})) . (36)

Since

{\hat{𝒩}}_{k s} | 0_{F} ⟩ = {\hat{a}}_{k s}^{†} {\hat{a}}_{k s} | 0_{F} ⟩ = 0, (37)

we have

T_{1} = T_{4} = 0. (38)

Now we will turn to the evaluation of the term T₃ in Eq. (31). The first step is to analyse the Hamiltonian ${\hat{H}}_{k s}^{I}$ (t₁). We may immediately rewrite Eq. (11) in the interaction picture as

{\hat{H}}_{k s}^{I} (t_{1}) = - \sum_{α = 1}^{N_{mol}} {\hat{μ}}_{α}^{I} (t_{1}) \cdot (i N_{k} p_{k s} {\hat{a}}_{k s}^{† I} (t_{1}) e^{ik \cdot r_{α}} - i N_{k} p_{k s}^{*} {\hat{a}}_{k s}^{† I} (t_{1}) e^{-ik \cdot r_{α}}) . (39)

The transition dipole moment operator is an operator of the system subspace. Therefore, when switching to the interaction representation using Eq. (25), $\hat{μ}$ _α is not being acted upon by the operator e^(i/^ħ⁾^Ĥ^F^t¹. Similarly, the creation and annihilation operators of the field commute with e^(i/^ħ⁾^Ĥ^S^t¹. Hence, we may switch to the interaction picture using the formulas

{\hat{μ}}_{α}^{I} (t_{1}) = e^{(i/ ℏ) {\hat{H}}_{S} t_{1}} {\hat{μ}}_{α} e^{- (i/ ℏ) {\hat{H}}_{S} t_{1}}, (40)

and

{\hat{a}}_{ks}^{I} (t_{1}) = e^{(i / ℏ) {\hat{H}}_{F} t_{1}} {\hat{a}}_{ks} e^{- (i / ℏ) {\hat{H}}_{F} t_{1}}, (41)

{\hat{a}}_{k s}^{† I} (t_{1}) = e^{(i / ℏ) {\hat{H}}_{F} t_{1}} {\hat{a}}^{†}_{k s} e^{- (i / ℏ) {\hat{H}}_{F} t_{1}}, (42)

which remain compatible with the definition in Eq. (25). Equation (40) will not be elaborated any further in the present derivation since that requires selecting a particular model for the molecular system. On the other hand, the quantities ${\hat{a}}_{ks}^{I} (t_{1})$ and ${\hat{a}}_{k s}^{† I} (t_{1})$ may be readily calculated, as the forms of the field Hamiltonian and the creation and annihilation operators are universal. Since the field Hamiltonian commutes with itself, we have

\begin{array}{l} {\hat{H}}_{F} = {\hat{H}}_{F}^{I} (t_{1}) \\ = \sum_{k, s} ℏ ω_{k} ({\hat{a}}_{k s}^{† I} (t_{1}) {\hat{a}}_{k s}^{I} (t_{1}) + \frac{1}{2}), (43) \end{array}

and the Heisenberg equation for the evolution of ${\hat{a}}_{k s}^{I}$ (t₁) takes the simple form,

\begin{array}{l} \frac{d}{d t_{1}} {\hat{a}}_{k s}^{I} (t_{1}) = \frac{i}{ℏ} [{\hat{H}}_{F}, {\hat{a}}_{k s}^{I} (t_{1})] \\ = - i ω {\hat{a}}_{k s}^{I} (t_{1}) . (44) \end{array}

Hence

\begin{array}{l} {\hat{a}}_{k s}^{I} (t_{1}) = {\hat{a}}_{k s} e^{- i ω k^{t} 1}, \\ {\hat{a}}_{k s}^{† I} (t_{1}) = {\hat{a}}^{†}_{k s} e^{i ω k^{t} 1.} (45) \end{array}

Using these results, Eq. (39) becomes

\begin{array}{l} {\hat{H}}_{k s}^{I} (t_{1}) = - \sum_{α = 1}^{N_{mol}} {\hat{μ}}_{α}^{I} (t_{1}) \cdot (i N_{k} p_{k s} {\hat{a}}_{k s} e^{ik \cdot r_{α} - i ω_{k} t_{1}} \\ - i N_{k} p_{k s}^{*} {\hat{a}}_{k s}^{†} e^{- ik \cdot r_{α} +i ω_{k} t_{1}}) . (46) \end{array}

Substituting Eqs. (46) and (18) into Eq. (31), we now have for T₃,

\begin{array}{l} T_{3} = \sum_{α, β}^{N_{mol}} {(i N_{k})}^{2} {Tr}_{T} ({\hat{N}}_{k s} ({\hat{μ}}_{α}^{I} (0) \cdot p_{k s} {\hat{a}}_{k s} e^{ik \cdot r_{α}} - {\hat{μ}}_{α}^{I} (0) \cdot p_{k s}^{*} {\hat{a}}_{k s}^{†} e^{- ik \cdot r_{α}}) \\ \times | 0_{F} ⟩ ⟨ 0_{F} | (μ_{_{β}}^{I} (t_{1}) . p_{k s} {\hat{a}}_{k s} e^{ik \cdot r_{β} - i ω_{k} t_{1}} - {\hat{μ}}_{β}^{I} (t_{1}) \cdot p_{k s}^{*} {\hat{a}}_{k s}^{†} e^{- ik \cdot r_{β} +i ω_{k} t_{1}})) . (47) \end{array}

Once we expand the braces, we will have four terms under the trace, but only one of them will contribute. This can be seen in Fig. 1(c) from the direction of the wavy arrows that represent the interaction with the field. For the photon to be emitted, the initial field state |0_k_s⟩⟨0_k_s| has to change into |1_k_s⟩⟨1_k_s| via two interactions with the field which are represented by two outgoing wavy arrows. The creation operator ${\hat{a}}_{k s}^{†}$ should act from the left on the ket, and a conjugate interaction should take place, which corresponds to the annihilation operator â_k_s acting on the bra from the right. Therefore, out of four possible combinations of the directions of the two wavy arrows, only the one shown in Fig. 1(c) will contribute. In passing, we note that no combination of one-sided creation and annihilation of photons in Figs. 1(a, d) could lead to the |1_k_s⟩⟨1_k_s| state in those cases.

Expanding the braces, we obtain

\begin{aligned} T_{3} = N_{k}^{2} \sum_{α, β}^{N_{mol}} {Tr}_{T} ({\hat{N}}_{k s} ({\hat{μ}}_{α}^{I} (0) \cdot p_{k s}) \\ \times {\hat{a}}_{k s} | 0_{F} ⟩ ⟨ 0_{F} | {\hat{ρ}}_{S}^{e q} ({\hat{μ}}_{β}^{I} (t_{1}) \cdot p_{k s}) {\hat{a}}_{k s} e^{i k \cdot (r_{α} + r_{β}) - i ω_{k} t_{1}} \\ - {\hat{N}}_{k s} ({\hat{μ}}_{α}^{I} (0) \cdot p_{k s}) {\hat{a}}_{k s} | 0_{F} ⟩ ⟨ 0_{F} | {\hat{ρ}}_{S}^{e q} ({\hat{μ}}_{β}^{I} (t_{1}) \cdot p_{k s}^{*}) \\ \times {\hat{a}}_{k s}^{†} e^{i k \cdot (r_{α} - r_{β}) + i ω_{k} t_{1}} \\ - {\hat{N}}_{k s} ({\hat{μ}}_{α}^{I} (0) \cdot p_{k s}^{*}) {\hat{a}}_{k s}^{†} | 0_{F} ⟩ ⟨ 0_{F} | {\hat{ρ}}_{S}^{e q} ({\hat{μ}}_{β}^{I} (t_{1}) \cdot p_{k s}) \\ \times {\hat{a}}_{k s} e^{i k \cdot (- r_{α} + r_{β}) - i ω_{k} t_{1}} \\ + {\hat{N}}_{k s} ({\hat{μ}}_{α}^{I} (0) \cdot p_{k s}^{*}) {\hat{a}}_{k s}^{†} | 0_{F} ⟩ ⟨ 0_{F} | {\hat{ρ}}_{S}^{e q} ({\hat{μ}}_{β}^{I} (t_{1}) \cdot p_{k s}^{*}) \\ \times {\hat{a}}_{k s}^{†} e^{i k \cdot (- r_{α} - r_{β}) + i ω_{k^{'}} t}) (,, (48)) \end{aligned}

Since â_ks|0_F⟩ = 0 and ⟨0_F| ${\hat{a}}_{k s}^{†}$ = 0, we indeed notice that only the third term, containing ${\hat{a}}_{k s}^{†}$ |0_F⟩⟨0_F| â_ks, will contribute. In order to simplify the trace, we introduce a schematic notation

| 0_{F} ⟩ = | 0_{k s} ⟩ | 0_{F - k s} ⟩, (49)

meaning that the vacuum state of the field can be partitioned as a product of the vacuum state of the (k, s) mode and the vacuum state of the rest of the field, and

{Tr}_{F} {= Tr}_{k}_{s} {Tr}_{F–} {_{k}}_{s}, (50)

meaning that we can first perform the trace over all modes but (k, s) and only then trace over the (k, s) mode. Then,

\begin{aligned} {Tr}_{F} ({\hat{A}}_{k s} | 0_{F} ⟩ ⟨ 0_{F} | {\hat{B}}_{k s}) \\ = {Tr}_{k s} {Tr}_{F - k s} ({\hat{A}}_{k s} | 0_{k s} ⟩ ⟨ 0_{k s} | {\hat{B}}_{k s} | 0_{F - k s} ⟩ ⟨ 0_{F - k s} |) \\ = {Tr}_{k s} ({\hat{A}}_{k s} | 0_{k s} ⟩ ⟨ 0_{k s} | {\hat{B}}_{k s}) \times {Tr}_{F - k s} (| 0_{F - k s} ⟩ ⟨ 0_{F - k s} |) \\ = {Tr}_{k s} ({\hat{A}}_{k s} | 0_{k s} ⟩ ⟨ 0_{k s} | {\hat{B}}_{k s}) \\ = \sum_{m = 0}^{\infty} ⟨ m_{k s} | {\hat{A}}_{k s} | 0_{k s} ⟩ ⟨ 0_{k s} | {\hat{B}}_{k s} | m_{k s} ⟩ . (51) \end{aligned}

Here Â_k_s and B̂_k_s are arbitrary operators acting only on the (k, s) mode. Thus, we obtain from Eq. (48),

\begin{matrix} T_{3} = - N_{k}^{2} \sum_{α, β}^{N_{m o l}} {Tr}_{S} (\sum_{m = 1}^{\infty} ⟨ m_{k s} | {\hat{N}}_{k s} ({\hat{μ}}_{α}^{I} (0) \cdot p_{k s}^{*}) \\ \times {\hat{a}}_{k s}^{†} | 0_{k s} ⟩ ⟨ 0_{k s} | {\hat{ρ}}_{S}^{e q} ({\hat{μ}}_{β}^{I} (t_{1}) \cdot p_{k s}) \\ \times {\hat{a}}_{k s} e^{i k \cdot (- r_{α} + r_{β}) - i ω_{k} t_{1}} | m_{k s} ⟩ . (52) \end{matrix}

Inspecting the term ⟨0_ks| â_ks|m_ks⟩, we see that it will only contribute if m = 1. Hence,

\begin{aligned} T_{3} = - N_{k}^{2} \sum_{α, β}^{N_{mol}} {Tr}_{S} (⟨ 1_{k s} | {\hat{N}}_{k s} ({\hat{μ}}_{α}^{I} (0) \cdot p_{k s}^{*}) {\hat{a}}_{k s}^{†} | 0_{k s} ⟩ ⟨ 0_{k s} | \\ \times {\hat{ρ}}_{S}^{e q} ({\hat{μ}}_{β}^{I} (t_{1}) \cdot p_{k s}) {\hat{a}}_{k s} e^{i k \cdot (- r_{a} + r_{β}) - i ω_{k} t_{1}} | 1_{k s} ⟩) \\ = - N_{k}^{2} \sum_{α, β}^{N_{m o l}} e^{i k \cdot (- r_{a} + r_{β}) - i ω_{k} t_{1}} \\ \times {Tr}_{S} (({\hat{μ}}_{α}^{I} (0) \cdot p_{k s}^{*}) {\hat{ρ}}_{S}^{e q} ({\hat{μ}}_{β}^{I} (t_{1}) \cdot p_{k s})) (53) \end{aligned}

Using Eqs. (28), (35) and (38), we have

T = T_{2} + T_{3} = T_{3}^{*} + T_{3} = 2 ℜ (T_{3}) . (54)

Together with Eq. (53) this result allows us to write Eq. (24) as

\begin{aligned} f_{k s} & = 2 \frac{N_{k}^{2}}{ħ^{2}} \sum_{α, β}^{N_{m o l}} ℜ \int_{0}^{\infty} d t_{1} e^{i k \cdot (- r_{α} + r_{β})} e^{- i ω_{k} t_{1}} \\ \times {Tr}_{S} (({\hat{μ}}_{β}^{I} (t_{1}) \cdot p_{k s}) ({\hat{μ}}_{α}^{I} (0) \cdot p_{k s}^{*}) {\hat{ρ}}_{S}^{e q}) . (55) \end{aligned}

Since we are working in the dipole approximation and neglect all higher order effects, we may set e^{ik·(-r_α+r_β)} ≈ 1 [3]. Indeed, for the dipole–dipole interaction between the constituents of the molecular system to be significant, the distance Δr = |–r_α + r_β| between the dipoles usually must not exceed ~100Å. This is the case, for example, for chlorophylls in photosynthetic systems [25]. Meanwhile, the frequency of the electronic transitions is ω_k ≈ 3 fs^–1, so |k| = ω_k/c ≈ 10⁷ m^–1. Therefore, the upper limit for the product |k|Δr is 0.1 ≪ 1, which supports our claim. This lets us define the total system transition dipole moment operator as

\hat{μ} = \sum_{α = 1}^{N_{m o l}} {\hat{μ}}_{α}, (56)

leading to

\begin{aligned} f_{k s} = 2 \frac{N_{k}^{2}}{ħ^{2}} ℜ \int_{0}^{\infty} d t_{1} e^{- i w_{k} t_{1}} \\ \times {Tr}_{S} (({\hat{μ}}^{I} (t_{1}) \cdot p_{k s}) ({\hat{μ}}^{I} (0) \cdot p_{k s}^{*}) {\hat{ρ}}_{s}^{e q}) . (57) \end{aligned}

In order to obtain the total fluorescence spectrum, we will have to sum over the polarizations, and it is appropriate to do so now. The summation of the relevant part of Eq. (57) proceeds as follows,

\begin{aligned} \sum_{s} ({\hat{μ}}^{I} (t_{1}) \cdot p_{k s}) ({\hat{μ}}^{I} (0) \cdot p_{k s}^{*}) \\ = \sum_{s} \sum_{i, j = 1}^{3} {({\hat{μ}}^{I} (t_{1}))}_{i} {(p_{k s})}_{i} {({\hat{μ}}^{I} (0))}_{j} {(p_{k s}^{*})}_{j} \\ = \sum_{i, j = 1}^{3} {({\hat{μ}}^{I} (t_{1}))}_{i} {({\hat{μ}}^{I} (0))}_{j} (δ_{i j} - k_{0 i} k_{0 j}) \\ = {\hat{μ}}^{I} (t_{1}) \cdot {\hat{μ}}^{I} (0) - ({\hat{μ}}^{I} (t_{1}) \cdot k_{0}) ({\hat{μ}}^{I} (0) \cdot k_{0}) . (58) \end{aligned}

Here, notation such as (p_k_s)_i is to be understood as the ith component of the vector p_k_s; k₀ denotes the unit vector in the direction of k (which is the wave vector of the emitted wave), and we have used the relation Σ_s(p_k_s)_i( $P_{k s}^{*}$ )_j = (δ_ij – k₀_ik₀_j) [19]. Now we will account for the fact that the total system transition dipole moment may be arbitrarily oriented with respect to the wave vector of the emitted wave. Therefore, we perform orientational averaging over all possible orientations of the transition dipole moments and the wave vectors k,

\begin{aligned} {⟨ {\hat{μ}}^{I} (t_{1}) \cdot {\hat{μ}}^{I} (0) - ({\hat{μ}}^{I} (t_{1}) \cdot k_{0}) ({\hat{μ}}^{l} (0) \cdot k_{0}) ⟩}_{or} \\ = {\hat{μ}}^{I} (t_{1}) \cdot {\hat{μ}}^{I} (0) - {⟨ ({\hat{μ}}^{I} (t_{1}) \cdot k_{0}) ({\hat{μ}}^{I} (0) \cdot k_{0}) ⟩}_{or} \\ = {\hat{μ}}^{I} (t_{1}) \cdot {\hat{μ}}^{I} (0) - \frac{1}{3} {\hat{μ}}^{I} (t_{1}) \cdot {\hat{μ}}^{I} (0) \\ = \frac{2}{3} {\hat{μ}}^{I} (t_{1}) \cdot {\hat{μ}}^{I} (0) . (59) \end{aligned}

In the third line, we have used the results from Ref. [26]. Returning to Eq. (57), we will now have

\begin{aligned} {⟨ f_{k} ⟩}_{o r} \equiv \sum_{s} {⟨ f_{k s} ⟩}_{o r} \\ = \frac{4}{3} \frac{N_{k}^{2}}{ħ^{2}} ℜ \int_{0}^{\infty} d t_{1} e^{- i ω_{k} t_{1}} {Tr}_{s} ({\hat{μ}}^{I} (t_{1}) \cdot {\hat{μ}}^{I} (0) {\hat{ρ}}_{s}^{e q}) . (60) (60) \end{aligned}

If we wish to find the total number of photons in some frequency interval (ω₁, ω₂) that the system emits per unit time, we have to perform the summation over the corresponding wave vectors k,

\begin{aligned} \sum_{k} {⟨ f_{k} ⟩}_{o r} \\ = \frac{V}{(2 π)^{3}} \int_{k_{1}}^{k_{2}} d k \int_{0}^{π} d ϑ \int_{0}^{2 π} d φ k^{2} \sin ϑ {⟨ f_{k} ⟩}_{o r} \\ = \frac{V}{(2 π)^{3}} \int_{ω_{1}}^{ω_{2}} d ω_{k} \int_{0}^{π} d ϑ \int_{0}^{2 π} d φ \frac{ω_{k}^{2}}{c^{3}} \sin ϑ {⟨ f_{k} ⟩}_{o r} \\ = \int_{ω_{1}}^{ω_{2}} d ω_{k} \frac{V}{2 π^{2}} \frac{ω_{k}^{2}}{c^{3}} {⟨ f_{k} ⟩}_{o r . ((61))} \end{aligned}

In the second line, we changed the sum over k by an integral in the spherical coordinates using $\sum_{k} \to \frac{V}{{(2 π)}^{3}} \int d^{3} k$ in the limit of large V and then denoting |k| = k. In the third line, we changed the integration variable, k = ω_k/c. Finally, we recognized that ⟨f_k⟩_or does not depend on the direction of k (see Eq. (60)), which let us perform the integration over the angles with the result given in the fourth line.

The integrand in the last line in Eq. (61) represents the total emission rate of photons whose frequencies lie in the interval (ω_k, ω_k + dω_k) and therefore represents the fluorescence spectrum f(ω_k) of the system. Using Eq. (60) and dropping the index “k” on ω_k, we obtain

\begin{aligned} f (ω) = \frac{ω^{3}}{3 π^{2} ħ ε_{0} c^{3}} ℜ & {\int_{0}^{\infty} d t_{1} e^{- i ω t_{1}} \\ \times {Tr}_{S} ({\hat{μ}}^{I} (t_{1}) \cdot {\hat{μ}}^{I} (0) {\hat{ρ}}_{S}^{e q}), (62) \end{aligned}

completing the derivation.

The trace in Eq. (62) denotes the dipole–dipole correlation function,

C_{μ μ} (t) = {Tr}_{s} ({\hat{μ}}^{I} (t) \cdot {\hat{μ}}^{I} (0) {\hat{ρ}}_{S}^{e q}), (63)

calculated to the zeroth order in the system–field interaction. Consequently, as we show in the following section, the spectra calculated using Eq. (62) do not exhibit the radiative broadening. Nevertheless, the accuracy of the obtained formula turns out to be sufficient for simulating molecular systems interacting with the environment, since the broadening of the spectral lines due to the system–environment interaction is not neglected in Eq. (62), and it is several orders of magnitude greater than the radiative broadening [8].

In general, calculation of the dipole–dipole correlation function requires application of dynamical theories for open quantum systems. However, if we focus on some simple model system, we may calculate C_μμ(t) without any advanced theories. As an example, let us consider an isolated two-level system.

4. Emission rate of a two-level system

The system Hamiltonian of a general two-level system is given by

{\hat{H}}_{S} = (\begin{array}{ll} ε_{1} & 0 \\ 0 & ε_{2} \end{array}), (64)

where ε_n is the energy of the nth energy level, and we set ε₁ > ε₂. The transition dipole moment matrix may be taken in the form

\hat{μ} = (\begin{array}{ll} 0 & μ_{12} \\ μ_{12}^{*} & 0 \end{array}), (65)

where the zero-valued diagonal elements express the absence of a permanent static dipole moment in the system. This is a valid simplification, as the static dipole moment has no impact on the spontaneous emission process. Finally, the state corresponding to the excited state equilibrium is described by the density matrix

{\hat{ρ}}_{S}^{e q} = (\begin{array}{ll} 1 & 0 \\ 0 & 0 \end{array}), (66)

which shows that the system occupies the higher energy level with probability 1. Equations (64–66) completely define the system, and the dipole–dipole correlation function may be readily calculated from Eq. (63). As we show in Appendix C, the result is

C_{μ μ} (t_{1}) = {| μ_{12} |}^{2} e^{i ω_{12} t_{1}}, (67)

and the integration in Eq. (62) yields

f (ω) = \frac{ω_{12}^{3}}{3 π ħ ε_{0} c^{3}} {| μ_{12} |}^{2} δ (ω - ω_{12}) . (68)

As expected from the application of Eq. (62), the obtained spectral line is not of a Lorentzian profile with a finite natural linewidth, but is an infinitely narrow line instead, which is a deficiency of this level of approximation. However, the coefficient in the obtained emission rate formula is precisely the well-known Einstein A coefficient,

A_{12} = \frac{ω_{12}^{3}}{3 π ħ ε_{0} c^{3}} {| μ_{12} |}^{2}, (69)

that was originally used to phenomenologically describe the process of spontaneous emission in a two-level system [19]. Thus, the derived stationary fluorescence spectrum formula (62) reduces to the expected result in the limit of the simplest quantum system.

5. Conclusions

In this work, we have presented a systematic derivation of the relation between the fluorescence spectrum of a condensed matter system and the relevant dipole–dipole correlation function. This requires describing the electromagnetic field as a quantum object, which can be conveniently done using the second quantization notation. Although conceptually and mathematically rather intricate, such a description is highly advantageous because it allows one to operate directly with the mode occupation operator. Consequently, the rate of change of mode occupation also has a corresponding quantum operator, and the theory thus naturally accounts for possible creation and emission of photons. This process cannot be adequately described if the electromagnetic field is treated classically; the only way to account for it without using the concepts of quantum electrodynamics is to postulate that spontaneous relaxation followed by the emission of a photon is possible. By contrast, the present derivation clearly shows that emission of photons is a result of the interaction between the quantum system and the vacuum state of the electromagnetic field.

The obtained formula for the fluorescence spectrum is a rather general one since it does not rely on any specific model of the quantum system. This generality was achieved by employing the density matrix formalism, which is mandatory when working with open quantum systems. In fact, the density matrix formalism proved to be equally convenient when studying a simple two-level system. As we have shown, the emission rate formula reduces to the Einstein A coefficient in that special case.

Acnowledgements

This work was supported by the Research Council of Lithuania (LMT Grant No. S-MIP-20-44). Y.B. acknowledges funding by the European Social Fund under Measure No. 09.3.3-LMT-K-712 ‘Development of Competences of Scientists, other Researchers and Students through Practical Research Activities’.

Appendices

A. Trace identities

In this section, we prove several trace identities used in the main text. Throughout this work, we will be using extensively the cyclic invariance of the trace,

Tr (\hat{A} \hat{B} \hat{C}) = Tr (\hat{C} \hat{A} \hat{B}) = Tr (\hat{B} \hat{C} \hat{A}), (A.1)

where $\hat{A}$ , $\hat{B}$ and $\hat{C}$ are arbitrary, generally non-commuting matrices.

The first identity is

Tr (\hat{A} [\hat{B}, \hat{C}]) = Tr (- \hat{B} [\hat{A}, \hat{C}]) . (A.2)

This is trivially proved by expanding the commutator,

\begin{aligned} Tr (\hat{A} [\hat{B}, \hat{C}]) = Tr (\hat{A} \hat{B} \hat{C} - \hat{A} \hat{C} \hat{B}) \\ = Tr (\hat{C} \hat{A} \hat{B} - \hat{A} \hat{C} \hat{B}) \\ = Tr ([\hat{C}, \hat{A}] \hat{B}) \\ = Tr (- \hat{B} [\hat{A}, \hat{C}]) . (A.3) \end{aligned}

This property of the trace is similar to the corresponding property of the scalar triple product of vectors, which is cyclically invariant as well.

Next, we prove that

Tr (\hat{L} \hat{M} \hat{N} \hat{P}) = (Tr (\hat{P} \hat{N} \hat{M} \hat{L}))^{*}, (A.4)

where $\hat{L}$ , $\hat{M}$ , $\hat{N}$ and $\hat{P}$ are arbitrary Hermitian operators. The proof proceeds as follows,

\begin{aligned} Tr (\hat{L} \hat{M} \hat{N} \hat{P}) = Tr ({\hat{L}}^{†} {\hat{M}}^{†} {\hat{N}}^{†} {\hat{P}}^{†}) \\ = Tr ((\hat{P} \hat{N} \hat{M} \hat{L})^{†}) \\ = Tr ((\hat{P} \hat{N} \hat{M} \hat{L})^{*}) \\ = (Tr (\hat{P} \hat{N} \hat{M} \hat{L}))^{*} . (A.5) \end{aligned}

In the third line, we recognized that Hermitian conjugate of an operator is the transpose of its complex conjugate, but the trace, being the sum of the diagonal elements, is invariant with respect to transposition.

Finally, we will show that

Tr (\hat{L} [\hat{M}, [\hat{N}, \hat{P}]) = (Tr (\hat{P} [\hat{N}, [\hat{M}, \hat{L}]))^{*} . (A.6)

Starting by expanding all the commutators, we have

\begin{aligned} Tr (\hat{L} [\hat{M}, [\hat{N}, \hat{P}]) = (Tr (\hat{L} \hat{M} \hat{N} \hat{P} - \hat{L} \hat{M} \hat{P} \hat{N} \\ - (\hat{L} \hat{N} \hat{P} \hat{M} - \hat{L} \hat{P} \hat{N} \hat{M})) \\ = (Tr (\hat{P} \hat{N} \hat{M} \hat{L} - \hat{N} \hat{P} \hat{M} \hat{L} - (\hat{M} \hat{P} \hat{N} \hat{L} - \hat{M} \hat{N} \hat{P} \hat{L})))^{*} \\ = (Tr (\hat{P} \hat{N} \hat{M} \hat{L} - \hat{P} \hat{M} \hat{L} \hat{N} - (\hat{P} \hat{N} \hat{L} \hat{M} - \hat{P} \hat{L} \hat{M} \hat{N}))) \\ = (Tr (\hat{P} \hat{N} \hat{M} \hat{L} - \hat{P} \hat{N} \hat{L} \hat{M} - (\hat{P} \hat{M} \hat{L} \hat{N} - \hat{P} \hat{L} \hat{M} \hat{N}))) \\ = (Tr (\hat{P} \hat{N} [\hat{M}, \hat{L}] - \hat{P} [\hat{M}, \hat{L}] \hat{N}))^{*} = \end{aligned}

= {(Tr (\hat{P} [\hat{N}, [\hat{M}, \hat{L}]))}^{*} . (A .7)

In the third line, we used the identity (A.4); in the fourth line, we used the cyclic invariance to make operator $\hat{P}$ appear first in every term; in the fifth line, we simply swapped the second and the third terms.

B. First order solution for the density matrix

Below, we provide the steps leading to the formal solution for the first order term of the density matrix, as quoted in Eq. (22).

Our aim is to find the solution of Eq. (21). We may use the result of Eq. (20) and write

\frac{d}{d t} {\hat{ρ}}_{T}^{(1)} (t) = - \frac{i}{ħ} [{\hat{H}}_{0}, {\hat{ρ}}_{T}^{(1)} (t)] - \frac{i}{ħ} [{\hat{H}}_{S F}, {\hat{ρ}}_{T}^{e q}] . (B.1)

We will apply the method of variation of constants. The general solution of the corresponding homogeneous equation,

\frac{d}{d t} {\hat{\tilde{ρ}}}_{T}^{(1)} (t) = - \frac{i}{ħ} [{\hat{H}}_{0}, {\hat{\tilde{ρ}}}_{T}^{(1)} (t)], (B.2)

is easily seen to be

{\hat{\tilde{ρ}}}_{T}^{(1)} (t) = e^{- (i / ħ) {\hat{H}}_{0} t} \hat{C} e^{(i / ħ) {\hat{H}}_{0} t}, (B.3)

where $\hat{C}$ is some constant matrix. Next, we promote that matrix to a time-dependent one, thus seeking for the solution of (B.1) in the form

{\hat{ρ}}_{T}^{(1)} (t) = e^{- (i / ℏ) {\hat{H}}_{0} t} \hat{C} (t) e^{(i / ℏ) {\hat{H}}_{0} t} . (B .4)

Substituting into Eq. (B.1), we obtain

e^{- (i / ħ) {\hat{H}}_{0} t} \frac{d \hat{C} (t)}{d t} e^{(i / ħ) {\hat{H}}_{0} t} = - \frac{i}{ħ} [{\hat{H}}_{S F}, {\hat{ρ}}_{T}^{e q}], (B.5)

and the formal solution is

\hat{C} (t) = \hat{C} (t_{0}) - \frac{i}{ħ} \int_{t_{0}}^{t} d t^{'} e^{(i / ħ) {\hat{H}}_{0} t^{'}} [{\hat{H}}_{S F}, {\hat{ρ}}_{T}^{e q}] e^{- (i / ħ) {\hat{H}}_{0} t^{'}} . (B .6)

In the equilibrium – at t = t₀, when the perturbation is absent – the system density matrix is fully described by the zeroth order term, ${\hat{ρ}}_{T}^{(0)}$ , and all higher order terms are equal to zero. Hence, ${\hat{ρ}}_{T}^{(1)}$ = 0, which implies $\hat{C}$ (t₀) = 0. Taking this into account and substituting Eq. (B.6) into Eq. (B.4), we finally arrive at

\begin{aligned} {\hat{ρ}}_{T}^{(1)} (t) \\ = - \frac{i}{ħ} \int_{t_{0}}^{t} d t^{'} e^{- (i / ħ) {\hat{H}}_{0} (t - t^{'})} [{\hat{H}}_{S F}, {\hat{ρ}}_{T}^{e q}] e^{(i / ħ) {\hat{H}}_{0} (t - t^{'})} \end{aligned} (., (B.7))

It will prove to be more convenient to work with a slightly modified expression, which we obtain by performing a change of variables, t–t' = t₁,

{\hat{ρ}}_{T}^{(1)} (t) = - \frac{i}{ħ} \int_{0}^{t - t_{0}} d t_{1} e^{- (i / ħ) {\hat{H}}_{0} t_{1}} [{\hat{H}}_{S F}, {\hat{ρ}}_{T}^{e q}] e^{(i / ħ) {\hat{H}}_{0} t_{1}} \cdot (B .8)

C. Detailed calculation of the emission rate of a two-level system

Here, we show at full length how one may calculate the emission rate of a two-level system defined by Eqs. (64–66) starting from Eq. (62).

The first step is to find the dipole moment operator in the interaction picture,

{\hat{μ}}^{I} (t_{1}) = e^{(i / ħ) {\hat{H}}_{S} t_{1}} \hat{μ} e^{- (i / ħ) {\hat{H}}_{s} t_{1}} . (C.1)

Using the expansion

e^{\hat{A}} = \sum_{n} e^{A_{n}} | n ⟩ ⟨ n |, (C.2)

where $\hat{A}$ is an arbitrary operator, A_n is the nth eigenvalue and |n⟩ is the nth eigenvector, we have

\begin{aligned} {\hat{μ}}^{I} (t_{1}) & = (\begin{array}{cc} e^{(i / ħ) ε_{1} t_{1}} & 0 \\ 0 & e^{(i / ħ) ε_{2} t_{1}} \end{array}) (\begin{array}{lc} 0 & μ_{12} \\ μ_{12}^{*} & 0 \end{array}) (\begin{array}{cc} e^{- (i / ħ) ε_{1} t_{1}} & 0 \\ 0 & e^{- (i / ħ) ε_{2} t_{1}} \end{array}) \\ = (\begin{array}{cc} 0 & μ_{12} e^{i ω_{22} t_{1}} \\ μ_{12}^{*} e^{- i ω_{12} t_{1}} & 0 \end{array}), ((C.3)) \end{aligned}

where we have defined the transition frequency, ω₁₂ = (ε₁–ε₂)/ћ. We note that it is a well-known result that each element (m, n) of an arbitrary time-independent operator simply acquires a factor of eⁱ^ωmnt in the Heisenberg picture [21]. (When viewing the calculation of Eq. (C.1) as an isolated problem, we may consider Ĥ_S to be the total Hamiltonian of the system, so that the interaction picture coincides with the Heisenberg picture.)

Noting that $\hat{μ}$ = ${\hat{μ}}_{}^{I}$ (0), performing the matrix multiplication and taking the trace, we obtain the dipole–dipole correlation function,

\begin{aligned} C_{μ μ} (t_{1}) & \equiv {Tr}_{S} ({\hat{μ}}^{I} (t_{1}) \cdot {\hat{μ}}^{I} (0) {\hat{ρ}}_{S}^{e q}) \\ = {| μ_{12} |}^{2} e^{i ω_{12} t_{1}} \end{aligned} . ((C.4))

Calculation of the integral in Eq. (62) then proceeds as follows,

\begin{aligned} ℜ \int_{0}^{\infty} d t_{1} e^{- i ω t_{1}} C_{μ μ} (t_{1}) \\ = {| μ_{12} |}^{2} R \int_{0}^{\infty} d t_{1} e^{- i (ω - ω_{12}) t_{1}} \\ = {| μ_{12} |}^{2} \int_{0}^{\infty} d t_{1} \cos ((ω - ω_{12}) t_{1}) \\ = {| μ_{12} |}^{2} \frac{1}{2} \int_{- \infty}^{\infty} d t_{1} \cos ((ω - ω_{12}) t_{1}) \\ = {| μ_{12} |}^{2} π δ (ω - ω_{12}) . ((C.5)) \end{aligned}

In the last line, we used the integral representation of the Dirac delta function. Substituting this result into Eq. (62), we arrive at the final expression for the emission rate,

f (ω) = \frac{ω_{12}^{3}}{3 π ħ ε_{0} c^{3}} {| μ_{12} |}^{2} δ (ω - ω_{12}), ((C.6))

where the presence of the delta function allowed us to change 𝜔³ to 𝜔³₁₂.

It can be straightforwardly verified that including arbitrary static dipole moments μ₁₁ and μ₂₂ leads to the same results.

DERIVATION OF THE STATIONARY FLUORESCENCE SPECTRUM FORMULA FOR MOLECULAR SYSTEMS FROM THE PERSPECTIVE OF QUANTUM ELECTRODYNAMICS

1. Introduction

2. Preliminaries

3. Derivation

4. Emission rate of a two-level system

5. Conclusions

Acnowledgements

Appendices

A. Trace identities

B. First order solution for the density matrix

C. Detailed calculation of the emission rate of a two-level system

References

MOLEKULINIŲ SISTEMŲ STACIONARIOSIOS FLUORESCENCIJOS SPEKTRO FORMULĖS IŠVEDIMAS NAUDOJANT KVANTINĘ ELEKTRODINAMIKĄ