Quantum process tomography of molecular dimers from two-dimensional electronic spectroscopy I: General theory and application to homodimers

Is it possible to infer the time evolving quantum state of a multichromophoric system from a sequence of two-dimensional electronic spectra (2D-ES) as a function of waiting time? Here we provide a positive answer for a tractable model system: a coupled dimer. After exhaustively enumerating the Liouville pathways associated to each peak in the 2D-ES, we argue that by judiciously combining the information from a series of experiments varying the polarization and frequency components of the pulses, detailed information at the amplitude level about the input and output quantum states at the waiting time can be obtained. This possibility yields a quantum process tomography (QPT) of the single-exciton manifold, which completely characterizes the open quantum system dynamics through the reconstruction of the process matrix. This is the first of a series of two articles. In this manuscript, we specialize our results to the case of a homodimer, where we prove that signals stemming from coherence to population transfer and viceversa vanish upon isotropic averaging, and therefore, only a partial QPT is possible in this case. However, this fact simplifies the spectra, and it follows that only two polarization controlled experiments (and no pulse-shaping requirements) suffice to yield the elements of the process matrix which survive under isotropic averaging. The angle between the two site transition dipole moments is self-consistently obtained from the 2D-ES. Model calculations are presented, as well as an error analysis in terms of the angle between the dipoles and peak overlap. In the second article accompanying this study, we numerically exemplify the theory for heterodimers and carry out a detailed error analysis for such case. This investigation provides an important benchmark for more complex proposals of quantum process tomography (QPT) via multidimensional spectroscopic experiments.

to the case of a homodimer, where we prove that signals stemming from coherence to population transfer and viceversa vanish upon isotropic averaging, and therefore, only a partial QPT is possible in this case.
However, this fact simplifies the spectra, and it follows that only two polarization controlled experiments (and no pulse-shaping requirements) suffice to yield the elements of the process matrix which survive under isotropic averaging.The angle between the two site transition dipole moments is self-consistently obtained from the 2D-ES.Model calculations are presented, as well as an error analysis in terms of the angle between the dipoles and peak overlap.In the second article accompanying this study, we numerically exemplify the theory for heterodimers and carry out a detailed error analysis for such case.This investigation provides an important benchmark for more complex proposals of quantum process tomography (QPT) via multidimensional spectroscopic experiments.
Multidimensional optical spectroscopies (MDOS) provide very powerful tools to study excited state dynamics of multichromophoric systems in condensed phases.These techniques distribute spectral features along several dimensions, uncluttering data which would otherwise appear obscured in linear spectroscopies and simultaneously yielding novel information on the dynamics of the probed system [1].Possibilities in multidimensional techniques include decongesting spectral lineshapes, differentiating between homogeneous and inhomogeneous broadening mechanisms, providing unambiguous signatures about couplings between chromophores, and yielding signatures of coherent and incoherent processes involving excited states at the amplitude level [2,3].Although MDOS have been historically inspired by their NMR analogues, the timescales of the physical and chemical processes studied through MDOS are quite different from the ones in NMR [4][5][6][7].The characteristic timescales of NMR are milliseconds, a resolution that does not allow for the observation of a wide variety of chemical dynamics in condensed phases ocurring in the orders of femto and picoseconds.On the other hand, femtosecond timescales can be easily accessed with ultrafast optical techniques.Examples of phenomena studied via MDOS are vast and include molecular reorientation processes and solvation dynamics [8,9], electron transfer [10], vibrational coherences in organometallic complexes [11][12][13] or halogens in rare gas matrices [14,15], phonon dynamics in carbon nanotubes [16], protein unfolding kinetics [17], excitonic dynamics in light-harvesting systems [18][19][20][21][22][23] and organic polymers [24,25], as well as many-body physics in semiconductor quantum wells [26][27][28] and quantum dots [29].
Traditionally, the spectroscopy of condensed phases is formulated as a response problem: The molecular system is perturbed with a sequence of short laser pulses and the coherent polarization response due to this set of perturbations (nonlinear polarization) is subsequently measured [1].Informally, we can describe the exercise as 'kicking' the quantum black box (molecular system) and 'listening to the whispers' (measuring the response) due to the kicks, from which some properties of the box can be inferred.This description of spectroscopy is reminiscent to an idea stemming from the quantum optics and quantum information processing (QIP) communities, namely, quantum process tomography (QPT) [30][31][32][33].Broadly speaking, QPT is a systematic procedure to characterize a quantum black box by sending a set of inputs, measuring their outputs, and analyzing the functional relationships between them.With the increasing effort of quantum engineering of gates and devices, QPT constitutes a cornerstone of QIP theory and experiment, as it provides a necessary check on the performance of the respective quantum black boxes.A natural question arises from the comparison of the two aforementioned concepts: Can the spectroscopy of condensed I. RELEVANT CONCEPTS OF QUANTUM PROCESS TOMOGRAPHY AND GEN-ERAL DEFINITIONS Consider an arbitrary quantum system (quantum black box) interacting with an environment.
We are interested in its evolution as a function of time T in the form of a reduced density matrix ρ(T ).Very generally, this evolution is a linear transformation on the initial quantum state [35,36]: (1) χ(T ) is the central object of this article, and shall be called process matrix.Eq. (1) can be regarded as an integrated equation of motion for every T [101].Eq. ( 1) can be expressed in terms of a basis for the Liouville space of the system: For purposes of this article, we present two useful definitions.Consider the Liouville space L of the system, and classify the vectors of L in proper and improper density matrices.A state or a density matrix is proper if it satisfies all the conditions of a physical quantum state; namely, this is Hermitian, positive semidefinite, and has trace one.An improper state is any matrix that lives in the Liouville space but is not a proper density matrix.Clearly, any improper density matrix in the same Liouville space may be written as a unique linear combination of proper density matrices.
In principle, Eq. (2), being a physical equation of motion, is restricted to the domain of proper density matrices ρ(0).However, by linearity, its extension to any linear combination of proper states is well defined, so its validity for improper density matrices is not under question.
The meaning of the process matrix χ(T ) is easy to grasp: Conditional on the initially state being prepared at ρ(0) = |c d|, χ abcd (T ) is the value of the entry ab of the quantum state after time T , ρ(T ), i.e. χ abcd (T ) = a|ρ(T )|b .Therefore, χ abcd (T ) denotes a state to state transfer amplitude.Note that ρ(0) = |c d| is an improper density matrix if c = d (coherences on their own are not valid quantum states).However, improper states are not necessarily unphysical as one expects at a first glance.Most of our intuition for nonlinear spectroscopies in the perturbative regime stems from the consideration of how a perturbative amplitude created at a certain entry |c d| of the total (proper) density matrix is transferred to other entries due to free evolution, as time progresses [1].It is not the evolution of the total density matrix (which to leading order is unperturbed, mostly in its ground state, and not yielding a time-dependent dipole) what is effectively monitored in the phase-matched signal, but the evolution of an effective density matrix, such as |c d|, which can be improper.Terms such as transfer from population to population, coherence to coherence, population to coherence, and coherence to population are all ubiquitous in the jargon of MDOS.However, the monitoring of the latter is often ambiguous, incomplete, and in most cases, qualitative.Obtaining quantitative information about these events amounts to finding each of the elements of χ(T ).
The transformation in Eq. ( 2) is limited by two classes of restrictions for the process matrix associated with Hermiticity and trace preservation: We derive these conditions in Appendix A, but their content is intuitive: If ρ(0) is a proper density matrix, ρ(t) remains as a valid quantum state as T evolves if these two requirements are preserved.
In particular, elements of the form χ aabb (T ), which denote population transfers, are purely real as one expects, whereas the other elements are in general complex.As a comment to our previous discussion, by linearity, these conditions must also be satisfied even if ρ(0) is improper (notice that χ(T ) does not depend on ρ(0)).
Equations (1) and (2) are remarkable because they guarantee that, in principle, if χ(T ) is known, the quantum black box described by ρ(T ) is perfectly understood, as it predicts by linearity the evolution of an arbitrary initial state of L. Although ρ(t) describes an open quantum system, details about the environment evolution need not be included explicitly, but only in an averaged sense in the elements of χ(T ).We shall operationally define QPT as any procedure to reconstruct χ(T ).A possible QPT is the following: (a) Prepare a linearly independent set of states ρ(0) that spans L; (b) for each of the prepared states, wait for a free evolution time T and determine the density matrix at that time.Any protocol for determining a density matrix for a system is called Quantum State Tomography (QST) [37][38][39][40][41].In essence, QPT can be carried out for any system if both a selective preparation of initial states and QST can be achieved.Variants of this methodology exist although all of them operate within the same spirit [30][31][32][33].QPT has been successfully implemented in a wide variety of experimental scenarios, including nuclear magnetic resonance [42][43][44], ion traps [45], single photons [46,47], solid state qubits [48], optical lattices [49], and Josephson junctions [50].In this article, we show how to perform QPT for a model coupled heterodimer using two-color polarization controlled heterodyne photon-echo experiments, extending the domain of application of QPT to systems of chemical and biophysical interest.

II. MODEL SYSTEM: COUPLED DIMER
Consider a molecular dimer described by the effective Hamiltonian [3,51,52]: where a + i and a i are creation and anhilation operators for a single Frenkel exciton in the site i ∈ {A, B}, ω A , ω B are the first and second site energies, and J is the coupling between the chromophores.
The standard diagonalization of this Hamiltonian, which is effectively a two-level system for the single-exciton manifold, follows from defining some convenient parameters: The average of the site , and the mixing angle θ = 1 2 arctan J ∆ .By introducing the operators: the Hamiltonian in Equation ( 5) can be readily written as: where the eigenvalues ω α and ω β of the single excitons are: , also plays a role in our study, as it is resonantly accesed through excited state absorption (ESA) after several pulses.Notice that the exciton binding or repulsion terms, so the energy level of the biexciton is just the sum of the two exciton energies, . Defining ω ij ≡ ω i − ω j , the following relations hold: Since we are concerned with the interaction of the chromophores with electromagnetic radiation, we make some remarks on the geometry of the transition dipoles (see Fig. 1).Let µ ij = i| μ|j .
Assume that the transition dipole moments from the ground to the single excitons in the site basis are µ gA = µ Ag = d A and µ gB = µ Bg = d B , respectively.It follows that the dipole moments µ ij for i, j ∈ {α, β, f } are located in the same plane, but in general have different magnitudes and directions: In this model, we shall consider µ ij = µ ji .As enumerated in our model, dipole mediated transitions only couple the ground state to the single excitons, and the single excitons to the biexciton.

III. PHOTON-ECHO EXPERIMENT AS QUANTUM PROCESS TOMOGRAPHY
Consider a four-wave mixing experiment where an ensemble of identical dimers interacts with a series of three ultrashort laser pulses.The perturbation due to these pulses is given by: where λ is the intensity of the electric field, which is assumed to be weak, μ is the dipole operator, e i , t i , k i , ω i denote the polarization [102], time center, wavector, and carrier frequency of the i − th pulse, and r is the position of the center of mass of the molecule.E(t) is the slowly varying in time pulse envelope, which we choose as a Gaussian with width σ, or full-width half-maximum where k s = lk 1 + mk 2 + nk 3 and l, n, m are integer numbers.Notice that k s equals to a linear combination of the wavevectors associated with each pulse.As expected, the action of a pulse on each molecule attaches a spatial phase to its quantum state, so the total phase accumulated by it equals e iks•r for each combination of perturbations.Each improper density matrix ρ s (t) corresponds to one of these phases, and can be calculated by keeping track of the actions of the pulses in the bra and the ket of the system using double-sided Feynman diagrams.Eq. (12) implies that the optical polarization induced on the molecule can also be Fourier decomposed into different components [53][54][55]: P (r, t) = T r( μ(r)ρ(r, t)) = s P s (t)e iks•r , where μ(r) denotes the dipole operator of the molecule located at r.The experimental setting we describe is analogous to the one of an array of dipole antennas which are spatially phased in a grating with respect to each other and oscillate in time.Classical electromagnetism predicts that the induced macroscopic polarization of this array emits radiation which is precisely concentrated along the vectors k s .This condition, which reflects conservation of momentum of the fields, is known as phase-matching [56].
A fourth pulse of the same wavevector as one of the k s , known as the local oscillator, is allowed to interfere with the radiation along that direction.By varying the phases of this fourth field, two heterodyne detections can be carried out to extract the real and imaginary components of , where e 4 is the polarization of the local oscillator [1].
In this article, we are interested in the signal along k P E = −k 1 + k 2 + k 3 , the so called photonecho (PE) direction [57].The frequency components of the pulses lie within the optical regime, so they can induce the transitions enumerated in the previous section.Traditionally, in the MDOS literature, the intervals between the time centers of the pulses are called coherence τ = t 2 − t 1 , waiting T = t 3 − t 2 , and echo t = t 4 − t 3 times, respectively.Here, t 4 is the time of detection of the signal [58].We shall only consider rephasing photon-echo signals, where where the inhomogeneous broadening is rephased [59].Due to these explicit interval dependences, the collected signal can be expressed as P P E (τ , T, t).As explained in our previous study [34], we may regard the PE experiment as a QPT of the single-exciton manifold dynamics of the dimer as a function of T .In fact, the polarization signal may be expressed as a linear combination of elements of the process matrix χ(T ): where, and, The coefficients C p ω i for p ∈ {α, β} are frequency amplitudes of the laser pulse which is centered at ω i , evaluated at the transition energy ω pg : and is the propagator of the optical coherences |i j| in the coherence and echo times, which, has been taken to be the product of a coherent oscillatory term beating at a frequency ω ij and an exponential decay with dephasing rate Γ ij .This propagator is defined only for τ > 0 via the step function Θ(τ ).
The frequencies of the coherences in the coherence and echo intervals have opposite signs, reflecting the rephasing character of the signal.In optical PE experiments, it is customary to assume that the free-induction decay characterized by the evolution of optical coherences in the coherence and echo times is well characterized, and given by expressions of the form (17).The reason is that the characteristic energetic scales of the vibrational degrees of freedom are much lower than the optical gap, so the only nonunitary dynamics they induce in the optical coherence is, to a good approximation, restricted to pure dephasing Γ ij which can be inferred from the polarization signal [103].The dynamics in the waiting time is more complex, consisting of small frequencies due to excitonic superpositions which are strongly influenced by the bath.It is the latter interval where QPT will prove most useful.
The polarization signal yields a linear combination of elements χ abcd (T ) weighted by the probability amplitude to prepare a state |c d| with the first two pulses and detect |a b| with the third pulse and the fourth heterodyning pulse.These probability amplitudes can be controlled by manipulating the polarization of the pulses e i and the frequency amplitudes for the resonant In essence, state preparation and QST are implicit in the coherence and echo times (see [34]).In a different context, Gelin and Kosov had previously hinted at a similar idea by identifying these times as "doorway" and "window" intervals [60].By conducting several experiments varying these control knobs and collecting the signal from each of these settings, a system of linear equations can be established whereby the elements of χ(T ) can be inverted, and therefore QPT is achieved.This statement is correct provided that besides the free-induction decay rates Γ ij , the parameters ω αg , ω βg , µ αg , µ αg , µ αg , and µ αg are all known or can be obtained selfconsistently during the experiment.We will elaborate on these points for the case of a homodimer in section V.
Notice that Eqs. ( 13), (14), and ( 15) monitor all the 12 real valued paramenters involving χ abcd (T ) for a, b, c, d ∈ {α, β}, so that they allow for the QPT of the single exciton manifold, which is an effective quantum bit (qubit) system.However, we have also kept track of the elements χ ggcd (T ) c, d ∈ {α, β}, that is, the possibility of amplitude leakage errors from the single-exciton dynamics channel to |g g|.It is known that whereas the excitonic dynamics occurs in femtosecond timescales, exciton recombination happens in the order of nanoseconds.Therefore, these decay channels could be potentially ignored in many experimental systems.We shall keep them in our theoretical analysis as they do not increase the complexity of the problem by much, although in situations where this could be problematic, we could accordingly disregard them.

IV. QPT FROM 2D SPECTRUM OF PE
As mentioned, QPT can be carried out from data resulting from a series of experiments varying colors and polarizations of the pulses.The necessary information can in principle be obtained by collecting a single point for a fixed pair of τ and t points for each of the experiments.Often, however, the PE signal is collected across many τ, T, t points, and conveniently processed into a 2D correlation spectrum in the conjugate frequency variables ω τ and ω t : which still evolves in the T coordinate [104].By performing the integrals of Eq. ( 18) using Eqs.
The peaks are centered about ω = ω mg and have a width parameter Γ mg .The difference in signs for the Fourier transform in Eq. ( 18) guarantees that all the resonances appear in the first quadrant of both frequency axes.The expressions for the amplitudes S mn (T ), associated with peaks centered at (ω τ , ω t ) = (ω mg , ω ng ), are given by [106]: Typically, the probed samples are in solution, so the molecules in the ensemble are isotropically distributed.The isotropic average ) is given by [61]: e 1 e 2 e 3 e 4 ;m where e i and m i are the polarizations of the pulses in the lab and the molecular frame, respectively.
The isotropic average consists of a sum of molecular frame products weighted by the isotropically invariant tensor I e 1 e 2 e 3 e 4 ;m 1 m 2 m 3 m 4 .Since the information in Eqs. ( 13), ( 14), ( 15) is in principle contained in Eqs. ( 13), (14), and (15), several conclusions from our previous study are immediately transferable: The elements of χ(T ) contained in Eqs. ( 22), ( 23), (24), and ( 25) can be all be extracted by repeating a number of experiments with different polarization configurations for the fields and two different waveforms for the pulses.Under different motivations, theoretical proposals for manipulating 2D-ES using pulseshaping capabilities have been previously reported [62,63].An extensive study of this possibility for a heterodimer will be presented in the second article accompanying this study.
Eqs. ( 22), ( 23), (24), and (25) can also be derived by book-keeping the double-sided Feynman diagrams that oscillate at the particular frequencies for the coherence and waiting times in each of the four resonances (we refer the reader to Fig. 1).In analyzing the possible pathways in Liouville space, we make use of the rotating-wave approximation (RWA): Perturbations which are proportional to e −ik i •r+iω i (t−t i ) μ•e i can excite the ket and de-excite the bra, whereas the ones proportional to e ik i •r−iω i (t−t i ) μ • e i can deexcite the ket and excite the bra.As an illustration, consider the signal S αβ (T ), which arises from diagrams oscillating with frequency ω gα at the coherence time and ω βg at the echo time (Fig. 1   Therefore, in principle, there are 4 × 3 = 12 possibilities for χ abcd (T ) which can be detected in at S αβ (T ).However, we assume that the state |g g| does not evolve to other states due to the bath: This assumption is quite reasonable, as we are ignoring processes where phonons can induce optical excitations from |g g|.This condition is present in Eqs.(14), and (15) in the δ−function terms and in Eqs. ( 22), ( 23), (24), and (25) in the "-1" terms, which correspond to −χ gggg (T ).This leaves 12 − 4 = 8 possibilities for χ abcd (T ).
To be more explicit, consider the pathways in S αβ (T ) that monitor the population to coherence process χ βααα (T ).These are displayed in Fig. (2).The pathway on the left involves represents an excited state absorption (ESA) from the single-exciton manifold, and is proportional to ), an expression which can be immediately read out from the diagram: Each interaction with the field picks up a factor corresponding to the amplitude of the transition, which depends on the alignment of the corresponding dipole with the polarization of the pulse, as well as the frequency amplitude of the pulse at the given transition.A minus sign is included if the perturbation is on the bra.Similarly, the pathway on the right involves stimulated emission (SE) and is proportional to ).The rest of the diagrams for all the peaks can be systematically analyzed in the way described above.In general, the pathways we need to consider can be classified in ESA, SE, and ground state bleaching (GSB) processes.GSB processes are the ones that take |g g| at the end of the waiting time to a dipole active coherence involving an excited state.ESA pathways, which are proportional to dipole transitions involving the excited state, differ in sign from SE and GSB pathways, as can be easily seen by inspection.
Fig. (3) provides a mnemonic device to keep track of the Liouville pathways that each peak in the 2D electronic spectrum monitors, and therefore, also provides a scheme of the QPT protocol.
The ω τ axis can be associated with a particular state preparation whereas the ω t axis with a particular detection.Each peak reflects a nontrivial number of processes in Liouville space.As an illustration (see Fig. ( 4)), we consider the ideal case where the bath does not interact with the system, in which case, a very simple picture is recovered: The off-diagonal peaks beat at the coherence frequency and the diagonals remain static.This case can be easily derived from Eqs. (22), ( 23), (24), and (25) by substituting χ abcd (T ) = δ ab δ cd + δ ac δ bd e −iω ab T , that is, populations remain static whereas coherences beat at difference frequencies.

V. THE CASE OF THE HOMODIMER
To gain insights into the described QPT protocol, we specialize the results above to a coupled homodimer.In the following subsections, we discuss, for this particular case, (A) the Hamiltonian and the transition dipole moments involved in the experiments, (B) properties of the spectroscopic signals under isotropic averaging, (C) stability of the numerical inversion, (D) analytical expressions of the elements of χ(T ) in terms of the peak amplitudes of the spectra, (E) a procedure to extract the angle φ between the dipoles, (F) a summary of the QPT procedure, and (G) a numerical example with a model system.A similar study focused on the heterodimer will be presented in the second part of this series.

A. Hamiltonian and transition dipole moments
In the homodimer, the two sites are identical chromophores with energies ω A = ω B = ω, and the Hamiltonian in Eq. ( 5) and ( 7) is given by: which we have diagonalized with the symmetric a + β |g and antisymmetric a + α |g single exciton states given by: The splitting between the two delocalized states is just 2J.Using Eq. ( 10), the transition dipoles take the simple forms: , Interestingly, these expressions are independent of the coupling J. Also, notice that µ αg and µ f α are perpendicular to µ βg and µ f β (see Fig. 5).Denoting the norm of each site dipole by the following relationships follow: As expected, in the degenerate limit that φ = 0 or π (the site dipoles are parallel or antiparallel), one of the delocalized excitons becomes dark and there is only one bright transition from the ground state.If this is not the case, in general, the two transitions are bright and their dipoles perpendicular to each other.Furthermore, as a difference with the heterodimer case, there are only three (instead of four) different transition dipoles in the homodimer, and two of them are just negatives of each other.The degenerate case will be discussed as a limit of the more general one in the next paragraphs.

An important observation regarding isotropic averaging follows:
Claim.-Upon isotropic averaging, signals stemming from coherence to population or population to coherence transfer cannot be monitored in the 2D-PE spectrum of a homodimer.
Proof.-For simplicity, align µ βg and µ αg in the y and z directions in the frame of the molecule and use Eq. ( 26) to argue that the only terms in the sum that contribute to an isotropic averaging are the ones where only two or four polarizations of the field are the same.In Eqs. ( 22), ( 23), (24), and ( 25), all the dipole-polarization terms corresponding to coherence to population and the opposite processes involve three dipoles of the same kind and a perpendicular one.Therefore, they vanish under isotropic averaging.As an example, consider the the terms associated with χ βααα (T ) in S αβ (T ): The claim above allows for a considerable simplification of Eqs. ( 22), ( 23), (24), and (25).Each of the peaks in the 2D spectrum keeps track of less elements of the process matrix χ(T ) upon isotropic averaging: Only population to population and coherence to coherence transfers can be monitored.The results for the peaks of an isotropically averaged 2D spectra are presented below.
We have taken the shortcut notation • e 1 e 2 e 3 e 4 , which denotes the isotropically averaged signal stemming from the two pulse polarization configurations (e 1 , e 2 , e 3 , e 4 ) = (z, z, z, z), (z, z, x, x), so that the terms S mn (T ) e 1 e 2 e 3 e 4 and S(ω τ , T, ω t ) e 1 e 2 e 3 e 4 have the obvious meanings.TABLE 1. Isotropically averaged 2D-ES peak amplitudes for the zzzz configuration TABLE 2. Isotropically averaged 2D-ES peak amplitudes for the zzxx configuration We focus our attention on experiments with short pulses that are broadband enough to create either exciton |α or |β with the same amplitude, that is, C p ω i = C, for a purely imaginary constant C, for all p and ω i .This condition can be easily relaxed, but we shall proceed with it to analyze our QPT protocol with more detail.Using the condition of Eq. ( 3), we can eliminate the variables χ ggαα (T ) and χ ggββ (T ) for χ αααα (T ), χ ββαα (T ), χ ααββ (T ), χ ββββ (T ).Also, taking advantage of Eq. ( 4), we discard χ βαβα (T ) and χ αββα (T ) and keep χ αβαβ (T ) and χ βααβ (T ).From the left hand side (LHS) of the real and imaginary parts of the spectrum (see Tables 1 and 2), we derive the following real valued matrix equation: Similarly, the right hand side (RHS) of the spectra yields: Inverting Eqs. ( 34) and ( 35) yields most of the elements of χ(T ) involving the single-exciton manifold.Whereas the presented QPT for the homodimer is partial, no complicated pulse shaping efforts need to be carried out.Instead, the requirement is standard pulse polarization control achievable with current experimental capabilities [11,12,[64][65][66][67][68][69].
The transition dipole moments must be well characterized in order to construct the matrices in Eqs.(34) and (35).This requirement is self-consistently fulfilled by collecting the spectra in the collinear and cross-polarized configurations.Notice that χ αααα (T ) and χ ββαα (T ) are exclusively monitored in the LHS, and χ ββββ (T ) and χ ααββ (T ) only detected in the RHS.However, coherence transfer terms χ αβαβ (T ) and χ βααβ (T ) are repeatedly monitored in different peaks in both sides of the spectra.Due to this repetition, there are redudant equations that allow for the self-consistent extraction of the angle φ without compromising the inversion of the elements of χ(T ).Details about this parameter extraction are developed in subsection E. For the time being, we assume that the information about the transition dipoles is previously known.

C. Stability of the Quantum Process Tomography protocol for a homodimer
In order to characterize the stability of inversion of χ(T ), we can arrange Eqs.(34) and (35) into a single matrix equation M χ(T ) = S(T ), where M is a 16 × 8 matrix of dipole moments, χ(T ) is a vector of 8 unknowns, and S(T ) is a vector of 16 real valued amplitudes extracted from the signal.Denoting with • the spectral norm of a vector or a matrix [70], we obtain a bound on the relative error of the inverted vector χ(T ) which yields the QPT: where ∆ χ(T ) and ∆S(T ) denote errors in χ(T ) and S(T ) upon inversion, respectively, and the condition number κ is given by: The lowest possible value for a condition number is κ = 1.In Fig. (6), a plot of κ vs. φ (red) is displayed.As expected, very large values of κ, which denote unstable inversions, are expected for systems where the site dipoles are aligned or antialigned, when one of the eigenstates of H S becomes dark.We also indicate a range of angles where κ is below a reasonable threshold, say κ ≤ 15 (blue horizontal line), which consists of angles in the range 0.3π ≤ φ ≤ 0.7π.Also, note that κ is symmetric about the minimum of κ at φ = π/2, where the best inversion is carried out with κ ∼ 3.9.

D. Analytical expressions for χ(T )
To gain insights into the QPT protocol, we shall derive explicit expressions for the elements of χ(T ) in terms of the amplitudes of the spectra and the angle φ between the dipoles.First, we substitute Eqs. ( 32) and ( 33) into the LHS of the spectra through Eq. (34).The following expressions are obtained after inverting the resulting matrix equation: ' Similarly, the RHS of the spectra through Eq. ( 35) yield: ' As expected, the RHS equations above can be obtained from the LHS equations upon the substitutions α → β, β → α, and φ → π − φ (see Eqs. ( 33)).
From Eqs. ( 43) and ( 50), we notice that the imaginary parts of S αβ (T ) zzzz , S αβ (T ) zzxx , S βα (T ) zzxx , S βα (T ) zzxx are all exclusively proportional to A careful study of the precise bath-induced dynamics in this system is beyond the scope of this study and shall be addressed in future work in collaboration with an experimental realization.
We consider a harmonic bath with an Ohmic spectral density: J(ω) = λ ωc ωe −ω/ωc , with λ = 100 cm −1 , ω c = 150 cm −1 at a temperature T = 273 K. Identical baths are assumed to be diagonally and linearly coupled to each of the sites.We closely follow the calculation reported in [72] and adapted in [34].The dynamics of the total excitonic system, which is a proper density matrix, is governed by the following equation of motion: where R denotes the time-independent sparse dissipative superoperator containing only a few non-zero elements listed in Table 3.Since ρ(T ) only depends on ρ(T ) and not on the value of the quantum state at previous times, the simulated dynamics are Markovian.
It is well known that the secular Redfield equations guarantee thermal equilibrium since the population transfer rates satisfy R ααββ /R ββαα = e −ω αβ /k B T , where k B is the Boltzmann constant.
Also, R f gf g will not be relevant for the calculations, as coherences between the ground state and the biexciton are never created in the PE experiment.The free-induction decay rates for the coherence and echo intervals will be taken for simplicity to be the same, Γ αg ≈ Γ βg ≈ 1 2 (R αgαg + R βgβg ).This restriction is by no means necessary, but will simplify the simulations below.
The picture of the secular Redfield equations is very simple and provides transparent means for understanding the QPT protocol for the homodimer: The evolution of populations and the coherences independently satisfy standard first-order kinetic equations, leading to multiexponential integrated dynamics.
In Fig. 7, we display the calculated 2D-ES of this model system.We consider the three pulses to be identical, centered about From the simulated spectra S(ω τ , T, ω t ) e 1 e 2 e 3 e 4 , the extraction of the terms S mn (T ) e 1 e 2 e 3 e 4 is achieved with high fidelity (>99%) by a nonlinear optimization routine based on the simplex  search method with bound constraints [73].The signals are fitted to a sum of four different resonances as in the isotropically averaged version of Eq. ( 19).The parameters ω mg , ω ng , Γ, and S mn (T ) e 1 e 2 e 3 e 4 are reconstructed from 2D-ES with a grid spacing of ∆ω τ = ∆ω t = 1 cm −1 and a grid size of 1050 cm −1 for every axis.We present the results of this calculation in Figure 8.Notice that the imaginary parts of the diagonal peaks are zero since no terms of the form χ βααβ (T ) = (χ αββα (T )) * are considered in the secular Redfield theory, and population transfer terms are purely real.
Eq. ( 52) is solved at every waiting time T ∈ [T c /2, 10T c ] yielding the roots ξ = 1.000, 0.4059 for every T , without variance after the fourth decimal digit, indicating its robustness for the inversion of φ.The same exercise with Eq. ( 53) gives ξ = −0.4059,0.4059.These values of ξ imply φ = 90 0 , 65 0 , 65 0 , 86.3i 0 , respectively.The last value can be discarded for it is not even a real number.The value of φ must be a root of both equations.Therefore, we can also discard 90 0 , since it is only solution of the first equation, but not of the second.The result φ = 65 0 follows unambiguously, as expected.
Finally, the terms S mn (T ) e 1 e 2 e 3 e 4 and the angle φ allow for the evaluation of the elements of χ(T ) which are extractable for the homodimer.Figure 9 shows that this reconstruction coincides with the analytical expressions presented in Table 4.The population decay terms χ αααα (T ), χ ββββ (T ) both start at 1 and reach 0 exponentially, the second faster than the first, since |β is the excitonic state of higher energy.The population transfer terms χ ααββ (T ), χ ββαα (T ) are complementary to the former ones, with the transfer from |β to |α being faster for the same reasons just mentioned.The coherence term decays exponentially, with real and imaginary parts π/2 phase shifted one from another.The calculated timescale of this decay (hundreds of femtoseconds) is similar to the one inferred from the experiment reported by Lee and coworkers, where a superposition of excitons in the bacteriopheophytin and bacteriochlorophyll sites in the reaction center of purple bacteria is monitored indirectly through a two-color experiment [74].

VI. DISCUSSION
In the present article, we have outlined a general theory for carrying out a QPT for a molecular dimer using the information contained in various frequency and polarization controlled 2D-ES.We started by providing the basic concepts of QPT, and operationally defined a QPT as a protocol to extract the process matrix χ(T ), which in principle completely characterizes a quantum black box, in our case, the box being the single-exciton manifold of the dimer.After reviewing the model Hamiltonian as well as the transition dipole moments of an excitonic dimer, we adapted the QPT theory presented in our previous work, where the nonlinear polarization was analyzed in real time for single values of τ and t times (see Eqs. ( 13), (14), and ( 15)) [34], to the more standard and visual Fourier transformed signal collected along several values of these interval times.The central result of this exercise was Eqs. ( 22), ( 23), (24), and ( 25), which from a purist stanpoint completes the QPT effort: The peaks in a heterodyne-detected 2D-ES can be expressed as linear combinations of elements of the process matrix χ(T ).This information can be distilled by carrying out several experiments alternating the frequency components of the pulses as well as their polarization.By setting up a system of linear equations with this data, a linear algebraic routine yields the inversion of χ(T ) for every waiting time T .In order to get a more intuitive picture of this procedure, the particular case of a homodimer was studied in detail.The degeneracies of this system yield a perpendicular set of transition dipole moments which considerably simplify the theory (see Fig. 5).It was shown that under isotropic average of the signal, no population to coherence processes or viceversa can be monitored, impeding a full QPT for the single-exciton manifold of this system.However, the partially achievable QPT is very simple, robust with respect to transition dipole moment parameters as long as they are not aligned or antialigned, and readily implemented without pulse shaping.The only requirement is the collection of two polarization controlled 2D-ES.Numerical examples with a model homodimer validated the presented theory.
The possibilities that QPT opens for the study of excited state dynamics in condensed molecular systems are as vast as the information acquired at the amplitude level of the evolving quantum state of the probed system.On the one hand, with the peaks in the 2D spectrum indicating a plethora of pathways in Liouville space, understanding of the dynamics is undoubtedly enhanced by the dissection of these peaks into processes described by the χ(T ) matrix.On the other hand, a plethora of questions can be addressed with this information, for instance: Is a Markovian description accurate?[75][76][77][78] If not, what is the degree of non-Markovianity of the dynamics?[79,80] If it is Markovian, is the secular Redfield equation appropriate or are non-secular processes important?[81].Is there any degree of entanglement in the quantum states produced in the single exciton manifold upon photoexcitation?[82,83] What is the rate of decoherence of a quantum superposition between excitonic states?[84][85][86][87][88] A few aspects have not been fully addressed with respect to the implementation of QPT of a molecular dimer.These issues will be carefully studied in future publications in collaboration with experimental groups.The role of static disorder in the eigenenergies of the system as well as in the distribution of the angle φ will necessarily yield an inhomogenously averaged signal from which the relevant information must be carefully extracted.We anticipate this feature to add another step of parameter fitting, but not change the results of our theory dramatically.
Furthermore, we have ignored the possibility of resolving the vibronic structure accompanying each of the four resonances in the 2D-ES.If this were to happen, it might be wiser to take the approach of Cina and coworkers [89][90][91] to consider the evolution of the nuclear wavepackets for a few modes strongly coupled to the system, and maybe regard the rest of the modes as a bosonic bath.This possibility would require an exponential increase of experimental resources [33], so either partial or compressed sensing approaches [92,93] would be necessary.Alternatively, by going back to the time-domain picture provided by the authors in their previous work [34], and applying novel concepts of QPT for initially correlated states [94][95][96], a coarse grained and consistent tomographic protocol could be designed to address this problem.Finally, it might be worth considering additional nonlinear optical spectroscopic techniques such as considering the analysis of both rephasing and non-rephasing signals [59,97], transient grating [24], pump probe [11], or phase cycling of multipulse induced fluorescence [98] to investigate if they provide additional information for a more robust QPT.
Albeit this article not being exhaustive, we hope to have convinced the reader that the QPT approach follows the spirit of MDOS in a very natural way.By systematically studying excited state dynamics as a quantum black box, an intriguing perspective on MDOS emerges that allows the use of tools designed in the QIP community in order to study excited state dynamics of condensed molecular systems.These possibilities will be the subject of future studies.

VII. ACKNOWLEDGMENTS
We wish to acknowledge valuable discussions with Dr. Jacob Krich and Dr. Paul Peng on nonlinear spectroscopy, Seungwan Ryu on time-bandwidth issues, and Patrick Wen, Dylan Arias, Gabriela Schlau-Cohen, and Eleonora de Re for insights on experimental realizations of QPT.This work was funded by DARPA-QUBE as well as the Center for Excitonics.
APPENDIX: Derivation of Eqs. ( 6), (4), and (3) Proof of Eq. ( 6).-Consider a system S interacting with a bath B. The total density matrix of the composite object is ρ total , whereas the reduced one for the system and the bath are ρ and ρ B , respectively.Suppose that the total initial state is a tensor product of the form: Where ρ B (0) is assumed to be fixed at: with p β ≥ 0, for every initial state ρ(0) of the system.
At time T , the state of the composite object is simply a rotation of the initial state (it is a closed system): ρ total (T ) = U(T )ρ total (0)U + (T ).
Here, U(T ) = T (e −i ´T 0 H total (t ′ )dt ′ ) is the propagator for the entire object, where T is the timeordering operator, and H total is given by: where H S , H B , H SB are terms in the Hamiltonian that depend only on S, on B, or on degrees of freedom of both, respectively.Taking the trace of Eq. ( 57) with respect to the states of B yields ρ(T ): where: is a Kraus operator and Eq. ( 59) is often known as the operator sum representation [35,36].By identifying: we have shown the equivalence between Eq. ( 59) and (6).
Proof of Eq. ( 3).-As explained in Section II and in the previous proof, Eq. ( 6) holds for any proper or improper state ρ(0) in the Liouville space of the system.With this in mind, Eq.
(3) is straightforward to prove.We want to enforce that Eq. ( 6) preserves trace throughout the Proof of Eq. ( 4).-Again, we exploit the fact that Eq. ( 6) is true for any ρ(0).

)
Denoting |g as the molecular ground state or the excitonic vaccuum, |A = a + A |g and |B = a + β |g are the excitons at each site, whereas |α = a + α |g , |β = a + β |g are the delocalized excitons.The biexcitonic state, expressed by |f (b)).The two states at the coherence time which can oscillate at ω gα are |g α| or |β f |, but the latter cannot be produced by a single action of the dipole operator on the initial ground state |g g|.Therefore, |g α| is the only possible state for the coherence interval, and is produced by acting the first pulse on the bra of the ground state: |g g| → |g α|.Similar considerations imply that the state at the echo time must be |β g| or |f α|.Given these constraints, we are ready to enumerate the possible initial and final states for the waiting time interval which are compatible with these restrictions.By exciting the ket or deexciting the bra of |g α| with the second pulse, the following initial states |c d| for the quantum channel can be produced: |c d| ∈ {|α α|, |β α|, |g g|}.The final states |a b| ∈ {|α α|, |β α|, |β β|, |g g|}

Figure 1 :
Figure 1: Liouville space pathways corresponding to each of the four resonances in the 2D-ES of a coupleddimer.The amplitude S mn (T ) corresponds to the peak located at (ω τ , ω t ) = (ω mg , ω ng ), which are the values of the optical frequencies at the coherence and echo time intervals τ and t, respectively.These amplitudes provide information on the coherent and incoherent excitonic processes at the waiting time T by enumerating the possible initial and final states at the waiting time T which satisfy the PE phase matching condition for the pulses acting at times t 1 , t 2 , t 3 , t 4 .The information contained in the amplitudes S mn (T ) can be distilled to reconstruct the process matrix for the single-exciton manifold of the dimer, thus allowing for a QPT.
can all give rise to |β g| or |f α| by exciting the ket or deexciting the bra with the third pulse.

Figure 2 :
Figure 2: A more detailed view on the Liouville space pathways corresponding to the monitoring of the population to coherence process χ βααα (T ) in the peak at (ω τ , ω t ) = (ω αg , ω βg ) of the 2D-ES.These diagrams belong to the amplitude S αβ (T ) and can be easy constructed by taking into account the PE phase-matching and the resonant conditions.

Figure 3 :
Figure 3: Summary of QPT for a coupled dimer in the 2D-ES.The Liouville pathways depicted in Fig. 1 can be condensed into this diagram.The horizontal axis for the coherence frequency ω τ is associated with a state preparation, whereas the vertical axis for the echo frequency ω t corresponds to a detection.The four resonances labeled as (m, n) correspond to peaks located at (ω τ , ω t ) = (ω mg , ω ng ).Their amplitudes contain information on χ abcd (T ), where cd is the state prepared at the beginning of the waiting time interval, and ab the state detected at the end of the same interval.For instance, the peak at (ω α , ω β ) keeps track of the elements χ abcd (T ) where |c d| ∈ {|α α|, |β α|, |g g|} and |a b| ∈ {|β α|, |β β|, |α α|, |g g|}.

Figure 4 :
Figure 4: 2D-ES of a coupled dimer in the absence of interactions with a bath.Under unitary dynamics of the excitonic system, each of the four resonances keep track of the elements of χ(T ) indicated in the diagram.Notice that diagonal peaks do not oscillate as a function of waiting time T , whereas off-diagonals beat at the frequency equal to the difference in energies of the single exciton eigenstates.

Figure 5 :
Figure 5: Transition dipole moments of a homodimer.Diagrams for (a) sites and (b) eigenstates.

Figure 6 :
Figure 6: Stability of the homodimer QPT protocol.As described in the text, the QPT protocol depends on the inversion of a matrix which is a function of transition dipole moments.The condition number of the matrix (κ) vs. the angle between site dipoles (φ) is plotted in red.The blue line plots the constant value of κ = 15, as a reference to indicate that for the range of angles 0.3π ≤ φ ≤ 0.7π, the condition number is below that value.
of F W HM = 20 f s, i.e., σ = 8.49 f s, which amount to an equal excitation amplitude C for both |α and |β .Rabi oscillations for a coherent superposition between |α and |β occur with a period T c = 47.5 f s.We present several snapshots of the real and imaginary parts of the spectra at values of waiting time T corresponding to multiples of T c /2, skipping T c = 0, as our theory has avoided pulse overlap effects.In principle, as Eqs.(23) and (25) indicate, the excitonic quantum beats associated with the term χ αβαβ (T ) can be monitored by looking at either cross peak of the spectra.This feature is subtly manifested in every column of the figure, but more easily perceived in the real part the zzxx spectrum, where the peak at (ω β , ω α ) changes from red to yellow/green every interval T c /2 before the bosonic bath has washed out significant portion of the coherent dynamics at about T = 5T c .Also, incoherent population transfer primarily from |β to |α (downhill) manifests as a decrease in amplitude of the peak at (ω β , ω β ) due to ESA and an increase in (ω β , ω α ) due toSE [2].This effect can also be more obviously seen in the real part of zzxx spectrum.

Figure 7 :
Figure 7: 2D-ES for coupled porphyrin homodimer with secular Redfield model.From left to right, we show the real and imaginary parts of the spectrum with zzzz polarization setting (first and second columns), and with zzxx setting (third and fourth columns).Each row represents a particular waiting time T , corresponding to (a) T c /2, (b) T c , (c) 3T c /2, (d) 2T c , (e) 9T c /2, (f) 5T c , where T c = 47.5 f s is the period for one Rabi oscillation between |α and |β .The colormap is such that red is associated with positive numbers, green with values about zero, and blue with negative numbers.

First, we consider
a population as the initial state, i.e. ρ(0) = |k k| for any k.Note that ρ(0) is a proper density matrix, and hence Hermitian.Then, ρ ab (T ) = cd χ abcd (T )δ ck δ dk = χ abkk (T ), 2 /(2σ 2 ) .The pulses are sent to the sample in a non-collinear fashion to the sample, generating a time-dependent dipole in each of the molecules.Since the characteristic size of a molecule is much smaller than the wavelength of the radiation, 2π/|k i |, each molecule only experiences a potential that changes in time but is uniform in space, in consistency with the I e 1 e 2 e 3 e 4 ;m 1 m 2 m 3 m 4 e 1 e 2 δ e 3 e 4 δ e 1 e 3 δ e 2 e 4 δ e 1 e 4 δ e 2 e 3 ]

TABLE 3 .
Values (in f s −1 ) of non-zero rates of the secular Redfield tensor