Stochastic and Quantum Thermodynamics of Driven RLC Networks

TL;DR

This paper develops a comprehensive stochastic thermodynamics framework for driven RLC networks, incorporating topology, nonisothermal conditions, and quantum effects, enabling the design of circuit-based thermal machines.

Contribution

It introduces a novel thermodynamics theory for complex RLC networks, including quantum effects and exact energy current formulas under time-dependent driving.

Findings

Designed a circuit-based refrigerator as an example application.

Derived exact low-temperature quantum results for electrical noise.

Generalized the Landauer-Buttiker formula for energy currents with driving protocols.

Abstract

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.

Equations

$$B_{i,j}=\{\begin{array}{l}1\text{ if edge j is in loop i with = orientation}\\ -1\text{ if edge j is in loop i with = orientation }\\ 0\text{ otherwise}\end{array}$$

$$B=\left[\begin{array}[]{cc}B_{\text{twig}}&\mathds{1}_{b_{l}}\\ \end{array}\right],$$

$$B_{\text{twig}}=\begin{array}[]{c|cccccc}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}V}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}C_{1}}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}C_{2}}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}R_{1}}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}R_{2}}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}R_{4}}\\ \hline\cr R_{3}&-1&0&1&0&0&-1\\ L&0&1&0&1&0&-1\\ I&1&1&-1&1&1&0\end{array},$$

$$Q_{i,j}=\{\begin{array}{l}1\text{ if edge j is in cut-set i with = orientation}\\ -1\text{ if edge j is in cut-set i with = orientation }\\ 0\text{ otherwise}\end{array}$$

$$Q=\left[\begin{array}[]{cc}\mathds{1}_{b_{t}}&Q_{\text{link}}\\ \end{array}\right].$$

$$Q_{\text{link}}=\begin{array}[]{c|ccc}&R_{3}&L&I\\ \hline\cr{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}V}&1&0&-1\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}C_{1}}&0&-1&-1\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}C_{2}}&-1&0&1\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}R_{1}}&0&-1&-1\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}R_{2}}&0&0&-1\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}R_{4}}&1&1&0\end{array},$$

$$\text{KVL:}Bv$$

$$\text{KCL:}Qj$$

$$v_l$$

$$j_t$$

$$j=\left[\begin{array}[]{c}j_{t}\\ j_{l}\end{array}\right]\qquad v=\left[\begin{array}[]{c}v_{t}\\ v_{l}\end{array}\right].$$

$$x_{t}=\left[\begin{array}[]{c}x_{E}\\ x_{C}\\ x_{R_{t}}\end{array}\right]\qquad x_{l}=\left[\begin{array}[]{c}x_{R_{l}}\\ x_{L}\\ x_{I}\end{array}\right],$$

$$Q_{\text{link}}=\left[\begin{array}[]{ccc}Q_{\text{ER}}&Q_{\text{EL}}&Q_{\text{EI}}\\ Q_{\text{CR}}&Q_{\text{CL}}&Q_{\text{CI}}\\ Q_{\text{RR}}&Q_{\text{RL}}&Q_{\text{RI}}\\ \end{array}\right].$$

$$v_{R_t}=R_t{j}_{R_t}{v}_{R_l}=R_l{j}_{R_l},$$

$$\frac{d}{dt}(C{v}_C)=j_C\frac{d}{dt}(L{j}_L)=v_L.$$

$$\frac{d}{dt}\left[\begin{array}[]{c}C\>v_{C}\\ L\>j_{L}\end{array}\right]=\left[\begin{array}[]{c}j_{C}\\ v_{L}\end{array}\right].$$

$$-j_C$$

$$v_L$$

$$\begin{split}&\frac{d}{dt}\left[\begin{array}[]{c}C\>v_{C}\\ L\>j_{L}\end{array}\right]=\overbrace{\left[\begin{array}[]{cc}&-Q_{\text{CL}}\\ Q_{\text{CL}}^{T}&\end{array}\right]}^{\mathcal{M}_{c}}\left[\begin{array}[]{c}v_{C}\\ j_{L}\end{array}\right]+\\ &\underbrace{\left[\begin{array}[]{cc}&-Q_{\text{CI}}\\ Q_{\text{EL}}^{T}&\end{array}\right]}_{\mathcal{M}_{s}}\left[\begin{array}[]{c}v_{E}\\ j_{I}\end{array}\right]+\underbrace{\left[\begin{array}[]{cc}-Q_{\text{CR}}&\\ &Q_{\text{RL}}^{T}\end{array}\right]}_{\mathcal{M}_{d}^{T}}\left[\begin{array}[]{c}j_{R_{l}}\\ v_{R_{t}}\end{array}\right],\end{split}$$

$$-R_t^{-1}{v}_{R_t}$$

$$R_l{j}_{R_l}$$

$$\begin{split}\overbrace{\left[\begin{array}[]{cc}R_{l}&-Q_{\text{RR}}^{T}\\ Q_{\text{RR}}&R_{t}^{-1}\end{array}\right]}^{{\alpha^{-1}}}\left[\begin{array}[]{c}j_{R_{l}}\\ v_{R_{t}}\end{array}\right]&=\overbrace{\left[\begin{array}[]{cc}Q_{\text{ER}}^{T}&\\ &-Q_{\text{RI}}\end{array}\right]}^{-\mathcal{M}_{sd}}\left[\begin{array}[]{c}v_{E}\\ j_{I}\end{array}\right]\\ &+\left[\begin{array}[]{cc}Q_{\text{CR}}^{T}&\\ &-Q_{\text{RL}}\end{array}\right]\left[\begin{array}[]{c}v_{C}\\ j_{L}\end{array}\right],\end{split}$$

$$\frac{d}{dt}\left[\begin{array}[]{c}C\>v_{C}\\ L\>j_{L}\end{array}\right]=\mathcal{A}(t)\left[\begin{array}[]{c}v_{C}\\ j_{L}\end{array}\right]+\mathcal{B}(t)\left[\begin{array}[]{c}v_{E}\\ j_{I}\end{array}\right],$$

$$\mathcal{A}(t)={\mathcal{M}}_c-{\mathcal{M}}_d^T\alpha (t){\mathcal{M}}_d,$$

$$\mathcal{B}(t)={\mathcal{M}}_s-{\mathcal{M}}_d^T\alpha (t){\mathcal{M}}_{sd}.$$

$$x=\left[\begin{array}[]{c}q\\ \phi\end{array}\right]=\underbrace{\left[\begin{array}[]{cc}C&\\ &L\end{array}\right]}_{\mathcal{H}^{-1}}\left[\begin{array}[]{c}v_{C}\\ j_{L}\end{array}\right],$$

$$\frac{dx}{dt}=\mathcal{A}(t)\mathcal{H}(t)x+\mathcal{B}(t)s(t),$$

$$E(x,t)=\frac{1}{2}{x}^T\mathcal{H}(t)x.$$

$$\frac{dx}{dt}=\mathcal{A}(t)\nabla E(x,t)+\mathcal{B}(t)s(t).$$

$$\begin{array}{rl}\frac{d}{dt}E(x,t)& =\nabla E(x,t{)}^T\frac{dx}{dt}+\frac{\partial }{\partial t}E(x,t)\\ & =\underset{\text{dissipation}}{\underbrace{\nabla E^T\mathcal{A}(t)\nabla E}}+\underset{\text{forcing}}{\underbrace{\nabla E^T\mathcal{B}(t)s(t)}}+\underset{\text{driving}}{\underbrace{\frac{\partial E}{\partial t}}}.\end{array}$$

$${\dot{E}}_{\text{diss}}=\nabla E^T\mathcal{A}\nabla E=-\nabla E^T{\mathcal{M}}_d^T(\alpha {)}_s{\mathcal{M}}_d\nabla E,$$

$$\dot{Q}=\sum_{r}\dot{Q}_{r}=\left[\begin{array}[]{cc}j_{R_{l}}^{T}&v_{R_{t}}^{T}\end{array}\right]\overbrace{\left[\begin{array}[]{cc}R_{l}&\\ &R_{t}^{-1}\end{array}\right]}^{R}\left[\begin{array}[]{c}j_{R_{l}}\\ v_{R_{t}}\end{array}\right]\geq 0.$$

Full text

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 kish2002 . 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 seifert2012 ; rao2018 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 van2004 ; garnier2005 ; ciliberto2013 ; pekola2015 , 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 graph-theoretic description of the network dynamics in terms of capacitor charges and inductor currents developed in electrical engineering balabanian1969 ; desoer2010 . 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 are 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. The reason is that, depending on their topology, some circuits might display diverging heat currents under the white noise idealization. This fact might be missed by an analysis of the global energy budget alone, and is related to the anomalous thermodynamic behaviour of overdamped models, that are however valid from a purely dynamical point of view celani2012 ; polettini2013 ; bo2014 ; murashita2016 . We clarify this issue and obtain a sufficient and necessary condition on the topology of a circuit for it to be thermodynamically consistent.

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 clerk2010 ; macklin2015 . In fact, the possible use of these circuits as cooling devices has been pointed out before niskanen2007 ; bergeal2010 . 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 schmid1982 . 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 an important 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 freitas2017 .

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 $u_{1},\cdots,u_{n}$ and the edge currents $j_{1},\cdots,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 conditions111 A given circuit can always be made to satisfy conditions (i) and (ii) by adding small stray inductances or capacitances.: (i) The circuit graph has no loops222A loop or cycle is a sequence of edges forming a path so that the first and last node coincide. formed entirely by capacitors and voltage sources, and (ii) The circuit graph has no cut-sets333A cut-set or cocycle is a subset of edges such that their removal splits the graph into at least two disconnected parts. As a loop, it can be oriented, with the orientation indicating a prefered direction from one of the disconnected parts to the other. 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).

The mathematical description of the circuit is easily constructed based on a tree of the circuit graph. A tree is a fully connected subgraph with no loops (see Figure 1-(b)). Under conditions (i) and (ii) it is always possible to find a tree 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. We will call such a tree a normal tree and we will base our dynamical description of the circuit on it. Following the terminology of balabanian1969 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). Thus, all capacitors and voltage sources are twigs, while all inductors and current sources are links. The edges corresponding to resistors can be split in two groups according to whether or not they are part of the normal tree. Every tree has $b_{t}=n-1$ twigs and therefore $b_{l}=b-(n-1)$ links.

II.1 Loops and cut-sets matrices

Given a set of oriented loops we define the loop matrix $B$ as follows:

[TABLE]

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:

[TABLE]

where $\mathds{1}_{k}$ is the $k\times k$ identity matrix. For the normal tree indicated in Figure 1-(b), the matrix $B_{\text{twig}}$ is:

[TABLE]

and specifies which twigs are involved in the loops corresponding to each link.

Given a set of oriented cut-sets (or cocycles) we define the cut-set matrix $Q$ as:

[TABLE]

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:

[TABLE]

With the same ordering as before, the matrix $Q_{\text{link}}$ for the circuit of Figure 1 is:

[TABLE]

and as we see specifies which links belong to the cut-set corresponding to a given twig. More details about the construction of the loop and cutset matrices can be found in Appendix A.

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_{\text{twig}}=-Q_{\text{link}}^{T}$. We also mention that loops and cut-sets can be identified algebraically from the incidence matrix of the graph, and that the loop and cut-set matrices satisfy additional algebraic relations polettini2015 .

II.2 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:

[TABLE]

These 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:

[TABLE]

where $v^{T}=[v_{t}^{T},v_{l}^{T}]$ and $j^{T}=[j_{t}^{T},j_{l}^{T}]$. 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$.

II.3 Block structure of $j$, $v$ and $Q_{t}$

In the following we will assume the previous splitting of the vectors $j$ and $v$:

[TABLE]

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:

[TABLE]

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_{\text{link}}$:

[TABLE]

Thus, for example, each column of the block $Q_{\text{CL}}$ correspond to an inductive edge (which is a link) and its rows specify which capacitive edges (twigs) are involved in the cut-set corresponding to that inductive edge. For the matrix of Eq. (6), we have $Q_{\text{CL}}=\left[\begin{smallmatrix}-1\\ 0\end{smallmatrix}\right]$.

II.4 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,

[TABLE]

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:

[TABLE]

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 and/or describe non-linearities.

III Derivation of dynamical equations

We want to obtain an 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:

[TABLE]

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 parts of Eq. (7):

[TABLE]

and we then obtain:

[TABLE]

where we have defined the matrices $\mathcal{M}_{c}$, $\mathcal{M}_{s}$ and $\mathcal{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_{R_{t}}$, since $v_{E}$ and $j_{I}$ are given. For this we use the constitutive relations of Eq. (12) and two additional Kirchhoff relationships:

[TABLE]

which can be rewritten as:

[TABLE]

where we have defined the additional matrices $\alpha$ and $\mathcal{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:

[TABLE]

where the matrix coefficients are given by:

[TABLE]

and

[TABLE]

Note that ${\alpha}$ might depend on time if resistances do. Some comments about the structure of the matrix $\mathcal{A}$ are in order. The first term in Eq. (20), $\mathcal{M}_{c}$, describes the conservative interchange of energy between capacitors and inductors. In some cases it can be 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 $\mathcal{M}_{d}^{T}{\alpha}(t)\mathcal{M}_{d}$ describes the interchange of energy between capacitors and inductors that is allowed by resistive channels, and the symmetric part describes the loss of energy of those elements (see Eq. 27). Importantly, both the symmetric and antisymmetric part of $\mathcal{M}_{d}^{T}{\alpha}(t)\mathcal{M}_{d}$ respect the block structure of $\mathcal{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 C.

III.1 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 $\phi=L\>j_{L}$, respectively. We will group them in a column vector $x$. Thus, we can write

[TABLE]

where we have defined the matrix $\mathcal{H}$. With these definitions, the dynamical state equation now reads:

[TABLE]

where $s^{T}=[v_{E}^{T},j_{I}^{T}]$ is a vector grouping the voltage and currents of the sources, that can depend on time.

III.2 Linear energy storage elements

If we consider the particular case in which the matrices $C$ and $L$, and thus also $\mathcal{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:

[TABLE]

Then, we can write the dynamical state equation as:

[TABLE]

Note that Eq. (25) still allows for non-linear resistive relations, and in that case the matrices $\mathcal{A}$ and $\mathcal{B}$ have a tacit dependence on $x$ (through ${\alpha}$).

III.3 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:

[TABLE]

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 $\mathcal{A}$ plays a role in the expression $\nabla E^{T}\mathcal{A}\>\nabla E$, and from the definition of $\mathcal{A}$ in Eq. (20), we see that it can only be different from zero if there are resistors in the circuit. Explicitly,

[TABLE]

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 is $\dot{Q}_{r}=j_{r}v_{r}$, and the total rate of heat production reads:

[TABLE]

Using Eq. (18) to eliminate the resistor currents and voltages we find:

[TABLE]

As we will see next, the first term of the previous expression is exactly $-\dot{E}_{\text{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 $-\dot{E}_{\text{diss}}$ and the first term of Eq. (29) can be established from the following property of the matrix ${\alpha}$:

[TABLE]

that can be proven from the definition of ${\alpha}$ via $2\times 2$ block matrix inversion. Using this, we can write the following expression for the energy balance of the circuit:

[TABLE]

where the total work rate $\dot{W}$ is the sum of the rates of work performed by the sources, $\dot{W}_{s}$, and by an external agent that drives the circuit parameters, $\dot{W}_{d}$. They are given by:

[TABLE]

and

[TABLE]

As expected, using the fact that $j_{C}^{T}v_{C}+j_{L}^{T}v_{L}=\frac{dE}{dt}-\frac{\partial E}{\partial t}$ and Tellegen’s theorem, we find that

[TABLE]

III.4 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 $\dot{\Sigma}_{r}=\dot{Q}_{r}/T_{r}$, where $T_{r}$ is the temperature of the considered resistor. Therefore, we can write the total entropy production as

[TABLE]

where $\beta$ is a diagonal matrix with the inverse temperatures of the resistors as non-zero elements. As with $\dot{Q}$, Eq. (18) can be employed to express $\dot{\Sigma}$ 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 ciliberto2013 ; ciliberto2013b , 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 delvenne2014 . 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.

IV Stochastic dynamics

IV.1 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 johnson1928 ; nyquist1928 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 $\Delta v$ of this source has zero mean, $\left\langle\Delta v(t)\right\rangle=0$, and for high temperatures is modeled as delta correlated:

[TABLE]

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(\omega)=Rk_{b}T/\pi$. For low temperatures or high frequencies, one should instead consider the ‘symmetric’ schmid1982 ; clerk2010 quantum noise spectrum:

[TABLE]

where $N(\omega)=(\exp(\hbar\omega/(k_{b}T))-1)^{-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.4. 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.

IV.2 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. (15-b) becomes

[TABLE]

and similar modifications to Eqs. (17) lead to the following generalization of Eq. (18):

[TABLE]

where $\Delta v_{R_{l}}$ and $\Delta v_{R_{t}}$ 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:

[TABLE]

where the vector $\eta(t)$ groups the random voltages in the following way:

[TABLE]

In what follows, for simplicity, we will omit the explicit dependence on time of $\mathcal{A}$, $\mathcal{B}$ and $\alpha$, which are anyway constant if the resistances are also 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 $\Delta v_{R_{l}}$. Also, for high temperatures, taking into account the form of the spectrum for each random voltage we can write:

[TABLE]

where $R$ is the matrix defined in Eq. (28), $\beta$ is the matrix of inverse temperatures defined in Eq. (35), and $\xi(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:

[TABLE]

where $\mathcal{A}$ and $\mathcal{B}$ are the same as in the deterministic case. The last sum is over every resistor in the circuit, and $\mathcal{C}_{r}$ is given by

[TABLE]

where $\Pi_{r}$ is a projector over the one-dimensional subspace corresponding to resistor $r$. Note that the matrices $\mathcal{C}_{r}$ depend on time only if resistances do. Also, they always satisfy $\mathcal{C}_{r}\mathcal{C}_{r^{\prime}}^{T}=0$ for $r\neq r^{\prime}$. However the other possible product $\mathcal{C}_{r}^{T}\mathcal{C}_{r^{\prime}}$ does not vanish in general. Eq. (43) is the stochastic generalization of the usual state equation for electrical circuits.

IV.3 Mean values and covariance matrix dynamics

For linear circuits, the mean values $\left\langle x\right\rangle$ of the charges and fluxes evolve according to the deterministic equation of motion:

[TABLE]

The evolution of the covariance matrix $\sigma(t)=\left\langle(x-\left\langle x\right\rangle)(x-\left\langle x\right\rangle)^{T}\right\rangle$ can be easily derived from the Langevin equation and reads:

[TABLE]

It is important to note that for linear circuits the evolution of the covariance matrix is completely independent from 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 nor their correlations.

V Stochastic thermodynamics

V.1 Fluctuation-dissipation relation

We now consider the particular case in which the circuit is not driven (its parameters are time independent), all the resistors are at the same temperature $T$, and no voltages or currents are applied ($s=0$). Then, according to equilibrium thermodynamics the system must attain a thermal state for long times, with a distribution:

[TABLE]

This is a Gaussian state with covariance matrix

[TABLE]

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 $\sigma_{\text{st}}$ must satisfy

[TABLE]

Also, using the definition of the matrices $\mathcal{C}_{r}$ (Eq. (44)), we see that:

[TABLE]

where in the second equality we used the previously mentioned property ${\alpha}R{\alpha}^{T}=({\alpha})_{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 $\sigma_{\text{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:

[TABLE]

where the mean values $\left\langle x\right\rangle$ are given by the stationary solution of Eq. (45). However, $p_{\text{disp}}$ is not an equilibrium state since, at variance with $p_{\text{th}}$, it has a non-vanishing entropy production given by Eq. (35).

V.2 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 $\left\langle x\right\rangle$ and the covariance matrix $\sigma$: $\left\langle E(x,t)\right\rangle=E(\left\langle x\right\rangle,t)+\operatorname{Tr}[\mathcal{H}\sigma]/2$. Therefore, we have

[TABLE]

where $E(x,t)$ is the energy function in Eq. (24). Thus, the energy balance takes the same form as before

[TABLE]

with the following expressions for the driving work and heat terms:

[TABLE]

and

[TABLE]

where $\dot{Q}(x,t)$ is the heat rate of Eq. (29). Finally, the work performed by the sources is equal to its deterministic value in Eq. (33) (for noiseless sources):

[TABLE]

V.3 Local heat currents

Using Eq. (46) for $d\sigma/dt$ and the FD relation of Eq. (50), we can rewrite Eq. (55) for the total heat dissipation rate as:

[TABLE]

where for clarity we omitted the explicit time dependence of $\mathcal{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 $\langle\dot{Q}_{r}\rangle$ (the rate of energy dissipation in resistor $r$):

[TABLE]

This heuristic definition of the local heat currents is sometimes considered in the literature (see, for example, parrondo1996 ), but it is not always correct. In fact, if we add to the quantities $\langle\dot{Q}_{r}\rangle$ defined in the previous equation a term of the form $\sum_{r^{\prime}}\Delta\dot{Q}_{r,r^{\prime}}$, for any antisymmetric tensor $\Delta\dot{Q}_{r,r^{\prime}}$, we find that the total heat rate $\langle\dot{Q}\rangle=\sum_{r}\langle\dot{Q}_{r}\rangle$ 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

[TABLE]

where $\Delta v_{r}$ is the random voltage associated to each resistor. That quantity, however, is found to be divergent in the general case. 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. There are two different reasons why this might happen. On one hand, the initial description of the circuit might be missing important degrees of freedom that are relevant for the thermodynamics, while still being valid from a purely dynamically point of view. This was recently discussed in celani2012 ; polettini2013 ; bo2014 ; murashita2016 , in connection to the overdamped approximation in stochastic dynamics, even though not in the specific context of electrical circuits. In particular, in murashita2016 it was shown that when a particle is in simultaneous contact with multiple reservoirs at different temperatures there is a transfer of heat associated to its momentum degree of freedom at a rate that is inversely proportional to the relaxation timescale, which is taken to be infinitesimal in the overdamped limit (Figure 2-(b) is actually an example of this situation, as explained in the caption). On the other hand, heat currents might also be ill-defined independently of whether or not some kind of underdamped approximation was made in the description of the circuit. In those cases one should consider a more realistic model for the noise in the resistors, going beyond the white noise idealization. This is the more general way to address this issue, which is discussed in the next section.

From this discussion and the examples of Figure 2, we expect that the possibility to define finite heat currents in a given circuit is related to its topology. Indeed, this intuition is confirmed in the next section, where it is shown that the quantities $\langle\dot{Q}_{r}\rangle$, as defined in Eq. (59), 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_{\text{RR}}=0$ in Eq. (11)). In that case, Eqs. (58) and (59) coincide.

V.4 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 is

[TABLE]

where $v_{r}=R_{r}j_{r}$ is the voltage drop in the ideal resistor and $\Delta v_{r}$ is the random voltage. We want to express the quantity $\dot{Q}_{r}$ in terms of the state $x$ of this circuit. For this, we first write it as

[TABLE]

where $\Pi_{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_{R_{t}}$. 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:

[TABLE]

Notice that the last term in this expression is quadratic in the noise variables. As a consequence, it will give a divergent contribution to $\langle\dot{Q}_{r}\rangle$ in the white noise limit, proportional to $\delta(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 rate $\dot{Q}$ 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}(\omega)$ vanishing for large frequencies, such that its noise spectrum is given by $S_{r}(\omega)=(Rk_{b}T_{r}/\pi)J_{r}(\omega)$. However, this is equivalent to appropriately ‘dressing’ a white-noise resistor with inductors and/or capacitors that can be considered part of the circuit444 This is fully analogous to well known ‘Markovian embedding’ techniques. (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 fact, we see that the quadratic terms in the noise vanish if the matrix $\alpha^{t}\Pi_{r}R\alpha$ is block diagonal. In turn, from the definition of $\alpha$ we can see that this happens if $Q_{\text{RR}}=0$. Thus, $Q_{\text{RR}}=0$ is a sufficient condition for the local heat currents $\dot{Q}_{r}$ to be well behaved. In Appendix B it is shown that this condition is also necessary. We also show in the same Appendix that this condition can be restated in a way that does not make reference to any tree. To finish this discussion, we note that since $\sum_{r}\Pi_{r}=\mathds{1}$ and $\alpha^{t}R\alpha$ is always block diagonal (recall Eq. (30)), we see that the total heat rate $\dot{Q}=\sum_{r}\dot{Q}_{r}$ is always well behaved in the white noise limit, even if $Q_{\text{RR}}\neq 0$.

Thus, assuming $Q_{\text{RR}}=0$, the local heat currents are simplified to:

[TABLE]

where $\eta$ is the vector of random voltages and currents defined in Eq. (42) and we used the fact that for $Q_{\text{RR}}=0$ we have $\alpha=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 gardiner2009 (Chapter 8), in the white noise limit the corresponding stochastic dynamics can be described by the following set of Ito differential equations:

[TABLE]

and

[TABLE]

where $dW$ is a vector of independent Wiener processes differentials. Taking the mean value of Eq. (65) we recover Eq. (58) for $\langle\dot{Q}_{r}\rangle$. The Fokker-Planck equation for the joint probability distribution $P(x,Q_{r})$ corresponding to the previous Ito differential equations reads:

[TABLE]

where $\nabla$ is the nabla operator with respect to the variables $x$. The corresponding equation for the reduced probability distribution $p(x)=\int dQ_{r}\>P(x,Q_{r})$ is just

[TABLE]

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 murashita2016 ; jakvsic2017 ; damak2019 , valid for finite time protocols or asymptotic steady states. However, in this article we will focus only on the behaviour of the mean values.

V.5 Entropy production

We consider the continuous Shannon entropy

[TABLE]

associated to the distribution $p(x)$ and the total entropy production rate

[TABLE]

Using Eq. (67) one can show that $\dot{\Sigma}$ is non-negative:

[TABLE]

where $j_{r}(x,t)=\mathcal{H}(t)x+k_{b}T_{r}\nabla\log(p(x,t))$. Eq. (70) is the Second Law of thermodynamics for the circuit.

This total entropy production can be decomposed as a sum of adiabatic and non-adiabatic contributionsesposito2010 ; van2010 ; esposito2010pre ; spinney2012 , as shown in Appendix F. The adiabatic contribution is positive definite and for time independent circuits it is the only non vanishing contribution for large times. In general the non-adiabatic contribution can have any sign, but for overdamped circuits (for example, circuits with no capacitors or no inductors) it is also positive definite. If the circuit is time independent, then this contribution equals $-k_{b}$ times the time derivative of the relative entropy $H(p|p_{\text{st}})$ between the instantaneous state $p(x,t)$ and the stationary one $p_{\text{st}}(t)$ (which for linear circuits is unique, although it might depend on the initial conditions). Thus, for time independent circuits $H(p|p_{\text{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_{\text{st}})$ due to a conservative flow in phase space, and vanish identically in isothermal conditions (when $p_{\text{st}}$ is an equilibrium state). These findings are analogous to the results of spinney2012 .

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 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 ciliberto2013 ; ciliberto2013b ; karimi2016 .

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},\phi)^{T}$, where $q_{i}$ is the charge in the capacitor $C_{i}$ and $\phi$ 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

[TABLE]

and from this we can construct the matrix $\mathcal{M}_{c}$ and $\mathcal{M}_{d}$:

[TABLE]

Also, in this case $\alpha=R^{-1}=\text{diag}(1/R_{1},1/R_{2})$ and therefore:

[TABLE]

Finally, we have $\mathcal{H}=\text{diag}(C_{1},C_{2},L)^{-1}$, $\mathcal{C}_{1}\mathcal{C}_{1}^{T}=\text{diag}(R_{1}^{-1},0,0)$ and $\mathcal{C}_{2}\mathcal{C}_{2}^{T}=\text{diag}(0,R_{2}^{-1},0)$.

VI.1 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

[TABLE]

where $\bar{T}=(T_{1}+T_{2})/2$ and $\Delta T=T_{1}-T_{2}$. Similar results can be found in ciliberto2013 ; ciliberto2013b . Thus, the first term in the previous expression is just the equilibrium covariance matrix corresponding to the mean temperature $\bar{T}$. 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:

[TABLE]

where we have introduced the two characteristic timescales of the circuit, $\tau_{d}=RC$ and $\tau_{0}=\sqrt{LC}$, respectively associated to the dissipation rate and the free oscillations period. The entropy production can be computed from the above heat rates and reads:

[TABLE]

Thus, in this case, the only dissipation in this circuit corresponds to static heat conduction from the hot to the cold reservoir.

VI.2 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+\Delta C\cos(\omega_{d}t)$ and $C_{2}(t)=C+\Delta C\cos(\omega_{d}t+\theta)$. Thus, both capacitances are driven with the same angular frequency $\omega_{d}$ and the same amplitude $\Delta C$, but with some fixed phase difference $\theta\in(-\pi,\pi)$. The matrices describing the circuit are the same as before except for the energy matrix $\mathcal{H}(t)$ that now depends on time. Since it is periodic, it can be decomposed as a Fourier series:

[TABLE]

To lower order in $\Delta C$ we have only three Fourier components: $\mathcal{H}_{0}=\text{diag}(C,C,L)^{-1}$ and $\mathcal{H}_{\pm 1}=-\Delta C/(2C^{2})\>\text{diag}(1,e^{\pm i\theta},0)$. We focus in regimes when a stable stationary state is reached. We note that this is not always the case due to the phenomenon of parametric resonance poulin2008 . 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

[TABLE]

Inserting the previous two Fourier decompositions into Eq. (46), we obtain an algebraic equation from which it is possible to obtain the coefficients $\sigma_{k,k^{\prime}}^{2}$ in terms of $\mathcal{H}_{0}$ and $\mathcal{H}_{\pm 1}$. This technique is explained in Appendix E. Once the Fourier components $\sigma_{k,k^{\prime}}^{2}$ have been determined, we compute the average heat and work rates per cycle. Explicitly, we consider the quantity

[TABLE]

where $X$ stands for $Q_{1}$, $Q_{2}$, $W$ or $E$. We note that $\langle\dot{E}\rangle_{c}=0$ since the asymptotic state of the system is periodic. Thus, averaging the balance of energy (Eq. (31)) during a cycle we obtain

[TABLE]

where $\langle\dot{W}\rangle_{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 E, to lower order in $\Delta C$ and to second order in $\omega_{d}$, $\langle\dot{W}\rangle_{c}$ is given by

[TABLE]

while the average heat currents are

[TABLE]

We then see that the rate of heat pumping from one resistor to the other (the first term in Eq. 82) is proportional to $\omega_{d}$, while the rate at which work is performed by the driving (or equivalently, the total dissipated heat rate) scales as $\omega_{d}^{2}$. Also, the pumping of heat is maximized and the dissipated work minimized for $\theta=\pm\pi/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

[TABLE]

We note that in this isothermal case the CoP is not bounded by the Second law, and in fact diverges in the quasistatic limit $\omega_{d}\to 0$. There is a maximum driving frequency such that cooling is not possible above it (in the considered regime of low $\omega_{d}$ and $\Delta C$). It corresponds to $\text{CoP}=0$ and for $\theta\in[0,\pi]$ reads

[TABLE]

There is also an optimal frequency $\omega_{d}^{\text{opt}}$, in the sense that the heat extracted from one of the resistors is maximized. We can obtain it by optimizing $\langle\dot{Q}_{1}\rangle_{c}$ in Eq. (82) with respect to $\omega_{d}$, and in this way we find that $\omega_{d}^{\text{opt}}=\omega_{d}^{\text{max}}/2$.

An intuitive understanding of how the cooling effect is achieved can be obtained by analyzing Eq. (58) for the heat currents. We note that the stochastic contribution is proportional to the difference between the energy contained in the circuit elements connected to a given resistor, and the energy they would have if they were in thermal equilibrium at the temperature of that resistor. For example, for the circuit of Figure 3, we have $\langle\dot{Q}_{1}\rangle=[\left\langle q_{1}^{2}\right\rangle\!(t)/(2C_{1})-k_{b}T_{1}/2]/(R_{1}C_{1}/2)$. Thus, we see that $\langle\dot{Q}_{1}\rangle$ will be positive whenever the average energy contained in $C_{1}$, $\left\langle q_{1}^{2}\right\rangle\!(t)/(2C_{1})$, is larger that the corresponding to equilibrium, $k_{b}T_{1}/2$. Therefore, to achieve cooling of $R_{1}$ we require the variance $\left\langle q_{1}^{2}\right\rangle\!(t)$ to be, on average during a cycle, lower than its equilibrium value $k_{b}T_{1}C_{1}$. This reduction in the variance in one degree of freedom with respect to its equilibrium value is the classical analogue to the well known concept of quantum squeezing rugar1991 ; natarajan1995 . In quantum electronic and quantum optical setups, squeezing is a useful resource for metrology and it is usually achieved by means of some form of parametric driving yariv1967 ; anisimov2010 ; pirkkalainen2015 ; wollman2015 . Here we see its additional role as a thermodynamic resource, which supports the already mentioned observation that parametric amplifiers can also be employed as refrigerators bergeal2010 . Thus, an optimal cooling strategy is one in which a highly squeezed state is created and maintained in a dissipative and non-equilibrium environment. Finally, we mention that the methods of Appendix E can be employed to numerically optimize thermal cycles.

VI.3 Non-isothermal case

If the resistor temperatures are different the heat currents are

[TABLE]

In contrast to the isothermal case, the term of second order in $\omega_{d}^{2}$ is too involved to be shown 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, $\Delta 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 $\theta=\pi/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 $\omega_{d}>\omega_{d}^{\text{min}}$. For $\theta=\pi/2$ this minimum cooling frequency $\omega_{d}^{\text{min}}$ reads

[TABLE]

and we can write the heat rate $\langle\dot{Q}_{1}\rangle$ as:

[TABLE]

The previous considerations do not take into account the terms of second order in $\omega_{d}$. From the expression of the heat currents in the isothermal case, Eq. (82), we know that these corrections correspond to heating and establish a maximum driving frequency $\omega_{d}^{\text{max}}$ such that cooling is not possible above it, Eq. (84). Thus, for cooling to be possible at all we need that $\omega_{d}^{\text{min}}<\omega_{d}^{\text{max}}$, 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 $\omega_{d}$ it is given by the following expression:

[TABLE]

Note that that if $\Delta T\neq 0$ then to lower order in $\omega_{d}$ the average work rate can be positive or negative, depending on the value of $\theta$. This two cases correspond to the device working as a refrigerator or a (non-autonomous) heat engine, respectively. From Eqs. (87) and (88) it follows that the cooling efficiency in this case, to lower order in $\omega_{d}$, is:

[TABLE]

which is of course bounded by the Carnot efficiency $\text{CoP}_{\text{Carnot}}=T_{1}/(T_{2}-T_{1})$.

VI.4 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. For the numerical evaluation we consider $\tau_{d}=\tau_{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: $\Delta C/C=1/2$, $\omega_{d}=10^{-2}\>2\pi/\tau_{d}$ and the optimal phase difference of $\theta=\pi/2$. The analytical and numerical results are compared in Figure 4-(b) for a fixed driving amplitude ($\Delta 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 $\omega_{d}^{\text{max}}$ such that both heat currents are positive if $\omega_{d}>\omega_{d}^{\text{max}}$. However, the analytical results overestimate the value of $\omega_{d}^{\text{max}}$. Analogously, we show in Figure 4-(c) the heat currents for fixed driving frequency ($\omega_{d}/(2\pi)=10^{-2}/\tau_{d}$) and increasing driving amplitude. Again, we see that for low driving amplitude the analytical expressions correctly describe the numerical results.

Finally, in Figure 5 we show the Coefficient of Performance as a function of the driving frequency for different values of $\Delta T$. We see that cooling is possible only in a clearly defined range of driving frequencies. The lower limit of this range increases with $|\Delta T|$, in accordance with Eq. (86), while the upper limit displays a weaker dependence on the same parameter. The points in each curve correspond to maximum cooling power, which is attained in all the cases at a frequency close to the optimal driving frequency of the isothermal case, $\omega_{d}^{\text{opt}}$. In all cases the maximum CoP is only about $1\%$ of the ideal value, $\text{CoP}_{\text{Carnot}}$. Thus, we see that this simple cooling scheme can only withstand small thermal gradients and operate at low efficiencies. More complex schemes are expected to improve these figures of merit, although they will probably share some of the properties of this elementary example, like the existence of minimum and maximum driving frequencies for cooling.

VII Quantum Johnson-Nyquist noise

We now turn to the low temperature regime of our theory. 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 $\hbar\omega$ at the relevant frequencies $\omega$, 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 caldeira1983 . 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 devoret1995 ; solgun2015 ; clerk2010 ; parra2019 ; burkard2004 ; girvin2011 ). 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 devoret1995 . As an example we can take the overdamped circuit of Fig. 2-(c). This circuit has well behaved classical 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 parra2019 ). To do that, one must give a more detailed description specifying stray inductances. Of course, any actual component in a real circuit is characterized by an impedance that is never purely resistive, capacitive or inductive. At a fundamental level we could consider elementary components which are always a combination of a single inductor and a single capacitor (resistors are then modeled as infinite arrays of them). Any circuit constructed in this way can be directly quantized, and based on this quantization the low temperature behaviour can be studied. Circuits like the one in Fig. 2-(c) result from having disregarded the inductive component in the impedance of the actual capacitors, which anyway might be an excellent approximation for practical purposes. However, this procedure of ignoring some degrees of freedom in the circuit is problematic, as we can see by considering the simple example of a series RLC circuit. The relevant frequency scales in this case are given by the dissipation rate $\gamma=R/L$, the oscillation frequency $\omega_{0}=1/\sqrt{LC}$, and the thermal frequency $\omega_{\text{th}}=k_{b}T/\hbar$. Under the high temperature condition $\omega_{\text{th}}\gg\omega_{0},\gamma$, the equilibrium charge variance is $\left\langle q^{2}\right\rangle\simeq k_{b}TC$ and does not depend on $L$. This justifies the use of an overdamped RC model, where the inductance $L$ is neglected compared to $R^{2}C$, in the high temperature regime. In contrast, if we take the overdamped limit $L/(R^{2}C)\to 0$ for low temperatures we obtain that the equilibrium charge variance diverges as $\left\langle q^{2}\right\rangle\propto L^{-1}$ while the flux variance vanishes as $\left\langle\phi^{2}\right\rangle\propto L$ (these scalings are obtained by keeping $\omega_{0}$ constant, otherwise we have $\left\langle q^{2}\right\rangle\propto L^{-1/2}$ and $\left\langle\phi^{2}\right\rangle\propto L^{1/2}$). Thus, in general the covariance matrix of the circuit state is not well defined in this limit, and to compute it one needs to have precise information about the value of $L$. In spite of this, we will see in the following that the heat currents are actually well defined in this kind of overdamped limits, and can be directly evaluated from the overdamped description of circuits even at low temperatures.

To show 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 $\xi(t)$ are not white anymore and display a quantum spectrum. Then, if $\xi_{r}$ is the adimensional noise variable associated with the $r$-th resistor, we must consider the following power spectrum:

[TABLE]

where $N_{r}(\omega)=(e^{\hbar\omega/(k_{b}T_{r})}-1)^{-1}$ is the Planck’s distribution at temperature $T_{r}$. Inverting the previous equation, we can compute the correlation functions:

[TABLE]

where $\Lambda$ 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). However, it is possible to show that for linear circuits that can be directly quantized, the results obtained in this way are fully equivalent to those obtained by a full quantization (under the Markovian approximation) schmid1982 ; freitas2017 .

VII.1 Covariance matrix and heat currents for quantum noise

Since $\left\langle\xi\right\rangle(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 $\xi(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 gardiner2009 ). 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-\left\langle x\right\rangle$:

[TABLE]

Given the initial value $y(0)$, the solution to this equation can be written as

[TABLE]

where the retarded Green’s function $G(t,t^{\prime})$ is defined as the solution of

[TABLE]

with $G(t,t^{\prime})=0$ for $t<t^{\prime}$ (from this it follows that $G(t^{\prime},t^{\prime})=\mathds{1}$). Then, the covariance matrix can be expressed as

[TABLE]

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, $\left\langle y(0)\xi(\tau)^{T}\right\rangle=0$, and therefore the second term in the previous equation vanish. We will also assume, for simplicity, that the resistances, and thus the matrices $\mathcal{A}$ and $\mathcal{C}_{r}$, are constant. Then, taking the time derivative of Eq. (95) and using Eq. (94), we obtain the following differential equation:

[TABLE]

where $\mathcal{I}_{r}(t)$ is the convolution between the Green’s function $G(t,t^{\prime})$ and the correlation function of resistor $r$:

[TABLE]

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 $\langle\xi_{r}(0)\xi_{r}(\tau)\rangle\to\delta(\tau)$, and therefore $\mathcal{I}_{r}(t)\to G(t,t)/2=\mathds{1}/2$.

Based on these results, we can now derive an expression for the local heat currents that, in contrast to Eq. (58), is exact and valid for arbitrary temperatures. In the quantum case the total heat rate is also given by Eq. (55): $\dot{\left\langle Q\right\rangle}=\dot{Q}(\left\langle x\right\rangle,t)-\frac{1}{2}\operatorname{Tr}\left[\mathcal{H}\frac{d\sigma}{dt}\right]$. However, this time we should replace $d\sigma/dt$ by Eq. (96). Using this and the FD relation, we can write:

[TABLE]

In analogy with the classical case, under the condition $Q_{\text{RR}}=0$, we can identify the local heat currents as:

[TABLE]

This expression can be evaluated using the above equations for $\sigma$ and $\mathcal{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.

VII.2 Asymptotic covariance matrix and heat currents

If the system is asymptotically stable, i.e, if $G(t,t^{\prime})\to 0$ for $|t-t^{\prime}|\to\infty$, then the first two terms in Eq. (95) can be neglected for sufficiently long times. By expressing the correlations $\left\langle\xi(\tau)\xi(\tau^{\prime})^{T}\right\rangle$ in terms of the power spectrum $\mathcal{S}_{r,r^{\prime}}(\omega)$ via inversion of Eq. (90), we can rewrite the last term in Eq. (95) as:

[TABLE]

where we have defined the following partial transform of the Green’s function:

[TABLE]

For circuits that can be directly quantized, Eq. (100) is equivalent to what one obtains from a fully quantum model of the network and its environment under the Markovian approximation freitas2017 .

If the circuit parameters are periodically driven, $\hat{G}(t,\omega)$ has the useful property of being asymptotically periodic in time with the same period as the driving, as shown in Appendix D. It also trivially satisfies $\hat{G}(t,\omega)^{*}=\hat{G}(t,-\omega)$, a property that is sometimes used implicitly in the derivations below. The convolution integral $\mathcal{I}_{r}(t)$ can also be expressed in terms of $\hat{G}(t,\omega)$:

[TABLE]

Finally, we note that $\hat{G}(t,\omega)$ can be directly obtained by solving its own evolution equation, that can be derived from Eq. (94) and reads

[TABLE]

with the initial condition $\hat{G}(t=0,\omega)=0$.

Introducing Eqs. (100) and (102) for $\sigma$ and $\mathcal{I}_{r}$ into Eq. (99), we can write the local heat currents as:

[TABLE]

where we introduced the shorthand definition $\mathcal{D}_{r}=\mathcal{C}_{r}\mathcal{C}_{r}^{T}$, that we will employ in the following to simplify the notation. The first term inside the trace, that is quadratic in $\hat{G}(t,\omega)$, is actually closely related to the second one, which is linear in $\hat{G}(t,\omega)$. We can see this by employing Eq. (103) to compute the derivative of $\hat{G}^{\dagger}\mathcal{H}\hat{G}$:

[TABLE]

Using this relationship and Eq. (50), it is possible to rewrite Eq. (104) in the following compact way:

[TABLE]

where $f_{r,r^{\prime}}(t,\omega)$ is a transfer function, specifying how the temperature of resistor $r^{\prime}$ affects the heat current of resistor $r$. For $r\neq r^{\prime}$ it is always positive and is given by:

[TABLE]

while the diagonal elements $f_{r,r}(t,\omega)$ are determined by the following expression for the sum over the first index:

[TABLE]

Equation (106) is the central result of this article. It is a fully general expression for the local heat currents valid for arbitrary temperatures and driving protocols. Although it has been derived based on Eq. (78), which is in principle only valid for circuit descriptions that can be quantized, nothing prevents the evaluation of Eqs. (106), (107) and (108) for general, overdamped circuits. In fact, as we show analytically in Appendix G, the overdamped limit of the transfer function $f_{r,r^{\prime}}(\omega)$ for a underdamped circuit correctly matches the transfer function directly obtained from the corresponding overdamped circuit. Later this is also verified numerically for the cooling scheme of Section VI.

To clarify the physical interpretation of the previous expressions we analyze first the particular case of time independent circuits.

VII.3 Undriven circuits

If the matrix $\mathcal{H}$ is time independent, then for sufficiently long times we have, from Eq. (103), $\hat{G}(t,\omega)=\hat{G}_{0}(\omega)=(i\omega-\mathcal{A}\mathcal{H})^{-1}$ (this is just the Laplace’s transform of the Green’s function evaluated at $i\omega$). Then, asymptotically, the transfer functions $f_{r,r^{\prime}}(\omega)$ are time independent and $\bar{f}_{r^{\prime}}(\omega)=0$, so that $f_{r,r}(\omega)=-\sum_{r^{\prime}\neq r}f_{r^{\prime},r}(\omega)$. Also, in this case $f_{r,r^{\prime}}(\omega)$ is symmetric under interchange of the indexes $r$ and $r^{\prime}$. To see this it is necessary to consider the block structure of the matrix $\mathcal{A}$, that is inherited by the matrix $\hat{G}_{0}(\omega)$, and of the matrices $\mathcal{D}_{r}$. This is discussed in detail in Appendix C, in connection with the invariant nature of the dissipation upon time inversion. Thus, using these properties we recover the usual Landauer-Büttiker expression for the heat currents rego1998 ; dhar2006 ; yamamoto2006 :

[TABLE]

From this equation, the quantity $f_{r,r^{\prime}}(\omega)d\omega$ can be naturally interpreted as the rate at which an excitation with frequency between $\omega$ and $\omega+d\omega$ is transported from resistor $r^{\prime}$ to $r$. We note that the symmetry of the transfer function in the undriven case causes the heat currents to only depend on the differences $(N_{r^{\prime}}(\omega)+1/2)-(N_{r}(\omega)+1/2)$. 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^{\prime}}(\omega)$ is not symmetric in general. In that case, the ground state fluctuations represented by the $1/2$ term 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 as 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 where the external driving is periodic. They are useful for the numerical study of thermal cycles.

VII.4 Periodically driven circuits

We now consider the situation where the matrix $\mathcal{H}(t)$ is a periodic function of time, and therefore can be decomposed as a Fourier series:

[TABLE]

where $\omega_{d}$ is the angular frequency of the driving. In this case, assuming stable dynamics and long times, the function $\hat{G}(t,\omega)$ is also periodic with the same period of the driving, as shown in Appendix D. Thus, the following decomposition holds asymptotically

[TABLE]

Then, the asymptotic covariance matrix of the system and the heat currents are also periodic, with period $\tau=2\pi/\omega_{d}$. We thus consider the average values of the heat currents during a driving period, that we denote $\langle\dot{Q}_{r}\rangle_{c}$, as we did in the example of Section VI. They are given by

[TABLE]

where $F_{r,r^{\prime}}(\omega)$ is the asymptotic average of $f_{r,r^{\prime}}(t,\omega)$ during a driving period, that can be expressed in terms of the Fourier components $\mathcal{H}_{k}$ and $\hat{G}_{j}(\omega)$:

[TABLE]

for $r^{\prime}\neq r$, and

[TABLE]

Finally we note that, given the Fourier components $\mathcal{H}_{k}$ of the external driving, the Fourier components $\hat{G}_{j}(\omega)$ of the Green’s function can be found by solving the following infinite set of algebraic equations

[TABLE]

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 D.

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 freitas2017 . 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.

VII.5 Quantum limits for cooling

In this section we explore the low temperature behaviour of the heat currents for the cooling scheme of Section VI. We show how the quantum corrections to the heat currents impose a minimum temperature below which it is not possible to extract heat. Also, this will serve us to show that the overdamped description of circuits can be directly employed to compute the heat currents in the low temperature regime.

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 $\Delta C$, frequency $\omega_{d}$, and phase difference $\theta$ such that for high temperatures heat is extracted from resistor $R_{1}$ (thus, $\langle\dot{Q}_{1}\rangle_{c}<0$), while it is dumped in resistor $R_{2}$ ($\langle\dot{Q}_{2}\rangle_{c}>0$). The heat currents are computed by evaluating Eqs. (112), (113) and (114), based on the overdamped circuit of Figure 3. In Figure 6-(a) 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 $\Delta C/C=10^{-1}$, $\omega_{d}=10^{-2}(2\pi/\tau_{0})$ and $\tau_{d}=2\tau_{0}$ (recall that $\tau_{0}=\sqrt{LC}$ and $\tau_{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. 6-(a)), $\langle\dot{Q}_{1}\rangle_{c}$ becomes positive and cooling stops. This is shown in more detail in Figure 6-(b), where it is also clear that the value of $T^{*}$ decreases with increasing $\tau_{d}/\tau_{0}$, or equivalently with decreasing dissipation rate (In this case $\omega_{d}$ is selected as the optimal driving frequency for each value of $\tau_{d}/\tau_{0}$). This breakdown of the cooling effect is a strong coupling result that cannot be captured with usual approaches based on master equations (see karimi2016 , for example), as discussed in detail in freitas2017 . All these results are independent of the cutoff frequency $\Lambda$ for large $\Lambda$.

In figure 6-(b) we also show the results obtained based on the circuit of Figure 7, for $L^{\prime}/L=10^{-3}$ (indicated by crosses). This circuit can be considered an extension of the one in Figure 3 in which the additional stray inductance $L^{\prime}$ was specified. As shown in more detail in figure 6-(c), the heat currents obtained from the two descriptions match as $L^{\prime}/L\to 0$. Thus, while overdamped descriptions of circuits are not enough to construct a quantum model and compute the quantum state, they are sufficient to compute the quantum corrections to the heat currents.

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 gnesotto2019 ; mura2018 . 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 using RLC circuits.

IX Acknowledgments

We acknowledge funding from the European Research Council project NanoThermo (ERC-2015-CoG Agreement No. 681456).

Appendix A Construction of the loop and cutset matrices

In this section we give a detailed explanation of how the loop and cutset matrices are constructed for the example of Figure 1-(a). The fist step is to identify a normal tree of the circuit graph (one possible normal tree is indicated in blue in Figure 1-(b)). By definition, all the capacitors and voltage sources should be part of the normal tree, while all inductors and current sources should be out of it. Also, as any tree, it cannot contain loops, and all the nodes should be connected. Thus, the edges $V$, $C_{2}$ and $C_{1}$ should be part of the normal tree, while the edges $L$ and $I$ should be out of it. To identify the edges $R_{1}$, $R_{2}$, $R_{3}$, and $R_{4}$, we see that $R_{1}$ should be part of the tree, since otherwise the node at which $L$ arrives would be disconnected. $R_{2}$ should also be part of the tree, since otherwise the node from which $I$ departs would remain disconnected. Thus, $R_{1}$ and $R_{2}$ are twigs. Finally, $R_{3}$ and $R_{4}$ cannot be both part of the tree, since in that case they would form a loop with $V$ and $C_{2}$. They also cannot be both outside the tree, since in that case the central node from which $L$ departs would be disconnected. Thus, one of them should be a twig and the other a link. We are free to choose between the two possible options, and in the example of Figure 1-(b) we have choosen $R_{4}$ as a twig.

Having identified the normal tree, we now proceed to construct the matrices of fundamental loops and cutsets associated to it. We begin with the loops. We can associate a loop to each link. For example, if we add link $L$ to the normal tree, then a loop if formed by the edges $L$, $R_{1}$, $C_{1}$ and $R_{4}$ (See figure 8).

The orientation of the loop is defined to be that of the added link, in this case $L$, which matches the orientation of $R_{1}$ and $C_{1}$, but is contrary to that of $R_{4}$. The row corresponding to this loop in the full loop matrix is then:

[TABLE]

The rows corresponding to the links $R_{3}$ and $I$ are constructed analogously. By adding link $I$ to the normal tree a loop is formed involving the twigs $R_{2}$, $V$, $C_{2}$, $C_{1}$ and $R_{1}$. Adding link $R_{3}$ we form a loop with twigs $V$, $C_{2}$ and $R_{4}$. All the loops are shown in Figure 8. Taking into account their orientations, the final full loop matrix is:

[TABLE]

We see that the first $3\times 3$ block is just the identity, and the rest of the matrix is what is denoted by $B_{\text{twig}}$ in the main text.

The construction of the cutset matrix is done as follows. A cutset can be associated to each twig. For example, if we remove twig $C_{1}$, then the start and end nodes of $R_{1}$ are disconnected from the others (in the tree, not in the full graph). So we consider the splitting of the set of nodes in two subsets: the two nodes connected to $R_{1}$, and all the rest. The edges going from one subset to the other are $I$, $L$, and of course $C_{1}$ (see Figure 9-(a)).

We give to this cutset the same orientation as $C_{1}$, which happens to be contrary to the orientation of $I$ and $L$. Then, the row corresponding to this particular cutset in the full cutset matrix is:

[TABLE]

The other rows corresponding to the twigs $V$, $C_{2}$, $R_{1}$, $R_{2}$ and $R_{4}$ are constructed in the same way. As an additional example, we note that removing twig $R_{4}$ from the normal tree we isolate the central node, and the corresponding cutset is formed by $R_{4}$, $R_{3}$ and $L$, which share the same orientation, as shown in Figure 9-(b). Adding the corresponding row to the previous one we have:

[TABLE]

Filling in the remaining 4 rows to obtain the full cutset matrix, we find that the first $6\times 6$ block is the identity, and the remaining block is the one denoted by $Q_{\text{link}}$ in the main text.

Appendix B Necessary and sufficient condition to have well defined local heat currents for white noise

resistors

The matrix $\alpha^{T}\Pi_{r}R\alpha$ appearing in the last line of Eq. (62) has the following block structure:

[TABLE]

The explicit form of each block can be derived from the definition of $\alpha$ and using $2\times 2$ block matrix inversion. Using this, we see that the mean value of the last term in Eq. (62) can be written as:

[TABLE]

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_{\text{RR}}$ and $Q_{\text{RR}}^{T}B_{r}$ should have null diagonals. The explicit form of $B_{r}$ is $B_{r}=(R_{l}+Q_{\text{RR}}^{T}R_{t}Q_{\text{RR}})^{-1}Q_{\text{RR}}^{T}\pi_{r}(R_{t}^{-1}+Q_{\text{RR}}R_{l}^{-1}Q_{\text{RR}}^{T})^{-1}$ if the index $r$ correspond to a twig resistor, or $B_{r}=-(R_{l}+Q_{\text{RR}}^{T}R_{t}Q_{\text{RR}})^{-1}\pi_{r}Q_{\text{RR}}^{T}(R_{t}^{-1}+Q_{\text{RR}}R_{l}^{-1}Q_{\text{RR}}^{T})^{-1}$ if $r$ correspond to a link resistor, where $\pi_{r}$ is the reduction of the projector $\Pi_{r}$ to the appropriate twigs or links subspace. In any case, it is easy to see that for the matrix $B_{r}Q_{\text{RR}}$ to have null diagonal for any $r$ and for arbitrary values of the resistances, all the components of the matrix $Q_{\text{RR}}$ should vanish.

The condition $Q_{\text{RR}}=0$ means that there are no fundamental cutsets associated to the normal tree containing simultaneously a link resistor and a twig resistor. However, since it is ultimately related to the possibility of defining finite heat currents, its validity should be independent of the choice of normal tree. In fact, it is possible to obtain the following equivalent condition that does not make any reference to a particular tree of the graph:

• •

Given the full graph of the circuit, remove the I and L edges.

• •

Merge the connected components formed by E and C edges into a single node (this will remove internal edges, that might be resistors).

• •

The heat currents are finite if and only if the resulting network of only R edges is acyclic.

The equivalence with the condition $Q_{\text{RR}}=0$ is given by the following two observations: i) if there is a cycle in the final R network, then any normal tree will have a fundamental cutset involving a resistor in the tree and a resistor outside it, ii) if there is no cycle in the final R network, then all the edges in it will be part of any normal tree, and the fundamental cutsets associated to them will not involve any of the previously eliminated resistors (if any was eliminated).

Appendix C Block structure of the matrix $\mathcal{A}$ and the invariance of

dissipation upon time inversion

In this section, we show explicitly how the block structure of the matrix $\mathcal{A}$ is related to the invariance of the dissipation upon time reversal and to the symmetry property of the transfer function $f_{r,r^{\prime}}(\omega)$ involved in the expression for the heat currents for undriven circuits. We first recall the necessary definitions:

[TABLE]

[TABLE]

and

[TABLE]

Using $2\times 2$ block matrix inversion, we see that $\alpha$ has also a block structure like its inverse. Then it is straightforward to see that the matrix $\mathcal{A}$ can be written as:

[TABLE]

where $\bm{s_{1}}$ and $\bm{s_{2}}$ are symmetric matrices with dimensions $N_{C}\times N_{C}$ and $N_{L}\times N_{L}$, respectively ($N_{C}$ is the number of capacitors and $N_{L}$ the number of inductors). Therefore, the rate of energy dissipation is

[TABLE]

where in the second line we introduced the block structure of the matrix $\mathcal{H}=\text{diag}(C^{-1},L^{-1})$ and the state vector $s=(q,\phi)^{T}$. Thus, we see that due to the block structure of the symmetric part of $\mathcal{A}$, the rate of energy dissipation has no cross terms coupling charges $q$ and fluxes $\phi$. Therefore it is an even quantity under time reversal.

We now turn to discuss the symmetry of the transfer function $f_{r,r^{\prime}}(\omega)$ for undriven circuits upon the interchange of $r$ and $r^{\prime}$. This function is given by (Eq. (107) in the main text):

[TABLE]

with $\hat{G}(\omega)=(i\omega\mathds{1}-\mathcal{A}\mathcal{H})^{-1}$. We first note that since $\mathcal{H}$ is a positive definite matrix we can write the product $\mathcal{H}\hat{G}(\omega)$ in the more symmetric form $K(\omega)=\mathcal{H}\hat{G}(\omega)=\sqrt{\mathcal{H}}\left(i\omega\mathds{1}-\sqrt{\mathcal{H}}\mathcal{A}\sqrt{\mathcal{H}}\right)^{-1}\sqrt{\mathcal{H}}$, and thus we see that if $\mathcal{A}$ were symmetric or antisymmetric, then the matrix $K(\omega)$ would inherit that property. In any of those cases the function

[TABLE]

would be trivially symmetric upon $r\leftrightarrow r^{\prime}$. However, $\mathcal{A}$ has no definite symmetry. It has, nevertheless, a block structure that is also inherited by $K(\omega)$:

[TABLE]

where the dimensions of $\bm{s^{\prime}_{1}}(\omega)$, $\bm{s^{\prime}_{2}}(\omega)$ and $\bm{a^{\prime}}(\omega)$ are the same as $\bm{s_{1}}$, $\bm{s_{2}}$ and $\bm{a}$ in Eq. (125). If $K_{s}=\text{diag}(\bm{s^{\prime}_{1}},\bm{s^{\prime}_{2}})$ and $K_{a}=K-K_{s}$ are the symmetric and antisymmetric part of $K$, we have

[TABLE]

were we omitted the dependences in $\omega$. Finally from the definitions of the matrices $\mathcal{C}_{r}$ it can be seen that their products $\mathcal{D}_{r}=\mathcal{C}_{r}\mathcal{C}_{r}^{T}$ are block diagonal whenever $Q_{RR}=0$, which is the condition for the heat currents to be properly defined in the white noise or high temperatures limit. From this it follows that the two last traces in the previous equation are zero, since their arguments have null diagonals. The remaining terms are easily shown to be invariant upon $r\leftrightarrow r^{\prime}$.

Appendix D Periodic driving

In this section we show some useful properties of the Green’s function of periodically driven and stable circuits. We begin with the expression for the transform $\hat{G}(t,\omega)$ given in Eq. (101):

[TABLE]

Here, $G(t,t^{\prime})$ is the solution to

[TABLE]

with $G(t,t^{\prime})=0$ for $t<t^{\prime}$. If $\mathcal{H}(t)$ is a periodic function with period $\tau_{d}=2\pi/\omega_{d}$, and the solution to the previous differential equation is unique, then we have that $G(t,t^{\prime})=G(t+\tau_{d},t^{\prime}+\tau_{d})$. Then:

[TABLE]

where in the first step we just employed a change of variables ($\tau\to\tau-\tau_{d}$), and in the last step we assumed that the system is stable, in the sense that $G(t,t^{\prime})\to 0$ for $|t-t^{\prime}|\to\infty$. 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,\omega)$ is asymptotically periodic, with period $\tau_{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 $\hat{G}(t,\omega)$:

[TABLE]

Then, inverting Eq. (131) and using Eq. (132) (or equivalently, transforming Eq. (132) to obtain Eq. (103) in the main text), we can find the following set of algebraic equations for the coefficients $\hat{G}_{j}(\omega)$:

[TABLE]

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 $|\mathcal{H}_{k}|\ll|\mathcal{H}_{0}|$ for all $k\neq 0$. Then, to first order in $\mathcal{H}_{k\neq 0}$ we have $G_{0}(\omega)\simeq(i\omega\mathds{1}-\mathcal{A}\mathcal{H}_{0})^{-1}$, that is just the transform of the Green’s function on the undriven circuit, and:

[TABLE]

Higher orders in $\mathcal{H}_{k\neq 0}$ can be easily computed. From this solution, we see that for weak driving the range of relevant Fourier components in $\hat{G}(t,\omega)$ is restricted by that of $\mathcal{H}(t)$. Thus, another non-perturbative method to solve Eq. (135) is just to truncate the Fourier space to some maximum number of components given by $|k|\leq k_{\text{max}}$ and $|j|\leq j_{\text{max}}$ and then numerically solve the resulting finite system of linear equations for each value of $\omega$.

Appendix E 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 $\sigma$ for large times can be obtained as the solution of the following Lyapunov equation:

[TABLE]

For driven circuits there is no time-independent asymptotic state and one must solve the dynamical equation:

[TABLE]

However, if the function $\mathcal{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 $\mathcal{H}$). Then we consider $\mathcal{H}(t)=\sum_{k=-\infty}^{k=+\infty}\mathcal{H}_{k}e^{ik\omega_{d}t}$ and introduce the following decomposition for $\sigma(t)$:

[TABLE]

Thus, the problem is now to find the coefficients $\sigma_{k,k^{\prime}}$ in terms of $\mathcal{H}_{k}$. Note that, at variance with a regular Fourier decomposition, the previous expression involves a double summation and as a consequence the coefficients $\sigma_{k,k^{\prime}}$ 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:

[TABLE]

[TABLE]

and

[TABLE]

Then, it can be seen that the coefficients $\sigma_{k,k^{\prime}}$ that give a solution to Eq. (138) can be obtained by solving the following generalized Lyapunov equation:

[TABLE]

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 $\mathcal{H}$ has only $3$ Fourier components, and therefore to lower order in $\mathcal{H}_{\pm 1}$ we can truncate the matrix $S$ to $3$ blocks in each direction. After doing, this we solved symbolically the Lyapunov equation of Eq. 143 using Mathematica.

Appendix F Adiabatic and non-adiabatic decomposition of the entropy production

The Fokker-Planck equation for the circuit state can be cast as

[TABLE]

where the total probability current $J(x,t)$ was split into conservative and dissipative parts, that are respectively given by

[TABLE]

and

[TABLE]

where $\mathcal{A}_{s}$ and $\mathcal{A}_{a}$ are the symmetric and antisymmetric parts of $\mathcal{A}$. The total entropy production rate is:

[TABLE]

with $j_{r}(x,t)=\mathcal{H}(t)x+k_{b}T_{r}\nabla\log(p(x,t))$. Using the FD relation the dissipative probability current $J_{d}$ can be expressed in terms of the functions $j_{r}$ as follows:

[TABLE]

For a given instantaneous value of $\mathcal{H}(t)$, we define the corresponding steady state distribution $p_{\text{st}}(x,t)$ as the one for which $\nabla^{T}J_{\text{st}}(x,t)=0$. It is the probability distribution to which the circuit would eventually relax if the parameters are frozen at the values given by $\mathcal{H}(t)$. The functions $j_{r}(x,t)$ corresponding to $p_{\text{st}}(x,t)$ are denoted as $j_{r}^{\text{st}}(x,t)$. Then, by replacing $j_{r}\to j_{r}-j_{r}^{\text{st}}+j_{r}^{\text{st}}$ in Eq. (147), we obtain the following decomposition of the total entropy production:

[TABLE]

where

[TABLE]

[TABLE]

with $\Delta j_{r}=j_{r}-j_{r}^{\text{st}}$, and finally

[TABLE]

The terms $\dot{\Sigma}_{\text{nad}}$ and $\dot{\Sigma}^{{}^{\prime}}_{\text{nad}}$ vanish in the adiabatic limit of infinitely slow driving, since the state is always the stationary one and $\Delta j_{r}\to 0$. For the same reason they vanish for long times if there is no driving. In contrast, $\dot{\Sigma}_{\text{ad}}$ converges to the entropy production in the stationary state, which is $\sum_{r}\langle\dot{Q}_{r}\rangle/T_{r}$. The non-adiabatic (nad) terms can be interpreted in terms of the relative entropy of the instantaneous state with respect to the stationary state. Indeed, $\dot{\Sigma}_{\text{nad}}$ can be expressed as:

[TABLE]

and therefore in the absence of external driving ($\mathcal{H}(t)$ and therefore $p_{\text{st}}(x,t)$ are constants) it equals $-k_{b}$ times the derivative of the relative entropy $H(p|p_{\text{st}})$. The term $\dot{\Sigma}^{{}^{\prime}}_{\text{nad}}$ also accepts the following expressions:

[TABLE]

and can therefore be interpreted as a the change in the relative entropy $H(p|p_{\text{st}})$ due to the conservative flow $J_{c}$ in phase space. At variance with $\dot{\Sigma}_{\text{nad}}$, $\dot{\Sigma}^{{}^{\prime}}_{\text{nad}}$ is not always positive. It vanishes identically in the following cases: i) for circuits with no inductors or no capacitors (in that case we can consider $\mathcal{A}_{a}=0$ and the dynamics is always overdamped), and ii) in isothermal conditions, since the steady state satisfies $\nabla p_{\text{st}}/p_{\text{st}}\propto Hx$ and therefore $\dot{\Sigma}^{{}^{\prime}}_{\text{nad}}\propto\operatorname{Tr}\left[\mathcal{A}_{a}\mathcal{H}\sigma\mathcal{H}\right]=0$. Finally, adding the last two equations and using Eq. (144), we can see that the sum of the non-adiabatic terms only depends on the dissipative current:

[TABLE]

We note that this quantity is not necesarily positive definite. It is positive for overdamped circuits, since we have seen that in that case $\dot{\Sigma}_{\text{nad}}=0$ and $\Sigma_{\text{nad}}\geq 0$.

Appendix G Underdamped limit for static heat conduction

In this section we illustrate with an example how the heat transfer function of an overdamped circuit is recovered by taking the appropiate limit in the heat transfer function of a underdamped circuit. We consider the circuit in Figure 10-(a), that in the limit $L/R^{2}C\to 0$ reduces to the circuit of Figure 2-(c) in the main text.

We consider the heat transfer function given by Eq. (107) in the main text (only the non-diagonal element is necessary in this case of stationary heat conduction):

[TABLE]

where

[TABLE]

The matrices $\mathcal{A}$, $\mathcal{H}$ and $\mathcal{D}_{1/2}$ for each circuit are to be constructed according to the procedure of Section II in the main text. In this way, we obtain the following transfer function for the circuit of Figure 10-(a):

[TABLE]

and for the circuit of Figure 10-(b):

[TABLE]

Thus, a simple calculation shows that in fact

[TABLE]