text:
# Dynamical Evolution: Spirals and Bars

## 1.  Bar formation and destruction

### 1.1. Rapid evolution of disks

Galaxy disks are far from stationary. They are unstable, with time-scales that can be very short according to the radius. This ranges from Myr in the centers to Gyr in the outer parts, According to the environment, companions, mass accretion, etc.. normal modes can be excited. With the mass distribution of normal nearby galaxies, the most rapidly growing mode is the m=2 (spirals and bars). Trailing waves transfer angular momentum to the outer parts, and concentrate the mass. This mass concentration itself changes the development condition of the modes, and weakens them. Alternatively, the disk is heated by waves and the perturbation fades away, until the next excitation.

### 1.2.  Self-regulation

The evolution of the disk is then controlled by recurrent waves, through regulation and feedback processes. One of these is related to gravitational instabilities, suppressed or favored by heating and gas cooling. When the disk is cold ( Toomre Q parameter $`<`$ 1), and therefore unstable to spiral and bar waves, it develops non-axisymmetric perturbations and gravity torques, that transfer the angular momentum outwards. The disk is progressively heated by the waves until the Q threshold is reached. A disk with only old stars cannot cool and will then remain stable, while a galaxy rich in gas will continue to be unstable by gas cooling. Young stars formed out of the gas reform a cool unstable stellar disk.

The gravity torques that drive the gas inwards play the role of a viscosity; the galaxy disk can then be considered as an accretion disk. An exponential stellar disk can then be built, through star formation, if there is approximate agreement between the two time-scales, viscosity and star-formation (Lin & Pringle 1987). This is the case, since the two processes depend exactly on the same gravitational instabilities.

The radial inflow of gas towards the center produces a mass accumulation, that can destroy the bar. This occurs when about 5% of the mass of the disk has sunk inside the inner Lindblad resonance (Hasan & Norman 1990, Pfenniger & Norman 1990, Hasan et al 1993). In the mean time, such a mass accretion in the center has created a nuclear disk, which may decouple kinematically from the rest of the disk; there can develop either a secondary bar (Friedli & Martinet 1993) or nuclear spirals (Barth et al. 1995, Regan & Mulchaey 1999). Alternatively, the decoupling may not occur, and the central pattern rotates at the same velocity as the extended one (Englmaier & Shlosman 2000). This depends on the equivalent viscosity and star-formation rates (Combes 1994, Sheth et al. 2000). It is possible that the decoupling involves also a warp (Schinnerer et al. 2000a,b).

After a bar has dissolved, the axisymmetric disk is devoid of torques, and the new gas accreted stays in the extended disk. If there is significant accretion with respect to the mass of the disk, the latter becomes unstable and another bar forms. External interactions and mergers can intervene to destroy the disk, but also to replenish it. The main agents to drive matter to the center in galaxy interactions and mergers, are also the bars triggered by the tidal forces.

### 1.3. Is the bar responsible for the fueling?

There have been several observational works revealing a correlation between nuclear activity and bars (Dahari 1984, Simkin et al 1980, Moles et al 1995). But the correlation is weak and depends on the definitions of the samples, the completion and other subtle effects. Near-infrared images have often revealed bars in galaxies previously classified unbarred, certainly due to gas and dust effects. However, Seyfert galaxies observed in NIR do not statistically have more bars nor more interactions than a control sample (McLeod & Rieke 1995, Mulchaey & Regan 1997).

Peletier et al (1999) have recently re-visited the question, and took a lot of care with their active and control samples. Their Seyfert and control samples are different at 2.5 $`\sigma `$, in the sense that Seyferts are more barred. They also measure the bar strength by the observed axial ratio in the images. In Seyferts, the fraction of strong bars is lower than in the control sample (Shlosman et al. 2000). Although a surprising result a priori, this is not unexpected, if bars are believed to be destroyed by central mass concentrations (cf section 1.2.).

Regan & Mulchaey (1999) have studied 12 Seyfert galaxies with HST-NICMOS. Out of the 12, only 3 have nuclear bars but a majority show nuclear spirals. However their criterium for nuclear bars is that there exist leading dust lanes along this nuclear bar. This is not a required characteristic, since these secondary bars in general are not expected to have ILRs themselves.

## 2. Lopsidedness, m=1 perturbations

The $`m=1`$ perturbation is present in most galactic disks (Richter & Sancisi 1994), often superposed to the $`m=2`$ ones. These perturbations can be of different nature and origin, according to their scale (nuclear or extended disk), or whether they involve the gaseous or stellar disks.

Linear analysis, supported by numerical calculations, show that gaseous disks rotating around a central mass, are unstable for low value of the central mass (Heemskerk et al 1992). The instability disappears when the central mass equals the disk mass.

Density waves for $`m=2`$ and higher generally use the corotation amplifier, together with reflections at resonances or boundaries, for a steady mode to grow. For $`m=1`$ in nuclear disk, the pattern speed is so low that the corotation falls outside the disk. But, due to the special character of $`m=1`$ waves, there exists another amplifier: the indirect term (Adams et al. 1989, Shu et al. 1990). This is due to the off-centering of the center of mass. Long lasting oscillations of a massive nucleus have been observed (Miller & Smith 1992, Taga & Iye 1998).

For extended disks, off-centered in an extended dark halo, lopsidedness could survive, if the disk remains in the region of constant density (or constant $`\mathrm{\Omega }`$(r)) of the halo (Levine & Sparke 1998).

In the special case of a stellar nuclear disk around a massive black-hole, it is possible that self-gravity is sufficient to compensate for the differential precession of the nearly keplerian orbits, and that a long-lasting $`m=1`$ mode develops, or is maintained long after an external excitation (Bacon et al. 2000). This could be the explanation of the double nucleus observed for a long time in M31 (Bacon et al. 1994, Kormendy & Bender 1999), and for which an eccentric disk model has been proposed (Tremaine 1995, Statler 1999).

In this $`m=1`$ mode, the maximum density is obtained at the apocenter of the aligned elongated orbits (see fig 1). The pattern speed, equal to the orbital frequency of the barycentre of the stellar disk, is slow (3km/s/pc, fig 2), with respect to the orbital frequency of the stars themselves (250km/s/pc in the middle of the disk). The excitation of the $`m=1`$ perturbation can then last more than 3000 rotation periods (fig 3).

## 3. Counter-Rotating Components

The phenomenon of counter-rotating components is a tracer of galaxy interactions, mass accretion or mergers. It has been first discovered in ellipticals with kinematically decoupled cores (likely to be merger remnants, e.g. Barnes & Hernquist 1992). It has been observed also in many spirals; either two components of stars are counter-rotating, or the gas with respect to the stars, or even two components of gas, in different regions of the galaxies (Galletta 1987, Bertola et al 1992). In the special case of NGC 4550 (Rubin et al. 1992), two almost identical counter-rotating stellar disks are superposed along the line of sight.

These systems pose a number of questions, first from their formation scenario, but also about their stability, their life-time, etc.. Do special waves and instabilities develop in counter-rotating disks, and does this favor central gas accretion?

### 3.1. Stability

First, it appears that the counter-rotation (CR) can bring more stability. Even a small fraction of CR stars has a stabilising influence with respect to bar formation ($`m=2`$), since the disk has then more velocity dispersion (Kalnajs 1977). But a one-arm instability is triggered, for a comparable quantities of CR and normal stars. This comes from the two-stream instability in flat disks, similar to that in CR plasmas (Lovelace et al. 1997). There develop two $`m=1`$ modes in the two components, with energies of opposite signs: the negative-E mode can grow by feeding energy in the positive-E mode.

A quasi-stationary one-arm structure forms, and lasts for about 1 to 5 periods (Comins et al. 1997). The structure is first leading, than trailing, and disappears. The formation of massive CR disks in spirals has been studied by Thakar & Ryden (1996, 1998).

### 3.2. CR with gas

The presence of two streams of gas in the same plane will be very transient: strong shocks will occur, producing heating and rapid dissipation. The gas is then driven quickly to the center. But the two streams of gas could be in inclined planes, or at different radii. This is the case in polar rings, discovered in 0.5 % of all galaxies (Whitmore et al. 1990). After correction for selection effects (non optimal viewing, dimming, etc..) 5% of all S0 would have polar rings.

If there is only one gas stream, the problem is more similar to the two-stream instabilities of stellar disks mentioned above. However the gas is cooling, and is not easy to stabilise against $`m=2`$ components. Both $`m=1`$ and $`m=2`$ may be present simultaneously in these systems. This is the case of the galaxy NGC 3593, composed of two CR stellar disks (Bertola et al 1996): in an extended stellar disk, is embedded a CR nuclear disk, possessing co-rotating gas. The molecular component associated with this nuclear disk reveals both a nuclear ring and a one-arm spiral structure outside of the ring (Garcia-Burillo et al. 2000). N-body simulations have shown that both structures can be explained by the superposition of $`m=1`$ and $`m=2`$ in the gas component, the ring being formed at the ILR of the bar (fig. 4). In the $`m=2`$ pattern, two counter-rotating bars develop (fig 5).

## 4.  Conclusions

Galactic disks are not stationary, they evolve greatly on time-scales of the order of a few dynamical times. In normal disks, there is a preponderance of $`m=2`$ spirals and bars. Spirals are the most transient features, but bars also come and go, and several bar episodes could occur in a Hubble time.

According to the mass concentration in the disk, a secondary bar can decouple and rotate with a higher pattern speed. This is important for the ultimate fueling of gas to the nucleus. Observational evidence of correlation between nuclear activity and bars is still controversial. It might be a question of spatial resolution, or the fact that bars are destroyed by massive central black holes.

Lopsidedness is also wide-spread in galactic disks, but much less studied. The nature and origin of the $`m=1`$ waves are diverse, either particular to a purely gaseous disk, or a nuclear stellar disk. They could be slowly damped modes (excited by companions), or two-stream instabilities in counter-rotating systems. There is also the special case of rapid off-centring instabilities (oscillations of the center of mass), or nearly keplerian eccentric disks near a central massive black hole.

Non-axisymmetric perturbations are transferring angular momentum outwards, and result in rapid evolution of galactic disks. It is likely that continuous disk reformation occurs all over the Hubble time, through mass accretion.

#### Acknowledgments.

I thank the organisers for a splendid conference, and my collaborators for letting me use common work prior to publication.

## References

Adams, F. C., Ruden, S. P., Shu, F. H.: 1989, ApJ 347, 959

Bacon R., Emsellem E., Monnet, G., Nieto, J. L.: 1994, A&A 281, 691

Bacon R., Emsellem E., Combes F. et al.: 2000, A&A submitted (B00)

Barnes J.E., Hernquist L.: 1992, ARAA 30, 705

Barth A., Ho L. C., Filippenko, A. V., Sargent, W. L.: 1995, AJ 110, 1009

Bertola F., Buson, L. M., Zeilinger, W. W.: 1992, ApJ 401, L79

Bertola, F., Cinzano, P., Corsini, E.M., et al.: 1996, ApJ, 458, L67

Combes F.: 1994, in “Mass-transfer induced activity in galaxies” ed. I. Shlosman, Cambridge Univ. Press, p. 170

Comins, N. F., Lovelace, R. V. E., Zeltwanger, T., Shorey, P.: 1997, ApJ 484, L33

Dahari, O.: 1984, AJ 89, 966

Englmaier, P., Shlosman, I.: 2000, ApJ 528, 677

Friedli, D., & Martinet, L. 1993, A&A, 277, 27

Galletta, G.: 1987,ApJ 318, 531

Garcia-Burillo S., Sempere M.J., Combes F., Hunt L.K., Neri R.: 2000, A&A, submitted (G00)

Hasan, H., & Norman, C.A. 1990, ApJ, 361, 69

Hasan H., Pfenniger D., Norman C: 1993, ApJ 409, 91

Heemskerk, M. H. M., Papaloizou, J. C., Savonije, G.: 1992, A&A 260, 161

Kalnajs, A.J.: 1977 ApJ 212, 637

Kormendy J., Bender R.: 1999, ApJ 522, 772

Levine S.E., Sparke L.S.: 1998, ApJ 496, L13

Lin D.N.C., Pringle J.E.: 1987, ApJ 320, L87

Lovelace, R. V. E., Jore, K. P., Haynes, M. P.: 1997, ApJ 475, 83

Mcleod, K. K., Rieke, G. H.: 1995, ApJ 441, 96

Miller R.H., Smith B.F.: 1992, ApJ 393, 508

Moles, M., Marquez, I., Perez, E.: 1995, ApJ 438, 604

Mulchaey J.S., Regan M.W.: 1997, ApJ 482, L135

Peletier, R. F., Knapen, J. H., Shlosman, I., et al.: 1999, ApJS 125, 363

Pfenniger, D., & Norman, C. 1990, ApJ, 363, 391

Regan, M. W., Mulchaey, J. S.: 1999 AJ 117, 2676

Richter O-G., Sancisi R.: 1994, A&A 290, L9

Rubin, V. C., Graham, J. A., Kenney, J. D. P.: 1992, ApJ 394, L9

Shinnerer, E., Eckart, A., Tacconi, L.: 2000a, ApJ 533, 826

Shinnerer, E., Eckart, A., Tacconi, L. et al.: 2000b, ApJ 533, 850

Sheth, K., Regan, M. W., Vogel, S. N., Teuben, P. J.: 2000 ApJ 532, 221

Shlosman, I., Peletier, R. F., Knapen, J. H. : 2000, ApJ 535, L83

Shu, F. H., Tremaine, S., Adams, F. C., Ruden, S. P.: 1990, ApJ 358, 495

Simkin S.M., Su H.J., Schwarz M.P.: 1980, ApJ, 237, 404

Statler T.S.: 1999, ApJ 524, L87

Taga M., Iye M.: 1998, MNRAS 299, 111

Thakar A.R., Ryden B.S.: 1996, ApJ 461, 55

Thakar A.R., Ryden B.S.: 1998, ApJ 506, 93

Tremaine S.: 1995, AJ 110, 628

Whitmore, B. C., Lucas, R. A., McElroy, D. B., et al.: 1990 AJ 100, 1489


format: text

text:
# Local minimal energy landscapes in river networks

## I Introduction

Understanding the development of a landscape in terms of fundamental mechanical principles is a formidable task . In spite of this high complexity, recent theoretical studies have resulted in considerable progresses by considering the issue from a viewpoint analogous to the one taken in conventional critical phenomena  where simple models are exploited to identify universality classes. The main idea of this approach is indeed to focus on few fundamental ingredients which, in the spirit of critical phenomena, are expected to provide a reasonable description of the large-scales properties of the network.

The remarkable properties of river networks have been known for some time and they can be condensed in few phenomenological scaling laws which have been confirmed to hold in the observational data . While these laws do not explain the underlying physical mechanisms, they nevertheless provide guidelines for their search. Hence any newly proposed model for river networks ought to be tested against these laws. In the language of critical phenomena, those scaling laws can be used to derive critical exponents and thus discriminate among different universality classes. Among such laws a fundamental role in the physics of the erosion of a landscape due to the flow of water over it, is played by the slope-area law (see Ref.  for a review).

A few Langevin-like equations have been recently proposed for describing the evolution of the landscape under the effect of erosional processes. In Ref. , reparametrization invariance arguments  were used to derive a dynamical equation which yields the slope-area law as a stationary state. The same equation was obtained by Somafai and Sander  using Landau arguments. Other proposals have also been advanced ,.

The equation proposed in Ref. , was studied both analytically in 1+1 dimensions (one spatial and one temporal) and numerically using a self-consistent (SC) solution in 2+1 dimensions, and it was shown to predict a fairly reasonable stationary state quite different from the starting network acting as initial condition. This type of analysis hinges on the observation that the network representing the river has an evolution mirroring that of the evolution profile. Other methods have also been devised which directly address the topological properties of the network itself, the best known being the optimal channel networks (OCN) . This is a lattice model where a functional describing the dissipated energy is minimized in order to find the optimal configuration and it is based on the idea that, presumably, the erosional process taking place on a landscape is driven by a strive for optimality. A simulating annealing procedure  has been implemented to this aim both at zero  and finite temperature . While the latter is aimed to find the absolute minimum, the former is expected to display local minima only. Using exact bounds and finite-size scaling, Maritan et al () showed that the absolute minimum belongs to the mean-field universality class, that, in turn, means that the corresponding network has a highly symmetric pattern with small rivers draining into bigger rivers in a predictable way (this network is akin to the Peano basin ). However this minimum is not easily reachable in the space of all configurations, and one is then led to suspect that real rivers are better described by configurations related to local minima. We further note that stationary states of the aforementioned dynamical equations, are also expected to be associated with local minima for the same reason.

In view of the above discussion, it is apparent that a more complete description of the stability and the scaling behavior associated with these local minima would be desirable. The present work is an attempt in this direction. Our aims are threefold. Our first goal is the full characterization of the possible universality class associated with local minima. Using the SC numerical procedure envisaged in Ref. , we extend previous numerical results both in size and statistics and we use a finite size procedure for a more accurate estimate of the exponents. We then compare these exponents with those calculated using a zero-temperature simulation in the OCN framework with similar sizes and statistics. Our second goal is the generalization of the dynamical equation to include a tunable parameter which turns out to be related to the $`0\gamma 1`$ parameter considered in OCNs and which is responsible of the aggregation mechanism. In Ref.  it was shown that the absolute minimum is insensitive to a variation of $`\gamma `$ in the interval $`1/2\gamma <1`$ . We find that although the final patterns display marked differences as a function of $`\gamma `$, critical exponents show a much smaller discrepancy, which our results indicate to be marginal (at the edge of error bars), both in the SC and OCN case. Finally we test the stability of our results against a variation in the initial conditions. We find that even in this case a marginal difference appears in the effective critical exponents.

The paper is organized as follows. In Sec.II, we review the dynamical equation whereas in Sec.III scaling laws and critical exponents are briefly recalled. Sec.IV contains the definition of OCN and our results are presented in Sec.V. Finally Sec.VI contains some concluding remarks and future perspectives.

## II The dynamical equation

A rather general dynamical equation consistent with general principles, and capturing the physics of erosional processes occurring in river basin, is ,:

$`{\displaystyle \frac{h(x,t)}{t}}`$ $`=`$ $`u+\nu ^2h(x,t)Q^\alpha (x,t)\left[h(x,t)\right]^2+\eta (x,t)`$ (1)
Here $`h(x,t)`$ and $`Q(x,t)`$ are the height of the landscape and the flow rate at position $`x`$ and time $`t`$, $`\eta (x,t)`$ is a noise term and $`\alpha `$ is a parameter which will be related to $`\gamma `$ (see below). The constant term $`u`$ on the right-hand-side of Eq. (1) mimics the so-called geological uplift which is known to originate by tectonic forces. The second term represents a local diffusion term of strength $`\nu `$, condensing both smoothing (rounding of hilltops) and sedimentation (filling of valleys) processes. The third term is non-local and corresponds to an erosion driven by the waterflow on a surface with $`Q(x,t)`$ representing the flow through site $`x`$. The exponent $`\alpha `$ is a parameter of the model which was assumed to be equal to $`1`$ in Ref. and Ref.. The last noise term is added to account for small-scale stochastic processes such as rainfall fluctuations. In the attempt of extracting the basic features associated with the third term, one can ignore both the diffusion and noise terms (the accuracy of this approximation has been discussed in Ref. ). In this case, a stationary state is obtained when the uplift term balances the flow dependent erosional term, that is when

$`|h(x,t)|`$ $``$ $`Q^{\gamma 1}(x,t)`$ (2)
where, for later convenience, we define $`\gamma =1\alpha /2`$. Most of previous work was performed with $`\alpha =1`$ ($`\gamma =1/2`$) (there are few exceptions but with aims and methodologies different from ours, see e.g. Ref. ). In this case Eq. (2) is known as slope-area law , which is a robust empirical law verified by real field observational data. We shall see later on how $`\alpha `$ is related to the dissipated energy functional. Here we note that if we assume (as we shall do hereafter) uniform rainfall (and no ground water), then $`Q(x,t)a(x,t)`$ where $`a(x,t)`$ is the area of the basin draining into point $`x`$ at time $`t`$. The basic result of this dynamical equation is that the non-linear term is able to account for the correct relationship between local slope and water flowrate.

## III Critical exponents and scaling laws

River networks are a remarkable example of systems governed by scaling laws. Although the concept of power and scaling laws has been known to the hydrologists for half a century, it was not until recently that this concept was put into a well defined framework , in analogy with conventional non-equilibrium critical phenomena where power laws are associated with critical exponents. For the sake of completeness, we shall briefly review here a few of the central laws appearing in river networks which we regards as the most fundamental.

Hack observed that the total area $`a`$ draining into a given point and the upstream length $`l`$ going from that point to the source though the path of maximum water flow, were not independent but related by

$`la^h`$ (3)
where $`h`$ is often referred as Hack’s exponent and it ranges in the interval $`[0.50.6]`$ in real rivers. The distributions of drainage areas and upstream lengths also follow a power law. Within the context of finite-size scaling, these may be written as

$`p(a,L)`$ $`=`$ $`a^\tau f(a/L^\varphi )`$ (4)
for areas and

$`\pi (l,L)`$ $`=`$ $`l^\psi g(l/L^{d_l})`$ (5)
for lengths. In Eq. (4) $`p(a,L)`$ is the distribution of drainage areas on a basin of size $`L`$ and $`f(x)`$ is a finite size function. It defines two critical exponents $`\tau `$ and $`\varphi `$ which are not independent . Similarly $`\pi (l,L)`$ is the distribution of the upstream lengths on a basin of size $`L`$, $`g(x)`$ is a finite size function and there exists a relation between $`\psi `$ and $`d_l`$. Note that $`\varphi `$ and $`d_l`$ define the “typical” area ($`a_{\text{typ}}L^\varphi `$) and length ($`l_{\text{typ}}L^{d_l}`$). Finally we remark that one usually distinguishes between self-affine basins ($`d_l=1`$, $`\varphi <2`$) and self-similar basins ($`d_l>1`$, $`\varphi =2`$)  and in both cases only one exponent out of four is independent. They are related via scaling relations as , ,

$`h={\displaystyle \frac{d_l}{\varphi }}\tau =2h\psi ={\displaystyle \frac{1}{h}}`$ (6)
## IV Zero-temperature OCN: Local minima

Within the context of the OCNs, the “dissipated energy” of a river network at a given time, can be written, in continuum notations, as :

$`E(t)`$ $`=`$ $`{\displaystyle 𝑑x|h(x,t)|Q(x,t)}`$ (7)
The local gradient $`|h(x,t)|`$ and the water flow rate $`Q(x,t)`$ are expected to satisfy Eq. (2). Hence one gets:

$`E(t)`$ $`=`$ $`{\displaystyle 𝑑xQ^\gamma (x,t)}`$ (8)
We note that although the energy expression (8) has been derived in the context of OCN , in principle it could be defined in the evolution equation (1) as well, thus providing a way of monitoring the rate of approach to steady state.

The basic assumption of the OCN is that there is a tendency of any real river basin, to assume a configuration which minimizes Eq. (8). A natural question arising is the characterization of the local and absolute minima of $`E`$. It was shown , that the absolute minimum of $`E`$ is insensitive to a variation of the value of $`\gamma `$ in the interval $`1/2\gamma <1`$  and then $`h=1/2`$, $`\tau =3/2`$ and $`d_l=1`$. The analogous issue in the case of local minima is much less clear. In fact, when $`\gamma =0.5`$, numerical work  suggests that there exists a set of local minima which presumably corresponds to a new universality class that is the relevant one for real rivers. Indeed only exceptional events are able to radically modify river courses and so local minima of dissipated energy trap the system with high probability and therefore dominate the statistics.

It is then a vital issue to discriminate whether or not local minima configurations are indeed related to a well defined and robust universality class independent on $`\gamma `$.

## V Results

In this section, we describe numerical procedures and results in detail for each case. All calculations are carried out on a $`L\times L`$ square lattice with periodic boundary conditions in one direction which we identify as the transverse direction. Multiple outlets are allowed in the outflowing longitudinal direction (West side in the figures) whereas an infinite wall is set up on the opposite side (East side). All averages considered in the statistics are carried out only over the river with the largest flow. It is worth stressing that the choice of considering only the maximum river in a multiple outlets environment corresponds to consider the statistical behavior of rivers that are in competition to drain a given region. On the other hand, a single river within a given region is more appropriate if geological constraints are known to exist.

Typically 8 nearest-neighbors (nn) - 4 associated with the square lattice and 4 associated with the two diagonal directions - are allowed. A somewhat more restricted choice considers only the 4 natural nn associated with the square lattice structure. In view of universality one would expect that the details of the lattice structure should not matter. This second choice has the considerable advantage of being less time consuming for numerical purposes. For this reason, although all following figures display patterns obtained with 8 nn, the results reported in this work are obtained with statistics based on 4 nn. In one example, we have explicitly checked that the outcomes using the two choices are consistent within the statistical errors.

Finally, in all our networks, the drainage area $`a(x,t)`$ is computed, at each time step, in a standard way according to :

$`a(x,t)`$ $`=`$ $`{\displaystyle \underset{y(x)}{}}a(y,t)+1`$ (9)
where $`y(x)`$ denotes all sites $`y`$ which drain into $`x`$ and the last term on the right hand side represents a unit rainfall input on each site at all times.

### A Initial condition

There are many possible initial conditions that can be used in numerical analysis of river networks. A popular one among hydrologists is a deterministic comb-like structure . However this initial condition suffers of various drawbacks  and its use in the multiple outlets case, such the one treated here, appears to be somewhat inconvenient. Hence we consider here two other physically reasonable choices, namely a spanning tree (ST) and a Scheidegger network (SN). Although STs are well known in statistical physics mainly due to their relation with the $`q0`$ limit of the Potts model , only recently their topological properties have been studied in details . A suitable variation of spanning trees has even been proposed as a topological model for river networks ,. We have generated STs with multiple outlets by using an adapted Broder’s algorithm (see Ref.  and references therein). We have considered sizes ranging from $`L=32`$ to $`L=512`$ and statistics based on a number of configurations ranging from $`500`$ to $`100`$ respectively. We have then performed a finite size analysis described in the next subsection to extract the most reliable values of the exponents. Our best result for the exponents are $`\tau =1.378\pm 0.002`$, $`\psi =1.596\pm 0.003`$, $`h=0.633\pm 0.003`$, $`d_l=1.25\pm 0.01`$, $`\varphi =2.00\pm 0.01`$ in excellent agreement with the exact results $`\tau =1.375`$, $`\psi =1.6`$, $`h=0.625`$, $`d_l=1.25`$ . This also provides a good test on the quality of our data analysis. Note that the exact value of $`\varphi `$, albeit not known, can be derived from the knowledge of $`h`$ using Eq. (6). This predicts $`\varphi =2.0`$ in perfect agreement with our computed value. We also note that the obtained results are nearly identical to the one predicted and numerically generated for a single outlet .

A ST is a self-similar network in that it is undirected and isotropic. A quite different choice (directed and anisotropic and hence self-affine) is a SN. This was again proposed as a topological model for river networks (see e.g. ) and the $`\tau `$ exponent has been exactly determined via a mapping to a one-dimensional model for mass aggregation . A similar mapping to a diffusion-reaction model also provides the solution in general dimensionality . Our best results for the critical exponents with the same statistics as above are $`\tau =1.337\pm 0.003`$, $`\psi =1.49\pm 0.02`$, $`h=0.69\pm 0.02`$, $`d_l=0.96\pm 0.01`$, $`\varphi =1.53\pm 0.02`$, which are in excellent agreement with the expected values ($`\tau `$ being exact, the others determined via Eq. (6)), that is $`\tau =4/3`$, $`\psi =1.5`$, $`h=0.75`$, $`d_l=1`$, $`\varphi =1.5`$.

In Figs.1 and 2, we depict typical patterns for a ST and a SN respectively. The presence of multiple outlets magnifies the main difference between ST and SN. A ST is constructed in such a way that total freedom is given to the meandering of streams, the only constraint being that they all terminate on the same line (West side in the Figures). Typically this allows the formation of a main big river of size considerably larger with respect to the others. For SNs, the East-West preferred direction prevents the forming of such big river and many smaller rivers are usually present as shown in Fig.2.

### B The Dynamical Equation: A Self-Consistent solution

As we discussed in Sec. II, we seek the stationary states of the following simplified equation:

$`{\displaystyle \frac{h(x,t)}{t}}`$ $`=`$ $`uQ^\alpha (x,t)\left[h(x,t)\right]^2+\eta (x,t)`$ (10)
where $`\alpha =2(1\gamma )`$. The stationary averaged states of Eq. (10) are expected to conform to Eq. (2). Despite the apparent simplicity of the equation, an explicit numerical solution proves to be rather slow . The reason for this can be traced back to the particular form of the erosion term. According to Eq. (10) only sites with a non-negligible combination $`Q^\alpha (x,t)\left[h(x,t)\right]^2`$ will affect the change of the pattern. As it turns out, this yields a long time transient during which the elevation profile evolves very slowly.

Since we are mainly interested in stationary states, there is a way out from this situation . The main idea is to start with an arbitrary network (e.g. a ST or a SN) and recursively construct the heights starting from the outlets with the aid of Eq. (2). From the derived landscape, a new network (in general different from the original one) can then be obtained by assuming that at each site the outward direction is along the steepest descent path and using Eq.(9). The noise term in Eq. (10) is mimicked by the unity term in Eq. (9). The procedure can then be iterated until self-consistency is finally achieved. The final configuration is, by definition, a stationary state of Eq. (10).

The convergence of the above procedure as a function of the number of iterations is reported in Fig. 3 for various values of $`\gamma `$. Two remarks are in order. First we note that the dissipated energy as obtained from the above SC procedure does not have any physical meaning in the transient state. In other words the definition of “time” for the independent variable appearing in Fig. 3 should be considered only as a short-hand notation for “number of iterations”. Second it is apparent a the value $`\gamma =0.5`$ seems to be the one with the slowest convergence ratio. This is probably due to the particular role played by the value $`\gamma =0.5`$ as it will become clear shortly.

Figs.4-8 depict typical patterns obtained on changing the parameter $`\gamma `$ in the interval $`[0,1]`$. The effect of the value of $`\gamma `$ on the aggregation pattern is evident. As $`\gamma `$ increases from $`0`$ to $`1`$, single big rivers draining the entire basin in a snake-like form are less and less favored. The pattern for $`\gamma =0`$ has a strong memory of the original initial tree. The overall effect of the SC procedure is to disfavor long meandering of the streams thus providing a self-affine character to the final tree. The main river becomes rather straight for $`\gamma =0.25`$ and the whole pattern appears to be more symmetric with respect to the case of $`\gamma =0`$. A noteworthy feature is that this tendency is inverted as $`\gamma 0.5`$ and returns back for $`\gamma =0.75`$. As $`\gamma 1`$ rivers become very directed as one would expect on the ground that this limit corresponds to the Scheidegger network . The fact that $`\gamma =0.5`$ most closely resembles real rivers is a reflection of the natural selection by the erosional processes of this value of $`\gamma `$ in terms of Eq. (2) as it has already been well documented in the literature (see e.g. ).

In our numerical estimates of the exponents, for simplicity, we shall restrict our attention to the half region $`[0,0.5]`$ . We also note that this is the region inaccessible to the analytical scheme of Ref. .

For a more accurate evaluation of the critical exponents $`\tau `$ and $`\psi `$ it proves convenient to introduce the integrated probabilities:

$`P(a,L)`$ $`=`$ $`{\displaystyle _0^a}𝑑a^{}p(a^{},L)=a^{1\tau }F({\displaystyle \frac{a}{L^\varphi }})`$ (11)
and

$`\mathrm{\Pi }(l,L)`$ $`=`$ $`{\displaystyle _0^l}𝑑l^{}\pi (l^{},L)=l^{1\psi }G({\displaystyle \frac{l}{L^{d_l}}})`$ (12)
where $`F(x)`$ and $`G(x)`$ are related to $`f(x)`$ and $`g(x)`$ defined in Eqs. (4) and (5), in an obvious way. An efficient way of computing the exponents is through the so-called “effective” (sometime also referred to as “running”) exponents. In the present case they are defined as:

$`\tau (a)`$ $`=`$ $`1{\displaystyle \frac{\mathrm{log}P(a,L)}{a}}`$ (13)
and

$`\psi (l)`$ $`=`$ $`1{\displaystyle \frac{\mathrm{log}\mathrm{\Pi }(l,L)}{l}}.`$ (14)
One then obtains an effective exponent for each value of the independent variable ($`a`$ or $`l`$).

Fig. 9 shows one typical result on a $`256\times 256`$ lattice. We can divide the obtained values roughly in four regions. The first region ($`1a<10`$) corresponds to the region of no scaling. Small rivers belonging to the second region ($`10a<100`$) have an exponent close to the absolute minimum value $`\tau =1.5`$. This is consistent with the picture of typical rivers (see Fig. 6) where small rivers display a marked straightness similar to the one of the absolute minimum  and it means that they quickly assume configurations consistent with the absolute minimum (independent of the initial conditions). Larger areas are associated with larger rivers which have longer memory of the initial condition (ST in the present case). Hence the corresponding exponent is sensibly smaller (closer to the ST value $`1.38`$). The last region corresponds to the finite size cut-off and must be discarded.

After discarding the first and forth regions the obtained values can then be grouped into local bins and a local average exponent can be associated with each of them. Statistical fluctuations within each box then yield an estimate of error bars. This provides our best estimate of the exponent for each value of $`L`$ and a simple $`1/L`$ extrapolation is then carried out to extract the final values. This is depicted in Fig. 10 and the corresponding best estimates of this method are reported in Table I. An analogous procedure leads to the best estimates for $`\psi `$ as reported again in Table I. Sizes and statistics are identical to those considered for the initial conditions and hence these simulations are rather time consuming.

One can notice a weak dependence on $`\gamma `$ for both $`\tau `$ and $`\psi `$. One way of computing $`\varphi `$ and $`d_l`$ is through the collapse plots of the probabilities $`P(a,L)`$ and $`\mathrm{\Pi }(l,L)`$. However we have found that a satisfactory collapse can achieve only within a limited range of the appropriate variable ($`a/L^\varphi `$ and $`l/L^{d_l}`$). Hence we have opted for an alternative scheme hinging on the calculation of the following ratios:

$`M_a^q(L)`$ $``$ $`{\displaystyle \frac{a^{q+1}_a}{a^q_a}}L^\varphi q=1,2,\mathrm{}`$ (15)
where averages are over the probability densities $`p(a,L)`$ (of the maximum river) and the $`L`$ dependence is straightforward . A similar relationship holds for the lengths:

$`M_l^q(L)`$ $``$ $`{\displaystyle \frac{l^{q+1}_l}{l^q_l}}L^{d_l}q=1,2,\mathrm{}`$ (16)
The results for $`q=1,2,3`$ are reported in Table III for the exponent $`\varphi `$ and in Table IV for the exponent $`d_l`$. Surprisingly the values for $`\varphi `$ are consistently larger that the space-filling value $`\varphi =2`$ which is expected to be the upper bound for this exponent. This is probably due to a deficiency of this procedure during a cross-over from a self-similar regime (the initial ST) to a self-affine pattern (the final configuration). Indeed we shall see later on that this feature is not present when one starts with a self-affine network (e.g. a SN) from the outset. A seemingly large value of $`d_l`$ is appearing probably due to the same reason. Overall, one can notice a rather weak (if any) dependence on $`\gamma `$ and almost no dependence on $`q`$.

Finally Hack’s exponent $`h`$ was computed from its definition Eq. (3) using an effective exponent method to be described below.

Let us assume a power-law dependence for a generic function as given by

$`f(x)`$ $``$ $`x^\theta \theta >0`$ (17)
On integrating $`f(x)`$ (between lower and upper limits, say $`x_0`$ and $`x`$), we find

$`\theta `$ $`=`$ $`{\displaystyle \frac{\left[xf(x)x_0f(x_0)\right]}{_{x_0}^x𝑑zf(z)}}1.`$ (18)
The effective exponent obtained using Eq.(18) in Eq.(3) with $`f(x)a`$, $`xl`$, and $`\theta 1/h`$, can then be analyzed with the same procedure (local average plus $`1/L`$ extrapolation) outlined before. We note that for $`a`$ we have used an averaged value over all areas corresponding to the same length. The final results are reported in Table I and one can see that there is an overall good agreement with scaling laws.

### C OCN

The minimization of the energy functional Eq. (8) goes through an algorithm akin to the one exploited by Sun et. al . It is based on the following steps:

1) An initial configuration (ST or SN) is generated and its dissipated energy is computed according to Eq. (8). By definition, this initial network has no loops.

2) A link is randomly selected and its local outflow is also randomly chosen.

3) This new configuration is tested for loop creation. If a loop has been created, the configuration is rejected and step 2) is repeated.

4) The energy of the new candidate configuration is computed. If it is smaller than the previous one, the new configuration is accepted, otherwise it is rejected.

5) Steps 2)-4) are repeated until the energy does not change within a given tolerance.

The final configuration is regarded as a local minimum. This scheme is patterned after a standard Metropolis algorithm at zero temperature, and averages are over many different configurations (ranging from $`500`$ at $`L=32`$ to $`100`$ at $`L=256`$). Fig. 11 depicts the dissipated energy per unit of length as a function of the convergence “time” (i.e. the number of total iteration of the algorithm). The similarity with patterns obtained from the SC procedure is evident. Here, however, the typical convergence time is much longer than before.

Critical exponents are computed with the same prescription given in the previous subsection. In Fig. 12 we show the resulting $`1/L`$ extrapolation. The final best estimates are reported in Table I.

Once more, a weak dependence on $`\gamma `$ ($`\tau `$ increases as $`\gamma `$ increases) can be noticed. All other exponents are consistent with scaling relations Eq. (6). It is worth stressing that exponents are nearly consistent (at the edge of statistical errors) with those previously obtained from the dynamical equation.

Regarding the exponents $`\varphi `$ and $`d_l`$, they can be found in Tables III and IV respectively. Here too the same feature, discussed in V B in connection with the exponent values of $`\varphi `$ and $`d_l`$, applies.

### D Independence of the initial condition

Our final task is a test of the sensitivity of critical exponents to the initial conditions. To this aim we have changed the initial condition from a ST to a SN for both the dynamical equation and the OCN. The main difference between the two initial conditions is that while the latter is a directed network, the former is not. Table II reports the comparison for the case $`\gamma =0.5`$, and Fig. 13 depicts the network resulting from the SC scheme, for a typical final state. Despite the obvious memory of the initial network (see Fig. 2), the typical distance among big rivers is larger than the initial one. This in turn yields a higher value for $`H`$ and hence non-trivial exponents.

It is remarkable that although these patterns appear to be considerably different from those obtained when starting with a ST, the two sets of critical exponents are nevertheless consistent within the statistical errors.

As a final comment, we find in this case values of $`\varphi `$ and $`d_l`$ in very good agreement with scaling predictions. The results for different ratios $`q`$ in the case $`\gamma =0.5`$ for both the SC and the OCN, along with the comparison with the corresponding values stemming from STs are reported in Tables V and VI respectively.

## VI Conclusions

In this paper we have addressed the issue of the existence and robustness of the universality class associated with landscapes corresponding to local minimal energies.

To this aim we have first extended previous studies for both the SC solution of a Langevin equation and the OCNs variational methods. Higher sizes and statistics have been exploited in both cases and, (to our knowledge) for the first time, the most physical procedure of basing the statistics only on the river with largest flow, has been implemented.

Second, we have monitored the dependence of critical exponents on a parameter $`\gamma `$ associated with the slope-area law in one case (SC) and with the dissipated energy in the other (OCN). Our results give compatible critical exponents between SC and OCN within the error bars, but a weak and similar dependence on $`\gamma `$ appears in both models.

Finally we have tested the stability of the obtained results for both the SC and the OCN with respect to changes in the initial conditions. Although the obtained final patterns display a dependence on initial conditions, critical exponents appear to be insensitive to this dependence.

As a by-product of our investigation, we have found that the SC solution of the dynamical equation is a very powerful method to investigate river networks as it is capable of providing useful informations on the stationary state in a simple and physical way. Another interesting point, from a numerical point of view, is that this procedure typically achieves convergence much faster with respect to OCN scheme.

In view of our results we can now summarize the arguments favoring and disfavoring the appearance of a new universality class associated with local minima. As we mentioned in our discussion of Fig. 10 the typical evolution of a network appears to depend on the considered length scale. Small rivers very quickly settle to their final state whereas much longer time is required to large rivers to forget their initial conditions. This is also reflected by the difficulty in collapsing the distribution probabilities of both $`a`$ and $`l`$ into a single plot for a reasonably extended range of the corresponding variable. It is then possible that although only two universality classes (one associated with the initial condition and the other to the absolute minimum) are present, an intermediate universality class sets on due to both the averaging over different regimes and the difficulty of reaching the absolute minimum.

On the other hand, in support of the existence of a new universality class, we cite the fact that critical exponents are different from those of both ST and SN – on starting with these initial conditions, critical exponents in the final configurations clearly deviate from their initial values. Furthermore all exponents are found to be robust and obey to scaling relations summarized in Sec.III.

Overall we believe that the evidence suggested by our results favors this second possibility which has also been hinted into a different context . In this respect the weak $`\gamma `$ dependence of critical exponents remains unexplained.

It would be interesting to generalize the present calculation in two aspects. For the dynamical equation, it would be instructive to tackle the problem of the explicit solution of Eq. (1). This route has the advantage that Eq. (2) (the key relation between erosional process and network topology) can be then derived rather than assumed as in the self-consistent procedure. Similarly, a parallel calculation could also be implemented in the OCN framework, upon starting with the more general expression given in Eq. (7).

###### Acknowledgements.

It is a pleasure to thank Jayanth Banavar, Amos Maritan and Andrea Rinaldo for many enlighting discussions, innumerable suggestions and a critical reading of the manuscript. Financial support from the Italian MURST (Ministero dell’Universitá e della Ricerca Scientifica e Tecnologica) and INFM (Istituto Nazionale di Fisica della Materia) is gratefully acknowledged.


text:
# 1 Introduction

## 1 Introduction

In statistical physics textbooks, one often encounters systems which consist of macroscopic numbers of identical small parts which do not interact with each other. This is not too unrealistic since there are many physical systems (such as certain spin systems, including nuclear spin systems) which are well approximated by such non-interacting models in some ranges of temperature and time.

The question that we wish to examine here is the following: Does a system which consists of many non-interacting pieces behave as a “healthy” thermodynamic system? It is evident from exercises in statistical physics that the answer is “yes” when only equilibrium properties are concerned. If one focuses on certain non-equilibrium aspects, however, the situation may be different. In fact Sato, Sekimoto, Hondou, and Takagi recently proved that, in a system which consists of many non-interacting small parts, a simple quasi-static process involving contacts with two heat baths can never be reversible in general. It is remarkable that such a system fails to provide us with reversible processes, which are among the building blocks of conventional thermodynamics222 The comparison hypothesis as in  is also violated, if we include a contact with a heat bath (whose temperature is chosen carefully so that no net energy is exchanged) at the end of each “adiabatic process.” . Sato, Sekimoto, Hondou, and Takagi then raise an interesting question whether one can develop a new thermodynamic framework which is capable of describing these unavoidable irreversible processes.

The basic idea of Sato, Sekimoto, Hondou, and Takagi is indeed quite simple. Suppose that a small piece (whose identical copies form the whole system) is a quantum system with $`n`$ energy levels $`\epsilon _1,\mathrm{},\epsilon _n`$. We first assume that the whole system is in equilibrium with a heat bath at inverse temperature $`\beta `$. Then the probability of finding the small system in the $`i`$-th state is $`p_i=e^{\beta \epsilon _i}/z(\beta )`$ where $`z(\beta )`$ is the partition function for the small system. We then gently decouple the system from the bath in such a way that the small system is still described by the same probability $`p_i`$. Then we change a parameter in the model Hamiltonian very slowly, modifying the energy levels to $`\epsilon _1^{},\mathrm{},\epsilon _n^{}`$. The adiabatic theorem tells us that, if the parameter change is sufficiently slow, the probability of finding the small system in the $`i`$-th state is still given by the same333 We number the states so that $`\epsilon _i\epsilon _{i+1}`$ and $`\epsilon _i^{}\epsilon _{i+1}^{}`$.  $`p_i`$. But this $`p_i`$ cannot be represented in the Gibbsian form $`p_i=e^{\beta ^{}\epsilon _i^{}}/z^{}(\beta ^{})`$ for any $`\beta ^{}`$ unless the two energy levels $`\epsilon _1,\mathrm{},\epsilon _n`$ and $`\epsilon _1^{},\mathrm{},\epsilon _n^{}`$ satisfy a special condition. (See (2.14).) Therefore if the system is brought into a contact with a heat bath after the operation, something nontrivial (which turns out to be a jump in the entropy) must take place, no matter how carefully one chooses the temperature of the bath.

Of course such a mechanism applies to any quantum systems, including truly macroscopic ones. But it is wellknown that this kind of deviation from a certain equilibrium distribution is irrelevant in thermodynamic limits. A crucial observation of Sato, Sekimoto, Hondou, and Takagi’s is that the smallness of each non-interacting piece makes this deviation thermodynamically observable.

If one regards non-interacting systems as models of more realistic weakly-interacting systems, then the above irreversibility of Sato, Sekimoto, Hondou, and Takagi should be observed only for operations which take place in a time interval which is not long enough for the system to redistribute its energy between the small parts. If the operation is slow enough and such a redistribution takes place, then the probability distribution for each part at the end of operation becomes essentially Gibbsian, and the process may be reversible.

## 2 Setup and main results

We consider a quantum system which consists of $`N`$ identical “small” quantum systems which do not interact with each other. A single small system has an $`n`$-dimensional Hilbert space, and the whole system has an $`nN`$-dimensional Hilbert space.

The initial Hamiltonian of the system is

$$H=\underset{i=1}{\overset{N}{}}h_i,$$
(2.1)
where, for all $`i`$, $`h_i`$ is identical to a certain Hamiltonian $`h`$ for the small system. Initially the system is in touch with a heat bath, and hence in the Gibbs state

$$\mathrm{}_1^{\mathrm{Gibbs}}=\mathrm{Tr}[(\mathrm{})\rho _1^{\mathrm{Gibbs}}],$$
(2.2)
with

$$\rho _1^{\mathrm{Gibbs}}=\frac{e^{\beta H}}{\mathrm{Tr}[e^{\beta H}]}.$$
(2.3)
We gradually decouple the system from the heat bath, making it thermally isolated. We assume that the state of the system is unchanged from (2.2).

We then perform a “classical” operation on the system. This corresponds to an adiabatic444 Note that the word “adiabatic” means that no “heat” is exchanged between the system and the outside world, and does not mean anything like “slow and gentle.”  process in thermodynamics. We model the operation by a time dependent Hamiltonian

$$H(t)=\underset{i=1}{\overset{N}{}}h_i(t),$$
(2.4)
where all $`h_i(t)`$ are identical. After the operation, the Hamiltonian settles to

$$H^{}=\underset{i=1}{\overset{N}{}}h_i^{},$$
(2.5)
where all $`h_i^{}`$ are identical to $`h^{}`$. We denote by $`U`$ the unitary operator which describes the quantum mechanical time evolution during the whole operation555 The unitary operator $`U`$ is formally defined as follows. Suppose that the operation is done during $`t=t_0`$ and $`t_1`$. Let $`U(t)`$ be the solution of the Schrödinger equation $`iU(t)/t=H(t)U(t)`$ with the initial condition $`U(t_0)=\mathrm{𝟏}`$. Then $`U=U(t_1)`$. In what follows, however, we only use the fact that $`U`$ is unitary. . Since the system is thermally isolated, its state after the operation is simply given by

$$\mathrm{}_2=\mathrm{Tr}[(\mathrm{})\rho _2],$$
(2.6)
with

$$\rho _2=U\rho _1^{\mathrm{Gibbs}}U^1.$$
(2.7)
Note that in general this is not a Gibbs state.

Next we put the system into a contact with another heat bath with the inverse temperature $`\beta ^{}`$. After equilibration, the state of the system becomes another Gibbs state

$$\mathrm{}_3^{\mathrm{Gibbs}}=\mathrm{Tr}[(\mathrm{})\rho _3^{\mathrm{Gibbs}}],$$
(2.8)
with

$$\rho _3^{\mathrm{Gibbs}}=\frac{e^{\beta ^{}H^{}}}{\mathrm{Tr}[e^{\beta ^{}H^{}}]}.$$
(2.9)
Here the inverse temperature is (carefully) chosen so that

$$H^{}_2=H^{}_3^{\mathrm{Gibbs}}$$
(2.10)
holds. This means that there is no net exchange of heat between the system and the second heat bath666 Because the Hamiltonian (2.5) is a sum of mutually commuting pieces, the fluctuation of the energy is $`O(\sqrt{N})`$ in the states $`\mathrm{}_2`$ and $`\mathrm{}_3^{\mathrm{Gibbs}}`$. Since the expectation values of the energy itself is $`O(N)`$, this means that the fluctuation of energy is simply negligible for large enough $`N`$. . This thermal contact thus looks quite innocent. We stress that for a normal thermodynamic system, such a thermal contact without energy exchange produces no observable effects. It is a peculiar nature of systems which consist of non-interacting small pieces that such a contact inevitably leads to irreversibility.

Finally we repeat the same process backwards. We first gently decouple the system from the heat bath (without changing the state of the system), and then make an arbitrary operation which brings back the Hamiltonian from $`H^{}`$ to the initial $`H`$. By denoting $`U^{}`$ the unitary operator for this time evolution, we can write the final (generally non-Gibbsian) state as

$$\mathrm{}_4=\mathrm{Tr}[(\mathrm{})\rho _4],$$
(2.11)
with

$$\rho _4=U^{}\rho _3^{\mathrm{Gibbs}}U^1.$$
(2.12)
From the energy conservation law, one sees that the difference $`H_4H_1^{\mathrm{Gibbs}}`$ is the total energy supplied to the system during the operation. In sufficiently large systems where the energy exchange between the second heat bath and the system is negligible, this difference is the total work done to the system by the outside agent (who performs the operation). Note that, from a thermodynamic point of view, the two operations together form a cycle because the Hamiltonian finally returns to the initial one. One of the basic assumptions in conventional thermodynamics is that the work during such a cycle can be made as small as possible by performing the operation slowly. The limiting operation in which the total work is vanishing is usually called a “reversible” cycle.

But the following theorem due to Sato, Sekimoto, Hondou, and Takagi clearly shows that a reversible cycle is never realized in a system which consists of non-interacting small pieces.

Theorem (Sato, Sekimoto, Hondou, and Takagi) For any unitary operators (which are consistent with the given initial and the final Hamiltonians777 See footnote 5. ) $`U`$ and $`U^{}`$, one has

$$H_4H_1^{\mathrm{Gibbs}}N\mathrm{\Delta }\epsilon (\beta ;h,h^{}),$$
(2.13)
where $`\epsilon (\beta ;h,h^{})`$ depends only on $`\beta `$ and the energy levels $`\epsilon _1,\mathrm{},\epsilon _n`$ and $`\epsilon _1^{},\mathrm{},\epsilon _n^{}`$ of the Hamiltonians $`h`$ and $`h^{}`$, respectively. The quantity $`\epsilon (\beta ;h,h^{})`$ is nonnegative in general, and is strictly positive unless

$$\epsilon _i\epsilon _1=A(\epsilon _i^{}\epsilon _1^{}),$$
(2.14)
for any $`i`$ with a constant $`A`$.

For an explicit form of $`\mathrm{\Delta }\epsilon (\beta ,h,h^{})`$, see (4.40).

## 3 Example

Simple (and probably realistic) examples are paramagnetic spin system. Note, however, that we do not expect the irreversibility due to the “smallness” in a paramagnetic $`S=1/2`$ system, where the small system has only two levels.

The simplest example is the $`S=1`$ model under external magnetic field $`\overline{H}`$ with a crystal field anisotropy $`D`$, whose Hamiltonian is

$$H=\underset{i=1}{\overset{N}{}}\left\{\overline{H}S_i^{(z)}D(S_i^{(z)})^2\right\}.$$
(3.1)
Here $`S_i^{(z)}`$ is the spin operator with three eigenvalues $`+1`$, 0, and $`1`$. Thus each small system has three energy levels $`\overline{H}D`$, 0, $`\overline{H}D`$.

As in the standard adiabatic cooling experiments, we thermally isolate the spin system, and change the external magnetic field $`\overline{H}`$. We find that the condition (2.14) is satisfied when $`D=0`$. Thus we expect that the irreversibility due to the smallness of non-interacting pieces is observable when $`|D/\overline{H}|`$ is sufficiently large.

## 4 Proof

### 4.1 Proof of the main theorem

We prove the main theorem by applying the ideas of Sato, Sekimoto, Hondou, and Takagi to quantum models. We stress that our argument is a straightforward implementation of that of Sato, Sekimoto, Hondou, and Takagi’s. The main points in the proof are the use of Gibbs-like entropy (4.2), (4.3), and the representation (4.13) of the entropy increase as a relative entropy.

Throughout the present section, we only treat a special (and unphysical) case with $`N=1`$, i.e., a single small system. Since the small systems do not interact with each other, results for $`N=1`$ immediately implies the corresponding results for general $`N`$.

For $`i=1,\mathrm{},n`$, let $`|\phi _i`$ and $`|\phi _i^{}`$ be the normalized eigenstates of $`h`$ and $`h^{}`$, respectively. We denote by $`\epsilon _i`$ and $`\epsilon _i^{}`$ the corresponding energy eigenvalues. We number the eigenstates so that $`\epsilon _i\epsilon _{i+1}`$ and $`\epsilon _i^{}\epsilon _{i+1}^{}`$ hold.

Let us define a unitary operator $`U`$ by

$$U|\phi _i=|\phi _i^{},$$
(4.1)
for all $`i=1,\mathrm{},n`$. The adiabatic theorem in quantum mechanics guarantees that the time evolution is described by this unitary operator when the operation is executed infinitely slowly (and if one tunes the phases of the basis states properly). We note that, in usual macroscopic quantum systems, an operation must be unphysically slow for the adiabatic theorem to be relevant. In the present case of systems consisting of non-interacting small systems, however, the slowness required by the theorem may be physically realized since each small system evolve independently.

Once we prove the desired bounds for the special unitary transformation $`U`$, the same bounds for arbitrary $`U`$ (and $`U^{}`$) follow immediately. This is because  for given $`h`$, $`h^{}`$, and $`\beta `$, the energy $`h^{}_2`$ after the first operation for arbitrary $`U`$ (consistent with the Hamiltonians $`h`$ and $`h^{}`$) does not exceed the same quantity obtained from the unitary operation (4.1). (This is nothing but the minimum work principle.) The same is true for the second operation which brings back the Hamiltonian from $`h^{}`$ to $`h`$. Therefore we shall only concentrate on $`U`$ defined by (4.1).

Let us define the Gibbs-like entropy with respect to Hamiltonian $`h`$ as888 Note that the present definition of entropy depends on the choice of Hamiltonian $`h`$, while the Gibbs entropy (and the von Neuman entropy) does not depend on Hamiltonians. In this sense, the present definition is closer (in spirit) to the Boltzmann entropy, which depends on the way of characterizing the system from a macroscopic view point (e.g., by using energy).

$$S_h[\rho ]=k\underset{i=1}{\overset{n}{}}\phi _i|\rho |\phi _i\mathrm{log}\phi _i|\rho \phi _i,$$
(4.2)
and the entropy with respect to $`h^{}`$ as

$$S_h^{}[\rho ]=k\underset{i=1}{\overset{n}{}}\phi _i^{}|\rho |\phi _i^{}\mathrm{log}\phi _i^{}|\rho \phi _i^{}.$$
(4.3)
Since the initial state $`\mathrm{}_1^{\mathrm{Gibbs}}`$ is Gibbsian, we write down the density matrix explicitly as

$$\rho _1^{\mathrm{Gibbs}}=\frac{e^{\beta h}}{\mathrm{Tr}[e^{\beta h}]}=\underset{i=1}{\overset{n}{}}|\phi _ip_i\phi _i|,$$
(4.4)
where

$$p_i=\frac{e^{\beta \epsilon _i}}{z(\beta )},z(\beta )=\underset{i=1}{\overset{n}{}}e^{\beta \epsilon _i}.$$
(4.5)
Thus, from the definition (4.2) of entropy, we find

$$S_h[\rho _1^{\mathrm{Gibbs}}]=k\underset{i=1}{\overset{n}{}}p_i\mathrm{log}p_i=k\beta h_1+k\mathrm{log}z(\beta ).$$
(4.6)
By using (4.4) and (4.1), we find that the density matrix after the first operation is simply written as

$$\rho _2=U\rho _1^{\mathrm{Gibbs}}U^1=\underset{i=1}{\overset{n}{}}|\phi _i^{}p_i\phi _i^{}|,$$
(4.7)
with the same $`p_i`$ as in (4.5). We therefore find from (4.3) that

$$S_h^{}[\rho _2]=k\underset{i=1}{\overset{n}{}}p_i\mathrm{log}p_i=S_h[\rho _1^{\mathrm{Gibbs}}].$$
(4.8)
Our entropy is invariant under slow operation.

We then examine the next Gibbs state $`\mathrm{}_3^{\mathrm{Gibbs}}`$ obtained by letting the system interact with the second heat bath. Its density matrix is simply

$$\rho _3^{\mathrm{Gibbs}}=\frac{e^{\beta ^{}h^{}}}{\mathrm{Tr}[e^{\beta ^{}h^{}}]}=\underset{i=1}{\overset{n}{}}|\phi _i^{}p_i^{}\phi _i^{}|,$$
(4.9)
where

$$p_i^{}=\frac{e^{\beta ^{}\epsilon _i^{}}}{z^{}(\beta ^{})},z^{}(\beta ^{})=\underset{i=1}{\overset{n}{}}e^{\beta ^{}\epsilon _i^{}}.$$
(4.10)
Recall that the inverse temperature $`\beta ^{}`$ is determined from the condition

$$h^{}_2=h^{}_3^{\mathrm{Gibbs}}.$$
(4.11)
From the definition (4.3) of the entropy, we have

$`S_h^{}[\rho _3^{\mathrm{Gibbs}}]`$ $`=`$ $`k{\displaystyle \underset{i=1}{\overset{n}{}}}p_i^{}\mathrm{log}p_i^{}`$ (4.12)
$`=`$ $`k\beta h^{}_3^{\mathrm{Gibbs}}+k\mathrm{log}z^{}(\beta ^{})`$
$`=`$ $`k\beta h^{}_2+k\mathrm{log}z^{}(\beta ^{})`$
$`=`$ $`k\beta {\displaystyle \underset{i=1}{\overset{n}{}}}p_i\{\epsilon _i^{}+k\mathrm{log}z^{}(\beta ^{})\}`$
$`=`$ $`k{\displaystyle \underset{i=1}{\overset{n}{}}}p_i\mathrm{log}p_i^{},`$
where we used (4.11) and (4.7). Combining this with (4.6), we find that the difference of the two entropies can be represented in a remarkable manner as

$$S_h^{}[\rho _3^{\mathrm{Gibbs}}]S_h[\rho _1^{\mathrm{Gibbs}}]=k\underset{i=1}{\overset{n}{}}p_i\mathrm{log}\frac{p_i^{}}{p_i},$$
(4.13)
where the right-hand side is nothing but the relative entropy. From the wellknown property of relative entropies, we find999 Proof: Since $`\mathrm{log}xx1`$, one finds $`(\text{L.H.S.})k_{i=1}^np_i\{(p_i^{}/p_i)1\}=k_{i=1}^n\{p_ip_i^{}\}=0`$. To show that the equality holds only when $`p_i=p_i^{}`$ for all $`i`$, it suffices to note that $`\mathrm{log}x<x1`$ whenever $`x1`$.

$$S_h^{}[\rho _3^{\mathrm{Gibbs}}]S_h[\rho _1^{\mathrm{Gibbs}}]0.$$
(4.14)
Here the equality holds only when $`p_i=p_i^{}`$ for all $`i`$. By examining the explicit formulas (4.5) and (4.10) of $`p_i`$ and $`p_i^{}`$, respectively, one finds that this is possible only when the energy levels satisfy the condition (2.14).

As for the second operation (which brings back the Hamiltonian from $`h^{}`$ to $`h`$), we assume that the time evolution is described by the unitary operator $`U^1`$. We repeat the same argument to show

$$S_h[\rho _4]=S_h^{}[\rho _3^{\mathrm{Gibbs}}].$$
(4.15)
As a theoretical reference, we define yet another Gibbs state $`\mathrm{}_5^{\mathrm{Gibbs}}`$ with the density matrix

$$\rho _5^{\mathrm{Gibbs}}=\frac{e^{\beta ^{\prime \prime }h}}{\mathrm{Tr}[e^{\beta ^{\prime \prime }h}]},$$
(4.16)
where $`\beta ^{\prime \prime }`$ is determined by the condition

$$h_5^{\mathrm{Gibbs}}=h_4.$$
(4.17)
Then exactly the same argument as before shows

$$S_h[\rho _5^{\mathrm{Gibbs}}]S_h^{}[\rho _3^{\mathrm{Gibbs}}]0.$$
(4.18)
From (4.14) and (4.18), we finally get

$$S_h[\rho _5^{\mathrm{Gibbs}}]S_h[\rho _1^{\mathrm{Gibbs}}]0,$$
(4.19)
where the equality holds only when the condition (2.14) is satisfied. Since $`\mathrm{}_1^{\mathrm{Gibbs}}`$ and $`\mathrm{}_5^{\mathrm{Gibbs}}`$ are the Gibbs states for the same Hamiltonian $`h`$ with the inverse temperatures $`\beta `$ and $`\beta ^{\prime \prime }`$, respectively, (4.19) implies that

$$\beta ^{\prime \prime }\beta ,$$
(4.20)
and hence

$$h_4=h_5^{\mathrm{Gibbs}}h_1^{\mathrm{Gibbs}}.$$
(4.21)
Again the equality holds only when the condition (2.14) is satisfied101010 It can be shown easily that the equality indeed holds when we have the condition (2.14) and the time evolution is described by the special $`U`$ and $`U^1`$ defined by (4.1). For a general time evolution, the equality may be violated for any Hamiltonians, but this is of course the usual irreversibility. . This proves the desired theorem since, for the present choice of unitary operator $`U`$, the quantity $`h_4h_1^{\mathrm{Gibbs}}`$ clearly depends only on $`\beta `$, $`\epsilon _1,\mathrm{},\epsilon _n`$, and $`\epsilon _1^{},\mathrm{},\epsilon _n^{}`$, and has been shown to be strictly positive unless (2.14) holds111111 We note that the above proof can be directly applied to classical systems. To do this, we fix a constant $`\delta >0`$, and (as in ) decompose the whole phase space into slices with the same volume $`\delta `$. The slices are determined so that the points in the $`i`$-th slice have initial energies between $`\epsilon _i`$ and $`\epsilon _{i+1}`$. We prepare similar decomposition for the final energy as well. Then we let $`p_i`$ be the probability that the initial state is found in the $`i`$-th slice. If we perform the operation slowly, one finds that, in the final state, the probability that the state is in the $`i`$-slice (in the decomposition according to the final energy) is still given by $`p_i`$. (We only use the adiabatic theorem.) Then we can simply repeat the above proof. By taking the limit $`\delta 0`$ finally, the desired proof for classical systems follow..

### 4.2 Explicit lower bound of the work

We stress that the argument of Sato, Sekimoto, Hondou, and Takagi presented in the previous subsection indeed proves the strict inequality

$$h_4h_1^{\mathrm{Gibbs}}>0,$$
(4.22)
for any $`h`$ and $`h^{}`$ which do not satisfy the condition (2.14). But it might be desirable to have an explicit positive lower bound for the required work $`h_4h_1^{\mathrm{Gibbs}}`$. Let us construct a simple lower bound here.

We start from the representation (4.13) of the entropy difference, and rewrite it as

$`S_h^{}[\rho _3^{\mathrm{Gibbs}}]S_h[\rho _1^{\mathrm{Gibbs}}]`$ $`=`$ $`k{\displaystyle \underset{i=1}{\overset{n}{}}}p_i\mathrm{log}{\displaystyle \frac{p_i^{}}{p_i}}`$ (4.23)
$`=`$ $`k{\displaystyle \underset{i=1}{\overset{n}{}}}p_i\left\{\left(1{\displaystyle \frac{p_i^{}}{p_i}}\right)+{\displaystyle \frac{p_i^{}}{p_i}}1\mathrm{log}{\displaystyle \frac{p_i^{}}{p_i}}\right\}`$
$`=`$ $`k{\displaystyle \underset{i=1}{\overset{n}{}}}p_ig\left({\displaystyle \frac{p_i^{}}{p_i}}\right),`$
where

$$g(x)=x1\mathrm{log}x.$$
(4.24)
Note that $`g(x)0`$ for any $`x>0`$ and $`g(x)>0`$ if $`x1`$.

Let us first treat the simplest nontrivial case with $`n=3`$, i.e., a three level system. Let

$$d_1=\epsilon _3\epsilon _1,d_2=\epsilon _3\epsilon _2,d_1^{}=\epsilon _3^{}\epsilon _1^{},d_2^{}=\epsilon _3^{}\epsilon _2^{}.$$
(4.25)
Then for $`p_i`$ and $`p_i^{}`$ as in (4.5) and (4.10), we have

$$\frac{p_3}{p_3^{}}\frac{p_1^{}}{p_1}=\mathrm{exp}[\beta d_1+\beta ^{}d_1^{}],\frac{p_3}{p_3^{}}\frac{p_2^{}}{p_2}=\mathrm{exp}[\beta d_2+\beta ^{}d_2^{}].$$
(4.26)
Introducing vectors $`𝐝=(d_1,d_2)`$ and $`𝐝^{}=(d_1^{},d_2^{})`$, we observe that

$$\underset{\beta ^{}}{\mathrm{min}}\left|\beta 𝐝\beta ^{}𝐝^{}\right|=\beta D(𝐝,𝐝^{}),$$
(4.27)
where

$$D(𝐝,𝐝^{})=|𝐝|\sqrt{1\frac{(𝐝𝐝^{})^2}{|𝐝|^2|𝐝^{}|^2}}.$$
(4.28)
Note that $`D(𝐝,𝐝^{})`$ is vanishing if $`𝐝`$ and $`𝐝^{}`$ are proportional, and is strictly positive otherwise. From (4.27) we see that

$$|\beta d_i\beta ^{}d_i^{}|\frac{\beta D(𝐝,𝐝^{})}{\sqrt{2}}$$
(4.29)
for at least one of $`i=1`$ or 2. From this bound we can show that, at least for one of $`i=1,2,3`$, we have

$$\left|\mathrm{log}\frac{p_i^{}}{p_i}\right|\frac{\beta }{2\sqrt{2}}D(𝐝,𝐝^{}).$$
(4.30)
To see this, suppose the contrary, i.e., $`|\mathrm{log}(p_i^{}/p_i)|<\beta D(𝐝,𝐝^{})/(2\sqrt{2})`$ for all $`i`$. Then for $`i`$ in (4.29)

$`\left|\mathrm{log}{\displaystyle \frac{p_i^{}}{p_i}}\right|`$ $`=`$ $`\left|\mathrm{log}{\displaystyle \frac{p_3}{p_3^{}}}{\displaystyle \frac{p_i^{}}{p_i}}+\mathrm{log}{\displaystyle \frac{p_3^{}}{p_3}}\right|`$ (4.31)
$`=`$ $`\left|\beta d_i+\beta ^{}d_i^{}+\mathrm{log}{\displaystyle \frac{p_3^{}}{p_3}}\right|`$
$``$ $`|\beta d_i\beta ^{}d_i^{}|\left|\mathrm{log}{\displaystyle \frac{p_3^{}}{p_3}}\right|`$
$``$ $`{\displaystyle \frac{\beta }{2\sqrt{2}}}D(𝐝,𝐝^{}),`$
which is a contradiction. We have here used (4.26), (4.29), the triangular inequality, as well as the assumption (4.30).

Let $`i`$ be such that (4.30) is valid. Since $`g(x)`$ is increasing in $`x`$ for $`x>1`$ and decreasing in $`x`$ for $`x<1`$, we find that $`g(p_i^{}/p_i)`$ is bounded from below by the smaller of $`g(\mathrm{exp}[\beta D(𝐝,𝐝^{})/(2\sqrt{2})])`$ and $`g(\mathrm{exp}[\beta D(𝐝,𝐝^{})/(2\sqrt{2})])`$. Since the latter is smaller, we have

$$g\left(\frac{p_i^{}}{p_i}\right)\frac{\beta }{2\sqrt{2}}D(𝐝,𝐝^{})+\mathrm{exp}\left[\frac{\beta }{2\sqrt{2}}D(𝐝,𝐝^{})\right]1>0.$$
(4.32)
By using the representation (4.23), we can bound the total entropy difference as

$$S_h^{}[\rho _3^{\mathrm{Gibbs}}]S_h[\rho _1^{\mathrm{Gibbs}}]kp_ig\left(\frac{p_i^{}}{p_i}\right)kp_3g\left(\frac{p_i^{}}{p_i}\right).$$
(4.33)
By combining this with (4.18) and121212 Although we can construct a lower bound for $`S_h[\rho _5^{\mathrm{Gibbs}}]S_h^{}[\rho _3^{\mathrm{Gibbs}}]`$ in the same manner, it will depend on $`\beta ^{}`$ and won’t be an explicit bound as it is. Therefore we simply use the weak inequality (4.18) for the second operation.  using (4.32), we get an explicit bound for the entropy increase131313If one prefers a bound which include only simple functions, one may use
$$S_h[\rho _5^{\mathrm{Gibbs}}]S_h[\rho _1^{\mathrm{Gibbs}}]\frac{ke^{\beta \epsilon _3}}{z(\beta )}\times \{\begin{array}{cc}\frac{\beta ^2}{24}\{D(𝐝,𝐝^{})\}^2\hfill & \text{if }\beta D(𝐝,𝐝^{})/(2\sqrt{2})(4/3)\hfill \\ \frac{\beta }{2\sqrt{2}}D(𝐝,𝐝^{})1\hfill & \text{if }\beta D(𝐝,𝐝^{})/(2\sqrt{2})(4/3)\text{,}\hfill \end{array}$$
which is obtained from (4.34) through elementary estimate.

$$S_h[\rho _5^{\mathrm{Gibbs}}]S_h[\rho _1^{\mathrm{Gibbs}}]\frac{ke^{\beta \epsilon _3}}{z(\beta )}\left\{\frac{\beta }{2\sqrt{2}}D(𝐝,𝐝^{})+\mathrm{exp}\left[\frac{\beta }{2\sqrt{2}}D(𝐝,𝐝^{})\right]1\right\}>0.$$
(4.34)
We can use the same estimate to treat models with general $`n`$. Let $`M`$ be the smallest integer which does not exceed $`n/3`$. For each $`m=1,\mathrm{},M`$, we choose three integers $`1\mathrm{}_1(m)<\mathrm{}_2(m)<\mathrm{}_3(m)n`$, in such a way that each $`i=1,\mathrm{},n`$ is chosen at most once141414 A simple choice is $`\mathrm{}_i(m)=3(m1)+i`$. . Then by using the previous estimate for the levels $`\epsilon _i`$, $`\epsilon _i^{}`$ with $`i=\mathrm{}_1(m),\mathrm{}_2(m),\mathrm{}_3(m)`$, we get

$$S_h[\rho _5^{\mathrm{Gibbs}}]S_h[\rho _1^{\mathrm{Gibbs}}]\underset{m=1}{\overset{M}{}}\frac{ke^{\beta \epsilon _{\mathrm{}_3(m)}}}{z(\beta )}\left\{\frac{\beta }{2\sqrt{2}}D(𝐝_m,𝐝_m^{})+\mathrm{exp}\left[\frac{\beta }{2\sqrt{2}}D(𝐝_m,𝐝_m^{})\right]1\right\}>0,$$
(4.35)
where

$$𝐝_m=(\epsilon _{\mathrm{}_3(m)}\epsilon _{\mathrm{}_1(m)},\epsilon _{\mathrm{}_3(m)}\epsilon _{\mathrm{}_2(m)}),$$
(4.36)
and

$$𝐝_m^{}=(\epsilon _{\mathrm{}_3(m)}^{}\epsilon _{\mathrm{}_1(m)}^{},\epsilon _{\mathrm{}_3(m)}^{}\epsilon _{\mathrm{}_2(m)}^{}).$$
(4.37)
It is clear that whenever the energy levels $`\epsilon _i`$, $`\epsilon _i^{}`$ do not satisfy the condition (2.14), one can choose $`\mathrm{}_i(m)`$ in such a way that $`D(𝐝_m,𝐝_m^{})>0`$ for at least one $`m`$.

Finally, we note that the entropy and the energy in the Gibbs state for the Hamiltonian $`h`$ at a general inverse temperature $`\stackrel{~}{\beta }`$ satisfies the wellknown relation

$$\frac{}{\stackrel{~}{\beta }}S_h(\rho _{\stackrel{~}{\beta },h}^{\mathrm{Gibbs}})=k\stackrel{~}{\beta }\frac{}{\stackrel{~}{\beta }}h_{\stackrel{~}{\beta },h}^{\mathrm{Gibbs}}.$$
(4.38)
Dividing this relation by $`k\stackrel{~}{\beta }`$, and integrating it from $`\beta ^{\prime \prime }`$ to $`\beta `$, one finds that

$`h_5^{\mathrm{Gibbs}}h_1^{\mathrm{Gibbs}}`$ $`=`$ $`{\displaystyle _{\beta ^{\prime \prime }}^\beta }𝑑\stackrel{~}{\beta }{\displaystyle \frac{1}{k\stackrel{~}{\beta }}}\left\{{\displaystyle \frac{}{\stackrel{~}{\beta }}}S_h(\rho _{\stackrel{~}{\beta },h}^{\mathrm{Gibbs}})\right\}`$ (4.39)
$``$ $`{\displaystyle \frac{1}{k\beta }}{\displaystyle _{\beta ^{\prime \prime }}^\beta }𝑑\stackrel{~}{\beta }\left\{{\displaystyle \frac{}{\stackrel{~}{\beta }}}S_h(\rho _{\stackrel{~}{\beta },h}^{\mathrm{Gibbs}})\right\}`$
$`=`$ $`{\displaystyle \frac{1}{k\beta }}\{S_h[\rho _5^{\mathrm{Gibbs}}]S_h[\rho _1^{\mathrm{Gibbs}}]\}.`$
By substituting (4.35) and recalling (4.17), we get our final estimate for the irreversible work

$$h_4h_1^{\mathrm{Gibbs}}\underset{m=1}{\overset{M}{}}\frac{e^{\beta \epsilon _{\mathrm{}_3(m)}}}{\beta z(\beta )}\left\{\frac{\beta }{2\sqrt{2}}D(𝐝_m,𝐝_m^{})+\mathrm{exp}\left[\frac{\beta }{2\sqrt{2}}D(𝐝_m,𝐝_m^{})\right]1\right\}.$$
(4.40)
We stress that the bound (4.40) is far from optimal. In the crucial estimate (4.27), for example, we are simply looking for $`\beta ^{}`$ which minimizes the distance between $`\beta 𝐝`$ and $`\beta ^{}𝐝^{}`$. But the true $`\beta ^{}`$ is determined from the condition (4.11) about the balance of energy, which condition is not taken into account in our bound.

It is a pleasure to thank Ken Sekimoto for useful discussions and for letting me know of the results in  before making them public.