# Nonlinear response theory for lossy superconducting quantum circuits

V. Vadimov<EMAIL_ADDRESS>QCD Labs, QTF Centre of Excellence,
Department of Applied Physics, Aalto University, P.O. Box 15100, FI-00076
Aalto, Finland  M. Xu  Institute for Complex Quantum Systems and IQST, Ulm
University - Albert-Einstein-Allee 11, D-89069 Ulm, Germany  J. T. Stockburger
Institute for Complex Quantum Systems and IQST, Ulm University - Albert-
Einstein-Allee 11, D-89069 Ulm, Germany  J. Ankerhold  Institute for Complex
Quantum Systems and IQST, Ulm University - Albert-Einstein-Allee 11, D-89069
Ulm, Germany  M. Möttönen  QCD Labs, QTF Centre of Excellence, Department of
Applied Physics, Aalto University, P.O. Box 15100, FI-00076 Aalto, Finland
VTT Technical Research Centre of Finland Ltd, QTF Centre of Excellence, P.O.
Box 1000, FI-02044 VTT, Finland

###### Abstract

We introduce a numerically exact and yet computationally feasible nonlinear
response theory developed for lossy superconducting quantum circuits based on
a framework of quantum dissipation in a minimally extended state space.
Starting from the Feynman–Vernon path integral formalism for open quantum
systems with the system degrees of freedom being the nonlinear elements of the
circuit, we eliminate the temporally non-local influence functional of all
linear elements by introducing auxiliary harmonic modes with complex-valued
frequencies coupled to the non-linear degrees of freedom of the circuit. In
our work, we propose a concept of time-averaged observables, inspired by
experiment, and provide an explicit formula for producing their
quasiprobability distribution. Furthermore, we systematically derive a weak-
coupling approximation in the presence of a drive, and demonstrate the
applicability of our formalism through a study on the dispersive readout of a
superconducting qubit. The developed framework enables a comprehensive fully
quantum-mechanical treatment of nonlinear quantum circuits coupled to their
environment, without the limitations of typical approaches to weak
dissipation, high temperature, and weak drive. Furthermore, we discuss the
implications of our findings to the quantum measurement theory.

††preprint: APS/123-QED

##  I Introduction

Quantum systems inherently interact with their environments, typically
consisting of a macroscopic number of degrees of freedom [1, 2, 3]. These
interactions give rise to a diverse range of phenomena, including decoherence,
dissipation-induced phase transitions [4, 5], and emergence of classical
physics in quantum systems [6] among many others [3]. Understanding and
accurately describing the dynamics of open quantum systems is of fundamental
importance in various fields, ranging from quantum information science and
quantum computing [7, 8, 9, 10] to condensed matter physics [11] and quantum
optics [1, 12]. Recent technological advances [13, 14, 15, 16, 17, 18, 19, 20]
have opened up new possibilities for experimental studies of open quantum
systems, necessitating the development of accurate and computationally
efficient theoretical models.

Traditional Markovian approaches [1, 21, 22, 23] for treating open quantum
systems rely on assumptions of weak coupling and time-scale separation of the
system and the bath, or substantially high temperature, which may not always
hold in experimentally relevant setups leading to non-Markovian dynamics [24,
25, 26, 27, 28, 29, 30]. Furthermore, these approaches may have limited
applicability in the presence of structured reservoirs with gaps or other
singularities in the spectral density. Non-perturbative treatments often
involve discretization of the bath [31, 32] or employ the Feynman-Vernon path
integral formalism [33]. The latter includes techniques such as path integral
Monte Carlo (PIMC) [34, 35, 36, 37], unraveling the time-nonlocal influence
functional using auxiliary stochastic fields [38, 39], introducing a set of
auxiliary density operators governed by a system of time-local hierarchical
equations of motion (HEOM) [40, 41], and others [42, 43, 44, 45, 46].
Recently, a broad variety of numerically exact methods has been unified within
the framework of quantum dissipation in a minimally extended state space (QD-
MESS) [47, 48] which provides high accuracy with relatively moderate
computational resources.

In addition to the inevitable interaction with the natural environment, the
study of open quantum systems is also motivated by the need to carry out
measurements on real quantum systems [49]. In such cases, the measurement
apparatus itself acts as an environment, leading to controllable decoherence.
In this setup standard von Neumann measurements [50] can be performed only on
the macroscopic measurement device or its part [49], hence a single
measurement gives a very limited information about the system of interest.
This has two consequences: first, no real measurement is instantaneous since
the result usually comes from integration of continuous measurements, and
second, even continuous measurements performed with a specific device cannot
provide access to arbitrary observables of the system of interest.

In many experimental setups it is common to probe the electromagnetic field
[51] emitted from the measured system. To this end, the input-output (IO)
theory [52, 53, 54] for quantum optical systems and superconducting quantum
circuits was developed to calculate the properties of the emitted far-field
radiation and its integral characteristics. A major advantage of the IO theory
is the possibility to generalize it for quantum networks [55], where the
output from one part of the system may serve as the input for another part.
These theories are quite accurate in quantum optical systems but rely on the
Born–Markov approximation, which may break down, e.g., in the strong coupling
regime. In fact, non-Markovian effects may be more pronounced in microwave
superconducting circuits [29, 30] which are among the leading platforms for
quantum computing applications [51, 56, 57]. Effects of retardation and
correlations between the system and environment are critical for the design
and optimization of superconducting qubits, quantum gates, initialization, and
readout protocols, especially in circuits with tunable dissipative elements
[16, 58, 59, 60, 61] and distributed-element qubits [62]. Thus, a rigorous and
computationally efficient IO theory that goes beyond the limitations of Markov
approximations is essential for accurate characterizing and harnessing the
dynamics of superconducting circuits.

The first non-Markovian IO theory was presented in Ref. [63] based on
temporally non-local Heisenberg equations of motion for the degrees of freedom
of the system and the field operators in the case of a one-dimensional chiral
field propagating in one direction. However, solving these equations can be
technically involved for nonlinear systems. Later, this approach was
generalized for quantum networks [64] and squeezed input fields [65]. In other
studies [66, 67], the IO formalism was developed for atoms in a cavity which
is coupled to a Markovian bath.

Due to the approximation of the chiral field, it is infeasible to properly
account for thermal effects. This approximation introduces negative-frequency
modes, which do not have a thermal state at a positive temperature. It is a
reasonable approximation in quantum optics, since the characteristic frequency
scale is much higher than the temperature scale, especially if the back-
scattering can be neglected. This is not the case for the superconducting
quantum circuits, where thermal effects are important and reflections from the
qubits and resonators are usually strong.

The main goal of this paper is to develop a general formalism for dissipative
superconducting quantum circuits which goes beyond the Born–Markov
approximation, utilizing modern techniques for open quantum systems. We refer
to this formalism as nonlinear response (NLR) theory. In superconducting
circuits, the sources and measurement devices of the microwave field are
typically connected to other components by dispersion-free coaxial or coplanar
waveguides. These transmission lines act as dissipation channels, in addition
to dielectric losses [68, 69] or quasiparticles in the superconductors [70,
71, 72, 73, 74].

In our formalism, we utilize the QD-MESS approach and account for the non-
Markovian environment which is represented by finite number of damped bosonic
modes [75, 47] coupled to the nonlinear degrees of freedom of the circuit. By
incorporating these auxiliary modes, we can capture the full non-Markovian
dynamics induced by the system-environment interaction, including memory
effects and correlations. Our formalism provides a systematic framework for
studying the input-output characteristics of the circuit, taking into account
the interplay between the dissipation, the nonlinearity of the circuit, and
the external driving fields. The suggested approach is similar to the black-
box quantization technique [76, 77], however it is free from the limitation of
weak dissipation [76] and the need of bath discretization [77].

The structure of the paper is as follows: In Sec. II, we provide a brief
introduction to the Caldeira–Leggett model, which describes a broad class of
dissipative quantum systems, and establish the connections between classical
Langevin equations and the influence functional components in Keldysh space.
We focus the remainder of the paper on the study of superconducting quantum
circuits. Section III contains the central result of our paper, presenting a
detailed derivation of the NLR theory. Starting from the classical Kirchhoff’s
law for the circuit and the insights gained from Sec. II, we derive the action
on the Keldysh contour for the circuit. We explicitly incorporate classical
driving and introduce a quantum source or counting field to obtain the
generating functional of the output field. We then proceed to eliminate the
linear degrees of freedom of the circuit and unravel the resulting time-
nonlocal influence functional by introducing auxiliary bosonic modes. This
leads to a time-local QD-MESS equation that governs the dynamics of the open
quantum system. We define observables as time averages of the output field and
provide an explicit formula for their quasiprobability distribution through
the generating functional. In Sec. IV, we focus on the systematic derivation
of the weak-coupling regime in the presence of possibly strong driving fields,
which are non-trivial to accommodate exactly in master equations. The coupling
to off-resonant auxiliary modes can be eliminated through the use of a
dynamical Schrieffer–Wolf transformation [78, 79]. In Sec. V, we illustrate
our theory by applying it to the problem of dispersive readout of a
superconducting transmon qubit. Finally, we summarize our findings and provide
our conclusions in Sec. VI.

##  II Caldeira–Leggett model

Throughout the paper we use the following notation: Lowercase bold characters
denote vectors, uppercase bold characters denote matrices. A bold capital
character with subscripts denotes a certain block of the matrix denoted by the
corresponding character. Functions of time and frequency denoted by the same
character are related to each other by Fourier transformation:

$\displaystyle
X(\omega)=\int\limits_{-\infty}^{+\infty}X(t)\mathrm{e}^{\imath\omega
t}\;\mathrm{d}t,$ (1) $\displaystyle
X(t)=\frac{1}{2\pi}\int\limits_{-\infty}^{+\infty}X(\omega)\mathrm{e}^{-\imath\omega
t}\;\mathrm{d}t.$ (2)

The operators in the Hilbert space are emphasized by $\hat{\cdot}$ accent, the
operators in the extended Liouville space are denoted by $\check{\cdot}$
accents.

We begin our analysis by considering the paradigmatic Caldeira–Leggett model
[80] for dissipative systems. The classical Lagrangian of this model is given
by

$\mathcal{L}(\dot{x},\dot{q}_{k},x,q_{k},t)=\frac{m\dot{x}^{2}}{2}-V(x,t)\\\
+\sum\limits_{k}\left[\frac{m_{k}\dot{q}_{k}^{2}}{2}-\frac{m_{k}\omega_{k}^{2}}{2}\left(q_{k}-\frac{c_{k}}{m_{k}\omega_{k}^{2}}x\right)^{2}\right],$
(3)

which describes a particle of mass $m$ in a possibly time-dependent potential
$V(x,t)$ interacting with a plethora of environmental harmonic oscillators
each with mass $m_{k}$, bare angular frequency $\omega_{k}$, and coupling
constant $c_{k}$ which provides the coupling between the system and the $k$:th
oscillator. The Feynman–Vernon path integral approach [33] is the conventional
quantum-mechanical treatment of this system. The time evolution of the full
density matrix of the system $W(\bm{x}_{+},\bm{x}_{-},t)$ from time
$t^{\mathrm{i}}$ to $t^{\mathrm{f}}>t^{\mathrm{i}}$ can be found as

$W\left(\bm{x}_{+}^{\mathrm{f}},\bm{x}_{-}^{\mathrm{f}},t^{\mathrm{f}}\right)=\int
K\left(\bm{x}_{+}^{\mathrm{f}},\bm{x}_{-}^{\mathrm{f}},t^{\mathrm{f}};\bm{x}_{+}^{\mathrm{i}},\bm{x}_{-}^{\mathrm{i}},t^{\mathrm{i}}\right)\\\
W\left(\bm{x}_{+}^{\mathrm{i}},\bm{x}_{-}^{\mathrm{i}},t^{\mathrm{i}}\right)\;\mathrm{d}\bm{x}_{+}^{\mathrm{i}}\;\mathrm{d}\bm{x}_{-}^{\mathrm{i}},$
(4)

where the subscripts “$+$” and “$-$” stand for the first and the second
arguments of the density matrix, respectively,
$\bm{x}_{\pm}=\begin{bmatrix}x_{\pm}&q_{1,\pm}&q_{2,\pm}&\ldots\end{bmatrix}^{\mathsf{T}}$
are vectors of all coordinates of all the degrees of freedom, and the
propagator $K(\ldots)$ is given by the following path integral:

$K(\ldots)=\int\mathcal{D}[\bm{x}_{+},\bm{x}_{-}]\exp\left\\{\frac{\imath}{\hbar}\left(S[\bm{x}_{+}]-S[\bm{x}_{-}]\right)\right\\}.$
(5)

Here, the path integration is carried out over the fixed-boundary trajectories
$\bm{x}_{\pm}\left(t^{\mathrm{i}/\mathrm{f}}\right)=\bm{x}_{\pm}^{\mathrm{i}/\mathrm{f}}$
and the functional $S[\bm{x}]$ is a classical action

$S[\bm{x}]=\int\limits_{t^{\mathrm{i}}}^{t^{\mathrm{f}}}\mathcal{L}\left(\dot{\bm{x}},\bm{x},t\right)\;\mathrm{d}t.$
(6)

It is convenient to carry out a Keldysh rotation to so-called classical and
quantum trajectories [11]

$\bm{x}_{\pm}=\bm{x}_{\mathrm{c}}\pm\frac{\bm{x}_{\mathrm{q}}}{2}.$ (7)

The choice of the Keldysh formulation for the path integral is justified by
its simplicity in the coupling to the external fields: classical fields couple
to quantum degrees of freedom of the system, and quantum counting fields
couple to classical degrees of freedom. We utilize this substitution in the
path integral (5), consider the limit $t^{\mathrm{i}}\to-\infty$ and
$t^{\mathrm{f}}\to+\infty$, and assume that the full system was initially in a
thermal state at temperature $T$. The drive, if there is any, is assumed to be
off in the past $\dot{V}(x,t)\to 0$ at $t\to-\infty$. We calculate the trace
over all degrees of freedom and obtain unity

$1=\int\mathcal{D}[\bm{x}_{\mathrm{c}},\bm{x}_{\mathrm{q}}]\exp\left\\{\frac{\imath}{\hbar}S[\bm{x}_{\mathrm{c}},\bm{x}_{\mathrm{q}}]\right\\},$
(8)

where the Keldysh action is given by [11]

$S[\bm{x}_{\mathrm{c}},\bm{x}_{\mathrm{q}}]=\int\limits_{-\infty}^{+\infty}\left\\{\frac{m}{2}\begin{bmatrix}x_{\mathrm{c}}&x_{\mathrm{q}}\end{bmatrix}\bm{D}^{-1}\left(\imath\partial_{t}\right)\begin{bmatrix}x_{\mathrm{c}}\\\
x_{\mathrm{q}}\end{bmatrix}\right.\\\
\left.-V\left(x_{\mathrm{c}}+\frac{x_{\mathrm{q}}}{2},t\right)+V\left(x_{\mathrm{c}}-\frac{x_{\mathrm{q}}}{2},t\right)\right.\\\
\left.+\sum\limits_{k}\left[\frac{m_{k}}{2}\begin{bmatrix}q_{k,\mathrm{c}}&q_{k,\mathrm{q}}\end{bmatrix}\bm{D}^{-1}\left(\imath\partial_{t}\right)\begin{bmatrix}q_{k,\mathrm{c}}\\\
q_{k,\mathrm{q}}\end{bmatrix}\right.\right.\\\
\left.\left.-\omega_{k}^{2}\left(q_{k,\mathrm{c}}-\frac{c_{k}x_{\mathrm{c}}}{m_{k}\omega_{k}^{2}}\right)\left(q_{k,\mathrm{q}}-\frac{c_{k}x_{\mathrm{q}}}{m_{k}\omega_{k}^{2}}\right)\right]\right\\}.$
(9)

Here, we introduced an inverse Green’s function

$\bm{D}^{-1}(\omega)=\begin{bmatrix}0&(\omega-\imath 0^{+})^{2}\\\
(\omega+\imath 0^{+})^{2}&2\imath
0^{+}\omega\coth\left(\frac{\hbar\omega}{2k_{\mathrm{B}}T}\right)\end{bmatrix},$
(10)

where $\bm{D}^{-1}(\imath\partial_{t})$ stands for the function of the
differential operator. The symbol $0^{+}$ should be treated as an
infinitesimal positive frequency introduced to provide proper causality. The
quantum-classical and classical-quantum components are the retarded and
advanced inverse Green’s functions of a free classical particle with unit
mass, respectively. The quantum-quantum or Keldysh component is a time non-
local regularization which arises from the thermal initial condition in the
action.

The action is quadratic in the bath degrees of freedom, and hence the path
integral over them can be evaluated explicitly. We obtain an effective action
for the system coordinate as

$S_{\mathrm{s}}\left[x_{\mathrm{c}},x_{\mathrm{q}}\right]\\\
=-\imath\hbar\ln\left(\int\mathcal{D}\left[q_{k,\mathrm{c}},q_{k,\mathrm{q}}\right]\exp\left\\{\frac{\imath}{\hbar}S\left[\bm{x}_{\mathrm{c}},\bm{x}_{\mathrm{q}}\right]\right\\}\right).$
(11)

This effective action is given by

$S_{\mathrm{s}}\left[x_{\mathrm{c}},x_{\mathrm{q}}\right]=\int\limits_{-\infty}^{+\infty}\left[m\dot{x}_{\mathrm{c}}\dot{x}_{\mathrm{q}}\phantom{\frac{x_{\mathrm{q}}}{2}}\right.\\\
\left.-V\left(x_{c}+\frac{x_{\mathrm{q}}}{2},t\right)+V\left(x_{\mathrm{c}}-\frac{x_{\mathrm{q}}}{2},t\right)\right]\;\mathrm{d}t\\\
-\frac{1}{2}\int\limits_{-\infty}^{+\infty}\begin{bmatrix}x_{\mathrm{c}}(t)&x_{\mathrm{q}}(t)\end{bmatrix}\bm{\Sigma}(t-t^{\prime})\begin{bmatrix}x_{\mathrm{c}}(t^{\prime})\\\
x_{\mathrm{q}}(t^{\prime})\end{bmatrix}\;\mathrm{d}t\;\mathrm{d}t^{\prime},\;$
(12)

where the integration over bath degrees of freedom results in the emergence of
a finite time non-local term called influence functional [33]. We refer to the
kernel of this functional $\bm{\Sigma}(t-t^{\prime})$ as self-energy in
analogy with the self-energy from condensed matter physics. The self-energy
has a similar structure to that of the inverse Green’s function, namely,

$\bm{\Sigma}(\omega)=\begin{bmatrix}0&\Sigma^{\mathrm{A}}(\omega)\\\
\Sigma^{\mathrm{R}}(\omega)&\Sigma^{\mathrm{K}}(\omega)\end{bmatrix}$ (13)

and $\bm{\Sigma}(t-t^{\prime})$ is its Fourier image

$\bm{\Sigma}(t-t^{\prime})=\frac{1}{2\pi}\int\limits_{-\infty}^{+\infty}\bm{\Sigma}(\omega)\mathrm{e}^{-\imath\omega(t-t^{\prime})}\;\mathrm{d}\omega.$
(14)

The retarded component $\Sigma^{\mathrm{R}}(\omega)$ is a classical force-
coordinate response function of the bath, the advanced component, extended to
the whole complex frequency plane,
$\Sigma^{\mathrm{A}}(\omega)=\left[\Sigma^{\mathrm{R}}(\omega^{\ast})\right]^{\ast}$
is its conjugate, and the Keldysh component is given by the fluctuation
dissipation theorem as

$\Sigma^{\mathrm{K}}(\omega)=\frac{1}{2}\left[\Sigma^{\mathrm{R}}(\omega)-\Sigma^{\mathrm{A}}(\omega)\right]\coth\left(\frac{\hbar\omega}{2k_{\mathrm{B}}T}\right).$
(15)

By applying the Hubbard–Stratonovich transformation to the Keldysh part of the
influence functional (12) and employing a saddle-point approximation, we
obtain the quasiclassical Langevin equation for the system as [11, 81]

$-m\ddot{x}_{\mathrm{c}}-\int\limits_{-\infty}^{+\infty}\Sigma^{\mathrm{R}}(t-t^{\prime})x_{\mathrm{c}}(t^{\prime})\;\mathrm{d}t^{\prime}-\frac{\partial
V(x_{\mathrm{c}},t)}{\partial x_{\mathrm{c}}}=\xi(t),$ (16)

where $\xi(t)$ is a Gaussian-noise term with the correlation function

$\langle\xi(t)\xi(t^{\prime})\rangle=\frac{\imath}{2\pi}\int\limits_{-\infty}^{+\infty}\Sigma^{\mathrm{K}}(\omega)\mathrm{e}^{-\imath\omega(t-t^{\prime})}\;\mathrm{d}\omega.$
(17)

A quantum counterpart of this equation is a Heisenberg–Langevin equation [82,
83, 84] which is written for coordinate operators and has an operator-valued
noise term with Gaussian statistics. However, its direct application beyond
linear systems is complicated due to the lack of a closed system of equations
for correlation functions.

The main conclusion from this section is that in order to construct the
Feynman–Vernon influence functional, it is enough to know only the classical
response function of the bath and its temperature. If we do not have a
microscopic model for the bath, but know its classical response function, we
can introduce the corresponding terms phenomenologically into the Keldysh
action.

##  III Nonlinear response theory for quantum circuits

Although many parts of this theory are general and agnostic to the physical
realization, we motivate our studies by a particular class of open quantum
systems, namely, superconducting quantum circuits. We consider circuits formed
by lumped elements with linear elements such as capacitors, inductors, Ohmic
resistors, and nonlinear elements such as Josephson junctions. Due to their
inherent nonlinearity, the latter are the key components of every
superconducting qubit. The control and the measurement of the circuit is
provided by coupling to semi-infinite transmission lines, which inevitably
introduce additional dissipation channels in the system.

Such systems indeed fall into the Caldeira–Legett class of systems. However,
from the experimental point of view, not all of its observables can be probed
directly. The only experimentally accessible quantities are the output fields
in the transmission lines, which belong rather to the environment. As we
further discuss below, the whole division between the system and the
environment may be ambiguous. Therefore, there is a need of a theory which
connects the classical drive field with the experimentally observable output
field.

###  III.1 Classical equations of motion

Figure 1: Node $k$ with flux $\varphi_{k}$, connected to a transmission line
and other nodes via capacitors, resistors, inductors, and Josephson junctions.
The output voltage in the transmission line is determined by the input voltage
and the voltage given by $\dot{\varphi}_{k}$.

Using conclusions from the previous section, we begin with writing classical
equations of motion for the electric circuit [85]. We define the flux
$\varphi_{k}$ for each node $k$ of the circuit and express Kirchhoff’s law for
it as

$-{\sum\limits_{k^{\prime}}C_{kk^{\prime}}\left(\ddot{\varphi}_{k}-\ddot{\varphi}_{k^{\prime}}\right)}-{\sum\limits_{k^{\prime}}\frac{\dot{\varphi}_{k}-\dot{\varphi}_{k^{\prime}}}{R_{kk^{\prime}}}}\\\
-{\sum\limits_{k^{\prime}}\frac{\varphi_{k}-\varphi_{k^{\prime}}}{L_{kk^{\prime}}}}-{\sum\limits_{k^{\prime}}i_{\mathrm{J},kk^{\prime}}}-{\sum\limits_{j}\frac{\dot{\varphi}_{k}-2v_{\mathrm{i},(kj)}}{Z_{(kj)}}}=0,$
(18)

where the first four sums describe the currents coming into the node $k$ from
capacitors $C_{kk^{\prime}}$, resistors $R_{kk^{\prime}}$, inductors
$L_{kk^{\prime}}$, and Josephson junctions which connect it to the other nodes
$k^{\prime}$, whereas the last sum corresponds to the current coming from the
transmission lines with characteristic impedances $Z_{(kj)}$ coupled to the
node $k$ as illustrated in Fig. 1. Here, index $j$ in a pair $(kj)$ enumerates
all the transmission lines connected to the node $k$. The Josephson current is
given by

$i_{\mathrm{J},kk^{\prime}}=\frac{2\pi
E_{\mathrm{J},kk^{\prime}}}{\Phi_{0}}\sin\left[\frac{2\pi}{\Phi_{0}}\left(\varphi_{k}-\varphi_{k^{\prime}}-\Phi_{kk^{\prime}}\right)\right],$
(19)

where $E_{\mathrm{J},kk^{\prime}}$ is the Josephson energy of the junction
between the nodes $k$ and $k^{\prime}$, $\Phi_{kk^{\prime}}$ is the flux bias
of the junction, and $\Phi_{0}$ is the superconducting flux quantum. If a pair
of nodes $k$ and $k^{\prime}$ is not connected by a capacitor then we put
$C_{kk^{\prime}}=0$. For the same cases related to resistors, inductors, and
Josephson junctions we put $R_{kk^{\prime}}\to\infty$,
$L_{kk^{\prime}}\to\infty$, and $E_{\mathrm{J},kk^{\prime}}=0$. If there is a
ground node in the circuit, we put the corresponding flux
$\varphi_{\mathrm{ground}}=0$. The output $v_{\mathrm{o},(kj)}$ voltage in the
transmission line $j$ connected to the node $k$ can be expressed through the
input field and the voltage at the node as

$v_{\mathrm{o},(kj)}=\dot{\varphi}_{k}-v_{\mathrm{i},(kj)}$ (20)

since the voltage at the node
$v_{\mathrm{o},(kj)}+v_{\mathrm{i},(kj)}=\dot{\varphi}_{k}$ must be single-
valued for all transmission lines $j$ connected to the node $k$.

We rewrite Eq. (18) in the matrix form, introducing notations $\bm{\varphi}$
as vector of the fluxes at all the nodes, $\bm{v}_{\mathrm{i/o}}$ as vectors
of input and output voltages in the transmission lines, and
$\bm{i}_{\mathrm{J}}$ as vector of currents through all the Josephson
junctions:

$\bm{K}^{\mathrm{R}}(\imath\partial_{t})\begin{bmatrix}\bm{\varphi}\\\
\bm{v}_{\mathrm{i}}\end{bmatrix}=\begin{bmatrix}\bm{\nabla}^{\mathsf{T}}\bm{i}_{\mathrm{J}}\\\
\bm{v}_{\mathrm{o}}\end{bmatrix},$ (21)

where the frequency dependent matrix $\bm{K}^{\mathrm{R}}(\omega)$ describes
the causal response of the linear elements, and matrix
$\bm{\nabla}^{\mathsf{T}}$ is a discrete divergence matrix for the Josephson
currents. The explicit expression for the elements of the matrix
$\bm{K}^{\mathrm{R}}(\omega)$ can be obtained by applying a Fourier
transformation to Kirchhoff’s law and the input-output relations. The matrix
$\bm{K}^{\mathrm{R}}(\omega)$ has a block structure, where blocks correspond
to the linear coupling between dynamical degrees of freedom and input and
output fields:

$\bm{K}^{\mathrm{R}}(\omega)=\begin{bmatrix}\bm{K}^{\mathrm{R}}_{\varphi\varphi}\left(\omega\right)&\bm{K}^{\mathrm{R}}_{\varphi\mathrm{i}}\left(\omega\right)\\\\[5.0pt]
\bm{K}^{\mathrm{R}}_{\mathrm{o}\varphi}\left(\omega\right)&\bm{K}^{\mathrm{R}}_{\mathrm{o}\mathrm{i}}\left(\omega\right)\\\
\end{bmatrix}.$ (22)

The first block $\bm{K}^{\mathrm{R}}_{\varphi\varphi}\left(\omega\right)$
describes the linear coupling between the nodes:

$K^{\mathrm{R}}_{\varphi\varphi;kk^{\prime}}\left(\omega\right)=\sum\limits_{j}\frac{\imath\omega}{Z_{(kj)}}\delta_{kk^{\prime}}\\\
-(-1)^{\delta_{kk^{\prime}}}\left(\omega^{2}C_{kk^{\prime}}+\frac{\imath\omega}{R_{kk^{\prime}}}-\frac{1}{L_{kk^{\prime}}}\right).$
(23)

The block $\bm{K}^{\mathrm{R}}_{\varphi\mathrm{i}}(\omega)$ describes the
coupling to the driving field:

$K^{\mathrm{R}}_{\varphi\varphi;k(k^{\prime}j)}=\frac{2}{Z_{kj}}\delta_{kk^{\prime}}.$
(24)

The latter blocks provide the expression of the output field through the node
fluxes and input voltages:

$\begin{gathered}K^{\mathrm{R}}_{\mathrm{o}\varphi;(kj)k^{\prime}}(\omega)=-\imath\omega\delta_{kk^{\prime}},\\\
K^{\mathrm{R}}_{\mathrm{o}\mathrm{i};(kj)(k^{\prime}j^{\prime})}(\omega)=-\delta_{kk^{\prime}}\delta_{jj^{\prime}}.\end{gathered}$
(25)

Note that the matrix $\bm{K}^{\mathrm{R}}(\omega)$ is nothing but the retarded
component of the self-energy coming from the dissipation channels, combined
with the contributions from the linear reactive elements.

###  III.2 Keldysh action

Let us write down the Keldysh action for the considered circuit, assuming that
all the dissipation channels are described by a uniform temperature $T$. The
full information on the dynamics of the circuit and its observables is stored
in the corresponding generating functional [11]. This can be expressed through
a path integral

$\mathcal{Z}[\bm{v}_{\mathrm{i}},\bm{\eta}_{\mathrm{o}}]=\int\mathcal{D}[\bm{\varphi}_{\mathrm{c}},\bm{\varphi}_{\mathrm{q}}]\exp\left\\{\frac{\imath}{\hbar}S_{\mathrm{circ}}[\bm{\phi}]\right\\},$
(26)

where
$\bm{\phi}=\begin{bmatrix}\bm{\varphi}_{\mathrm{c}}^{\mathsf{T}}&\bm{v}_{\mathrm{i}}^{\mathsf{T}}&\bm{\varphi}_{\mathrm{q}}^{\mathsf{T}}&\bm{\eta}_{\mathrm{o}}^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}$,
$\bm{\varphi}_{\mathrm{c}}$ and $\bm{\varphi}_{\mathrm{q}}$ stand for
classical and quantum trajectories of the dynamical degrees of freedom
$\bm{\varphi}$, and $\bm{\eta}_{\mathrm{o}}$ is an auxiliary quantum field
which is an argument of the generating functional. The latter plays the role
of a counting field used to generate the statistics of observables [2, 11] The
action of the circuit $S_{\mathrm{circ}}[\bm{\phi}]$ is given by

$S_{\mathrm{circ}}[\bm{\phi}]=S_{\mathrm{J}}[\bm{\varphi}_{\mathrm{c}},\bm{\varphi}_{\mathrm{q}}]+\frac{1}{4\pi}\int\limits_{-\infty}^{+\infty}\bm{\phi}^{\dagger}(\omega)\bm{G}^{-1}(\omega)\bm{\phi}(\omega)\;\mathrm{d}\omega,$
(27)

where

$S_{\mathrm{J}}[\bm{\varphi}_{\mathrm{c}},\bm{\varphi}_{\mathrm{q}}]=-\int\limits_{-\infty}^{+\infty}V_{\mathrm{J}}[\bm{\varphi}_{\mathrm{c}}(t),\bm{\varphi}_{\mathrm{q}}(t)]\;\mathrm{d}t$
(28)

and

$V_{\mathrm{J}}(\bm{\varphi}_{\mathrm{c}},\bm{\varphi}_{\mathrm{q}})\\\
=2\sum\limits_{k<k^{\prime}}E_{\mathrm{J},kk^{\prime}}\sin\left[\frac{2\pi\left(\psi_{\mathrm{c},kk^{\prime}}-\Phi_{kk^{\prime}}\right)}{\Phi_{0}}\right]\sin\left(\frac{\pi\psi_{\mathrm{q},kk^{\prime}}}{\Phi_{0}}\right).$
(29)

is the action for Josephson junctions and
$\psi_{\mathrm{c}/\mathrm{q},kk^{\prime}}=\varphi_{\mathrm{c}/\mathrm{q},k}-\varphi_{\mathrm{c}/\mathrm{q},k^{\prime}}$
are the classical and quantum variables for the flux difference across the
junction between the nodes $k$ and $k^{\prime}$. The inverse bare Green’s
function assumes the form

$\bm{G}^{-1}(\omega)=\begin{bmatrix}\bm{0}&\bm{K}^{\mathrm{A}}(\omega)\\\
\bm{K}^{\mathrm{R}}(\omega)&\bm{K}^{\mathrm{K}}(\omega)\end{bmatrix}.$ (30)

The retarded component is the same matrix $\bm{K}^{\mathrm{R}}(\omega)$
introduced in the previous subsection, the advanced component is its hermitian
conjugate

$\bm{K}^{\mathrm{A}}(\omega)=\left[\bm{K}^{\mathrm{R}}(\omega^{\ast})\right]^{\dagger},$
(31)

and the Keldysh component is determined by the fluctuation dissipation theorem

$\bm{K}^{\mathrm{K}}(\omega)=\frac{1}{2}\coth\left(\frac{\hbar\omega}{2k_{\mathrm{B}}T}\right)\begin{bmatrix}\bm{K}_{\varphi\varphi}^{\mathrm{R}}(\omega)-\bm{K}_{\varphi\varphi}^{\mathrm{A}}(\omega)&0\\\
0&0\end{bmatrix}.$ (32)

In addition to the classical driving field $\bm{v}_{\mathrm{i}}$, we have
introduced an auxiliary quantum source $\bm{\eta}_{\mathrm{o}}$ which couples
to the output field in the transmission lines. The saddle-point approximation
with respect to the quantum degrees of freedom and the quantum source
$\bm{\eta}_{\mathrm{o}}$ yields the classical equations of motion (21).

It is convenient to use vectors $\bm{\psi}_{\mathrm{c}/\mathrm{q}}$ of all
classical and quantum flux differences across all junctions as dynamical
variables. By doing so, we move from a node description of the electric
circuit to an element description. For simplicity, we assume that there are no
Josephson junction loops since otherwise we need to pick a minimal set of
linearly independent variables which fully parametrize the Josephson energy of
the circuit. They should be complemented by dynamical variables
$\bm{\psi}^{\prime}_{\mathrm{c}/\mathrm{q}}$ such that a linear transformation
of chosen real-valued matrices $\bm{Q}_{\varphi\psi}$ and
$\bm{Q}_{\varphi\psi^{\prime}}$

$\bm{\varphi}_{\mathrm{c}/\mathrm{q}}=\begin{bmatrix}\bm{Q}_{\varphi\psi}&\bm{Q}_{\varphi\psi^{\prime}}\end{bmatrix}\begin{bmatrix}\bm{\psi}_{\mathrm{c}/\mathrm{q}}\\\
\bm{\psi}^{\prime}_{\mathrm{c}/\mathrm{q}}\end{bmatrix}$ (33)

is reversible. Note that the action becomes quadratic with respect to the
$\bm{\psi}^{\prime}$ degrees of freedom, and hence they can be integrated out.
By doing this, we shift the boundary between the system and environment and
consider the linear subcircuit as part of the environment. The subcircuit of
remaining nonlinear elements we refer to as core system. The effective action
arising from the Gaussian integration over these degrees of freedom is given
by

$S_{\mathrm{core}}[\bm{\phi}_{\mathrm{core}}]=S_{\mathrm{J}}[\bm{\psi}_{\mathrm{c}},\bm{\psi}_{\mathrm{q}}]\\\
-\frac{1}{4\pi}\int\limits_{-\infty}^{+\infty}\bm{\phi}_{\mathrm{core}}^{\dagger}(\omega)\bm{\Sigma}_{\mathrm{core}}(\omega)\bm{\phi}_{\mathrm{core}}(\omega),$
(34)

where
$\bm{\phi}_{\mathrm{core}}=\begin{bmatrix}\bm{\psi}_{\mathrm{c}}^{\mathsf{T}}&\bm{v}_{\mathrm{i}}^{\mathsf{T}}&\bm{\psi}_{\mathrm{q}}^{\mathsf{T}}&\bm{\eta}_{\mathrm{o}}^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}$
and the self-energy becomes

$\bm{\Sigma}_{\mathrm{core}}(\omega)=\begin{bmatrix}\bm{0}&\bm{\Sigma}_{\mathrm{core}}^{\mathrm{A}}(\omega)\\\\[5.0pt]
\bm{\Sigma}_{\mathrm{core}}^{\mathrm{R}}(\omega)&\bm{\Sigma}_{\mathrm{core}}^{\mathrm{K}}(\omega)\end{bmatrix}\\\
=\tilde{\bm{G}}^{-1}_{\tilde{\bm{\phi}}\bm{\psi}^{\prime}}(\omega)\tilde{\bm{G}}_{\bm{\psi}^{\prime}\bm{\psi}^{\prime}}(\omega)\tilde{\bm{G}}^{-1}_{\bm{\psi}^{\prime}\tilde{\bm{\phi}}}(\omega)-\tilde{\bm{G}}^{-1}_{\tilde{\bm{\phi}}\tilde{\bm{\phi}}}(\omega),$
(35)

with
$\tilde{\bm{G}}^{-1}(\omega)={\bm{Q}}^{\mathsf{T}}\bm{G}^{-1}(\omega){\bm{Q}}$
and

${\bm{Q}}=\begin{bmatrix}\bm{Q}_{\varphi\psi}&\bm{Q}_{\varphi\psi^{\prime}}&0&0&0&0\\\
0&0&1&0&0&0\\\ 0&0&0&\bm{Q}_{\varphi\psi}&\bm{Q}_{\varphi\psi^{\prime}}&0\\\
0&0&0&0&0&1\end{bmatrix}.$ (36)

###  III.3 Analytical properties of the self-energy

The self-energy is a meromorphic matrix function of complex-valued frequency.
We can classify its poles into two categories. The first category of poles
corresponds to the frequencies $\omega_{k}^{\mathrm{d}}$ where the matrix
$\tilde{\bm{G}}_{\bm{\psi}^{\prime}\bm{\psi}^{\prime}}(\omega)$ is singular.
They have a purely classical origin and we refer to them as to dynamical
poles. Since these poles come from a matrix inverse, the corresponding
residues have unit rank. The other category of poles arise from the
Bose–Einstein distribution function $\coth[\hbar\omega/(2k_{\mathrm{B}}T)]$
and corresponds to bosonic Matsubara frequencies
$\imath\omega_{n}^{\mathrm{m}}=2\pi\imath nk_{\mathrm{B}}T/\hbar$, $n=\pm
1,\pm 2,\ldots$. Since all dynamical variables and source fields are real-
valued in the time domain, the self-energy has the symmetry property
$\bm{\Sigma}_{\mathrm{core}}(t)=\bm{\Sigma}^{\mathsf{T}}_{\mathrm{core}}(-t)$.
Hence the action $S_{\mathrm{core}}[{\bm{\phi}}_{\mathrm{core}}]$ can be
expressed through a modified self-energy with a retarded causality. Due to
this symmetry property

$\int\limits_{-\infty}^{+\infty}\bm{y}^{\dagger}(\omega)\bm{\Sigma}^{\mathrm{R}}_{\mathrm{core}}(\omega)\bm{x}(\omega)\;\mathrm{d}\omega\\\
=\int\limits_{-\infty}^{+\infty}\bm{x}^{\dagger}(\omega)\bm{\Sigma}^{\mathrm{A}}_{\mathrm{core}}(\omega)\bm{y}(\omega)\;\mathrm{d}\omega,$
(37)

where $\bm{x}(\omega)$ and $\bm{y}(\omega)$ are Fourier images of arbitrary
real-valued vector functions, the advanced component of self-energy can be
replaced with the retarded one. We apply a Hilbert transformation to the
Keldysh component

$\tilde{\bm{\Sigma}}_{\mathrm{core}}^{\mathrm{K}}(\omega)=\frac{\imath}{2\pi}\int\limits_{-\infty}^{+\infty}\frac{\bm{\Sigma}_{\mathrm{core}}^{\mathrm{K}}(\omega-\omega^{\prime})}{\omega^{\prime}+\imath
0^{+}}\;\mathrm{d}\omega^{\prime}.$ (38)

The resulted $\tilde{\bm{\Sigma}}_{\mathrm{core}}^{\mathrm{K}}(\omega)$ has a
retarded causality. It means that all of the poles of
$\tilde{\bm{\Sigma}}_{\mathrm{core}}^{\mathrm{K}}(\omega)$ are located in the
lower complex half-plane and coincide with the poles of
$\bm{\Sigma}_{\mathrm{core}}^{\mathrm{K}}(\omega)$. Then we construct a
rectangular self-energy matrix

$\tilde{\bm{\Sigma}}_{\mathrm{core}}(\omega)=\begin{bmatrix}\bm{\Sigma}^{\mathrm{R}}_{\mathrm{core}}(\omega)&\tilde{\bm{\Sigma}}_{\mathrm{core}}^{\mathrm{K}}(\omega)\end{bmatrix}$
(39)

which satisfies

$\frac{1}{2}\int\limits_{-\infty}^{+\infty}\begin{bmatrix}\bm{x}^{\dagger}(\omega)&\bm{y}^{\dagger}(\omega)\end{bmatrix}\bm{\Sigma}_{\mathrm{core}}(\omega)\begin{bmatrix}\bm{x}(\omega)\\\
\bm{y}(\omega)\end{bmatrix}\;\mathrm{d}\omega\\\
=\int\limits_{-\infty}^{+\infty}\bm{y}^{\dagger}(\omega)\tilde{\bm{\Sigma}}_{\mathrm{core}}(\omega)\begin{bmatrix}\bm{x}(\omega)\\\
\bm{y}(\omega)\end{bmatrix}\;\mathrm{d}\omega.$ (40)

A meromorphic function allows the Mittag–Leffler expansion which has the
following form for the self-energy:

$\tilde{\bm{\Sigma}}_{\mathrm{core}}(\omega)=\bm{P}(\omega)+\sum\limits_{k}\frac{\bm{R}_{k}^{\mathrm{d}}}{\omega-\omega_{k}^{\mathrm{d}}}+\sum\limits_{n}\frac{\begin{bmatrix}0&\bm{R}_{n}^{\mathrm{m}}\end{bmatrix}}{\omega+\imath\omega_{n}^{\mathrm{m}}},$
(41)

where $\bm{P}(\omega)=\bm{P}^{(0)}+\bm{P}^{(1)}\omega+\bm{P}^{(2)}\omega^{2}$
is a polynomial matrix, $\bm{R}_{k}^{\mathrm{d}}$ is a unit rank matrix
residue at the dynamical pole $k$, and $\bm{R}_{n}^{\mathrm{m}}$ is a matrix
residue at the Matsubara pole $n$ which contributes only to the Keldysh
component of the self-energy.

In connection to the above expansion, we note that a sharp cut-off of the
Matsubara series gives a quite poor approximation for the cotangent function
present in the Bose–Einstein distribution function. In practice it is more
efficient to use a rational approximation in the form

$\coth\left(\frac{\hbar\omega}{2k_{\mathrm{B}}T}\right)\approx\omega\tau+\frac{2k_{\mathrm{B}}T}{\hbar\omega}+\sum\limits_{n=1}^{N}\frac{2\omega
r_{n}^{\mathrm{gm}}}{\omega^{2}+(\omega_{n}^{\mathrm{gm}})^{2}},$ (42)

where $\omega_{n}^{\mathrm{gm}}$ are the poles of rational approximations,
$r_{n}^{\mathrm{gm}}$ are corresponding residues, and $\tau$ is a constant
parameter. Such an approximation can be obtained, e.g., with the symmetrized
version of AAA algorithm [86, 75, 47] which turns out to be very accurate in a
finite range of frequencies even with very few poles taken into account. We
refer to the frequencies $\omega_{n}^{\mathrm{gm}}$ as generalized Matsubara
frequencies.

Finally, the action for the circuit reads as

$S_{\mathrm{core}}[\bm{\phi}_{\mathrm{core}}]=S_{\mathrm{J}}[\bm{\psi}_{\mathrm{c}},\bm{\psi}_{\mathrm{q}}]\\\
+\int\limits_{-\infty}^{+\infty}\sum\limits_{jj^{\prime}=0}^{1}\frac{\mathrm{d}^{j}\bm{\phi}_{\mathrm{core,q}}^{\mathsf{T}}(t)}{\mathrm{d}t^{j}}\bm{A}^{(jj^{\prime})}\frac{\mathrm{d}^{j^{\prime}}\bm{\phi}_{\mathrm{core}}(t)}{\mathrm{d}t^{j^{\prime}}}\;\mathrm{d}t\\\
-\frac{1}{2\pi}\int\limits_{-\infty}^{+\infty}\sum\limits_{k}\bm{\phi}_{\mathrm{core,q}}^{\dagger}(\omega)\frac{\bm{R}_{k}^{\mathrm{d}}}{\omega-\omega_{k}^{\mathrm{d}}}\bm{\phi}_{\mathrm{core}}\;\mathrm{d}\omega\\\
-\frac{1}{2\pi}\int\limits_{-\infty}^{+\infty}\sum\limits_{n}\bm{\phi}_{\mathrm{core,q}}^{\dagger}(\omega)\frac{\bm{R}_{n}^{\mathrm{gm}}}{\omega+i\omega_{n}^{\mathrm{gm}}}\bm{\phi}_{\mathrm{core,q}}(\omega)\;\mathrm{d}\omega,$
(43)

where the influence functional in Eq. (34) has been replaced by the three
terms following
$S_{\mathrm{J}}[\bm{\psi}_{\mathrm{c}},\bm{\psi}_{\mathrm{q}}]$. Here

$\begin{gathered}\bm{A}^{(00)}=-\bm{P}^{(0)},\\\
\bm{A}^{(01)}-\bm{A}^{(10)}=-\imath\bm{P}^{(1)},\\\
\bm{A}^{(11)}=-\bm{P}^{(2)},\end{gathered}$ (44)

sum $\bm{A}^{(01)}+\bm{A}^{(10)}$ can be arbitrary, and
$\bm{R}^{\mathrm{gm}}_{n}$ are the matrix residues at generalized Matsubara
poles. We have denoted the quantum components of the vector
$\bm{\phi}_{\mathrm{core}}$ as
$\bm{\phi}_{\mathrm{core,q}}=\begin{bmatrix}\bm{\psi}_{\mathrm{q}}^{\mathsf{T}}&\bm{\eta}_{\mathrm{o}}^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}$.

###  III.4 Auxiliary modes

We treat each non-local term in the action separately introducing auxiliary
bosonic degrees of freedom via the Hubbard–Stratonovich transformation [87,
88, 47] as

$\exp\left[-\frac{\imath}{2\pi\hbar}\int\limits_{-\infty}^{+\infty}\bm{y}^{\dagger}(\omega)\frac{\bm{R}}{\omega-\omega_{\ast}}\bm{x}(\omega)\;\mathrm{d}\omega\right]\\\
=\int\mathcal{D}[\bm{a},\bm{a}^{\dagger}]\exp\left\\{\frac{\imath}{\hbar}\int\limits_{-\infty}^{+\infty}\left[\imath\hbar\bm{a}^{\dagger}(t)\dot{\bm{a}}(t)-\hbar\omega_{\ast}|\bm{a}(t)|^{2}\right.\right.\\\
\left.\phantom{\int}\left.+\bm{a}^{\dagger}(t)\bm{V}^{\mathsf{T}}\bm{x}(t)+\bm{y}^{\mathsf{T}}(t)\bm{U}\bm{a}(t)\right]\;\mathrm{d}t\right\\},$
(45)

where $\bm{x}(t)$ and $\bm{y}(t)$ are arbitrary real valued-vector functions,
$\bm{R}$ is an arbitrary matrix which is factorized as
$\bm{R}=\hbar\bm{U}\bm{V}^{\mathsf{T}}$, and $\bm{a}$ is an auxiliary complex-
valued vector degree of freedom. Such a factorization is favorable since the
matrix $\bm{R}$ is not full-rank, and can be obtained, e.g., using the
singular value decomposition. The dimension of the auxiliary vector degree of
freedom $\bm{a}$ is given by the rank of the matrix $\bm{R}$. Consequently, we
obtain a time-local action of the QD-MESS equation as

$S_{\text{QD-MESS}}\left[\bm{\phi}_{\mathrm{core}},\bm{a},\bm{b}\right]\\\
=\int\limits_{-\infty}^{+\infty}\mathcal{L}_{\text{QD-
MESS}}\left(\dot{\bm{\phi}}_{\mathrm{core}},\dot{\bm{a}},\dot{\bm{b}},\bm{\phi}_{\mathrm{core}},\bm{a},\bm{b}\right)\;\mathrm{d}t,$
(46)

where the QD-MESS Lagrangian is given by

$\mathcal{L}_{\text{QD-
MESS}}\left(\dot{\bm{\phi}}_{\mathrm{core}},\dot{\bm{a}},\dot{\bm{b}},\bm{\phi}_{\mathrm{core}},\bm{a},\bm{b}\right)\\\
=-V_{\mathrm{J}}(\bm{\psi}_{\mathrm{c}},\bm{\psi}_{\mathrm{q}})+\sum\limits_{jj^{\prime}=0}^{1}\frac{\mathrm{d}^{j}\bm{\phi}_{\mathrm{core,q}}^{\mathsf{T}}}{\mathrm{d}t^{j}}\bm{A}^{(jj^{\prime})}\frac{\mathrm{d}^{j^{\prime}}\bm{\phi}_{\mathrm{core}}}{\mathrm{d}t^{j^{\prime}}}\\\
+\sum\limits_{k}\left(\imath\hbar
a_{k}^{\ast}\dot{a}_{k}-\hbar\omega_{k}^{\mathrm{d}}|a_{k}|^{2}\right.\\\
\left.+\bm{\phi}_{\mathrm{core,q}}^{\mathsf{T}}\bm{u}_{k}^{\mathrm{d}}a_{k}+a_{k}^{\ast}\bm{v}_{k}^{\mathrm{d}\mathsf{T}}\bm{\phi}_{\mathrm{core}}\right)\\\
+\sum\limits_{n}\left(\imath\hbar\bm{b}_{n}^{\dagger}\dot{\bm{b}}_{n}+\imath\omega_{n}^{\mathrm{gm}}|\bm{b}_{n}|^{2}\right.\\\
\left.+\bm{\phi}_{\mathrm{core,q}}^{\mathsf{T}}\bm{U}_{n}^{\mathrm{gm}}\bm{b}_{n}+\bm{b}_{n}^{\dagger}\bm{V}_{n}^{\mathrm{gm}\mathsf{T}}\bm{\phi}_{\mathrm{core,q}}\right),$
(47)

where
$\bm{R}_{k}^{\mathrm{d}}=\hbar\bm{u}_{k}^{\mathrm{d}}\bm{v}_{k}^{\mathrm{d}\mathsf{T}}$,
$\bm{R}_{n}^{\mathrm{gm}}=\hbar\bm{U}_{n}^{\mathrm{gm}}\bm{V}_{n}^{\mathrm{gm}\mathsf{T}}$,
$a_{k}$ and $\bm{b}_{n}$ are auxiliary degrees of freedom corresponding to
dynamical and Matsubara poles of the self-energy, and all the dynamical
variables and sources are given in the time domain. In this Lagrangian all the
auxiliary bosonic modes are coupled to the flux degrees of freedom
$\bm{\psi}_{\mathrm{c}}$ and $\bm{\psi}_{\mathrm{q}}$, hence we refer to the
QD-MESS equation produced by this Lagrangian as flux QD-MESS. However, this
form is not unique. Instead of applying a Hubbard–Stratonovich transformation
as in Eq. (45), we may expand in the integrand

$\frac{1}{\omega-\omega_{\ast}}=-\frac{1}{\omega_{\ast}}+\frac{\omega}{\omega_{\ast}(\omega-\omega_{\ast})}=\\\
-\frac{1}{\omega_{\ast}}-\frac{\omega}{\omega_{\ast}^{2}}+\frac{\omega^{2}}{\omega_{\ast}^{2}(\omega-\omega_{\ast})}.$
(48)

Multiplication by $\omega$ of a Fourier image corresponds to applying the
$\imath\partial_{t}$ operator in time domain. Hence we can express the time
non-local term as

$\exp\left[-\frac{\imath}{2\pi\hbar}\int\limits_{-\infty}^{+\infty}\bm{y}^{\dagger}(\omega)\frac{\bm{R}}{\omega-\omega_{\ast}}\bm{x}(\omega)\;\mathrm{d}\omega\right]\\\
=\int\mathcal{D}[\bm{a},\bm{a}^{\dagger}]\exp\left\\{\frac{\imath}{\hbar}\int\limits_{-\infty}^{+\infty}\left[\imath\hbar\bm{a}^{\dagger}(t)\dot{\bm{a}}(t)\right.\right.\\\
\left.\left.+\hbar\omega_{\ast}|\bm{a}(t)|^{2}+\bm{y}^{\mathsf{T}}(t)\bm{R}\left(\frac{1}{\omega_{\ast}}+\frac{i\partial_{t}}{\omega_{\ast}^{2}}\right)\bm{x}(t)\right.\right.\\\
\left.\phantom{\int}\left.+\bm{a}^{\dagger}(t)\bm{V}^{\mathsf{T}}\dot{\bm{x}}(t)+\dot{\bm{y}}^{\mathsf{T}}(t)\bm{U}\bm{a}(t)\right]\;\mathrm{d}t\right\\}.$
(49)

This way, with a proper renormalization of the time-local terms in the action
(43), one can obtain a QD-MESS Lagrangian where the auxiliary degrees of
freedom couple to the voltages $\dot{\bm{\psi}}_{\mathrm{c}}$ and
$\dot{\bm{\psi}}_{\mathrm{q}}$ for some of the modes (see Appendix A for
details). This formulation may be beneficial for practical calculations, but
for the sake of simplicity of the presentation, we proceed with the flux QD-
MESS Lagrangian.

###  III.5 QD-MESS equation

To derive the quantum equation which corresponds to the Lagrangian (47), we
proceed with a standard transition to the Schrödinger picture. First, we
introduce the canonical conjugate variables to the fluxes, i.e., the charges

$\begin{gathered}\bm{q}_{\mathrm{c}}=\frac{\partial\mathcal{L}_{\text{QD-
MESS}}}{\partial\dot{\bm{\psi}}_{\mathrm{q}}^{\mathsf{T}}},\\\
\bm{q}_{\mathrm{q}}^{\mathsf{T}}=\frac{\partial\mathcal{L}_{\text{QD-
MESS}}}{\partial\dot{\bm{\psi}}_{\mathrm{c}}}\end{gathered}$ (50)

Then, we apply the Legendre transformation associated with Eq. (50) and obtain
the Liouvillian

$\mathfrak{L}=\bm{q}_{\mathrm{q}}^{\mathsf{T}}\dot{\bm{\psi}}_{\mathrm{c}}+\dot{\bm{\psi}}_{\mathrm{q}}^{\mathsf{T}}\bm{q}_{\mathrm{c}}+\imath\hbar\sum\limits_{k}a_{k}^{\ast}\dot{a}_{k}+\imath\hbar\sum\limits_{n}\bm{b}_{n}^{\dagger}\dot{\bm{b}}_{n}-\mathcal{L}_{\mathrm{QD-
MESS}}.$ (51)

This Liouvillian, when written as a function of canonical variables, can be
interpreted as a quantum Liouvillian by imposing the canonical commutation
relations

$\begin{gathered}\left[\check{\psi}_{\mathrm{c},\alpha},\check{q}_{\mathrm{q},\alpha^{\prime}}\right]=\imath\hbar\delta_{\alpha\alpha^{\prime}},~{}\left[\check{\psi}_{\mathrm{q},\alpha},\check{q}_{\mathrm{c},\alpha^{\prime}}\right]=\imath\hbar\delta_{\alpha\alpha^{\prime}},\\\
\left[\check{a}_{k},\check{a}_{k^{\prime}}^{\dagger}\right]=\delta_{kk^{\prime}},~{}\left[\check{b}_{n,\alpha},\check{b}_{n^{\prime},\alpha^{\prime}}^{\dagger}\right]=\delta_{nn^{\prime}}\delta_{\alpha\alpha^{\prime}},\end{gathered}$
(52)

with all other single-operator commutators vanishing. Finally, the quantum
Liouvillian, which governs the dynamics of the superconducting quantum
circuit, reads

$\check{\mathfrak{L}}=\check{\bm{x}}_{\mathrm{q}}^{\mathsf{T}}\bm{L}_{\mathrm{quad}}\check{\bm{x}}+V_{\mathrm{J}}\left(\check{\bm{\psi}}_{\mathrm{c}},\check{\bm{\psi}}_{\mathrm{q}}\right)\\\
+\sum\limits_{j}\left(\hbar\omega_{k}^{\mathrm{d}}\check{a}_{k}^{\dagger}\check{a}_{k}-\check{a}_{k}^{\dagger}\bm{v}_{k}^{\mathrm{d}\mathsf{T}}\check{\bm{\phi}}_{\mathrm{core}}-\check{\bm{\phi}}_{\mathrm{core,q}}^{\mathsf{T}}\bm{u}_{k}^{\mathrm{d}}\check{a}_{k}\right)\\\
-\sum\limits_{n}\left(\imath\hbar\omega_{n}^{\mathrm{gm}}\check{\bm{b}}^{\dagger}_{n}\check{\bm{b}}_{n}+\check{\bm{b}}_{n}\bm{V}_{n}^{\mathrm{gm}\mathsf{T}}\check{\bm{\phi}}_{\mathrm{core,q}}+\check{\bm{\phi}}_{\mathrm{core,q}}^{\mathsf{T}}\bm{U}_{n}^{\mathrm{gm}}\check{\bm{b}}_{n}\right),$
(53)

where $\bm{L}_{\mathrm{quad}}$ is a constant matrix, which describes the
coupling between the classical and the quantum charge and flux operators and
their coupling to the sources,
$\check{\bm{x}}=\begin{bmatrix}\check{\bm{q}}_{\mathrm{c}}^{\mathsf{T}}&\check{\bm{\psi}}_{\mathrm{c}}^{\mathsf{T}}&\dot{\bm{v}}_{\mathrm{i}}^{\mathsf{T}}&\bm{v}_{\mathrm{i}}^{\mathsf{T}}&\check{\bm{q}}_{\mathrm{q}}^{\mathsf{T}}&\check{\bm{\psi}}_{\mathrm{q}}^{\mathsf{T}}&\dot{\bm{\eta}}_{\mathrm{o}}^{\mathsf{T}}&\bm{\eta}_{\mathrm{o}}^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}$
is a vector of all charges, fluxes, and sources, while
$\check{\bm{x}}_{\mathrm{q}}=\begin{bmatrix}\check{\bm{q}}_{\mathrm{q}}^{\mathsf{T}}&\check{\bm{\psi}}_{\mathrm{q}}^{\mathsf{T}}&\dot{\bm{\eta}}_{\mathrm{o}}^{\mathsf{T}}&\bm{\eta}_{\mathrm{o}}^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}$
is its quantum component,
$\check{\bm{\phi}}_{\mathrm{core}}=\begin{bmatrix}\check{\bm{\psi}}_{\mathrm{c}}^{\mathsf{T}}&\bm{v}_{\mathrm{i}}^{\mathsf{T}}&\check{\bm{\psi}}_{\mathrm{q}}^{\mathsf{T}}&\bm{\eta}_{\mathrm{o}}^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}$
is a vector of all fluxes and sources, and
$\check{\bm{\phi}}_{\mathrm{core,q}}=\begin{bmatrix}\check{\bm{\psi}}_{\mathrm{q}}^{\mathsf{T}}&\bm{\eta}_{\mathrm{o}}^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}$
is its quantum component. Note that the ordering of the operators is such that
the operators corresponding to the quantum degrees of freedom lie at the left,
similarly to the creation operators of the auxiliary bosonic modes. The
obtained Liouvillian by its form resembles a Hamiltonian of a bosonic impurity
problem, where a nonlinear impurity is coupled to a continuum of harmonic
modes. The most studied case in this context is the paradigmatic spin-boson
problem [4]. An important difference to the situation is though that here the
harmonic modes are discrete and have complex frequencies and couplings.

The QD-MESS equation, which governs the dynamics of the dissipative quantum-
electric circuit, reads as

$\imath\hbar\frac{\mathrm{d}}{\mathrm{d}t}|W\rangle=\left.\check{\mathfrak{L}}\right|_{\bm{\eta}_{\mathrm{o}}=0}|W\rangle,$
(54)

where $|W\rangle$ is a multidimensional vector in an extended Liouville space,
which completely describes a state of both system and its environment.
Equation (54) has a conservation law

$\langle 1|\frac{\mathrm{d}}{\mathrm{d}t}|W\rangle=0,$ (55)

where $\langle 1|$ is a left zero eigenvector of the quantum operators of all
dynamical variables $\check{\bm{\psi}}_{\mathrm{q}}$ and
$\check{\bm{q}}_{\mathrm{q}}$, and of the creation operators
$\check{\bm{a}}^{\dagger}$ and $\check{\bm{b}}^{\dagger}$ of all the auxiliary
modes, i.e. it corresponds to zero occupation number of each of these modes.
Consequently, the proper ordering of the operators in the quantum Liouvillian
(53) results in

$\langle 1|\left.\check{\mathfrak{L}}\right|_{\bm{\eta}_{\mathrm{o}}=0}=0.$
(56)

The product $\langle 1|W\rangle$ is equivalent to the trace operation of the
full system-environment density operator, hence Eq. (55) is equivalent to the
conservation of probability for the standard Liouville–von Neumann equation.
Thus, it is natural to choose a normalization of both $|W\rangle$ and $\langle
1|$ such that $\langle 1|W\rangle=1$. The thermal equilibrium state
corresponds to a right zero eigenvector of the Liouvillian in the absence of
classical and quantum sources fields

$\left.\check{\mathfrak{L}}\right|_{\bm{\varphi}_{\mathrm{i}}=0,~{}\bm{\eta}_{\mathrm{o}}=0}|W_{\mathrm{eq}}\rangle=0.$
(57)

The physical meaning of the components of $|W\rangle$ may be non-trivial to
identify, especially given that the choice of the degrees of freedom included
in the system and in the environment may be arbitrary. This raises two issues:
the preparation of the initial state for the simulations and the calculation
of the probability distribution of observables. The first issue can be
overcome by an explicit simulation of the initialization pulse, assuming the
initial state to be the equilibrium state $|W_{\mathrm{eq}}\rangle$. Below, we
introduce a method how to calculate mutual quasiprobability distributions of
various observables.

###  III.6 Observables

Let us focus on a case where we apply an input field and measure the output
field within a time interval $(t^{\mathrm{i}},t^{\mathrm{f}})$. Before
$t^{\mathrm{i}}$, the open system is in the thermal state
$|W_{\mathrm{eq}}\rangle$. The result of a measurement is a time trace of the
output field $\bm{v}_{\mathrm{o}}(t)$ sampled from a stochastic process. The
observables $\bm{o}$ we are interested in are weighted averages of these
fields over the above-mentioned time interval, and hence given by

$\bm{o}=\int\limits_{t^{\mathrm{i}}}^{t^{\mathrm{f}}}{\bm{F}}(t)\bm{v}_{\mathrm{o}}(t)\;\mathrm{d}t,$
(58)

where $\bm{F}(t)$ is a matrix-valued weight function. Since the output field
is stochastic, these observables $\bm{o}$ are random variables. Their mutual
quasiprobability distribution can be formally found as [11]

$p(\bm{o})=\left.\left\\{\delta\left[-\imath\hbar\int\limits_{t^{\mathrm{i}}}^{t^{\mathrm{f}}}\bm{F}(t)\frac{\delta}{\delta\bm{\eta}_{\mathrm{o}}^{\mathsf{T}}(t)}\;\mathrm{d}t-\bm{o}\right]\mathcal{Z}[\bm{v}_{\mathrm{i}},\bm{\eta}_{\mathrm{o}}]\right\\}\right|_{\bm{\eta}_{\mathrm{o}}=0},$
(59)

where the argument of the delta function is a linear functional acting on the
generating functional
$\mathcal{Z}[\bm{v}_{\mathrm{i}},\bm{\eta}_{\mathrm{o}}]$. The latter itself
can be calculated by solving the QD-MESS equation (54) with a non-zero quantum
source

$\mathcal{Z}[\bm{v}_{\mathrm{i}},\bm{\eta}_{\mathrm{o}}]=\langle
1|\mathbb{T}\exp\left[-\frac{\imath}{\hbar}\int\limits_{t^{\mathrm{i}}}^{t^{\mathrm{f}}}\check{\mathfrak{L}}(t)\;\mathrm{d}t\right]|W_{\mathrm{eq}}\rangle,$
(60)

where $\mathbb{T}$ is a time-ordering symbol. By employing a Fourier
transformation of a multidimensional delta function, we find that the
quasiprobability density can be calculated as

$p(\bm{o})=\int\mathrm{e}^{\imath\bm{s}^{\mathsf{T}}\bm{o}}\mathcal{Z}[\bm{v}_{\mathrm{i}},\hbar\bm{F}^{\mathsf{T}}(t)\bm{s}]\;\frac{\mathrm{d}\bm{s}}{(2\pi)^{\operatorname{dim}(\bm{o})}}.$
(61)

Since $\mathcal{Z}[\bm{v}_{\mathrm{i}},0]=1$ for an arbitrary drive field
$\bm{v}_{\mathrm{i}}$, the quasiprobability density $p(\bm{o})$ integrates to
unity. However, it is not clear at this point whether it satisfies the non-
negativity requirement for a probability density. For example, the Wigner
quasiprobability function, which may have negative values, falls into this
class of distribution functions if one chooses $\bm{q}_{\mathrm{c}}$ and
$\bm{\psi}_{\mathrm{c}}$ as observables instead of the time-averaged output
fields.

##  IV Weak-coupling regime

Even though it is possible to write a time-local equation of motion which
fully describes the dynamics of the dissipative quantum circuit, solving this
equation remains challenging due to the exponentially high dimension of the
state vector $|W\rangle$. This issue is pronounced in the ultra-low
temperature domain, where the generalized Matsubara frequencies accumulate to
zero.

Curiously, the QD-MESS equation written explicitly in Fock basis of the
auxiliary modes, is identical to HEOM [47, 75] which can be solved using
tensor network algorithms [89, 90]. Here, we suggest a systematic approach to
eliminate weakly coupled auxiliary modes, which results in a substantial
reduction of dimensionality required for the Liouville space of linear degrees
of freedom. The general form of the QD-MESS Liouvillian reads as

$\check{\mathfrak{L}}(t)=\check{\mathfrak{L}}^{(0)}(t)+\\\
\sum\limits_{k}\left\\{\left[\check{x}_{k}^{\prime}+\alpha_{k}^{\prime}(t)\right]\check{a}_{k}+\left[\check{x}_{k}^{\prime\prime}+\alpha_{k}^{\prime\prime}(t)\right]\check{a}_{k}^{\dagger}+\hbar\omega_{k}\check{a}_{k}^{\dagger}\check{a}_{k}\right\\},$
(62)

where we denote the Liouvillian of all non-linear degrees of freedom as
$\check{\mathfrak{L}}^{(0)}(t)$, unify the dynamical and generalized Matsubara
modes, denote all the coupling operators as $\check{x}_{k}^{\prime}$ and
$\check{x}_{k}^{\prime\prime}$, and introduce the drive variables
$\alpha_{k}^{\prime}(t)$ and $\alpha_{k}^{\prime\prime}(t)$ which may
originate from both, classical $\bm{v}_{\mathrm{i}}$ and quantum
$\bm{\eta}_{\mathrm{o}}$ sources. Note that this form is not unique since one
can incorporate some of the auxiliary modes to
$\check{\mathfrak{L}}^{(0)}(t)$.

First, we eliminate time-dependent source terms by carrying out a shift
transformation

$|W\rangle=e^{\check{T}(t)}|W_{\mathrm{shift}}\rangle,$ (63)

where
$\check{T}(t)=\sum_{k}\left[\beta_{k}^{\prime}(t)\check{a}_{k}+\beta_{k}^{\prime\prime}(t)\check{a}_{k}^{\dagger}\right]$
and the time-dependent shifts $\beta_{k}^{\prime}(t)$ and
$\beta_{k}^{\prime\prime}(t)$ are required to satisfy the differential
equations

$\begin{gathered}\imath\hbar\dot{\beta}_{k}^{\prime}(t)=-\hbar\omega_{k}\beta_{k}^{\prime}(t)+\alpha_{k}^{\prime}(t),\\\
\imath\hbar\dot{\beta}_{k}^{\prime\prime}(t)=\hbar\omega_{k}\beta_{k}^{\prime\prime}(t)+\alpha_{k}^{\prime\prime}(t).\\\
\end{gathered}$ (64)

Thus, the dynamics of $|W_{\mathrm{shift}}\rangle$ is governed by the
following Liouvillian

$\check{\mathfrak{L}}_{\mathrm{shift}}(t)=\check{\mathfrak{L}}^{(0)}_{\mathrm{shift}}(t)+\\\
\sum\limits_{k}\left(\check{x}_{k}^{\prime}\check{a}_{k}+\check{x}_{k}^{\prime\prime}\check{a}_{k}^{\dagger}+\hbar\omega_{k}\check{a}_{k}^{\dagger}\check{a}_{k}\right),$
(65)
$\check{\mathfrak{L}}^{(0)}_{\mathrm{shift}}(t)=\check{\mathfrak{L}}^{(0)}(t)+\sum\limits_{k}\left[\check{x}_{k}^{\prime}(t)\beta_{k}^{\prime\prime}(t)-\check{x}_{k}^{\prime\prime}(t)\beta_{k}^{\prime}(t)\right]\\\
+\frac{1}{2}\sum\limits_{k}\left[\alpha_{k}^{\prime}(t)\beta_{k}^{\prime\prime}(t)-\alpha_{k}^{\prime\prime}(t)\beta_{k}^{\prime}(t)\right].$
(66)

Up to this point, we have not done any approximations to arrive at the above
Liouvillian since the shift transformation is exact. Next, we proceed with
eliminating the transverse coupling to the auxiliary modes by a dynamical
Schrieffer–Wolf transformation [78, 79]

$|W_{\mathrm{shift}}\rangle=e^{\check{S}(t)}|W_{\mathrm{SW}}\rangle,$ (67)

where
$\check{S}(t)=\sum_{k}\left[\check{p}_{k}^{\prime}(t)\check{a}_{k}+\check{p}_{k}^{\prime\prime}(t)\check{a}_{k}^{\dagger}\right]$
and operators $\check{p}_{k}^{\prime}(t)$ and
$\check{p}_{k}^{\prime\prime}(t)$ act in the space of the Liouvillian
$\check{\mathfrak{L}}_{\mathrm{shift}}^{(0)}(t)$. We emphasize that this
Schrieffer–Wolf transformation is performed at the Liouvillian level and is
non-unitary. The first-order expansion of the exponents with respect to
operators $\check{p}_{k}^{\prime}(t)$ and $\check{p}_{k}^{\prime\prime}(t)$
reveals that the transverse coupling may be eliminated if
$\check{p}_{k}^{\prime}(t)$ and $\check{p}_{k}^{\prime\prime}(t)$ satisfy the
following equations:

$\begin{gathered}\imath\hbar\frac{\mathrm{d}\check{p}_{k}^{\prime}(t)}{\mathrm{d}t}=\left[\check{\mathfrak{L}}_{\mathrm{shift}}^{(0)}(t),\check{p}_{k}^{\prime}(t)\right]-\hbar\omega_{k}\check{p}_{k}^{\prime}(t)+\check{x}_{k}^{\prime},\\\
\imath\hbar\frac{\mathrm{d}\check{p}_{k}^{\prime\prime}(t)}{\mathrm{d}t}=\left[\check{\mathfrak{L}}_{\mathrm{shift}}^{(0)}(t),\check{p}_{k}^{\prime\prime}(t)\right]+\hbar\omega_{k}\check{p}_{k}^{\prime\prime}(t)+\check{x}_{k}^{\prime\prime}.\end{gathered}$
(68)

The second-order expansion yields the Schrieffer–Wolf correction to the
Liouvillian

$\check{\mathfrak{L}}_{\mathrm{SW}}\approx\check{\mathfrak{L}}^{(0)}_{\mathrm{shift}}(t)+\\\
\sum\limits_{k}\left\\{\hbar\omega_{k}\check{a}_{k}^{\dagger}\check{a}_{k}+\frac{1}{2}\left[\check{x}_{k}^{\prime}\check{a}_{k}+\check{x}_{k}^{\prime\prime}\check{a}_{k}^{\dagger},\check{S}(t)\right]\right\\}.$
(69)

Next, we apply a partial rotating-wave approximation for the auxiliary degrees
of freedom, keeping only the diagonal terms
$\check{a}_{k}^{\dagger}\check{a}_{k}$. This results in

$\begin{gathered}\check{\mathfrak{L}}_{\mathrm{SW}}(t)\approx\check{\mathfrak{L}}_{\mathrm{SW}}^{(0)}(t)+\sum\limits_{k}\left\\{\hbar\omega_{k}+\left[\check{x}_{k}^{\prime},\check{p}_{k}^{\prime\prime}(t)\right]\right\\}\check{a}_{k}^{\dagger}\check{a}_{k},\\\
\check{\mathfrak{L}}_{\mathrm{SW}}^{(0)}(t)=\check{\mathfrak{L}}_{\mathrm{shift}}^{(0)}(t)+\frac{1}{2}\sum\limits_{k}\left[\check{x}_{k}^{\prime}\check{p}_{k}^{\prime\prime}(t)-\check{p}_{k}^{\prime}(t)\check{x}_{k}^{\prime\prime}\right].\end{gathered}$
(70)

The last step is valid only if there are no pairs $(k,k^{\prime})$ for which
$\omega_{k}\pm\omega_{k^{\prime}}$ is small. Such a delicate case appears, for
example, for two weakly decaying dynamical modes with frequencies
$\omega_{k}=-\omega_{k^{\prime}}^{\ast}$. These pairs should be treated
carefully: One may either keep the off-diagonal terms
$\check{a}_{k}\check{a}_{k^{\prime}}$ and
$\check{a}_{k}^{\dagger}\check{a}_{k^{\prime}}^{\dagger}$ in the transformed
Liouvillian, or take these modes into account exactly by incorporating them
into $\check{\mathfrak{L}}^{(0)}(t)$ in the original Liouvillian (62).

The dynamical Schrieffer–Wolf transformation is perturbative, hence it can be
applied only to systems weakly coupled to the auxiliary modes. This
perturbative treatment is similar to standard Born–Markov approximations
widely used for open quantum systems in a weak coupling regime. However, the
presented approach provides a systematic way to derive dissipative terms in
the Liouvillian in case of complicated circuits, for which employing the
standard IO theory [52, 53, 54] becomes hindered. On top of that, resonances
can be treated accurately by incorporating weakly decaying auxiliary modes
into the Liouvillian $\check{\mathfrak{L}}^{(0)}(t)$.

In addition to the weak coupling requirement, one of the most restrictive
limitations of this approach is that it can be used only for sufficiently
small systems, for which Eq. (68) can be solved exactly. This makes it
challenging for circuits with multiple Josephson junctions or many weakly
decaying dynamical modes. However, accurate solutions of even small circuits
may be of great interest, for example, from the point of view of single-qubit
operation fidelities.

##  V Dispersive readout of a transmon qubit

Figure 2: Schematic of a dispersive readout setup for a transmon qubit. The circuit consists of the qubit, capacitively coupled to a resonator, which in turn capacitively coupled to two semi-infinite transmission lines. Table 1: Parameters of the circuit shown in Fig. 2 used for numerical calculations. $C_{\mathrm{t}}$ | $E_{\mathrm{J}}/h$ |  $C_{\mathrm{r}}$ |  $L_{\mathrm{r}}$ |  $C_{\mathrm{c}}$ |  $C_{\mathrm{o}}$ |  $Z$
---|---|---|---|---|---|---
80 fF | 12.5 GHz |  570 fF |  0.7 nH |  40 fF |  15 fF |  50 $\Omega$

Here, we demonstrate our theory in an experimentally feasible scenario by
applying it to a model of the dispersive readout of a superconducting qubit
[51]. As an example system shown in Fig. 2, we employ a typical transmon qubit
[91] which is capacitively coupled to an off-resonant readout resonator. The
resonator is capacitively coupled to two semi-infinite transmission lines, one
of which is used to drive the readout mode and the other to measure the output
field. No input is applied on the transmission line for the measurement and
the output field to the drive line is disregarded.

Kirchhoff’s law for the circuit shown in Fig. 2 yields

$\begin{gathered}-\left(C_{\mathrm{t}}+C_{\mathrm{c}}\right)\ddot{\varphi}_{\mathrm{t}}+C_{\mathrm{c}}\ddot{\varphi}_{\mathrm{r}}-\frac{2\pi
E_{\mathrm{J}}}{\Phi_{0}}\sin\left(\frac{2\pi\varphi_{\mathrm{t}}}{\Phi_{0}}\right)=0,\\\
-\left(C_{\mathrm{r}}+C_{\mathrm{c}}+C_{\mathrm{o}}\right)\ddot{\varphi}_{\mathrm{r}}+C_{\mathrm{c}}\ddot{\varphi}_{\mathrm{t}}+C_{\mathrm{o}}\ddot{\varphi}_{\mathrm{o}}-\frac{\varphi_{\mathrm{r}}}{L_{\mathrm{r}}}=0,\\\
-C_{\mathrm{o}}\left(\ddot{\varphi}_{\mathrm{o}}-\ddot{\varphi}_{\mathrm{r}}\right)-\frac{2}{Z}\left(\dot{\varphi}_{\mathrm{o}}-v_{\mathrm{i}}\right)=0,\\\
v_{\mathrm{o}}=\dot{\varphi}_{\mathrm{o}}.\end{gathered}$ (71)

The blocks of inverse Green’s functions (30) for this circuit are given by

$\begin{gathered}\bm{K}^{\mathrm{R}}(\omega)=\begin{bmatrix}\omega^{2}\tilde{C}_{\mathrm{t}}&-\omega^{2}C_{\mathrm{c}}&0&0\\\
-\omega^{2}C_{\mathrm{c}}&\omega^{2}\tilde{C}_{\mathrm{r}}-L_{\mathrm{r}}^{-1}&-\omega^{2}C_{\mathrm{o}}&0\\\
0&-\omega^{2}C_{\mathrm{o}}&\omega^{2}C_{\mathrm{o}}+2\imath\omega
Z^{-1}&2Z^{-1}\\\ 0&0&-\imath\omega&0\end{bmatrix},\\\
\bm{K}^{\mathrm{A}}(\omega)=\left[\bm{K}^{\mathrm{R}}(\omega^{\ast})\right]^{\dagger},\\\
\bm{K}^{\mathrm{K}}(\omega)=\begin{bmatrix}0&0&0&0\\\ 0&0&0&0\\\
0&0&2\imath\omega
Z^{-1}\coth\left(\frac{\hbar\omega}{2k_{\mathrm{B}}T}\right)&0\\\
0&0&0&0\end{bmatrix},\end{gathered}$ (72)

where $\tilde{C}_{\mathrm{t}}=C_{\mathrm{t}}+C_{\mathrm{c}}$ and
$\tilde{C}_{\mathrm{r}}=C_{\mathrm{r}}+C_{\mathrm{c}}+C_{\mathrm{o}}$. The
dynamical degree of freedom of the core system is represented by transmon flux
$\psi=\varphi_{\mathrm{t}}$. The rest of the circuit can be described by three
dynamical auxiliary modes, two of them have frequencies close to $\pm
1/\sqrt{L_{\mathrm{r}}\tilde{C}_{\mathrm{r}}}$ with small negative imaginary
parts and one with purely imaginary frequency around
$-2\imath/(ZC_{\mathrm{o}})$. For our calculations we use parameters of the
circuit shown in Tab. 1.

It appears to be more convenient to use coupling to the auxiliary modes
through voltage $\dot{\psi}$ instead of flux coupling shown in Eq. (47). This
leads to the coupling through the charge of the transmon island to the
auxiliary modes in Liouvillian. More details can be found in Appendix A. Since
the coupling is weak, we may apply the dynamical Schrieffer–Wolf approach
developed in Sec. IV for all of the generalized Matsubara modes and the
dynamical mode with imaginary frequency. We restrict the linear space of
$\check{\mathfrak{L}}^{(0)}(t)$ and the coupling operators in the Liouvillian
(62) to the four lowest levels for the isolated transmon and six levels for
the weakly dissipative auxiliary modes. See more technical details on the
construction of the Liouvillian in Appendix B.

###  V.1 Transmission coefficient

Figure 3: Absolute value of the transmission coefficient $S_{21}(\omega)$ for (a) cold transmon at temperature $T=0.01$ K and (b) hot transmon at temperature $T=0.1$ K. The parameters of the circuit are given in Tab. 1. Frequencies and linewidths of the qubit and the resonator, qubit anharmonicity and dispersive shift are shown in Tab. 2. Table 2: Frequency $\omega_{\mathrm{t}}$, linewidth $\gamma_{\mathrm{t}}$, and anharmonicity $\alpha$ of the transmon, frequency $\omega_{\mathrm{r}}$ and linewidth $\gamma_{\mathrm{r}}$ of the resonator, and the dispersive shift $\chi$, obtained from the diagonalization of the system Liouvillian. The qubit linewidth is given at low (top cell) and high (bottom cell) temperatures. All the quantities are given in GHz. $\omega_{\mathrm{t}}/(2\pi)$ |  $\gamma_{\mathrm{t}}/(2\pi)$ | $\alpha/(2\pi)$ | $\omega_{\mathrm{r}}/(2\pi)$ |  $\gamma_{\mathrm{r}}/(2\pi)$ |  $\chi/(2\pi)$
---|---|---|---|---|---
3.834 |  $10^{-6}$ (0.01 K) | $-0.175$ | 7.719 |  $1.71\cdot 10^{-3}$ |  $-6.4\cdot 10^{-3}$
$6.9\cdot 10^{-5}$ (0.1 K) |  |

As the first step we probe the spectrum of the device by calculating the
transmission coefficient $S_{21}(\omega)$. It is defined as $\langle
v_{\mathrm{o}}(\omega)\rangle=S_{21}(\omega)v_{\mathrm{i}}(\omega)$ for
infinitesimal drive $v_{\mathrm{i}}$, where the average is taken over the
realizations of the output field. It can be evaluated as [11]

$S_{21}(\omega)=\left.-\imath\hbar\frac{\delta^{2}Z[v_{\mathrm{i}},\eta_{\mathrm{o}}]}{\delta
v_{\mathrm{i}}(\omega)\delta\eta_{\mathrm{o}}^{\ast}(\omega)}\right|_{v_{\mathrm{i}}=0,~{}\eta_{\mathrm{o}}=0}.$
(73)

In order to calculate it, we apply a Schrieffer–Wolf transformation for
$v_{\mathrm{i}}=0$ and $\eta_{\mathrm{o}}=0$ to the static Liouvillian and to
the drive terms. The resulting Liouvillian assumes the form

$\check{\mathfrak{L}}_{\mathrm{SW}}=\check{\mathfrak{L}}_{\mathrm{SW}}^{(\mathrm{static})}\\\
+\begin{bmatrix}\eta_{\mathrm{o}}&\dot{\eta}_{\mathrm{o}}&1\end{bmatrix}\begin{bmatrix}L^{(00)}_{\mathrm{oi}}&L^{(01)}_{\mathrm{oi}}&\check{V}_{\mathrm{c}}^{(0)}\\\
L^{(10)}_{\mathrm{oi}}&L^{(11)}_{\mathrm{oi}}&\check{V}_{\mathrm{c}}^{(1)}\\\
\check{V}_{\mathrm{q}}^{(0)}&\check{V}_{\mathrm{q}}^{(1)}&0\end{bmatrix}\begin{bmatrix}v_{\mathrm{i}}\\\
\dot{v}_{\mathrm{i}}\\\ 1\end{bmatrix}.$ (74)

Here the static Liouvillian is equal to

$\check{\mathfrak{L}}_{\mathrm{SW}}^{(\mathrm{static})}=\left.\check{\mathfrak{L}}_{\mathrm{SW}}\right|_{v_{\mathrm{i}}=0,\eta_{\mathrm{o}}=0},$
(75)

operators which couple system degrees of freedom with the input and output
field are given by

$\begin{gathered}\check{V}_{\mathrm{c}}^{(0)}=\left.\left(\frac{\partial\check{\mathfrak{L}}}{\partial\eta_{\mathrm{o}}}+\left[\frac{\partial\check{\mathfrak{L}}}{\partial\eta_{\mathrm{o}}},\check{S}\right]\right)\right|_{v_{\mathrm{i}}=0,\eta_{\mathrm{o}}=0},\\\
\check{V}_{\mathrm{c}}^{(1)}=\left.\left(\frac{\partial\check{\mathfrak{L}}}{\partial\dot{\eta}_{\mathrm{o}}}+\left[\frac{\partial\check{\mathfrak{L}}}{\partial\dot{\eta}_{\mathrm{o}}},\check{S}\right]\right)\right|_{v_{\mathrm{i}}=0,\eta_{\mathrm{o}}=0},\\\
\check{V}_{\mathrm{q}}^{(0)}=\left.\left(\frac{\partial\check{\mathfrak{L}}}{\partial
v_{\mathrm{i}}}+\left[\frac{\partial\check{\mathfrak{L}}}{\partial
v_{\mathrm{i}}},\check{S}\right]\right)\right|_{v_{\mathrm{i}}=0,\eta_{\mathrm{o}}=0},\\\
\check{V}_{\mathrm{q}}^{(1)}=\left.\left(\frac{\partial\check{\mathfrak{L}}}{\partial\dot{v}_{\mathrm{i}}}+\left[\frac{\partial\check{\mathfrak{L}}}{\partial\dot{v}_{\mathrm{i}}},\check{S}\right]\right)\right|_{v_{\mathrm{i}}=0,\eta_{\mathrm{o}}=0},\end{gathered}$
(76)

and direct coupling between the input and output field is provided by

$\begin{gathered}L_{\mathrm{oi}}^{(00)}=\left.\frac{\partial^{2}\check{\mathfrak{L}}}{\partial\eta_{\mathrm{o}}\partial
v_{\mathrm{i}}}\right|_{v_{\mathrm{i}}=0,\eta_{\mathrm{o}}=0},~{}L_{\mathrm{oi}}^{(01)}=\left.\frac{\partial^{2}\check{\mathfrak{L}}}{\partial\eta_{\mathrm{o}}\partial\dot{v}_{\mathrm{i}}}\right|_{v_{\mathrm{i}}=0,\eta_{\mathrm{o}}=0},\\\
L_{\mathrm{oi}}^{(10)}=\left.\frac{\partial^{2}\check{\mathfrak{L}}}{\partial\dot{\eta}_{\mathrm{o}}\partial
v_{\mathrm{i}}}\right|_{v_{\mathrm{i}}=0,\eta_{\mathrm{o}}=0},~{}L_{\mathrm{oi}}^{(11)}=\left.\frac{\partial^{2}\check{\mathfrak{L}}}{\partial\dot{\eta}_{\mathrm{o}}\partial\dot{v}_{\mathrm{i}}}\right|_{v_{\mathrm{i}}=0,\eta_{\mathrm{o}}=0}.\end{gathered}$
(77)

The latter are c-numbers since the both classical and quantum fields couple
linearly to the system and auxiliary degrees of freedom in the Liouvillian
(53). Here, we drop the terms quadratic in $\eta_{\mathrm{o}}$, since they
contribute only to the fluctuations. We emphasize that the generator of the
Schrieffer–Wolf transformation $\check{S}$ should be assumed for the absent
drive $v_{\mathrm{i}}=0$ and $\eta_{\mathrm{o}}=0$.

Then, the transmission coefficient can be expressed as

$S_{21}(\omega)=\sum\limits_{jj^{\prime}=0}^{1}\imath^{j-j^{\prime}}\omega^{j+j^{\prime}}\langle
1|\check{V}_{\mathrm{c}}^{(j)}\check{\mathfrak{G}}(\omega)\check{V}_{\mathrm{q}}^{(j^{\prime})}|W_{\mathrm{eq}}\rangle\\\
+\sum\limits_{jj^{\prime}=0}^{1}\imath^{j-j^{\prime}}\omega^{j+j^{\prime}}L_{\mathrm{oi}}^{(jj^{\prime})},$
(78)

where
$\check{\mathfrak{G}}(\omega)=\left[\omega-\check{\mathfrak{L}}_{\mathrm{SW}}^{(\mathrm{static})}\right]^{-1}$
which can be efficiently calculated using numerical diagonalization of the
static Liouvillian.

The frequency dependence of the absolute value of $S_{21}(\omega)$ is shown in
Fig. 3 at a low temperature of $T=0.01$ K and a high temperature of $T=0.1$ K.
At the low temperature, the circuit is close to its ground state and there are
only two transitions visible. There is a strong feature at the resonator
frequency $\omega_{\mathrm{r}}$ and a weak signature at the transmon frequency
$\omega_{\mathrm{t}}$. At the high temperature, both qubit and resonator
become thermally populated, and consequently more transitions become visible.
In Fig. 3(b), we can clearly observe lines which correspond to transitions
between different transmon levels as well as the splitting of the resonator
feature caused by the dispersive shift (see the inset). A spurious feature
below the resonator peak arises due to the hard cut-off of the transmon levels
which results in the increased dispersive shift for the highest considered
transmon level. Positions of the peaks and their widths correspond to the
eigenvalues of the Liouvillian
$\check{\mathfrak{L}}_{\mathrm{SW}}^{(\mathrm{static})}$. The results are
shown in Tab. 2. Note that the linewidth of the qubit transition is increased
at high temperature which would result in the lower $T_{2}$ time of the qubit.

###  V.2 Readout of a hot transmon

Figure 4: Quasiprobability density of normalized in-phase and quadrature
components of the output field. Temperature is $T=0.1$ K, drive frequency is
$\omega_{\mathrm{d}}=2\pi\times 7.7146$ GHz, drive amplitude is
$V_{\mathrm{d}}=0.2$ $\mu$V, integration time is $t_{\mathrm{int}}=5.18$
$\mu$s.

Next, we proceed to the simulation of the state of the output field during
transmon readout. For the sake of simplicity, we do not consider the
preparation of the qubit state and the simulation of its subsequent dynamics
with respect to the readout pulse. Instead, we consider a continuous
monochromatic drive
$v_{\mathrm{i}}(t)=V_{\mathrm{d}}\cos\left(\omega_{\mathrm{d}}t\right)$ and
the qubit-resonator system at temperature $T$ relaxed to the steady state with
respect to the drive. Details of the calculation can be found in Appendix C.

The quantity of interest is a homodyne demodulated signal averaged over a
finite time interval of length $t_{\mathrm{int}}$:

$S_{\mathrm{i}}+\imath
S_{\mathrm{q}}=\frac{2}{t_{\mathrm{int}}V_{\mathrm{d}}}\int\limits_{0}^{t_{\mathrm{int}}}v_{\mathrm{o}}(t)e^{\imath\omega_{\mathrm{d}}t}\;\mathrm{d}t.$
(79)

The real and imaginary components $S_{\mathrm{i}}$ and $S_{\mathrm{q}}$
correspond to in-phase and quadrature components of the output field,
respectively, normalized by drive amplitude. These components are measured
also in the single-shot readout experiments although typically using a
weighted time average over the duration of the readout pulse. We calculate the
quasiprobability density $p(S_{\mathrm{i}},S_{\mathrm{q}})$ using Eq. (61) and
show the results and simulation parameters in Fig. 4. We can clearly
distinguish four disk-shaped clouds which correspond to the different transmon
states. The intensities of the clouds are proportional to the thermal
populations of the qubit states. In addition there are faint lines which
connect different clouds. These lines correspond to the spontaneous
transitions between the qubit states which may occur during the integration
time [92]. Similar connecting lines have indeed been observed in experiments
[93].

The value of $t_{\mathrm{int}}$ is pivotal for the accurate readout protocol.
For long integration times, the spontaneous transitions become more and more
probable and in the limit $t_{\mathrm{int}}\to\infty$, the quasiprobability
distribution collapses into a delta function centered at the real and
imaginary parts of the transmission coefficient for the finite-amplitude
drive. In the opposite limit, with the decrease of $t_{\mathrm{int}}$, the
clouds begin to grow and, at some point, overlap. Therefore, there is an
optimal integration time, not too short such that the clouds can still be
distinguishable, and not too long such that probability of spontaneous
transition during the measurement is small.

###  V.3 Connection to measurement theory

The choice of the qubit readout problem, in addition to the purely
illustrative purpose, is also motivated by the problem of the quantum
measurement. The qubit readout falls into a class of continuous inefficient
measurements [49] since the state of the qubit is not probed directly, but
through the averaging of the output field over a finite time interval. In a
contrast to the standard presentation of measurements [49], we do not consider
an initial product state of the system and the environment or the measurement
device, but a thermal initial state. Moreover, we do not separate the system
and the environment and study the open quantum system as a whole. Such an
approach may be favorable outside of the weak coupling regime, where product
initial states gives rise to non-exponential short-time decoherence [94, 95].

##  VI Conclusions

We presented an exact nonlinear response theory, free of the limiting
requirements of weak coupling, low temperature, and weak drive. The
dissipation and the driving field introduced by the coupled transmission lines
is reduced to the coupling of the non-linear system degrees of freedom to a
discrete set of auxiliary bosonic modes. We split the modes into two classes:
dynamical and generalized Matsubara modes. The first ones provide coupling to
classical field, while latter are responsible for thermal and quantum
fluctuations in the circuit. This results in the emergence of a time-local QD-
MESS equation which fully describes the open-quantum-system dynamics of the
circuit.

We introduced observables as weighted time averages of the output fields
emitted to the transmission lines and derived an analytical expression of
their mutual quasiprobability distribution based on the Fourier transformation
of the generating functional. This formalism allows to calculate homo- and
heterodyne demodulated signal as well as calculate broadband characteristics
of the output field.

In the weak coupling regime, we employed a dynamical Schrieffer–Wolf
transformation, which provides a consistent way of calculating the dissipators
for arbitrary circuits and driving fields. Finally, we illustrated our
findings by studying a dispersive-readout setup of a superconducting transmon
qubit.

The developed NLR theory provides a powerful framework for the theoretical
analysis of open quantum systems such as superconducting quantum circuits,
including extraction of experimentally observable quantities. Going beyond the
weak-coupling regime can be achievable by solving the QD-MESS equation with
advanced numerical techniques such as tensor networks. In principle, the
developed formalism can be applied not only to lumped-element circuits but
also to distributed elements. The necessary self-energies can be computed
explicitly using high-frequency simulation software for linear components of
the distributed element circuits. This would pave the way for exact quantum
simulations of practical devices.

In the future, the presented formalism can also be modified for the
calculation of thermodynamic properties of the circuit. A non-equilibrium case
can be simulated by introducing generalized Matsubara modes corresponding to
different temperatures. Here, the counting field corresponding to the heat
flow couples to the classical expression for Joule losses which are quadratic
in terms of the microwave field in the circuit. Such a generalization would
allow to study heat transport in non-equilibrium systems and implement
accurate simulations of quantum heat engines.

Among our findings there is a treatment of open quantum systems free from
explicit system-environment separation and from an initial product state
condition. Such a description is valuable for the development of the
measurement theory in the regime of intermediate to strong coupling between
the system and the measurement device.

###### Acknowledgements.

V.V. thanks Kalle Kansanen and Riya Baruah for fruitful discussions. This work
has been financially supported by the Academy of Finland Centre of Excellence
program (project no. 336810) and THEPOW (project no. 349594), the European
Research Council under Advanced Grant no. 101053801 (ConceptQ), by Horizon
Europe programme HORIZON-CL4-2022-QUANTUM-01-SGA via the project 101113946
OpenSuperQPlus100, the German Science Foundation (DFG) under AN336/12-1
(For2724), the State of Baden-Wüttemberg under KQCBW/SiQuRe, and the BMBF
within the QSolid project. We acknowledge the computational resources provided
by the Aalto Science-IT project.

##  Appendix A Charge QD-MESS equation

In this appendix, we provide the charge QD-MESS equation in contrast to the
flux version considered in the main article. We begin with the classical QD-
MESS Lagrangian with the voltage coupling to the auxiliary degrees of freedom,
given by

$\tilde{\mathcal{L}}_{\text{QD-
MESS}}\left(\dot{\bm{\phi}}_{\mathrm{core}},\dot{\bm{a}},\dot{\bm{b}},\bm{\phi}_{\mathrm{core}},\bm{a},\bm{b}\right)=-V_{\mathrm{j}}(\bm{\psi}_{\mathrm{c}},\bm{\psi}_{\mathrm{q}})\\\
+\sum\limits_{jj^{\prime}=0}^{1}\frac{\mathrm{d}^{j}\bm{\phi}_{\mathrm{core,q}}^{\mathsf{T}}}{\mathrm{d}t^{j}}\tilde{\bm{A}}^{(jj^{\prime})}\frac{\mathrm{d}^{j^{\prime}}\bm{\phi}_{\mathrm{core}}}{\mathrm{d}t^{j^{\prime}}}\\\
+\sum\limits_{k}\left(\imath\hbar
a_{k}^{\ast}\dot{a}_{k}-\hbar\omega_{k}^{\mathrm{d}}|a_{k}|^{2}\right.\\\
\left.+\tilde{\bm{\phi}}_{\mathrm{core,q}}^{\mathsf{T}}\tilde{\bm{u}}_{k}^{\mathrm{d}}a_{k}+a_{k}^{\ast}\tilde{\bm{v}}_{k}^{\mathrm{d}\mathsf{T}}\tilde{\bm{\phi}}_{\mathrm{core}}\right)\\\
+\sum\limits_{n}\left(\imath\hbar\bm{b}_{n}^{\dagger}\dot{\bm{b}}_{n}+\imath\omega_{n}^{\mathrm{gm}}|\bm{b}_{n}|^{2}\right.\\\
\left.+\tilde{\bm{\phi}}_{\mathrm{core,q}}^{\mathsf{T}}\tilde{\bm{U}}_{n}^{\mathrm{gm}}\bm{b}_{n}+\bm{b}_{n}^{\dagger}\tilde{\bm{V}}_{n}^{\mathrm{gm}\mathsf{T}}\tilde{\bm{\phi}}_{\mathrm{core,q}}\right),$
(80)

where
$\tilde{\bm{\phi}}_{\mathrm{core}}=\begin{bmatrix}\dot{\bm{\psi}}_{\mathrm{c}}^{\mathsf{T}}&\bm{v}_{\mathrm{i}}^{\mathsf{T}}&\dot{\bm{\psi}}_{\mathrm{q}}^{\mathsf{T}}&\bm{\eta}_{\mathrm{o}}^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}$
and
$\tilde{\bm{\phi}}_{\mathrm{core,q}}=\begin{bmatrix}\dot{\bm{\psi}}_{\mathrm{q}}^{\mathsf{T}}&\bm{\eta}_{\mathrm{o}}^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}$.
The expressions for the classical and quantum charges (50), Liouvillian (51)
and the commutation relations (52) remain valid but the part of the
Liouvillian which describes the dynamics of the auxiliary modes becomes non-
diagonal. The form of Eq. (62) can be obtained by applying the Bogoliubov
transformation to the operators of the auxiliary modes. The difference from
the flux QD-MESS Liouvillian is that the operators of the auxiliary modes are
coupled to the charge degrees of freedom.

The main advantage of using the charge QD-MESS equation for the transmon
readout circuit shown in Fig. 2 is that it evades spurious inductive terms in
the Liouvillian which would break quasi-charge conservation of the transmon
island.

##  Appendix B Discrete approximation of QD-MESS for the transmon

Direct numerical simulation of the QD-MESS equation is impossible, since even
the Liouvillian $\check{\mathfrak{L}}^{(0)}(t)$ in Eq. (62) is infinite-
dimensional and one needs to come up with its finite-dimensional approximation
in some suitable basis. For a general circuit with intermediate or strong
dissipation it may be a challenging problem. For the transmon readout circuit
we use, the coupling of the qubit to the resonator and to the transmission
lines is weak, and hence we can use a finite number of bare transmon basis
states for our calculations. The dissipation-free part of the transmon
Liouvillian reads as

$\check{\mathfrak{L}}_{\mathrm{t}}=\frac{\check{q}_{\mathrm{q}}\check{q}_{\mathrm{c}}}{C_{\mathrm{eff}}}+2E_{\mathrm{j}}\sin\left(\frac{2\pi\check{\psi}_{\mathrm{c}}}{\Phi_{0}}\right)\sin\left(\frac{\pi\check{\psi}_{\mathrm{q}}}{\Phi_{0}}\right).$
(81)

We make a substitution

$\begin{gathered}\check{\psi}_{\mathrm{c}}=\frac{\check{\psi}_{+}+\check{\psi}_{-}}{2},~{}\check{q}_{\mathrm{c}}=\frac{\check{q}_{+}+\check{q}_{-}}{2},\\\
\check{\psi}_{\mathrm{q}}=\check{\psi}_{+}-\check{\psi}_{-},~{}\check{q}_{\mathrm{q}}=\check{q}_{+}-\check{q}_{-},\end{gathered}$
(82)

which is a transformation from the classical and quantum variables to the
forward and backward-path variables. In terms of these variables, the
Liouvillian assumes the form

$\begin{gathered}\check{\mathfrak{L}}_{\mathrm{t}}=\check{H}_{+}-\check{H}_{-},\\\
\check{H}_{\pm}=\frac{\check{q}_{\pm}^{2}}{2C_{\mathrm{eff}}}-E_{\mathrm{j}}\cos\left(\frac{2\pi\check{\psi}_{\pm}}{\Phi_{0}}\right).\end{gathered}$
(83)

This Liouvillian is nothing but a commutator with a standard transmon
Hamiltonian

$\begin{gathered}\check{\mathfrak{L}}_{\mathrm{t}}\operatorname{vec}\left(\hat{\rho}\right)=\operatorname{vec}\left(\left[\hat{H},\hat{\rho}\right]\right),\\\
\hat{H}=\frac{\hat{q}^{2}}{2C_{\mathrm{eff}}}-E_{\mathrm{j}}\cos\left(\frac{2\pi\hat{\psi}}{\Phi_{0}}\right),\end{gathered}$
(84)

where $\operatorname{vec}\left(\hat{\rho}\right)$ stands for the vectorization
of the density operator. For the latter, we can construct a finite-difference
approximation, diagonalize it numerically and use a few lowest transmon levels
$|n\rangle$ with energies $E_{n}$ as a basis. Consequently, an arbitrary
operator $\check{A}_{\pm}$ acting in the Liouville space can be expressed
through its counterpart $\hat{A}$ acting in the Hilbert space as

$\check{A}_{+}\tilde{\rightarrow}\hat{A}\otimes\hat{1},~{}\check{A}_{-}\tilde{\rightarrow}\hat{1}\otimes\hat{A}^{\mathsf{T}},$
(85)

where $\tilde{\rightarrow}$ stands for isomorphism relation, $\otimes$ denotes
Kronecker product, and $\hat{1}$ is the identity operator in the Hilbert
space. For the bosonic operators of the auxiliary modes, which we do not
eliminate with dynamical Schrieffer–Wolf transformation, we cut the Liouville
space with a few lowest levels. For our calculations, we use four levels of
the transmon and six levels of the auxiliary modes.

##  Appendix C QD-MESS equation for the periodically driven system

For the transmon readout we solve the QD-MESS equation with the following
classical and quantum sources:

$\begin{gathered}v_{\mathrm{i}}(t)=V_{\mathrm{d}}\cos(\omega_{\mathrm{d}}t),\\\
\eta_{\mathrm{o}}(s_{1},s_{2},t)=s_{1}\cos(\omega_{\mathrm{d}}t)+s_{2}\sin(\omega_{\mathrm{d}}t),\end{gathered}$
(86)

for different values of the parameters $s_{1}$ and $s_{2}$. The total
Liouvillian is periodic in time, and hence Floquet theory can be applied to
such a system. We consider the evolution operator over a single period as

$\check{U}(s_{1},s_{2})=\mathbb{T}\exp\left[-\frac{\imath}{\hbar}\int\limits_{0}^{{2\pi}/{\omega_{\mathrm{d}}}}\check{\mathfrak{L}}(t)\;\mathrm{d}t\right].$
(87)

The steady state of the system corresponds to the unit eigenstate of this
operator in the absence of the quantum drive

$\check{U}(0)|W_{\mathrm{ss}}\rangle=|W_{\mathrm{ss}}\rangle.$ (88)

We consider that the integration time $t_{\mathrm{int}}$ is an integer
multiple $N$ of the drive period $2\pi/\omega_{\mathrm{d}}$. Thus the
characteristic function can be evaluated as

$\mathcal{Z}(s_{1},s_{2})=\langle
1|\left[\check{U}\left(\frac{s_{1}}{N},\frac{s_{2}}{N}\right)\right]^{N}|W_{\mathrm{ss}}\rangle.$
(89)

These quantities can be efficiently evaluated using the dynamical
Schrieffer–Wolf approach for weakly coupled auxiliary modes.

## References

  * Gardiner and Zoller [2004] C. Gardiner and P. Zoller,  _Quantum noise: a handbook of Markovian and non-Markovian quantum stochastic methods with applications to quantum optics_ (Springer Science & Business Media, 2004).
  * Breuer  _et al._ [2002] H.-P. Breuer, F. Petruccione, _et al._ ,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Weiss [2012] U. Weiss,  _Quantum dissipative systems_ (World Scientific, 2012).
  * Leggett  _et al._ [1987] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. Fisher, A. Garg, and W. Zwerger, Dynamics of the dissipative two-state system, Reviews of Modern Physics 59, 1 (1987).
  * Schön and Zaikin [1990] G. Schön and A. D. Zaikin, Quantum coherent effects, phase transitions, and the dissipative dynamics of ultra small tunnel junctions, Physics Reports 198, 237 (1990).
  * Zurek [2003] W. H. Zurek, Decoherence, einselection, and the quantum origins of the classical, Reviews of modern physics 75, 715 (2003).
  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang,  _Quantum computation and quantum information_ (Cambridge university press, 2010).
  * Averin  _et al._ [2012] D. V. Averin, B. Ruggiero, and P. Silvestrini,  _Macroscopic quantum coherence and quantum computing_ (Springer Science & Business Media, 2012).
  * Beige  _et al._ [2000] A. Beige, D. Braun, B. Tregenna, and P. L. Knight, Quantum computing using dissipation to remain in a decoherence-free subspace, Physical review letters 85, 1762 (2000).
  * Verstraete  _et al._ [2009] F. Verstraete, M. M. Wolf, and J. Ignacio Cirac, Quantum computation and quantum-state engineering driven by dissipation, Nature physics 5, 633 (2009).
  * Kamenev [2023] A. Kamenev,  _Field theory of non-equilibrium systems_ (Cambridge University Press, 2023).
  * Scully and Zubairy [1999] M. O. Scully and M. S. Zubairy, Quantum optics (1999).
  * Georgescu [2020] I. Georgescu, 25 years of quantum error correction, Nature Reviews Physics 2, 519 (2020).
  * Preskill [2018] J. Preskill, Quantum computing in the nisq era and beyond, Quantum 2, 79 (2018).
  * Reiter  _et al._ [2017] F. Reiter, A. S. Sørensen, P. Zoller, and C. Muschik, Dissipative quantum error correction and application to quantum sensing with trapped ions, Nature communications 8, 1822 (2017).
  * Tan  _et al._ [2017] K. Y. Tan, M. Partanen, R. E. Lake, J. Govenius, S. Masuda, and M. Möttönen, Quantum-circuit refrigerator, Nature communications 8, 1 (2017).
  * Ilias  _et al._ [2022] T. Ilias, D. Yang, S. F. Huelga, and M. B. Plenio, Criticality-enhanced quantum sensing via continuous measurement, PRX Quantum 3, 010354 (2022).
  * Ronzani  _et al._ [2018] A. Ronzani, B. Karimi, J. Senior, Y.-C. Chang, J. T. Peltonen, C. Chen, and J. P. Pekola, Tunable photonic heat transport in a quantum heat valve, Nature Physics 14, 991 (2018).
  * Hernández-Gómez  _et al._ [2022] S. Hernández-Gómez, S. Gherardini, N. Staudenmaier, F. Poggiali, M. Campisi, A. Trombettoni, F. Cataliotti, P. Cappellaro, and N. Fabbri, Autonomous dissipative maxwell’s demon in a diamond spin qutrit, PRX Quantum 3, 020329 (2022).
  * Bouton  _et al._ [2021] Q. Bouton, J. Nettersheim, S. Burgardt, D. Adam, E. Lutz, and A. Widera, A quantum heat engine driven by atomic collisions, Nature Communications 12, 2063 (2021).
  * Redfield [1965] A. Redfield, The theory of relaxation processes, in  _Advances in Magnetic and Optical Resonance_ , Vol. 1 (Elsevier, 1965) pp. 1–32.
  * Lindblad [1976] G. Lindblad, On the generators of quantum dynamical semigroups, Communications in Mathematical Physics 48, 119 (1976).
  * Gorini  _et al._ [1976] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Completely positive dynamical semigroups of n-level systems, Journal of Mathematical Physics 17, 821 (1976).
  * Liu  _et al._ [2011] B.-H. Liu, L. Li, Y.-F. Huang, C.-F. Li, G.-C. Guo, E.-M. Laine, H.-P. Breuer, and J. Piilo, Experimental control of the transition from markovian to non-markovian dynamics of open quantum systems, Nature Physics 7, 931 (2011).
  * Smirne  _et al._ [2011] A. Smirne, D. Brivio, S. Cialdi, B. Vacchini, and M. G. Paris, Experimental investigation of initial system-environment correlations via trace-distance evolution, Physical Review A 84, 032112 (2011).
  * Cialdi  _et al._ [2014] S. Cialdi, A. Smirne, M. G. Paris, S. Olivares, and B. Vacchini, Two-step procedure to discriminate discordant from classical correlated or factorized states, Physical Review A 90, 050301 (2014).
  * Gessner  _et al._ [2014] M. Gessner, M. Ramm, T. Pruttivarasin, A. Buchleitner, H.-P. Breuer, and H. Häffner, Local detection of quantum correlations with a single trapped ion, Nature Physics 10, 105 (2014).
  * Tang  _et al._ [2012] J.-S. Tang, C.-F. Li, Y.-L. Li, X.-B. Zou, G.-C. Guo, H.-P. Breuer, E.-M. Laine, and J. Piilo, Measuring non-markovianity of processes with controllable system-environment interaction, Europhysics Letters 97, 10002 (2012).
  * White  _et al._ [2020] G. A. White, C. D. Hill, F. A. Pollock, L. C. Hollenberg, and K. Modi, Demonstration of non-markovian process characterisation and control on a quantum processor, Nature Communications 11, 6301 (2020).
  * Andersson  _et al._ [2019] G. Andersson, B. Suri, L. Guo, T. Aref, and P. Delsing, Non-exponential decay of a giant artificial atom, Nature Physics 15, 1123 (2019).
  * Bulla  _et al._ [2005] R. Bulla, H.-J. Lee, N.-H. Tong, and M. Vojta, Numerical renormalization group for quantum impurities in a bosonic bath, Physical Review B 71, 045122 (2005).
  * Bulla  _et al._ [2008] R. Bulla, T. A. Costi, and T. Pruschke, Numerical renormalization group method for quantum impurity systems, Reviews of Modern Physics 80, 395 (2008).
  * Feynman and Vernon Jr [2000] R. P. Feynman and F. Vernon Jr, The theory of a general quantum system interacting with a linear dissipative system, Annals of physics 281, 547 (2000).
  * Suzuki [1993] M. Suzuki,  _Quantum Monte Carlo methods in condensed matter physics_ (World scientific, 1993).
  * Egger and Mak [1994] R. Egger and C. Mak, Low-temperature dynamical simulation of spin-boson systems, Physical Review B 50, 15210 (1994).
  * Egger  _et al._ [2000] R. Egger, L. Mühlbacher, and C. Mak, Path-integral monte carlo simulations without the sign problem: Multilevel blocking approach for effective actions, Physical Review E 61, 5961 (2000).
  * Mühlbacher  _et al._ [2005] L. Mühlbacher, J. Ankerhold, and A. Komnik, Nonequilibrium dynamics of correlated electron transfer in molecular chains, Physical review letters 95, 220404 (2005).
  * Stockburger and Mak [1999] J. T. Stockburger and C. H. Mak, Stochastic liouvillian algorithm to simulate dissipative quantum dynamics with arbitrary precision, The Journal of chemical physics 110, 4983 (1999).
  * Stockburger and Grabert [2002] J. T. Stockburger and H. Grabert, Exact c-number representation of non-markovian quantum dissipation, Physical review letters 88, 170407 (2002).
  * Tanimura and Kubo [1989] Y. Tanimura and R. Kubo, Time evolution of a quantum system in contact with a nearly gaussian-markoffian noise bath, Journal of the Physical Society of Japan 58, 101 (1989).
  * Tanimura [2020] Y. Tanimura, Numerically “exact” approach to open quantum dynamics: The hierarchical equations of motion (heom), The Journal of chemical physics 153 (2020).
  * Beck  _et al._ [2000] M. H. Beck, A. Jäckle, G. A. Worth, and H.-D. Meyer, The multiconfiguration time-dependent hartree (mctdh) method: a highly efficient algorithm for propagating wavepackets, Physics reports 324, 1 (2000).
  * Wang and Thoss [2003] H. Wang and M. Thoss, Multilayer formulation of the multiconfiguration time-dependent hartree theory, The Journal of chemical physics 119, 1289 (2003).
  * Makri and Makarov [1995] N. Makri and D. E. Makarov, Tensor propagator for iterative quantum time evolution of reduced density matrices. i. theory, The Journal of chemical physics 102, 4600 (1995).
  * Makri [1998] N. Makri, Quantum dissipative dynamics: A numerically exact methodology, The Journal of Physical Chemistry A 102, 4414 (1998).
  * Gull  _et al._ [2011] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Continuous-time monte carlo methods for quantum impurity models, Reviews of Modern Physics 83, 349 (2011).
  * Xu  _et al._ [2023] M. Xu, V. Vadimov, M. Krug, J. T. Stockburger, and J. Ankerhold, A universal framework for quantum dissipation:minimally extended state space and exact time-local dynamics (2023), arXiv:2307.16790 [quant-ph] .
  * Vilkoviskiy and Abanin [2023] I. Vilkoviskiy and D. A. Abanin, A bound on approximating non-markovian dynamics by tensor networks in the time domain (2023), arXiv:2307.15592 [quant-ph] .
  * Jacobs [2014] K. Jacobs,  _Quantum measurement theory and its applications_ (Cambridge University Press, 2014).
  * Landau and Lifshitz [2013] L. D. Landau and E. M. Lifshitz,  _Quantum mechanics: non-relativistic theory_ , Vol. 3 (Elsevier, 2013).
  * Blais  _et al._ [2021] A. Blais, A. L. Grimsmo, S. Girvin, and A. Wallraff, Circuit quantum electrodynamics, Reviews of Modern Physics 93, 025005 (2021).
  * Collett and Gardiner [1984] M. Collett and C. Gardiner, Squeezing of intracavity and traveling-wave light fields produced in parametric amplification, Physical Review A 30, 1386 (1984).
  * Gardiner and Collett [1985] C. W. Gardiner and M. J. Collett, Input and output in damped quantum systems: Quantum stochastic differential equations and the master equation, Physical Review A 31, 3761 (1985).
  * Yurke and Denker [1984] B. Yurke and J. S. Denker, Quantum network theory, Phys. Rev. A 29, 1419 (1984).
  * Combes  _et al._ [2017] J. Combes, J. Kerckhoff, and M. Sarovar, The SLH framework for modeling quantum input-output networks, Advances in Physics: X 2, 784 (2017).
  * Kjaergaard  _et al._ [2020] M. Kjaergaard, M. E. Schwartz, J. Braumüller, P. Krantz, J. I.-J. Wang, S. Gustavsson, and W. D. Oliver, Superconducting qubits: Current state of play, Annual Review of Condensed Matter Physics 11, 369–395 (2020), https://doi.org/10.1146/annurev-conmatphys-031119-050605 .
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell, _et al._ , Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019).
  * Vadimov  _et al._ [2022] V. Vadimov, A. Viitanen, T. Mörstedt, T. Ala-Nissila, and M. Möttönen, Single-junction quantum-circuit refrigerator, AIP Advances 12 (2022).
  * Silveri  _et al._ [2019] M. Silveri, S. Masuda, V. Sevriuk, K. Y. Tan, M. Jenei, E. Hyyppä, F. Hassler, M. Partanen, J. Goetz, R. E. Lake, _et al._ , Broadband lamb shift in an engineered quantum system, Nature Physics 15, 533 (2019).
  * Estève  _et al._ [2018] J. Estève, M. Aprili, and J. Gabelli, Quantum dynamics of a microwave resonator strongly coupled to a tunnel junction (2018), arXiv:1807.02364 [cond-mat.mes-hall] .
  * Aiello  _et al._ [2022] G. Aiello, M. Féchant, A. Morvan, J. Basset, M. Aprili, J. Gabelli, and J. Estève, Quantum bath engineering of a high impedance microwave mode through quasiparticle tunneling, Nature Communications 13, 7146 (2022).
  * Hyyppä  _et al._ [2022] E. Hyyppä, S. Kundu, C. F. Chan, A. Gunyhó, J. Hotari, D. Janzso, K. Juliusson, O. Kiuru, J. Kotilahti, A. Landra, _et al._ , Unimon qubit, Nature Communications 13, 6895 (2022).
  * Diósi [2012] L. Diósi, Non-markovian open quantum systems: Input-output fields, memory, and monitoring, Phys. Rev. A 85, 034101 (2012).
  * Zhang  _et al._ [2013] J. Zhang, Y.-x. Liu, R.-B. Wu, K. Jacobs, and F. Nori, Non-markovian quantum input-output networks, Phys. Rev. A 87, 032117 (2013).
  * Gross  _et al._ [2022] J. A. Gross, B. Q. Baragiola, T. M. Stace, and J. Combes, Master equations and quantum trajectories for squeezed wave packets, Phys. Rev. A 105, 023721 (2022).
  * Link  _et al._ [2022] V. Link, K. Müller, R. G. Lena, K. Luoma, F. Damanet, W. T. Strunz, and A. J. Daley, Non-markovian quantum dynamics in strongly coupled multimode cavities conditioned on continuous measurement, PRX Quantum 3, 020348 (2022).
  * Propp [2022] T. B. Propp, How to describe collective decay of uncoupled modes in the input–output formalism, Journal of the Optical Society of America B 39, 3128 (2022).
  * Oliver and Welander [2013] W. D. Oliver and P. B. Welander, Materials in superconducting quantum bits, MRS bulletin 38, 816 (2013).
  * Martinis  _et al._ [2005] J. M. Martinis, K. B. Cooper, R. McDermott, M. Steffen, M. Ansmann, K. Osborn, K. Cicak, S. Oh, D. P. Pappas, R. W. Simmonds, _et al._ , Decoherence in josephson qubits from dielectric loss, Physical review letters 95, 210503 (2005).
  * Gordon  _et al._ [2022] R. Gordon, C. Murray, C. Kurter, M. Sandberg, S. Hall, K. Balakrishnan, R. Shelby, B. Wacaser, A. Stabile, J. Sleight, _et al._ , Environmental radiation impact on lifetimes and quasiparticle tunneling rates of fixed-frequency transmon qubits, Applied Physics Letters 120 (2022).
  * Wang  _et al._ [2014] C. Wang, Y. Y. Gao, I. M. Pop, U. Vool, C. Axline, T. Brecht, R. W. Heeres, L. Frunzio, M. H. Devoret, G. Catelani, _et al._ , Measurement and control of quasiparticle dynamics in a superconducting qubit, Nature communications 5, 5836 (2014).
  * Barends  _et al._ [2011] R. Barends, J. Wenner, M. Lenander, Y. Chen, R. C. Bialczak, J. Kelly, E. Lucero, P. O’Malley, M. Mariantoni, D. Sank, _et al._ , Minimizing quasiparticle generation from stray infrared light in superconducting quantum circuits, Applied Physics Letters 99 (2011).
  * Catelani and Pekola [2022] G. Catelani and J. P. Pekola, Using materials for quasiparticle engineering, Materials for Quantum Technology 2, 013001 (2022).
  * Rainis and Loss [2012] D. Rainis and D. Loss, Majorana qubit decoherence by quasiparticle poisoning, Physical Review B 85, 174533 (2012).
  * Xu  _et al._ [2022] M. Xu, Y. Yan, Q. Shi, J. Ankerhold, and J. Stockburger, Taming quantum noise for efficient low temperature simulations of open quantum systems, Physical Review Letters 129, 230601 (2022).
  * Nigg  _et al._ [2012] S. E. Nigg, H. Paik, B. Vlastakis, G. Kirchmair, S. Shankar, L. Frunzio, M. H. Devoret, R. J. Schoelkopf, and S. M. Girvin, Black-box superconducting circuit quantization, Phys. Rev. Lett. 108, 240502 (2012).
  * Solgun  _et al._ [2014] F. Solgun, D. W. Abraham, and D. P. DiVincenzo, Blackbox quantization of superconducting circuits using exact impedance synthesis, Phys. Rev. B 90, 134504 (2014).
  * Schrieffer and Wolff [1966] J. R. Schrieffer and P. A. Wolff, Relation between the anderson and kondo hamiltonians, Physical Review 149, 491 (1966).
  * Bravyi  _et al._ [2011] S. Bravyi, D. P. DiVincenzo, and D. Loss, Schrieffer–wolff transformation for quantum many-body systems, Annals of physics 326, 2793 (2011).
  * Caldeira and Leggett [1983] A. O. Caldeira and A. J. Leggett, Path integral approach to quantum brownian motion, Physica A: Statistical mechanics and its Applications 121, 587 (1983).
  * Schmid [1982] A. Schmid, On a quasiclassical langevin equation, Journal of Low Temperature Physics 49, 609 (1982).
  * Benguria and Kac [1981] R. Benguria and M. Kac, Quantum langevin equation, Physical review letters 46, 1 (1981).
  * Ford and Kac [1987] G. Ford and M. Kac, On the quantum langevin equation, Journal of statistical physics 46, 803 (1987).
  * Ford  _et al._ [1988] G. W. Ford, J. T. Lewis, and R. O’Connell, Quantum langevin equation, Physical Review A 37, 4419 (1988).
  * Vool and Devoret [2017] U. Vool and M. Devoret, Introduction to quantum electromagnetic circuits, International Journal of Circuit Theory and Applications 45, 897 (2017).
  * Nakatsukasa  _et al._ [2018] Y. Nakatsukasa, O. Sète, and L. N. Trefethen, The aaa algorithm for rational approximation, SIAM Journal on Scientific Computing 40, A1494 (2018).
  * Hubbard [1959] J. Hubbard, Calculation of partition functions, Physical Review Letters 3, 77 (1959).
  * Stratonovich [1957] R. Stratonovich, On a method of calculating quantum distribution functions, in  _Soviet Physics Doklady_ , Vol. 2 (1957) p. 416.
  * Shi  _et al._ [2018] Q. Shi, Y. Xu, Y. Yan, and M. Xu, Efficient propagation of the hierarchical equations of motion using the matrix product state method, The Journal of Chemical Physics 148, 174102 (2018).
  * Borrelli [2019] R. Borrelli, Density matrix dynamics in twin-formulation: An efficient methodology based on tensor-train representation of reduced equations of motion, The Journal of Chemical Physics 150, 234102 (2019).
  * Koch  _et al._ [2007] J. Koch, T. M. Yu, J. Gambetta, A. A. Houck, D. I. Schuster, J. Majer, A. Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Charge-insensitive qubit design derived from the cooper pair box, Phys. Rev. A 76, 042319 (2007).
  * Gambetta  _et al._ [2007] J. Gambetta, W. Braff, A. Wallraff, S. Girvin, and R. Schoelkopf, Protocols for optimal readout of qubits using a continuous quantum nondemolition measurement, Physical Review A 76, 012325 (2007).
  * Gunyhó  _et al._ [2023] A. M. Gunyhó, S. Kundu, J. Ma, W. Liu, S. Niemelä, G. Catto, V. Vadimov, V. Vesterinen, P. Singh, Q. Chen, and M. Möttönen, Single-shot readout of a superconducting qubit using a thermal detector (2023), arXiv:2303.03668 [quant-ph] .
  * Braun  _et al._ [2001] D. Braun, F. Haake, and W. T. Strunz, Universality of decoherence, Physical review letters 86, 2913 (2001).
  * Tuorila  _et al._ [2019] J. Tuorila, J. Stockburger, T. Ala-Nissila, J. Ankerhold, and M. Möttönen, System-environment correlations in qubit initialization and control, Phys. Rev. Res. 1, 013004 (2019).