Engineering dissipation with resistive elements in circuit quantum electrodynamics

The importance of dissipation engineering ranges from universal quantum computation to non-equilibrium quantum thermodynamics. In recent years, more and more theoretical and experimental studies have shown the relevance of this topic for circuit quantum electrodynamics, one of the major platforms in the race for a quantum computer. This article discusses how to simulate thermal baths by inserting resistive elements in networks of superconducting qubits. Apart from pedagogically reviewing the phenomenological and microscopic models of a resistor as thermal bath with Johnson-Nyquist noise, the paper introduces some new results in the weak coupling limit, showing that the most common examples of open quantum systems can be simulated through capacitively coupled superconducting qubits and resistors. The aim of the manuscript, written with a broad audience in mind, is to be both an instructive tutorial about how to derive and characterize the Hamiltonian of general dissipative superconducting circuits with capacitive coupling, and a review of the most relevant and topical theoretical and experimental works focused on resistive elements and dissipation engineering.


I. INTRODUCTION
Circuit quantum electrodynamics (circuit QED) is nowadays at the forefront of experimental quantum physics, offering an experimental platform for topical research fields such as quantum information processing [1][2][3], quantum simulation [4,5], and quantum thermodynamics [6][7][8]. By employing few circuit elements, namely capacitors, inductors, Josephson junctions and transmission line resonators, the experimentalist using this platform is able to initialize, measure and control the dynamics of one or more qubits, interacting or not with an external field [9,10].
The quantization of the above-mentioned circuit components is quite straightforward and well-understood, but including resistors in the formalism of circuit QED is not a trivial task. This is due to the fact that, while capacitors, inductors, junctions and finite transmission line resonators are non-dissipative elements modeleded by Hamiltonians, resistors induce energy loss in the circuits. This problem can be solved by adapting the phenomenological Caldeira-Leggett model of dissipation employing an infinite number of bosonic modes to describe resistive elements in circuit QED [11]. When quantized through this method, a resistor behaves as a thermal bath [12] with an Ohmic spectral density, reproducing the well-known Johnson-Nyquist (or thermal) noise [13,14]. This idea is fundamental for the study of engineered dissipation in circuit QED. Our aim is to provide a broad introduction to this method, together with the necessary tools for its application in general scenarios, which may guide the reader through the study of networks of qubits and resistors in circuit QED.
We start by presenting a pedagogical, step-by-step derivation of the Hamiltonian of the Caldeira-Leggett model and of the spectral density of the voltage fluctuations across a resistor. We do not assume for our presentation any previous knowledge about the topic, and close familiarity with quantum electromagnetic circuits is not strictly necessary throughout the paper, which we have written with a broad audience of theoretical and experimental quantum physicists in mind.
Having understood the principles of the Caldeira-Leggett model and how to apply them, it is important to clarify the physical intuition of how and why this phenomenological model works at all. Indeed, the phenomenological models for resistors do not depend on the specific materials used in their fabrication or on specific microscopic models; rather, they are constrained only by the general principles of thermodynamics, as first pointed out by Nyquist already in 1928 [14]. In order to understand the microscopic origin of fluctuations in resistors, we introduce the Landauer-Büttiker method for the study of scattering processes in a conductor as presented in the original papers [15][16][17][18] and in a well-known review on the topic [19]. We show that a description of the resistor based on this method provides a microscopic model that correctly reproduces the Johnson-Nyquist spectrum.
After introducing the phenomenological and microscopic models of resistors in circuit QED, we focus on how to employ the latter as thermal baths to simulate open quantum systems [12]. In order to do so, it is crucial to understand how the coupling between them and any circuit element (such as a transmon qubit [20]) works, and how to tune it. We present the derivation of the general Hamiltonian of a circuit with capacitively coupled qubits and resistors, and discuss the master equations driving the open qubit evolution. Obtaining the final circuit Hamiltonian is not an easy task, and we rely on a recent seminal paper that has paved the way for the treatment of any complex circuit network of coupled systems of this kind [21]. Furthermore, we introduce some new results about the derivation in the weak coupling limit, leading to Markovian master equations. In particular, we show how the reservoir engineering allows for the simulation of the open system dynamics of one or two qubits coupled to common or separate baths. The same procedure trivially applies for the study of harmonic oscillators (LC circuits in the present framework) immersed in thermal reservoirs. Moreover, we briefly show how one can introduce a LC circuit acting as a filter, in order to shape the spectral density of the dissipation on the qubits.
To conclude the discussion, we review some of the major experimental methods to insert dissipative elements into superconducting circuits, with an eye to reservoir engineering. Fabricating resistive components in circuit QED is relatively easy, since they usually consist of simple normal-metal elements in a circuit made of superconducting materials [22]. The dissipation induced by these elements is nowadays well-understood, as well as their thermodynamic properties, including the role played by the electron-phonon interaction [23,24]. Resistors as independent elements of quantum circuits have been recently employed in disparate fields, with a particular focus on quantum thermodynamics [6], and experiments relying on them have been presented in the past few years [22,[25][26][27]. Besides resistors, transmission lines and lossy cavities have been used in circuit QED either as 3D or coplanar microwave components, and they can be modeled by similar techniques [28][29][30]. In addition, we outline a few of the recent inroads into modern approaches on the active control of dissipation. While reducing decoherence is nowadays a major effort on the road to quantum computing, here we show, in contradistinction, that dissipation can be engineered to achieve stabilization of certain states in superconducting three-level systems and in resonators. A student in microwave engineering or rf circuits discovers quite early one little puzzle of this type: in a series RLC circuit, the quality factor is (1/R) L/C, in agreement with our naive intuition that an increased resistance leads to a worse quality factor for the oscillator. But for a parallel RLC circuit, the quality factor is R C/L, increasing as R gets larger -the reason being, of course, that for large R the current will flow less through the shunt resistor and more into the ideal LC oscillator itself. The end of our paper highlights such situations where dissipation is beneficial, reviewing recent microwave engineering methods based on filtering, tunability, and microwave-assisted couplings to lossy transmission lines. Having at our disposal a controllable dissipative reservoir acting on superconducting qubits would be of great interest for both quantum engineers and quantum physicists [31]. The goal is to use dissipation as a resource for creating, in conjunction with coherent driving [32], novel quantum states of matter and novel forms of quantum control. The benefit is stabilization and robustness, which comes for free. For digital quantum computing, the usefulness of dissipation is mostly related to the realization of precise and fast initial states, and to some extend, if a modularbased approach is to be pursued, that of quantum memories [33]. However, alternative directions for quantum computing based on control of dissipation can be envisioned, where the steady-state negociated by the environment and coherent control encodes the result of a computation [34]. These are forms of analog quantum processing of information [5,35]. For simulating manybody systems, superconducting circuits provide time-resolved measurement data for each lattice site, allowing the extraction of higher-order correlations, which enable the investigation of thermalization processes and dynamics of phase transitions.
The paper is organized as follows: in Sec. II we introduce the phenomenological and microscopic model of resistors in circuit QED. After presenting the concept of impedance and its extension for physical signals, we discuss the Caldeira-Leggett model for resistors and calculate the spectral density, showing its relation with the Johnson-Nyquist noise. Then, we introduce the Landauer-Büttiker method to sketch a microscopic model of the resistor reproducing the Johnson-Nyquist noise. After reviewing the capacitive coupling in Sec. III, in Sec. IV we address how to employ it to connect resistors with other circuit elements, focusing in particular on the cases of one or two qubits, and discussing how to add a spectral filter. The weak coupling limit plays here a major role, while the general derivation for strong couplings is discussed in Appendix D. Sec. V addresses the main available control methods of dissipation, and reviews some experiments that make use of them. Finally, in Sec. VI we briefly draw some concluding remarks and perspectives.
In summary, our purpose is to establish a toolbox of theoretical models for dissipation engineering, with a focus on the comprehension of the Caldeira-Leggett model and its application in networks of coupled qubits.

II. A HAMILTONIAN FORMALISM FOR RESISTORS IN CIRCUIT QED
In this section we discuss how to describe a resistor in the time-independent Hamiltonian formalism for the quantization of quantum circuits, which the reader can find in Appendix A with an application to the LC circuit. We introduce basic notions of circuit quantization and for further details we refer the reader to many excellent reviews about this topic [1,2,4,5,31,[36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51]. After introducing the definition of impedance (or equivalently admittance) and its generalization for physical signals, we focus on the Caldeira-Leggett model of a resistor presented in one of the very first lectures about quantized electromagnetic circuits [11], and we follow the lines of this work to pedagogically review the model derivation. Next, we calculate the spectral density of the voltage fluctuations across the resistor, and we show that it is equal to the one of a microscopic physical model based on the Landauer-Büttiker method.

A. Impedance and admittance
In the language of electrical engineering, the impedance of a circuit element is defined when dealing with AC signals. Already at this initial step in the journey, the astute and diligent scholar will notice the discrepancy between the conventions used by physicists and engineers when it comes to defining Fourier transforms, phasors, and positive/negative frequencies. Throughout this paper, we adopt the definition employed in Ref. [52], and we give below a number of considerations on how to relate this to the standard notations used in microwave engineering. Given a time-dependent operator or complex function F (t), its Fourier transform is (1) with inverse If F is an operator, the direct and inverse Fourier transform of the adjoint are, by the definition above, and similarly (with complex conjugation instead of the dagger symbol) for a complex function. The following useful identities regarding a change of sign in frequency, hold true: One notices that these definitions are generally at loggerheads with the ones used in elementary signal analysis. We provide additional clarifications about these conventions in Appendix B.
We define the complex impedance Z(ω) of a branch of a circuit by The same definition can be given in terms of the phasors V e −iφ V and Since V (t) and I(t) are real functions, using the properties of the Fourier transform above we also have As mentioned already above, this relation allows us to translate our results with ease into the engineering convention if needed. The impedance can be interpreted as the opposition to the passage of alternate current through the branch when a certain voltage is applied. Equivalently, we can use the admittance Y as a measure of how easily the current runs under the driving of an applied voltage: We can easily compute the impedance of all the passive circuit elements. Indeed, since Ohm's law for a resistor with resistance R reads V = RI, we have being characterized by a constant value all over the range of frequencies. When dealing with a capacitor with capacitance C, we employ the relation I = CdV /dt and obtain: Finally, the constitutive equation for an inductor with inductance L reads V = LdI/dt, therefore From Eq. (4) we observe that the impedance of a series of circuit elements is given by the sum of the impedances of the elements. On the contrary, the inverse of the impedance of some circuit elements in parallel is given by the sum of the inverses of the impedances of the elements. Hence, the impedance of a parallel LC circuit reads which, after few steps, can be rewritten as where ω LC = 1/ √ LC is the resonant frequency of the circuit. When discussing quantum electromagnetic circuits, we do not want to restrict ourselves to the AC circuit regime, but, on the contrary, we are interested in dealing with generic signals. Due to the linearity of the newtork (at least in the dissipative part), such signals, namely V (t) and I(t), can be written as a sum of the AC signals over all the frequencies, i.e. as a Fourier transform: and equivalently for I(t). In Eq. (12) we consider negative frequencies also, since V (t) and I(t) are real, physical signals, and the negative-frequency contribution deletes the imaginary part. Since V (t) and I(t) are indeed physical signals, they properly vanish at t = ±∞, therefore their Fourier transforms are well-defined. Starting from Eq. (4), we write the relation between voltage and current by making use of a linear response function: where Z(t) is the Fourier transform of the complex impedance Z(ω). Some divergence issue may arise in the last step; it turns out that, to avoid any problem, we need to employ a generalized impedance function, defined as: Using Eq. (14), the relation between current and voltage reads In Appendix C we carefully discuss the derivation of Eq. (13) and the meaning of Eq. (14). Here, we just point out that such definition is necessary to introduce the linear response function of the LC circuit, given that Eq. (11) does not admit Fourier transform because of the poles at ω = ±ω LC . In this case, employing Eqs. (14) and (15) corresponds to shifting the poles in the lower half plane. Moreover, notice that the negative exponential in Eq. (14) does not cause any divergence issue, since for t < 0 causality requires Z(t) = 0 [53]. Finally, the definition of a linear response function as in Eq. (13) implies a couple of relations between the real and imaginary part of Z(ω), known as Kramers-Kronig relations [53]: where P denotes the principal value integral.

B. Caldeira-Leggett phenomenological model of a resistor
A phenomenological model for a resistor can be obtained by adapting the seminal concepts introduced by Caldeira and Leggett in 1983 to describe dissipation in quantum mechanics [54]. The key idea is to describe the environment by a collection of harmonic oscillators -which, in the case of electrical circuits, can be taken as a collection of LC oscillators. Such a model was formulated with clarity and in modern notations and electrical-engineering terminology by Michel Devoret in the lecture notes of Les Houches summer school in 1995 [11], and further refined in Ref. [45]. The model is based on a method of network synthesis used to represent, under certain assumptions, a generic purely imaginary impedance as a collection of LC circuits. Different ways of realizing such synthesis are available, and here we focus on the ones called first and second Foster's form [55,56], extending the method to dissipative impedances, characterized by a real part also. We mention here that a renewed interest in the description of a generic impedance (or admittance) in circuit QED has emerged in the last twenty years, and has led to the development of novel theoretical methods. For instance, a couple of seminal papers by Burkard  times of superconducting flux and charge qubits in the presence of an environment, where the latter is described as a collection of harmonic oscillators in a Caldeira-Leggett fashion, considering weak coupling and possible leakages to higher qubit levels [57,58]. Black box quantization, presented in Ref. [59], is a method of obtaining the effective low energy Hamiltonian of one (or more) Josephson junctions in parallel with a generic environment, by describing the latter through the Foster's first form.
A refined version of this method, relying on Brune's synthesis instead of Foster's one, has been presented in Ref. [60] and later extended to multiport impedances [61]. The well-established input-output theory to compute dissipation [62] has been broadly employed since the very early stages of circuit QED, with a special attention on the analysis of open transmission line resonators with an infinite number of modes [63]. Quantization of generic networks in circuit quantum electrodynamics has been extensively analysed by Parra-Rodriguez et al. [64], considering both reciprocal [21] and non-reciprocal [65] devices. Finally, a recent paper proposes a new method to derive the master equation of a system of superconducting qubits coupled to a low admittance, corresponding to weak coupling to the environment, by extracting only a few relevant degrees of freedom of the circuit and treating the rest in a Lindbladian way [66].

Building the model
Let us start with the Foster's first form to describe an impedance, say Z(ω), by employing capacitors and inductors only. We consider a network composed of a collection of parallel LC circuits in series, as in Fig. 1. If the number of parallel LC circuits is finite, the overall impedance of the network must be purely imaginary. On the contrary, if we consider an infinite number of them, the real part of the impedance appears. Indeed, by making use of Eq. (11) and of the law on the impedance of a series of elements, the generalized impedance of the network in Fig. 1, which we name as Z ∞ (ω), reads where C j and ω j are respectively the capacitance and the resonant frequency of the jth LC circuit. Now, we carefully choose the values of each C j and ω j , such that: , where Z(ω) is the impedance whose action we want to mimic and ∆ω is very small. The corresponding inductance L j of each LC circuit is readily obtained as: If ∆ω is small enough to be considered as infinitesimal, the infinite summation in Eq. (18) can be seen as an integral over the (now continuous) variable ω j , and after inserting the chosen values we obtain: We can now apply the Plemelj-Sokhotski theorem for the real line[67] to the above integral, splitting the latter into a real and an imaginary part. Notice that the presence of a real part is due to the infinite number of LC circuits characterized by dense values of the resonant frequency. Indeed, a non-zero real part, which will be proportional to the resistance R, indicates the presence of dissipation in the circuit. From a physical perspective, its origin can be traced back to the infinite number of LC circuits: intuitively, an input signal can be sent through the chain of LC circuits and, since it will never reach the "end" of the chain, it will never "bounce back" and return to the start. This means that input excitations can travel "forever" through the LC chain without any recurrence time, that is, energy at the input is dissipated. This heuristic picture shows how the infinite-dimensional Foster's decomposition is able to simulate physical dissipation by employing non-dissipative elements only, provided we use an infinite number of them. Physically, the poles of the impedance at the resonance frequencies ω j correspond to the resonant excitation of the L j C j oscillator when the circuit is open-ended. By studying separately the case ω > 0 and ω < 0, we obtain that the real and imaginary parts of Z ∞ (ω) are given by: We have been successful in reproducing the real part of the impedance Z(ω) by employing an infinite network of passive capacitors and inductors only. The relation between real and imaginary part given by Eq. (21) expresses the Kramers-Kronig relations in Eq. (16), stressing that voltage and current depend on each other through a linear response function. A suitable expression for the impedance of a resistor fulfilling Eq. (21) is given by: where ω C is a suitable cutoff frequency, and the above equation is valid in the range ω ω C . For a discussion about the cutoff frequency in circuit quantum electrodynamics, we refer the reader to Refs. [68,69]. Eq. (22) can be derived through the classical Lorentz-Drude dispersion model [53], and is usually referred to as the Ohmic spectrum of a resistor. The reader can verify that it fulfils the condition of Eq. (21) by substituting the value of Z(ω) in Eq. (21) with Z eff R (ω) in Eq. (22), and solving the integrals in Eq. (21b).
In the derivation of Eq. (21) we have discussed the cases ω > 0 and ω < 0, but not the value ω = 0. Indeed, from Eq. (17) we observe that, if ω = 0, then the impedance vanishes. This is a requirement of the Foster's first form, and it can be interpreted as the fact that, for a DC signal, the capacitors in Fig. 1 do not allow for the passage of current. Therefore, the latter flows through each inductor, which does not oppose any resistance, i.e. the impedance is zero. This is a pathological case appearing only for the singularity ω = 0, therefore, if we deal with a physical signal whose contribution given by frequency ω = 0 (i.e. DC current) is negligible, we should not bother about it. However, if we want to treat the response of the system to DC current, we need to add an initial capacitor or inductor, depending on the behavior we want to mimic [11,45].

Quantization of the resistor
Let us now address the quantization of the resistor. We consider here the open circuit of Fig. 1, i.e. the two terminals of the resistor are grounded. In this case, we can parametrize the system as shown in Fig. 1, where we term the flux associated to the voltage difference between the jth LC circuit as φ j . Using this convention, we can see the LC circuits as a collection of decoupled harmonic oscillators, and we quantize the system through the standard quantization rules of circuit QED discussed in Appendix A. The resistor Hamiltonian is readily written: where Q J = C jφj , a † j and a j are the standard creation and annihilation operators [70] and the values of C j , L j , ω j are given by Eq. (18) and Eq. (19), while the resistor impedance thereof can be found in Eq. (22). Note that, using the standard convention of quantum optics, we have neglected an infinite constant term in the right-hand side of Eq. (23).
To conclude the exposition of the Caldeira-Leggett model for a resistor, note that all the previous results can be obtained for an admittance instead of an impedance, by describing the equivalent network using the Foster's second form [45]. All the equations of the present section still hold, switching voltage with current and capacitance with inductance. In this case, the poles of the admittance at ω j correspond to the resonant excitation of the corresponding L j C j oscillator, as the input of the line is shortcut. The only difference consists in the fact that Foster's second form requires a vanishing admittance at ω = 0, that is to say, infinite impedance.

Thermal noise
We are now able to calculate the spectrum of the voltage fluctuations over the resistor, if the latter is in a thermal state at temperature T = 1/βk B . The thermal state is defined as: where ρ B (β, ω j ) describes the thermal state of each mode: where |n j are the eigenstates of the number operator of the mode j, n j = a † j a j , and n(ω j ) = n j ρ B (β,ωj ) = 1/(e ωj β − 1) is the Bose-Einstein distribution. Fore more details about the thermal state we refer the reader to Ref. [70].
Let us now calculate the correlation function of the voltage across the jth LC circuit. Such function is defined as where the evolution of any operator is computed according to the Hamiltonian Eq. (23), and we are exploiting its time-translation invariance to define the correlation function as depending only on t. From standard Hamilton's equations, we observe thaṫ φ j (t) =φ j (0) cos ω j t − ω j φ j (0) sin ω j t. We recall that φ j (0) andφ j (0) are related to the creation and annihilation operators through: Therefore, where we have used Eq. (25) and the commutator [a j , a † j ] = 1. Analogously, we have Finally, we calculate the correlation function: The operator describing the voltage difference across the resistor is nothing but the sum of the voltages over each LC circuit, V = jφ j . Since all the harmonic oscillators in Eq. (23) are uncorrelated, after substituting the value of C j given by Eq. (18), the correlation function of the voltage operator reads where we have replaced the summation over all the modes with an integral, and we have introduced the notation N (ω j ) = (n(ω j ) + 1)Θ(ω j ) − n(−ω j )Θ(−ω j ). Here n(ω) = (e β ω − 1) −1 and Θ(ω) is the Heaviside theta function. From the result Eq. (30) above it is clear that positive frequencies correspond to absorption of energy and negative frequencies to emission. This can be seen by considering T = 0 (vacuum state). In this case only the positive-frequency part remains (indicating that the system is still able to absorb energy), while the negative-energy part is zero (indicating that it is not possible to emit the energy of vacuum). This is reflected in the definition of the quantity N (ω), which is N (ω) = n(ω) + 1 for absorptive processes (ω > 0) and N (ω) = −n(−ω) for emission (ω < 0). The spectral density of an operator O is defined as the Fourier transform of its correlation function [52]: At this point it is useful to mention the connection with the Wiener-Khinchin theorem in signal analysis. Given a random process X(t), the power spectral density (PSD) is defined as being the windowed Fourier transform and ..., ... the ensemble average -i.e. expectation value with respect to the probability distribution function (PDF) of the process. The quantity |X t W (ω)| 2 /t W is referred to as periodogram or sample spectrum. Suppose we record the signal with an oscilloscope in a time-window t W . We can take the Fourier tranform of the data and calculate |X t W (ω)| 2 /t W , which is an unbiased estimator of the true power spectral density Eq. (33) obtained by averaging over an ensemble of such measurements. For wide-sense stationary processes, the autocorrelation function depends only on the lag between the two times: Then, the Wiener-Khinchin theorem says that Hence, the power spectral density is the Fourier transform of the auto-correlation function, and following the reasoning above Eq. (32) can be interpreted as a mathematical result instead of a definition. For a detailed proof of this result we refer to Ref. [71] Finally, the spectral density of the fluctuations of the voltage operator across a resistor, reads Eq. (37) represents the quantum limit of the well-known Johnson-Nyquist noise induced by thermal agitation [13,14]; we will see throughout the paper how this expression can be interpreted as the standard fluctuation-dissipation theorem for electrical circuits. For low frequencies or high temperatures we find the classical result S V V (ω) = 2k B T Re[Z(ω)], with a factor 2 instead of a factor 4 because of the employed definition of spectral density including negative frequencies also (see the discussion below). This expression is referred to as the double-sided spectral density, since ω can take both positive and negative values.
The zero-voltage noise due to the thermal agitation of electric charges in a conductor was first observed experimentally by John B. Johnson [13], and then described theoretically by Harry Nyquist [14] during 1926-1928, while both scientists were working at Bell Labs. Nyquist's derivation is a beautiful demonstration of the heuristic power of fundamental concepts from thermodynamics and statistical mechanics. In his seminal paper from 1928, he considered two resistors with resistance R at the same temperature T , connected by conductive wires. Since the resistors are in thermal equilibrium and have the same resistance, the electrical power due to thermal agitation flowing towards a resistor or the other must be same. This must hold true not only for the total power, but also the power exchanged within any frequency bandwidth. If this were not the case, by inserting a filter or a resonant circuit in-between the resistors one obtains immediately a contradiction with the second law of thermodynamics, as such a circuit would allow extracting work from a single reservoir at temperature T . It then follows that the voltage across the resistor must be a function of frequency, resistance and temperature only. Next, to derive the final formula, Nyquist analysed a circuit where two resistors are connected by a transmission line of length l. Let us first calculate the power transmitted from one resistor to the other. If V is the voltage generated by thermal noise corresponding say to the left resistor, then this creates a voltage V R/(R + R) = V /2 on the right resistor, and the corresponding power transmitted through the line to be dissipated in that resistor is P = V 2 /(4R). Next, we recall Parseval's theorem, showing that the energy content of the signal is the same irrespective to whether it is expressed in time-domain or in frequency domain.
Further, let us define the time-average power of the signal V (t) by We can now use Parseval's theorem in the relation above and, further, we can move the limit inside the integral by taking the ensemble average. But since we assume that the process is ergodic, the ensemble average should give the same result as the temporal average. Thus we obtain the mean squared voltage where S V V (ω) = lim t W →∞ |V (ω)| 2 /t W is the voltage power spectral density, with ..., ... denoting ensemble average. This means that, in the frequency space, the power is distributed over frequencies and therefore the quantity represents the average voltage power in an interval dω at a frequency ω. This can be understood also as the result of filtering the voltage V (t) with a bandpass filter centered at ω with bandwidth dω.
The same result can be obtained by applying the Wiener-Khinchin theorem. From Eq. (36) we have R XX (τ ) = ∞ −∞ dω/(2π)S XX (ω)e −iω and for τ = 0 the autocorrelation function becomes simply the power, R XX (t) = X 2 (t) . Identifying X = V we obtain We can also see from the above expression that, when restricting the integrals to only positive frequencies, where the quantity is called the symmetrized single-sided spectral density, and has the meaning of voltage power per bandwidth at positive frequencies. Typical laboratory instruments such as spectrum analysers measure directly S V V (ω), since they use a local oscillator with positive frequencies to mix, downconvert, and power-detect the signal. Returning now to our circuit, let us consider the positive-frequency case ω > 0 as well. Since V (t) is real, as we have seen in Next, we can write the average power dissipated in the resistor in a frequency interval dω as 1 2π Therefore, the energy dE in this bandwidth can be obtained from the power above, transmitted through the transmission line during the time l/v (the time needed for photons to traverse the line) by the two sources (the left and right resistors) On the other hand, the energy dE can be obtained by the following reasoning. After the equilibrium has been established, the transmission line is suddenly disconnected and shorted at both ends. This ensures that the energy is stored in modes satisfying l = nλ/2, where λ = 2πv/ω is the wavelength, v is the speed of light in the line, and ω/(2π) is the frequency. Therefore the number of modes in the frequency interval dω/(2π) is dn = ldω/(πv). From the equipartition theorem, we have to ascribe an energy k B T per mode (each mode has two degrees of freedom, magnetic and electric). As a result, which is valid in the large-temperature limit k B T ω. Nyquist also showed how to generalize this formula for any complex impedance, which amounts in our model to the replacement of R by Re[Z(ω)]. To cover also the large-frequency limit, he proposed to use Planck's distribution instead of the equipartition theorem from classical thermodynamics, which amounts to replacing k B T with ω/[exp( βω) − 1] in Eq. (46). Moreover, we should add a factor of 1/2 coming from vacuum fluctuations, resulting finally in Note that we still recover Eq. (47) in the limit βω → 0 since coth(β ω/2) → 2.
Comparing to our previous results, we find that this is the symmetrized single-sided spectral density For ease of comparison with the full quantum results derived above, note that n(ω) = 1 2 [coth(β ω/2) − 1] and n(ω) being the number of photons. We also have the relation n(−ω) = −n(ω)− 1. Also we note that the difference of these spectra is Finally, we would like to highlight a similarity between the phenomenological derivation of the Caldeira-Leggett model and Nyquist's derivation: the first Foster's form of a resistor is nothing but the normal mode decomposition of the electrical oscillations in the resistive circuit, as given in Eq. (23). Postulating that the resistor is in a bosonic thermal state means that, to each mode with frequency ω j , we assign an energy contribution of ω j n(ω j ) given by the Bose-Einstein distribution (see Eq. (25)). This is the quantum generalization of the standard equipartition principle. Analogously, Nyquist took into account all the modes of short-ended transmission line and assigned to each of them an energy contribution as in Eq. (46), according to the classical energy equipartition principle. The generalization of this to the quantum case requires employing the Bose-Einstein distribution, which again leads to the factor ω j n(ω j ) as above. We therefore see how both phenomenological derivations rely on similar arguments about mode decomposition and quantum energy distribution.
Both the Nyquist and the Caldeira-Leggett approaches are material-independent, that is to say, they do not consider the internal microscopic structure of the resistor, while being based on purely thermodynamic principles. The discussion of a microscopic model for the resistor and corresponding voltage fluctuations are the topic of the following section.

C. Microscopic description of quantum fluctuations through a resistor
In this section we employ the Landauer-Büttiker method to describe the quantum fluctuations across a resistive element in a mesoscopic circuit. For simplicity, we apply the method to a simple case, namely a DC voltage across a NIN junction, and then provide some more specialized references for the NIS junction or AC voltage.
Our aim is to write down the operator describing the current passage through a mesoscopic NIN junction and to compute its spectral density. We follow the lines of the well-known review on the topic by Blanter and Büttiker [19]. We describe the passage of current as a scattering process across the insulator barrier, schematized as in Fig. 2: we assume that the electrons can run freely in the left (L) and right (R) leads, so that the total wavefunction (far from the barrier) can be written as a sum of free particles wavefunctions. We do not need to describe in details the scattering process inside and at the borders of the insulating barrier, while, for our purposes, it is sufficient to relate the output modes to the input modes through the scattering matrix. Specifically, using a second-quantization formalism, we call α L1 , . . . , α LN (α R1 , . . . , α RN ) the input modes in the left (right) lead, while β L1 , . . . , β LN (β R1 , . . . , β RN ) are the output modes, as shown in Fig. 2. α Ln,Rm and β Ln,Rm are fermionic annihilation operators, and for simplicity we have taken an equal number of modes ("channels") in the left and in the right lead. Let us insert them in two vectors as α = (α L1 , . . . , α LN , α R1 , . . . , α RN ) T , β = (β L1 , . . . , β LN , β R1 , . . . , β RN ) T . Then, by looking at the characteristics of the scattering process, we can write the output modes as functions of the input modes: where S is the scattering matrix. We can write the latter as where r, r are block matrices describing the reflection coefficients in respectively the left and right lead, while t and t describe the transmission ones. The total wavefunction for the left lead for example can be written as: FIG. 2: Scattering channels in a NIN junction, where the left and the right leads are coupled to reservoirs with temperatures respectively T L and T R and chemical potentials respectively µ L and µ R .
where χ Ln (x, y) are transverse wavefunctions, k Ln (E) = 2m(E − E Ln )/ is the free-particle momentum of the nth channel, and v Ln (E) = k Ln (E)/m is the carrier velocity and m the carrier mass. By making use of Eq. (52), we can write the current operator in the left lead as: After some tedious algebra, and under the assumption that the carrier velocity varies slowly as a function of energy, we can rewrite Eq. (53) as [19]: with where Starting from Eq. (54), we can calculate all the non-equilibrium features of the current across the mesoscopic NIN junction. In particular, let us assume that the left lead is connected to a reservoir with temperature T L and chemical potential µ L as shown in Fig. 2 (and correspondingly for the right lead). Then, we can compute the spectral density of the operator I L (t) in Eq. (54), defined as: For the zero frequency case, current conservation immediately implies S I R I R (0) = S I L I L (0), therefore without losing generality we can focus on the latter. After some long but trivial manipulations of Eq. (54) in the definition of Eq. (56), it can be shown [19] that, for zero frequency and zero applied voltage, the spectral density has the following form: where t n (E) are the eigenvalues of the matrix t † t, while f L (E) is the Fermi-Dirac distribution in the left lead at energy E.
In the common case where t n (E) varies slowly as a function of energy, we can replace it with its value at the Fermi energy, t n (E F ). Integrating the Fermi-Dirac distributions in the equilibrium scenario (T L = T R ), we obtain: where G = e 2 π N n=1 t n (E F ) is the conductance according to Landauer's formula [72]. Now we can compare Eq. (58) with Eq. (37) (the latter in the low-frequency regime): they both display the Johnson-Nyquist thermal noise at equilibrium (the former for current while the latter for voltage), although they have been obtained through completely different approaches. The spectrum of Eq. (37) has been derived following a phenomenological model, based on the impedance behavior of the resistor in the lumped element model of the circuit, and makes use of an effective temperature of the resistor. On the contrary, Eq. (58) has been derived through a microscopic physical model describing the scattering of charge carriers in a NIN junction, and the temperature appearing in the noise formula is the temperature of the left and right leads of the junction. Similar results can be found in the case of non-zero frequency spectral density [17] or NIS junctions [18,73], leading to Eq. (37). We can therefore conclude that the microscopic model discussed in this section provides us with a good intuitive physical motivation of the Caldeira-Leggett model used throughout the paper.
In general terms, the Caldeira-Leggett model fails to reliably describe the circuit dynamics when the so-called "lumped element model" breaks down, i.e., when the microscopic behavior of the carriers in the circuit cannot be described anymore by well-separated capacitive, inductive, and resistive elements. This approximation is widely employed in the study of circuit quantum electrodynamics [50]. A detailed discussion about the validity of the lumped element model in electrical engineering is provided in Ref. [74], while some related comments for superconducting circuits can be found, for instance, in Refs. [43,51].

III. CAPACITIVE COUPLING IN CIRCUIT QED
In order to understand how to couple superconducting qubits with resistors, which is the subject of Sec. IV, we first have to get familiar with the capacitive coupling. Indeed, although several ways of coupling elements in circuit QED are available [38], the capacitive coupling is the natural one when we want to connect transmon qubits and resistors.
First of all, let us discuss how the capacitive coupling can be realized between two LC circuits [43]. The coupling scheme is depicted in Fig. 3: two LC circuits, named A and B, are connected via a capacitor with capacitance C g . Following the standard quantization procedure for quantum circuits, we first of all have to write the Lagrangian of the system and then find the Hamiltonian through Legendre transformation. We refer the reader to Appendix A for a detailed presentation of this procedure. The system Lagrangian reads: where we have introduced the flux vector φ = (φ A , φ B ) T and the capacitance matrix defined as Note that, for convenience, we express the potential energies of the inductors as V (φ i ) = 1 2Li φ 2 i . After having obtained the momenta Q i = j=A,B C ijφj , collected in the vector Q = (Q A , Q B ) T , we write the Hamiltonian of the system as and we can proceed with the quantization of the appropriate conjugate variables: Note that, unlike in the case of a simple LC circuit discussed in Appendix A, each momentum Q j is not proportional to the voltage across the jth LC circuitφ j , but is a more involved function of both the voltage differences. In particular, {φ j ,φ k } = C −1 kj , and C −1 reads: We observe that the voltage across each LC circuit can still be considered as the conjugate variable with respect to the node flux, up to a constant with the units of the inverse of capacitance, if the coupling capacitance is low, i.e. under the weak coupling condition C g C A , C B . In this limit, we assume that C A and C B are of the same order, i.e. O(C B /C A ) = O(1); therefore, without losing generality, we choose to employ the perturbation order O(C g /C A ). For instance, the weak coupling condition . Note that the above description is valid for any potential energy V (φ j ).
The off-diagonal elements in C −1 represent couplings between the quantized modes of the system. For instance, with V (φ j ) = 1 2Li φ 2 i , we can see that the off-diagonal elements correspond to a coupling between the p quadratures [75] of two quantum harmonic oscillators, in the form p A p B .
The description relying on the flux vector and capacitance matrix is valid for a general number of capacitively coupled LC circuits. The reader can verify that, for instance, if we deal with three LC circuits coupled capacitively in a line as in Fig. 4, the Hamiltonian of the system includes coupling between all the three modes, but in the weak coupling limit C L , C R C A , C B , C C we find only first-neighbor coupling of the order of O(C L,R /C A,B,C ), while the second-neighbor coupling is of the order of O((C L,R /C A,B,C ) 2 ), and therefore it can be neglected. More generally, we can summarize the action of the capacitive coupling as: when we connect LC circuits in a line through capacitors, the voltage fluctuations across each circuit are coupled only to the first neighbours. This is reflected in the matrix C appearing in the Langrangian. The conjugate variables, however, are not the voltage fluctuations but the momenta, which are all coupled between each other as described by the matrix C −1 . Therefore, in general, the Hamiltonian displays a long-range coupling between the modes of the system. In the weak coupling limit, i.e. when all the coupling capacitors are far weaker than the capacitors of each LC circuit, the long-range interaction is broken and the Hamiltonian describes a first-neighbour coupling as well.
Trivially, the statement above holds not only for LC circuit, but also for superconducting qubits, i.e. when the inductive energy V (φ j ) has a contribution given by a Josephson junction. In the next section, we will observe a very similar behaviour when discussing the coupling between qubits and resistors.
FIG. 4: Three parallel LC circuits in a line, with two coupling capacitors C L and C R . Although from a "diagrammatic point of view" we may say that the LC circuits are connected exclusively to the first neighbors in the line, in the Hamiltonian this is true only in the weak coupling limit, while in general a strong coupling between the three modes of the circuit emerges.

IV. COUPLING QUBITS TO RESISTORS AS THERMAL BATHS
In this section we introduce some relevant examples of resistors acting as thermal baths coupled to superconducting qubits. We focus on the most basic scenarios with one or two qubits, in particular addressing the case of one qubit coupled to a bath, of two qubits coupled to a common bath, with and without a direct coupling between them, and finally of two coupled qubits interacting with separate baths.
For each case, we present the Hamiltonian of the system and, when possible, the associated master equation describing the dynamics of the qubits. We do not provide here the derivation of each Hamiltonian, which is quite involved and can be found in Appendix D. Note that all the cases discussed in this section assume weak coupling between qubits and baths, and for this reason we focus here on Markovian master equations [12]. As we will see, the weak coupling limit is controlled by the ratio between the qubit-resistor coupling capacitance and the internal qubit capacitance. The treatment of the strong coupling limit under the Caldeira-Leggett model, emerging when the qubit-resistor coupling capacitance is of the same order as the internal qubit capacitance, requires a numerical derivation of both the circuit Hamiltonian and of the qubit dynamics -see for instance the discussion in Ref. [76]. The experimental realization of a strong coupling is possible as well, being simply related to the value of the coupling capacitor. In this regime, however, the final form of the Hamiltonian is more involved, and sometimes it may not even be possible to express it in a closed form through the Caldeira-Leggett method, as we show in Appendix D. The open system dynamics in the presence of strong coupling with the bath typically displays non-Markovian dynamics [77], revivals, and has crucial consequences also for the study of quantum thermodynamics [78][79][80]. This regime has been numerically tested in the case where the open system is a single superconducting qubit, for instance, in Refs. [76,81].
Finally, note that, throughout the paper, we represent a superconducting qubit as a capacitor in parallel with an inductive "potential energy" V . Not writing the explicit form of V is quite convenient, because, in this way, we can describe different circuit elements by varying V , such as a qubit or a LC circuit, and at the same time the method of deriving the interaction Hamiltonian does not depend on V . More details about it can be found in Appendix A and in Ref. [11]. For each scenario, we provide the Hamiltonian keeping the abstract notation with the potential energy V , and then we rewrite it in the specific scenario with superconducting transmon qubits. Therefore, our results may be trivially extended to the case of harmonic oscillators coupled to thermal baths, by replacing the transmon qubits with LC circuits.
A. Single qubit Fig. 5 depicts the circuit representing a single qubit, which we name as A, capacitively coupled to a resistor with resistance R. The capacitance of the coupling capacitor is C g . The coordinates of the systems are the flux φ A associated to qubit A as in Fig. 5, and the "fictional" infinite non-interacting internal variables of the resistor. We term as p A the momentum corresponding to φ A . In the Hamiltonian, we describe the internal variables of the resistor as modes labelled by α and created by the operator a † α . The details of the calculation to obtain the final Hamiltonian can be found in Appendix D 1. We now consider the weak coupling limit, which, eluding some subtleties discussed in Appendix D 1 a, reads C g C A . Neglecting all the terms of the order of O((C g /C A ) 2 ), the circuit Hamiltonian is given by: where C Σ A = C A + C g , ∆ω is the infinitesimal frequency gap introduced in the Caldeira-Leggett model of a resistor in Eq. (18), which will be absorbed when computing the bath spectral density, and ω α = α∆ω are the frequencies of the internal modes of the resistor, forming a continuum in the limit ∆ω → 0. If we choose as V (φ A ) the inductive energy of a transmon qubit [20] with frequency ω A , we can rewrite Eq. (63) as: where we have replaced the spectrum of the resistor with Eq. (22) and we have introduced a summation over the "dense" index k, as usually done in quantum optics. λ A is a constant with the units of charge, depending on the features of the transmon qubit [82], which allows us to pass from p A to σ y A through p A ≈ λ A σ y A [20]. We observe that Eq. (64) describes the standard dissipative spin-boson model with an Ohmic spectral density of the bath [12], defined as: where we have used the infinitesimal frequency gap ∆ω to transform the summation into an integral, and we have introduced: If the resistor is in a thermal state with temperature T , in the weak coupling limit the master equation describing the state of the qubit at time t, ρ A (t), is [12]: H LS is the Lamb-Shift Hamiltonian, which only induces a frequency-shift [12]. Γ ↓ (β) and Γ ↑ (β) are respectively the emission and absorption rates, which depend on the bath temperature. By looking at the qubit-bath interaction Hamiltonian in Eq. (64), we observe that it can be rewritten as where V is the operator describing the overall voltage across the resistor, i.e. the one whose correlation function we have computed in Eq. (30). The form of the correlation function, and therefore of the spectral density, does not change in the weak coupling limit [83]. Hence, the dissipation induced on the qubit is driven by the fluctuations of the voltage difference (see for instance the derivation in Refs. [22,83,84]), and we can think of the qubit as a quantum spectrum analyser [52]; the decaying rates are readily obtained: where S V V (ω) is the temperature-dependent spectral density in Eq. (37). Note that Γ ↓ (β) = Γ(n(ω) + 1) and Γ ↑ (β) = Γn(ω), where Γ = πJ(ω A ) is the zero-temperature decay rate, Γ = Γ ↓ (β = ∞). We also have Γ ↑ (β) = exp(−β ω)Γ ↓ (β), reflecting the detailed balance principle from thermodynamics.

Adding an LC filter
The spectral density in Eq. (65) can be shaped by inserting a LC circuit between qubit and resistor as shown in Fig. 6. The LC circuit, with frequency ω f = 1/ C f L f , capacitance C f and inductance L f , is directly coupled to the transmon qubit through the capacitor C g . Once again, we assume the weak coupling limit, i.e. C g C A , C f , where C A is the total capacitance of the qubit. Following the discussion for the previous case, under this assumption the qubit operator σ y A gets coupled to the operator describing the voltage V f at the LC circuit node (see Fig. 6) in the final Hamiltonian. The bath spectral density of this dissipative coupling can be obtained, according to Eqs. (63), (65) and (68), by evaluating the autocorrelation function of the voltage fluctuations V f (t)V f (0) and then the corresponding spectral density S V f V f (ω). The latter can be derived by analysing the circuit in Fig. 6 where, following the convention in electrical engineering, the noisy resistor can be decomposed as a noiseless resistor with resistance R in series with a voltage noise source with suitable spectral density. In the weak coupling limit, the current passing through the small capacitor C g is negligible, therefore we can forget about the transmon qubit and focus on the circuit on the right-hand side of C g ; then, we can find the value of V f (t) as a function of V (t), where the latter is the voltage produced by the noise source, whose spectral density Eq. (37) is known. If Z f (ω) is the impedance of the parallel LC circuit (see Eq. (11)), then current conservation under the weak coupling limit yields: FIG. 7: Two qubits A and B coupled to a common resistor. The resistor is described according to the Foster's first form, and the node fluxes of the internal degrees of freedom (in the present picture not shown) are the ones depicted in Fig. 1.
where h f (ω) is the transfer function. Therefore, , and the open dynamics of the transmon qubit is described, once again, by the spin-boson model, with the prescription that the bath spectral density is modified as J(ω) → |h f (ω)| 2 J(ω). Then, the decay rates can be found through the same formula Eq. (68). A simple calculation gives where |h f (ω)| 2 is a Lorentzian distribution resulting in band-pass filtering with bandwidth ∆ω = 1/RC f and quality factor Several other schemes of RLC filters with different transfer functions can be engineered in order to obtain the desired spectral density [85]. They can also be studied mathematically by means of the Vernon transform [86]. This filtering scheme can be really useful, for instance, to avoid leakages to higher transmon levels, see e.g. the similar strategies for the theoretical proposals in Refs. [83,84] and for the experiments in Refs. [22,87]. Several other schemes that use filters are discussed in Section V.
B. Two qubits, common bath

No direct coupling
Two superconducting qubits coupled to a common resistor through two distinct capacitive couplings, and without a direct coupling between them, are shown in Fig. 7. For simplicity, we have chosen the same value C g for both the coupling capacitors, but other choices are possible, as discussed in Appendix D 2. The system coordinates are the fluxes φ A and φ B associated respectively to qubit A and qubit B as in Fig. 7, and, as always, the "fictional" infinite non-interacting internal variables of the resistor. We term p A and p B the momenta corresponding respectively to φ A and φ B . In the Hamiltonian, we describe the internal variables of the resistor as modes labelled by α and created by the operator a † α . In Appendix D 2 we extensively discuss how to obtain the correct momenta, the modes of the resistor and the final Hamiltonian of the system. In the weak coupling limit, expressed by the condition C g C A , C B , such Hamiltonian reads: where C Σj = C j + C g , and ω α = α∆ω according to Eq. (18). The weight of each mode in the interaction Hamiltonian leads to an Ohmic spectral density as in Eq. (65).
As for the single qubit case, we express the circuit Hamiltonian by choosing as V (φ A ) and V (φ B ) the inductive energy of two transmon qubits [20] of frequency ω A and ω B : where λ A and λ B are constants with the units of charge, depending on the intrinsic values of the energies of the transmon qubits [88], given by the substitution of p j by σ y j , through p j ≈ λ j σ y j . Eq. (72) describes the standard Hamiltonian of two qubits interacting through a common bath, which is at the basis of several collective quantum phenomena [89], such as superradiance [90], entanglement generation [91] and quantum synchronization [92].
The master equation driving the evolution of the two-qubit state ρ AB is given by [93]: where H S is the free Hamiltonian of the two qubits, H LS is the Lamb-Shift Hamiltonian [12] and Γ ↓ jk (β), Γ ↑ jk (β) are the coefficients of the master equation, which depend on the spectral density Eq. (65), on the qubit frequencies and on temperature of the bath, and they fulfil the condition Γ ↓↑ jk (β) = Γ ↓↑ kj (β) * . We refer the interested reader to Ref. [93] for their explicit form. Note that the behaviour of the capacitive coupling between qubits and resistor is analogous to the one observed in Sec. III: in the weak coupling limit C g C A , C B we do not have a direct long-range interaction between the two qubits, while, as proven in Appendix D 2, for strong coupling this is not true. In many situations, we are actually interested in having a direct qubit-qubit coupling in the Hamiltonian. In the next section we show how to obtain it by inserting an additional coupling capacitor.

Direct coupling
Let us consider the circuit analysed in the previous section and add an additional (weak or strong) coupling capacitor C c as in Fig. 8. In the same weak coupling limit, i.e. assuming C g C A , C B , it can be proven that such capacitor induces a direct qubit-qubit coupling (see Appendix D 2 for details). Using the same notation as in the previous section, the Hamiltonian, keeping the general form of the inductive potentials V (φ j ), reads: where S is the matrix The Hamiltonian can be rewritten for the case of two transmon qubits with frequency ω A and ω B , as done in the previous sections: FIG. 8: Two qubits A and B coupled to a common resistor, with an additional coupling capacitor C c directly connecting them. The latter is necessary to induce a direct coupling σ y A σ y B between the qubits in the weak coupling limit. The resistor is described according to the Foster's first form, and the node fluxes of the internal degrees of freedom (in the present picture not shown) are the ones depicted in Fig. 1.
where, as before, λ A and λ B are constants depending on the intrinsic values of the energies of the transmon qubits [20], given by the substitution of p j by σ y j , through p j ≈ λ j σ y j . In this scenario, the Hamiltonian is not diagonal in the canonical basis of the qubits A and B, and the master equation must be derived by finding the eigenmodes of the system. This makes the problem and the final expression of the master equation more complex. We refer the interested reader to Ref. [93] for its specific structure.

C. Two qubits, separate baths
Finally, let us analyse the case in which we couple two different resistors to two qubits, and then couple the qubits at the center of the circuit as depicted in Fig. 9. The qubit-bath coupling capacitors have the value C g , while the qubit-qubit coupling capacitor C c . We find the Hamiltonian of the system in the weak coupling limit, under which all the coupling capacitors are weak compared to the capacitances of the qubits: C g , C c C A , C B . Under this condition, the circuit in Fig. 9 may be seen as the chain of capacitively coupled LC circuits discussed in Sec. III. In analogy with that scenario (see the discussion in Appendix D 3), we can thus write the Hamiltonian in the weak coupling limit as: where S is the matrix and we have neglected any contribution beyond the first order of the weak coupling approximation. The bath frequencies are given by ω α (ω β ) = ∆ωα(∆ωβ), and taking the limit ∆ω → 0 we recover an Ohmic spectral density as in Eq. (65).
FIG. 9: Two qubits A and B coupled through a capacitor with capacitance C c , and independently coupled to two resistors R L and R R , playing the role of separate baths, at least if the coupling capacitances C g are small. The resistors are described according to the Foster's first form, and the node fluxes of the internal degrees of freedom (in the present picture not shown) are the ones depicted in Fig. 1.
Considering A and B as transmon qubits with frequency ω A and ω B , we have: where λ A and λ B are constants depending on the intrinsic values of the energies of the transmon qubits [20], given by the substitution of p j by σ y j , through p j ≈ λ j σ y j (see the discussion for the previous examples). In the weak coupling limit (77) is the Hamiltonian of two coupled qubits interacting with separate baths. This is a fundamental model, for instance, for the study of heat transport in quantum thermodynamics [94]. The requirement C c C A,B induces a weak coupling between the qubits (in the circuit Hamiltonian, S −1 12 S −1 11 , S −1 22 ). This allows us to solve the system dynamics by employing a local master equation [93,95], in which the dissipation is computed locally on each qubit, and the interaction between them appears only in the unitary part of the master equation. If ρ AB (t) is the two-qubit state at time t, its dynamics is given by: where is the system Hamiltonian containing the unitary interaction between the qubits, while H LS is a frequency-shifting Lamb-shift term [12,93]. The local coefficients Γ ↓↑ j (β) depend on the temperature of the bath coupled to the qubit j, and on the corresponding Ohmic spectral density J j (ω) obtained according to Eq. (65); analogously to Eq. (68), they are given by:

V. EXPERIMENTAL METHODS FOR CONTROL OF DISSIPATION
For real-life applications, the simplest dissipative circuit elements can be realized either by using resistive components or lossy transmission lines. In both cases the experimentalist has the choice of integrating them on-chip or to realize separately the dissipative unit from components. In the first case, metal evaporation and e-beam lithography in the same process as the qubits are the methods of choice. Also several types of resistors are commercially available, and in general they can be utilized at low temperatures. Commercial types of resistors include wirewound resistors, metal-foil resistors, thin-film (tantalum, nickelchromium, carbon-boron), carbon-based resistors, and ceramic-metal powder resistors. A semi-infinite transmission line is another way of realizing a resistive element. This type of model can be written down based on the standard procedure of quantization of transmission lines, and the spectral power density can be calculated. Alternatively, the spectral power density can be obtained directly from our formalism by the following argument: a semi-infinite transmission line can be constructedboth in theory and in practice -by terminating the line with a matched impedance. If the line has impedance Z 0 = L /C , where L and C are the inductance and capacitance per unit length, then we add a load of the same value Z 0 (typically 50 Ω). From the elementary classical treatment of transmission lines, we know that in this case the impedance seen at the input of the line is also Z 0 . By applying Eq. (37) we find S V V = 2 ωZ 0 N (ω), which is the same as the result obtained from the quantization procedure. It is also possible to create filtered spectral spectral power densities using transmission lines. For example, consider a transmission line of length l = λ/2 with a shortcut at one end. An analysis of the impedance as seen from the other end shows that the device can be approximated around the resonance as an RLC series circuit with L = lL /2, C = lC /2, and R = lR /2, where L , C and R are respectively the inductance, capacitance, and resistance per unit length. Here the factor 1/2 is due to the formation of standing waves in the line, with voltages and currents varying from zero to the maximum value.
However, fabrication of a standard resistor or the transmission line are not the only methods to produce, and especially to control dissipation. Indeed, to date various forms of dissipation engineering have been proposed or realized in experiments, involving a variety of techniques, from passive methods to active techniques employing tunable devices and microwave-assisted dissipation. Irrespective of the physical realization (typical examples included below are low-Q coplanar waveguide resonators, Josephson junctions, and metallic resistors fabricated on-chip), it is important to understand that the dynamics of the electromagnetic collective degrees of freedom is, generically, the same. For applications such as quantum computing, at first sight the use of any dissipative circuit element is something to avoid, as one aims at high coherences. However, if controlled, dissipation can serve as a way for high-speed preparation of states that are useful for quantum computing, for example for initializing the qubits in the ground state |00000.. at the beginning of a quantum algorithm, or for dumping the excess photons in the resonator used for readout. The latter situation is useful to eliminate the dephasing of the qubit caused by the parasitic photons from the resonator and to reduce the measurement time and the lag between successive measurements. In quantum computing, high-fidelity and fast measurements are essential in running error-correction codes -and more generally for feedback experiments (where the results on measurements performed on ancillas are used to correct the state of the logical qubits).

A. Passive control and circuit design
In mesoscopic physics, dissipative elements have traditionally played a role mostly in realizing filters and ensuring good thermalization. But recently they took a more prominent role [6] as thermal reservoirs in quantum thermodynamics -the study of fluctuations and small-scale thermal engines at the single-quantum level. For instance, they have been key ingredients for some proposals and realizations of calorimeters in circuit QED [25,96,97], with applications to the measurement of work in a quantum system [98,99]. Moreover, they play a major role in the study of quantum heat transport [100,101], in the proposal for a quantum heat switch [102] and an Otto refrigerator [103], in an experimental demonstration of quantum-limited heat conduction [26], and in the creation of a heat sink for quantum circuits [104]. In the experiments with a qubit operated as a heat valve [22] as well as in the demonstration of heat flux rectification [87] and of electrostatic control of radiative heat transfer [105], the reservoirs were fabricated as Cu-metal resistors on the same chip. The resistors can be Joule-heated, thus providing the temperature difference needed to study heat transport phenomena, and, in the future, to build mesoscopic-scale thermal engines. The slow electron-phonon relaxation is an advantage for the use of resistors as thermal baths, since they do not require excessive Here h(ω) is the transfer function -the Fourier transform of the impulse function h(t), see Fig. 10. One such very useful linear time-invariant device is a band-pass filter, which allows us to shape the spectrum of the noise in very useful ways. We have seen already an example of this in Section 4.1.1, where we analysed in detail the effect of an LC filter. An ideal bandpass filter with low-pass frequency ω L > 0 and high-pass frequency ω H > 0 has the transfer function Since we are interested in voltage power spectral density, the shaping action of a filter with transfer function h(ω) is given by We assume here that the quantum system at the output of the filter, which acts as a load, has a very large impedance in the bandwidth region, therefore its effect is neglected. It is important to notice that the action of the filter applies both to negative and positive frequencies ω; that is, the filter works in the same way both when the system absorbs and emit energy. The bandwidth ω H − ω L of the filter can be designed narrow enough such that it ensures selective coupling with the quantum system of interest (with only one of the transitions for a transmon qubit or with the resonant frequency if it is a superconducting resonator). Two immediate applications of this scheme are the realization of Purcell filters and the filtering of thermal reservoirs in quantum thermodyamics. The readout resonator used in circuit QED to measure the qubits has to satisfy contradictory requirements: it needs to have a high quality factor so that the qubit is not damped and at the same time it needs to allow photons to pass in and out such that the measurement is performed in an as short measurement time as possible. To protect the qubit, it is easy to design circuits such that the qubit frequency is placed on the tail of the response curve of the dissipative element, as far as possible from the maximum. However, additional filtering offers an even better protection, and also enables simultaneous measurement of several qubits by multiplexing. For example, single-pole bandpass filtering of the signal from the resonator is a good compromise for the requirements above, preserving a high T 1 > 10 µs and allowing a high-fidelity (99.8%) and at the same time fast (120 ns) measurement of the qubit [106]. In this case, the filter (with quality factor of about 30) was realized on-chip as a single λ/4 coplanar waveguide resonator. The readout resonators are placed in the passband, thus preventing the dissipative modes to dampen the qubit and at the same time allowing for a significant outcoupling for the modes around the resonator frequencies. Another design uses two λ/4 resonators [107] to achieve a zero real admittance environment at the qubit frequency, thus preventing the spontaneous decay through the outcoupled resonator modes. In this case, an improvement in T 1 by a factor of 50 was observed at a qubit frequency of 6.7 GHz -at the same frequency at which the transmission through the resonator followed by the Purcell filter showed a 30 dB drop. Finally, in quantum thermodynamics the use of resonators placed between resistors and the qubit and acting as spectral filters is a necessary requirement in order to achieve selective thermalization of the qubit with either one of the reservoirs [22,87].
b. Tuning of dissipation Tunability is an important property of a dissipative device, opening up the way for many applications. The easiest way to achieve this is to incorporate a SQUID into the device, and change the critical current by the use of an external magnetic field. Proposals along this line include the dissipator as a frequency-tunable Josephson junction with small quality factor (≈ 100) [108]. A design where two resonators, one with high quality factor and another one tunable and coupled to the ground by a resistor, has been studied in Refs. [83,104]. It was observed that the loaded quality factor of the high-Q cavity (10 5 ) can be tuned down by a factor of almost 100. In this experiment the resonators were made of Nb, while the resistor was made of Cu (R=375 Ω). Another scheme for creating wideband tunable dissipative elements has been studied [109]; the device consists of a double chain of SQUIDs capacitively coupled to each other, which sustains modes with wide dispersion. Dissipation in this circuit is realized by the capacitive coupling to a transmission line of both of the SQUID chains. These devices can be employed to efficiently initialize the state of a superconducting qubit [84], a well-known fundamental task in circuit QED. Further experiments could use photon-assisted tunneling in an NIS junction to cool a device such as a qubit or a resonator [27]. Present practical methods used in more advanced quantum processors [110] for initialization also make use of tuning of dissipation, based on the observation that it is more convenient experimentally to tune the qubit rather than the bath. In a standard transmon qubit coupled to the readout coplanar waveguide cavity, the qubit (with a higher frequency than the cavity) is adiabatically tuned down towards resonance with the cavity, which swaps the excess population from the excited states of the qubit into the cavity. The qubit is then held at a frequency lower than the cavity. Meanwhile, the excitations transferred to the cavity are fast dissipated into the environment. Note that the cavity, with decay rate κ ≈ 1/(45 ns), relaxes much faster than the qubit (T 1 ≈ 14 µs), therefore it truly plays the role of a dissipative environment. After the cavity is sufficiently depleted (typically 200 ns), the qubit is returned diabatically to the operating frequency. During this last step, it crosses again the cavity, but it will not get populated since the cavity is empty. The whole process lasts only 280 ns and achieves a ground state with error below 0.5%.
Finally, certain tasks in quantum thermodynamics, for example related to the operation of proposed Stirling engines [111], require a fast-tunable coupling between a qubit and a reservoir. For this, various designs for tunable coupling have been developed in the context of making two-qubit gates in quantum computing [112,113], which can be easily adapted to coupling a qubit to a bath.

B. Active control of dissipation by coherent drives
An interesting recent development is the realization of schemes where dissipation is controlled by coherent microwave drives combined with filtering.
For example, a reset method using a Josephson photomultiplier as a photon counter was demonstrated for a transmon qubit [114]. The Josephson photomultiplier consists of a capture cavity with relatively low quality factor (≈ 1300) capacitively coupled to a Josephson junction placed in a rf-SQUID and biased near the phase-slip critical flux. In this experiment, the dissipating device was fabricated as a standalone device from Al on high-resistivity Si, packaged in an Al box, and connected by a coaxial wire to the measuring resonator of the qubit. With this device, a deterministic reset of the measuring cavity can be realized by a two-pulse sequence. The first is a fast flux pulse that resets the photomultiplier itself by putting it in a single-well potential configuration. The second pulse puts the qubit cavity in resonance with the capture cavity (and nearly-resonant with the qubit cavity), with the result of depleting both cavities of spurious photons. Ramsey interference fringes show a clear increase in visibility as the depletion time is increased. Another method based on two microwave drives was demonstrated in Ref. [115].
In the same vein, a scheme for initialization or reset can be constructed by using a parametric modulation of the qubit frequency ω(t) = ω + A cos(νt), which creates sidebands [116]. When the n's sideband becomes resonant with the readout resonator ω d (typically lower than ω), nν + ω d = ω, transfer of population can occur. The microwave-induced coupling [117] between the qubit and the resonator results in the interaction Hamiltonian between the qubit and the lossy cavity (written in the rotating frame) where σ ge = |g e|, d is the mode of the lossy cavity (in this case the readout cavity), g is the bare (unmodulated) qubit-resonator interaction, and J n is the Bessel function of first kind. Swapping the population between the qubit and the resonator, due to the results in the suppression of the residual excited population to 0.08 % in only 34 ns [116]. This type of interaction has multiple application in circuit QED, such as entanglement production [118], realization of cross-resonant qubit-qubit gates [119,120], motional averaging effects [121], and Landau-Zener-Stückelberg-Majorana interference [122,123]. Selective coupling to dissipative elements is not limited to two-level systems. Transmon circuits have higher energy levels that can be accessed experimentally. In Fig. 10 (b) we show such a configuration, where a transmon is truncated to the first three-levels |g , |e , and |f . The e − f transition is coupled to a dissipative element via a filter ensuring that the g − e transition is off-resonant, therefore not affected. If the g − f transition is pumped resonantly at a frequency of half the sum of the g − e and e − f transition frequencies, two-photon processes promote the population to state f , from which it decays fast to state e. This results in population inversion in the {|g , |e } manifold, stabilizing the excited state |e . That this is the case can be seen by considering the truncation of the Hamiltonian in the {|e , |f } manifold, which in the rotating-wave approximation is where σ ef = |e f |, σ + ef = |f e| and g is the coupling with the dump cavity d. Since d is a highly dissipative mode, the stabilized state |ss will necessarily have the cavity in the vacuum state |ss |0 d . Enforcing that this is an eigenvalue of the above Hamiltonian results in the condition σ ef |ss = 0 or |ss = |e . Note also that the two-photon process is sufficiently off-resonant with both the g − e and e − f transitions. Another scheme that uses three levels and is based on a microwave drive coupling the states |f, 0 and |g, 1 of a transmon coupled to a dissipative broadband resonator with κ = 56.52 MHz has achieved a fast reset (500 ns) resulting in 0.2% residual population [124]. Subsequently, it was shown that the read-out resonator can serve as the dissipative dump resonator, therefore reducing considerably the number of circuit elements needed [125]. This concept can also be generalized to multiple qubits. For example, similar ideas [126] were used to create arbitrary superpositions of the ground and excited state [28], to autonomously stabilize two transmon qubits in a three-dimensional cavity in a Bell state [29], and to produce a stable bosonic Mott insulator state in an array of qubits [127].
Dissipation, when balanced with amplification, gives rise to a generalization of standard unitary evolution known as PTsymmetric quantum mechanics [128]. Experiments on exceptional points and broken PT -symmetry have been realized on many experimental platforms. In circuit QED, a digital simulation of single-qubit and two-qubit non-Hermitian dynamics has been realized on the IBM superconducting processor [129]. It is possible to emulate the physics of these Hamiltonians in a dedicated transmon in a three-dimensional cavity setup [130], following the same principles of noise filtering as presented above, but with the interesting difference that now the dissipative channel is coupled selectively to the g − e transition, see Fig.  10 (c). Experimentally, this is done by increasing the decay rate Γ e through the insertion of a mismatching element between the cavity and the parametric amplifier. This produces standing waves in the transmission line, resulting in different density of modes at different frequencies and therefore different Purcell rates. Since the transmon is tunable, it is possible to match the g − e transition frequency to a frequency region with large density of states resulting in a dissipation rate Γ e = 5.25 µs −1 , more than one order of magnitude larger than Γ f = 0.25 µs −1 . The condition Γ e Γ f has as consequence an effective dynamics on the {|e , |f } manifold governed by where ∆ is the detuning of the drive from the e − f transition and Ω is the coupling strength. Note that in this case we are not interested in the steady-state imposed by dissipation (which, because of the conditions above, is trivially just the ground state for the transmon), but on the dynamics on a timescale shorter than 1/Γ e . One recognizes here the paradigmatic non-Hermitian Hamiltonian for qubits, with the dissipative term on the diagonal and a σ x -coupling in the off-diagonal elements. Both in the case of digital simulation [129] and in the case of emulation by a threelevel system, postselection plays a key role: each experimental run ends with a measurement and only the data corresponding to the correct manifold is retained. In the digital simulation [129] this is flagged by the ancilla qubit being projected in the state 1, while in the three-level system emulation this is enabled by the single-shot quantum tomography, where the occurrences of projection to the ground state are eliminated.
With continuous variables, such an active method can be employed for the stabilization of highly squeezed states (exceeding the standard 3 dB achieved by parametric drive with a single pump) in a cavity, by using an additional dump cavity [131]. The concept is similar to the case presented above in Fig. 10, with the difference that now we have a harmonic oscillator instead of a transmon. In addition, in this experiment the two cavities were coupled to each other by a Josephson ring modulator, which enables the modulation of the couplings by applied microwave tones. The cavity had resonant frequency ω c and decay rate κ c , while the dump cavity had ω d > ω c and the decay rate κ d κ c (κ d /κ c = 200 in the experiment). By pumping the Josephson ring modulator simultaneously at the difference frequency ω − = ω d − ω c and at the sum frequency ω + = ω d + ω c one obtains by standard methods involving the rotating-wave approximation (RWA) is an effective beam-splitter Hamiltonian and H + = g + (c † d † + cd) is an effective parametric down-conversion Hamiltonian. The truly interesting feature emerges when the high damping of mode d is considered. Indeed, arranging the phases and amplitudes of the pumps such that g ± are positive and real with g − > g + , the full effective Hamiltonian can be rewritten as a beam-splitter Hamiltonian between the dump mode and a Bogoliubov mode where G = g 2 − − g 2 + and B is the Bogoliubov mode of the cavity c, defined by B = c cosh r + c † sinh r with squeezing coefficient r = tanh −1 (g + /g − ). Now, if κ d κ c , dissipation forces the dump resonator into the vacuum state |0 d . As a result, by inspecting the Hamiltonian Eq. (85), we conclude that the steady-state |ss c of the cavity c must be the vacuum state of the Bogoliubov mode β|ss c = 0, which is a squeezed state with squeezing parameter determined by the ratio g + /g − . In principle, one can increase r to very large values by having g + nearly equal to g − . This however reduces the coupling G, so in practice an optimal value has to be found. In the experiment [131], an optimal squeezing slightly above 8 dB was obtained. Another way of seeing this is by analysing the time it takes to reach the steady state, which is of the order of κ d /G 2 . This time increases as we aim at higher squeezing factors by increasing g + to approach from below g − , because G decreases as well in this situation.
Along the same lines, an important application of these concepts is the realization of a single-photon detector [132]. This device uses two cavities (a buffer cavity where the photon to be detected is first stored, and a second highly-dissipative waste cavity), with a coupling transmon qubit in-between. Denoting the buffer and waste annihilation operators by b and w, and the qubit lowering operator by σ, it is possible to derive an effective Hamiltonian where g 3 is an effective three-wave mixing rate, activated by the parametric pumping of the transmon. This Hamiltonian converts a photon incident on the buffer into a qubit excitation and a waste photon. To realize the effective Hamiltonian Eq. (86), the transmon is pumped at a frequency ω p such that ω p =ω ge + ω e w − ω g b , whereω ge is the ac-Stark shifted qubit frequency, ω e w is the waste-cavity frequency with the qubit being in the excited state, and ω g b is the buffer cavity frequency with the qubit in the ground state. When the buffer contains one photon, the pump provides the missing energy of transforming that photon into an excitation on the qubit and a photon in the waste. Now, if the dissipation rate k w of the waste is such that |g 3 | k w , then the waste photon is eliminated fast into the environment, while the qubit remains in the excited state. By adiabatically eliminating the waste, a non-local and non-linear loss operator 2|g 3 |k −1/2 w bσ † is obtained, which shows that the qubit goes into the excited state conditioned on the presence of a photon in the buffer. Finally, the excited state of the qubit is detected by a standard singleshot dispersive technique. This detector has a detection efficiency of 58% and a dark count rate of 1.4 photons per millisecond.

VI. CONCLUSIONS AND PERSPECTIVES
Modeling of electrical circuits comprising resistors and the other dissipative elements has experienced a new momentum recently, driven by burgeoning applications in circuit quantum electrodynamics. Here we have reviewed some of the standard models in the field, and we outlined the main techniques for coupling resistors with qubits and controlling dissipation. This opens a number of interesting perspectives for future works, as many theoretical proposals rely on engineering thermal baths acting on quantum devices.
In quantum information and quantum computing, as we have already seen, engineered reservoirs provide an advantage with respect to adiabatic state preparation, which is slow -since it should be slower compared to the energy gap. Also, in dissipation engineering of quantum states the final state is protected from perturbations -implementing in fact a form of self-correction. Tunable Markovian dissipation provides a scheme for universal quantum computation [34,133] and for the characterization of information-preserving steady states [134,135]. Moreover, non-Markovian effects may provide new insights into the exchange of information between a system and the environment that surrounds it [136][137][138]. Looking forward to applications in circuit QED, ideas exploring quantum optimal control and dissipative networks could allow for novel reservoir-engineering reset schemes of superconducting qubits [139][140][141]. Furthermore, resistors acting as thermal bath can be used to study how non-Markovianity effects affect the work done by a driving force on a superconducting qubit [142], and to address different collective effects in a pair of transmon qubits immersed in a common bath [89]. Besides, they have been used to theoretically explore quantum phase transition with dissipative frustration [143]. More information can be found also in Ref. [144]. Finally, a very recent proposal of the experimental observation of the quantum trajectory of a single microwave photon detection is based on a Cu-metal resistor that plays again the role of a thermal bath [145].
In superconducting circuit simulators, access to each site will enable further experiments on many-body localization, impurity (defect) dynamics and non-equilibrium thermodynamics such as pre-thermalization, allowing for further manipulation (driving the system adiabatically from incompressible to compressive phases) of the many-body photonic state. There are theoretical proposals for realizing a photonic Bose-Einstein condensate in strongly driven microcavity arrays using two interacting superconducting qubits each coupled to a neighbouring cavity, where the dissipation here would be realized by the spontaneous decay of the antisymmetric mode [146], which produces scattering of photons in low-momentum states thus stimulating condensation in the zero-momentum state. Exotic state preparation using engineered dissipation is possible in circuit QED [147], drawing on ideas of observing photon blockade and the Mott-insulator-to-superfluid transition with polaritons [148,149]. Other possibilities include dissipative preparation of topological superconductors [150], parametric thermalization [151] and creation of compressible phases and exotic nonequilibrium processes associated with generation of entropy [152].
Finally, in quantum thermodynamics, proposals to realize heat engines in circuit QED will require filtering of reservoir noise together with precise control of the operation point of the qubit [111]. Dissipation-by-measurement is also an emerging, fundamentally new approach [153].
To conclude, in the near future, tailored dissipation will be more and more employed as a resource for the preparation of otherwise inaccessible or unstable many-body states and for targeting and purifying states with high degree of correlations. Engineering a thermal bath acting on quantum devices is of great importance for both fundamental physics and technological applications. We also notice that some of the techniques we have presented, related to control of dissipation by the application of external fields, are not limited to the superconducting circuits. These concepts can be implemented in platforms such as trapped ions or ultracold gases, bringing in new physics and additional opportunities for experimental control.
We discuss here the original method to quantize the degrees of freedom of electromagnetic circuits, relying on the quantization of the branch fluxes of the circuital network [45,62] (for a different approach, see [154]). The flux of a certain branch of the circuit is defined as the integral of the voltage difference across it: Straightforwardly, we haveφ(t) = V (t). Once we have defined the fluxes of all the circuit branches, we can identify the degrees of the freedom of the system by using the standard method of nodes [45]. This is why the fluxes of the circuit are usually called "node fluxes": the reduction given by the node method allows us to interpret the remaining relevant fluxes as associated with each node of the network, once we have defined a proper ground point. This can be observed in Fig. 11, representing a parallel LC circuit: we only have one degree of freedom given by the node flux φ, equal to the fluxes of both the capacitive and inductive branches under a suitable choice of sign. More complex networks require a more involved analysis (see [45] for a deeper discussion).
Let us now focus on the LC circuit shown in Fig. 11. Assuming the lumped element model is valid, all the inductive elements of the circuit are described by the inductor in the left branch of circuit. Therefore, making use of Faraday's law, we observe that φ is, up to a choice of sign, the magnetic flux through the inductor, related to the current I running in the circuit by φ = LI, where L is the inductance. We recall that the electromagnetic energies contained in a perfect capacitor and inductor respectively read U C = 1 2 CV 2 and U L = 1 2L φ 2 , where C is the capacitance, L the inductance, V the voltage across the capacitor and φ the flux through the inductor. Then, we can consider U C as the kinetic energy and U L as the potential energy of the system, and we can introduce the Lagrangian of the parallel LC circuit as: One can verify the validity of this Lagrangian by writing the Euler-Lagrange equations and observing that they correspond to the proper Kirchoff's circuit laws. We can now obtain the Hamiltonian by finding the momentum and applying the Legendre transformation: with Poisson bracket {φ, Q} = 1. Note that Q represents the physical charge on the island associated to the node flux φ. The validity of the Lagrangian method explained above is not restricted to the parallel LC circuit: it holds for any network made of inductive and capacitive elements, since it correctly describes each term we would usually write in the Kirchhoff's laws. However, it does not consider the last category of passive circuit elements, i.e. resistors, since it would include a cumbersome and unwieldy expression to describe the energy dissipated by the resistor at time t. How to include resistors in the Hamiltonian formalism is the subject of Sec. II.
FIG. 12: Effective diagram to describe a Josephson junction, or, more generally, a superconducting qubit. The "potential energy" V (φ) is the inductive energy of the junction, namely V (φ) = −E j cos ϕ = −E j cos 2π φ φ0 , where E j is an energy constant and φ 0 is the magnetic flux quantum. Now that we have identified the Poisson brackets of the degrees of freedom of the system and written its Hamiltonian, we are ready to quantize the circuit. However, first of all we should explain the meaning of such a procedure for electromagnetic circuits; indeed, it is not trivial to understand when and why a mesoscopic circuit behaves quantum-mechanically, i.e. when all of its constituents are not well described anymore by a classical approximation. The requirements for the quantum regime basically reads [43]: the circuit must be made of superconducting materials, and must be kept at a temperature (usually around 20 mK) such that the thermal excitations are not comparable with the superconducting energy gap necessary to break a Cooper's pair. Moreover, additional excitations such as plasmons are usually of a frequency way higher than the relevant frequencies of the circuit (typically of a few GHz). Therefore, the electrons in the circuit are bulked in a superconducting collective state that cannot be broken by any excitation, and a quantum-mechanical description is meaningful and necessary. We refer the reader to Refs. [37,43,45] for an extensive discussion on this topic.
Let us go back to the Hamiltonian of the parallel LC circuit in Eq. (A3). If the conditions discussed above are fulfilled, we can quantize the macroscopic flux and charge of the circuit and write: From now on, we drop the hat sign over the operators. The quantum circuit Hamiltonian is now the Hamiltonian of a quantum harmonic oscillator with frequency ω = 1 √ LC , which is the resonant frequency of the LC circuit. We can define the standard creation and annihilation operators as: and the Hamiltonian reads A necessary circuit element to create superconducting qubits is the Josephson junction, i.e. two superconducting electrodes connected by a thin insulating barrier [155]. The first ever superconducting qubit, the so-called charge qubit, essentially consisted of a Josephson junction shunted by a voltage source [9]. As depicted in Fig. 12, the Josephson junction can be described by a capacitor in parallel with a non-linear inductance, giving rise to a "potential energy" of the form U L = −E j cos ϕ, where E j is an energy constant depending on the junction and ϕ is a variable proportional to the flux φ through it. By employing this non-linear effect, it is possible to isolate the two lowest quantum states of the junction Hamiltonian and to fabricate a qubit (see Refs. [1,36,43] for a review). For the purpose of Hamiltonian analysis of a quantum circuit, we describe superconducting qubits, or Josephson junctions, as circuit elements made of a capacitor in parallel with a source of (inductive) potential energy V (φ), as shown in Fig. 12.

Appendix B: Fourier Transform in quantum mechanics versus signal analysis
We give here some clarifications on the opposite sign choices in the definition of Fourier transform in quantum mechanics and signal analysis.
Consider for example a quantum harmonic oscillator with frequency ω 0 and Hamiltonian H = ω 0 (a † a + 1/2). The Heisenberg equations of motion yield a(t) = e −iω0t a(0) and a † (t) = e iω0t a † (0), and the position operator is where m is the mass. We can see that both components appear, one corresponding to the annihilation operator and rotating clockwise in the complex plane, and the other corresponding to the creation operator and rotating counterclockwise. None of them is "more fundamental" and since they are non-Hermitian operators they are not observables. The observable is the Hermitian position operator x(t). However, it is convenient to write the equations of motion in terms of the annihilation operator a, to avoid carrying along the dagger sign. Thus, for the Fourier transform we would like to have a convention that attributes the Fourier component at frequency ω 0 to the annihilation operator. The definition Eq. (1) clearly achieves this, since a(ω) = 2πa(0)δ(ω − ω 0 ), but we stress that in the Hermitian x(t) both Fourier components are present, Now, let us examine the conventions from signal analysis and electrical-circuits engineering. Consider a generic harmonic signal x(t) = X cos(ω 0 t + φ), which can be also written as with Fourier transform Again it is clear that both frequency components are present but we have to make a choice. This choice in electrical engineering is motivated by the representation of complex numbers as phasors. Phasors can be summed by the usual rules of voltage addition, thus, several signals at the same frequency can be analysed in a time-independent framework. For the signal above, the standard definition of its phasor is Xe iφ , with the property that angles are measured counterclockwise in the complex plane, and the real signal x(t) is then obtained as x(t) = Re[Xe iφ e iω0t ]. However, we can see that there is nothing fundamental in this choice. Equally well, we could designate Xe −iφ and follow its evolution, and finally write x(t) = Re[Xe −iφ e −iω0t ]. Doing so however changes the standard textbook definition of the complex ac impedances and admittances for standard components such as inductors and capacitors, as used in electrical engineering. If we wish to recover these definitions, we simply have to replace ω → −ω in all the results of the paper.

Appendix C: Linear response function and generalized impedance
Let us start with the definition of complex impedance Eq. (4): Provided it is well-defined, we can take the Fourier transform of the above equation and obtain: and assuming that the integration order can be switched, which coincides with Eq. (13). Both the above assumptions, however, are quite strong, and for instance the impedance function of a parallel LC circuit does not fulfil them. We clearly understand that this is not a physical problem, but a mathematical issue due to the divergence of the impedance function in a zero-measure set. We can therefore try to avoid this pathological point of the frequency range by slightly changing the definition of impedance. The trick consists in introducing a real parameter η in Eq. (C3), and then taking its value to zero so that it cannot change any physical result [156], as done in Eq. (15). The negative exponential does not cause any trouble because Z(t) vanishes for negative times [53]. At this point, we can follow the inverse path that led us to Eq. (C3), and write: Thanks to the negative exponential, the Fourier transform of Z(t) is now well-defined and the integrals can be interchanged. Hence, we have finally obtained the desired relation Eq. (C1) with the new generalized impedance defined in Eq. (14).

Single qubit
In what follows we discuss the coupling of one qubit to a resistor via a capacitor with capacitance C g , as shown in Fig. 5. We describe the resistor using the Caldeira-Leggett model presented in Sec. II B, thus considering as internal degrees of freedom the voltage differences across each LC circuit of the Foster's first form, associated to the fluxes φ j with j = 1, 2, . . .. From now on, we use greek letters to refer to such internal degrees of freedom, therefore we write φ α instead of φ j .
The Lagrangian of the circuit in Fig. 5 is readily written: where we have introduced the vector φ = (φ A , φ 1 , φ 2 , . . .) T and the capacitance matrix where e α = (1, 1, . . .) T is an infinite vector filled with ones, and C α = diag(C 1 , C 2 , . . .) contains the values of the capacitances of the internal degrees of freedom of the resistor. We have also introduced the inductance matrix, defined as where L −1 α = diag(L −1 1 , L −1 2 , . . .). Our aim is to apply the Legendre transformation to the above Lagrangian and obtain the Hamiltonian. Moreover, we want such Hamiltonian to be in a readable form, that is to say, we want it to be divided into three parts: the free Hamiltonian of the qubit, the free Hamiltonian of the resistor as a thermal bath, and the interaction Hamiltonian coupling the two elements. This corresponds to writing, in the final Hamiltonian, a capacitance matrix C −1 such that the (infinite) α × α submatrix associated to the degrees of freedom of the resistor is diagonal, and the same for the inductance matrix. In order to do so, we rely on a method first employed by Paladino et al. [157], and very recently improved in Ref. [21].
The idea consists in searching for a coordinate transformation (a point transformation in the language of Hamiltonian mechanics) such that the new coordinates read z = Zφ, where Z is a symmetric matrix, and the new associated capacitance (inductance) matrix is given by . We choose to parametrize such point transformation as with M α = C α + ξe α e T α , where ξ is a parameter we can freely choose, and M 0 is a free constant with the units of capacitance, whose value does not induce any physical effect. The new capacitance matrix reads where for convenience we have introduced the notation C Σ A = C A + C g , and f α = M 1/2 0 M −1/2 α e α . In order to obtain the Hamiltonian as a function of z and their conjugate momenta, we now have to find C −1 z . This inversion is obtained through Eq. (E4) proven in Appendix E, in turn derived from a formula presented in Ref. [21]. The inverse matrix reads: with coefficients (D7) In order to remove any interaction between the modes of the resistor, from Eq. (D6) it is clear that we need to set δ = 0, i.e. we choose the parameter ξ = C g − C 2 g /C Σ A , and finally we obtain |f α | 2 is the square norm of an infinite vector, hence we need to worry about a possible divergence. Fortunately, as proven in Ref. [21] this norm always converges, since by using the Sherman-Morrison inversion formula we get therefore which is finite even in the pathological case e T α C −1 α e α → ∞. Anyway, if we choose the values of C α to represent a resistor as in Eq. (18) and Eq. (22), this quantity remains finite as well.
If p are the momenta associated to the transformed coordinates z, the Hamiltonian is now given by We have finally managed to decouple the internal modes of the resistor in the capacitance matrix. Unfortunately, this is not true for the inductance matrix, which now reads Therefore, in order to get the Hamiltonian into the desired form, we need to find the unitary transformation that diagonalizes L −1 z , by calculating the eigenvalues and eigenvectors of Eq. (D12). This requires finding M −1/2 α , and none of these tasks can lead to an explicit formula for a generic system. After knowing the value of all the relevant circuit parameters, we would most likely need to cut off a certain number of modes in order to deal with finite matrices, and then to perform the diagonalization numerically.
It is worth pointing out that a closed-form expression for the Hamiltonian can, in general, be found by following a different procedure, which relies on a description of the dissipative part of the network as a generic transmission line. This leads to solving a boundary condition differential equation singular value problem. The essential difference lies in the fact that the internal modes of the resistive element are there treated as a continuum, in contrast with the Caldeira-Leggett method we employ in this paper, where the modes are discrete and obtained through a Foster's decomposition. This alternative approach has been differently introduced and developed in Refs. [21,65,158,159], and an extensive review can be found in the Ph.D. thesis by Adrián Parra-Rodriguez [64]. As already said, in this paper we do not discuss this approach, and we just focus on the Caldeira-Leggett method, so as to recover the standard description of a thermal bath as a collection of infinite, discrete bosonic modes. We will later see that a closed-form expression for the Hamiltonian can be obtained also in the weak coupling limit of the Caldeira-Leggett model.
The suitable unitary transformation U takes the form and the variables transform according to z = U z, p = U p. U is a canonical transformation, i.e. it preserves the Poisson brackets. L −1 z is diagonalized through L z −1 = U L z −1 U T , while the capacitance matrix is rotated according to C z −1 = U C z −1 U T , which only affects the coupling vector f α , transformed into f α = U α f α .
Finally, the Hamiltonian H = 1 2 p C z −1 p + 1 2 z L z −1 z + V (φ A ) can be rewritten as [160] H = (C z −1 ) 11 We observe that, by capacitively coupling a qubit to a resistor, the interaction Hamiltonian displays a qubit-bath coupling through the qubit momentum p A with the units of charge (through σ y in the case of a transmon qubit), in accordance with the discussion on the capacitive coupling in Sec. III. Therefore, the capacitance matrix inversion and the search for a diagonalizing unitary U are necessary only to obtain the weight of the coupling to each bath mode, described by the vector f α in Eq. (D14), which also defines the spectral density. We can absorb the constant M 0 by transforming the bath quadratures into annihilation and creation operators of each mode (respectively a α and a † α ) [70]: ω α = (L −1 z ) αα /M 0 are the frequencies of each resistor mode, forming a continuum, while the constant M 0 in the interaction Hamiltonian is going to be absorbed by the coefficients f α , according to their definition in Eq. (D5). Following the convention of quantum optics, we drop the infinite constant term in the resistor Hamiltonian. Finally, note that the qubit Hamiltonian H q = 1 2 (C z −1 ) 11 p A + V (φ A ) has been renormalized by the interaction with the resistor, since in general (C z −1 ) 11 = C −1 A . This is an important effect of circuit QED, which may be negligible in the weak coupling limit, but may also display non-trivial effects in broader scenarios [157].

a. Weak coupling limit
Despite a closed-form expression for the Hamiltonian is in general not available, i.e. we cannot know a priori all the coefficients of Eq. (D14) as functions of the circuit parameters, in the weak coupling limit we are able to recover a closed-form expression. Defining the weak coupling limit in the presence of a resistor is less trivial than in Sec. III. Indeed, while there we could set C g C A , C B , the internal degrees of freedom of the resistor are related to an infinite collection of capacitances stored in the diagonal matrix C α , expressed in Eq. (18). We set the weak coupling limit as C g C A , C g (C α ) jj . If this is the case [161], we can neglect the effects of the order of C g /(C α ) jj in the point transformation and in the renormalization of the internal modes of the resistor, and the qubit-bath coupling is weak compared to the energy of the qubit (with a perturbation parameter of the order of O(C g /C A )). In this case, recalling Eq. (D9) with ξ = C g − C 2 g /C Σ A , at the zeroth order we have M −1 α ≈ C −1 α , and we obtain the final Hamiltonian: Note that requiring C g (C α ) jj is sufficient to obtain a closed-form expression for the interaction Hamiltonian, while the additional condition C g C A introduces the standard weak coupling limit in the language of open quantum systems, for which the system-bath interaction is a perturbation that can be treated through the Born-Markov approximation, leading to the master equation (67).

Two qubits, common bath
The two-qubit scenario is similar to the single qubit one presented in Sec. IV A, although there are some subtleties that are worth discussing. Let us start with the case where the two qubits, named A and B, are capacitively coupled to the same resistor, but not directly coupled together, as depicted in Fig. 7. This corresponds to the well-known "common bath" scenario in the language of open quantum systems. For simplicity, we consider the balanced case in which the two coupling capacitors have equal capacitances; the derivation of the Hamiltonian in the unbalanced case can be obtained by following the same method we discuss here.
The circuit Lagrangian is: where we have used the same notation of the previous section in Eq. (D1), with φ = (φ A , φ B , φ 1 , φ 2 , . . .) T . The capacitance matrix reads: where a = (1, 1) T , e α = (1, 1, . . .) T and S is the system matrix, describing the qubits: Note that any possible imbalance between the qubit-bath couplings, e.g. when the coupling capacitors have different capacitances, is reflected in the elements of the vector a, which in this case would not be equal anymore. The inductance matrix L −1 is the same as in Eq. (D3), with an additional zero row and column because of the second qubit. Following the path of the previous section, we introduce new coordinates through the point transformation Eq. (D4) (once again with additional zero row and column), and obtain the new capacitance matrix given by with D = M 0 + |f α | 2 (2C g − ξ) − C 2 g a T S −1 a . In order to remove the capacitive interaction between different modes, we choose ξ = 2C g − Tr S −1 C 2 g , so that the submatrix describing the resistor modes (block (2, 2) of Eq. (D21)) is proportional to the identity matrix and D = M 0 . Therefore, the inverse capacitance matrix is now given by where s = S −1 a = (C −1 Σ A , C −1 Σ B ) T , and We now face the same issue as in the previous section: while there is no capacitive interaction between the resistor modes, the inductive matrix is not diagonal. We need to numerically find the suitable canonical transformation Eq. (D13) in order to completely decouple the modes of the bath, and find the spectral density of their interaction with the qubit. Anyway, although we are not able to find a priori a closed-form expression for the Hamiltonian, we can already observe that the qubits are not only coupled to the same bath, but a direct qubit-qubit interaction term appears in the Hamiltonian, proportional to σ y A σ y B in the case of two transmon qubits. This is the effect of the non-diagonal matrix N appearing in the block (1, 1) of Eq. (D22), and is in accordance with the result found when introducing the capacitive coupling in Sec. III. As in the case of a single qubit, we are still able to find a closed-form expression for the Hamiltonian in the weak coupling limit.
a. Weak coupling limit The weak coupling limit applies as in the case of a single qubit, as explained in Sec. D 1 a: we assume C g (C α ) jj to get a closed-form expression for the interaction Hamiltonian, and C g C A , C B to remove any long-range coupling between the qubits and to recover the weak coupling limit of the system-bath interaction. With these prescriptions, if p are the momenta associated to the coordinates z, the circuit Hamiltonian now reads: where following the discussion for the single qubit scenario we have introduced the creation and annihiliation operators of each bath mode (respectively a † α and a α ). Note that, because of the assumption C g C A , C B , any direct interaction between the two qubits is removed. This is due to the fact that in the block (1, 1) of Eq. (D22) the matrix N, inducing the σ y A σ y B interaction, appears in a term of the order of O((C g /C A,B ) 2 ), while we are only keeping terms of the first order. This mimics the behaviour of the capacitive coupling in chains of LC circuits discussed in Sec. III. How can we re-introduce a direct qubit-qubit coupling without dropping the weak coupling approximation? The circuit scheme presented in Fig. 8 does the job: an additional coupling capacitor (weak or strong) directly connecting the qubits modifies the off-diagonal terms of the system matrix S in Eq. (D19). Since the inverse circuit matrix still has the form of Eq. (D21), we obtain the required interaction. Note that, in this scenario, even if the direct coupling between the qubit B and bath were switched off (i.e. a = (1, 0) T ), we would still obtain an interaction between them in the circuit Hamiltonian, through the vector s = S −1 a = (S −1 11 , S −1 12 ) T . A weak qubit-qubit coupling capacitor (C c C A,B ) would eliminate this direct interaction, since S −1 12 is of the first order in C c .

Two qubits, separate baths
Obtaining the general Hamiltonian in the separate baths scenario for a strong coupling is a cumbersome task, since it requires considering two sets of infinite internal variables of the resistors, and finding the inverse of a Lagrangian as a function of all these variables. This would lead to the presence of a direct interaction between all the elements of the circuit (and between the two sets of internal variables as well), greatly complicating the problem. We do not address the general case here, although a solution can still be obtained, by following the path presented in the previous sections and using the tools of Appendix E and Ref. [21]. Anyway, we can still find a closed-form expression for the Hamiltonian under the weak coupling approximation: considering the circuit in Fig. 9, the low capacitances C g and C c remove all the "long range" interactions, in analogy with the case of Sec. III, and we can obtain the desired Hamiltonian that is presented in Sec. IV C. Note that, unlike in the case of a direct coupling and a common bath, in this scenario we must consider a weak direct qubit-qubit coupling as well, otherwise the qubit B, for instance, may be directly coupled to the bath on the left in the interaction Hamiltonian. This coupling would be mediated by the vector s as in Eq. (D22), and it would vanish in the limit C c C A,B , as explained at the end of the weak coupling discussion in Appendix D 2.