Relation between topology and heat currents in multilevel absorption machines

The steady state heat currents of continuous absorption machines can be decomposed into thermodynamically consistent contributions, each of them associated with a circuit in the graph representing the master equation of the thermal device. We employ this tool to study the functioning of absorption refrigerators and heat transformers with an increasing number of active levels. Interestingly, such an analysis is independent of the particular physical implementation (classical or quantum) of the device. We provide new insights into the understanding of scaling up thermal devices concerning both the performance and the magnitude of the heat currents. Indeed, it is shown that the performance of a multilevel machine is smaller or equal than the corresponding to the largest circuit contribution. Besides, the magnitude of the heat currents is well-described by a purely topological parameter which in general increases with the connectivity of the graph. Therefore, we conclude that for a fixed number of levels, the best of all different constructions of absorption machines is the one whose associated graph is as connected as possible, with the condition that the performance of all the contributing circuits is equal.


Introduction
Continuous quantum absorption machines [1] are multilevel systems connected to several thermal baths at different temperatures. Their autonomous functioning can be rigorously described by using the theory of open quantum systems [2]. Some basic models such as the three-level [3,4], the two-qubit [5] and the three-qubit [4,6] absorption refrigerators have been widely employed in establishing fundamental relations in quantum thermodynamics [7]. Besides, several experimental proposals have been put forward, for example those based on nano-mechanical oscillators or atoms interacting with optical resonators [8,9], atoms interacting with nonequilibrium electromagnetic fields [10], superconducting quantum interference devices [11,12], and quantum dots [13]. Further, an experimental realization of a quantum absorption refrigerator has been recently reported [14].
The dynamics of quantum machines is described by a master equation when the coupling with the baths is weak enough. Along this paper we consider in addition systems for which two states with the same energy cannot be connected to a third one through the same bath. This assumption greatly simplifies the quantum master equation as the population and coherence dynamics are decoupled in the system energy eigenbasis [2], and will be referred in the following as the PCD condition. It guarantees the thermodynamic consistency of the models [15,16], which may be broken when some uncontrolled approximations are introduced [17]. Under this assumption coherences decay with time and are irrelevant in the steady state functioning of the device, contrary, for example, to externally driven devices [18] and systems including matter currents [19], where they may play an important role on the thermodynamic properties. When the PCD condition holds, the populations follow a continuous time Markov master equation [20], given in terms of the rates for the transitions between states, which are always allowed in both directions. In this case a thermal device implemented in a quantum system can be described within the framework of stochastic thermodynamics [21,22,23].
The analysis of the thermodynamic quantities can be realized at different levels of description [24]. From a macroscopic point of view, where the relevant quantities are the bath temperatures and the physical (total) heat currents, to a microscopic description that considers in addition the device structure. This latter perspective is more and more relevant as the advance of the experimental techniques allows for the design and the manipulation of the device. A prominent tool for this microscopic analysis is graph theory, where the stochastic master equation for the populations is represented by a graph. Schnakenberg theory [25] is a popular approach that gives a decomposition of the total entropy production based on a set of fundamental circuits in the graph. Basically, Schnakenberg applies Kirchhoff's current laws to reduce the number of terms appearing in the entropy production, which may be highly beneficial for optimization procedures. It has been used for example in linear irreversible thermodynamics [26] and in the study of steady-state fluctuation theorems [27,28]. This method does not intend to associate an entropy production with each circuit. In particular, the attempt to interpret individually each term in the decomposition may lead to apparent negative entropy productions, although this problem can be avoided by a convenient choice of the fundamental circuits [29,30,31]. However, it has been shown that the diagnosis of the machine performance greatly benefits from considering the thermodynamic analysis of not only the fundamental but all the possible circuits in the graph [32,33,34]. A convenient approach is then Hill theory [35]. Schnakenberg and Hill theory assign the same affinity to each circuit, but the latter considers all the possible circuits and leads to thermodynamically consistent entropy productions. Both methods coincide when the fundamental set of circuits contains all the possible ones.
In this paper we use Hill theory to fully characterize the two relevant quantities in the study of continuous absorption devices: the steady state heat currents and performance. Our aim is to find out under what conditions these quantities are as large as possible, i.e. what is the best construction of multilevel machines. Graph theory allows us to answer this key question from a very general perspective, looking only at the topological structure of the graph. Although we are motivated by the study of quantum models and in the following we will assume the PCD condition, our analysis also applies to classical stochastic models, including mesoscopic systems where the relevant degrees of freedom are identified by a coarse graining procedure [36,37]. In fact, the main advantage of this approach is that many properties of a device can be inferred from its graph representation irrespectively of its underlying, microscopic or mesoscopic, quantum or classical, realization.
It has been shown that systems with degenerate energy levels and driven by an external field may present a linear increment of the heat currents with the number of states [38,39]. Furthermore, two-stroke models in the quasi-equilibrium regime show an improvement in the performance with the number of levels [40]. However, using a particular construction of continuous absorption devices by merging three-level systems, Correa [41] found no changes in the performance and a fast saturation in the magnitude of the heat currents as the number of levels increases. Thus the question arises whether this limitation may be overcome by different designs of the absorption device.
We are interested in continuous machines that either extract energy from the coldest bath (refrigerators) or inject energy to the hottest bath (heat transformers). We do not consider devices designed for complicated tasks involving more than one target bath, although our procedure could also be applied to such systems. The best refrigerators and heat transformers should generally provide the largest possible heat currents and be also able of reaching the reversible limit for a particular set of the parameters. In order to identify them, we first justify that a machine coupled to three baths is capable of achieving the same currents than more complicated devices which consider additional heat reservoirs. As multilevel machines are composed by multiple circuits, our next step is to identify the optimal circuit to be used as building block. In general the magnitude of the heat currents increases with the transitions rates for any circuit. Hence, to elucidate the role of the circuit structure in the currents we set the rates to fixed values. Moreover, this condition avoids processes which prevents the machine from reaching the reversible limit when considering multiple circuits. The following step is to determine the graph structure leading to the largest heat currents considering optimal blocks. Finally, we relax the condition of fixed rates to improve the scaling of the currents with the number of levels without introducing harmful processes as heat leaks.
The paper is organized as follows: in section 2 we motivate the generic nature of our work by describing two different models of absorption devices which are represented by the same four-state graph. The master equation for all the quantum models used as illustration of the general results can be obtained using Appendix A with the Hamiltonians provided in Appendix B. We also introduce in section 2 the essential concepts of graph theory needed to characterize the heat currents associated with a circuit inside a general graph. This result allows us to relate each circuit to a thermodynamic mechanism and classify it attending to its contribution to the overall functioning of a device coupled to three baths. Although we have used previously the circuit decomposition in a different context, the analysis of the irreversible mechanisms arising in thermal devices indirectly connected to environments [34], we provide now a derivation of it using Hill theory in Appendix C. The differences between Hill and Schnakenberg decompositions are discussed and worked out for the four-state graph in Appendix D and Appendix E. In section 3 we analyze multilevel machines represented by a graph circuit. Explicit expressions for the scaling of the heat currents with the number of levels in the high and low temperature limits are provided in Appendix F. Machines represented by graphs with multiple circuits are studied in section 4. A simple example to illustrate the relation between the heat currents and the graph connectivity is presented in Appendix G. We draw our conclusions in section 5.

Motivation and background
We will motivate our approach by first considering two different models of absorption devices, both connected to three reservoirs at temperatures T c (cold), T h (hot) and T w (referred in the following as the temperature of the work bath, in analogy with devices driven by an external field), with T c < T h < T w . Depending on internal parameters, the devices can either work as heat transformers, transferring energy from the hot to the work bath, or as refrigerators, extracting energy from the cold bath assisted by the work bath. The first model is the two-qubit device [5] shown in figure 1(a). Each qubit is connected to a bosonic heat bath at temperatures T c and T h . The interaction between them is mediated by another bath at temperature T w . The state of this machine can be expanded in the product state basis |1 ≡ |0 h 0 c , |2 ≡ |0 h 1 c , |3 ≡ |1 h 0 c and |4 ≡ |1 h 1 c , with energies E 1 = 0, E 2 = ω c , E 3 = ω h and E 4 = (ω c + ω h ). When the system is weakly coupled to the reservoirs, the dynamics of the populations is described by a master equation where p i are the populations, W α ij ≥ 0 the transition rates associated with the coupling with the bath α, and W α ii = − j =i W α ji . For our purpose now is only important that the non-zero rates associated with the cold bath correspond to transitions 1 ↔ 2, 3 ↔ 4, with the hot bath to 1 ↔ 3, 2 ↔ 4, and with the work bath to 2 ↔ 3.
The second model is a photoelectric device [32,42,43,44] composed of two singlelevel, spinless quantum dots with energies E c and E h that can be both occupied at the same time. Each dot is connected to a metal electrode with chemical potential µ < E c < E h and temperatures T c and T h , see figure 1(b). We choose the same chemical potential in order to avoid introducing mechanical work and assume that neither the temperatures nor the chemical potential are modified by the interchange of electrons through the quantum dots. The system states are where now 0 α and 1 α are the number of electrons in the dot α. Transitions between the two dots are supported by an additional radiation source (for example the Sun in photovoltaic models [32,44]) at an effective temperature T w . Considering a weak coupling with the electrodes and a negligible line broadening of the energy levels, the device dynamics can be described by an stochastic master equation [45] in the form (1), but with transition rates determined by the particularities of the physical model under consideration. In the case of the absorption device with bosonic baths the rates are proportional to Planck distributions, while for the photoelectric device they are proportional to Fermi functions.
The relevant point for our analysis is that since the master equations have the same structure, the devices share several thermodynamic properties that stem directly from it. The different physics involved in each case is only reflected in the particular values of the transition rates. The master equation may be represented by a network, a weighted and labeled multi-digraph. However, as transition between states are always allowed in both directions, we will use a simpler representation consisting in a labeled graph, with vertices associated with the system states and undirected edges with the transitions [25,35]. When necessary, an arbitrary orientation can be assigned to the graph and a weight to each edge, given by the corresponding transition rate. For example, the representation of (1), denoted in the following by G 4 , is shown in figure 1(c). The circuits of G 4 , defined as a cyclic sequence of distinct edges, are displayed in figure 1(d), (e) and (f). Circuits C 1 and C 2 participate in different processes depending on their two possible orientations, referred as cycles. For example the cycle C 1 ≡ {1, 2, 3, 1} absorbs energy from the cold and work baths that is rejected to the hot bath, whereas the opposite cycle − C 1 ≡ {1, 3, 2, 1} absorbs energy from the hot bath and rejects it to the cold and work baths. In both processes there is a net exchange of energy with the three baths. The circuit C 3 involves only two baths (cold and hot) and in our models does not lead to any net exchange of energy. This is a consequence of having the same energy gap for transitions assisted by the same bath. However, in more general setups with different transition energies, E 34 = E 12 + ∆ and E 24 = E 13 + ∆ with E ji = E i − E j , there is a heat leak which increases with the energy shift ∆ [33].
The overall physical heat currentsQ α and the performance of the device are then the result of the interplay of the different mechanisms related to each circuit. In spite of the simplicity of the previous qualitative interpretation, the microscopic currentṡ q α (C ν ) corresponding to each circuit in the graph are not straightforwardly obtained from the physical currents. We introduce below the concepts of graph theory needed to characterize them.

Graph, circuits and steady state heat currents
For simplicity we consider systems with N states of energies E i , 1 ≤ i ≤ N, represented by a connected graph and coupled with thermal baths. The generalization for systems exchanging particles without involving any mechanical work, as the absorption device of figure 1(b), is straightforward. The system transitions may be coupled to one or several independent heat baths, each one in equilibrium at temperature T α , 1 ≤ α ≤ R. The system evolution is described by a master equation where p i is the normalized probability distribution to be in the state i, W α ij ≥ 0 is the transition rate from the state j to the state i due to the coupling with the bath α, and The transition matrix W, with elements W ij = R α=1 W α ij , is singular, which guarantees the existence of a non-trivial steady state solution of (2) and the conservation of the normalization. In addition, we assume that where k B is the Boltzmann constant. If the transition rates W α ij for j > i are known, the remaining rates can be determined by using (3) and (4). The master equation (2) is represented by a graph G(N, U) composed of N vertices and U undirected edges. Let x e be an edge in the graph, 1 ≤ e ≤ U. In the following x e will denote an edge oriented from vertex i e to j e , whereas − x e connects j e to i e , in both cases due to the coupling with the bath α e . Oriented edges are related to rate coefficients by W ( x e ) = W αe jeie and W (− x e ) = W αe ieje . An algebraic value A may be assigned to any oriented subgraph G s of G, composed of s ≤ U oriented edges x e [25], where, if the subgraph involves edges associated with the bath α, with e∈s,α the product over all the directed edges of G s corresponding to this bath, and otherwise A α ( G s ) = 1. Both A( G s ) and A α ( G s ) are positive real numbers. A maximal tree T µ , 1 ≤ µ ≤ N T , is a subgraph of G containing N − 1 edges without forming any closed path. The oriented subgraph T µ i is a maximal tree in which all the edges are directed towards the vertex i. A chord of a maximal tree is one of the U − N + 1 edges that are not part of it. The subgraph obtained when a chord is added to a maximal tree has only a circuit C ν , 1 ≤ ν ≤ N C . When removing the circuit from the previous subgraph, a collection of edges remains. Orienting them towards the circuit, a forest F β ν is found. The index β indicates that for a given circuit different forests can be found, resulting from different maximal tress. The number of maximal trees (N T ), circuits (N C ) and forests depend on the topological structure of the graph G. Each circuit C ν may be oriented in one of the two possible directions, leading to the cycles C ν and − C ν . Some examples are shown in figure 2. In Appendix C we use Hill theory to show that the steady state heat current associated with a circuit is given bẏ The factor D is calculated using The quantity D(G) > 0 increases with the complexity (both the number of vertices and edges). It is a factor which reduces the population in a circuit and therefore the corresponding heat currents when considering machines with an increasing number of them. The matrix W is obtained from the transition matrix W by replacing the elements of an arbitrary row by ones, whereas the matrix (−W|C ν ) is obtained by removing from −W all the rows and columns corresponding to the vertices of the circuit. Indeed, det(−W|C ν ) is the sum of the forests of C ν and can be thought of as an "injection of population" through edges not belonging to it. We have also introduced the cycle affinity associated with the bath α, and then the total cycle affinity is The quantity −T α X α ( C ν ) is just the net amount of energy interchanged between the bath α and the system when performing the cycle C ν . Notice that X α (− C ν ) = −X α ( C ν ) and hence each cycle is related to a process where some energy is either absorbed from or rejected to the bath. The circuit heat current (7) can be viewed as the result of the competition between the two cycles, described by As a consequence of (4) reflecting that the net energy exchanged by the system with the baths along a complete cycle is zero. Using it the following relation is found, and since the only contribution to the steady state entropy production is due to finiterate heat transfer effects, the circuit entropy production iṡ where the inequality is shown in Appendix C. These two last equations assure the consistency of the circuit heat currents and the entropy production with the first and second laws of thermodynamics. Finally, the total entropy production is given bẏ S = N C ν=1ṡ (C ν ) and the physical heat currents byQ α = N C ν=1q α (C ν ). They can be directly obtained from the transition rates by using (7), without determining the steady state populations. As an illustration of the circuit decomposition, the heat currents for the graph G 4 are worked out in Appendix D. Let us remind that other decompositions ofṠ are possible and we briefly discuss them in Appendix E.
The circuit heat currents (7) are homogeneous functions of degree 1 with respect to the transition rates, that is with σ > 0. Therefore, the currents can be always modified by changing the rates, provided that the assumptions to obtain the master equation remain valid. This property emphasizes the importance of the graph topology.

Classification of circuits
The contribution of each circuit C ν to the physical heat currents can be classified attending to their non-zero affinities X α : (i) X α ( C ν ) = 0 for all the baths. These circuits will be referred as trivial circuits, as they do not contribute neither to the steady state heat currents nor to the entropy production.
(ii) Condition (11) prevents any circuit from having only a non-zero affinity X α .
(iii) X α ( C ν ) = 0 only for two baths, α = α 1 , α 2 . Then there is only a net energy transfer between them, although other baths could participate in the cycle. Using (12) and (13), the following condition is founḋ Taking T α 1 < T α 2 , the heat currents verifyq α 2 (C ν ) > 0 andq α 1 (C ν ) < 0. Therefore the net heat current associated with these circuits always flows from the higher temperature bath to the lower temperature one. In the context of refrigerators and heat transformers these circuits are related to heat leaks that decrease the performance [33,34].
(iv) X α ( C ν ) = 0 for three baths, α = α 1 , α 2 , α 3 . They will be referred as three-bath circuits in the following. Equation (11) implies that, given a circuit orientation, two of the affinities and their corresponding heat currents must have the same sign.
Considering sgn(X α 1 ) = sgn(X α 2 ) = −sgn(X α 3 ) and using again (12) and (13), we obtainq The formalism applies also to circuits with non-zero affinities associated with more than three baths, but they are not relevant for our analysis.

Circuits in refrigerators and heat transformers
For simplicity we discuss now refrigerators, but the results are also valid for heat transformers. In general, the environment may be composed by the target coldest bath, a collection of sink baths with temperatures {T h,i } (where the surplus energy is rejected) and work baths with temperatures {T w,i } (supplying energy to complete the cycles). Let {X α,i } be the affinities of a particular circuit. Equation (7) implies that we can always find a hot and a work bath with temperatures and affinities given by , such that tuning their rate values (14) we obtain the same or larger heat currents than in the original system. Therefore we focus in the following on circuits and thermal machines coupled to three thermal baths with temperatures T c < T h < T w .
In the construction of the device we do not consider circuits with two edges associated with different baths connecting the same vertices, as it would lead directly to heat leaks (iii). To perform useful tasks we must include three-bath circuits (iv), which can be classified as: In cases (a) and (b) heat is simply transferred from the work to the cold bath, whereas the hot bath absorbs or gives up some energy. In (c) two different directions for the heat currents are possible: , which correspond to the conditions for the heat currents in heat transformers and refrigerators respectively. Therefore equation (17) settles the condition for the affinities in useful circuits. The particular working mode will depend on the system parameters.

Thermal machines represented by a circuit graph
In this section we analyze thermal machines that are represented by a circuit graph, G = C N , with N ≥ 3 states (vertices) and U = N undirected edges. We consider useful three-bath circuits for which (17) holds. Along this section we shall make explicit the circuit length (the number of states or edges) by the superscript N. In this case the physical and circuit heat currents coincide. From (10) we obtain , and then the physical heat currents are given byQ where Notice that the dependence on the arrangement of the edges in the circuit is contained in R( C N ) and the currents vanish for X( C N ) = 0. Using (11), the circuit affinity is rewritten as The device operating mode depends only on the parameter which is independent of the particular circuit orientation, and for X( C N ) = 0 results in When x < x r , the device operates as an absorption refrigerator whose coefficient of performance is The coefficient of performance reaches the Carnot value when x approaches to x r from below but at vanishing heat currents (X( C N ) = 0). When x > x r , the machine operates as a heat transformer with efficiency ] when x approaches to x r from above. In consequence, the device performance depends only on the circuit affinities X α ( C N ), irrespective of the value R( C N ), and they may be suitably tuned to reach the reversible limit for any graph circuit.

Circuit structure, performance and heat currents
In the following and without loss of generality, we choose a circuit orientation such that X( C N ) > 0. The affinities and the algebraic value A( C N ) depend only on the number of edges and their associated transitions rates. In particular, A( C N ) is the product of N transition rates. However, the factor D depends also on the arrangement of the edges through the oriented maximal trees in (8). The N T = N maximal trees are obtained by removing in each case one of the edges in the circuit. We denote by T j the maximal tree obtained by removing the edge starting in the state j. The term D(C N ) is the sum of N 2 terms A( T j i ), each one composed of the product of N − 1 transition rates.
Then, for a given circuit orientation, a path from i − 1 to i + 1 consists in two jumps either absorbing energy from or rejecting energy to the bath. In both cases the transition rate W α is the same. (b) When considering more general graphs, the PCD condition implies that the maximum number of edges connecting a state is six in a machine connected to three baths.
3.1.1. Dependence on the transition rates From the previous results for A and D and after a straightforward calculation, the heat currents are bounded by where W m is the minimum rate in A( C N ). As intuitively expected, increasing the lowest rates may result in larger heat currents for any circuit. The remaining question is then what kind of circuit shows the largest heat currents for a set of fixed transition rates.
In order to answer it, we assume in the following that the available resource in the machine design is a set of three undirected edges with fixed transitions rates, W α and W −α , associated the first with energy transfer to and the second with energy absorption from the bath α = c, w, h. This construction can always overcome complicated ones with more that three edges using a proper scaling of the rates, see (14). Besides, it implies fixed energy gaps |E ij | ≡ E α for transitions assisted by the same bath and: (i) When two edges, x i−1 and x i , connecting the state i are associated with the same bath, then W ( x i−1 ) = W ( x i ) for any of the two cycles as a consequence of the PCD condition, see figure 3(a).
(ii) The minimum number of edges required to construct a useful three-bath circuit is three, therefore E c + E w = E h as a result of (11) and (17). For simplicity we take Considering these points, any circuit must be constructed adding either two-edge sets {αα}, with α = c, w, h, or three-edge sets {cwh} to guarantee that the change of energy of the system in a complete cycle is zero. We denote by m = m c + m w + m h the number of two-edge sets in a circuit. Each one of them contributes with the product W −α W α to the algebraic value A( C N ), independently of the circuit orientation. Circuits constructed only by adding two-edge sets (N = 2m) are trivial circuits, X α ( C 2m ) = 0. The n = n + + n − three-edge sets {cwh} in a circuit contribute either with the product . Notice that when changing the circuit orientation to − C, n + and n − are interchanged. The smallest useful circuit is a triangle denoted by C 3 , see for example figure 1(d) and (e). Large circuits C N with N = 3n + 2m states are obtained adding additional two and three-edges sets to C 3 . Their smallest instances are shown in figure 4(a) and (d).  3.1.2. Circuit affinities and R( C N ) The circuit affinities are given by X α ( C N ) = (n + − n − )X α ( C 3 ) for a proper choice of the cycles. Both C N and C 3 have the same value of the parameter x, provided that n + − n − = 0, and then the performance of the circuit C N , given by (22) or (23), is necessarily equal to the performance of C 3 . In other words, for a fixed set of transition rates the circuit performance is independent of the number of edges.
The remaining question is whether larger circuits result in an increment of the magnitude of the heat currents with respect to C 3 . As X α ( C N ) increases at most linearly with N, the term depending on the affinities in (18) increases at most as N 2 , but only when NX( C 3 )/k B remains small. However, the increment of the affinities with the number of states is compensated by the factor R( C N ). As the number of terms in D grows quadratically with N, one would expect that in most cases R( C N ) decreases when adding new states and edges to the circuit. In fact, numerical evidence indicates that when adding two and three-edge sets to a circuit, R( C N ) decreases equal or faster than N −z , with z ≥ 1 for large enough N, see for example figures 4(b) and (e). Notice that we do not claim that R( C N ′ ) < R( C N ) for arbitrary values of N and N ′ subjected to the condition N ′ > N. Our statement only applies to the construction where the circuit C N ′ is obtained by adding two-edge and three-edge sets to C 3 , while keeping the edges and the orientation of the latter. Explicit expressions for R( C N ) in the high and low temperature limits supporting this result are given in Appendix F.
We have shown that typically the affinity term in (18) depends linearly on N whereas R( C N ) decreases faster than N −1 , and therefore in most cases the heat currents will decrease when adding additional edges to C 3 , see figures 4(c) and (f). An increment in the heat currents may be obtained at some extent by adding some three-edge sets when z < 2 while the circuit affinity remains small enough to grow quadratically with the number of states, as shown in figure 4(f) for intermediate temperatures. This improvement, although modest, may be relevant in situations where the heat currents are intrinsically small. Intuitively, increasing the circuit size implies the addition of states with larger energies and small populations except for specific values of the parameters. This small population makes harder closing the cycles and then effectively reduces the heat currents. Thus, the triangle C 3 is in general the optimal choice as building block for multilevel devices. This machine is one of the reference models used in quantum thermodynamics and it has been studied in both the cwh [1,3,4] and wch [23] configurations. Since X α ( C 3 cwh ) = X α ( C 3 wch ) for a proper orientation, the circuits show the same thermodynamic performance. Notice that A( C 3 cwh ) = A( C 3 wch ) but D(C 3 cwh ) = D(C 3 wch ). Using (18), the heat currents are related bẏ Q wch For high temperatures, y α ≡ exp[−E α /(k B T α )] ≈ 1, the arrangement of the edges in the circuit is irrelevant andQ wch α /Q cwh α ≈ 1. A different picture appears at low temperatures, For W c < W w , the ratioQ wch α /Q cwh α < 1 and for W c > W w ,Q wch α /Q cwh α > 1. Then the most favorable configuration corresponds to the lowest transition rate being associated with transitions from the ground state, by far the most populated in the low temperature limit.

Thermal machines represented by a graph with multiple circuits
We study now multilevel absorption machines with multiple circuits. We start by analyzing the relation between the heat currents and the performance of a circuit C ν in an arbitrary graph G, and the corresponding quantities for the (isolated) graph circuit C iso ν . To this end, we rewrite (7) aṡ Using this expression we find: The first result indicates that the magnitude of the heat currents associated with a circuit in a graph is always smaller than the one corresponding to the isolated circuit. It follows from (8) by noticing that the product between a term in the forest det(−W|C ν ) and a term of D(C iso ν ) gives the algebraic value of one of the oriented maximal trees of G. Therefore det(−W|C ν )D(C iso ν ) = N ′ T µ=1 i∈ν A( T µ i ), with i∈ν the summation over all the vertices of C ν and being the number of maximal trees involved N ′ T ≤ N T . The second result derives directly from (22) and (23) and indicates that the circuit performance is not modified when the circuit is included in an arbitrary graph.

General bound for the performance
A consequence of (ii) is that the device performance cannot exceed the corresponding to the circuit with the best performance. For example, let us consider a device working as an absorption refrigerator,Q c andQ w > 0. The coefficient of performance is given by whereq c (C ν ) is positive for the N ′ C circuits contributing to the cooling cycle, and negative for the N ′′ C − N ′ C "counter-contributing" circuits, corresponding for example to heat leaks and circuits with finite counter-currents which flow in directions against the operation mode [33,34]. The N C −N ′′ C trivial circuits are irrelevant in this discussion. In consequence, denoting by ε(C ν ) max the largest performance of a circuit in the graph, and the equality, ε = ε(C ν ) max , is reached when N ′′ C − N ′ C = 0 and ε(C ν ) = ε(C ν ) max for all the circuits. In particular, ε = ε C only if all of them achieve the Carnot performance for the same value of the affinity. A similar analysis applies to the device working as a heat transformer. Therefore, with regard to the performance, optimal multilevel machines are represented by graphs without "counter-contributing" circuits. We will impose this condition in the design of the optimal graph.

Graph topology and heat currents
The magnitude of the physical heat currents is determined by the graph topology and the value of the transition rates. We first explore the graph topology of an arbitrary graph with the only restriction that two vertices can be connected by just one edge (see section 2.3).
In generalQ α = N C ν=1q α (C ν ) increases with the number of positive contributing circuits N ′ C ≤ N C , which operate in the same way as the entire device. However, this increment may be hindered by the unavoidable decrease in D(G) −1 det(−W|C ν ) when adding new states and edges to a graph. The factor D is the sum of NN T terms. For circuits of length L, det(−W|C L ν ) is the sum of det( A|C L ν ) terms. In this expression the submatrix ( A|C L ν ) is obtained by removing from A all the rows and columns corresponding to the vertices of the circuit C L ν . The matrix A is calculated by replacing the diagonal elements a ii of the adjacency matrix A (see for example [46]) by the vertex degree of the corresponding state i. The non diagonal elements are a ij = 1 when states i and j are connected by an edge, and a ij = 0 otherwise. Therefore, when attending to the number of terms, the magnitude of the heat currents resulting from the positive contributions of N L ≤ N ′ C circuits of length L is related to the topological parameter where λ(C L ν ) ≡ det( A|C L ν )/N T , with N −1 T ≤ λ(C L ν ) < 1. The ratio λ may depend on the position of the circuit in the graph and in general λ(C L ′ ν ) < λ(C L ν ) when L ′ > L. Although τ L can be readily calculated, we found that the upper bound τ b L incorporates the relevant information about the graph topology. In particular, it makes clear the relevance of the graph connectivity: favorable graphs consist in as many small positive contributing circuits as possible (that is, avoiding heat leaks and another negative contributions), built with the smallest possible number of states, implying a large graph connectivity. This dependence on the graph topology is weighted by the transition rates. For high temperatures all circuits participate in the heat currents. However, only small circuits including the ground state will contribute significantly in the low temperature regime, independently of the total number of circuits in the graph.

Graphs constructed by merging triangles
The optimal choice for the building block is the triangle, the smallest possible contributing circuit as described before. We consider that all the triangles have fixed energy gaps for transitions assisted by the same bath. This is a necessary condition to achieve the maximal possible connectivity because otherwise adjacent triangles cannot share any edge. In order to analyze the dependence on the graph topology we consider now the more restrictive condition of fixed transition rates for each bath. This assumption will be relaxed latter. Moreover, we assume the PCD condition, that implies now that the maximum vertex degree in the graph is six, i.e. each state may be connected at most to another six ones, see figure  3(b). As a consequence, all the constructed graphs are planar and τ b 3 incorporates the relevant topological information. The number triangles is easily accessible by using the adjacency matrix, N 3 = Tr{A 3 }/6, where Tr{} denotes the trace. The graph with the largest connectivity compatible with our restrictions is denoted by G B 4 . It is constructed using B units of two triangles sharing one edge, for example one associated with the work bath, see figure 5(a). We consider square graphs with 1, 4, 9, . . . units, being the smallest instance G B=1 4 ≡ G 4 . By construction, the two configurations of the triangle, cwh and wch, are present. Besides, many other circuits can be identified. For example {(0, 0), (0, 1), (1, 1), (2, 0), (1, 0), (0, 0)} is a circuit C 3+2m h with m h = 1, and {(0, 0), (0, 1), (0, 2), (1, 1), (2, 0), (1, 0), (0, 0)} a circuit C 3+3n + with n + = 1. All of them follow the same operation mode. There are also many trivial circuits, for example {(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)}. This construction is optimal with respect to the performance because it can attain the reversible limit as there are not "countercontributing" circuit, see the discussion for C 3 in section 2.
We also consider two subgraphs of G B 4 for comparison purposes. The first one is a row of this units, denoted by G B 4L , which represents the absorption device studied in [41]. The second one is obtained considering only a row and removing the upper hot edges. We use this graph, denoted by G B 3 , to compare τ b L with other measure of the graph connectivity in Appendix G. Figure 5(b) shows the parameter τ b 3 for G B 4 , G B 4L and G B 3 , considering only complete units in each case. For a given number of states, larger values of τ b 3 correspond to larger number of circuits and therefore to a larger connectivity. When the number of states increases the parameter τ b 3 saturates to a different constant value in each case. This is reflected in the physical heat currents shown in figure 5(c) for different bath temperatures. This saturation is due to the difficulty of exploring big circuits or those which are distant from the ground state in complex graphs. The simple picture based on the parameter τ b 3 is weighted by the transition rates. For decreasing bath temperatures, all the currents converge to the same result, independently of the number of circuits, since only the triangle including the ground state contributes significantly to them.
In summary, given a set of transition rates and a number of levels, the best topology corresponds to the most connected planar graph G B 4 . This construction only contains trivial and positive contributing circuits and provides in general the largest heat currents for fixed rates.

Transition rates and heat currents
We have shown that for fixed transition rates the heat currents saturate to a constant value when increasing the number of states. To overcome this limitation, we now consider a graph with the optimal topology given by G B 4 and relax the condition on the rates but keeping fixed energy gaps. The circuit affinities and then the performance are not modified. Considering (4), all the transition rates must be taken as sW ±α , with s ≥ 1, and W ±α the smallest rate. As discussed for circuit graphs, increasing s will lead to larger heat currents.
In particular, we analyze the construction shown in figure 6(a), denoted by G B HO , which has a simple physical implementation as discussed below. Seeking a measure of the graph connectivity when the transition rates increase with s, and in analogy with the adjacency matrix, we define A ′ with elements a ′ ij = s when states i and j are adjacent with transition rates sW ±α , and we denote its spectral radius as ρ(A ′ ), see Appendix G. When incorporating additional building units into the graph, the spectral radius defined in this way increases nearly linearly with the number of states, see figure 6(b).
The graph G B HO , allowing an infinity number of building blocks, represents the master equation of a device composed of two harmonic oscillators [5]. Each oscillator is connected to a thermal bath at temperatures T c and T h . The coupling operators areŜ c =â c andŜ h =â h (see Appendix A), beingâ α the annihilation operator of the oscillator coupled to the bath α. A third bath at temperature T w is coupled to the system through the operatorŜ w =â † câ h . For simplicity we assume a very large value of T w , a regime for which the heat currents can be easily calculated. Figure 6(c) shows the heat currents as a function of n c n h , which gives a rough estimation of the number of states populated and then of the effective graph size, calculated for increasing bath temperatures. In this expression n α is the average number of excitations in the oscillator α = c, h. When the temperature increases, larger areas of the graph are populated involving a larger number of circuits, which results in an increment of the magnitude of the heat currents. This example illustrates that given a machine with the optimal topology, the rates can always be carefully designed to achieve larger currents without diminishing the performance.

Conclusions
We have determined the steady state heat currents associated with all possible circuits in the graph representing the master equation of multilevel continuous absorption machines. Each circuit is related to a thermodynamically consistent mechanism in the device functioning. Although the number of circuits may be very large when increasingly complex graphs are considered, efficient standard algorithms, which scale as N C (N +2U) [47], can be used for determining them. For example, in the graphs studied in previous sections U increases linearly and N C quadratically with the number of states and the computational cost scales as N 3 . The main result of the decomposition is an equation for the circuit heat currents depending only on the transition rates, without any prior knowledge of the steady state populations. This expression allows us to analyze the two relevant quantities for refrigerators and heat transformer, the magnitude of the physical heat currents and the performance. We focus on devices coupled to three baths, since they can provide the same currents than more complicated setups.
In order to elucidate the role of the graph topology in the thermodynamic properties, we have analyzed machines constructed by a fixed set of transition rates. In devices represented by a single graph circuit, the performance depends only on the circuit affinities, which can be tuned to reach the reversible limit, and the magnitude of the heat currents decreases in general with the number of states. Then the simplest graph, a triangle, leads to the largest heat currents in most cases and is the proper building block for optimal multilevel machines.
When considering generic devices, we have found that the performance of the device cannot exceed the corresponding to the circuit with maximum performance. Besides the magnitude of the heat currents is described by a topological parameter that increases with the graph connectivity. As a consequence, if the construction of larger graphs including additional circuits presents a limited connectivity, then the magnitude of the resulting physical heat currents saturates to a constant value, which is different for different constructions. We use triangles with fixed energy gaps for transitions assisted by the same bath to construct the graph with the largest possible connectivity, denoted by G B 4 . This is a planar graph containing neither heat leaks nor "counter-contributing" circuits.
The assumption of a fixed set of transition rates can be relaxed. We give the necessary condition to improve the currents without modifying the performance. We provide an example using a system of harmonic oscillators. In this case the magnitude of the heat currents increases almost linearly with the effective size of the graph, determined by the achievable range of temperatures. An interesting question is whether there are other physical feasible implementations leading to a faster than linear dependence of the currents on the number of states.
The circuit decomposition could be employed in other different scenarios, from the study of heat transport through quantum wires to the analysis of machines designed for complicated tasks involving more than three baths. Besides, our formalism also applies to the case of reservoirs exchanging both energy and particles with the system, and even to periodically driven machines. The only condition required is that the population and coherence dynamics are decoupled in a certain basis. However, this is not always possible, as for example in weakly driven systems. Finally, it is worth to mention that the study of four-stroke many-particle thermal machines has been recently addressed in [48], the analysis of their continuous counterparts is another interesting issue we can explore in the future by using the circuit decomposition. We expect these findings will help in the experimental design of absorption devices.

Acknowledgments
We thank L. Correa and A. Ruiz for useful discussions.
Appendix A. Transition rates for quantum systems weakly coupled with thermal baths In this appendix we describe how to calculate the transition rates W α ij in the master equation (2) for a quantum system with HamiltonianĤ S = N i=1 ω i |i i|, and coupled with R bosonic baths at temperatures T α . We assume that the PCD condition holds. The total Hamiltonian readŝ whereĤ α are the bath Hamiltonians and the coupling terms are given bŷ withŜ α andB α a system and a bath operator respectively. The rates γ α determine the time scale of the system relaxation dynamics. Finally, the system operators in the coupling terms arê We consider the following assumptions: the system is weakly coupled with the environments, Then the Born-Markov and the rotating wave approximation applies and the master equation for the populations of the eigenstates ofĤ S [2] is given by (2) with transition rates (i < j) The functions Γ α only depend on bath operators whereB α (t) = exp(iĤ α t/ )B α exp(−iĤ α t/ ) andρ α denotes the bath thermal state. We will consider bosonic baths of physical dimensions d α and coupling operators The summation is over all the bath modes of frequencies ω µ and annhilation operatorsb µ . With this choice the rates Γ α ±ω are [2] The frequency ω 0 depends on the physical realization of the coupling with the bath. The condition (3) derives now directly from the conservation of the normalization of the system density matrix. Besides, the Kubo-Martin-Schwinger relation in (A.6) implies (4).

Appendix B. Quantum implementation of the graphs
We introduce here a possible quantum physical realization of the graphs described in the main text by specifying their Hamiltonians and coupling operators. Considering bosonic heat baths, the results of Appendix A can be used to obtain the corresponding transition rates. In all cases ω c + ω w = ω h .
Appendix C. Hill theory and the steady state heat currents We apply Hill theory [35] to obtain (7). The starting point is the steady state probability of finding the system in the state i [25,35] with D given by (8). Introducing the steady state fluxes along a directed edge and the corresponding affinities the total steady state entropy production is given by [22,25] where the orientation of each edge is arbitrary. Introducing the populations in the product between fluxes and affinities where µ∈Me denotes the summation only over the maximal trees for which x e is a chord, since otherwise the term between brackets is zero. The product W αe jeie A( T µ ie ) is no more than the algebraic value A of the oriented subgraph T µ ie + x e , composed of the maximal tree T µ ie and its chord x e . Then the entropy production (C.4) can be written asṠ Each term between brackets is only related to a circuit oriented in the two possible directions, C ν and − C ν , associated with T µ ie + x e and T µ je − x e respectively. When removing these two cycles from the corresponding subgraphs, the same forest F β ν remains, see for example figure C1(a) and (b). Using this result and the properties of A, each term in (C.6) can be written as The number of such terms with the same forest F β ν equals the number of edges of the circuit C ν , as shown in figure  C1(c). Next we introduce the cycle affinity (10), X( C ν ) = e∈ν X( x e ) with e∈ν the summation over all edges of C ν , to obtaiṅ where β∈ν denotes the summation over all the different forests associated with C ν . This expression can be further simplified applying the matrix-tree theorem [49], β∈ν A( F β ν ) = det(−W|C ν ). The flux associated with each cycle is The cycle affinities associated with each bath (9) are X c ( C 1 ) = E 21 /T c , X w ( C 1 ) = E 32 /T w , and X h ( C 1 ) = E 13 /T h , where E ij = E j − E i . The contribution of the forests is det(−W|C 1 ) = W h 24 + W c 34 , from which the cycle flux is given by where D(G 4 ) is determined by using (8). Then the circuit heat currents areq c (C 1 ) = E 12 I( C 1 ),q w (C 1 ) = E 23 I( C 1 ) andq h (C 1 ) = E 31 I( C 1 ). The consistency of the circuit currents with the first lawq c (C 1 ) +q w (C 1 ) +q h (C 1 ) = 0 follows from E 12 + E 23 + E 31 = 0. The cycle affinity (10) is X( C 1 ) = E 21 /T c + E 32 /T w + E 13 /T h from which the circuit entropy production can be determined with (C.9). A similar procedure can be used in order to obtain the quantities associated with the circuit C 2 .
For the circuit C 3 we denote by C 3 the cycle {1, 2, 4, 3, 1}. Now A c ( C 3 ) = W c 21 W c 34 , A w ( C 3 ) = 1 (there is not any edge associated with the work bath) and A h ( The cycle affinities associated with each bath are X c ( C 3 ) = (E 34 − E 12 )/T c , X w ( C 3 ) = 0 and X h ( C 3 ) = (E 13 − E 24 )/T h . Notice that (E 13 − E 24 ) = −(E 34 − E 12 ). When the transition energies are equal, E 34 = E 12 and E 24 = E 13 , all the affinities are zero. The circuit C 3 involves all the graph vertices and therefore there is not any forest associated with it. Then (−W|C 3 ) is an empty matrix and det(−W|C 3 ) = 1. The cycle flux is given by and the circuit heat currents byq c ( .

Appendix E. Other decompositions of the entropy production
There are several possible decompositions of the steady state entropy production in terms of circuits. Schnakenberg [25] designed a method based on the identification of a set of U − N + 1 fundamental circuits. The circuits are determined by choosing an arbitrary maximal tree and adding each one of its chords. Taking a particular orientation for the circuits, a set of fundamental cycles is found. The total steady state entropy production is thenṠ , where x ν is the chord giving the circuit C ν and J( x ν ) the corresponding flux. The previous decomposition is simple and specially relevant when U − N is small. However, it is not unique, since it depends on the choice of the maximal tree, and some terms in the sum may be not positive definite, which discards a possible consistent thermodynamic interpretation of each circuit contribution. Besides the evaluation of J( x ν ) requires the calculation of the steady state populations.
As an example we apply Schnakenberg method to the four-state model. The procedure requires an arbitrary set of fundamental circuits of the graph G 4 . We choose the maximal tree shown in figure 2(a), which has two chords, {1, 2} and {2, 4}. By adding the chord {1, 2} the circuit C 1 is obtained. Choosing an arbitrary orientation, for example C 1 as in figure 2(d), the directed chord x 1 goes from state 1 to state 2. In this decomposition the flux associated with each cycle is taken as the corresponding to the directed chord (C.2), J( x 1 ) = W c 21 p s 1 − W c 12 p s 2 . The cycle affinity is defined by (10) and was calculated in Appendix D, X( C 1 ) = E 21 /T c + E 32 /T w + E 13 /T h . When adding the chord {2, 4} we obtain the circuit C 2 . Choosing the orientation {2, 3, 4, 2}, the directed chord x 2 goes from state 4 to state 2. The flux is J( and the affinity X( C 2 ) = E 43 /T c + E 32 /T w + E 24 /T h . Then the entropy production iṡ S = J( x 1 )X( C 1 ) + J( x 2 )X( C 2 ) . (E.1) The cycles { C 1 , C 2 } are the elements of one of the possible fundamental sets of G 4 . Notice that for our choice the circuit C 3 is not involved. A related decomposition is obtained by the algorithm of Kalpazidou [29,30]. For systems showing dynamical reversibility the algorithm leads to a Schnakenberg decomposition with a clever choice of the fundamental set of cycles, such that all the terms in the sum are positive. Therefore a positive entropy production can be assigned to each cycle, which is required in many applications [50,51]. The algorithm is based on choosing an orientation for the graph such that all the fluxes (C.2) for the directed edges are positive. Next a cycle is identified and the entropy production J min ( x ν )X( C ν ) > 0 is assigned to it, where J min ( x ν ) is the smallest flux associated with an edge of C ν . Then J min ( x ν ) is subtracted to each flux in the cycle to obtain a new flux field and the process is repeated for new cycles until a fundamental set is completed [50,51]. For example, let us assume parameter values for which the two triangles of the four-state model work as refrigerators. Then the fluxes along x 1 , x 2 , x 3 (from 3 to 4), x 4 (from 3 to 1) and x 5 (from 2 to 3) are positive. With this orientation only the cycles C 1 and C 2 appear in the directed graph and the entropy production can be written as (E.1), where the two terms are guaranteed to be positive. If we modify the system parameters such that the circuit C 2 works as a heat transformer but the overall device remains working as a refrigerator, the total entropy production can still be determined using (E.1), but the positivity of each term is not guaranteed. Now the fluxes are positive along the edges x 1 , − x 2 , − x 3 , x 4 and x 5 . Only the cycles C 1 and C 3 remain with this graph orientation and the algorithm of Kalpazidou leads tȯ S = J( x 5 )X( C 1 ) + J(− x 2 )X( C 3 ) . (E.2) However, in this expression the contribution of each mechanism, refrigerator (C 1 ), heat transformer (C 2 ) and heat leak (C 3 ) could not be isolated.
Appendix F. R( C N ) in the high and low temperatures limits In the high temperature limit, y α ≡ exp[−E α /(k B T α )] ≈ 1, the transition rates satisfy W −α ≈ W α , leading to vanishing affinities and heat currents. Now A( T j i ) is in a good approximation independent of the orientation, what considerably facilitates the calculations to obtain with r α = n + m α and r c + r w + r h = N. In this limit R( C N ) decreases quadratically (z = 2) with N, except when one or two of the terms r α /W α are much larger than the others and r α remains constant when increasing N, which can only happens adding two-edge sets. In this limit R( C N ) decreases as N −1 for small enough values of N.
In the low temperature limit, y α ≪ 1 and W −α ≪ W α . Again this limit implies vanishing heat currents. The cycle algebraic value is proportional to the small factors y α , A( C N ) ∝ α y uα+u ′ α α , where u α and u ′ α are the number of W −α transitions before and after the highest-energy state respectively. Besides, the largest contributions to D comes from two terms that include the lowest number of rates W −α , A( T h−1 Using this result the physical heat currents are given bẏ where 1 ≤ K ≤ 2 and K = 3N C /(2N C + 1) in the high temperature limit. For this graph λ(C ν ) = 1 3 and τ 3 = τ b 3 /3. The parameter τ b 3 as a function of the number of edges is shown in figure G1(b). The spectral radius ρ(A), defined as the largest eigenvalue of the adjacency matrix of the unweighted graph [52], is also shown. The spectral radius is a measure of the graph connectivity which increases monotonically with the number of edges. However, it does not reflect the decrease in the heat currents each time a pendant edge is added to the graph, see figure G1(c). An increment in the total heat currents is only found when a new triangle is completed, saturating to a constant value when the addition of new circuits does not improve significantly the graph connectivity. This behavior is well described by τ b 3 .