Assessment of weak-coupling approximations on a driven two-level system under dissipation

The standard weak-coupling approximations associated to open quantum systems have been extensively used in the description of a two-level quantum system, qubit, subjected to relatively weak dissipation compared with the qubit frequency. However, recent progress in the experimental implementations of controlled quantum systems with increased levels of on-demand engineered dissipation has motivated precision studies in parameter regimes that question the validity of the approximations, especially in the presence of time-dependent drive fields. In this paper, we address the precision of weak-coupling approximations by studying a driven qubit through the numerically exact and non-perturbative method known as the stochastic Liouville-von Neumann equation with dissipation. By considering weak drive fields and a cold Ohmic environment with a high cutoff frequency, we use the Markovian Lindblad master equation as a point of comparison for the SLED method and study the influence of the bath-induced energy shift on the qubit dynamics. We also propose a metric that may be used in experiments to map the regime of validity of the Lindblad equation in predicting the steady state of the driven qubit. In addition, we study signatures of the well-known Mollow triplet and observe its meltdown owing to dissipation in an experimentally feasible parameter regime of circuit electrodynamics. Besides shedding light on the practical limitations of the Lindblad equation, we expect our results to inspire future experimental research on engineered open quantum systems, the accurate modeling of which may benefit from non-perturbative methods.


Introduction
Driven quantum systems are ubiquitous in quantum technologies. They appear, for example, in the control and measurement protocols as well as in the studies of non-equilibrium dynamics [1,2]. One of the simplest paradigmatic examples encompasses a two-level quantum system, a qubit, subjected to a classical drive field which promotes population dynamics in the eigenbasis of the bare qubit. Despite its simplicity, such a model has been applied in many contexts ranging from the coherent control in quantum computing to the simulation of a number of important photochemical reactions [3][4][5][6][7][8][9][10][11][12]. Moreover, its properties have been investigated through different descriptions such as the dressed and Floquet state formalisms [13][14][15][16][17][18], also being associated with various physical phenomena, such as coherent suppression of tunneling [14,19,20] and interference between successive Landau-Zener transitions [21][22][23][24][25].
Another well-known example of drive-induced quantum phenomena is attributed to the work by Mollow in Ref. [26]. Remarkably, Mollow theoretically showed that the fluorescence spectrum of a driven qubit may turn into a triplet in the presence of weak dissipation. If the Rabi frequency of the classical field well exceeds the dissipation rate, two sidebands emerge in the spectrum with an offset equal to the Rabi frequency from the center peak at the drive frequency. A more sophisticated explanation of such a phenomenon was later provided using a quantum treatment also for the drive field [27][28][29]. In this so-called dressed-state picture, the energy levels of the composite qubit-field system are split due to the dynamic Stark effect promoted by the drive. The Mollow triplet has been verified experimentally in many different physical scenarios [30][31][32][33][34][35][36][37][38].
The approach for solving the open-quantum-system dynamics in Mollow's study [26] and in the follow-up work in Refs. [27][28][29] assumes a weak coupling between the system, i.e., the qubit, and its bath of quantized bosonic modes, thus motivating a perturbative treatment of the dissipation [39][40][41]. Such an approach is guided by the so-called Born-Markov approximations, where one assumes a stationary bath and fast decay of bath correlations in the typical timescales of the system evolution. Furthermore, upon the elimination of fastoscillating terms, the typically non-unitary evolution of the system is usually expressed by a Lindblad master equation (LME) [42,43] with positive decay, excitation, and dephasing rates, in addition to which the environment introduces a rescaling of the transition frequency of the free system. However, one of the possible drawbacks of such a form of the LME is that the interplay between the drive and dissipation is not fully contemplated. As a consequence, the presence of a time-dependent drive field may rise questions on the validity of the above-mentioned approximations. In the literature, a vast amount of strategies have been presented to approach such a scenario, each with their own range of applicability and assumptions motivated by the details of the system under study. For example, still within and non-perturbative method has been proposed to capture more general initial systembath states, including correlated ones [89]. In this paper, however, we restrict our studies to the case of factorized initial states in such a way that the use of SLED is well justified. In principle, one has the choice to prepare such a factorized state in the beginning of the dynamics.
In this context, the goal of this work is the following: Firstly, we investigate the main characteristics introduced by SLED on the properties of the driven qubit by considering the usual LME [39,40] as a point of comparison. To this end, we show that the often overlooked bath-induced energy shift term in the LME may give an important contribution to the dynamics and greatly alter the actual steady state of the driven qubit, in stark contrast to non-driven systems with otherwise matching parameter values. Even though the relevance of the bath-induced energy shift has already been pointed out in other contexts, e.g. [90,91], our study illustrates the regimes where such effects are more pronounced, implying that they should be carefully considered in the related experiments. Taking advantage of the nonperturbative characteristics of SLED, we also propose a scheme to experimentally witness the failure of the asymptotic predictions of the LME. Secondly, we study signatures of the Mollow triplet using the SLED in a pump-probe spectroscopy configuration with parameter regimes that may be implemented in the framework of cQED.
This manuscript is organized as follows. In section 2, we present the theoretical and numerical models, emphasizing the main differences between the Lindblad and SLED master equations for a driven system. In section 3, we apply both methods in the case of a monochromatic transverse drive field and show that the overlap fidelity between the SLED and Lindblad solutions is drastically reduced if the bath-induced energy shift term is not taken into account in the LME. Furthermore, we analyse the steady-state properties of the driven system and describe a protocol to witness the failure of the LME. The concern here is not to prove that the approximate treatment with the LME fails, but instead to show the practical relevance of non-perturbative approaches such as the one we utilize in this work. In section 4, we numerically study the meltdown of the Mollow triplet within the SLED and Lindblad formalisms for experimentally feasible parameters for cQED. Section 5 concludes this work.
Some particular choices ofĤ d (t) have been historically used in the study of driven quantum systems and in the discovery of novel physical phenomena. For example, with f (t) = Ω 0 + Ω d cos(ω d t), the dynamics of the expectation value ofσ x may be frozen in the fast drive regime ω d ω q for selected values of Ω d and ω d even when Ω 0 = 0. This phenomenon is usually referred to as coherent destruction of tunneling [14]. For f (t) ∝ t, the transition probability between the lower eigenstate of the instantaneous Hamiltonian H S +Ĥ d (t) at t → −∞ and its excited eigenstate at t → +∞ can be found analytically, a case generally referred to as Landau-Zener transitions [21]. In contexts where the drive field couples weakly to the qubit, such as in laser-atom interactions, and its frequency lies near the bare qubit frequency ω q , one may apply the rotating-wave approximation and writeĤ d (t) ≈ Ω d (|0 1|e iω d t + |1 0|e −iω d t )/2, as in the original investigation of the qubit fluorescence spectrum in Mollow's work [26].
For the background theory considered in this section, we do not consider any specific form of f (t). In sections 3 and 4, we consider the case of oscillating transverse fields without a constant term. The only assumption made here is that such a classical field is a good approximation of a coherent quantum field over the time scales of interest [92][93][94]. This restricts the subsequent analysis to a system described by a two-dimensional Hilbert space, thus reducing computational time of the numerically exact protocol.
We choose the qubit to interact linearly with a dissipative bosonic bath which is modeled by an infinite set of quantum harmonic oscillators. The j:th oscillator has creation and annihilation operatorsb † j andb j , respectively, and an angular frequency ω j . The Hamiltonian of the bare bath can thus be written asĤ B = j ω jb † jb j , and the system-bath interaction Hamiltonian asĤ SB = σ x j g j (b j +b † j ), with {g j } being the corresponding coupling strengths. The total Hamiltonian is given bŷ Conveniently, the system-bath interaction can be fully characterized by the spectral density function where δ(ω − ω j ) is Dirac delta function. In the continuum limit, the spectral density becomes a smooth function of ω that approaches zero as ω → ∞. This so-called ultraviolet cutoff is physically motivated in cQED, for example, by the finite bandwidth of the transmission lines coupled to the qubit. Specifically, we use an Ohmic environment with a quartic Drude cutoff characterized by the spectral density [76,77] J(ω) = 2ηω where η is an effective dimensionless coupling constant and ω c is the bath cutoff frequency, which is assumed to be much higher than the qubit bare frequency. The quartic cutoff has been chosen over the typical quadratic Drude and exponential cutoffs in order to speed up the convergence of the numerical method since it produces a conveniently narrow spectrum for a given cutoff frequency ω c . Nevertheless, this spectral density can capture the relevant physics of a tunable resistor coupled to a superconducting transmon qubit [76] shown in figure 1.
In this work, we assume that the total density operator is initially factorized, i.e., being the Gibbs state of the bath at temperature T = (k B β) −1 and with mean excitation numbern(ω) = (e βω − 1) −1 . Consequently, the internal dynamics of the heat bath can be represented by its autocorrelation function whereξ = j g j (b j +b † j ) and S(ω) = J(ω)[n(ω) + 1] is the power spectrum of the noise. The behavior of L(τ ) also plays an important role in the reduced dynamics of the qubit, with its real and imaginary parts, denoted hereafter by L r (τ ) and L i (τ ), respectively.
Motivated by experiments, for example those in nuclear magnetic resonance [95][96][97], a perturbative description of the reduced qubit dynamics under the Hamiltonian in equation (1) and the choice ofρ(0) has been extensively studied in the weak-coupling case, i.e., for |g j /ω q | 1 [39,40,[98][99][100]. Such an approach relies on assuming that the system-bath correlations created dynamically are negligible for the dynamics of the system so that in the reduced master equation of the system, we may use the thermal equilibrium state of the environment, i.e.,ρ B (t) ≈ρ B (0). In addition, memory effects on the reduced dynamics arising from the finite decay time of L(τ ) are usually neglected. These assumptions constitute the Born-Markov approximations, and along with the elimination of quickly oscillating terms, produce a time-local master equation for the reduced density operatorρ S (t) = Tr B [ρ(t)] of the qubit which can be cast into a Lindblad form as (see Appendix A for details) where the superoperators D ij (ρ) = |i j|ρ|j i| − {|j j|,ρ}/2 express incoherent transitions between the eigenstates ofĤ S with positive rates Γ(±ω q ) = 2Re[ ∞ 0 dτ e ±iωqτ L(τ )] and, excluding constant terms,Ĥ s = − ∆ sσz /2 is the bath-induced energy shift of the qubit characterized by the correction to its bare frequency ω q , with Λ(±ω q ) = Im[ ∞ 0 dτ e ±iωqτ L(τ )]. Note that such a correction comprises both the vacuum and thermal contributions of the bath, traditionally referred to as Lamb and Stark shifts, respectively. For a non-driven system characterized byĤ d (t) = 0, equation (5) describes the thermalization process of the qubit with the heat bath, the superoperator of which commutes with the unitary part of the master equation. Consequently, the specific value of the shift ∆ s does not affect the steady-state quantities since the coherences in the eigenbasis ofĤ S vanish in the limit t → ∞. For a driven system in contrast, neglecting the shift may lead to incorrect predictions both dynamically and in the steady state as we show in section 3.
The assumption of a stationary thermal bath implies the drive field to contribute only to the unitary part of the master equation (Appendix A).
As discussed in the introduction, the limitations imposed by the above-mentioned approximations can be handled through different strategies. Here we use the framework of stochastic Liouville equations, where the non-perturbative treatment of the system-bath coupling relies on a stochastic unraveling of the reduced density operatorρ S (t). In this procedure, one resorts to the path integral description of quantum mechanics [101] to establish a numerically exact model of the open-quantum-system dynamics with the help of a classical stochastic process [72,102]. For the Ohmic spectral density in equation (3), and in the limit ω c → ∞, a single trajectory for the state of the system,ρ S (t), is given by the so-called stochastic Liouville equation with dissipation, SLED, as [72,76] d dtρ where ξ(t) is a real-valued classical random variable with null mean and autocorrelation with E[.] denoting the ensemble average over the noise trajectories (see Appendix B for details). Upon a suitable choice of a real-valued kernel, the quantity ξ(t) can be generated from a delta-correlated Gaussian noise, indeed providing for it a fully stochastic interpretation (see Appendix B.1). This term in equation (7) encodes the quantum fluctuations neglected when one treats the bath classically, as in the high-temperature limit through the Caldeira-Leggett master equation [40,72,103]. The actual reduced density operator of the system is obtained asρ S (t) = E[ρ S (t)] in the limit of infinite trajectories. Note that for a single generation of ξ(t), equation (7) is deterministic, so that parallelization and usual numerical methods for quantum evolution can be combined to speed up the SLED calculations. In order to simplify the comparisons between the SLED and Lindblad approaches in the next sections, we establish a connection between the transition rates in equation (5) and equation (7). This is achieved by writing Γ(ω q ) = γ[n(ω q ) + 1] and Γ(−ω q ) = γn(ω q ), with γ = 2ηω q being an effective qubit dissipation rate calculated in the limit ω c ω q . Consequently, η = γ/(2ω q ) in equation (7).
Below, we investigate the features promoted by the numerically exact SLED on the reduced state of the qubit for a nearly resonant drive fieldĤ d (t) and for different dissipation rates which are produced by the tunable environment of the qubit. We compare the exact predictions of such a method against the approximate LME. To this end, we define the fidelity of a density operator of an approximate solutionρ 1 against the numerically exact SLED solutionρ 2 for qubits as [104] The infidelity, 1 − F, describes the distance between the two states, and hence provides information on the amount of error introduced by the approximations used in the inexact approaches. In particular, we compare the fidelities of the Lindblad solutions obtained with and without the inclusion of the bath-induced energy shift defined in equation (6). In the following, the latter case is referred to as LME-nES.

Monochromatic periodic field
In this section, we compare the Lindblad and SLED approaches for the dissipative dynamics of a driven qubit. Namely, we solve the Lindblad equation (5) and compare the results with 250 GHz bath temperature T 48 mK those of SLED (7) using a sufficiently large number of noise realizations ξ(t) to obtain the reduced system density operator as an average over individual trajectories. Both LME and SLED are solved in the Liouville space [105] upon second-order Magnus expansions of the discretized time propagator associated to the corresponding Liouvillians. Particularly in the SLED, each noise realization ξ(t) is generated before the temporal evolution as in the recipe in equation (B.20). Here, we focus on the case in which the qubit is driven by a monochromatic and periodic transverse field at angular frequency ω d as where Ω d is referred to as the drive Rabi frequency. The corresponding full microscopic Hamiltonian is given by equation (1).
In figure 2a, we show the dynamics of the Bloch vector components ) in a frame rotating with the drive frequency ω d and for a qubit initially prepared in the excited state |1 . The drive frequency is set at the bare qubit frequency (ω d = ω q ) and the other parameters are chosen in compliance with the current state of art of cQED implementations, as shown in table 1. Note that the temperature of the bath is chosen such that βω q = 5, corresponding ton(ω q ) ≈ 6.8 × 10 −3 . For a typical superconducting-qubit frequency of ω q /(2π) = 5.0 GHz, such temperature corresponds to approximately 48 mK, lying within the achievable temperatures using dilution refrigerators in typical circuit QED setups. Moreover, the Rabi frequency is kept fixed at 1% of the qubit frequency, Ω d /(2π) = 50 MHz, and the environment cutoff frequency at ω c /(2π) = 250 GHz, which is extended to agree with the high-cutoff approximation ω c → ∞ in the SLED formalism. We simulate the qubit dynamics for a broad range of effective qubit dissipation rates γ, comprising values of 2π × 2.5 MHz (0.05% of ω q ) up to 2π × 500 MHz (10% of ω q ). Such tunability has been demonstrated in similar physical setups in recent protocols for engineered environments [78][79][80][81][82][83]. For simplicity of comparison between Lindblad and SLED in the present model, we assume that the intrinsic dephasing and decay rates of the qubit are low compared to Ω d and γ, such that they have a negligible effect on the qubit dynamics.
We highlight that the parameters associated to the superconducting resonator in table 1 are chosen as reference and do not enter in the numerical simulations of this work. They are only considered in the phenomenological description of the measurement setup implemented for both qubit readout and pump-probe spectroscopy as it will be detailed in section 4. The large detuning between the bare qubit and resonator frequencies has been extensively used in the dispersive readout of superconducting qubits [106].
In figure 2a, we find a very good agreement between the SLED and Lindblad methods if the bath-induced energy shift is not neglected and if we use a relatively small value of γ = 5ω q × 10 −3 . Thus the Born-Markov approximations in this weak-coupling regime seem valid. Such an agreement is manifested by the considerably high value of the corresponding temporally averaged fidelityF shown in figure 2b. However, neglecting the bath-induced energy shift alters the dynamics significantly, reducingF by approximately 10%. The reduction is even more pronounced in the case γ = ω q × 10 −2 , where the dissipation rate and the Rabi frequency are equal in magnitude. Therefore, such a deviation is maximum at the critical ratio c ≡ γ/Ω d ≈ 1, which can be analytically obtained from the asymptotic fidelities between the Lindblad solutions with and without the bath-induced energy shift, see Appendix A.1. Taking such an energy shift into account is naturally addressed beyond the LME [91]. However, we observe quantitatively that it can be relevant for the dynamics of a dissipative driven system even within the weak-coupling approximations. For on-demand dissipation, such a shift needs to be carefully considered in the experiments.
In addition to giving rise to faster stabilization time scales, the progressive increase of dissipation over the Rabi frequency (γ/Ω d > 1) attenuates the relevance of the bath-induced energy shift in the Lindblad dynamics, as shown in figure 2a for γ = ω q ×10 −1 , and in figure 2b for several different dissipation rates. We attribute this intriguing behavior of the dissipative dynamics of the driven qubit to the amount of coherent superposition between the eigenstates of the bare qubit promoted by the drive, which is increased for a resonant drive but inhibited in the strongly dissipative regime. Thus, for the chosen drive frequency ω d = ω q in figure 2, the inclusion of the bath-induced energy shift in the LME renders the drive nonresonant with the actual qubit frequency modified by the bath. This in turn leads to a non-vanishing x component of the Bloch vector in the rotating frame. With increasingly strong environmental coupling however, the amount of coherence promoted by the drive field decreases, even in the resonant case, and the decay towards the thermal steady state is favored. This state commutes withσ z , and hence is dynamically unaffected by the bath-induced energy shift.
Despite the negligible effect of the bath-induced energy shift on the steady state of the LME in the dissipation-dominated regime γ Ω d , the overall validity of the LME is compromised. Namely, figure 2b shows a progressive reduction of the fidelityF as a function of increasing γ Ω d . Manifestations of this breakdown are shown for the z component of the Bloch vector in Figure 2c, where we observe non-exponential short-time dynamics and a shift in the steady-state values for the SLED method, both being not contemplated by the Born-Markov approximations [76]. Note that owing to the low bath temperature, the effect of the bath-induced energy shift on the thermal populations in SLED solution is expected to be negligible here.

Steady state: bath-induced energy shift and failure of the Lindblad master equation
Let us further detail the importance of the bath-induced energy shift on the Lindblad description of the open dynamics by inspecting the steady-state properties of the driven qubit using equation (5). This serves as a guide for a measure of the inadequacy of the LME for asymptotic predictions, as we introduce below. The use of the SLED is aimed here at simulating the qubit dynamics in a typical cQED experiment.
For convenience, we denote the components of the steady state of the system in the rotating-frame for an arbitrary detuning ∆ q = ω q + ∆ s − ω d as where i = x, y, z, ∆ s is the bath-induced energy shift, andρ ss S =ρ S (t → ∞) is the asymptotic density operator of the system in Schrödinger's picture. In addition, in the case where the qubit is driven by a weak field (Ω d ω q ), the drive Hamiltonian can be written in the rotating-wave approximation such thatĤ Consequently, the asymptotic components σ ssf i (∆ q ) can be found analytically from the LME (5) for an arbitrary detuning ∆ q . Based on this result, we express how a nonresonant drive modifies the steady state of the qubit in comparison to the resonant-drive case ∆ q = 0 by defining the difference where the superscript 'L' highlights that such quantities are obtained through the Lindblad equation. As shown in Appendix A.1, we find where γ β = γ[2n(ω q ) + 1]. Therefore, the quantities ∆σ ss i indicate the LME predictions for the qubit sensitivity on the frequency change of a weak and monochromatic transverse drive.
In figure 3a, we show the dependence of ∆σ ss i on the qubit decay rate γ and on the angular frequency detuning of the drive ∆ q for a low-temperature environment (γ β ≈ γ).   (13), as a function of the dissipation rate γ and the detuning of the drive angular frequency from resonance ∆ q . Assuming that the drive is set to the angular frequency of the bare qubit, the bath-induced energy shift induces a finite detuning ∆ q as shown for SLED (blue markers) and for the LME (green markers). (b) The x component of the Bloch vector in the Schödinger picture as a function of time for SLED (solid line) and LME (dash-dotted line) simulations without drive (Ω d = 0). The qubit is initially prepared in the eigenstate ofσ x with eigenvalue 1. We fit an exponentially damped cosine function (dashed line) to the SLED solution. The obtained oscillation frequency corresponds to the qubit frequency modified by the bath which is used in (a) for the SLED data. (c) Difference ∆σ ss z , defined in equation (14), with the choice ∆ q = ∆ s for SLED (circles), numerically solved Lindblad (squares), and Lindblad with rotating-wave approximation (crosses) as a function of the qubit dissipation rate, γ. The colored region corresponding to ∆σ ss z < 0 indicates the failure of LME as expected from the definition in equation (14). All qubit parameters not given here are fixed according to table 1.
A non-resonant drive on the qubit (∆ q = 0) affects the different components of the Bloch vector in different ways. The difference in the x component changes its sign with that of the detuning and is clearly pronounced in the region of moderate dissipation, tending to vanish at |∆ q /Ω d | 1. On the other hand, the detuning barely affects the y component for γ/Ω d 1, and in this regime, ∆σ ss z saturates to a high value for a sufficiently large |∆ q /Ω d |. In general, the difference in all components of the Bloch vector decreases and becomes independent of the detuning with increasing γ/Ω d 1, which is another manifestation of dissipation dominating over the drive dynamics.
Provided that the drive frequency is set to ω d = ω q , we have ∆ q = ∆ s . For this choice of ∆ q , the markers in figure 3a show the bath induced energy-shift as function of the qubit decay rate. Whereas the perturbative approach of the dissipative dynamics allows one to obtain ∆ s directly from equation (6), we obtain it for the SLED method by fitting an exponentially damped cosine function to the early decay of the qubit coherence as illustrated in figure 3b. For the chosen parameters, the relation between ∆ s and γ is well approximated by a linear fit in both methods. As shown in [76], this dependence ceases to be linear for strong enough system-bath coupling strength.
Interestingly, equations (13) along with figure 3a show that ∆σ ss y and ∆σ ss z are symmetric with respect to ∆ q . Focusing our attention to the z component, we observe that ∆σ ss z ≥ 0 in equations (13) for any choice of parameters. Based on this result, we introduce the measure where σ ssf,L z (0) is the steady-state z component of the Bloch vector given by the LME at resonance and σ ssf z (∆ q ) is obtained by our method of choice or even experimentally. Note that if σ ssf z (∆ q ) is also obtained from the LME, equation (14) reduces to ∆σ ss z defined in equations (13), which is always positive. Thus, this measure can be used to identify regimes where the perturbative approach encoded in the LME is not sufficient to correctly predict the steady state of the weakly driven qubit. In an experiment, σ ssf,L On the other hand, σ ssf z (∆ q ) can be obtained through usual steady-state readout of the qubit driven out of resonance. The negativity of ∆σ ss z violates the lower bound imposed by equation (13), being a sufficient condition for the failure of the time-local LME. Figure 3c shows ∆σ ss z for selected values of γ/Ω d at ∆ q = ∆ s . The values of σ ssf z (∆ q ) are calculated from long-time solutions of the LME and SLED, which are intended to simulate the dynamics of the qubit in an experiment. The good agreement between equation (13) and the numerical results from the Lindblad equation highlights the validity of the RWA in the drive HamiltonianĤ d (t). In these cases, as expected from equation (13), ∆σ ss z is positive for all decay rates and achieves its maximum for γ ≈ Ω d . However, the detuned asymptotic z component given by SLED produces ∆σ ss z < 0 for dissipation rates γ > 0.05 × ω q , or in terms of the parameters of table 1, for γ > 2π ×250 MHz. This is a clear evidence of incompatibility with the used weak-coupling approximations.
As pointed out in Ref. [76] for a non-driven qubit, the shift of σ z (t) given by the SLED compared to that by the Lindblad equation cannot be fully attributed to a bathinduced energy shift since the correlations between the qubit and the bath created during the dynamics contributes as well. By turning on a very weak drive field, one is potentially able to study threshold conditions where ∆σ ss z = 0, indicating that such correlations may be significant. From a different perspective, our measure serves as a fine benchmark of the time-local Lindblad equation for a weakly driven qubit, thus shedding light on the validity limits of weak-coupling assumptions on the open dynamics. However, the threshold is still relaxed in the sense that possible deviations owing to the weak-coupling assumptions are not directly detected if ∆σ ss z > 0 even though they tend to increase with γ in our particular case as shown in figure 3c. This observation clearly illustrates the practical relevance of using non-perturbative approaches for accurate predictions.

Pump-probe spectroscopy
The second example of a driven dissipative system presented in this work consists of a qubit driven by a bichromatic field of the form which is a sum of the monochromatic drive field of equation (10), referred to as the primary drive, and a probe field with angular frequency ω p and an associated Rabi angular frequency Ω p . Below, we employ the SLED and compare it with the Lindblad formalism to study the signatures of the qubit fluorescence spectrum by means of a pump-probe approach [32]. Specifically, assuming Ω p Ω d so that the probe field acts as a weak perturbation to the driven qubit, information about the spectrum under the primary drive is obtained from the response of the system to the probe as the angular frequency ω p is swept.
Rather than monitoring the radiation spectrum of the qubit, we study the response of the system through the temporally averaged z component of the Bloch vector where σ z (t) = Tr[σ zρS (t)] as above, n p t p is the length of the integration interval, and t f is the final time chosen such that the initial transient dynamics has a negligible effect onσ z . If the probe and the drive are out of resonance, σ z (t) tends to oscillate in time with an amplitude h z and frequency |∆ p | ≈ |ω p −ω d |, so that the average in equation (16) is calculated over multiple integers n p of its period t p = 2π/|∆ p |. At resonance ∆ p = 0, the temporal dependence of σ z (t) is negligible due to the single oscillation frequency in equation (15) and the considerably small chosen values of Ω d and Ω p compared to the qubit angular frequency (see table 1). In a typical cQED experiment, one can relateσ z with the field transmitted through a readout resonator dispersively coupled to the qubit [106]. In the semiclassical approximation [107], a phenomenological inclusion of the resonator yields for the asymptotic field amplitude  transmitted from the readout resonator to its output port (see Appendix C) where Ω m is the amplitude of a weak measurement drive continuously applied on the input port of the resonator, κ is the resonator energy decay rate that is assumed to be dominated by leakage to the output port, and χ is the so-called dispersive shift associated to the qubitresonator coupling. Figure 4 showsσ z and A as functions of the probe frequency for various dissipation rates γ. Similar to section 3, the parameters are chosen according to table 1 unless otherwise stated. Here, the drive frequency ω d is adjusted to the resonance with the frequency of the qubit including any bath-induced energy shifts for each γ. The bath-induced frequency shift is calculated as in section 3.1, that is, through equation (6) for the Lindblad master equation and through a fit to the damped decay of the qubit coherence for the SLED.
We show in figure 4 that for very weak dissipation there is a good agreement between the two methods forσ z and A. As γ increases, the solutions given by the two approaches tend to separate, indicating that the steady state of the Lindblad master equation significantly deviates from the one given by SLED. Despite numerical fluctuations caused by the finite number of noise trajectories used for SLED, we observe that some resonance-like features tend to be preserved even for the dissipation rate of the order of the Rabi angular frequency of the drive. The main differences arise from the probe-frequency-independent shift of the response. Similar to section 3, this effect is caused by a shift of σ z (t) given by SLED as compared to that produced by the Lindblad master equation, becoming more pronounced as γ increases. In contrast to figure 2, in which the drive frequency is fixed at the bare qubit frequency in both methods (ω d = ω q ), the shift appears in figure 4 for drive frequencies ω d matching the qubit transition frequency shifted by the bath (ω d = ω q + ∆ s ). This suggests that such a phenomenon is not primarily caused by the renormalization of the qubit frequency, or any unitary effect, but rather it may be attributed to the appearance of asymptotic system-bath correlations as the strong coupling is approached.
A qualitative analysis may connect the presented results with the actual qubit fluorescence spectrum predicted by Mollow in Ref. [26]. Typically, the radiation spectrum is proportional to the Fourier transform of a two-time correlation function of the system evaluated at its steady state, e.g., R(ω) = ∞ −∞ dτ e −iωτ Ĉ † (τ )Ĉ(0) ss , withĈ = |0 1| being an example in the case of a qubit. For a weak environmental coupling, the calculation of such correlation functions is usually obtained through the quantum regression theorem [39,40], where one resorts to the Born-Markov approximations. Within this approach, for γ/Ω d 1, the fluorescence spectrum of a dissipative qubit that is driven by a resonant field of the form of equation (10) in the RWA presents three peaks centered at frequencies ω 0 = ω d and ω ± = ω d ± Ω d [26]. The sideband peaks have a Lorentzian shape that becomes broadened and flattened as the ratio γ/Ω d increases. These features are present in the top and middle panels of figure 4a, where we fit the data provided by the SLED with Lorenztian functions peaked roughly at ω − = 0.99 × ω d and ω + = 1.01 × ω d . In the bottom panels, however, the sideband peaks are absent due to the high qubit dissipation rate. The small oscillations in this case, which are noticeable in both the SLED and Lindblad data, may be attributed to the different number of periods n p used in the numerical integration of equation (16).
In contrast to the studies of Ref. [26], no central peak at ω 0 = ω d is observed in figure 4, since the probe field does not excite the qubit at resonance, according to the definition of f (t) in equation (15). This can also be checked by writing the driven qubit Hamiltonian in the frame rotating at the primary drive frequency. Since the pump-probe approach renders the long-time behavior of σ z (t) an indicator of the presence of the probe field, the latter is not perceived by the qubit when ∆ p = 0. Alternatively, this and the above-mentioned qualitative features of the qubit fluorescence spectrum can be checked in figure 5, where we show the amplitude h z of the probe-induced oscillations of σ z (t) as a function of the probe frequency ω p and a broad range of qubit dissipation rates γ. We clearly observe that the regions of high amplitude indicate the sideband peaks of the Mollow triplet, these being pronounced and narrow at small γ/Ω d . These peaks become flat and broad at high qubit dissipation rates, eventually disappearing at γ/Ω d 1. Similar damped oscillations as shown here have also been observed in different physical setups, for instance, in two coupled degenerate resonators with significantly different leakage rates [79]. In addition, the radiation spectrum of a qubit under a bichromatic field in the RWA has also been obtained through the quantum regression theorem and presents a rich variety of phenomena depending on the choice of parameters [108]. However, the features of the qubit spectrum under the drive field of equation (10) are preserved assuming that it is much stronger than the probe field, the case considered in this work.

Conclusions
We assessed the precision of weak-coupling assumptions of a driven qubit interacting with a bosonic environment through examples where the non-perturbative stochastic Liouvillevon Neumann equation, or SLED, is appropriate. We focused our attention on the case where the qubit interacts with weak and nearly resonant transverse fields along with a cold Ohmic bath with low dissipation rates compared with the bare system frequency. Such a scenario is typically addressed by the Lindblad master equation, or LME, and it is of practical relevance in state-of-art implementations of engineered environments in circuit quantum electrodynamics aimed, for example, at optimized initialization protocols of the system. Thus, our investigation complements the recent studies published in [76] and [77] on the benchmark of SLED over perturbative master equations.
We carried out a quantitative comparison of SLED with the LME and showed that the often overlooked bath-induced energy shift in the LME becomes less relevant for the dynamics with the strength of the dissipation increasing well beyond the drive Rabi frequency. However, new effects arising from the failure of the weak-coupling assumptions emerge in these regimes, being captured by the non-perturbative treatment of the drive and dissipation given by the SLED. In addition, we proposed a measure based on the sensitivity of the qubit population to the drive frequency. As a consequence, we identified regimes where the SLED yields for the steady state of the qubit dynamics distinctive and quantitatively measurable differences to the results of the Lindblad equation. Moreover, we have used the SLED and Lindblad approaches to study the signatures of the qubit fluorescence spectrum for different dissipation rates that may be produced by tunable environments in cQED.
In conclusion, our results may guide future experiments to probe driven open quantum systems and the validity of the weak-coupling approximations in describing their dynamics. This potentially allows for the exploration of undiscovered frontiers which are not well captured by the weak-coupling Markovian dynamics. In particular, our work may motivate further investigations on the validity of other perturbative master equations, such as the Floquet-Born-Markov equation, where the driving has been more accurately taken into account in the derivation of the master equation. However, despite improvements arising from deriving the dissipators in the dressed state basis, such equations are nevertheless perturbative. As a consequence, a precise proposal for experiments and corresponding parameters in scenarios of increasing dissipation as presented in this work calls for a model contemplating both the drive-dissipation interplay and high-order corrections to the systembath correlations, as given by the SLED. A more transparent comparison of the SLED with other non-perturbative approaches in the context of circuit quantum electrodynamics emerges as natural future line of research.

Appendix A. Lindblad master equation
Below, we review the derivation of the Lindblad master equation for the driven dissipative qubit described in section 2 of the main text. Similar derivations have been reported in the existing literature [39,40,99], and hence the discussion in this Appendix is given mainly for the sake of completeness of our notation.
We begin by employing the interaction picture with respect to the free Hamiltonian H S +Ĥ B , such that the exact temporal evolution of the total density operatorρ i (t) is given by the Liouville-von Neumann equation and the superscript i stands for the interaction picture. A recursive integration up to the second order and a trace over the bath degrees of freedom yield is the reduced density operator of the system at time instant t and we have assumed an initially factorized state as mentioned in the main text. Since the first moments of observables calculated in the stateρ i B (0) may be chosen to vanish, one naturally obtains Tr B [Ĥ i SB (t),ρ i (0)] = 0, which is hence not visible in equation (A.1). In the absence of the drive [f (t) = 0], only the third term on the right side of equation (A.1) remains.
The structure of equation (A.1) is rather complicated as it is an integro-differential exact equation. However, it can be simplified under a series of assumptions. First, one writes the total density operator asρ is present only when system-bath correlations are created during the dynamics. Naturally, Tr[ω i (t)] = 0 in order to preserve the normalization. The so-called Born approximation physically asserts that the influence of the system on the bath dynamics is small so that it essentially stays in the Gibbs state throughout the interaction. Consequently, one neglects the correlation termω i (t ) when consideringρ i (t ) in equation (A.1), and hence we obtain Note that the Born approximation implemented inside the the double commutators of equation (A.1) does not guarantee the conservation of the system entropy as in a unitary evolution. We note that in equation (A.2), the two last terms in equation (A.1) have been dropped. This is a consequence of the Born approximation, i.e.,ω i (t) → 0. Therefore, the Born approximation does not account for the drive-dissipation interplay promoted by such terms in the case of linear interaction of the system with a heat bath. The overlooking of such an interplay also allows one to combine the first two terms of equation (A.2) into the single commutator −i[Ĥ i d (t),ρ i S (t)]/ so that the drive contributes only to the unitary dynamics of the qubit.
The last term of equation (A.2) is here simplified by assuming that the system dynamics is memoryless and by coarse-graining in time (Markov approximation). The joint effect of these approximations allows one to neglect the dependence of the state of the system on its past history such thatρ i S (t ) →ρ i S (t) during the time integration, and to extend the integration limit up to infinity. This is usually justified as long as the bath autocorrelation function (4) arising from the double commutator decays faster than the relaxation time of the system. Upon the change of variable t → t − τ , these approximations lead to the master equation In general, the double commutator in equation Here, we show the analytical expressions for the components of the steady-state Bloch vector in the rotating frame according to the Lindblad equation, σ ssf,L i (∆ q ), which can be found following the procedure described in section 3.1. They read . (A.4) Using equations (A.4) in the definition for ∆σ ss i in equation (12), yields equations (13) of the main text.
One may also be interested in the fidelity between the steady states of the Lindblad master equation for the qubit driven at the bare frequency (∆ q = ∆ s ) and driven at the frequency shifted by the system-bath interactions (∆ q = 0). Such a fidelity can be obtained analytically by writing the steady states in terms of the Bloch vector components in equation (A.4) and using the definition in equation (9). Assuming γ β ≈ γ, we find .
As shown in figure 3a and from a direct evaluation of equation (6) where the integration measure with c being a normalization factor, goes along all possible paths q(t ) between q i = q(0) and q f = q(t). Here, the temporal dependence ofĤ(t ) is manifested by the explicit dependence on t in the classical Lagrangian L(q, t ) = H(q, t ) − V(q, t ), with H(q, t ) and V(q, t ) being the kinetic-like and potential energy of the system, respectively. Suppose that the considered system is multipartite and one is just interested in the reduced dynamics of a subpartition with position and momentum operatorsq andp, respectively, and the corresponding eigenstates |q and |p . Moreover, suppose that the temporal dependence ofĤ(t) arises strictly from the free Hamiltonian of this main subsystem, denoted here byĤ S (t). By assuming that the state of the main subsystem is initially factorized from the rest, the path integral formalism allows one to express the reduced density operator in position representation, ρ S (q f , q f ) = q f |ρ S (t)|q f , as where classical action S S [q, t] is associated withĤ S (t). All the dynamical effects of the secondary subsystems on the main one are encoded in the so-called influence functional F [q, q ], which equals unity in absence of interaction.
In this work, one assumes that the secondary subsystems form a thermal bosonic bath and its interaction with the main subsystem is linear through the position coordinates as described by the Caldeira-Leggett model [103]. If the main subsystem is a driven qubit as described in section 2,Ĥ S (t) =Ĥ S +Ĥ d (t), the Caldeira-Leggett model reduces to the Hamiltonian in equation (1), and the influence functional can be cast into the form is a phase functional with real and imaginary parts [72,109] where we defined new integration path variables u = q − q and v = q + q . Despite the exact expression for the influence functional F [u, v], its calculation is nontrivial since it is nonlocal in time. For the choice of the spectral density J(ω) in equation (3) and in the limit of high cutoff frequency (ω c → ∞), such temporal nonlocality can be only attributed to the finite temperature of the bath. In this case, one can rewrite the phase functional as the sum of temporally nonlocal and temporally local phases, i.e., In equation (B.6) one has defined as the white noise deducted from the real part of L(τ ).
Here, one can establish a numerically exact correspondence of F [u, v] with a temporally local averaged functional E{F ξ [u, v]} arising from a classical stochastic process. This process is described by the real-valued classical random variable ξ(t) with null mean and autocorrelation Placing equation (B.9) into equation (B.6) and using a Hubbard-Stratonovich transformation [110,111], one can write the influence functional F [u, v] as and consequently equation (B.3) reduces to ρ S (q f , q f ) = E[ρ S,ξ (q f , q f )]. Namely, the actual density operator ρ S (q f , q f ) can be regarded as the initial state ρ S (q i , q i ) evolving according to the stochastic influence functional F ξ [u, v] and averaged over a large number of noise trajectories. By returning to the operator representation and making the replacementŝ q →σ x ,p → ω qσy , the evolution corresponding to a single realization of ξ(t) is given by equation (7) of the main text. Such equation is deterministic and a single realization of ξ(t) can be generated from an arbitrary Gaussian random variable (see Appendix B.1). Interestingly, equation (7) has the form of the Caldeira-Leggett master equation [40] with −ξ(t)σ x added to the unitary part. Physically, such a term is responsible to account for the quantum fluctuations neglected in the classical treatment of the dissipative environment [72]. It has been shown [102,109] that the full stochastic unraveling of ρ S (q f , q f ) for an arbitrary spectral density J(ω) requires the inclusion of two complex-valued random variables, therefore making the numerical convergence slower than that of SLED. A thorough analysis involving such an extended method is out of the scope of this work.

Appendix B.1. Noise generation in SLED
Here, we describe the procedure for the generation of the stochastic noise ξ(t) which appears in equation (7). The initial point is to consider a Gaussian random variable r(t), the autocorrelation function of which is given by (B.11) Then one can define a real-valued convolution kernel G(t) in such a way that the noise ξ(t) is written as By using the definition (B.12) and the property (B.11), the autocorrelation function of the noise ξ(t) becomes Therefore, ξ(t) can be regarded as the inverse transform ofG(ω)r(ω), withG(ω) being obtained through equation (B.18) andr(ω) being generated from a Gaussian random variable r(t). In this work, we have used Python built-in functions for generating r(t) and for calculating the Fourier and inverse transforms.
the tunable resistor, the latter thus being the main source of dissipation for the qubit within the considered times scales. On the other hand, the dissipative dynamics of the resonator is caused by its finite quality factor so that it behaves as a lossy cavity where photons can leak out incoherently. These features can be represented by a master equation of the form d dtρ qr (t) = − i Ĥ qr (t),ρ qr (t) + L γ [ρ qr (t)] + L κ [ρ qr (t)] (C. 2) where the commutator describes the unitary dynamics determined by the Hamiltonian in equation (C.1), and L γ/κ [ρ qr (t)] describe the non-unitary dynamics promoted by the tunable resistor and the lossy resonator, respectively. Hence, L γ/κ [ρ qr (t)] contains only operators either in the qubit or the resonator subspace and already includes all possible bath-induced energy shifts and dissipative effects. Based on the arguments presented above, equation (C.2) allows one to separate the nonunitary effects produced by each local bath and describe the mutually induced frequency shifts in the qubit-resonator system by the effective HamiltonianĤ qr (t). While the average number of photons in the resonator contributes to the shift of the qubit frequency, the qubitdependent frequency of the resonator changes the behavior of the transmitted field providing an indirect measurement of the driven qubit spectrum as a function of the probe frequency ω p . In order to visualize such a phenomenon, we study the temporal evolution of the expectation value ofâ in a frame rotating at ω m , a(t) = Tr[âe iωmtâ †âρ qr (t)e −iωmtâ †â ]. First, we assume that L κ [ρ qr (t)] is phenomenologically described in the Lindblad form L κ [ρ qr (t)] = κ âρ qr (t)â † − 1 2 â †â ,ρ qr (t) , (C.3) where for shortness of notation we omitted the energy shift term, assuming that it is already incorporated to the definition of ω r . Using equation (C.2), we can write the dynamical equation for the expectation value of the annihilation operator aṡ where we have defined ∆ rm = ω r − ω m and a z (t) = Tr[σ zâ e iωmtâ †âρ qr (t)e −iωmtâ †â ]. In the semiclassical approximation [107], one can neglect the influence of qubit-resonator entanglement on the temporal evolution of a z (t), in such a way that a z (t) ≈ a(t)σ z (t), with σ z (t) = Tr[σ zρqr (t)]. Consequently, equation (C.4) can be rewritten aṡ a(t) = −i Ω m 2 − i∆ rm + iχσ z (t) + κ 2 a(t). (C.5) Here, we assume that the measurement field Ω m is turned on at a sufficiently long time after the initial transient dynamics of the dissipative driven qubit. Consequently, provided that the probe is much weaker than the drive, i.e., Ω p Ω d , we can write the solution to equation (C.5) as a(t) = a(0) − iΩ m e [i(∆rm+χσz)+κ/2]t − 1 κ + 2i (∆ rm + χσ z ) e −[i(∆rm+χσz)+κ/2]t , (C. 6) whereσ z is the temporal average of σ z (t). For the choice of t = n p t p as in section 4,σ z may be written as in equation (16). In our approach, we solve the dissipative dynamics of the driven qubit and feed the solution of a(t) with the values ofσ z . As explained in the main text, in this work σ z (t) is obtained either from the Lindblad master equation (5) or by the average solution of the SLED in equation (7). Regardless on the method of solution of the qubit dynamics, equation (C.6) clearly shows that a(t) contains information about the qubit population and, therefore, contains information about its spectrum. One can access the features of the spectrum, for instance, through the amplitude, phase, or Fourier transform of the transmitted field. Defining the field quadratures as I(t) = Re[a(t)] and Q(t) = Im[a(t)], the amplitude A(t) and phase φ(t) of the transmitted signal can be expressed as [106] A(t) = I 2 (t) + Q 2 (t), φ(t) = Arg[a(t)]. (C.7) Finally, by setting ∆ rm = 0, assuming κt/2 = κn p t p /2 1, and using the definition for A(t) in equation (C.7), one finds the transmitted field amplitude as defined in equation (17).