Theory of single molecule emission spectroscopy.

THE JOURNAL OF CHEMICAL PHYSICS 142, 174104 (2015)

Theory of single molecule emission spectroscopy Golan Bel1,2,3,a) and Frank L. H. Brown3,b)

1

Department of Solar Energy and Environmental Physics, Jacob Blaustein Institutes for Desert Research, Ben-Gurion University of the Negev, Sede Boqer Campus 84990, Israel 2 Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA 3 Department of Chemistry and Biochemistry and Department of Physics, University of California, Santa Barbara, California 93106, USA

(Received 21 January 2015; accepted 8 April 2015; published online 4 May 2015) A general theory and calculation framework for the prediction of frequency-resolved single molecule photon counting statistics is presented. Expressions for the generating function of photon counts are derived, both for the case of naive “detection” based solely on photon emission from the molecule and also for experimentally realizable detection of emitted photons, and are used to explicitly calculate low-order photon-counting moments. The two cases of naive detection versus physical detection are compared to one another and it is demonstrated that the physical detection scheme resolves certain inconsistencies predicted via the naive detection approach. Applications to two different models for molecular dynamics are considered: a simple two-level system and a two-level absorber subject to spectral diffusion. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4918709]

I. INTRODUCTION

Single molecule spectroscopy (SMS) is a widely used and versatile experimental technique in physics, chemistry, and biology.1–12 Often, SMS studies (experimental and theoretical)9,12–50 focus on number statistics for the photons emitted from a molecule under laser irradiation. These statistics may literally count individual photons, or may examine fluctuations in time-dependent intensities. In either case, significantly more information is obtained than would be possible from a simple absorption experiment on a bulk sample. In the majority of SMS studies, photons are detected indiscriminately, without frequency (or polarization or propagation direction) resolution. However, SMS experiments that resolve emitted photons by color have the potential to reveal even more information from spectroscopically active systems than do their broadband counterparts. Perhaps the best known example of multi-color detection schemes in SMS involves the exploitation of resonance energy transfer as a molecular ruler,8,26,27,51 where two differently colored dyes are attached to a molecule or super molecular assembly to monitor system dynamics via the simultaneous detection of photons emitted from both dyes. Similar studies can also be performed with semiconductor quantum dots52 and more elaborate detection schemes have been developed to extend two-channel photon detection to multiple channels53 with an eye toward the study of complex condensed phase and biologically relevant systems. Yet, even the simplest physical systems can yield surprisingly complex structure in their emission spectrum, indicating the utility of achieving frequency resolution in the emission spectroscopy for any SMS system. For example, the emission spectrum from a simple two-level absorber was predicted by Mollow,54 and subsequently experimentally observed in various systems,55–60 a)Electronic address: [email protected] b)Electronic address: [email protected]

0021-9606/2015/142(17)/174104/42/$30.00

to display a triplet structure when driven sufficiently strongly or off-resonance. Theoretical approaches building upon the pioneering work of Mollow have extended predictions for the emission spectrum of a two-level absorber to the calculation of frequency dependent intensity fluctuations and correlations.61–64 However, these studies have focused primarily on analytical approaches that are confined to certain limiting regimes and are applicable only to specific observables. Much of this work has been motivated by the beautiful experimental work of Aspect et al.,58 where strong temporal correlations were observed between opposite sidebands in the off-resonance excitation of strontium atoms. As important as this pioneering work is, the experiment was designed to allow for the simplest possible theoretical description-based upon a two-level emitter and widely separated emission peaks that can be unambiguously assigned to individual detectors. It is not immediately clear how the theoretical tools developed for this idealized situation could be carried over to more complex condensed phase systems that are typically studied via SMS. In Ref. 37, we suggested an alternate approach to calculating photon emission statistics that is not limited to two level systems. However, that approach is limited to situations where all emitted photons are associated with individual spectral transitions that are sufficiently wellseparated to be unambiguously associated with corresponding detector channels. The approach can be viewed as an extension of generating function methods65 developed to interpret singlepair fluorescence resonance energy transfer (FRET) measurements,26,27 but with the ability to include quantum dynamics of the system and multiple detection channels. The methods discussed above are not applicable to situations involving poorly resolved spectral lines and a general treatment of photon emission statistics should be able treat such a case. In a preliminary report,66 we presented a somewhat general scheme for the calculation of frequency resolved photon counting statistics that is capable of dealing with poorly

142, 174104-1

© 2015 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Fri, 08 May 2015 07:26:14

174104-2

G. Bel and F. L. H. Brown

J. Chem. Phys. 142, 174104 (2015)

resolved spectral lines and multi-state dynamics. Indeed, the only approximations of this approach are the assumption of Markovian dynamics for the driven chromophore system and the assumption that moments of the photon number operator for various field modes correspond to observable physical quantities. In this work, we will refer to the methodology of Ref. 66 as involving the “emission based” approach to photon statistics. The purpose of this paper is three-fold. First, we present details and an extended discussion of our emission based treatment introduced in Ref. 66. Second, we demonstrate that the emission based approach can lead to negative photon intensities at short times, which is nonsensical in the context of experimental measurements. Third, we generalize the emission based approach to explicitly account for the response of physically realizable detectors, thereby removing the problem of negative intensities and providing concrete predictions for direct comparison to experiments. The paper is organized as follows. In Sec. II, we introduce the generating function for photons emitted into discrete frequency bins and the formal relationship between this construct and various photon counting moments. In Sec. III, we introduce a Hamiltonian for the combined chromophore system, semi-classical continuous wave (CW) exciting field, and quantum radiation field describing the emitted photons. The dynamics of this composite system are solved to present the generating function and moments of Sec. II in terms of timeordered integrals of correlation functions of the chromophore system operators. Section IV presents an explicit scheme for the calculation of these correlation functions under the Markov approximation for chromophore dynamics. Though the correlation functions themselves are relatively simple to express in terms of simple matrix operations, the time-ordering operations inherent to the moments are tedious and the final expressions are derived and presented in Appendix A. The generalization of the foregoing results to account for photon detection (as opposed to photon emission) is presented in Sec. V; these results are closely related to the emission based results and involve the same correlation functions of chromophore operators. However, the required time-integrations are distinct from the emission based scheme and the final expressions are presented in Appendix B. Sections VI and VII present applications of the theory to model two-level and four-level chromophore systems, respectively. In Sec. VIII we conclude.

II. THE GENERATING FUNCTION

Our starting point is the Heisenberg representation of the → − number operator for photons with wave vector k and polari→ − zation ε , which is defined using the creation and annihilation operators as † (t)  Nk,ε (t) =  ak,ε ak,ε (t) .

(1)

† Here,  ak,ε , ak,ε are the creation and annihilation operators, → − respectively, for photons with wave vector k and polarization → −ε and t is the time.

The moment generating function is defined in a manner analogous to the usual classical definition, namely,

G (⃗s,t) ≡





N (t) sk,εk, ε

.

(2)

k,ε

The auxiliary variables sk,ε are introduced to facilitate extraction of photon statistics and the vector ⃗s ≡ {sk,ε }, on the left hand side of Eq. (2), is simply a shorthand notation to indicate the collective set of auxiliary variables associated with all possible photon wave vectors and polarizations. The averaging operation denoted by the angular brackets indicates a quantum mechanical average over the initial density matrix of the combined molecule/radiation field system that will be further specified in the following section, i.e., ⟨. . .⟩ ≡ Trace {. . . ρ (0)}. From the definition of the generating function, it follows that the factorial moments of the number operator Nk,ε (t) are obtainable by differentiating the generating function and setting ⃗s = 1. To be more explicit, the n th (n is a positive integer) factorial moment is given by  (n) 

Nk,ε (t) ≡ Nk,ε (t) (Nk,ε (t) − 1) . . . (Nk,ε (t) − n + 1) ∂ n G (⃗s,t) (3) = . n ∂sk,ε ⃗s =1 Similarly, the multivariate factorial moments are given by 

 ∂ n+m G (⃗s,t) (n) , (t) Nk(m) (t) Nk,ε = ′ε ′ n ∂skm′ε ′ ∂sk,ε ⃗s =1

(4)

(n and m are positive integers) with extension to moments involving more than two photon modes following immediately by introducing more derivatives. Making use of the identity x y = e y ln x , Eq. (2) can be recast as    * + † . /  G (⃗s,t) = exp ak,ε (t) ak,ε (t) ln (sk,ε ) , k,ε    * + †  = N exp . ak,ε (t) ak,ε (t) (sk,ε − 1)/ . , k,ε -

(5)

The second equality in Eq. (5) follows from a standard operator identity67 and employs the normal ordering operator, N . The normal ordering operator acts on all creation and annihilation operators appearing to its right and rearranges them by placing the creation operators to the left of the annihilation operators. When the generating function is written in this form, the derivatives indicated in Eq. (4) immediately lead to the rather simple expression for the factorial moments  (n)  Nk,ε (t) Nk(m) ′ε ′ (t) )n( † ( † )m ( ) n =  ak,ε (t)  ak′ε′ (t) ( ak′ε ′ (t))m  ak,ε (t) . (6) The formalism presented in this section is completely general and applies to all photons in the system, regardless of their origin (e.g., photons from vacuum fluctuations, photons emitted from chromophores, photons from an applied field, etc.). Our primary concern in this work is with photons emitted from a driven chromophore; Sec. III details the model used to study these photons and the manipulations necessary to cast the equations of this section into a form suitable for practical calculations.


174104-3


J. Chem. Phys. 142, 174104 (2015)

)  → − −µ ( D + + D − , µ =→ 0

III. MODEL SYSTEM AND ITS DYNAMICS

The content of Eqs. (7)–(20) below should be viewed as the natural extension of well-known results in quantum optics for two-level chromophores to multi-level chromophore systems. These results are presented in concise fashion, primarily to introduce notation and precisely define our system; more detailed discussions (in the context of two-level chromophores) are available in Refs. 62, 68, and 69. The Hamiltonian,  H(t), of the systems considered in this work is of the form  =H sys + H I e (t) + H I + H R . H(t)

(7)

sys is the unperturbed Hamiltonian of the molecule-environH ment system that is being irradiated by a laser source and I e (t) stands for the is thus driven into emitting photons. H interaction between the system and the applied laser field; this interaction is treated semi-classically leading to an explicit I describes the time dependence within the Hamiltonian. H interaction of the atom/molecule with the quantized radiation R is the Hamiltonian of the quantized radiation field. field. H We will focus on model systems consisting of two electronic levels (ground |g⟩ and excited |e⟩), so that sys = |g⟩ Hg ⟨g| + |e⟩ He ⟨e| H b + V ≡H CH + H b + V , +H

(8)

(12)

+/ D − are the raising/lowering operators between chrowhere D mophore excited and ground states, with matrix elements + m D = ⟨me ng , e ;n g (13)

n− ;m = ng |me ⟩ . D g e ± operators live entirely within the CH We stress that the D portion of the total system Hilbert space. We also introduced the Rabi frequency → − → E 0 · −µ 0 . (14) Ω0 = ~ We stress that the semi-classical treatment for the applied field described above should not be viewed as an approximation. Equation (11) can be derived from a fully quantum-mechanical treatment of the incident field without approximation.70,71 The semi-classical representation of the interaction is very convenient for the analysis that follows. The system is also coupled to the quantum radiation field. This interaction is written as → −  − I = −→ H µ · E q (0) ( ) ( ) (→ † −µ · → −) + + D − = −i D Ek  ak,ε −  ak,ε (15) 0 ε . k,ε

Hg and He are the Hamiltonians for nuclear motion within the ground and excited states of the chromophore (CH), respectively, with eigenfunctions and eigenvalues specified by Hg ng = ~ω n g ng , (9) He |me ⟩ = ~ω m e |me ⟩ for ng = 1, . . . Nmaxg and me = 1, . . . Nmaxe . We will be concerned only with a finite number of states Nmaxg and Nmaxe . b is the Hamiltonian of the thermal bath (environment) and V  H describes the interaction of the chromophore with the environment. Note that we use the identifier “sys” to specify that part of the total Hamiltonian associated with the “system,” as opposed to the quantum radiation field and applied perturbing field. The designation “CH” reflects that sub-portion of the system (the chromophore itself) that is actually responsible for interaction with radiation. Those portions of the system not associated with the chromophore itself constitute a reservoir/ bath/environment that may be capable of inducing relaxation and thermal fluctuation among the various chromophore states. The system is driven by a CW plane-polarized external → − laser field characterized by wave vector k L and frequency → − ω L = c| k L | such that ) → → − (→ − →− →− E e −r ,t = E 0ei k L · r cos (ω Lt) . (10) −r = 0 and treating Assuming that the molecule is located at → the interaction within the dipole approximation, we are led to ) → − (−  − I e = −→ H µ · Ee → r = 0,t ( ) + + D − cos (ω Lt) . = −~Ω0 D (11) On the right hand side of the above equation, we have used the dipole moment operator, calculated in the Condon approximation, such that



In the above, Ek = 2ε Vk , where ε 0 is the vacuum permit0 0 tivity and V0 is the volume of the cubic box used to quantize the field (in what follows, the limit V0 → ∞ will be taken such that V0 does not appear in any final results). The radiation Hamiltonian is ( )  1 † R =  H ~ωk  ak,ε ak,ε + . (16) 2 k,ε ~ω

A. The dynamics of the creation and annihilation operators and their relation to the system dipole moment operators

In what follows, it is convenient to work with the creation and annihilation operators in the interaction picture ak,ε (t) ≡  ak,ε (t) eiω k t , † † (t) ≡  (t) e−iω k t , ak,ε ak,ε

(17)

→ − where ωk = c k , c is the speed of light, and the operators † (t) are the Heisenberg operators associated  ak,ε (t) and  ak,ε † , respectively. The number operator of Eq. (1) with  ak,ε and  ak,ε † (t) ak,ε (t) in the interacconveniently remains Nk,ε (t) = ak,ε tion representation. Similarly, any operator involving only † pairs of  ak,ε and  ak,ε (e.g., Eq. (6)) retains its functional form when moving to the interaction representation due the cancellation of the positive and negative complex exponentials. Further, it is convenient to introduce the slowly varying rotating frame operators for the system operators ± . D ± ≡ e∓iω L t D

(18)

In Sec. II, we expressed the generating function for the number of photons in terms of photon creation and annihilation operators. Our goal is to express all photon counting moments


174104-4


J. Chem. Phys. 142, 174104 (2015)

in terms of the system operators alone. For this purpose, we derive the dynamics of the field operators in relation to the dipole moment operators of the system. The dynamics of the creation and annihilation operators are given by  dak,ε (t) 1   (t) + iωk ak,ε (t) = ak,ε (t) , H dt i~ (− → ) ) Ek → µ 0 · −ε ( = D +ei(ω k +ω L )t + D −ei(ω k −ω L )t , ~ †  (19) dak,ε (t) 1  †  (t) − iωk a† (t) = ak,ε (t) , H k,ε dt i~ (− → ) ) Ek → µ 0 · −ε ( = D +ei(−ω k +ω L )t + D −e−i(ω k +ω L )t , ~ which follow immediately by applying the Heisenberg equations of motion to Eq. (17). Applying the Rotating Wave Approximation (RWA) to the equations above (i.e., neglecting terms oscillating with the laser frequency, or higher) yields the formal solutions (− → ) t Ek → µ 0 · −ε  D − (τ) eiω k Lτ dτ, ak,ε (t) = ak,ε (0) + ~ 0

† (t) ak,ε

=

† (0) ak,ε

+

(− → ) t Ek → µ 0 · −ε  ~

(20) D + (τ) e−iω k Lτ dτ.

0

In the above equations, ωk L ≡ ωk − ω L .

B. The generating function in terms of system operators

In principle, the results of Eq. (20) could be substituted directly into the expressions of Sec. II to derive expressions for photon counting statistics. Unfortunately, the field operators at time zero do not commute with the system operators at positive times and this naive substitution results in cumbersome theoretical expressions involving both field operators at time zero and system operators over all positive times entangled together. We can simplify matters considerably by recognizing that Eq. (5) can be rewritten as   ∂  ∂ G (⃗s,t) = exp  (sk,ε − 1) ∂pk,ε ∂r k,ε   k,ε      † p a (t) × e k, ε k, ε k, ε e k, ε rk, ε ak, ε (t) →− →− , (21) p = r =0 where the variables pk,ε and r k,ε have been introduced as a mathematical device to extract the required number of creation and annihilation operators via differentiation; they do not appear in the generating function as they are set to zero following the indicated differentiations. Writing the generating function in this way enforces the normal ordering operation explicitly. Looking at Eq. (20), we see that the time derivatives of ak,ε (t) commute with ak,ε (t) since the derivative involves only the system operator D −(t) and no field components; similarly for the creation operators and associated time derivative. This implies that we may write

(→ −µ · → −ε ) ,   E † k 0 ∂ k, ε pk, ε ak,† ε (t) p a (t) e = e k, ε k, ε k, ε *. D +(t)e−iω k L t +/ , pk,ε ∂t ~ , k,ε ) (→ −µ · → −ε  E   k 0 ∂ k, ε rk, ε ak, ε (t) * e = . r k,ε D −(t)e+iω k L t +/ e k, ε rk, ε ak, ε (t) ∂t ~ , k,ε which can be solved in terms of time-ordered exponentials of system operators (− → )  t    Ek → µ 0 · −ε † † → − p a (t) p a (0) * D+(t 1)e−iω k L t1+/ , e k, ε k, ε k, ε = e k, ε k, ε k, ε T exp . dt 1 pk,ε ~ 0 k,ε ) , (→ (22)  − → −ε t  E µ ·   k ← − 0 * r k, ε a k, ε (0) − +iω k L t 1+ r k, ε a k, ε (t) k, ε k, ε /e . e D (t 1)e = T exp . dt 1 r k,ε ~ 0 k,ε , ← − The operator T serves to rearrange the operators appearing to its right in standard time order, with later times appearing to the → − left of earlier times. The operator T rearranges the operators in the reversed time order, with later times appearing to the right of earlier times. Assuming that the composite system/field density matrix factorizes at t = 0 and that the quantum field contains zero photons at this time (i.e., ρ(0) = σsys(0) |0⟩⟨0| with σsys(0) the reduced system density matrix), the disentangled expressions of Eq. (22) lead immediately to (− → )  t      Ek → µ 0 · −ε  → − * + −iω k L t 1+   T exp .  / dt p D (t )e 1 k,ε 1      k,ε ~    † 0 p a (t) r a (t) , - . ( ) =  e k, ε k, ε k, ε e k, ε k, ε k, ε (23)  → − → −  t    Ek µ 0 · ε   ← − * + × T exp . dt 1 r k,ε D −(t 1)e+iω k L t1/    k,ε ~   0 , -  Substituting this expression in Eq. (21) yields our final expression for the generating function (− →   t t − )2 2 → *. Ek µ 0 · ε +/ + − −iω k L (u−v) (s (u) (v) G (⃗s,t) = TN exp . − 1) D D e dudv (24) kε / . ~2 k,ε 0 0 , -

⟨

⟩


174104-5


J. Chem. Phys. 142, 174104 (2015)

The operator TN acts on all operators to the right of it by first arranging all the “+” operators to the left of all the “−” operators and subsequently placing all “−” operators in standard time order (latest times at the left) and all “+” operators in reversed time order (latest times at the right). Equation (24) is fully general and opens up the possibility to calculate a seemingly limitless array of statistics related to photon emission, once the dynamics for the radiating system are established. In the interest of both concreteness and simplicity, we focus further attention on statistical quantities related to the number of photons emitted within a given frequency window, irrespective of propagation direction and polarization state. The number operator for photons emitted within such a window is defined as  †   ak,ε ak,ε , (25) N(ω,∆) ≡ ∆ ∆ k,ε:(ω− 2 ) ≤ω k ≤ (ω+ 2 ) where the indicated sum over k, ε is understood to include all photon modes consistent with the indicated frequency window centered at ω with width ∆. Our previously introduced generating functions (Eqs. (2), (5), and (24)) are readily converted to this frequency resolved picture by consolidating all photon modes associated with (ω, ∆) together under the common auxiliary variable sω . This leads to   (t) N G (⃗s,t) = sω (ω,∆) ω   * + = N exp N(ω,∆)(t) (sω − 1) ,ω   ω+∆/2 * Γ0  (sω − 1) = TN exp .. dω ′ 2π ω ω−∆/2 ,  t t +/ −i[ω ′−ω L ](u−v) − + dudv / . D (u) D (v) e × 0

0

-

(26)

In the final expression, we have carried out the sum over wave vector directions and polarizations by moving to the continuum limit of V0 → ∞ where the sums become integrals.68 Here, ω3 | µ ⃗ |2

eg 0 Γ0 ≡ 3πε 3 and ω eg is the frequency corresponding to the 0~c energy gap between the excited and ground electronic states of the molecule, i.e., ωeg = ω m e − ω n g where it is assumed that the electronic energy scale overwhelms the nuclear energetics and ω m e − ω n g is effectively constant for all values of ng and me . This result approximates a factor of ω ′3 within the frequency

3 integral (sum) by the constant factor ωeg , which is analogous to the same approximation introduced in the Wigner-Weisskopf treatment of spontaneous emission.68 The constant Γ0 is, in fact, the spontaneous emission rate for a two-level absorber with frequency splitting ωeg . It is critical to recognize that the sum (product) over frequencies appearing in Eq. (26), with windows defined by the summation in Eq. (25), assumes that the chosen frequency windows are completely non-overlapping. Overlapping windows would invalidate the second and third lines of Eq. (26) as these results rely upon the operator identity introduced in Eq. (5), which assumes that the creation/annihilation operators associated with different auxiliary variables commute with one another. We will also assume for convenience, though it is not essential, that every frequency window considered −s in shares the same width, ∆. We also stress that the vector → Eq. (26) now represents the set of frequency associated variables {sω }, in contrast to the usage introduced previously. This should not lead to any confusion below. The case of broadband statistics, where all emitted photons are consolidated under a single auxiliary variable s, is realized by taking the limit ∆ → ∞ in Eq. (26). The resulting frequency integral reduces to a delta function in this limit, leaving     t ′ + ′ − ′ G(s,t) = TN exp (s − 1)Γ0 dt D (t )D (t ) . (27) 0

This expression is equivalent to other formulations of broadband photon emission statistics in resonance fluorescence.34,72 In particular, broadband emission statistics calculated with Eq. (27) reproduce the predictions of simplified “generating function” calculations which count all photons indiscriminately via tracking instantaneous spontaneous emission events from the chromophore.33,34,37,65 However, it is worth emphasizing that the simplified approach used in the calculation of broadband statistics cannot be extended to the general case. Decay of an electronic excitation into a particular field mode or narrow band of modes is a non-Markovian process and there is no simple protocol to relate spontaneous emission events to photon emission with frequency resolution. C. Expressions for the moments in terms of system operators

Performing differentiations on Eq. (26) with respect to the sω variables results in factorial moments analogous to the mode resolved expressions previously discussed. For example,

n



 ( Γ0 ) n+m (n) (m) (t) N(ω N(ω,∆) ′,∆) (t) = 2π

⟨

ω+ ∆   t t  *.  2 +    TN  .. dω1 D + (u) D − (v) e−iω1L (u−v)dudv ///    ∆ 0 0  , ω− 2 m , ∆ ′ ω +2  t t    *.  +/    + − −iω 2L (u−v) × .. dω2 D (u) D (v) e dudv //      ∆ 0 0 ′ ,ω − 2 - 

⟩

(28)

which is general enough to cover all expressions we explicitly compute below. This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Fri, 08 May 2015 07:26:14

174104-6


J. Chem. Phys. 142, 174104 (2015)

In what follows, we will focus on a few quantities of special interest. The average number of photons radiated into a given frequency window is provided by ω+ ∆2



Γ0 N(ω,∆) (t) = 2π

t t

dω1 0

ω− ∆2

D + (u) D − (v) × e−iω1L (u−v)dudv.

(29)

0

Another quantity of interest is the second multivariate moment, i.e., C2(ω,ω′,∆) (t) ≡ N(ω,∆) (t) N(ω′,∆) (t) − δω,ω′

Γ2 N(ω,∆) (t) = 02 4π

ω+ ∆2

t t t t



dω1 ω− ∆2

+

+

ω ′+ ∆2



dω2

ω ′− ∆2

0

0

0

0

× TN D (u) D (v) D (w) D ( y) eiω1L (u−w)eiω2L (v−y)dudv dwd y.

−

−

(30)

Carrying out both the normal and time ordering, the expression becomes Γ2 C2(ω,ω′,∆) (t) = 02 4π

ω ′+ ∆2

ω+ ∆2





dω1 ω− ∆2

ω ′− ∆2

t t t t

dω2 0

0

0

D + (u) D + (v)

0

iω 1L (u−w) iω 2L (v−y)

e e *. iω1L (u−y) iω2L (v−w)+/ +e e × D − (w) D − ( y) ... iω1L (v−y) iω2L (u−w)/// Θ (v − u) Θ (w − y) dudv dwdy, +e e / . iω (v−w) iω 2L (u−y) e ,+e 1L where Θ (t) is the Heaviside theta function. From Eqs. (29) and (30), we introduce a normalized correlation function as

N(ω,∆) (t) N(ω′,∆) (t) − N(ω,∆) (t) N(ω′,∆) (t) ′ C∆ (ω,ω ,t) = − δω,ω′. 

N(ω,∆) (t) N(ω′,∆) (t)

(31)

(32)

The chosen normalization ensures that, in the specific case of ω = ω ′, the correlation function reduces to Mandel’s Q parameter for the number of photons emitted within a given frequency range, defined as

2 N(ω,∆) (t) N(ω,∆) (t) − N(ω,∆) (t)

− 1. (33) Q∆ (ω,t) = C∆ (ω,ω,t) = N(ω,∆) (t) IV. CALCULATION OF THE MOMENTS FOR THE MODEL SYSTEM OF SEC. III A. Evaluating correlation functions of chromophore operators

To make practical use of the expressions appearing in Sec. III C, one must specify how the multipoint correlation functions, central to these formulae, are to be evaluated. To this end, we introduce a Markov approximation for evolution of chromophore dynamics; the relevant correlation functions involve only D ± operators acting in the CH space of the total system Hilbert space. The Markov approximation involves the traditional set of assumptions routinely introduced in the derivation of the optical Bloch equations/optical master equation and Redfield theory; namely, (1) a weak interaction between the chromophore and both the quantum radiation field and I and V  appearing in Eqs. (7) system bath (i.e., that the H and (8) are “small”); (2) much faster correlation times within the bath and quantum field than the time-scales associated with chromophore dynamics; and (3) negligible influence of the chromophore upon the bath and quantum field such that the density matrices associated with bath and quantum field

remain thermalized over all times for the purposes of calculating the dynamics of the chromophore. It is well known that the above assumptions will lead to an elementary dynamics for the reduced chromophore density matrix in the rotating frame σC H (t) ≡ Traceb,R{ ρ(t)}, σ ˙ C H (t) = LσC H (t), σC H (t) = eL(t−t1)σC H (t 1).

(34)

Explicit expressions for the matrix representation of the timeindependent rotating frame reduced Liouville super-operator L are presented in Sec. IV A. For present purposes, we simply note that σC H (t) contains (Nmaxg + Nmaxe )2 elements; the matrix representation of L has dimensions (Nmaxg + Nmaxe )2 × (Nmaxg + Nmaxe )2. The TN operator ensures that all of the correlation functions related to photon counting observables are of the form

D + (u1) D + (u2) . . . D + (um ) D − (vm ) D − (vm−1) . . . D − (v1) − = Trace D (vm ) D − (vm−1) . . . D − (v1) ρ(0) × D + (u1) D + (u2) . . . D + (um ) (35)


174104-7


J. Chem. Phys. 142, 174104 (2015)

with um ≥ um−1 ≥ · · · ≥ u1 and vm ≥ vm−1 ≥ · · · ≥ v1. Multipoint correlation functions with the time-ordering characteristics displayed in Eq. (35) can be calculated through a straightforward generalization of the considerations leading to Eq. (34).73 The result requires the introduction of some new notation. First, one has to order all the times ui , vi into a single set of times τ1 . . . τ2m such that τ2m ≥ τ2m−1 . . . ≥ τ1. Then, the super-operator D(i) is defined by its action upon an arbitrary quantum mechanical operator A depending upon whether time τi was originally associated with a D +(u j ) or a D −(vk ) operator D(i) A = D− A = D − A D(i) A = D+ A = AD +

if if

τi ⇐⇒ vk , τi ⇐⇒ u j .

Using this formalism, Eq. (37) may be written as

+ D (u1) D+ (u2) . . . D + (um ) D − (vm ) D − (vm−1) . . . D − (v1) = (I|D(2m)eL(τ2m −τ2m−1)D(2m−1)eL(τ2m−1−τ2m−2) . . . D(2)

(36)

(The operators D − and D + appearing above are Schrödinger picture operators (in the rotating frame) and carry no time dependence.) Then,73

+ D (u1) D + (u2) . . . D + (um ) D − (vm ) D − (vm−1) . . . D − (v1)  = TraceC H D(2m)eL(τ2m −τ2m−1)D(2m−1)  × eL(τ2m−1−τ2m−2) . . . D(2)eL(τ2−τ1)D(1)eLτ1σC H (0) . (37) For explicit calculations it is most convenient to reformulate Eq. (37) within the Liouville space formalism.74 Recall that Liouville space is the vector space composed of all linear operators acting upon the state vectors of a quantum mechanical Hilbert space; this includes both Hermitian operators and non-Hermitian operators. Adopting the bra-ket-like notation of Blum,74 a standard quantum mechanical operator in the Schrödinger representation A is associated with an element | A) in Liouville space. The inner product of Liouville space is defined as (A|B) = Trace A† B , (38) and the super-operators introduced above, which act on quantum operators to yield new quantum operators, are just simple operators within Liouville space acting on Liouville space elements to yield new Liouville space elements. For example,

we have various equivalent ways to express the non-equilibrium expectation value of the operator A, ⟨A(t)⟩ = Trace { Aρ(t)} = Trace AeLt ρ(0)   † = Trace [eLt ]† A† ρ(0) = A†| ρ(t) = A†|eLt ρ(0) = A† eLt | ρ(0)) = [eLt ]† A†| ρ(0) . (39)

× eL(τ2−τ1)D(1)eLτ1|σC H (0)),

(40)

where I is the identity operator. We now consider Liouville state elements that are eigenfunctions of L. These normalized elements and their associated eigenvalues are defined by L |α) = λ α |α) , (α| L = (α| λ α ,

(41)

(α| β) = δα, β , and obey the completeness relation  |α)(α| = I

(42)

α

with identity super-operator I. It is important to recognize that L is not Hermitian and that the left and right eigenvectors are not expected to be simply related to one another. However, each distinct eigenvalue is associated with a pair of left and right eigenvectors and these eigenvectors obey the orthonormality condition indicated above. The zero eigenvalue λ 0 = 0 is of particular importance; it is associated with the steady state chromophore density matrix as a right eigenvector |0) = | ρ ss ) and the identity operator as a left eigenvector (0| = (I|. The remaining eigenvalues are expected to have negative real parts indicating a decay to the steady state. By repeated use of Eqs. (42) and (40) becomes

D + (u1) D + (u2) . . . D + (um ) D − (vm ) D − (vm−1) . . . D − (v1)

(2m) ) (2m−1) (τ −τ ) λ α (τ −τ λα  *.Dα1,α2m e 2m 2m 2m−1 Dα2m,α2m−1e 2m−1 2m−1 2m · · ·+/ // (α1|σC H (0)). .. = / . λ α 2(τ 2−τ 1) (1) λ α 1τ 1 α 1 ···α 2m × D(2) e D e α 3,α 2 α 2,α 1 , (2m) λ α (τ −τ ) (2m−1) (τ λα −τ )  *.D0,α2m e 2m 2m 2m−1 Dα2m,α2m−1e 2m−1 2m−1 2m · · ·+/ // .. = / . λ α 2(τ 2−τ 1) (1) α 2 ···α 2m D (if σ (0) = σ ). × D(2) e CH ss α 3,α 2 α 2,0 ,

(43)

Dα, β ≡ (α|D| β).

As indicated, the third line holds true if the system begins in steady state at t = 0, which will be assumed for the cases studied in this work. Equation (43) presents an explicit expression for the multipoint chromophore correlation functions pre-

sented in this work. The primary burdens associated with evaluating this expression involve diagonalization of L, transforming the D(i)’s to the basis of L eigenfunctions, and carefully accounting for the interleaving of u, v times to ensure that the


174104-8


J. Chem. Phys. 142, 174104 (2015)

individual D(i)’s are associated with the proper D ± character depending upon the set of time arguments. Before proceeding to specific cases in Subsections IV B and IV C, we briefly comment on an apparent paradox associated with the Markov approximation invoked in this section. A key assumption of the Markov approximation is that the quantum field remains in thermal equilibrium at all times, for the purposes of calculating chromophore dynamics (see above). For the optical frequencies of interest in this work and typical laboratory temperatures, this corresponds to a field effectively devoid of photons. However, the entire point of expression (43) is to use it to count photons emitted into the quantum radiation field. The resolution to this apparent paradox is to recognize that, although the chromophore emits photons which are measurable, these photons are few in number and are not expected to significantly influence the dynamics of the chromophore system. Photon emission from the chromophore is driven entirely by the applied laser field. Stated differently, the Markov approximation assumes that photons emitted from the chromophore have no effect on future chromophore dynamics and future emissions, although the emission of these photons into the field is explicitly tracked and quantified by Eq. (26) and derived quantities. If we assume that all the emitted photons are immediately detected, and thereby removed from the system, this assumption is quite reasonable. B. Calculation of the first moment

Equation (29) can be written in a manifestly time-ordered form ω+ ∆2

Γ0 2 Re N(ω,∆) (t) = 2π



t t

dω1 0

ω− ∆2

D + (u) D − (v)

0

× e−iω1L (u−v)Θ (u − v) dudv.

(44)

Applying Eq. (43) to the correlation function in the expression above enables us to write it as 

+ D (u) D − (v) Θ (u − v) = D+0,α D−α,0e λ α (u−v)Θ (u − v) , α

(45) where we have assumed σC H (0) = σ ss . It then follows that   Γ0 N(ω,∆) (t) = 2 Re  D+0,α D−α,0nα (t,ω, ∆) , 2π  α  where

ω+ ∆2



nα (t,ω, ∆) ≡

t t e(λ α −iω1L )(u−v)Θ (u − v) dudv

dω1 0

ω− ∆2 ω+ ∆2



= ω− ∆2

(46)

dω1

0

e t(λ α −iω1L ) − 1 − t(λ α − iω1L ) . (λ α − iω1L )2

(47)

The first moment of the number of emitted photons is proportional to t in the long time limit and the time derivative

of average emitted photon number saturates to constant at long times. We thus define the intensity of emitted photons as   dN(ω,∆) (t) , (48) I(ω,∆) (t) ≡ dt which becomes I(ω,∆) (t) =

  Γ0 2 Re  D+0,α D−α,0 x α (t,ω, ∆) , 2π  α 

(49)

with x α (t,ω, ∆) given by ω+ ∆2



x α (t,ω, ∆) =

t e(λ α −iω1L )(t−v)dv

dω1 0

ω− ∆2 ω+ ∆2



=

dω1

1 − e(λ α −iω1L )t . iω1L − λ α

(50)

ω− ∆2

The long time limit is conventionally decomposed into a coherent intensity associated with elastic scattering and an incoherent part via I(ω,∆) = I(ω,∆)coh + I(ω,∆)incoh with

I(ω,∆)coh = Γ0 D + D − δω,ω L = Γ0D+0,0D−0,0δω,ω L (51) and   ω+ ∆2     ′ Γ0 1 , 2 Re  dω1 I(ω,∆)incoh = D+0,α D−α,0 2π iω1L − λ α   α  ω− ∆2  (52) where the prime denotes the fact that we omit the 0 eigenvalue from the summation. It is to be understood that the Kronecker delta function appearing in Eq. (51) assumes the value of one when ω L lies within the ∆-wide frequency window centered at ω and is otherwise zero. The formalism detailed in this section allows the calculation of the emission spectrum for any detuning frequency (within the validity regime of the RWA) and for any system of the form described in Sec. III (as evaluated within the Markov approximation detailed in Sec. IV A). It is important to note that the spectrum defined by Eqs. (48)-(50) is the quantum equivalent to the Page75 and Lampard76 spectra. Thus, for short times, it can yield apparently non-physical results.77 In particular, I(ω,∆) (t) can assume negative values at a given time, which would correspond to the disappearance of previously emitted photons from the field. While such a process (reabsorption of photons previously emitted by the chromophore) is perfectly allowable within quantum mechanics, the possibility of negative “intensities” raises serious questions about how one intends to actually detect photons experimentally— certainly photons that have been registered on a measuring apparatus cannot subsequently be unmeasured. The question of the experimentally measured spectrum and higher order statistics will be discussed in Sec. V.


174104-9


J. Chem. Phys. 142, 174104 (2015)

C. Calculation of the second multivariate moment

The second multivariate moment is described in Eqs. (30) and (31). The nature of the multivariate second moment requires that it be symmetric under ω ω ′ exchange. In order to use this symmetry, we write it as C2(ω,ω′,∆) (t) =

Γ02 4π 2

ω+ ∆2

ω ′+ ∆2





dω1 ω− ∆2

dω2

× [ f (t,ω1,ω2) + f (t,ω2,ω1)] ,

(53)

where t t t t dudv dwd y 0 +

0

0

0

× D (u) D + (v) D− (w) D − ( y) eiω1L (u−w)eiω2L (v−y) × * iω1L (u−y) iω2L (v−w)+ e ,+ e × Θ (v − u) Θ (w − y) .

(54)

In order to calculate function of Eq. (32), we the correlation express N(ω,∆) (t) N(ω′,∆) (t) in a way that will enable us to subtract it from the second moment. Using Eq. (29), we write it as

Γ02 4π 2

N(ω,∆) (t) N(ω′,∆) (t) =

ω+ ∆2

ω ′+ ∆2





dω1 ω− ∆2

dω2

ω ′− ∆2

× [h (t,ω1,ω2) + h (t,ω2,ω1)] , (55) where 1 h (t,ω1,ω2) = 2

t t t t

0

0

0

D + (u) D − (v)

0

+

× D (w) D − ( y) e−iω1L (u−v)

× e−iω2L (w−y)dudv dwdy.

ω ′+ ∆2

ω+ ∆2





dω1 ω− ∆2

dω2

ω ′− ∆2

× [M (t,ω1,ω2) + M (t,ω2,ω1)] ,

In the standard theory of photon detection,68,69,78 one typically focuses interest on correlation functions of the quantized electric field. These quantities are closely related to, but different from, those introduced in Sec. IV. To establish the relationship between these pictures, we first review the case of broadband photon detection. In contrast to Sec. IV, we imagine counting photons that are actually measured by a detector, as opposed to counting photons that have been emitted from the irradiated system. Photodetection typically relies on technology that directly measures the intensity of the measured light incident upon the detector, e.g., as in a photomultiplier. We imagine an idealized photodetector with the following properties. (1) The detector completely covers a spherical shell surrounding the system, which is located at the center of the shell. (2) The distance from the system to the detector, r, is in the far field of the emitting system. (3) The detector displays no frequency or polarization dependence and has a quantum efficiency of one; detector response is instantaneous and there is no dead time. (4) The presence of the detector does not affect the planewave modes of the cavity used to quantize the field; i.e., the photons are unaffected by the detector until the time at which they impinge upon it and are absorbed. The Heisenberg electric field operators for the quantum radiation field are written62,68,69 → − → − → − E (r,t) = E +(r,t) + E −(r,t),  → −+ → −ε E  ik·r E (r,t) = i , (58) k ak,ε (t)e k,ε

(56)

Subtracting the expression of Eq. (55) from the expression of Eq. (53), we can formally write the numerator of the correlation function C∆ (ω,ω ′,t) (Eq. (32)) as

C2(ω,ω′,∆) (t) − N(ω,∆) (t) N(ω ′,∆) (t) Γ2 = 02 4π

V. CALCULATION OF THE MOMENTS WITH EXPLICIT CONSIDERATION DETECTOR EFFECTS A. Correspondence between the results of Sec. IV and the standard theory of broadband photon detection

ω ′− ∆2

f (t,ω1,ω2) =

(Eq. (32)) as a combination of terms that are explicitly evaluated in Appendix A. Numerical values of the correlation function (Eq. (57)) are obtained via numerical integration over the frequencies ω1 and ω2.

(57)

with M (t,ω1,ω2) = f (t,ω1,ω2) − h (t,ω1,ω2). As written, this expression is of little practical use. To obtain numerical values for the correlation function, M (t,ω1,ω2) must be expressed in terms of elementary pieces that are simple to calculate. This tedious procedure is detailed in Appendix A. The function M (t,ω1,ω2), together with the expression for the first moment in Eq. (46), yields the correlation function C∆ (ω,ω ′,t)

 → −− † → −ε E  −ik·r E (r,t) = −i . k ak,ε (t)e k,ε

The creation and annihilation operators appearing above can be written as time integrals over the system operators as shown in Eqs. (17) and (20). Making this substitution and performing the sums over k, ε, by taking the V0 → ∞ limit and evaluating as an integral, leads to ( c rr ) −  → −+ µ0 e−iω L t E (r,t) = 1 − 2 · → 2 r 8π ε 0r  ∞ × dk k 2 eik r − e−ik r 0



t

× 0

=

→ − dt 1 D −(t 1)eiω k L (t1−t) + E +f (r,t)

 ω2 rr ) → eg −µ · e−iω L (t−r /c) D −(t − r/c) 0 2 r 4πε 0c2r → − → − → − + E +f (r,t) ≡ E +S (r,t) + E +f (r,t), (59) (

1−


174104-10


J. Chem. Phys. 142, 174104 (2015)

where the second equality is obtained by extending the lower limit of the k integral to −∞ and approximating k ∼ ωeg /c as discussed under Eq. (26). The remaining k (equivalently, ω) integral reduces to a Dirac delta function, which explains the → − remarkable simplification seen in this step. E +f (r,t) represents the “free” contributions to the field that would be present even in the absence of the radiating system (i.e., the contribution → − from the a(0) pieces in Eq. (20)) and E +S (r,t) the contribution due to radiation from the system. The explicit form provided → − for E +S (r,t) presumes that the detector is in the far field; contributions of order 1/r 2 and smaller have been intentionally dis→ − carded. The similar calculation for E −(r,t) yields 2 ( rr ) −  ωeg → −− µ0 eiω L (t−r /c) E (r,t) = 1 − 2 · → r 4πε 0c2r → − → − → − × D +(t − r/c) + E −f (r,t) ≡ E −S (r,t) + E −f (r,t). (60) The intensity operator for light incident on a given point of our far-field detector is69 → − → − → − I(⃗r ,t) = 2ε 0c E −S (r,t) · E +S (r,t) + terms including E ±f . (61) Dividing this quantity by the energy per photon, ~ωeg provides an operator corresponding to the local photon flux into the detector at position r. We integrate this quantity over the detector surface to provide the total rate of photon detection operator (all photons hitting the detector are registered due to the assumed perfect efficiency) and further integrate over time to yield an operator corresponding to the total number of detected photons in the interval [0,t],  d d N(t) 2ε 0cr 2 → − → − = dΩ E −(r,t) · E +(r,t) dt ~ωeg → − + terms including E ±f , (62)   2ε 0cr 2 t ′ → − → − d N(t) = dt dΩ E −(r,t ′) · E +(r,t ′) ~ωeg 0 → − + terms including E ±f

(63)

with the Ω integral extending over the full 4π solid angle of the detector. Carrying out the angular integration explicitly yields  t d N(t) = Γ0 dt ′D +(t ′ − r/c)D −(t ′ − r/c) 0

→ − + terms including E ±f .

into the expression for the generating function for detected  photons d G(s,t) = pn (t)s n yields d

where the T: ordering operator arranges all E − fields (irrespective of S or f character) to the left of the E + fields, and timeorders the (−) fields with latest times to the right, and timeorders the (+) fields with latest times to the left. Inserting this

(66)

Although not obvious on first inspection, it can be → − shown71,80 that the terms involving E ±f make no contribution to Eq. (66), which is why we have not been careful in explicitly tracking these contributions above. This is due to the particular time ordering associated with T: and the fact that causality dictates that the free field at time τ must commute with all system operators at times t < τ. Any term involving a free field contribution can be written such that the free field contribution(s) appear at the leftmost and rightmost edges of the expression to be averaged. The averaging operation then renders these terms zero due to the assumed absence of photons in the field at t = 0. The only terms that remain in the generating function are those comprised solely of system operators and we may write d

G(s,t)     t ′ + ′ − ′ = TN exp (s − 1)Γ0 dt D (t −r/c)D (t − r/c) . 0

(67) This expression is nearly identical to Eq. (27), which corresponds to the monitoring of all emitted photons, irrespective of frequency, propagation direction, and polarization. The only difference appearing between the two expressions is the translation in time of the system operators in the current expression by −r/c, reflecting the finite time it takes for the emitted photons to reach the detector. With this single (trivial) exception, the two methods for calculating photon statistics are identical as applied to broadband monitoring/detection schemes. We emphasize that the correspondence between Eqs. (26) and (67) follows from our completely idealized treatment of the detector, which ensures the detection of all emitted photons. Certain aspects of a non-ideal detector cause only minor changes to this expression. If the detector is assumed to occupy only some fraction of the spherical surface surrounding the emitting system, this affects only the constant pre-factor multiplying the time integral in Eq. (67) due to an incomplete integration over Ω in Eq. (63). Similarly, a detector with quantum efficiency η < 1 requires further scaling of d N(t) by this factor. A slightly more general version of Eq. (67) can therefore be written

(64)

The probability of actually observing n photons at the detector over the time interval [0,t] is a classic problem in quantum optics,72,78,79 which we shall not re-derive here. Using the present notation, the result is   d d [ N(t)]n e− N (t) p(n,t) = T: , (65) n!

G(s,t) = T: exp (s − 1)d N(t) .

 d





G(s,t) = TN exp (s − 1)α d Γ0

t ′

+ ′

 − ′

dt D (t )D (t )

,

(68)

0

where the constant α d < 1 accounts for the inability of the detector to capture all emitted photons. We have also translated the origin of time by r/c in this expression relative to Eq. (67) in order to make the correspondence with Eq. (27) manifest. B. Detection of frequency filtered photons

In Sec. IV B, we briefly commented that frequencyresolved statistics based on photon emission, as opposed to


174104-11


J. Chem. Phys. 142, 174104 (2015)

photon detection, are problematic at short times and can give rise to unphysical results. Some deficiencies of the emission approach will be explicitly demonstrated in the calculations of Sec VI. To remedy these shortcomings, it is necessary to consider the detection process explicitly in the theoretical description. Unlike the case of broadband statistics, there is no clear way to establish correspondence between experimental photon counting measurements at finite times and the theoretical approach of Sec. II; the approach of Sec. V A must be extended to account for the physical process of spectral filtering as realized experimentally. The procedure we will adopt to model the effect of spectral filtering on photon detection was first discussed by Eberly and Wodkiewicz77 in the context of classical fields. They suggested that considering the response function of the physical device used to filter the photons prior to detection is necessary to compute the spectrum at finite times and that such an approach removes the unphysical artifacts associated with schemes related to that presented in Sec. IV B. Subsequent work63,81 has confirmed that the procedure suggested by Eberly and Wodkiewicz77 can be rigorously justified within the context of a quantized radiation field. Both the classical approach and practical implementations of the quantum approach assume that the presence of spectral filters and the detectors themselves do not affect the dynamics of the emitting system-i.e., any radiation scattered off the filters and detectors is ignored for the purpose of calculating the evolution of the system, which is assumed to evolve (as discussed above) in a quantum radiation field driven by a semiclassical CW excitation at frequency ω L . We assume a filtering/detection scheme as diagrammed in Figure 1. Fluoresced light emitted from the system is intercepted by a filter immediately prior to impinging upon the detector. The filter is introduced because of its spectral filtering properties on the light incident upon it. For our purposes, this filtering is most conveniently expressed by the effect on the electric field

in the time domain77  → −+ E D (r,t) = → −− E D (r,t) =

t

→ − dt ′H(t − t ′) E +(r,t ′),

−∞ t



′

∗

−− ′→

(69) ′

dt H (t − t ) E (r,t ). −∞ → −± Here, E D (r,t) represents the operators associated with the electric field that reaches the detector after passage through the → − filter, as compared to the bare fields E ±(r,t), from Sec. V, that would be present in absence of the filter. The response function H(t) reflects the influence of the filtering device. For practical purposes, we assume a functional form for H(t) motivated by the action of a lossless and highly reflective Fabry-Perot interferometer with central frequency ω1 and bandwidth ∆, H(t) = ∆e−(∆+iω1)t .

(70)

This yields detectable fields of 2 ( rr ) −  ∆ωeg → −+ µ0 E D (r,t) = 1 − 2 · → r 4πε 0c2r  t ′ ′ → − × dt ′e−(∆+iω1)(t−t )e−iω L t D −(t ′) + E +f (r,t), −∞

→ −− E D (r,t) =

( 1− 

t

× −∞

 ∆ω2 rr ) → eg −µ · 0 r2 4πε 0c2r

(71)

′ ′ → − dt ′e−(∆−iω1)(t−t )e+iω L t D +(t ′) + E −f (r,t),

where we have invoked the same translation in time by r/c here as was introduced in Eq. (68). Though this filter response can be generalized in principle, the computational scheme we adopt below largely depends on the exponential form of H(t). → − → − Substituting these results for E ±D (r,t) in place of E ±(r,t) in the equations of Sec. V A (i.e., Eqs. (63) and (66)) leads to an expression for the generating function of spectrally filtered photons

  t t1 t1 *. 2 + G (s, ∆,t) = TN exp .∆ α d Γ0 (s − 1) D + (u) D− (v) e−iω1L (u−v)e∆(u+v−2t1)dudv dt 1// , 0 −∞ −∞ , 

d

(72)

where we have again introduced the α d factor to account for the fraction of solid angle covered by the detector and the detector quantum yield. Considering multiple detectors arranged around the sample, each with a different central frequency, but identical bandwidths and α d factors, this result generalizes to   t t1 t1 + *. 2 G (⃗s, ∆,t) = TN exp . ∆ α d Γ0 sω1 − 1 D + (u) D − (v) e−iω1L (u−v)e∆(u+v−2t1)dudv dt 1// , ω 0 −∞ −∞ , 1 

d

where each sω1 is the auxiliary variable associated with the counts at a particular detector. In what follows, we will always make the assumption that α d = 1,

(74)

in order to make direct comparison between the results of two theoretical schemes we have discussed (based on emitted photons and detected photons, respectively). However,

(73)

it should be stressed that this situation is not physically realizable, even in an idealized gedanken experiment. It is not possible to place multiple α d = 1 filtering detectors around the sample since each of these perfect detectors requires a full 4π solid angle around emitter and the action of one filter would necessarily interfere with other filters placed behind it. Notwithstanding this point, the α d factors are simple constants that contribute to the theoretical predictions


174104-12


J. Chem. Phys. 142, 174104 (2015)

for the photon counting moments in a transparent way. In Subsections V C and V D, we will focus on the first two moments of detected photon number. These expressions

are readily generalized to account for arbitrary values of α d simply by replacing Γ0 → Γ0α d wherever factors of Γ0 appear.

C. Calculation of the first moment

The first moment is given by t1

t ⟨ N (ω1, ∆,t)⟩ = ∆ Γ0 d

2

dt 1e

−2∆t 1

du

−∞

0

t1

dv⟨D + (u) D − (v)⟩e−iω1L (u−v)e∆(u+v)

−∞

  t1 u t  (λ α −iω 1L )(u−v) ∆(u+v) −2∆t 1 2 + −   du dve = 2∆ Γ0 Re  D0,α Dα,0 dt 1e e  α  −∞ −∞ 0   1 . = t∆2Γ0 Re  D+0,α D−α,0 (75) 2∆(∆ − λ α + iω1L )   α  The second line follows from similar considerations as introduced in Sec. IV B and the third has simply evaluated the integrals. The intensity follows immediately and is stationary in time d

  D+0,α D−α,0 . I (ω1, ∆) = 2∆2Γ0 Re   α 2∆(∆ − λ α + iω1L ) 

(76)

D. Calculation of the second multivariate moment

The correlation function is defined as

d d

C∆ (ω1,ω2,t) ≡

N(ω1,∆) (t) d N(ω2,∆) (t) − d N(ω1,∆) (t) d N(ω2,∆) (t) − δω1,ω2, 

d dN N(ω2,∆) (t) (ω 1,∆) (t)

(77)

with the generalized Q parameter following as  d

d

Q∆ (ω,t) ≡ dC∆ (ω,ω,t) =



2 2 (t) − d N(ω,∆) (t) N(ω,∆)

d − 1. N(ω,∆) (t)

(78)

Using similar motivation to the approach outlined in Sec. IV C, the second multivariate moment follows from Eq. (73), as t d

C2(ω1,ω2,∆) (t) ≡

d

N(ω1,∆) (t) N(ω2,∆) (t) − δω1,ω2 d

d

N(ω1,∆) (t) =

Γ02∆4

e 0

t2 ×

t1

t1 dw

−∞

t2 e

dt 2 0

−2∆t 1

dt 1

du

−∞

d y TN D + (u) D + (v) D − (w) D − ( y) eiω1L (u−w)eiω2L (v−y)e∆(u+v+w+y).

dv

−∞

t −2∆t 2

(79)

−∞

To calculate the numerator of the correlation function, the product of first moments, t ⟨ N (ω1, ∆,t)⟩ ⟨ N (ω2, ∆,t)⟩ = d

d

t

∆4Γ04 0

0

t1 ×

t1 du

−∞

dt 2e−2∆(t1+t2)

dt 1

−∞

t2 dv −∞

t2 dw

dy⟨D + (u) D − (v)⟩ ⟨D + (w) D − ( y)⟩e−iω1L (u−v)−iω2L (w−y)e∆(u+v+w+y),

−∞

(80) must be subtracted. In Appendix B, we show how the correlation functions of Eqs. (79) and (80) can be expressed as sums of simple terms. The explicit forms of the individual terms are provided in Appendix B, allowing calculation of the detected correlation functions. This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Fri, 08 May 2015 07:26:14

174104-13


J. Chem. Phys. 142, 174104 (2015)

FIG. 1. Schematic for experimental photon detection. Frequency dependent filters are placed just in front of broadband photon detectors. Photons emitted from the source are filtered prior to detection and the apparatus counts only those photons that survive the filtering process. Multiple filter-detector pairs can be oriented around the sample to capture correlations between photons of different frequency.

VI. RESULTS FOR A TWO-LEVEL CHROMOPHORE

The two-level system serves as the prototypical example of an irradiating chromophore in quantum optics and, indeed, most theoretical studies to date concentrate focus solely on this simplest possible emitting system. Adopting the notation of Sec. III, the system and chromophore of the two-level-system are the same; there is no bath coupled to the optically active chromophore. That is, the terms Vˆ and Hˆ b are absent from Eq. (8) and 1 Hˆ sys = Hˆ CH = ~ω0 (|e⟩ ⟨e| − |g⟩ ⟨g|) . 2

(81)

The dipole moment operator is given simply by −µ D+ + D − ; D ≡→ 0

D + = |e⟩ ⟨g| ;

D − = |g⟩ ⟨e| .

(82)

The reduced chromophore density matrix σC H (t) is 2 × 2 in dimension and corresponds exactly to the density matrix associated with the optical Bloch equations. The Liouville superoperator L from Sec. IV A is well known for this case with the matrix representation −Γ0 *. .. . Γ0 L = ... ..−i Ω0 .. 2 Ω0 ,+i 2

0 0 Ω0 2 Ω0 −i 2 +i

Ω0 2 Ω0 −i 2

Ω0 2 Ω0 +i 2

+i

−i

−iδ L − 0

Γ0 2

+/ // // // , // 0 / Γ0 / iδ L − 2

(83)

using our previously introduced notations for the spontaneous emission rate, Γ0 and Rabi frequency, Ω0. δ L = ω L − ω0 is the detuning frequency of the applied field. This matrix is assumed to act on a vector representation of the density matrix ordered as σC H = (σee , σg g , σg e , σeg )†. In this same basis, the superoperators associated with the dipole operator are

0 *. 0 D+ = ... .0 ,1

0 0 0 0

0 1 0 0

0 + 0// /, 0// 0-

0 *. 0 D− = ... .1 ,0

0 0 0 0

0 0 0 0

0 + 1// /. 0// 0-

(84)

The three matrices in Eqs. (83) and (84), when supplemented with the parameter ∆ quantifying tracked/detected bandwidth, completely specify the photon counting statistics for both emission and detection schemes presented in Secs. IV and V. Numerical diagonalization of L allows the L, D+, and D− matrices to be transformed to the basis of L eigenstates, which then allows for immediate (if tedious) calculations via the formulae derived and presented in Appendices A and B. Thus, for the two-level chromophore, four physical parameters completely determine the photon counting statistics: the Rabi frequency Ω0, the spontaneous emission rate Γ0, the detuning frequency δ L , and the observation bandwidth ∆. For brevity of notation, we will henceforth set Γ0 ≡ Γ and Ω0 ≡ Ω when discussing the two level chromophore; there are only a single possible spontaneous emission rate and a single possible Rabi frequency in the absence of nuclear motions coupled to the electronic transition.

A. Statistics in the long observation time limit

Mollow54 first showed that, in the case of strong laser field (Ω > Γ) or off-resonant excitation (δ L > Γ, Ω), the emission spectrum of a driven two-level chromophore shows a triplet shape, while for low intensity near resonant excitation, the spectrum has a single peak. We consider both of these regimes and compare the emission statistics of Sec. IV and detection statistics of Sec. V. In this section, we consider only the limit of long observation times and explicitly avoid the more problematic artifacts of the emission based scheme. Short time observations will be discussed in Sec. VI B. The case of longtime emission statistics was considered in our preliminary report66 and the primary focus of this section is to contrast the results obtained via the emission and detection approaches. In order to demonstrate the effect√of the detector, we show in Fig. 2 single peak results (Ω = Γ/ 2) for both the emission and detection schemes; both the spectrum itself (intensity) and Q parameter spectrum are displayed. In the left panels, we used a detector of width ∆ = Γ/10, comparable to the width of the peak. In the right panels, we used ∆ = Γ/104, which is orders of magnitude narrower than the peak. In Fig. 3, we present similar plots√ for the case of a triplet spectrum, corresponding to Ω = 5Γ/ 2. Figures 2 and 3 display clear differences between the spectra as obtained via the two different prescriptions (emission vs. detection). The primary effect at play here is signal “leakage” between adjacent detection intervals. The emission scheme of Sec. IV precisely assigns photons to a single


174104-14


J. Chem. Phys. 142, 174104 (2015)

FIG. 2. The singlet spectrum and Q parameter spectrum for finite frequency resolution. √The Rabi frequency is given by Ω = Γ/ 2. The laser frequency is on resonance with the system, ω L − ω 0 = 0. In the left panels, we show the intensity (top) and the Q parameter (bottom) as measured by a detector of width ∆ = Γ/10. The triangles (red) correspond to calculations including the response function of the detector, and the circles (blue) correspond to the naive calculations based solely on photon emission. In all the calculations here, we took the limit t → ∞ analytically. In the right panels, we show the same quantities, but for a detector of width ∆ = Γ/104. In order to show the fine structure in these plots, we truncated the y axis, and the insets show the full data (non-truncated). In the insets, the x axis shows the frequency (the same units as in the main figure), and the y axis shows the intensity (not scaled).

frequency bin, however, the Lorentzian response function associated with physical detection is inherently less discriminate. A physical detector can register photons falling directly within its ∆-wide frequency window and also photons with frequencies considerably more distant from the window center. This convolution between properties of the physical field and response of the detector leads to an “over-counting” in the detected intensity spectrum as multiple windows can all lay

claim to the same photon. Hence, the detection spectrum is shifted to substantially greater magnitude than the corresponding emission spectra. The case of the Q parameter spectrum is more subtle, although the disparity between emission and detection cases derives from the same origin involving convolution between the field and detector response. As discussed in our original report,66 binning the Q spectrum can result in complete

FIG. 3. The triplet spectrum and Q parameter for finite frequency resolution. The √ Rabi frequency is given by Ω = 5Γ/ 2. All other parameters are the same as in Fig. 2.


174104-15


reversal of sign, proceeding from broadband measurement with only a single frequency bin to multiple narrow frequency bins. The Q parameter is a non-linear function of the physical field and there is no requirement that the detector response lead to a simple “overcounting” type effect. Indeed, the physically detectable Q parameter spectrum is seen to include both regions that are higher than and lower than the emission Q spectrum based solely on emission. Although one could argue that the emission and detection plots in Figs. 2 and 3 are qualitatively similar, they contain obvious differences and speak to the importance of using the proper theory including detector response in the modeling of any future experiments. Fig. 4 explicitly shows that the signal expected from detected photons is as much a function of the detector width as of the underlying physics of the physical field itself. Here, the two physical regimes from Figs. 2 and 3 are considered on the same plot. Only the middle frequency bin centered on ω0 is considered and the signals are plotted as a function of detector width δ. The singlet and triplet signals show a reversal in dominance at intermediate bin sizes; this is true for both the intensity and Q parameter. So, while the standard broadband Q parameter of the singlet regime is considerably more negative than the corresponding triplet regime quantity, at small bin sizes there is exactly the opposite behavior. Note that the plateau seen at large ∆ reflects convergence to broadband regime. As discussed in Secs. III–V, both emission and detection schemes converge to the same results in this particular limit. In Fig. 5, we show the correlations between photons emitted at different frequencies. The correlations were calculated using the naive emission-based approach (Eq. (32)) and the long time limit was taken analytically using the results of Appendix A. The two panels correspond to the singlet √ √ (Ω = Γ/ 2) and the triplet (Ω = 5Γ/ 2) parameters. The observation width was set to ∆ = Γ/10, which is narrow enough to resolve frequencies within the individual emission peaks. Both singlet and triplet cases display positive correlations along the ω, −ω diagonal reflecting the importance of energy conservation;58,66 fluctuations with an enhanced number of +ω emissions are correlated with an enhanced

J. Chem. Phys. 142, 174104 (2015)

number of −ω emissions. As we have previously noted,66 the narrow band Q parameter centered at the excitation frequency, corresponding to the ω1L = ω2L pixel at the origin of the two panels (n.b. the definition in Eq. (33)), is strongly positive and reflects bunching statistics. Anti-correlations between the numbers of photons at other frequencies are the origin of the anti-bunching observed in broad band photon counting statistics.72 In Fig. 6, we plot the correlations between photons detected at different frequencies accounting for the detectors’ response functions (see Eq. (77)). These plots parallel the analogous calculations presented in Fig. 5, which neglect the detectors. Although the general features of these plots are similar to those determined via consideration of photon emission, there are obvious differences. Most apparent are the differences in the vicinity of the ω1L = 0 and ω2L = 0 axes. Unlike the emission scheme of Fig. 5, the detector response allows the large contribution from the elastic or “coherent”54,62 delta-function at ω1L = 0 to bleed over into nearby frequency regions. The mixture of coherent and incoherent contributions within all detector windows extends the positive and negative correlations to the whole range of the spectrum. The appreciable tails of the detector response also weaken the correlations along the off-diagonal as the energy conservation condition is smeared out relative to the emission-based statistics. We note that, while statistics based on emission are not directly experimentally observable, the physical picture afforded by the emission calculations is somewhat more transparent than that obtained via direct calculation of detected statistics. The convolution of the instrument response can only obscure the behavior of the radiated field. The emission and detection based schemes are complementary to one another in the longtime limit. Calculations based on photon detection provide a direct means to calculate experimentally observable phenomena, while the emission scheme may provide the clearest window into the underlying physics. B. Statistics at short times

In Sec. IV, a possible problem with the naive emissionbased approach to calculating photon statistics was

FIG. 4. The intensity (left panel) and the Q parameter (right panel) as measured by a detector centered at ω = ω L as a function of the detector width ∆. The solid line (blue) correspond √ to the singlet spectrum Ω = Γ/ 2, and the dashed line (red) correspond to the √ triplet spectrum Ω = 5Γ/ 2. In both cases, the laser is on resonance with the two-level system (ω L − ω 0 = 0), and the infinite time limit (t → ∞) was taken analytically. We see that the intensity saturates as the width of the detector becomes larger than the width of the emission spectrum. For the Q parameter, a clear transition between positive and negative values is seen.


174104-16


J. Chem. Phys. 142, 174104 (2015)

FIG. 5. The correlation function, C ∆=Γ/10 (ω 1, ω 2, t → ∞), calculated according to the naive approach. ∆ = Γ/10; each pixel represents a single calculation and √ is ∆ × ∆ in size, reflecting the actual size and location of the frequency windows. The left panel corresponds to the singlet spectrum Ω = Γ/ 2, and the right √ panel corresponds to the triplet spectrum Ω = 5Γ/ 2. In both cases, the laser is on resonance with the two-level system (ω L − ω 0 = 0), and the infinite time limit (t → ∞) was taken analytically.

anticipated. The formalism presented there for the time derivative of the first moment of photon number (i.e., the intensity) has the appearance of a quantum generalization of the time-dependent spectrum introduced by Page75 and Lampard76 and it is thus expected that the expression for I∆(ω,t) is not restricted to be always positive.77 It is readily verified that the emission-based intensity does take on negative values at short times. In Fig. 7, we plot the time evolution of both the singlet spectrum and the triplet spectrum from Sec. VI A under resonant excitation conditions. Although the spectra converge to the fully positive results of Sec. VI A at sufficiently long times, the early time results are highly oscillatory as a function of detuning frequency and include negative regions. This behavior

is restricted to the emission-based calculation. The detected intensity, d I∆(ω,t) is independent of time and is always nonnegative. Whether one views negative values of I∆(ω,t) as being entirely unphysical, or simply unmeasurable, is really a matter of interpretation. If one insists on associating I∆(ω,t) with an intensity of the electric field, then it is true that this quantity should be non-negative. However, returning to our derivation of Sec. IV, we see that I∆(ω,t) was never associated with a field intensity, but rather was derived as the time derivative of the number of photons emitted into the field by the chromophore. Insofar as some emitted photons may be subsequently reabsorbed, it is not expected that I∆(ω,t) must remain strictly

FIG. 6. The correlation function, dC ∆=Γ/10 (ω 1, ω 2, t → ∞), as measured by detectors of width ∆ = Γ/10. The left panel corresponds to the singlet spectrum √ √ Ω = Γ/ 2, and the right panel corresponds to the triplet spectrum Ω = 5Γ/ 2. In both cases, the laser is on resonance with the two-level system (ω L − ω 0 = 0), and the infinite time limit (t → ∞) was taken analytically.


174104-17


J. Chem. Phys. 142, 174104 (2015)

FIG. 7. The time evolution of the singlet (left) and triplet (right) spectra for frequency resolution of√width ∆ = Γ/102 using the naive √ approach based on photon emission without the detector response function. For the singlet, we set the Rabi frequency Ω = Γ/ 2 and for the triplet Ω = 3Γ/ 2. In both cases, the laser is on resonance with the two-level system (ω L − ω 0 = 0). At short times, negative “intensities” are obtained.

positive. The real problem with I∆(ω,t) is that there would not appear to be any means to directly measure it. Physical detectors measure true intensities, which are non-negative (registered photons are removed from the system by the detector). The scheme presented in Sec. V for calculation of the detected intensity, d I∆(ω,t) is thus to be preferred for comparison to experiment, not only for the intensity itself, but for all higher order statistics as well. Therefore, in what follows, we only present the results of the scheme accounting for the response function of the detector. Fig. 8 displays the time evolution of the Q parameter for a very broad detector, ∆ = 1000Γ. A detector of this width yields (with insignificant discrepancies) the traditional Mandel Q parameter, in which all photons are counted; the results presented here are indistinguishable from the traditional calculations for broad band detection72 and serve as a validation of the computational approach. For the singlet parameters, the Q parameter monotonically converges to its long time limit while for the triplet parameters the convergence is not monotonic and the Q parameter oscillates before it converges. For the very broad detector used in this figure, the oscillations

reflect the system parameters showing a frequency of ≈Γ. The inset displays the short-time behavior of the main figure to emphasize the oscillations. The time evolution of the generalized Q parameter, as would be measured by detectors of width ∆ = Γ/10, is shown in Fig. 9. At short times, the generalized Q is small at all frequencies reflecting nearly Poissonian statistics for small average photon number. At later times, Q converges to the long-time limit seen in Fig. 3. The figure also suggests that convergence requires times t >> 1/∆ = 10/Γ, rather than the shorter time scale introduced by the system (∼1/Γ). The fact that observation times must exceed both 1/∆ and intrinsic system time scales to achieve convergence is seen more clearly in Fig. 10. Observation times shorter than (or comparable to) either 1/Γ or 1/∆ lead to statistics significantly different from the long-time results (as indicated by t = 1000/Γ in the figure). When the observation time is shorter than system time scales, the statistics are different from long observation times regardless of detector width. For times exceeding system time scales, convergence to the long-time asymptotics can be achieved, but only if the time exceeds 1/∆. The

FIG. 8. The time evolution of the Q parameter as recorded by a detector of width ∆ = 1000Γ centered at the laser √ frequency. The Rabi frequency corresponding to √ the singlet is Ω = Γ/ 2 (solid blue) and for the triplet spectrum, Ω = 5Γ/ 2 (dashed red). The intensity is constant in time and equal to 0.25 for the singlet and ≈0.48 for the triplet. The laser is on resonance with the two-level system (ω L − ω 0 = 0).

FIG. 9. The Q parameter as recorded by a detector of width ∆ = Γ/10. The different symbols correspond to different times as denoted in the figure. We have √used the Rabi frequency corresponding to the triplet spectrum Ω = 5Γ/ 2. The laser is on resonance with the two-level system (ω L − ω 0 = 0).


174104-18


J. Chem. Phys. 142, 174104 (2015)

FIG. 10. The Q parameter vs. the width of the detector at different times. The left panel corresponds to the singlet spec√ trum (Ω = Γ/ 2) and the right√panel to the triplet spectrum (Ω = 5Γ/ 2). For both panels the laser is on resonance with the two-level system (ω L − ω 0 = 0).

figure also shows that, for a narrow band detector the Q parameter is positive due to the correlations imposed by the detector, while for broad band detector it is negative reflecting the sub-Poissonian statistics of the photons from a single emitter.

Fig. 11 displays the time evolution of the correlation function (Eq. (77)) for the case of the triplet spectrum. As was the case for the single frequency observables discussed above, the correlations grow into the long-time results of Fig. 6 but require observation times longer than 1/∆ to achieve convergence.

FIG. 11. The time evolution √ of the correlation function as recorded by a detector of width ∆ = Γ/10. We have used the Rabi frequency corresponding to the triplet spectrum Ω = 5Γ/ 2. The laser is on resonance with the two-level system (ω L − ω 0 = 0). (a) shows the correlation function at t = 5/Γ, (b) at t = 10/Γ, (c) at t = 100/Γ, and (d) at the limit t → ∞ corresponding to the results of Fig. 6.


174104-19


C. Off-resonance excitation

The results presented above have all assumed resonant excitation by the laser field with ω L = ω0. However, the methods presented in this work are not limited to this case and can be also used to calculate the photon statistics for offresonance excitation. To demonstrate the case of off-resonant excitation, we chose the system parameters similar to the pioneering experiment of Aspect et al.58 In this experiment, atoms were excited far from resonance such that the detuning frequency was much larger than all system inverse time scales. As mentioned earlier, in this case the spectrum shows a triplet–three peaks corresponding to the laser frequency and two side bands at frequencies ω L ± ω L0. It was found that there is a strong correlation between the photons emitted at the two side bands. In our calculations, the Rabi frequency is set to Ω = 400Γ, and the detuning frequency is given by ω L0 = 20 000Γ. The limit t → ∞ was taken analytically. In Fig. 12, we show the intensity and Q parameters for detectors of width ∆ = 50Γ and ∆ = 500Γ. These widths were chosen to ensure that all the photons emitted from a side band are detected by a detector centered at the side band frequency. For these parameters, the coherent peak is 104 times larger than side bands’ peaks. Therefore, in the top panels, showing the intensity, the y-axes were truncated such that the two side bands are visible. The insets of these panels show the full triplet spectra. The bottom panels show the corresponding Q parameter which is very small but exhibits significant peaks at

J. Chem. Phys. 142, 174104 (2015)

the side bands. Effectively, the narrow-band photon statistics is Poissonian at all frequencies. In Fig. 13, the correlation function for this system is displayed. The strong temporal correlations between emitted photons in the opposing side bands58 are reflected in the correlation function, while there is no correlation between the photons in the central peak and those at the side-bands. These results agree with the measurements of Aspect et al.58 and reflect the energy conservation between adsorbed and emitted photons; if a photon comes out of the atom with significantly higher energy than the incoming photons, that emission will be accompanied by another emission with a photon at lower energy to enforce energy conservation. The effect can be seen if the detectors are simply in the vicinity of the two side-bands and show some overlap, but is strongest if the detectors are centered exactly at ω = ω L ± ω L0.

VII. RESULTS FOR A FOUR-STATE CHROMOPHORE

To demonstrate the applicability of our method to more complex systems, we consider the case of a chromophore coupled to a two-level system (TLS). This model is of both theoretical and practical interest. A two-state chromophore coupled to a single TLS yields a four-state quantum system, which is considerably more complex than a simple two-level chromophore, but remains simple enough to study numerically without difficulty. Physically, TLS’s are known to play an

FIG. 12. The spectrum and the Q parameter for off-resonance excitation. We set the parameters to have the same ratios as in Ref. 58. The Rabi frequency is set to Ω = 400Γ, and the detuning frequency is given by ω L0 = 20 000Γ. The limit t → ∞ was taken analytically. In the left panels, we show the intensity (top) and the Q parameter (bottom) as measured by a detector of width ∆ = 50Γ. In the right panels, we show the same observables for a detector of width ∆ = 500Γ. For these parameters, the coherent peak of the intensity is ∼104 times larger than the side bands peaks, and therefore, we truncated the intensity axes in the main panels. In the insets, we show the full (un-truncated) versions.


174104-20


J. Chem. Phys. 142, 174104 (2015)

excited electronic states. Within a traditional stochastic lineshape model,86,87 this picture would correspond to a chromophore undergoing two-state hop spectral diffusion as driven by phonon-assisted tunneling of the TLS. The full quantum mechanical treatment used here is somewhat more general than this due to the mixing between chromophore and TLS states within the Hamiltonian, but it is convenient to think of this model as a proper quantum treatment of spectral diffusion. The chromophore density matrix for this system, σC H (t) is 4 × 4, corresponding to super-operators of dimension 16 × 16. The super-operator L governing dynamics of the chromophore density matrix (see Sec. IV A) is comprised of several additive pieces FIG. 13. The correlation function for off-resonance excitation corresponding to the same system as in Fig. 12. The detector’s width is ∆ = 500Γ. Correlations between photons in the two opposing side bands are clearly seen (as was measured in Ref. 58), while photons at different frequencies (including the coherent photons) are not correlated with others. The inset shows a 3D version of the main contour plot.

important role in both the thermal physics and spectroscopy of low temperature glassy systems82–84 and the model we present can be viewed as representing the behavior of a single dye molecule embedded in a glassy matrix, which we emphasized in our original report.66 The mathematical model we adopt here for a coupled two-level chromophore and TLS has been repeatedly described in some detail in the literature.37,84,85 For this reason, we will present only the equations necessary to specify the pertinent super-operators L, D+, and D− necessary to calculate photon statistics. Further details behind the derivation of these equations can be found in Ref. 37, where we considered exactly the same model in the context of broadband emission statistics and a much simplified and approximate treatment of spectrally resolved emission statistics. Physically, the model generalizes that of the two-level chromophore considered in Sec. VI by introducing two nuclear states within both the excited and ground state manifolds of the chromophore (see Fig. 14) and coupling these four states to a thermal bath representing the long-wavelength phonons of the solid matrix. The states |a⟩ and |b⟩ are associated with the ground electronic manifold and the states |c⟩ and |d⟩ with the excited electronic manifold. In the notation of Sec. III, a and b are associated with the ng index for ground nuclear states, whereas c and d are associated with the ne index for

L = W + R + L E + LΓ .

(85)

The operator W is associated with the unperturbed dynamics of the “chromophore” component of the total system, which in this case is comprised of the coupled bare chromophore CH is fully and TLS pieces. Using the language of Sec. III, H 84,85 specified by ) ( A α ~ω0 σ TLS + + − Hg = − 2 2 4r 3 z ( ) ~ω0 A α He = + + + σ TLS + 2 2 4r 3 z

J TLS σ , 2 x (86) J TLS σ . 2 x

Here, A and J are the asymmetry and tunneling matrix elements, respectively, for the TLS, and σ zTLS, σ TLS are Pauli x matrices in the basis of the TLS’s localized “minima” states. ω0 is the chromophore electronic transition frequency in the absence of interactions. α is a chromophore-TLS coupling constant and the strain-field mediated interaction between the two is of a dipolar nature—leading to the displayed 1/r 3 dependence in the physical distance between the two. DiagonalCH leads to the four eigenstates |g⟩|a⟩, |g⟩|b⟩, |e⟩|c⟩, izing H and |e⟩|d⟩, increasing in energy alphabetically. In this basis, Eq. (86) can be written as Hg = ~ωa |a⟩⟨a| + ~ωb |b⟩⟨b|, He = ~ωc |c⟩⟨c| + ~ω d |d⟩⟨d|,

(87)

with the frequencies

FIG. 14. Schematic description of the energy levels of the composite chromophore-two level system. The solid lines with arrows represent all the possible radiative transitions while the wavy curves represent all the possible thermal transitions.

1 1 ω a = − ω0 − 2 2~



J 2 + (A − P)2,

1 1 ω b = − ω0 + 2 2~



J 2 + (A − P)2,

1 1 ω c = + ω0 − 2 2~



J 2 + (A + P)2,

1 1 ω d = + ω0 + 2 2~



J 2 + (A + P)2,

(88)

and we have set P ≡ 2rα3 . Written in the |a⟩, |b⟩, |c⟩, |d⟩ basis (in super-operator elements specified below, indices with subscript g take the values a or b and indices with subscript e


174104-21


J. Chem. Phys. 142, 174104 (2015)

take the values c or d), W is diagonal, with elements

with

Wn g m g ;n g m g = −iω n g m g , Wn e m e ;n e m e = −iω n e m e , Wn e m g ;n e m g = −i(ω n e m g − ω L ),

Γ0 ≡ (89)

Wn g m e ;n g m e = −i(ω n g m e + ω L ). As introduced in Secs. III and IV A, σ CH (t) is a rotating frame density matrix within the RWA approximation. This is why ω L enters into W. The dynamics due to the driving laser field follow immediately from the form of H I e in Sec. III to yield L E , with nonzero elements given by

i LnEg m g ;k e l g = −LkEe l g ;n g m g = Ω0 ng | k e δ m g l g , 2

i LnEg m g ;k g l e = −LkEg l e ;n g m g = − Ω0 l e | mg δ n g k g , 2 (90)

i E E Ln e m e ;k e l g = −Lk e l g ;n e m e = − Ω0 l g | me δ n e k e , 2

i LnEe m e ;k g l e = −LkEg l e ;n e m e = Ω0 ne | k g δ m e l e . 2 Γ

R and L contain the non-Hamiltonian dissipative contributions to the dynamics of σ CH (t) obtained by integrating out the acoustic phonons of the bath and photons of the quantum radiation field, respectively. The elements of R are calculated using the full Redfield formalism74 described previously in Refs. 37 and 85. For example, a sampling of non-zero elements is provided by 1 , 1 − e−β~ω e 1 , = e β~ω g Rbb;aa = Cωg J 2 1 − e−β~ω g   ωg 1 ωe Rca;db = Rcc;dd + Raa;bb , 2 ωg ωe   ωg 1 ωe Rdb;ca = Rdd;cc + Rbb;aa , 2 ωg ωe ωg ≡ ωb − ωa , ωe ≡ ω d − ωc ,

Rcc;dd = e β~ω e Rdd;cc = Cωe J 2 Raa;bb

(91)

however, there are many more non-zero contributions, which we do not explicitly list here. Rcc;dd and Raa;bb are the thermal transition rates from d to c and from b to a, respectively. The rates of the reverse transitions are seen to obey detailed balance. C is a collection of constants incorporating the coupling strength between the TLS and the bath, which is typically taken as a parameter used to fit an experiment rather than estimated from first principles.88 Of course, the top two lines simply express the phonon-assisted transition rates from state d to c and b to a, as expected. Other elements follow similarly. LΓ involves state-to-state spontaneous emission rates from states ne = c or d to states mg = a or b, 2 Γn e m g = Γ0 ⟨ne | mg ⟩ ,

(92)

− 2 3 → ωeg µ 0 3πε 0~c3

(93)

as previously defined. The non-zero elements of LΓ are given by  LΓn e n e ;n e n e = − Γn e m g , mg

LΓm g m g ;n e n e = Γn e m g ,

(94)

 1  Γik + Γ j k +/ , LΓi j;i j = − *. 2 k,i k, j , with all other elements being zero. The i, j, k indices in the third line run over all states irrespective of ground or excited affiliation, however, Γi j = 0 unless the i and j indices have ne and mg character, respectively. The elements of D± follow immediately from the formalism specified in Sec. IV A. The non-vanishing elements are given by D+in g ;i m e = ⟨me | ng ⟩, D−n g i;m e i = ⟨ng | me ⟩,

(95)

where the index i runs over states of both excited and ground character. As a first example, we consider the case of a slow TLS modulation, and strong coupling between the TLS and chromophore. (The case of weak TLS-chromophore coupling is closely related to a purely stochastic description and is somewhat less interesting.37) To orient ourselves with respect to the relevant physics, we first calculate the absorption spectrum and traditional broadband Mandel Q parameter as a function of detuning frequency in Fig. 15. We emphasize that these quantities involve photon counting statistics without any resolution of the emitted/detected photons and the results can be obtained either via the ∆ → ∞ limit of the equations presented in this work, or via the simpler methods we have discussed previously.34,37,65 Analytically, both approaches are equivalent and we have verified that the numerical results are identical as well, though the two routes to the numerical results are quite different. The model parameters are Γ0 = 40 Ms−1, T = 1.7 K, A = 0.006 K, J = 0.008 K, α = 3.75 × 1011 nm3 s−1, C = 3.9 × 108 K−3 s−1, and r = 5.72 nm. These parameters, which were previously considered in Ref. 37, correspond to thermal transition rates (from state a to b or c to d) of the order of ≈ 4 × 104 s−1. Figure 15 clearly shows 4 separate molecular transitions; depending on the detuning frequency, the laser is exciting one of the following transitions: a → c, a → d, b → c, or b → d. The relative magnitude of the line shape peaks results primarily from the probability splitting to excite “diagonal” vs. “off-diagonal” TLS transitions as the electronic transition from g to e is made. Looking at Eq. (88), the diagonal transitions that would be allowed in the absence of TLS-chromophore coupling (P = 0) correspond to a → c and b → d, with the remaining two being off-diagonal. First, consider the Ω0 = 0.1 Ms−1 results. Under these weak


174104-22


J. Chem. Phys. 142, 174104 (2015)

FIG. 15. The absorption spectrum (left panel) and Q parameter (right panel) for a chromophore coupled to a two-level system. The parameters, which are specified in the text, correspond to a slow modulation of the two-level system-C = 3.9 × 108. In the left panel, the blue line with square symbols corresponds to the intensity for Rabi frequency Ω0 = 105, and the black line with diamond symbols corresponds to the intensity divided by 104 for Ω0 = 4 × 107. In the right panel, the red line with circle symbols corresponds to the Mandel Q parameter multiplied by 103 for Ω0 = 105, and the purple line with star symbols corresponds to the Mandel Q parameter for Ω0 = 4 × 107.

excitation conditions, one expects the strength of the transition to be proportional to the square of the Rabi frequency associated with that transition.69 The ratio of peak heights for the a → c vs. a → d transitions in the line shape, for example,

should thus look like |⟨a | c⟩|2/|⟨a | d⟩|2 (∼ 1.9 for present parameters), which is in good agreement with the results. The mirror symmetry about zero detuning is due to the small (relative to kT) energy difference between the states a and

FIG. 16. The emission intensity (top panels), Q parameter (medium panels), and correlation function (bottom panels) for a chromophore coupled to a two-level system. The width of the detector is ∆ = Γ0. The system parameters, which are specified in the text, correspond to a slow modulation of the two-level systemC = 3.9 × 108 and a relatively small Rabi frequency-Ω0 = 105. Each column corresponds to a different excitation frequency as specified in the figure. The y-axis of the middle row and the color map of the bottom panels correspond to the Q parameter and correlation function multiplied by 103, with the exception of the middle panels in these rows where the y-axis and colormap correspond to the values multiplied by 106 (the additional 103 factor is indicated within these panels).


174104-23


b—there is no significant preference to find the system in one TLS state over the other and the scenario just described for transitions out of a is repeated for transitions out of b. With respect to the Q parameter, detuning to the favored diagonal transitions is qualitatively different from detuning to offdiagonal transitions. A system excited via a diagonal transition is more likely to relax (via spontaneous emission) back to its starting point than to end up in the other possible ground state. This is because the spontaneous emission rates for a given transition are proportional to the same squared matrix elements as are the excitation probabilities just discussed (see Eqs. (92) and (90)). The opposite holds true for off-diagonal excitations; these excitations are more likely to end up the other available ground state than where they began. This qualitative difference in preferred pathways manifests itself through a reversal of the sign of the Q parameter: diagonal detuning yields a bunching behavior (Q > 0) whereas off-diagonal excitations yield anti-bunching. Note however, that the effect is not very strong—the Q parameter values are all quite close to zero. When the system is pushed out of the weak excitation limit by increasing Ω0, the most obvious effects are the drastic increase in the rate of photon emission and power-broadening of the spectrum. At these powers, however, the individual transitions are still clearly recognizable. The Q parameter does qualitatively change, with a strong increase in magnitude and positive values observed over nearly all detuning frequencies.

J. Chem. Phys. 142, 174104 (2015)

In Figs. 16 and 17, we show the detected emission spectrum (upper panels), the detected Q parameter spectrum (middle panels), and the emission correlation function (bottom panels) for five different excitation frequencies (corresponding to the different columns). The two figures correspond to the two different Rabi frequencies studied in Fig. 15 and all other physical constants are the same as in Fig. 15. The emission spectrum is more complicated and strongly depends on the detuning frequency. For the weaker excitation case (Fig. 16, Ω0 = 105) and applied frequencies on resonance with a particular transition, the intensity spectrum consists of two peaks corresponding to emission from the state excited by the laser to the two ground states. TLS flipping is significantly slower than photon emission for these parameters and emissions out of the other (not directly excited) excited state are very rare. For instance, in the leftmost column, when the laser excites the population from state b to state c, there is significant intensity at frequencies corresponding to the transitions c → a and c → b. The peaks at the other possible transitions (i.e., the transitions d → b and d → a) are orders of magnitude smaller and therefore are not visible. For off-resonant excitation (ω L = ω0, middle column), the intensity is much lower and there are three peaks corresponding to the coherent part of the intensity (at the laser frequency) and sidebands associated with the nearest resonant transition with peaks at ω − ω L = ±(ω L − ωbc ). The spectrum of the Q parameter shows that the coherent photons

FIG. 17. Same as Fig. 16 but with a stronger excitation – Ω0 = 4 × 107. The color map of the bottom panels corresponds to the correlation function.


174104-24


J. Chem. Phys. 142, 174104 (2015)

FIG. 18. The absorption spectrum (left panel) and Q parameter (right panel) for a chromophore coupled to a two-level system. The parameters, which are specified in the text, correspond to a fast modulation of the two-level system–C = 3.9 × 1012. In the left panel, the blue line with square symbols corresponds to the intensity for Rabi frequency Ω0 = 105, and the black line with diamond symbols corresponds to the intensity divided by 105 for Ω0 = 4 × 107. In the right panel, the red line with circle symbols corresponds to the Mandel Q parameter multiplied by 105 for Ω0 = 105, and the purple line with star symbols corresponds to the Mandel Q parameter for Ω0 = 4 × 107.

(emitted at the frequency of the laser) exhibit bunching. The photons emitted in other frequencies show nontrivial statistics due to the off-resonance excitation of the corresponding transitions. The cross correlation panels show that for excitation

in resonance with one of the possible transitions, the only correlations are between the numbers of photons in the same frequency (the Q parameter). For the off-resonance excitation we see, in addition, correlations corresponding to two-photon

FIG. 19. The emission intensity (top panels), Q parameter (medium panels), and correlation function (bottom panels) for a chromophore coupled to a two-level system. The width of the detector is ∆ = Γ0. The system parameters, which are specified in the text, correspond to a fast modulation of the two-level system – C = 3.9 × 1012 and a relatively small Rabi frequency – Ω0 = 105. Each column corresponds to a different excitation frequency as specified in the figure. The color map of the bottom panels corresponds to the correlation function multiplied by 108.


174104-25


processes (the off-diagonal which corresponds to frequencies whose sum is equal to 2ω L ). It should be stressed that, although the second and third rows of Fig. 16 display some visually striking structures, the absolute magnitude is everywhere small in these plots and the statistics are close to Poissonian with limited correlations. Fig. 17 provides similar plots to Fig. 16, but for much stronger laser excitation: Ω0 = 4 × 107 s−1. Due to the higher Rabi frequency (stronger laser excitation), the spectra for both on-resonance and off-resonance excitations show three peaks rather than the two seen for the weaker laser. The three peaks are due to the same effect yielding the Mollow triplet in simple TLS. Note that for the four-level chromophore, the height of the peaks is determined by the properties of the chromophore and the excitation frequency. The Q parameter shows significant super-Poissonian statistics of the resonant photons while the incoherent photons show only small deviation from Poissonian statistics. Excitation and emission are now both significantly faster than TLS flipping and it is now expected to see multiple resonant photon emissions before a TLS flip event; this leads to bunching. The bunching effect is particularly pronounced for excitation of the diagonal transitions which favor repeated identical emission events the most strongly. The cross-correlation panels show significant correlations between the number of resonant photons and the number of off-resonant photons sharing a common

J. Chem. Phys. 142, 174104 (2015)

excited state with the resonant transition (for example, in the left most panel, the excited transition is the c → b and it shows strong correlations with the number of photons with frequency corresponding to the transition c → a). This makes sense as the emission of such a non-resonant photon likely coincided with recent resonant excitations as the system cycled up and down. When the coupling between TLS and environment is stronger, TLS dynamics increase and the phenomenon of mutual narrowing is observed in the absorption spectrum. Before turning to frequency resolved photon counting statistics for this case, the absorption spectrum and the broadband Q parameter are presented in Fig. 18. The figure shows that due to the strong coupling, the separate transitions are not distinguishable in the absorption spectrum or Mandel Q parameter. The results of Fig. 18 are similar to those presented in Ref. 37. For relatively weak laser excitation (Fig. 19), Ω0 = 105, the emission lineshape is almost identical to the absorption lineshape with the main difference being the additional coherent peak at the laser frequency in the emission lineshape. Note that the shape of the emission spectrum (with the exception of the coherent peak) is insensitive to the excitation frequency while it is within the width of the lineshape (the shift between the panels is because we measure the frequency relative to the laser frequency and this frequency varies between the

FIG. 20. Same as Fig. 19 but with a stronger excitation – Ω0 = 4 × 107. The color map of the bottom panels corresponds to the correlation function multiplied by 103.


174104-26


panels). The generalized Q parameter is practically zero everywhere, reflecting a Poissonian statistics. Similarly, there are no significant cross correlations due to the weak laser excitation and the large lag between successive emissions of photons. (Although the structure looks striking in the lower two rows of Fig. 19, the magnitude of the numbers is very small.) In Fig. 20, we show the same quantities for the case of stronger laser excitation, Ω0 = 4 × 107. The lineshape is almost the same as the lineshape for the case of weaker excitation (of course the amplitude here is much higher). The generalized Q parameter is still very small showing only minor deviations from Poissonian statistics and, similarly, the cross-correlations are weak. The stronger coupling to the thermal bath allows for rapid transfer of energy between the thermal and optical degrees of freedom; thereby weakening the cross-correlations due to energy conservation and multi-photon processes. An interesting question to pose is “how do the results of the detected frequency-resolved statistics compare to the results derived in Ref. 37, using a method that separately counts photons from different radiative transitions?” As we will see, the answer to this question is that the transition based approach can range from either a good to poor approximation of the proper detection based calculation, depending upon exactly what is being studied. In Fig. 21, the results of Fig. 5 from Ref. 37 are reproduced and overlaid with results for the detected emission spectra and Q parameter spectra; the

J. Chem. Phys. 142, 174104 (2015)

physical system under study is a TLS-chromophore system identical to that introduced above and five different exciting laser frequencies are considered (marked in the figure using black arrows). The excitation frequencies, from left to right, are ω0 + ωcb , ω0 + ωca , ω0, ω0 + ω db , and ω0 + ω da . The spontaneous emission rate and the Rabi frequency are Γ0 = 40 Ms−1 and Ω0 = 4 Ms−1, respectively. The model parameters are as follows: T = 1.7 K, A = 0.006 K, J = 0.008 K, α = 3.75 × 1011 nm3 s−1, C = 3.9 × 108 K−3 s−1, and r = 5.72 nm. In the case of the detected results, detectors of width ∆ = Γ0, centered at the frequencies of the allowed radiative transitions, are assumed. The figure demonstrates good agreement between the transition based scheme and measurable statistics, when the system is excited on resonance with one of the electronic excitation frequencies for the system. Indeed, for the intensity, the two schemes yield indistinguishable results at the resolution of these plots. Slight deviations are seen in the Q parameter spectrum, but the features are semi-quantitatively reproduced. However, when the excitation is off resonance, the two methods yield completely different results (see the middle column). For off-resonant excitation conditions, the coherent contribution to the emission is not aligned with any radiative transition—for the case studied in the middle column, the coherent emission is centered at ω0. Most of the photons are being emitted in the coherent peak and detectors placed at the transition frequencies only capture a small fraction of these

FIG. 21. A comparison of the intensity and the Q parameter as calculated using a generating function resolving specific atomic transition (denoted as Old, red diamonds)37 and using the frequency-resolved method (denoted as Detector, blue circles) presented in this paper. The parameters are specified in the text.


174104-27


J. Chem. Phys. 142, 174104 (2015)

photons in the tail of their response function. The detected intensity summed over all the “expected” transition frequencies is orders of magnitude lower than the intensity inferred from the transition based scheme. All emitted photons can be mathematically associated with individual transitions and the sum over all transitions is equivalent to broadband photon statistics. The coherent peak at ω0, which is largely missed by physical detectors placed at the transition frequencies, is captured by the monitoring of the transitions themselves and this explains the vastly higher intensities seen in the transition scheme as opposed to the detection scheme in the middle panel. Similar and related disparities are seen in the results for the Q spectrum. The transition based scheme of Ref. 37 can provide a good approximation to physical detection under the right circumstances. When all photons are being emitted at frequencies associated with the individual molecular transitions, the scheme works well (albeit with quantitative discrepancies in the higher order statistics), however, the method fails miserably when this criteria is not met. It should be emphasized that there would not appear to be any way to measure transition based statistics directly—for comparisons to experiment the approach presented in this work is strongly preferred, with the transition-based method serving as a semi-quantitative approximation when introduced in appropriate circumstances.

VIII. CONCLUSION

Single molecule spectroscopy has grown to become an indispensable tool in many fields of scientific research. The technique will see its utility enhanced further still when it becomes possible to routinely perform photon counting experiments with spectral resolution. As experimental efforts in this direction take shape, it is important to develop theoretical tools that will help in the interpretation of these measurements. The present work has extended the traditional analytical approaches of quantum optics to numerical computation schemes applicable to multi-state chromophores with coupling to a thermal bath. This makes it possible to directly calculate frequency-resolved photon counting observables for a broad class of models commonly employed in the study of condensed-phase spectroscopy.89 It is immediately clear upon inspection of Appendices A and B that the theoretical predictions of this work are exceedingly complex and must be evaluated numerically for all but the simplest possible cases (i.e., broadband statistics for a two level chromophore, etc.). However, it should be stressed that the underlying physical considerations contained in these equations are really no more complex than the conventional optical Bloch equations. Section IV provides the prescription for calculation of the system correlation functions that underly all the terms detailed in the appendices (as well as higher order statistics not explicitly considered here). These

correlation functions are evaluated as elementary matrix operations involving only the Liouville super-operator for density matrix evolution and super-operator representations for leftand right-multiplication by system dipole moment operators. Each of these correlation functions is thus rather simple— it is only the various time and frequency integrals (and the associated time ordering permutations within the correlation functions) that lead to the complex results. So, in principle, it is possible to calculate frequency-resolved photon counting statistics for any chromophore system with a known dynamics in the Markov approximation. The appendices explicitly lay out the prescription to calculate these statistics up to second order. Two different theoretical approaches have been presented in this work based upon analysis of emitted photon statistics and detected photon statistics, respectively. The former approach (adopted in our preliminary report66) has the advantage of theoretical simplicity, however, the resulting expressions are not directly comparable to realizable experimental measurements. Explicitly accounting for the physical process of spectral filtering at the detector leads to expressions that can be compared directly to experiment and the detection based approach is preferred for this purpose. That said, emission based statistics lead to an elementary interpretation of the generating function (Eq. (2)) and the long-time statistics generated from these results are qualitatively similar to the detection based scheme—the primary differences between the two being attributable to the long tails in the detector response function. The two approaches are complementary to one another as the emission based scheme may pick up on fine structure which is lost due to convolution with the detector response; comparing the two approaches makes it clear when such a “washing out” is due to the detector wings and not the underlying physics of the chromophore (e.g., as in Figs. 5 and 6). ACKNOWLEDGMENTS

This research was supported in part by the National Science Foundation (Grant Nos. CHE-1153096 and CHE1465162) and by a Grant from the GIF, the German-Israeli Foundation for Scientific Research and Development (2259/2010). APPENDIX A: DETAILS OF THE SECOND MOMENT CALCULATION IN SEC. IV C

In order to calculate the correlation function of Eq. (32), it is convenient to write the dipole moment operators using their deviation from their average value, similar to the approach used in the separation of the intensity into coherent and incoherent parts. Doing so and combining terms according to the degree of the correlation function they include, we rewrite the function f (t,ω1,ω2) as

f (t,ω1,ω2) = f 0 (t,ω1,ω2) + f 2 (t,ω1,ω2) + f 3 (t,ω1,ω2) + f 4 (t,ω1,ω2) .

(A1)

The functions introduced above are defined as follows: f 0 (t,ω1,ω2) = D

+ 2

D

− 2

t t t t eiω1L (u−w)eiω2L (v−y)Θ (v − u) dudv dwdy. 0

0

0

0

(A2)


174104-28


J. Chem. Phys. 142, 174104 (2015)

The second term consists of three contributions detailed below, f 2 (t,ω1,ω2) = f 2a (t,ω1,ω2) + f 2b (t,ω1,ω2) + f 2c (t,ω1,ω2) ,

(A3)

where f 2a (t,ω1,ω2)   t t t t

+ 2

− iω (y−v) iω (w−u)  − 1L 2L = 2 Re  D δD (w) δD ( y) e e Θ (w − y) dudv dwdy  ,   0 0 0 0

f 2b (t,ω1,ω2) = D

+

t t t t

D

−

δD + (u) δD − (w) eiω1L (y−u)eiω2L (w−v)dudv dwd y,

(A5)

δD + (u) δD− (w) eiω1L (y−v)eiω2L (w−u)dudv dwdy.

(A6)

0

0

0

(A4)

0

and f 2c (t,ω1,ω2) = D

t t t t

+

D

−

0

0

0

0

In the equations above, we introduced the notation δD ± = D ± − ⟨D ±⟩. f 3 (t,ω1,ω2) is given by   t t t t

+

+   D δD (v) δD − (w) δD − ( y)  . 2 Re   (  0 0 0 0 ) × eiω1L (y−v)eiω2L (w−u) + eiω1L (y−u)eiω2L (w−v) Θ (w − y) dudv dwd y 

(A7)

The last portion of f (t,ω1,ω2) to be specified is t t t t f 4 (t,ω1,ω2) = 0

0

0

dudv dwd y δD + (u) δD + (v) δD− (w) δD − ( y)

0

) ( × eiω1L (y−v)eiω2L (w−u) + eiω1L (y−u)eiω2L (w−v) Θ (v − u) Θ (w − y) . (A8)

In order to calculate the numerator of the correlation function, Eq. (32), we express N(ω,∆) (t) N(ω′,∆) (t) in a way that will enable us to subtract it from the second multivariate moment. Using Eq. (29), we write it as

N(ω,∆) (t)

ω+ ∆2

ω ′+ ∆2



Γ2 N(ω′,∆) (t) = 02 4π



dω1 ω− ∆2

dω2 [h (t,ω1,ω2) + h (t,ω2,ω1)] .

(A9)

ω ′− ∆2

To proceed, we again group terms with the same order of dipole operators correlation functions and express h (t,ω1,ω2) as h (t,ω1,ω2) = h0 (t,ω1,ω2) + h2 (t,ω1,ω2) + h22 (t,ω1,ω2) ,

(A10)

where h0 (t,ω1,ω2) = D

+ 2

D

− 2

t t t t 0

h2 (t,ω1,ω2) = D

+

0

0

t t t t D

−

0

0

0

eiω1L (u−w)eiω2L (v−y)Θ (v − u) dudv dwdy,

(A11)

δD + (u) δD − (w) eiω1L (y−v)eiω2L (w−u)dudv dwdy,

(A12)

0

0

and 1 h22 (t,ω1,ω2) = 2

t t t t

0

0

0

0

δD + (u) δD − ( y) δD + (v) δD − (w) eiω1L (y−u)eiω2L (w−v)dudv dwdy.

(A13)


174104-29


J. Chem. Phys. 142, 174104 (2015)

Combining the various terms, the numerator of the correlation function C∆ (ω,ω ′,t) (Eqs. (32) and (57)) can be written as Γ2

C2(ω,ω′,∆) (t) − N(ω,∆) (t) N(ω′,∆) (t) = 02 4π

ω+ ∆2

ω ′+ ∆2





dω1 ω− ∆2

dω2 [M (t,ω1,ω2) + M (t,ω2,ω1)] ,

(A14)

ω ′− ∆2

where we defined M (t,ω1,ω2) = [ f 4 (t,ω1,ω2) + f 3 (t,ω1,ω2) + f 2a (t,ω1,ω2) + f 2b (t,ω1,ω2) − h22 (t,ω1,ω2)] .

(A15)

As mentioned in Sec. IV C, the expression above for the numerator of the correlation function is of little use from a practical point of view. In order to carry out the actual calculation and obtain numeric values for the correlation function for different parameters, we write the function M (t,ω1,ω2) as the sum of terms simple to calculate. The dipole operators’ correlation functions may only be calculated if they are expressed in a fully time-ordered manner. Then, one can use the formalism given at the beginning of Sec. IV to carry out the calculation. We start with f 4 (t,ω1,ω2). For this term, there are six possible time orderings which, after some algebra, can be reduced to the real parts of three terms, f 4 (t,ω1,ω2) = 2 Re [A (t,ω1,ω2) + K (t,ω1,ω2) + G (t,ω1,ω2)] ,

(A16)

where  t t4 t3 t2 A (t,ω1,ω2) = t 4=0 t 3=0 t 2=0 t 1=0

eiω1L (t1−t4)eiω2L (t2−t3)

dt 1dt 2dt 3dt 4 δD + (t 3) δD + (t 4) δD − (t 2) δD − (t 1) * iω1L (t1−t3) iω2L (t2−t4)+ . e ,+e -

(A17)

K (t,ω1,ω2) consists of two terms and is written as K (t,ω1,ω2) = K1 (t,ω1,ω2) + K2 (t,ω1,ω2) ,

(A18)

where  t t4 t3 t2 K1 (t,ω1,ω2) =

dt 1dt 2dt 3dt 4 δD + (t 2) δD + (t 4) δD− (t 3) δD − (t 1) eiω1L (t1−t4)eiω2L (t3−t2)

(A19)

dt 1dt 2dt 3dt 4 δD + (t 2) δD+ (t 4) δD − (t 3) δD − (t 1) eiω1L (t1−t2)eiω2L (t3−t4).

(A20)

t 4=0 t 3=0 t 2=0 t 1=0

and  t t4 t3 t2 K2 (t,ω1,ω2) = t 4=0 t 3=0 t 2=0 t 1=0

Similarly, we write G (t,ω1,ω2) as G (t,ω1,ω2) = G1 (t,ω1,ω2) + G2 (t,ω1,ω2) ,

(A21)

where  t t4 t3 t2 G1 (t,ω1,ω2) =

dt 1dt 2dt 3dt 4 δD + (t 1) δD + (t 4) δD − (t 3) δD − (t 2) eiω1L (t2−t4)eiω2L (t3−t1)

(A22)

dt 1dt 2dt 3dt 4 δD + (t 1) δD+ (t 4) δD − (t 3) δD − (t 2) eiω1L (t2−t1)eiω2L (t3−t4).

(A23)

t 4=0 t 3=0 t 2=0 t 1=0

and  t t4 t3 t2 G2 (t,ω1,ω2) = t 4=0 t 3=0 t 2=0 t 1=0

From the formalism for the calculation of the correlation functions in the Markov approximation (see Sec. IV), we know that the evolution operator U may be separated into the two components U0 and U1. As explained above (Sec. IV), the action of U0 yields the average of the terms to its right. Taking into account the fact that, by definition, ⟨D±⟩ = 0 allows us to express the correlation function using only U1 as the effective evolution operator. For this purpose, we use the notation ⟨{}⟩ for the correlation function, with the change of the full evolution operator to U1 only. For example,

+ D (t 2) D− (t 1) = Trace U1 (t 2 − t 1) D −σi ss D + , (A24) where we assumed t 2 ≥ t 1. Applying this notation to A, K , G allows us to write them as follows: A (t,ω1,ω2) = A a (t,ω1,ω2) + A b (t,ω1,ω2) ,

(A25)


174104-30


J. Chem. Phys. 142, 174104 (2015)

where  t t4 t3 t2 A a (t,ω1,ω2) =

dt 1dt 2dt 3dt 4

δD + (t 3) δD + (t 4) δD − (t 2) δD− (t 1)

t 4=0 t 3=0 t 2=0 t 1=0

eiω1L (t1−t4)eiω2L (t2−t3) * + iω 1L (t 1−t 3) iω 2L (t 2−t 4) e ,+e -

(A26)

and  t t4 t3 t2 A b (t,ω1,ω2) =

dt 1dt 2dt 3dt 4

δD − (t 2) δD− (t 1)

δD + (t 3) δD + (t 4)

t 4=0 t 3=0 t 2=0 t 1=0

eiω1L (t1−t4)eiω2L (t2−t3) * + iω 1L (t 1−t 3) iω 2L (t 2−t 4) . +e e , -

(A27)

Similarly, K1 (t,ω1,ω2) = K1a (t,ω1,ω2) + K1b (t,ω1,ω2) ,

(A28)

where  t t4 t3 t2 K1a (t,ω1,ω2) =

dt 1dt 2dt 3dt 4

δD + (t 2) δD + (t 4) δD− (t 3) δD − (t 1)

eiω1L (t1−t4)eiω2L (t3−t2)

(A29)

t 4=0 t 3=0 t 2=0 t 1=0

and  t t4 t3 t2 K1b (t,ω1,ω2) =

dt 1dt 2dt 3dt 4

+ δD + (t 2) δD − (t 1) δD (t 4) δD − (t 3) eiω1L (t1−t4)eiω2L (t3−t2).

(A30)

t 4=0 t 3=0 t 2=0 t 1=0

K2 (t,ω1,ω2) = K2a (t,ω1,ω2) + K2b (t,ω1,ω2) ,

(A31)

where  t t4 t3 t2 K2a (t,ω1,ω2) =

dt 1dt 2dt 3dt 4

δD + (t 2) δD+ (t 4) δD − (t 3) δD − (t 1) eiω1L (t1−t2)eiω2L (t3−t4)

(A32)

t 4=0 t 3=0 t 2=0 t 1=0

and  t t4 t3 t2 K2b (t,ω1,ω2) =

dt 1dt 2dt 3dt 4

δD + (t 2) δD− (t 1)

δD + (t 4) δD − (t 3)

eiω1L (t1−t2)eiω2L (t3−t4).

(A33)

t 4=0 t 3=0 t 2=0 t 1=0

And the last term of f 4 (t,ω1,ω2) is written as G1 (t,ω1,ω2) = G1a (t,ω1,ω2) + G1b (t,ω1,ω2) ,

(A34)

where  t t4 t3 t2 G1a (t,ω1,ω2) =

dt 1dt 2dt 3dt 4 δD + (t 1) δD + (t 4) δD − (t 3) δD − (t 2) eiω1L (t2−t4)eiω2L (t3−t1)

(A35)

dt 1dt 2dt 3dt 4 δD + (t 1) δD + (t 4) δD− (t 3) δD − (t 2) eiω1L (t2−t4)eiω2L (t3−t1).

(A36)

t 4=0 t 3=0 t 2=0 t 1=0

and  t t4 t3 t2 G1b (t,ω1,ω2) = t 4=0 t 3=0 t 2=0 t 1=0

f 3 (t,ω1,ω2) is a bit simpler and can be written as f 3 (t,ω1,ω2) = 2 Re [S (t,ω1,ω2) + P (t,ω1,ω2) + B (t,ω1,ω2)] ,

(A37)

where

+

×

S (t,ω1,ω2) = D

 t  t t3 t2

dt 1dt 2dt 3du

u=0 t 3=0 t 2=0 t 1=0

D+ (t 3) D − (t 2) D − (t 1)

(

) eiω1L (t1−t3)eiω2L (t2−u) + eiω1L (t1−u)eiω2L (t2−t3) ,

(A38)


174104-31


J. Chem. Phys. 142, 174104 (2015)

P (t,ω1,ω2) = D

+

 t  t t3 t2

dt 1dt 2dt 3du

D + (t 2) D − (t 3) D − (t 1)

u=0 t 3=0 t 2=0 t 1=0

(

× eiω1L (t1−t2)eiω2L (t3−u) + eiω1L (t1−u)eiω2L (t3−t2)

)

(A39)

and B (t,ω1,ω2) = D

+

 t  t t3 t2

dt 1dt 2dt 3du

) ( D + (t 1) D − (t 3) D − (t 2) eiω1L (t2−t1)eiω2L (t3−u) + eiω1L (t2−u)eiω2L (t3−t1) .

u=0 t 3=0 t 2=0 t 1=0

(A40) In the above expressions, we used the fact that the three-time averaged correlation function of the operators δD ± is equal to the modified correlation function with only U1 instead of the full evolution operator. Next, we write f 2a (t,ω1,ω2) as f 2a (t,ω1,ω2) = 2 Re [ψ (t,ω1,ω2)], where ψ (t,ω1,ω2) = D

+ 2

 t  t  t t2

0

δD − (t 2) δD − (t 1)

eiω1L (t1−v)eiω2L (t2−u)dt 1dt 2dv du.

(A41)

0 t 2=0 t 1=0

Similarly, we write f 2b (t,ω1,ω2) = 2Re [Ψ (t,ω1,ω2)], where Ψ (t,ω1,ω2) = D

+

 t  t  t t2 D

−

0

δD + (t 2) δD − (t 1)

eiω1L (y−t2)eiω2L (t1−v)dt 1dt 2d y dv.

(A42)

0 t 2=0 t 1=0

The last part of M (t,ω1,ω2) remaining to be written in a fully time-ordered manner is h22 (t,ω1,ω2), which we express as h22 (t,ω1,ω2) = 2 Re [ϕ (t,ω1,ω2) + φ (t,ω1,ω2)] ,

(A43)

where  t t4 t3 t2 ϕ (t,ω1,ω2) =

dt 1dt 2dt 3dt 4 t 4=0 t 3=0 t 2=0 t 1=0

 D + (t 4) D − (t 2) D + (t 3) D − (t 1) eiω1L (t1−t3)eiω2L (t2−t4) 

  + + D (t 4) D − (t 2) D + (t 1) D − (t 3) eiω1L (t2−t4)eiω2L (t3−t1)  iω (t −t ) iω (t −t ) + ×  + − − + D (t 1) D (t 4) D (t 3) D (t 2) e 1L 4 1 e 2L 2 3  + D + (t 4) D − (t 1) D + (t 3) D − (t 2) eiω1L (t2−t3)eiω2L (t1−t4)

(A44)

and  t t4 t3 t2 φ (t,ω1,ω2) =

dt 1dt 2dt 3dt 4 t 4=0 t 3=0 t 2=0 t 1=0

 D + (t ) D − (t ) D + (t ) D − (t ) eiω1L (t2−t1)eiω2L (t3−t4)  1 2 4 3 .

×  + + D (t 1) D − (t 2) D + (t 3) D − (t 4) eiω1L (t2−t1)eiω2L (t4−t3) Combining all the expressions above, we write M (t,ω1,ω2) as A a (t,ω1,ω2) + K1a (t,ω1,ω2) + K2a (t,ω1,ω2) + G1a (t,ω1,ω2)   +G2a (t,ω1,ω2) + S (t,ω1,ω2) + P (t,ω1,ω2) + B (t,ω1,ω2)   . M (t,ω1,ω2) = 2 Re   +A b (t,ω1,ω2) + K1b (t,ω1,ω2) + G1b (t,ω1,ω2) + ψ (t,ω1,ω2)    +Ψ (t,ω1,ω2) − ϕ (t,ω1,ω2)

(A45)

(A46)

Our main advantage so far is that we already subtracted all terms that could be subtracted before the actual calculation of the integrals and the correlation functions. In what follows, we will express the evolution operator in the correlation functions on the diagonal basis, so that we will be able to carry out all the integrals immediately. The price we have to pay for that is the summation of many terms. The details of the transformation from correlation functions to the sum of terms are given below. It is important to note that the solutions of the integrals given below are not valid for the singular cases (i.e., ω1 = ω2 and/or terms in This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Fri, 08 May 2015 07:26:14

174104-32


J. Chem. Phys. 142, 174104 (2015)

which two or more of the summation indices α, β, γ are equal). In those cases, the limit has to be taken properly. The final results have no singularities. We start with the terms including four-point correlation functions. A a (t,ω1,ω2) is written as ′ A a (t,ω1,ω2) = Aα;β;γ a1α;β;γ (t,ω1,ω2) , (A47) α;β;γ

where + − − Aα;β;γ = D+0,γ Dγ, β D β,α Dα,0

(A48)

and  t t4 t3 t2 a1α;β;γ (t,ω1,ω2) = t 4=0 t 3=0 t 2=0 t 1=0

eiω1L (t1−t4)eiω2L (t2−t3) dt 1dt 2dt 3dt 4e λγ (t4−t3)e λ β (t3−t2)e λ α (t2−t1) * iω1L (t1−t3) iω2L (t2−t4)+ . e ,+e -

(A49)

The solution of the integrals is a1α;β;γ (t,ω1,ω2) =

t(ω1L + ω2L + 2iλ γ ) (λ α − iω1L )(ω1L + iλ γ )(ω2L + iλ γ )(iλ β + ω1L + ω2L ) + e t λγ − e t λγ

(ω2L + iλ γ

)2(−λ

e−itω2L β + iω1L + λ γ )(λ α − iω1L + iω2L − λ γ )

e−itω1L (λ γ − λ α )(ω1L + iλ γ )2(−λ β + iω2L + λ γ )

+

e t(λ β −i(ω1L +ω2L ))(2λ β − i(ω1L + ω2L ) − 2λ γ ) (−λ α + λ β − iω2L )(λ β − i(ω1L + ω2L ))2(λ β − iω1L − λ γ )(λ β − iω2L − λ γ )

+

e t(λ α −iω1L )(2λ α − iω1L + iω2L − 2λ γ ) (λ α − iω1L )2(λ α − λ γ )(λ α − λ β + iω2L )(λ α − iω1L + iω2L − λ γ )

+

(ω1L + iλ γ )(iλ α + iλ β + 2ω1L + ω2L ) + (ω1L + iλ α )(iλ β + ω1L + ω2L ) (λ α − iω1L )2(ω1L + iλ γ )2(iλ β + ω1L + ω2L )2

+

(ω2L + iλ γ )(iλ α + iλ β + 2ω1L + ω2L ) + (ω1L + iλ α )(iλ β + ω1L + ω2L ) . (λ α − iω1L )2(ω2L + iλ γ )2(iλ β + ω1L + ω2L )2

(A50)

The leading order term in the long time limit a1α;β;γ (t,ω1,ω2) ∼

t(ω1L + ω2L + 2iλ γ ) . (λ α − iω1L )(ω1L + iλ γ )(ω2L + iλ γ )(iλ β + ω1L + ω2L )

K1a (t,ω1,ω2) and K2a (t,ω1,ω2) are combined together and written as ′ K1a (t,ω1,ω2) + K2a (t,ω1,ω2) = Kα;β;γ k1α;β;γ (t,ω1,ω2) ,

(A51)

(A52)

α;β;γ

where − + − Kα;β;γ = D+0,γ Dγ, β D β,α Dα,0

(A53)

eiω1L (t1−t4)eiω2L (t3−t2) dt 1dt 2dt 3dt 4e λγ (t4−t3)e λ β (t3−t2)e λ α (t2−t1) * iω1L (t1−t2) iω2L (t3−t4)+ . e ,+e

(A54)

and  t t4 t3 t2 k1α;β;γ (t,ω1,ω2) = t 4=0 t 3=0 t 2=0 t 1=0

The solution of the integrals is e

λβt

k1α;β;γ (t,ω1,ω2) =

( 1 λ 2β (−λ α +λ β +iω 1L )

+

e −i t (ω 1L −ω 2L ) (λ α −λ β −iω 2L )(i λ β +ω 1L −ω 2L )2

)

λ β + iω2L − λ γ − +

i(λ α +λ β −iω 1L ) λ 2β (ω 2L +i λγ )

+

λ α −iω 1L λ β (ω 2L +i λγ )2

+

λ α +λ β −2iω 1L +iω 2L (λγ −iω 1L )(i λ β +ω 1L −ω 2L )2

+

ω 1L +i λ α (ω 1L +i λγ )2(i λ β +ω 1L −ω 2L )

(λ α − iω1L )2


174104-33


J. Chem. Phys. 142, 174104 (2015)

e t(λ α −iω1L )

(

+

+

1 (λγ −λ α )(−λ α +λ β +iω 2L )

)

(λ α − iω1L )2 e t λγ

+

(

e −i t ω 2L (ω 2L +i λγ )2(−λ α +i(ω 1L −ω 2L )+λγ )

+

e −i t ω 1L (λγ −λ α )(ω 1L +i λγ )2

)

λ β + iω2L − λ γ ( t

+

1 (−λ α +λ β +iω 1L )(−λ α +i(ω 1L −ω 2L )+λγ )

1 (ω 1L +i λγ )(i λ β +ω 1L −ω 2L )

The leading order term in the long time limit is

−

i λ β ω 2L +i λ β λγ

λ α − iω1L ( t

k1α;β;γ (t,ω1,ω2) ∼

) .

1 (ω 1L +i λγ )(i λ β +ω 1L −ω 2L )

(A55)

−

i λ β ω 2L +i λ β λγ

)

λ α − iω1L

.

(A56)

In a similar way, we write G1a (t,ω1,ω2) + G2a (t,ω1,ω2) =

′ Gα;β;γ g1α;β;γ (t,ω1,ω2) ,

(A57)

α;β;γ

where − − + Gα;β;γ = D+0,γ Dγ, β D β,α Dα,0,

(A58)

eiω1L (t2−t4)eiω2L (t3−t1) dt 1dt 2dt 3dt 4e λγ (t4−t3)e λ β (t3−t2)e λ α (t2−t1) * iω1L (t2−t1) iω2L (t3−t4)+ . e ,+e

(A59)

and  t t4 t3 t2 g1α;β;γ (t,ω1,ω2) = t 4=0 t 3=0 t 2=0 t 1=0

The solution of the integrals is ( g1α;β;γ (t,ω1,ω2) = t −

+

+

1 λ β (ω1L − iλ α )(ω2L + iλ γ ) (λ α + iω2L )(ω1L + iλ γ )(−iλ β − ω1L + ω2L ) ( ) e λ β t λ 2β e−it(ω1L −ω2L ) + (λ β − i(ω1L − ω2L ))2 1

)

−

λ 2β (−λ α + λ β − iω1L )(λ β − i(ω1L − ω2L ))2(λ β + iω2L − λ γ ) ( ) itω itω e λ α t (λ e +iω1L )2 + (λ e +iω2L )2 α

1L

α

2L

(λ α − λ β + iω1L )(λ α + i(ω1L + ω2L + iλ γ ))

e t(λγ −i(ω1L +ω2L )) eitω1L (ω1L + iλ γ )2 + eitω2L (ω2L + iλ γ )2 + (ω1L + iλ γ )2(λ γ − iω2L )2(λ β + iω2L − λ γ )(λ α + i(ω1L + ω2L + iλ γ )) +

λ α + λ β − i(ω1L − 2ω2L ) 1 + (λ α + iω2L )2(λ γ − iω1L )(iλ β + ω1L − ω2L )2 (ω2L − iλ α )(ω1L + iλ γ )2(iλ β + ω1L − ω2L )

+

ω1L − i(λ α + λ β ) 1 + 2 . 2 λ β (λ α + iω1L )(ω2L + iλ γ ) λ β (λ α + iω1L )2(ω2L + iλ γ )

The leading order term in the long time limit is ( ) 1 1 g1α;β;γ (t,ω1,ω2) ∼ t − − . λ β (ω1L − iλ α )(ω2L + iλ γ ) (λ α + iω2L )(ω1L + iλ γ )(−iλ β − ω1L + ω2L ) Next, we write the terms including the three-time correlation functions ′ S (t,ω1,ω2) = Aα;β;0 s1α;β (t,ω1,ω2) ,

(A60)

(A61)

(A62)

α;β

where  t  t t3 t2 s1α;β (t,ω1,ω2) = u=0 t =0 t =0 t =0

eiω1L (t1−t3)eiω2L (t2−u) e λ β (t3−t2)e λ α (t2−t1) * iω1L (t1−u) iω2L (t2−t3)+ dt 1dt 2dt 3du. e ,+e -

(A63)

3 2 1 This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Fri, 08 May 2015 07:26:14

174104-34


J. Chem. Phys. 142, 174104 (2015)

The solution of the integrals is e−itω1L eitω1L + 2 2 λ αω1L (λ β − iω2L ) ω1L (ω1L + iλ α )(iλ β + ω1L + ω2L )

s1α;β (t,ω1,ω2) = −

eitω2L e−itω2L + 2 + iλ α )(iλ β + ω1L + ω2L ) ω2L (ω1L + iλ β )(iλ α + ω1L − ω2L ) e t(λ α −iω1L ) λ α ω1L −1 + eitω2L + ω2L −1 + eitω1L − iω2L −1 + eitω1L (ω1L − ω2L ) + λ αω1Lω2L (λ α − iω1L )(λ α − λ β + iω2L )(iλ α + ω1L − ω2L ) e t(λ β −i(ω1L +ω2L )) λ β ω1L −1 + eitω2L + ω2L −1 + eitω1L − iω1Lω2L eitω1L + eitω2L − 2 + ω1Lω2L (λ β − iω1L )(λ β − iω2L )(−λ α + λ β − iω2L )(iλ β + ω1L + ω2L ) 2 2 2 2 2 λ α ω2L ω1L + ω1Lω2L + 2ω2L + iλ β ω1L + 2ω2L + ω1Lω2L (λ β − i(ω1L + ω2L )) + 2 2 λ αω1Lω2L (λ α − iω1L )(λ β − iω2L )(iλ β + ω1L + ω2L )

+

+

2 ω2L (ω1L

2 ω2L (λ β

1 . − iω1L )(λ α − i(ω1L − ω2L )) ′ P (t,ω1,ω2) = Kα;β;0 p1α;β (t,ω1,ω2) ,

(A64) (A65)

α;β

where  t  t t3 t2 p1α, β (t,ω1,ω2) = u=0 t 3=0 t 2=0 t 1=0

eiω1L (t1−t2)eiω2L (t3−u) e λ β (t3−t2)e λ α (t2−t1) * iω1L (t1−u) iω2L (t3−t2)+ dt 1dt 2dt 3du. e ,+e

(A66)

The solution of the integrals is p1α, β (t,ω1,ω2) =

1 − e−itω1L −1 + eitω1L 1 − e−itω2L + + 2 2 2 λ αω1L (λ β + iω2L ) ω1L (ω1L + iλ α )(iλ β + ω1L − ω2L ) ω2L (ω2L − iλ β )(iλ α + ω1L − ω2L ) ie λ α t 1 − e−itω1L −1 + eitω2L − − 2 λ β ω2L (λ α − iω1L ) λ αω1L (λ α − iω1L )(λ α − λ β − iω2L ) −1 + eitω1L e t(λ β −i(ω1L −ω2L )) ω1L (λ β + iω2L )(−λ α + λ β + iω2L )(iλ β + ω1L − ω2L ) −1 + eitω2L e t(λ α −iω1L ) e λ β t −1 + eitω2L + + . ω2L (λ α − iω1L )(λ α − λ β − iω1L )(iλ α + ω1L − ω2L ) λ β ω2L (ω2L − iλ β )(λ α − λ β − iω1L )

+

(A67) ′ B (t,ω1,ω2) = Gα;β;0b1α, β (t,ω1,ω2) ,

(A68)

α, β

where  t  t t3 t2 b1α, β (t,ω1,ω2) = u=0 t 3=0 t 2=0 t 1=0

eiω1L (t2−t1)eiω2L (t3−u) e λ β (t3−t2)e λ α (t2−t1) * iω1L (t2−u) iω2L (t3−t1)+ dt 1dt 2dt 3du. e ,+e -

(A69)

The solution of the integrals is b1α, β (t,ω1,ω2) =

1 − e−itω1L 2 ω1L (λ β + iω2L )(λ α + i(ω1L + ω2L )) −

2 ω1L (ω2L

−1 + eitω1L 1 − e−itω2L − 2 − iλ α )(iλ β + ω1L − ω2L ) ω2L (ω2L − iλ β )(−iλ α + ω1L + ω2L )

e λ β t −1 + eitω2L −1 + eitω2L − + 2 λ β ω2L (λ α + iω1L ) λ β ω2L (λ β + iω2L )(−iλ α + iλ β + ω1L ) This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Fri, 08 May 2015 07:26:14

174104-35


J. Chem. Phys. 142, 174104 (2015)

−1 + eitω1L e t(λ β −i(ω1L −ω2L )) + ω1L (ω2L − iλ β )(−iλ α + iλ β + ω1L )(iλ β + ω1L − ω2L ) −1 + eitω1L e t(λ α +iω2L ) + ω1L (ω2L − iλ α )(−iλ α + iλ β + ω1L )(−iλ α + ω1L + ω2L ) i −1 + eitω2L e t(λ α +iω1L ) − . ω2L (λ α + iω1L )(λ α − λ β + iω1L )(λ α + i(ω1L + ω2L ))

(A70)

In the long time limit, the terms involving the three-point correlation functions do not contribute to the generalized Mandel parameter because they have no terms proportional to t. The terms including the two-time correlation functions are expressed below.  A b (t,ω1,ω2) = Aα;0;β η α, β (t,ω1,ω2) , (A71) α, β

where  t t4 t3 t2 η α, β (t,ω1,ω2) = t 4=0 t 3=0 t 2=0 t 1=0

eiω1L (t1−t4)eiω2L (t2−t3) e λ β (t4−t3)e λ α (t2−t1) * iω1L (t1−t3) iω2L (t2−t4)+ dt 1dt 2dt 3dt 4. e ,+e

(A72)

The solution of the integrals is η α, β (t,ω1,ω2) =−

( ) 2 2 + 8ω1Lω2L + 3ω2L iλ α −2iλ 3β − 5λ 2β (ω1L + ω2L ) + iλ β 3ω1L (λ α − iω1L )2(ω1L + iλ β )2(λ β − iω2L )2(ω1L + ω2L )2

) ( 2 2 i −2λ 3β (2ω1L + ω2L ) + iλ 2β (ω1L + ω2L )(8ω1L + 3ω2L ) iλ α ω1L + ω1Lω2L + ω2L − − (λ α − iω1L )2(ω1L + iλ β )2(λ β − iω2L )2(ω1L + ω2L ) (λ α − iω1L )2(ω1L + iλ β )2(λ β − iω2L )2(ω1L + ω2L )2 3 2 2 3 2 2 iλ β 4ω1L + 13ω1L ω2L + 8ω1Lω2L + ω2L ω1L ω1L + 2ω1Lω2L + 2ω2L − − (λ α − iω1L )2(ω1L + iλ β )2(λ β − iω2L )2(ω1L + ω2L )2 (λ α − iω1L )2(ω1L + iλ β )2(λ β − iω2L )2(ω1L + ω2L ) +

e t(λ α −iω1L )(2λ α − 2λ β − i(ω1L − ω2L )) i(ω1L + ω2L − 2iλ β )e−it(ω1L +ω2L ) + (λ α − λ β )(λ α − iω1L )2(λ α + iω2L )(λ α − λ β − i(ω1L − ω2L )) (ω2L − iλ α )(λ β + iω1L )(ω2L − iλ β )(ω1L + ω2L )2

+

t(2iλ β + ω1L + ω2L ) (λ α − iω1L )(ω1L + iλ β )(ω2L + iλ β )(ω1L + ω2L )

+

e λ β t−itω1L e λ β t−itω2L − . 2 (λ α − λ β )(λ β − iω1L ) (λ β + iω2L ) (λ β + iω1L )(λ β − iω2L )2(λ α − λ β − i(ω1L − ω2L ))

(A73)

The leading order term in the long time limit is η α, β (t,ω1,ω2) ∼

t(2iλ β + ω1L + ω2L ) . (λ α − iω1L )(ω1L + iλ β )(ω2L + iλ β )(ω1L + ω2L )

K1b (t,ω1,ω2) =



Kα;0;β k2α, β (t,ω1,ω2) ,

(A74)

(A75)

α, β

where  t t4 t3 t2 k2α, β (t,ω1,ω2) =

e λ β (t4−t3)e λ α (t2−t1)eiω1L (t1−t4)eiω2L (t3−t2)dt 1dt 2dt 3dt 4.

(A76)

t 4=0 t 3=0 t 2=0 t 1=0


174104-36


J. Chem. Phys. 142, 174104 (2015)

The solution of the integrals is k2α, β (t,ω1,ω2) =

t −1 + e t(λ α −iω1L ) + 2 (ω ) (λ − ω − iω1L ) ω1L + iλ β λ α − λ β (λ α − iω1L ) (λ α − iω2L ) 1L 2L α i 1 − e−it(ω1L −ω2L ) −1 + e t ( λ β −iω1L ) + . 2 + (ω1l − ω2l )2 (λ α − iω2l ) ω2l + iλ β λ β − λ α λ β − iω1L λ β − iω2L

(A77)

The leading order term in the long time limit is k2α, β (t,ω1,ω2) ∼

t . (ω1L − ω2L ) (λ α − iω1L ) ω1L + iλ β

G1b (t,ω1,ω2) =



(A78)

Gα;0;β g2α, β (t,ω1,ω2) ,

(A79)

α, β

where  t t4 t3 t2 g2α, β (t,ω1,ω2) =

e λ β (t4−t3)e λ α (t2−t1)eiω1L (t2−t4)eiω2L (t3−t1)dt 1dt 2dt 3dt 4.

(A80)

t 4=0 t 3=0 t 2=0 t 1=0

The solution of the integrals is g2α, β (t,ω1,ω2) =

λ β − iω1 +

2

t −1 + e t ( λ β −iω1) + (ω ) (λ − ω + iω2) ω1 + iλ β 1 2 α λ β − iω2 −λ α + λ β − i (ω1 + ω2)

i 1 − e−it(ω1−ω2) (ω1 − ω2)2 (λ α + iω1) ω2 + iλ β

+

−1 + e t(λ α +iω2) (λ α + iω1) (λ α + iω2)2 λ α − λ β + i (ω1 + ω2)

.

(A81)

The leading order term in the long time limit is g2α, β (t,ω1,ω2) ∼

t . (ω1 − ω2) (λ α + iω2) ω1 + iλ β

(A82)

The next term to be written is ψ (t,ω1,ω2), ψ (t,ω1,ω2) =



Aα;0;0ψ1α (t,ω1,ω2) ,

(A83)

α

where  t  t  t t2 ψ1α (t,ω1,ω2) = 0

e λ α (t2−t1)eiω1L (t1−v)eiω2L (t2−u)dt 1dt 2dv du.

(A84)

0 t 2=0 t 1=0

The solution of the integrals is 1 − e−itω1L 1 − e−itω2L ψ1α (t,ω1,ω2) = ω1Lω2L (λ α − iω1L )(λ α + iω2L )(ω1L + ω2L ) ( ( ) ) × −iλ α −1 + eit(ω1L +ω2L ) − (ω1L + ω2L )e t(λ α +iω2L ) + ω2L eit(ω1L +ω2L ) + ω1L , Ψ (t,ω1,ω2) =



Kα;0;0ψ1α (t,ω2, −ω1) .

(A85) (A86)

α

The last function to be expressed is ϕ (t,ω1,ω2),  ϕ (t,ω1,ω2) = Kα;0;β ϕ1α, β (t,ω1,ω2) + G β;0;α ϕ1α, β (t,ω2, −ω1) ,

(A87)

α, β

where  t t4 t3 t2  (iω1L −λ β )(t2−t4) (iω2L −λ α )(t1−t3)  e  e  ϕ1α, β (t,ω1,ω2) =  (iω2L −λ β )(t1−t4) −(iω1L −λ α )(t3−t2) dt 1dt 2dt 3dt 4. e + e  t =0 t =0 t =0 t =0

(A88)

4 3 This article is copyrighted as indicated in the article. Reuse of 2AIP 1content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Fri, 08 May 2015 07:26:14

174104-37


J. Chem. Phys. 142, 174104 (2015)

The solution of the integrals is i −1 + e t(λ β −iω1L ) −1 + e t(λ β −iω1L ) ϕ1α, β (t,ω1,ω2) = + (λ α + iω2L )(λ β − iω1L )2(ω1L − ω2L ) (λ α + iω1L )(λ β − iω1L )(ω2L + iλ β )(ω1L − ω2L ) −

−1 + e t(λ β −iω1L ) (λ α + iω1L )(λ α + iω2L )(λ β − iω1L )(λ α − λ β + i(ω1L + ω2L ))

+

t −1 + e t(λ α +iω2L ) + (λ α + iω1L )(λ α + iω2L )2(λ α − λ β + i(ω1L + ω2L )) (λ α + iω2L )(ω1L + iλ β )(ω1L − ω2L )

i 1 − e−it(ω1L −ω2L ) . + (λ α + iω1L )(ω2L + iλ β )(ω1L − ω2L )2

(A89)

The leading order term in the long time limit is t . (λ α + iω2L )(ω1L + iλ β )(ω1L − ω2L ) The function M (t,ω1,ω2) can now be expressed as the sum of simple terms as  ′   Aα;β;γ a1α;β;γ (t,ω1,ω2) + Kα;β;γ k1α;β;γ (t,ω1,ω2)   α;β;γ ′ ′   + Gα;β;γ g1α;β;γ (t,ω1,ω2) + Aα;0;β η α;β (t,ω1,ω2)      α;β;γ α;β  + ′ A (t,ω ) (t,ω ) (t,ω ) s1 ,ω + K p1 ,ω + G b1 ,ω α;β;0 α;β 1 2 α;β;0 α;β 1 2 α;β;0 α;β 1 2    α;β ′  . M (t,ω1,ω2) = 2 Re  + Kα;0;β k2α;β (t,ω1,ω2) − ϕ1α;β (t,ω1,ω2)     α;β ′   (t,ω ) (t,ω ) + G g2 ,ω − ϕ1 , −ω α;0;β α;β 1 2 β;α 2 1   α;β    ′   + Kα;0;0ψ1α (t,ω2, −ω1) + Aα;0;0ψ1α (t,ω1,ω2)   ϕ1α, β (t,ω1,ω2) ∼

(A90)

(A91)

α

The values of the moments and the correlation functions presented in the manuscript were obtained by numerically calculating the integrals over the frequencies ω1 and ω2.

APPENDIX B: DETAILS OF THE SECOND MOMENT CALCULATION IN SEC. V

The goal of this appendix is to provide the details of the calculation of the second multivariate moment of the number of photons, including the effect of the response function of the finite width detector. Starting with the generating function of Eq. (73), the second multivariate moment is given by t t

d

d d 2 4 −2∆t 2 N(ω,∆) (t) N(ω′,∆) (t) − δω,ω′ N(ω,∆) (t) = Γ0 ∆ dt 2 e−2∆t1dt 1 e 0

0

t2

t2 ×

du −∞

t1 dw

−∞

t1 dv

−∞

d y TN D + (u) D + (v) D − (w) D − ( y) eiω1L (u−w)

−∞

× eiω2L (v−y)e∆(u+v+w+y).

(B1)

We follow a similar procedure to the one presented in Appendix A and express the second multivariate moment as ⟨N (ω1, ∆,t) N (ω2, ∆,t)⟩ − δω1,ω2⟨N (ω1, ∆,t)⟩ =

∆6Γ02 [Υ (ω1,ω2, ∆,t) + Υ (ω2,ω1, ∆,t)] , π2

(B2)

where t Υ (t,ω1,ω2, ∆) =

t dt 2e

0

−2∆t 2

t1 dt 1e

0

−2∆t 1 −∞

t1 dv

t2 dw

−∞

−∞

t2 du

dyΘ (u − v) Θ (w − y)

−∞

× D + (v) D + (u) D − (w) D − ( y) eiω1L (u−y) eiω2L (v−w)e∆(u+v+w+y)


174104-38


J. Chem. Phys. 142, 174104 (2015)

+

t1

t

t

dt 1e−2∆t1

dt 2e−2∆t2

dv

−∞

0

0

t2

t2 dw

−∞

t1 du

−∞

dyΘ (u − v) Θ (w − y)

−∞

× D + (v) D + (u) D − (w) D − ( y) eiω1L (u−w) eiω2L (v−y)e∆(u+v+w+y).

(B3)

As was the case for the implicit treatment of the detector, in order to actually calculate the correlation functions, we express them in a fully time-ordered manner. Therefore, we rewrite the equation above as t

t Υ (t,ω1,ω2, ∆) = 2 Re

dt 2e

−2∆t 2

dt 1e

−2∆t 1

−∞

0

0

t1 t2 t2 t1

[ dτ4 dτ3 dτ2 dτ1 D + (τ3) D+ (τ4) D − (τ2) D − (τ1) −∞

−∞

−∞

× eiω1L (τ4−τ1)eiω2L (τ3−τ2)e∆(τ4+τ3+τ2+τ1)Θ (τ4 − τ3) Θ (τ3 − τ2) Θ (τ2 − τ1) t2 +

t1 dτ4

−∞

t2 dτ3

−∞

t1 dτ2

−∞

dτ1 D + (τ3) D + (τ4) D − (τ2) D − (τ1)

−∞


t2 dτ4

−∞

t2 dτ3

−∞

t1 dτ2

−∞

dτ1 D + (τ1) D + (τ3) D − (τ4) D− (τ2)

−∞


t2 dτ4

−∞

t1 dτ3

−∞

t1 dτ2

−∞

dτ1 D+ (τ1) D + (τ3) D − (τ4) D− (τ2)

−∞


t2 dτ4

−∞

t1 dτ3

−∞

t1 dτ2

−∞

dτ1 D+ (τ1) D + (τ4) D − (τ3) D− (τ2)

−∞


t1 dτ4

−∞

t2 dτ3

−∞

t1 dτ2

−∞

dτ1 D+ (τ1) D + (τ4) D − (τ3) D− (τ2)

−∞

× eiω1L (τ4−τ2)eiω2L (τ1−τ3)e∆(τ4+τ2+τ3+τ1)Θ (τ4 − τ3) Θ (τ3 − τ2) Θ (τ2 − τ1)].

(B4)

In order to calculate the numerator, we also need the multiplication of the first moments at two different frequencies. This can be written as ⟨d N (ω1, ∆,t)⟩ ⟨d N (ω2, ∆,t)⟩ = 2 Re [Π (ω1,ω2, ∆,t) + Π (ω2,ω1, ∆,t)] .

(B5)

Here, we introduced t Π (ω1,ω2, ∆,t) =

t dt 1e

−2∆t 1

t1 dt 2e

−2∆t 2

dτ4

−∞

0

0

t1

t2 dτ3

−∞

t2 dτ2

−∞

dτ1 D + (τ4) D − (τ3)

−∞

× D + (τ2) D − (τ1) eiω1L (τ2−τ1)eiω2L (τ4−τ3)Θ (τ4 − τ3) Θ (τ3 − τ2) Θ (τ2 − τ1) e∆(τ1+τ2+τ3+τ4)

t +

t dt 1e

0

−2∆t 1

t1 dt 2e

−2∆t 2 −∞

0

t1 dτ4 −∞

t2 dτ1 −∞

t2 dτ3

dτ2 D + (τ4) D − (τ1)

−∞



174104-39


J. Chem. Phys. 142, 174104 (2015)

t1

t

t +

dt 1e

−2∆t 1

dt 2e

−2∆t 2

dτ3

−∞

0

0

t1

t2 dτ1

−∞

t2 dτ4

−∞

dτ2 D+ (τ3) D − (τ1)

−∞


t1

t

t +

dt 1e

−2∆t 1

dt 2e

−2∆t 2

dτ4

−∞

0

0

t1

t2 dτ3

−∞

t2 dτ1

−∞

dτ2 D+ (τ4) D − (τ3)

−∞


t1

t

t +

dt 1e

−2∆t 1

dt 2e

−2∆t 2

dτ4

−∞

0

0

t1

t2 dτ1

−∞

t2 dτ2

−∞

dτ3 D + (τ4) D − (τ1)

−∞


t

t +

dt 1e

t1

−2∆t 1

dt 2e

−2∆t 2

dτ3

−∞

0

0

t1

t2 dτ1

−∞

t2 dτ2

−∞

dτ4 D + (τ3) D − (τ1)

−∞

× D + (τ2) D − (τ4) eiω1L (τ2−τ4)eiω2L (τ3−τ1)Θ (τ4 − τ3) Θ (τ3 − τ2) Θ (τ2 − τ1) e∆(τ1+τ2+τ3+τ4).

(B6)

The correlation functions are then calculated according to the procedure described in Sec. IV. In a way similar to what was done for the case of no implicit detector effect, we express the integrals over the multi-point correlation function by a sum of simple terms. The numerator of the correlation function (Eq. (32)) is written as Γ02∆4 d M (t,ω1,ω2, ∆) + d M (t,ω2,ω1, ∆) , (B7) where    Aα;β;γ d k α;β;γ (t, −ω1,ω2, ∆) + d k α;β;γ (t,ω1,ω2, ∆)   α;β;γ    d d  (t, (t,ω + K a −ω ,ω , ∆) + G k , −ω , ∆) α;β;γ α;β;γ 2 1 α;β;γ α;β;γ 1 2   α;β;γ    d d   . (t, (t, + G g −ω ,ω , ∆) + K g −ω , −ω , ∆) d α;β;γ α;β;γ 1 2 α;β;γ α;β;γ 1 2 M (t,ω1,ω2, ∆) = 2 Re  (B8)  α;β,0;γ   d  d   − Kα;0;β ϑα;β (t,ω1,ω2, ∆) + ϱα;β (t,ω1,ω2, ∆)   α;β  d   d (t,ω (t,ω − G ϑ , −ω , ∆) + ϱ , −ω , ∆) α;0;β α;β 1 2 α;β 1 2   α;β The additional quantities we defined are detailed below. t d

aα;β;γ (t,ω1,ω2, ∆) =

t1

t dt 2e

−2∆t 2

0

dt 1e 0

−2∆t 1 −∞

t2 dτ4 −∞

t2 dτ3 −∞

t1 dτ2

dτ1Θ (τ4 − τ3)

−∞

× Θ (τ3 − τ2) Θ (τ2 − τ1) e λ α (τ2−τ1)e λ β (τ3−τ2)e λγ (τ4−τ3)e−iω1L (τ4−τ1)e−iω2L (τ3−τ2)e∆(τ4+τ3+τ2+τ1).

(B9)

The solution of the integrals is d

aα;β;γ (t,ω1,ω2, ∆) =

8∆3(−λ

e−2∆t (−λ γ + ∆ + iω1L ) α + ∆ + iω1L )(ω1L + i(λ γ − 3∆))(λ γ + ∆ − iω1L )(i(λ β − 2∆) + ω1L + ω2L )

−

it 4∆2(−λ α + ∆ + iω1L )(−λ γ + ∆ + iω1L )(i(λ β − 2∆) + ω1L + ω2L )

−

e t(λγ −∆−iω1L ) 1 (−λ α + ∆ + iω1L )(ω1L + i(λ γ − 3∆))(−λ γ + ∆ + iω1L )2(λ γ + ∆ − iω1L ) (i(λ β − 2∆) + ω1L + ω2L )

+

7i∆2 − 4∆(ω1L + iλ γ ) − i(ω1L + iλ γ )2 . 8∆3(ω1L + i(λ α − ∆))(ω1L + i(λ γ − 3∆))(ω1L + i(λ γ − ∆))2(i(λ β − 2∆) + ω1L + ω2L )

(B10)


174104-40


J. Chem. Phys. 142, 174104 (2015)

The leading order term in the long time limit is d

aα;β;γ (t,ω1,ω2, ∆) ∼ − t

d

4∆2(−λ

t

k α;β;γ (t,ω1,ω2, ∆) =

dt 2e−2∆t2 0

it , α + ∆ + iω1L )(−λ γ + ∆ + iω1L )(i(λ β − 2∆) + ω1L + ω2L )

t2 dt 1e−2∆t1

t1 dτ4

−∞

0

t2 dτ3

−∞

t1 dτ1Θ (τ4 − τ3)

dτ2

−∞

(B11)

−∞

× Θ (τ3 − τ2) Θ (τ2 − τ1) e λ α (τ2−τ1)e λ β (τ3−τ2)e λγ (τ4−τ3)eiω1L (τ4−τ2)eiω2L (τ3−τ1)e∆(τ4+τ3+τ2+τ1).

(B12)


k α;β;γ (t,ω1,ω2, ∆) = −

ie−2∆t (λ γ − ∆ + iω1L ) 8∆3(λ γ − 3∆ + iω1L )(λ γ + ∆ + iω1L )(−λ α + ∆ − iω2L )(−i(λ β − 2∆) + ω1L + ω2L ) t

−

4∆2(−λ

γ

+ ∆ − iω1L )(ω2L − i(λ α − ∆))(−i(λ β − 2∆) + ω1L + ω2L )

1 ie t(λγ −∆+iω1L ) (λ γ − 3∆ + iω1L )(λ γ − ∆ + iω1L )2(λ γ + ∆ + iω1L )(−λ α + ∆ − iω2L ) −i(λ β − 2∆) + ω1L + ω2L

−

+

8∆3(ω

−7∆2 + 4∆(λ γ + iω1L ) + (ω1L − iλ γ )2 . 2 1L − i(λ γ − 3∆))(ω1L − i(λ γ − ∆)) (−λ α + ∆ − iω2L )(−i(λ β − 2∆) + ω1L + ω2L )

(B13)

The leading order in the long time limit is d

k α;β;γ (t,ω1,ω2, ∆) ∼ − t

d

t , 4∆2(−λ γ + ∆ − iω1L )(ω2L − i(λ α − ∆))(−i(λ β − 2∆) + ω1L + ω2L )

t

gα;β;γ (t,ω1,ω2, ∆) =

dt 2e

−2∆t 2

0

t2 dt 1e

−2∆t 1

t2 dτ4

−∞

0

−∞

t1 dτ3 −∞

(B14)

t1 dτ2

dτ1

−∞

× Θ (τ4 − τ3) Θ (τ3 − τ2) Θ (τ2 − τ1) e(λ α −iω2L )(τ2−τ1)e λ β (τ3−τ2)e( λγ −iω1L )(τ4−τ3)e∆(τ4+τ3+τ2+τ1).

(B15)


gα;β;γ (t,ω1,ω2, ∆) =

eλ βt λ 2β

(

λ 2β

−

4∆2

)

(λ α − ∆ − iω2L )(λ β − λ γ + ∆ + iω1L )

+

t 4λ β ∆2(−λ γ + ∆ + iω1L )(λ α − ∆ − iω2L )

e−2∆t λ β (λ γ − ∆) + 4∆2 − iλ β ω1L + 3 8∆ (2∆ − λ β )(λ β + 2∆)(λ γ + ∆ − iω1L )(−λ γ + 3∆ + iω1L )(−λ α + ∆ + iω2L ) +

+

−

+

1 e t(λγ −∆)−itω1L (ω1L + i(λ γ − 3∆))(ω1L + i(λ γ − ∆))2(ω1L + i(λ γ + ∆))(−λ α + ∆ + iω2L ) λ β − λ γ + ∆ + iω1L ( ) 2 12∆4 + ∆2 7λ 2β + 12λ β λ γ + 4λ γ2 − 4ω1L − 4iω1L (3λ β + 2λ γ ) − λ 2β (ω1L + iλ γ )2 8λ 2β ∆3(λ β − 2∆)(−λ γ + ∆ + iω1L )2(−λ γ + 3∆ + iω1L )(−λ α + ∆ + iω2L ) (9λ β + 8λ γ − 8iω1L ) 4λ 2β (λ β

4λ β

− 2∆)(−λ γ + ∆ + iω1L )2(−λ γ + 3∆ + iω1L )(−λ α + ∆ + iω2L )

∆2(λ

(ω1L + iλ γ )(ω1L + i(2λ β + λ γ )) . 2 β − 2∆)(−λ γ + ∆ + iω1L ) (−λ γ + 3∆ + iω1L )(−λ α + ∆ + iω2L )

(B16)


174104-41


J. Chem. Phys. 142, 174104 (2015)


ϱα;β (t,ω1,ω2, ∆) =

t1

t

t d

gα;β;γ (t,ω1,ω2, ∆) ∼

dt 2e

−2∆t 2

dt 1e

−2∆t 1

t1 dτ4

−∞

0

0

t , 4λ β ∆2(−λ γ + ∆ + iω1L )(λ α − ∆ − iω2L ) t2 dτ1

−∞

t2 dτ3

−∞

(B17)

dτ2

−∞

× Θ (τ4 − τ3) Θ (τ3 − τ2) Θ (τ2 − τ1) e λ β (τ4−τ1)e λ α (τ3−τ2)e−iω1L (τ3−τ2)e−iω2L (τ4−τ1)e∆(τ4+τ3+τ2+τ1).

(B18)


ϱα;β (t,ω1,ω2, ∆) =

8∆3(λ

ie−2∆t β + ∆ − iω2L )(−λ β + 3∆ + iω2L )(i(λ α + λ β − 2∆) + ω1L + ω2L )

+

it 4∆2(ω2L + i(λ β − ∆))2(i(λ α + λ β − 2∆) + ω1L + ω2L )

+

e t(λ β −∆−iω2L ) (ω2L + i(λ β − 3∆))(ω2L + i(λ β − ∆))3(ω2L + i(λ β + ∆))(i(λ α + λ β − 2∆) + ω1L + ω2L )

+

7i∆2 − 4∆(ω2L + iλ β ) − i(ω2L + iλ β )2 . 8∆3(ω2L + i(λ β − 3∆))(ω2L + i(λ β − ∆))3(i(λ α + λ β − 2∆) + ω1L + ω2L )

(B19)


ϱα;β (t,ω1,ω2, ∆) ∼ t

d

ϑα;β (t,ω1,ω2, ∆) =

it , 4∆2(ω2L + i(λ β − ∆))2(i(λ α + λ β − 2∆) + ω1L + ω2L )

t dt 2e

−2∆t 2

0

t1 dt 1e

0

−2∆t 1

t1 dτ4

−∞

−∞

t2 dτ1 −∞

(B20)

t2 dτ3

dτ2

−∞

× Θ (τ4 − τ3) Θ (τ3 − τ2) Θ (τ2 − τ1) e λ β (τ4−τ1)e λ α (τ3−τ2)e−iω1L (τ4−τ2)e−iω2L (τ3−τ1)e∆(τ4+τ3+τ2+τ1).

(B21)


ϑα;β (t,ω1,ω2, ∆) =− +

e−2∆t (−λ β + ∆ + iω1L ) 8∆3(λ β − 3∆ − iω1L )(λ β + ∆ − iω1L )(ω2L + i(λ β − ∆))(i(λ α + λ β − 2∆) + ω1L + ω2L ) 4∆2(ω

t 1L + i(λ β − ∆))(λ β − ∆ − iω2L )(i(λ α + λ β − 2∆) + ω1L + ω2L )

−

ie t(λ β −∆−iω1L ) 1 (ω1L + i(λ β − 3∆))(ω1L + i(λ β − ∆))2(ω1L + i(λ β + ∆)) (λ β − ∆ − iω2L )(i(λ α + λ β − 2∆) + ω1L + ω2L )

+

−7∆2 + 4∆(λ β − iω1L ) + (ω1L + iλ β )2 1 . 3 2 8∆ (λ β − 3∆ − iω1L )(−λ β + ∆ + iω1L ) (ω2L + i(λ β − ∆)) i(λ α + λ β − 2∆) + ω1L + ω2L

(B22)


ϑα;β (t,ω1,ω2, ∆) ∼

4∆2(ω

t . 1L + i(λ β − ∆))(λ β − ∆ − iω2L )(i(λ α + λ β − 2∆) + ω1L + ω2L )

(B23)

This appendix, although very formal and cumbersome, details the exact procedure required for the actual evaluation of numerical values for the correlation function of Eq. (32). This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Fri, 08 May 2015 07:26:14

174104-42 1W.


E. Moerner and L. Kador, Phys. Rev. Lett. 62, 2535 (1989). Orrit and J. Bernard, Phys. Rev. Lett. 65, 2716 (1990). 3W. E. Moerner and M. Orrit, Science 283, 1670 (1999). 4T. Plakhotnik, E. A. Donley, and U. P. Wild, Annu. Rev. Phys. Chem. 48, 181 (1997). 5H. P. Lu, L. Xun, and X. S. Xie, Science 282, 1877 (1998). 6X. S. Xie, J. Chem. Phys. 117, 11024 (2002). 7X. Zhuang, L. E. Bartley, H. P. Babcock, R. Russell, T. Ha, D. Herschlag, and S. Chu, Science 288, 2048 (2000). 8S. Weiss, Science 283, 1676 (1999). 9Y. Jung, E. Barkai, and R. J. Silbey, J. Chem. Phys. 117, 10980 (2002). 10M. Orrit, J. Chem. Phys. 117, 10938 (2002). 11M. Bohmer and J. Enderlein, ChemPhysChem 4, 792 (2003). 12M. Lippitz, F. Kulzer, and M. Orrit, ChemPhysChem 6, 770 (2005). 13E. Geva and J. L. Skinner, J. Phys. Chem. B 101, 8920 (1997). 14J. Wang and P. Wolynes, Phys. Rev. Lett. 74, 4317 (1995). 15J. Wang and P. Wolynes, J. Chem. Phys. 110, 4812 (1999). 16J. J. Portman and P. G. Wolynes, J. Phys. Chem. A 103, 10602 (1999). 17N. Agmon, J. Phys. Chem. B 104, 7830 (2000). 18J. Cao, Chem. Phys. Lett. 327, 38 (2000). 19S. Yang and J. Cao, J. Phys. Chem. B 105, 6536 (2001). 20V. Barsegov, V. Chernyak, and S. Mukamel, J. Chem. Phys. 116, 4240 (2002). 21V. Barsegov and S. Mukamel, J. Phys. Chem. A 108, 15 (2004). 22E. Barkai, Y. Jung, and R. Silbey, Phys. Rev. Lett. 87, 207403 (2001). 23Y. Jung, E. Barkai, and R. Silbey, Chem. Phys. 284, 181 (2002). 24F. L. H. Brown, Phys. Rev. Lett. 90, 028302 (2003). 25L. Fleury, J.-M. Segura, G. Zumofen, B. Hecht, and U. P. Wild, Phys. Rev. Lett. 84, 1148 (2000). 26I. V. Gopich and A. Szabo, J. Chem. Phys. 118, 454 (2003). 27I. V. Gopich and A. Szabo, J. Chem. Phys. 122, 014707 (2005). 28O. Flomenbom, J. Klafter, and A. Szabo, Biophys. J. 88, 3780 (2005). 29D. E. Makarov and H. Metiu, J. Chem. Phys. 115, 5989 (2001). 30G. K. Schenter, H. P. Lu, and X. S. Xie, J. Phys. Chem. A 103, 10477 (1999). 31A. M. Berezhkovskii, A. Szabo, and G. H. Weiss, J. Phys. Chem. B 104, 3776 (2000). 32S. Mukamel, Phys. Rev. A 68, 063821 (2003). 33Y. Zheng and F. L. H. Brown, Phys. Rev. Lett. 90, 238305 (2003). 34Y. Zheng and F. L. H. Brown, J. Chem. Phys. 119, 11814 (2003). 35Y. Zheng and F. L. H. Brown, J. Chem. Phys. 121, 3238 (2004). 36Y. Zheng and F. L. H. Brown, J. Chem. Phys. 121, 7914 (2004). 37G. Bel, Y. Zheng, and F. L. H. Brown, J. Phys. Chem. B 110, 19066 (2006). 38M. B. Plenio and P. L. Knight, Rev. Mod. Phys. 70, 101 (1998). 39I. S. Osad’ko, J. Exp. Theor. Phys. 86, 875 (1998). 40R. Verberk and M. Orrit, J. Chem. Phys. 119, 2214 (2003). 41Y. Jung, E. Barkai, and R. J. Silbey, Adv. Chem. Phys. 123, 199 (2002). 42Y. He and E. Barkai, Phys. Rev. Lett. 93, 068302 (2004). 43Y. He and E. Barkai, J. Chem. Phys. 122, 184703 (2005). 44A. M. Berezhkovski, A. Szabo, and G. H. Weiss, J. Chem. Phys. 110, 9145 (1999). 45I. V. Gopich and A. Szabo, J. Phys. Chem. B 114, 15221 (2010). 46W. Dong-Sheng and Z. Yu-Jun, Chin. Phys. B 19, 083202 (2010). 47Y. Peng, D. Wang, Y. Zheng, and S. Xie, Physica E 42, 2242 (2010). 48E. A. Bloemsma and J. Knoester, J. Chem. Phys. 136, 224507 (2012). 49Y. Zheng and F. L. H. Brown, J. Chem. Phys. 139, 164120 (2013). 50B. Han, Y. Pan, C. Han, F. Hu, and Y. Zheng, J. Phys. B: At., Mol. Opt. Phys. 47, 025502 (2014). 51T. Ha, T. Enderle, D. F. Ogletree, D. S. Chemla, P. R. Selvin, and S. Weiss, 2M.

J. Chem. Phys. 142, 174104 (2015) Proc. Natl. Acad. Sci. U. S. A. 93, 6264 (1996). 52E. Moreau, I. Robert, L. Manin, V. Thierry-Mieg, J. M. Gerard, and I. Abram,

Phys. Rev. Lett. 87, 183601 (2001). K. Luong, C. C. Gradinaru, D. W. Chandler, and C. C. Hayden, J. Phys. Chem. B 109, 15691 (2005). 54B. R. Mollow, Phys. Rev. 188, 1969 (1969). 55F. Schuda, C. R. Stroud, Jr., and M. Hercher, J. Phys. B: At. Mol. Phys. 7, L198 (1974). 56R. E. Grove, F. Y. Wu, and S. Ezekiel, Phys. Rev. A 15, 227 (1977). 57G. Wrigge, I. Gerhardt, J. Hwang, G. Zumofen, and V. Sandoghdar, Nat. Phys. 4, 60 (2008). 58A. Aspect, G. Roger, S. Reynaud, J. Dalibard, and C. Cohen-Tannoudji, Phys. Rev. Lett. 45, 617 (1980). 59O. Astafiev, A. M. Zagoskin, A. A. Abdumalikov Jr., Yu. A. Pashkin, T. Yamamoto, K. Inomata, Y. Nakamura, and J. S. Tsai, Science 327, 840 (2010). 60C. Matthiesen, A. N. Vamivakas, and M. Atatüre, Phys. Rev. Lett. 108, 093602 (2012). 61G. Nienhuis, Phys. Rev. A 47, 510 (1993). 62C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg, Atom-Photon Interactions (Wiley-Interscience, New York, 1992). 63K. Joosten and G. Nienhuis, J. Opt. B: Quantum Semiclassical Opt. 2, 158 (2000). 64M. Alexanian and S. K. Bose, Phys. Rev. A 74, 063418 (2006). 65F. L. H. Brown, Acc. Chem. Res. 39, 363 (2006). 66G. Bel and F. L. H. Brown, Phys. Rev. Lett. 102, 018303 (2009). 67R. M. Wilcox, J. Math. Phys. 8, 962 (1967). 68M. O. Scully and M. S. Zubairy, Quantum Optics (Cambridge University Press, Cambridge, UK, 1997). 69R. Loudon, The Quantum Theory of Light, 3rd ed. (Oxford, 2000). 70B. R. Mollow, Phys. Rev. A 12, 1919 (1975). 71B. R. Mollow, J. Phys. A: Math. Gen. 8, L130 (1975). 72L. Mandel, Opt. Lett. 4, 205 (1979). 73C. W. Gardiner, Quantum Noise (Springer-Verlag, Berlin, 1991). 74K. Blum, Density Matrix Theory and Applications, 2nd ed. (Plenum Press, New York, 1981). 75C. H. Page, J. Appl. Phys. 23, 103 (1952). 76D. G. Lampard, J. Appl. Phys. 25, 802 (1954). 77J. H. Eberly and K. Wodkiewicz, J. Opt. Soc. Am. 67, 1252 (1977). 78R. J. Glauber, Quantum Optics and Electrodynamics, Les Houches 1964, edited by C. de Witt, A. Blandin, and C. Cohen-Tannoudji (New York: Gordon and Breach Science Publishers, Inc., 1965), p. 63. 79P. L. Kelley and W. H. Kleiner, Phys. Rev. 136, A316 (1964). 80W. Vogel and D.-G. Welsch, Quantum Optics, 3rd ed. (Wiley-VCH, Weinheim, Germany, 2006). 81L. Knoll, W. Vogel, and D. G. Welsch, Phys. Rev. A 36, 3803 (1987). 82P. W. Anderson, B. I. Halperin, and C. M. Varma, Philos. Mag. 25, 1 (1972). 83W. A. Phillips, J. Low Temp. Phys. 7, 351 (1972). 84R. Silbey, Relaxation Processes in Molecular Excited States (Kluwer Academic Publishers, 1989), pp. 243–276. 85F. L. H. Brown and R. J. Silbey, J. Chem. Phys. 108, 7434 (1998). 86R. Kubo, Fluctuation, Relaxation and Resonance in Magnetic Systems, edited by D. TerHaar (Oliver and Boyd, Edinburgh, 1962). 87R. Kubo, M. Toda, and N. Hashitsume, Statistical Physics II, 2nd ed. (Springer-Verlag, Berlin, 1991). 88E. Geva and J. L. Skinner, Chem. Phys. Lett. 288, 225 (1998). 89S. Mukamel, Principles of Nonlinear Optical Spectroscopy (Oxford, 1995). 53A.


Journal of Chemical Physics is copyrighted by AIP Publishing LLC (AIP). Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. For more information, see http://publishing.aip.org/authors/rights-and-permissions.

Single Molecule Spectroscopy of Monomeric LHCII: Experiment and Theory.

Fast single-molecule FRET spectroscopy: theory and experiment.

Single molecule laser spectroscopy.

Single-molecule optical spectroscopy.

Stimulated emission properties of fluorophores by CW-STED single molecule spectroscopy.

Bloch Surface Wave-Coupled Emission from Quantum Dots by Ensemble and Single Molecule Spectroscopy.

Testing Landscape Theory for Biomolecular Processes with Single Molecule Fluorescence Spectroscopy.

Single molecule microscopy and spectroscopy: concluding remarks.

Ultra-nonlocality in density functional theory for photo-emission spectroscopy.

Single-molecule spectroscopy of unfolded proteins and chaperonin action.

Single-molecule spectroscopy reveals chaperone-mediated expansion of substrate protein.

Dynamic single-molecule force spectroscopy of rhodopsin in native membranes.

Unraveling the gymnastics of reverse transcription through single molecule spectroscopy.

Fluorescence-excitation and Emission Spectroscopy on Single FMO Complexes.

Investigating single molecule adhesion by atomic force spectroscopy.

Correction to molecular approaches to chromatography using single molecule spectroscopy.

Single-molecule surface-enhanced Raman spectroscopy with nanowatt excitation.

Radio frequency scanning tunneling spectroscopy for single-molecule spin resonance.

Silicon photon-counting avalanche diodes for single-molecule fluorescence spectroscopy.

Directly measuring single-molecule heterogeneity using force spectroscopy.

Second harmonic generation correlation spectroscopy for single molecule experiments.

Single-molecule spectroscopy and imaging over the decades.

Improved single molecule force spectroscopy using micromachined cantilevers.

Mechanisms of small molecule-DNA interactions probed by single-molecule force spectroscopy.