Perturbative theoretical model of electronic transient circular dichroism spectroscopy of molecular aggregates

A chiral analog of transient absorption spectroscopy, transient circular dichroism (TCD) spectroscopy is an emerging time-resolved method. Both spectroscopic methods can probe the electronic transitions of a sample, and TCD is additionally sensitive to the dynamic aspects of chirality, such as those induced by molecular excitons. Here, we develop a theoretical description of TCD for electronic multi-level models in which the pump pulse is linearly polarized and probe pulse is alternately left-and right-circularly polarized. We derive effective response functions analogous to those often used to describe other four-wave mixing methods and then simulate and analyze TCD spectra for three representative multi-level electronic model systems. We elaborate on the presence and detection of the spectral signatures of electronic coherences


I. INTRODUCTION
Transient circular dichroism (TCD) spectroscopy is an emerging time-resolved spectroscopic technique-analogous to transient absorption (TA) spectroscopy-that aims to resolve macromolecular structural evolution and the dynamics of molecular excitons, among other processes.Multiple groups have developed fast and ultrafast TCD spectrometers for performing laboratory measurements of the dynamics of photodissociation of molecular enantiomers, 1 photo-initiated molecular ring-opening reactions, 2 peptide and protein structural changes, [3][4][5][6] excited states in molecules, [7][8][9] electronic and structural relaxation of molecular macromolecules and nanostructures, [10][11][12] energy transfer in photosynthetic proteins, 13 exciton states in inorganic halide perovskites, 14 and excitons and phonons in polymers. 15,16Interestingly, the analysis and interpretation of the measured spectra in these studies used either a phenomenological treatment of the measured signal or no underlying model at all, even though other groups have developed and presented microscopic, quantum-mechanical theoretical descriptions of TCD and related time-resolved chiral spectroscopic techniques.For example, in a 2003 pioneering report, Cho adapted the responsefunction formalism to study chiral nonlinear-optical signals and focused on the case of a circularly polarized pump pulse and linearly polarized probe. 17Abramavicius and Mukamel used the nonlinear exciton equations method to calculate the second-order and third-order nonlinear-optical responses of anharmonic oscillators to model peptides. 18,19More recently, Holdaway et al. used a theoretical doorway-window approach to study exciton-coherence signals generated in simulations of TCD measurements. 20Very few TCD studies have analyzed and interpreted measured spectra using an underlying microscopic theory. 21ne possible reason for the separation between theory and measurement is the lack of simulations based on model systems.Therefore, here we present a theoretical description of thirdorder TCD spectroscopy for the common laboratory conditions of a linearly polarized pump and circularly polarized probe using the adapted response-function approach. 22We use the theoretical method to simulate TCD spectra for a variety of multi-level The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp3][24] In addition, a multi-level electronic model is an ideal starting point to describe a molecular aggregate, where electronic coupling between or among chromophores leads to delocalized excitation states commonly known as excitons. 25,268][29] A steady-state CD spectrum is an incisive probe of the chirality induced by electronic coupling in these types of molecular aggregate samples, and our theory and simulations reveal the information contained in TCD spectra can yield insights about the dynamics induced by these couplings.Developing a microscopic theoretical description of the TCD spectra of molecular aggregates is challenging for at least four reasons.First, CD signals do not fall within the ubiquitous electric dipole approximation; the magnetic dipole and electric quadrupole interactions must be included. 30These interactions pose a challenge for nonlinear-optical spectroscopy methods, where the electric dipole approximation is invoked in nearly all instances. 22,31As a simplification, we do assume that the intrinsic magnetic dipole and the intrinsic electric quadrupole of each molecule are zero.The second complication is that CD techniques require polarized optical signals, and therefore the theoretical analysis must track the full vector representing the optical polarization of each beam.In contrast, one can often obtain significant insight into TA spectroscopy signals even when the vectors are replaced with scalars.Tracking the full vectors rather than scalars approximately triples the number of terms.The third complication is that TCD is a differencespectroscopy technique, where signals arise due to incomplete cancellation between two polarization cases.This approximately doubles the number of terms required to model a given signal compared to conventional TA spectroscopy.The final complication is that the signals require a potentially complicated basis transformation between the molecular sites and the delocalized electronic states of the aggregate sample.Accomplishing this transformation correctly required algebraic manipulations using the strict rules of tensor-sum methodology.In Sec.II, we satisfy these complications to develop the general expressions by adapting the familiar response-function formalism.The general expressions generate some physical insights.Then, in Sec.III, we generate further insights using simulations of representative model multi-level electronic systems.

II. THEORETICAL RESULTS AND DISCUSSION
In this section, we first derive the source term for the signal field tracking both the macroscopic polarization and magnetization fields.We then develop the perturbation theory of the density operator for the electromagnetic perturbation.Using these results, we then work out the first-order and third-order responses both generally and specifically for molecular aggregates.

A. Source term of signal field
The interaction between the electromagnetic field and the sample requires us to account for the radiation caused by the macroscopic magnetization dipole (M) and the symmetric rank-2 tensor for the macroscopic polarization quadrupole (Q ), in addition to the conventional polarization dipole (P) of the material.(To avoid confusion with multiple meanings of variables i and k, we reserve the variables {a, b, c} as generic indices of the three spatial coordinates, {x, y, z}.)When all of these terms are included in Maxwell's equations, the wave equation for the ath component of the electric field becomes [32][33][34] where c is the speed of light, we have assumed a transverse electric field (E) with ∇ ⋅ E = 0, and the notation x 1,2,3 in the partial derivatives indices the set {x, y, z}.We choose to work with the traceless quadrupole moment, 33,34 ∑ a Q aa = 0. We can write the polarization dipole and quadrupole and magnetization dipole, respectively, as where Ps(t), Q s (t), and Ms(t) represent the time-dependent envelope functions and where the frequency ωs and wavevector k ′ s are determined by sums of plus or minus the frequencies and wavevectors of the incident fields that created the polarization and magnetization, These are the phase-matching relationships that ensure both energy and momentum are conserved.The detected signal field takes the form where the signal wavevector (ks) is not necessarily equal to the polarization and magnetization wavevector.We insert these forms of the polarization, magnetization, and signal fields into the wave equation, Eq. ( 1), and then we assume that the propagation is entirely in the z-direction through a thin, weakly absorbing sample of length l.Finally, we invoke the slowly varying envelope approximation to yield The Journal of Chemical Physics

ARTICLE scitation.org/journal/jcp
where n is the index of refraction of the signal in the medium and Δk = |k ′ s − ks| is the phase-matching condition.We can then expand the standard relationship E sig ∝ iP to include the polarization quadrupole and magnetization dipole, where ω and k are the signal frequency and wavevector, respectively.CD spectroscopy relies on the subtraction of two intensitylevel measurements, and therefore it is important to clarify that for heterodyne detection, one interferes the signal field, E sig , with a relatively intense reference field, E ref .
Here, we assume that the reference and signal fields have the same handedness, such that As in conventional TA spectroscopy, we will subtract a control measurement of I ref to cancel the |E ref | 2 term, and we will assume that |E sig | 2 is weak.Using the form of the signal field from Eq. ( 6), the measurement becomes where we have used some vector identities and the relationship between the electric and magnetic fields in an electromagnetic wave in the last step.

B. Perturbation theory of density operator
Having identified the contributions to the source term, we must describe their microscopic origins in the molecular response.To begin, we assume the total Hamiltonian is partitioned into a timeindependent component, Ĥ0 , and an electromagnetic perturbation, V(t), given by where μ, q , and m are the electric dipole, electric quadrupole, and magnetic dipole operators, respectively, and E and B are the electric and magnetic fields carried by a laser beam.The quadrupole operator can take one of several forms; hence, we chose to define the matrix elements of q for a molecule consisting of N charges with charge q (p) located at positions We assume ρ 0 = ρ(−∞) is given by a thermal distribution of the eigenstates of Ĥ0 .For electronic spectroscopy, it is generally accurate to assume that all molecules are in the ground electronic states, ρ 0 = ρ gg .Following convention, 22 we expand the density operator in orders of the perturbation such that ρ ≡ ρ (0) + ρ (1) + ρ (2) + ⋅ ⋅ ⋅ , (11)   where ρ (n) represents the nth-order contribution.In the interaction picture, we define the density operator and the perturbation, respectively, as where Û0 is the time-evolution operator given by Û0 (t, t 0 ) = e −i Ĥ 0 (t−t 0 )/ ̵ h .Each order of the density operator is given by ρ (n) where the τn are the interaction time intervals given by τn = t n+1 − tn.This density operator can be used to calculate the signal field using the perturbative contributions to the polarization dipole, P(t), polarization quadrupole matrix elements, Q ab (t), and magnetization dipole, M(t), which can be written as In Sections II C and II D, we use the above expressions to address the first-order and third-order responses for molecular aggregate systems.Pedagogically, the first-order case is a useful midpoint that assists in developing and comprehending the third-order case, where the number of terms is larger.

C. First-order response
To calculate the first-order response of a molecular aggregate, we need the first-order density operator, Similar to the conventional analysis, 22 we insert the perturbation expression and find that the polarization dipole elements are given by The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp the polarization quadrupole matrix elements are given by and the magnetization dipole vector is given by We have suppressed the interaction-picture notation subscripts, and in Eqs. ( 17) and ( 18),we have ignored any term having more than one combined interaction with the electric quadrupole or magnetic dipole because they will produce a signal that is much weaker than the retained terms.
We next assess the polarization components of the fields.We use Jones vectors to describe the polarization states of the fields.
Although not necessary at first order, we elaborate on the field components and the phase-matching conditions because this level of detail will be required at third order.Each semi-impulsive field can be decomposed into its time envelope, polarization, wavevector, and frequency components as where ϵ L/R,± are unit vectors that identify the positive and negativemomentum electric-field polarization for left-(L) and right-(R) circularly polarized fields and where the parameter T accounts for the relative delay between laser pulses necessary in the third-order calculations.In the first-order calculations, we are free to choose T = 0. We refer to each component by the sign of the momentum contribution. 23The phase-matching condition for absorption is k (1) sig = k 0 , which means that we must select the same momentum components for the input and signal fields.The complex conjugate symbol ( * ) indicates to use the negative-wavevector component of the indicated field.While the field amplitudes, E 0 , are real-valued and therefore the conjugate symbol seems unnecessary, we found it helpful to retain the symbol as a bookkeeping reminder of the sign of the wavevector component and also acknowledge that there could be a phase difference.
We assume the fields propagate in the z-direction so that the polarization vectors are given by where ϵ x/y are unit polarization vectors in the x and y direction, respectively.The resulting magnetic fields (B = 1 ω k × E) are given by with unit vectors where k = |k|ϵz.
The time-dependent polarization dipole, polarization quadrupole, and magnetization dipole elements due to a left-or right-circularly polarized incident field can then be written as, respectively, and where the spatial dependence is omitted.A steady-state CD measurement is given by the difference of intensity-level measurements having left-(L) and right-(R) circularly polarized fields, measured in the frequency domain.Therefore, the signal is given by where E L,R 0 is the incident light in each measurement and we have applied a Fourier transformation from time variable t to detection frequency variable ω to the input and signal fields as The Journal of Chemical Physics ARTICLE scitation.org/journal/jcpand We follow convention and assume that all fields are spectrally white and can therefore be represented by δ-functions in the time domain. 22,23,31his yields The first-order CD signal for a fixed molecular orientation can then be written as In the case of a sample composed of a large number of arbitrarily oriented chromophores, the signal must be averaged over all molecular orientations.It is convenient to do the averaging in a coordinate system fixed in the molecular frame and averaging over all orientations of the lab frame.Using the orbit symbol ( ) to indicate a coordinate-independent orientational average, we need to compute (31)   where, for example, ϵ x a represents the projection of the unit vector ϵ x in the lab frame onto the ath axis of the molecular frame and we have assumed that the molecular frame is fixed so that we can remove the operators from the orientational average.The general results for the orientational averages of these tensors are established. 33For the second-rank tensor, 33 the general result is where δ is the Kronecker delta and i and j are any of {x, y, z} in the lab frame.For the third rank tensor, 33 the general result is where {i, j, k} are cyclic permutations of {x, y, z} and for {i, j} ∈ {x, y, z}.This yields (32)   where we have assumed the vector operators can be separated into their constituent orientational and operator components: The Journal of Chemical Physics ARTICLE scitation.org/journal/jcpμ(0) = μ(0)n(0) or m(0) = m(0)η(0), where ns and ηs correspond to electric dipole or magnetic dipole orientations, respectively.Electric quadrupole contributions average to zero because the quadrupole moment is symmetric, giving ∑ bc ε abc qbc = 0.
There is insight to be gleaned from an analysis of Eq. ( 32).First, unlike the analogous steady-state absorption signal, 22 which contains a single commutator involving two electric dipole interactions, the steady-state CD signal contains two commutators that each involve one magnetic dipole interaction and one electric dipole interaction.Second, each commutator can be cast into the responsefunction formalism. 22For example, we can label the terms arising from the expansion of the commutators as where we use the † symbol to denote raising operators.We incorporate the dot product of the dipole orientations into the response functions and use the subscript "CD" and the orientational-average notation to distinguish these effective response functions from the more traditional tensor response functions.Although we do not do so, identifying the response functions would allow one to draw double-sided Feynman diagrams akin to those commonly used in nonlinear-optical spectroscopy by making one adjustment, which is to track each interaction as a magnetic or an electric dipole.Holdaway et al. present such diagrams. 20o make a strong connection to literature 22,31 and to set up the labeling scheme required at third order, we write the orientationally averaged first-order CD signal as (34)  where, using H.c. to indicate Hermitian conjugate, we can define In this work, we consider a molecular aggregate sample composed of multiple molecules each having at least two electronic states, and we evaluate the two first-order response functions individually.Working in the delocalized excitation basis {|G⟩, |E 1 ⟩, . . ., |En⟩}, we evaluate the trace, set ρ 0 = |G⟩⟨G|, and insert a complete set of states to yield (36)  Next, we use the result of Appendix A to convert the magnetic dipole operator and orientation vector into its equivalent electric dipole operator and orientation vector.This yields (37)   where we have used the fact that μ † = μ because the electric dipole operator is real-valued and where we have used the vector identity a ⋅ (b × c) = b ⋅ (c × a).We then make the time dependence explicit and use an identify derived in Appendix B to transform the matrix elements and the orientation vectors from the delocalized excitation basis to the site basis.This yields (38)   where Ω a,b = ω a,b − iγ a,b is the complex frequency containing phenomenological dephasing γ a,b for the optical coherence between delocalized excitation states |Ea⟩ and |E b ⟩.Here, Roman characters index the delocalized states and Greek characters index the sites.The constraint α ≠ β on the sums over these two indices is a reminder that when α = β, the cross product of orientation vectors will be zero.
A similar analysis for the second response function yields (39)  where the only notable distinction is that there were two sign changes that resulted in no net change to the sign.
Subtracting the two response functions and simplifying yields (40)  where R αβ = Rα − R β .We can then insert this result into Eq.( 34) to derive the general linear CD result for molecular aggregates, (41)   This expression yields multiple physical insights.(1) The CD spectrum will have a Lorentzian peak at each delocalized excitation frequency, ωj,G.(2) The sign and amplitude of each peak is given by the coefficients of the linear combination of sites that give rise to the specific delocalized excitation state, the c j α/β terms.(3) Also affecting the sign and amplitude of the peaks are the relative distance and orientation of the dipoles R αβ ⋅ (n β × nα).These are independent of the specific delocalized excitation state.(4) The matrix elements indicate an excitation event on site α and a de-excitation event on site β, which can occur only when the two sites are coupled.The coupling is contained in the coefficients c j α/β .

D. Third-order response
Following the methodology of the first-order case, we next introduce the appropriate perturbative contribution to the density The Journal of Chemical Physics ARTICLE scitation.org/journal/jcpoperator and field polarization and magnetization contributions to the signal under the appropriate laboratory conditions.There are three significant modifications from the first-order signal calculation.First, the third-order signal is subject in the phase-matching condition of the pump-probe geometry, k sig = k pump − k pump + k 0 , where the probe field is indicated by k 0 because it also serves as the reference.A second key difference from the first-order case is that the pump fields are not necessarily circularly polarized.Prior works have shown that any linear polarized pump pulse is adequate for a time-resolved third-order CD measurement that uses circularly polarized probe pulses. 35Third, the electric quadrupole terms do not average to zero in general for isotropic samples in thirdorder experiments as they do in first-order experiments.The electric quadrupole contribution is dependent on the relative angle between the pump and probe beams; Cho showed 17 that it is suppressed for a crossing angle of θ ∠ = arctan( 1 √ 2 ) ≈ 35.26 ○ .We assume a circularly polarized probe pulse propagating along the z axis and a horizontally polarized pump pulse that crosses the z axis at angle θ ∠ in the (x, z)-plane and assume again the semi-impulse limit, and where kpu = k(sin θ ∠ ϵx + cos θ ∠ ϵz), the H polarization unit vectors are given by ϵ H = cos θ ∠ ϵx − sin θ ∠ ϵz, and  H = k k × ϵ H = ϵy.We choose to have the probe propagate along the z axis with kpr = kprϵz so that the polarization vectors for the pump match the prior definitions.We assume that the pump pulse envelope peaks at t = −T and the probe at t = 0.

General third-order response
We consider the general third-order response in which we include electric dipoles and quadrupoles and magnetic dipoles.In this case, we need the third-order density operator, ρ (3) We insert the perturbation and find that the polarization and magnetization are, respectively, and and the elements of the polarization quadrupole are The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp These are the nonlinear versions of Eqs. ( 16)- (18), respectively.Each term can arise from an interaction with the electric dipole, electric quadrupole, or magnetic dipole of the molecule.The terms can be grouped according to the field interactions.Terms having zero or one interactions with either the magnetic dipole or electric quadrupole will dominate the response, and therefore we neglect terms having two or more magnetic dipole or electric quadrupole interactions from the remainder of the work.The terms that survive this approximation are + P (3) (3) where the subscripts indicate from right to left the time-ordered sequence of interaction.For example, P In general, all of these terms can contribute.However, after subtracting the right-circularly polarized signal from the left-circularly polarized signal, P μμμ (t) will cancel for an isotropic sample.Furthermore, the specific choice of geometry described above 17 eliminates contributions from P ab,μμμ (t) after orientational averaging.This leaves six contributions to the signal, each with eight response functions.After evaluating the τ integrals for the semiimpulsive fields, the time-dependent polarization and magnetization terms can be written as and The Journal of Chemical Physics

ARTICLE scitation.org/journal/jcp
where the spatial dependence is omitted.These expressions assume that-under the rotating-wave approximation-the pump will contribute one positive-momentum interaction and one negativemomentum interaction.The sign(k 1(2) ) tracks the sign of the momentum of the quadrupole moment interaction with the field.These expressions also assume a positive-momentum interaction with the pump field and a negative-momentum reference field.
Analogous to the first-order CD signal, the third-order CD signal is given in the frequency domain by where E L,R * pr (ω) = E * pr ϵ L/R,− as in the first-order case, and we can write the ath component of the third-order signal as For bookkeeping, we sort the contributions to the total CD spectrum by the order of the interactions with a fourth label to identify the interaction leading to the emission of the signal field, where the subscripts again identify from right to left the sequence of interactions.For the third-order calculation, none of these contributions can be neglected in general.The choice to work at the suppression angle does suppress the last two in an isotropic medium, 17 .We consider the terms individually, first expanding into tensorsum notation to ease the orientational averaging.For example, where a, b, c, d again reference indices corresponding to the {x, y, z} components of the vectors.However, the summed results are valid in any Cartesian coordinate system, and therefore we may choose those to represent a coordinate system {x ′ , y ′ , z ′ } fixed in the molecular frame and average over all orientations of the {x, y, z} lab frame.This is convenient because the laser polarization and propagation vectors are perpendicular to each other.The orientational averaged contribution of this term to the CD signal in an isotropic sample becomes (53)   where we used the fourth-rank tensor orientational average, 33,36 (54)

ARTICLE scitation.org/journal/jcp
The Kronecker δ-functions arising from the polarizations create dot products between the dipole operators, which are straightforward to evaluate when considering individual response functions.
The remaining terms work out similarly, assuming the pump and probe cross at the suppression angle, where cos θ ∠ = √ 2 3 .Terms involving the quadrupole moment require the fifth-rank tensor isotropic average. 33The final results for an isotropic sample assuming the pump and probe cross at θ ∠ are (55a) Each of the six orientationally averaged TCD signal contributions in Eq. (55) has 6 unique response functions (plus 6 Hermitian conjugate response functions) for a total of 36 unique response functions.We do expand the nested commutators and define tensor response functions for each of the 36 contributions to the total signal.We define R 1 → R 6 to consistently represent the same type of pathway labeled SE (stimulated emission), GSB (ground state bleach), and ESA (excited-state absorption) and Reph (rephasing) or NR (nonrephasing) next to each response function. 31The total signal will include the Hermitian conjugate terms that appear at negative signal frequencies.Appendix C contains the full list of 36 response functions.These response functions naturally separate into nine groups of four.For each of SE, GSB, and ESA type, there are three unique contributions: One has a magnetic dipole interaction with the pump pulse, one has a magnetic dipole interaction during the probe pulse and signal emission, and one has an electric quadrupole interaction with the pump pulse.In general, it would also be necessary to consider electric quadrupole interactions during the probe, but they are suppressed by the choice of geometry.
Each response function can be evaluated in a manner analogous to the first-order case, albeit with more terms.For example, consider the rephasing stimulated-emission term which yields a contribution to the total CD signal of (57) where we use μ † α = ⟨eα|μ † α |gα⟩ to represent the transition dipole moments.

ARTICLE scitation.org/journal/jcp
We assume the coefficients and electric dipole matrix elements are real and incorporate the model from Appendix A to write the signal contribution expression in terms of the magnitude of the matrix elements and a sum over the dipole orientations, (58) where we have introduced a phenomenological population decay rate (Γ) during the pump-probe delay time.

Summary of terms
After many mathematical manipulations, we are left with only five nonzero expressions.There are two for stimulated emission, These five expressions are all that is required to simulate a TCD spectrum for a multi-level electronic system, and they are the primary result of this work.

ARTICLE scitation.org/journal/jcp
While deriving these expressions, several key steps generated physical insights.(1) For terms in which the pump pulse contributed the magnetic dipole moment or electric quadrupole moment interaction, the rephasing and nonrephasing contributions have a relative sign change.For the magnetic dipole, this happens because the swap from rephasing to nonrephasing flips m ↔ m † .For the electric quadrupole interaction, this happens because the sign of the momentum of the electric field involved in the interaction flips.This cancels all terms that leave the density matrix in a population after the pump pulse interactions.The implication is that terms that include a coherence (|E j ′ ⟩⟨E j|) during the delay time do not cancel, and therefore these terms are a key to revealing the contribution of coherences in TCD signals.Furthermore, our model assumes that each site is a two-level electronic system with no vibrational modes, and therefore the coherences are purely electronic, cannot be created on the ground state, and the GSB contributions with a magnetic dipole or electric quadrupole interaction during the pump pulse will sum to zero, leaving seven nonzero terms.(2) In stark contrast, terms in which the probe pulse contributed the magnetic dipole interaction yield both population and coherence terms.(3) We reduce the initial seven nonzero terms to five by summing SE and ESA contributions that involve a magnetic dipole or electric quadrupole interaction during the pump pulse.Separately, these terms retain a dependence on the choice of origin in the molecular frame, but the summed contribution depends only on relative positions of the molecular sites in an aggregate.(4) The ESA terms required an additional transformation between matrix elements in the delocalized exciton basis and the site basis that involve transitions to the doubly excited delocalized exciton states, |F k ⟩.We derive this relationship in Appendix B. This means that the ESA signals-if they can be unambiguously distinguished in the TCD spectrum-could be exquisitely sensitive to the properties of the doubly excited delocalized exciton states.

III. SIMULATION RESULTS AND DISCUSSION
Although the general equations derived in Sec.II provided some meaningful physical insights, simulations of multi-level electronic systems can generate more intuition regarding CD and TCD signals arising from molecular aggregates.In this section, we detail features of spectra simulated for a homodimer, homotrimer, and homotetramer.We use Eqs.( 41) and ( 59)-( 63) to simulate steadystate CD and TCD spectra, respectively.To produce the maximum number of peaks having nonzero amplitude, we evaluate a consistent set of arbitrarily located and oriented vectors.Table I lists the dipole coordinates and Fig. 1 displays their geometric arrangement in a Cartesian coordinate system.Studies of homodimers have generated a significant amount of intuition about spectroscopic signals arising from molecular aggregates. 28In this model, the transition dipoles of a pair of two-level systems are coupled, and this coupling yields two delocalized excitation states, having energies above and below the site energies.The Hamiltonian is given by where we have partitioned the diagonal elements, Eι, from the off-diagonal elements, J ι,κ , and where sums over ι and κ run over the set of {α, β, γ, δ} as appropriate for each aggregate.The diagonal elements are the site energies and the off-diagonal elements are the site-site couplings.We compute the site-site coupling using the expression where which is appropriate in the dipole-dipole approximation. 26,28,37We invoke this approximation in our simulations, and, further, we compute the two-exciton couplings in a similar manner.For example, the coupling between two-exciton states |α, β⟩ and |α, γ⟩ is J βγ .
In the simulations, each site has a frequency of 500 THz and a transition-dipole magnitude of 1 a.u.microscopic model of electronic-vibrational interactions. 22,26While one might expect that the coherence dephasing time should be the same as the linewidth, 38 here we found that doubling the dephasing time to 250 fs was helpful to make each simulated coherence signal distinguishable.Finally, in all simulations we scale each coupling value by an arbitrary factor of 40 to resolve the transitions more clearly.From Table I, the dimer simulation uses sites α and β, the trimer simulation includes also site γ, and the tetramer simulation uses all four sites.In the dimer simulation, the relative angle between dipoles is about 65 ○ .All simulations used a frequency step of 0.5 THz and a time step of 2.5 fs, which required about one second of computation time on a personal computer using scripts written in Python and appropriate libraries.The simulated spectra in Fig. 2 contain several notable features.First, the steady-state CD and TCD spectra show the same number of main peaks for each system: The dimer simulations have two peaks, the trimer simulations have three main peaks, and the tetramer simulations have four main peaks.The number of one-exciton states in both simulations is N, the number of sites, and the number of two-exciton states, N 2x , used in the TCD simulations is given by Although the number of main peaks is conserved between the steady-state CD and TCD spectra for a given system, their relative signs are not.This is not unexpected because comparing Eq. ( 41) to Eqs. ( 59)-( 63) reveals a number of mathematical components of the signal expressions that could lead to a sign change.For example, the terms that account for the relative orientation of the dipoles are distinct.
Second, in the dimer simulation, the rightmost panel shows that the peaks oscillate at about four cycles in 0.5 ps; in other words, the oscillation frequency appears to be about 8 THz.This frequency corresponds to the difference between the exciton eigenfrequencies, 496 and 504 THz.Due to interference effects between overlapping peaks of opposing signs, the peak maxima are shifted slightly to 493 and 505 THz when the linewidth is broad.The distinction between exciton eigenfrequencies and the peak maxima is observable in the steady-state CD spectrum of the dimer, where the peak maxima in the broadband simulation are further apart than the peak maxima in the narrowband simulation.At first glance, the quantum beats arising from the exciton-exciton coherences at the two peaks appear to have a π phase shift.However, the rightmost panels in Fig. 2 are normalized, which means the negative-amplitude (blue) trace was divided by a factor of −1.Hence, our TCD simulations are consistent with transient absorption spectroscopy, in which quantum beats arising from exciton-exciton coherences will oscillate with the same phase. 23,39Holdaway et al. examined several relative phase differences for coherent oscillations using a distinct theoretical treatment for a vibronic dimer system. 20These phase shifts may help distinguish electronic coherences from vibrational or vibronic coherences in future studies.To make an additional connection to literature, the TCD measurement by Trifonov et al. conducted on merocyanine dye helical aggregates showed dimerlike steady-steady CD and TCD spectra.However, those authors did not display the early time, <1-ps, dynamics, and hence we cannot compare to any coherent oscillations in that spectra.The trimer and tetramer TCD simulations in Fig. 2 contain analogous, albeit more complicated, TCD spectra and coherent dynamics.In these systems, the number of transitions and number of possible exciton-exciton coherences increase significantly.These coherences generate a linear combination of oscillation frequencies that lead to highly structured time-domain dynamics.
Third, the effects of the two-exciton states appear to be weak relative to the transitions involving only one-exciton states.In fact, there are no signatures of two-excitons in the dimer simulations because there is only a single two-exciton state and-because we neglect biexciton interaction contributions 26 -its frequency is exactly double that of the site energy.This means there is perfect destructive interference between the possible pathways that involve the two-exciton state.In the trimer and tetramer simulations, however, the simulations reveal small shoulders and extra peaks that arise from transitions involving the two-exciton states.The peak at about 450 THz in the tetramer TCD simulation is one of the more visible peaks involving a two-exciton transition.The lowestfrequency one-exciton peak is at 471 THz in this simulation.The peak at 450 THz therefore must involve a two-exciton state, and we observe that it also contains oscillating quantum beats due to exciton-exciton coherences.
One advantage of a theoretical model that is unavailable to laboratory measurements is the ability to decompose the total signal into its constituent components for further study.In the left panels of Fig. 3, we present the five components of the total TCD signal for the dimer case, wherein each panel is individually normalized so that weak components are visible.The color bar bounds indicate the relative normalization values.The SE1, SE2, GSB, ESA1, and ESA2 contributions arise from Eqs. ( 59)-(63), respectively.For the dimer, the SE1, GSB, and ESA1 terms contain no coherent oscillations.Instead, these pathways are dominated by signals arising from population decays.These three pathways all have the same overall amplitude, and the largely two-peak structure of the total signal shown in Fig. 2 mainly arises from these three pathways.In contrast, the SE2 and ESA2 terms contain no population decays and only have coherent oscillations.These coherences have an intriguing phase profile, and there is a relative phase shift of π between the two pathways.These two pathways are weaker, at ∼16% and 63% amplitude relative to the three other pathways.A primary conclusion from these separated components of the signal is that-in stark contrast to transient absorption spectroscopy where there is a difference in relative sign-there is no clear way to distinguish excited-state absorption pathways from bleach or stimulated-emission pathways in TCD spectroscopy measurements of dimers.Finally, we examine the separate components of the TCD signal for the case of the trimer.The right panels of Fig. 3 show these simulated signals.Like the dimer case, the GSB signal remains strong (normalized amplitude of 1), and its dynamics are dominated by population decays.The SE1 and ESA1 signals also contain population decay pathways but also now contain coherences, and their relative amplitudes are only 26% and 20% of the GSB amplitude, respectively.This is distinct from the dimer case wherein these two signals had no coherences and had the same amplitude as the GSB signal.The coherences contained in the SE1 and ESA1 signals are relatively uncomplicated; the one-exciton eigenfrequencies are 480, 504, and 516 THz, meaning the highest possible frequency is about 36 THz.By contrast, the SE2 and ESA2 signals contain coherences that appear to be of higher frequencies.However, the two-exciton eigenfrequencies are 980, 1004, and 1016 THz, and hence the highest possible frequency is also 36 THz.Indeed, the complicated peak pattern arises from interference effects due to sign changes among the coherences.A Fourier transform analysis (not shown) supported this conclusion.Like the dimer case, the SE1 and ESA2 signals at 13% and 29%, respectively, are weaker than the other signals.In the case of the trimer, transitions involving the two-exciton states are no longer automatically suppressed due to destructive interference as in the dimer case.Instead, the ESA1 and ESA2 components contain signals that extend to nearly 550 THz, arising from the difference between the two-exciton state having a frequency of 1016 THz and the one-exciton state having a frequency of 480 THz.These signals do not appear in the SE1, SE2, and GSB components.Hence, for trimers and higher-order aggregates, it may be possible to distinguish excited-state absorption signals in TCD spectra by identifying peaks that appear at frequencies higher or lower than any peak present in the steady-state CD or absorption spectra.

IV. CONCLUSIONS
Transient circular dichroism spectroscopy is a chiral analog to the ubiquitous method of transient absorption spectroscopy.This nascent method is slowly developing on both the measurement and theoretical fronts.Here, we have developed effective response functions for steady-state CD and TCD signals.This required including an electromagnetic perturbation that incorporated the magnetic dipole and the electrical quadrupole moments, applying orientational averaging, and accounting for the full vector fields rather than scalars.We used the effective response functions to simulate spectra for an arbitrary electronic dimer, trimer, and tetramer.These are multi-level electronic models of molecular aggregate samples.The simulations revealed several key spectral signatures regarding the peak structure and their dynamics for analysis of coherences and excited-state absorption components.
We anticipate that-because achiral molecules produce no CD or TCD signals-TCD will be a more incisive probe of the excited-state dynamics of molecular aggregates, which can be mixtures of chiral and achiral species. 40This method may be useful for interpreting measurements of natural molecular aggregates such as biological pigment-protein complexes, and better resolve the excitonic structure of such systems.Future work could transcend the static-dipole approximation and lead to interesting studies of the dynamics of excimers. 41The theoretical methods presented here could be expanded and adapted to simulate chiral coherent multidimensional spectra as well.Another interesting avenue to explore is adding line shape theory to include the effects on the electronic spectra of intramolecular and intermolecular vibrational modes. 42,43This would be an important advance because the pigment-pigment and pigment-vibration interactions are of approximately equal strength in some natural biological molecular aggregates.

ARTICLE
scitation.org/journal/jcpdistribution along the length of the rod-is bound at the ends and has zero width in the other two dimensions.The spatial distribution of the current density has the form where J 0 (y − y 0 ) is any distribution that is bound at ±d/2 and zero at the end points, and where we denote r = xî + y + z k.The charge density and current density can be related through the continuity equation The electric dipole moment is given by μ = ∫ rρ(r)d 3 r.For this example, the x and z components integrate to zero and the y component is given by where the factor of i indicates a phase shift.The magnetic dipole moment is given by Finally, comparing Eqs.(A4) and (A5) and defining μ = μ yields The relationship derived in this last step arose from an extrapolation of the example for the conditions indicated in Fig. 4.However, for a single dipole, we should have the freedom to place and orient the dipole in any location relative to the origin.Therefore, Eq. (A6) should be general for all orientations and locations, and indeed it is equivalent to previous results. 44e can use the classical result in Eq. (A6) to generate the corresponding dipole operators, using the orientation vectors η and n defined in Sec.II.

Electric quadrupole
In this work, electric quadrupole interactions with the pump pulse can contribute to the TCD signal, and therefore we use a model to relate the electric quadrupole operator to the electric dipole operator and a displacement.This model is appropriate for molecules with zero net charge and a nonzero electric dipole moment.Similar to the magnetic dipole moment, the electric quadrupole moment is dependent on the choice of origin.We begin with the definition of the electric quadrupole moment in two coordinate systems where the origin of the r ′ coordinate system is shifted by R 0 from the origin of the r coordinate system such that r (p) ′ = r (p) − R 0 as shown in Fig. 5, giving where we have assumed that the net charge on each molecule is zero and used the definition of the electric dipole moment.For rod-like molecules with a nonzero electric dipole, a center of dipole exists where the quadrupole moment is zero. 33We choose the origin of the r ′ coordinate system to sit at that this point, so that the quadrupole moment can be found in the r coordinate system in terms of the dipole operator and the displacement,

Coordinate system independence
While the electric quadrupole moment and magnetic dipole moment are both dependent on the choice of coordinate system, the radiated signal field must be independent of the arbitrary choice of origin in the molecular frame.To this order in the multipole expansion, the Hamiltonian describing the laser-molecule interaction includes both the electric quadrupole and magnetic dipole interactions.When both contributions are summed, the resulting Hamiltonian is independent of the choice of coordinate system up to a term that results from the phase of the laser field.In the calculated signals, this contributes a dependence only on relative positions of the molecular sites in an aggregate.Therefore, all final measurable contributions are independent of the choice of molecular origin.We confirmed this by performing additional simulations to those shown in Sec.III with a common origin offset and found that the signals were identical regardless of the choice of origin.

APPENDIX B: SITE AND DELOCALIZED EXCITATION BASIS TRANSFORMATIONS
In this work, we treat each molecule-known as a site-as a two-level electronic system including its ground state, |g⟩, and The Journal of Chemical Physics ARTICLE scitation.org/journal/jcpits excited state, |e⟩.The electric transition-dipole operator for site ι is given by μι = μι|e⟩⟨g| + μι|g⟩⟨e|, and the total electric transitiondipole operator in the site basis is given by the tensor sum 45,46 of the electric transition-dipole operator for each site, Because this notation may be unfamiliar, we explicitly write the three-site system operator, The tensor-sum notation allows us to construct a Hilbert space that includes all possible permutations of the states of the sites.
where we used lowercase Greek symbols to index sites 1 through η.The interpretation of the single-excitation expression is that the jth delocalized excitation state is a linear combination of the site-basis states for which only the γth site is excited, weighted by the coefficients c j γ .The inequality in the last expression accounts for the fact that any individual two-level system cannot be doubly excited and that the permutation of the two excitations is irrelevant.
Given these definitions, the matrix elements in the delocalized excitation basis can be written as in the site basis.To arrive at the last line, we used ⟨gi|μi|gi⟩ = 0.The tensor-sum notation can accommodate coupling among an arbitrary number of sites having arbitrary numbers of excited states.Analysis of excited-state absorption pathways involve matrix elements such as ⟨F k |μ|Ej⟩, which can be written in the site basis as We next identify that only certain combinations of bras, operators, and kets will be nonzero.The bra will always have two sites in the excited state, whereas the ket will only have one site in the excited state.The index of the site that remains in the excited state has complete freedom.The index of the site that transitioned from the excited state to the ground state needs to match the index of the transition-dipole operator, thereby linking the ζ and γ indices as well as the ι and ϵ indices, or vice versa.Hence, instead of 4 free indices, there are really only 2. Importantly, because we assume the sites are two-level systems, γ ≠ ϵ.Simplification yields where, in the second step, to aid the eye, we removed all notation for the sites that remained in the ground state, and then in the third step, we suppressed the notation for the site that remained in the excited state.Because of the permutation symmetry of the doubly excited states, ci,j = cj,i, and hence the expression can be reduced to In this appendix, we list all of the third-order response functions for TCD.The first set involves those arising from an interaction between the probe field and the magnetic dipole, R The second set contains terms that involve an interaction between the pump field and the magnetic dipole, R  The third set contains terms arising from an interaction between the pump field and the electric quadrupole, R ρ

FIG. 2 .
FIG. 2.Simulation results for dimer (top row), trimer (middle row), and tetramer (bottom row).Steady-state CD spectra are shown in the left column.Dashed lines indicate simulations using γ = 0.2 THz to highlight the exciton eigenfrequencies.The middle column contains the full TCD simulated spectra, which are individually normalized.The rightmost column contains normalized dynamics of the brightest positive (red) and negative (blue) amplitude peaks.

FIG. 3 .
FIG. 3. Components of the dimer (left) and trimer (right) TCD signal.The signals have relative maximum absolute value amplitudes as indicated by the color bars.

FIG. 5 .
FIG. 5.Illustration of field lines from (left) an electric (μ) dipole and (right) an origin independent electrical quadrupole (q ab ) displaced by R 0 arbitrarily along the x axis.While generally the fields are distinct, the two moments produce similar fields at the origin.

TABLE I .
Dipole coordinates and orientations used in simulations.
Table II lists the resulting exciton and two-exciton frequencies for each aggregate.In all steadystate CD simulations, we set |k||E 0 | 2

TABLE II .
Frequencies of exciton and two-exciton states in each aggregate.
This notation-where state vectors are an ordered direct-product basis indicating the state of each site, for example |e 1 ⟩ ⊗ |g 2 ⟩ = |e; g⟩-becomes unwieldy for larger aggregates.Therefore, we will also work in the delocalized excitation basis, often known as the exciton basis.The delocalized basis state vectors are {|G⟩, |E 1 ⟩, . . ., |En⟩, |F 1 ⟩, . . ., |F k ⟩}.The electric dipole operator, the ground state, and the delocalized excitation states can be written in the site basis, respectively, as