Laser Control of Dissipative Two-Exciton Dynamics in Molecular Aggregates

There are two types of two-photon transitions in molecular aggregates, that is, non-local excitations of two monomers and local double excitations to some higher excited intra-monomer electronic state. As a consequence of the inter-monomer Coulomb interaction these different excitation states are coupled to each other. Higher excited intra-monomer states are rather short-lived due to efficient internal conversion of electronic into vibrational energy. Combining both processes leads to the annihilation of an electronic excitation state, which is a major loss channel for establishing high excitation densities in molecular aggregates. Applying theoretical pulse optimization techniques to a Frenkel exciton model it is shown that the dynamics of two-exciton states in linear aggregates (dimer to tetramer) can be influenced by ultrafast shaped laser pulses. In particular, it is studied to what extent the decay of the two-exciton population by inter-band transitions can be transiently suppressed. Intra-band dynamics is described by a dissipative hierarchy equation approach, which takes into account strong exciton-vibrational coupling in the non-Markovian regime.


Introduction
Molecular aggregates continue to provide inspiration and challenges to experiment and theory [1,2,3]. In terms of a microscopic understanding of the dynamics of elementary excitations recent advances due to time-resolved nonlinear spectroscopy have been tremendous [4,5]. However, despite the success of laser pulse control in rather different areas of research [6], relatively little emphasis has been put on exciton dynamics. On the experimental side, feedback laser control has been applied to manipulate the branching ratio between internal conversion (IC) and energy flow in the light-harvesting antenna of purple bacteria [7]. There has been a number of simulations for light-harvesting systems by May and Brüggemann et al. where the focus was put on the transient generation of a single localized excitation within the aggregate [8,9,10,11]; non-biological aggregates have received no attention so far.
The quantum dynamics of excitons is well founded within the Frenkel exciton approach, which starts from a classification of aggregate's excitation states in terms of the number of simultaneously excited monomers [12,13]. Under weak irradiation conditions only a single excitation will be present, but in nonlinear optical experiments or in the presence of strong irradiation multiple excited states play an important role. The organic building blocks of molecular aggregates have a multitude of electronically excited states, i.e. there is always a state S n such that the S 0 -S 1 excitation energy approximately matches that of a S 1 -S n transition. As a consequence one needs to distinguish between local double excitations (LDE) and nonlocal double excitations (NDE) where the two excitations are localized at different monomers. Note that the NDE should not be confused with bi-exciton states, i.e. bound states formed by two excitons in the presence of different permanent dipoles in the involved electronic states [13]. LDE and NDE can couple via the Coulomb interaction. The manifestation of this coupling in nonlinear spectroscopy has been investigated in Refs. [14,15,16]. The presence of LDEs has important consequences for the exciton dynamics since it leads to exciton-exciton annihilation (EEA) by virtue of an intramolecular IC triggered by nonadiabatic electronic transitions. The presence of this process in molecular aggregates is well established [17,18,19] and various theoretical approaches exist [20,21,22,23,24], although it must be stated that a first principles determination of the respective IC rates is still out of reach.
Understanding exciton dynamics in aggregates is impossible without taking into account the effect of exciton-vibrational coupling [3,12,25]. Here one can distinguish between approaches which either account for all vibrations on the same footing, i.e. in terms of a heat bath, or treat a few active vibrational degrees of freedom explicitly, but coupled to the heat bath of the remaining modes. Needless to say, the latter approach is more demanding as the dimension of the density matrix increases rapidly with the number of explicit modes. Early investigations along these lines therefore have been restricted to molecular dimer models with one active vibrational coordinate per monomer [26,27,28,29]. More recent approaches include a Green's function description of exciton dynamics [30], a multi-configuration time-dependent Hartree simulation [31] or the Frenkel excitonic polaron treatment [32]. These methods cannot treat finite temperature effects due to the coupling to a further heat bath in a microscopic manner.
The coupling of the electronic degrees of freedom to a single or site-specific heat bath is usually treated using dissipation theory [3]. Recently, the so-called hierarchy equation (HE) method [33,34,35,36,37,38,39] has enjoyed great popularity as it promises a nonperturbative and non-Markovian description of exciton dynamics, which is numerically equivalent, e.g. to the influence functional approach within a path integral formulation [40]. In the present context one should note that a distinction between different types of modes, i.e. low-frequency solvent modes versus high-frequency intramolecular vibrations, can be introduced via the spectral density. This is of particular relevance since the chromophores used to assemble artificial molecular aggregates usually show vibronic side bands in their absorption spectra pointing to the prominent role of intramolecular high-frequency vibrations [41].
The focus of the present contribution is on the laser control of dissipative twoexciton dynamics in linear aggregates. Thereby we take into account the effect of EEA and of a coupling to a heat bath being composed of an effective high-frequency mode as well as of a continuous distribution of low-frequency modes. For the solution of the density matrix equations of motion we will use a HE approach. Two-exciton populations are usually quenched by EEA. Since the latter is a local process we will ask the question whether an excitonic wave packet can be prepared such as to transiently suppress EEA by its composition in terms of LDE states. In the following Section 2 we will outline the density matrix approach and its interface to optimal control theory (OCT). Section 3 will start with a discussion of the field free dynamics focussing on the comparison between the Markovian and non-Markovian limits. Subsequently, laser driven dynamics is discussed. The results are summarized in Sec. 4.

Frenkel Exciton Hamiltonian
Consider an aggregate composed of N monomers, each carrying three adiabatic electronic states (a = g, e, f ) corresponding to the S 0 , S 1 , and some S n state. Thus we have the adiabatic electronic basis of local states |m, a , m = 1, . . . , N . The electronic states of the aggregate can be classified as the zero excitation (ground) state [3] the singly-excited state |m = |m, e Π n =m |n, g , the doubly-excited (LDE and NDE) states |mn = |m, e |n, e Π k =m,n |k, g .
Restricting ourselves to these excitation the completeness relation reads The electronic Hamiltonian can be written aŝ with the bare exciton Hamiltonian, being composed of the ground state part where E 0 is the electronic ground state energy, the single-exciton Hamiltonian and the two-exciton Hamiltonian Here, E m,e and E m,f are the energies of electronic excitation at site m for S 0 -S 1 and S 0 -S n , respectively. Further, J mn is the coupling between site m and site n leading to single exciton transfer and J (ef) mn is the Coulomb coupling responsible for the two-exciton dynamics.
It is customary to introduce Frenkel exciton eigenstates which follow from the diagonalization of the exciton Hamiltonian The eigenstates can be decomposed in terms of local excitation states, |a (a = 0, m, mm, mn), as follows The eigenstates separate into M-exciton manifolds (here M = 0, 1, 2). Frequently, we will also use a notation where the states in the M-exciton manifolds are counted according to increasing energy, i.e. |M k . In Eq. (6) the aggregate is assumed to interact with an external laser field treated in dipole approximation, i.e.
where the field E(t) is oriented parallel to the transition dipoles. In the local Frenkel exciton basis the dipole operator for the aggregate readŝ µ = m µ m,e |m 0| + µ m,f |mm m| + n =m µ n,e |mn m| + h.c.. (14) Here, µ m,e (µ m,f ) is the transition dipole for S 0 -S 1 (S 1 -S n ) transition of monomer m. Finally, we consider the IC process described byĤ IC in Eq. (6). IC has its origin in the breakdown of the Born-Oppenheimer approximation, e.g. at conical intersections. The respective non-adiabaticity operator triggering transitions between adiabatic electronic states is proportional to the momentum operatorP m,ξ of the involved nuclear degrees of freedom (coordinatesQ m,ξ ) counted by the index ξ at monomer m. The non-radiative life time of the S 1 state is usually in the nanosecond range. Therefore we will restrict ourselves to the IC between adiabatic states |m, f and |m, e . Hence the coupling becomeŝ m,ξP m,ξ , In the following we will assume that the coupling strength is site-independent, i.e. g (IC) m,ξ = g (IC) ξ

Exciton-Vibrational Coupling: System-Bath Model
Exciton-vibrational coupling leading to intra-band phase and energy relaxation is introduced in the spirit of the system-bath approach, i.e. we havê The bath modes are treated in harmonic approximation, i.e.Ĥ vib is the standard harmonic oscillator Hamiltonian. Concerning the system-bath coupling we will distinguish between two types of modes (see also Ref. [42]). First, local high-frequency modes, {q m,ξ }, which usually give rise to a vibrational progression in the monomer absorption spectrum [41]. Second, global low-frequency modes, {x ξ }, contributing via a continuous spectrum. Hence we havê For the coupling to the local high-frequency modes we will use the model Hamiltonian m,ξ,f |mm mm| m,ξ,e |mn mn| .
In the following we will assume that the coupling is the same for all monomers, i.e. g For the coupling to the low-frequency modes we will invoke the same approximations and arrive at the Hamiltonian 2.3. Dissipation Models 2.3.1. Internal Conversion The IC rate between S n and S 1 states is rather large for typical chromophores, reaching time scales on the order of about 100 fs [44]. On this time scale the S n population decays into the S 1 state, thereby passing many electronic states (for the case of perylene bisimides n would be of the order of about 30 [41]). Hence it is justified to assume that the memory time associated with the IC process is even shorter than 100 fs and EEA can be treated in Markovian approximation. Treatinĝ H IC in Eq. (15) as a perturbation of the bare exciton Hamiltonian within second order perturbation theory one can write the contribution to the Quantum Master Equation for the reduced exciton density operatorρ in terms of the Lie operator R (IC) which operates onρ as [3] i HereΞ m is defined aŝ and α m,ξP m,ξ of Eq. (15) (note that we assume that the momenta at different sites are uncorrelated). By using the completeness relation for the eigenstates, one haŝ Here, J

Separation of Time Scales in the Response Function of the Oscillator Bath
In general the influence of the bath degrees of freedom due to exciton-vibrational coupling cannot be treated in Markovian approximation. It is fully characterized by the response function, α(t), of the bath, which in turn is determined by the spectral density function J(ω) Here, g ξ , µ ξ , and ω ξ are the coupling constant, reduced mass, and frequency, respectively, of the ξth oscillator and β = 1/k B T . In the following we will specify the spectral density to the coupling models discussed in Section 2.2. First, we consider the case of the local high-frequency mode, Eq. (18) for which we assume a damped Brownian oscillator model to hold. This can be described by the following spectral density [45] Below we will assume the underdamped limit, where the central mode frequency ω 0 is much larger than the cutoff Λ H . For this spectral density the response function can be expanded into a sum of exponentials as follows where is the kth Matsubara frequency of the bath, and N H is the smallest integer satisfying ν N H +1 Λ H . As discussed in Ref. [33] the response function can be split into two parts: a long memory contribution (the first two terms in the last line of Eq. (27)) and a short memory part α (short) H (t). The latter one contains all the terms with short memory in the Matsubara summation and will be treated in Markovian approximation within the hierarchy equations to be discussed below.
The low-frequency bath will be described by the Debye spectral density where η L is the coupling strength and Λ L is the frequency cutoff of the bath. The response function in Eq. (25) then can be expressed as follows where ν 0 = Λ L , ν k = 2πk/( β)(k > 0) is the kth Matsubara frequency of the bath, and N L is the smallest integer satisfying ν N L +1 Λ L . The same time scale separation has been introduced as for the high-frequency bath.

Hierarchy Equations for the Dissipative Multi-exciton Dynamics
In order to describe the non-Markovian and non-perturbative exciton dynamics we will utilize a HE approach to propagate the reduced exciton density matrix in eigenstate representation. Specifically we have employed the stochastic decoupling procedure due to Shao and coworkers which is briefly sketched in the Appendix [37,38]. Since the HE approach starts from the influence functional [33,46], it is based on the fact that the response function for bath is written as a sum of exponentials (see previous section).
The resulting HEs of motion for the dissipative exciton dynamics of the aggregate are given by where A is an N × 2 integer matrix, with its matrix element A mj corresponding the contribution from Ω j in m-th molecule (Eq. (27)). B is an N × N H integer matrix with its matrix element B mj counting the j-th Matsubara frequency in the m-th molecule (Eq. (27)). And V is an N L + 1 dimensional integer vector. V 0 refers to the cutoff frequency and V j (j = 0) to the j-th Matsubara frequency in the low frequency bath (Eq. (29)). The operator " + / − " is defined as M ± mj nl = M nl ± δ mn δ jl for the matrix and V ± j k = V k ± δ kj for the vector. R is the Redfield super-operator accounting for the exciton-exciton annihilation and the short memory effects of the high and low frequency environments, [·, ·] denotes the commutator and {·, ·} denotes the anti-commutator. For the definition of α k,± and c (H) 0 , see Appendix.
In the above equations, the first termρ 000 is the reduced density matrix for the exciton and the other elements reflect the finite memory effect due to the excitonvibration coupling. The initial condition for the hierarchy isρ 000 (0) =ρ(0) and ρ ABV (0) = 0 if any of the elements of A, B and V is nonzero.
In the Markovian limit, the HEs are reduced to the Quantum Master Equation where R is the Redfield superoperator Here R (IC) is the superoperator accounting for the IC as defined in Eq. (21).Ξ (tot) andΞ m are the dissipative operators defined in the similar way as Eq. (22) for the low frequency and high frequency bath in m-th molecule, respectively.

Optimal Control Theory
For the design of laser fields for driving the exciton dynamics we will employ optimal control theory (see, eq. Ref. [47]; the present implementation follows Ref. [3]). Here the goal is to find a laser field E(t) such as to optimize a target at a certain time t = T . This can be cast into an optimization problem for the functional whereÔ is the projection operator onto the target state, λ is the penalty factor for strong fields. In this work λ is treated as a Lagrangian multiplier such as to keep the integrated intensity close to I 0 . To simplify matters we will not employ the HE approach at this point but resort to the Markovian approximation, which leads to a set of two coupled equations.
Optimizing the functional J with respect to the field one gets the expression whereσ(t) is an auxiliary operator propagating backward in time, starting fromÔ at t = T with the following differential equation Eqs. (31) and (35)   The fields obtained by the OCT equations will be characterized by their XFROG (cross-correlation frequency-resolved optical gating [48]) trace defined as where G(t) is the gate function having the form In the above equation τ is the width of the gate and ∆ is the width of the shoulder.

Parameters
We will apply the formalism outlined in Section 2 to a model mimicking the situation in perylene bisimide aggregates.
The photophysical properties of PBI monomers have been studied in Ref. [41]. As far as the present model is concerned the following properties will be used: The S 0 -S 1 transition energy is E m,e = 2.13 eV. The monomer S 0 -S 1 transition dipole moment is 3.34 ea 0 . This transition is coupled to an effective vibrational mode of frequency 1415 cm −1 with a Huang-Rhys factor of 0.44. Further we used κ = 1 in Eq. (19). While the effective mode captured the vibrational side band observed in the experiment, it could not account for the general line broadening, which amounts to 1110 cm −1 thus indicating substantial coupling to a continuous distribution of low-frequency bath modes. Hence, for the local high-frequency mode we take in the spectral density, Eq. The situation is more complicated for the monomeric excited state absorption (ESA). On the experimental side, the ESA spectrum has been obtained from a pumpprobe spectrum, recorded after equilibration in the S 1 state. Since the procedure requires knowledge about ground state bleach and stimulated emission spectra at a detail which cannot be obtained, the extracted ESA spectrum may contain spurious features. On the theoretical side calculating highly excited states for a molecule as large as PBI is a challenge, in particular because standard time-dependent density functional theory will fail to describe transitions of double excitation character. In Ref. [41] we applied the multireference variant of density functional theory to a reduced PBI system. Combining experimental and theoretical data it can be stated that (i) the excited state responsible for the S 1 -S n transition is at n > 20. In other words, the internal conversion down to the S 1 state proceeds via a dense manifold of electronic states, whose microscopic description is out of reach. This suggests using a description in terms of an effective internal conversion rate as outlined in Section 2. For the actual value we have used an inverse rate of 100 fs. Within the spectral density model, Eq. (24), this amounts to choosing Λ IC =1000 cm −1 and η IC = 0.4. (ii) Since reliable information on the transition dipole moment µ m,f are not available we will use the value obtained within the harmonic oscillator picture, i.e. µ m,f = √ 2µ m,e [15]. (iii) Concerning the anharmonicity ∆ m = E m,f − 2E n,e we first note that in the range where E m,f ≈ 2E n,e two transitions have been observed. Since the magnitude of ∆ m decides about the mixing between local and nonlocal double excitation states [15] we will consider two cases. In case A the anharmonicity is large, ∆ m = −0.26 eV, whereas in case B it is small, ∆ m = −0.04 eV.
In Ref. [41] a PBI monomer has been investigated only. The Coulomb coupling between adjacent monomers in the PBI aggregate has been estimated from the experimentally observed absorption line shifts upon aggregation. It has also been calculated in Ref. [50]. Below we will use the calculated value of J mn = −515 cm −1 and only nearest neighbor coupling will be considered. The coupling between excited state transitions is not known and we will again use the harmonic oscillator approximation giving J  Figure 2. Population dynamics in the field-free case after initial population of a specific two-exciton eigenstate. Shown is the total population of the two-exciton manifold in comparison with the monomer case.  Figure 3. Population dynamics of the dimer in the field-free case after initial population of the highest state of the two-exciton band. Shown are the populations of the one-and two-exciton eigenstates. N A/N B). The one-and two-exciton eigenstates are shown in Fig. 1

Population Dynamics in Field-Free Case
In Fig. 2 we compare the internal conversion dynamics for the different models as a function of the initial preparation of the two-exciton manifold. Also shown is the decay of the monomer LDE state, which is always faster than that for the aggregate. Furthermore, it is found that the decay in case B is always slower as compared to case A, with the difference becoming more apparent upon increasing the aggregate size. This observation can be traced to the fact, that the mixing of LDE and NDE states in case B leads to a different intra-band dynamics and in particular to a slow-down of the interband dynamics as can be seen from Fig. 3. From this figure we notice that even though the decay of the initial state population is faster in case B, the population gets trapped for a longer time in state |2 1 , i.e. at the bottom of the two-exciton band. In case A the state |2 1 is of more local character as compared with case B and therefore inter-band relaxation is faster. Similar arguments apply to the aggregates with N = 3 and 4. So far we have presented results from the HE simulation. In Fig. 4   HE and the Markovian dynamics if the system is prepared in the uppermost two-exciton state. In order to scrutinize this effect we have plotted the population dynamics of the two-exciton eigenstates in Fig. 5. First, we notice that in general the populations from the HE are showing an oscillatory behavior, which is not present in the Markovian approximation. Despite this difference in details, the total two-exciton population is rather similar in the two limits. Only for case A and for initial preparation in state |2 3 we notice a marked deviation for the dimer considered in Fig. 5. In the HE case state |2 3 relaxes rapidly and states |2 2 and |2 1 become populated before the overall decay due to IC sets in. Notice that the energy gap between states |2 3 and |2 2 amounts to 2518 cm −1 . The spectral density for the system-bath interaction is composed of a lowfrequency part with cut-off at 100 cm −1 and a high-frequency contribution at 1415 cm −1 . Hence relaxation within second order perturbation theory as implied in the Markovian, i.e. truncated HE, approach will be very inefficient due to the smallness of the spectral density. The HE, on the other hand, accounts for higher order effects and yields a rapid  energy relaxation. Notice that in case of model B this energy gap is only 1203 cm −1 , i.e. within the frequency range covered by the spectral density. According to the argument given one would not expect a marked difference between HE and Markovian dynamics in case B, what is in line with the results shown in Fig. 4. Inspecting Fig. 1 we note that a similar energy gap exists also for larger aggregates what leads to a corresponding behavior of the population dynamics in Fig. 4. However, upon further increase of the aggregate size this gap closes and it can be expected that HE and Markovian dynamics behave rather similar at least from the total population point of view. In passing we note that ultrafast inter-band relaxation suppresses the effect of excitonic-polaron formation, which, however, will play a role in the one-exciton manifold (see, Ref. [39]).

Laser Control of Two-Exciton Dynamics
In the following we will investigate the possibility to control two-exciton dynamics with ultrashort shaped laser pulses. In particular we will ask the question whether one can transiently suppress EEA. Since EEA is a local process one might argue that a two-exciton wave packet, which is prepared in a way such that the contribution coming from LDE is minimized will lead to a slow down of EEA since the latter is mediated by intra-band relaxation processes mixing LDE and NDE states. For this purpose we will use the target state |1N within the OCT scheme, i.e. the state corresponding to the situation where the two NDEs have the largest separation in real space. The fastness of the IC process puts some restriction on the pulse duration which has been set to T = 50 fs. As mentioned in Sec. 2.5 the field-optimization will be performed in the Markovian limit, starting with a very broad band pulse. The field is then used to generate the dynamics using the HE approach.
In Fig. 6 we show the optimized fields (XFROG) for cases A and B together with the population dynamics of target states as well as the exciton eigenstates for the dimer model. Since dipole transitions are possible between adjacent exciton bands only, the dynamics necessarily involves a two-step process. First, the one-exciton band is excited (transition frequency ω 1 1 ,0 = 2.06 eV) and this process is followed by a one-to twoexciton transition. In case A this way a superposition of states |2 1 and |2 3 is prepared (transition frequencies ω 2 1 ,1 1 = 1.88 eV, ω 2 3 ,1 1 = 2.06 eV), which, however, rapidly decays due to intra-and inter-band relaxation processes. For case B, where the LDE and NDE states are strongly mixed, the compromise found by the OCT equations has been to prepare just state |2 1 which has a 42 % overlap with the target state. Although state |2 3 would have had a larger overlap, its transition dipole matrix element is about a factor of 13 smaller for that case (transition frequencies: ω 2 1 ,1 1 = 2.04 eV, ω 2 3 ,1 1 = 2.30 eV). Similar to the scenario of the field free dynamics in Fig. 3, the overall decay of the two-exciton population is slower in case B as compared to case A. Finally, we notice that the maxima of the XFROG traces do not match with the bare exciton transitions in all cases, since superposition states are prepared by the broadband excitation. Further, exciton-vibrational coupling and the dynamic Stark shift will modify the bare excitonic resonances.
One might argue that in the dimer EEA will always be very effective since the NDE are localized at neighboring sites. This situation might change for larger aggregates. As an example we consider the tetramer model in Fig. 7. The convergence of the OCT algorithm is rather slow in this case and for the given constraints only a small population of the target state can be reached. In both cases the pulse initially excites the lowest transition of the one-exciton band ( ω 1 1 ,0 = 2.03 eV). For case A the subsequent oneto two-exciton transition populates mainly state |2 5 ( ω 2 5 ,1 1 = 2.12 eV) which has the largest overlap with the target state (32 %). In case B state |2 1 is dominantly populated ( ω 2 1 ,1 1 = 2.01 eV) although it has only a small overlap with the target state (3%). However, as compared with state |2 4 , which has a 37 % overlap with the target state, the transition dipole moment to state |2 1 is larger by a factor of about 5. Overall, the comparison between cases A and B resembles again Fig. 6.
Comparing the dimer and tetramer cases one notices that the initial hypothesis that the longer the aggregate the longer the time scale during which one can maintain a two-exciton population does not hold in general. Instead the difference between cases    A and B points to the importance of mixing between LDE and NDE states within the two-exciton band. However, focussing on case B only there is the anticipated increase of the time-scale of transient two-exciton state population. In this case the optimized pulse populates the lowest state of the two-exciton manifold whose overlap with the target state decreases with increasing system sizes due to the "dilution" of the zeroorder excitation states. For the same reason, however, the overlap of state |2 1 with LDE states decreases yielding a longer time scale for the two-exciton decay.
Next we comment on the use of pulses derived by using the OCT equations and Markovian dynamics within the HE scheme. Inspection of the different cases shows that the resulting population of the target state is comparable in the two cases, i.e. this procedure does not lead to a degradation of control for the given pulse shape. Needless to say that in line with the discussion of Fig. 4  case A case B Figure 8. Comparison of laser-driven dynamics for the dimer models using the pulses obtained from the OCT equations (red) and a fit of these pulses to two Gaussians (2G, black). The lower panels show the resulting target and two-exciton manifold populations.
HE might result in different pulses. However, in view of the difficulties arising from the short life time of the two-exciton population this will have little relevance for the present discussion. We need to emphasize that our model is limited to the two-exciton space. In principle for the given field intensities it is likely that higher exciton manifolds could be excited as well (for a study of intensity-dependent multi-exciton dynamics within a Bloch model, see Ref. [52]). The expected rapid decay of these states would lead to an additional channel for populating the two-exciton manifold, what could in principle influence the details of the dynamics. Needless to say, that this effect could be suppressed by reducing the field intensity, but at the expense of the population of the two-exciton manifold.
Finally, the question arises whether all the details of the control pulses play a role or in other words how robust is the achieved population control and can the pulse be simplified. Exemplarily, in Fig. 8 we show a comparison of the dynamics obtained for the OCT derived pulses (cf. Fig. 6) and their fits to two Gaussian shaped pulses for the dimer model. For case A we observe a strong sensitivity, i.e. the two-exciton population drops down considerably when using two Gaussian pulses. Apparently, the preparation of the |2 1 and |2 3 superposition state depends on the details of the control field. In contrast, in case B where dominantly the state |2 1 is prepared a simple Gaussian-shaped field is sufficient to achieve a control comparable to the optimized pulse.

Summary
Two-exciton dynamics has been studied on the basis of a non-perturbative and non-Markovian hierarchy equation approach. The dissipative dynamics has been combined with optimal control theory for obtaining laser fields designed such as to trigger longlived two-exciton state populations. Specific results have been obtained for short aggregates made of J-aggregate forming perylene bisimide dyes for which Frenkel exciton parameters are available [41]. Establishing a two-exciton population one has to compete with the very efficient internal conversion, which is a consequence of the breakdown of the Born-Oppenheimer approximation. For the monomer this process takes place on a time scale of about 100 fs, which therefore sets the upper bound for laser control. The fact that two-exciton populations can be maintained on a longer time scale is due to the mixing between local and non-local double excitations of the aggregate; the latter do not decay on an ultrafast time scale. The considered dyes support two possible double excitation states in the relevant energy range. Therefore, we considered two scenarios corresponding to small and substantial mixing between the two type of aggregate excitations. It turned out that the case of strong mixing allows for maintaining a twoexciton population on a longer time scale (in the present cases about 1 ps as compared with 0.5 ps for the weak mixing case). This effect should be observable in a pump-probe experiment. At this point it should be noted that two-exciton populations in aggregates made of organic dyes like PBI have been observed for the first time only very recently [53], what demonstrates feasibility of preparation and spectral identification.
The present investigation highlights the importance of zero-order state mixing within the two-exciton band for the inter-band transitions. The more diluted the zeroorder states the more the decay is slowed down. However, this situation might change if static disorder is taken into account which will lead to a localization of the two excitations on different parts of the aggregate.

Appendix: Derivation of the Hierarchy Equations of Motion
For the Caldeira-Leggett model of dissipation, Eq. (19), withĤ S−B = f (ŝ)g(b), one can employ a stochastic procedure to derive the equation of the motion for the reduced density matrix of the relevant systemρ if the whole system is prepared as a factorized state, i.e.ρ tot (0) =ρ(0) exp(−βĤ B )/Tr[exp(−βĤ B )] [37,38]. Within this approach, the influence of the bath is completely characterized by its response function and its effect is described by a random forceḡ(t) in the Itô stochastic differential equation where w 1,t , w 2,t are two independent complex-valued Wiener processes and w * 1,t , w * 2,t are their complex conjugates. In Eq. (A.1)ḡ(t), which plays the same role as the influence functional in the path integral treatment, can be regarded as a stochastic field fully characterizing the influence of the environment. It is defined as For simplicity, here it is assumed that the response function of the bath is where Ω 1 = Ω * 2 (cf. Section 2.3.2). Note that the extension to multiple exponentials is straightforward. For a more general form of the response function consult Ref. [38,33].  = |b 1 b 2 |. Notice that using a scaling like in the above equation has been suggested by YiJing Yan and co-workers [36,51]. The present prefactor is slightly different from that suggested in Refs. [36,51], where the numerator was set equal to one. This choice will keep the terms in the same order size consistent in cases where the decay constants in the response functions are rather different.
Carrying out the stochastic average for the random variables in Eq. (A.1), elementary stochastic calculus yields the differential equation forρ mn (t) This equation needs to be solved for the initial conditionsρ 00 (0) =ρ(0) andρ m,n (0) = 0 (m = 0 or n = 0). Application to the present system-bath model yields the equations of motion, Eq. (30).