The excitation of molecular systems by natural incoherent light relevant, for example, to photosynthetic light-harvesting is examined. We show that the result of linear excitation with natural incoherent light can be obtained using incident light described in terms of transform limited pulses, as opposed to conventional classical representations with explicit random character. The derived expressions allow for computations to be done directly for any thermal light spectrum using a simple wave function formalism and provide a route to the experimental determination of natural incoherent excitation using pulsed laser techniques. Pulses associated with solar and cosmic microwave background radiation are provided as examples.

“

One of the most important attributes of a stationary random process is its spectrum.”L. Mandel and E. Wolf

^{1}

## I. INTRODUCTION

Light incident on matter can initiate fundamental dynamics in systems as diverse as semiconductors, quantum dots, and living organisms. Of particular interest is the interaction of sunlight with molecules such as photosynthetic light-harvesting complexes,^{2} the foremost step responsible for primary energy production on earth, or sunlight interacting with visual pigments, initiating the first steps in vision.^{3}

While the development of advanced experiments such as non-linear spectroscopy^{4} allows characterizing a system’s dynamical response upon photoexcitation, the molecular dynamics measured in this way differs significantly from that resulting from natural excitation conditions,^{5–10} a consequence of the differences in the character of the incident radiation.^{5,11} Although it is possible to measure molecular dynamics directly from incoherent excitation,^{12} technical difficulties are such that, in practice, experimental data are extremely limited, and the spectrum is not necessarily that of natural light (i.e., sunlight and moonlight). At present, the dynamics of a system excited by thermal light can be best inferred from appropriate theoretical models,^{13–15} which require a proper description of thermal excitation. Two different approaches are available:^{16} (i) a full quantum treatment, using a quantized radiation field (e.g., the *P*-representation, with the light represented by a direct product of incoherent mixtures of coherent states of a single frequency—typically plane waves), or (ii) a semi-classical approach based on a classical description of the radiation field. It was recently shown that the quantum state of light, in first order, can be composed of a mixture of single pulses.^{17} We focus here on the semi-classical formalism of thermal excitation and provide a description of the associated molecular response.

As a proper mixture, thermal light can, in principle, only be represented statistically, i.e., by a probability distribution function or, alternatively, as a set of all possible realizations. The most conventional classical representation, the collision-broadening model,^{16} relies on realizations with stochastic phase jumps to reproduce the random character of thermal light, with the average time interval between jumps reproducing the (short) coherence time. Producing such a set of realizations in the laboratory, so as to reproduce the spectrum of natural light, is difficult. Here we show, in the relevant domain of *linear* excitation, i.e., weak excitation appropriate to natural incoherent light, that an ensemble of well defined *deterministic* realizations can reproduce the result of thermal-light excitation. Indeed, even realizations comprised of transform limited pulses are possible, forging a link with modern pulsed-laser laboratory techniques. The result provides a route to exploring incoherent excitation with a wide variety of possible spectral profiles.

This issue was motivated by ongoing interest in excitation of, for example, photosynthetic pigment-protein complexes such as the Fenna-Matthew-Olson (FMO) complex^{18} or PC645,^{19} or the LH2 light-harvesting complex.^{20,21} Since our focus is on the molecular response to the light, the relevance and interpretation of the results are independent of the effects of the system’s environment. Thus, we deal with a closed molecular system. Further, due the weakness of the interaction of sunlight with such systems, first order perturbation theory provides an essentially exact treatment.

## II. FORMALISM

Consider the Hamiltonian

where the system Hamiltonian (first term on the r.h.s.) is written on the basis of molecular eigenstates |{*α*}〉_{} with eigenenergies ε_{α}, and the second term describes the interaction with the field in the semi-classical approach. Here, $d= \u2211 \alpha \mu \alpha |\alpha \u3009 \u3008g | n \u02c6 \alpha +h.c.$, with |*g*〉_{} the molecular ground state, *μ*_{α} is the transition dipole moment operator, oriented along the unit vector $ n \u02c6 \alpha $, and **E**^{(γ)}(*t*) denotes the positive frequency part of one individual field realization. Specifically, the superscript *γ* is used throughout the paper to label a particular realization, either of the field or of the corresponding molecular excited state arising from **E**^{(γ)}(*t*) incident on the molecule. Hence, this Hamiltonian characterizes the molecular dynamics upon interaction with one of the field realizations.

In the rotating wave approximation and in the limit of weak-field interaction, first-order perturbation theory provides a solution for the excited state wave function resulting from the |*g*〉_{} state excitation as a superposition of eigenstates,

where *N*_{α} is the number of molecular eigenstates and *ω*_{αg} ≡ (*ϵ*_{α} − *ϵ _{g}*)/

*ħ*denotes the eigenfrequency,

*ϵ*being the ground state energy. Here, we assume the excited-state level splitting to be larger than their spontaneous decay widths and neglect effects from radiative decay.

_{g}^{14}For the sake of clarity, we consider a scalar field and assume it lies along the orientation of the transition dipole moment $ n \u02c6 \alpha $, taken identical for all eigenstates.

Considering that the molecule’s response to a stochastic field is comprised of an ensemble of individual excitation events, the principle quantity of interest is the excited state density matrix,

where $ \u2026 \Gamma $ denotes averaging over the ensemble of realizations. From Eq. (2), the elements of this density matrix are

with *ω*_{αβ} = *ω*_{α} − *ω*_{β}. Here and below, we follow the literature (e.g., Ref. 22 and references therein) in assuming that the system encounters the radiation field suddenly at *t* = 0. Work on the general effect of slow turn-on on realistic systems is in progress.^{23,24}

Requiring that this density matrix be representative of the result of thermal-light excitation imposes restrictions on the fields *E*^{(γ)}(*t*) and ensemble $ \u2026 \Gamma $. For the system considered in Eq. (1), interaction with thermal light would lead to the excited state described by a density matrix *ρ*^{tot} with elements

where the first-order correlation function between the *i* and *j* Cartesian components of the field, for thermal light, is^{25}

with *n*(*ω*) = (*e*^{ħω/(kBT)} − 1)^{−1} being the average photon number at temperature *T* and *k _{B}* the Boltzmann’s constant.

Requiring that Eq. (4) describes the same molecular dynamics as that resulting from linear thermal-light excitation then implies that

To calculate the left-hand-side, i.e., ensemble-averaged auto-correlation function of the field, we use the Fourier representation for the individual realizations given by

where $ E \u0303 ( \gamma ) ( \omega ) $ denotes the Fourier components. We choose all individual fields with identical spectral distribution $| E \u0303 ( \gamma ) ( \omega ) |$ and only allow differences in the phase factor. This is a sufficient but not necessary condition, and we could well reproduce the excitation from different lineshapes.^{17} The Fourier coefficients $ E \u0303 ( \gamma ) ( \omega ) $ can therefore be written as the product of a real function *f*(*ω*) defining the spectral distribution (independent of *γ*) and a phase factor *ϕ*^{(γ)}(*ω*) (dependent on *γ*),

The field auto-correlation function averaged over the ensemble of realizations is then

By choosing the phase factor with a linear dependence on frequency and interpreting the ensemble averaging as an integral over the phases, i.e.,

where *t _{s}* is a scaling constant with unit of time that can be chosen arbitrarily without affecting the total molecular response. Note that this is only one of many possible representations; any collection of fields that would recover a delta function

*δ*(

*ω*−

*ω*′) in the last line of Eq. (10), thus reducing the two frequency integrals to a single one, would suffice. Requiring that Eq. (10) matches the thermal-light first-order correlation function [Eq. (6)] then restricts the lineshape

*f*(

*ω*) so that the electric field for each individual realization Eq. (8) becomes

Using this field as the individual field realizations, the individual excitation events [Eq. (2)] become

where the complex constant *K* and coefficients *C*_{α}(*t*, *γ*) are provided in the Appendix.

The state |*ψ*^{(γ)}(*t*)〉_{}〈*ψ*^{(γ)}(*t*)|_{} obtained from Eq. (13) represents one contribution to the total density matrix following excitation from a single realization of the field. The total density matrix *ρ*^{tot} is then obtained by averaging over the ensemble according to Eq. (11b).

It is easy to analytically show that, with this definition of the field realizations in Eq. (12) and using Eq. (11b), the density matrix *ρ*^{tot} built from the ensemble of individual excitation events recovers that obtained directly from thermal light excitation (5), i.e.,

Note that the formulation outlined above is applicable to alternate incoherent spectral densities with differing *n*(*ω*), e.g., that associated with a wide range of temperatures (e.g., sunlight at 5777 K to cosmic microwave background (CMB) at 2.7 K^{22}), perturbed underwater sunlight^{26} or excitation from deep-sea hydrothermal vents.^{27} Below, we focus on excitation from sunlight at T = 5777 K and comment on the CMB case.

## III. COMPUTATIONAL RESULTS

We compute the molecular response from the total density operator, constructed as average over the ensemble of individual responses (15). Here and below, the Hamiltonian is given by Eq. (1) with parameters that are representative of a photosynthetic light-harvesting complex: *ϵ*_{α} = 1 eV, *ϵ*_{β} = 1.1 eV, *ϵ _{g}* = 0,

*μ*

_{α}=

*μ*

_{β}= 1 e.Å. A sample comparison of these results, obtained via realizations and those obtained from direct numerical integration of (5) with the use of thermal correlation function (6) is shown in Fig. 1, where perfect agreement is seen. Convergence to

*ρ*

^{tot}(

*t*) is obtained using only 20

*t*realizations, where

*t*is in units of fs. As in Refs. 22, 28, and 29, where the incoherent radiation is modeled as white-noise, we find that sudden turn-on incoherent excitation leads to coherent excitation, with the amplitude of the coherence inversely proportional to the energy splitting between eigenstates. In this approach, the coherence survives destructive interference between the individual excitation events with primary contribution from realizations with 0 ≤

*γ*≤

*t*. Although non-zero, the coherences are of constant amplitude and quickly, i.e., after a few femtoseconds, become negligible compared to the linearly growing populations.

Of particular interest is the nature of the realizations of the field. As expected from the linear phase relationship between the different modes imposed by (11a), field (8) takes the form of a widely used class of pulses; see, e.g., Ref. 30. This also appears clearly in Fig. 2, showing that the linear dependence of the phase on frequency localizes the individual realizations, forming pulses. In addition, the *γω* phase factor acts such that the field realizations only differ from one another by a time translation, evident from Eq. (12) and the pictorial representation in Fig. 2. Thus, the averaging process over all realizations Eq. (11b) becomes equivalent to averaging over all times, which recovers the classical time average. Further, the spectral density of each realization Eq. (12), and hence the coherence time (here ≈1.3 fs for 5777 K sunlight), is identical to that of Planck’s radiation law because of the matched correlation functions [Eq. (7)] and the fact that realizations differ only through the phase factor *ϕ*^{(γ)}(*ω*).

Figure 3 shows several individual excitation events corresponding to the particular field realizations sampled in Fig. 2. Each event corresponds to a response from a pulsed driving force, i.e., with coherent evolution and constant population after excitation.^{13} Although the field realizations are all identical to within a time factor, the individual molecular responses are more complex. Specifically, the field realizations differ by a phase factor, reflected in a time shift of the individual diagonal elements contributing to the total populations $ \rho \alpha \alpha tot ( t ) $. On the other hand, the individual molecular coherence contributions $ \rho \alpha \beta tot ( t ) $ present, in addition to the time translation, variations in amplitude. Consequently, although the individual molecular responses are not phase related (because of different onsets), for sudden turn-on the total contribution can still exhibit some non-zero oscillatory coherences that survive destructive interferences.

The above results are representative of pulses associated with thermal excitation with sunlight. However, a wide variety of cases may be considered, from ultra-hot thermal sources, to excitation with thermalheat vents, to studies of the excitation from the CMB, all represented merely by differing photon number distribution *n*(*ω*), and hence different spectra and associated pulses. For example, studies of atoms and molecules in CMB radiation are difficult since the radiation is extremely weak, (5777/2.7)^{4} = 2.1 × 10^{13} weaker than sunlight. However, the linear nature of the dynamics allows us to study the light-molecule interaction in the laboratory by scaling up the field strength. Current CMB amplifiers^{22,31} do not permit such large scale amplification, but the approach advocated here does, and amplifying the pulses intensity by 10^{A} would correspondingly amplify the population by the same factor, allowing one to reach measurable quantities (cf. Ref. 22 for the population of state 65s of Rb, for a possible application). For example, for 90-dB amplification, such pulse would exhibit a peak power of 1/2 *cϵ*_{0}|*E*^{(γ)}(*γ*)|^{2} ∼ 0.6 W.cm^{−2} per surface area. The power per surface area of the ensemble of pulses would then be $ E ( \gamma ) ( t ) \u2217 E ( \gamma ) ( t ) \Gamma \u223c3.3kW. m \u2212 2 $, which represents, as expected, 10^{9} times that at the surface of a black-body radiation at CMB temperature, given by Stefan-Boltzmann law. Specifically, a useful laboratory pulse for CMB radiation is shown in Fig. 4. Utilizing this pulse for, e.g., excitation of Rydberg atoms would allow for studies of, for example, deviations from Markovian behavior.^{14}

## IV. DISCUSSION AND SUMMARY

In general, ensembles can be represented by a wide variety of realizations, the so-called issue of “unraveling” the ensemble.^{32} The particular approach advocated here has yielded, for the case of thermal light, a representation of the light consisting of a set of identical deterministic transform-limited pulses, translated in time from one another. This result has formal, computational, and practical (experimental) advantages. Formally, it demonstrates that each individual realization of the field can be treated deterministically, thus eliminating any randomness in each realization. Beyond the fundamental implications, this result also simplifies computational realizations. Experimentally, using incident transform limited pulses allows for a direct connection to well-established pulsed laser techniques. Interestingly, each successive realization utilizes the same pulse, shifted in time, considerably simplifying preparation of the incident pulses.

Carrying out the proposed procedure experimentally requires irradiation of different molecules, each with a different *γ* value, corresponding, as shown above, to irradiation at different starting times. The challenge is that different molecules do not possess a “clock” that makes them aware of these different starting times. Meeting this requirement is, however, relatively straightforward. That is, only a small fraction of the molecular sample is affected by any given pulse. Therefore, a second pulse is unlikely to excite the same molecule that was affected by the first pulse. If one now excites a sample at consecutive times *and collects the entire signal resulting from all excitation events with a single detector*, the detector, in viewing all the excited molecules collectively, will automatically carry out the required averaging over realizations and do so with an awareness that the molecules have been excited at different times.

It is noteworthy that in the particular approach adopted here, the ensemble average can be performed in many ways. That is, similar results, fulfilling Eq. (14), can be obtained with individual realizations that differ from Eq. (12), e.g., non-pulsed light could be used instead. In particular, the only requirement to achieve the thermal result is to match the thermal-light first-order correlation function, which requires (i) matching the spectral distribution (and hence the coherence time), and (ii) that the ensemble of field realizations fulfills $ e i \varphi ( \gamma ) ( \omega \u2032 ) \u2212 \varphi ( \gamma ) ( \omega ) \Gamma \u221d\delta ( \omega \u2212 \omega \u2032 ) $. Satisfying these conditions will reproduce the stochastic character of the first-order correlation of thermal-light. The latter requirement resembles the condition that the Fourier components of a stationary random process belonging to different frequencies are uncorrelated.^{1} While the Fourier components of the derived fields in Eq. (12) are correlated at the level of the individual realizations, this correlation is lost upon averaging over a collection of independent realizations. This is in accordance with the fact that, for a stationary and ergodic ensemble of random functions, one can replace the auto-correlation function defined by the time average by that defined as an ensemble average.^{33} Our results demonstrate that this not only holds at the level of the field but also applies to describe linear excitation of molecules.

Three additional comments are in order. First, note specifically that here *coherent* light is used to generate the *incoherent* result. This is distinctly different from previous work, wherein *spectrally incoherent* (i.e., non-transform limited) light is used to study *coherent* excitation (see Appendix 10B of Ref. 34 and associated references therein). Second, our results apply to *linear* interaction, appropriate for natural incoherent light. In general, however, one cannot represent the state of thermal light (i.e., the full density matrix) with an ensemble of single broadband coherent states (representative of pulses), because such a representation already fails for the second-order correlation describing broadband delayed excitation.^{11} Third, we emphasize that the individual realizations have no physical meaning with respect to excitation with incoherent light. Rather, it is only the ensemble average that describes the response to incoherent excitation.

In summary, the implications of these results are two-fold for the examination of the interaction of matter with natural incoherent light. Conceptually, linear excitation by thermal light can be represented by an ensemble average of deterministic pulse realizations—all effects of incoherence are recast in the post-excitation averaging. Practically, Eq. (13) provides a simplified description of the molecular excitation that allows studies of natural incoherent excitation using molecules excited by coherent transform limited pulses at the wave function level; explicit and complicated random fields (e.g., Refs. 10 and 28) need not be used. In addition, these results motivate both experimental and computational studies of natural incoherent excitation using pulsed laser techniques.

## Acknowledgments

We thank Professors Jennifer Ogilvie, Amar Vutha, and Junrong Zheng for discussions on the experimental protocol and Dr. T. V. Tscherbul for comments on the manuscript. A.C. acknowledges funding from the Swiss National Science Foundation. P.B.’s research was supported by U.S. Air Force Office of Scientific Research under Contract No. FA 9550-13-1-0005 and the Natural Sciences and Engineering Research Council of Canada.

### APPENDIX: DERIVATION AND DETAILS OF THE FORMALISM

In this appendix, we provide details of the wave function dynamics after excitation by one realization of the field, Eq. (2) in the main text.

Taking, as exciting field, the individual realizations derived in (12) yields

Here, we have assumed the field to lie along the orientation of the transition dipole moment $ n \u02c6 \alpha $, and we have defined

Expanding $ n ( \omega ) \u2261 e \u2212 \u0127 \omega 2 k B T ( 1 \u2212 e \u2212 \u0127 \omega k B T ) \u2212 1 / 2 $ in terms of power of exponentials yields to a series of terms for which the integral can be evaluated, i.e.,

Note that $ p \nu ( t ) = ( \u2212 1 / i ) \nu \u2202 \nu \u2202 t \nu p 0 ( t ) $. We therefore take the time derivative to obtain

We add the integration in time and define

with $F ( x ) \u2261 e \u2212 x 2 \u222b 0 x e y 2 dy$ is the Dawson function and

Although Eq. (A2) involves an infinite number of terms, the sum converges quickly. For example, at room temperature, results are typically within 2% accuracy with the first three terms only.

Taking

the dynamics of system reduces to the expression given in the main text [Eq. (13)].