Perturbation approach for computing frequency- and time-resolved photon correlation functions

# Perturbation approach for computing frequency- and time-resolved photon correlation functions

## Abstract

We propose an alternative formulation of the sensor method presented in [Phys. Rev. Lett 109, 183601 (2012)] for the calculation of frequency-filtered and time-resolved photon correlations. Our approach is based on an algebraic expansion of the joint steady state of quantum emitter and sensors with respect to the emitter-sensor coupling parameter . This allows us to express photon correlations in terms of the open quantum dynamics of the emitting system only and ensures that computation of correlations are independent on the choice of a small value of . Moreover, using time-dependent perturbation theory, we are able to express the frequency- and time-resolved second-order photon correlation as the addition of three components, each of which gives insight into the physical processes dominating the correlation at different time scales. We consider a bio-inspired vibronic dimer model to illustrate the agreement between the original formulation and our approach.

###### pacs:
42.50 Ar , 03.65 yz, 87.80.Nj, 87.15.hj
\externaldocument

[SI-]SI_new.tex

## I Introduction

Single-photon coincidence measurements have been recognised as a fundamental theoretical and experimental methodology to characterize quantum properties both of light Glauber (1963, 2006); Grangier et al. (1986); Lounis and Moerner (2000) as well as those of the emitting source Olaya-Castro et al. (2001); Michler et al. (2000); Moreau et al. (2001). Particular focus has been placed on investigation of the second-order photon correlation function as the lowest order of correlations capable of probing non-classical phenomena Glauber (2006). Formally, such normally-ordered two-photon correlation function is defined as Vogel and Welsch (2006)

 g(2)(t1,t2)=⟨T−[^A†(t1)^A†(t2)]T+[^A(t2)^A(t1)]⟩⟨^A†(t1)^A(t1)⟩⟨^A†(t2)^A(t2)⟩, (1)

with being the field operator and and the time-ordering and antiordering superoperators necessary for a consistent physical description Vogel and Welsch (2006). Here increases time arguments to the right in products of creation operators, while increases time arguments to the left in products of annihilation operators.

In the context of photon counting experiments it has also become clearer that spectral filtering of optical signals –and its associated trade off between frequency and time resolution– opens up the door for the investigation of variety of phenomena in quantum optics Vogel and Welsch (2006); Sallen et al. (2010); Peiris et al. (2015); Grünwald et al. (2015); Silva et al. (2016). The energy-time Fourier uncertainty relation imposes a constraint on the precision with which arrival time and frequency of a photon can be measured Eberly and Wódkiewicz (1977); Brenner and Wódkiewicz (1982). Rather than being a limitation, this uncertainty has shown to offer a potential for novel investigations of quantum phenomena ranging from the identification and manipulation of new types of photon quantum correlations del Valle et al. (2012, 2016); Gonzalez-Tudela et al. (2013); González-Tudela et al. (2015); Peiris et al. (2015) to the development of new protocols for the preparation and readout of entangled photons Flayac and Savona (2014); Peiris et al. (2017). It has also recently been discussed how frequency- and time-resolved photon correlation measurements can provide insights into the emitter dynamics which are complementary to the information provided by coherent multi-dimensional spectroscopy Brixner et al. (2005) –the later being the ultrafast, non-linear technique capable of probing of quantum coherence dynamics in a variety of biomolecular and chemical systems (for a review see Scholes et al. (2017)).

Filter-dependent correlation functions are defined in terms of filtered emission operators and with the one sided Fourier transform of the frequency filter function. The filtered two-time correlation function can be written as and is defined identically to Eq. (1), but with the substitutions and with the time and space filter functions for each detector Vogel and Welsch (2006); Knoll and Weber (1986); Cresser (1987); Bel and Brown (2009); Kamide et al. (2015); Shatokhin and Kilin (2016). Due to the convoluted definition of , calculating involves computing a four dimensional integral with the time ordering applying within this set of integrals, thereby making such a calculation non-trivial. Higher-order correlations are defined in a similar way, although their theoretical computation becomes more difficult. Thus, a full theoretical understanding of the effects of such filters in the photon statistics has only recently been possible with the development of methods that can overcome the computational complexity del Valle et al. (2012, 2016).

In particular, Ref.del Valle et al. (2012) has put forward an efficient sensor method for calculating these frequency and time-resolved correlation functions, which avoids the need to explicitly compute the multidimensional integral set. The methodology involves accounting, in a quantum mechanical manner, for a weak coupling between the quantum emitter and a set of sensors, each of which is represented as a two-level system. In the limit of vanishing system-sensor coupling, the sensor population correlations are shown to quantify the photon correlations of interest. As originally proposed, this method relies on the explicit use of a numerically small value for the system-sensor coupling. Hence, accurate determination of the photon-correlations function demands testing for convergence. The need for a small parameter may also lead to numerical calculations exhibiting instabilities. Moreover, the method requires solving the quantum dynamics of the joint emitter-sensors state and therefore the dimensionality of the density matrix can become a problem for quantum systems of large dimensions.

In this paper, we report an alternative formulation of the sensor method that addresses the above issues. We propose a formalism that allows us to express photon correlations fully in terms of the quantum dynamics of the emitting system while at the same time eliminating the dependence on a specific value for the small parameter. In our formalism the small parameter vanishes algebraically in all the expressions for multi-photon correlations. Furthermore, using time-dependent perturbation theory, we are able to express the second-order photon correlation as the addition of three components, which provide insight into the physical processes dominating the emission at different time scales.

The paper is organised as follows: Sec. II summarises the original presentation of the sensor method and motivates the development of an alternative formulation. Sec. III presents the approach to derive the steady state system and photon correlations at zero-time delay. Sec. IV explains the derivation of time-dependent correlation functions for finite detection delays, Sec. V illustrates the agreement between our approach and the original sensor method for a bio-inspired vibronic dimer model, and Sec. VI concludes.

## Ii Motivation

As proposed in del Valle et al. (2012), the sensor method for calculating the photon correlation function involves simulating the dynamics of a quantum emitter with Hamiltonian , weakly coupled to sensors represented by two-level systems, labelled with , and ground and excited states and , respectively. Each sensor has an associated Hamiltonian with annihilation operator and transition frequency set to match the emission frequency to be measured. The interaction Hamiltonian between the quantum emitter and the th sensor is given by , with the coupling strength being small enough to neglect back action. For generality, we have considered that the emission operators coupled to each sensor can be different. This is the case whenever local resolution is achievable in a multipartite quantum emitter or when emitting transitions can be distinguished via fluorescence polarization detection as it happens, for instance, in single light-harvesting complexes Tubasum et al. (2011). In such a scenario, the frequency filters illustrated in the envisioned experimental setup (see Fig. 1(a)) will also be polarizing filters.

Considering Markovian relaxation channels for both the emitter and the sensors, the joint emitter-sensors density matrix satisfies the master equation , with the Liouvillian conveniently split as ()

 L(^ρ)=L0(^ρ)+M∑m=1(Lm(^ρ)−i[He,m,^ρ]), (2)

with

 L0(^ρ) Extra open brace or missing close brace (3) Lm(^ρ) =−i[Hm,^ρ]+12ΓmLςm(^ρ). (4)

The superoperators on the right hand side of Eq. (3) have the Lindblad form i.e. for a system jump operator and a relaxation process at rate . Same holds for in Eq. (4) describing the decay of the th sensor with jump operator at a rate . In the limit of satisfying with the smallest transition rate within the emitter dynamics, and sensor populations satisfying , intensity-intensity correlations of the form are directly related to the th order photon correlation functions del Valle et al. (2012, 2016):

 g(M)Γ1...ΓM(ω1,T1;...;ωM,TM)=limϵ1,...,ϵM→0⟨:n1(T1)...nM(TM):⟩⟨n1(T1)⟩...⟨nM(TM)⟩ (5)
 ⟨:n1(T1)...nM(TM):⟩=ϵ21...ϵ2MΓ1...ΓM(2π)MS(M)Γ1...ΓM(ω1,T1;...;ωM,TM) (6)

with Vogel and Welsch (2006)

 S(M)Γ1...ΓM(ω1,⋯,ωM;T1,⋯,TM)==∫∞−∞dt′1∫∞−∞dt′M+1F∗1(T1−t′1)F1(T1−t′M+1)...∫∞−∞dt′M∫∞−∞dt′2MF∗M(TM−t′M)FM(TM−t′2M)⟨T−[a(+)1(t′1)...a(+)M(t′M)]T+[a(−)1(t′M+1)...a(−)M(t′2M)]⟩. (7)

The filter functions correspond to a Cauchy–Lorentz distribution i.e. with the Heaviside function, and can be realised, for instance, via a Fabry-Perot interferometer when the reflection coefficient tends to unity Eberly and Wódkiewicz (1977). The experimental setup we envision is sketched in Fig. 1(a) and the theoretical calculation of frequency-filtered photon correlations through the sensor method is illustrated in Fig. 1(b).

The original presentation of Eq. (5) in del Valle et al. (2012) omitted the normal order. Without the normal order this function yields unphysical results for a finite delay time. In an Erratum del Valle et al. (2016) the authors clarified that normal order is implied through the proof of Eq. (5), though it turns out to be unnecessary for zero time delay. Since Eq. (5) is the departing point of our work, we have carried out a consistency check of its proof as discussed in Appendix VII.

The method as proposed in del Valle et al. (2012) is conceptually clear yet, in practice, its computation involves tackling some numerical challenges. Assuming that all the sensor couplings are identical, , the numerical calculations of photon correlations rely on the choice of a system-sensor coupling that is numerically small, but not so small that adding or subtracting terms of order to or from terms of order causes problems within double precision arithmetic. The procedure then involves checking convergence and stability of the numerical results for different values of . Most importantly, computation of photon correlations at zero-time delay requires to numerically finding the zero eigenvalue of the Liouvillian superoperator associated to the joint emitter system plus sensors. This means that computing when , involves calculating the eigenvector with a zero eigenvalue of a matrix times larger than that of the quantum emitter alone Gonzalez-Tudela et al. (2013). Similarly, for time-resolved correlations, the calculation involves time propagation in the joint state space of the system and sensors. Evidently, as the dimensionality of the system is larger these numerical challenges become more demanding.

We were therefore motivated to find an approach that would allow us to avoid the issues above mentioned. In what follows we show that by expanding algebraically in one can propose an approach that eliminates the explicit numerical dependance on while at the same time reducing the dimensionality of the Hilbert space needed for computation.

## Iii Frequency-filtered spectrum and photon correlations at zero delay time

### iii.1 M=1: power spectrum

We begin by demonstrating the basics of our derivation by considering the emitter system coupled to only one sensor. Let us denote the steady state of the joint emitter-plus-sensor system. From Eq. (6) we can calculate the power spectrum as:

 S(1)Γ1(ω1)=Γ12πϵ2⟨n1⟩=Γ12πϵ2Tr[n1^ρss]. (8)

Considering the identity operator in the sensor Hilbert space i.e. , we can write the full steady state as

 ^ρss=1s1^ρss1s1=∑j1,j′1=0,1^ρj′1j1,⊗|j1⟩⟨j′1|, (9)

where the matrices are therefore only related to the degrees of freedom of the quantum emitter. Hermitian conjugates are obtained by swapping the upper and lower indices. Notice that each matrix is thus of order . With this definition the power spectrum given in Eq. (8) becomes

 S(1)Γ1(ω1)=Γ12πϵ2Tr[^ρ11] (10)

To find the steady state, and in particular the matrix , we solve for the combined quantum emitter plus sensor system . We consider the action of the Liouvillian given in Eq. (2) () on every term in Eq. (9):

 L(^ρ00⊗|0⟩⟨0|)= L0(^ρ00)⊗|0⟩⟨0| (11) −iϵ(a1^ρ00⊗|1⟩⟨0|−^ρ00a†1⊗|0⟩⟨1|), L(^ρ01⊗|1⟩⟨0|)= (L0−Γ1/2−iω1)(^ρ01⊗|1⟩⟨0| (12) −iϵ(a†1^ρ01⊗|0⟩⟨0|−^ρ01a†1⊗|1⟩⟨1|), L(^ρ11⊗|1⟩⟨1|)= (L0−Γ1)(^ρ11⊗|1⟩⟨1|) (13) +Γ1^ρ11⊗|0⟩⟨0| −iϵ(a†1^ρ11⊗|0⟩⟨1|−^ρ11a1⊗|1⟩⟨0|),

and the expression for is the complex conjugate of Eq. (III.1). We can rewrite the sum of these expressions, in a similar way to Eq. (9), by grouping together terms related to populations or coherences of the sensor:

 L(^ρss) =∑j1,j′1=0,1^Bj′1j1,⊗|j1⟩⟨j′1|=0, (14a) ^Bj′1j1, =0,For all j1,j′1, (14b)

In this way we can see our problem reduces to solving the set of coupled equations for such that the operators of the system are null matrices (zero at every element). Notice that Eq. (11) has one term contributing to , which is zeroth order in , and terms contributing to and , which are of linear order in . Similarly, Eq. (III.1) contributes terms to , as well as to and via a term proportional to . Hence, the equation to be solve for , for instance, becomes

 ^B00=L0(^ρ00)−iϵ(a†1^ρ01−^ρ10a1)+Γ1ρ11=0. (15)

For an arbitrary value of , the set of coupled equations given by Eq. (14b) does not have a simple solution. However, in the limit of weak coupling where and , we can neglect terms of the order of . For instance, for (Eq. (15) ) the terms and are of the order of , so will be neglected. Likewise, in this weak coupling limit we have in the equation for . These approximations can be generalised to a concept of ignoring down coupling, that is, the prefactor matrix for has only contributions from terms with . This is equivalent to a formal expansion in as all the matrices are of order . Using these approximations, we can write out the equations governing the steady state as

 L0(^ρ00)∼0 (16a) L0(^ρ01)−(Γ1/2+iω1)^ρ01−iϵa1^ρ00∼0 (16b) L0(^ρ10)−(Γ1/2−iω1)^ρ10+iϵ^ρ00a†1∼0 (16c) L0(^ρ11)−Γ1^ρ11−iϵ(a1^ρ10−^ρ01a†1)=0. (16d)

We can solve these equations in a chain from top to bottom, starting with . In practice, we need not solve for as it is equal to . Numerically we formulate the problem in Liouville space, such that is the zero eigenvector of the (square) matrix given in Eq. (3). The remaining equations can be solved as

 |^ρ01⟩⟩ ∼iϵa1|^ρ00⟩⟩L0−(Γ1/2+iω1)1 (17a) |^ρ11⟩⟩ Missing or unrecognized delimiter for \right (17b)

and are written in the Liouville space form, and is the identity operator in the emitter Hilbert space. Notice Eq. (17)(b) has an equality as for that case no term is discarded. In the above equations has a prefactor of and has a prefactor of . Therefore the dependence of the power spectrum (Eq. (10)) on vanishes algebraically. The numerical calculation of the matrices given by Eqs. (17) can, in principle, be done by carrying over a small value for . This procedure, however, could lead to numerical instabilities due to the smallness of . With our method, such instabilities are prevented by computing the re-scaled matrices (such that and ), which are now -independent system operators. From the trace of the -independent matrix we can calculate the sensor count rate as:

 SΓ1(ω1)=Γ12πTr[~^ρ11]. (18)

### iii.2 M≥2 zero-delay correlations

The normalised second-order () photon correlation at zero delay time can be written as

 g(2)Γ1Γ2(ω1,ω2)=S(2)Γ1,Γ2(ω1,ω2)S(1)Γ1(ω1)S(1)Γ2(ω2) (19)

where and are the mean count rates for the two sensors, as given in Eq. (8), and:

 S(2)Γ1,Γ2(ω1,ω2)=Γ1Γ2(2π)2ϵ4⟨:n1n2:⟩ (20)

Since time-independent sensor number operators commute, normal order in Eq. (20) is unnecessary. Following the same procedure as before, we write our steady state density matrix, with two sensors included, as:

 ^ρss=∑j1,j2,j′1,j′2=0,1^ρj′1,j′2j1,j2⊗|j1⟩⟨j′1|⊗|j2⟩⟨j′2| (21)

where and are counters over the states of sensor 1 and sensor 2, respectively. As before, the matrices are defined in the Hilbert space of the quantum emitter alone and can be re-scaled as . With this definition the second-order photon coincidence becomes

 S(2)Γ1,Γ2(ω1,ω2)=Γ1Γ2(2π)2Tr[~^ρ1,11,1], (22)

with the power spectrum re-defined as

 S(1)Γ1(ω1)=Γ12πTr[~^ρ1,01,0],S(1)Γ2(ω2)=Γ22πTr[~^ρ0,10,1]. (23)

To compute the matrices we solve for the steady state with two sensors and by ignoring down coupling terms, that is, the matrix prefactor for has only contributions from terms satisfying the condition . The resultant full set of linearly independent equations (besides those which are Hermitian conjugates of others) are then given by

 L0(~^ρ0,00,0)∼0 (24a) [L0−Γ1/2−iω1](~^ρ0,01,0)∼ia1~^ρ0,00,0 (24b) [L0−Γ2/2−iω2](~^ρ0,00,1)∼ia2~^ρ0,00,0 (24c) [L0−Γ1](~^ρ1,01,0)∼i(a1~^ρ1,00,0−~^ρ0,01,0a†1) (24d) [L0−Γ2](~^ρ0,10,1)∼i(a2~^ρ0,10,0−~^ρ0,00,1a†2) (24e) [L0−(Γ2+Γ1)/2−i(ω1+ω2)](~^ρ0,01,1)∼ i(a1~^ρ0,00,1+a2~^ρ0,01,0) (24f) [L0−(Γ2+Γ1)/2−i(ω1−ω2)](~^ρ0,11,0)∼ i(a1~^ρ0,10,0−~^ρ0,01,0a†2) (24g) [L0−(Γ1/2+Γ2)−iω1](~^ρ0,11,1)∼ i(a1~^ρ0,10,1+a2~^ρ0,11,0−~^ρ0,01,1a†2) (24h) [L0−(Γ2/2+Γ1)−iω2](~^ρ1,01,1)∼ i(a1~^ρ1,00,1−~^ρ0,01,1a†1+a2~^ρ1,01,0) (24i) [L0−(Γ1+Γ2)](~ρ1,11,1)= i(a1~^ρ1,10,1−~^ρ0,11,1a†1+a2~^ρ1,11,0−~^ρ1,01,1a†2). (24j)

The solutions, in analogy to those in Eq. (17), are the following:

 |~^ρ0,01,0⟩⟩∼ia1|~^ρ0,00,0⟩⟩L0−(iω1+Γ1/2)1, (25a) |~^ρ0,00,1⟩⟩∼ia2|~^ρ0,00,0⟩⟩L0−(iω2+Γ2/2)1, (25b) |~^ρ1,01,0⟩⟩∼i(a1|~^ρ1,00,0⟩⟩−|~^ρ0,01,0⟩⟩a†1)L0−Γ11, (25c) |~^ρ0,10,1⟩⟩∼i(a2|~^ρ0,10,0⟩⟩−|~^ρ0,00,1⟩⟩a†2)L0−Γ21, (25d) |~^ρ0,01,1⟩⟩∼i(a1|~^ρ0,00,1⟩⟩−a2|~^ρ0,01,0⟩⟩)L0−(iω1+iω2+Γ1/2+Γ2/2)1, (25e) |~^ρ0,11,0⟩⟩∼i(a1|~^ρ0,10,0⟩⟩−|~^ρ0,01,0⟩⟩a†2)L0−(iω1−iω2+Γ1/2+Γ2/2)1, (25f) |~^ρ0,11,1⟩⟩∼i(a1|~^ρ0,10,1⟩⟩+a2|~^ρ0,11,0⟩⟩−|~^ρ0,01,1⟩⟩a†2)L0−(iω1+Γ1/2+Γ2)1, (25g) |~^ρ1,01,1⟩⟩∼i(a1|~^ρ1,00,1⟩⟩−|~^ρ0,01,1⟩⟩a†1+a2|~^ρ1,01,0⟩⟩)L0−(iω2+Γ2/2+Γ1)1, (25h) |~^ρ1,11,1⟩⟩=i(a1|~^ρ1,10,1⟩⟩−|~^ρ0,11,1⟩⟩a†1+a2|~^ρ1,11,0⟩⟩−|~^ρ1,01,1⟩⟩a†2)L0−(Γ1+Γ2)1. (25i)

Generalisation of this formalism for the -order frequency-resolved correlation at is straightforward and requires writing out the general steady state for the emitter and sensors in the form analogous to Eqs. (9) and (21):

 ^ρss=∑j1,j′1,…,jM,j′M=0,1^ρj′1…j′m…j′Mj1…jm…jM⊗|j1⟩⟨j′1|⊗…⊗|jM⟩⟨j′M|. (26)

where are counters over the state of sensor and . We define the re-scaled matrices . The th order photon-coincidence at zero-delay time is given in terms of the trace of matrix with for all

 S(M)Γ1...ΓM(ω1,...ωm,...ωM)=Γ1...Γm...Γ