Non-Markovianity in the optimal control of an open quantum system described by hierarchical equations of motion

Optimal control theory is implemented with fully converged hierarchical equations of motion (HEOM) describing the time evolution of an open system density matrix strongly coupled to the bath in a spin-boson model. The populations of the two-level sub-system are taken as control objectives; namely, their revivals or exchange when switching off the field. We, in parallel, analyze how the optimal electric field consequently modifies the information back flow from the environment through different non-Markovian witnesses. Although the control field has a dipole interaction with the central sub-system only, its indirect influence on the bath collective mode dynamics is probed through HEOM auxiliary matrices, revealing a strong correlation between control and dissipation during a non-Markovian process. A heterojunction is taken as an illustrative example for modeling in a realistic way the two-level sub-system parameters and its spectral density function leading to a non-perturbative strong coupling regime with the bath. Although, due to strong system-bath couplings, control performances remain rather modest, the most important result is a noticeable increase of the non-Markovian bath response induced by the optimally driven processes.


I. INTRODUCTION
Open quantum systems (OQS) are ubiquitous in physics and chemistry and have many uses from setting quantum technology in condensed phase to exploring long-lived coherence in biological systems [1][2][3][4][5][6]. They consist in selecting a given partitioning into a central quantum system and a statistical surrounding bath. The reduced system dynamics is nonunitary and can be called Markovian or non-Markovian according to the importance of memory effects [2]. The comparison of system and bath typical timescales is a relevant qualitative measure to separate both situations : if the timescale characterizing the bath is shorter than the one of the system, dynamics can be said Markovian, non-Markovian if not. For a two-level system, this characteristic time is the Rabi period whereas the bath dynamics can be estimated from the time decay of the two-time correlation function of the system bath coupling related to the Fourier transform of the bath spectral density. A nearly delta correlated bath leads to a Markovian behavior usually described by Lindblad [7] or Redfield [2,5,8] approaches involving unidirectional relaxation. Non-Markovianity is described by strong quantum memory effects leading to temporary information back flow from the environment to the system. Several measures of non-Markovianity have been proposed and compared recently in the literature [3,4]. Among them one can mention the distinguishability of quantum states estimated by their trace distance that can transitively decrease during the relaxation, as opposite to a Markovian evolution in which it continuously increases [9,10]. Other non-Markovianity signatures refer to a re-amplification of the volume of accessible states during the decay process [11], the detection of a negative canonical decay rate [12,13], or a non-monotonous time evolution of the system von Neumann entropy [14].
Our main purpose is to take advantage of the back flow of information from the surrounding bath, characterizing non-Markovianity, to enforce the control of the central system physical observables, protecting them against decoherence. At that respect, the present paper is a second one of a series of three [29,30] where an optimal control scheme is worked out, still acting on the central system alone, aiming at some protection against decoherence (population revivals, or robust and efficient transfers) and subsequently examine its consequences in terms of the bath non-Markovian response. More precisely, we analyze non-Markovianity during an ultra-short field pulse optimized by quantum control [32][33][34] in a spin-boson (SB) model [1,2,35] where the active sub-system strongly interacts with the bath. The controlled dynamics ends before the complete decay of the volume of accessible states in the Bloch sphere [11], i.e. before the decay of the bath correlation function which means before quantum memory (or non-Markovian) effects are expected to vanish. The control is also shorter than the full relaxation time of the state populations towards equilibrium. The interaction of the two-level system with the bath is described by the standard spin boson Hamiltonian (SB) which can used in many different situations ranging from qubit in quantum dots to exciton or charge transfer. In the present work, it is built and calibrated to simulate a charge transfer between donor and acceptor electronic states in a heterojunction [36,37]. The model addresses ultra-short control of electronic dynamics in a complex system strongly coupled to the nuclear vibrational motion [5]. Similar coherent control of excitation energy transfers in photosynthetic systems has already been investigated, but in weak coupling regimes, referring to Markovian approaches [39,40]. Here we analyze a nonperturbative situation, described through hierarchical equations of motion (HEOM) [41][42][43][44][45][46].
We focus on early dynamics and we investigate the extent to which optimal control field enhances non-Markovianity during control. The canonical decoherence rates and the von Neumann entropy are taken as signatures of non-Markovianity. In a recent work, the enhancement of non-Markovianity during laser driven dynamics has been studied with simple periodic fields in a SB model with a smooth Lorentzian spectral density [25]. This example shows an enhancement of non-Markovianity signatures but for weak coupling only. On the contrary, in the present work, we obtain non Markovian behaviours even in the strongly coupled case.
Optimal control theory (OCT) is implemented here together with the HEOM method. The Rabitz monotonous algorithm in Liouville space we are referring to [47][48][49][50], requires the forward and the backward propagations of the master equation. The memory kernel occurring in a time non-local master equation with a final condition has been discussed in different works. It has been implemented at second order level keeping the memory kernel [48,51,52] and by the auxiliary matrix method leading to time local coupled equations [50,53]. We generalize here this methodology with HEOM equations at higher order. The HEOM master equation can be rewritten as a time dependent Lindblad superoperator with time dependent canonical rates to get a witness of non-Markovianity [12,13]. This interesting Krauss decomposition [54,55] has already been suggested to analyze the control in ref. [20]. In a first attempt, we do not impose any constraint on the field area so that the optimal field is not necessarily an optical one with zero area [56][57][58][59]. Such constraint could be added in a second step, but this issue would go beyond the scope of this paper. The electric field is assumed to have a dipole interaction with the central system only. However, since the memory kernel depends on the external field through the system Hamiltonian, this latter has an influence on the bath dynamics so that control and dissipation are strongly correlated. The modification of the bath dynamics is probed here from the HEOM formalism by analyzing the first moment of the bath collective mode [61].
The paper is organized as follows. Section II describes the SB model calibrated from data simulating a charge transfer in a heterojunction. The HEOM equations, the signatures of non-Markovianity and the optimal control theory in dissipative system are presented in Section III. Section IV gives the results for three ultra-short control cases, two for which the target is the initial state itself (a revival), and one for which the control enforces a transition between the two levels. Finally, some perspectives are presented in Section V.

II. THE MODEL
The spin-boson model is a two-level quantum system linearly coupled to a bosonic bath of harmonic oscillators at thermal equilibrium. The Hamiltonian reads in mass weighted coordinates and H SB = S k c k q k . Atomic units are used with = 1. The system operator is S = σ z with σ i operators taken as Pauli matrices. The control field E(t) only acts on the two-level system and is assumed to be linearly polarized. The µ matrix is the matrix of the corresponding component of the dipole operator. In the context of a charge transfer between a donor and an acceptor in a heterojunction, H S (t = 0) corresponds to the diabatic representation of the two electronic states for which the parameters are estimated at the equilibrium geometry. The diabatic parameters δ and W are taken from a model heterojunction between oligothiophene and fullerene [36,37]. The inter fragment distance is fixed to R = 3Å leading to δ = 0.21 eV and W = 0.13 eV. The corresponding Rabi period is 12.3 fs and the eigenenergy gap is 0.33 eV. The dipole matrices are not calibrated from ab initio calculations and different dipole models have been used to discuss the stability of the observed behaviors. In this electron transfer framework, the bath is formed by all the normal modes of the two fragments (here 264). The harmonic frequencies are assumed to be the same in both electronic states but the equilibrium geometries differ by a distance d k .
Taking the origin of bath coordinates at a middle position between these equilibrium points, the vibronic coupling coefficients are c k = ω k 2 d k /2 .
The bath is fully characterized by the spectral density leading to the two-time correlation function where is the Boltzmann equilibrium density matrix of the bath and β = 1/k B T . Spectral density and correlation functions (real, imaginary parts and modulus) of this heterojunction model are displayed in Fig.1. In this example, the Rabi period (12.3 fs) is smaller than the correlation time (25 fs) so that non-Markovian dynamics is expected.
As displayed in Fig.1, the spectral density J(ω) is fitted by four four-pole functions where Λ l, (1,2) Cauchy's residue theorem is used to compute the integral of Eq. and an infinity of poles on the imaginary axis ∀j ∈ N * / 0, ν j = 2π β j called the Matsubara frequencies. In practice, the number of Matsubara terms is limited ensuring convergence for a given temperature.

A. HEOM equations
The system density matrix is the partial trace of the full density matrix Ξ(t) over the bath degrees of freedom ρ(t) = T r B [Ξ(t)]. The initial condition is assumed to be factorized The hierarchical equations of motion have been established from the path integral method [44] or from the stochastic Liouville equation [41][42][43]. The non-Markovian master equationρ is solved by a time local system of coupled equations among auxiliary matrices arranged in a hierarchical structure. The algorithm requires a particular parametrization of the correlation function as a sum of n cor exponential terms, written as: Analytical expressions for the α k and γ k parameters can be derived when the spectral density is fitted by a sum of two-poles [60] or four-pole Lorentzian functions leading to an Ohmic or super Ohmic behavior at low frequencies [38]. The complex conjugate of the correlation function can be expressed by keeping the same coefficients γ k in the exponential functions but using modified coefficientsα k according to: k being a collective index such that,α l,1 = α * l,2 ,α l,2 = α * l,1 ,α l,3 = α * l,4 ,α l,4 = α * l,3 and α j,matsu = α j,matsu , where α l,m ,α l,m with m = 1, 4 are related to the four poles of each Lorentzian l [62].
The level L of the hierarchy corresponds to an order 2L in the perturbation expansion of the initial non-Markovian equation. Auxiliary matrices are labeled by a collective index n = {n 1 , · · · , n ncor } specifying the number of occupation of each artificial mode associated with one of n cor decaying components. The system density matrix ρ(t) has the index n = {0, · · · , 0}. The first level L = 1 contains n cor auxiliary matrices with a single excitation only k n k = 1. The HEOM coupled differential equations are given by : with n + k = {n 1 , · · · , n k + 1, . . . , n ncor } and n − k = {n 1 , · · · , n k − 1, . . . , n ncor }. Each matrix is coupled only to the superior and inferior levels in the hierarchy. The level of the hierarchy is chosen until convergence is reached for the system density matrix.
The HEOM formalism allows one to get insight into the correlated system-bath dynamics by probing the different moments [61]. In particular, the expectation value of B in each state is given by the diagonal elements of the X (1) (t) operator given by the sum of the first level auxiliary matrices where the sum runs over all index vectors n = {n 1 , · · · , n ncor } with l n l = 1. Recursive formula for higher orders can be found in ref. [61]. This first moment already provides a signature of the induced correlated system-bath dynamics. As discussed in [61], the master equation can be recast to emphasize the role of X (1) (t) in the system dynamics by writinġ B. Non-Markovian witnesses.
Signature of non-Markovianity is discussed here through the volume of the accessible states [11] and through the canonical decoherence rates of a time-dependent Lindblad form [12,13]. In the two-level case, the dynamical map The volume of accessible states may be obtained from the matrix representation of the This volume may also be expressed as a function of the decoherence canonical rates. The master equation is then recast in a canonical Lindblad form but with time dependent rates associated with time-dependent decay channels. Details can be found in refs [12,13]. The master equation is reformulated aṡ In order to describe the decrease of the Bloch volume independently of the translation of its center, the contribution of the unity operator is separated by gathering terms containing coefficients a j0 . One then defines an operator O = a 00 /2d + where D jk (t) is the decoherence matrix. Its diagonalization provides the decoherence canonical rates g k (t) and the decay channels C k (t). Eq. (15) becomeṡ It is worthwhile noting that the occurrence of negative canonical decoherence rates g k (t) yields another characterization of non-Markovianity [12]. The rates are linked to the volume of accessible states through the relation with The criterion based on the volume can be considered as an average measure since it depends of the sum of the rates only. Thus, it can be considered as a less stringent witness of non-Markovianity than a negative canonical decoherence rate g k (t).
A possible numerical strategy to compute the decoherence matrix D ij (t) has been discussed in [12] and is given by with Besides the analysis of the decoherence canonical rates, we also compute the von Neumann entropy of the system that should vary monotonously in a Markovian evolution [14] S where λ k are the eigenvalues of the system density matrix.

C. Optimal Control Theory
We use optimal control theory in the Liouville space [47,48,50] to optimize the field driven state-to-state transfer at the end of the pulse of total duration t f . The cost functional is built from a chosen performance index from state to state at time t f , here T r ρ † (t f )ρ target with a contraint on the field intensity and one assuring the respect of the mater equation at any time. The corresponding Lagrange multipliers are the scalar α 0 and the density matrix χ(t) respectively. We do not enforce here the constraint on the pulse area which is required for a purely optical field [59]. The optimal field is obtained from the system matrix density propagated by the master equation with initial condition ρ(t = 0) = ρ ini and from the Lagrange multiplier propagated with a final condition χ(t = t f ) = ρ target .
The corresponding master equations with initial and final conditions take the form with When the master equation is solved by the HEOM algorithm, the operational equations for the Lagrange multiplier can be derived by using Eqs. (7) and (8) In practice Eq. (24) is solved backwards starting from χ {0,0,..,0} (t = t f ) = ρ target with all the auxiliary matrices set equal to zero. The field at iteration k is obtained by

A. Field-controlled dynamics
We consider three control objectives defined by the populations of the system. In the two first strategies that are denoted C1-1 and C2-2, the target is the revival of initial zeroorder state, either state 1 or 2, at the end of the control. A third control denoted C1-2, enforces the fast decay from state 1 to state 2 (a fast switch from 1 to 2). We compare the control without or with dissipation and analyze both the system and bath responses (memory effects) during the corresponding field-driven dynamics.
The field free and field driven populations in the initial state during the three control strate- Obviously, control performances remain rather modest. This point can be explained both by limitations of the control parameters (rather low field amplitudes and short pulse duration), and more importantly, by the way the strong system-bath interactions inherent to the specific molecular situation at hand interplays with the control. This can be numerically rationalized through the analysis of the first moment of the bath collective mode in each state given by the diagonal elements of the 2 × 2 X (1) matrix as depicted in Eq. (10). Although the control field is explicitly introduced only to act on the system Hamiltonian, it affects the overall dynamics through the memory terms included in the right-hand-side of Eq.(9).
Control fields indirectly modify the bath response leading to a strong correlation between control and dissipation. This is illustrated in Fig.3 which displays the first moment of the bath collective mode, in terms of the diagonal elements X 2,2 , starting from the time when the gap is at its maximum value, i.e., close to 6 fs for C1-1 and 4 fs for C1-2 control strategies. Actually, when dealing with these two strategies, the gap is decreasing during the first femtoseconds, such that the system internal transition frequencies better match bath resonant phonons transitions. As a consequence, the amplitudes of collective modes oscillations are expected to increase. For initial state 2 and the corresponding C2-2 control strategy, early Stark shifts have an opposite sign leading to increasing gaps, preventing bath resonant processes from occuring. Discrepancy from the field-free situation occurs at the very beginning of the control process. Such observations on the first moment X (1) can be considered as additional insight for a comprehensive rationalization of control strategies as they evolve in time. Actually, it turns out that control fields take advantage from two simultaneous mechanisms: (i) population transfer improved by modifying the Rabi frequency, through the Stark shift directly affecting the central system; (ii) dynamical decoupling effects, through indirect process in the bath, preventing overall decoherence. It is worthwhile noting that, we have previously reported similar mechanisms with single cycle or dc fields [30]. As a final remark, these mechanisms being dynamically mixed, a non-Markovian diagnostic cannot merely be inferred from their analysis. This motivates the need to resort to other non-Markovian witnesses as is done hereafter. During the field-free evolution, the volume of accessible states illustrated in Fig.4 decreases very fast, in about 30 fs with a smooth monotonous decreasing profile. Nevertheless the decay is not exponential as it should be in a Markovian process. The duration of the control is fixed to 20 fs, i.e. less than the time for a complete decay of the volume. The resulting behaviors are displayed with the three control strategies C1-1, C2-2 and C1-2. Basically, after 5 fs, the decay is slightly faster than the field-free case and, more importantly, one observes some bumps, considered as clear signatures of non-Markovianity. Actually, the bumps arise at times close to 12 fs (for C1-1) or 17 fs (for C2-2) which could be associated with the maxima of the Stark shifts affecting the system energy gaps as displayed in Fig.2.  These analyses conclude that the field-dressed dynamics during the optimal control is more non-Markovian than the field free evolution. Moreover, one may question about the particular role of the quantum channel with the negative rate that should correspond to some backward flow. In order to observe the role of the different decoherence channels (Eq. (16)) during the evolution of a given initial state, we compute the weight of the three quantum channels as: Note that the operator G 0 (corresponding to the unity matrix) is not involved in the computation of coherence matrices so that the initial sum of |c k (t = 0)| 2 is equal to 0.5 and this sum is not conserved during the evolution since the decoherence matrix only describes the volume decrease and not its translation in the Bloch sphere. The upper panels of Fig.6 show the three canonical rates during the field-free and field-controlled evolution. The rates are given in increasing order so that channel k = 1 corresponds to the negative rate, which may to the field free case at the end of the C2-2 strategy; (iii) The highest positive rates are increased by the control fields, but more importantly their weights may decrease, for instance in the range 10-15 fs during control C2-2. As a consequence, the effective decay rate rate is basically affected by the combination of these effects. The increase of non-Markovianity during control does not necessarily imply that the channel with the negative rate plays the most significant role. In other words, an efficient control strategy for enforcing the bumps in the volume evolution, cannot merely be the tracking at each time of the channel with the negative partial rate, as it could be expected.
The volume of reachable states is a global property of the system. It is built from the dynamical map so that when it exhibits non Markovian witness, it is expected that similar signatures could be found in properties related to the evolution of a particular initial state.
As already discussed [14], non-Markovian witness can be seen in the system entropy (Eq. (21)) shown in Fig.7. The Markovian evolution of the entropy when the initial state is a pure state should be a monotonous evolution towards the value associated to the final Boltzmann mixture. In the present case, due to the energy gap, the final state is nearly the ground eigenstate so that the entropy profile should be a monotonous bell shape function. The non-Markovian signature is linked to any local decrease in the entropy which corresponds to a similar local bump in the purity T r(ρ 2 (t)) and therefore to an enhancement of the coherence. For instance, such a non-Markovian information back flow occurs between 11 and 18 fs in the field free evolution of state 1 and between 6 and 12 fs for state 2 (black curves). One observes that the dressed dynamics enhances this effect and more interestingly reduces the maximum entropy in a given time interval as during the control 1 to 2 (red dots in Fig.7).

V. CONCLUSION
This work is devoted to a detailed analysis of external field control versus dissipation in non-Markovian strongly coupled open quantum systems. A heterojunction is taken as an illustrative example with its specific parameters and spectral density, building up a spin-boson type Hamiltonian. With respect to methodology, the originality relies on a complete implementation of an optimal control scheme, together with a fully converged HEOM treatment of the master equation describing the time evolution of the two-level sub-system density matrix beyond a perturbative regime.
We put the emphasis on control scenarios aiming at producing physically relevant processes within the two-level sub-system interacting with its environmental bath. The ultimate goal is to protect against decoherence, the sub-system (such as a qubit), the control taking advantage from memory effects to draw back some information content from the bath to the sub-system. As a first attempt, we consider two targets, namely, the revival of an initial state |i (i = 1, 2) or a transition between the two states of the sub-system. The optimal control is precisely concerned with these goals through the populations of these states given in terms of the diagonal elements ρ 11 (t) and ρ 22 (t) of the sub-system density matrix. Once such control fields have been found, we address the consequences on the bath memory responses. Basically, we observe that non-Markovianity is increased during the optimally driven process. This is actually quantified through some typical signatures: time-dependent behavior of the volume of accessible states displaying bumps during its monotonic decay or the time-dependent behavior of the entropy exhibiting transitory decreases. At that point we have shown that a control aiming at the protection against decay of the subsystem characteristics provides, as a consequence, higher non-Markovian response of the bath. However, one of the main conclusions is that the mechanism does not necessarily increase the component on the quantum decay channel with the negative rate. We observe in most of the cases a decrease of the weight of the channel with the largest decay rate.
Similar behaviors have been obtained for other targets such as the one inducing relaxation towards the ground system eigenstate. The control performances remain however rather modest. The main reasons are the strong system-bath coupling and the limited range of our flash field amplitudes, in relation with their experimental feasibility. To go beyond such limitations, we have to refer to ultra short and intense laser pulses. This requires the introduction of an additional constraint in the optimal control scheme to correct the time integrated pulse area that, following Maxwell equations should be zero [56][57][58][59]. Finally, even more realistic calculations should be conducted with ab-initio transition dipoles, resulting from quantum chemistry codes.
As mid-term perspectives, future works should deal with exerting control directly on bath dynamics, in such a way to decrease decoherence of the sub-system, or in other words, achieve appropriate control of non-Markovianity to better protect sub-system characteristics. To that end, different strategies can be proposed: (i) Additional control of the environment through the introduction of a transition dipole among bath normal modes; (ii) Extraction of a collective mode from the bath so as to deal with a control involving an augmented active system, as has already been done in field-free heterojunction [38] or in a SQUID model [24].
We are actively pursuing research of these topics.