Efficient exploration of Hamiltonian parameter space for optimal control of non-Markovian open quantum systems
Gerald E. Fux, Eoin Butler, Paul R. Eastham, Brendon W. Lovett, Jonathan Keeling
EEfficient exploration of Hamiltonian parameter space for optimal control ofnon-Markovian open quantum systems
Gerald E. Fux, Eoin Butler, Paul R. Eastham, Brendon W. Lovett, and Jonathan Keeling SUPA, School of Physics and Astronomy, University of St Andrews, St Andrews, KY16 9SS, United Kingdom School of Physics, Trinity College Dublin, College Green, Dublin 2, Ireland (Dated: January 11, 2021)We present a general method to efficiently design quantum control procedures for non-Markovianproblems and illustrate it by optimizing the shape of a laser pulse to prepare a quantum dot ina specific state. The optimization of open quantum system dynamics with strong coupling tostructured environments—where time-local descriptions fail—is a computationally challenging task.We modify the numerically exact time evolving matrix product operator (TEMPO) method, suchthat it allows the repeated computation of the time evolution of the reduced system density matrixfor various sets of control parameters at very low computational cost. This method is potentiallyuseful for studying numerous quantum optimal control problems, in particular in solid state quantumdevices where the coupling to vibrational modes is typically strong.
One of the main challenges in current quantum tech-nology is to avoid or mitigate decoherence. Quantumoptimal control theory seeks to address this challenge bysearching for protocols that achieve the highest fidelityof a process for a given experimental setup [1–5]. Forthis to be successful, it is necessary to be able to accu-rately compute the dynamics of the system under theinfluence of the environment. The majority of researchon optimal control for open quantum systems has beencarried out in the Markovian limit, where one assumes aweak system-environment coupling and environment cor-relations that are short compared to the timescale of thesystem evolution [6, 7]. However, in many solid-state de-vices and other systems this assumption breaks down [8–12] so one cannot use simple time-local density matrixequations of motion. In addition to non-Markovianitybeing common, it can be desirable [13–24]: it has beenshown that non-Markovianity of open quantum systemscan lead to higher fidelity of quantum operations due tothe possible information backflow from the environmentto the system. The simulation of general non-Markovianopen quantum systems is, however, a computationallychallenging task, which hampers progress on the designof optimal control procedures.There exist a number of available methods applica-ble to simulating specific non-Markovian scenarios [25–37] (see [38] for a review of some of these). One ap-proach is to extend the notion of the system, in caseswhere the environment can be modeled as extra “noisequbits” which couple strongly to the system and weaklyto some Markovian environment [13, 17, 22]. This isdone systematically in the time evolving density oper-ator with orthogonal polynomials (TEDOPA) [28, 29]method, which maps the environment to a chain of cou-pled sites. Instead of augmenting the system space,one can write coupled hierarchical equations of motion(HEOM) [19, 20, 25, 37]; this performs well for spectraldensities that are well-approximated by a small numberof Lorentzians. Most relevant to this Letter are methods based on an augmented density tensor, capturing the sys-tem history, specifically the quasi adiabatic path integral(QUAPI) [26, 27]. When considering optimal control ap-plications, a major impediment to all of these methodsis that they are computationally intensive and the entirecalculation needs to be repeated for each set of controlparameters. This makes numerical optimization unfeasi-bly costly.In this Letter we present a general method to efficientlydesign quantum control procedures for non-Markovianproblems. A crucial step is a recasting of the TEMPOmethod [39, 40] within the framework of process tensorsintroduced by Pollock et al. in Ref. [41]. With the result-ing process tensor TEMPO (PT-TEMPO) method [42]one can perform the bulk of the computation indepen-dent of the system control parameters. This enables us torepeatedly find the system density matrix time evolutionfor various sets of control parameters at very low compu-tational cost. We can use this to optimize quantum con-trol procedures with respect to any chosen aspect of thesystem evolution, taking full account of non-Markovianeffects. To demonstrate this method we apply it to aquantum dot that is driven by a shaped laser pulse andstrongly coupled to a super-Ohmic phonon environment(see Fig. 1a for a sketch of the experimental setup wemodel). We show that our method can explore a thirty-five dimensional space of control parameters, and findoptimized pulses for state preparation in an ensemble offive qubits of differing energies.
TEMPO and process tensors — The most general sce-nario that we consider is a small system coupled to abosonic environment with a total Hamiltonian of the formˆ H = ˆ H S ( t, { p n } ) + ∞ (cid:88) k =0 (cid:104) ˆ O S (cid:16) g k ˆ b † k + g ∗ k ˆ b k (cid:17) + ω k ˆ b † k ˆ b k (cid:105) , (1)where ˆ H S is the system Hamiltonian, ˆ O S is the sys-tem coupling operator, g k are the coupling constants, ω k are the environment mode frequencies, and ˆ b † k and a r X i v : . [ qu a n t - ph ] J a n ˆ b k are the bosonic environment raising and loweringoperators. The coupling constants and environmentmode frequencies are encoded in the spectral density J ( ω ) = (cid:80) ∞ k =0 | g k | δ ( ω − ω k ). The system Hamiltonianˆ H S can be a function of time and of a set of control pa-rameters { p n } which parametrise the set of experimen-tally accessible system Hamiltonians. We assume that atsome initial time t the total state can be written in aproduct state ρ ( t ) = ρ S ( t ) ⊗ ρ E ( t ), where ρ E ( t ) is aGaussian state of the bosonic environment (for example athermal state). We make no assumption on the couplingstrength or the total state at any later time.The TEMPO method is numerically exact and basedon Feynman-Vernon influence functionals [43]. LikeQUAPI, the TEMPO method utilizes an augmented den-sity tensor (ADT) to encode the system’s history andits auto-correlations over time. It employs tensor net-work methods to compress this ADT in the form of amatrix product state (MPS) [44, 45], which allows it toinclude the history over hundreds of time steps. Fig-ure 1b exemplifies the TEMPO tensor network for fourtime steps. Each node of this network represents a ten-sor and each edge (called leg ) corresponds to an index.When a leg connects two tensors it signals a summationbetween them.The tensor network underlying the TEMPO methodworks in Liouville space, which means that densitymatrices are represented as vectors and the maps be-tween them as matrices (so called super-operators ).One needs to choose a small enough time step δt such that a Trotterization between the system Hamil-tonian ˆ H S ( t, { p n } ) and the rest is valid. The bluecircle in Fig. 1b is the vectorized initial system state ρ S ( t ) and the green circles are the system propaga-tors for a single time step U m = exp ( L S ( t m + δt/ δt )at times t m = t + mδt with the system Liouvillian ... pulsedlasergrating QDotphonon environment S L M system propagatorsinitial system stateprocess tensor (b)(a) (c) FIG. 1. (a) Sketch of the experimental setup to drive aquantum dot (QDot) with a shaped laser pulse. The pulseform can be modified with a spatial light modulator (SLM).(b) The TEMPO tensor network for four time steps. (c) Con-traction scheme to obtain the process tensor in MPS form. super-operator L S ( t ) = − i [ ˆ H S ( t ) , · ]. The blue squaresin Figs. 1b and 1c are the Feynman-Vernon influencefunctionals I n , which can be directly constructed formthe coupling operator ˆ O S , the spectral density J ( ω ) andthe initial environment state ρ E ( t ). An influence func-tional I n quantifies how the system evolution at sometime t m is influenced by the state of the system at theearlier time t m − n , thus allowing for a non-Markovian dy-namics of the system. Because the individual Feynman-Vernon influence functionals in the network are derivedfrom the corresponding time intervals in the environmentauto-correlation function, the TEMPO method performsbest when this function is smooth and decays to zerowithin some finite time. For more details on this method,see Refs. [39, 40].The crucial point we make use of is that since most ofthe TEMPO tensor network consists solely of influencefunctionals, which do not depend on the system Hamilto-nian or the initial system state (see red dashed region inFig. 1b), we may perform the bulk of the computation—contraction of the tensor network—before specifying thesystem Hamiltonian. This provides an efficient methodenabling optimization over Hamiltonian parameters.As realized by Jørgensen and Pollock the TEMPO net-work can be contracted to yield a process tensor [41, 46].The process tensor framework takes an operational ap-proach to characterize non-Markovian open quantum sys-tems. It has been used to resolve common misconceptionson implications of completely positive divisibility [47] andit gives a natural definition of quantum non-Markovianitythat coincides with the classical definition in the classi-cal limit [48]. Here we show that in addition to theseconceptional advantages it also has computational bene-fits. By constructing the process tensor without the sys-tem Hamiltonian, the formalism supposes a finite set oftime slots to perform system control operations. Thesecan be any completely positive map, for example, thepreparation of a particular state or a projective measure-ment. The process tensor has two legs for each timeslot and encodes the outcome of every possible sequenceof control operations. The red dashed area in Fig. 1bcan be identified as such a process tensor with respectto a total Hamiltonian excluding the pure system partˆ H (cid:48) = ˆ H − ˆ H S ( t, { p n } ). When we choose the time steps δt small enough such that a Trotterization between thesystem part and the rest is valid, it is also justified toput the pure system dynamics into the control sequenceas unitary evolutions. Thus, in the language of this for-malism, the red dashed area in Fig. 1b is a process tensorand the system propagators are control operations.To optimize the control of a non-Markovian open quan-tum system, we propose to make use of the ideas aboveand split the computation into two steps. First, we con-tract the influence functionals column by column as sug-gested in [46] and depicted schematically in Fig. 1c. Theresult of this computation yields the process tensor inMPS form, which we save. Then, we perform a searchby repeatedly inserting different system propagators as-sociated with various sets of control parameters. Thisallows us to compute the reduced system dynamics atvery low computational cost and thus enables us to findthe control parameters that optimize the fidelity of theprocess. Application to a quantum dot — We demonstrate theperformance of the PT-TEMPO approach by applying itto a quantum dot that is strongly coupled to its phononenvironment and driven by a configurable laser pulse.We aim to drive the quantum dot into a coherent su-perposition within a few picoseconds. Modelling this ischallenging because the evolution timescale is compara-ble to the memory time, so non-Markovian effects playan important role.Figure 1a shows a sketch of the experimental setup wemodel for this purpose. The laser pulse shape is modifiedwith a pair of diffraction gratings, lenses, and a spatiallight modulator (SLM). We consider the ground state andthe exciton state of the quantum dot and denote themwith |↓(cid:105) and |↑(cid:105) respectively. Under the rotating waveapproximation the system Hamiltonian (with ¯ h = 1) isˆ H S ( t ) = ω ↑↓ σ z + Ω( t )2 e − iω t ˆ σ + + Ω ∗ ( t )2 e iω t ˆ σ − , (2)where Ω( t ) is the positive frequency part of the classi-cal electrical field amplitude, ω is the laser carrier fre-quency and ω ↑↓ is the exciton energy. Also, σ z is thePauli matrix, σ + = |↑(cid:105)(cid:104)↓| , and σ − = |↓(cid:105)(cid:104)↑| . In addition,the quantum dot couples strongly to its phonon environ-ment with the coupling operator ˆ O S = σ z / J ( ω ) = 2 αω ω − c exp( − ω /ω c ),with the unit-less coupling constant α = 0 .
126 and thecut-off frequency ω c = 3 .
04 ps − [10, 12, 49]. The initialstate is assumed to be the product of the quantum dotground state |↓(cid:105) and the thermal state of the environ-ment at T = 1 K. We note that the environment auto-correlation function dies off only after about 2 . H S ( t ) = E ( t )2 ˆ σ + + E ∗ ( t )2 ˆ σ − , (3)where E ( t ) = Ω( t ) exp( − i ∆ t ) is the positive frequencypart of the electric field in the rotating frame, and∆ = ω − ω ↑↓ is the detuning of the carrier frequency ofthe input pulse with respect to resonance. The inputpulse (before it enters the pulse shaper) is assumed tobe Gaussian, i.e. E in ( t ) ∝ τ − exp (cid:0) − t /τ (cid:1) exp ( − i ∆ t ),with the input pulse duration τ assumed to be between30 fs and 300 fs. The pair of diffraction gratings and ap-propriate lenses enable the spatial separation of the fre- quency components of the input pulse, with an approx-imately linear relationship between frequency and posi-tion at the SLM. Therefore each SLM pixel modifies aparticular frequency range of the pulse. We further as-sume that the pulse also has a finite spatial width witha Gaussian profile, which results in a finite spot size foreach frequency part at the SLM. Given that each SLMpixel can induce a phase shift φ n to its correspondingfrequency Ω n the output pulse E ( t ) ∝ ( h ∗ E in ) ( t ) is aconvolution with the SLM’s impulse response function h ( t ) ∝ sinc (cid:18) δ Ω p t (cid:19) e − δ Ω2 st (cid:88) n ∈ pixels e i (Ω n t + φ n ) , (4)where δ Ω p is the pixel width and δ Ω s is the spot sizein terms of their corresponding frequency range [50, 51].We assume 512 SLM pixels centered at the pulse carrierfrequency and evenly spaced over a frequency range of2 π × . − . Also, we assume that the spot size ofthe pulse covers about two pixels, i.e. δ Ω s = 2 . × δ Ω p .The setup described here leads to 515 experimentalparameters to modify the laser pulse, namely, the initialpulse length τ , the initial pulse detuning ∆, the pulsearea Θ, and one phase shift φ n for each of the 512 SLMpixels. Instead of directly using the 512 parameters onthe SLM, we use a continuous phase mask function f ( x )on the interval x ∈ [ − , − φ n of pixel n is φ n = 2 π frac( f ( x ( n )) / π ) ∈ [0 , π ), where x ( n ) = ( n − /
256 and frac( y ) = y − (cid:98) y (cid:99) denotes thefractional part.To study the non-Markovian dynamics of the quan-tum dot as a function of these experimental parameterswe employ the PT-TEMPO method for which we firstcompute the process tensor. Similar to the conventionalTEMPO method, the accuracy of the result as well asthe necessary computation time depends on the choice ofthe simulation parameters. We choose time steps of 10 fs,a memory time of 2 . − . relative to the largest valueduring the contraction. With this, the computation ofthe process tensor takes approximately 167 s, which onlyneeds to be calculated once. The application of a systemHamiltonian to this process tensor takes only 1 . |↓(cid:105) to the | y + (cid:105) = ( |↑(cid:105) + i |↓(cid:105) ) / √ τ = 50 fs and the pulse area toΘ = 10 π . We also fix the shape of the phase mask func-tion to a downward facing parabola f ( x ) = Φ − x with a central phase shift Φ. This parabola results in abroadened and chirped output laser pulse that starts bluedetuned and ends red detuned with respect to its car-rier frequency. The central phase shift induces an overallphase which rotates the x-y coordinate system. Applyingour method we can easily map out the trace distance ofthe final state to the | y + (cid:105) target state, as a function of thetwo open parameters ∆ ∈ [ − ,
50] ps − and Φ ∈ [ − π, π ].Figure 2a shows the results of 201 ×
81 full non-Markoviansimulations corresponding to the different parameter sets.Employing the PT-TEMPO method the entire computa-tion takes less than 8 hours on a single core of an Intel i7processor, while it would take approximately 1040 hoursor 43 days employing the conventional TEMPO method.We find two local minima on this landscape which aremarked by a star and a diamond in Fig. 2a. The laserpulse that corresponds to the star parameter set is achirped pulse that starts strongly detuned and finisheson resonance. This can be thought of as an interruptedadiabatic rapid passage, which has the advantage to bealmost independent of an inaccurate pulse area, but hasthe disadvantage to be sensitive towards the detuning ofthe pulse. The laser pulse that corresponds to the dia-mond parameters, on the other hand, starts on resonanceand ends strongly red detuned. In this case, like a dy-namical gate, the fidelity of the protocol is sensitivelydependent on the pulse area, but tolerant towards de-tuning inaccuracies, similar to a simple π/ -40 -20 0 20 40 detuning /ps -202 c en t r a l pha s e s h i ft (a) -101 ( t ) (b) -0.5 0.0 0.5 1.0 time t /ps -50050 ( t ) / p s (c) (d) -0.5 0.0 0.5 1.0 time t /ps (e) x ( t ) y ( t ) z ( t )Re ( t )| ( t )| FIG. 2. The dynamics of a quantum dot as a function ofthe detuning and overall phase of a chirped laser pulse. (a) Aheat map indicating the trace distance of the final state tothe target state | y + (cid:105) . (b-e) Dynamics of the quantum dotand the electric field for the laser pulse parameters markedwith the symbols (cid:5) and (cid:63) in (a) respectively. to simultaneously drive them to the equator of the Blochsphere. The detunings of the quantum dots relative tothe middle dot are chosen to be [ − , − , , ,
10] ps − .We perform a global optimization search employing a dif-ferential evolution algorithm on 35 pulse parameters. Weparameterize the phase mask function by splitting it into32 segments and assigning one parameter to the slope ofeach segment. In addition to these 32 parameters we alsooptimize over all three input pulse parameters τ , ∆ andΘ. To avoid oscillations in time, the phase mask functionis smoothed out with a 3rd order spline. We expect thata short π/ π/ laserpulse fi vedetunedQDots FIG. 3. Optimization of a laser pulse driving an ensembleof five quantum dots (QDots). (a) A sketch of the ensembletaking the place of the single QDot in the setup from Fig. 1a.(b) The phase mask function for the π/ π/ • and ∗ respectively. The pulse length of both pulsesprior to the pulse shaper is 245 fs, the pulse areas for the ini-tial pulse and the optimized pulse are 0 . × π and 7 . × π respectively. The plots in (c) and (e) show the expectationvalues (cid:104) σ xy ( t ) (cid:105) = (cid:112) (cid:104) σ x ( t ) (cid:105) + (cid:104) σ y ( t ) (cid:105) for all five QDots. (RMS) distance to the equator of the Bloch sphere of0.10, which is significantly better than the performanceof a π/ τ = 245 fs (see Figs. 3c and 3d). Also, it performs slightlybetter than the shortest π/ τ = 30 fs, whichyields a RMS distance of 0.12. However, we note that,unlike the π/ Conclusion — We have shown that the PT-TEMPOmethod makes optimal control of non-Markovian openquantum systems a feasible task. It is applicable tosmall systems that couple to a structured bosonic en-vironment and it performs well if the environment cor-relation function is smooth and decays to zero withinsome finite time. The key idea of the method is to mod-ify the contraction order of the TEMPO tensor networksuch that the result of the bulk of the computation—corresponding to the contraction of the Feynman-Vernoninfluence functionals—can be stored and reused for eachnew trial system Hamiltonian. Finally, we note that thisidea can be also applied to other tensor network meth-ods [28, 29, 33, 34, 36], opening the way to a family ofefficient methods to design quantum control proceduresfor non-Markovian open quantum systems.We thank F. Pollock for insightful discussions on theprocess tensor framework. G.E.F. acknowledges sup-port from the EPSRC under grant no. EP/L015110/1.B.W.L. and J.K. acknowledge support from EPSRC(EP/T014032/1). E.B. acknowledges support from theIrish Research Council (GOIPG/2019/1871). P.R.E.acknowledges support from Science Foundation Ireland(15/IACA/3402). [1] M. Shapiro and P. Brumer,
Quantum Control of Molecu-lar Processes: Second Edition (Wiley-VCH Verlag GmbH& Co. KGaA, Weinheim, Germany, 2011).[2] E. Torrontegui, S. Ib´a˜nez, S. Mart´ınez-Garaot, M. Mod-ugno, A. del Campo, D. Gu´ery-Odelin, A. Ruschhaupt,X. Chen, and J. G. Muga, Chapter 2 - shortcuts to adia-baticity, in
Advances in Atomic, Molecular, and OpticalPhysics , Vol. 62, edited by E. Arimondo, P. R. Berman,and C. C. Lin (Academic Press, 2013) pp. 117 – 169.[3] S. J. Glaser, U. Boscain, T. Calarco, C. P. Koch,W. K¨ockenberger, R. Kosloff, I. Kuprov, B. Luy,S. Schirmer, T. Schulte-Herbr¨uggen, D. Sugny, and F. K.Wilhelm, Training Schr¨odinger’s cat: Quantum optimalcontrol: Strategic report on current status, visions andgoals for research in Europe, Eur. Phys. J. D , 279(2015).[4] C. P. Koch, Controlling open quantum systems: Tools,achievements, and limitations, J. Phys. Condens. Matter , 213001 (2016).[5] C. P. Koch, M. Lemeshko, and D. Sugny, Quantum con-trol of molecular rotation, Rev. Mod. Phys. , 035005(2019). [6] H.-P. Breuer and F. Petruccione, The Theory of OpenQuantum Systems (Oxford University Press, 2002).[7] L. Chirolli and G. Burkard, Decoherence in solid-statequbits, Adv. Phys. , 225 (2008).[8] I. Wilson-Rae and A. Imamo˘glu, Quantum dot cavity-QED in the presence of strong electron-phonon interac-tions, Phys. Rev. B , 235311 (2002).[9] C. Galland, A. H¨ogele, H. E. T¨ureci, and A. Imamo˘glu,Non-Markovian Decoherence of Localized Nanotube Ex-citons by Acoustic Phonons, Phys. Rev. Lett. ,067402 (2008).[10] A. J. Ramsay, T. M. Godden, S. J. Boyle, E. M. Gauger,A. Nazir, B. W. Lovett, A. M. Fox, and M. S. Skolnick,Phonon-induced Rabi-frequency renormalization of opti-cally driven single InGaAs/GaAs quantum dots, Phys.Rev. Lett. , 177402 (2010).[11] C. Roy and S. Hughes, Phonon-Dressed Mollow Tripletin the Regime of Cavity Quantum Electrodynam-ics: Excitation-Induced Dephasing and NonperturbativeCavity Feeding Effects, Phys. Rev. Lett. , 247403(2011).[12] S. L¨uker, K. Gawarecki, D. E. Reiter, A. Grodecka-Grad,V. M. Axt, P. Machnikowski, and T. Kuhn, Influenceof acoustic phonons on the optical control of quantumdots driven by adiabatic rapid passage, Phys. Rev. B ,121302 (2012).[13] P. Rebentrost, I. Serban, T. Schulte-Herbr¨uggen, andF. K. Wilhelm, Optimal control of a qubit coupled toa non-Markovian environment, Phys. Rev. Lett. ,090401 (2009).[14] R. Schmidt, A. Negretti, J. Ankerhold, T. Calarco, andJ. T. Stockburger, Optimal control of open quantumsystems: Cooperative effects of driving and dissipation,Phys. Rev. Lett. , 130404 (2011).[15] B. Hwang and H. S. Goan, Optimal control for non-Markovian open quantum systems, Phys. Rev. A ,032321 (2012).[16] F. F. Floether, P. De Fouquieres, and S. G. Schirmer,Robust quantum gates for open systems via optimal con-trol: Markovian versus non-Markovian dynamics, New J.Phys. , 073023 (2012).[17] D. M. Reich, N. Katz, and C. P. Koch, Exploiting Non-Markovianity for Quantum Control, Sci. Rep. , 12430(2015).[18] C. Addis, E. M. Laine, C. Gneiting, and S. Manis-calco, Problem of coherent control in non-Markovianopen quantum systems, Phys. Rev. A , 052117 (2016).[19] R. Puthumpally-Joseph, E. Mangaud, V. Chevet,M. Desouter-Lecomte, D. Sugny, and O. Atabek, Basicmechanisms in the laser control of non-Markovian dy-namics, Phys. Rev. A , 033411 (2018).[20] E. Mangaud, R. Puthumpally-Joseph, D. Sugny,C. Meier, O. Atabek, and M. Desouter-Lecomte, Non-markovianity in the optimal control of an open quan-tum system described by hierarchical equations of mo-tion, New J. Phys. , 043050 (2018).[21] M. H. Goerz and K. Jacobs, Efficient optimization ofstate preparation in quantum networks using quantumtrajectories, Quantum Sci. Technol. , 045005 (2018).[22] J. Fischer, D. Basilewitsch, C. P. Koch, and D. Sugny,Time-optimal control of the purification of a qubit incontact with a structured environment, Phys. Rev. A ,033410 (2019).[23] N. Mirkin, P. Poggi, and D. Wisniacki, Entangling proto- cols due to non-Markovian dynamics, Phys. Rev. A ,020301 (2019).[24] S. Alipour, A. Chenu, A. T. Rezakhani, and A. delCampo, Shortcuts to Adiabaticity in Driven Open Quan-tum Systems: Balanced Gain and Loss and Non-Markovian Evolution, Quantum , 336 (2020).[25] Y. Tanimura and R. Kubo, Time Evolution of a QuantumSystem in Contact with a Nearly Gaussian-MarkoffianNoise Bath, J. Phys. Soc. Japan , 101 (1989).[26] N. Makri and D. E. Makarov, Tensor propagator for itera-tive quantum time evolution of reduced density matrices.I. Theory, J. Chem. Phys. , 4600 (1995).[27] N. Makri and D. E. Makarov, Tensor propagator for it-erative quantum time evolution of reduced density ma-trices. II. Numerical methodology, J. Chem. Phys. ,4611 (1995).[28] J. Prior, A. W. Chin, S. F. Huelga, and M. B. Plenio,Efficient simulation of strong system-environment inter-actions, Phys. Rev. Lett. , 050404 (2010).[29] A. W. Chin, ´A. Rivas, S. F. Huelga, and M. B. Plenio, Ex-act mapping between system-reservoir quantum modelsand semi-infinite discrete chains using orthogonal poly-nomials, J. Math. Phys. , 092109 (2010).[30] J. Cerrillo and J. Cao, Non-Markovian dynamical maps:Numerical processing of open quantum trajectories,Phys. Rev. Lett. , 110401 (2014).[31] J. Iles-Smith, N. Lambert, and A. Nazir, Environmentaldynamics, correlations, and the emergence of noncanon-ical equilibrium states in open quantum systems, Phys.Rev. A , 032114 (2014).[32] D. Tamascelli, A. Smirne, S. F. Huelga, and M. B. Plenio,Nonperturbative Treatment of non-Markovian Dynamicsof Open Quantum Systems, Phys. Rev. Lett. , 30402(2018).[33] D. Tamascelli, A. Smirne, J. Lim, S. F. Huelga, and M. B.Plenio, Efficient Simulation of Finite-Temperature OpenQuantum Systems, Phys. Rev. Lett. , 90402 (2019).[34] A. D. Somoza, O. Marty, J. Lim, S. F. Huelga, and M. B.Plenio, Dissipation-Assisted Matrix Product Factoriza-tion, Phys. Rev. Lett. , 100502 (2019).[35] F. Mascherpa, A. Smirne, A. D. Somoza, P. Fern´andez-Acebal, S. Donadi, D. Tamascelli, S. F. Huelga, and M. B.Plenio, Optimized auxiliary oscillators for the simulationof general open quantum systems, Phys. Rev. A ,052108 (2020).[36] M. Brenes, J. J. Mendoza-Arenas, A. Purkayastha, M. T.Mitchison, S. R. Clark, and J. Goold, Tensor-NetworkMethod to Simulate Strongly Interacting Quantum Ther-mal Machines, Phys. Rev. X , 031040 (2020).[37] Y. Tanimura, Numerically “exact” approach to openquantum dynamics: The hierarchical equations of mo-tion (HEOM), J. Chem. Phys. , 020901 (2020).[38] I. De Vega and D. Alonso, Dynamics of non-Markovianopen quantum systems, Rev. Mod. Phys. , 015001(2017).[39] A. Strathearn, P. Kirton, D. Kilda, J. Keeling, and B. W.Lovett, Efficient non-Markovian quantum dynamics us-ing time-evolving matrix product operators, Nat. Com-mun. , 3322 (2018).[40] A. Strathearn, Modelling Non-Markovian Quantum Sys-tems Using Tensor Networks , Springer Theses (SpringerInternational Publishing, Cham, 2020).[41] F. A. Pollock, C. Rodr´ıguez-Rosario, T. Frauenheim, M. Paternostro, and K. Modi, Non-Markovian quantumprocesses: Complete framework and efficient characteri-zation, Phys. Rev. A , 012127 (2018).[42] The TEMPO collaboration, TimeEvolvingMPO: APython 3 package to efficiently compute non-Markovianopen quantum systems. (2020).[43] R. P. Feynman and F. L. Vernon, The theory of a gen-eral quantum system interacting with a linear dissipativesystem, Ann. of Phys. , 118 (1963).[44] R. Or´us, A practical introduction to tensor networks:Matrix product states and projected entangled pairstates, Ann. of Phys. , 117 (2014).[45] I. Cirac, D. Perez-Garcia, N. Schuch, and F. Verstraete,Matrix product states and projected entangled pairstates: Concepts, symmetries, and theorems, 2011.12127(2020), preprint.[46] M. R. Jørgensen and F. A. Pollock, A discrete memory-kernel for multi-time correlations in non-Markovianquantum processes, Phys. Rev. A , 052206 (2020).[47] S. Milz, M. S. Kim, F. A. Pollock, and K. Modi, Com-pletely Positive Divisibility Does Not Mean Markovian-ity, Phys. Rev. Lett. , 040401 (2019).[48] F. A. Pollock, C. Rodr´ıguez-Rosario, T. Frauenheim,M. Paternostro, and K. Modi, Operational Markov Con-dition for Quantum Processes, Phys. Rev. Lett. ,040405 (2018).[49] P. R. Eastham, A. O. Spracklen, and J. Keeling, Lind-blad theory of dynamical decoherence of quantum-dotexcitons, Phys. Rev. B , 195306 (2013).[50] A. M. Weiner, D. E. Leaird, J. S. Patel, and J. R. Wullert,Programmable Shaping of Femtosecond Optical Pulsesby Use of 128-Element Liquid Crystal Phase Modulator,IEEE J. Quantum Electron. , 908 (1992).[51] A. M. Weiner, Femtosecond pulse shaping using spatiallight modulators, Rev. Sci. Instrum.71