Stochastic and Quantum Thermodynamics of Driven RLC Networks

We develop a general stochastic thermodynamics of RLC electrical networks built on top of a graph-theoretical representation of the dynamics commonly used by engineers. The network is: open, as it contains resistors and current and voltage sources, nonisothermal as resistors may be at different temperatures, and driven, as circuit elements may be subjected to external parametric driving. The proper description of the heat dissipated in each resistor requires care within the white noise idealization as it depends on the network topology. Our theory provides the basis to design circuits-based thermal machines, as we illustrate by designing a refrigerator using a simple driven circuit. We also derive exact results for the low temperature regime in which the quantum nature of the electrical noise must be taken into account. We do so using a semiclassical approach which can be shown to coincide with a fully quantum treatment of linear circuits for which canonical quantization is possible. We use it to generalize the Landauer-Buttiker formula for energy currents to arbitrary time-dependent driving protocols.


I. INTRODUCTION
Electronic circuits are versatile dynamical systems that can be designed and built with a high degree of precision in order to perform a great variety of tasks, from simple filtering and transmission of analog signals to the complex processing of digital information in modern computers. Traditional trends in the miniaturization of electronic components and in the increase of operation frequencies are nowadays facing serious challenges related to the production of heat and to the detrimental effects of thermal noise in the reliability of logical operations [1]. Thus, thermodynamical considerations are central in the search for new information processing technologies or improvements on the actual ones. In light of this, it might be surprising that a general thermodynamical description of electrical circuits is not available.
Perhaps one of the reasons for the absence of such general theory is the fact that a satisfactory understanding of thermodynamics and fluctuations in systems out of equilibrium has only been achieved in recent years. Stochastic thermodynamics [2,3] is now emerging as a comprehensive framework in which it is possible to describe and study thermodynamical processes arbitrarily away from thermal equilibrium, and to obtain different 'fluctuation theorems' clarifying and constraining the role of fluctuations. Simple electrical circuits have already been employed to experimentally study non-equilibrium processes and to confirm the validity of fluctuation theorems [4,5], however a general treatment is still lacking.
In this work we start by putting forward a general stochastic thermodynamic description of electrical circuits composed of resistors, capacitors and inductors, as well as current and voltage sources. The circuit components may parametrically depend on time due to an external controller. Our theory makes use of the graphtheoretic description of the network dynamics in terms of capacitor charges and inductor currents developed in electrical engineering. Noise is introduced via random voltage sources associated to resistors which may lie at different temperatures. For high temperatures, we use the standard Johnson-Nyquist noise which can be considered white. The first and second law of thermodynamics is formulated based on an underdamped Fokker-Planck description of the stochastic dynamics. A proper definition of local heat currents (the rates at which energy is dissipated in each resistor) turns out to require some care for reasons which depend on the topology of the network and which we relate to recent observations that simple nonisothermal overdamped models may be thermodynamically inconsistent [6][7][8][9].
After having established the high temperature theory, we proceed by demonstrating how it can be used to design thermodynamic machines made of electrical circuits. We do so by considering a simple electrical circuit where a resistor is cooled by driving two capacitors connected to it in a periodic manner. Such schemes may be employed to design new cooling strategies within electronic circuits. This is particularly interesting in the quantum domain of low temperatures, where parametrically driven circuits are commonly employed as low-noise amplifiers for the detection of small signals down to the regime of single quantum excitations [10,11]. In fact, the possible use of these circuits as cooling devices has been pointed out before [11,12]. Thus, we generalize our theory to the low temperature quantum regime, where the spectrum of the Johnson-Nyquist noise is not flat anymore and is given by the Planck distribution. Our approach does not involve the quantization of the degrees of freedom of the circuit, but can be shown to be equivalent to an exact quantum treatment of linear circuits that have a direct quantum analogue via canonical quantization. In this context we derive a generalization of the Landauer formula for transport in non-driven systems that is valid for arbitrary driving protocols. This is the central result of this article as it provides an efficient computational tool that can be applied to arbitrary circuits with any number of resistors at arbitrary temperatures, and subjected to arbitrary driving protocols. Finally, we show numerically that this formalism is able to capture the quantum limits for cooling recently identified in [13]. This article is organized as follows. In Section II we quickly review the basic graph-theoretical concepts involved in the description of electrical circuits. This will serve also to introduce notation and to define the basic objects to be employed later. In Section III we derive the deterministic equation of motion for a given circuit and analyze the energy balance and entropy production. In Section IV we expand the deterministic description to consider the noise associated to each resistive element. We provide an expression for the stochastic correction to the local heat currents, which carry on to the energy balance and the entropy production, and construct the Fokker-Planck equation describing the stochastic evolution of the circuit state and heat currents. In Section VI we show how to apply our formalism to study a simple example of an electrical heat pump. Finally, in Section VII we generalize our results to the low temperature quantum regime.

II. DESCRIPTION OF RLC CIRCUITS
We consider circuits composed of two-terminal devices connected with each other forming a network. A given circuit is mapped to a connected graph in which each two-terminal device is represented by an oriented edge or branch (see Figure 1). The orientation of each edge serves as a reference to indicate voltages drops and currents in the standard way. The state of the circuit is specified by the nodes voltages v 1 , · · · , v n and the edge currents j 1 , · · · , j b , where n is the number of nodes and b is the number of edges. We consider five types of devices: voltage sources, current sources, capacitors, inductors and resistors. The charges of all capacitors and the magnetic fluxes of all inductors are the dynamical variables of the circuit (alternatively, the voltages and currents, respectively). In order for these variables to be truly independent, we consider circuits fulfilling the following two conditions: (i) The circuit graph has no loops formed entirely by capacitors and voltage sources, and (ii) The circuit graph has no cut-sets formed entirely by inductors and current sources. If condition (i) is fulfilled, then the voltages of all capacitors and voltage sources are independent variables, i.e, they cannot be directly related via the Kirchhoff's voltage law (KVL). Analogously, if condition (ii) is fulfilled, then the currents of all inductors and current sources are independent variables since they are not directly constrained via the Kirchhoff's current law (KCL). Thus, it can be seen that under conditions (i) and (ii) it is always possible to find a tree of the graph (a fully connected subgraph with no loops) for which all the edges corresponding to capacitors and voltage sources are part of the tree, and all the edges corresponding to inductors and current sources are out of it 1 . We will call such a tree a normal tree and we will base our dynamical description of the circuit on it. The edges corresponding to resistors can be split in two groups according to whether or not they are part of the normal tree. Following the terminology of [14] we will refer to edges in the normal tree as twigs and edges outside it as links (also known as co-chords and chords, respectively). Every tree has b t = n − 1 twigs and therefore b l = b − (n − 1) links.

A. Loops and cut-sets matrices
Given a set of oriented loops (or cycles) we define the loop matrix B as follows: (1) A normal tree can be used to define a set of independent oriented loops in the following way: take the tree subgraph and add a link edge to it, then a loop will be formed (otherwise that link should have been a twig). The orientation of the loop is chosen to coincide with that of the added link. This can be done for every link. If we order the edges by counting first the twigs and then the links, the loop matrix thus obtained has the following structure: where 1 k is the k × k identity matrix. For the normal tree indicated in Figure 1-(b), and ordering the edges as {V, C 1 , C 2 , R 1 , R 2 , R 4 , R 3 , L, I}, the matrix B twig is: Given a set of oriented cut-sets (or cocycles) we define the cut-set matrix Q as: (4) As before, a normal tree can be used to define a set of independent cut-sets. The procedure is as follows. Take the tree subgraph and remove a twig, then the graph is split into two disconnected subgraphs. Consider all the edges going from one subgraph to the other, including the removed twig. These edges then form a cut-set which is oriented as the removed twig. This is then repeated for every twig, obtaining the following cut-set matrix: With the same ordering as before, the matrix Q link for the circuit of Figure 1 is: The matrices B and Q are orthogonal in the sense that BQ T = 0. This property follows naturally from their definitions and implies that B twig = −Q T link . We also mention that loops and cut-sets can be identified from the incidence matrix of the graph, and that the loop and cut-set matrices satisfy additional relations [15].

B. Kirchhoff 's Laws and Tellegen's theorem
If j and v are column vectors with the edge current and voltage drops as components, respectively, then the Kirchhoff's laws can be expressed in terms of the loop and cut-set matrices as follows: This form a set of b t + b l = b independent algebraic equations that can be used to eliminate half of the 2b variables contained in j and v. In particular, if we order the components of j and v such that twig variables appear first we can employ Eqs (2) and (5) and write: . We see then that it is enough to give the twigs voltages and the links currents to determine the rest of the variables.
From the orthogonality of B and Q and the two Kirchhoff's laws it is possible to prove Tellegen's theorem: any vector v compatible with the KVL and any vector j compatible with the KCL are orthogonal, i.e, j T v = 0.
C. Block structure of j, v and Qt In the following we will assume the previous splitting of the vectors j and v: Also, since all voltage sources (E) and capacitors (C) are twigs, while all current sources (I) and inductors (L) are links, we can further split the vectors j and v as follows: where x stands for j or v and R t and R l indicate the resistive edges which are respectively twigs or links. This partitioning induces the following block structure in the matrix Q link :

D. Constitutive relations
Each two-terminal device in the circuit is characterized by a particular relation between the electric potential difference across its terminals and the current through it. Voltage sources just fix a definite value for the potential difference, regardless of the current, and current sources fix a current value regardless of the voltage. Resistors are described by an algebraic relation between voltage and current. We can write, where R t and R l are diagonal matrices with the resistances of the twigs and links resistors as non-zero elements, respectively. R t and R l can be time-dependent. Also, for non-linear resistors like diodes these matrices can also be functions of the voltages or currents. Finally, capacitors and inductors are described by the following set of differential equations: Here, C and L are symmetric matrices describing the capacitances and inductances of the circuit. C is usually diagonal, but L could account for cross-couplings between different inductors. As with resistors, they can also depend on time or describe non-linearities.

III. DERIVATION OF DYNAMICAL EQUATIONS
We want to obtain a equation of motion for the circuit describing the evolution of the voltages of all capacitors and the currents in all inductors. We begin with the constitutive equations for them: The task now is to use the Kirchhoff's laws and the algebraic constitutive equations for the resistors to express the variables j C and v L in terms of the dynamical ones v C and j L . First, we use the following KCL and KVL equations which are part of Eq. (7): and we then obtain: where we have defined the matrices M c , M s and M d (the subindices stand for conservative, sources and dissipation, as will be justified in the following). We now only need to eliminate j R l and v Rt , since v E and j I are given. For this we use the constitutive relations of Eq. (12) and two additional Kirchhoff laws: which can be rewritten as: where we have defined the additional matrices α and M sd (in this case the subindex stands for source dissipation).
Inserting this relation in Eq. (16), a closed dynamical equation for the variables v C and j L is obtained. We can express it in the following concise form: where the matrix coefficients are given by: and Note that α might depend on time if resistances do. Some comments about the structure of the matrix A are in order. The first term in Eq. (20), M c , describes the conservative interchange of energy between capacitors and inductors. In some cases it can interpreted as analogous to the symplectic matrix in Hamiltonian mechanics. It has a block structure that stems from the separation of variables according to its behaviour under time reversal (voltages are even under time reversal while currents are odd). The second term takes into account the effect of the resistances in the dynamics. It is twofold: the antisymmetric part of M T d α(t)M d describes the interchange of energy between capacitors and inductors that is allowed by resistive channels, and the symmetric part describes the energy loss associated to it (see Eq. 27). Importantly, both the symmetric and antisymmetric part of M T d α(t)M d respect the block structure of M c . The meaning of this property is that energy dissipation is bound to be an invariant quantity upon time reversal, and has important consequences regarding the non-equilibrium thermodynamic behaviour of the circuit, as discussed in more detail in Appendix B.

A. Charge and flux variables
Instead of working with v C and j L as dynamical variables, from a physical point of view it is more natural to work with the charges in the capacitors and the magnetic fluxes in the inductors. They are defined as q = C v C and φ = L j L , respectively. We will group them in a column vector x. Thus, we can write where we have defined the matrix H. With these definitions, the dynamical state equation now reads: where s T = [v T E , j T I ] is a vector grouping the voltage and currents of the sources.

B. Linear energy storage elements
If we consider the particular case in which the matrices C and L, and thus also H, do not depend on the state vector x, we can express the energy contained in the circuit at a given time as a quadratic function of the circuit state: Then, we can write the dynamical state equation as: Note that Eq. (25) still allows for non-linear resistive relations, and in that case the matrices A and B have a tacit dependence on x (through α).

C. Energy balance
We now analyze the balance of energy between the different elements of the circuit and the entropy produced during its operation. We consider linear storage elements since in this case we have a simple notion of energy associated to a given state x of the circuit, which is given by Eq. (24). We begin by writing down the variation in time of the circuit energy: We see that it naturally splits into tree distinct terms that account for different mechanisms via which the energy stored in capacitors and inductors can change. The first one describes how this energy is dissipated into the resistors. Note that only the symmetric part of A plays a role in the expression ∇E T A∇E, and from the definition of A in Eq. (20) we see that it can only be different from zero if there are resistors in the circuit. Explicitly, where (X) s indicates the symmetric part of X. Secondly, voltage and current sources in the circuit can give energy to the capacitors and inductors, and this is described by the second term. Finally, external changes in the capacitances or inductances can also contribute to the circuit energy. This is clearly identified as work performed on the circuit by an external agent.
To complete the understanding of the energy balance of the circuit we analyze the total rate of energy dissipation in the resistors (i.e, Joule heating). Thus, the instantaneous dissipated power in each resistor isQ r = j r v r , and the total rate of heat production reads: (28) Using Eq. (18) to eliminate the resistor currents and voltages we find: As we will see next, the first term of the previous expression is exactly −Ė diss defined in Eq. (27), and thus represents the part of the energy dissipated into the resistors that is lost by capacitors and inductors. The second term corresponds to the part of the energy that is dissipated into the resistors directly by the voltage or current sources. The identity between −Ė diss and the first term of Eq. (29) can be established from the following property of the matrix α: that can be proven from the definition of α via 2×2 block matrix inversion. Using this, we can write the following expression for the energy balance of the circuit: where the total work rateẆ is the sum of the rates of work performed by the sources,Ẇ s , and by an external agent that drives the circuit parameters,Ẇ d . They are given by: As expected, using the fact that and Tellegen's theorem, we find that

D. Entropy production
At the level of description considered so far, the only aspect of entropy production that can be accounted for is the one related to the generation of heat in the resistors due to non-vanishing net currents. Thus, to every resistor r we associate a instantaneous entropy production equal toΣ r =Q r /T r , where T r is the temperature of the considered resistor. Therefore, we can write the total entropy production aṡ where β is a diagonal matrix with the inverse temperatures of the resistors as non-zero elements. As witḣ Q, Eq. (18) can be employed to expressΣ in terms of known quantities. However, we know in advance that the previous expression for entropy production is incomplete, since it doesn't take into account the effects of fluctuations. There are three main ways or mechanisms via which the fluctuations can play a role. In general, if there are fluctuations then the state of the circuit and its evolution become stochastic. Therefore an entropy can be associated to the probability distribution of the state x at any time, which will in general change and evolve in a non-trivial way. Of course, if the circuit is operating in a steady state this internal contribution to the total entropy production will vanish. But even in that situation fluctuations can still play a role, in combination with other two conditions. First, if the resistors are at different temperatures then fluctuations alone can transport heat from a hot resistor to a cold one [4], and this non-equilibrium process has an associated entropy production. Secondly, even if all the resistors are at the same temperature, heat transport can be induced by the driving of the circuit parameters [16]. For example, it is possible to devise a cooling cycle in which two capacitors are used as a working medium to extract heat from a resistor and dump it into another resistor (see Section VI).
In linear systems, this can be done only through fluctuations, without affecting the mean values of current and voltages and even in the case in which they vanish at all times. All these aspects of the full entropy production in linear circuits will be explored using stochastic thermodynamics in the next sections.

A. Johnson-Nyquist noise
In this section we extend the previous description of RLC circuits to take into account the effects of the Johnson-Nyquist noise [17,18] originating in the resistors. To do so we model a real resistor as a noiseless resistor connected in series with a random voltage source. The voltage ∆v of this source has zero mean, ∆v(t) = 0, and for high temperatures is modeled as delta correlated: where R is the resistance of the considered resistor, T is its temperature, and k b is the Boltzmann constant. This corresponds to a flat noise spectrum S(ω) = Rk b T /π. For low temperatures or high frequencies, one should in-stead consider the 'symmetric' [10] quantum noise spectrum: where N (ω) = 1 e ω/(k b T ) −1 is the Planck's distribution. The 1/2 term added to the Planck's distribution takes into account the environmental ground state fluctuations. It has no consequences regarding the thermodynamics of static circuits, even for non-equilibrium conditions. However, it does have important consequences in the ultralow temperature regime in the case of driven circuits, as discussed in detail in Section VII D. In the following we will just consider the classical case of high temperatures and delta correlated noise, where simpler results hold and the usual machinery of stochastic calculus can be employed. Then we present in Section VII the full treatment of quantum noise using Green's functions techniques.

B. Langevin dynamics
Thus, if we want to consider Johnson-Nyquist noise in our general description of circuits we need to introduce a random voltage source with the previous characteristics for every resistor in the given circuit. This can be done directly by considering each voltage source as a new element in the circuit graph, or indirectly by taking into account the presence of the new voltage sources in the Kirchhoff's laws of the original circuit. The latter approach has the advantage of not modifying the graph of the circuit and is the one we are going to use in the following. Thus, Eq.
and similar modifications to Eqs. (17) lead to the following generalization of Eq. (18): where ∆v R l and ∆v Rt are column vectors with the random voltage associated with links or twigs resistors, respectively. These modifications to the Kirchhoff's laws are propagated straightforwardly to the final equation of motion in Eq. (25), which now reads: where the vector η(t) groups the random voltages in the following way: (In what follows, for simplicity, we will omit the explicit dependence on time of A, B and α, which is anyway trivial if the resistances are constant). The previous expression for the equation of motion can be cast into a slightly more symmetric form. For this we note that since the statistics of the noise variables are invariant under sign inversion we can neglect the minus sign in front of ∆v R l . Also, for high temperatures, taking into account the form of the spectrum for each random voltage we can write: where R is the matrix defined in Eq. (28), β is the matrix of inverse temperatures defined in Eq. (35), and ξ(t) is a column vector of unit-variance white noise variables (as many as there are resistors). Then, the equation of motion for the circuit takes the following Langevin form: where A and B are the same as in the deterministic case. The last sum is over every resistor in the circuit, and C r is given by where Π r is a projector over the one-dimensional subspace corresponding to resistor r. Note that the matrices C r depend on time only if resistances do. Also, they always satisfy C r C T r = 0 for r = r . However the other possible product C T r C r does not vanish in general. Eq. (43) is the stochastic generalization of the usual state equation for electrical circuits.

C. Mean values and covariance matrix dynamics
For linear circuits, the mean values x of the charges and fluxes evolve according to the deterministic equation of motion: The evolution of the covariance matrix σ(t) = (x − x )(x − x ) T can be easily derived from the Langevin equation and reads: (46) It is important to note that for linear circuits the evolution of the covariance matrix is completely independent of the forcing function s(t) (if it is noiseless, as we considered). Thus, deterministic forcing of the circuit via voltage or current sources can only change the mean values of the charges and fluxes in the circuit but not the fluctuations around them, or their correlations.

A. Fluctuation-dissipation relation
We now consider the particular case in which the circuit is not driven (its parameters are time independent), all the resitors are at the same temperature T , and no voltages or currents are applied. Then, we expect the system to attain a thermal state for long times, with a distribution: This is a Gaussian state with covariance matrix It is instructive to check that in fact this result is obtained from the previous dynamical description of the circuit. Thus, if there is an asymptotic stationary state we see from Eq. (46) that in isothermal conditions its covariance matrix σ st must satisfy Also, using the definition of the matrices C r (Eq. (44)), we see that: where in the second equality we used the previously mentioned property αRα T = (α) s . Eq. (50) is nothing else that the fluctuation-dissipation (FD) relation for general RLC circuits. Using the FD relation we can easily verify that the covariance matrix σ th is in fact a solution of Eq. (49). We finally note that, due to the linearity of the circuit, if the circuit is forced by constant voltages and/or current sources, then the asymptotic state is given by a displaced thermal state: where the mean values x are given by the stationary solution of Eq. (45). However, p disp is not an equillibrium state since, at variance with p th , it has a non-vanishing entropy production given by Eq. (35).

B. Energy balance
To write down the energy balance for the stochastic description we begin by noticing that since the energy is a quadratic function its mean value can be easily expressed in terms of the mean values x and the covariance matrix σ: E(x, t) = E( x , t) + Tr[Hσ]/2. Therefore, we have (52) Thus, the energy balance takes the same form as before with the following expressions for the heat and work terms:˙ and˙ while the work performed by the sources is equal to its deterministic value (for noiseless sources): C. Local heat currents Using Eq. (46) for dσ/dt and the FD relation of Eq. (50) we can rewrite Eq. (54) for the total heat dissipation rate as: where for clarity we omitted the explicit time dependence of H. We see that each term in the sum of the previous equation can be associated to a particular resistor in the circuit. Then, they are a sensible definition for the local heat currents Q r (the rate of energy dissipation in resistor r): This heuristic definition of the local heat currents is often considered in the literature (see, for example, [19]), but it is not always correct. In fact, if we add to the quantities Q r defined in the previous equation a term of the form r ∆Q r,r , for any antisymmetric tensor ∆Q r,r , we find that the total heat rate Q = r Q r remains unchanged. Then, it is in general not possible to derive local heat currents via a decomposition of the global one.
From another perspective, in the context of electrical circuits the local heat currents can be naturally defined as where ∆v r is the random voltage associated to each resistor. That quantity, however, is found to be divergent in the general case, as we explain in detail in the next section. The physical origin for this divergence is the fact that for some circuits thermal fluctuations of arbitrarily high frequencies originating in a resistor (which are always present in the model due to the white-noise idealization) can be dissipated into another, as illustrated in Figure 2 with a simple example. Another way to think about this problem is to consider that the initial description of the circuit, which might be a valid one with respect to the dynamics, is missing inertial degrees of freedom relevant for the thermodynamics. This was recently discussed in [6][7][8][9], although not in the specific context of electrical circuits. In particular, in [9] it was shown that inertial degrees of freedom transfer heat at a rate that is inversely proportional to the their relaxation timescale, which is taken to be infinitesimal in the overdamped limit. This is an alternative way of understanding the divergence of the heat currents.
In the next section it is shown that the quantities Q r as defined in the last equation are always well behaved if and only if, given a normal tree of the graph of the circuit, there are no fundamental cut-sets associated to it containing link resistors and twig resistors simultaneously (i.e, if Q RR = 0 in Eq. (11)). In that case, Eqs. (58) and (59) coincide.

D. Fokker-Planck equation for the circuit state and heat currents
Each resistor was modeled as an ideal resistor in series with a random voltage noise. If j r is the instantaneous current flowing through the resistor, then the rate of energy dissipation in it iṡ where v r = R r j r is the voltage drop in the ideal resistor and ∆v r is the random voltage. We want to express the quantityQ r in terms of the state x of this circuit. For this we first write it aṡ where Π r is the projector associated with the r-th resistor appearing in Eq. (44), and R is the matrix of resistances defined in Eq. (28). We can now use Eq. (39) to eliminate the variables j R l and v Rt . For ease of notation, in the following we will not consider the terms associated to the sources of the circuit, as we know that they can only affect the mean values of the current and voltages and we are only interested in the stochastic contributions to the heat currents. In this way, after some manipulations, we obtain the following expression: Notice that the last term in this expression is quadratic in the noise variables. As a consequence, it will give a divergent contribution to Q r in the white noise limit, proportional to δ(0). As we already mentioned, the physical origin of this divergence is the possibility of direct heat transport at arbitrarily high frequencies between resistors at different temperatures. We stress that this problem only arises with regard to the definition of local heat currents, while the state of the circuit x and the total heat currentQ are always well behaved quantities (see below). A solution to this problem would be to give a more realistic description of the resistive thermal noise, associating to each resistor a spectral density J r (ω) vanishing for large frequencies, such that its noise spectrum is given by S r (ω) = (Rk b T r /π)J r (ω). However, this is equivalent to appropriately 'dressing' a white-noise resistor with inductors and/or capacitors that can be considered part of the circuit 2 (a capacitor in parallel or a inductor in series to a given resistor being the most simple options to filter out high frequencies, see Figure 2). This observation hints at a relationship between the topology of the circuit and the possibility of defining well behaved local heat currents in the white noise limit. In fact, we see that the quadratic terms in the noise vanish if the matrix α t Π r Rα is block diagonal. In turn, from the definition of α we can see that this happens if Q RR = 0. Thus, Q RR = 0 is a sufficient condition for the local heat currentsQ r to be well behaved. In Appendix A it is shown that this condition is also necessary. To finish this discussion we note that since r Π r = 1 and α t Rα is always block diagonal (recall Eq. (30)), we see that the total heat rateQ = rQ r is always well behaved in the white noise limit, even if Q RR = 0. Thus, assuming Q RR = 0, the local heat currents are simplified to: where η is the vector of random voltages and currents defined in Eq. (42) and we used the fact that for Q RR = 0 we have α = R −1 . This equation is a Langevin equation for the heat Q r that is of course coupled to the Langevin equation in Eq. (43) for the circuit state x. Their integration is to be performed according to the Stratonovich procedure. As explained in [20] (Chapter 8), in the white noise limit the corresponding stochastic dynamics can be described by the following set of Ito differential equations: and where dW is a vector of independent Wienner processes differentials. Taking the mean value of Eq. (65) we recover Eq. (58) for Q r . The Fokker-Planck equation for the joint probability distribution P (x, Q r ) corresponding to the previous Ito differential equations reads: where ∇ is the nabla operator with respect to the variables x. The corresponding equation for the reduced probability distribution p(x) = dQ r P (x, Q r ) is just Equation (66) allows to analyze the full statistics of the heat currents. Indeed, different integrated and detailed fluctuation theorems can be derived for this kind of linear and in general underdamped stochastic systems [9,21,22], valid for finite time protocols or asymptotic steady states. However, in this article we will focus only on the behaviour of the mean values.

E. Entropy production
We consider the continuous Shannon entropy associated to the distribution p(x) and the total entropy production rateΣ Using Eq. (67) it is possible to show thatΣ accepts the following explicitly positive expression: , t)). This establish the validity of the second law of thermodynamics, but also offers a way to compute the full entropy production (and thus, to characterize the irreversibility) of arbitrary non-equilibrium processes. This total entropy production can be decomposed as a sum of adiabatic and non-adiabatic contributions [23][24][25][26], as shown in Appendix E. The adiabatic contribution is positive definite and for time independent circuits it is the only non vanishing contribution for large times. For overdamped circuits (for example, circuits with no capacitors or no inductors) the non-adiabatic contribution is also positive definite, and for time independent circuits it equals −k b times the time derivative of the relative entropy H(p|p st ) between the instantaneous state p(x, t) and the stationary one p st (t) (which for linear circuits is unique, although it might depend on the initial conditions). Thus, for time independent circuits H(p|p st ) is a always decreasing Lyapunov function. For underdamped circuits a third non-adiabatic term appears which, at variance with the previous two, is not positive definite. It is related to the change in H(p|p st ) due to a conservative flow in phase space, and vanish identically in isothermal conditions (when p st is an equilibrium state). These findings are analogous to the results of [26].

VI. A SIMPLE CIRCUIT-BASED MACHINE
We now illustrate how our formalism can be used to design thermodynamic machines made of RLC circuits. External driving on the circuit allows to implement thermodynamical cycles that might extract work from a thermal gradient (non-autonomous heat engine) or to extract heat from some resistors (non-autonomous refrigerator). We illustrate the basic techniques by considering a minimal circuit which can work both as an engine or a refrigerator. The circuit is shown in Figure 3 and consists of two parallel RC circuits coupled by an inductor. The capacitances in each RC circuit can be driven externally, for example by changing in time the distances between the plates of each capacitor, or, more practically, using varicap diodes. The circuit has no loops consisting only of capacitors or cutsets of all inductors, and therefore the previous formalism can be directly applied. We first analyze the simplest case of regular heat conduction for constant parameters and different temperatures. Then we show that the capacitances in the circuit can be driven in time in order to cool (i.e, extract heat) from one of the resistors, while dumping the extracted energy into the other one. Similar circuits were analyzed before in [4,27].
We begin by describing the circuit by the procedure of Section II. The state of the circuit is encoded in the vector x = (q 1 , q 2 , φ) T , where q i is the charge in the capacitor C i and φ is the magnetic flux in L. According to the edge orientations of Fig. 3-(b), the cutset matrix associated to the normal tree is specified by and from this we can construct the matrix M c and M d : Also, in this case α = R −1 = diag(1/R 1 , 1/R 2 ) and therefore: Finally, we have H = diag(C 1 , C 2 , L) −1 , C 1 C T 1 = diag(R −1 1 , 0, 0) and C 2 C T 2 = diag(0, R −1 2 , 0).

A. Heat conduction
Given the above matrices, we can readily solve Eq. (49) to find the stationary covariance matrix (for time independent parameters). The solution is particularly simple in the symmetric case in which R 1 = R 2 = R and C 1 = C 2 = C. It reads whereT = (T 1 +T 2 )/2 and ∆T = T 1 −T 2 . Thus, the first term in the previous expression is just the equilibrium covariance matrix corresponding to the mean tempera-tureT . By examining the second term, we see that a temperature bias will establish correlations between the capacitors and the inductor, but not between the capacitors themselves. Since this circuit is not forced by voltage or current sources, the mean values of voltages and currents in any branch will vanish in the stationary state. Then, introducing the previous expression for the covariance matrix in Eq. (58) for the rate of heat dissipation in each resistor, we obtain: The entropy production can be computed from the above heat rates and reads:

B. Cooling cycle
We now turn to analyze the more interesting situation in which the two resistors are at the same temperature T 1 = T 2 = T and the capacitors are driven periodically in time such that C 1 (t) = C +∆C cos(ω d t) and C 2 (t) = C + ∆C cos(ω d t+θ). Thus, both capacitances are driven with the same angular frequency ω d and the same amplitude ∆C but with some fixed phase difference θ ∈ (−π, π). The matrices describing the circuit are the same as before except for the energy matrix H(t) that now depends on time. It is periodic so it can be decomposed as a Fourier series: To lower order in ∆C we have only three Fourier components: H 0 = diag(C, C, L) −1 and H ±1 = −∆C/(2C 2 ) diag(1, e ±iθ , 0). We are interested in settings in which there is a stable stationary state. We note that this is not always the case due to the phenomenon of parametric resonance. However, if there is a stable stationary state it will be such that the mean values of voltages and currents vanish, while the covariance matrix is periodic with the same period as the driving. Thus, we can decompose it as Inserting the previous two Fourier decompositions into Eq. (46) we obtain an algebraic equation from which it is possible to obtain the coefficients σ 2 k,k in terms of H 0 and H ±1 [28]. This technique is explained in Appendix D. Once the Fourier components σ 2 k,k have been determined, we compute the average heat and work rates per cycle. Explicitly, we consider the quantity where X stands for Q 1 , Q 2 , W or E. We note that Ė c = 0 since the asymptotic state of the system is periodic. Thus, averaging the balance of energy (Eq. (31)) during a cycle we obtain where Ẇ c is the average rate of work corresponding to the external driving, which is the only source of work in this case. As shown in Appendix D, to lower order in ∆C and to second order in ω d , Ẇ c is given by while the average heat currents are We then see that the rate of heat pumping from one resistor to the other (the first term in Eq. 82) is proportional to ω d , while the rate at which work performed by the driving (or equivalently, the total dissipated heat rate) scales as ω 2 d . Also, the pumping of heat is maximized and the dissipated work minimized for θ = ±π/2. This is natural since in this case the left/right asymmetry induced by the driving is maximum. The cooling efficiency or 'Coefficient of Performance' is which for θ = π/2 and in the limit L → 0 is just There is a maximum driving frequency such that cooling is not possible above it (in the considered regime of low ω d and ∆C). It corresponds to CoP = 0 and for θ ∈ [0, π] reads There is also an optimal frequency ω opt d , in the sense that the heat extracted from one of the resistors is maximized. We can obtain it by optimizing Q 1 c in Eq. (82) with respect to ω d , and in this way we find that ω opt If the resistor temperatures are different the heat currents are In contrast to the isothermal case, the term of second order in ω 2 d is too involved to show it here. The first term corresponds to regular heat conduction in response to the thermal gradient. The second term is a correction to the regular heat conduction, due to the driving, while the third term describes the pumping of heat. We consider the case in which T 1 < T 2 (then, ∆T = T 1 − T 2 < 0) and analyze the conditions under which it is possible to extract heat from R 1 . The pumping of heat out of R 1 is, as before, optimized for θ = π/2. From the previous equation we see that in general the driving frequency must be above a minimum value in order for the heat pumping to overcome the heat conduction imposed by the thermal gradient. Thus, we will have effective cooling of R 1 only if ω d > ω min d . For θ = π/2 this minimum cooling frequency ω min and we can write the heat rate Q 1 as: The previous considerations do not take into account the terms of second order in ω d . From the expression of the heat currents in the isothermal case, Eq. (82), we know that this corrections correspond to heating and establish a maximum driving frequency ω max d such that cooling is not possible above it, Eq. (85). Thus, for cooling to be possible at all we need that ω min d < ω max d , which imposes a condition on the temperature difference.
We now turn to analyze the total heat rate, or work rate. Up to second order in ω d it is given by the following expression: Note that that if ∆T = 0 then to lower order in ω d the average work rate can be positive or negative, depending on the value of θ. This two cases correspond to the device working as a refrigerator or a (non-autonomous) heat engine, respectively.

D. Exact numerical results
The previous analytical results for the cooling protocol are limited to low driving amplitude and frequency. In order to assess their validity in that regime and to study the behaviour of the system away from it, we numerically compute the heat currents. For this we integrate the differential equation for the time evolution of the covariance matrix (Eq. (46)). Then we compute the instantaneous expected values for the heat currents via Eq. (58), and obtain their averages during a cycle for sufficiently long times. The circuit has two characteristic time scales: τ 0 = √ LC and τ d = RC, respectively related to the oscillations period and the dissipation rate. For the numerical evaluation we consider τ d = τ 0 and take this quantity as the unit of time. As an example we show in Figure 4-(a) the long time oscillations of the heat currents, as well as their averages, for an isothermal setting and the following driving parameters: ∆C/C = 1/2, ω d = 10 −2 2π/τ d and the optimal phase difference of θ = π/2.
The analytical and numerical results are compared in Figure 4-(b) for a fixed driving amplitude (∆C/C = 0.5) and increasing driving frequency, while the temperatures are the same and the phase difference is the optimal. We see that the analytical expressions indeed match the numerical results in the low driving frequency regime. We also see that there is, as expected from the theoretical analysis, a maximum driving frequency ω max Analogously, we show in Figure 4-(c) the heat currents for fixed driving frequency (ω d /(2π) = 10 −2 /τ d ) and increasing driving amplitude. Again, we see that for low driving amplitude the analytical expressions correctly describe the numerical results.

VII. QUANTUM JOHNSON-NYQUIST NOISE
If the typical temperature in the circuit is low enough that the thermal energy k b T starts to be comparable to the quantum of energy ω at the relevant frequencies ω, then the quantum nature of the noise in each resistor must be taken into account. One approach to work in this regime is to construct a quantum model of the circuit, assigning quantum mechanical operators to the charge and flux degrees of freedom in it (or a proper combination of them), satisfying the usual commutation relations. The dissipation and diffusion effects induced by the resistors are typically introduced using the Caldeira-Legget model for quantum Brownian motion [29]. This is certainly the way to go if one is interested in having access to the full quantum state of the circuit (for modern treatments, see [10,[30][31][32][33][34]). However, this cannot be directly done for any circuit, since the canonical quantization procedure requires the detailed specification of stray or parasitic capacitances and inductances [30]. As an example we can take the circuit of Fig. 2-(c). This circuit has well behaved dynamics and heat currents, but cannot be directly quantized since it is missing inertial degrees of freedom. In other words, it is not possible to define canonical 'momentum' coordinates that are conjugate of the capacitor charges (there are no kinetic energy terms [32]). To do that one must give a more detailed description specifying stray inductances.
Thus, we would like to have a tool to study the lowtemperature behaviour of circuits like the one in Fig. 2-(c) in a direct way, without having to worry about the explicit description of irrelevant degrees of freedom. For this we put forward a semiclassical approach that is based on the same stochastic equations of motion of Eq. (43), where now the noise variables ξ(t) are not white anymore and display a quantum spectrum. Then, if ξ r is the adimensional noise variable associated with the r-th resistor we must consider the following power spectrum: where we use the shorthand N r (ω) = (e ω/(k b Tr) − 1) −1 for the Planck's distribution at temperature T r . Inverting the previous equation we can compute the correlation functions: (91) where Λ is a high-frequency cutoff, that must be large compared to any other frequency scale of the problem.
The variables x specifying the circuit state remain classical (they are not promoted to quantum mechanical operators). In spite of this, it is possible to show that for linear circuits the results obtained in this way are fully equivalent to those obtained by a full quantization, if such a quantization is possible.

A. Covariance matrix and heat currents for quantum noise
Since ξ (t) = 0, for linear circuits the equation of motion for the mean values is still the fully deterministic one given in Eq. (45). However, in the quantum low temperature regime the differential equation for the covariance matrix, Eq. (46, must be modified. The reason is that the circuit state x(t) at time t will be in general correlated to ξ(t) (i.e, it will not be a non-anticipating function, and therefore the usual assumptions of stochastic calculus underlying the derivation of Eq. (46) are not valid [20]). For stable systems it is still possible to obtain a simple expression for the covariance matrix at large times. To see this explicitly it is convenient to employ techniques based on the Green's function of the circuit, like it is done in fully quantum mechanical models. We start with the equation of motion for y = x − x : Given the initial value y(0), the solution to this equation can be written as where the retarded Green's function G(t, t ) is defined as the solution of with G(t, t ) = 0 for t < t (from this it follows that G(t , t ) = 1). Then, the covariance matrix can be expressed as The first term in this expression is just the deterministic evolution of the fluctuations present in the initial state. The second term takes into account the effect of the correlations between the circuit initial state and the environmental noise. The last term, which for stable systems dominates the long time behaviour, represents the diffusion induced by the environment. We will assume in the following that the initial state is not correlated in any way with the environmental noise, y(0)ξ(τ ) T = 0, and therefore the second term in the previous equation vanish. We will also assume, for simplicity, that the resistances, and thus the matrices A and C r , are constant. Then, taking the time derivative of Eq. (95) and using Eq. (94), we obtain the following differential equation: where I r (t) is the convolution between the Green's function G(t, t ) and the correlation function of resistor r: Eq. (96) is the generalization for quantum noise of Eq. (46), which is recovered in the limit of high temperatures. To see this, we note that for high temperatures ξ r (0)ξ r (τ ) → δ(τ ), and therefore I r (t) → G(t, t)/2 = 1/2. Based on these results we can now derive a expression for the local heat currents that, in contrast with Eq.
(58), is exact and valid for arbitrary temperatures. In the quantum case the total heat rate is also given by Eq.
However this time we should replace dσ/dt by Eq. (96). Using this and the FD relation we can write: Q = r j r v r + Tr (Hσ(t)H − 2k b T r HI r (t))C r C T r .
(98) In analogy with the classical case, under the condition Q RR = 0, we can identify the local heat currents as: (99) This expression can be already evaluated given the above equations for σ and I r . However, if we are only interested in the asymptotic heat currents under the assumption that the dynamics of the system is stable, we can express them as frequency integrals that might be easier to compute, and that also have a clear physical interpretation in terms of elementary transport processes.

B. Asymptotic covariance matrix and heat currents
If the system is asymptotically stable, i.e, if G(t, t ) → 0 for |t−t | → ∞, then the first two terms in Eq. (95) can be neglected for sufficiently long times. By expressing the correlations ξ(τ )ξ(τ ) T in terms of the power spectrum S r,r (ω) via inversion of Eq. (90), we can rewrite the last term in Eq. (95) as: where we have defined the following partial transform of the Green's function: For circuits with periodically driven parameters,Ĝ(t, ω) has the useful property of being asymptotically periodic in time with the same period as the driving, as shown in Appendix C. It also trivially satisfiesĜ(t, ω) * = G(t, −ω), a property that is sometimes used implicitly in the derivations below. The convolution integral I r (t) can also be expressed in terms ofĜ(t, ω): Finally, we note thatĜ(t, ω) can be directly obtained by solving its own evolution equation, that can be derived from Eq. (94) and reads with the initial conditionĜ(t = 0, ω) = 0. Introducing Eqs. (100) and (102) for σ and I r into Eq. (99), we can write the local heat currents as: where we introduced the shorthand definition D r = C r C T r , that we will employ in the following to simplify notation. The first term inside the trace, that is quadratic inĜ(t, ω), is actually closely related to the second one, which is linear inĜ(t, ω). We can see this by employing Eq. (103) to compute the derivative ofĜ † HĜ: (105) Using this relationship and Eq. (50) it is possible to rewrite Eq. (104) in the following compact way: dω ω f r,r (t, ω) (N r (ω) + 1/2) (106) where f r,r (t, ω) is a transfer function, specifying how the temperature of resistor r affects the heat current of resistor r. For r = r is always positive and is given by: (107) while the diagonal elements f r,r (t, ω) are determined by the following expression for the sum over the first index: Equation (106) is the central result of this article. It is a fully general expression for the local heat currents that can be applied to arbitrary RLC circuits, with any number of resistors at arbitrary temperatures, and is valid for arbitrary driving protocols.
To clarify the physical interpretation of the previous expressions we analyze first the particular case of time independent circuits.

C. Undriven circuits
If the matrix H is time independent then for sufficiently long times we have, from Eq. (103),Ĝ(t, ω) = G 0 (ω) = (iω − AH) −1 (this is just the Laplace's transform of the Green's function evaluated at iω). Then, asymptotically, the transfer functions f r,r (ω) are time independent andf r (ω) = 0, from where it follows that f r,r (ω) = − r =r f r ,r (ω). Also, in this case f r,r (ω) is symmetric under interchange of the indexes r and r . To see this it is necessary to consider the block structure of the matrix A, that is inherited by the matrix G 0 (ω), and of the matrices D r . This is discussed in detail in Appendix B, in connection with the invariant nature of the dissipation upon time inversion. Thus, using the mentioned properties we can obtain the usual Landauer-Büttiker expression for the heat currents [35][36][37]: From this equation, the quantity f r,r (ω)dω can be naturally interpreted as the rate at which an excitation with frequency between ω and ω + dω is transported from resistor r to r. We note that the symmetry of the transfer function in the undriven case makes the heat currents to depend only on the differences (N r (ω) + 1/2) − (N r (ω) + 1/2), and the 1/2 terms added to each Planck's distribution cancel each other. However, this is not the case for driven circuits, where f r,r (ω) is not symmetric in general. In that case, the ground state fluctuations represented by the 1/2 are responsible for the dissipation of heat into the resistors due to the parametric driving even if all the temperatures are zero.
Thus, Eq. (106) can be considered the generalization to arbitrary driving protocols of the Landauer-Büttiker formula for the static case, Eq. (109). In the following section we provide simplified expressions for the transfer functions in the case that the external driving is periodic, which are useful for the numerical study of thermal cycles.

D. Periodically driven circuits
We now consider the situation in which the matrix H(t) is a periodic function of time, and therefore can be decomposed as a Fourier series: where ω d is the angular frequency of the driving. In this case, assuming stable dynamics and long times, the functionĜ(t, ω) is also periodic with the same period of the driving, as shown in Appendix C. Thus, the following decomposition holds asymptoticallŷ Then, the asymptotic covariance matrix of the system and the heat currents are also periodic, with period τ = 2π/ω d . We thus consider the average values of the heat currents during a driving period, that we denote Q r c , as we did in the example of Section VI. They are given by where F r,r (ω) is the asymptotic average of f r,r (t, ω) during a driving period, that can be expressed in terms of the Fourier components H k andĜ j (ω): for r = r, and F r (t, ω) = r F r,r (t, ω) Finally, we note that given the Fourier components H k of the external driving, the Fourier componentsĜ j (ω) of the Green's function can be obtained by solving the following infinite set of algebraic equations that is obtained by introducing the decompositions of Eqs. (110) and (111) into Eq. (103). Some methods to solve this equation are discussed in Appendix C. The interpretation of the previous expressions in terms of elementary transport processes is not as straightforward as in the regular Landauer-Büttiker formula for the undriven case. For open mechanical systems composed of quantum harmonic oscillators, a physically clear decomposition of the local heat currents in terms of assisted transport and pair creation of excitations was obtained recently [13]. In particular, the pair creation mechanism was shown to be dominant at low temperatures and to be responsible for the ultimate limit for cooling in those systems. In the next section, this quantum limit for cooling is illustrated numerically for the cooling scheme introduced in Section VI.

E. Quantum limits for cooling
In this section we show how the quantum corrections to the heat currents impose a minimum temperature below which it is not possible to extract heat. We consider the particular cooling scheme of Section VI in isothermal conditions. Thus, both resistors are at the same temperature T , and we choose a driving amplitude ∆C, frequency ω d , and phase difference θ such that for high temperatures heat is extracted from resistor R 1 (thus, Q 1 c < 0), while it is dumped in resistor R 2 ( Q 2 c > 0). In Figure 5 we compare the classical (i.e, high temperature) heat currents, obtained by averaging Eq. (58) over a driving period, with the quantum heat currents according to Eq. (112). The parameters are ∆C/C = 10 −1 , ω d = 10 −2 (2π/τ 0 ) and τ d = 2τ 0 (recall that τ 0 = √ LC and τ d = RC are the oscillation and dissipation time scales, respectively). As expected, the quantum heat currents approach the classical ones for increasing temperature. However, below a given value of T (indicated as T * in Fig. 5), Q 1 c becomes positive and cooling stops. This is shown in more detail in Figure 6, where it is also clear that the value of T * decreases with increasing τ d /τ 0 , or equivalently with decreasing dissipation rate. This breakdown of the cooling effect is a strong coupling result that cannot be captured with usual approaches based on master equations (see [27], for example), as discussed in detail in [13].

VIII. CONCLUSIONS
We presented a general study of the non-equilibrium thermodynamics of driven electrical circuits. We derived the stochastic evolution of the circuit state and of the heat currents dissipated in each resistor. A relation between the topology of the circuit and the possibility of defining finite heat currents under the white noise idealization was established. As a first and simple example of application, we showed how to use our formalism to study the transport and pumping of heat in a minimal circuit of two driven RC circuits coupled by an inductor.
The initial classical treatment was then generalized in order to consider the effects of quantum low-temperature noise. We considered a semiclassical treatment in which the classical equations of motion are driven by noise with a quantum spectrum. In contrast with treatments based on the full quantization of the degrees of freedom in the circuit, our method has the advantage of being directly applicable to circuits that cannot be quantized without the additional specification of stray inductances or capacitances, but are however detailed enough to properly describe the dynamics and also the thermodynamics. Based on these results we expressed the heat currents for static circuits in terms of the familiar Landauer-Büttiker formula, and also obtained the generalization of this expression for arbitrary driving protocols.
Our results offer a general formalism to study and design thermodynamical processes in electrical systems from a first principles perspective and working in strongly non-equilibrium conditions, and also away from the adiabatic and weak coupling regimes. A direct application of the expressions provided and of the techniques illustrated in this article is the automatic optimization of thermal cycles in complex and large electrical circuits. This is particularly straightforward in the regime of high temperatures, where optimal cycles can be obtained that could be later refined to take into account quantum effects.
Finally, nontrivial networks with stochastic linear dynamics are commonly used to describe various kinds of complex systems including biological ones [38,39]. In principle any such network can be emulated by a suitable RLC circuit. This means that not only our results may apply to a very broad class of systems, but also that experimental studies of those models could be carried out experimentally using RLC circuit.
where T l and T r are diagonal matrices with the temperatures of the link and twig resistors, respectively. For the previous expression to vanish for arbitrary temperatures, the matrices B r Q RR and Q T RR B r should have null diagonals. The explicit of B r is B r = (R l + Q T RR R t Q RR ) −1 Q T RR π r (R −1 t + Q RR R −1 l Q T RR ) −1 if the index r correspond to a twig resistor, or B r = −(R l + Q T RR R t Q RR ) −1 π r Q T RR (R −1 t + Q RR R −1 l Q T RR ) −1 if r correspond to a link resistor, where π r is the reduction of the projector Π r to the appropriate twigs or links subspace. In any case it is easy to see that for the matrix B r Q RR to have null diagonal for any r and for arbitrary values of the resistances, all the components of the matrix Q RR should vanish.
with G(t, t ) = 0 for t < t . If H(t) is a periodic function with period τ d = 2π/ω d , and the solution to the previous differential equation is unique, then we have that G(t, t ) = G(t + τ d , t + τ d ). Then: where in the first step we just employed a change of variables (τ → τ − τ d ), and in the last step we assumed that the system is stable, in the sense that G(t, t ) → 0 for |t − t | → ∞. If that condition holds then for sufficiently large t we can neglect the contribution of the first part of the integration domain. Thus, we have shown that under this condition the function G(t, ω) is asymptotically periodic, with period τ d . We note that the stability condition does not always hold, even in the presence of strong dissipation, since it is possible for the circuit to continuously absorb energy from the driving and have a divergent dynamics. This is the phenomenon of parametric resonance, that we exclude from our analysis.
Therefore, for long times t we can give the following Fourier decomposition of the functionĜ(t, ω): G(t, ω) = k=+∞ k=−∞Ĝ j (ω)e ijω d t . (C4) Then, inverting Eq. (C1) and using Eq. (C2) (or equivalently, transforming Eq. (C2) to obtain Eq. (103) in the main text), we can find the following set of algebraic equations for the coefficientsĜ j (ω): i(ω + jω d )Ĝ j (ω) = 1δ j,0 + A k=+∞ k=−∞ H kĜj−k (ω), (C5) A simple method to solve these equations is to use a perturbative approach in which the strength of the driving is considered small, i.e, we assume |H k | |H 0 | for all k = 0. Then, to first order in H k =0 we have G 0 (ω) (iω1 − AH 0 ) −1 , that is just the transform of the Green's function on the undriven circuit, and: G j (ω) Ĝ 0 (ω + jω d )AH jĜ0 (ω) for j = 0. (C6) Higher orders in H k =0 can be easily computed. From this solution we see that for weak driving the range of relevant Fourier components inĜ(t, ω) is restricted by that of H(t). Thus, another non-perturbative method to solve Eq. (C5) is just to truncate the Fourier space to some maximum number of components given by |k| ≤ k max and |j| ≤ j max and then numerically solve the resulting finite system of linear equations for each value of ω.

Appendix D: Generalized Lyapunov Equation
In this section we introduce a generalization of the Lyapunov equation that is useful to compute the asymptotic state of periodically driven linear systems subjected to white noise. For undriven circuits, the covariance matrix σ for large times can be obtained as the solution of the following Lyapunov equation: For driven circuits there is not a time independent asymptotic state and one must solve the dynamical equation: (D2) However, if the function H(t) is periodic, then we known from the results of the previous section that the system state for large times will also be periodic (with the same period of H). Then we consider H(t) = k=+∞ k=−∞ H k e ikω d t and introduce the following decomposition for σ(t): Thus, the problem is now to find the coefficients σ k,k in terms of H k . Note that, at variance with a regular Fourier decomposition, the previous expression involves a double summation and as a consequence the coefficients σ k,k are not uniquely defined. This choice, however, allows to cast our problem as an extended Lyapunov equation. Indeed, it is useful to introduce the following definitions: Then, it can be seen that the coefficients σ k,k that give a solution to Eq. (D2) can be obtained by solving the following generalized Lyapunov equation: Of course, to numerically solve this problem we need to truncate the dimensions of the matrices A and S. As we saw in the previous section, this is justified for sufficiently weak driving. In the example given in Section VI the function H has only 3 Fourier components, and therefore to lower order in H ±1 we can truncate the matrix S to 3 blocks in each direction. After doing this we solved symbolically the Lyapunov equation of Eq. D7 using Mathematica.