Funneling dynamics is a key mechanism in the framework of artificial light-harvesting processes, biological antennas, or organic photovoltaic devices^[1–3]. Transfer from an initial donor to a target acceptor proceeds through an excited state gradient, leading to a final localization from which the photonic energy may be captured. In complex molecular systems, such energy transport mechanisms are driven by ultra-fast non-adiabatic transitions through conical intersections (CIs)^[4–9]. Investigation of the coherent preparation of the donor state from the initial ground state and the control of the subsequent dynamics using intense, ultra-short laser pulses are major issues for modifying the energy transport mechanisms^[10–12]. In this context, we have recently studied a rather complex molecular system, namely, polyphenylene ethynylene dendrimer, with a control objective aiming at the coherent preparation of two donor states involved in the dynamics, or their symmetric versus asymmetric superposition^[13]. In these systems, the acceptor state is radiatively coupled to the ground state. It is only by discarding this direct excitation that we could build control strategies to increase the coherence (asymmetry) lifetime to populate the acceptor through the vibrational baths of the donors.

In the present work, we address the coherent control of the pyrazine (C4H4N2) photophysics and more precisely its S2→S1 internal conversion leading to challenging issues, both with respect to its ultra-fast non-adiabatic CI transfer (internal conversion time scale of 20 fs^[14–18]) and to its 24-dimensional (24D) vibrational excitation dynamics. Some previous works in the literature are also concerned by the coherent control of population transfer from S0 to S2 and S1 electronic states, putting the emphasis on a significant delay of the internal conversion (or, in general, radiationless transition), which is considered as an undesired phenomenon^[19–21]. This is in relation with the fact that a successfully populated excited state (S2) would be readily available for applications (such as the observation of desired photoproducts, mode selective chemistry, light-triggered molecular rotors, separation of racemic mixtures) for only a brief period of time. In this respect, Refs. [20,21] focus specifically on an optimal control scheme to suppress radiationless transitions on nanosecond (ns) time scales after the external control is over. As for Ref. [19], the authors proceed to a detailed study of the S0→S2→S1 process, putting the emphasis on maximal S2 population during the control scheme (about 50 fs duration). Still other works^[22,23] refer to more general schemes by using overlapping resonances. The main computational limitation is identified as the reduced number of vibrational modes, usually four rather than 24. However, within some approximations, model extensions to 24 dimensions have been achieved, making use of a so-called QP (Q denotes the small subspace of interest and P the rest) partitioning^[24]. Motivated by light-harvesting systems, we pursue an opposite goal. Our purpose is to maximize the population of the energetically lower and thus more stable excited state S1. The basic idea is to deposit the exciton energy from the bright donor state S2 to the dark acceptor state S1, through a funneling mechanism triggered by a CI and non-adiabatic internal coupling, taking advantage of a decreasing electronic gradient from S2 to S1.

Control objectives being opposite, the strategies we are developing also differ completely from the ones referred to in previous works. More precisely, when maximizing the S2 population, the emphasis is put on stable, electronically localized eigenstates among the very dense manifold of vibronically coupled levels. References [20,21] are concerned by an optically coherent superposition of a specific eigenstate localized in the S2 potential well. Reference [19] also deals with high-lying S2 vibrational states with significantly slower decay rates into S1. A complete study of excitation frequencies is undertaken to resonantly reach the superposition of such excited vibrational states, using rather long duration pulses (180 fs for Refs. [20,21] and 40 to 100 fs for Ref. [19]). Some other works deal either with the dynamic Stark effect used as a basic mechanism to shift the CI away from the Franck–Condon (FC) region^[25–27], or with optimal control theory^[28], but they are limited to three or four modes of vibronic couplings. We are adopting different control strategies based on two mechanisms: either interference or sudden excitation kicks, as extended to high-dimensional complex molecular systems. This is conducted by using ultra-short broad-band pulses, without a specific resonant excitation frequency. We note that such strategies remain in the spirit of what has been suggested in the preliminary work of Ref. [29]. Our control schemes are built so as to decrease the ground state S0 population and simultaneously reach the optimal contrast favoring S1 population as compared to the one of S2. As a signature of a stable deposit of electronic energy, we are interested in the post-pulse S1 population, taking into account the dephasing processes of the full 24-mode dynamical description of pyrazine. On the experimental side, it is worth mentioning that several works are addressing the control of the internal conversion to the energy transfer ratio, using in particular, open learning loop setups and evolutionary algorithms^[30–33]. In particular, the multi-pulse excitation (a degenerate four-wave-mixing sequence) has its counterpart in the present work, in terms of ultra-short pulse trains to propel the wave packet towards the CI.

The manuscript is organized as follows. Section 2 is devoted to the Frenkel vibronic Hamiltonian involving the 24-mode model of pyrazine, with all relevant parameters taken from Ref. [34]. The methodology adopted to treat this model is the multi-configurational time-dependent Hartree (MCTDH)^[34–39], which is shortly reviewed with respect to its numerical convergence criteria in the CI branching space. The post-pulse field-free dynamics and its recursive patterns in terms of wavepacket vibrational periods on both S2 and S1 potential energy curves are analyzed. A detailed interpretation is provided when going from 2D to four-dimensional (4D) and 24D models. Control strategies are presented in Section 3, putting the emphasis on control observations and parameters on the one hand and interference and kick mechanisms from weak to strong field regimes on the other hand. Results are presented in Section 4 for both the mechanisms with their interpretation. Conclusions and some perspectives are given in Section 5.

As in many previous works, we are considering three electronic states of pyrazine, namely, the ground state S0(A1g), together with the two lowest excited electronic states, S1(B13u) and S2(B12u). These states are well-separated from others, close in energy, and vibronically coupled to each other^[15]. Figure 1 displays a schematic view of the corresponding states and their couplings.

Following the usual model of displaced harmonic oscillators in open quantum systems, we work with a linear vibronic Hamiltonian expressed in a diabatic representation. In this section, we examine the field-free vibronic Frenkel Hamiltonian, the MCTDH method retained for solving it, and the generic field-free dynamics after launching part of the ground state wavepacket on the excited states. That is more precisely conducted with successive 2D, 4D, and 24D models for a clear understanding of the role of dimensionality.

We briefly recall the generic field-free Frenkel Hamiltonian, which reads^[40]HS(Q)=∑n|n〉Hnn(Q)〈n|+∑n≠m|n〉Hnm(Q)〈m|,(1)where n and m=0,1,2 denote the electronic states and Q collectively the nuclear coordinates. Using the linear vibronic coupling of Ref. [34] in mass weighted coordinates, we have Hnn(Q)=ϵ˜(n)+∑αNα[Pα2+ωα(n)2(Qα−dα(n))2],(2)and Hnm(Q)=∑αNαλα(nm)Qα,(3)where Qα and Pα=−iℏ∂∂Qα are the coordinate and momentum associated with the mode α. There are Nα such modes in the electronic state n with frequency ωα(n) and displacement dα(n). Reorganizing Eq. (2), one gets Hnn(Q)=ϵ(n)+∑αNα(Pα2+ωα(n)2Qα2)+∑αNακα(n)Qα,(4)where the gradient at the reference point (Qα=0) is given by κα(n)=−2ωα(n)2dα(n), and the gradient of the interstate coupling is noted as λα(nm)=W12. Finally, the renormalized site energy is ϵ(n)=ϵ˜(n)+ωα(n)2dα(n)2. Table 1 of Ref. [34] collects all the parameters that are relevant for the system Hamiltonian given by Eq. (1). We have to position S2 with respect to S0, and, for this, we take the value 4.84 eV from Refs. [18,41].

It is to be noted that the two normal modes that span the CI branching space, using the symmetry group notations, are Q6a for state S1 and Q10a for state S2^[40]. As mentioned before, the three models will be referred to as 2D (Nα=2) involving Q6a, Q10a, 4D (Nα=4) when enlarging the branching space to Q9a and Q1, and the full 24D (Nα=24) description involving all normal modes.

The MCTDH approach for multi-electronic states and multi-mode nuclear dynamics has been presented in detail in the literature. For pyrazine, with its three electronic states in consideration, the multi-state vibrational wave function is taken as a three-dimensional column vector (Ψ1,Ψ2,Ψ3)T, of single-state nuclear wave functions Ψn(n=1,2,3), T being the transpose. The total electro-nuclear eigenvector is then written as |ΨTot(Q,t)〉=∑n=13Ψn(Q,t)|n〉.(5)

n labels the electronic states, and the unknown nuclear wave functions are solutions of the following close-coupled equations: −iℏ∂∂tΨn(Q,t)=HnnΨn(Q,t)+∑m≠nHnmΨm(Q,t).(6)

A standard multi-configuration approach consists of expanding these functions on a time-independent basis set ΦJ,n, with time-dependent coefficients AJ,n, as Ψn(Q,t)=∑JNJAJ,n(t)ΦJ,n(Q).(7)

The major improvement brought by MCTDH is that not only the coefficients but also the basis functions are taken as time-dependent. The challenge is that, by adapting the basis set functions to the temporal evolution, one must reduce the total number of these functions (i.e., NJ) for a given convergence criterion. To proceed along this line, the now time-dependent basis set functions Φ˜J,n(Q,t) are given as a tensorial product of time-dependent single-particle functions (SPFs) φjα(α)(Qα,t) describing a given nuclear degree of freedom α. The complete expansion reads as Ψn(Q,t)=∑j1=1n1…∑jf=1nfAj1,…jf(t)∏α=1fφjα(α)(Qα,t),(8)where f is the number of nuclear degrees of freedom. The index jα stands for one of the nα possible single-particle functions for the αth nuclear degree of freedom. The number of configurations is thus given by the product n1…nf. The single-particle functions are ultimately expressed in time-independent, so-called primitive basis set functions χik(α)(Qk) as φjα(α)(Qα,t)=∑iαMαciα(α,jα)(t)χiα(α)(Qα).(9)

Interestingly enough, there is no prescription for SPFs to depend on a single coordinate, and they may depend on several coordinates. To further reduce the memory and the numerical effort, several physical coordinates can be regrouped into a so-called combined mode^[34]: Q˜κ≡(Qκ,1,Qκ,2,…,Qκ,d),(10)φj(κ)(Q˜κ,t)=φj(κ)(Qκ,1,Qκ,2,…,Qκ,d,t),(11)where d is the number of nuclear degrees of freedom in the combined mode Q˜κ. The MCTDH wavefunction is now expanded as Ψ(Q˜1,…,Q˜p,t)=∑j1⋯jpn1⋯npAj1⋯jp(t)∏κ=1pφj(κ)(Q˜κ,t),(12)where p is the number of combined modes, and the SPFs themselves are expanded as φj(κ)(Q˜κ,t)=∑i1⋯idMκ1⋯MκdCi1⋯id(κ,j)(t)χ(κ,1)(Qκ,1)⋯χ(κ,d)(Qκ,d).(13)

Here, we use p=8 combined modes Q˜j (j=1,…8). More precisely, Q˜1 results from a combination of Q6a and Q10a. For the others, the combination scheme is indicated in parenthesis, i.e., Q˜2(1,9a,8a), Q˜3(2,4,5), Q˜4(6b,3,8b), Q˜5(7b,16a,17a), Q˜6(12,18a), Q˜7(19a,13,18b,14), and Q˜8(19b,20b,16b,11).

As already said, we use the normal coordinates of pyrazine in its electronic ground state for the dynamics (f=24). For the primitive basis set [χ(κ,j)(Qκ,j) correspond to the primitive basis functions for the nuclear coordinate Qκ,j and Mκj their number], we have used 22 harmonic oscillator discrete variable representation (DVR) functions for the nuclear (normal) coordinates Q6a and Q10a, whereas only 12 and 11 such functions are retained for the normal coordinates Q9a and Q1 normal modes, respectively. All other nuclear coordinates, not playing a crucial part in the dynamics, are merely described by four DVR functions. Numerical convergence is obtained by a much smaller basis set (n1⋯n8) of single-particle functions, which are no more one-dimensional (φj(κ)(Q˜κ,t)), as given in Table 1.

The equations of motion for the coefficients A of Eq. (8) result from the time-dependent evolution equation involving the total Hamiltonian, i.e., the one of the system, together with the molecule field coupling. The close-coupled system of differential equations [Eq. (6)] is solved by projecting on the basis function of the combined modes. The initial condition, at time t=0, being taken as one for the vibrationless ground state n=1, and zero for all other states, the nuclear wave functions are built following Eq. (8). Finally, the time-dependent population in each electronic state is given as Pn(t)=∫|Ψn(Q,t)|2dQ.(14)

In order to have a typical generic overview of the S2→S1 population transfer, we now proceed to three dynamical calculations extending over a period of 500 fs. More precisely, we successively analyze 2D and 4D approximate models and relate them with the full 24D model. The minimal 2D model merely involves the two normal modes α=6a and α=10a building up the CI branching space. The population transfer dynamics are induced by some Franck–Condon vertical launching of the initial vibrationless (vα=0 for all α) S0 wavepacket on the excited states. To fix the ideas, this is practically done referring to a low intensity (5×1012W/cm2), short duration (14 fs), resonant laser pulse (carrier-wave frequency of 4.8 eV), which is symbolized by the vertical blue arrow of Fig. 2 positioned at time t=0. Note that several time-resolved photoelectron spectroscopy experimental studies have been reported using pump pulses around 4.7 eV: they do not report any dipole-allowed resonant (1 + 1) transition to a higher lying electronic state but a slow decay (of around 20 ps) to the electronic ground state; see the introduction of Ref. [42].

It is worthwhile noting that the precise parameters of this excitation are an illustrative example of the kind of sudden excitation we refer to initiate the transfer dynamics. During the 14 fs S0→S2 resonant excitation process, it is basically the S2 state that is populated up to P2=0.37, whereas P1 does not exceed 0.03. The short pulse here is just to initiate a typical population transfer from the initial ground S0 to the bright donor state S2, without any control purpose. Actually, what we call free-field dynamics is to be understood as the wave packet dynamics following the switch-off of this pulse.

The post-pulse field-free dynamics starts with this configuration and proceeds towards the descending energy gradient from S2 to S1 electronic potential energy surfaces coupled by W12. In particular, due to the CI, an important amount of population is transferred to S1. After a delay corresponding to the S1 vibrational period (T1=56fs with the parameters of our model), the wavepacket reflected on the outer right turning point returns back to the CI region with a partial back transfer to S2 and then evolves with another oscillation with the vibrational period (T2=44fs) of the state S2. The combination of these two oscillation modes in the vibrational baths of S1 and S2 gives rise to a recurrent pattern in the population dynamics that is then periodically repeated. Figure 2 provides a complete illustration of such recursive patterns. Calculations based on longer propagation times show that the revival structures last for more than 1 ps. They are progressively attenuated, and their periodicity is slightly altered due to wavepacket broadening. We now proceed by including higher degrees of freedom through additional normal modes. Figure 3 displays such calculations over 500 fs of propagation time, with the same excitation initiating the population transfer from S0. It is to be noted that the additional degrees of freedom in multi-dimensional calculations dramatically increase the density of vibrational states. As a consequence, the specificity of the T1 and T2 vibrational periods of the α=6a and α=10a normal modes, building up the 2D model, is partially lost. The wavepacket is spreading over all the additional modes through all the vibrational states that form a quasi-continuum. Typically after 250 fs, the revival structures are completely suppressed. The P1 and P2 populations are moderately oscillating and progressively stabilizing in time.

The pyrazine molecule is studied in planar geometry and assumed to be oriented in a plane Oyz orthogonal to the propagation direction Ox of the electromagnetic field E(t). The laser is supposed to be linearly polarized along the Oz axis. The time-dependent total Hamiltonian is written in the length gauge and within the dipole approximation as H(Q;t)=HS(Q)+V(Q;t),(15)HS being the molecular Hamiltonian taken from Eq. (1) and V the radiative coupling: V(Q;t)=−μ(Q)E(t).(16)

The transition dipole matrix elements between the ground S0 and excited S1, S2 states are, respectively, noted as μ01 and μ02. Their explicit spatial expansions in terms of the normal mode coordinates are given by^[26]μ01(Q)=ξ10a(01)Q10a,(17)μ02(Q)=μ02(0)+∑α≠10aξα(02)Qα+12ρ10a(02)Q10a2.(18)

It is worthwhile noting that for symmetry arguments, μ02 is much larger than μ01 due to its permanent dipole component μ02(0). Numerical values we are using are extracted from Ref. [27]. It is also to be noted that we are working in such field conditions where referring to the polarizability is not necessary, as opposed to Ref. [27], concerned by very strong intensities and excitation conditions far from resonance. We will now examine the relevance of some control observables, together with two strategies to reach them in an optimal way, by exploiting two different mechanisms; namely, interference and kicks.

The primary objective of this study is to optimize the population P1 of the acceptor state S1. But for pyrazine, the acceptor state S1 is a quasi-dark one, as is clear from Eq. (17), practically without the possibility of a direct radiative excitation from the ground state. The population transfer process first proceeds with the radiative excitation of the ground S0 state to the bright excited S2, which then acts as a donor. The second step is monitored by the interstate coupling W12, leading to a CI-induced strong non-adiabatic transfer from S2 to S1. A rather intuitive control objective could be taken as the maximization of the ratio r(t)=P1(t)/P2(t). But, this goes with the difficulty of getting very high ratio r, between actually very small excited states populations, that is without any practical interest. To properly take into account the two steps, excitation and funneling, we built a population contrast type control observable combining two requirements: (i) maximization of the excited states populations, which actually amounts to the minimization of the residual population P0 of the ground state S0 and (ii), among the two excited states, bring the maximum of population on the dark acceptor state S1 or, in other words, increase the contrast between excited states population. Such a time-dependent contrast is given by p˜(t)=P1(t)−P2(t)P0(t).(19)

It is worth noting that, in terms of r, Eq. (19) can be written as p˜(t)=r(t)−11/P2(t)−[r(t)+1].(20)

Maximizing p˜ requires simultaneously maximizing the numerator, that is r>1, and minimizing the (always positive) denominator, which is achieved with P2(t)<1/[1+r(t)]. Finally, combining these two requirements, in particular limiting the P2 increase, our control objective goes much beyond the passive observation of the acceptor state population P1 increase (i.e., increasing the donor state population P2 and waiting for its CI-mediated transfer towards the acceptor).

On numerical grounds, in weak fields, P0 would be considered not less than 0.5, meaning that no more than half of the initial population is assumed to be transferred to the excited states. As a consequence, even in the case of the best contrast that could ideally be reached between the donor and acceptor states (that is P1=0.5 and P2=0), the maximum possible value for p˜ should not exceed one. In other words, it is only in the strong field regime that we could expect a contrast p˜>1.

Finally, instead of the time-dependent contrast given by Eq. (19), we choose a scalar control observable by taking time averages of p˜(t): p=1tmax−tmin∫tmintmaxp˜(t)dt.(21)

This involves the definition of time windows over which the average is performed. A local contrast can be defined over a time window displaying the maximum of contrast, which happens during the first vibrational pattern, lasting over 45 fs, as illustrated in Fig. 2, and leading to the choice [tmin,tmax]=[45fs,90fs]. Although, during this 45 fs time window, some experimental observables are reachable, more interestingly, an asymptotic contrast is adopted in the following defining a semi-finite time interval above 300 fs. Incidentally, it is to be noticed that 300 fs corresponds to typical electronic coherence times for such molecular systems [the dissipation comes from the intramolecular vibrational redistribution (IVR)] and an upper limit of 500 fs is taken as our final propagation time, leading to [tmin,tmax]=[300fs,500fs].

In the following, we refer to two laser control strategies based on two mechanisms, namely, interference and kicks. In both cases, the control field is made up of a train of N time-delayed individual pulses. Two such pulses are illustrated in Fig. 4. The general expression for the field involved in Eq. (16) is given by E(t)=∑i=1NIEi(t)sin(ωt),(22)where I is the laser leading intensity and ω its frequency. Ei(t) is the pulse envelope, taken as Ei(t)=sin2[πT(t−ti)]H(t−ti)H(ti+T−τ),(23)where T is the pulse duration. H is the Heaviside function, being zero or one depending on whether its argument is negative or positive. The time intervals ti are given as ti=(i−1)τ, where τ is the time delay between the two successive pulses. Moreover, we wish to disentangle the roles of strong field and interference effects. When referring to our control observable, a relevant comparison would require excitation conditions where the same total electromagnetic energy is deposited in the system by the pulse train over its full duration and given by A=∫0∞E2(t)dt.(24)

More precisely, we compare pulse trains with the same total energy, namely, the one provided by the single pulse. To fix the ideas, for two identical pulses E1=E2=ϵ, A can be written as A=a∫0∞ϵ2(t)dt.(25)

Referring to Fig. 4, three cases are examined: (i) τ≥T, no overlapping pulses, resulting in a=2; (ii) τ=0, full overlapping, with a=4; and (iii) 0<τ<T, partial overlapping, leading to a delay-dependent a (2<a<4), to be calculated according to Eqs. (24) and (25). Our control fields amplitudes are normalized so as to provide the same total electromagnetic energy. This is practically done by dividing the nominal intensities I by a, or equivalently, the electric field amplitudes by a.

We also wish to emphasize that when developing these strategies, we always have in mind the experimental feasibility with respect to pulse duration, intensities, and maximum number of pulses in the train. For the interference scenario, the laser leading frequency ω is taken to be approximately resonant with the excited states vibrational levels at the Franck–Condon vertical region of the ground state. The optimal value that is adopted is ω=4.77eV, i.e., a wavelength of λ=260nm in the UV region. The pulse duration T is calculated such that the corresponding energy band broadening covers enough vibrational levels taking part in the funneling process. For the kicks scenario, the resonance condition is no more a relevant requirement, since the ultra-short pulses involved in this strategy would lead to excitation with even broader energy bands. All ultra-short pulses are taken as identical and typical values, which are retained (T=14fs or even T=10fs) and correspond to band broadening ranging from 2382cm−1 up to 3335cm−1. Following a few attempts to roughly optimize the frequency ω and the duration T, we are finally left with two main control parameters, namely, the field intensity I and the time delay τ.

The first strategy we are referring to relies on vibrational wave packets interference as a basic mechanism, envisioned either in a pump-probe or a pump-pump process. More precisely, as illustrated in Fig. 1, the pump-probe process involves two routes to reach S1, starting from the ground state S0. Route 1 (pump pulse) proceeds through the intermediate state S2 with the lowest-order transition amplitude given by^[43]T1=V+VG0V,(26)where G0 is the lowest-order Green’s function of the system Hamiltonian involving the diabatic electronic states, and V stands for interstate (W12) and radiative (μ20E1) couplings. More precisely, for Route 1, in the absence of direct radiative coupling S0→S1, one has T1=〈S1|W121ES0−ES2+ℏωμ20E1|S0〉,(27)where ESi are the energies of the corresponding electronic states, and E1 is the pump pulse electric field. Route 2 (probe pulse) is the direct transition S0→S1, with a transition amplitude: T2=〈S1|μ01E2|S0〉,(28)with E2 the probe pulse electric field. When applying the two pulses with a time delay τ, the pump-probe strategy finally results in an S1 population given by P1=|T1+T2e−iEτ/ℏ|2,(29)with E=ES0+ℏω. Due to very low transition dipole μ01 as compared to μ02, the transition amplitude T2 is actually negligible unless very strong field amplitudes E2 are used with possible ionization or dissociation damage on the molecule. This is why we would rather focus on a pump-pump process retaining only Route 1 and resulting in P1=|T1+T1e−iEτ/ℏ|2.(30)

Finally, the interference scheme we are considering through Eq. (30) is between the transition amplitudes T1, with a controlled delay τ. Such interference mechanisms have already been recently exploited in the control of rotational anisotropy^[44]. A physical understanding of the mechanism can be obtained by considering the vibrational wavepacket back and forth oscillations in the excited states potentials. Following a first laser pulse leading to a vertical Franck–Condon launching from S0 to S2, the vibrational wavepacket oscillates during its early dynamics t<70fs in the excited S2 and S1 states harmonic potentials. A second pulse is then applied with a controlled delay τ launching a second wavepacket from the ground state S0, which can interfere with the first one, as they may overlap if the delay is adequately chosen. A constructive interference would produce an amplitude enhancement close to the CI, and thus an efficient population transfer from S2 to S1. Later on (t>70fs), due to wavepacket spatial broadening and also its dispatching over additional degrees of freedom, the control efficiency is expected to be lost, as a consequence of less important successive wavepacket overlapping.

The interference mechanism can be extended and complemented by a kick mechanism, in the spirit of the one already discussed previously in the context of molecular alignment or orientation processes^[45]. At this respect, we are considering a train of N=5 identical individual pulses in Eq. (22). The kick mechanism assumes sudden excitation of the system. We are therefore addressing even shorter duration T=7fs to 10 fs pulses without the specific need for a resonance condition. The molecule can be kicked every time when P1 starts decreasing.

As has previously been stated, the control observable given by Eq. (19) highly depends on the field intensity. In particular, it is only for strong field regimes that it can take the values exceeding 0.5. In order to fix the intensity regimes for the present system described within a 4D model, we proceed with dynamical calculations using both a single pulse or two (equal fluence) pulses with a typical delay of 13 fs and for a series of increasing intensities. Figure 5 plots the bright donor state population P2 at the final time (t=50fs) as a function of the field leading intensity. As is clear from the figure, for the single excitation, the weak field regime extends up to intensities I=1013W/cm2, for which P2 increases rather linearly. The strong field regime shows non-linear behaviors: first, a saturation in S2 population, for about I=2×1013W/cm2, followed by a decrease down to I=7×1013W/cm2, corresponding to partial population trapping in the ground state S0, and then again an increase. We checked that the limit between weak (linear behavior) and strong fields (non-linear behaviors) remains practically unchanged when applying two pulses of equal fluence with a typical delay of τ=13fs, as will basically be the case for the following control issues.

We will successively examine the efficiency and robustness of the interference and kick mechanisms as implemented in their respective control strategies, when going from a reduced 4D to a full 24D model. The guiding principle is to define the parameters of the control field in a 4D realistic and tractable model. This field, once obtained, is later used in a full-dimensional dynamic calculation including all 24 degrees of freedom. The challenge is to discover the efficiency and robustness of such a control field that persists when confronted with the presence of the other numerous degrees of freedom of the molecule.

Having in mind the general post-pulse evolution of the excited states populations as illustrated by Fig. 2, two possible interference schemes could be envisioned. One concerns the nuclear wavepacket early dynamics for times typically less than t<70fs within the first vibrational periods, the second corresponding to longer times, typically t>70fs when the wavepacket revisits the FC region, following the revival patterns that are observed in Fig. 2. The asymptotic contrast, as defined by Eq. (21), is calculated in a 4D model as a function of the delay τ between the two pulses of fixed duration T=14fs. All calculations are conducted within the appropriate field amplitude renormalization condition [cf. Eq. (24)] leading to the same fluence as the one of the single pulse. The results are displayed in Fig. 6. Following the analysis of Fig. 5, we consider three laser intensities: two pertaining to the weak field regime I=5×1012W/cm2 and I=1013W/cm2 and one to the strong field regime I=8.8×1013W/cm2. Weak fields, apart from being easier to realize experimentally, have the advantage of better supporting the first-order perturbation approach of Eq. (28). However, the rather low contrast obtained for I=5×1012W/cm2 makes us favor the choice of I=1013W/cm2 for a typical illustration of the weak field regime. In this regime, the optimal contrast p=0.8 is obtained for a delay of τ=25fs, which, as discussed before, denotes a rather good efficiency, taking into account that less than 20% of the initial state population is transferred to the donor state S2, as seen from Fig. 5.

As a clear indication of the efficiency of the optimal control strategy, we compare the two population transfer dynamics resulting from either a single pulse or two pulses of equal fluence with the optimal delay of τ=25fs (Fig. 7). The single pulse depletes the ground state up to P0=0.65. The second pulse in constructive interference, producing a depletion up to P0=0.38. The contrast is increased by a factor 3.3 from p=0.24 up to p=0.8. The final population, which is deposited on state S1, is P1=0.5, meaning that about 50% of the ground state population is in the acceptor state, while only 16% remains in the donor.

The strong field regime offers much better contrasts by depleting more of the ground state population P0. The results displayed in Fig. 6 show that the maximum is reached for a delay of τ=6fs and leads to p=1.55 as a contrast. Moreover, fast oscillations as a function of the delay show the sensitivity of the observable to this delay, which is another signature of interference effects. Four characteristic delays are considered in the 4D reduced model for an intensity I=8.8×1013W/cm2. The single pulse (corresponding to τ=0fs, after equal fluence renormalization) and τ=9fs (corresponding to a local minimum in Fig. 6) result in comparable dynamics leading to small asymptotic contrasts (typically p<0.3). The third choice for the delay τ=6fs corresponds to the optimal value reached during the early dynamics p=1.55, whereas the fourth choice τ=75fs stands for the optimal value p=0.9 obtained during the late dynamics. In an a posteriori way, we also checked the sensitivity and optimality with respect to the resonance condition. The first rows of Table 2 for the interference scenario collect the asymptotic contrasts for two off-resonance conditions. The frequency ω=4.56eV below the resonance is better fitted for a population transfer towards state S1. But, as such transfers cannot operate in a direct way, the efficiency in terms of contrast is not satisfactory. The frequency ω=5.15eV above the resonance, although opening the possibility for some interesting superposition of vibrational levels mixing the two excited states, turns out to be less efficient. This study validates the choice made for the resonance frequency, as the presumably optimal one.

Finally, we introduce the optimal control parameter values obtained from the reduced 4D dynamics into a full 24D model describing population evolution. The results are displayed in Fig. 8 for the early dynamics and Fig. 9 for the late dynamics. It is worthwhile noting that such full-dimensional calculations are done once the laser parameters are optimized within the context of the reduced 4D model. As for the 24D calculations, they necessitate a central processing unit (CPU) time of about 32 h using the OpenMP parallelization scheme with 16 processors on an Intel E5-2665 computer. The most striking observation is that the population evolution in the 24D full dynamics closely follows the one of the 4D calculations. Both end with very comparable asymptotic contrast exceeding p≃1.6 for τ=6fs and p≃0.9 for τ=75fs. Following this observation, late dynamics control turns out to be less efficient. Actually, better interference schemes are operating for early dynamics where the vibrational wavepacket dynamics is more accurately periodic, with better marked revival structures, better adapted for the interference mechanism. The asymptotic contrasts obtained in the 24D model are comparable with the ones of the reduced 4D model, even very slightly better. This is presumably due to the fact that the wavepacket is more efficiently and rapidly dispatching towards the additional degrees of freedom, resulting in a faster population stabilization between the excited states. As a conclusion, the optimal result obtained for an early dynamics control scheme gives an acceptor dark state population of about P1=0.6, which means that 60% of the ground state initial population is transferred to the acceptor S1, while only 16% remains on the donor S2.

The kick mechanism operates with ultra-short (broad-band) pulses, leading to successive sudden, non-resonant momentum transfer to the molecule^[45]. Figure 10 displays the asymptotic contrast as a function of a constant single delay τ between the kicks imparted to the system from a train of five identical pulses, with fixed intensity I=5×1012W/cm2. Here, we deal with two control parameters, namely, the pulse duration T and the delay τ. It is important to note that a precise value for the excitation frequency is not very relevant, since the broad-band character of the individual pulses renders the strategy essentially non-resonant. As with the previously discussed mechanism, once the optimal pulse parameters are determined on the reduced more tractable 4D model, they are transposed to the full 24D dynamics. More precisely, we obtain from Fig. 10 an optimal asymptotic contrast p=1.71, for individual pulse duration T=10fs and delay τ=13fs. Here, also we checked the sensitivity of the results with respect to the excitation frequency ω. The last rows of Table 2, corresponding to the kick scenario, collect the results for three off-resonant frequencies. As expected, the contrast is not very sensitive to frequencies within a reasonable window covering ω=4.77eV to ω=4.91eV, leading to almost the same p=1.7. It is only far from these values that the results differ but with a less efficient control.

Figure 11 displays the corresponding time evolution of the electronic states populations. As can be seen from the step-by-step decreasing evolution of P0, at each kick, the momentum imparted to the system induces an increasing population towards the excited states. Our primary aim was to increase regularly the population P1. Such a control could presumably be possible with different delays between the kicks. But, for easier experimental requirements, we assume in this control scheme a maximum of equally delayed, five identical pulses pertaining to a weak field regime I=5×1012W/cm2. With such restrictions, the population P1 is not regularly increasing, even during the time when the kicks are applied. Another strategy could be to apply kicks every time P1 starts to decrease, as has been previously done in Ref. [45]. Targeting experimental feasibility, within this restrictive parameters sampling space, we get very encouraging results, with contrasts up to p=1.75. We also observe from the full-dimensional calculation of Fig. 11 that 60% of the population of the ground state is transferred to the dark acceptor, while only 16% remains on the donor. About 54 CPU hours are needed for this calculation with the same computer as mentioned before. It is worthwhile noting that short duration and low intensity five-pulse trains can do as well, even better than the two interfering long duration (resonant) intense pulses, as can be observed from Figs. 8 and 11.

In this work, we consider the intense laser control of funneling dynamics between electronic states with decreasing energy gradient and passing through a CI. The pyrazine molecule offers such a framework involving two excited electronic states, presenting a strong non-adiabatic coupling in relation with their CI. Moreover, the initial ground state is radiatively coupled to the highest energy excited state, which is considered as a bright donor (D). We are concerned by the population transfer dynamics from this state, to the lowest energy dark acceptor (A) state, motivated by the possibility of long term deposit in (A), like in light-harvesting systems. We propose some efficient and robust external control schemes to achieve stable optimal final population in (A). Facing an indirect population transfer process, our control observable should incorporate both a maximum depletion of the ground state and a maximum final population in (A), when sharing between (D) and (A). This is done by defining an asymptotic contrast in terms of the ratio of the population difference between (A) and (D) to the remaining ground state population. In order to maximize this contrast, and inspired by a thorough understanding of the 2D (CI branching space) post-pulse field-free population evolution, we refer to two basic coherent control mechanisms, namely, pump-pump interference and kicks. The search for the optimal control parameters for an electromagnetic field in terms of a train of ultra-short laser pulses is conducted on a reduced but still realistic 4D model of pyrazine. More precisely, we optimize the pulse leading frequency and intensity, duration, and the inter-pulse delay. Once obtained and fixed, these parameters are used in the full 24D dynamics to calculate the time evolution of the electronic states populations. It is very gratifying to see that this extension of control mechanisms, often established for over simplified systems involving only a small number of levels, can survive and be effective when applied to systems with a large number of degrees of freedom. This is probably the most relevant message of this study.

Having in mind the limitations of realistically achievable experimental conditions, the excited states population contrasts we obtained can be considered as efficient ones for pyrazine CI-mediated funneling dynamics. Referring to either interference or kick mechanisms, about 60% of the ground state population is deposited in the acceptor state, while about 16% remains in the donor state. As a perspective, we can mention other control strategies that could be checked in terms of their possible extension to multi-dimensional systems. One is the strong field pump-probe interference we already mentioned [cf. Eq. (27)]; the other being the frequency chirp. The density of levels for the vibrational baths being important, addressing in an optimal way some specific levels playing a major role in the funneling dynamics control may deserve interest. Even more importantly, pyrazine should be considered as an example of a modestly complex system. As a conclusion and consequence of our findings, we anticipate the possibility of potential transposition of such control mechanisms to other larger biological systems. This would certainly be conditioned by having a small number of active normal modes among all the others, which are spectators in the dynamics we wish to control. Then, similar to what we have observed for pyrazine, external control fields once defined and fixed in low-dimensional models should be effective at higher dimensions.