Emergence of energy-avoiding and energy-seeking behaviours in nonequilibrium dissipative quantum systems

A longstanding challenge in nonequilibrium thermodynamics is to predict the emergence of self-organized behaviours and functionalities typical of living matter. Despite the progress with classical complex systems, it remains far from obvious how to extrapolate these results down to the quantum scale. Here, we employ the paradigmatic master equation framework to establish that some lifelike behaviours and functionalities can indeed emerge in elementary dissipative quantum systems driven out of equilibrium. Specifically, we find both energy-avoiding (low steady dissipation) and energy-seeking behaviours (high steady dissipation), as well as self-adaptive shifts between these modes, in generic few-level systems. We also find emergent functionalities, namely, a self-organized thermal gradient in the system's environment (in the energy-seeking mode) and an active equilibration against thermal gradients (in the energy-avoiding mode). Finally, we discuss the possibility that our results could be related to the concept of dissipative adaptation.


II. INTRODUCTION
When matter is driven far from thermal equilibrium and dynamically dissipates energy, self-organized states and processes may emerge. Examples are as diverse as tornadoes and living cells. The kind of self-organization that has led to the transition from inanimate to living matter may look distinctive perhaps in part due to the emergence of seemingly purposeful functionalities and finely-tuned adaptive behaviours [1,2]. We are thus left with the challenge of looking for general principles behind the emergence of exceptional behaviours that look like those found in living organisms [3]. From a fundamental viewpoint, this can be significant for a better understanding of life as we know it. For instance, because the notion of natural selection does not apply to situations prior to the existence of the first living cells, one may turn to nonequilibrium statistical mechanics [4] and thermodynamics [5] as attempts to generalize Darwinian evolution. From a practical viewpoint, searching for physical prin- * Electronic address: valente.daniel@gmail.com ciples may allow us to extend to diverse systems the selforganized functionalities typically found in living matter. Studies with classical many-body systems have recently made progress in this direction [2,3,[6][7][8][9][10][11][12][13][14][15][16], showing, in particular, the emergence of energy-seeking [7,10] and of energy-avoiding behaviours [10,12,15], of dynamical self-healing [7,9,12], and of self-adaptation under certain changes in the environment [6,9,12,14]. While, on the one hand, we may expect emergent lifelike properties to require some degree of complexity and many-body systems, on the other hand, in our search for principles we usually try to simplify things as far as possible, keeping the essential ingredients only. Finding just the right balance is all but a trivial task: quoting Schwille, how simple could life be? [1].
Here, we look for the simplest possible scenarios where those lifelike self-organized behaviours and functions can emerge, with special interest in extending them to the quantum scale. We consider the dissipative dynamics of generic driven quantum systems, as described by the wellestablished Markovian quantum master equation. We show self-organized energy-avoiding behaviour (leading to low dissipation of energy in the steady state) and energy-seeking behaviour (leading to high dissipation) in three-level systems. We also show a self-adaptation between these two modes in a four-level system, depending on a change of stimulus. We find that, when the fourlevel system is in its energy-seeking mode, it can build a self-organized thermal gradient in its own environment. In its energy-avoiding mode, by contrast, the four-level system is shown to act against thermal gradients, actively seeking to equilibrate its own environment. Finally, we sketch a plausible link between our results and the notion of dissipative history, as conveyed by the concept of dissipative adaptation [2,3,17]. Throughout the text, we present biologically-inspired processes in order to support our hypothesis that some lifelike behaviours can emerge in far from equilibrium matter even in the absence of natural selection.

A. Energy-avoiding behaviour
Living organisms have evolved to keep themselves structurally stable and functional. To maintain stability, one way used by living systems is to avoid absorbing en-arXiv:2201.04496v1 [quant-ph] 12 Jan 2022 ergy of specific types. Bones, shells and cell walls are specially adapted to guarantee mechanical stability against wounds, just as photosynthetic pigments are well-suited to avoid the damaging absorption of ultraviolet light. From a broad thermodynamic perspective [2,3], energyavoiding structures may self-assemble when a physical system finds stable states achieved by means of a history of nonequilibrium work absorption and dissipation. These stable states have, in turn, very low susceptibility to further absorption from that specific nonequilibrium source. The system then becomes finely tuned to avoid absorbing more of that specific nonequilibrium energy input that drove the transition in the first place [10,12,15]. This is the energy-avoiding mechanism, which can also be deemed a negative feedback loop [6,8].
Here, we look for the emergence of this type of energyavoiding mechanism, as described above, in the dynamics of an elementary physical system. Specifically, we search for an open system that departs from thermal equilibrium and, once driven by a nonequilibrium source, transiently absorbs enough energy so as to finally stabilize in a nonequilibrium steady state characterized by a strong reduction in the amount of energy steadily absorbed and dissipated from that given drive. We find this energyavoiding behaviour emerging in a lambda (Λ) three-level system. Let us label the energy eigenstates as |a , |b and |e , in increasing order, as illustrated in Fig.(1)(a). We suppose that, initially (at times t < 0), the Λ system is in thermal equilibrium with an environment at temperature T . A nonequilibrium energy source, more specifically, a classical resonant monochromatic light beam, is suddenly turned on and kept constant (at t ≥ 0). As a result, detailed balance is broken, and the system undergoes a dissipative dynamics towards a self-organized nonequilibrium quantum state. A self-organized state here means an atypical (exceptional) stationary quantum state (as compared with the thermal equilibrium state of the same system), achieved autonomously (i.e, with no pre-determined goal encoded in the energy source). To emphasize the energetic exchanges driving the dynamics, we recast the master equation as where P is the absorbed power (at the origin of broken detailed balance), E em = E e − E m is the energy gap between the states |e and |m = |a , |b , J em is the dissipation current in the transition from |e to |m , and p e (t) ≡ ρ ee (t) is the excitation probability (with ρ ij ≡ i|ρ|j , and ρ being the density operator of the system). In turn, the state of this nonlinear, bistable system strongly affects back how energy will be absorbed. That is, the power absorption, depends on the state of the system, specially on the imaginary part of the quantum coherence ρ ae between the states |a and |e , whereas the dissipation rates, arise from stochastic transitions, involving the populations of all states, namely, |a , |b , and |e . In Eq. (2), Ω is the system-field coupling strength (aka the Rabi frequency, proportional to the amplitude of the incoming field). In Eq.(3), κ is the spontaneous emission rate, and n em describes a Bose-Einstein distribution. This closes a nonlinear (state-dependent) feedback loop, where energy consumption leads to a dynamical change, which strongly affects back how much energy is absorbed. At all times, J em provides a measure of broken detailed balance (see further details in the Methods).
After the transient absorption peak, the system achieves the (nonequilibrium) stationary state with equality p b (∞) = 1 holding at T = 0, where p b (t) ≡ ρ bb (t) is the population of state |b . We call this a self-organized state, since we are considering E b > E a (which means that Eq.(4) is quite far from the thermalequilibrium values for the populations, namely, p b → 0 and p a → 1, in the T → 0 limit). See Fig.(1)(b). In the stationary regime, the self-organized nonequilibrium state makes both P and the total dissipation |J|, where J ≡ m J em , to converge towards very low values, as we see in Fig.(1)(c). In other words, the system seems specially adapted to avoid absorbing (hence to avoid dissipating) more of that incoming light that drove the transition in the first place. We note that, at precisely T = 0, detailed balance is reestablished in the steadystate, P (∞) = J ea (∞) = J eb (∞) = 0, although a small power consumption, P (∞) = −J(∞) ≈ E ea κn eb , still holds at very low temperatures (n eb → 0). Fig.(1)(c) also shows that higher coupling strengths Ω accelerate the system's adaptation towards this energy-avoiding behaviour of low steady P . In Fig.(1)(b) and (c), we compare the zero-temperature limit, T = 0, with the finite temperature T = 0.1E ea /k B , where k B is the Boltzmann constant. This shows that the energy-avoiding behaviour is homeostatically robust against thermal fluctuations, at low-enough temperatures.

B. Energy-seeking behaviour
Living cells have also developed a sophisticated network of processes, called metabolism, characterized by the capacity of the cell to absorb specific types of nonequilibrium energy available in its environment, and then use it to maintain its (nonequilibrium) living state. Coming back to our broad thermodynamic perspective, a generic fluctuating physical system driven far from equilibrium may also dynamically (dissipatively) self-adjust towards an exceptional stationary state that seems specially adapted to maintain high rates of energy consumption from the specific drive that led the system to that Energy-avoiding behaviour in a threelevel system. (a) Λ system, with eigenstates |e , |b and |a (black horizontal bars), with corresponding eigenenergies Ee E b > Ea. The spontaneous emission rate of each transition is κ (dashed arrows), and the resonant drive with the coupling strength Ω induces transitions between states |a and |e (full arrow). Energy exchange rates are characterized as power consumption P and dissipation J. (b) Departing from thermal equilibrium (at times t < 0), the pump is suddenly turned on and kept constant (at t ≥ 0), with Ω = 0.5κ, so the state of the system asymptotically self-organizes towards a nonequilibrium state, with populations p b (∞) → 1 p b (0) and pa(∞) → 0 pa(0). In the meantime, the system undergoes a finite excitation probability, pe(t) > 0. (c) In the transient regime, both the absorption P (t) (black) and the dissipation J(t) (blue) are appreciable (both in units of Eeaκ, where Eea = Ee−Ea). Asymptotically, the stationary self-organized nonequilibrium state lowers the power consumption and the absolute value of the dissipation, P (∞) = |J(∞)| → 0. In (b) and (c), we compare T = 0 (solid curves) with the finite (low) temperature T = 0.1Eea/kB (dashed curves).
nonequilibrium state in the first place [7,10]. Such an emergent energy-seeking behaviour can also be regarded as a positive feedback loop [6,8].
We find this energy-seeking behaviour in a three-level quantum system, but in V configuration. We label the energy eigenstates as |a , |b , and |g , in decreasing order (see Fig.(2)(a)). Let us consider that, at times t < 0, the V system is at thermal equilibrium. By pumping one of its transitions with a classical resonant monochromatic light at t ≥ 0, the V system undergoes a transition to a self-organized state, which seems to be exceptionally adapted to keep absorbing more energy from the nonequilibrium source. That is, the population p a (∞) ≡ ρ aa (∞) increases with respect to its equilibrium value (see Fig.2(b)), and a quantum coherence between the states |a and |g , namely ρ ag (∞), is built up (see Methods), as if the system was purpose-fully trying to keep high stationary power consumption. In turn, the system dissipates high amounts of energy in a steady manner, maintaining its nonequilibrium state. Figure (2)(c) shows the dynamics of the absorption power P and the dissipation J as functions of time, at distinct coupling strengths Ω. The higher the Ω, the more power is steadily absorbed from light and the higher is the dissipation. Similarly to the energy-avoiding mode, the energy-seeking behaviour here is reasonably resilient against thermal perturbations, as shown in the comparison between T = 0 and T = 0.3E ag /k B . Eg. The spontaneous emission rate of each transition is κ (dashed arrows) and the resonant drive induces the coupling between states |a and |g with strength Ω (full arrow). Energy exchange rates are characterized as power consumption P and dissipation J. (b) Departing from thermal equilibrium (at times t < 0), the pump is turned on and kept constant (at t ≥ 0), with Ω = 0.5κ, so the state of the system asymptotically self-organizes towards a nonequilibrium state where the populations satisfy pa(∞) > pa(0) and pg(∞) < pg(0). (c) After a transient regime, the stationary self-organized state maintains high power consumption and dissipation, P (∞) = |J(∞)| > 0 (both in units of Eagκ, where Eag = Ea − Eg), which increase with Ω. In (b) and (c), we compare T = 0 (solid curves) with the finite (low) temperature T = 0.3Eag/kB (dashed curves).

C. Self-adaptative shifts
Can a few-level quantum system present a versatile behaviour where it autonomously shifts from the energyseeking to the energy-avoiding mode and back? Selfadaptive shifts in a system's behaviour [9,12,14] represent a close analogy with the highly selective, switch-like response of receptors in living organisms [6]. We implement such a self-adaptive shifting behaviour in a four-level system with a diamond energy structure (♦). Our motivation is that the ♦ system integrates, in some sense, features from the Λ and the V energy-level structures. We set out to analyze the change in the system's response under two distinct stimuli. To that end, we compare the case where the power source drives the transition between states |a and |e (Λ-type) with the case where it drives the transition between states |g and |a (V -type). Here, we have labeled the energy levels in increasing order as |g , |a , |b and |e (see Fig.(3)).
We find that, as expected, when the Λ-type transition is driven, the energy-avoiding behaviour emerges (see Fig.(3), panel (a)), whereas driving the V -type transition leads to the energy-seeking behaviour (as in Fig.(3), panel (b)). A remark must be made. In contrast to the results shown in Fig.(1), we see directly from the blue curves in Fig.(3)(a) that the absorbed power stabilizes at finite values, instead of dropping to (or close to) zero. This can be seen in Fig.(3)(a), where we plot P as a function of time. The reason is the finite lifetime of state |b in the ♦ configuration. By decaying towards |g , the system opens a pathway for a steady non-vanishing dissipation to the environment. To better clarify that this can indeed be characterized as an energy-avoiding behaviour, we show, in addition to the very low values of power absorption, that the population of state |a decreases with respect to its thermal equilibrium value. As we can see from the inset of panel (a) in Fig.(3), stronger drives tend to empty state |a more strongly, meaning that the system behaves as if it was purposefully avoiding to absorb energy from the driving light. Also, we see that the increase in Ω ae leads to a slight increase in P . In the energy-seeking mode, Fig.(3)(b), we find the opposite behaviour, namely, a stationary increase in the population of state |a with respect to the equilibrium value. This is shown in the inset of panel (b). Finally, panel (b) evidences that not only P is much higher in the energy-seeking mode, but also that the increase in Ω ga leads to a visible increase in the absorption power.

D. Self-organized thermal gradients and active thermalization
Temperature gradients can be crucial to nonequilibrium self-organization. For instance, they can play a significant role in the emergence of complex chemistry [9,[18][19][20][21][22], with possible implications to the problem of the origin of life on earth [20][21][22]. Temperature gradients may also result in broken detailed balance, thus competing with the effects of other nonequilibrium sources (the driving light, in our case) [14]. This competition makes it more challenging to achieve homeostatic resilience, as compared to thermal-equilibrium fluctuations. This alone could be reason enough for us to report on the effects of thermal gradients in this paper. To our surprise, Ee. The spontaneous emission rate of each transition is κ (dashed arrows) and the resonant drive induces the coupling between states |a and |e with strength Ωae (full arrow). P stands for the power consumption (units of Eeaκ, where Eea = Ee − Ea). At times t < 0, the diamond system is at thermal equilibrium. At t ≥ 0 the Λ-type transition undergoes stimulation (Ωae > 0), making the energy-avoiding behaviour to emerge, so that the stationary power absorption P is low and almost insensitive to the increase from Ωae = 0.4κ (dashed blue) to 0.5κ (solid blue). (Inset) the stationary population of state |a decreases, as if the system was intentionally avoiding absorbing higher powers. (b) Same ♦ system as in (a), but with the resonant drive coupling states |g and |a with strength Ωga (full arrow). At times t < 0, the diamond system is at thermal equilibrium. At t ≥ 0 the Vtype transition undergoes stimulation (Ωga > 0), making the energy-seeking behaviour to emerge, so that the stationary power absorption P is high and quite sensitive to the increase from Ωae = 0.4κ (dashed blue) to 0.5κ (solid blue). (Inset) the stationary population of state |a increases, as if the system was intentionally seeking to absorb higher powers. In all plots, we consider the temperature T = 0.5Eea/kB. however, we have found that the emergent behaviours analyzed above imply self-organized functionalities related to thermal gradients: a four-level system can autonomously act upon the temperature gradients in its own environment and modify it. This reminds us of the apparently purposeful, end-directed actions taken by living organisms when altering their own surroundings, as well as of the measurable thermal gradients within the boundaries of single living cells [23].
Specifically, we find that the energy-seeking mode results in a self-organized thermal gradient in the environment, as explained below. Let us suppose that the ♦ system is coupled to two distinct environments, namely, left (L) and right (R), both initially at thermal equilibrium. Environment L is the one coupled to the |g -|a -|e transitions, whereas R is coupled to the |g -|b -|e transitions. The above assumption is inspired by models of energy transport through small quantum chains and their applications in biophysical scenarios [24][25][26][27]. When the (V -type) transition |g -|a is stimulated, we find that the energy-seeking mode emerges. The system thus preferentially dissipates towards L, that is, |J L | > |J R | (here, J L ≡ J ea + J ag and J R ≡ J eb + J bg ). The key point is that, if the heat capacities of both L and R are finite and comparable, the left temperature will raise faster, implying that T L > T R at some point. Once L is warmer, the system will keep favoring dissipation towards L, therefore tending to increase the disequilibrium even further. We also find that the system preferentially dissipates towards environment L even if the temperature gradient is initially (externally) set against the natural tendency of the self-organized gradient, that is, when T L is initially lower than T R . Figures (4) (a) and (b) characterize this preferential dissipation towards environment L. In (a), we set T L < T R . At times t < 0, the stimulating field is off (so that P = 0), and the heat currents follow standard steady thermal conduction (J R = −J L > 0), thermally breaking detailed balance. At t ≥ 0, the stimulating field is turned on and kept constant, so that high power consumption is achieved. Because this power consumption leads to increased jumps between states |a and |g , the dissipation J ag becomes predominant, hence |J L | > |J R |. In (b), we reverse the thermal gradient, setting T L > T R . At times t < 0, the absence of the light field makes for J L = −J R > 0, as expected. At t ≥ 0, the high power consumption still leads to higher J ag . Therefore the net dissipation is again more favorable towards L, i.e., |J L | > |J R |.
We find the opposite behaviour when the (Λ-type) transition |a -|e is stimulated. In that case, the energyavoiding mode emerges. If the environments have initially equal temperatures, the dissipation towards R will be slightly higher (not shown). Higher dissipation towards R leads to an increase in T R . However, as soon as T R becomes sufficiently higher than T L , dissipation towards environment L turns out becoming higher than that towards R (as further discussed in the following paragraph). That is, the dissipation rates become reversed, |J L | > |J R | (in comparison to the T L = T R case, in which |J L | < |J R |). Due to this higher dissipation towards L, the temperature T L is now led to increase faster. As a consequence, this will eventually lead to T L > T R . When T L > T R , the dissipation rates will readjust themselves again, turning back to |J R | > |J L | (as also discussed in the following paragraph). To sum up, we find that, whenever there is a thermal gradient, dissipation Ea < E b Ee. Spontaneous emission transitions (dashed arrows) are considered with rates κ. The coupling strength with the resonant drive (full arrow) is set to Ωga = 0.5κ (at times t ≥ 0). Power consumption is characterized by P . Heat currents are characterized by JL and JR, such that J L(R) is coming from the left (right) environment, at temperature T L(R) . We plot P (t) (black), JL(t) (purple), and JR(t) (blue), in units of Eeaκ (where Eea = Ee −Ea). We also set TL = 0.2Eea/kB and TR = 0.4Eea/kB. (b) Same as in (a), but with TL = 0.4Eea/kB and TR = 0.2Eea/kB. In both (a) and (b), dissipation is higher towards environment L, i.e., |JL(∞)| > |JR(∞)|, thus enabling a self-organized thermal gradient (towards TL > TR), and consuming high powers.
will always be favored towards the coldest environment, in the energy-avoiding mode. In other words, the system actively seeks to reestablish thermal equilibrium in its own environment. Figures (5) (a) and (b) characterize this preferential dissipation towards the coldest environment, as stated in the previous paragraph. In Fig.(5)(a), we set T L < T R . At times t < 0, the stimulating field is off (so that P = 0), and the heat currents follow standard steady thermal conduction (J R = −J L > 0). At t ≥ 0, the stimulating field is turned on and kept constant, so that low power consumption takes place (see inset). Because this power consumption ends up by decreasing the popula-tion of state |a and increasing the population of state |b , the emission towards environment R increases, so that the (positive) difference between the absorbed and the emitted heat decreases, resulting in |J L | > |J R |. In Fig.(5)(b), we set T L > T R . At times t < 0, the nonequilibrium field is absent, making for J L = −J R > 0, as expected from standard thermal conduction. At t ≥ 0, the energy-avoiding behaviour again leads to low power consumption (see inset). Once more, the increased population transfer from |a to |b increases the emission towards environment R, but now resulting in |J R | > |J L |, since the heat emission towards R was higher than the absorption from R (the colder environment) in the absence of the drive. In summary, the plots show that, in the energy-avoiding mode, dissipation is always favored towards the coldest environment, be it L or R.

IV. DISCUSSION
We have established emergent lifelike behaviours and functionalities in elementary, generic open quantum systems. This represents a step further in the direction of looking for general thermodynamic principles underlying far-from-equilibrium self-organization and the onset of lifelike exceptional behaviours, ultimately allowing us to imitate the transition from inanimate to living matter in diverse systems. The ubiquity of the master equation framework used here suggests that the disclosed effects can, perhaps, be found in a much wider variety of quantum and classical systems undergoing stochastic dynamics, provided that the systems are sufficiently nonlinear and multistable, therefore giving rise to statedependent feedback loops, and that the nonequilibrium environments enable transient broken detailed balance.
We mention the time-dependent nature of broken detailed balance so as to suggest the relation between the dissipative history and the dynamic responses considered here. By dissipative history we mean the total (timeintegrated) amount of energy dissipated by the system during the process departing from thermal equilibrium and reaching a nonequilibrium steady state. We hypothesize that, in the context of trying to find a thermodynamic principle characterizing the nonequilibrium stationary state of an arbitrary open physical system, the total dissipated heat may be a relevant quantity (as opposed to the stationary dissipation rate, as analyzed by Kondepudi et al. [7], for instance). We base this hypothesis on two main facts. The first is that the stationary heat dissipation rate here cannot be a candidate for unifying principle. To see that, we compare the Λ and the V systems. When subject to the same environmental conditions (equal temperatures and equal-intensity resonant drives), the Λ and the V systems are clear examples of opposed behaviours in terms of steady power consumption and dissipation. The second is that the significance of the time-integrated dissipated heat as a promising candidate has already been Active thermalization in the energyavoiding mode. (a) ♦ system with eigenstates |g , |a , |b and |e (black horizontal bars), and eigenenergies Eg Ea < E b Ee. Spontaneous emission transitions (dashed arrows) are considered with rates κ. The coupling strength with the resonant drive (full arrow) is set to Ωae = 0.5κ (at times t ≥ 0). Power consumption is characterized by P . Heat currents are characterized by JL and JR, such that J L(R) is coming from the left (right) environment, at temperature T L(R) . We plot P (t) (black), JL(t) (purple), and JR(t) (blue), in units of Eeaκ (where Eea = Ee − Ea). We also set TL = 0.2Eea/kB and TR = 0.4Eea/kB. (b) Same as in (a), but with TL = 0.4Eea/kB and TR = 0.2Eea/kB. In both (a) and (b) (see insets), the dissipation at long times is higher towards the coldest environment, thus inhibiting thermal gradients, and consuming low powers. discussed in the literature [2,3,17,28], giving rise in particular to the concept of dissipative adaptation [2,3,17]. Mathematically, the dissipative adaptation has been conceived from Crooks' microscopically reversible condition [29], and establishes a fluctuation theorem where the nonequilibrium work absorption and the heat dissipation along dynamical (transient) trajectories of a system driven far from equilibrium provide a general thermodynamic mechanism explaining driven self-organization (from self-assembly, in particular, to biological adaptation, in general). The dissipative adaptation paradigm has been recently extended to the quantum realm [17], where the zero-temperature divergences of the classical theorem have been solved. The models analyzed here are closely related to that quantum dissipative adaptation [17], since the same types of self-organized nonequilibrium quantum states have been obtained (specially with the Λ system studied in both papers), but here we have bridged from zero to finite temperatures, comprising thermal gradients, have considered continuous classical light sources (instead of the single-photon pulses), and have achieved self-organized behaviours and functionalities (going beyond self-organized quantum states). The dissipative histories of the driven Λ systems at zero temperature, as analyzed both here and in Valente et al. [17], have been crucial for understanding their self-organized quantum states (see Methods for more details). Additionally, we highlight that Cook and Endres [28] have also put forward the relevance of the transient over the steady-state dissipation to the stability of chemical nonequilibrium systems described by stochastic dynamical equations.
As a perspective, we would like to further investigate how the dissipative history could be a thermodynamic principle of self-organization across quantum and classical regimes. Two quantum properties have been relevant in our modeling in this paper, namely, the discreteness of states and the quantum coherences underlying the broken detailed balance. We believe that quantum coherences, along with quantum-coherent broken detailed balance, will be limited, or even extinguished, by additional environmental noise. However, incoherent nonequilibrium sources may still enable broken detailed balance in systems described by discrete states, even in the presence of dephasing noises. In those incoherent scenarios (i.e., those where broken detailed balance does not require quantum superpositions), our results would be equivalent to self-organization in classical multistable systems containing a finite number of local minima. For that, it suffices that each local minimum is modeled as a discrete, coarse-grained state. A study in this direction could mean an integration of self-organization across quantum and classical open systems, possibly bridging our results to chemical self-organization and thermophoresis [21,22].

Master Equation
The Hamiltonian of the system is H S = i E i |i i|. In the interaction picture with respect to H S , the dynamics of the system under the Born-Markov approximations is governed by the Lindblad master equation for its density operator ρ, H L describes the interaction Hamiltonian between the system and the nonequilibrium drive. The Lindblad su-peroperator describes the thermal processes, where κ ij are the spontaneous emission rates, σ ij ≡ |i j| are the jump operators, and each thermal reservoir is at temperature T ij , with mean number of excitations given by In the models analyzed in this paper, we consider classical resonant light fields, in the rotating-wave approximation, so that H L = i>j Ω ji (σ ji + σ ij ), where Ω ji are the system-field coupling strengths stimulating the transition between states |j and |i (also known as Rabi frequencies).

Energy exchanges
We assume that the average energy of the system is given by E ≡ Tr (H S ρ) in the interaction picture. Energy exchanges are obtained by using the continuity equation (describing energy conservation), namely ∂ t E = Tr (H S ∂ t ρ), where we have used that H S is timeindependent. By applying Eq.(5), we define the input power as the unitary contribution [30,31], and the dissipation rates (heat currents) as the nonunitary contribution [24,30,31], By construction, they satisfy ∂ t E = P + J, where J ≡ i>j J ij .

Feedback loops
We emphasize here how feedback loops play a key role in our results: the power supply first drives the nonlinear system towards nonequilibrium steady-states which, in turn, strongly affect back how much power is absorbed from the drive. This strong dependence of power absorption on the system's state can be regarded as a kind of state-dependent nonlinear susceptibility. To make this idea more explicit, we recast Eq.(5), more precisely, the populations ρ nn ≡ n|ρ|n , in terms of power consumption and dissipation rates as In turn, the power consumption and the dissipation rates depend on the state of the system, that is, so that P = i>j P ij , and where we have defined Note that Γ ij consists of a time-dependent quantitative measure of broken detailed balance.

Λ system
We consider i = a, b, e, so that E a < E b E e . The coupling strengths are such that Ω ae = Ω and Ω ab = Ω be = 0. The spontaneous emission rates are κ ea = κ eb = κ, and κ ba = 0. Similarly, the temperatures are T ea = T eb = T . At T = 0 and Ω = 0, the steady-state is not uniquely defined, hence introducing a kind of multistability. At T = 0 and Ω > 0, the system self-organizes to the (unique) asymptotic state ρ(∞) = |b b|, implying that P = J ea = J eb = 0 (hence preserving detailed balance in the stationary regime, whereas breaking it in the transient regime).
To show the equivalence between our results here (coherent drive, at T = 0) and the single-photon drive analyzed in Valente et al. [17], we combine Eqs.(9)-(12) and find that where the total incoming work is W ≡ ∞ 0 P (t )dt and the dynamics is calculated with the initial state |a (as in Valente et al. [17]). This means that the asymptotic transition probability from |a to |b is maximized along with the total work W performed by the coherent field on the system. Similarly, the total dissipated heat |Q| (where Q ≡ ∞ 0 J(t )dt ) is maximized conditional to the same asymptotic transition probability, namely, ρ bb (∞)| ρaa(0)=1 = (|Q| + E a )/(2E ea − E ba ) (again equivalent to the results in Valente et al. [17]). This highlights how the dissipative history underlies the dynamic transition undergone by the Λ system at zero temperature, both in this semiclassical and in the fully-quantum [17] regimes.
In the numerical calculations, we have set E eb = 0.99E ea .

V system
We consider i = g, b, a, so that E g E b < E a . The coupling strengths are such that Ω ga = Ω and Ω ba = Ω gb = 0. The decay rates are κ ag = κ bg = κ, and κ ab = 0. Similarly, the temperatures are T ag = T bg = T . The steady-state of the V system is uniquely defined for any set of parameters. At T = 0, the system maintains a nonequilibrium coherence ρ ag = −i2κΩ/(κ 2 +8Ω 2 ) and a nonequilibrium population ρ aa = 4Ω 2 /(κ 2 + 8Ω 2 ), implying that P = E ag 4κΩ 2 /(κ 2 + 8Ω 2 ) = −J ag , and J bg = 0. In the numerical calculations, we have set E bg = 0.99E ag .

♦ system
We consider i = g, a, b, e, so that E g E a < E b E e . The coupling strenghts are such that Ω gb = Ω be = Ω ab = Ω ge = 0. For the energy-seeking mode, Ω ga > 0 and Ω ae = 0. For the energy-avoiding mode, Ω ga = 0 and Ω ae > 0. The decay rates are κ ea = κ ag = κ eb = κ bg = κ, and κ ba = κ eg = 0. In the numerical calculations, we have set E ag = 0.9E ea , E bg = 1.1E ea , and E eb = 0.8E ea .
In the case of a thermal equilibrium environment, we consider T ea = T ag = T eb = T bg = T . At T = 0, the stationary power and the dissipation rates exactly coincide with those in the Λ-system (for the energy-avoiding mode) and the V -system (for the energy-seeking mode).
In the case of thermal gradients, we consider T L ≡ T ea = T ag and T R ≡ T eb = T bg , as well as J L ≡ J ea + J ag and J R ≡ J eb + J bg .

VI. DATA AVAILABILITY
Data sharing not applicable to this article as no datasets were generated or analysed during the current study.

VII. CODE AVAILABILITY
The code used to produce the figures in this article is available from the corresponding author upon request.

IX. COMPETING INTERESTS
The authors declare no competing interests.