Modeling the unphysical pseudomode model with physical ensembles: simulation, mitigation, and restructuring of non-Markovian quantum noise

The influence of a Gaussian environment on a quantum system can be described by effectively replacing the continuum with a discrete set of ancillary quantum and classical degrees of freedom. This defines a pseudomode model which can be used to classically simulate the reduced system dynamics. Here, we consider an alternative point of view and analyze the potential benefits of an analog or digital quantum simulation of the pseudomode model itself. Superficially, such a direct experimental implementation is, in general, impossible due to the unphysical properties of the effective degrees of freedom involved. However, we show that the effects of the unphysical pseudomode model can still be reproduced using measurement results over an ensemble of physical systems involving ancillary harmonic modes and an optional stochastic driving field. This is done by introducing an extrapolation technique whose efficiency is limited by stability against imprecision in the measurement data. We examine how such a simulation would allow us to (i) perform accurate quantum simulation of the effects of complex non-perturbative and non-Markovian environments in regimes that are challenging for classical simulation, (ii) conversely, mitigate potential unwanted non-Markovian noise present in quantum devices, and (iii) restructure some of some of the properties of a given physical bath, such as its temperature.


I. INTRODUCTION
The definition of a physical system requires a distinction between internal and external degrees of freedom.The interaction between these constituents causes the internal, or closed, system dynamics to become open, i.e., affected by the external environment.For example, information stored in the system can propagate towards a measurement device or it can simply be lost, or dissipated, in the continuum of a thermal bath.Interestingly, this distinction between internal and external degrees of freedom becomes more ambiguous as the coupling between the system and the environment increases.In this case, information can be coherently exchanged between the system and the bath before the dynamics eventually stabilizes into a state characterizing hybridization properties of the whole system-environment.
An interesting characteristic of the pseudomode model (PM) is that, while its ancillary degrees of freedom superficially resemble simple physical harmonic modes and driving fields, in reality they can take on unphysical features and exhibit unphysical dynamics.This combination of superficial simplicity and unintuitive unphysicality is a consequence of the effective nature of the model which focuses on reproducing the reduced system dynamics after averaging over the ancillary degrees of freedom.This point of view can be more formally described as a map E phys → PM can between a physical environment E phys and its corresponding canonical pseudomode model PM can , whose unphysical parameters allow us to reduce the number of pseudomodes, and hence computational resources, needed to represent a given environment.In other words, the collection of these models directly correspond to real physical environments.In this context, pseudomodes have primarily been seen as a numerical tool that allows us to simulate non-Markovian and nonperturbative environments [91,92] in a transparent and simple way.Here, we enlarge the domain of applicabil-ity of the method by considering how pseudomodes can be used as ancilla degrees of freedom in experiments as a means to simulate, mitigate, or restructure physical environments.
To develop a framework encompassing all these applications, it is useful to consider pseudomode models PM with fully unconstrained parameters, i.e., corresponding to a map E PM → PM which generalizes the canonical case, PM can , to simulate a potential unphysical pseudoenvironment E PM .
To give meaning to this analysis, it is necessary to analyze what applications would be made available given the possibility to engineer these generalized environments in experiments (either via analog or digital techniques).
As an intuitive example, it is useful to consider what happens when such environments are coupled to the system alongside some originally present environmental bath E phys .We can then restrict the analysis to the situation in which the resulting 'total' environment E ′ phys is physical, i.e., original bath "restructured" (unphysical) (physical) (physical) (1) as exemplified in Fig. 1.This simple setting can be used to show different benefits of a possible simulation of the unphysical pseudomode environment E PM (henceforth dubbed pseudo-environment).Specifically, we will analyze three different applications.
• Quantum simulation.First, we note that, by imposing E phys = 0 in the equation above, we recover the canonical pseudomode model PM can in which the pseudo-environment directly correspond to a physical one.As a consequence, a protocol able to engineer the canonical pseudo-environment would correspond to performing a quantum simulation of non-Markovian quantum noise by using physical analogs or digital simulations [93][94][95][96][97][98][99][100][101][102].
• Mitigation of non-Markovian noise.Second, by imposing E ′ phys = 0, the equation above tells us that the pseudo-environment E PM is the one which exactly counteracts all the effects of the original E phys , resulting in the possibility to mitigate the effects non-Markovian noise generated by the physical environment E phys .As a consequence, in this case the simulation of such a pseudo-environment can be interpreted as an extrapolation/error-mitigation procedure on the lines of [103][104][105][106][107][108][109][110][111][112][113][114][115][116].The generality of these mitigation techniques usually rely on an underlying error model.For example, in [108], non-Markovian noise is analyzed in terms of Lindblad equations with negative rates, whose explicit derivation can be involved unless pertubative assumptions are made.In this context, the technique presented here can FIG. 1. Logical structure of the presented framework.The parameters of a pseudomode model, PM, can be extended to describe environments, EPM, which are physical, E ′ phys , when coupled to a system alongside another physical environment E phys as in (a).In the most general case, this pseudoenvironment EPM corresponds to an abstract unconstrained pseudomode model (b) which can be realized with a physical ensemble leading to a "restructuring" of the original bath E phys (c).More specific applications are found (blue ellipses) by imposing E phys = 0 (a) which lead to a simulation of the environment E ′ phys (b) through the canonical pseudomode model (c).Similarly (red ellipses), imposing E ′ phys = 0 (a), corresponds to a mitigation of the effects of the bath E phys (b) on the system through what we call "anti-modes" (c).
be interpreted as a non-Markovian mitigation protocol whose underlying error model is defined by the pseudomode mapping.On one hand, this has the advantage of being an explicit non-perturbative model but its applicability is, on the other hand, limited to noise originating from Gaussian Bosonic baths.
• General bath restructuring.Third, more generally, when E ′ phys is a modified version of E phys , the formalism can be interpreted as a way to restructure some environmental properties without assuming them to be experimentally accessible.As an example, E ′ phys could specify a new bath having the very same properties as the original E ′ phys apart from a different temperature.This restructuring can be realized when the pseudo-environment defined by Eq. ( 1) is coupled to the system alongside the original bath E phys .
The opening of the route towards these applications (for the simulation, mitigation, or reshaping of the properties of non-Markovian noise) is conditioned on the possibility to physically reproduce the unconstrained pseudo-environment.In this article we describe a protocol to enable the experimental realization of such unphysical pseudomode models.This is achieved by ana-lytical continuation of measurement results over an ensemble of physical systems [117], see Fig. 2. By using bounds provided by the theory of polynomial approximation and extrapolation, we describe the limitations of this technique in terms of an interplay between bias and stability errors.Mainly, while the order of the extrapolating polynomial allows one to increase the predictive power of the method under perfect conditions (reducing its bias), it also increases the sensitivity to imprecision in the initial data (observable outcomes over a physical ensemble).This leads to a limitation in the complexity of the effects of the environment on the system which can be recovered using this method.
This article is organized as follows.In section II we start by reviewing open quantum systems to introduce the pseudomode model in section III.In section IV, we describe the main result of this article, i.e., a protocol to simulate the pseudomode model using an ensemble of physical systems.In section V, we present an overview of the possible use of these results to simulate, mitigate, and restructure non-Markovian noise.In section VI, we provide specific numerical examples for each of these usecases for an environment described by a underdamped Brownian spectral density.In section VII, we further provide an error analysis of the algorithm.We finish presenting our conclusions in section VIII.

II. OPEN QUANTUM SYSTEMS
In this article, we focus on Gaussian bosonic environments, whose effects on a system can be fully characterized by the correlation of the operator through which they couple to the system.In this section, tildes are used over some of the parameters to explicitly distinguish them from conceptually related ones used to describe the effective baths introduced in the next sections.We consider a system S interacting with a physical bosonic environment E with a total Hamiltonian (ℏ = 1) where H S is the Hamiltonian of the system and H E = k ω k ã † k ãk is the free Hamiltonian of the environment with ãk the annihilation operator for the environmental mode k with frequency ω k.The interaction term is a function of a generic system operator ŝ and it has a linear dependency on the environment, i.e., where we introduced the coupling strengths λk ∈ R. The reduced dynamics of the system is where ρ S+E is the density matrix of the full sys-tem+environment.In general, the reduced dynamics of FIG. 2. Graphical representation of the protocol to simulate unphysical pseudo-environments.In (a), a system S interacts with a continuum of bosonic degrees of freedom.In (b), the reduced dynamics of the system is approximated by the one in which the system interacts with a discrete set of pseudomodes or stochastic fields, see Eq. (10).Some of these modes/fields can be unphysical (represented by the color red).In (c), an ensemble of physical pseudomodes/fields (represented by the color blue) leads to the same reduced system dynamics as the original pseudomode model in (b) and the original open quantum system in (a) once the parameter Λ parametrizing the ensemble is analytically continued to the critical value Λc.This prescription for the physical implementation of unconstrained pseudomode model opens an avaneue to the applications presented in Fig. 1 the system depends on all the n-point free correlation functions involving the interaction operator X(t) (where the time dependence indicates the Heisenberg picture).However, for Gaussian environments, all this information is encoded, through Wick's theorem, in the two-point correlation function where the average is taken over the equilibrium distribution ρ eq E of the environment.When the environment is initially in a thermal equilibrium at inverse temperature β, the correlation function is translational invariant and takes the form C E (t 2 , t 1 ) = C E (t 2 − t 1 ), where in terms of the spectral density function In the next section, we introduce the pseudomode model as a map between the reduced dynamics for the open quantum system described here and the one computed by averaging over the effects of a discrete set of ancillary harmonic quantum modes and classical stochastic fields.

III. PSEUDOMODE MODEL
The pseudomode method consists in replacing the original continuum of environmental modes with a discrete set of dissipative harmonic modes and stochastic driving fields.The main purpose of these ancillary degrees of freedom is to reproduce the correlation function characterizing the original Gaussian environment and, ultimately, to reproduce the original reduced system dynamics in Eq. ( 4).We refer to [84,85,90] and to Appendix A and B for more details.
The term "pseudomode model", PM, will be used to indicate a, possibly stochastic, operator M acting linearly on the density matrix ρ S-PM in a space composed by the system and N PM ancillary Fock spaces.The influence of M can thereby be interpreted as the effect of an artificial environment made out of N PM dissipative harmonic modes (pseudomodes) and an optional stochastic driving field acting on the system.Perhaps the most important feature of the model is the existence of a non-pertubative relation between these "pseudo-environments" and physical, Gaussian bosonic environments as described in section II.This "pseudomode mapping" defines the dependency of M on a set of parameters G PM which directly depend on the spectral density in Eq. ( 7), hence on the properties of the physical bath.This mapping is defined by matching the coupling statistics of the pseudomode model and the original bath [encoded in the correlation in Eq. ( 6)].
These definitions take a practical form in specifying the dynamics of the density matrix ρ S-PM of the pseudomode model as The parameters in G PM are optimized so that the reduced system dynamics is equivalent to the one in the original open quantum system, i.e., that Here, E indicates the expected value over the potentially, but not necessarily, non-trivial, stochastic properties of M.
It is useful to stop for a moment in order to introduce a few examples describing how to qualitatively different environmental effects can be modeled using the formalism above.

A. Reproducing environmental effects
The pseudomode model is designed to exactly reproduce all the effects of a Gassian bosonic bath linearly coupled to a system.For concreteness, here we analyze how to model specific effects.
Markovian effects.Markovian effects arise as a consequence of memory-less contributions C δ (t) = Γδ(t) to the correlation in Eq. ( 6), so that Γ can be interpreted as a decay rate.In this limit, the operator M takes a Lindblad form, i.e., it can be written as the sum of dissipators where Non-Markovian Classical effects.This case, which includes the previous one, amounts in modeling effects which can be ascribed to a symmetric contribution C class (t) = C class (|t|) to the correlation in Eq. ( 6).The corresponding model can be written in terms of a stochastic drive where the coefficients {c n }, n = 1, • • • , N ξ , define the spectral representation of a Gaussian field with zero bias and correlation matching C class (t), i.e., which can always be achieved by increasing the cut-off parameter N ξ .Non-Markovian Quantum effects.This case, which includes both previous ones, amounts in modeling general effects related to a non-classical contribution C Q (t) to the correlation in Eq. ( 6).The corresponding model can be written in terms of a sum of terms taking the form written in terms of a pseudomode operator a and dependent on its frequency Ω, the coupling strength λ, the decay rate Γ, and the distribution n.This form implies the effects due to N PM pseuomodes to be encoded in a correlation contribution written as where The pseudomode mapping consists in finding the optimal value for these parameters such that, in general, By considering all these different effects together, the canonical pseudomode model associated with the environment E is determined by the collection where k = 1, • • • , N PM , n = 0, • • • , N ξ , characterizing the environmental correlation in Eq. ( 17) in terms of pseudomodes and fields (where we omitted the Markovian component as its effects can always be included in the other degrees of freedom).
In [85], it was shown that, the correlation of environments characterized by a Brownian spectral density in the underdamped regime, can be well approximated by a sum of purely quantum contributions, even in the zerotemperature case.In this case, some of the modes were shown to require imaginary couplings λ to the system.In [90], it was further shown that such an environment can always be modeled by a single, physical, quantum contribution with n = 0 alongside a classical one representing a single, imaginary, stochastic field, for any temperature of the original bath.These results imply that, in order to correctly reproduce the reduced system dynamics, the artificial environment representing the pseudomode model might require to have "unphysical" parameters.
While preventing a direct physical interpretation, this enlarged parameter domain also implies the possibility for a more optimized description of the original environment.At the same time, it is worth to explicitly point out that, since the form of the dynamical equation does not depend on the physicality of the parameters, its unphysical solutions can be interpreted as the analytical continuation of the physical solutions (considered as functions of the variables in the set G PM ).For example, non-Hermitian contributions to the unitary dynamics are not compensated by taking the Hermitian conjugate on the Hamiltonian when acting on the right of the density matrix [85].This makes the procedure qualitatively different than the orthodox concept of non-Hermitian quantum mechanics [118,119].
In order to give a more formal terminology, we define a parameter η ∈ G PM as physical if η ∈ R + and unphysical otherwise.This language sets the basis for the generalizations which we are going to analyze in this article.

B. Beyond the canonical pseuodomode model
The pseudomodel described above is labeled "canonical" because, despite its possible unphysicality, all the parameters are constrained to mimic the effects of physical environments on the system through Eq. (10).Here, we are interested in a more general case whose parameters are unconstrained to lie outside the canonical model, thereby corresponding to unphysical environments.Given a set of parameters G PM of an unconstrained pseudomode model, we define E GPM as its corresponding pseudo-environment which is made out of the (potentially unphysical) modes and fields characterized by the parameters in G PM , through Eq. ( 18).We will show that, while unphysical, these environments can be used to effectively change, or "restructure" the properties of a physical bath.Given this general setting, our goal is to describe how to reproduce the effects of such an unconstrained pseudomode model from measurement results over physical ensembles.We will achieve this by defining an analytical continuation protocol on the unphysical parameters to allow observable outcomes in a physical ensemble to reproduce the effects of any general unphysical pseudomode model.In turn, this will lead to a physical representation of models associated to physical environments (thereby defining a tool for their simulation) and models corresponding to more general effects, such as noise mitigation or the restructuring of some environmental properties such as temperature.

IV. UNPHYSICAL PSEUDOMODES WITH PHYSICAL ENSEMBLES
In this section, we analyze a physical characterization of the unconstrained pseudomode model.The possibility to do so stems from the fact that the analytical continuation considered here can intuitively be interpreted as a Wick's rotation (used to map time to an imaginary temperature-like quantity) on the parameter Λ.In this context, [117] (see also [120]) have shown that observable results on an ensemble of quantum systems can be used to perform a Wick rotation thereby opening an avenue to analyze classical statistical systems from quantum measurements and vice-versa.
Here, we implement a similar strategy to analytically continue observables extracted from a physical pseudomode model into observables of the corresponding unphysical pseudomode model.To start, we denote the unphysical/physical parameters as η unphys j /η phys j , where k = 1, . . ., N phys and j = 1, . . ., N unphys , specifying their total number.These definitions can be used to write the set G PM characterizing the model as We then consider a set of N unphys functions Ξ j : Here, the first set of conditions allow to interpret the values Ξ j (Λ) as physical parameters on the domain D phys which, here, is arbitrarily chosen to be The second condition requires that the unphysical parameters η unphys j can be recovered by analytical continuation of the functions Ξ j (Λ) for Λ → Λ c ∈ C. While not necessary, it is convenient to further impose to reflect bounds on accessible physical regimes.The parameter Λ is the one we will, eventually analytically continue.We now define a "physical regularization" of the set G PM as which uses Λ to parametrize an ensemble of pseudomodes models which are physical in the domain Λ ∈ [−1, 1].
Using the second condition in Eq. ( 21), it is possible to directly verify that This shows that, indeed, the analytical continuation of the physically regularized model recovers the unconstrained pseuodomode model.For concreteness, it is possible to define the function Ξ[Λ] as which, together with fulfills the constraints in Eq. ( 21) and Eq. ( 22) through the definition where and where is defined in terms of the step function θ(x)=1 for x ≥ 1 and zero otherwise.While other choices are possible, this will be the default one used throughout this article.
In summary, we formally defined physical versions of the pseudomodes parametrized by a single parameter Λ whose analytical continuation to Λ → Λ c reproduces the unphysical pseudomodes model as shown in Eq. (24).Our attention now shifts towards the implementation of this analytical continuation using measurement results over the physical ensemble G phys PM (Λ) corresponding to sweeping over different values of Λ.

A. Analytical continuation
In the previous sections we showed that, given a pseudomode model defined by a set G PM (Λ) of unconstrained parameters, there exists a regularized version G phys PM (Λ) whose parameters are physical for Λ ∈ D phys .This implies that the corresponding models M phys PM (Λ) can be interpreted as a physical ensemble parametrized by Λ ∈ D phys .Here, the word "physical" takes an operative meaning [121], as each of these models can, in principle, be experimentally realized by coupling the system to a set of ancillary resonators and by driving it with a classical stochastic field.This means that, for such Λs, and for a given time t, we can always perform a sequence of measurements to reconstruct the density matrix which simply corresponds to the dynamics of the model as in Eq. ( 8).The functional dependence of the density matrix in the parameter Λ can then be analytical continued over to the value Λ c to finally achieve our goal, i.e., the simulation of the unphysical pseudomdode model corresponding to the density matrix This result can also be rephrased in terms of equivalent pseudo-environments.In fact, generalizing Eq. ( 19), we can define E phys PM (Λ) as the physical pseudo-environment corresponding to the regularized set G phys PM (Λ), i.e., In this way, Eq. ( 30) can be interpreted as an analytical continuation of the corresponding pseudo-environments, i.e., While this can be considered just as an equivalent way to interpret Eq. ( 30), it also allows us to grasp its meaning from a different point of view.In fact, Eq. ( 32) can be intuitive interpreted as manifesting the possibility to define a parametrization for physical environments such that, upon analytical continuation, the effects on the system are equivalent to those generated by a generalized, unphysical pseudomode model.This different prospective paves the way for the applications which will be analyzed in section V.
In practice, we can only assume the physical ensemble E phys PM (Λ) to be realized for a discrete grid of points, which is the case we are going to analyze in the next section.

B. Polynomial extrapolation
In this section, we analyze the analytical continuation of the expectation values of a system observable over the physical ensemble described by a set of We note that this is a weaker version of the case considered in the previous section.In fact, the analytical continuation of the full density matrix can be considered as the limiting, and more resource-consuming, task where we analytically continue not just one, but a complete set of observables.We now introduce the following practical procedure for analytical continuation.We first define the function p M (Λ) as the order M polynomial which minimizes the distance This least square fitting can be used to reconstruct the observable as The ability to reconstruct the full pseudomode model from an ensemble of pseudomodes constrained to physical dynamics opens up the possibility of several interesting applications that would be, otherwise, physically inaccessible.In the next section, we are going to analyze these applications, in terms of simulation, restructuring, and mitigation of non-Markovian quantum noise.

V. APPLICATIONS
In this section we describe different possible applications of our protocol.Mainly, we show that, by driving a system with classical noise and by coupling it to ancillary harmonic modes, it is possible to analogically simulate, mitigate, or restructure the effects of a Gaussian Bosonic environment.

A. Simulation
We start by defining a procedure to simulate a physical environment E ′ , which we assume to be well described by a pseudomode model characterized by a parameter set G ′ PM .We note that this corresponds to imposing E phys = 0 in Eq. ( 1), which then translates to In other words, a simulation of the open quantum system S +E ′ can be achieved by coupling the (closed) system to the pseudo-environment E PM corresponding to the pseudomode model with parameters G ′ PM .As noted in our previous discussion, the pseudo-environment E PM is, in general, not directly physically realizable.To circumvent this problem, we follow the analysis in section IV, and define a set of regularized parameters G phys PM (Λ) corresponding to a collection of physical pseudo-environments E phys PM (Λ).By measuring observables in this ensemble, it is then possible to proceed with the analytical continuation protocol illustrated in section IV, i.e., which corresponds to simulating the effects of the original environment E ′ on the system.

B. Mitigation
Another interesting application is the use of the pseudomode model to mitigate the effects of a non-Markovian environment E, which we assume to be well described by a pseudomode model characterized by a parameter set The mitigation procedure correspond to imposing E ′ phys = 0 in Eq. ( 1), which brings it to the form This equation describes the (unphysical) pseudoenvironment E PM which, once coupled to the system, completely cancels all effects of the environment E. Since we assume the knowledge of the pseudomode model describing E, the parameters defining this "antimode" environment E PM can be explicitly written as To derive this equation we simply noted that, if the original set G PM describes a correlation C E (t), then the set since Eq.(A12) and Eq. ( 14) are linear in c n and g 2 k , respectively.In other words, the compound effect of the pseudomode model and its "antimode" version does not have any influence on the system since the two correlations sum up to zero: A system simultaneously in contact with a Gaussian environment and its unphysical "mirror" or anti-environment should evolve as if no environment were present at all, i.e., noise-free.Because of its effective noise-cancelling action, the antimode environment is necessarily unphysical.However, following section IV, we can regularize its parameters to define the corresponding physical ensemble E phys PM (Λ) parametrized by Λ, from which gathering all the information needed for the analytical continuation to achieve mitigation using measurements on a physical ensemble.

C. Restructuring
In its most general form, Eq. ( 1) reads and can, for example, be interpreted as a "restructuring" of a given physical environment E phys (which we assume associated to a pseudo-environment phys ) characterized by different physical properties (such as temperature or system-bath coupling).This is achieved by coupling the system to a pseudoenvironment E restructuring PM whose parameters must be a function of those in E PM and E PM ′ encoding the specific "changes" to be imposed on the bath.For the most interesting cases, the resulting environment is going to be unphysical so that, as done in the previous cases, we will need to resort to the analytical continuation procedure outlined in section IV for the physical implementation of its effects.
Here, we note that, while more optimized versions might be possible (as shown in the example in section VI B), in the worse case, the restructuring can always be defined by using both the simulation and mitigation techniques presented in the previous sections.In fact, thanks to our assumptions and the definition in Eq. (38) we can directly check that satisfies Eq. ( 42).In other words, in order to restructure a bath E phys into E ′ phys it is always possible to couple the system to both the anti-environment relative to E phys and an additional one simulating E ′ phys .As mentioned, this worse case scenario can be optimized depending on the specific requirements of the restructuring as we will show in the example in section VI B.
In the following, we are going to present explicit numerical examples for each of these applications.

VI. NUMERICAL EXAMPLES
To show a practical numerical implementation [122][123][124][125] of the applications presented in the previous section, here we present three examples: the mitigation and simulation of the effects of an environment characterized by a underdamped Brownian spectral density at zero and finite temperature and the "restructuring" of the finite temperature case into the zero temperature one.We note that, for exposition clarity, the order in which these examples are presented is different with respect to the previously reported one.
For concreteness, we will focus on a Gaussian Bosonic environment characterized by the spectral density function written in terms of a resonance frequency ω 0 , a frequencywidth γ, and a (frequency) 3/2 strength λ.It describes a Ohmic behavior at low frequency, i.e., and it has a polynomial cut-off at high-frequencies, i.e., in terms of the adimensional scale α = λ 2 γ/ω 4 0 .We further restrict to the underdamped regime which requires By inserting the spectral density J B (t) in Eq. ( 6) the correlation function can be computed as [90] where , is the symmetric, or "classical" contribution which contains all temperature effects which can, thereby, be modeled within the statistics of a single temperature dependent, classical field ξ β (t), defined explicitly in appendix A, and such that where t = t 2 − t 1 .We can label the remaining contribution as "quantum", since it does not have specific symmetries, but it has the form of Eq. (A10), i.e., it can be reproduced using a single harmonic mode initially at zero temperature.
To be specific, we will consider the open system to be made out of a two-level system coupled to the Brownian environment B described above.Specifically, given a full density matrix ρ S+B (t) in the system+environment space whose dynamics is described by a full Hamiltonian H S+B as we want to compute the reduced dynamics encoded in By defining the system Hamiltonian as and assuming that the operator ŝ = σ x mediates the interaction to the environment, the reduced dynamics can be computed by using the pseudomode model, i.e., by computing which requires to solve the differential equation Here, ρ S+PM is the density matrix in the sys-tem+pseudomode space.This dynamics is driven by a single field ξ β (t) whose statistics depends on the inverse temperature β to reproduce the Matsubara contribution to the correlation in Eq. ( 64) and a single pseudomode a res to describe the resonant properties of the spectral density in Eq. (44).Following the standard pseudomode mapping [84,85,90], we have in terms of the Hamiltonian the zero-temperature dissipator written for a generic operator a , and the coefficients The initial state is ρ S+PM (0) = ρ S (0) ⊗ |0⟩ ⟨0| in terms of the vacuum |0⟩ annihilated by the harmonic mode a.
An interesting feature of this hybrid model associated to the spectral density in Eq. ( 44) is that [90]: • all the parameters associated with the resonant mode a res are physical, i.e., λ res , Ω res , Γ res are real and positive.
• all temperature effects are encoded in the statistics of the field ξ β (t); i.e., in its autocorrelation function.
Given this specific environment, we are now going to analyze three specific applications in terms of mitigation, restructuring, and analog simulation.

A. Noise mitigation
In this section, we are going to analyze the mitigation of the effects of the environment specified by the spectral density in Eq. ( 44) both at zero and finite temperature to ultimately compute observables corresponding to the noise-free dynamics To do this, we add stochastic driving and coupling to ancillary quantum modes to the original S + B open quantum system such that, after analytical continuation of a single parameter, their correlation is exactly the opposite of the original one C B (t).We can achieve this using the formalism described in section V B. First, we want to define the unphysical "antimode" model PM whose correlation satisfies Using Eq. ( 39) this can be done by introducing a resonant "anti-mode" āres and an "anti-field" ξβ , whose parameters are the same as those associated with a res and ξ β in Eq. ( 56), except for an additional complex-rotation in the interaction to the system.Specifically, where the field coefficients c n are explicitly defined in Eq.A14.
To continue, we need to define a physical ensemble whose effects are the same as the antimode model after analytical continuation.Before doing this, it is worth differentiating between the zero and finite temperature case.

Zero Temperature
In the zero-temperature limit, it is interesting to note that a simplification occurs in the classical expression for the correlation function as [85,90]: This zero-temperature "Matsubara" contribution is negative at zero time, hinting to the fact that imaginary fields are necessary to ensure the correct sign is reproduced in Eq. ( 50).This is, in fact, the case, and the field ξ β=∞ (t) is purely imaginary; which, using Eq. ( 63), corresponds to a real anti-field ξβ=∞ (t).Therefore, the stochastic part of the model does not require introducing any additional analytical continuation procedure [see Eq. ( 27)].Thus, the analytical continuation needs only to be performed upon a single unphysical parameter (the coupling λ res between the system and the resonant mode).Given this construction of the antimode model, we can now define its physically-regularized version, i.e., we want to couple the system to a physical ensemble which can be analytically continued to the antimode model above.Specifically, given the Hamiltonian H S+B in the original system+bath space, we couple the system to a physical regularized-antimode āphys res and drive it by a physical field ξphys β=∞ (t), such that the dynamics is described by Here, the physicality of the field ξphys β=∞ (t) = ξI β=∞ (t) explicitly corresponds to having an autocorrelation function of the form In parallel, the parameters characterizing the physicallyregularized antimode āphys res are the same as those for the antimode āres in Eq. ( 56), except for the presence of the additional parameter Λ [introduced through the function Ξ defined in Eq. ( 25)], renormalizing the coupling to the system into the physical domain.We assume such a coupling to be constrained in the range [0, λ res ], corresponding to Λ ∈ [−1, 1].The correlation function of this mode is then given by which, by construction, gives rise to a minus sign for Λ = Λ c = −1 + 2i using Eq.(25).In this way, Eq. ( 66) and Eq. ( 67) fulfill the identity in Eq. ( 62) defining the antimode model.In other words, by performing the analytical continuation Λ → Λ c , the model in Eq. ( 65) adds a correlation −C B res (t) to the open quantum system in the (S + B) space in Eq. ( 52), i.e., it completely counteracts, in principle, all environmental noise ultimately leading to the free dynamics in Eq. (61).
We give a graphical exemplification of this proceudure in Fig. 3.It is important to note that, in order to compute the dynamics plotted in this figure, we did not solve the differential equation in Eq. ( 65) which involves the original environmental continuum.Instead, the physically-regularized antimode model is introduced on top of the deterministic pseudomode model characterizing such a continuum.Specifically, we considered which corresponds to using the antimode model alongside the deterministic pseudomode model in [85] [used to describe the continuum B in Eq. ( 65)].We refer to Appendix B 3 for a brief overview.The solutions of this differential equation for different values of Λ in Λ n ∈ [−1, 1] are used to compute the observable ⟨ ÔS ⟩ phys (t; Λ).This is done by interpolating the values ⟨ ÔS ⟩ phys (t; Λ n ) with an order M polynomial p M (t; Λ) which is then analytically continued to Λ → Λ c to compute ⟨ ÔS ⟩(Λ c ; t) = ⟨ ÔS ⟩ free (t) corresponding to Eq. ( 61).To improve the accuracy of the reconstruction, we also perform a regularization procedure by analytically continuing a complete set of observables (in this case, the three Pauli matrices) and imposing physicality on the extrapolated result.Specifically, the regularization we considered here corresponds to writing where In Fig. 3(a) we show the original dissipative dynamics alongside the ones determined by coupling additional noise (i.e., the regularized antimode environment) to the system.The values of these observables are then used to recover the noise-free dynamics.In Fig. 3(b) the underlying analytical continuation procedure is shown explicitly in the complex plane for a specific time of the dynamics.
As shown in Fig. 3(a), the effect of the further regularization in Eq. ( 69) helps in reproducing the correct values for the noise-free dynamics.
We present a further exemplification in Fig. 4, where we plot the fidelity of a simple single-qubit gate, consisting of a π/2 rotation around the y-axis in the Bloch sphere.The qubit is in contact with a zero-temperature bath which lowers the quality of the gate as the gatetime is increased.We show that the effects of the bath on this single-qubit operation can be mitigated through analytical continuation.

Finite Temperature
As mentioned in the introduction to this section, the effects of a bath with spectral density given by Eq. ( 44) in the underdamped regime can be modeled using a single resonant harmonic mode and a classical stochastic field even at finite temperature.In fact, all temperature effects of the bath are encoded in the statistics of this driving field.As a consequence, in order to mitigate temperature effects, we only need to update the way we handle this classical stochastic process with respect to the (zero temperature) analysis in the previous section.Such a case was, in fact, rather special as the stochastic driving field in the antimode model was real, thereby not requiring regularization and analytical continuation.This is no longer the case at finite temperature where, in order to satisfy its defining Eq. ( 50), the field ξ β (t) is no longer restricted to imaginary values but must, in general, be written as in terms of real and imaginary parts.Here, the notation simply follows from using the definition ξβ (t) = ξR β (t)+i ξI β (t) for the antifield in Eq. ( 63).The mitigation of the imaginary part follows the same procedure as in the zero-temperature case, i.e., it only requires the introduction of a corresponding real anti-field.However, the mitigation of the real part now requires to be regularized and analytically continued, similarly to what was done before for the system-pseudomode coupling.In turn, this leads to the following modification to the regularized antimode model in Eq. ( 65) as in terms of the physically regularized version of the antifield which, following Eq.( 27), reads These equations define the ensemble whose physical observables can be used to mitigate the noise of the environment B at finite temperature.In Fig. 5, we analyze a specific example to highlight this model.Restructuring of a finite temperature environment into a zero temperature one.The black/blue curves correspond to the dissipative dynamics when the system is in contact with a finite/zero temperature environment.For the dashed red curve, a stochastic field is added to effectively approximate the zero-temperature dynamics after analytical continuation.The specific parameters used here are: ωz = ω0, ∆ = 0, λ = 0.1ω 3/2 0 , γ = 0.3ω0, and β = 1/ω0.

B. Restructuring the environment
In this section, we consider a scenario in which one might wish to modify, i.e., restructure, some properties of an environment without assuming them to be directly accessible.We specifically focus on the case in which the temperature of an environment is algorithmically reduced by using the analytical continuation procedure presented in the previous sections.
Interestingly, in the underdamped Brownian bath considered here, the modeling of this situation does not require the introduction of any additional quantum degree of freedom.Indeed, it is possible to algorithmically modify the temperature of the original bath by simply performing analytical continuation over the intensity of a classical stochastic drive.In fact, we can explicitly write the finite temperature stochastic pseudomdomode model as whose field satisfies Eq. (50).Therefore, the differential equation for a pseudomode model at a different temperature β ′ can be obtained by simply adding an extra field ξ ∆β (t) to the previous equation with Given these considerations, the temperature β of a bath B acting on a system S can be algorithmically modified The blue curve shows the true dissipative dynamics, in dashed-red the simulated one obtained by coupling the system to a harmonic mode and by analytically continuing a classical stochastic drive.The dashed black curve shows, the simulation done without the additional field, which is necessary to predict the correct hybridization to the environment, when the system is in contact with a zero temperature environment.The dashed red curve shows the dynamics after analytical continuation is performed upon an additional stochastic field.In the inset, a restricted scale of the main plot is shown to better highlight the effects of the simulation.The specific parameters used here are: ωs = ω0, ∆ = 0, λ = 0.3ω to a new temperature β ′ by driving the system with a regularized stochastic field To be more explicit, the dynamics in the full (S + B) space becomes so that, assuming the original dynamics to be well described by the model in Eq. ( 74), the analytical continuation to Λ c effectively changes the classical correlation function to C B class (t; β ′ ); thereby achieving the mentioned algorithmic change in temperature.We give an example of this procedure in Fig. 6.

C. Simulation
Perhaps the most immediate application of the formalism presented here is in the analog simulation of the non-Makovian effects of an environment.This can be implemented in a rather direct way since the resonant pseudomode is physical.Given a closed system with Hamiltonian H S , the effects of a Brownian environment at inverse Classical spectrum corresponding to an underdamped Brownian spectral density as a function of frequency.At low temperatures, the correlation is characterized by the Matsubara contribution leading to a negative spectrum.In this regime, imaginary fields are required for the simulation of this enviornment.At higher temperatures, resonant terms (vertical lines highlight ω = ±ω0) dominate, leading to a positive spectrum.In this regime, the classical contribution to the correlation can be simulated using physical (real) fields.In the inset, we show, in arbitrary units, the integral of the spectrum as a function of inverse temperature β to highlight the cross-over (as the curve intersects zero) between the two regimes.Other parameters used here are: γ = 0.1ω0, and temperature β can be reconstructed from the physical ensemble obtained from Eq. ( 74) with the replacement It is interesting to note that this expression is the "dual" of Eq. ( 73).In fact, while Eq. ( 73) is used to mitigate the classical effects of the environment, Eq. ( 78) is used to simulate them.We present a specific example of this analog procedure in Fig. 7, where we simulate the effects of a zero temperature environment.As mentioned above, the only obstacle preventing from a direct simulation of an environment is the imaginary component of the field ξ(t).It is then worth asking whether there are regimes where such component is negligible.To do this, we can recall that the stochastic fields depend on the square-root of the parameters c n defining the spectral decomposition of the classical contribution to the bath correlation function, see Eq. (A13).Assuming the time dynamics to be the largest time-scale in the model, these coefficients become which can be interpreted as a continuum version of Eq. (A14) and written in terms of the "classical" spectrum of the bath In the limit, the question about the possibility to simulate a bath using physical fields then becomes equivalent to check whether this quantity is positive, i.e., S class [ω] ∈ R + for all ω.For the Brownian spectral density considered here, the expression for the classical contribution to the correlation is given explicit in Eq. (B58) in terms of a sum of decaying exponentials, corresponding to a spectrum represented as an infinite sum of Lorentzians functions.In Fig. 8, we plot the spectrum for a specific set of parameters.From this figure we can recover the previously analyzed case: at low temperatures, the classical correlation is determined by the Matsubara contribution to the correlation which is negative, corresponding to imaginary fields.However, at higher temperatures, resonant terms in the correlation become dominant leading to a positive spectrum.This argument can be used (see Appendix B 3) to derive an estimate for the inverse temperature locating the cross-over between the two regimes as For β < β * , it is possible to produce a direct analogical simulation of non-Markovian effects at higher temperatures, such as in the case of the excitonic energy transfer in molecular dimer systems [64,[126][127][128][129].In Fig. 9, we simulate one of the environments analyzed in [129].

VII. ERROR ANALYSIS
The limitations of the method we presented rely on the interplay between imperfections in the experimental data and the complexity of the environmental effects to be mitigated.
It is possible to make this intuitive consideration more precise by estimating the error in the computed observables after analytical continuation.In Appendix C, we follow [130] and use the properties of Chebyshev polynomials (see, for example, [131][132][133]) to present precise upper bounds on the expected error.In general, the expected difference between the true observable following a noise-free dynamics and the one reconstructed through analytical continuation can be written as Err = Err bias + Err stability . ( This expression depends on the sum of two qualitatively different terms characterizing the stability and the bias of the algorithm.In the case of perfect knowledge on the measured observables, the stability term is zero.In this case, the error is only due to the bias, i.e., by how precisely the M -order fitting polynomial is able to reconstruct the analytical continuation.In general, the bias becomes exponentially accurate by increasing M .Unfortunately, sensitivity to imperfections in the initial data does also increase exponentially with the order of the polynomial.This is encoded in the term Err stability , which is also inversely proportional to the square-root of the number N exp of data available. .Dynamics of a two-level system corresponding to a model for excitonic energy transfer dynamics in a molecular dimer system [129] at room temperature.The blue curve shows the results simulated by the HEOM and the dashed red curve shows the ones arising from a quantum simulation protocol in which the system is coupled to a single harmonic mode at zero temperature and a classical field.The high-temperature regime allows this field to be physical, thereby not requiring any analytical continuation procedure.The parameters were chosen to be on the same order as the ones used in Fig. 5 in [129] and they explicitly are: ω0 = 220 cm −1 , ωs = 100 cm −1 , ∆ = 100 cm −1 , γ = 20 cm −1 , α = 2λ 2 /ω 2 0 = 80 cm −1 , T = 300 K [corresponding to 1/β = kBT = (kBT )λ0/(2πℏc)ℏω0 = 0.94ℏω0, in terms of the Boltzmann constant kB, the speed of light c, the wavelength λ0 = (1/220) cm] in units such that ℏ = 2πc = 1.For consistency, in this simulation, the coupling operator to the environment is set to ŝ = σz.
This interplay between stability and bias is a common feature to extrapolation algorithm and it can be interpreted as the balance between underfitting and overfitting.For example, in Fig. 10 we computed the error as a function of the degree of the interpolating polynomial showing the origin of an optimal range for the degree M .

VIII. CONCLUSIONS
We presented an extrapolation technique to mitigate, restructure, and simulate the effects of Gaussian non-Markovian environments on a quantum system.The method relies on the interaction of the system with ancillary leaky modes and a stochastic driving field to define a physical ensemble parameterized by a single parameter.Measurements results over this ensemble can be used to define an analytical continuation procedure, allowing us to use these ancilla modes to perform one of the following tasks: • Simulating the effects of a non-Markovian environment.
• Mitigating non-Markovian noise affecting the system.Reconstruction error as a function of the order M of the interpolating polynomial.The error refers to the open quantum system analyzed in Fig. 3.More precisely, we defined Err = |⟨σz(t)⟩ free − ⟨σz⟩ reg (Λc; t)| where the averages are taken with respect to ρ free (t) and ρ reg (t), respectively, and the time corresponds to the grey vertical line in Fig. 3.For small values of M the error is biased, i.e., the polynomial is not able to encode the complexity of the environmental effects we want to mitigate.For high values of M , the algorithm is unstable with respect to errors in the initial data.In this example, the artificially injected error is modeled as a normal variable with zero mean and standard deviation σ = 10 −5 .As the variance of the error decreases, the optimal parameter range of the algorithm (shaded blue) extends towards the right, allowing to model more subtle features.
• Restructuring some of the properties of a given environment (already interacting with the system) without directly accessing it.
We presented the details and proof of the general formalism, as well as several numerical examples to showcase the flexibility of the algorithm to adapt to different environments (such as zero and finite temperature ones, and a physically motivated example related to simulating excitonic energy transfer in a molecular dimer systems).These examples also demonstrated how the above range of applications can be used in practise (such as noise mitigation for the dynamics of quantum observables or quantum gate operations, simulation of zero-temperature environments, and restructuring of a finite temperature bath to a zero-temperature one).The algorithm used in all these tasks is based on polynomial extrapolation and it is thereby ultimately limited by instability against imperfections in experimental data.
As an outlook, the protocol presented here could also be adapted to specify or optimize digital quantum algorithms for the mitigation of non-Markovian noise in quantum computing tasks, or for the simulation of zerotemperature environments for ground-state engineering, see [134].In this context, it would be relevant to relax any requirement on prior knowledge about the environment by preceding the presented protocols with noise characterization techniques such as in [116].We further note that the regularization procedure used in here could be extended to bigger system sizes by analytically continuing a reconstructed version of the full state using classical shadow tomography [135][136][137].In addition, the possibility of restructuring a given environment could, possibly, be helpful to push experimental setups past the limits imposed by the presence of physical baths elusive of direct manipulation.a. Introduce an open quantum system S +E ′ in terms of an alternative environment E ′ , such that the interaction operator with the system has a correlation given by , with the C k PM (t) defined in Eq. (A10).The environment E ′ is made by (N PM + 1) pseudomodes and (N + 1) bosonic baths for each of them.The reduced system dynamics of this model is equivalent to that of the original one, provided b. Show that reduced dynamics in the Sys-tem+Pseudomode space of the open quantum system S + E ′ above is equivalent to the one originating from the Lindblad Eq. (A1).
c.As a consistency check, we show that this Lindblad model has correlations given by Eq. (A9).
The next three sections, labelled a, b, c will prove the three points itemized above.

a. Pseudomodes as an open quantum system
Let us consider an open quantum system in which a subsystem of interest S interacts with an alternative environment E ′ made out of a set of (N + 1) pseudomodes a k , each coupled to its own bosonic bath (whose modes are denoted as b kα ) as where the Hamiltonian of the environment H B is with (B3) We also defined the interaction operator We assume that the environment of each pseudomode is populated by modes with both positive and negative frequencies, i.e., ω kα ∈ (−∞, +∞).We further suppose that the spectral density associated with the bath of the kth pseudomode is constant, i.e., and that each of the pseudomodes a k and the corresponding environmental modes b kα are in an initial state ρ 0 /Z kα respectively (with Z k and Z kα imposing unit trace), together with the consistency condition which, in the continuum limit, reads We explicitly highlight the slight abuse in notation as β k (ω) denotes the inverse temperature associated with environmental mode at frequency ω for the pseudomode k, while β k denotes the inverse temperature associated with the pseudomode k).We note that these unorthodox conditions define a state for the full environment which is not a thermal state but it is the closest quantum idealization of classical white noise (see [2], pag.164).
Here, we defined and Z 0 as a constant to impose unit trace.The free correlation function of the interaction operator can be obtained as where X PM (t) = e iH B t X PM e −iH B t .The unorthodox definition of the equilibrium state of the bath is designed to allow which, in turn, makes the correlation translational invariant in time since (B12) Note that this is non-trivial because the Hamiltonian H B involves a Jaynes-Cummings interaction between the pseudomodes and the modes of their baths.Because of this, the previous relation does not hold in the presence of a true thermal equilibrium [2].In order to make progress evaluating C PM (t), we can first compute the formal solution for the Heisemberg equation of motion ḃkα = i[H B , b kα ], and use it in the Heisemberg equation of motion for the pseudomodes ȧk = i[H B , a k ].This leads to the following result for the Laplace transforms xk , pk of the quadratures where (B14) Now, using the expression for the spectral density in Eq. (B5) we find which leads to We can now use this result in the expression for the translational invariant correlation Eq. (B12) to obtain where we defined ⟨•⟩ ≡ Tr PM (• k ρ k ), the trace being over the pseudomodes space.We note that, in the above derivation, translational invariance in time (derived thanks to the condition in Eq. (B6)) was essential as it allowed a great simplification through the iden-  The dynamics in the System+Pseudomodes space can be explicitly written in terms of the following influence functional expression where T is the time-ordering operator and where the influence superoperator is where (using the shorthand (B21) where we defined B23) where 2n k (ω) + 1 = coth β k (ω)ω/2 and, importantly, we used the condition in Eq. (B6) to obtain 2n k (ω) + 1 = coth β k Ω k /2, effectively implying Using Eq. (B23) into Eq.(B20), we get B25) where we used t 0 dt ′ δ(t−t ′ ) = 1/2 (see Eq. 5.3.12 in [2]).Going back to the Shroedinger picture, using Eq.(B19), and taking a time derivative, we get where and This is Eq.(A1) in the main text.

c. Correlations for the dissipative pseudomodes
As a consistency check, we now show that the correlation function computed in Eq. (B47) can, equivalently, also be computed using the Lindblad model developed in the previous section.Specifically, we want to show that where with (B33) which implies Eq. (B29).Here, we prove that where = 0 , (B36) where, in the last step, we used Eq.(B6).

A hybrid pseudomode model
As shown in [90], it is possible replace some of the quantum degrees of freedom present in the fully deterministic pseudomode model presented in the previous section with a classical stochastic colored noise ξ(t) which is stationary, Gaussian, and with zero mean.In fact, the addition of a driving term ξ(t)ŝ in the Hamiltonian in Eq. (B2) effectively adds, after averaging over the noise, a term where t 2 − t 1 = t and which accounts for the statistics of the field.In turn, the dynamics of this hybrid model can be explicitly written by adding the stochastic noise in Eq. (B26) to write the Lindblad operator as While we refer to [90] for more details, here we note that the addition of noise in the model corresponds to a decomposition of the effects of the original environment into a classical and a quantum part, i.e., in which the second term is symmetric under time reversal (so that it can be modeled using classical resources) while the first term is more general and requires ancillary quantum degrees of freedom.Interestingly, this decomposition is not unique allowing for further possibility of optimization in the model.For example, it is possible to chose the classical contribution C class (t) [hence the field ξ(t)] in such a way that all the pseudomodes modeling C Q (t) are, initially, at zero temperature.This is the choice which is used throughout this article.

Underdamped Brownian Spectral Density
In this section, we describe the explicit form of the pseudomode model for the spectral density J B (t) in Eq. (44).At zero temperature, the deterministic pseudomode mode can be defined using a single resonant mode at N mats zero-frequency Matsubara modes.This follows from the decomposition in terms of the resonant and Matsubara contributions .
(B41) It is possible to model these terms by introducing a resonant and N Mats Matsubara harmonic modes a res , a k (k = 1, • • • , N Mats ) and define a pseudomode model through the Lindblad operator where , with These parameters are explicitly given by The remaining contribution to the Lindbladian can be written as in terms of with These parameters can be estimated by fitting the corresponding correlation (B47) to M (t), see [85].However, because of the absence of oscillatory behavior in M (t), it is possible to directly impose It is interesting to note that, in general, the number of Matsubara modes needed to reproduce the effects of M (t) depends on the simulation time.In fact, we note that the domain of integration in the integral in Eq. (B41) is effectively restricted to values x ≪ 1/t by the presence of the exponential [138].Therefore, This shows that M (t) is characterized by an asymptotic (∼ 1/t 2 ) polynomial decay which contrasts with the exponential decay present in the Matsubara correlation, see Eq. (B47).At finite temperatures, the decomposition of the correlation becomes where In order to reproduce it, we can introduce three "resonant" modes a j;res , j = 1, 2, 3 alongside N mats "Matsubara" modes a k , k = 1, • • • N mats characterized by the Lindbladian The Matsuabara contribution is formally the same as in the zero temperature case, but its parameters are defined to fit M (t; β) instead of its zero-temperature limit in Eq. (B41).However, the resonant part has to be updated and defined by where B54) as a function of the parameters B55) characterizing the resonant modes through the definitions We note that other choices which only use two resonant modes initially prepared with complex-temperature values are possible [139].
The hybrid pseudomode model relies on writing the correlation function of the bath in the following alternative form where B58) Since the "classical" contribution C B class (t) is symmetric under time reversal, it can be reproduced by a single classical stochastic process ξ B (t) with autocorrelation function The remaining "quantum" contribution can be reproduced with a single "resonant" pseudomode characterized by the same L res Lindblad as the one defined in the zero temperature case.Despite the presence of a series dependent on the Matsubara frequencies, the spectrum S[ω] can be computed in a closed form.In fact, we can write (B60) and in terms of the Fourier transforms of the quantum and classical contributions to the correlation Assuming continuity and J(0) = 0, the zero-frequency limit becomes where J * (ω) ≡ J(ω)/ω which we assume to be defined at ω = 0 by continuity.For the Brownian spectral density above, we have which allows to write the classical spectrum in a closed form as where We can further write the zero-frequency limit in the following explicit form We now present an intuitive estimation for the temperature above which this spectrum is positive.To do this, we note that, for low temperatures, the classical spectrum is negative (due to the fact that the classical correlation is dominated by the Matsubara contribution) and peaked at zero frequency while, at high temperature, the spectrum becomes positive over all frequencies.As a consequence of this consideration, there must be a cross-over temperature β * such that S class [0] = 0. Using the results above, this temperature can be computed explicitly.In fact, we have By imposing S B class [0] = 0, the value β * can be written as available data f exp n .As explained in section IV B, we are ultimately interested in reconstructing the expectation ⟨ ÔS ⟩ reconstructed corresponding to an analytical continuation at a specific unphysical value denoted by Λ c ∈ C, i.e., In the following, we further limit the analysis to the case N ≥ M , so that the number of data points available is bigger than the order of the extrapolating polynomial.The goal of the following section is to provide an estimate for the extrapolation error of the quantity in Eq. (C4), i.e., where ⟨ ÔS ⟩(Λ c ) is the "true" analytically-continued value which we ultimately want to compute.Importantly, the average E exp takes into account the uncertainties on the experimental data encoded in the random variables ϵ n .
To analytically approach the estimation of this error, in the next section we follow [130] and take advantage of an explicit extrapolation form, expressed in terms of Chebyshev polynomials.

Chebyshev Polynomials
This subsection reviews the definitions and properties of the Chebyshev polynomials of first kind, adapting the references [131][132][133].
The Chebyshev polynomials of first kind are defined as for m ≥ 0 and for x ∈ [−1, 1], such that x = cos θ, which immediately implies the bound While not immediately apparent from the definition, these functions satisfy the recurrence relation which, together with the initial conditions T 0 (x) = 1 and T 1 (x) = x imply that each T m (x) is a polynomial of order m.It is further possible to analytically continue these polynomials and to write them as C9) where z ∈ C. Given a parameter ρ > 0, it is possible to define the so-called Bernstein ellipse Given a point z ∈ C is very useful to define the following function i.e., such that z ∈ Γ ρz or, in other words, the Bernstein ellipse with parameter ρ z passes through the point z.
In the following, we will denote by Ẽρ the open Bernstein ellipse, i.e., the region of the complex plane inside E ρ (i.e., the one containing the origin).This constitutes a rather interesting identity because it implies that, on the Bernstein ellipse E ρ the Chebyshev polynomials satisfy the bounds and, since ρ > 1, the cleaner, but less strict one, The Chebyshev polynomials T (x) are a basis for functions in the interval [−1, 1].This means that every function g(x) can be written [131], in this interval, in terms of the Chebyshev series It is further possible to follow an elegant geometric construction [131], to show that an analytical continuation of g(x) in Ẽρ corresponds to analytically continue G(z) inside the annulus ρ −1 ≤ |z| ≤ ρ.Intuitively, this is a consequence of the fact that G(z) can be interpreted as the pull-back of g(x) on the unit circle by the Joukoasky map J(z) = (z + z −1 )/2, z ∈ C (which projects the unit circle onto the real axis, i.e., J(z) = Re[z] for |z| = 1).In other words, F (z) = g(x).In turn, this allows to take advantage of the Laurent series in Eq. (C15) to write [131]  (C19) We note that this bound is valid for all parameters ρ > 1 defining an open Bernstein ellipse Ẽρ inside which g(x) admits an analytical continuation.
In the following subsection, we show how to use these results in order to provide an estimate for the error made in a polynomial extrapolation.

Error analysis using Chebyshev Polynomials
In this subsection, we adapt the results elegantly presented in [130] to analyze in more detail the expression for the extrapolation error in Eq. (C5) which averages the difference between the "true" extrapolation value f (Λ c ) = ⟨ ÔS ⟩(Λ c ) and its M -order polynomial estimate p exp M (Λ c ).Assuming that the function f (Λ) admits analytical continuation, we can always write it in terms of a Chebyshev series as In this way, the coefficients a m (m = 0, . . ., M ) can be written explicitly as a m ≡ ⃗ a m , where Here, we introduced the vector which, as can be noted by comparison, is not equivalent to the one in Eq. (C23) since it relies on the truncated series rather than the "true" value.With this notation, we have (C30) Our goal is now to compute the quantities in the expressions above.To begin, using Eq.(D6), we have with the intention to find an analytical bound.To make progress, we follow [130], as done throughout the whole section.It is possible to define a projector P = T (T † T ) −1 T † onto the range of T which has the property (T † T ) −1 T † P = (T † T ) −1 T † .(C47) This means that we can write C48) where, in the second step, we used Eq.(D8), and where we defined We also used E exp [X] 2 ≤ E exp [X 2 ] for a generic random variable X (special case of Jensen's inequality), in the last step.Now, since P projects onto the range R T of T , we can write exp E exp || ⃗ v∈R T 1 = (M + 1)σ 2 exp , (C50) where ⃗ v ∈ R T are orthonormal vectors and where we used a notation borrowed from quantum mechanics to indicate vectors and their duals as kets and bras, respectively.Note that, in the last step, we assumed that M < N , i.e., that the number of measurements is bigger than the degree of the interpolating polynomial, as it is the regime relevant for us.In that case, since T is a (N +1)×(M +1) matrix, its rank cannot be bigger than (M +1), justifying the last step.Using Eq. (C50) into Eq.(C48), we obtain The definitions above regarding the vector spaces C n can be used to define induced norms in the vector space of operators A : C q → C p as ||A|| n = sup It is interesting to explicitly compute ||A|| 2 .To do this, we notice that A can be represented by a p × q matrix A pq which, using singular value decomposition, can be written as A p×q = V p×p Σ p×q U q×q in terms of unitary matrices V and U and an upper diagonal Σ.We have where |A| = √ A † A. For future reference, we also write here the following identities: with |r| < 1 for the second identity to hold.

Least Squares
We want to define the best estimate for a vector ⃗ c ∈ C q such that its image under a linear mapping T : C q → C p is as close as possible to a given vector ⃗ f ∈ C p under the || • || 2 norm.Equivalently, we want to find the vector ⃗ c which minimizes the function We have which allows to write The equations ∂ cα L = 0 and ∂ c * α L = 0 are equivalent and lead to When T † T is invertible, this leads to which is the explicit expression for the vector ⃗ c minimizing the least squares distance L in Eq. (D12).

FIG. 3 .
FIG. 3.Mitigation of an underdamped Brownian environment at zero temperature.(a) Dynamics of the observable σz.In blue, the original dissipative dynamics, solving Eq.(56).A physical version of the antimode model is added, to generate the different violet curves corresponding to different Λ ∈ [−1, 1] in Eq.(68).In green, partial and full mitigation as an analytical continuation on the function σz(t;Λ) for Λ ∈ [−1, −1 − 2i].In red, a regularization is applied to the analytically continuation of the full density matrix.(b) Analytical continuation at the specific time described by the grey line in (a).Violet points corresponds to the measured observables for different real Λ.The green points correspond to the partially and fully implemented analytical continuation.In blue, the original value in the solely presence of the bath B. The overall grey surface is the polynomial pM (t; Λ) used to define the analytical continuation as an extrapolation.The specific parameters used are: α = 0.02, Γ = 0.3ω0, ωs = ∆ = ω0, M = 10, Nexp = 12, N ξ = 100.
FIG. 6.Restructuring of a finite temperature environment into a zero temperature one.The black/blue curves correspond to the dissipative dynamics when the system is in contact with a finite/zero temperature environment.For the dashed red curve, a stochastic field is added to effectively approximate the zero-temperature dynamics after analytical continuation.The specific parameters used here are: ωz = ω0, ∆ = 0, λ = 0.1ω

8 FIG. 7 .
FIG.7.Analog simulation of a zero temperature environment.The blue curve shows the true dissipative dynamics, in dashed-red the simulated one obtained by coupling the system to a harmonic mode and by analytically continuing a classical stochastic drive.The dashed black curve shows, the simulation done without the additional field, which is necessary to predict the correct hybridization to the environment, when the system is in contact with a zero temperature environment.The dashed red curve shows the dynamics after analytical continuation is performed upon an additional stochastic field.In the inset, a restricted scale of the main plot is shown to better highlight the effects of the simulation.The specific parameters used here are: ωs = ω0, ∆ = 0, λ = 0.3ω
FIG.9.Dynamics of a two-level system corresponding to a model for excitonic energy transfer dynamics in a molecular dimer system[129] at room temperature.The blue curve shows the results simulated by the HEOM and the dashed red curve shows the ones arising from a quantum simulation protocol in which the system is coupled to a single harmonic mode at zero temperature and a classical field.The high-temperature regime allows this field to be physical, thereby not requiring any analytical continuation procedure.The parameters were chosen to be on the same order as the ones used in Fig.5in[129] and they explicitly are: ω0 = 220 cm −1 , ωs = 100 cm −1 , ∆ = 100 cm −1 , γ = 20 cm −1 , α = 2λ 2 /ω 2 0 = 80 cm −1 , T = 300 K [corresponding to 1/β = kBT = (kBT )λ0/(2πℏc)ℏω0 = 0.94ℏω0, in terms of the Boltzmann constant kB, the speed of light c, the wavelength λ0 = (1/220) cm] in units such that ℏ = 2πc = 1.For consistency, in this simulation, the coupling operator to the environment is set to ŝ = σz.
FIG. 10.Reconstruction error as a function of the order M of the interpolating polynomial.The error refers to the open quantum system analyzed in Fig.3.More precisely, we defined Err = |⟨σz(t)⟩ free − ⟨σz⟩ reg (Λc; t)| where the averages are taken with respect to ρ free (t) and ρ reg (t), respectively, and the time corresponds to the grey vertical line in Fig.3.For small values of M the error is biased, i.e., the polynomial is not able to encode the complexity of the environmental effects we want to mitigate.For high values of M , the algorithm is unstable with respect to errors in the initial data.In this example, the artificially injected error is modeled as a normal variable with zero mean and standard deviation σ = 10 −5 .As the variance of the error decreases, the optimal parameter range of the algorithm (shaded blue) extends towards the right, allowing to model more subtle features.
. acknowledges support from NSFC (Grant No. 11935012) and NSAF (Grant No. U2330401).F.N. is supported in part by: Nippon Telegraph and Telephone Corporation (NTT) Research, the Japan Science and Technology Agency (JST) [via the Quantum Leap Flagship Program (Q-LEAP), and the Moonshot R&D Grant Number JPMJMS2061], and the Asian Office of Aerospace Research and Development (AOARD) (via Grant No. FA2386-20-1-4069).
B18) as in Eq. (A10) in the main text upon the definition g k = λ 2 k /2Ω k .Since, by hypothesis, this is a Gaussian open quantum system, the reduced dynamics of the system is fully determined by the functional form of C PM (t) through Dyson equation.As a consequence, this model reproduces the same dynamics as the original model as long as C E (t) = C PM (t).
C16)where the factor 1/2 should not appear on the left-hand side for m = 0.By further assuming |g(z)| ≤ Q ρ for x ∈ Ẽρ the expression above implies|g 0 | ≤ Q ρ , and |g m | ≤ 2Q ρ ρ −m for m ≥ 1 ,(C17)see theorem 8.1 in[131].We can now use these bounds to analyze the error made by truncating the Chebyshev series in Eq. (C15).To do this, we can define g M (x) as the truncated seriesg M (x) = M m=0 g m T m (x) .(C18)Using Eq. (C7) and (C17), the error made by this truncation can be quantified as sup x∈[−1,1] |g(x) − g M (x)| = sup