Dynamics of photosynthetic light harvesting systems interacting with N-photon Fock states

We model the lowest two accessible electronic states of each chromophore as a two level system, corresponding to the Q_y transition [6]. The dipole-dipole coupling between chromophores, under the rotating wave approximation, does not change the number of electronic excitations. Therefore the chromophoric Hamiltonian will commute with the operator that counts the number of excitons in the system, and the Hamiltonian is block diagonal, with the j^th block corresponding to the j-excitation subspace, $\mathcal{H}_{j}$. The excitation energy of a typical chromophore is $\sim 15,000\,\text{cm}^{-1}$, about 75 times larger than $k_{B}T\approx 200\text{ cm}^{-1}$ at room temperature. Therefore, to a very good approximation, we may regard the initial thermal state of the isolated chromophoric system as the absolute ground state of $\mathcal{H}_{0}$, denoted as $|g\rangle$, which is defined to have zero energy. Due to the weak system-field coupling in a natural photosynthetic system, the probability to have multiple excitations in the system is much smaller than the probability to have a single excitation, so we will only consider the $(1+N)$-dimensional subspace $\mathcal{H}_{0}\bigoplus\mathcal{H}_{1}$. We denote the singly excited state where site $j$ is excited and all other sites are in their ground states as $|j\rangle$. The chromophoric system Hamiltonian is written as

$$H_{\text{sys}}=\sum^{N}_{j=1}\epsilon_{j}|j\rangle\langle j|+\sum_{j\neq k}J_{jk}|j\rangle\langle k|,$$

(1)

where $\epsilon_{j}$ is the excited state energy of site $j$, and $J_{jk}$ is the dipole-dipole interaction between chromophores $j$ and $k$. The effects of the phonon environment on the site energies are discussed in Section 2.8. The numerical values for the $\epsilon_{j}$ and $J_{jk}$ parameters in an LHCII monomer system are taken from [18]. For simplicity, we shall refer to the chromophore system as the system, unless otherwise noted.

The quantized electromagnetic field in 3-dimensional space is described by harmonic oscillators indexed with a 3-dimensional wavevector and a polarization index [19, 20]. However, it is useful to decompose the electric field into 1-dimensional fields, so that the input-output formalism can be applied [11, 21]. We will first consider the incoming N-photon Fock state in a fixed paraxial beam mode. The electric field at the center of the beam is written as

$$\mathbf{E}_{\text{para}}(t)=\int_{0}^{\infty}d\omega\,\sqrt{\frac{\hbar\omega}{4\pi\epsilon_{0}c}}\tilde{\mathbf{u}}(\mathbf{0})(ia(\omega)e^{-i\omega t}+\text{h.c.}),$$

(2)

(see derivation in appendix A or Ref. [22]), where $\tilde{\mathbf{u}}(\mathbf{x})$ is the normalized spatial mode function such that the integral over all transverse area is unity, i.e., $\int dA_{T}\,|\tilde{\mathbf{u}}|^{2}=1$. Here $a(\omega)$ is the annihilation operator for the frequency $\omega$ component of the paraxial mode. Note that the field is now indexed by the 1-dimensional parameter $\omega$.

To capture the spontaneous emission into other modes, we must also consider field modes other than the incoming mode. One way to do this is to partition the $4\pi$ solid angle of the 3-dimensional wavevector into finite numbers of small solid angle sections, indexed by $m$, and write the electric field at position $\mathbf{x}=\mathbf{0}$ as

	$$\mathbf{E}(t)=\sum_{m,\lambda}\mathbf{E}_{m,\lambda}(t),$$		(3a)
where
	$$\mathbf{E}_{m,\lambda}(t)=\int_{0}^{\infty}d\omega\,\sqrt{\frac{\hbar\omega^{3}\Delta\Omega_{m}}{16\pi^{3}\epsilon_{0}c^{3}}}(ia_{m,\lambda}(\omega)e^{-i\omega t}\hat{\epsilon}_{m,\lambda}+\text{h.c.})$$		(3b)

(see appendix B). Here $\lambda$ indexes the two polarizations in a small solid angle section, $\Delta\Omega_{m}$ is the amount of solid angle (in steradian units) of section $m$, and $\hat{\epsilon}_{m,\lambda}$ is the polarization vector corresponding to $m$ and $\lambda$. One can show that in the case of the simplest transverse electromagnetic mode TEM₀₀,if we think of the boundary of the beam as the location where the beam intensity is $1/e^{4}\approx 2\%$ of the intensity at the center, then Eq. (2) matches Eq. (3b) for a particular $(m,\lambda)$. Therefore we can treat the incoming TEM₀₀ paraxial mode as one of the $(m,\lambda)$ modes in the small angle decomposition (Eq. (3b)), as illustrated in Figure (1).

In Eq. (3b), the electric field is decomposed into a finite number of 1-D fields, which is not enough to describe all degrees of freedom in the 3-D electromagnetic field. However, appendix C shows that one can describe the 3-D electromagnetic field by a set of countably infinite 1-D fields. Then using the small solid angle decomposition (see Eq. (3b) and appendix B), we can decompose the 3-D electromagnetic field into a finite number of small solid angle 1-D fields plus a countably infinite number of 1-D fields that are needed to describe all degrees of freedom in the 3-D field. The countably infinite number of 1-D fields are chosen to be orthogonal to each other and to the small solid angle 1-D fields in the sense described in appendix C. Under this decomposition scheme, the total system+field Hamiltonian under the dipole-electric field coupling ($-\mathbf{d}\cdot\mathbf{E}$) can be written as

$$H_{\text{sys+field}}=H_{\text{sys}}-\mathbf{d}\cdot\sum_{m,\lambda}\mathbf{E}_{m,\lambda}+\sum_{m,\lambda}\int^{\infty}_{0}d\omega\,\hbar\omega a_{m,\lambda}^{\dagger}(\omega)a_{m,\lambda}(\omega)+\sum_{s}^{\infty}\int^{\infty}_{0}d\omega\,\hbar\omega a_{s}^{\dagger}(\omega)a_{s}(\omega),$$

(4)

where $\sum_{m,\lambda}$ sums over the finite number of small solid angle modes, and the $\sum_{s}^{\infty}$ sums over the remaining countably infinite number of 1-D field modes. Since the last term in Eq. (4) involving the infinite sum is decoupled from the rest of the Hamiltonian, we can describe the system-light interaction as system interacting with just a finite number of 1-D fields.

An alternative way to decompose the electric field into a finite sum of 1-D fields is to write it as a sum of the x, y, and z polarization components, i.e. $\mathbf{E}(t)=E_{x}(t)\hat{x}+E_{y}(t)\hat{y}+E_{z}(t)\hat{z}$. $E_{x}$ then takes the form

$$E_{x}(t)=\int^{\infty}_{0}d\omega\,\sqrt{\frac{\hbar\omega^{3}}{6\pi^{2}\epsilon_{0}c^{3}}}(ia_{x}(\omega)e^{-i\omega t}+\text{h.c.}),$$

(5)

with $E_{y}$ and $E_{z}$ defined similarly (see appendix D).

To put the system-light interaction in the language of input-output formalism, we follow the procedure in [12], with addition of some details specific to our modeling of light harvesting systems. Since the size of light harvesting complexes is much smaller than the wavelength of the light they interact with, we treat all chromophores as being in the same location. The dipole operator $\mathbf{d}$ of the system is written as

$$\mathbf{d}=\sum_{j=1}^{N}(\tilde{\mathbf{d}}_{j}\sigma^{(j)}_{-}+\text{h.c.}),$$

(6)

where $\tilde{\mathbf{d}}_{j}=\langle g|\mathbf{d}|j\rangle$ is the transition dipole moment of site $j$, and $\sigma^{(j)}_{-}\equiv|g\rangle\langle j|$ is the lowering operator of site $j$. The total system+field Hamiltonian is written as

$$H_{\text{sys+field}}=H_{\text{sys}}+\frac{\hbar}{\sqrt{2\pi}}\sum_{l}\int^{\infty}_{0}d\omega\,\big{(}-ia_{l}(\omega)L_{l}^{\dagger}(\omega)+\text{h.c.}\big{)}+\sum_{l}\int^{\infty}_{0}d\omega\,\hbar\omega a_{l}^{\dagger}(\omega)a_{l}(\omega),$$

(7)

where

$$L_{l}(\omega)=\sqrt{\frac{d_{0}^{2}\omega^{3}\Delta\Omega}{8\pi^{2}\epsilon_{0}\hbar c^{3}}}\sum_{j=1}^{N}\mathbf{d}_{j}\cdot\hat{\epsilon}\,\sigma^{(j)}_{-},$$

(8)

for a small solid angle mode (Eq. (3b)). For notational simplicity, we denote the mode indices $(m,\lambda)$ by $l$. The dipole moment $\tilde{\mathbf{d}}_{j}$ is expressed as a scalar unit dipole $d_{0}$ multiplying a dimensionless vector dipole moment $\mathbf{d}_{j}$. In Eq. (7), we have used a rotating wave approximation to drop terms that do not conserve excitation numbers, i.e. $a^{\dagger}(\omega)L^{\dagger}(\omega)$ and $a(\omega)L(\omega)$. The constant factor $\frac{\hbar}{\sqrt{2\pi}}$ is chosen for later convenience. For a paraxial mode (Eq. (2)), the prefactor in Eq. (8) is $\sqrt{\frac{d_{0}^{2}\omega}{2\epsilon_{0}\hbar c}}\tilde{u}(\mathbf{0})$. For a polarization mode (Eq. (5)), the prefactor is $\sqrt{\frac{d_{0}^{2}\omega^{3}}{3\pi\epsilon_{0}\hbar c^{3}}}$. In the current work one of the field modes indexed by $l$ in Eq. (7) is the incoming paraxial mode containing the N-photon Fock state pulse. In the following, we shall denote quantities pertaining to the incoming mode by replacing the subscript “$l$” with subscript “inc”.

Since the system is resonant with the field in a narrow frequency range near the pulse center frequency $\omega_{0}$, only the frequency near $\omega_{0}$ contributes significantly to the overall dynamics. We make the approximation $L_{l}(\omega)\approx L_{l}(\omega_{0})$, and simply denote it as $L_{l}$ hereafter. In order to have a unified approach to describe different ways to decompose the electric field, it is useful to express the prefactor in Eq. (8) as $\sqrt{\Gamma_{0}\eta}$, where $\Gamma_{0}=\frac{d_{0}^{2}\omega_{0}^{3}}{3\pi\epsilon_{0}\hbar c^{3}}$ is the spontaneous emission rate of a unit dipole $d_{0}$ with transition frequency $\omega_{0}$ in vacuum, and $\eta$ is a geometric factor depending on the type of the field mode. For a paraxial mode (Eq. (2)), $\eta=\frac{3\pi c^{2}}{2\omega_{0}^{2}}|\tilde{u}(\mathbf{0})|^{2}$; for a small solid angle mode (Eq. (3b)), $\eta=\frac{3}{8\pi}\Delta\Omega$; for a polarization mode (Eq. (5)), $\eta=1$.

To summarize, $L_{l}$ is now written as

$$L_{l}=\sqrt{\Gamma_{0}\eta_{l}}\sum_{j}\mathbf{d}_{j}\cdot\hat{\epsilon}|g\rangle\langle j|.$$

(9)

To account for the dielectric medium, we will replace the spontaneous emission rate in vacuum with the spontaneous emission rate in the empty cavity model with index of refraction $n$, given by

$$\Gamma_{0}=n\left(\frac{3n^{2}}{2n^{2}+1}\right)^{2}\frac{d_{0}^{2}\omega_{0}^{3}}{3\pi\epsilon_{0}\hbar c^{3}}$$

(10)

[23, 24]. We take $n=1.33$ as the index of refraction of water for the numerical simulations.

$L_{l}$ can be simplified further by writing

$$L_{l}=\sqrt{\Gamma_{l}}|g\rangle\langle B_{l}|,$$

(11)

with $\Gamma_{l}$ the effective coupling constant

$$\Gamma_{l}=\Gamma_{0}\eta_{l}\sum_{j}|\mathbf{d}_{j}\cdot\hat{\epsilon}|^{2}$$

(12)

and $|B_{l}\rangle$ the normalized bright state

$$|B_{l}\rangle=(\sum_{j}|\mathbf{d}_{j}\cdot\hat{\epsilon}|^{2})^{-1/2}\sum_{j}\mathbf{d}_{j}\cdot\hat{\epsilon}|j\rangle.$$

(13)

As a physical motivation for this notation, we note that the excited state $|B_{l}\rangle$ spontaneously emits into mode $l$ at the rate $\Gamma_{l}$. A more detailed discussion of the emission rate, or more generally, the photon flux, is given in Section 2.7.

Going into an interaction frame where we rotate out the zeroth order non-interacting Hamiltonian

$$H_{0}=\sum^{N}_{j=1}\hbar\omega_{0}|j\rangle\langle j|+\sum_{l}\int^{\infty}_{0}d\omega\,\hbar\omega a_{l}^{\dagger}(\omega)a_{l}(\omega),$$

(14)

the interaction frame Hamiltonian becomes

$$H_{\text{int}}(t)=\underbrace{\bigg{(}\sum^{N}_{j=1}\Delta_{j}|j\rangle\langle j|+\sum_{j\neq k}J_{jk}|j\rangle\langle k|\bigg{)}}_{H}+\frac{\hbar}{\sqrt{2\pi}}\sum_{l}\int^{\infty}_{0}d\omega\,\big{(}-ia_{l}(\omega)e^{-i(\omega-\omega_{0})t}L_{l}^{\dagger}+\text{h.c.}\big{)},$$

(15)

where $\Delta_{j}\equiv\epsilon_{j}-\hbar\omega_{0}$ is the site energy detuning (see Eq. (1)). The term inside the first parenthesis acts on the system only and does not depend on time. We shall denote this simply as $H$.

Changing to the relative frequency variable $\omega^{\prime}=\omega-\omega_{0}$ and re-indexing the frequency components $a(\omega)\rightarrow a(\omega-\omega_{0})$, Eq. (15) becomes

$$H_{\text{int}}(t)=H+\frac{\hbar}{\sqrt{2\pi}}\sum_{l}\int^{\infty}_{-\infty}d\omega^{\prime}\,\big{(}-ia_{l}(\omega^{\prime})e^{-i\omega^{\prime}t}L_{l}^{\dagger}+\text{h.c.}\big{)}.$$

(16)

Since the dynamics is dominated by the field modes in a narrow bandwidth near resonance ($\omega^{\prime}\approx 0$) and $\omega_{0}$ is much larger than this bandwidth, we let the lower bound of the integral, $-\omega_{0}$, go to $-\infty$.

Finally, we define the quantum white noise operator

$$a_{l}(t)\equiv\sqrt{\frac{1}{2\pi}}\int^{\infty}_{-\infty}d\omega^{\prime}\,a_{l}(\omega^{\prime})e^{i\omega^{\prime}t},$$

(17)

a central object in the input-output formalism. With this definition, the field operators $a_{l}(t)$ satisfy the bosonic commutation relations: $[a_{l}(t),a_{l^{\prime}}^{\dagger}(t^{\prime})]=\delta_{l,l^{\prime}}\delta(t-t^{\prime})$ and $[a_{l}(t),a_{l^{\prime}}(t^{\prime})]=[a_{l}^{\dagger}(t),a_{l^{\prime}}^{\dagger}(t^{\prime})]=0$. The system-light interaction (Eq. (16)) can now be written in a simple-looking form as

$$H_{\text{int}}(t)=H+\sum_{l}\left(-i\hbar a_{l}(t)L_{l}^{\dagger}+i\hbar a_{l}^{\dagger}(t)L_{l}\right).$$

(18)

An N-photon Fock state pulse in mode $l$ with temporal profile $\xi(t)$ is defined as

$$|N_{\xi}\rangle_{l}=\frac{1}{\sqrt{N!}}\bigg{(}\int d\tau\,\xi(\tau)a_{l}^{\dagger}(\tau)\bigg{)}^{N}|\phi\rangle,$$

(19)

where $\xi(t)$ is normalized according to $\int d\tau\,|\xi(\tau)|^{2}=1$ and $|\phi\rangle$ is the vacuum state of the field. More general N-photon Fock states are possible, as discussed in Ref. [9], but we shall restrict ourselves to this single mode product form here.

Another important object in the input-output formalism is the unitary $U(t)$ defined by

$$\frac{dU(t)}{dt}=-\frac{i}{\hbar}H_{\text{int}}(t)U(t)$$

(20)

with initial condition $U(0)=1$. The system+field unitary in the Schrodinger picture, $e^{-\frac{i}{\hbar}H_{\text{sys+field}}t}$, is related to the interaction picture unitary $U(t)$ by

$$e^{-\frac{i}{\hbar}H_{\text{sys+field}}t}=e^{-\frac{i}{\hbar}H_{0}t}U(t)$$

(21)

(see Eq. (14)). We will set $\hbar=1$ from now on.

This section summarizes some results in the input-output formalism that will be used later. The operators $a(t)$ and $U(t)$ satisfy the so-called input-output relation [11]:

$$U(t^{\prime})a_{l}(t)U^{\dagger}(t^{\prime})=\begin{cases}a_{l}(t)\quad t^{\prime}<t\\ a_{l}(t)+\frac{1}{2}L_{l}(t)\quad t^{\prime}=t\\ a_{l}(t)+L_{l}(t)\quad t^{\prime}>t,\end{cases}$$

(22)

with $L_{l}(t)\equiv U^{\dagger}(t)L_{l}U(t)$ (see appendix E). At time $t^{\prime}<t$, $a_{l}(t)$ has not interacted with the system, so time evolution has no effect on it (i.e. $U(t^{\prime})a_{l}(t)U^{\dagger}(t^{\prime})=a_{l}(t)$). For this reason, $a_{l}(t)$ is called the input field. At time $t^{\prime}>t$, $a(t)$ has interacted with the system and will not interact with the system anymore. For this reason, the time-evolved operator $U(t^{\prime})a_{l}(t)U^{\dagger}(t^{\prime})=a_{l}(t)+L_{l}(t)\equiv a_{l,\text{out}}(t)$ is called the output field. Using the Markov property that $a(t)$ only interacts with the system at time $t$, we have thereby expressed the time-evolved field operator, which usually mixes the system and field degrees of freedom in some complicated way, as a simple sum of a pure field operator and a time-evolved system operator.

Left-multiplying Eq. (22) by $U(t)$, we then obtain the following commutation relation at $t=t^{\prime}$,

$$[a_{l}(t),U(t)]=\frac{1}{2}L_{l}U(t).$$

(23)

Using this commutation relation, we can rewrite Eq. (20) (also see Eq. (18)) in normal-ordered form as

$$\frac{d}{dt}U(t)=\big{(}-iH-\frac{1}{2}\sum_{l}L_{l}^{\dagger}L_{l}\big{)}U(t)-\sum_{l}L_{l}^{\dagger}U(t)a_{l}(t)+\sum_{l}a_{l}^{\dagger}(t)L_{l}U(t).$$

(24)

In the following sections, we will see that the commutation relation and the normal-ordered form of the time evolution derivative ensure that we only need to consider the action of normal-ordered field operators acting on field states, greatly simplifying the calculations.

Input-output theory is traditionally presented in terms quantum stochastic differential equations (QSDE) where the input field operators are taken to be differentials, e.g. $dA_{l,t}=a_{l}(t+dt)-a_{l}(t)$, analogous to the Wiener process increment $dW_{t}$ in classical stochastic calculus. As in classical stochastic calculus, integrals over the quantum stochastic differential $\int X_{t}dA_{t}$ can be interpreted either as Stratonivich integrals or as Ito integrals, resulting in different values [25, 26]. The Stratonovich integral takes a “mid-point” time approximation for the integrand, while the Ito integral takes the initial time approximation for the integrand. To make connections to the QSDE approach, we note that Eq. (24) can be written as

$$dU_{t}=-iHU_{t}dt-\sum_{l}L_{l}^{\dagger}U_{t}\circ dA_{l,t}+\sum_{l}L_{l}U_{t}\circ dA_{l,t}^{\dagger},$$

(25)

where the symbol $\circ$ indicates that the differential expression is to be integrated in the Stratonovich sense, or it can also be written as

$$dU_{t}=\big{(}-iH-\frac{1}{2}\sum_{l}L_{l}^{\dagger}L_{l}\big{)}U_{t}dt-\sum_{l}L_{l}^{\dagger}U_{t}\cdot dA_{l,t}+\sum_{l}L_{l}U_{t}\cdot dA_{l,t}^{\dagger},$$

(26)

where the symbol $\cdot$ indicates that the differential expression is to be integrated in the Ito sense.

It was shown in [27] that a usual time-ordered product of iterated time integrals that include products of the field operators $a_{l}(t)$ and $a_{l}^{\dagger}(t)$ should be interpreted as a quantum Stratonovich integral. Conversely, if the iterated time integrals were instead written in normal-order then that corresponds to a quantum Ito integral. Thus by writing Eq. (24) in normal order, the solution will have the same mathematical properties as a quantum Ito integral (Eq. (26)), without delving into the mathematics of quantum Ito integration and its modified rules of calculus.

Let $|N_{\xi}\rangle$ be the field state where a single photon pulse with temporal profile $\xi(t)$ is in the incoming mode ($l=\text{inc}$), and all other field modes are in the vacuum state. The system density matrix in the interaction picture, $\tilde{\rho}_{\text{sys}}(t)$, is obtained by evolving the initial state then partially tracing over the field, i.e.,

$$\tilde{\rho}_{\text{sys}}(t)=\text{Tr}_{\text{field}}\Big{(}U(t)\rho_{0}\otimes|N_{\xi}\rangle\langle N_{\xi}|U^{\dagger}(t)\Big{)}.$$

(27)

We assume at $t=0$, the system+field state is a factorizable state, $\rho_{0}\otimes|N_{\xi}\rangle\langle N_{\xi}|$. Using Eqs. (23) and (24), it can be shown that the reduced system dynamics follow a hierarchy of equations [9] (see appendix F)

$$\begin{split}\frac{d\rho_{m,n}}{dt}=&-i[H,\rho_{m,n}]+\sum_{l}\mathcal{D}[L_{l}](\rho_{m,n})\\ &+\sqrt{m}\xi(t)[\rho_{m-1,n},L_{\text{inc}}^{\dagger}]+\sqrt{n}\xi^{*}(t)[L_{\text{inc}},\rho_{m,n-1}],\end{split}$$

(28)

where $\mathcal{D}[L]$ is the Lindblad superoperator defined as $\mathcal{D}[L](\rho)\equiv-\frac{1}{2}L^{\dagger}L\rho-\frac{1}{2}\rho L^{\dagger}L+L\rho L^{\dagger}$ and the density matrices $\rho_{m,n}$ are defined as

$$\rho_{m,n}(t)\equiv\text{Tr}_{\text{field}}\Big{(}U(t)\rho_{0}\otimes|m_{\xi}\rangle\langle n_{\xi}|U^{\dagger}(t)\Big{)}.$$

(29)

Given an N-photon input, $\rho_{N,N}(t)=\tilde{\rho}_{\text{sys}}(t)$ is the physical density matrix that describes the system state. This couples to auxiliary density matrices $\rho_{m,n}$ with lower indices, where $0\leq m,n\leq N$. When solving the equations without phonons, we can take advantage of the fact that $\rho_{m,n}=\rho_{n,m}^{\dagger}$. At $t=0$, the initial condition is $\rho_{m,n}=\delta_{m,n}\rho_{0}$. The Fock state master equation (Eq. (28)) was originally derived by Baragiola et. al. [9] using the more mathematical quantum stochastic differential equations (QSDE). In appendix F, we give a more accessible alternative derivation based on ordinary differential equations (ODE).

In the special case of a single photon Fock state input, the Fock state master equations (Eq. (28)) can be solved analytically. However, rather than solving Eq. (28) directly, a more physically intuitive approach is to write the system+field pure state $|\psi(t)\rangle$ as

$$|\psi(t)\rangle=|\beta(t)\rangle|\text{vac}\rangle+|g\rangle\sum_{l}\int^{\infty}_{-\infty}dt_{r}\,\phi_{l}(t,t_{r})a_{l}^{\dagger}(t_{r})|\text{vac}\rangle,$$

(30)

where $|\beta(t)\rangle$ is an unnormalized system state in the excited subspace, and $\phi_{l}(t,t_{r})$ is an unnormalized “wavefunction” of the single photon field state in mode $l$ at time $t$. This is in fact the most general form of the system+field pure state when there is exactly one excitation in the system and the field. We can solve for $|\beta(t)\rangle$ and $\phi_{l}(t,t_{r})$ to obtain (see appendix G)

	$\displaystyle|\beta(t)\rangle$	$\displaystyle=-\int^{t}_{0}d\tau\,\xi(\tau)e^{(-iH-\frac{1}{2}\sum_{l}L^{\dagger}_{l}L_{l})(t-\tau)}L^{\dagger}_{\text{inc}}|g\rangle$		(31a)
	$\displaystyle\phi_{l}(t,t_{r})$	$\displaystyle=\begin{cases}\delta_{l,\text{inc}}\xi(t_{r})\quad,t<t_{r}\\ \delta_{l,\text{inc}}\xi(t_{r})+\frac{1}{2}\langle g|L_{l}|\beta(t_{r})\rangle\quad,t=t_{r}\\ \delta_{l,\text{inc}}\xi(t_{r})+\langle g|L_{l}|\beta(t_{r})\rangle\quad,t>t_{r}\end{cases}$		(31b)

Since $\langle\psi(t)|\psi(t)\rangle=1$, we have the property

$$\langle\beta(t)|\beta(t)\rangle+\sum_{l}\int dt_{r}\,|\phi_{l}(t,t_{r})|^{2}=1.$$

(32)

Using Eq. (29), we find that the solution to the single photon Fock state master equation (Eq. (28)) is

$$\begin{split}&\rho_{1,1}(t)=|\beta(t)\rangle\langle\beta(t)|+|g\rangle\langle g|\sum_{l}\int^{\infty}_{-\infty}dt_{r}\,|\phi_{l}(t,t_{r})|^{2}\\ &\rho_{1,0}(t)=|\beta(t)\rangle\langle g|\\ &\rho_{0,0}(t)=|g\rangle\langle g|.\end{split}$$

(33)

The total probability of being in the excited state is $\langle\beta(t)|\beta(t)\rangle$. However, since the system energy scale $||H||$ is much larger than the spontaneous emission rate $||\sum_{l}L^{\dagger}_{l}L_{l}||$, if we are interested in the system behavior at short times we will have $||\sum_{l}L^{\dagger}_{l}L_{l}||\,t\ll 1$. It is then useful to drop the spontaneous emission terms in Eq. (31a) and approximate $|\beta(t)\rangle$ as

$$|\beta_{\xi}^{\prime}(t)\rangle\equiv-\int^{t}_{0}d\tau\,\xi(\tau)e^{-iH(t-\tau)}L^{\dagger}_{\text{inc}}|g\rangle,$$

(34)

where for later convenience, we added the subscript $\xi$ to emphasize the dependence of the excited state $|\beta^{\prime}(t)\rangle$ on the temporal profile $\xi(t)$.

The photon flux of a field mode is defined as the rate at which photons pass through the mode and has the dimension of 1/[time]. At time $t$, the photon flux immediately downstream of the system is given by $a^{\dagger}_{\text{out}}(t)a_{\text{out}}(t)$ [9, 19].

The expectation value of the photon flux $F_{l}$ in spatial mode $l$ is the trace over the initial state:

$$F_{l}=\text{Tr}\Big{(}a^{\dagger}_{l,\text{out}}(t)a_{l,\text{out}}(t)\rho_{0}\otimes|N_{\xi}\rangle\langle N_{\xi}|\Big{)}.$$

(35)

Substituting $a_{l,\text{out}}(t)=a_{l}(t)+L_{l}(t)$ in Eq. (22) and following a similar procedure as in appendix F, we obtain the following explicit expression for the photon flux into mode $l$

$$F_{l}=\begin{cases}N|\xi(t)|^{2}+\big{(}\sqrt{N}\xi^{*}(t)\langle L_{l}\rangle_{N,N-1}+\text{c.c.}\big{)}+\langle L_{l}^{\dagger}L_{l}\rangle_{N,N}\quad,\,l=\text{inc}\\ \langle L_{l}^{\dagger}L_{l}\rangle_{N,N}\quad,\,\text{otherwise},\end{cases}$$

(36)

where $\langle X\rangle_{n,m}\equiv\text{Tr}(X\rho_{n,m})$. In the expression on the first line for the incoming mode, the first term represents the transmission of the incoming photon and the second term arises from the coherent dynamics between the system and the field. A negative value of the latter term represents the absorption of photons. A positive value for this term corresponds to stimulated emission. It is interesting to note that the second term can be positive also in the case of a single photon input, meaning that a single photon can stimulate its own emission. The final term on the first line for the incoming mode has the same form as the flux expression for the non-incoming modes and represents the spontaneous emission into the particular mode $l$.

We model the interaction with the phonon using the shifted harmonic oscillator model, where the excited state potential energy surface consists of harmonic potentials shifted from the ground state potential energy surface [28, 1]. The overall Hamiltonian can be written as

$$H_{\text{sys+vib}}=\sum_{j}\epsilon_{j}|j\rangle\langle j|+\sum_{j\neq k}J_{jk}|j\rangle\langle k|+H_{\text{vib}}+\sum_{j}|j\rangle\langle j|u_{j},$$

(37)

where

$$H_{\text{vib}}=\sum_{j,\xi}\frac{p^{2}_{j,\xi}}{2}+\frac{\omega^{2}_{j,\xi}q^{2}_{j,\xi}}{2}$$

(38)

is the vibrational Hamiltonian in the electronic ground state, and

$$u_{j}=\sum_{\xi}c_{j,\xi}q_{j,\xi}.$$

(39)

$\xi$ indexes the phonon modes, and $p_{j,\xi}$ and $q_{j,\xi}$ are the momentum and position of the phonon mode $\xi$ on site $j$. $u_{j}$ is a collective phonon coordinate characterized by the spectral density $J(\omega)$, defined as

$$J(\omega)=\frac{\pi}{2}\sum_{\xi}\frac{c_{j,\xi}^{2}}{\omega_{j,\xi}}\delta(\omega-\omega_{j,\xi}).$$

(40)

Since the system-vibration coupling term does not involve the system ground state, to a very good approximation, the initial thermal state $\propto e^{-\beta H_{\text{sys+vib}}}$ is a product state between the system ground state $|g\rangle\langle g|$ and the vibrational thermal state $\propto e^{-\beta H_{\text{vib}}}$. For the rest of Section 2, we analyze the effect of vibration using two different methods. The first method considers an initial vibrational pure state, and solves for the pure state analogous to Section 2.6. Then by averaging the dynamics starting from a thermal distribution of vibrational pure states, we can analyze the behavior given an initial vibrational thermal state. This method can be applied numerically only for a small number of discrete vibration modes, but it gives us useful analytical expressions in the continuum limit. The second method uses the HEOM formalism to trace out the vibration degrees of freedom and represent the effect of a continuum of vibration modes by a set of auxiliary density matrices containing only the system degrees of freedom.

Generalizing the analysis presented in Section 2.6, one can write down a system+field+vibration pure state as a function of time, given the system+vibration state initialized in the pure product state $|g\rangle\otimes|v\rangle$, where $|v\rangle$ is an arbitrary vibration state. In Section 4, we shall take $|v\rangle$ to be the eigenstate of $H_{\text{vib}}$, the vibrational Hamiltonian in the electronic ground state with energy $E_{v}$, so that $H_{\text{vib}}|v\rangle=E_{v}|v\rangle$. The overall pure state $|\psi(t)\rangle$ is written as

$$|\psi(t)\rangle=|\gamma(t)\rangle|\text{vac}\rangle+|g\rangle\sum_{l}\int^{\infty}_{-\infty}dt_{r}\,|\chi_{l}(t,t_{r})\rangle a^{\dagger}_{l}(t_{r})|\text{vac}\rangle,$$

(41)

with

	$\displaystyle|\gamma(t)\rangle$	$\displaystyle=-\int^{t}_{0}d\tau\,\xi(\tau)e^{(-iH_{\text{sys+vib}}-\frac{1}{2}\sum_{l}L^{\dagger}_{l}L_{l})(t-\tau)}L^{\dagger}_{\text{inc}}|g\rangle e^{-iH_{\text{vib}}\tau}|v\rangle$		(42a)
	$\displaystyle\chi_{l}(t,t_{r})$	$\displaystyle=\begin{cases}\delta_{l,\text{inc}}\xi(t_{r})e^{-iH_{\text{vib}}t}|v\rangle\quad,t<t_{r}\\ \delta_{l,\text{inc}}\xi(t_{r})e^{-iH_{\text{vib}}t}|v\rangle+\frac{1}{2}e^{-iH_{\text{vib}}(t-t_{r})}\langle g|L_{l}|\gamma(t_{r})\rangle\quad,t=t_{r}\\ \delta_{l,\text{inc}}\xi(t_{r})e^{-iH_{\text{vib}}t}|v\rangle+e^{-iH_{\text{vib}}(t-t_{r})}\langle g|L_{l}|\gamma(t_{r})\rangle\quad,t>t_{r}.\end{cases}$		(42b)

Here $|\gamma(t)\rangle$ is an unnormalized system+vibration state with the system being in the excited subspace, and $\chi_{l}(t,t_{r})$ is an unnormalized vibration state at time $t$. One can check that

$$\begin{split}&\rho_{1,1}(t)=|\gamma(t)\rangle\langle\gamma(t)|+|g\rangle\langle g|\sum_{l}\int^{\infty}_{-\infty}dt_{r}\,|\chi_{l}(t,t_{r})\rangle\langle\chi_{l}(t,t_{r})|\\ &\rho_{1,0}(t)=|\gamma(t)\rangle\langle g|\langle v|e^{iH_{\text{vib}}t}\\ &\rho_{0,0}(t)=e^{-iH_{\text{vib}}t}|v\rangle|g\rangle\langle g|\langle v|e^{iH_{\text{vib}}t}\end{split}$$

(43)

solves the single photon Fock state master equation (Eq. (28)) with vibration. For later convenience, we drop the spontaneous emission terms and define

$$|\gamma_{\xi,v}^{\prime}(t)\rangle\equiv-\int^{t}_{0}d\tau\,\xi(\tau)e^{-iH_{\text{sys+vib}}(t-\tau)}L^{\dagger}_{\text{inc}}|g\rangle e^{-iH_{\text{vib}}\tau}|v\rangle$$

(44)

to emphasize the dependence on the temporal profile $\xi(t)$ and the initial vibrational state $|v\rangle$.

Another way to treat the vibration effects is to use the HEOM formalism. We let each chromophore couple to a phonon bath with a Drude-Lorentz spectral density

$$J(\omega)=\frac{2\lambda\gamma\omega}{\omega^{2}+\gamma^{2}}.$$

(45)

The numerical values of the parameters in the spectral density are taken from [29]. Specifically, $\lambda$ is the reorganization energy, taken to be 37 $\text{cm}^{-1}$ for all sites, and $\gamma$, with physical dimension of frequency, is the decay rate of the phonon correlation function (see appendix H), which characterizes the time scale of vibration-induced fluctuations in the electronic excitation energy. For chlorophyll a, we take $\gamma=30\,\text{cm}^{-1}$; for chlorophyll b, we take $\gamma=48\,\text{cm}^{-1}$. Under high temperatures (characterized by $\hbar\gamma/k_{B}T\ll 1$, where $k_{B}$ and $T$ are the Boltzmann constant and temperature, respectively), the interaction picture system density matrix (i.e., with $\sum_{j}\omega_{0}|j\rangle\langle j|$ rotated out) follows the hierarchical equations of motion [13, 1]

$$\frac{d}{dt}\rho^{\vec{n}}(t)=-iH^{\times}\rho^{\vec{n}}-(\sum_{j}n_{j}\gamma_{j})\rho^{\vec{n}}-\sum_{j}\lambda_{j}P_{j}^{\times}\rho^{\vec{n}+\hat{e}_{j}}+n_{j}(2k_{B}TP_{j}^{\times}-i\gamma_{j}P_{j}^{o})\rho^{\vec{n}-\hat{e}_{j}},$$

(46)

where $P_{j}\equiv|j\rangle\langle j|$, $A^{o}B\equiv\{A,B\}$ is the anticommutator superoperator, and $A^{\times}B\equiv[A,B]$ is the commutator superoperator. $\vec{n}=(n_{1},\cdots,n_{N})$ is a vector of N integers, where the element $n_{j}$ is the hierarchical level of the j^th site. $\hat{e}_{j}\equiv(0,\cdots,0,1,0,\cdots,0)$ is the “unit vector” with the j^th element equals to 1 and all other elements equal to 0. $\rho^{(\vec{0})}$ is the physical density matrix, and the other $\rho^{(\vec{n})}$’s are non-physical auxiliary density matrices that capture the non-Markovian effects of the phonon environment. The initial contition is

$$\rho^{\vec{n}}(0)=\begin{cases}\rho_{\text{sys}}(0)\quad,\vec{n}=\vec{0}\\ 0\quad,\vec{n}\neq\vec{0}.\end{cases}$$

(47)

Numerically, a cutoff level $N_{\text{cutoff}}$ has to be introduced so that only a finite number of auxiliary density matrices with $\sum_{j}n_{j}\leq N_{\text{cutoff}}$ are solved. Given the cutoff, the total number of auxiliary density matrices is $\binom{N+N_{\text{cutoff}}}{N_{\text{cutoff}}}$ [30]. The auxiliary density matrices having $\sum_{j}n_{j}=N_{\text{cutoff}}$ are called the terminators. To capture the effects of one-higher-level set of auxiliary density matrices $\rho^{(\vec{n})}$, we have developed the following modified equations of motion for the terminators (appendix I):

$$\begin{split}\frac{d}{dt}\rho^{\vec{n}}(t)=&-iH^{\times}\rho^{\vec{n}}-(\sum_{j}n_{j}\gamma_{j})\rho^{\vec{n}}+\sum_{j}n_{j}(2k_{B}TP_{j}^{\times}-i\gamma_{j}P_{j}^{o})\rho^{\vec{n}-\hat{e}_{j}}\\ &-\sum_{j,k}\lambda_{j}\frac{n_{k}+\delta_{j,k}}{\gamma_{j}+\sum_{l}n_{l}\gamma_{l}}P_{j}^{\times}(2k_{B}TP_{k}^{\times}-i\gamma_{k}P_{k}^{o})\rho^{\vec{n}+\hat{e}_{j}-\hat{e}_{k}}.\end{split}$$

(48)

To simultaneously study the effects of the single photon and the phonon bath on the excitonic system, we now combine the input-output and HEOM formalisms. Two different formalisms are needed to treat the effects of these two bosonic baths because of their different properties. The input-output formalism is based on the frequency-independent coupling and the wide-band approximation (see Section 2.3), which has been used extensively in quantum optics to treat the interaction of matter with the photon field. The coupling to phonons, on the other hand, is frequency-dependent, and therefore cannot be treated with the assumptions of the input-output formalism. The input-output formalism allows us to explicitly calculate properties of the outgoing photon field (see Section 2.7), while the HEOM formalism traces out the bath degrees of freedom. The HEOM formalism is well-suited to treat the coupling to phonons, since the phonon correlation function is Gaussian (see appendix H), while the correlation function (e.g. $\langle a^{\dagger}(t_{2})a^{\dagger}(t_{1})\rangle$, $\langle a^{\dagger}(t_{2})a(t_{1})\rangle$, etc.) of an N-photon Fock state is not Gaussian. As an aside, we note that in contrast to a Fock state, for a multimode coherent state the correlation function is Gaussian, and for a coherent state input one can in fact treat the interaction with photons using the HEOM formalism. In addition, because the second cumulant $\langle a(t_{2})a^{\dagger}(t_{1})\rangle-\langle a(t_{2})\rangle\langle a^{\dagger}(t_{1})\rangle$ is proportional to a delta function for a coherent state, the resulting reduced system dynamics is Markovian and does not involve auxiliary density matrices (see Section 3.1) [10, 7].

To combine the input-output and HEOM formalisms, we use Eqs. (7) and (37) to write the full Hamiltonian as

$$H_{\text{total}}=H_{\text{field}}+H_{\text{vib}}+\sum_{j}\omega_{0}|j\rangle\langle j|+\underbrace{(H_{\text{sys}}-\sum_{j}\omega_{0}|j\rangle\langle j|)}_{H}+H_{\text{sys-field}}+H_{\text{sys-vib}}.$$

(49)

Here the term inside the parenthesis is the Hamiltonian appearing in the Fock state master equation, and will be denoted simply as $H$. Moving into the interaction picture where we now rotate out $H_{\text{field}}+H_{\text{vib}}+\sum_{j}\omega_{0}|j\rangle\langle j|$, the full Hamiltonian becomes

$$H_{\text{total}}(t)=H+\sum_{j}|j\rangle\langle j|u_{j}(t)+\sum_{l}(-ia_{l}(t)L_{l}^{\dagger}+\text{h.c.}).$$

(50)

This is to be compared with the interaction picture Hamiltonian for the system+field state alone, Eq. (18). Given a multimode Fock state photon in one spatial mode as the input field, the reduced dynamics in the system+vibration degrees of freedom is then given by the Fock state master equation Eq. (28), with $H$ replaced by $H+\sum_{j}|j\rangle\langle j|u_{j}(t)$. To apply the HEOM formalism to this, we rewrite the Fock state master equation in the block matrix form

$$\begin{split}&\frac{d}{dt}\begin{pmatrix}\rho_{N,N}\\ \vdots\\ \rho_{m,n}\\ \vdots\\ \rho_{0,0}\end{pmatrix}=\\ &\begin{pmatrix}-i\sum_{j}(P_{j}u_{j}(t))^{\times}\rho_{N,N}+\big{(}-iH^{\times}+\sum_{l}\mathcal{D}[L_{l}]\big{)}\rho_{N,N}-\sqrt{N}\xi(t)L^{\dagger\times}_{\text{inc}}\rho_{N-1,N}+\sqrt{N}\xi^{*}(t)L^{\times}_{\text{inc}}\rho_{N,N-1}\\ \vdots\\ -i\sum_{j}(P_{j}u_{j}(t))^{\times}\rho_{m,n}+\big{(}-iH^{\times}+\sum_{l}\mathcal{D}[L_{l}]\big{)}\rho_{m,n}-\sqrt{m}\xi(t)L^{\dagger\times}_{\text{inc}}\rho_{m-1,n}+\sqrt{n}\xi^{*}(t)L^{\times}_{\text{inc}}\rho_{m,n-1}\\ \vdots\\ -i\sum_{j}(P_{j}u_{j}(t))^{\times}\rho_{0,0}+\big{(}-iH^{\times}+\sum_{l}\mathcal{D}[L_{l}]\big{)}\rho_{0,0}\end{pmatrix}.\end{split}$$

(51)

Eq. (51) can be written in a more compact notation as

$$\frac{d}{dt}\Xi(t)=(\mathcal{V}(t)+\mathcal{W}(t))\Xi(t),$$

(52)

where

$$\Xi(t)=\begin{pmatrix}\rho_{N,N}\\ \vdots\\ \rho_{0,0}\end{pmatrix},$$

(53)

and $\mathcal{V}(t)$ and $\mathcal{W}(t)$ are linear operators on $\Xi$. $\mathcal{V}(t)$ is the operator that acts nontrivially on the vibrational degrees of freedom. Its effect on $\Xi$ is given by

$$\mathcal{V}(t)\Xi=\begin{pmatrix}-i\sum_{j}(P_{j}u_{j}(t))^{\times}\rho_{N,N}\\ \vdots\\ -i\sum_{j}(P_{j}u_{j}(t))^{\times}\rho_{0,0}\end{pmatrix}.$$

(54)

The effect of $\mathcal{W}(t)$ on $\Xi$ is to produce the rest of the terms in Eq. (51). We note that $\mathcal{W}(t)$ acts trivially on the vibrational degrees of freedom. Now, from Eq. (52), we can write the vector $\chi(t)$ of reduced Fock state auxiliary density matrices on the system, i.e.,

$$\chi(t)=\text{Tr}_{\text{vib}}(\Xi(t))=\begin{pmatrix}\text{Tr}_{\text{vib}}(\rho_{N,N})\\ \vdots\\ \text{Tr}_{\text{vib}}(\rho_{0,0})\end{pmatrix}$$

(55)

formally as a time-ordered exponential

$$\chi(t)=\text{Tr}_{\text{vib}}\Big{(}\hat{T}\exp\big{(}\int^{t}_{0}d\tau\,(\mathcal{V}(\tau)+\mathcal{W}(\tau)\big{)}\rho_{\text{vib,thermal}}\Big{)}\chi(0),$$

(56)

where we have used the fact $\Xi(0)=\chi(0)\otimes\rho_{\text{vib,thermal}}$ to pull out $\chi(0)$ from the partial trace. Using the Gaussian property of the vibration correlation function (see appendix H), we perform a generalized cumulant expansion [31] on the time-ordered exponential and obtain

$$\chi(t)=\hat{T}\mathcal{Z}\chi(0),$$

(57)

where $\mathcal{Z}$ is defined as

$$\mathcal{Z}=\exp(\int^{t}_{0}dt_{1}\,\mathcal{W}(t_{1})-\sum_{j}\int^{t}_{0}dt_{2}\int^{t_{2}}_{0}dt_{1}\,\lambda_{j}e^{-\gamma_{j}(t_{2}-t_{1})}P_{j}^{\times}(t_{2})(2k_{B}TP_{j}^{\times}(t_{1})-i\gamma_{j}P_{j}^{o}(t_{1}))).$$

(58)

The integrand of the double integral is now understood as an operator that applies to every block matrix $\text{Tr}_{\text{vib}}(\rho_{m,n})$ of $\chi$. Following a similar procedure as in appendix H, we then obtain the Fock state + HEOM master equation

$$\frac{d}{dt}\chi^{\vec{n}}(t)=\mathcal{W}(t)\chi^{\vec{n}}-(\sum_{j}n_{j}\gamma_{j})\chi^{\vec{n}}-\sum_{j}\lambda_{j}P_{j}^{\times}\chi^{\vec{n}+\hat{e}_{j}}+n_{j}(2k_{B}TP_{j}^{\times}-i\gamma_{j}P_{j}^{o})\chi^{\vec{n}-\hat{e}_{j}}.$$

(59)

Written in terms of individual auxiliary density matrices, this is equivalent to

$$\begin{split}\frac{d}{dt}\rho^{\vec{n}}_{m,n}=&(-iH^{\times}+\sum_{l}\mathcal{D}[L_{l}])\rho^{\vec{n}}_{m,n}-\sqrt{m}\xi(t)L^{\dagger\times}_{\text{inc}}\rho^{\vec{n}}_{m-1,n}+\sqrt{n}\xi^{*}(t)L^{\times}_{\text{inc}}\rho^{\vec{n}}_{m,n-1}\\ &-(\sum_{j}n_{j}\gamma_{j})\rho^{\vec{n}}_{m,n}-\sum_{j}\lambda_{j}P_{j}^{\times}\rho^{\vec{n}+\hat{e}_{j}}_{m,n}+n_{j}(2k_{B}TP_{j}^{\times}-i\gamma_{j}P_{j}^{o})\rho^{\vec{n}-\hat{e}_{j}}_{m,n},\end{split}$$

(60)

with the initial condition

$$\rho^{\vec{n}}_{m,n}(0)=\begin{cases}\rho_{\text{sys}}(0)\quad,\vec{n}=\vec{0}\text{ and }m=n\\ 0\quad,\text{otherwise}.\end{cases}$$

(61)

Note that $\rho^{\vec{0}}_{N,N}$ is the only physical density matrix.