Reconstruction of the dynamics (quantum process tomography) of the single-exciton manifold in energy transfer systems is proposed here on the basis of two-dimensional fluorescence spectroscopy (2D-FS) with phase-modulation. The quantum-process-tomography protocol introduced here benefits from, e.g., the sensitivity enhancement ascribed to 2D-FS. Although the isotropically averaged spectroscopic signals depend on the quantum yield parameter Γ of the doubly excited-exciton manifold, it is shown that the reconstruction of the dynamics is insensitive to this parameter. Applications to foundational and applied problems, as well as further extensions, are discussed.

## I. INTRODUCTION

Since Chuang and Nielsen’s 1996 seminal proposal to experimentally reconstruct the evolution operator of a quantum black box,^{1} a variety of experiments have been proposed^{2–4} and some of these have been implemented.^{5–7} Although most of these ideas involved relatively “clean” optical systems in the context of quantum information processing (QIP), recently there has been growing interest in the application of quantum process tomography to study electronically coupled molecular systems.^{8–10} Such combined theoretical and experimental studies on excitation energy transfer serves to bring together the QIP and physical chemistry communities.

The current version of Quantum Process Tomography (QPT) for excitonic systems^{8,9,11} is based on the method of 2D Photon Echo Spectroscopy (2D-PES).^{12,13} It therefore relies on the wave-vector phase-matching condition, which works for macroscopic systems with many chromophores (see Refs. 11–13 for details). For this reason, the proposal developed in Refs. 8–10 is not suitable for single-molecule QPT, and it is hence desirable to implement the phase-cycling technique (PCT) or phase-modulation techniques (PMTs).^{14,15} In most 2D-PES experiments, the system interacts with three non-collinear ultrafast laser pulses, which gives rise to a third-order polarization that propagates in the wave-vector matched direction.^{11–13} The transmitted signal must be separated from background laser light, which is inadvertently scattered by the sample in the same direction as the third-order signal. The presence of background scattering is a limiting factor to the sensitivity of 2D-PES experiments. In two-dimensional fluorescence spectroscopy (2D-FS), the system interacts with four collinear laser pulses, and the signal is detected by monitoring nonlinear contributions to the ensuing fluorescence signal.^{15,16} The red-shifted fluorescence can be easily separated from background scattered laser light by using long-pass spectral filters. Specific nonlinear contributions to the fluorescence signal are isolated according to the phase-modulation schemes described in the previous work.^{15,16} Because the 2D-FS method is based on the detection of incoherent fluorescence signals, it may be applied to systems of small numbers of chromophores, quantum dots, or thin film materials.

Despite the high-sensitivity advantages afforded to fluorescence-detected PMTs, which are useful for studies of biological molecules and molecular aggregates,^{17–19} these methods have been less commonly practiced than four-wave mixing approaches to 2D-PES. However, recent theoretical^{15} and experimental^{16,20} progress with classical light have enabled PMTs for a variety of complex molecular systems relevant to exciton dynamics. These recent developments and the general theory of open quantum systems—i.e., quantum systems coupled to the environment—are combined here to formulate a self-consistent theory of QPT that is based on collinear PMT with synchronous detection.

For the ideal situation, when nonradiative processes are neglected in the doubly excited-exciton manifold, the quantum yield parameter of this manifold is set to Γ = 2. Under this circumstance, it was shown that 2D-FS coincides with 2D-PES.^{15,16,20} It is shown below that this equivalence also holds at the level of quantum process tomography, i.e., the protocol introduced here generalizes the protocol in Refs. 8 and 9 to the more realistic situation 0 ≤ Γ < 2.

## II. INITIAL CONSIDERATIONS ON QPT, SYSTEM MODEL, AND 2D-FS

Before introducing the reconstruction of the dynamics, it is necessary to state some remarks on the basics of process tomography, the system-of-interest model, and 2D-FS.

*Quantum Process Tomography Tensor*—In quantum mechanics, the state of a physicochemical system S is described by a density operator $ \rho \u02c6 $. Time evolution of quantum states is governed by the Schrödinger equation, which is linear in the state of the system. This linearity allows for a description of the system’s dynamics in terms of a linear map, $ \chi \u02c6 t : \rho \u02c6 0 \u21a6 \rho \u02c6 t $. After projecting onto a complete orthonormal basis {|*n*〉}, the map reads as

where *χ*_{nmνμ}(*t*) stands for the process tomography tensor.^{2,8,21} For Hamiltonian dynamics with $ H \u02c6 |n\u3009= E n |n\u3009$, *χ*_{nm,νμ}(*t*) = e^{−i(Em−En) t/ħ}*δ*_{nν}*δ*_{mμ}.^{22–24} Thus, population-to-coherence [*χ*_{nmνν}(*t*)] and the reverse [*χ*_{nnνμ}(*t*)] process are prevented by the Kronecker deltas *δ*_{nν}*δ*_{mμ}. Clearly, this restriction is not present if driving fields are present or if the system of interest is coupled to its environment.^{22–24}

In the general case of open quantum systems, the functional form of Eq. (1) remains valid under some conditions: (i) If the coupling to the bath is weak, Eq. (1) holds for Markovian and non-Markovian processes and the process tensor is independent of the initial state (see, e.g., Refs. 8 and 9 and references therein). (ii) If the coupling to the bath is strong, and the initial system-environment correlations cannot be neglected, Eq. (1) holds after including those initial correlations in *χ*_{nmνμ}(*t*) (see Refs. 22–24 for details). (iii) Because the initial bath correlations vanish at high temperature, even for strong coupling,^{25} *χ*_{nmνμ}(*t*) can be defined independently of the initial state in the strong coupling regime entered at high temperatures.

After identifying the conditions under which Eq. (1) holds, it is relevant to consider some of the main properties of the QPT tensor,^{8} namely,

where *z* is any complex valued vector. Equation (2) ensures the Hermitian character of the density operator, $ \rho \u02c6 = \rho \u02c6 \u2020 $, while Eq. (3) guaranties probability conservation, $tr\u2009 \rho \u02c6 ( t ) =1$. The last property is a consequence of the fact that $ \rho \u02c6 ( t ) $ remains positive-semidefinite under unitary operations.

The objective of QPT is the experimental reconstruction of the process tomography tensor *χ*_{nmνμ}(*t*).

*Model*—Consider an excitonic dimer described by $ H \u02c6 S $ and given by

where $ a \u02c6 i \u2020 $ and $ a \u02c6 i $ are the creation and annihilation operators for site *i*, and *ϖ*_{1} ≠ *ϖ*_{2} are the site energies while *J* ≠ 0 is the Coulombic coupling between chromophores. By defining the average frequency $\varpi = 1 2 ( \varpi 1 + \varpi 2 ) $, the half-difference $\Delta = 1 2 ( \varpi 1 \u2212 \varpi 2 ) $, and the mixing angle $\theta = 1 2 arctan ( J / \Delta ) $, it is possible to introduce the creation and annihilation operators, $ c \u02c6 p =cos\theta a \u02c6 1 +sin\theta a \u02c6 2 $ and $ c \u02c6 p \u2020 =sin\theta a \u02c6 1 \u2020 +cos\theta a \u02c6 2 \u2020 $, of the *p*th delocalized exciton state with energy *ϖ _{p}* =

*ϖ*± Δsec2

*θ*and

*p*∈ {

*e*,

*e*′}. Starting from the ground state |

*g*〉, the single-exciton states are conveniently defined as $|\u2009e\u2009\u3009= c \u02c6 e \u2020 |\u2009g\u2009\u3009$ and $|\u2009 e \u2032 \u2009\u3009= c \u02c6 e \u2032 \u2020 |\u2009g\u2009\u3009$, while the biexciton state as $|\u2009f\u2009\u3009= a \u02c6 1 \u2020 a \u02c6 2 \u2020 |\u2009g\u2009\u3009= c \u02c6 e \u2020 c \u02c6 e \u2032 \u2020 |\u2009g\u2009\u3009$ with

*ϖ*=

_{f}*ϖ*

_{1}+

*ϖ*

_{2}=

*ϖ*+

_{e}*ϖ*

_{e′}. The dipole vectors at each site are set to

**d**

_{1}=

*d*

_{1}

**e**

_{z}and

**d**

_{2}=

*d*

_{2}cos(

*ϕ*)

**e**

_{z}+

*d*

_{2}sin(

*ϕ*)

**e**

_{x}. So that,

*μ*_{eg}=

*d*

_{2}sin

*θ*sin

*ϕ*

**e**

_{x}+ (

*d*

_{1}cos

*θ*+

*d*

_{2}sin

*θ*cos

*ϕ*)

**e**

_{z},

*μ*_{e′g}=

*d*

_{2}cos

*θ*sin

*ϕ*

**e**

_{x}+ (−

*d*

_{1}sin

*θ*+

*d*

_{2}cos

*θ*cos

*ϕ*)

**e**

_{z},

*μ*_{fe}=

*d*

_{2}cos

*θ*sin

*ϕ*

**e**

_{x}+ (

*d*

_{1}sin

*θ*+

*d*

_{2}cos

*θ*cos

*ϕ*)

**e**

_{z}, and $ \mu f e \u2032 =\u2212 d 2 sin\theta sin\varphi \u2009 e x + d 1 cos \theta \u2212 d 2 sin \theta cos \varphi e z $.

Although exciton-exciton binding or repulsion terms are not included here, it is considered that each excitonic manifold contributes to the spectroscopic signal with a weight given by their fluorescence quantum yield coefficients Γ_{ν}. Specifically, it is assumed that the quantum yield of the two singly excitonic states is the same and equal to 1, while for the doubly excitonic manifold, it is assumed that Γ_{f} = Γ with 0 ≤ Γ ≤ 2. In the ideal case in which two photons are emitted via the path | *f* 〉 → | *e*, *e*′ 〉 → | *g* 〉, Γ = 2. It is possible that singlet-singlet annihilation would convert a doubly excited excition into a singly excited exciton,^{26} which in the absence of non-radiative decay would result in Γ = 1. However, because of the abundance of non-radiative relaxation pathways for highly excited states, the quantum yield of the doubly excitonic manifold is expected to be smaller than that of the singly excitonic manifold, so that values smaller than unity are expected. For example, for membrane-supported self-assembled porphyrin dimers, it was found that Γ = 0.31.^{16,20}

For convenience, the dimer Hamiltonian can be written as $ H \u02c6 S = \u2211 \nu = { g , e , e \u2032 , f} \omega \nu \nu \nu $. To account for the influence of the local vibrational environment in the excitonic dimer, coupling to a thermally equilibrated phonon bath at inverse temperature *β* is considered next. Specifically, the Hamiltonian of the environment is given by $ H \u02c6 E = \u2211 p = e , e \u2032 \u2211 n \Omega n , p ( b \u02c6 n , p \u2020 b \u02c6 n , p + 1 / 2 ) $, where Ω_{n,p} denotes the frequency of the environment modes. The interaction is described by $ H \u02c6 SE = E \u02c6 e e e + E \u02c6 e \u2032 e \u2032 e \u2032 + E \u02c6 e + E \u02c6 e \u2032 e \u2032 e \u2032 $ with $ E \u02c6 p = \u2211 n \lambda n , p ( b \u02c6 n , p \u2020 + b \u02c6 n , p ) $. $ b \u02c6 n , p \u2020 $ and $ b \u02c6 n , p $ are the creation and annihilation bosonic operators of the *n*th mode of the vibrational environment in the *p*-site. *λ*_{n,p} measures the interaction strength between the *n*th mode of the environment and the *p*th site. The net effect of the local environment is encoded in the spectral density $ J n = \u2211 n \Omega n , p 2 \lambda n , p 2 \delta ( \omega \u2212 \omega n ) $.

*2D-FS*—The main difference between the QPT scheme introduced below and the previous QPT proposals is the spectroscopic technique, 2D-FS, which the present proposal is based on. It is therefore relevant to discuss the main differences and advantages that 2D-FS has over, e.g., 2D-PES. The 2D-FS method isolates the nonlinear optical response of a material system by monitoring fluorescence signals. Because fluorescence can be efficiently separated from background scattered light, the 2D-FS approach can be used to perform experiments that require very high detection sensitivity.^{15,16,20} Moreover, the collinear beam geometry used in 2D-FS has the advantage that every illuminated molecule experiences the same optical phase condition at every instant in time. Since the incoherent fluorescence signal is emitted isotropically, very small numbers of molecules may be studied in this way.^{15,16,20}

The 2D-FS observable is proportional to the fourth-order excited populations,

with $ A \u02c6 = \u2211 \nu = { e , e \u2032 , f} \Gamma \nu \nu \nu $, generated by the action of the operator $ V \u02c6 ( t \u2032 ) $ that comprises the excitation by four weak non-overlapping laser pulses,

Here, *λ* denotes the maximum intensity of the pulses’ electric field, and $ \mu \u02c6 $ denotes the dipole operator. **e**_{i}, *t _{i}*,

*ω*, and

_{i}*ϕ*stand for the polarization vector, time center, frequency, and phase of the

_{i}*i*th laser pulse, respectively. The pulse envelope

*E*(

*t*) is chosen to be Gaussian with fixed width

*σ*, i.e.,

*E*(

*t*) = e

^{−t2/2σ2}. In the model under consideration, the only optically allowed transitions are between states differing by one excitation. Hence, the only non-vanishing dipole transition matrix elements are

*μ*_{ij}=

*μ*_{ji}with

*ij*= {

*eg*,

*e*′

*g*,

*fe*,

*fe*′}. Details about the derivation and the explicit functional form of the fourth-order density matrix can be found in Appendices A and B.

For the purpose of extracting the QPT tensor from the 2D-FS experimental signals, only the rephasing signals with global phase *ϕ*_{reph} = − *ϕ*_{1} + *ϕ*_{2} + *ϕ*_{3} − *ϕ*_{4} will be considered below (see Fig. 1). Thus, assuming that the rotating wave approximation (RWA) holds, the interactions with the electromagnetic fields are characterized by

where $ \mu \u02c6 < = \u2211 \omega p < \omega q \mu p q |p\u3009\u3008q|$ promotes emissions from the ket and absorptions on the bra, and $ \mu \u02c6 > = ( \mu \u02c6 < ) \u2020 $ induces the opposite processes. For this particular selection of the global phase, *ϕ*_{reph} = − *ϕ*_{1} + *ϕ*_{2} + *ϕ*_{3} − *ϕ*_{4}, the 2D-FS signals are equivalent to the rephasing spectroscopic signals in the photon-echo direction **k**_{PE} = − **k**_{1} + **k**_{2} + **k**_{3} when Γ = 2.^{15,16,20} It is shown below that the present QPT protocol reduces to the protocol in Refs. 8 and 9 when Γ = 2 as well.

## III. 2D-FS QPT

As stated above, the main goal of QPT is the reconstruction of the dynamics of the density operator. In doing so, it is assumed that the structural parameters of the model, namely, the transition frequencies *ϖ _{ij}* =

*ϖ*−

_{i}*ϖ*and the electric dipole transition matrix element

_{j}

*μ*_{ij}are all known. This prerequisite is not an issue because information about the transition frequencies is routinely obtained from linear absorption spectra, and the transition dipole directions can be inferred from structural measurements and polarization spectroscopy.

^{8}

Once the structural parameters are defined, the reconstruction of the dynamics comprises three main parts: (i) initial state preparation, (ii) evolution, and (iii) final state detection. In describing these stages, it is useful to introduce the standard time intervals {*τ*, *T*, *t*} instead of the time center *t _{i}* of each pulse.

^{12,13}The time difference between the second and the first pulse defines the

*coherence time*interval

*τ*=

*t*

_{2}−

*t*

_{1}. The time interval between the third and the second pulse,

*T*=

*t*

_{3}−

*t*

_{2}, is known as the

*waiting time*, which defines the quantum channel to be characterized by the QPT scheme. Finally, the difference between the fourth and the third pulse,

*t*=

*t*

_{4}−

*t*

_{3}, denotes the

*echo time*.

*Initial State Preparation*—The excitonic system, before any electromagnetic perturbation, is assumed to be in the ground state, $ \rho \u02c6 ( \u2212 \u221e ) = g g $. Thus, the basic idea is to make use of the first two pulses to prepare the effective initial density matrix at *T* = 0, $ \rho \u02c6 e 1 , e 2 \omega 1 , \omega 2 ( T = 0 ) $, and use the last two pulses to read out the state.

After applying second order perturbation theory in *λ*, and under the assumption that the RWA holds in this case (see Appendix A for details), the effective initial state reads as

where $ G i j ( \tau ) $ is the propagator of the optical coherence $ i j $. For simplicity, it can be assumed as $ G i j ( \tau ) =\Theta ( \tau ) exp [ ( \u2212 i \varpi i j \u2212 \Gamma i j ) \tau ] $ begin Γ_{ij} dephasing rates, and the Heaviside function Θ(*τ*) ensures causality. The coefficients $ C \omega i p $ are purely imaginary and given by $ C \omega i p =i\lambda 2 \pi \sigma 2 e \u2212 \sigma 2 ( \varpi p g \u2212 \omega i ) $. Because at this level there is no influence of the doubly excited exciton manifold, the effective initial state in Eq. (12) coincides with the effective initial state prepared by 2D-PES in Ref. 8.

The contributions to $ \rho \u02c6 e 1 , e 2 \omega 1 , \omega 2 ( 0 ) $ in Eq. (12) are clearly depicted in the lower panel of Fig. 1 (see Figs. A–D). Starting from the system the ground state $ g g $, in the RWA and selecting only those contributions with −*ϕ*_{1} + *ϕ*_{2}, the first pulse can only excite the bra and then creates an optical coherence $ g p $ with probability $ C \omega 1 p $. This coherence $ g p $ evolves under the action of $ G g p ( \tau ) $ for a time *τ* when the second pulse prepares the state $ q p $ with probability $ C \omega 2 q $ or a hole $\u2212 g g $ with probability $ C \omega 2 p $.

As in the case of QPT based on 2D-PES,^{8,9} to prepare the set of four linearly independent states in Eq. (12) [see also Fig. 1(A–D)], it suffices to consider a pulse toolbox of two waveforms with carrier frequencies {*ω*_{+}, *ω*_{−}} that create | *e* 〉 and | *e*′ 〉 with different amplitudes. Of course, the discrimination in the preparation of | *e* 〉 and | *e*′ 〉 depends on how close to resonance the carrier frequencies are. For an extensive and detailed analysis on this respect, see Refs. 8 and 9.

*Evolution*—Once the initial state $ \rho \u02c6 e 1 , e 2 \omega 1 , \omega 2 ( 0 ) $ is effectively prepared, i.e., after the action of the first-two pulses *T* ≳ 3*σ*, the system evolves over a time *T* under the action of the superoperator $ \chi \u02c6 ( T ) $, according to

To avoid contamination of the initial state by terms proportional to a hole every time there is a single-exciton population $ p p $, it is assumed that $\u3008ab| \chi \u02c6 ( T ) |gg\u3009= \chi a b g g ( T ) = \delta a g \delta b g $, which is equivalent to neglect processes where phonons induce upward optical transitions and spontaneous excitation from the single to the double exciton manifolds.^{8} Up to this condition, $ \chi \u02c6 ( T ) $ in Eq. (13) describes the dynamics induced by any bath model and accounts for any system-bath coupling.

*Final State Detection*—The very nature of the fluorescence detection in 2D-FS suggests considering contributions from excitation configurations that lead to populations only. In the upper panel of Fig. 1, those contributions are schematically displayed and grouped according to the probability $ C \omega 3 p C \omega 4 q $ that they occur.

In contrast to QPT, which is based on 2D-PES, there are here 14 possibilities for the final state instead of 10. Thus, 20 independent experiments are needed instead of 16. This comes at the expense of the different role that the fourth pulse has in each technique, namely, heterodyne detection in 2D-PES and the generation of populations in 2D-FS. However, the signals that lead to population of the doubly excited exciton manifold $ f \u3009 \u3008 f $ from the coherence $ f \u3009 \u3008 e $ must be summed up to the signal that lead to the population $ e \u3009 \u3008 e $. In the summation, the process $ f \u3009 \u3008 e \u2192 f \u3009 \u3008 f $ is weighted with a factor Γ, while the process $ f \u3009 \u3008 e \u2192 e \u3009 \u3008 e $ has weight −1. The same procedure applies to the process $ f \u3009 \u3008 e \u2032 \u2192 f \u3009 \u3008 f $ and $ f \u3009 \u3008 e \u2032 \u2192 e \u2032 \u3009 \u3008 e \u2032 $. This procedure leads to the 16 independent signals that are needed to reconstruct the 16 elements of the process tensor *χ*_{nmνμ}(*T*).

By using the same toolbox as in the preparation state and following closely the notation in Ref. 8, the total signal is given by

with

and

Analogous expressions hold for $ P e 1 , e 2 , e 3 , e 4 p , q , e \u2032 , e \u2032 $ and $ P e 1 , e 2 , e 3 , e 4 p , q , e \u2032 , e $ after interchanging *e*↔*e*′. The 2D-FS signals $ [ S FS ] e 1 , e 2 , e 3 , e 4 \omega 1 , \omega 2 , \omega 3 , \omega 4 ( \tau , T , t ) $ in Eq. (14) with (15) and (16) are the main result of this article. They allow for the reconstruction of the dynamics of excitonic systems based on 2D-FS that is an attractive approach to reach QPT at the level of single molecules. Remarkably, the appealing form of Eqs. (14)–(16) allows for an immediate connection with the protocol derived in Refs. 8 and 9. Specifically, up to a global minus sign that is consistent with previous investigations,^{16,20} results in Refs. 8 and 9 are obtained by simply setting Γ = 2 in Eqs. (14)–(16).

Because the probed sample is an ensemble of isotropically distributed molecules in solution, an isotropic average of $ \mu a \u22c5 e 1 \mu b \u22c5 e 2 \mu c \u22c5 e 3 \mu d \u22c5 e 4 $ is needed. In doing so, standard procedures are followed (see, e.g., Chap. 11 in Ref. 27 or Sec. 3.3 in Ref. 13). Specifically, the isotropic average is given by

where **e**_{i} and **m**_{i} denote the polarization of the pulses in the laboratory and molecule-fixed frames, respectively. *e _{i}* = {

*e*,

_{xi}*e*,

_{yi}*e*} and

_{xi}*m*= {

_{i}*m*,

_{xi}*m*,

_{yi}*m*} are the components of the polarization vectors

_{xi}**e**

_{i}and

**m**

_{i}, respectively. The isotropically invariant tensor I

^{(4)}is given by

The explicit expression for the relevant case of interest in the collinear configuration used in 2D-FS, **e _{1}** =

**e**=

_{2}**e**=

_{3}**e**=

_{4}**z**, can be found in Appendix C. Thus, after isotropically averaging,

Because in 2D-FS the laser pulses are collinear, it is possible to set **e**_{j} = **z** at this point. Additionally, it is assumed below that *τ* = 0 and *t* = 0 so that only the signals $ P e 1 , e 2 , e 3 , e 4 p , q , r , s ( 0 , T , 0 ) $ are considered. However, following a similar procedure as in Ref. 8, this restriction can be relaxed and arbitrary *τ* and *t* can be considered. For the sake of generality, the polarization vectors are denoted independently by **e**_{j} and *τ*, and *t* are not set to zero in the above expressions.

The extraction procedure of the matrix elements of $ \chi \u02c6 $ from $\u3008 [ S FS ] e 1 , e 2 , e 3 , e 4 \omega 1 , \omega 2 , \omega 3 , \omega 4 ( \tau , T , t ) \u3009 iso $ follows from Eqs. (19), (15), and (16). Note that in doing so, the 16 2D-FS signals $\u3008 [ S FS ] e 1 , e 2 , e 3 , e 4 \omega 1 , \omega 2 , \omega 3 , \omega 4 ( \tau , T , t ) \u3009 iso $ and the 16 auxiliary signals $\u3008 P e 1 , e 2 , e 3 , e 4 p , q , r , s ( \tau , T , t ) \u3009 iso $ can be grouped into the 16-dimensional vectors 〈[S_{FS}](*τ*, *T*, *t*)〉_{iso} and 〈P(*τ*, *T*, *t*)〉_{iso}, respectively. This allows writing Eq. (19) as 〈[S_{FS}](*τ*, *T*, *t*)〉_{iso} = C〈P(*τ*, *T*, *t*)〉_{iso}, where the matrix elements of C contains the probabilities $ C \omega 1 p C \omega 2 q C \omega 3 r C \omega 4 s $. Then, the first step in the extraction procedure is to invert the matrix C so that the signals $\u3008 P e 1 , e 2 , e 3 , e 4 p , q , r , s \u3009 iso $ can be extracted from the measured signals $\u3008 [ S FS ] e 1 , e 2 , e 3 , e 4 \omega 1 , \omega 2 , \omega 3 , \omega 4 ( \tau , T , t ) \u3009 iso $, i.e., 〈P(*τ*, *T*, *t*)〉_{iso} = C^{−1}〈[S_{FS}](*τ*, *T*, *t*)〉_{iso}. The second step comprises the extraction of the 16 elements of the process tensor $ \chi \u02c6 ( T ) $ from the isotropically averaged version of Eqs. (15) and (16). This process can be accomplished by conveniently defining a 16-dimensional vector **χ**(*T*), such that **χ**(*T*) = M^{−1}〈P(0, *T*, 0)〉_{iso}. See the Appendix for further details.

## IV. NUMERICAL EXAMPLE

As a concrete example, consider parameters of relevance in the context of light-harvesting systems.^{28,29} Specifically, to compare with previous results,^{8} consider *ϖ*_{1} = 12 881 cm^{−1}, *ϖ*_{2} = 12 719 cm^{−1}, *J* = 120 cm^{−1}, *d*_{2}/*d*_{1} = 2, and *ϕ* = 0.3. The two-waveform toolbox is assumed to have frequencies *ω*_{+} = 13 480 cm^{−1} and *ω*_{−} = 12 130 cm^{−1} so that *ϖ _{i}* = {

*ω*

_{+},

*ω*

_{−}}, ∀

*i*and pulse width

*σ*= 40 fs. To simulate the signals, the spectral density

*J*(

_{n}*ω*) of the local vibrational environments is assumed identical and given by

*J*(

_{n}*ω*) = (

*λ*

_{r.e.}/

*ω*

_{c})

*ω*exp(−

*ω*/

*ω*

_{c}) where the cut-off frequency is set as

*ω*

_{c}= 120 cm

^{−1}while the reorganization energy is chosen as

*λ*

_{r.e.}= 30 cm

^{−1}. These set of parameters are relevant for light-harvesting systems and were used in Ref. 8.

In the simulations below, an inhomogeneously broadened ensemble of 10^{4} dimers with diagonal disorder is considered. Specifically, it is assumed that the site energies $ \varpi 1 \u2032 $ and $ \varpi 2 \u2032 $ in the ensemble follow a Gaussian distribution centered at *ϖ*_{1} and *ϖ*_{2} with standard deviation *σ*_{inh} = 40 cm^{−1}. The dynamics are solved at the level of the secular Redfield master equation at room temperature and for *T* ≥ 3*σ*.

Figure 2 depicts the nonvanishing real parts of $\u3008 P z , z , z , z p , q , r , s \u3009 iso $ for a variety of values of the quantum yields parameter Γ. From the functional dependence on Γ of the signals $\u3008 P z , z , z , z p , q , r , s \u3009 iso $ in Eqs. (15) and (16), three cases are of interest: (i) For Γ = 0, the contribution from the excited state absorption (ESA) pathways has the same sign as the stimulated emission (SE) and ground-state bleach (GSB) contributions (see the double-sided Feynman’s diagrams in Fig. 1 or the discussion in Ref. 16). Thus, the amplitude of the signals is the largest possible. (ii) For Γ = 1, the ESA pathways do not contribute to the signal and the amplitude of the signals is expected to be smaller than in the case of Γ = 0. (iii) For Γ = 2, the contribution from the ESA pathways has opposite sign to the SE and GSB contributions, so that the amplitude is expected to be smaller than in the previous case Γ = 1. These expectations are confirmed by simulations in Fig. 2. For completeness, the intermediate cases Γ = 0.5 and Γ = 1.5 were also depicted in Fig. 2.

Based on the signals $\u3008 P z , z , z , z p , q , r , s \u3009 iso $ obtained above, the QPT tensor is reconstructed in Fig. 3. The reconstruction appears to be insensitive to the value of the quantum yield parameter Γ. This unexpected result can be understood after noticing that the QPT tensor *χ*_{nmνμ}(*T*) is a characteristic of the singly exited exciton manifold and Γ is a function of the doubly excited exciton manifold. Thus, QPT of the singly exited exciton manifold by 2D-FS is robust against nonradiative processes of the doubly excited exciton manifold and benefits from the quality of the signals discussed above.

Results in Fig. 3 agree with the secular Redfield tensor used to simulate the 2D-FS signals. If the signal were obtained from experimental data, a careful analysis of the propagation of errors would be in order. In particular, it is necessary to include fluctuations in the laser intensity at each time *T* at which the signals are collected, and to pay attention to the stability conditions imposed by the invertibility of the matrices C and M.^{8} In this article, interest was in providing a proof-of-principle for the scheme derived above, so that Fig. 3 is aimed to depict the type of information that can be extracted from the protocol.

Specifically, (i) if for a particular photochemical system, non-negligible, non-secular terms emerge during the reconstruction of the process tensor $ \chi \u02c6 $ from experimental data that would imply, e.g., that coherent control schemes assisted by the environment^{23,30} may be applied in that particular system. (ii) If the decay of the tensor elements associated to the coherences of the density matrix, *χ _{nmnm}* with

*n*≠

*m*, is non-exponential, it may indicate the presence of non-Markovian dynamics.

^{31}The deviation from exponential decay behavior may even be considered as a measure of the non-Markovian character of the dynamics—a relevant topic in the context of open quantum systems. (iii) Although in multilevel systems the decay rate of the elements

*χ*cannot be directly associated to the decay rate of the coherences $\u3008n| \rho \u02c6 |m\u3009$ of the density matrix, the decay rate of

_{nmnm}*χ*provides information about the lifetime of particular transfer and coherent mechanisms.

_{nmnm}## V. DISCUSSION

Having experimental access to the process tomography tensor *χ*_{nmνμ}(*t*) is fundamental to revealing energy pathways in exciton dynamics, and in designing control strategies to increase transport efficiency. Specifically, applications of QPT to photosynthetic light-harvesting systems can (i) rule out certain transfer mechanisms proposed in the literature and (ii) address the question about the quantum/classical nature of the energy transport in certain biological systems from an experimental viewpoint. In addressing these issues, a complete analysis of the classical/quantum correlations encoded in the process tomography tensor, as well as an analysis of the main contributing elements to energy transport, is required and will be discussed elsewhere.

A variety of applied and foundational problems can be addressed once the process tomography tensor is reconstructed. From a foundational viewpoint, if the process tomography tensor is translated into the phase-space representation of quantum mechanics, it reduces to the propagator of the Wigner function.^{32–34} Based on this object, it is possible to *experimentally* reconstruct signatures of quantum chaos such as scars with *sub-Planckian resolution*.^{32} Phase-space resolution below ħ can be achieved here because the process tomography tensor, or equivalently the propagator of the Wigner function, is not a physical state, and therefore, it is not restricted by the uncertainty principle.^{32}

In the same way that 2D-PES was extended to study chemical exchange to obtain reaction rates under well controlled conditions (see, e.g., Chap. 10 in Ref. 13), a straightforward extension of QPT is the accurate measurement of concentration of different species in chemical reactions. This has been considered very recently in the literature.^{35} In this context, interest is in the population dynamics *χ*_{nnνν}(*T*) of the different chemical species, which under Markovian dynamics are in accordance to detailed balance and Onsager’s regression hypothesis.^{25,36,37} In this respect, because ultrafast spectroscopy allows for the study of chemical exchange with no need of pressure, temperature, pH, nor concentration jumps, it is expected that the proposed approach provides experimental evidence for the failure of the Onsager’s regression hypothesis induced by non-Markovian dynamics at the quantum level.^{25,36,37}

The QPT scheme introduced here can be readily implemented at the experimental level and constitutes a first step toward the formulation of QPT at the single-molecule level. Such a scheme would certainly incorporate quantum aspects of the electromagnetic radiation such as the use of time-energy-entangled photons.^{38} This is already under development in our laboratories. Finally, based on present non-linear optical activity spectroscopy (see, e.g., Chap. 16 in Ref. 13), by introducing time-polarization-entangled photons, instead of time-energy-entangled photons,^{38} single-molecule QPT may be extended to study optically active materials at the single-molecule level. These materials exhibit unique optical properties and are constantly finding applications in science and industry.

## Acknowledgments

Discussions with Keith Nelson and Joel Yuen-Zhou are acknowledged with pleasure. This work was supported by the Center for Excitonics, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, and Office of Basic Energy Sciences, under Award No. DE-SC0001088, by *Comité para el Desarrollo de la Investigación*—CODI—of Universidad de Antioquia, Colombia under the *Estrategia de Sostenibilidad* 2015-2016, by the *Colombian Institute for the Science and Technology Development*—COLCIENCIAS—under the Contract No. 111556934912 and by the National Science Foundation, Chemistry of Life Processes Program (No. CHE-1307272- to A.H.M).

### APPENDIX A: INITIAL STATE PREPARATION

Because the effective initial state $ \rho \u02c6 e 1 , e 1 \omega 1 , \omega 2 ( t + T ) $ is prepared by the first two pulses, it is of second order in *λ*. Thus, after applying second order perturbation theory to the time evolution of the system density matrix (see, e.g., Chap. 5 in Ref. 12), the effective initial state reads as

where it was assumed that $ \rho \u02c6 ( \u2212 \u221e ) =|g\u3009\u3008g|$. Symbols in calligraphic font denote superoperators; in particular, $ V \u0303 = \u2211 i = 1 4 V \u02c6 i ( t ) ,$ where $ V \u0303 i = [ V \u02c6 i , \u22c5 ] $. Assuming that the RWA holds, and that the rephasing signal is synchronously phase detected at *ϕ* = − *ϕ*_{1} + *ϕ*_{2} + *ϕ*_{3} − *ϕ*_{4}, the interaction with the electromagnetic radiation is conveniently described by

where $ \mu \u02c6 < = \u2211 \varpi p < \varpi q \mu p q |p\u3009\u3008q|$ promotes emissions from the ket and absorptions on the bra, and $ \mu \u02c6 > = ( \mu \u02c6 < ) \u2020 $ induces the opposite processes.

In the following, it is considered that the first and the second pulses are well separated, i.e., *τ* = *t*_{2} − *t*_{1} > 3*σ*. This allows for the substitutions 𝒱 (*t*″)=𝒱_{2} (*t*″) and 𝒱 (*t*″)=𝒱_{1} (*t*′),

If the duration of the pulse *σ* is much sorter than the dynamics induced by the environment characterized by Γ_{nm}, i.e., if $\sigma \u226a \Gamma n m \u2212 1 $, then decoherering contributions can be neglected so that $ G g p ( t 1 \u2212 t \u2032 ) \u2248exp i \varpi p g ( t 1 \u2212 t \u2032 ) $, $ G q p ( t 2 \u2212 t \u2033 ) G g p ( t \u2033 \u2212 t 2 ) \u2248exp \u2212 i \varpi q p ( t 2 \u2212 t \u2033 ) exp i \varpi p g ( t 2 \u2212 t \u2033 ) =exp \u2212 i \varpi q g ( t 2 \u2212 t \u2033 ) $, $ G g p ( t \u2033 \u2212 t 2 ) \u2248exp \u2212 i \varpi p g ( t \u2033 \u2212 t 2 ) $, and $ G g g ( t \u2033 \u2212 t 2 ) \u2248exp \u2212 i \varpi g g ( t \u2033 \u2212 t 2 ) $. Moreover, $ G g p ( t \u2033 \u2212 t \u2032 ) \u2248 G g p ( t \u2033 \u2212 t 2 ) G g p ( t 1 \u2212 t \u2032 ) G g p ( t 2 \u2032 \u2212 t 1 ) $. After some manipulations,

with

As mentioned in the main text, this effective initial state coincides with the one prepared by 2D-PES in Ref. 8.

### APPENDIX B: FINAL STATE DETECTION

To derive the explicit form of the density operator after the action of the four pulses, it is assumed that the third and the fourth pulses are well separated as well. To account for the action of the third pulse and a subsequent period of free evolution, perturbation theory is applied once more, so that the density operator of the system reads

where

Finally, the fourth pulse prepares the system in the state,

Each contribution can be easily associated to the double-sided Feynman diagrams in Fig. 1 in the main text.

Once the state of the system is obtained, the spectroscopy signals $ [ S FS ] e 1 , e 2 , e 3 , e 4 \omega 1 , \omega 2 , \omega 3 , \omega 4 ( \tau , T , t ) $, synchronously detected at *ϕ* = − *ϕ*_{1} + *ϕ*_{2} + *ϕ*_{3} − *ϕ*_{4}, follow from the calculation of $\u3008 A \u02c6 ( \tau + T + t ) \u3009=tr\u2009 A \u02c6 \rho \u02c6 e 1 , e 2 , e 3 , e 4 \omega 1 , \omega 2 , \omega 3 , \omega 4 ( \tau + T + t ) $ with $ A \u02c6 = \u2211 \nu = { e , e \u2032 , f} \Gamma \nu \nu \nu $. Specifically,

After replacing the explicit functional form of the density matrix elements $\u3008\nu \rho \u02c6 e 1 , e 2 , e 3 \omega 1 , \omega 2 , \omega 3 ( \tau + T + t ) \nu \u3009$ in Eq. (B4) and after conveniently collecting terms, Eq. (B4) leads to the 2D-FS signals in Eq. (14), which are the main result of this article.

### APPENDIX C: ISOTROPIC AVERAGES

Before proceeding to the calculation of the isotropic average, it is necessary to express the dipole transition operators in the molecular frame. In doing so, take as reference the transition dipole operator *μ*_{eg} = *μ _{eg}*

**m**

_{z}. Hence,

*μ*_{e′g}=

*μ*

_{e′g}cos(

*θ*

_{e′g})

**m**

_{z}+

*μ*

_{e′g}sin(

*θ*

_{e′g})

**m**

_{x},

*μ*_{fe}=

*μ*cos(

_{fe}*θ*)

_{fe}**m**

_{z}+

*μ*sin(

_{fe}*θ*)

_{fe}**m**

_{x}, and

*μ*_{fe′}=

*μ*

_{fe′}cos(

*θ*

_{fe′})

**m**

_{z}+

*μ*

_{fe′}sin(

*θ*

_{fe′})

**m**

_{x}. The angle between the different transition dipole moments is given by

The isotropically averaged signals can then be written in the compact form

with

In defining the vector **χ**(*T*), besides its properties in Eqs. (2)-(4), the following relations were instrumental, *χ _{ggee}*(

*T*) − 1 = −

*χ*(

_{eeee}*T*) −

*χ*

_{e′e′ee}(

*T*),

*χ*

_{gge′e′}(

*T*) − 1 = −

*χ*

_{eee′e′}(

*T*) −

*χ*

_{e′e′e′e′}(

*T*), and

*χ*

_{ggee′}(

*T*) = −

*χ*

_{eeee′}(

*T*) −

*χ*

_{e′e′ee′}(

*T*).