###### Abstract

We theoretically investigate polarization-filtered two-photon correlations for
the light emitted by a multichromophoric system undergoing excitation
transport under realistic exciton-phonon interactions, and subject to
continuous incoherent illumination. We show that for a biomolecular aggregate,
such as the Fenna-Matthews Olson (FMO) photosynthetic complex, time-
asymmetries in the cross-correlations of photons corresponding to different
polarizations can be exploited to probe both quantum coherent transport
mechanisms and steady-state coherence properties, which are not witnessed by
zero-delay correlations. A classical bound on correlation asymmetry is
obtained, which FMO is shown to violate using exact numerical calculations.
Our analysis indicates that the dominant contributions to time-asymmetry in
such photon cross-correlations are population to coherence transfer for
Frenkel-Exciton models. Our results therefore put forward photon correlation
asymmetry as a promising approach to investigate coherent contributions to
excited-stated dynamics in molecular aggregates and other many-site quantum
emitters.

Introduction:– The direct observation of quantum coherent phenomena in complex
systems is a fundamental interdisciplinary scientific challenge that aims to
advance our understanding of quantum processes in solid-state Stotz  _et al._
(2005); Tran  _et al._ (2016); Reigue  _et al._ (2017), biomolecular, and
chemical systems Scholes  _et al._ (2011); Huelga and Plenio (2013); Scholes
_et al._ (2017); Polisseni  _et al._ (2016); Grandi  _et al._ (2019); Mattioni
_et al._ (2024). In these scenarios, the rich interaction of the system of
interest with phonon environments makes the task of characterising quantum
behaviour challenging. New approaches to probing the influence of coherent
mechanisms in the dynamics of such systems, under relevant phonon
environments, then promise to open new pathways for the development of robust
quantum technologies Toninelli  _et al._ (2021).

Crucial in each of the aforementioned fields is the development of
experimental approaches which may probe effects of quantum coherence in energy
transport Hanschke  _et al._ (2018); Bluvstein  _et al._ (2022); Fassioli and
Olaya-Castro (2010); Scholes  _et al._ (2011); Fassioli  _et al._ (2012); Chin
_et al._ (2013); O’Reilly and Olaya-Castro (2014); Romero  _et al._ (2014,
2017); Scholes  _et al._ (2017); Kim  _et al._ (2021); Nation  _et al._
(2023). In particular, verifying quantum coherence in photosynthetic
excitation transport has been a topic of considerable debate Brixner  _et al._
(2005); Olaya-Castro  _et al._ (2008); Collini  _et al._ (2010); Wilkins and
Dattani (2015); Holdaway  _et al._ (2016, 2017); Duan  _et al._ (2017);
Scholes  _et al._ (2017); Knee  _et al._ (2018); Kim  _et al._ (2021); Werren
_et al._ (2023), which was first given experimental credence by optical
multidimensional spectroscopic observations of long-lived ‘coherence beats’ in
the two-dimensional electronic spectra of the Fenna-Matthews-Olson (FMO)
complex, first at cryogenic temperatures Engel  _et al._ (2007), and later at
room temperature Panitchayangkoon  _et al._ (2010). Such nonlinear
spectroscopic methods have indeed become the cornerstone of experimental
procedures for the investigation of ultrafast excitation dynamics in organic
molecules and materials 2Dp (2023). They however rely on coherent pulsed laser
light that, while giving access to a wealth of information, probe molecular
systems in regimes outside of their natural conditions, such as weak
incoherent irradiation Fassioli  _et al._ (2012). Furthermore, such methods
rely on the optical response of an ensemble of molecules whose properties vary
from one to another, making it difficult to discern the quantum processes that
are relevant at the single-molecule level.

Since the seminal work of Glauber Glauber (1963), photon correlations have
become the defining measurements to characterise the quantum nature of light,
and are equally pivotal to assess the quantum nature of emitters Carmichael ;
Aspect  _et al._ (1980); Schrama  _et al._ (1992); Campisi  _et al._ (2011);
Holdaway  _et al._ (2018); del Valle  _et al._ (2012); Anderson  _et al._
(2018); Velez  _et al._ (2019); Sánchez Munõz and Schlawin (2020); Miller  _et
al._ (2021); Landi  _et al._ (2023); Cygorek  _et al._ (2023); Downing  _et
al._ (2023). More recently photon-counting experiments have been performed on
larger and more complex systems Hollars  _et al._ (2003); Hedley  _et al._
(2021), probing even single large supramolecular complexes Wientjes  _et al._
(2014); Kim  _et al._ (2019); Li  _et al._ (2023). The study of photon
correlations thus promises new insights into quantum features of single
multichromophoric systems, as experimental progress has made near-term study
of coherence in energy transport through such methods viable. What is still
uncertain are the fundamental signatures of quantum behaviour in multi-
chromophoric systems which may be extracted from such experiments.

Recent studies Iles-Smith  _et al._ (2017); Holdaway  _et al._ (2018); Cygorek
_et al._ (2023); Li  _et al._ (2023); Humphries  _et al._ (2023); Nation  _et
al._ (2024) have indeed illustrated the use of second-order photon
correlations for the investigation of quantum coherent phenomena. Photon
correlations can be bounded under classical assumptions Loudon (1980);
Auffèves  _et al._ (2011); Sánchez Munõz and Schlawin (2020), and thus a
violation of such bounds demonstrates clearly the presence of coherence.
However, theoretical investigations of photon-correlations that account for
the complexity of the electron-phonon interactions characteristic of the
systems of interest have just started to emerge Humphries  _et al._ (2023);
Nation  _et al._ (2024). We have previously shown Nation  _et al._ (2024) that
frequency-filtered and time-resolved photon cross-correlations, and
particularly their asymmetry on interchange of measurement order, can be
exploited to probe vibronic couplings in a composite emitter system, and the
associated coherent dynamical process.

In this work we investigate polarization-filtered photon correlations of the
light emitted by a multichromophoric system of several interacting chromophore
emitters, each locally coupled to realistic phonon environments. To do this,
we deploy the hierarchical equations of motion (HEOM) formalism to open
quantum system dynamics and include non-additive incoherent pump and decay
channels. We systematically investigate the cases of two- and three-site
emitters, and the FMO complex, and show that for small size systems, the
violation of classical bounds for the zero-delay photon correlations is
remarkably robust to environmental influence. However, since the classical
bound grows with the number of coupled emitters, zero-delay photon
correlations are not sensitive enough to witness coherence in large systems
such as FMO. We then present our main result: time-asymmetry of photon cross-
correlations for a multichromophoric system interacting with phonon
environments is a sensitive measure to witnesses both coherence properties of
the non-steady state as well as the underlying dynamical coherent processes,
for which we provide analytical arguments in the supplemental material (SM) sm
. Notably, there is no single known quantum measure able to capture both
steady-state and dynamically induced coherences in quantum processes Streltsov
_et al._ (2017). The polarization filtered correlations here investigated
capture several of the relevant properties of such a measure when applied to
quantum transport in multichromophoric systems.

Set-up:–

Figure 1: Zero delay photon correlations for different relative filter angles
$\phi_{1}$ and $\phi_{2}$. a)-b) Show the dimer model with intersite coupling
of $J=-46.8$ cm-1 (see SM sm  for additional details) with a) intermediate
environmental coupling $\lambda=35$ cm-1 and b) strong environmental coupling
$\lambda=100$ cm-1 (PC645 value). c), d) Show LH2 subunit with c)
$\max(J)=140$ cm-1 = $\lambda_{1}=\lambda_{2}$ and d) $\max(J)=339$ cm-1 (LH2
subunit value). f) Zero delay correlations for the FMO complex. e) Shows the
autocorrelations for equal polarization angles of FMO (diagonal elements of
f)). Additional model details are given in SM sm .

Collective electronic excitations of molecular aggregates are well described
by a Frenkel exciton Hamiltonian Mukamel (1995); Cho , with $M$ chromophores
(sites) and a maximal number of $N_{ex}=2$ excitations,

$\displaystyle H_{S}=\sum_{m}^{M}\epsilon_{m}|m\rangle\langle
m|+\sum_{\begin{subarray}{c}m,n\\\
m>n\end{subarray}}^{M}J_{mn}(|m\rangle\langle n|+|n\rangle\langle m|)$ (1)
$\displaystyle+\sum_{\begin{subarray}{c}m,n\\\
m>n\end{subarray}}^{M}(\epsilon_{m}+\epsilon_{n})|m,n\rangle\langle
m,n|+\sum_{\begin{subarray}{c}m,n,l\\\ (m>n)\neq
l\end{subarray}}^{M}J_{mn}(|n,l\rangle\langle m,l|+|m,l\rangle\langle n,l|),$

with the state $|m\rangle$ denoting a single excitation on site $m$, and
$|m,n\rangle$ denoting excitations on sites $m$ and $n$. The environment is
modelled by a set of local harmonic oscillators with creation (annihilation)
operator $b^{\dagger}_{\xi,m}\,(b_{\xi,m})$ for mode $\xi$ on site $m$,
$H_{B}=\sum_{m,\xi}\omega_{\xi,m}b^{\dagger}_{\xi,m}b_{\xi,m},$ each linearly
coupled to respective sites via $H_{SB}=\sum_{m}|m\rangle\langle
m|\sum_{\xi}g_{\xi,m}(b_{\xi,m}+b^{\dagger}_{\xi,m})+\sum_{n}\sum_{m\neq
n}|m,n\rangle\langle m,n|\sum_{\xi}g_{\xi,m}(b_{\xi,m}+b^{\dagger}_{\xi,m}).$
We have included the contribution of doubly excited electronic states, without
which polarization filtered two-photon correlations trivially vanish at zero-
delay. Details of the definition of the second excited manifold for arbitrary
operators are given in the SM sm .

We capture the influence of the environment via a Drude-Lorentz spectral
density Mukamel (1995) for each site,
$J_{m}(\omega)=\Theta(\omega)2\lambda_{m}\Omega_{m}\frac{\omega}{\omega^{2}+\Omega_{m}^{2}},$
where $\Theta(\omega)$ is the Heaviside step function, $\Omega_{m}$ is the
cutoff frequency, and the reorganisation energy is
$\lambda_{m}:=\frac{1}{\pi}\int_{-\infty}^{\infty}d\omega
J_{m}(\omega)/\omega$. On-site energies are given by
$\epsilon_{m}=\epsilon_{m}^{0}+\lambda_{m}$ with $\epsilon_{m}^{0}$ the bare
electronic excitation energy for site $m$. Exact open system dynamics with the
Drude-Lorentz environment is obtained via the HEOM method Ishizaki and
Tanimura (2006); Jin  _et al._ (2008); Ishizaki and Fleming (2009); Tanimura
(2014); Ke and Zhao (2017).

Interaction with the electromagnetic environment
$H_{EM}=\sum_{m}\sum_{p}\omega_{m,p}a^{\dagger}_{m,p}a_{m,p}$ can support
photonic excitations of three polarizations $p\in\\{x,y,z\\}$, and each couple
in the dipole approximation via
$H_{I}=\sum_{i}\sum_{p}\sum_{m}^{N}g_{i,m}(\hat{\mu}_{m,p}a^{\dagger}_{i,p}+\hat{\mu}_{m,p}^{\dagger}a_{i,\lambda})$.
The the transition-dipole annihilation operators for polarization $p$ are
given by $\hat{\mu}_{p}=\sum_{m}|0\rangle\langle
m|\sigma_{m}+\sum_{m>n}(\mu_{a,n}|m\rangle\langle
mn|+\mu_{a,m}|n\rangle\langle mn|))$ , and may be calculated from the
transition-dipole moments $\mu_{m,x},\,\mu_{m,y},\,\mu_{m,z}$ on each site. We
assume incoherent pumping of the highest energy excitonic state (eigenstate of
$H_{S}$ on the singly excited subspace) at rate $P_{X_{1}}$, and an incoherent
decay of each exciton at rate $\gamma$, unless otherwise stated these are
taken to be $\gamma=1$ ns-1 and $P_{X_{1}}=0.1$ ns-1. We show in the SM sm
(Eq. (77)) a consistent approach to combining in a non-additive manner these
incoherent environmental contributions with the HEOM approach used for the
local vibrational environments, in a similar framework to Ref. Fay and Limmer
(2022).

We will analyse three models of the form given above. These are a heterodimer
($M=2$) model, for which we choose parameters to resemble the central dimer in
the photosynthetic antenna phycocyanin 645 Mirkovic  _et al._ (2007); Blau
_et al._ (2018), with $\lambda=\Omega=100$ cm-1 unless otherwise specified, as
well as a subunit of ($M=3$) of the Light-harvesting 2 (LH2) complex Cogdell
_et al._ (2006); Cupellini  _et al._ (2020) found in purple bacteria. The
latter has site-dependent reorganisation energies with the higher energy two
sites having $\lambda_{1}=\lambda_{2}=140$ cm-1, and the lower energy site
having $\lambda_{3}=40$ cm -1. These simple models will be useful to
characterise bounds on zero-delay correlations, where we will see that system
size $M$ is of key importance. Our main results are then shown for the full
Fenna-Matthews-Olson (FMO) complex Kell  _et al._ (2016); Hein  _et al._
(2012a); Kramer  _et al._ (2018); Sarovar  _et al._ (2010), with $M=8$,
reorganization energy $\lambda=35$ cm-1 and cutoff $\Omega=106.1$ cm-1.
Further model details are shown in the SM. All calculations are performed at
300 K.

Photon correlation functions give information on the conditional probability
of time-delayed photon coincidence events Scully and Zubairy (1997). That is,
given the detection of a photon with polarization $a$ at $t=0$, the studied
correlations describe the conditional probability of detection of a photon of
polarization $b$ at time $t=\tau$. The second-order photon correlation
function reads

$\displaystyle
g^{(2)}_{ab}(\tau)=\frac{\langle\hat{\mu}_{a}^{\dagger}\hat{\mu}_{b}^{\dagger}(\tau)\hat{\mu}_{b}(\tau)\hat{\mu}_{a}\rangle}{\langle\hat{\mu}_{a}^{\dagger}\hat{\mu}_{a}\rangle\langle\hat{\mu}_{b}^{\dagger}\hat{\mu}_{b}\rangle},$
(2)

where $\langle\cdots\rangle:=\operatorname{Tr}[\cdots\rho]$, with $\rho$ the
steady-state of the system. Assuming the quantum regression theorem, we
replace
$\hat{\mu}_{b}^{\dagger}(\tau)\hat{\mu}_{b}(\tau)\to(\hat{\mu}_{b}^{\dagger}\hat{\mu}_{b})(\tau)$,
where $O(\tau)=e^{\mathcal{L}^{\dagger}\tau}O$ is the Heisenberg time
evolution of the operator $O$ under the dual map $\mathcal{L}^{\dagger}$,
which is calculated via HEOM as described in the SM. Negative time photon
correlations are defined as the reversal of the order of measurements on the
system, that is, swapping of the indices $a\leftrightarrow b$. The dynamical
map $\mathcal{L}$ is the generator of evolution of the system state, with
steady-state $\mathcal{L}[\rho]=0$. We note the inclusion of the double
excited subspace greatly increases the complexity of the HEOM calculations, as
the effective state dimension including the ground state is
$d=\dim(H)=M(M-1)/2+M+1$. We thus exploit multiple speedups, such as
auxiliary-operator filtering Shi  _et al._ (2009), and an efficient steady-
state solver Ye  _et al._ (2016) to overcome prohibitive scaling of numerical
complexity.

In simulating the emission of molecular systems Mukamel (1995); Adolphs and
Renger (2006); Ishizaki and Tanimura (2006); Hein  _et al._ (2012b); Gelzinis
_et al._ (2015), in Eq. (2) the operators $\mu_{a}$ are the transition-dipole
operators. In the SM we show how these may be obtained for an arbitrary
polarization filter and transition-dipole moments. We show in the SM that over
the nanosecond timescale of decay to the electronic ground state, photon
correlations indeed tend to their uncorrelated value
$g^{(2)}_{ab}(\tau\to\infty)=1$. Here we focus on picosecond timescales of
interaction with the vibrational environment.

We note that current experiments are able to obtain delay time-resolution of a
few picoseconds Korzh  _et al._ (2020), whilst this may limit a complete
characterization of time-resolved correlations, the total time-asymmetry is
more readily accessible. This is defined as
$\mathcal{A}[f(t)]=||f(t)-f(-t)||,$ where $||f(t)||=\sqrt{h\sum_{k}f(t_{k})}$,
$h=\Delta t/N_{points}$, with $\Delta t$ and $N_{points}$ the range in and
number of time values, respectively. We observe that this total asymmetry is
related to both steady-state coherences and dynamical evolution of coherences
in the system, and investigate this looking at the asymmetry over differing
timescales. In order to capture the behaviour of steady-state coherences in
the system, in the following we compare the asymmetry to the basis independent
coherence Ma  _et al._ (2019); Le and Olaya-Castro (2020), defined as,

$\displaystyle C_{1}(\rho)=S(\rho||\mathbb{1}/d)=\log_{2}(d)-S(\rho),$ (3)

where $S(\rho)=-\operatorname{Tr}[\rho\log_{2}(\rho)]$ is the von Neumann
entropy, and $S(\rho||\sigma)=-\operatorname{Tr}[\rho\log(\sigma)]-S(\rho)$ is
the relative entropy; $C_{1}(\rho)$ hence measures a distance from the
infinite temperature state, and behaves similarly to the purity of the state.

Zero-delay coherence witness:– Zero delay photon-correlations may be exploited
to witness the existence of steady-state quantum coherence via the violation
of a bound Sánchez Munõz and Schlawin (2020): $g^{(2)}_{ab}(0)\leq
2-\frac{2}{M}$, which we label the zero-delay bound. Here, in Figure 1 we show
zero-delay polarization filtered correlations for each model. Notably, the
structure of bunched/anti-bunched regions is dominated by the molecule-
detector orientation. The regions where the zero-delay bound is violated are
designated by the inside of the black solid lines, indicating the presence of
steady-state excitonic coherence.

Figure 1a), b) shows the results for the dimer model for two reorganisation
energies, with in Figure 1b) $\lambda=100$ cm-1, much larger than the
electronic coupling of $-46.8$ cm-1. We note an increase of the region of
violation of this bound when strong environmental coupling is present, which
is consistent with the increased excitonic coherences observed as $\lambda$ is
increased, shown in SM. In the SM we show that strong Markovian pure dephasing
can remove violation of the bound under different pumping conditions as in
Ref. Sánchez Munõz and Schlawin (2020). For $M=2$, violation of the zero-delay
bound is robust to the presence of strong environmental coupling, leading to
the potential to witness coherence in, for example, fluorescent proteins Curti
_et al._ (2023); Sánchez-Pedreño Jiménez  _et al._ (2023) with zero-delay
photon cross-correlations.

We note that measurement of the zero-delay photon correlations is
experimentally challenging in the presence of fast environmental dephasing
Auffèves  _et al._ (2011), as the required time resolution is lower than the
dephasing time, typically from hundreds of femtoseconds to a few picoseconds.
More easily observed in experiment is then a dephased quasi-equilibrium value,
which is typically lower than $g^{(2)}(0)$ Nation  _et al._ (2024), as
observed below.

Figure 1c), d) show the results for the LH2 subunit Cupellini  _et al._
(2020), looking at two cases with altered coupling between the most strongly
coupled sites: Figure 1c) $\max(J)\approx\lambda$ and Figure 1d)
$\max(J)\gg\lambda$ (note that the latter is the biologically relevant value).
We see a region of zero-delay bound violation when $\max(J_{mn})>\lambda$,
however for reduced couplings the bound is no longer violated, despite the
presence of steady-state coherences (see Figure 2).

In Figure 1e), f) we show zero delay photon cross-correlations for the FMO
complex, and observe that the bound is not violated for any combination of
polarization filters, even in the absence of the vibrational environment (see
SM). This is due to the fact that the zero-delay bound increases towards the
value 2 for large $M$, and thus for multichromophoric systems such as FMO,
additional bunching is required for the bound to be violated. We also note
that with knowledge of the dipole moments one can strengthen the bound, and
that in this case the bound is similarly not violated (see SM). Thus, as $M$
is increased the zero-delay bound is less sensitive to the presence of steady-
state coherence (see Figure 3h)).

Figure 2: Time-dependence a), b), and total asymmetry c), e) of photon
correlations and their comparison to basis independent coherence d), f). e)
Shows integration times close to (orange squares, right axis) and larger than
(blue circles) the system-environment coherence time
$\tau_{S-E}\sim\frac{1}{\sqrt{\Omega\lambda}}$. Dependence on reorganization
energy $\lambda$ and coupling strength scaling parameter $\alpha$ shown. HEOM
truncation = 20. Figure 3: Time dependence of polarization filtered photon
correlations for different polarizations a)-c) and their asymmetry for d) for
the FMO complex. Dashed yellow line shows the maximal value of the classical
bound on correlation asymmetry, which is that of the $x,z$ cross correlation.
e), f) Show pump power dependence of the total asymmetry of photon
correlations out to 50 ps, and $C_{1}(\rho)$ (Eq. (3)), respectively. g), h)
Show the same for dependence on reorganization energy $\lambda$. HEOM
truncation is 4, temperature = 300K. We note over nanosecond timescales not
shown here each $g^{(2)}_{ab}(\tau)\to 1$.

Asymmetric photon correlations and quantum coherent transport:– Our central
study is the relation of the behaviour of photon cross-correlations to quantum
coherent exciton transfer in photosynthetic biomolecular systems. In the SM
(sections II-IV) we provide analytical arguments under the assumption that
photon emission processes occur much slower than exciton transfer via an
approximation assuming dipole transitions occur predominantly between
excitonic states. We show that the dominant sources leading to time-
asymmetrical photon cross-correlations are due to excitonic coherences in the
steady-state or those generated by early-time Liouvillian dynamics, which we
refer to as dynamical coherences. Dynamical coherences are additionally a
relevant quantum feature of excitation transport that we wish to probe the
influence of in the system. In the following we will not attempt do
distinguish contributions to correlation asymmetry from dynamical or steady-
state coherence, rather, to probe the presence of either form of coherence.

We recover symmetric correlations for classical multi-level emitters obeying
detailed balance, however, we observe asymmetries may occur in a semi-
classical description of coupled two-level systems. These asymmetries can be
bounded from the assumptions of an incoherent steady-state, combined with the
assumption that excitonic coherences do not contribute to dynamics (see SM,
section II), from which we obtain

$\displaystyle\mathcal{A}[g^{(2)}_{ab}(\tau)]\leq\frac{\sum_{\alpha,\alpha^{\prime}}|(|\hat{\mu}_{b\alpha^{\prime}}|^{2}|\hat{\mu}_{a\alpha}|^{2}-|\hat{\mu}_{a\alpha^{\prime}}|^{2}|\hat{\mu}_{b\alpha}|^{2})|}{S^{(1)}_{a}S^{(1)}_{b}}p^{(2)},$
(4)

where $S^{(1)}_{a}=\operatorname{Tr}[\mu^{\dagger}_{a}\mu_{a}\rho]$. We note
this bound is proportional to the population of the doubly excited subspace
$p^{(2)}$, and thus expected to be very small. We then analyse the
contributions to the Liouvillian of a Frenkel-Exciton model with local pure-
dephasing, which may induce asymmetric cross-correlations. We show that in a
non-Hermitian Hamiltonian approach the asymmetry inducing processes are those
coupling excitonic coherences and populations. Including the influence of
quantum jumps, we show that only a single term describing exciton population-
population transfer can contribute to photon cross-correlation asymmetry,
which appears as a contribution due to excitonic interference enabled by
delocalisation. Additional terms which lead to an asymmetry of photon
correlations either arise from coherence-coherence or population-coherence
coupling terms in the master equation, or coherences in the steady-state,
suggesting steady-state and/or dynamical coherence is the dominant source of
correlation asymmetry.

Photon correlation asymmetry can be seen to behave as a witness of coherent
transport in Figure (2), where we study the dependence of the total asymmetry
on both electronic couplings and reorganization energy $\lambda$ in the LH2
subunit. We study the dependence on electronic couplings by scaling each
coupling strength by a unitless parameter $\alpha$, sending $J_{mn}\to\alpha
J_{mn}$ in Eq. (S3). We see that when altering $\alpha$, both $C_{1}(\rho)$
and total asymmetry each have very similar dependence. When altering the
reorganization energy, however, there more subtle behaviour owing to
environmentally induced dynamical coherence at short times. We first note that
the excitonic coherences, which increase with $\lambda$, and $C_{1}(\rho)$,
have an opposite dependence on $\lambda$. This is due to the fact that
increasing $\lambda$ leads to localisation of the steady-state in the site
basis, increasing coherences in the excitonic basis, whereas the basis
independent measure captures the decrease in entropy due to this localisation,
leading to decreased $C_{1}(\rho)$. The total asymmetry, on the other hand,
has a dependence on $\lambda$ that depends on the timescale of integration. We
define a system-environment coherence time as
$\tau_{S-E}\sim\frac{1}{\sqrt{\Omega\lambda}}$, over which the system-
environment correlations decay O’Reilly and Olaya-Castro (2014). If we choose
an integration measure much larger than this coherence timescale (blue circles
in Figure (2)e)), we see that the total asymmetry follows a similar behaviour
to $C_{1}(\rho)$. Close to this coherence time (orange squares in Figure
(2)e)), we observe that total asymmetry increases with $\lambda$. We thus see
that correlation asymmetry at early times is able to capture the increased
coherent environmentally induced transport, and thus has concomitant behaviour
to pertinent coherence measures in relevant regimes.

In Figure 3a)-d) we show the second-order correlations and their asymmetries
for multiple polarization angles for the FMO complex, and observe asymmetries
over the expected exciton transport timescales of a few picoseconds in each
case Engel  _et al._ (2007); Kramer  _et al._ (2018) that violate the
classical bound of Eq. (4). For simplicity we use the $x,y$ and $z$ component
of the dipoles. In Figure 3e), g) we show how the total asymmetry over a
timescale of 50 ps depends on the incoherent pump power and reorganization
energy, respectively, noting that pump power is controllable in the laboratory
and thus experimentally accessible. We compare these to the basis independent
coherence in Figure 3f), h), and observe the same trend in each case as the
total asymmetry. The dependence for $C_{1}(\rho)$ we attribute to the
increasing entropy of the state for larger pump power. We thus observe that
close to biological conditions photon correlation asymmetry has similar
behaviour to $C_{1}(\rho)$. Indeed, in the SM we show that assuming a maximal
entropy steady-state and incoherent dynamics, the total asymmetry is zero, as
is $C_{1}(\rho)$.

Conclusions:– In this work we have shown that photon correlation asymmetry
acts as a sensitive probe of quantum coherent energy transport in
multichromophoric systems. We have used exact numerical calculations of a
prototypical photosynthetic light-harvesting system, the FMO complex, to show
that, under realistic electron-phonon interactions, zero delay polarisation-
filtered photon correlations may not suffice to witness quantum coherence in
non-equilibrium steady-states of such photoexcited system. We have shown in
that case correlation time-asymmetries are more reliable reporters of
coherence: asymmetry is present over timescales comparable to exciton transfer
times at room temperature, and violates a bound obtainable from assuming
completely classical dynamics and an incoherent steady state. Further, the
total asymmetry is observed to closely follows the behaviour of the basis
independent coherence.

Analytical arguments based on a simplified model show that whilst correlation
asymmetry may be induced by population-population transfer terms in a general
quantum master equation, for Frenkel-exciton Hamiltonian models with locally
coupled environments the dominant sources of asymmetry are excitonic
coherences.

Our results thus show that time-resolved photon cross-correlations can be used
to probe the quantum coherent nature of transport in single-molecule quantum
optical experiments, which can be realised by deploying current on-chip
technology Davanco  _et al._ (2017). This will not only facilitate the much
sought after experimental insight into the effect that such coherences may
have in the function of photo-active supramolecular complexes, but also the
exploration of such complex quantum emitters (and their phonon environments)
in the development of next-generation quantum photonic technologies Toninelli
_et al._ (2021).

Acknowledgements:– We thank the Gordon and Betty Moore Foundation grant
GBMF8820 and the Engineering and Physical Sciences Research Council (EPSRC UK)
Grant EP/V049011/1 for financial support. We would like to thank Luca
Sapienza, Stefano Vezzoli, Thomas Fay, Graham Fleming, Quanwei Li, K. Birgitta
Whaley and Robert L Cook for insightful discussions. The authors acknowledge
use of the UCL Myriad High Performance Computing Facility.

## References

  * Stotz  _et al._ (2005) J. A. H. Stotz, R. Hey, P. V. Santos, and K. H. Ploog, Nature Materials 4, 585 (2005).
  * Tran  _et al._ (2016) T. T. Tran, K. Bray, M. J. Ford, M. Toth, and I. Aharonovich, Nature Nanotechnology 11, 37 (2016).
  * Reigue  _et al._ (2017) A. Reigue, J. Iles-Smith, F. Lux, L. Monniello, M. Bernard, F. Margaillan, A. Lemaitre, A. Martinez, D. P. S. McCutcheon, J. Mørk, R. Hostein, and V. Voliotis, Phys. Rev. Lett. 118, 233602 (2017).
  * Scholes  _et al._ (2011) G. D. Scholes, G. R. Fleming, A. Olaya-Castro, and R. van Grondelle, Nature Chemistry 3, 763 (2011).
  * Huelga and Plenio (2013) S. Huelga and M. Plenio, Contemporary Physics 54, 181 (2013).
  * Scholes  _et al._ (2017) G. D. Scholes, G. R. Fleming, L. X. Chen, A. Aspuru-Guzik, A. Buchleitner, D. F. Coker, G. S. Engel, R. van Grondelle, A. Ishizaki, D. M. Jonas, J. S. Lundeen, J. K. McCusker, S. Mukamel, J. P. Ogilvie, A. Olaya-Castro, M. A. Ratner, F. C. Spano, K. B. Whaley, and X. Zhu, Nature 543, 647 (2017).
  * Polisseni  _et al._ (2016) C. Polisseni, K. D. Major, S. Boissier, S. Grandi, A. S. Clark, and E. A. Hinds, Opt. Express 24, 5615 (2016).
  * Grandi  _et al._ (2019) S. Grandi, M. P. Nielsen, J. Cambiasso, S. Boissier, K. D. Major, C. Reardon, T. F. Krauss, R. F. Oulton, E. A. Hinds, and A. S. Clark, APL Photonics 4, 086101 (2019).
  * Mattioni  _et al._ (2024) A. Mattioni, J. K. Staab, W. J. A. Blackmore, D. Reta, J. Iles-Smith, A. Nazir, and N. F. Chilton, Nature Communications 15, 485 (2024).
  * Toninelli  _et al._ (2021) C. Toninelli, I. Gerhardt, A. S. Clark, A. Reserbat-Plantey, S. Götzinger, Z. Ristanović, M. Colautti, P. Lombardi, K. D. Major, I. Deperasińska, W. H. Pernice, F. H. L. Koppens, B. Kozankiewicz, A. Gourdon, V. Sandoghdar, and M. Orrit, Nature Materials 20, 1615 (2021).
  * Hanschke  _et al._ (2018) L. Hanschke, K. A. Fischer, S. Appel, D. Lukin, J. Wierzbowski, S. Sun, R. Trivedi, J. Vučković, J. J. Finley, and K. Müller, npj Quantum Information 4, 43 (2018).
  * Bluvstein  _et al._ (2022) D. Bluvstein, H. Levine, G. Semeghini, T. T. Wang, S. Ebadi, M. Kalinowski, A. Keesling, N. Maskara, H. Pichler, M. Greiner, V. Vuletić, and M. D. Lukin, Nature 604, 451 (2022).
  * Fassioli and Olaya-Castro (2010) F. Fassioli and A. Olaya-Castro, New Journal of Physics 12, 085006 (2010).
  * Fassioli  _et al._ (2012) F. Fassioli, A. Olaya-Castro, and G. D. Scholes, The Journal of Physical Chemistry Letters 3, 3136 (2012).
  * Chin  _et al._ (2013) A. W. Chin, J. Prior, R. Rosenbach, F. Caycedo-Soler, S. F. Huelga, and M. B. Plenio, Nature Physics 9, 113 (2013).
  * O’Reilly and Olaya-Castro (2014) E. J. O’Reilly and A. Olaya-Castro, Nature Communications 5 (2014), 10.1038/ncomms4012.
  * Romero  _et al._ (2014) E. Romero, R. Augulis, V. I. Novoderezhkin, M. Ferretti, J. Thieme, D. Zigmantas, and R. van Grondelle, Nature Physics 10, 676 (2014).
  * Romero  _et al._ (2017) E. Romero, V. I. Novoderezhkin, and R. van Grondelle, Nature 543, 355 (2017).
  * Kim  _et al._ (2021) Y. Kim, F. Bertagna, E. M. D’Souza, D. J. Heyes, L. O. Johannissen, E. T. Nery, A. Pantelias, A. Sanchez-Pedreño Jimenez, L. Slocombe, M. G. Spencer, J. Al-Khalili, G. S. Engel, S. Hay, S. M. Hingley-Wilson, K. Jeevaratnam, A. R. Jones, D. R. Kattnig, R. Lewis, M. Sacchi, N. S. Scrutton, S. R. P. Silva, and J. McFadden, Quantum Reports 3, 80 (2021).
  * Nation  _et al._ (2023) C. Nation, V. Notararigo, and A. Olaya-Castro, Phys. Rev. A 108, 063712 (2023).
  * Brixner  _et al._ (2005) T. Brixner, J. Stenger, H. M. Vaswani, M. Cho, R. E. Blankenship, and G. R. Fleming, Nature 434, 625 (2005).
  * Olaya-Castro  _et al._ (2008) A. Olaya-Castro, C. F. Lee, F. F. Olsen, and N. F. Johnson, Phys. Rev. B 78, 085115 (2008).
  * Collini  _et al._ (2010) E. Collini, C. Y. Wong, K. E. Wilk, P. M. G. Curmi, P. Brumer, and G. D. Scholes, Nature 463, 644 (2010).
  * Wilkins and Dattani (2015) D. M. Wilkins and N. S. Dattani, Journal of Chemical Theory and Computation 11, 3411 (2015).
  * Holdaway  _et al._ (2016) D. I. H. Holdaway, E. Collini, and A. Olaya-Castro, The Journal of Chemical Physics 144, 194112 (2016).
  * Holdaway  _et al._ (2017) D. I. H. Holdaway, E. Collini, and A. Olaya-Castro, Opt. Express 25, 6383 (2017).
  * Duan  _et al._ (2017) H. G. Duan, V. I. Prokhorenko, R. J. Cogdell, K. Ashraf, A. L. Stevens, M. Thorwart, and R. J. Miller, Proceedings of the National Academy of Sciences of the United States of America 114, 8493 (2017).
  * Knee  _et al._ (2018) G. C. Knee, M. Marcus, L. D. Smith, and A. Datta, Phys. Rev. A 98, 052328 (2018).
  * Werren  _et al._ (2023) N. Werren, W. Brown, and E. M. Gauger, PRX Energy 2, 013002 (2023).
  * Engel  _et al._ (2007) G. S. Engel, T. R. Calhoun, E. L. Read, T.-K. Ahn, T. Mančal, Y.-C. Cheng, R. E. Blankenship, and G. R. Fleming, Nature 446, 782 (2007).
  * Panitchayangkoon  _et al._ (2010) G. Panitchayangkoon, D. Hayes, K. A. Fransted, J. R. Caram, E. Harel, J. Wen, R. E. Blankenship, and G. S. Engel, Proceedings of the National Academy of Sciences 107, 12766 (2010).
  * 2Dp (2023) Nature Reviews Methods Primers 3, 83 (2023).
  * Glauber (1963) R. J. Glauber, Physical Review 130, 2529 (1963).
  * (34) H. Carmichael,  _Statistical Methods in Quantum Optics 2: Non-Classical Fields_.
  * Aspect  _et al._ (1980) A. Aspect, G. Roger, S. Reynaud, J. Dalibard, and C. Cohen-Tannoudji, Physical Review Letters 45, 617 (1980).
  * Schrama  _et al._ (1992) C. A. Schrama, G. Nienhuis, H. A. Dijkerman, C. Steijsiger, and H. G. M. Heideman, Physical Review A 45, 8045 (1992).
  * Campisi  _et al._ (2011) M. Campisi, P. Hänggi, and P. Talkner, Reviews of Modern Physics 83, 771 (2011).
  * Holdaway  _et al._ (2018) D. I. H. Holdaway, V. Notararigo, and A. Olaya-Castro, Phys. Rev. A 98, 063828 (2018).
  * del Valle  _et al._ (2012) E. del Valle, A. Gonzalez-Tudela, F. P. Laussy, C. Tejedor, and M. J. Hartmann, Phys. Rev. Lett. 109, 183601 (2012).
  * Anderson  _et al._ (2018) M. D. Anderson, S. Tarrago Velez, K. Seibold, H. Flayac, V. Savona, N. Sangouard, and C. Galland, Phys. Rev. Lett. 120, 233601 (2018).
  * Velez  _et al._ (2019) S. T. Velez, K. Seibold, N. Kipfer, M. D. Anderson, V. Sudhir, and C. Galland, Phys. Rev. X 9, 041007 (2019).
  * Sánchez Munõz and Schlawin (2020) C. Sánchez Munõz and F. Schlawin, Physical Review Letters 124, 203601 (2020).
  * Miller  _et al._ (2021) H. J. D. Miller, M. H. Mohammady, M. Perarnau-Llobet, and G. Guarnieri, Physical Review E 103, 052138 (2021).
  * Landi  _et al._ (2023) G. T. Landi, M. J. Kewming, M. T. Mitchison, and P. P. Potts, (2023), arXiv:2303.04270 .
  * Cygorek  _et al._ (2023) M. Cygorek, E. D. Scerri, T. S. Santana, Z. X. Koong, B. D. Gerardot, and E. M. Gauger, Physical Review A 107, 023718 (2023).
  * Downing  _et al._ (2023) C. A. Downing, E. del Valle, and A. I. Fernández-Domínguez, Physical Review A 107, 023717 (2023).
  * Hollars  _et al._ (2003) C. W. Hollars, S. M. Lane, and T. Huser, Chemical Physics Letters 370, 393 (2003).
  * Hedley  _et al._ (2021) G. J. Hedley, T. Schröder, F. Steiner, T. Eder, F. J. Hofmann, S. Bange, D. Laux, S. Höger, P. Tinnefeld, J. M. Lupton, and J. Vogelsang, Nature Communications 12, 1327 (2021).
  * Wientjes  _et al._ (2014) E. Wientjes, J. Renger, A. G. Curto, R. Cogdell, and N. F. Van Hulst, Nature Communications 5, 1 (2014).
  * Kim  _et al._ (2019) Y. Kim, H. L. Puhl, E. Chen, G. H. Taumoefolau, T. A. Nguyen, D. S. Kliger, P. S. Blank, and S. S. Vogel, Biophysical Journal 116, 1918 (2019).
  * Li  _et al._ (2023) Q. Li, K. Orcutt, R. L. Cook, J. Sabines-Chesterking, A. L. Tong, G. S. Schlau-Cohen, X. Zhang, G. R. Fleming, and K. B. Whaley, Nature (2023), 10.1038/s41586-023-06121-5.
  * Iles-Smith  _et al._ (2017) J. Iles-Smith, D. P. S. McCutcheon, J. Mørk, and A. Nazir, Phys. Rev. B 95, 201305 (2017).
  * Humphries  _et al._ (2023) B. S. Humphries, D. Green, M. O. Borgh, and G. A. Jones, Phys. Rev. Lett. 131, 143601 (2023).
  * Nation  _et al._ (2024) C. Nation, V. Notararigo, H. O. Gestsson, L. Sapienza, and A. Olaya-Castro, (2024), arXiv:2403.00857 .
  * Loudon (1980) R. Loudon, Reports on Progress in Physics 43, 913 (1980).
  * Auffèves  _et al._ (2011) A. Auffèves, D. Gerace, S. Portolan, A. Drezet, and M. F. Santos, New Journal of Physics 13, 093020 (2011).
  * (57) Supplemental Material .
  * Streltsov  _et al._ (2017) A. Streltsov, G. Adesso, and M. B. Plenio, Rev. Mod. Phys. 89, 041003 (2017).
  * Mukamel (1995) S. Mukamel,  _Principles of Nonlinear Optics and Spectroscopy_ (New York: Oxford University Press, 1995).
  * (60) M. . Cho,  _Two-Dimensional Optical Spectroscopy (1st ed.)_ (CRC Press.).
  * Ishizaki and Tanimura (2006) A. Ishizaki and Y. Tanimura, The Journal of Chemical Physics 125, 084501 (2006).
  * Jin  _et al._ (2008) J. Jin, X. Zheng, and Y. Yan, The Journal of Chemical Physics 128, 234703 (2008).
  * Ishizaki and Fleming (2009) A. Ishizaki and G. R. Fleming, The Journal of Chemical Physics 130, 234111 (2009).
  * Tanimura (2014) Y. Tanimura, The Journal of Chemical Physics 141, 044114 (2014).
  * Ke and Zhao (2017) Y. Ke and Y. Zhao, The Journal of Chemical Physics 146, 214105 (2017).
  * Fay and Limmer (2022) T. P. Fay and D. T. Limmer, The Journal of Chemical Physics 157, 174104 (2022).
  * Mirkovic  _et al._ (2007) T. Mirkovic, A. B. Doust, J. Kim, K. E. Wilk, C. Curutchet, B. Mennucci, R. Cammi, P. M. Curmi, and G. D. Scholes, Photochemical and Photobiological Sciences 6, 964 (2007).
  * Blau  _et al._ (2018) S. M. Blau, D. I. G. Bennett, C. Kreisbeck, G. D. Scholes, and A. Aspuru-Guzik, Proceedings of the National Academy of Sciences 115 (2018), 10.1073/pnas.1800370115.
  * Cogdell  _et al._ (2006) R. J. Cogdell, A. Gall, and J. Köhler, Quarterly Reviews of Biophysics 39, 227–324 (2006).
  * Cupellini  _et al._ (2020) L. Cupellini, F. Lipparini, and J. Cao, The Journal of Physical Chemistry B 124, 8610 (2020).
  * Kell  _et al._ (2016) A. Kell, R. E. Blankenship, and R. Jankowiak, The Journal of Physical Chemistry A 120, 6146 (2016).
  * Hein  _et al._ (2012a) B. Hein, C. Kreisbeck, T. Kramer, and M. Rodríguez, New Journal of Physics 14, 023018 (2012a).
  * Kramer  _et al._ (2018) T. Kramer, M. Noack, A. Reinefeld, M. Rodríguez, and Y. Zelinskyy, Journal of Computational Chemistry 39, 1779 (2018).
  * Sarovar  _et al._ (2010) M. Sarovar, A. Ishizaki, G. R. Fleming, and K. B. Whaley, Nature Physics 6, 462 (2010).
  * Scully and Zubairy (1997) M. O. Scully and M. S. Zubairy,  _Quantum Optics_ (Cambridge University Press, 1997).
  * Shi  _et al._ (2009) Q. Shi, L. Chen, G. Nan, R.-X. Xu, and Y. Yan, The Journal of Chemical Physics 130, 084105 (2009).
  * Ye  _et al._ (2016) L. Ye, X. Wang, D. Hou, R.-X. Xu, X. Zheng, and Y. Yan, WIREs Computational Molecular Science 6, 608 (2016).
  * Adolphs and Renger (2006) J. Adolphs and T. Renger, Biophysical Journal 91, 2778 (2006).
  * Hein  _et al._ (2012b) B. Hein, C. Kreisbeck, T. Kramer, and M. Rodríguez, New Journal of Physics 14, 023018 (2012b).
  * Gelzinis  _et al._ (2015) A. Gelzinis, D. Abramavicius, and L. Valkunas, The Journal of Chemical Physics 142, 154107 (2015).
  * Korzh  _et al._ (2020) B. Korzh, Q.-Y. Zhao, J. P. Allmaras, S. Frasca, T. M. Autry, E. A. Bersin, A. D. Beyer, R. M. Briggs, B. Bumble, M. Colangelo, G. M. Crouch, A. E. Dane, T. Gerrits, A. E. Lita, F. Marsili, G. Moody, C. Peña, E. Ramirez, J. D. Rezac, N. Sinclair, M. J. Stevens, A. E. Velasco, V. B. Verma, E. E. Wollman, S. Xie, D. Zhu, P. D. Hale, M. Spiropulu, K. L. Silverman, R. P. Mirin, S. W. Nam, A. G. Kozorezov, M. D. Shaw, and K. K. Berggren, Nature Photonics 14, 250 (2020).
  * Ma  _et al._ (2019) Z.-H. Ma, J. Cui, Z. Cao, S.-M. Fei, V. Vedral, T. Byrnes, and C. Radhakrishnan, Europhysics Letters 125, 50005 (2019).
  * Le and Olaya-Castro (2020) T. P. Le and A. Olaya-Castro, (2020), arXiv:2011.15016 .
  * Curti  _et al._ (2023) M. Curti, V. Maffeis, L. G. Teixeira Alves Duarte, S. Shareef, L. X. Hallado, C. Curutchet, and E. Romero, Protein Science 32, e4579 (2023), https://onlinelibrary.wiley.com/doi/pdf/10.1002/pro.4579 .
  * Sánchez-Pedreño Jiménez  _et al._ (2023) A. Sánchez-Pedreño Jiménez, H. L. Puhl, S. S. Vogel, and Y. Kim, Phys. Chem. Chem. Phys. 25, 19532 (2023).
  * Davanco  _et al._ (2017) M. Davanco, J. Liu, L. Sapienza, C.-Z. Zhang, J. V. De Miranda Cardoso, V. Verma, R. Mirin, S. W. Nam, L. Liu, and K. Srinivasan, Nature Communications 8, 889 (2017).
  * Chin  _et al._ (2010) A. W. Chin, A. Datta, F. Caruso, S. F. Huelga, and M. B. Plenio, New Journal of Physics 12, 065002 (2010).
  * Chin  _et al._ (2012) A. W. Chin, S. F. Huelga, and M. B. Plenio, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 370, 3638 (2012).
  * Ishizaki and Tanimura (2005) A. Ishizaki and Y. Tanimura, J. Phys. Soc. Jpn. 74, 3131 (2005).
  * Nakajima (1958) S. Nakajima, Progress of Theoretical Physics 20, 948 (1958).
  * Zwanzig (1960) R. Zwanzig, The Journal of Chemical Physics 33, 1338 (1960).
  * Tokuyama and Mori (1976) M. Tokuyama and H. Mori, Progress of Theoretical Physics 55, 411 (1976).
  * Breuer and Petruccione (2002) H. P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press, Great Clarendon Street, 2002).
  * Gardiner and Zoller (2000) C. W. Gardiner and P. Zoller,  _Quantum Noise_ , 2nd ed., edited by H. Haken (Springer, 2000).
  * Brian and Sun (2021) D. Brian and X. Sun, Chinese Journal of Chemical Physics 34, 497 (2021).

Supplemental Material:
Photon correlation time-asymmetry and dynamical coherence in multi-
chromophoric systems

Charlie Nation, Hallmann Óskar Gestsson, and Alexandra Olaya-Castro

1Department of Physics and Astronomy, University College London, Gower Street,
WC1E 6BT, London, United Kingdom

##  I Second excitation manifold of the Frenkel exciton model

In this section we describe the approach to define the doubly excited
manifold, for which we follow Ref. Hein  _et al._ (2012b) and references
therein. We define the matrix elements of the second excited subspace
Hamiltonian $H^{(2)}$ in terms of the singly excited subspace $H^{(1)}$ as,

$\displaystyle\langle mn|H^{(2)}|mn\rangle=\langle m|H^{(1)}|m\rangle+\langle
n|H^{(1)}|n\rangle$ (S1)

for diagonal elements, and

$\displaystyle\langle kl|H^{(2)}$
$\displaystyle|mn\rangle=\delta_{km}(1-\delta_{ln})\langle l|H^{(1)}|n\rangle$
(S2) $\displaystyle+\delta_{kn}(1-\delta_{lm})\langle
l|H^{(1)}|m\rangle+\delta_{lm}(1-\delta_{kn})\langle k|H^{(1)}|n\rangle$
$\displaystyle+\delta_{ln}(1-\delta_{km})\langle k|H^{(1)}|m\rangle.$

It will be illustrative below to initially analyse the case where doubly
excited states $|mm\rangle$ are included, which amounts to removing the
$(1-\delta_{mn})$ factors in each term above.

The full Hamiltonian including the ground state is then of the form

$\displaystyle H=\begin{pmatrix}0&0&0\\\ 0&H^{(1)}&0\\\
0&0&H^{(2)}\end{pmatrix}.$ (S3)

Then, for an operator $O$ in the second excited subspace we have

$\displaystyle O_{ij,ik}=\langle ij|O|ik\rangle=\langle j|O|k\rangle$ (S4)

for $j\neq k$. If we then decompose a dynamical map $M$ into a sum over
contributions from individual operators $O$ such that $M=\sum_{O}M_{O}$ with
$M_{O}[\cdot]=O_{L}\cdot O_{R}$, where each operator $O_{L/R}$ is defined on
the second excited state in the manner above. We then have that

$\displaystyle(M_{O})_{ij\to ik}$ $\displaystyle=\langle
ik|M_{O}[|ij\rangle\langle ij|]|ik\rangle$ $\displaystyle=\langle
ik|O_{L}|ij\rangle\langle ij|O_{R}|ik\rangle$ (S5) $\displaystyle=\langle
k|O_{L}|j\rangle\langle j|O_{R}|k\rangle$ $\displaystyle=(M_{O})_{j\to k}.$

##  II Photon correlations for classical maps

In this section we assert classicality by i) enforcing that the steady-state
has no coherences in some basis $|\alpha\rangle$, ii) the generator of
dynamics does not create coherences from populations, that is
$\mathcal{E}(t)[|\alpha\rangle\langle\alpha|]=\sum_{\beta}W_{\alpha\to\beta}(t)|\beta\rangle\langle\beta|$,
where $W_{\alpha\to\beta}$ is understood as the transition probability from
state $|\alpha\rangle$ to state $|\beta\rangle$. If this is true of a system,
then there exists a basis in which a classical map describes the complete
dynamics. Here this basis is taken as the excitonic basis diagonalising $H$,
in which we are interested in probing steady-state coherences, as well as
coherences induced by dynamics after photon emission from the steady-state.

In order to calculate the second-order photon correlation functions we must
express the dipole operators in terms of the two-excitation subspace of a
Frenkel-Exciton model. The first simplifying assumption we make in describing
photon emission processes is the following:

$\displaystyle\hat{\mu}_{a}$
$\displaystyle=\sum_{\alpha}(\mu_{a\alpha}|0\rangle\langle\alpha|+\sum_{\beta>\alpha}(\mu_{a\beta}|\alpha\rangle\langle\alpha\beta|+\mu_{a\alpha}|\beta\rangle\langle\alpha\beta|))$
$\displaystyle=\hat{\mu}_{a}^{(1)}+\hat{\mu}_{a}^{(2)},$ (S6)

where $\hat{\mu}_{a}^{(n)}$ acts on the $n$ excitation manifold to map onto
that of $n-1$ excitations for $n\in[1,2]$, and
$|\alpha\rangle=\sum_{m}c_{m\alpha}|m\rangle$ is the excitonic state. The
state $|\alpha\beta\rangle$ denotes an excitation on both excitons $\alpha$
and $\beta$. We note that in general, we have $\langle
0|\mu^{(1)}_{a}|\alpha\rangle=\sum_{m}\mu_{m,a}\langle m|\alpha\rangle$ and
$\langle\alpha|\mu^{(2)}|\beta\gamma\rangle=\sum_{m}\sum_{n\neq
m}\mu_{m}\langle\alpha|m\rangle\langle mn|\beta\gamma\rangle$, and thus the
approximation ignores emission processes from the doubly excited subspace of
the form $|\alpha,\beta\rangle\to|\gamma\rangle$ for
$\gamma\neq\alpha,\,\beta$.

The steady-state in the doubly excited space is

$\displaystyle\rho=\sum_{\begin{subarray}{c}\alpha\beta\alpha^{\prime}\beta^{\prime}=0\\\
\alpha>\beta,\alpha\prime>\beta^{\prime}\end{subarray}}\rho_{\alpha\beta\alpha^{\prime}\beta^{\prime}}|\alpha\beta\rangle\langle\alpha^{\prime}\beta^{\prime}|.$
(S7)

An incoherent state has the additional assertion that
$\alpha=\alpha^{\prime}$, $\beta=\beta^{\prime}$. We assume that $\rho$ is the
steady-state of the dynamical generator $\mathcal{E}(t)$, such that
$\mathcal{E}(t)[\rho]=\rho$. We may understand the second-order photon
correlation functions of Eq. (2) as a sequence of two maps
$\tilde{M}_{a}[\cdot]=\frac{1}{S^{(1)}_{a}}\hat{\mu}_{a}\cdot\hat{\mu}_{a}^{\dagger}=\frac{1}{S^{(1)}_{a}}M_{a}[\cdot]$
with $S^{(1)}_{a}=\langle\hat{\mu}^{\dagger}_{a}\hat{\mu}_{a}\rangle$. Up to a
normalisation factor of $\frac{1}{S^{(1)}_{a}S^{(1)}_{b}}$ the two-point
photon correlation function of Eq. (2) is obtained via

$\displaystyle
S_{ab}^{(2)}(\tau):=\operatorname{Tr}[M_{b}[\mathcal{E}(\tau)[M_{a}[\rho]]]].$
(S8)

We begin by writing an incoherent state in a convenient form:

$\displaystyle\rho_{inc}=p^{(0)}|0\rangle\langle
0|+\sum_{\alpha=1}^{N}p^{(1)}_{\alpha}|\alpha\rangle\langle\alpha|+\sum_{\begin{subarray}{c}\alpha\beta\\\
\alpha>\beta\end{subarray}}p^{(2)}_{\alpha\beta}|\alpha\beta\rangle\langle\alpha\beta|$
$\displaystyle=\rho^{(0)}+\rho^{(1)}+\rho^{(2)},$ (S9)

where we have split the contributions into the zero excitation, single
excitation, and doubly excited subspaces. We can simply ignore the
contribution of the zero excited state to $M_{a}[\rho]$, as this only has
terms $\propto\langle 0|\alpha\rangle$ for $\alpha\neq 0$, and thus
$,M_{a}[\rho^{(0)}]=0$. We then have

$\displaystyle M_{a}[\rho_{inc}]$
$\displaystyle=\sum_{\alpha}|\mu_{a\alpha}|^{2}p^{(1)}_{\alpha}|0\rangle\langle
0|$ (S10) $\displaystyle+\sum_{\begin{subarray}{c}\alpha\beta\\\
\alpha>\beta\end{subarray}}p^{(2)}_{\alpha\beta}(|\mu_{a\beta}|^{2}|\alpha\rangle\langle\alpha|+|\mu_{a\alpha}|^{2}|\beta\rangle\langle\beta|$
$\displaystyle\qquad\qquad+\mu_{a\beta}\mu_{a\alpha}^{*}|\alpha\rangle\langle\beta|+\mu_{a\alpha}\mu_{a\beta}^{*}|\beta\rangle\langle\alpha|).$

We now require the action of $\mathcal{E}$ onto $M_{a}[\rho_{inc}]$, which can
be obtained by defining $\mathcal{E}$ as an incoherent map and thus having the
properties:

$\displaystyle\mathcal{E}_{inc}(\tau)[|\alpha\beta\rangle\langle\alpha\beta|]=\sum_{\begin{subarray}{c}\alpha^{\prime}\beta^{\prime}\\\
\alpha^{\prime}>\beta^{\prime}\end{subarray}}W_{\alpha\beta\to\alpha^{\prime}\beta^{\prime}}(\tau)|\alpha^{\prime}\beta^{\prime}\rangle\langle\alpha^{\prime}\beta^{\prime}|$
(S11)
$\displaystyle\mathcal{E}_{inc}(\tau)[|\alpha\beta\rangle\langle\alpha^{\prime}\beta^{\prime}|]=0\quad\textrm{if}\quad\alpha,\beta\neq\alpha^{\prime},\beta^{\prime}.$
(S12)

In this sense we have asserted that both the steady-state and its dynamical
generator are ‘classical’, in that they do not contain or generate coherences.
We have thus effectively enforced dynamics under a Pauli master equation. We
then have

$\displaystyle\mathcal{E}_{inc}$
$\displaystyle(\tau)[M_{a}[\rho_{inc}]]=|\hat{\mu}_{a\alpha}|^{2}p^{(1)}_{\alpha}\mathcal{E}_{inc}(\tau)[|0\rangle\langle
0|]$ (S13) $\displaystyle+\sum_{\begin{subarray}{c}\alpha\beta\\\
\alpha>\beta\end{subarray}}p^{(2)}_{\alpha\beta}(|\hat{\mu}_{a\alpha}|^{2}\mathcal{E}_{inc}(\tau)[|\beta\rangle\langle\beta|]+|\hat{\mu}_{a\beta}|^{2}\mathcal{E}_{inc}(\tau)[|\alpha\rangle\langle\alpha|])$
$\displaystyle=\sum_{\alpha=1}^{N}\sum_{\begin{subarray}{c}\alpha^{\prime}\beta^{\prime}=0\\\
\alpha^{\prime}>\beta^{\prime}\end{subarray}}|\hat{\mu}_{a\alpha}|^{2}p^{(1)}_{\alpha}W_{00\to\alpha^{\prime}\beta^{\prime}}(\tau)|\alpha^{\prime}\beta^{\prime}\rangle\langle\alpha^{\prime}\beta^{\prime}|$
$\displaystyle+\sum_{\begin{subarray}{c}\alpha\beta=1\\\
\alpha>\beta\end{subarray}}\sum_{\begin{subarray}{c}\alpha^{\prime}\beta^{\prime}=0\\\
\alpha^{\prime}>\beta^{\prime}\end{subarray}}p^{(2)}_{\alpha\beta}(|\hat{\mu}_{a\alpha}|^{2}W_{\beta
0\to\alpha^{\prime}\beta^{\prime}}(\tau)$
$\displaystyle\qquad\qquad+|\hat{\mu}_{a\beta}|^{2}W_{\alpha
0\to\alpha^{\prime}\beta^{\prime}}(\tau))|\alpha^{\prime}\beta^{\prime}\rangle\langle\alpha^{\prime}\beta^{\prime}.|$

Here the first term corresponds to an initial de-excitation from the first
excited state to the ground state by emission of photon $a$ followed by
evolution from $|00\rangle\to|\alpha^{\prime}\beta^{\prime}\rangle$, including
the possibility of no change in the state. The second term, corresponds to
photon emission from the doubly excited state to the singly excited state,
followed by evolution under $\mathcal{E}_{inc}(\tau)$ to any other state.

In order to calculate the second-order correlation function, we then wish to
find the application of the map $M_{b}$ onto the above result, for which we
find it convenient to re-express the above result into 4 terms corresponding
to the transitions that are caused between manifolds by the action of
$\mathcal{E}_{inc}(\tau)$. We denote these as terms containing
$W^{(01)}_{00\to\alpha^{\prime}0}$, $W^{(11)}_{\alpha 0\to\alpha^{\prime}0}$,
$W^{(02)}_{00\to\alpha^{\prime}\beta^{\prime}}$, and $W^{(12)}_{\alpha
0\to\alpha^{\prime}\beta^{\prime}}$, corresponding to transitions from ground
to first excited states, evolution within the first excited state, ground to
second excited states, and transitions from the first excited states to the
second excited states, respectively.

We note that typically only the $W^{(11)}_{\alpha 0\to\alpha^{\prime}0}$ term
will contribute to early time asymmetries, as the transition rate between
states of differing excitation number is $\sim 1000$ times smaller than the
transition rate between states in the singly excited manifold, and we are
interested in the limit of weak incoherent pumping, where the relevant pumping
rates are significantly smaller than the decay rates. Thus, we ignore terms
that involve transitions between excited state manifolds, as the dominant
contribution is that of exciton transfer by $W^{(11)}_{\alpha
0\to\alpha^{\prime}0}$. We additionally exclude any terms that leave the state
in the ground state, as these are destroyed by the action of $M_{b}$. Thus,
the relevant terms are

$\displaystyle\mathcal{E}_{inc}(\tau)$ $\displaystyle[M_{a}[\rho_{inc}]]$
(S14) $\displaystyle\approx\sum_{\begin{subarray}{c}\alpha\beta=1\\\
\alpha>\beta\end{subarray}}\sum_{\alpha^{\prime}=1}p^{(2)}_{\alpha\beta}(|\hat{\mu}_{a\alpha}|^{2}W^{(11)}_{\beta
0\to\alpha^{\prime}0}(\tau)$
$\displaystyle\qquad\quad+|\hat{\mu}_{a\beta}|^{2}W^{(11)}_{\alpha
0\to\alpha^{\prime}0}(\tau))|\alpha^{\prime}\rangle\langle\alpha^{\prime}|$
$\displaystyle=\sum_{\begin{subarray}{c}\alpha\beta=1\end{subarray}}\sum_{\alpha^{\prime}=1}p^{(2)}_{\alpha\beta}|\hat{\mu}_{a\alpha}|^{2}W^{(11)}_{\beta
0\to\alpha^{\prime}0}(\tau)|\alpha^{\prime}\rangle\langle\alpha^{\prime}|$

where we have extended the definition
$p_{\alpha\beta}^{(2)}=p_{\beta\alpha}^{(2)}$ and $p_{\alpha\alpha}^{(2)}=0$.
The above form simplifies the calculation of the action of $M_{b}$, as we are
interested only in diagonal terms (due to the trace being taken after
$M_{b}$), which allows us to write

$\displaystyle\operatorname{Tr}[M_{b}[$
$\displaystyle\mathcal{E}_{inc}(\tau)[M_{a}[\rho_{inc}]]]]$
$\displaystyle\approx\operatorname{Tr}[M_{b}^{(1\to
0)}[\mathcal{E}_{inc}^{(1\to 1)}(\tau)[M_{a}^{(2\to 1)}[\rho_{inc}]]]]$ (S15)
$\displaystyle=\operatorname{Tr}[\hat{\mu}_{b}^{(1)}\sum_{ij}\sum_{f=1}p^{(2)}_{\alpha\beta}|\hat{\mu}_{a\alpha}|^{2}W^{(11)}_{\beta
0\to\alpha^{\prime}0}(\tau)|\alpha^{\prime}\rangle\langle\alpha^{\prime}|\hat{\mu}_{b}^{(1)\dagger}]$
$\displaystyle=\sum_{\alpha\beta}\sum_{\alpha^{\prime}}p_{\alpha\beta}^{(2)}|\hat{\mu}_{b\alpha^{\prime}}|^{2}|\hat{\mu}_{a\alpha}|^{2}W^{(11)}_{\beta\to\alpha^{\prime}}(\tau)$
$\displaystyle=\sum_{\alpha,\alpha^{\prime}}|\hat{\mu}_{b\alpha^{\prime}}|^{2}|\hat{\mu}_{a\alpha}|^{2}\sum_{\beta}p_{\alpha\beta}^{(2)}W_{\beta\to\alpha^{\prime}}(\tau),$

where the superscript $(n\to n^{\prime})$ is used on maps in the first line to
represent the terms of the map which take states from the $n$-excitation
subspace to the $n^{\prime}$ excitation subspace. This captures the key
approximation implied by the weak pumping and decay condition.

Notably, this expression allows for the total asymmetry to be bounded, as we
can write this asymmetry as

$\displaystyle A$
$\displaystyle\propto|\sum_{\alpha,\alpha^{\prime}}(|\hat{\mu}_{b\alpha^{\prime}}|^{2}|\hat{\mu}_{a\alpha}|^{2}-|\hat{\mu}_{a\alpha^{\prime}}|^{2}|\hat{\mu}_{b\alpha}|^{2})\sum_{\beta}p_{\alpha\beta}^{(2)}W_{\beta\to\alpha^{\prime}}(\tau)|$
$\displaystyle\leq\sum_{\alpha,\alpha^{\prime}}|(|\hat{\mu}_{b\alpha^{\prime}}|^{2}|\hat{\mu}_{a\alpha}|^{2}-|\hat{\mu}_{a\alpha^{\prime}}|^{2}|\hat{\mu}_{b\alpha}|^{2})|\sum_{\beta}p_{\alpha\beta}^{(2)}W_{\beta\to\alpha^{\prime}}(\tau),$
(S16)

as $p^{(2)}_{\alpha\beta}>0$ and $W_{\alpha\to\beta}>0$. Then using
$\sum_{\beta}W_{\beta\to\alpha}=\langle\alpha|e^{\mathcal{L^{\dagger}\tau}}[\sum_{\beta}|\beta\rangle\langle\beta|]\alpha\rangle=\langle\alpha|e^{\mathcal{L^{\dagger}\tau}}[\mathbb{1}]|\alpha\rangle=1$,
where we have used that the adjoint map is unital, which can be seen from the
trace preserving condition on the Kraus decomposition of the dynamical map.
The Cauchy-Schwarz inequality implies that
$\sum_{\beta}p_{\alpha\beta}^{(2)}W_{\beta\to\alpha^{\prime}}(\tau)\leq\sum_{\beta}p_{\alpha\beta}^{(2)}\leq
p^{(2)}$, where $p^{(2)}$ is the total occupation of the doubly excited state.
This leads us to,

$\displaystyle
A\leq\frac{\sum_{\alpha,\alpha^{\prime}}|(|\hat{\mu}_{b\alpha^{\prime}}|^{2}|\hat{\mu}_{a\alpha}|^{2}-|\hat{\mu}_{a\alpha^{\prime}}|^{2}|\hat{\mu}_{b\alpha}|^{2})|}{S^{(1)}_{a}S^{(1)}_{b}}p^{(2)}.$
(S17)

We thus have a classical bound on correlation asymmetry.

###  II.1 Examples of correlation symmetry

We here analyse the case of $N$ multi-level emitters - that is, a classical
model as above however where $p_{\alpha\alpha}\neq 0$. In this case, as a
consequence of the definition in Section I of the doubly excited subspace we
have that
$\sum_{\beta}p_{\alpha\beta}^{(2)}W_{\beta\to\alpha^{\prime}}(\tau)=\sum_{\beta}p_{\alpha\beta}^{(2)}\langle\alpha^{\prime}|\mathcal{E}_{inc}(\tau)[|\beta\rangle\langle\beta|]|\alpha^{\prime}\rangle=\langle\alpha^{\prime}|\mathcal{E}_{inc}(\tau)[\sum_{\beta}p_{\alpha\beta}^{(2)}|\beta\rangle\langle\beta|]|\alpha^{\prime}\rangle=p_{\alpha\alpha^{\prime}}$.
The final equality here is seen noting that
$p^{(2)}_{\alpha\beta}=p_{\alpha}p_{\beta}$ for a classical emitter, and thus
$\sum_{\beta}p_{\alpha\beta}^{(2)}|\beta\rangle\langle\beta|=p_{\alpha}\sum_{\beta}p_{\beta}|\beta\rangle\langle\beta|$,
which is proportional to the steady-state of the dynamical map. We then have

$\displaystyle\sum_{\alpha,\alpha^{\prime}}|\hat{\mu}_{b\alpha^{\prime}}|^{2}|\hat{\mu}_{a\alpha}|^{2}\sum_{\beta}p_{\alpha\beta}^{(2)}W_{\beta\to\alpha^{\prime}}(\tau)=\sum_{\alpha,\alpha^{\prime}}|\hat{\mu}_{b\alpha^{\prime}}|^{2}|\hat{\mu}_{a\alpha}|^{2}p_{\alpha\alpha^{\prime}}^{(2)},$
(S18)

which, crucially, is symmetric on exchange $a\leftrightarrow b$, and hence
$S_{ab}^{(2)}=S_{ba}^{(2)}$. We thus recover symmetry of correlations in the
case of classical multi-level emitters.

In this example we have seen that the relation $p^{(2)}_{\alpha\beta}\propto
p_{\alpha}p_{\beta}$ leads to the symmetry of photon correlations. This
relation is true for example for thermal equilibrium states
$\rho_{\beta}=Z^{-1}\sum_{\alpha}e^{-\beta
E_{\alpha}}|\alpha\rangle\langle\alpha|$, with $Z=\sum_{\alpha}e^{-\beta
E_{\alpha}}$.

The above considered case of a state which leads to correlation symmetry leads
us to note that for a fully classical system with the assumption of non-
invasive measurability, we trivially expect fully time-symmetric correlations
of emission. That is, if the initial measurement $M_{a}$ does not influence
the underlying probability distribution: $M_{a}[\rho]=m_{a}\rho$ for some
value $m_{a}$.

We can further see from the fact that for weak pumping and decay
$\sum_{\beta}W_{\beta\to\alpha}\approx 1$, if the state on the doubly excited
manifold is a maximal entropy state $\frac{1}{d^{(2)}}\mathbb{1}$, where
$d^{(2)}$ is the dimension on the doubly excited subspace, we recover once
more correlation symmetry. This explains why the asymmetry behaves in a
similar manner to the basis independent coherence in the limits of high
entropy where each are expected to be small.

###  II.2 Diagonal Map

In the above treatment of an incoherent map we completely ignored the presence
of coherences, which is a drastic simplification. We can see from Eq. (S10)
that even for an initial state with no coherences, the initial photon emission
may lead to a transient state with non-zero coherences. We thus extend the
above analysis to the case of a diagonal map, $\mathcal{E}_{D}(\tau)$, where
populations map to populations, and coherences to cohernces, such that
populations and coherences are completely decoupled. We further saw that the
only contribution of the dynamical map in the limit of weak pumping and decay
is to the single excited state manifold, thus, we write

$\displaystyle\mathcal{E}_{D}(\tau)[|\alpha\rangle\langle\alpha|]=\sum_{\beta}W_{\alpha\to\beta}(\tau)|\beta\rangle\langle\beta|$
(S19) $\displaystyle\mathcal{E}_{D}(\tau)[|\alpha\rangle\langle\beta|]\,=$
$\displaystyle\qquad\qquad\sum_{\alpha^{\prime}\beta^{\prime}}W^{(C-C)}_{\alpha\beta\to\alpha^{\prime}\beta^{\prime}}(\tau)|\alpha^{\prime}\rangle\langle\beta^{\prime}|.$

In the above section we found the term corresponding to the evolution of
populations under this map, here we see that the evolution of coherences may
also contribute. We obtain, using $M_{b}^{(1\to
0)}[|\alpha\rangle\langle\beta|]=\mu_{b\alpha}\mu_{b\beta}^{*}|0\rangle\langle
0|$, that the contribution of this coherence term to photon cross correlations
is

$\displaystyle\sum_{\begin{subarray}{c}\alpha\beta\\\
\alpha\neq\beta\end{subarray}}\sum_{\begin{subarray}{c}\alpha^{\prime}\beta^{\prime}\\\
\alpha^{\prime}\neq\beta^{\prime}\end{subarray}}p_{\alpha\beta}^{(2)}\mu_{a\alpha}\mu_{a\beta}^{*}\mu_{b\alpha^{\prime}}\mu_{b\beta^{\prime}}^{*}W^{(C-C)}_{\alpha\beta\to\alpha^{\prime}\beta^{\prime}},$
(S20)

and thus for a diagonal map we have:

$\displaystyle\operatorname{Tr}[M_{b}[$
$\displaystyle\mathcal{E}_{D}(\tau)[M_{a}[\rho_{inc}]]]]=\sum_{\alpha,\alpha^{\prime}}|\hat{\mu}_{b\alpha^{\prime}}|^{2}|\hat{\mu}_{a\alpha}|^{2}\sum_{\beta}p_{\alpha\beta}^{(2)}W_{\beta\to\alpha^{\prime}}(\tau)$
$\displaystyle+\sum_{\begin{subarray}{c}\alpha\beta\\\
\alpha\neq\beta\end{subarray}}\sum_{\begin{subarray}{c}\alpha^{\prime}\beta^{\prime}\\\
\alpha^{\prime}\neq\beta^{\prime}\end{subarray}}p_{\alpha\beta}^{(2)}\mu_{a\alpha}\mu_{a\beta}^{*}\mu_{b\alpha^{\prime}}\mu_{b\beta^{\prime}}^{*}W^{(C-C)}_{\alpha\beta\to\alpha^{\prime}\beta^{\prime}}.$
(S21)

Now, we can see that the coherence-coherence term contributes the asymmetry:

$\displaystyle\sum_{\begin{subarray}{c}\alpha\beta\\\
\alpha\neq\beta\end{subarray}}\sum_{\begin{subarray}{c}\alpha^{\prime}\beta^{\prime}\\\
\alpha^{\prime}\neq\beta^{\prime}\end{subarray}}p_{\alpha\beta}^{(2)}\mu_{a\alpha}\mu_{a\beta}^{*}\mu_{b\alpha^{\prime}}\mu_{b\beta^{\prime}}^{*}W^{(C-C)}_{\alpha\beta\to\alpha^{\prime}\beta^{\prime}}$
(S22) $\displaystyle\qquad\qquad-\sum_{\begin{subarray}{c}\alpha\beta\\\
\alpha\neq\beta\end{subarray}}\sum_{\begin{subarray}{c}\alpha^{\prime}\beta^{\prime}\\\
\alpha^{\prime}\neq\beta^{\prime}\end{subarray}}p_{\alpha\beta}^{(2)}\mu_{b\alpha}\mu_{b\beta}^{*}\mu_{a\alpha^{\prime}}\mu_{a\beta^{\prime}}^{*}W^{(C-C)}_{\alpha\beta\to\alpha^{\prime}\beta^{\prime}}$
$\displaystyle=\sum_{\begin{subarray}{c}\alpha\beta\\\
\alpha\neq\beta\end{subarray}}\sum_{\begin{subarray}{c}\alpha^{\prime}\beta^{\prime}\\\
\alpha^{\prime}\neq\beta^{\prime}\end{subarray}}p_{\alpha\beta}^{(2)}(\mu_{a\alpha}\mu_{a\beta}^{*}\mu_{b\alpha^{\prime}}\mu_{b\beta^{\prime}}^{*}-\mu_{b\alpha}\mu_{b\beta}^{*}\mu_{a\alpha^{\prime}}\mu_{a\beta^{\prime}}^{*})W^{(C-C)}_{\alpha\beta\to\alpha^{\prime}\beta^{\prime}}.$
(S23)

From which we recover that auto-correlations ($a=b$) are trivially time-
symmetric. The above can be re-written by instead relabelling the summed
indices:

$\displaystyle\sum_{\begin{subarray}{c}\alpha\beta\\\
\alpha\neq\beta\end{subarray}}\sum_{\begin{subarray}{c}\alpha^{\prime}\beta^{\prime}\\\
\alpha^{\prime}\neq\beta^{\prime}\end{subarray}}\mu_{a\alpha}\mu_{a\beta}^{*}\mu_{b\alpha^{\prime}}\mu_{b\beta^{\prime}}^{*}(p_{\alpha\beta}^{(2)}W^{(C-C)}_{\alpha\beta\to\alpha^{\prime}\beta^{\prime}}-p_{\alpha^{\prime}\beta^{\prime}}^{(2)}W^{(C-C)}_{\alpha^{\prime}\beta^{\prime}\to\alpha\beta}).$
(S24)

##  III Frenkel Exciton with time-dependent pure-dephasing

Here we show a specific example of an $N$ site Frenkel exciton model with a
pure-dephasing environment, as studied extensively in Chin  _et al._ (2010,
2012) in the context of environmentally induced transport processes.
Concretely, we have

$\displaystyle\mathcal{L}[\rho]=-i[H(t),\rho]+\sum_{i}\frac{\gamma_{i}(t)}{2}[2|i\rangle\langle
i|\rho|i\rangle\langle i|-\\{|i\rangle\langle i|,\rho\\}],$ (S25)

where $|i\rangle$ label the sites of the model. In the following, for
simplicity, we drop the explicit time dependence on dephasing rates
$\gamma_{i}$ and the Hamiltonian $H(t)$, however note that no assumptions
regarding their time dependence are made.

We first wish to write this in the excitonic basis, in order to explore the
role of excitonic coherences in the asymmetry of photon correlations. To do so
we define $C_{\alpha\beta}^{(i)}:=\langle\alpha|i\rangle\langle
i|\beta\rangle$, and the master equation becomes

$\displaystyle\mathcal{L}[\rho]=\bigg{(}$
$\displaystyle-i\sum_{\alpha\beta}(E_{\alpha}\rho_{\alpha\beta}-E_{\beta}\rho_{\alpha\beta})$
$\displaystyle-\sum_{i}\frac{\gamma_{i}}{2}\sum_{\alpha\beta\alpha^{\prime}}(C_{\alpha\alpha^{\prime}}^{(i)}\rho_{\alpha^{\prime}\beta}+C_{\alpha^{\prime}\beta}^{(i)}\rho_{\alpha\alpha^{\prime}})$
$\displaystyle+\sum_{i}\gamma_{i}\sum_{\alpha\beta\alpha^{\prime}\beta^{\prime}}C_{\alpha\beta^{\prime}}^{(i)}C_{\alpha^{\prime}\beta}^{(i)}\rho_{\beta^{\prime}\alpha^{\prime}}\bigg{)}|\alpha\rangle\langle\beta|.$
(S26)

The first two terms of the above master equation can be written in terms of a
non Hermitian Hamiltonian
$H_{NH}=H-i\sum_{i}\frac{\gamma_{i}}{2}|i\rangle\langle
i|=H-i\sum_{i}\frac{\gamma_{i}}{2}\sum_{\alpha\beta}C^{(i)}_{\alpha\beta}|\alpha\rangle\langle\beta|$,
by defining $\mathcal{L}_{NH}[\cdot]=-i(H_{NH}\cdot-\cdot H_{NH}^{\dagger})$.
The third term is that due to quantum jumps, which we write as
$\mathcal{L}_{J}$.

We have seen in the above that the asymmetry of photon correlations can be
analysed in four steps. 1) Initial measurement on steady-state
$\rho^{(a)}=M_{a}^{(2\to 1)}[\rho]$. 2) Time evolution under $\mathcal{L}$. 3)
Second measurement $M_{b}^{(1\to 0)}[\rho^{(a)}(t)]$. 4) Analysis of symmetry
under interchange of measurement order via that of indexes $a,b$. In order to
determine the symmetry, for simplicity, in step 2) rather than calculating the
full time evolution we can calculate the application of the dynamical
generator, or Liouvillian $\mathcal{L}$.

We first write the steady-state of the system, which we assume to be
incoherent, as

$\displaystyle\rho_{inc}$ $\displaystyle=p^{(0)}|0\rangle\langle
0|+\sum_{\alpha=1}^{N}p^{(1)}_{\alpha}|\alpha\rangle\langle\alpha|+\sum_{\begin{subarray}{c}\alpha\beta\\\
\alpha>\beta\end{subarray}}p^{(2)}_{\alpha\beta}|\alpha\beta\rangle\langle\alpha\beta|$
$\displaystyle=\rho^{(0)}+\rho^{(1)}+\rho^{(2)}.$ (S27)

Beginning with step 1), then, we apply the measurement via emission from the
doubly excited subspace to the single excited subspace:

$\displaystyle\rho^{(a)}_{inc}=\sum_{\begin{subarray}{c}\alpha\beta\\\
\alpha\neq\beta\end{subarray}}p^{(2)}_{\alpha\beta}(|\mu_{a\alpha}|^{2}|\beta\rangle\langle\beta|+\mu_{a\alpha}\mu_{a\beta}^{*}|\alpha\rangle\langle\beta|).$
(S28)

We thus observe that the initial emission process creates coherences from an
initial incoherent state.

For step 2) we then calculate the action of the Liouvillian $\mathcal{L}$,
which we split as above into two parts. First, the non-Hermitian Hamiltonian
part, which acts on matrix elements $|\alpha\rangle\langle\beta|$ to obtain

$\displaystyle\mathcal{L}_{NH}[|\alpha\rangle\langle\beta|]$
$\displaystyle=-i\sum_{\alpha^{\prime}}(|\alpha^{\prime}\rangle\langle\alpha^{\prime}|H_{NH}|\alpha\rangle\langle\beta|-|\alpha\rangle\langle\beta|H_{NH}^{\dagger}|\alpha^{\prime}\rangle\langle\alpha^{\prime}|)$
$\displaystyle=-i\sum_{\alpha^{\prime}}((E_{\alpha}\delta_{\alpha\alpha^{\prime}}-i\sum_{i}\frac{\gamma_{i}}{2}C^{(i)}_{\alpha^{\prime}\alpha})|\alpha^{\prime}\rangle\langle\beta|$
(S29)
$\displaystyle\qquad\qquad-(E_{\beta}\delta_{\beta\alpha^{\prime}}+i\sum_{i}\frac{\gamma_{i}}{2}C^{(i)}_{\beta\alpha^{\prime}})|\alpha\rangle\langle\alpha^{\prime}|)$
$\displaystyle=\sum_{\alpha^{\prime}}((E_{\alpha}\delta_{\alpha\alpha^{\prime}}-i\Gamma_{\alpha^{\prime}\alpha})|\alpha^{\prime}\rangle\langle\beta|$
$\displaystyle\qquad\qquad-(E_{\beta}\delta_{\beta\alpha^{\prime}}+i\Gamma_{\beta\alpha^{\prime}})|\alpha\rangle\langle\alpha^{\prime}|),$
(S30)

where we have defined
$\Gamma_{\alpha^{\prime}\alpha}=\sum_{i}\frac{\gamma_{i}}{2}C^{(i)}_{\alpha^{\prime}\alpha}$.
Acting on $\rho^{(a)}_{inc}$, we thus recover

$\displaystyle\mathcal{L}_{NH}[\rho^{(a)}_{inc}]=$
$\displaystyle-\sum_{\begin{subarray}{c}\alpha\beta\\\
\alpha\neq\beta\end{subarray}}p^{(2)}_{\alpha\beta}|\mu_{a\alpha}|^{2}\sum_{\alpha^{\prime}}(\Gamma_{\alpha^{\prime}\beta}|\alpha^{\prime}\rangle\langle\beta|-\Gamma_{\beta\alpha^{\prime}}|\beta\rangle\langle\alpha^{\prime}|)$
$\displaystyle-i\sum_{\begin{subarray}{c}\alpha\beta\\\
\alpha\neq\beta\end{subarray}}p^{(2)}_{\alpha\beta}\mu_{a\alpha}\mu_{a\beta}^{*}\sum_{\alpha^{\prime}}\bigg{[}(E_{\alpha}\delta_{\alpha\alpha^{\prime}}-i\Gamma_{\alpha^{\prime}\alpha})|\alpha^{\prime}\rangle\langle\beta|)$
$\displaystyle\qquad\qquad-(E_{\beta}\delta_{\beta\alpha^{\prime}}+i\Gamma_{\beta\alpha^{\prime}})|\alpha\rangle\langle\alpha^{\prime}|)\bigg{]}.$
(S31)

Similarly, for the jump term, we have

$\displaystyle\mathcal{L}_{J}[|\alpha\rangle\langle\beta|]=$
$\displaystyle\sum_{i}\gamma_{i}\sum_{\alpha^{\prime}\beta^{\prime}}C_{\alpha^{\prime}\alpha}^{(i)}C_{\beta\beta^{\prime}}^{(i)}|\alpha^{\prime}\rangle\langle\beta^{\prime}|,$
(S32)

such that

$\displaystyle\mathcal{L}_{J}[\rho^{(a)}_{inc}]=\sum_{\begin{subarray}{c}\alpha\beta\\\
\alpha\neq\beta\end{subarray}}p^{(2)}_{\alpha\beta}$
$\displaystyle\sum_{i}\gamma_{i}\sum_{\alpha^{\prime}\beta^{\prime}}(|\mu_{a\alpha}|^{2}C_{\alpha^{\prime}\beta}^{(i)}C_{\beta\beta^{\prime}}^{(i)}$
(S33)
$\displaystyle+\mu_{a\alpha}\mu_{a\beta}^{*}C_{\alpha^{\prime}\alpha}^{(i)}C_{\beta\beta^{\prime}}^{(i)})|\alpha^{\prime}\rangle\langle\beta^{\prime}|.$

Now, for step 3), we calculate the second measurement process, given by the
emission to the ground state via $M_{b}^{(1\to 0)}$. This can similarly be
seen to act on matrix elements on the single excited space to yield

$\displaystyle M_{b}^{(1\to
0)}[|\alpha\rangle\langle\beta|]=\mu_{b\alpha}\mu_{b\beta}^{*}|0\rangle\langle
0|.$ (S34)

We thus have, for the non-Hermitian contribution,

$\displaystyle\operatorname{Tr}[M_{b}^{(1\to 0)}[$
$\displaystyle\mathcal{L}_{NH}[\rho^{(a)}_{inc}]]]=$
$\displaystyle-\sum_{\begin{subarray}{c}\alpha\beta\\\
\alpha\neq\beta\end{subarray}}p^{(2)}_{\alpha\beta}|\mu_{a\alpha}|^{2}\sum_{\alpha^{\prime}}(\mu_{b\alpha^{\prime}}\mu_{b\beta}^{*}\Gamma_{\alpha^{\prime}\beta}-\mu_{b\beta}\mu_{b\alpha^{\prime}}^{*}\Gamma_{\beta\alpha^{\prime}})$
$\displaystyle-i\sum_{\begin{subarray}{c}\alpha\beta\\\
\alpha\neq\beta\end{subarray}}p^{(2)}_{\alpha\beta}\mu_{a\alpha}\mu_{a\beta}^{*}\sum_{\alpha^{\prime}}(E_{\alpha}\delta_{\alpha\alpha^{\prime}}-i\Gamma_{\alpha^{\prime}\alpha})\mu_{b\alpha^{\prime}}\mu_{b\beta}^{*})$
$\displaystyle\qquad\qquad-(E_{\beta}\delta_{\beta\alpha^{\prime}}+i\Gamma_{\beta\alpha^{\prime}})\mu_{b\alpha}\mu_{b\alpha^{\prime}}^{*})$
(S35)

and secondly, for the quantum jump term

$\displaystyle\operatorname{Tr}[M_{b}^{(1\to
0)}[\mathcal{L}_{J}[\rho^{(a)}_{inc}]]]=$
$\displaystyle\sum_{\begin{subarray}{c}\alpha\beta\\\
\alpha\neq\beta\end{subarray}}p^{(2)}_{\alpha\beta}\sum_{i}\gamma_{i}\sum_{\alpha^{\prime}\beta^{\prime}}(|\mu_{a\alpha}|^{2}C_{\alpha^{\prime}\beta}^{(i)}C_{\beta\beta^{\prime}}^{(i)}$
$\displaystyle+\mu_{a\alpha}\mu_{a\beta}^{*}C_{\alpha^{\prime}\alpha}^{(i)}C_{\beta\beta^{\prime}}^{(i)})\mu_{b\alpha^{\prime}}\mu_{b\beta^{\prime}}^{*}.$
(S36)

Finally, then, we can observe in each of these contributions which terms lead
to symmetric photon correlations, and which lead to an asymmetry on
interchange of indices $a,b$. We crucially wish to distinguish the role of
terms which couple populations to populations, those which couple coherences
to coherences, and those coupling populations and coherences.

First analysing the non-Hermitian term, Eq. (III), we can note that the first
term manifests from the populations in $\rho^{(a)}_{inc}$, and consists of
terms in the summation coupling these populations to both populations and
coherences. This can be seen via noting that

$\displaystyle
p^{(2)}_{\alpha\beta}|\mu_{a\alpha}|^{2}\mu_{b\alpha^{\prime}}\mu_{b\beta}^{*}\Gamma_{\alpha^{\prime}\beta}$
(S37)

arises from the matrix element
$p^{(2)}_{\alpha\beta}|\alpha\beta\rangle\langle\alpha\beta|$ of the initial
state being acted upon by the initial measurement $M_{a}^{(2\to 1)}$ to obtain
$p^{(2)}_{\alpha\beta}|\mu_{a\alpha}|^{2}|\beta\rangle\langle\beta|$. The non-
Hermitian part of the generator then acts on this state, the first term of
which includes
$\Gamma_{\alpha^{\prime}\beta}=\sum_{i}\frac{\gamma_{i}}{2}C^{(i)}_{\alpha^{\prime}\beta}$
which arises due to the non-Hermitian part of the effective Hamiltonian acting
from the left onto the population $\propto|\beta\rangle\langle\beta|$. This
term thus maps the population $|\beta\rangle\langle\beta|$ to itself for
$\alpha^{\prime}=\beta$, and otherwise maps populations to coherences in the
excitonic basis. We note that the term that maps the population to itself,
$p^{(2)}_{\alpha\beta}|\mu_{a\alpha}|^{2}|\mu_{b\beta}|^{2}\Gamma_{\beta\beta}$,
is symmetric under interchange of indices $a$ and $b$, and thus is symmetric
with respect to reversal of photon measurement order. The terms which map
populations to coherences in the dynamics contribute the asymmetry in Eq.
(S37).

The second term (lines 2 and 3 of the right hand side) of Eq. (III) arises
from the coherences generated by the initial measurement $M_{a}^{(2\to 1)}$,
and contains terms of the form

$\displaystyle
p^{(2)}_{\alpha\beta}\mu_{a\alpha}\mu_{a\beta}^{*}\sum_{\alpha^{\prime}}(E_{\alpha}\delta_{\alpha\alpha^{\prime}}-i\Gamma_{\alpha^{\prime}\alpha})\mu_{b\alpha^{\prime}}\mu_{b\beta}^{*}.$
(S38)

Here we see that the generated coherence $\propto|\alpha\rangle\langle\beta|$
with $\alpha\neq\beta$ under the initial measurement $M_{a}$ interacts (in
this case again from the left) with the non-Hermitian Hamiltonian. The
Hamiltonian part maps eigenstates to themselves, and thus contributes
$E_{\alpha}\delta_{\alpha\alpha^{\prime}}$, and thus this term maps coherences
to coherences. Notably the term
$p^{(2)}_{\alpha\beta}E_{\alpha}\mu_{a\alpha}\mu_{a\beta}^{*}\mu_{b\alpha}\mu_{b\beta}^{*}$
is symmetric on exchange of indices $a$ and $b$.

The second term in Eq. (S38) includes elements which map coherences to
coherences, and coherences to populations. The former include each term in Eq.
(S38) with $\alpha^{\prime}\neq\alpha$, which can be seen from the
corresponding term of Eq. (III). Notably, the term
$p^{(2)}_{\alpha\beta}\Gamma_{\alpha^{\prime}\alpha}\mu_{a\alpha}\mu_{a\beta}^{*}\mu_{b\alpha^{\prime}}\mu_{b\beta}^{*}$
is symmetric under interchange of indices $a$, $b$ only for
$\alpha^{\prime}=\alpha$ \- the term mapping coherences to coherences. Again,
those mapping populations to coherences manifest correlation asymmetry.

We thus see that terms in the non-Hermitian effective Hamiltonian that lead
coherences and populations to evolve separately do not lead to correlation
asymmetry, whereas terms coupling populations and coherences can lead to
asymmetric photon correlations.

We now turn to the quantum jumps term of the master equation. Taking as an
example the first term in (III), which describes the evolution and subsequent
emission event of the populations term in Eq. (S28), the factor
$C_{\alpha^{\prime}\beta}^{(i)}C_{\beta\beta^{\prime}}^{(i)}$ arises from the
jump operator acting on the population $|\beta\rangle\langle\beta|$, which is
mapped to $|\alpha^{\prime}\rangle\langle\beta^{\prime}|$, and thus this term
maps populations to populations when $\alpha^{\prime}=\beta^{\prime}$, and
otherwise maps populations to coherences. We note that when
$\alpha^{\prime}=\beta^{\prime}$ the relevant term is
$p^{(2)}_{\alpha\beta}\gamma_{i}|\mu_{a\alpha}|^{2}C_{\alpha^{\prime}\beta}^{(i)}C_{\beta\alpha^{\prime}}^{(i)}|\mu_{b\alpha^{\prime}}|^{2}$,
which is not in general symmetric under interchange of $a$ and $b$, meaning
that the quantum jumps term mapping populations to populations may indeed
generate asymmetries in photon correlations. Here symmetric correlations are
recovered in the case of no exciton delocalisation, when
$C_{\alpha^{\prime}\beta}^{(i)}=\delta_{\alpha^{\prime}\beta}$. We further
note that the term
$C_{\alpha\beta}^{(i)}C_{\beta\alpha}^{(i)}=\langle\beta|i\rangle\langle
i|\alpha\rangle\langle\alpha|i\rangle\langle i|\beta\rangle$ can be understood
in terms of interferences between excitons $\alpha$ and $\beta$ on site $i$.
For $\alpha=\beta$ the sum over each site corresponds to the inverse
participation ratio of the excitonic state, capturing the delocalisation of
the exciton over all sites.

The second term in Eq. (III) describes the evolution by quantum jumps of the
coherences formed via the initial measurement of $M_{a}^{(2\to 1)}$. This term
maps coherences to populations for $\alpha^{\prime}=\beta^{\prime}$, and
coherences to coherences otherwise. The relevant factor for this term is given
by
$p^{(2)}_{\alpha\beta}\gamma_{i}\mu_{a\alpha}\mu_{a\beta}^{*}C_{\alpha^{\prime}\alpha}^{(i)}C_{\beta\beta^{\prime}}^{(i)}\mu_{b\alpha^{\prime}}\mu_{b\beta^{\prime}}^{*}$,
which is symmetric for real $\mu$ in the $\alpha^{\prime}=\beta^{\prime}$
case, and otherwise only symmetric when this condition is again enforced in
the limit of zero exciton delocalisation. We note that for Hamiltonians
describing time-reversal invariant systems the eigenstates are real valued,
and thus $\mu$ is similarly real in the exciton basis. This occurs for the
Frenkel exciton Hamiltonian in the cases studied where all couplings are real
valued.

We thus see that the terms which may lead to a violation of the symmetry on
interchange of measurement order either manifest as a result of the mapping of
coherences to populations and vice-versa, or arise in the quantum jumps term
as a result of exciton delocalisation. This thus leads to the conclusion that
correlation asymmetry in Frenkel exciton models with environmental couplings
locally acting in the site basis is a direct result of manifestly quantum
behaviours of the system.

##  IV Models

In this section we describe in full the models used in the main text. Each
model consists of three components required to calculate the steady-state and
time evolution of photon cross correlations: i) the Hamiltonian, ii) the
environment as characterised by a spectral density and/or jump operators with
respective rates, iii) the dipole moments, which may be used to obtain their
projection onto a polarization direction as shown in Section V, and obtain the
dipole operators $\hat{\mu}$. All calculations are performed at room
temperature, 300 K.

###  IV.1 Dimer

We choose the dimer model parameters to resemble PC645 Blau  _et al._ (2018).
We thus choose site energies as $\epsilon^{0}_{1}=17317$ cm-1,
$\epsilon^{0}_{1}=15405$ cm-1, and the coupling strength as $J=-46.8$ cm-1
(when not specified), and with a Drude-Lorentz environment with
$\lambda=\Omega=100$ cm-1.

We take the dipole moment of the dimer model to be given by a single angle
$\theta$ dictating the relative orientation of the moments. Facing on in the
molecular $x$ direction, and the other an angle $\theta$ in the $x-y$ plane we
then have $\vec{\mu}_{1}=(\mu,\,0,\,0)$ and
$\vec{\mu}_{2}=(\mu\cos(\theta),\,\mu\sin(\theta),\,0)$. We take $\mu=1$ for
simplicity, and use $\theta=\pi/2$ in the calculations shown. The relative
angle of the molecular and lab frames is further specified by the angles
$\theta_{det},\,\phi_{det}$ (see Sec V below), which we take to be equal to
0.8 $\pi$ and $0.4\pi$, respectively. We note that it is these values as well
as the relative dipole orientations which dictate the overall shape of the
cross-correlation colour plots in e.g Figure 1. Thus we do not observe the
shape change significantly with Hamiltonian parameters, but rather the
magnitude of correlations.

We note that, as a naively applied dimer model we do not expect the direct
results of this model to in practice resemble PC645, which itself has 8 sites
and multiple environmental modes which play an important role in exciton
dynamics, which are not considered here. This model is chosen as a simple
model with similar parameters to those expected in photosynthetic systems.

In Figure S2 we use different parameters chosen to be identical to those in
Ref. Sánchez Munõz and Schlawin (2020), and add an additional pure dephasing
term to the environment via a Gorini, Kossakowski, Lindblad and Sudarshan
(GKSL) master equation with jump operators $L_{pd}=|m\rangle\langle m|$ acting
on each site with an identical rate $\gamma_{pd}$.

###  IV.2 LH2 subunit

The subunit for which we show results for zero-delay photon correlations in
Figures 1c), d) in the main text is chosen as a single $\alpha-\beta$ subunit
of LH2, consisting of a pair of chromophores from the B850 ring, which
consists of 18 total sites, and a single chromophore of the B800 ring, which
is made up of 9 sites. This subunit is thus repeated 9 times in a ring
structure in a single LH2 complex. The parameters we use for this model are
taken from Ref. Cupellini  _et al._ (2020), which we repeat here for clarity.

The Hamiltonian is given by

$\displaystyle H=\begin{pmatrix}12798.4&339.0&-15.6\\\ 339.0&12805.8&-9.7\\\
-15.6&-9.7&13021\end{pmatrix},$ (S39)

where the upper two diagonal elements are the (strongly coupled) $\alpha$ and
$\beta$ sites of the B850 ring, and the third is the B800 site, which is far
more weakly coupled. The Drude-Lorentz environments of the B850 chromophores,
sites 1 and 2, each have reorganization energy $\lambda_{1}=\lambda_{2}=140$
cm-1 and cutoffs $\Omega_{1}=\Omega_{2}=100$ cm-1 , whereas the B800
chromophore has a lower reorganization energy $\lambda_{3}=40$ cm-1, and an
identical cutoff $\Omega_{3}=100$ cm-1. In Figures 1c) we reduce the coupling
$J_{12}$ from 339.0 cm${}^{-1}\gg\lambda_{1}$ to 140 cm${}^{-1}=\lambda_{1}$.
The dipole moments are $\vec{\mu}_{1}=(8.7509,2.9632,1.1907)$,
$\vec{\mu}_{2}=(4.3143,-7.4567,0.842)$, and
$\vec{\mu}_{3}=(-5.6751,6.3565,1.2896)$, respectively. In Figures 1c), d) we
use $\theta_{det}=\pi/2$ and $\eta_{det}=\pi/20$, which are chosen to show the
area in which the zero-delay bound is violated in Figures 1d).

###  IV.3 FMO

The Hamiltonian used is ‘Model C’ of Ref. Kell  _et al._ (2016):

$H_{FMO}=\begin{pmatrix}12405.&-87&4.2&-5.2&5.5&-14.&-6.4&21.\\\
-87&12505.&28.&6.9&1.5&8.7&4.5&4.2\\\ 4.2&28.&12150.&-54.&-0.2&-7.6&1.2&0.6\\\
-5.2&6.9&-54.&12300.&-62.&-16.&-51.&-1.3\\\
5.5&1.5&-0.2&-62.&12470.&60.&1.7&3.3\\\
-14.&8.7&-7.6&-16.&60.&12575.&29.&-7.9\\\
-6.4&4.5&1.2&-51.&1.7&29.&12375.&-9.3\\\
21.&4.2&0.6&-1.3&3.3&-7.9&-9.3&12430\end{pmatrix}.$

Table 1: Components of dipole moments of FMO complex Kramer  _et al._ (2018) can be calculated from the protein data bank (code: 3ENI). | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---
x | -0.74100554 | -0.85714086 | -0.19712114 | -0.79924043 | -0.73692547 | -0.13501747 | -0.49511476 | -0.13838472
y | -0.56060174 | 0.5037757 | 0.95741018 | -0.53357321 | 0.6557619 | -0.879218 | -0.70834118 | 0.82141214
z | -0.36964371 | -0.10732938 | -0.21097155 | -0.27661231 | 0.16406458 | 0.45688729 | -0.50310451 | -0.55329175

As in Ref. Kell  _et al._ (2016) the employed form of the environment is not
amenable to description in terms of the HEOM formalism, we instead use a
Drude-Lorentz form with $\lambda=35$cm-1 and $\Omega=106.1$ cm-1 as in Refs.
Hein  _et al._ (2012a); Kramer  _et al._ (2018) which is also similar to that
used/observed in Sarovar  _et al._ (2010). The dipole moments are assumed to
lie in the direction of a line connecting the NB and ND atoms. The atomic
coordinates were taken from the crystal structure of the FMO trimer. Dipole
moments are shown in Table 1.

In our calculations for the time dependent correlations FMO complex we have
used the molecular $x,\,y,\,$ and $z$ components for simplicity, noting that
these are possible to select from combinations of polarization angle and
molecular orientation, though more importantly that qualitative features with
different selections of these parameters are not altered. For the zero delay
photon statistics, we use the same relative molecule-detector orientations as
the dimer model of $\theta_{det}=0.8\pi$ and $\phi_{det}=0.4\pi$.

##  V Projection of molecular dipole moment onto effective polarization-filter

Figure S1: Examples of relative molecule-detector orientations for obtaining
effective polarization-filtered dipole moment.

The dipole moment of site $i$ is defined as a three component vector
$\vec{\mu}_{m}=(\mu_{m,x},\,\mu_{m,y},\,\mu_{m,z})$. This vector can be
defined in terms of its magnitude
$\mu_{m}=|\vec{\mu}_{m}|=\sqrt{\mu_{m,x}^{2}+\mu_{m,y}^{2}+\mu_{m,z}^{2}}$,
and the azimuthal and polar angles, $\eta=\arccos(\frac{\mu_{z,i}}{\mu})$ and
$\theta=\arctan(\frac{\mu_{y,i}}{\mu_{x,i}})$, respectively. In the simplest
case where the detector is in the molecular $z$ direction, as in Figure S1a),
we have that the projection onto the polarization-filtered component of the
dipole is given by (dropping the site index $m$ for ease of notation)

$\displaystyle\mu_{eff}(\phi)$ $\displaystyle=P(\phi)\vec{\mu}$
$\displaystyle=\begin{pmatrix}\cos(\phi)&\sin(\phi)&0\end{pmatrix}\begin{pmatrix}\mu_{x}\\\
\mu_{y}\\\ \mu_{z}\end{pmatrix}$
$\displaystyle=\mu[\cos(\phi)\sin(\eta)\cos(\theta)+\sin(\phi)\sin(\eta)\sin(\theta)]$
$\displaystyle=\mu\sin(\eta)\cos(\theta-\phi).$ (S40)

We note that for $N$ sites, there is thus a list of $N$ dipole moments each
characterised by the above magnitude and relevant angles in the molecule
frame. The effective polarization-filtered dipole operator is then given by
$\hat{\mu}(\phi)=\sum_{i}^{N}\mu_{i}(\phi)|0\rangle\langle i|+H.c$ in the
single excitation manifold.

We can also define different molecule-detector orientations via the angles
$\theta_{det}$ and $\eta_{det}$ between the detector and molecular frames, as
in Figure S1. For finite $\theta_{det}$ in this case we may simply make the
substitution $\phi\to\phi+\theta_{det}$. In Figure S1b) we illustrate the case
of $\eta_{det}=\frac{\pi}{2},\,\theta_{det}=0$, corresponding to the detector
plane aligning with the molecular y axis. This such a scenario of finite
$\eta_{det}$ we have a rotation of the molecular basis describing the dipole
vector $\vec{\mu}=(\mu_{x},\,\mu_{y},\,\mu_{z})$ with respect to the lab
frame, via

$\displaystyle\vec{\mu}^{\prime}(\eta_{det},0)=P(\phi)R(\eta_{det})\vec{\mu}$
(S41)
$\displaystyle=\begin{pmatrix}\cos(\phi)&\sin(\phi)&0\end{pmatrix}\begin{pmatrix}1&0&0\\\
0&\cos(\eta_{det})&\sin(\eta_{det})\\\
0&-\sin(\eta_{det})&\cos(\eta_{det})\end{pmatrix}\begin{pmatrix}\mu_{x}\\\
\mu_{y}\\\ \mu_{z}\end{pmatrix}.$

Similarly, in S1c) we have both finite $\eta_{det}$ and $\theta_{det}$, such
that the relevant rotation becomes

$\displaystyle\vec{\mu}^{\prime}(\eta_{det},\theta_{det})$
$\displaystyle=P(\phi)R(\eta_{det})S(\theta_{det})\vec{\mu}$
$\displaystyle=\begin{pmatrix}\cos(\phi)&\sin(\phi)&0\end{pmatrix}\begin{pmatrix}1&0&0\\\
0&\cos(\eta_{det})&\sin(\eta_{det})\\\
0&-\sin(\eta_{det})&\cos(\eta_{det})\end{pmatrix}\begin{pmatrix}\cos(\theta_{det})&\sin(\theta_{det})&0\\\
-\sin(\theta_{det})&\cos(\theta_{det})&0\\\
0&0&1\end{pmatrix}\begin{pmatrix}\mu_{x}\\\ \mu_{y}\\\ \mu_{z}\end{pmatrix}.$
(S42)

In general, then, the effective dipole moment can be obtained by a two step
process of the required rotations $\vec{\mu}\to\vec{\mu}^{\prime}$ such that
the dipole aligns with the $z-$axis, followed by taking the projection in that
basis
$\vec{\mu}^{\prime}\to\mu^{\prime}_{x}\cos(\phi)+\mu^{\prime}_{y}\sin(\phi)$.
We thus define the dipole operator for site $m$ as
$\hat{\mu}_{\phi,m}=\vec{\mu}_{\phi}^{\prime}(\eta_{det},\theta_{det})|m\rangle\langle
0|+h.c$, and the total effective dipole operators as
$\hat{\mu}_{\phi}=\sum_{m}\hat{\mu}_{\phi,m}$.

Figure S2: zero delay photon correlations for dimer model of Ref. Sánchez
Munõz and Schlawin (2020) and the dependence on pure dephasing. a) Shows
$\gamma_{pd}=0$, b) shows $\gamma_{pd}=10\gamma$, c) shows
$\gamma_{pd}=1000\gamma$. We note that the latter is closest to the case of
photosynthetic systems, where the timescales of dephasing and decoherence due
to vibrational background are $\sim 1000\times$ faster (ps) than those of
excited state decay to the EM environment (ns).

Figure S3: a) Long time dependence of photon cross-correlations for LH2
subunit, and b) the corresponding asymmetry. c) Shows long time dependence of
photon cross-correlations of the FMO complex. Figure S4: Correlation asymmetry
dependence on coupling rate $J$ and pure dephasing rate $\gamma_{pd}$ for
Dimer model. a), d) show $\gamma_{pd}=0$ as in main text, b), e) show
$\gamma_{pd}=10\gamma$, and c), f) show $\gamma_{pd}=1000\gamma$. For each
plot $\phi_{1}=0,\,\phi_{2}=\frac{\pi}{2}$.

##  VI Additional numerical results

In this section we provide some additional numerical analysis to supplement
the conclusions of the main text. Firstly, for the Dimer model with identical
parameters to that studied in Ref. Sánchez Munõz and Schlawin (2020) with an
included Markovian environmental pure-dephasing, we show in Figure S2 that a
Markovian environment can remove the ability of zero delay correlations to
witness coherence in the steady-state for large values of pure dephasing
$\gamma_{pd}$.

In Figure S3 we show the time dependence of photon correlations for the LH2
subunit model and the FMO complex, showing indeed that at timescales longer
than the $\sim 1ns$ decay timescale of the electronic excited states photons
correlations tend to the uncorrelated limit of $g^{(2)}_{\phi_{1}\phi_{2}}=1$.
It is further notable that the correlation asymmetry similarly lasts for this
long time, which is not observed in the FMO complex (see Figure S6). This is
due to the larger system size $M$ of FMO, combined with the essentially random
projection of the polarization filter onto the exciton basis, leading to
decaying asymmetry over the picosecond timescales of the vibrational
environment. In the LH2 subunit, after this picosecond timescale some
asymmetry remains as the projection of the polarization filter onto the
excitonic basis biases strongly particular states. Future work will analyze in
more detail the role of the projection of the polarization filter onto
particular states onto the correlation asymmetry, which dictates the role of
the vibrational environment in the decay of asymmetries over short times.

In Figure S4 we similarly show how photon correlation asymmetry depends on the
pure dephasing rate $\gamma_{pd}$. We observe that coherent oscillations are
rapidly destroyed by the pure dephasing, and the observed asymmetries are
significantly reduced. Nonetheless, for each case in which coherences in the
steady-state are present, that is, for finite $V$, asymmetries in photon
correlations are similarly present. In the case of $V=0$ the dynamics is
effectively completely classical, with no coherences in the steady-state,
asymmetries are not observed in this case.

In Figure S5 we show the PC645 inspired dimer model cross-correlation time-
dependence, as well as the total asymmetries, excitonic coherences, and basis
independent coherence, and how each change with both electronic coupling $J$,
and reorganization energy $\lambda$. Here we see similar results to the LH2
subunit shown in Figure 2 of the main text, where the correlation asymmetry
closely resembles the coherence measures dependence on electronic coupling,
but deviates in its behaviour for reorganization energy, as whilst the steady-
state indeed has higher excitonic and basis independent coherences as a result
of the environmental coupling, the dynamics is more dominated by environmental
degrees of freedom and is thus more incoherent.

In Figure S5a), b) we show the photon correlations and their asymmetry for the
dimer model with $\lambda=\Omega=100$ cm-1, for various values of the inter-
site electronic coupling $J$. We see that electronic coupling induces an
asymmetry in photon correlations. In the $J=0$ case, the excitonic and site
bases are identical, and the environment thus does not enable exciton
transport or sustain coherences in the steady-state, leading to symmetric
correlations upon interchange of measurement order. We see that asymmetries
are introduced even for very weak inter-site couplings, and hence asymmetry
acts as a very sensitive witness of coherent excitation transfer.

In the dimer model we see a more complex dependence of the basis independent
coherence on reorganization energy, which can be explained as the singly
excited state, which dominates the change in the basis independent coherence,
only has two possible excitonic states. At weak pumping, there is no transfer,
and thus pumping the highest energy exciton leads to a steady state with only
population in this excitonic state, and thus a high $C_{1}(\rho)$ (low
entropy), as $\lambda$ is increased, exciton transfer is enabled, reducing the
excitonic coherence. For higher $\lambda$ values, however, transport to the
lowest energy exciton is stronger than pumping of the high energy, and thus
steady-state population is dominant in the low energy exciton, again
increasing the basis independent coherence due to the low entropy of the
state. Thus, we see that for the dimer model, $C_{1}(\rho)$ does not fully
capture the relevant coherent properties of excitation transport, which
decrease as the environment dominates this process. This is observed for the
correlation asymmetry, as the transport is more incoherent as $\lambda$ is
increased.

Figure S5: Photon correlations of the dimer model. Time dependence a), time-
asymmetry b), total asymmetry c), d) over $\tau\approx 13.25$ ps
$>>\tau_{S-E}$, excitonic coherence e), f) and basis independent coherence g),
h) and their dependence on reorganization energy $\lambda$ and coupling
strength $J$. $N_{trunc}=58$ (note that this is a very large value as the
basis independent coherence for the dimer model converges very slowly due to
its high sensitivity to small changes in populations and coherences in the
steady-state.)

In Figure S6 we show how photon correlations and their asymmetries depend on
the pump power in the FMO complex. We note that a change is observed in the
values of $g^{(2)}$ as the pump power is altered for all delay times. This is
due to the fact that the dynamical generator driving dynamics at early times
is dominated by the inter-exciton dynamics, and independent of pump power,
however the steady-state populations depend strongly on the pump power. Thus
the time evolution follows near identical trends, with differing absolute
values of the photon correlations. The data in Figure S6 are examples of those
used in the main text to obtain Figure 3g), in which we observe a steady
decrease in the asymmetry.

Figure S6: Pump power dependence of polarization filtered photon correlations
a)-c) and their asymmetry d)-f) for the FMO complex. We use polarization
filters corresponding to the $x,y,$ and $z$ components for simplicity,
observing asymmetries in each case. We see that asymmetries can be observed
over a few picoseconds, which is the typical timescale for exciton transport.

We additionally in Figure S7 show the Fourier transform (exploiting the
standard Blackman window filtering method) of the photon correlation asymmetry
results in Figure 3f) of the main text. Here we observe that, indeed, the
prominent early time oscillatory behaviours are of comparable frequency to the
exciton energy gaps.

Figure S7: Fourier transform of photon correlation asymmetry of FMO. Solid
vertical lines are all exciton energy gaps. We observe, indeed, that the
prominent oscillation frequencies occur in the same frequency range as that
expected in coherent inter-exciton transport.

In Figure S8 we show the asymmetry of polarization filtered correlation
functions for the FMO complex at early times, with the time integration of
asymmetry taken as 30 fs. We see that as the reorganization energy $\lambda$
is increased the asymmetry for some values indeed increases, which manifests
as a result of early time environmentally induced coherent transfer. This
effect is much smaller in FMO compared to the LH2 subunit shown in the main
text, which we associate to the multichromophoric system having many
contributions to the transport in different directions, leading to a washing
out of the induced asymmetries.

Finally, in Figure S9 we show that zero time photon correlations are robust in
their form to environmental changes in FMO by altering the Drude-Lorentz
reorganisation energy, observing only very slight changes in the zero-delay
correlations over a wide range of environmental parameters. We note that this
is consistent with the results observed for the dimer model, where whilst we
observe that the region of violation of the zero-delay bound indeed changes,
this is around regions where a very small change in the zero-delay values has
a large effect. We further show that a tightened form of the bound in Ref.
Sánchez Munõz and Schlawin (2020), written as

$\displaystyle g^{(2)}_{\phi_{1},\phi_{1}}(0)\leq
2-\frac{2\sum_{\alpha}|\mu_{\phi_{1},\alpha}|^{2}\langle\sigma_{\alpha}^{\dagger}\sigma_{\alpha}\rangle}{\sum_{\alpha,\beta}|\mu_{\phi_{1},\alpha}|^{2}|\mu_{\phi_{1},\beta}|^{2}\langle\sigma_{\alpha}^{\dagger}\sigma_{\alpha}\rangle\langle\sigma_{\beta}^{\dagger}\sigma_{\beta}\rangle},$
(43)

where $\mu_{\phi_{1},\alpha}$ is the component of the dipole moment on exciton
$\alpha$, and $\sigma_{\alpha}=|0\rangle\langle\alpha|$, is also not violated
for any of the shown angles in FMO.

Figure S8: Asymmetry of FMO at early times, with integration time of asymmetry
= 30 fs. Figure S9: Zero time polarization filtered correlations for varying
reorganization energy (solid line), and their comparison to the tightened
classical bound (dashed line). We see here that the zero-delay correlation
functions change very slowly with environmental parameters of the background
Drude-Lorentz mode. There is no violation of the zero-delay bound for any
value of the reorganization energy.

##  VII Combined Nakajima-Zwanzig Hierarchical equations of motion approach

###  VII.1 Set up

Here we apply the Nakajima-Zwanzig (NZ) formalism to the hierarchical
equations of motion (HEOM) in order to obtain a non-additive description of
both the interaction with the non-Markovian environment of the pigment protein
complex via HEOM, as well as the incoherent Markovian interaction with the
electromagnetic environment and pump. We begin by writing our Hamiltonian in
the form

$H=H_{S}+H_{E_{0}}+H_{SE_{0}}+\sum_{n=1}^{N_{heom}}(H_{E_{n}}+H_{SE_{n}}),$
(44)

where our goal is to treat the environmental contribution of $H_{E_{0}}$
perturbatively, and the other $N_{heom}$ environmental contributions with
HEOM. We thus redefine $H_{S^{\prime}}=H_{S}+H_{E_{0}}+H_{SE_{0}}$, which in
our case is the molecular electronic degrees of freedom $H_{S}$, and the
perturbative environment $E_{0}$ of the pump and decay into the
electromagnetic field. We indicate a density operator of the combined space of
$S$ and $E_{0}$ as $\sigma$. We write the interaction Hamiltonian between the
molecular electronic degrees of freedom and its vibrational environment, which
we wish to treat non-perturbatively via HEOM, as

$\displaystyle H_{E_{n}}=V_{n}\otimes B_{n},$ (45)

where $V_{n}$ is a system operator on the space of $H_{S}$, and $B_{n}$ a bath
operator on that of $H_{B,n}$.

HEOM is a formalism for non-perturbative simulations of open systems with
exponentially decaying bath correlation functions $C_{m}(t)$ on each site $m$
Ishizaki and Tanimura (2006); Ishizaki and Fleming (2009),

$C_{m}(t)=\sum_{k}^{K_{m}}c_{m,k}e^{-\gamma_{m,k}t}.$ (46)

We then have a set of coefficients $c_{m,k}$, as well as a corresponding set
of rates $\gamma_{m,k}$, and operators $B_{m,k}=B_{m}$. We relabel $k$ by
merging these indices into a single index $k\in[1,K]$, with $K=\sum_{m}K_{m}$.
The formalism then becomes identical to the case of a single bath, except that
each of the coefficients may now have different coupling operators.

Then the HEOM describe a set of coupled equations of motion for a set of
auxiliary density operators (ADOs) $\sigma_{\mathbf{n}}$, where the multi-
index $\mathbf{n}=(n_{1},\cdots,n_{K})$ of the ADO is a set of $K$ positive
integers, where $\sum_{k}^{K}\mathbf{n}_{k}=T_{\mathbf{n}}$ is the tier of the
ADO. The zero tier consists of one ADO only, which is the reduced density
operator of the system, and the coupling additional ADOs captures non-
perturbative dynamics of the environment. The ADOs are coupled between
adjacent tiers, and not within the same tier.

For the Drude-Lorentz environments treated here we have

$\displaystyle c_{0}$
$\displaystyle=\lambda\Omega\left(\cot\left(\frac{\beta\Omega}{2}\right)-i\right),$
(47) $\displaystyle c_{k}$
$\displaystyle=\frac{4\lambda\Omega}{\beta}\frac{\nu_{k}}{\nu_{k}^{2}-\Omega^{2}},$
(48) $\displaystyle\gamma_{0}$ $\displaystyle=\Omega,$ (49)
$\displaystyle\gamma_{k}$ $\displaystyle=\nu_{k},$ (50)

with $k=1,2,3\ldots$ and $\nu_{k}=\frac{2\pi k}{\beta}$, with
$\beta=(k_{B}T)^{-1}$ is the inverse temperature, with $k_{B}=1$ the Boltzmann
constant. $\nu_{k}$ are the Matsubara frequencies. We note that for high
temperatures treated here the expansion can be cut-off at simply the zero
order term $c_{0},\gamma_{0}$.

The HEOM for $H_{S^{\prime}}$ can be written as Ishizaki and Fleming (2009)

$\displaystyle\partial_{t}\sigma_{\mathbf{n}}(t)$
$\displaystyle=\mathfrak{L}^{D}_{\mathbf{n}}\sigma_{\mathbf{n}}(t)+\sum_{k}^{K}(\mathfrak{L}^{+}_{\mathbf{n}}\sigma_{\mathbf{n}_{k+}}(t)+\mathfrak{L}^{-}_{\mathbf{n}}\sigma_{\mathbf{n}_{k-}}(t))$
$\displaystyle=\sum_{\mathbf{n}}\mathfrak{L}_{\mathbf{n}}\sigma_{\mathbf{n}^{\prime}}(t)$
(51)
$\displaystyle=\sum_{\mathbf{n}^{\prime}\in\\{\mathbf{n}+{k-},\mathbf{n},\mathbf{n}_{k+}\\}}\mathfrak{L}_{\mathbf{n}}\sigma_{\mathbf{n}^{\prime}}(t),$

with
$\mathfrak{L}^{D}_{\mathbf{n}}[\sigma_{\mathbf{n}}(t)]=\mathfrak{L}_{0}[\sigma_{\mathbf{n}}(t)]-\sum_{n_{k}\in\mathbf{n}}n_{k}\nu_{k}\mathbb{1}_{E_{0}}$
describing the ‘diagonal’ contribution of the ADO, where
$\mathfrak{L}_{0}=-iH_{S^{\prime}}^{\times}$, and the inter-tier couplings are
$\mathfrak{L}_{m,k}^{-}=-in_{m,k}(\text{Re}(c_{m,k})V_{m}^{\times}+i\text{Im}(c_{m,k})V_{m}^{\circ})$,
$\mathfrak{L}_{m,k}^{+}=-iV_{m}^{\times}$, and $\mathbf{n}_{k\pm}$ differs
from $\mathbf{n}$ by $\pm 1$ in the $(k)$-th entry of the tuple $\mathbf{n}$,
$n_{k}$. Here we have used the notation $O^{\times}\cdot:=[O,\cdot]$ and
$O^{\circ}\cdot:=[O,\cdot]_{+}$.

Note that the hierarchical expansion that describes the environments
$E_{i}\,|\,i\geq 1$ has not altered the form of the environment $E_{0}$ or
it’s coupling to the system, however, the effect of this latter environment
alters the self interaction term $\mathfrak{L}_{\mathbf{n}}^{D}(t)$ at every
level of the hierarchy. We note that in our calculations we exploit the
Ishizaki-Tanimura truncation scheme Ishizaki and Tanimura (2005), which adds a
term to the diagonal $\mathfrak{L}^{D}_{\mathbf{n}}$ aiding in convergence of
the Hierarchy, which however does not alter otherwise the form of the HEOM.

Crucially, $E_{0}$ also effects the form of each ADO, as
$\rho^{\prime}_{\mathbf{n}}(0)=\rho_{\mathbf{n},S}(0)\otimes\rho_{\mathbf{n},E_{0}}(0)$.
If we wish to trace out the environment $E_{0}$, we must thus do so at every
tier in the hierarchy.

We note that in the absence of the environment $E_{0}$, we have the regular
HEOM acting on density operators $\rho\in\mathcal{H}$, which is similarly
written as

$\displaystyle\partial_{t}\rho_{\mathbf{n}}(t)$
$\displaystyle=\mathcal{L}^{D}_{\mathbf{n}}\rho_{\mathbf{n}}(t)+\sum_{k}^{K}(\mathcal{L}^{+}_{\mathbf{n}}\rho_{\mathbf{n}_{k+}}(t)+\mathcal{L}^{-}_{\mathbf{n}}\rho_{\mathbf{n}_{k-}}(t)).$
(52)

with
$\mathcal{L}^{D}_{\mathbf{n}}[\rho_{\mathbf{n}}(t)]=-i[H_{S},\rho_{\mathbf{n}}(t)]-\sum_{n_{k}\in\mathbf{n}}n_{k}\nu_{k}\mathbb{1}_{E_{0}}$,
and similarly to above, $\mathcal{L}_{0}=-iH_{S}^{\times}$,
$\mathcal{L}_{m,k}^{-}=-in_{m,k}(\text{Re}(c_{m,k})V_{m}^{\times}+i\text{Im}(c_{m,k})V_{m}^{\circ})$,
$\mathcal{L}_{m,k}^{+}=-iV_{m}^{\times}$. Similarly, if only the environment
$E_{0}$ was present, and we treated this as a Markovian environment of GKSL
form, we would write

$\partial_{t}\rho_{\mathbf{0}}(t)=\mathcal{L_{S}}\rho_{\mathbf{0}}(t)+\mathcal{D}_{E_{0}}\rho_{\mathbf{0}}(t)$
(53)

where $\mathcal{D}_{E_{0}}$ is the GKSL dissipator. Were we to treat these two
environments under an additive approximation, we could simply replace
$\mathcal{L}^{D}_{\mathbf{0}}\to\mathcal{L}^{D}_{\mathbf{0}}+\mathcal{D}_{E_{0}}$,
where the additional term acts only on the system density operator
$\rho_{\mathbf{0}}$. We will see, however, that the dissipation must act on
every tier of the hierarchy.

###  VII.2 Derivation of NZ-HEOM

We begin by writing total density operator of the entire HEOM space as the
total density operator as
$\sigma_{\Delta}(t)=\sum_{\mathbf{n}}\sigma_{\mathbf{n}}$, with

$\displaystyle\partial_{t}\sigma_{\Delta}(t)=\mathfrak{L}_{\Delta}[\sigma_{\Delta}(t)],$
(54)

in the following we use a subscripted ${\Delta}$ to refer to the entire
hierarchy of coupled ADOs - the space of these ADOs, either $\mathcal{H}_{S}$
or $\mathcal{H}_{S}\otimes\mathcal{H}_{E_{0}}$, is denoted by the notation
$\mathcal{L},\rho_{\mathbf{n}}$ or $\mathfrak{L},\sigma_{\mathbf{n}}$,
respectively. Our approach will be similar to that in Ref. Fay and Limmer
(2022), involving tracing out the $E_{0}$ environment with a projection
operator approach to obtain a hierarchical set of equations for
$\rho_{\mathbf{n}}$ from the exact HEOM of Eq. (54).

There are multiple formulations of projection operator based approaches to the
description of open quantum systems, the most well known of which are
Nakajima-Zwanzig Nakajima (1958); Zwanzig (1960) (NZ) and time-convolutionless
(TCL) Tokuyama and Mori (1976) master equations. Projection operator
techniques rely on the definition of projectors $\mathcal{P}$ onto the
‘relevant’ Hilbert space that defines the open system, and it’s complement
$\mathcal{Q}$, the ‘irrelevant’ part. Here we follow the NZ approach, which is
summarized in its more standard form in Section VIII below.

The total system Liouvillian can be written as
$\mathfrak{L}_{0}=\mathfrak{L}_{S}+\mathfrak{L}_{E_{0}}+\alpha\mathfrak{L}_{SE_{0}}$,
where we introduce the unitless coupling strength $\alpha$ for later
convenience in keeping track of the order in the interaction between the
system part $S$ and perturbative environment $E_{0}$. We then define the
projection operators $\mathcal{P}_{\mathbf{n}}$ by
$\mathcal{P}_{\mathbf{n}}\sigma_{\Delta}=\operatorname{Tr}_{E_{0}}[\sigma_{\mathbf{n}}]\otimes\rho_{E_{0}}^{\beta}=\rho_{\mathbf{n}}\otimes\rho_{E_{0}}^{\beta}$,
where $\rho_{E_{0}}^{\beta}$ is the thermal state of $E_{0}$ at inverse
temperature $\beta$ (note that in general this can be any reference state,
typically chosen to be a Gaussian state). The irrelevant part projector is
then
$\mathcal{Q}=\mathbb{1}-\mathcal{P}=\mathbb{1}-\sum_{\mathbf{n}}\mathcal{P}_{\mathbf{n}}=\sum_{\mathbf{n}}(\mathbb{1}_{\mathbf{n}}-\mathcal{P}_{\mathbf{n}})=\sum_{\mathbf{n}}\mathcal{Q}_{\mathbf{n}}$.

We can now follow the standard NZ projection operator technique, however in
this case, there are many ‘relevant’ projection operators that together form a
coupled set of equations

$\displaystyle\partial_{t}\mathcal{P}_{\mathbf{n}}\sigma_{\Delta}(t)$
$\displaystyle=\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{\Delta}\sigma_{\Delta}(t)$
(55)
$\displaystyle=\sum_{\mathbf{n}^{\prime}}\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{\Delta}\mathcal{P}_{\mathbf{n}^{\prime}}\sigma_{\Delta}(t)+\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{\Delta}\mathcal{Q}\sigma_{\Delta}(t)$

and for the irrelevant part

$\displaystyle\partial_{t}\mathcal{Q}\sigma_{\Delta}(t)$
$\displaystyle=\mathcal{Q}\mathfrak{L}_{\Delta}\sigma_{\Delta}(t)$ (56)
$\displaystyle=\sum_{\mathbf{n}^{\prime}}\mathcal{Q}\mathfrak{L}_{\Delta}\mathcal{P}_{\mathbf{n}^{\prime}}\sigma_{\Delta}(t)+\mathcal{Q}\mathfrak{L}_{\Delta}\mathcal{Q}\sigma_{\Delta}(t).$

We can then take the Laplace transform of each of the above two expressions to
find

$\displaystyle s\mathcal{P}_{\mathbf{n}}\tilde{\sigma}_{\Delta}(s)$
$\displaystyle-\mathcal{P}_{\mathbf{n}}\sigma_{\Delta}(0)$ (57)
$\displaystyle=\sum_{\mathbf{n}^{\prime}}\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{\Delta}\mathcal{P}_{\mathbf{n}^{\prime}}\tilde{\sigma}_{\Delta}(s)+\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{\Delta}\mathcal{Q}\tilde{\sigma}_{\Delta}(s)$

and

$\displaystyle s\mathcal{Q}\tilde{\sigma}_{\Delta}(s)$
$\displaystyle-\mathcal{Q}\sigma_{\Delta}(0)$ (58)
$\displaystyle=\sum_{\mathbf{n}^{\prime}}\mathcal{Q}\mathfrak{L}_{\Delta}\mathcal{P}_{\mathbf{n}^{\prime}}\tilde{\sigma}_{\Delta}(s)+\mathcal{Q}\mathfrak{L}_{\Delta}\mathcal{Q}\tilde{\sigma}_{\Delta}(s)$

where $\tilde{\sigma}(s)_{\Delta}=(s-\mathfrak{L}_{\Delta})\sigma(0)$. We can
rearrange Eq. (58) to obtain

$\displaystyle\mathcal{Q}\tilde{\sigma}_{\Delta}(s)=\frac{1}{s-\mathcal{Q}\mathfrak{L}_{\Delta}}(\mathcal{Q}\sigma_{\Delta}(0)+\sum_{\mathbf{n}^{\prime}}\mathcal{Q}\mathfrak{L}_{\Delta}\mathcal{P}_{\mathbf{n}^{\prime}}\tilde{\sigma}_{\Delta}(s))$
(59)

and then substitute this expression for the irrelevant part into that for the
relevant part, obtaining

$\displaystyle s\mathcal{P}_{\mathbf{n}}$
$\displaystyle\tilde{\sigma}_{\Delta}(s)-\mathcal{P}_{\mathbf{n}}\sigma_{\Delta}(0)$
$\displaystyle=\sum_{\mathbf{n}^{\prime}}\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{\Delta}\mathcal{P}_{\mathbf{n}^{\prime}}\tilde{\sigma}_{\Delta}(s)$
(60)
$\displaystyle+\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{\Delta}\frac{1}{s-\mathcal{Q}\mathfrak{L}_{\Delta}}(\mathcal{Q}\sigma_{\Delta}(0)+\sum_{\mathbf{n}^{\prime}}\mathcal{Q}\mathfrak{L}_{\Delta}\mathcal{P}_{\mathbf{n}^{\prime}}\tilde{\sigma}_{\Delta}(s)).$

The inverse Laplace transform then gives us

$\displaystyle\partial_{t}\mathcal{P}_{\mathbf{n}}\sigma_{\Delta}(t)=\sum_{\mathbf{n}^{\prime}}$
$\displaystyle\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{\Delta}\mathcal{P}_{\mathbf{n}^{\prime}}\sigma_{\Delta}(t)+\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{\Delta}e^{\mathcal{Q}\mathfrak{L}_{\Delta}t}\mathcal{Q}\sigma_{\Delta}(0)$
$\displaystyle+\sum_{\mathbf{n}^{\prime}}\int_{0}^{t}d\tau\mathcal{K}_{\mathbf{n},\mathbf{n}^{\prime}}(t-\tau)\sigma_{\Delta}(\tau),$
(61)

where we have defined the superoperator valued NZ kernel
$\mathcal{K}_{\mathbf{n},\mathbf{n}^{\prime}}(t)=\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{\Delta}e^{\mathcal{Q}\mathfrak{L}_{\Delta}t}\mathcal{Q}\mathfrak{L}_{\Delta}\mathcal{P}_{\mathbf{n}^{\prime}}$.
We thus see that the NZ kernel may, in general, mix tiers of the ADOs. In the
following we ignore the second term on the RHS of the above expression, the
so-called inhomogeneous term, which identically vanishes for product state
initial states.

###  VII.3 Weak coupling limit: Markov NZ-HEOM

The set of projections onto $\mathfrak{L}_{\Delta}$ define the elements of the
dynamical generator in the relevant subspace:

$\displaystyle\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{\Delta}$
$\displaystyle\mathcal{P}_{\mathbf{n}^{\prime}}=$ (62)
$\displaystyle\mathcal{P}_{\mathbf{n}}(\mathfrak{L}_{\mathbf{n}}^{D}{\delta}_{\mathbf{n}\mathbf{n}^{\prime}}+\sum_{k}^{K}(\mathfrak{L}^{+}_{\mathbf{n}}{\delta}_{\mathbf{n}\mathbf{n}^{\prime}_{k+}}+\mathfrak{L}^{-}_{\mathbf{n}}{\delta}_{\mathbf{n}\mathbf{n}^{\prime}_{k-}})\mathcal{P}_{\mathbf{n}^{\prime}}$

we then have projectors mixing relevant and irrelevant subspaces,

$\displaystyle\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{\Delta}\mathcal{Q}=\alpha\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{SE_{0}}\mathcal{Q}$
(63a)
$\displaystyle\mathcal{Q}\mathfrak{L}_{\Delta}\mathcal{P}_{\mathbf{n}}=\alpha\mathcal{Q}\mathfrak{L}_{SE_{0}}\mathcal{P}_{\mathbf{n}}$
(63b)

and finally the irrelevant subspace contributions, which we can see differ for
each ADO. Using that
$\mathcal{Q}=\mathbb{1}-\sum_{\mathbf{n}}\mathcal{P}_{\mathbf{n}}=\sum_{\mathbf{n}}(\mathbb{1}_{\mathbf{n}}-\mathcal{P}_{\mathbf{n}})$,
we can define
$\mathcal{Q}_{\mathbf{n}}=\mathbb{1}_{\mathbf{n}}-\mathcal{P}_{\mathbf{n}}$,
and write

$\displaystyle\mathcal{Q}\mathfrak{L}_{\Delta}\mathcal{Q}$
$\displaystyle=\sum_{\mathbf{n},\mathbf{n}^{\prime}}\mathcal{Q}_{\mathbf{n}}\mathfrak{L}_{\Delta}\mathcal{Q}_{\mathbf{n}}$
$\displaystyle=\sum_{\mathbf{n}}\mathcal{Q}_{\mathbf{n}}\mathfrak{L}_{\mathbf{n}}\mathcal{Q}_{\mathbf{n}}$
$\displaystyle=\sum_{\mathbf{n}}(\mathbb{1}_{\mathbf{n}}-\mathcal{P}_{\mathbf{n}})\mathfrak{L}_{\mathbf{n}}\mathcal{Q}_{\mathbf{n}}$
(64)
$\displaystyle=\sum_{\mathbf{n}}\mathfrak{L}_{\mathbf{n}}\mathcal{Q}_{\mathbf{n}}-\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{SE_{0},\mathbf{n}}\mathcal{Q}_{\mathbf{n}}$
$\displaystyle=\sum_{\mathbf{n}}\mathfrak{L}_{\mathbf{n}}\mathcal{Q}_{\mathbf{n}}-(\mathbb{1}_{\mathbf{n}}-\mathcal{Q}_{\mathbf{n}})\mathfrak{L}_{SE_{0},\mathbf{n}}\mathcal{Q}_{\mathbf{n}}$
$\displaystyle=\sum_{\mathbf{n}}(\mathfrak{L}_{S,\mathbf{n}}\mathcal{Q}+\mathfrak{L}_{E_{0},\mathbf{n}}\mathcal{Q}+\mathcal{Q}\mathfrak{L}_{SE_{0},\mathbf{n}}\mathcal{Q}).$

where in the second line we have used that the irrelevant subspace, that
including the interaction term with $E_{0}$, has no terms that mix tiers of
the hierarchy, and in the final line defined

$\mathfrak{L}_{S,\mathbf{n}}=\mathfrak{L}_{S}-\sum_{n_{k}\in\mathbf{n}}n_{k}\nu_{k}\mathbb{1}_{S}+\sum_{k}(\mathfrak{L}^{+}_{\mathbf{n}_{k+}}+\mathfrak{L}^{-}_{\mathbf{n}_{k-}}),$
(65)

and
$\mathfrak{L}_{E_{0},\mathbf{n}}=\mathfrak{L}_{E_{0}}-\sum_{n_{k}\in\mathbf{n}}n_{k}\nu_{k}\mathbb{1}_{E_{0}}$.
We have included the subscript $\mathbf{n}$ in the $SE_{0}$ coupling to
specify on which ADO they act, though note these terms are identical for each
ADO. We then wish to obtain the HEOM including the influence of $H_{E_{0}}$ up
to second-order in $\alpha$, and to do so, as above, take the Laplace
transform of the relevant part(s) of the full ADOs given in Eq. (55)

$\displaystyle s\mathcal{P}_{\mathbf{n}}$
$\displaystyle\tilde{\sigma}_{\Delta}(s)=\mathcal{P}_{\mathbf{n}}\tilde{\sigma}(0)+\sum_{\mathbf{n}^{\prime}}\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{\Delta}\mathcal{P}_{\mathbf{n}^{\prime}}\tilde{\sigma}_{\Delta}(s)+\alpha^{2}\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{\Delta}\frac{1}{s-\mathcal{Q}\mathfrak{L}_{\Delta}}\mathfrak{L}_{\Delta}\sum_{\mathbf{n}^{\prime}}\mathcal{P}_{\mathbf{n}^{\prime}}\tilde{\sigma}_{\Delta}(s).$
(66)

We note that this equation is exact up to the inhomogeneous term, which is
identically zero for an initial product state
$\sigma_{\mathbf{n}}(0)=\rho_{\mathbf{n}}(0)\otimes\rho_{E_{0}}^{\beta}$ for
each ADO. Using Eqs (63) and (VII.3) we obtain

$\displaystyle
s\mathcal{P}_{\mathbf{n}}\tilde{\sigma}_{\Delta}(s)=\mathcal{P}_{\mathbf{n}}$
$\displaystyle\tilde{\sigma}(0)+\sum_{\mathbf{n}^{\prime}}\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{\Delta}\mathcal{P}_{\mathbf{n}^{\prime}}\tilde{\sigma}_{\Delta}(s)$
(67)
$\displaystyle+\alpha^{2}\sum_{\mathbf{n}^{\prime}}\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{SE_{0},\mathbf{n}}\frac{1}{s-(\mathfrak{L}_{S,\mathbf{n}^{\prime}}\mathcal{Q}_{\mathbf{n}^{\prime}}+\mathfrak{L}_{E_{0},\mathbf{n}^{\prime}}\mathcal{Q}_{\mathbf{n}^{\prime}}+\alpha\mathcal{Q}_{\mathbf{n}^{\prime}}\mathfrak{L}_{SE_{0},\mathbf{n}^{\prime}}\mathcal{Q}_{\mathbf{n}^{\prime}})}\mathfrak{L}_{SE_{0},\mathbf{n}}\mathcal{P}_{\mathbf{n}^{\prime}}\tilde{\sigma}_{\Delta}(s).$

We now treat the interaction with $E_{0}$ to lowest order in $\alpha$,
obtaining

$\displaystyle
s\mathcal{P}_{\mathbf{n}}\tilde{\sigma}_{\Delta}(s)=\mathcal{P}_{\mathbf{n}}\tilde{\sigma}(0)+$
$\displaystyle\sum_{\mathbf{n}^{\prime}}\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{\Delta}\mathcal{P}_{\mathbf{n}^{\prime}}\tilde{\sigma}_{\Delta}(s)+\alpha^{2}\sum_{\mathbf{n}^{\prime}}\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{SE_{0}}\frac{1}{s-(\mathfrak{L}_{S,\mathbf{n}^{\prime}}+\mathfrak{L}_{E_{0},\mathbf{n}^{\prime}})}\mathfrak{L}_{SE_{0}}\mathcal{P}_{\mathbf{n}^{\prime}}\tilde{\sigma}_{\Delta}(s).$
(68)

Performing the inverse Laplace transform we recover

$\displaystyle\frac{d}{dt}\mathcal{P}_{\mathbf{n}}\sigma_{\Delta}(t)$
$\displaystyle=\sum_{\mathbf{n}^{\prime}}\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{\Delta}\mathcal{P}_{\mathbf{n}^{\prime}}\sigma_{{\Delta}(t)}+\alpha^{2}\sum_{\mathbf{n}^{\prime}}\int_{0}^{t}d\tau\mathcal{P}_{\mathbf{n}}\mathfrak{L}_{SE_{0}}e^{(\mathfrak{L}_{S,\mathbf{n}^{\prime}}+\mathfrak{L}_{E_{0},\mathbf{n}^{\prime}})(t-\tau)}\mathfrak{L}_{SE_{0}}\mathcal{P}_{\mathbf{n}^{\prime}}\sigma_{\Delta}(\tau)$
(69)
$\displaystyle=\sum_{\mathbf{n}^{\prime}}\operatorname{Tr}_{E_{0},\mathbf{n}}[\mathfrak{L}_{\Delta}\operatorname{Tr}_{E_{0},\mathbf{n}^{\prime}}[\sigma_{{\Delta}(t)}]\otimes\rho_{E_{0},\mathbf{n}^{\prime}}^{\beta}]\otimes\rho_{E_{0},\mathbf{n}}^{\beta}$
$\displaystyle\qquad\qquad\qquad+\alpha^{2}\sum_{\mathbf{n}^{\prime}}\int_{0}^{t}d\tau\operatorname{Tr}_{E_{0},\mathbf{n}}[\mathfrak{L}_{SE_{0}}e^{(\mathfrak{L}_{S,\mathbf{n}^{\prime}}+\mathfrak{L}_{E_{0},\mathbf{n}^{\prime}})(t-\tau)}\mathfrak{L}_{SE_{0}}\operatorname{Tr}_{E_{0},\mathbf{n}^{\prime}}[\sigma_{\Delta}(\tau)]\otimes\rho_{E_{0},\mathbf{n}^{\prime}}^{\beta}]\otimes\rho_{E_{0},\mathbf{n}}^{\beta}$

where we have used that
$\sum_{\mathbf{n}}\mathcal{P}_{\mathbf{n}}\sigma_{\Delta}=\sum_{\mathbf{n}}\operatorname{Tr}_{E_{0},\mathbf{n}}[\sigma_{\mathbf{n}}]\otimes\rho_{E_{0},\mathbf{n}}=\sum_{\mathbf{n}}\rho_{\mathbf{n}}\otimes\rho_{E_{0},\mathbf{n}}$
to obtain the second-order NZ/HEOM equation for each reduced ADO by tracing
over $E_{0}$ such that

$\displaystyle\frac{d}{dt}\rho_{\mathbf{n}}(t)=\mathcal{L}^{D}_{\mathbf{n}}\rho_{\mathbf{n}}(t)+\sum_{k}^{K}(\mathcal{L}^{+}_{\mathbf{n}}\rho_{\mathbf{n}_{k+}}(t)+\mathcal{L}^{-}_{\mathbf{n}}\rho_{\mathbf{n}_{K-}}(t))+\alpha^{2}\sum_{\mathbf{n}^{\prime}}\int_{0}^{t}d\tau\mathcal{K}^{(2)}_{\mathbf{n},\mathbf{n}^{\prime}}(t-\tau)\rho_{\mathbf{n}^{\prime}}(\tau).$
(70)

Here the first two terms on the RHS are the regular HEOM, the last term is the
second-order NZ correction for the influence of an additional environment,
with NZ memory Kernel

$\displaystyle\mathcal{K}^{(2)}_{\mathbf{n},\mathbf{n}^{\prime}}(\tau)=\operatorname{Tr}_{E_{0},\mathbf{n}}[\mathfrak{L}_{SE_{0},\mathbf{n}}e^{(\mathfrak{L}_{\mathbf{n}^{\prime}}+\mathfrak{L}_{E_{0},\mathbf{n}^{\prime}})\tau}\mathfrak{L}_{SE_{0},\mathbf{n}^{\prime}}\sigma_{E_{0}}^{\beta}].$
(71)

We note that this involves a summation over reduced ADOs from every tier,
however, the only term that mixes tiers is
$e^{\mathcal{L}_{\mathbf{n^{\prime}}}\tau}$, which acts, as one may expect, to
couple tiers of increasing distance as time evolves. For example, if we are to
expand this term as

$\displaystyle
e^{\mathfrak{L}_{\mathbf{n}}\tau}=\mathbf{1}_{\mathbf{n}}+\mathfrak{L}_{\mathbf{n}}\tau+\frac{1}{2}(\mathfrak{L}_{\mathbf{n}}\tau)^{2}+\cdots$
(72)

we see that each term couples tiers with an additional separation of a maximum
2 (one above, if available, and one below). In Ref. Fay and Limmer (2022) a
similar approach has been exploited, and a Redfield like term is obtained
making the Born-Markov approximation of the NZ kernel.

Now, we note that
$e^{\mathfrak{L}_{\mathbf{n}}\tau}\mathfrak{L}_{SE_{0},\mathbf{n}+\mathfrak{L}_{E_{0},\mathbf{n}^{\prime}}}=\tilde{\mathfrak{L}}_{SE_{0},\mathbf{n}}(\tau)$,
where the tilde denotes an interaction picture representation, where the
‘system’ is the uncoupled full hierarchy, and the environment is $E_{0}$.
Then,

$\displaystyle\mathcal{K}_{\mathbf{n},\mathbf{n}^{\prime}}(\tau)$
$\displaystyle=\operatorname{Tr}_{E_{0}}[\mathfrak{L}_{SE_{0},\mathbf{n}}\tilde{\mathfrak{L}}_{SE_{0},\mathbf{n}^{\prime}}(\tau)\sigma_{E_{0}}^{\beta}]$
$\displaystyle=-\operatorname{Tr}_{E_{0}}[[H_{SE_{0},\mathbf{n}},[\tilde{H}_{SE_{0},\mathbf{n}^{\prime}}(\tau),\sigma_{E_{0}}^{\beta}]]]$
(73)
$\displaystyle=-\operatorname{Tr}_{E_{0}}[[s_{E_{0},\mathbf{n}}B_{E_{0},\mathbf{n}},[\tilde{s}_{E_{0},\mathbf{n}^{\prime}}(\tau)\tilde{B}_{E_{0},\mathbf{n}^{\prime}}(\tau),\sigma_{E_{0}}^{\beta}]]].$

Then, defining the free environmental correlation functions
$C_{\mathbf{n},\mathbf{n}^{\prime}}^{(E_{0})}(\tau)=\operatorname{Tr}_{E_{0}}[\tilde{B}_{E_{0},\mathbf{n}^{\prime}}(\tau)B_{E_{0},\mathbf{n}^{\prime}}\sigma_{E_{0},\mathbf{n}^{\prime}}^{\beta}]$
we have for an interaction Hamiltonian $H_{SE,\mathbf{n}}=\alpha
s_{\mathbf{n}}B_{\mathbf{n}}$ (we note that the operators are identical for
each tier $\mathbf{n}$, though the label is nonetheless useful to keep track
of the ADO tiers on which they act)

$\displaystyle\mathcal{K}_{\mathbf{n},\mathbf{n}^{\prime}}(\tau)\rho_{\mathbf{n}^{\prime}}$
$\displaystyle=s_{\mathbf{n}}s_{\mathbf{n}^{\prime}}(t-\tau)\rho_{\mathbf{n}^{\prime}}(\tau)C_{\mathbf{n},\mathbf{n}^{\prime}}^{(E_{0})}(\tau-t)$
$\displaystyle-
s_{\mathbf{n}}\rho_{\mathbf{n}^{\prime}}(\tau)s_{\mathbf{n}^{\prime}}(t-\tau)C_{\mathbf{n}^{\prime},\mathbf{n}}^{(E_{0})}(t-\tau)$
(74) $\displaystyle-
s_{\mathbf{n}^{\prime}}(t-\tau)\rho_{\mathbf{n}^{\prime}}(\tau)s_{\mathbf{n}}C_{\mathbf{n},\mathbf{n}^{\prime}}^{(E_{0})}(\tau-t)$
$\displaystyle+\rho_{\mathbf{n}^{\prime}}(\tau)s_{\mathbf{n}^{\prime}}(t-\tau)s_{\mathbf{n}}C_{\mathbf{n}^{\prime},\mathbf{n}}^{(E_{0})}(t-\tau)$

We can then take the Markovian limit, where
$C_{\mathbf{n},\mathbf{n}^{\prime}}^{(E_{0})}(\tau)=\delta(\tau)\delta_{\mathbf{n},\mathbf{n}^{\prime}}$,
and extending upper limit of the time integral to infinity, we can write the
correction term
$\Gamma_{\mathbf{n},\mathbf{n}^{\prime}}=\int_{0}^{\infty}d\tau\mathcal{K}_{\mathbf{n},\mathbf{n}^{\prime}}(\tau)\rho_{\mathbf{n}^{\prime}}$
as

$\displaystyle\Gamma_{\mathbf{n},\mathbf{n}}(t)=s_{\mathbf{n}}s_{\mathbf{n}}\rho_{\mathbf{n}}(t)$
$\displaystyle-
s_{\mathbf{n}}\rho_{\mathbf{n}}(t)s_{\mathbf{n}}-s_{\mathbf{n}}\rho_{\mathbf{n}}(t)s_{\mathbf{n}}$
$\displaystyle+\rho_{\mathbf{n}}(t)s_{\mathbf{n}}s_{\mathbf{n}}$ (75)

Now, if $E_{0}$ is an optical environment coupled via, say
$s=\sigma_{x}=\sigma_{+}+\sigma_{-}$, we can take a rotating wave
approximation, and recover a GKSL type term for each tier:

$\displaystyle\Gamma_{\mathbf{n},\mathbf{n}}(t)=2\sigma_{-}\rho_{\mathbf{n}}(t)\sigma_{+}-\\{\sigma_{+}\sigma_{-},\rho_{\mathbf{n}}(t)\\}=:\mathcal{D}_{\mathbf{n}}[\rho_{\mathbf{n}}(t)]$
(76)

and we have the GKSL-HEOM of the form

$\displaystyle\partial_{t}\rho_{\mathbf{n}}(t)=(\mathcal{L}^{D}_{\mathbf{n}}$
$\displaystyle+\mathcal{D}_{\mathbf{n}})\rho_{\mathbf{n}}(t)$ (77)
$\displaystyle+\sum_{k}^{K}(\mathcal{L}^{+}_{\mathbf{n}}\rho_{\mathbf{n}_{k+}}(t)+\mathcal{L}^{-}_{\mathbf{n}}\rho_{\mathbf{n}_{k-}}(t))$

We thus see that the standard set of approximations that lead to a GKSL
equation yield non-additive dynamics with the HEOM - that is, the resulting
dissipators act on every tier of the ADOs.

We further note that in the most general form of the NZ-HEOM, Eq. (70), the
memory kernel may act to couple tiers of the hierarchy, and even ADOs within
the same tier not otherwise coupled via standard HEOM.

##  VIII Nakajima-Zwanzig Equation

Here, for completeness, we derive the NZ quantum master equation without the
presence of the non-perturbative environment described by HEOM. We will
largely follow references Breuer and Petruccione (2002); Gardiner and Zoller
(2000); Brian and Sun (2021), and refer the reader to Brian and Sun (2021) for
a more detailed review. We begin by writing the equations of motion for the
relevant and irrelevant parts of the total density matrix by applying the
projection operators directly to the LvN equation,

$\displaystyle\frac{d}{dt}\mathcal{P}\rho(t)=\mathcal{PLP}\rho(t)+\mathcal{PLQ}\rho(t)$
(78a)
$\displaystyle\frac{d}{dt}\mathcal{Q}\rho(t)=\mathcal{QLP}\rho(t)+\mathcal{QLQ}\rho(t)$
(78b)

where we have used $\mathcal{P+Q}=\mathbb{1}$, and that both projection
operators commute with the derivative. We than define the Laplace transform
$\tilde{f}(s)=\mathfrak{L}\\{f(t)\\}(s)=\int_{0}^{\infty}ste^{-st}f(t)$, which
has the useful properties,

$\displaystyle\mathfrak{L}\left\\{\frac{df(t)}{dt}\right\\}(s)=\int_{0}^{\infty}e^{-st}\frac{df(t)}{dt}dt=s\tilde{f}(s)-f(0)$
$\displaystyle\mathfrak{L}\\{g(t)\star f(t)\\}=\tilde{g}(s)\tilde{f}(s),$ (79)

where we have defined the convolution $g(t)\star f(t)=\int_{0}^{\infty}d\tau
g(t-\tau)f(\tau)$. Using the first of these properties we can take the Laplace
transforms of Eq. (78), obtaining

$\displaystyle
s\mathcal{P}\tilde{\rho}(s)-\mathcal{P}\rho(0)=\mathcal{PLP}\tilde{\rho}(s)+\mathcal{PLQ}\tilde{\rho}(s)$
(80a) $\displaystyle
s\mathcal{Q}\tilde{\rho}(s)-\mathcal{Q}\rho(0)=\mathcal{QLP}\tilde{\rho}(s)+\mathcal{QLQ}\tilde{\rho}(s).$
(80b)

Rearranging Eq. (80b) we obtain

$\displaystyle\mathcal{Q}\tilde{\rho}(s)=\frac{1}{s-\mathcal{QL}}\mathcal{Q}\left(\rho(0)+\mathcal{LP}\tilde{\rho}(s)\right)$
(81)

which we can substitute into Eq. (80a) to find

$\displaystyle s\mathcal{P}\tilde{\rho}(s)-\mathcal{P}$
$\displaystyle\rho(0)=\mathcal{PLP}\tilde{\rho}(s)$
$\displaystyle+\mathcal{PL}\frac{1}{s-\mathcal{QL}}\mathcal{Q}\left(\rho(0)+\mathcal{LP}\tilde{\rho}(s)\right).$
(82)

We then take the inverse Laplace transform of this expression, which for the
LHS and first term of the RHS is trivial (as they are themselves the result of
a Laplace transform above), and for the second term of the RHS requires use of
the convolution property in Eq. (VIII). We thus obtain

$\displaystyle\frac{d}{dt}\mathcal{P}\rho(t)$
$\displaystyle=\mathcal{PLP}\rho(t)+\mathcal{PL}e^{\mathcal{QL}t}\mathcal{Q}\rho(0)$
$\displaystyle\qquad+\int_{0}^{t}d\tau\mathcal{PL}e^{\mathcal{QL}(t-\tau)}\mathcal{QLP}\rho(\tau)$
(83)
$\displaystyle=\mathcal{PLP}\rho(t)+\mathcal{I}(t)+\int_{0}^{t}d\tau\mathcal{K}(t-\tau)\mathcal{P}\rho(\tau).$

Where we have defined the inhomogeneous term $\mathcal{I}(t)$, which vanishes
is $\mathcal{P}\rho(0)=\rho(0)$, and the memory kernel $\mathcal{K}(t)$.

We note that up to this point we have made no reference to the form of the
projection operators $\mathcal{P},\,\mathcal{Q}$, and thus the result is
equally applicable to the HEOM case. First, however, we simplify the result
using the more standard projection operators defined by
$\mathcal{P}\rho=\operatorname{Tr}_{E}[\rho]\otimes\rho_{E}$, and
$\mathcal{Q}=\mathbf{1}-\mathcal{P}$. In this case, we can further simplify
the problem,

$\displaystyle\frac{d}{dt}\rho_{S}(t)=\mathcal{L}_{S}\rho_{S}(t)+\int_{0}^{t}d\tau\mathcal{K}_{S}(t-\tau)\rho_{S}(\tau),$
(84)

where
$\mathcal{K}_{S}(t):=\operatorname{Tr}_{E}[\mathcal{PL}e^{\mathcal{QL}t}\mathcal{QLP}]$.

###  VIII.1 Weak coupling limit

We now wish to take the weak coupling and Markovian limits of Eq. (84), which
is straightforwardly obtained by once again looking at the Laplace transformed
equations for the relevant and irrelevant subspaces in Eq. (80). We then
redefine $\mathcal{L}=\mathcal{L}_{S}+\mathcal{L}_{E}+\alpha\mathcal{L}_{I}$,
such that the parameter $\alpha$ determines the coupling strength to the
environment. We can then use that
$\mathcal{PL}_{S}=\mathcal{L}_{S}\mathcal{P}$,
$\mathcal{PL}_{E}=\mathcal{L}_{E}\mathcal{P}=0$, and
$\mathcal{P}\mathcal{L}_{I}\mathcal{P}=0$ to obtain

$\displaystyle\mathcal{PLP}=\mathcal{P\mathcal{L}_{S}P}$ (85a)
$\displaystyle\mathcal{PLQ}=\alpha\mathcal{PL}_{I}\mathcal{Q}$ (85b)
$\displaystyle\mathcal{QLP}=\alpha\mathcal{QL}_{I}\mathcal{P}$ (85c)