Quantum control of electron localization in molecules driven by trains of half-cycle pulses

We present numerical simulations demonstrating efficient control of electron localization dynamics in small molecular systems by a train of half-cycle pulses using the H2+ molecule as a prototypical model system. Simulations are performed using a direct numerical integration of the Schrödinger equation on a grid for a two-degree of freedom system (one electronic, one nuclear coordinate) as well as for an expansion in terms of Born–Oppenheimer basis states. By varying the parameters of the driving field, we systematically explore the underlying control mechanism of this periodically approximately impulsively driven molecular system and demonstrate the robustness of the proposed control scheme.


Introduction
For more than two decades there has been an increasing desire to develop and implement protocols to actively control quantum systems, with possible applications reaching from physics and chemistry to biology [1,2]. Many different control schemes have been proposed and successfully applied to quantum and classical systems; for reviews see, e.g. [3]- [7]. Very often, control protocols involve electromagnetic pulses that are tailored for the weak-field regime and in this limit they may aid understanding of the dynamics of the unperturbed system (e.g. [8,9]). Other control schemes invoke strong fields, such as stimulated Raman adiabatic passage [10] and its relatives, Stark-control [11,12], selective population of dressed states [13] and optimization of the chirp [14]- [16]. As the dynamics in strong fields is often very complicated, an analysis of why a given scheme works may become increasingly difficult.
For ultrashort pulses on the order of a few femtoseconds, carrier-envelope phase (CEP) effects start to become important. Changing the CEP of a few-cycle pulse allows for an efficient control of molecular dynamics [17]- [19]. Recently, experimental control of electron localization dynamics in the dissociative ionization of D 2 has been demonstrated [20]. By changing the CEP, the atomic constituent of the dissociating D + 2 molecular ion on which the electron localizes could be controlled. The influence of other control parameters has not yet been examined systematically. Two-pulse setups, involving mainly a combination of a few-cycle infrared (IR) and a short extreme ultraviolet (XUV) pulse [18,21,22], introduce the delay time τ between the pulses as an additional control parameter.
An alternative tool for quantum control is (trains of) half-cycle pulses (HCPs). Simulations have shown that HCP of variable width are very efficient to control dissociation of OH [23] and of HDO [24], to control dynamics of isomerization reactions [17,25], and for orientation of rotational wave packets [26,27]. Trains of HCPs have also been applied, both theoretically as well as experimentally, to control the dynamics of Rydberg wave packets in atoms [28]- [31] and three-body Coulombic systems [32].
In this paper, we study trains of HCPs to control the coupled nuclear-electronic dynamics in the dissociation of small molecular systems. One salient feature of the HCPs for quantum control is their approximate unidirectionality. In the following, we exploit this breakdown of inversion symmetry to localize the electron at a preselected target site. We will demonstrate 3 that this mechanism of control is also quite robust to small changes in the pulse parameters. The HCPs we are using in the simulation can be experimentally synthesized and the suggested quantum control protocol should thus be realizable.
We consider HCPs (frequently named 'kicks') each with a duration of ≈800 attoseconds. The initial nuclear wave-packet mimics the ionization of H 2 by an ultra-short XUV pulse and evolves on the electronic ground state of H + 2 . The subsequent interaction with the pulse train excites the higher lying electronic state with ungerade symmetry, which is dissociative. During the interaction with the pulse train, a coherent superposition of electronic states is created, and the molecule dissociates. In the double-well picture, the wells (Coulombic potentials of H + ) separate and the barrier in between rises, eventually trapping the electron in one of the two wells. The role of the HCP train is now to drive the electron into the desired well. This process resembles the experiment on localization dynamics by Kling et al [20], but with the substantial difference that the peak fields we are using are much weaker, the fields are longer, and the driving by HCPs is approximately unidirectional. Another important difference is that in the present scenario, the ionization H 2 → H + 2 is induced by an XUV pulse rather than the control (HCP) pulse thereby resembling the two-pulse setup by He et al [21].
We highlight the critical role of the field strength of the HCP train in this control scheme. Depending on the field strength of the unidirectional pulse train, either strong localization in the direction of the 'kicks' or much weaker localization in the opposite direction is observed. We discuss the origin of this somewhat counterintuitive behavior.
The paper is organized as follows: section 2 describes briefly the model system of H + 2 we are using. In section 3, the numerical results for the exploration of the control landscape of electronic dynamics is studied. The physical mechanism underlying the localization is explored in section 4 by means of a simplified model. Finally, section 5 contains a brief summary and conclusions.

Model for H
As a model system, we consider the coupled nuclear-electronic dynamics of the hydrogen molecular ion, H + 2 . This choice is motivated by the recent experimental and theoretical studies [19,20,22,33,34,35,36] using the CEP of a few-cycle IR pulse as a control parameter and by the calculations for a two-pulse excitation scheme considering the time-delay between pulses as the control parameter [21,22], both for H + 2 and its isomers. We treat this problem in reduced dimension with one electronic degree of freedom (−∞ < x < ∞) and one nuclear degree of freedom, the internuclear distance (0 R < ∞). Accordingly, the model Hamiltonian is given bŷ The smoothing parameter ε eN (R) in the soft-core potentials is a function of the internuclear distance R in order to yield both the correct equilibrium distance (here R eq = 1.25 Å) and good approximations for the Born-Oppenheimer (BO) potential curves [37]. The full potential V (x, R) and representative cuts V (x, R) for fixed internuclear distances R are displayed in figure 1.
The interaction with the driving field F(t) of the pulse is treated in the length gauge of the dipole approximation. Certainly, this two-degree-of-freedom model can represent only a limited class of physical processes. For example, the dependence on the relative orientation between the internuclear axis and the field axis are beyond the scope of this model. However, several strong-field phenomena for which the polarization along the internuclear axis dominates the dynamics can be reasonably well represented by equations (1) and (2).
The time-dependent Schrödinger equation for the coupled electronic and nuclear coordinate reads as with ϕ g (x; R) the gerade electronic eigenstate and χ 0 (R) the initial nuclear wave packet. Our analysis is based on calculations performed on a two-dimensional grid. We use the split-operator technique [38], with a total number of 1024 points in R and in x, where the grid in R is defined from 0.1 to 64 a.u. and the corresponding grid in x from −100 to 100 a.u. The eigenstates are obtained by imaginary time propagation [39]. An optical potential cos 2 (R) absorbs the escaping wave packet 3 a.u. before the grid end, avoiding grid reflections. Note that with the large box used in the simulations, virtually the whole wave packet is staying on the grid for the times considered here. Alternatively, we have also performed calculations using an expansion in terms of the BO basis spanned up by the gerade and ungerade states. These calculations follow the nuclear Schrödinger equation: with the nuclear Hamiltonian The highest lying curve is the ionization potential curve of H 2+ 2 . Also shown is the initial nuclear wave packet.
In the above equation, m P is the proton mass, I is the unit matrix, and µ gu (R) is the transition dipole between the ground and upper electronic states. The electronic eigenenergies depend parametrically on R, yielding the potential curves 2 + g = V g (R) and 2 + u = V u (R), see figure 2. We use this simplified approach with the same grid in R which neglects the higher lying electronic states when we perform time-consuming scans in the multi-dimensional control parameter space. Its accuracy is tested by the comparison with the full two-dimensional solution.
To compare the BO results with the results from the full calculation, we expand the latter as Above, the ϕ i (x; R) are the electronic eigenstates, and the χ i (R) are the nuclear wave functions. The time-evolved wave packets (x, R, t) or χ(R, t) are characterized by several observables, in particular the populations ||χ g (t)|| = χ g (t)|χ g (t) and ||χ u || of the gerade and ungerade states (and correspondingly ||g(t)|| = ϕ g | (t) and ||u(t)||) and the expectation values of the bond length R(t) = (t)|R| (t) . Of particular interest in the present context of localization dynamics are projections onto coherent superpositions which converge for large R to the localized atomic wavefunctions in the right or left potential well, i.e. the two-dissociation limits H + + H and H + H + , see also figure 1.

Half-cycle pulses
An HCP features a time-dependent unipolar electric field of short duration and a single field maximum. Such pulses can be generated by conventional pulse generators on the (sub) nanosecond scales [28,40,41]. For much shorter pulses that represent propagating solutions of Maxwell's equation, the time integral over the pulse F(t)dt = 0 must vanish and truly unipolar pulses are ruled out. Nevertheless, approximately unipolar pulses can be generated with a strong peak field in one direction, F + (t), accompanied by a low-amplitude long-lasting off-set field in the opposite direction, F − (t). Employing semiconductor-based techniques, quasi-HCPs with durations as short as 500 fs [42] have been achieved. More recently, techniques based on harmonic generation have been suggested [43,44], promising the creation of HCP trains down to the attosecond regime. Such a HCP train can be described by the Fourier decomposition of the harmonics where f (t) is the normalized envelope, F 0 the overall strength and F n the amplitude of the harmonics with phase φ n . Thus an HCP train can be formed by mixing the fundamental frequency with several higher harmonics. To shape a pulse with preferred directionality a mixture of even and odd harmonics is necessary. Furthermore, for a nearly unipolar pulse train, the phases of the harmonics should be aligned as Here, k is an integer, and φ is related to a time delay, t = − φ/ω. This condition can be met by an appropriate combination of filters. To minimize fluctuations of the peak height, the amplitudes of the harmonics should ideally follow a Gaussian distribution, F n = exp(−(n/n 0 ) 2 ). Very recently, trains of HCPs have indeed been generated experimentally, (see e.g. [45]- [48]). A typical HCP trained used in our simulation is shown in figure 3. The unipolar strong field pulses in the positive direction F + are off-set by a weak-field pointing intermittently in the negative x direction, with strength The latter is an immediate consequence of the constraint F(t)dt = 0 for a propagating Maxwell field. For all our simulations, the width (FWHM) of the individual HCPs is 0.2 T p with T p the period of the train. The fields we consider are of moderate intensity (<5 × 10 13 W cm −2 ), ensuring that no significant amount of ionization occurs and a two-level description is applicable to a good degree of approximation. In none of the cases considered, we observe substantial ionization (< 5%) or population of higher lying electronic states, as checked by the full calculation for the highest intensity.

Numerical results
Starting point of our simulation is a vibrational wave packet (R, t = 0) in the electronic ground state of H + 2 formed, e.g. by Franck-Condon transition from the H 2 ground state after photoionization by an ultrashort XUV pulse. This initial wave packet is subject to the train of HCPs, which we optimize using a genetic algorithm to maximize the localization of the Optimal HCP train found by the genetic algorithm. The pulse is characterized by the amplitude F 0 , by the delay τ of the center of the HCP relative to the exciting XUV pulse at t = 0, the length of the train T (full width at half maximum (FWHM)), and the pulse CEP φ CEP (in this case the shift of the largest HCP relative to the peak of the envelope). Also shown is the envelope function f (t), (equation (9)). electronic wave packet, i.e. maximizing the asymmetry of the dissociation fragment distribution. We present results obtained from both the two-level BO expansion (equation (4)) and the solution of the fully coupled electronic-nuclear system. For computational efficiency reasons, the optimization of the control parameters are performed within the framework of the BO expansion.

Quantum control of the asymmetry in electron localization
We search for a pulse train that induces the highest possible degree of electron localization in the left potential well, corresponding to the direction of the force from the HCPs. The parameters of the HCP train modified in the genetic algorithm are field strength F 0 , F n , the CEP of the generating pulse φ CEP , the pulse length T , and the time delay τ of a driving pulse train relative to the ionizing XUV pulse (see figure 3). We thereby restrict the parameter space for field strength to F 0 < 0.04, avoiding ionization and pulse duration T <50 fs, keeping rotational dynamics frozen during the interaction. The envelope of the driving pulse was assumed to be of Gaussian shape and the maximum number of colors (harmonics) included in the HCP train was fixed to six. The fitness function for localization is defined as the asymmetry, The asymmetry A can take values between −1 and +1. The sign of A indicates in which well the electron is preferentially localized. Note that for the positive maximum of the field F, the classically preferred localization would be on the left potential well, i.e. A > 0.  The optimal pulse train found after less than 50 generations has a field strength of F 0 = 0.029 a.u., a length of T = 20 fs, centered at τ = 18 fs and a CEP of φ CEP = 0.12π , see figure 3. All of these values are well inside the parameter domain used for the genetic algorithm. The dynamics from the solution of the full coupled nuclear electronic Schrödinger equation starting with an initial nuclear wave packet in the electronic ground state, (x, R, 0) = ϕ g (x; R)χ 0 (R), R = 0.7 Å and driven by the optimal pulse train, features a partial population transfer to the upper electronic state, see figure 4. With increasing R, the population transfer becomes more efficient due to (a) an increased transition dipole moment µ gu (R), and (b) a reduced potential gap V (R) between the BO potential curves. The population transfer is very efficient with more than 40% of the population being transferred to the excited state. Population of the higher lying excited states is very small (<3%) and also the ionization is small (< 6%) justifying the twostate BO expansion (equation (4) figure 4, however, for a lower field strength F 0 = 0.018. Note that the localization is entirely different with preference for R rather than L .
pushes the electron to the left side. When one kick is over, and the field has a negative sign (F − (t)), the electron oscillates back. Upon the next kicks, this scenario is repeated. After about 12 fs, the degree of localization is frozen, and the following peaks-although having higher field strength-just keep the electron on this well until the dissociation limit is reached (>15 fs). A significant portion of the localization by 'kicking' occurs when either the excited state or both states are still above the barrier separating the double well (figure 1). Only after the nuclear wave packet has reached a distance of ≈ 5 Å, both the BO states are located below the barrier and are energetically degenerate. The electronic motion then freezes out and the electron is irreversibly localized in one potential well. The HCP-induced localization dynamics is remarkably efficient reaching A ≈ 0.85 in the optimal setting. The intuitive simple 'classical' picture of an electron being 'pushed' in one direction by the HCP sequence is somewhat deceiving. To illustrate the subtlety of the underlying quantum control mechanism we repeat the simulation, setting, however, 'by hand' the field strength of the train F 0 to be about half of the 'optimized' pulse train ( figure 5). Although the state L shows intermittently prominent peaks at early times, little localization results at the end of the pulse, which even prefers the opposite side (A = −0.27). Until around 10 fs, both cases behave qualitatively similarly; however, taking a closer look at t = 8 fs, we see that the kicks around this time seems to be of special importance. This observation will be explained in section 4. This counterintuitive result demonstrates the quantum properties of quantum control in this two-state system where quantum-classical correspondence is, a priori, not expected to hold. These figures illustrate that the relative phase between the electronic ground-state and the excited-state wave packets is key in the ensuing localization dynamics. The calculations show that for the highest degree of asymmetry possible, control over both the amplitudes and relative phases of the wave packet must be achieved. The profound differences between the different strengths of the driving field are highlighted in the evolution of the density of the electronic wave packet, as shown in figure 6. The ray towards negative x describes the localization to the left, i.e. L (dissociation H + H + ), whereas the ray towards positive x corresponds to R (or dissociation H + + H). The entirely different outcome in both directional emission of the dissociation products and the degree of dissociation for the two different field strengths F 0 is clearly visible. Also the oscillations in the density correlated with the temporal distribution of the 'kicks' is clearly recognizable.

Dependence of the asymmetry on the pulse parameters
In the preceding section, we highlighted the strong dependence of the localization asymmetry A (equation (12)) on the magnitude of the field in the HCP train. Before exploring the underlying physical processes, we explore the quantum control landscape [49,50] and systematically study the dependence on other control parameters. The following simulations scanning the different parameters are performed within the BO approach (following equation (4)). As a reference, we take the pulse shown before, with φ CEP = 0.123, F 0 = 0.0290 a.u. and λ = 800 nm, with a width T of 20 fs, containing six harmonics, and centered around τ = 18 fs relative to the ionizing XUV pulse. As a first step, we scan the field strength of the HCP train ( figure 7). As expected, the field strength influences the asymmetry A dramatically. The variation is, however, oscillatory rather than monotonic, pointing to quantum phase interference effects. Variation of the length of the pulse train T and of the delay relative to the XUV pulse τ (with ( T, τ ) = (0, 0) corresponding to the optimal pulse) displays pronounced variations in both τ and T (figure 8). Either variation changes the temporal position of the kicks relative to the dissociation dynamics induced by the ionizing XUV pulse initiating the vibrational wave packet propagation. Note that this dependence on T and τ results from the intrinsic timescales present determined by the time at which the nuclear wave packet reaches the region with R ≈ 5 Å where the dynamics freezes out. While the two-dimensional cut through the control landscape (figure 8) displays a clear maximum, its width is sufficiently broad such that the optimization algorithm yields a stable and robust maximum under realistic experimental conditions. Another important two-dimensional cut through the multi-dimensional control landscape is the (F 0 , φ CEP ) plane (figure 9). While the dependence of A on F 0 is pronounced, the CEP dependence is weak, attesting to the robustness of the quantum control protocol. Figure 9 suggests that the pronounced asymmetry can even be observed for HCP trains without stabilization of the CEP of the fundamental driving field. (Clearly the relative phases φ n still need to be well defined.) Comparison between figures 9(a) and (b) displays the dependence of the optimized peak field F 0 on the fundamental wavelength λ of the driver. This dependence emphasizes the role of the 'kicks', i.e. the momentum transfer during a single unidirectional pulse suggesting that the optimum asymmetry should appear at constant F λ and, in turn, F 0 λ, as is indeed observed. As generation of a large set of phase-aligned harmonics is technically challenging, we have also checked how well pronounced the character of a HCP train needs to be, i.e. how a few harmonics are sufficient in order to produce significant asymmetries ( figure 10). Obviously, the requirements are quite modest. Mixing only the fundamental frequency and its third harmonics leads to an asymmetry A oscillating around A = 0, with an efficiency comparable to the onecolor case only. As long as low-order even and odd harmonics are both present and the inversion symmetry of the pulse is clearly broken, pronounced asymmetries appear. Adding up to three harmonics leads already to well-converged results.

Model for electron localization
In order to unravel the underlying mechanism responsible for localization identified by the genetic algorithm, we resort to a simplified model that couples the expectation value of the internuclear distance R with the time t, R → R(t). Clearly, spreading, beats and revivals of the nuclear wave packet are thereby eliminated. Accordingly, we study the electronic dynamics parametrically dependent on the bond length expectation value. We choose R(t) as the internuclear distance expectation value (figure 3(c)) and consider the parametric dependence of the BO potential curves V g,u (R(t)) and the dipole matrix elements µ(R(t)) on R(t). The Schrödinger equation simplifies to . We have confirmed that the simplified model reproduces the asymmetry found in the exact calculation quite well, see figure 11. Within this model, the electronic evolution can be viewed as an adiabatic propagation in the slowly varying background field F − (t) together with non-adiabatic transitions caused by the HCPs.
In the absence of the periodic kicks, the electronic wavefunction would evolve quasiadiabatically along the field-dressed eigenstates φ 1,2 , i.e. eigenstates of the Hamiltonian matrix equation (15) with F − (t) instead of F(t), see figure 12(a). The adiabatic state φ 1 smoothly  (15)). In (a), we use the optimal HCP field (see figure 4) and in (b) the lower field strength ( figure 5).
connects the gerade state ϕ g at small R (when µ(R(t))F(t) is small compared to V (R)) with state R (corresponding to H + + H), the lower lying state in the field. Analogously, the fielddressed state φ 2 adiabatically connects ϕ u with L . The HCPs drive non-adiabatic transitions between the two-field dressed states φ 1 and φ 2 . Its strength is controlled by the transition amplitude ('boost operator') where p is the momentum transfer imparted by a single HCP as defined in equation (14). For small t (or R), the transition amplitude B 12 increases (see figure 12(b)) since the field strength increases. As we calculate B 1,2 between the dressed states (and not between the BO states) it vanishes again for large times. Here, the adiabatic eigenstates 1,2 are already closely resembling polarized atomic orbitals residing in separated potential wells. Further note the interference minimum at t = 10 fs, which corresponds to px = π. For effectively driving the transition, the off-diagonal elements in equation (15) should match the energy difference between the diagonal elements. The parametric dependences of |B 1,2 |, µF and V on t (and in turn on R(t)) show that the transfer between 1 and 2 for the optimal field is most effective near t ≈ 8 fs (see figures 12(b) and (c)), as is indeed seen in figure 11(a).
Clearly, such an impulsively driven two-state system does not necessarily comply with the classical intuition, which would suggest a monotonic dependence of the degree of electronic polarization on the strength of the unipolar driving field. The optimal pulse shape parameters found within the quantum control protocol can thus be viewed as the optimal placement in time (or internuclear distance) of the 'decisive' kick such that the transition amplitude from 1 to 2 is optimized. The subsidiary constraint is that subsequent kicks of the sequence come sufficiently late such that they are already ineffective to transfer the state back to 1 . The wave packet will then continue to adiabatically evolve on 2 leading to the polarized state L . For the suboptimal case, the last condition is not fulfilled. Here, the decrease of the transition amplitude starts at later t, i.e. px increases more slowly with t. Hence, two successive HCPs cause Transition amplitude B 1,2 , equation (16) corresponding to F − for the optimal field and the weaker field with F 0 = 0.018. (c) Matrix elements V (R(t)) and µ(R(t))F(t) used in equation (15).
strong transfer between 1 and 2 ( figure 11(b)). The second transfer cancels the effect of the first one, leading to zero net transfer between the adiabatic states. This is in agreement with the 'counterintuitive' result of localization in the direction opposite to the HCPs found in the full simulation.

Summary and conclusion
We have presented numerical simulations demonstrating the very efficient control of the electron localization dynamics in small molecular systems by applying a train of HCPs. Using H + 2 as an example, we find high degrees of electron localization during the dissociation in one potential well (≈ 85%). Compared to other control mechanisms recently proposed and realized in H + 2 to localize the electron in one potential well, the present approach is very robust to changing parameters. Parameters such as the CEP serve only for fine tuning of the absolute degree of localization. The efficiency and even direction of the final localization is, however, strongly and non-monotonically dependent on the field strength of the HCP train. We study a simplified two-level model considering the electronic states at the bond-length expectation value only and explain the latter finding as caused by a combination of adiabatic motion in the smooth off-set field opposite to the HCPs and non-adiabatic impulsive couplings by certain HCPs.
A low-order approximation of the HCPs consisting of the first and second harmonic frequencies of the fundamental driving frequency only is sufficient to achieve a substantial degree of asymmetry. This shows that the unidirectionality-the breaking of the inversion symmetry of the driving field-is the key to the electron localization by a half-cycle-pulse. An experimental verification of the localization including at least two colors would be highly desirable.
The results in the present contribution were obtained in a model including one electronic and one nuclear degree of freedom. On the timescales we are dealing with (t <50 fs), the rotational motion can be neglected. Since we are interested in the dynamics along the molecular axis using a linearly polarized field, we can extract the main features from a one-dimensional description. A full three-dimensional picture might change the absolute asymmetry but will not change the basic mechanism we are studying.
In future, we plan to apply HCP trains to larger chain-like molecules and develop coherent control schemes for electron dynamics along such chains.