Phonon Signatures in Photon Correlations

We show that the second-order, two-time correlation functions for phonons and photons emitted from a vibronic molecule in a thermal bath result in bunching and anti-bunching (a purely quantum effect), respectively. Signatures relating to phonon exchange with the environment are revealed in photon-photon correlations. We demonstrate that cross-correlation functions have a strong dependence on the order of detection giving insight into how phonon dynamics influences the emission of light. This work offers new opportunities to investigate quantum effects in condensed-phase molecular systems.

Correlation measurements of photon emission provide powerful tools for demonstrating quantum effects.Among the most striking examples is the experimental discovery of antibunching in the photon emission of fluorescing atoms [1,2], which provided the first direct demonstration of the quantum properties of light [3].Anti-bunching [4] is the phenomenon whereby the emission of a second photon immediately after a first is suppressed.The joint probability of detecting two photons at time t and t + τ is quantified by the second-order correlation function [5].While in classical emission, this function may have a maximum for τ = 0, second-order correlations falling off as τ → 0 can appear only as a purely quantum phenomenon.The second-order correlation function thus forms a powerful statistical tool that has been used to study fundamental properties of photon emission and photon-mediated interactions, for example bunching and anti-bunching in transmission through waveguides [6] and in emission from plasmonic nanojunctions [7], pattern formation in photoinduced nucleation [8], photon-blockade effects [9] in optical cavities [10][11][12][13][14] (including modified response at strong coupling [15]), atomic arrays [16,17], as well as superatom behavior in ensembles of quantum emitters [17,18].Higher-order correlations can further reveal two-photon blockade [19].
In all of these examples, the correlations are exhibited by emitted light.However, the usefulness of quantum correlation functions extends beyond the study of photons.Similar techniques have also been applied to describe phonon blockade in opto-mechanical [20,21] and spin-mechanical [22,23] systems.Intriguingly, second-order cross-correlation functions between non-identical particles and quasiparticles can reveal, for example, photon-magnon blockade in a ferrimagnetic material coupled to a microcavity [24], and photon-phonon bunching and anti-bunching in a qubit-phonon-plasmon system under strong coupling [25].
In recent years, there has been significant interest in the nature of electronic and vibrational coherence in condensedphase resonance energy transfer in molecular systems [26,27].However, its exact nature, whether classical or quantum, remains controversial [28,29].This is because experiments typically used to investigate this measure optical responses resulting from an induced macroscopic polarization, which is a classical property [30].Quantum-optical techniques, such as correlation measurements [31,32], offer an avenue to investi-gate genuine quantum effects in molecular systems directly.
Here we develop an open quantum system model for a vibronic molecule driven by a continuous monochromatic laser field (Fig. 1).We include both vibrational and electronic degrees of freedom as well as coupling to a thermal environment.We calculate photon-photon, phonon-phonon and photon-phonon correlations during cyclic pumping while the molecule undergoes vibrational relaxation (VR).The key result is that signatures relating to phonon exchange with the environment, which are revealed in phonon-phonon correlations, can be accessible through the measurement of photon emission.By examining photon-phonon cross-correlation functions [24,25,[32][33][34], we explain how phonon dynamics influences the emission of light.Measurements of such features could help elucidate the impact of vibrational excitations on the quantum nature of light-matter interaction processes in systems ranging from subwavelength molecular arrays [35] to large organic molecules [36].
We consider a model system consisting of a simple molecule with two electronic levels, each with a set of N vibrational states, coupled to an infinite ensemble of overdamped quantum harmonic oscillator modes [37][38][39].Similar models have previously been used to investigate the role of a vibrational environment on the open quantum system dynamics of molecules [40,41].
The total Hamiltonian is the formal sum which describes the system, the bath, and the system-bath interaction, as well as the coupling between the system and the photon field.( The nuclear Hamiltonians for the ground and excited electronic states are, respectively, defined by where ω 0 is the system mode frequency, b † and b are system phonon creation and annihilation operators, corresponding to the vibrational states of the molecule, and λ is the system reorganization energy [44].This model is constructed in a diabatic (D) basis that separates the vibrational levels from the electronic states, leading to explicit off-diagonal couplings.As a result, the electronic excited state, Eq. ( 4), appears displaced by ∆ = 2λω −1 0 relative to the electronic ground state [see also Fig. 1(b)].This displacement accounts for the change in the equilibrium geometry of the electronic excited state.Note that an increase in the system reorganization energy, λ, corresponds to an increased displacement, ∆.The excited potential also experiences an energy shift ℏω eg corresponding to the fundamental transition.
Instead of working in the diabatic basis, used above for conceptual clarity in the construction of the Hamiltonian, one can obtain adiabatic (A) eigenstates, |g, 0⟩ , |g, 1⟩ , . . ., |e, 0⟩ , |e, 1⟩ , . . .(see Fig. 1), by diagonalising Ĥtot through a unitary transformation ĤA tot = ( ÛAD ) † ĤD tot ÛAD [41,45].The presence of the system reorganization energy in Eq. ( 4) means that the energy eigenstates (the energy of the adiabatic states) are identical to the energies of the vibrational levels in the diabatic picture, as illustrated in Fig. 1(b).The adiabatic states correspond directly to the laboratory observables.
Vibrational relaxation occurs as a result of escape of system phonons to the environment, modelled as an infinite ensemble of harmonic oscillator modes.The system-bath coupling is then described by where Q = ( b + b † )/ √ 2, and m α , pα and xα are the mass, momentum and the coordinate of the environmental harmonic modes, which correspond to bath phonons.The coupling strength g α of the αth harmonic oscillator is determined by the spectral density, J(ω) = α g 2 α (2m α ) −1 ω α δ(ω − ω α ).This model for the surrounding environment is very general and allows in principle even the modelling of non-Markovian system-bath coupling [40].Here, however, we work in the Markovian limit and simplify the environment to an overdamped Brownian oscillator profile, J(ω) = 2ηωΛ(ω 2 +Λ 2 ) −1 , with bath reorganization energy η and dissipation rate Λ.The overdamped spectral density introduces stochastic, Gaussian fluctuations in the nuclear dynamics, representative of low frequency intermolecular modes from the interaction of the molecule with the solvent.The coupling of the system and bath nuclear coordinates leads to vibrational dephasing [47,48] and dissipation.Several different approaches exist for solving the generally computationally demanding equations of motion resulting from Eq. ( 1) and similar open quantum systems [7,13,31,[49][50][51].Here we employ the hierarchical equations of motion (HEOM) method [52,53] in the overdamped limit from Ref. [40] to evolve the vibronic molecule, equivalent to the Hamiltonian vibration model for a vibronic monomer in Ref. [41].
The hierarchical equations of motion simulations of the quantum dynamics allow us to numerically compute the correlation functions for the emission of photons and phonons from the molecule.In particular, the correlated emission of photons and phonons is quantified by the normalized secondorder correlation function When the operators are chosen such that ĉ1,2 = â = µ eg |g⟩⟨e|, the photon annihilation operator, we obtain the photon-photon correlation function g (2)  aa , which reflects the joint probability of a photon being emitted at time t + τ given that a photon was emitted at time t.By appropriately choosing ĉ1,2 from the photon and phonon operators â and b, respectively, we can correspondingly construct the phonon-phonon correlation function g (2)  bb and, notably, the photon-phonon and phonon-photon cross-correlation functions g (2)  ab and g (2)  ba in a manner similar to Refs.[24,25,[32][33][34].
We now employ the correlation functions ( 6) to probe quantum effects [54] in the molecular system in its thermal environment.In particular, we determine the presence of antibunching [4], defined as g (2)  c 1 c 2 (t, τ = 0) < g (2)  c 1 c 2 (t, τ > 0).This implies that the probability of a second emission event immediately following a first is suppressed.Note that this definition encompasses not only photon-photon or phonon-phonon correlation, but is also generalized [24,25] to include crosscorrelations where the two emission events consist of one photon and one phonon.Correspondingly, bunching is defined to occur when the probability of simultaneous emission is enhanced, g (2)  c 1 c 2 (t, τ = 0) > g (2)  c 1 c 2 (t, τ > 0).In the following we assume that all emitted photons and phonons are detected, regardless of scattering directions, e.g., by imagining the system enclosed by a detector and use the quantum regression theorem [5,42,[55][56][57] to compute second-order correlations as where ρ = ρ(t) is the density matrix at time t and L is the Liouvillian operator for the time evolution of the system.The molecule is initially equilibrated with its thermal environment, in the absence of the driving field, so that the density matrix of the system is correlated with the bath.This ensures that the vibrational states of the molecule are in the correct Boltzmann distribution.The overdamped hierarchy is then used to evolve the dynamics in the presence of the driving field over several optical cycles to find ρ(t).To compute g (2) c 1 c 2 (t, τ), one then takes ĉ2 ρ(t)ĉ † 2 as the initial state for a subsequent evolution in τ taking care to preserve continuity of the driving field.
Figure 2(a)-(c) show the two-time photon correlation function, g (2)  aa (t, τ), as a function of the photon-photon separation time τ for the molecular system defined in Eq. (1).For this and following simulation results, we assume (unless otherwise specified) η = 5 cm −1 , Λ = 200 cm −1 , ω 0 = 500 cm −1 , ∆ = 1.2 such that λ ≈ 260 cm −1 , ω eg = 10 4 cm −1 , E 0 = 10 7 NC −1 , and T = 298 K.These parameters ensure that coupling to the bath is weak, and we are operating in the Markovian limit.We truncate the number of vibrational levels at N = 10, which is sufficiently large for the results to be insensitive to the truncation.These parameters are comparable to real molecules with electronic and vibrational transition frequencies ∼ 10 4 cm −1 and ∼ 10 2 cm −1 , respectively [58,59].The weakly coupled Markovian bath parameters are typical of commonplace nonpolar solvents [40,41].
As population is initialised to a Boltzmann distribution, excitation results in a wavepacket which moves within the harmonic potential [60][61][62].When the system reorganization energy λ = 0, the effect of the monochromatic laser field is to drive population between the ground vibronic states and the equivalent excited states (|g, 0⟩ → |e, 0⟩ , |g, 1⟩ → |e, 1⟩ , ...).The resulting Rabi oscillations are reflected in g (2)  aa and show in photon anti-bunching.However when λ > 0, VR occurs and the excited state wavepacket population becomes different to that of the ground state.This results in the emergence of a minor oscillation in g (2)  aa , at the vibrational mode frequency, which implies that the experimentally measurable second-order photon correlation function contains an observable phonon signature despite the fact that phonons are not directly detected.This signature appears because of the change in Franck-Condon overlap (i.e., the overlap integral of the bound eigenstates of the electronic excited state with the ground state) [63].The Franck-Condon overlap of the fundamental transitions reduces with increasing λ, increasing the Rabi oscillation period and more population enters the vibronic |e, 0⟩ state via VR.
Increased bath reorganization energy η damps the Rabi oscillations as phonons dissipate into the bath, leading to the formation of a steady state [Fig.2(b) and (c)].We then evaluate g (2)  aa (t, τ) as a function of τ at a time t after reaching the steady state, in keeping with common convention.When no steady state forms [η = 0, Fig. 2(a)], but Rabi-like oscillations persist indefinitely [64], the choice of t, and therefore the normalisation of g (2)  aa as a function of τ [cf.Eq. ( 6)] is not obvious [43,45,65].We choose t such that the denominator in Eq. ( 6) corresponds to its value in the steady state when η > 0.
Figure 2(d)-(f) show the corresponding phonon-phonon correlation, g (2)  bb (τ).When λ = η = 0 the population moves resonantly between the ground and excited states and with no VR [66], resulting in a constant phonon-phonon g (2)  bb (τ).Despite lack of phonon dissipation (η = 0), a pronounced oscillation at the mode frequency appears for λ > 0 as the excited state displacement results in a non-stationary population distribution out of thermal equilibrium.For the same reason, g (2)  bb (τ) also tracks the Rabi oscillation when λ > 0 and phonon bunching is apparent.Additionally introducing coupling to the environment bath, η > 0, we find a rapid decay of g (2)   bb with τ due to the strong dissipation.
Both the g (2)  aa and g (2)  bb , by definition, have no dependence on the order of detection events since τ separates the detection of identical particles.In both cases there is a single source of vibrational character.For the g (2)  aa this is indirect, from the strong dependence of photons on the vibrational populations, aa (τ) photon-photon correlation function; (d)-(f) g (2)  bb (τ) phonon-phonon correlation function, scanning over bath (η) and system (λ) reorganization energies.whilst for g (2)  bb this is from direct measurement of the phonon number.In both cases, vibrational character is observed as oscillations at the vibrational mode frequency ω 0 .
We can understand the appearance of phonon signatures in the photon correlations from the cross-correlation functions g (2)  ab (τ) and g (2)  ba (τ) (Fig. 3), where the order of detection does matter.Specifically, the second detection event determines the dominant character of the cross-correlation function as a function of τ.In g (2)  ab , the phonon detection is second, and we observe the primary behaviour of the phonon correlation [cf.Fig. 2 (d)-(f)] with photon correlation-function characteristics superimposed.The first detection event can be thought of as an instantaneous measurement of photon number and contains no vibrational information.The second detection event-the phonon-occurs a time τ later, during which vibrational transitions may occur.However, because the fast phonon signatures are very small with respect to the electronic contributions their impact on the excited-state adiabatic population is minimal, i.e., there is no significant minor oscillation [45].Consequently, neither detection event incurs vibrational character.Figure 3(a)-(c) also show phonon bunching, i.e., a photon detection is likely to be immediately followed by another phonon, reflecting the non-equilibrium population distribution following photon emission.
By contrast, the phonon-photon correlation function g (2)  ba (τ) [Fig.3(d)-(f)] corresponds to the observation of a phonon followed by a photon.Since the photon detection is second, characteristics of g (2)  aa [Fig.2(a)-(c)] dominate, with phonon characteristics superimposed.The first detection event-the phonon-is an instantaneous measurement of phonon number and thus has intrinsic vibrational character at the molecule mode frequency.However, the second detection event-the photon-also introduces additional vibrational character due to vibrational transitions occurring between the detections.Consequently, there are two sources of vibrational character in g (2)  ba : 1. intrinsically from the first detection event, and 2. from phonon effects during the optical cycles leading to the photon emission.
For λ = 0, g (2)  ba (τ) remains at a small, non-zero constant value regardless of the bath coupling.However, increasing the system reorganization energy introduces strong oscillations and the correlation function may drop below the λ = 0 value.This is explained by a large proportion of the wavepacket population undergoing vibrational transitions.This effect persists for small τ even with strong bath dissipation, but it is destroyed for later τ by the influence of the environment.Note that similar to g (2)  aa (τ), g (2)  ba (τ) exhibits anti-bunching-like behavior , i.e., a photon is less likely to be emitted directly after a phonon.This is because photon emission at ω eg from higher vibrational levels is increasingly suppressed for larger λ due to decreasing Franck-Condon overlap.This means that phonon emission tends to inhibit subsequent photon emission when τ ≈ 0.
Dynamical impact of lattice phonons on quantum-dot emitters embedded in a solid-state system was recently theoretically and experimentally shown to result in a characteristic signature in the photon spectrum [67].System phonons are not defined in the quantum-dot model, but are integral to molecules and the proximate source of the signatures predicted here.Our results suggest that measurements of twophoton correlations could explicitly elucidate differences in the impact of system phonons and the phonon environment.
In conclusion, we have demonstrated theoretically photon anti-bunching in the fluorescence of a vibronic molecule under continuous laser drive and a thermal environment and that the photon-photon correlations exhibit signatures of the phonon interaction with the bath, suggesting that these are experimentally directly measurable.These appear as oscillations at the system-mode frequency on top of slower modulations associated with the electronic Rabi-like oscillations.Theoretically also considering phonon detection and photon-phonon cross-correlation functions, we have shown how vibrational contributions are understood as arising either directly, through phonon detection, or indirectly, through photon detection subsequent to phonon emission.As such, the order of particle detection can dramatically impact the behaviour of the correlation functions, which could in principle be exploited to investigate the phonon impact on photon emission.More immediately, these correlation functions present an opportunity to investigate phonon dynamics indirectly using existing quantum- where the weights are given by the Boltzmann distribution with canonical partition function Here ϵ k is the energy of the kth vibrational level in the electronic ground state.The diabatic density matrix is then calculated by rotating equation S16 Calculations are then run in the diabatic basis with a Hilbert space truncation of N = 10.

SIII. SIMULTANEOUS TIME CORRELATION FUNCTIONS AND NON-NORMALISED EQUIVALENTS
In order to calculate correlation functions we have assumed that the incident field is not detected as part of the scattered field.Experimentally can be achieved, for example, by addition of a thin wire as in the dark-ground imaging technique of reference (S1).
The detection probability for photons is given by, where ρ = ρ(t), â † is the associated creation operator, and τ = 0 corresponds to simultaneous time.Exchanging the creation/annihilation operators with phonon equivalents results in a corresponding phonon detection probability.In Fig. S1 (a) g (1)  a tracks the monomer excited state for the four different system reorganization energies.In all cases Rabi oscillations persist through the entire simulation, and the period of the oscillation is dependent on the displacement, ∆, and subsequent reduction of Franck-Condon overlap resulting in a hindered transition into the |e, 0⟩ state.Fig S1 (b) and (c) present the introduction of the bath through the bath reorganization energy, η.The addition of phonon dissipation results in the formation of a steady state in all cases apart from when the displacement of the excited state is zero.When λ = 0 the ground state Boltzmann distribution is able to move perfectly to an equivalent Boltzmann distribution in the excited state resulting in a Rabi oscillation.When λ > 0 and with a bath reorganization energy of η a steady state forms in approximately 5 ps, whereas for 2η, it takes approximately 2.5 ps.
Figure S1 (d), (e) and (f) are the equivalent simultaneous time detection probabilities for a phonon.It is clear that there is a definite minimum value of phonon transfer present within the monomer system which occurs when λ = 0. Any increase in λ (forcing the system away from equilibrium) results in increased VR and a linear increase in the g (1)  b .Excluding the λ = 0 case, all other values of system reorganization energy also exhibit the beating pattern of the electronic Rabi oscillation reflecting the departure from equilibrium.Just as in Fig. S1 (b), and (c), the results for (e) and (f) present a dramatic reduction in the formation time of the steady state when the bath is introduced.The non-normalised, second-order correlation function where ρ = ρ(t), represents the correlation of a photon followed by a phonon.The exchange of these operators results in a phonon-photon correlation G (2)  ba , and versions containing purely â( †) and b( †) operators result in photon-photon and phonon-phonon correlations, G (2)  aa and G (2)  bb , respectively.Figure S2 (a)-(c) demonstrate the underlying physics of the chosen model.Based on the chosen monomer system, at simultaneous time (t, τ = 0) it is impossible to have a double excitation as there is no electronic doubly excited state.As such, we must have a constant zero probability of two photon detection.In contrast, the system has many accessible vibrational states, so simultaneous detection of two phonons is possible [Fig.S2 (d)-(f)].
Similarly, Fig. S3 (a)-(e), are a demonstration of the nature of the second order cross-correlation functions.When at simultaneous time, the order of operations is irrelevant.As such the simultaneous time, τ = 0, cross-correlations are equivalent and a superimposition of the photon detection Rabi oscillations on the behaviour of the phonon detection probabilities.

SIII.A. Normalization
Figures S4-S7, present non-normalised second-order correlation functions at two different values of t.The 0.5 ps increment corresponds to one half of the period for the nondisplaced (λ = 0) Rabi oscillation cycle: moving from 3.5 to 4.0 ps goes from a trough in the photon detection probability to a peak.This is crucial for τ dependent results, as the starting point of the evolution dictates the starting amplitudes for all of the subsequent correlations, which can tend to favour or disfavour certain modes.
Within all of the results, there are two fundamental oscillatory modes: the major electronic Rabi oscillation, and the minor vibrational system mode frequency.The electronic Rabi oscillation is a result of the likelihood of photon emission from either the ground and excited states and is periodic due to the continuous driving by the laser field.The vibrational oscillation is a consequence of changes in the excited state wavepacket population with respect to the ground state Boltzmann distribution, corresponding to phonon transitions.
Based on the excited state population profile within Fig. S1  (a), at a specific time each system is at a different position within its electronic Rabi oscillation and this determines the relative strength of further evolution.In figs S4, and S6, t = 3.5ps and λ 2 , and λ are the closest to their maximum value, whereas, 2λ, and λ = 0 are closest to their minimum.This minimises the impact of the electronic Rabi oscillation within the latter which consequently improves the resolution of the much smaller vibrational contributions.This is especially clear in the appearance of a minor oscillation at the system mode frequency in the G (2)  ab for η = 0 with system reorganization of 2λ in Fig. S6 (a).The opposite is true (ob-scuring vibrational contributions) in figs S5, and S7, where t = 4.0 ps.This motivates our choice to normalise either at the steady state value (for η > 0), or at a time corresponding to an amplitude of one half, equivalent to the steady state average (for η = 0), resulting in a unique time value for each system reorganization energy.

FIG. 1 .
FIG. 1.(a) Schematic of the molecule coupled to bath modes and driven by laser field E I , resulting in the scattered field Êsc .Phonon movement between system and environment indicated by arrows.(b) Diabatic energy levels, with excited state displacement ∆, system reorganization energy λ, fundamental transition frequency ω eg , and system mode frequency ω 0 .Corresponding adiabatic levels on the far right.