domain: phys chim

text:
Amine Jaouadi

Etienne Mangaud

Michèle Desouter-Lecomte

Re-exploring Control Strategies in a Non-Markovian Open Quantum System by Reinforcement Learning

In this study, we reexamine a recent optimal control simulation targeting the preparation of a superposition of two excited electronic states in the UV range in a complex molecular system. We revisit this control from the perspective of reinforcement learning, offering an efficient alternative to conventional quantum control methods. The two excited states are addressable by orthogonal polarizations and their superposition corresponds to a right or left localization of the electronic density. The pulse duration spans tens of femtoseconds to prevent excitation of higher excited bright states what leads to a strong perturbation by the nuclear motions. We modify an open source software by L. Giannelli et al., Phys. Lett. A, 434, 128054 (2022) that implements reinforcement learning with Lindblad dynamics, to introduce non-Markovianity of the surrounding either by timedependent rates or more exactly by using the hierarchical equations of motion with the QuTiP-BoFiN package. This extension opens the way to wider applications for non-Markovian environments, in particular when the active system interacts with a highly structured noise.

I. INTRODUCTION

Since several decades, quantum-state manipulation with electromagnetic fields is a central problem in many areas of physics and chemistry to prepare particular initial states or realizing unitary transformations (quantum gates). The hardware systems and the spectral range are very different, from spins in nuclear magnetic resonance (NMR) [START_REF] Vandersypen | Nmr techniques for quantum control and computation[END_REF] to molecular systems [START_REF] Keefer | Pathways to new applications for quantum control[END_REF] or complex photosynthetic systems [START_REF] Caruso | Coherent optimal control of photosynthetic molecules[END_REF][START_REF] Hoyer | Realistic and verifiable coherent control of excitonic states in a lightharvesting complex[END_REF] and systems involved in emerging quantum technology [START_REF] Glaser | Training Schrödinger's cat: quantum optimal control[END_REF][START_REF] Koch | Quantum optimal control in quantum technologies. strategic report on current status, visions and goals for research in europe[END_REF] for instance, superconducting quantum interference device (SQUID) [START_REF] Yang | n-qubit-controlled phase gate with superconducting quantuminterference devices coupled to a resonator[END_REF], trapped ions [START_REF] Rangan | Control of trapped-ion quantum states with optical pulses[END_REF] or atoms [START_REF] Jones | Fast quantum state control of a single trapped neutral atom[END_REF], nitrogenvacancy diamond centers [START_REF] Avinadav | Time-optimal universal control of two-level systems under strong driving[END_REF] or quantum dots [START_REF] Räsänen | Optimal laser control of double quantum dots[END_REF]. Quantum control has developed through various theoretical strategies as pump-dump schemes [START_REF] Tannor | Control of selectivity of chemical reaction via control of wave packet evolution[END_REF][START_REF] Kosloff | Wavepacket dancing: Achieving chemical selectivity by shaping light pulses[END_REF], coherent control [START_REF] Brumer | Laser control of molecular processes[END_REF], adiabatic methods [START_REF] Vitanov | Stimulated raman adiabatic passage in physics, chemistry, and beyond[END_REF][START_REF] Bergmann | Roadmap on stirap applications[END_REF], local [START_REF] Engel | Local control theory: Recent applications to energy and particle transfer processes in molecules[END_REF][START_REF] Mališ | Local control theory for superconducting qubits[END_REF] or Lyapunov control [START_REF] Hou | Optimal lyapunov-based quantum control for quantum systems[END_REF], Pontryagin optimal control [START_REF] Boscain | Introduction to the pontryagin maximum principle for quantum optimal control[END_REF], adaptative tracking [START_REF] Zhu | Quantum control design via adaptive tracking[END_REF] and optimal control theory (OCT) [START_REF] Peirce | Optimal control of quantum-mechanical systems: Existence, numerical approximation, and ap-plications[END_REF] that involves a rich variety of optimization algorithms, for instance Ra- * michele.desouter-lecomte@universite-paris-saclay.fr bitz monotonous convergent algorithm [START_REF] Zhu | Rapidly convergent iteration methods for quantum optimal control of population[END_REF][START_REF] Maday | New formulations of monotonically convergent quantum control algorithms[END_REF] or Krotov method [START_REF] Krotov | Global methods in optimal control theory[END_REF][START_REF] Palao | Optimal control theory for unitary transformations[END_REF], gradient ascent pulse engineering GRAPE [START_REF] Khaneja | Optimal control of coupled spin dynamics: design of nmr pulse sequences by gradient ascent algorithms[END_REF][START_REF] De Fouquieres | Second order gradient ascent pulse engineering[END_REF], chopped random basis optimization CRAB [START_REF] Doria | Optimal control technique for many-body quantum dynamics[END_REF][START_REF] Riaz | Optimal control methods for quantum gate preparation: a comparative study[END_REF].

Recently, there has been significant interest in applying reinforcement learning (RL), a distinctive machine learning technology, to quantum control and this begs the questions : Does RL suggest new control strategies ? How does it compare with standard algorithms [START_REF] Zhang | When does reinforcement learning stand out in quantum control? a comparative study on state preparation[END_REF]? What is the efficiency for control in dissipative environment [START_REF] An | Quantum optimal control of multilevel dissipative quantum systems with reinforcement learning[END_REF]? RL has already been applied in control for state preparation [START_REF] Bukov | Reinforcement learning in different phases of quantum control[END_REF] and gate realization [START_REF] An | Deep reinforcement learning for quantum gate control[END_REF][START_REF] Niu | Universal quantum control through deep reinforcement learning[END_REF] and quantum compiling [START_REF] Moro | Quantum compiling by deep reinforcement learning[END_REF] in quantum technology. Recently RL has recovered the well-known counter-intuitive STIRAP (Stimulated Raman Adiabatic Passage) pulse sequence in a three-state system [START_REF] Vitanov | Stimulated raman adiabatic passage in physics, chemistry, and beyond[END_REF][START_REF] Bergmann | Roadmap on stirap applications[END_REF]. In reference [START_REF] Porotti | Coherent transport of quantum states by deep reinforcement learning[END_REF] the laser may be on or off leading to a so-called digital-STIRAP ensuring an efficient transfer without the constraint of adiabaticity conditions [START_REF] Paparelle | Digitally stimulated raman passage by deep reinforcement learning[END_REF]. Conversely, in references [START_REF] Brown | Reinforcement learningenhanced protocols for coherent populationtransfer in three-level quantum systems[END_REF][START_REF] Giannelli | A tutorial on optimal control and reinforcement learning methods for quantum technologies[END_REF], the pulses are continuous, while the RL algorithm optimizes either the laser detuning or the Rabi frequencies.

In this work, our objective is to revisit with RL an optimal control simulation recently performed in a molecular system (phenylene ethynylene dimer) with C 2v symmetry [START_REF] Jaouadi | Laser-controlled electronic symmetry breaking in a phenylene ethynylene dimer: Simulation by the hierarchical equations of motion and optimal control[END_REF]. This benchmark case is interesting for different reasons. The system consists of two quasidegenerate excited electronic states of different symmetries. They are addressable by orthogonal polarizations. The preparation of a superposition of the two states corresponds to a right or left localization of the electronic density in a way similar to the localization in a doublewell by superposing the two lowest eigenstates. Creating such a coherence in complex systems with orthogonal polarizations has been recently discussed as a prospective avenue for achieving coherent control over excitonic energy transfer [START_REF] Tomasi | Coherent and controllable enhancement of light-harvesting efficiency[END_REF]. In absence of dissipation, this is a V-type three-state system where two excited states are coupled to the ground state by two orthogonal transition dipoles. An analytical solution may be derived to prepare the superposition with equal weights [START_REF] Breuil | Funneling dynamics in a phenylacetylene trimer: Coherent excitation of donor excitonic states and their superposition[END_REF]. This is an important landmark to test the control. On the other hand, the ideal electronic V-type system strongly interacts with the surrounding leading to a non-Markovian non-perturbative open system. The nuclear vibrations form two baths called the tuning and the coupling baths making fluctuate the energy gaps (also called longitudinal noise) and the electronic coupling (transversal noise) respectively.

In our previous OCT simulation, we employed hierarchical equations of motion (HEOM), which represent an exact method for addressing non-perturbative and non-Markovian open systems with Gaussian statistics [START_REF] Tanimura | Numerically "exact" approach to open quantum dynamics: The hierarchical equations of motion (HEOM)[END_REF][START_REF] Shi | Efficient propagation of the hierarchical equations of motion using the matrix product state method[END_REF][START_REF] Yan | Efficient propagation of the hierarchical equations of motion using the tucker and hierarchical tucker tensors[END_REF][START_REF] Dan | Theoretical study of nonadiabatic hydrogen atom scattering dynamics on metal surfaces using the hierarchical equations of motion method[END_REF][START_REF] Borrelli | Expanding the range of hierarchical equations of motion by tensortrain implementation[END_REF][START_REF] Mangaud | Survey of the hierarchical equations of motion in tensor-train format for non-markovian quantum dynamics[END_REF]. In order to investigate the RL control, we start from the open source software [50] presented in reference [START_REF] Giannelli | A tutorial on optimal control and reinforcement learning methods for quantum technologies[END_REF]. This software uses the libraries QuTip [START_REF] Johansson | Qutip: An open-source python framework for the dynamics of open quantum systems[END_REF] and TensorFlow [START_REF] Abadi | TensorFlow: Large-scale machine learning on heterogeneous systems[END_REF]. The software already implements Lindblad dynamics with the QuTip collapse operators. In this work, we introduce non-Markovianity in different ways. In a simplified strategy, we first consider time-dependent rates [START_REF] Piilo | Non-markovian quantum jumps[END_REF][START_REF] Hall | Canonical form of master equations and characterization of non-markovianity[END_REF] by using time-dependent QuTip collapse operators. The rates are calibrated from the decoherence matrix [START_REF] Hall | Canonical form of master equations and characterization of non-markovianity[END_REF] extracted from the exact HEOM dynamics [START_REF] Mangaud | Analysis of the non-markovianity for electron transfer reactions in an oligothiophene-fullerene heterojunction[END_REF][START_REF] Mangaud | Non-markovianity in the optimal control of an open quantum system described by hierarchical equations of motion[END_REF]. We then address the exact non-Markovian dynamics with the QuTip HEOM-solver [START_REF] Lambert | Qutip-bofin: A bosonic and fermionic numerical hierarchical-equations-of-motion library with applications in light-harvesting, quantum control, and single-molecule electronics[END_REF].

The paper is organized as follows. In section II we describe the model treated as an isolated or an open system interacting with two Bosonic baths. We summarize the Lindblad and HEOM operational equations in section III. The control by RL or OCT is presented in section IV. The RL results are given in section V and a comparison with OCT is made in section VI before concluding in section VII.

II. MODEL

The V-type three-state model of the dimer (1,3-bis(phenylethynyl)benzene) is schematized in Fig. 1. It is calibrated from ab initio data computed by B. Lasorne et al. [START_REF] Jaouadi | Laser-controlled electronic symmetry breaking in a phenylene ethynylene dimer: Simulation by the hierarchical equations of motion and optimal control[END_REF][START_REF] Ho | Vibronic properties of para-polyphenylene ethynylenes: TD-DFT insights[END_REF][START_REF] Ho | Diabatic pseudofragmentation and nonadiabatic excitationenergy transfer in meta-substituted dendrimer building blocks[END_REF][START_REF] Galiana | On the unusual Stokes shift in the smallest PPE dendrimer building block: Role of the vibronic symmetry on the band origin?[END_REF] with the density-functional theory (DFT) for the ground state S 0 and the time-dependent density-functional theory (TDFT) for the two excited states (S 1 (B 2 ) and S 2 (A 1 )). The energies of the two excited states (E S1 = 4.43 eV and E S2 = 4.47 eV) are taken at the equilibrium geometry of the ground state (planar with C 2v symmetry). The two states are bright and coupled to the ground state by orthogonal transition dipoles. The axes are chosen so that z(A 1 ) is the C 2 rotation axis, y(B 2 ) lies within the molecular plane, and x(B 1 ) is orthogonal to it. The respective transition dipoles are ⃗ µ 01 = (0, µ y , 0) with µ y = 3.96ea 0 and ⃗ µ 02 = (0, 0, µ z ) with µ z = -1.83ea 0 .

A. Isolated system

The system is defined by selecting the ground and the first two excited electronic states at the equilibrium geometry of the ground state (vertical Franck-Condon transition). It is an ideal system assumed to be frozen at this geometry. The two excited states are coupled by non-adiabatic interactions via a conical intersection [START_REF] Ho | Diabatic pseudofragmentation and nonadiabatic excitationenergy transfer in meta-substituted dendrimer building blocks[END_REF][START_REF] Galiana | On the unusual Stokes shift in the smallest PPE dendrimer building block: Role of the vibronic symmetry on the band origin?[END_REF]. We choose a diabatic representation with states of A 1 or B 2 symmetry so that the electronic coupling vanishes at that reference position and the system Hamiltonian is simply:

H 0 S = 2 j=0 |j⟩ ⟨j| (1)
The electronic coupling (between the two diabatic excited states) becomes different from zero when any vibration of B 2 symmetry is active. The two bright states are coupled to the ground state only radiatively and the time-dependent Hamiltonian at the dipolar approximation is :

H S (t) = H 0 S - 2 j=1 ⃗ µ 0j ⃗ E j (t) |0⟩ ⟨j| + hc (2)
where hc designates the hermitian conjugate.

The two fields are linearly polarized with ⃗ E 1 = (0, E y , 0) and ⃗ E 2 = (0, 0, E z ). In interaction representation (I) and with the rotating wave approximation (RWA) [START_REF] Fujii | Introduction to the rotating wave approximation (rwa): Two coherent oscillations[END_REF] the Hamiltonian becomes :

H RW A S,I (t) = - ℏ 2   0 Ω y (t) Ω z (t) Ω y (t) -2∆ y 0 Ω z (t) 0 -2∆ z   (3)
where the Rabi frequencies are Ω y (t) = µ y E y (t)/ℏ, Ω z (t) = µ z E z (t)/ℏ and E j (t) (j = y, z) are the pulse envelopes. ∆ y and ∆ z are the field detunings. The two Rabi frequencies or their variations with respect to a guess field are the actions that will be optimized by the RL algorithm. They are represented in Fig. 2. In our application, the target is the superposition of the two excited states with equal weights

|0⟩ → 1 √ 2 (|1⟩ + |2⟩) . ( 4
)
This target is different from the superposition of the initial state and one excited state that may be prepared by fractional-STIRAP [START_REF] Vitanov | Stimulated raman adiabatic passage in physics, chemistry, and beyond[END_REF] or by a π/2 pulse. By imposing equal Rabi frequencies at all times, i.e. pulses of the same duration T with amplitudes in the inverse ratio of the dipole moments, the target transition is realized if each area is equal to π/ √ 2 [START_REF] Breuil | Funneling dynamics in a phenylacetylene trimer: Coherent excitation of donor excitonic states and their superposition[END_REF] T 0 with j = y, z. It is a generalization of the wellknown π rule for complete population transfer [START_REF] Thomas | Validity of the rosen-zener conjecture for gaussian-modulated pulses[END_REF] or π/2 for creating a superposition involving the initial state.

Ω j (t)dt = π/ √ 2 (5)

B. Open quantum system

According to the chosen partition, all the nuclear vibrations belong to the baths. The A 1 modes make fluctuate the energies and the B 2 vibrations modify the electronic coupling that becomes different from zero when the C 2v symmetry is broken. The generic Hamiltonian of the system-bath partition is then written where H B is the ensemble of the vibrational modes assumed to be harmonic. The normal modes are assumed to be the same in each electronic states but their equilibrium positions differ. H SB is the linear system-bath coupling. The two groups of A 1 or B 2 modes then constitute different baths that may be called the tuning baths coupled to the diagonal elements |1⟩ ⟨1| and |2⟩ ⟨2| of the system Hamiltonian and the coupling bath coupled to the off-diagonal elements |1⟩ ⟨2| and |2⟩ ⟨1| between the two excited states. This kind of partition in the case of a conical intersection has been discussed in different works applying HEOM dynamics [START_REF] Jaouadi | Laser-controlled electronic symmetry breaking in a phenylene ethynylene dimer: Simulation by the hierarchical equations of motion and optimal control[END_REF][START_REF] Breuil | Funneling dynamics in a phenylacetylene trimer: Coherent excitation of donor excitonic states and their superposition[END_REF][START_REF] Duan | Quantum mechanical wave packet dynamics at a conical intersection with strong vibrational dissipation[END_REF][START_REF] Mangaud | Statistical distributions of the tuning and coupling collective modes at a conical intersection using the hierarchical equations of motion[END_REF]. The partition leads to a strong system-bath coupling and a non-Markovian master equation. The analysis of the dimer model is given in our previous work [START_REF] Jaouadi | Laser-controlled electronic symmetry breaking in a phenylene ethynylene dimer: Simulation by the hierarchical equations of motion and optimal control[END_REF] where is explained how are obtained the continuous spectral densities J(ω) for the tuning and coupling baths from ab initio data, i.e. from the energy gradients and gradient of the electronic coupling at the reference position. The spectral densities give the strength of the coupling to the system for each energy range of the baths. We select the main part of the spectral densities consisting in very sharp peaks around 1700 and 2300 cm -1 . The two spectral densi-ties J tuning (ω) and J coupling (ω) are presented in Fig. 3(a). Figure 3(b) gives the corresponding correlation functions of the collective mode of each bath

H(t) = H S (t) + H SB + H B ( 6
C (t) = 1 π +∞ -∞ dω J (ω) e iω(t) e βω -1 ( 7
)
where β is the Boltzmann constant. Due to the peaked shape of the spectral densities, the effective collective modes are under-damped and decay in about 200 fs. We zoom on the early time range of 20 fs, which is the pulse duration chosen in our simulations. Indeed, the pulse duration must be longer than about 10 fs to correspond to a sufficiently narrow spectral range to avoid the excitation of higher bright states. It is noticeable that the collective bath modes undergo a full oscillation within the 20 fs timescale. This serves as an indicator of non-Markovian behavior, as we will delve into further in the subsequent discussion.

III. DYNAMICAL METHODS

The dynamics of open quantum systems [START_REF] Breuer | The Theory of Open Quantum Systems[END_REF] has been reviewed, as seen in references [START_REF] Vega | Dynamics of non-Markovian open quantum systems[END_REF], and [START_REF] Koch | Controlling open quantum systems: tools, achievements, and limitations[END_REF] specifically focusing on control aspects. The applications of RL control [START_REF] An | Quantum optimal control of multilevel dissipative quantum systems with reinforcement learning[END_REF][START_REF] Giannelli | A tutorial on optimal control and reinforcement learning methods for quantum technologies[END_REF] usually assume a weak coupling and a Markovian bath treated by a Lindblad master equation [START_REF] Lindblad | On the generators of quantum dynamical semigroups[END_REF][START_REF] Manzano | A short introduction to the Lindblad master equation[END_REF]. We first summarize the main relations to introduce time-dependent rates and then we recall the operational equations for HEOM.

A. Lindblad master equation

For a N -dimensional system coupled to M dissipative channels, the generic Lindblad operator reads:

D(t) = M k,l Γ kl L k ρ t L † l - 1 2 L † l L k , ρ t (8)
where {A, B} = AB + BA denotes the anticommutator and the rates Γ kl are constant. Of-FIG. 3: (a) Spectral densities giving the strength of the system-bath coupling to the diagonal (J tuning ) or off-diagonal (J coupling = J S1S2 ) element of the system Hamiltonian block related to the two excited states. J tuning is J S1 and J S2 is assumed to be similar with J S2 = 0.797J S1 [START_REF] Jaouadi | Laser-controlled electronic symmetry breaking in a phenylene ethynylene dimer: Simulation by the hierarchical equations of motion and optimal control[END_REF]. The tuning modes are of A 1 symmetry and the coupling ones of B 2 symmetry. (b) Real and imaginary parts of the normalized bath correlation function C(t) (Eq.( 7)) at room temperature (zoom on the early 20 fs timescale). C(t) decays in 200 fs.

ten, only the main dissipative processes are retained, as the radiative decay towards a sink [START_REF] Giannelli | A tutorial on optimal control and reinforcement learning methods for quantum technologies[END_REF], or dephasing processes affecting the diagonal elements of the density matrix or relaxation inducing population transfer due to inter-state coupling [START_REF] An | Quantum optimal control of multilevel dissipative quantum systems with reinforcement learning[END_REF]. In this work, we will consider these two processes induced by the tuning and coupling baths respectively leading to four operators

L 1 = |1⟩ ⟨1|, L 2 = |2⟩ ⟨2|, L 3 = |1⟩ ⟨2| and L 4 = |2⟩ ⟨1|
(they are the collapse operators in QuTip [START_REF] Johansson | Qutip: An open-source python framework for the dynamics of open quantum systems[END_REF]). Moreover, we want to account for non-Markovian baths, at least on an approximate way by introducing time-dependent rates [START_REF] Piilo | Non-markovian quantum jumps[END_REF][START_REF] Hall | Canonical form of master equations and characterization of non-markovianity[END_REF]. This has given rise to many fundamental analysis [START_REF] Hall | Canonical form of master equations and characterization of non-markovianity[END_REF][START_REF] Mangaud | Non-markovianity in the optimal control of an open quantum system described by hierarchical equations of motion[END_REF][START_REF] Ángel Rivas | Quantum non-markovianity: characterization, quantification and detection[END_REF][START_REF] Witt | Exploring complete positivity in hierarchy equations of motion[END_REF] concerning the nonmarkovianity signature or the positivity of the dynamical map [START_REF] Théret | Complete positivity, positivity, and long-time asymptotic behavior in a two-level open quantum system[END_REF]. Indeed, the master equation of non-Markovian dynamics can be recast in a Lindblad form with time-dependent rates that may be transitory negative. If the Lindblad dissipator is expressed with the orthogonal basis set of N 2 operators formed by the normalized identity [START_REF] Kimura | The Bloch vector for N-level systems[END_REF][START_REF] Aerts | The extended Bloch representation of quantum mechanics and the hidden-measurement solution to the measurement problem[END_REF], which are the Pauli matrices for N = 2 and the Gell-Mann matrices for N = 3, the corresponding rate matrix is also called the time dependent decoherence matrix D jk [START_REF] Hall | Canonical form of master equations and characterization of non-markovianity[END_REF] [START_REF] Hall | Canonical form of master equations and characterization of non-markovianity[END_REF][START_REF] Andersson | Finding the kraus decomposition from a master equation and vice versa[END_REF]. This requires (N 2 -1) propagations of the basis operators performed here with HEOM as shown in [START_REF] Mangaud | Analysis of the non-markovianity for electron transfer reactions in an oligothiophene-fullerene heterojunction[END_REF][START_REF] Mangaud | Non-markovianity in the optimal control of an open quantum system described by hierarchical equations of motion[END_REF]. Some elements of the decoherence matrix will be used to calibrate the time-dependent rates of the selected collapse operators. This is a low-cost way to introduce easily some non-Markovianity with time-dependent collapse operators.

G 0 = I/ √ N and the N 2 -1 generators of SU (N ), G i (i = 1, ..., N 2 -1)
D(t) = N 2 -1 j,k=1 D jk (t) G j ρ t G k - 1 2 {G k G j , ρ t } . ( 9
D ij (t) = N 2 -1 m=1 T r [G m G i Λ t [G m (t)] G j ] (10) where Λ t [.] denotes the map of the time lo- cal non-Markovian master equation Ġm (t) = Λ t [G m (t)]
However, its efficiency might be somewhat limited, as the rates are calibrated based on field-free dynamics and are subject to potential modification by the applied fields [START_REF] Mangaud | Non-markovianity in the optimal control of an open quantum system described by hierarchical equations of motion[END_REF].

B. HEOM

We now summarize the main operational equations of the HEOM method. The system density matrix is the partial trace of the full density matrix ρ tot (t) over the bath degrees of freedom ρ(t) = T r B [ρ tot (t)]. The initial condition is assumed to be factorized ρ tot (0) = ρ(0)ρ eq where ρ eq is the density matrix of the baths at Boltzmann equilibrium at a given temperature. The HEOM may be considered as a numerically exact method for non-Markovian dynamics with harmonic baths when convergence is achieved by a relevant truncation of the hierarchy. The method is abundantly described in the literature, see for instance references [START_REF] Tanimura | Numerically "exact" approach to open quantum dynamics: The hierarchical equations of motion (HEOM)[END_REF] for a recent review or [START_REF] Mangaud | Survey of the hierarchical equations of motion in tensor-train format for non-markovian quantum dynamics[END_REF] for a pedagogical survey, [START_REF] Shi | Efficient propagation of the hierarchical equations of motion using the matrix product state method[END_REF][START_REF] Yan | Efficient propagation of the hierarchical equations of motion using the tucker and hierarchical tucker tensors[END_REF][START_REF] Dan | Theoretical study of nonadiabatic hydrogen atom scattering dynamics on metal surfaces using the hierarchical equations of motion method[END_REF][START_REF] Borrelli | Expanding the range of hierarchical equations of motion by tensortrain implementation[END_REF] for applications with the tensor-train format and [START_REF] Lambert | Qutip-bofin: A bosonic and fermionic numerical hierarchical-equations-of-motion library with applications in light-harvesting, quantum control, and single-molecule electronics[END_REF] for a review of different softwares, in particular that implemented in QuTip. We briefly recall that the master equation is solved by a time local system of coupled equations among auxiliary density matrices or auxiliary density operators (ADOs) arranged in a hierarchical structure. The algorithm requires a particular fit of the correlation function C(t) as a sum of K damped oscillatory terms also called artificial decaying modes

C (t) = K k=1 α k e iγ k t (11)
and

C * (t) = K k=1
αk e iγ k t . Analytical expressions for the α k , αk and γ k parameters can be derived from Eq.( 7) [START_REF] Pomyalov | The importance of initial correlations in rate dynamics: A consistent non-Markovian master equation approach[END_REF] when the spectral density is fitted by a sum of two-poles Tannor-Meier Lorentzian functions [77]

J (ω) = n l l=1 p l ω (ω + Ω l ) 2 + Γ 2 l (ω -Ω l ) 2 + Γ 2 l . ( 12
)
The parameters fitting the spectral densities of Fig. 3 are given in reference [START_REF] Jaouadi | Laser-controlled electronic symmetry breaking in a phenylene ethynylene dimer: Simulation by the hierarchical equations of motion and optimal control[END_REF]. Each ADO is labelled by a collective index n = {n 1 , • • • , n K } specifying the occupation number of each artificial mode. The system density matrix has the index n = {0, • • • , 0}. The HEOM equations are:

ρn (t) = L S (t)ρ n (t) + i K k=1 n k γ k ρ n (t) -i S, K k=1 ρ n + k (t) -i K k=1 n k α k Sρ n - k (t) -αk ρ n - k (t)S ( 13
)
where L S (t) is the system Liouvillian and

n + k = {n 1 , • • • , n k + 1, . . . , n K }, and n - k = {n 1 , • • • , n k -1, . . . , n K }.
The OCT simulations make use of the HEOM (Eqs. ( 13)) via our in-house developed software [START_REF] Jaouadi | Laser-controlled electronic symmetry breaking in a phenylene ethynylene dimer: Simulation by the hierarchical equations of motion and optimal control[END_REF]. In RL simulations [START_REF] Giannelli | A tutorial on optimal control and reinforcement learning methods for quantum technologies[END_REF] based on QuTip software, we use the QuTip-BOFiN HEOMsolver [START_REF] Lambert | Qutip-bofin: A bosonic and fermionic numerical hierarchical-equations-of-motion library with applications in light-harvesting, quantum control, and single-molecule electronics[END_REF] that allows the description of the Bosonic baths by giving the real and imaginary parts of the correlation function

C(t) = C R (t)+ iC I (t). For each bath, they are parametrized by a combination of decaying terms C R (t) = N R k=1 c R k e -γ R k t and C I (t) = N I k=1 c I k e -γ I k t
where the c k and γ R,I k are complex. This is an alternative to expansion of Eq.( 11) already adopted for the second order time non-local [START_REF] Meier | Non-Markovian evolution of the density operator in the presence of strong laser fields[END_REF] or time-local non-Markovian equations [START_REF] Kleinekathöfer | Non-Markovian theories based on a decomposition of the spectral density[END_REF]. The HEOM equations adapted to this partition of the correlation function in real and imaginary part are given in Eq.( 11) of reference [START_REF] Lambert | Qutip-bofin: A bosonic and fermionic numerical hierarchical-equations-of-motion library with applications in light-harvesting, quantum control, and single-molecule electronics[END_REF]. The analytical expressions of the c k and γ R,I k parameters when the spectral density is fitted by the twopole Lorentzian functions (Eq.( 12)) are given in references [START_REF] Meier | Non-Markovian evolution of the density operator in the presence of strong laser fields[END_REF][START_REF] Kleinekathöfer | Non-Markovian theories based on a decomposition of the spectral density[END_REF].

IV. CONTROL

A. Reinforcement learning

The RL algorithm is summarized in many references, for instance [START_REF] Sutton | Reinforcement Learning, Second Edition : An Introduction[END_REF][START_REF] Carleo | Machine learning and the physical sciences[END_REF]. By using the generic vocabulary, the principle is as follows. A target must be reached in an environment. At each time, an agent makes an observation and gets information about its state. The agent then chooses an action according to a policy to modify the state. The agent obtains a reward that estimates the progress towards the target. In our application, we have thus to define the environment, the agent, the action, and the reward. The four points are represented in Fig. 4. The RL environment is the active system and its surrounding, i.e., the V-three-level system coupled to both tuning and coupling baths. The observation is the state of the system de-scribed by the reduced density matrix solution of a master equation. The agent is an algorithm called REINFORCE [START_REF] Sutton | Reinforcement Learning, Second Edition : An Introduction[END_REF]. It uses a neural network with three hidden layers in our application. The input layer contains all the density matrix elements (nine in our case). The output layer may provide discrete values or a continuous distribution. The outputs are discrete when only some actions are available, for instance if the laser may be only on or off, giving four outcomes in the two-pulse case [START_REF] Paparelle | Digitally stimulated raman passage by deep reinforcement learning[END_REF]. When the distribution is continuous, the output layer gives parameters of this distribution, for instance the mean of Gaussian distributions for the actions a y,z that provide the Rabi frequencies Ω y and Ω z of the two pulses. The reward is the control fidelity r t = T r(ρ † target ρ(t)). The policy is the conditional probability π (a t |s t ) that the agent takes action a when the system state is s. The action at a time t only depends on the state at that time and the process is called a Markov decision chain. Note that this does not mean that the dynamics of the system must be Markovian. The Markov decision chain means that when two states s and s ′ observed by the agent are the same, the probability to choose a is the same regardless of the history to reach the state s.

We employ the policy gradient method [START_REF] Williams | Simple statistical gradientfollowing algorithms for connectionist rein-forcement learning[END_REF], which is a technique used in reinforcement learning. It involves adjusting the policy's parameters aiming to maximize the cumulative reward over time. All the parameters of the network are represented by the global index θ. They are initialized at random. One then generates a batch of M episodes with the current policy π (a t |s t ). An episode or trajectory τ is divided in N time steps and lasts T = N δt. For each episode, one collects the N data triplets (s t , a t , r t ) where t = iδt for the i th time step. The performance of the agent is estimated by the so-called return R that depends on the network parameters θ and is the main tool to optimize the policy. For each episode, the return R(τ ) may be defined in more or less sophisticated ways by the simple sum of all the rewards r t or a weighted sum of these with a discount rate [START_REF] An | Quantum optimal control of multilevel dissipative quantum systems with reinforcement learning[END_REF][START_REF] Bukov | Reinforcement learning in different phases of quantum control[END_REF][START_REF] Martín-Guerrero | Reinforcement learning and physics[END_REF]. Here, R(τ ) is the sum of the r t FIG. 4: Representaion of a cycle of the RL. At each time, the reduced density matrix of the three-state system coupled to its surrounding is the input of the neural network. The policy π (a t |s t ) is optimized and provides two actions that are the pulse Rabi frequencies. The reward is the fidelity to reach the target, which is here the superposed state of the two excited states. and all the r t = 0 if t < T so that R(τ ) = r N , i.e. it is given by the final control fidelity. For the bunch of M episodes driven by the same policy, the return is the expectation value

E [R] = M τ =1 p θ (τ )R(τ ) ( 14
)
where p θ (τ ) is the probability of driving trajectory τ . Following reference [START_REF] Marquardt | Machine learning and quantum devices[END_REF] we summarize the main points of the policy optimization. The probability of each trajectory is different since the actions are taken at random in the current policy. It is a product for each time step of the probability p(s t+1 |a t , s t ) to have a transition from state s t to state s t+1 induced by action a t times the probability for the agent to choose action a t for state s t :

p θ (τ ) = N t p(s t+1 |a t , s t )π θ (a t , s t ). ( 15
)
p(s t+1 |a t , s t ) does not depend on the parameters θ but only on the system dynamics. Therefore, the gradient of the average return involves only the gradient of the policy

∇ θ p θ (τ ) = N t=1 N t ′ =1 p( s t ′ +1 |a t ′ , s t ′ )π θ (a t ′ , s t ′ ) × ∇ θ ln π θ (a t , s t ) = p θ (τ ) N t=1 ∇ θ ln π θ (a t , s t ). ( 16
)
The network parameters are optimized so that their gradient ∇θ is parallel to the gradient of the average return with a factor η called the learning rate ∇θ = +η∇ θ E [R] and by Eqs. ( 14) and ( 16) one has

∇θ = +η M τ =1 p θ (τ )R(τ ) N t=1 ∇ θ ln π θ (a t , s t ) = +ηE R(τ ) N t=1 ∇ θ ln π θ (a t , s t ) . ( 17
)
The parameters are updated according to the logarithmic gradient of the policy times the return and the learning rate that must be chosen not too fast and not too slow. Since the gradient contains the return, all the actions become more likely, the more the return is larger. An important point is that the optimization algorithm of the network parameters is independent of the underlying dynamical model. This is a difference with the standard optimization in OCT where the gradient of the final fidelity involves the system Hamiltonian. RL operates beyond static databases; it collects data during training.

B. Optimal control

The optimal field is built by iterations to maximize the cost functional that is the fidelity F = T r ρ † target ρ(T ) at the final time T with constraints to restrain the field intensity and to fulfill the master equation at any time. The optimization is performed here by Rabitz' monotonously convergent algorithm [START_REF] Zhu | Rapidly convergent iteration methods for quantum optimal control of population[END_REF][START_REF] Ohtsuki | Monotonically convergent algorithm for quantum optimal control with dissipation[END_REF] that involves forward and backward propagation of the system density matrices with initial condition ρ(0) and of an auxiliary system density matrix χ(t) with final condition χ(T ) = ρ target .

It is worth noting that a two-point boundary-value quantum control paradigm (TBQCP) has been presented in the literature as an accelerated convergent algorithm [START_REF] Ho | Accelerated monotonic convergence of optimal control over quantum dynamics[END_REF]. However for the sake of simplicity this method is not used in our simulations. Dynamics is driven with HEOM. The iterations begin with a guess field that strongly influences the final field. The operational relations for the backward propagation are given in references [START_REF] Jaouadi | Laser-controlled electronic symmetry breaking in a phenylene ethynylene dimer: Simulation by the hierarchical equations of motion and optimal control[END_REF][START_REF] Mangaud | Non-markovianity in the optimal control of an open quantum system described by hierarchical equations of motion[END_REF]. The field at each iteration k is obtained by ε (k) = ε (k-1) + ∆ε (k) where ∆ε (k) is estimated by

∆ε(t) = 1 α ℑm T r χ(t) p µ p , ρ(t) ( 18
)
where α is the intensity penalty factor. Note that we do not use RWA in this approach.

V. CONTROL BY RL

We choose a pulse duration of T = 20 fs. This is relevant to avoid a too large spectral band that would imply higher bright excited states not included in the model system. For RL simulations, the laser detunings are assumed to vanish, so the only optimized parameters are the two Rabi frequencies Ω y and Ω z at all times. In all the RL examples, the results are given in reduced units for the time (t/T ) and the Rabi frequency (T Ω). Conventional units are used in some OCT examples.

The three hidden network layers contain 100, 50 and 30 neurons. The learning rate has its standard value η = 10 -3 . Our investigation has confirmed the critical significance of these meta-parameters. Reducing the number of neurons results in a decelerated convergence rate, while elevating it prolongs computational duration without commensurate convergence enhancement. Increasing the learning rate by a factor 10 is not efficient to reach the desired target. Data are collected during a bunch of M = 10 episodes that are divided in 50 time steps.

A. RL in the isolated system

We first consider the ideal case without dissipation. In RL, the process begins by two Rabi frequencies chosen at random in a given range. The analytical solution (Eq.( 5)) is a landmark. By assuming a constant pulse envelope, the integrated Rabi frequency is T Ω and the best value should be π/ √ 2 = 2.22. We begin the RL simulation with an initial interval with T Ω min = 0 and T Ω max = 3 that would give an area larger than the best analytical value of 2.22. Figure 5(a) displays the first iteration with initial random Rabi frequencies, here 1.54 for Ω y and 1.55 for Ω z . The integrated frequencies are too small and the ground state is not completely depopulated. Figure 5(b) gives the outcome after 100 episodes. It reproduces the analytical result and is obtained after about 60 episodes as illustrated in Fig. 6. RL provides the good integrated frequencies but with very simple pulses since the envelopes are quasi constant. Figure 6 displays the return achieved during five simulations of 100 episodes. The random initial conditions differ from T Ω = 1.5 by about 10 percents giving a return close to 0.75. Notably, the convergence rate exhibits variability and doesn't follow the typical monotonic pattern observed in conventional OCT algorithms. However, from 60 iterations, the rate achieved its highest value.

To test the algorithm, we start with a larger initial interval with T Ω max = 9. As the initial frequencies are random in this range, RL converges towards different possibilities but it is worth noting that it always finds a solution close to the analytical result. Figure 7(a) illustrates a case where convergence occurs towards the expected value π/ √ 2 with a very good final coherence. In Fig. 7(b) one sees that according to the random initial values, optimization provides a solution with higher intensity and an area close to 3π/ √ 2 leading to a supplementary complete Rabi oscillation before the final coher-FIG. 5: Optimization of the superposed target state (Eq.( 4)) by RL in the isolated V-threelevel system and T Ω max = 3. The upper panels give the Rabi frequencies in reduced units (T Ω), the middle ones show the populations in each state and the lower ones, the modulus of the coherence ρ 12 between the two excited states. ence creation close to 0.5.

B. RL with Lindblad dynamics

In this section, we will undertake a comparative analysis between scenarios involving constant Lindblad rates and those incorporating time-dependent variations. The constant rates associated to the four selected Lindblad operators (QuTip collapse operators) L k defined in section III A are (in reduced units T Γ): T Γ 11 = 1 (L 1 ), T Γ 22 = 0.8 (L 2 ), T Γ 12 = 0.36 (L 3 ) and T Γ 21 = 0.16 (L 4 ). The operators L 1 and L 2 couple to the tuning baths in excited states S 1 and S 2 respectively. The ratio of the rates is approximated by J S2 /J S1 as elaborated upon in reference [START_REF] Jaouadi | Laser-controlled electronic symmetry breaking in a phenylene ethynylene dimer: Simulation by the hierarchical equations of motion and optimal control[END_REF]. The operators L 3 and L 4 induce non-adiabatic transitions. The reduced FIG. 6: Return (Eq.( 14)) during the optimization in the isolated system presented in figure 5 for five different simulations displayed in different colors. Each simulation runs 100 episodes. The green dotted line serves as an indicator, highlighting the target return value of 1. The random initial conditions differ from T Ω = 1.5 by about 10 percents giving a return close to 0.75. rate T Γ 12 and T Γ 21 are different as expected from the detailed balance. They are calibrated to roughly approximate the exact HEOM fieldfree dynamics at least during the early dynamics. Dynamics is performed with the QuTip mesolve solver [START_REF] Johansson | Qutip: An open-source python framework for the dynamics of open quantum systems[END_REF].

Non-Markovianity may be taken into account in an approximated way by time-dependent rates. The transitory negativity of the sum of the canonical rates that are the eigenvalues of the decoherence matrix (Eq.( 9)) is one of the signature. This sum obtained for the field-free dynamics is given in Fig. 8(a) in reduced units (T = 20 fs). It is obvious that its damped oscillation follows that of the bath correlation functions (see Fig. 3(b)). This illustrates that for this type of partition, the non-Markovianity is closely linked to the damped vibrational motion of the collective modes. Indeed, if the collective effective mode is underdamped, the nuclei oscillate and transitory return to the initial reference position. This tends to restore the system FIG. 7: Optimization of the superposed target state (Eq.( 4)) by RL in the isolated V-threelevel system and T Ω max = 9. The upper, middle and lower panels are as in in its initial condition. The decay towards equilibrium is not monotonous.

Figure 8(b) presents the time-dependent rates associated to the four Lindblad operators describing the energy tuning (L 1 , L 2 ) and the interstate transition (L 3 , L 4 ) induced by the coupling bath. Their shape are approximated from those of some elements of the decoherence matrix (Eq.( 9)) by considering the basis operator G k corresponding to the 1 -2 transition (analog of σ x in the two-state case) and one operator corresponding to an energy gap. The amplitudes are calibrated as in the constant rate case from the field-free HEOM dynamics. The functions are fitted by polynomials or by the product of a sine function times a decreasing exponential. These functions are introduced in the time-dependent collapse operators of QuTip by using the mesolve solver [START_REF] Johansson | Qutip: An open-source python framework for the dynamics of open quantum systems[END_REF].

Figure 9(a) shows the dynamics after 100 episodes with constant decay rates. The re-FIG. 8: (a) Sum of the canonical rates T Γ c k (in reduced units) of the field-free dynamics. They are the eigenvalues of decoherence matrix (Eq.( 9)). (b) Time-dependent rates T Γ associated to the four Lindblad operators describing the energy tuning (L 1 , L 2 with rates Γ 11 and Γ 22 ) and the interstate transition (L 3 , L 4 with rates Γ 12 and Γ 21 ) induced by the coupling bath.

turn is only 0.8 and it saturates after 60 iterations. Figure 9(b) presents the control with the time-dependent rates. The return is slightly improved. However, this is due mainly to the better depletion of the ground state and not to a better superposition. The optimal envelopes remain very simple and quasi constant in both cases. It is a bit disconcerting that RL behaves on a very similar way with constant or timedependent rates. We will compare these results with the OCT optimization in section VI.

C. RL with HEOM dynamics

For a 20 fs simulation, truncating at level 6 of the hierarchy proves to be satisfactory. However, it's worth noting that for longer dynamics, such as achieving a field-free asymptotic state, a higher level 9 becomes necessary. The implementation of the RL algorithm with HEOM requires some comments. (i) During the Markov FIG. 9: Optimization of the target state (Eq.( 4)) by RL with Lindblad or HEOM dynamics. The upper, middle and lower panels are as in Fig. 5. The first episode with random initial conditions is given in dashed lines. The optimized results are given after 100 episodes. In each case, the return saturates after about 60 episodes. decision chain the solver is called repetitively for each time step of the chain . . . s t → a t → r t → s t+1 → a t+1 → r t+1 . . .. All the ADOs describing the state of the surrounding must be saved for the following decision step so that each bath retains its configuration and does not restart with the initial conditions of the baths with ADOs equal to zero. This is a difficulty that does not concern the local Lindblad dynamics. (ii) In our application, each spectral density (see Fig. 3) is fitted by two Tannor-Meier Lorentzian functions [START_REF] Meier | Non-Markovian evolution of the density operator in the presence of strong laser fields[END_REF] leading to four decay modes for each bath. As the spectral densities are centered at high frequencies, we do not include Matsubara terms at room temperature. We use the description of the baths by the expansion of the real and imaginary parts of the correlation functions using the BosonicBath ap-plication of the QuTip HEOMSolver [START_REF] Lambert | Qutip-bofin: A bosonic and fermionic numerical hierarchical-equations-of-motion library with applications in light-harvesting, quantum control, and single-molecule electronics[END_REF]. When working with reduced units the real or imaginary parts of the c R,I k coefficients must be scaled by T 2 and the rates γ R,I k by T as usually. (iii) In the HEOM solver, the system-bath coupling operators are not written in interaction representation. Therefore, we have re-transformed the Hamiltonian in Schrödinger representation without the RWA approximation. We use 100 time steps in each episode which is enough to satisfy the Nyquist-Shannon sampling rule [START_REF] Marks | Handbook of Fourier Analysis & Its Applications[END_REF] for the fields in the Schrödinger representation.

The RL optimization with HEOM in Schrödinger representation without RWA is given in Fig. 9(c). The actions are more erratic due to the oscillation of the field in this representation. After 100 episodes, the envelopes have a slightly higher amplitude than in the Lindblad simulations. The populations and coherence behave on a similar way in each simulation. The return is only 0.53 primarily due to the less-than-optimal depletion of the ground state and the difference of population in the two excited states. Other examples with guess fields are given in Fig. 10.

Finally, we explore another strategy. We impose a guess field for the RL control by choosing the actions to be the variation δΩ of the Rabi frequencies with respect to the guess (see Fig. 2(b)). These trial fields are a sine square envelope or a constant satisfying the π/ √ 2 rule. Simulations are carried out in Schrödinger representation without RWA at level 6 of the hierarchy. The actions T δΩ are taken in an interval [-2,2] in reduced units. The Rabi frequencies of the guess fields and of the RL optimization after 100 episodes are given in Fig. 10. The envelopes are only very slightly modified during the optimization. The first action is always the largest and shifts the guess envelopes by adding a constant value. The further fluctuations remain of weak amplitude. Increasing the initial interval only modifies the initial shift. Only the area increases, which generally enhances the depletion of the ground state but not the target coherence. The sine square envelope is the best guess giving a return of 80%. The constant envelope provides only 70%.

VI. COMPARISON RL-OCT

The envelopes generated by RL consistently exhibit a high degree of simplicity, characterized by their quasi-constant profile when no guess is imposed. Given that the fields generated by standard OCT typically exhibit a higher degree of structure [START_REF] Jaouadi | Laser-controlled electronic symmetry breaking in a phenylene ethynylene dimer: Simulation by the hierarchical equations of motion and optimal control[END_REF], we compare in Fig. 11 the fields optimized by OCT with the same initial guess fields drawn in solid lines in Fig. 10. Simulations are performed utilizing our HEOM code [START_REF] Jaouadi | Laser-controlled electronic symmetry breaking in a phenylene ethynylene dimer: Simulation by the hierarchical equations of motion and optimal control[END_REF] in Schrödinger representation without RWA at level 6 of the hierarchy. Standard OCT optimizes the field amplitudes and not only the envelopes on a time grid. This offers more flexibility and may slightly modify the carrier frequency. The α penalty factor (Eq.18)) is fixed to 2 × 10 -4 . It influences the optimization rate. The field amplitudes increase regularly at each iteration. We take two snapshots to remain in the same order of magnitude as in the RL simulation. The results after 15 iterations are shown in Fig. 11 (a) and(b). We draw the fields times the dipole moment T µ y E y (t) and T µ z E z (t) in reduced units so that the envelopes may be compared with the reduced Rabi frequencies T Ω used in the RL optimization. OCT reshapes the envelopes more strongly than RL. In particular, the two envelopes do not remain similar. However, when the maximum field amplitude are in the same range, the return is similar around 80% for the sine square case and reaches 80% versus 70% in RL for the constant guess.

VII. SUMMARY AND CONCLUSION

Examining the potential of RL in quantum control has primarily been explored within the domain of quantum information [START_REF] An | Quantum optimal control of multilevel dissipative quantum systems with reinforcement learning[END_REF][START_REF] An | Deep reinforcement learning for quantum gate control[END_REF][START_REF] Niu | Universal quantum control through deep reinforcement learning[END_REF][START_REF] Moro | Quantum compiling by deep reinforcement learning[END_REF]. This analysis is particularly interesting in the context of retrieving the STIRAP scheme using either digital pulses [START_REF] Paparelle | Digitally stimulated raman passage by deep reinforcement learning[END_REF] or continuous ones [START_REF] Brown | Reinforcement learningenhanced protocols for coherent populationtransfer in three-level quantum systems[END_REF][START_REF] Giannelli | A tutorial on optimal control and reinforcement learning methods for quantum technologies[END_REF]. The outstanding property is the ability to propose strategies without any prior knowledge of the system leading to the denomination as a "model-free algorithm" [START_REF] Marquardt | Machine learning and quantum devices[END_REF]. The main question is to see whether RL will find new strategies in particular in presence of an environment. Most of the previous works have studied examples with dissipation treated by Lindblad master equation, i.e. for a Markovian dynamics [START_REF] An | Quantum optimal control of multilevel dissipative quantum systems with reinforcement learning[END_REF][START_REF] Giannelli | A tutorial on optimal control and reinforcement learning methods for quantum technologies[END_REF]. However, even if RL is built on a succession of decision Markov processes, non-Markovian noise could influence the system dynamics [START_REF] Porotti | Coherent transport of quantum states by deep reinforcement learning[END_REF].

In this work, we have revisited a control in a system with a strong non-Markovian dynamics due to the coupling to baths with highly structured spectral densities leading to long bath correlation times. The model is calibrated from ab initio data [START_REF] Jaouadi | Laser-controlled electronic symmetry breaking in a phenylene ethynylene dimer: Simulation by the hierarchical equations of motion and optimal control[END_REF][START_REF] Ho | Vibronic properties of para-polyphenylene ethynylenes: TD-DFT insights[END_REF][START_REF] Galiana | On the unusual Stokes shift in the smallest PPE dendrimer building block: Role of the vibronic symmetry on the band origin?[END_REF]. We have used an open source software [START_REF] Giannelli | A tutorial on optimal control and reinforcement learning methods for quantum technologies[END_REF]50] based on the policy gradient method for the optimization and on QuTip mesolve solver of the Lindblad master equation [START_REF] Johansson | Qutip: An open-source python framework for the dynamics of open quantum systems[END_REF]. We have enhanced its functionality to address non-Markovian dynamics. In a first approximate go-between step we have incorporated time-dependent rate constants derived from the field-free HEOM dynamics. Then we have interfaced the RL algorithm with the HEOM solver of the QuTip BoFiN package [START_REF] Lambert | Qutip-bofin: A bosonic and fermionic numerical hierarchical-equations-of-motion library with applications in light-harvesting, quantum control, and single-molecule electronics[END_REF].

An analytical solution exists to create the target superposition of two excited states addressable with orthogonal dipoles in an isolated Vsystem. It is of significant interest to assess the proficiency of RL to recover the expected solution from random initial conditions in a given interval for the Rabi frequencies. RL finds a very simple but efficient solution of straightforward quasi-constant envelopes satisfying the integrated Rabi frequency rule. The amplitudes display minimal variations, typically within a few percentage points.

For each level of complexity of the dynamics, we observe that the return saturates after about 60 iterations even if the target is not perfectly reached. The Rabi frequencies are always very simple, nearly constant when no guess is imposed. The proposal is basic and robust. Even if it is not completely satisfactory, RL does not go further. When dynamics is driven in Schrödinger representation without RWA, the process of optimizing using reinforcement learning (RL) exhibits increased complexity. This complexity is reflected in the erratic behavior of the envelopes from one step to another. However, it's worth noting that despite these fluctuations, a smoother average trajectory is observed with only minor fluctuations.

Lindblad dynamics even with some timedependent rates cannot take into account a possible influence of the field on the baths. On the contrary, the memory kernel of HEOM contains the time dependent Hamiltonian [START_REF] Meier | Non-Markovian evolution of the density operator in the presence of strong laser fields[END_REF] and this could in principle induce an effect on the bath dynamics [START_REF] Chenel | Control in a dissipative environment: The example of a cope rearrangement[END_REF]. Indeed, the decoherence matrix and thus the rates are modified by the field [START_REF] Mangaud | Non-markovianity in the optimal control of an open quantum system described by hierarchical equations of motion[END_REF]. However, in our application the behavior is qualitatively the same for the approximate non-Markovian approach or for exact HEOM. RL successfully captures the crucial condition regarding the integrated Rabi frequency, yet it does not discover novel strategies to fight dissipation. It would be powerful to increase the number of available actions, enabling optimization of detuning parameters as well or directly the amplitude of the fields and not only the envelopes. Another possibility would be to let the algorithm choose a guess field. Furthermore, the exploration of more sophisticated RL algorithms holds promise for future investigations [START_REF] Reuer | Realizing a deep reinforcement learning agent discovering realtime feedback control strategies for a quantum system[END_REF][START_REF] Schulman | Trust region policy optimization[END_REF].

The simplicity of the RL envelopes suggests to confront the standard OCT and to see if it can yield superior results. OCT exploits the system dynamics and may seem more flexible since it optimizes the field amplitude and the carrier frequency and possibly finds some chirp effect. However, in our example OCT is not more efficient to reach the target with dissipation. By imposing the same guess fields, RL and OCT provide different optimal fields ensuring similar return. The reshaping is stronger in OCT that proposes different envelopes for the two polarizations what RL does not do. When the envelope amplitudes are maintained in the same range as in RL, OCT slighly improves the depletion of the ground state but not really the preparation of the superposition with equal weights. Both strategies, RL and OCT depend on the guess fields and the optimal fields are different. However, they ensure similar final dynamics and the perfect target is not achieved neither by RL nor by OCT control due to the strong dissipation.

Our example is a complex system strongly coupled to a structured environment with laser pulses in the femtosecond range. RL seems more adapted to treat quantum information in another spectral range operating with very simple square box envelopes and weakly coupled Markovian noises [START_REF] Baum | Experimental deep reinforcement learning for error-robust gateset design on a superconducting quantum computer[END_REF].

VIII. DATA AVAILABILITY

The data are available upon request to the authors. The modified ThreeLS.py file of Luigi Giannelli's open-source software [START_REF] Giannelli | A tutorial on optimal control and reinforcement learning methods for quantum technologies[END_REF]50] allowing dynamics with HEOM in Schrödinger representation without RWA by using the QuTip BoFiN package [START_REF] Lambert | Qutip-bofin: A bosonic and fermionic numerical hierarchical-equations-of-motion library with applications in light-harvesting, quantum control, and single-molecule electronics[END_REF] is given in supplementary material [91].

FIG. 1 :

 1 FIG. 1: Schematic representation of the V-type model in 1,3-bis(phenylethynyl)benzene. The two excited states of symmetry A 1 and B 2 at the geometry of the ground state are delocalized over the two sites. Their superposition with equal weights corresponds to a localization on the left or right sites. The excited states are addressable by orthogonal dipole moments. The electronic sub-system is coupled to the vibrational baths. The tuning bath (longitudinal noise) making fluctuate the energy gap gathers the symmetric A 1 modes and the coupling bath (transversal noise) varying the electronic coupling contains the B 2 modes.

)FIG. 2 :

 2 FIG. 2: Actions optimized by the RL algorithm to prepare the target superposed state in the V-type system. (a) The actions are the Rabi frequencies Ω y and Ω z ; (b) the actions are the variations of the Rabi frequencies with respect to guess fields. The laser detuning ∆ y and ∆ z are assumed to vanish in this application.

  )The eigenvalues are the canonical decay rate Γ c k associated to the time-dependent decay channels. The decoherence matrix is given by[START_REF] Hall | Canonical form of master equations and characterization of non-markovianity[END_REF]

  FIG.5: Optimization of the superposed target state (Eq.(4)) by RL in the isolated V-threelevel system and T Ω max = 3. The upper panels give the Rabi frequencies in reduced units (T Ω), the middle ones show the populations in each state and the lower ones, the modulus of the coherence ρ 12 between the two excited states. (a) First iteration with initial random Rabi frequencies. The integrated frequencies are 1.54 for Ω y and 1.55 for Ω z . (b) After 100 episodes. Both areas are 2.22 (π/ √ 2), the best expected result.

  FIG. 7: Optimization of the superposed target state (Eq.(4)) by RL in the isolated V-threelevel system and T Ω max = 9. The upper, middle and lower panels are as in Fig. 5. (a) Convergence towards the best expected result (π/ √ 2). The integrated frequencies are 2.286 for Ω y and 2.254 for Ω z . (b) Convergence towards (3 π/ √ 2) leading to a supplementary Rabi oscillation. The integrated frequencies are 6.673 for Ω y and 6.825 for Ω z .

  FIG. 9: Optimization of the target state (Eq.(4)) by RL with Lindblad or HEOM dynamics. The upper, middle and lower panels are as in Fig. 5. The first episode with random initial conditions is given in dashed lines. The optimized results are given after 100 episodes. In each case, the return saturates after about 60 episodes. (a) Dynamics with constant Lindblad rates. The area of the Rabi frequencies are 2.35 and 2.48 respectively. (b) Dynamics with timedependent rates. The area are 2.38 and 2.39. (c) HEOM dynamics at level 6 of the hierarchy in Schödinger representation without RWA.

FIG. 10 :

 10 FIG. 10: RL optimization by RL with different guess fields. The actions are the variation of the Rabi frequencies δΩ in a range of reduced units [-2, 2]. Dynamics is computed by HEOM in Schrödinger representation without RWA at level 6 of the hierarchy. The upper panels give the Rabi frequencies Ω y and Ω z in reduced units, the middle ones, the populations and the lower ones, the modulus of the coherence between the two excited states. (a) guess fields with sine square envelopes of integrated Rabi frequencies π/ √ 2 in solid line. The corresponding RL fields after 100 episodes in dashed line (the areas are 3.77 and 3.96). (b) guess fields with constant envelope in solid line and the corresponding RL fields after 100 episodes in dashed line (the areas are 3.87 and 3.95).

FIG. 11 :

 11 FIG. 11: Optimization by standard OCT after 15 iterations. (a) Sine square guess fields, (b) constant envelope guess fields. The upper panels give the fields times the corresponding transition dipole in reduced units T µ y E y (t) and T µ z E z (t). Dynamics is computed by HEOM in Schrödinger representation without RWA at level 6 of hierarchy. The lower panels display the populations and the modulus of the coherence between the two excited states.

IX. ACKNOWLEDGEMENTS

The authors thanks Dr O. Atabek for stimulating discussions and encouragements to explore this project. Dr B. Lasorne and PhD J. Galiana are warmly acknowledged for providing the ab initio data calibrating the model.