Dynamics of energy fluctuations in equilibrating and driven-dissipative systems

We start by considering the approach of the combined system (composed of systems 1 and 2) to its steady state. If no driving is present (scenario 1), then this steady state is thermal equilibrium. We derive an expression for the evolution of the variance σ²₁ = 〈E₁²〉 − 〈E₁〉² during the entire equilibration process. The result is derived in section 3 by solving for the evolution of the first two moments of equation (1), and using the key relation, equation (2).

We find for the equilibrating systems that the variance is given by

$$\begin{equation} \sigma _{1,eq}^{2}\left( \left\langle E_{1}\right\rangle \right) = \sigma _{1_{0}}^{2}\frac{A_{1}^{2}\left( \left\langle E_{1}\right\rangle \right) }{A_{1}^{2}\left( \left\langle E_{1}\right\rangle _{0}\right) }+ 2A_{1}^{2}\left( \left\langle E_{1}\right\rangle \right) \int _{\left\langle E_{1}\right\rangle _{0}}^{\left\langle E_{1}\right\rangle }\frac{1}{A_{1} ^{2}\left( E^{\prime }\right) \left( \beta _{1}-\beta _{2}\right) } {\rm d}E^{\prime } . \end{equation} \tag{ 3 }$$

Here, 〈E₁〉₀ and $\sigma _{1_{0}}^{2}$ are 〈E₁〉 and σ²₁, respectively, at the initial time. Recall that E_total = E₁ + E₂ is held constant in this expression. It is easy to extend these results when E_total varies between experiments. It is interesting to note that when β₂ is set to zero, this expression is identical to that obtained for a single-driven isolated system [7]. This means that within this theory, driving a system is formally equivalent to attaching it to a bath with infinite temperature. It is straightforward to show, that when system 2 is a bath, so that β₂ can be taken to be a constant, the width σ²₁ approaches the equilibrium value: $\left( \partial _{E_{1}}\beta _{1}\right) _{E_{\rm eq}} ^{-1}=k_{B}T^{2}C$, where E_eq is the equilibrium value of 〈E₁〉 and C is the heat capacity (see e.g., [13]). To see this, note that at equilibrium A₁ must vanish, and β₁ = β₂. Therefore, the entire expression for σ²₁ is controlled by the final approach of E to E_eq where $\beta _{1} \simeq \beta _{2}+( \partial _{E_{1}}\beta _{2}) _{E_{\rm eq}}( E-E_{\rm eq})$, $A_{1}( E^{\prime }) =\frac{{\rm d}A_{1} }{{\rm d}E_{1}}\big|_{E_{\rm eq}}( E-E_{\rm eq})$ and the equilibrium expression follows. Note that away from the final equilibration regime A₁(E) need not be linear.

In the case of driven-dissipative systems (when system 2 is large), we obtain for the variance

$$\begin{equation} \sigma _{1,dd}^{2}\left( \left\langle E_{1}\right\rangle \right) =\sigma _{1_{0}}^{2}\frac{A_{1}^{2}\left( \left\langle E_{1}\right\rangle \right) }{A_{1}^{2}\left( E_{1_{0}}\right) }+2A_{1}^{2}\left( \left\langle E_{1}\right\rangle \right) \int _{E_{1_{0}}}^{\left\langle E_{1}\right\rangle }\frac{Z\left( E^{\prime }\right) }{A_{1}^{2}\left( E^{\prime }\right) }{\rm d}E^{\prime }\ , \end{equation} \tag{ 4 }$$

where Z = [1 − β₂A_F/(β₁A₁)]/(β₁ − β₂). Equations (3) and (4) are our main results for the approach to the steady state. They predict the size of fluctuations in E₁ around its average value. They depend only on the rates of energy injection into the system A_F (which is zero for equilibrating systems) and the rate of energy transfer to the bath A₂. In principle, both these quantities can be measured separately: as a consequence of the independence of the drive and interaction between the systems, A_F can be measured by the rate of energy absorption when system 1 is isolated, and A₂ in an equilibration experiment without the drive.

The framework described above can also be used to study fluctuations in the steady state of driven-dissipative systems, specifically fluctuations of E₁ around 〈E₁〉. Typical fluctuations around the steady state are expected to decay exponentially (see section 3), $\left\langle e_{1}\left( t_{1}\right) e_{1}\left( t_{2}\right) \right\rangle =\frac{B_{s}\tau }{2}{\rm e}^{-\left|t_{2} -t_{1}\right|/\tau }$ where e₁ ≡ E₁ − 〈E₁〉. τ is the relaxation time to the steady state. When A_F ≠ 0, namely for a driven-dissipative system, we find that

$$\begin{equation} \tau A_{F}=\frac{\beta _{1}}{\beta _{2}}\left( \beta _{2}-\beta _{1}\right) \sigma _{1}^{2} . \end{equation} \tag{ 5 }$$

This is our main result for the steady state of driven-dissipative systems. A_F is the rate of energy injected into the system from the drive. In the steady state, this energy is then dissipated into the bath. This expression therefore relates three central quantities characterizing the steady state: the size of the energy fluctuations σ²₁, the rate of energy dissipation A_F and τ which is the relaxation time in the steady state.

The above results can be easily generalized to analyzing other settings. For example, consider a setup where a system is attached to two baths, at temperatures β₂ and β₃, see figure 2. The system approaches a steady state. Here one can show that the size of the variance during the approach is given by equation (4) now with Z(E₁) given by

$$\begin{equation} Z\left( E_{1}\right) =\frac{A_{2}}{\beta _{1}-\beta _{2}}+\frac{A_{3}}{\beta _{3}-\beta _{1}}\ , \end{equation} \tag{ 6 }$$

where A₂(E₁)andA₃(E₁) are the average energy transfer rates and β₁(E₁) is the microcanonical temperature of the system at E₁. For the steady state, one obtains

$$\begin{equation} \tau A_{F}=\frac{\left( \beta _{1}-\beta _{2}\right) \left( \beta _{3} -\beta _{1}\right) }{\beta _{3}-\beta _{2}}\sigma _{1}^{2}\ , \end{equation} \tag{ 7 }$$

where A_F is the average rate of energy passing from one bath, through the system, to the other bath. Other more elaborate setups may be constructed similarly.

Zoom In Zoom Out Reset image size

At the steady state the probability distribution P_s(E₁) is independent of time. As ∂_t〈E₁〉 = 〈A〉 = A(〈E₁〉), at the steady state A(〈E₁〉) must vanish. Using the narrowness of the distribution (recall that the width scales at most as N^1/2), we expand A₁ and B₁₁ to lowest order in e₁ ≡ E₁ − E⁰₁:

$$\begin{equation} A_{1}=-\frac{1}{\tau }e_{1}\ ,\ \ \ B_{11}=B_{s}, \end{equation} \tag{ 13 }$$

where B_s and τ are constants. Equivalently, in this regime the FP equation describes the Brownian motion of the energy in a harmonic potential $\dot{e}_{1}=-e_{1}/\tau +\sqrt{B_{s}}\eta$, where the white noise η(t) satisfies 〈η(t)η(t')〉 = δ(t − t'). τ is then interpreted as the relaxation time, as can be seen from the two time correlation functions $\left\langle e_{1}\left( t_{1}\right) e_{1}\left( t_{2}\right) \right\rangle =\frac{B_{s}\tau }{2}{\rm e}^{-\left|t_{2}-t_{1}\right|/\tau }$. The variance of the energy fluctuations is given by σ²₁ = 〈e₁(t₁)²〉 = B_sτ/2.

When A_F = 0, equations (2) and (13) imply that$\ e_{1}=-\frac{\tau B_{s}}{2}\left( \beta _{1}-\beta _{2}\right)$. Then expanding β₁ around β₂ as done above, we find that $\sigma _{1}^{2}=B_{s}\tau /2={-}( \partial _{E_{1}^{0}}\beta _{1}) ^{-1}$which again reproduces the canonical distribution width. The present derivation gives a dynamic interpretation to this formula.

When A_F ≠ 0 namely for a driven-dissipative system we find, using A₁(E⁰₁) = 0 in equation (2), and σ²₁ = B_sτ/2, we obtain equation (5).

We now illustrate our main results on a gas of hard-sphere particles in a box, simulated by an event-driven molecular dynamics simulation [8]. The simulations demonstrate the validity of the theory for fully deterministic dynamics (apart from bath particles, where a bath is present). In addition, the simulations involve relatively small systems, with just few tens of particles, thus demonstrating the applicability of the predictions to small systems. The systems are completely mixed, making this setup more challenging.

The hard-sphere gas is composed of N₁ particles of mass m₁ and N₂ particles of mass m₂, all of equal size, corresponding to systems 1 and 2 respectively, see figure 3(a). Although the entropy of the two systems between collisions indeed factorizes, the collision process involves a strong interaction, which changes the velocities of the particles by a significant amount. A collision calculation shows that if the two masses are very different, then the energy transfer in each collision is small. In this case, energy transfer occurs over many collisions, fulfilling the assumption of timescale separation (see above). In what follows, we take m₁ = 10⁻⁴, m₂ = 1. (Throughout we use arbitrary units.) The large ratio of masses allows us to demonstrate the validity of the results on a single set of system parameters, for a wide range of driving strengths. The results still hold for a smaller mass ratios, albeit in a narrower range of validity. The box is a unit cube with reflecting boundary conditions, and the particles are taken to occupy a volume fraction of 0.05.

We first consider the approach to equilibrium of two systems in contact, figure 1(a), to be compared with the predictions of equation (3). We take N₁ = 30 for the first system and N₂ = 20 for the second system. N₁andN₂ are chosen to be relatively small in order to test the theory on a mesoscopic system. The initial velocities are sampled from a Maxwell–Boltzmann distribution with β₁ = 60 and β₂ = 3, corresponding to average energies per particle of 〈E₁〉/N₁ = 0.025 and 〈E₂〉/N₂ = 0.5. We start all runs from a fixed total energy E_total = 〈E₁〉 + 〈E₂〉, by performing a (small) rescaling of the m₂-particles' velocities. Gathering statistics over many runs, we calculate at each time the average energy 〈E₁〉(t) and the variance σ²₁(t). The function A₁(〈E₁〉) is obtained by plotting A₁(t) = d〈E₁〉/dt as a function of 〈E₁〉(t). Given A₁(〈E₁〉)¹, we use equation (3) to predict σ²₁(〈E₁〉), and find a good fit with the simulation results, see figures 4(a) and (b).

Zoom In Zoom Out Reset image size

We now turn to the driven-dissipative scenario and test equation (5) for the steady state. We run simulations on a system with N₁ = 10 and N₂ = 50 particles. System 1 is driven by applying short impulses to the m₁-particles, changing their velocity by Δv = Fδt/m₁, where F is a constant force and δt is the impulse duration. This impulse is applied at a constant rate. In order to mimic the behavior of a very large system 2, the velocity of the m₂ particles is changed upon reflection from the wall [9], so as to maintain a constant 〈E₂〉. The quantities τ, A_F, β₁ and σ²₁ are computed from the numerics. The energy of system 2 is maintained at 〈E₂〉/N₂ = 1/2, or β_M = 3. Figure 4(c) shows the results obtained for the two sides of equation (5) as a function of β₁/β₂ for different strengths of the drive. Good agreement is found over a wide range of drive strengths and temperature differences. In this range, A_F increases by a factor of 1000, σ²₁ by a factor of 17 and the relaxation time τ decreases by a factor of 3.

In the second example, we consider a one-dimensional XY model. In this model, a two-dimensional vector is associated with each site. The interaction between neighboring sites i and j is given by H_ij = 1 − cos (θ_i − θ_j), where θ_i is the angle of the vector at site i. The total Hamiltonian is H = ∑_{〈i, j〉}H_ij, running over all neighboring pairs; in the case of a linear spin chain with N spins, this becomes H = ∑^{N − 1}_{i = 1}cos z_i, where z_i = θ_{i + 1} − θ_i.

We attach the system to baths at both ends of the system, with inverse temperatures β₂ and β₃, and demonstrate the derivation of equation (7) in this case. A bath is modeled by adding an interaction to the leftmost (and similarly for the rightmost) spin θ₁, which applies random changes to the angle of that spin. The changes in the angle θ₁ are taken in the range −α ⩽ δθ₁ ⩽ α, with α ≪ 1, and at rates

$$\begin{equation*} w\left( \delta \theta _{1}\right) =w_{0}\frac{{\rm e}^{-\beta _{2}\delta e\left( \delta \theta _{1}\right) /2}}{{\rm e}^{-\beta _{2}\delta e\left( \delta \theta _{1}\right) /2}+{\rm e}^{-\beta _{2}\delta e\left( -\delta \theta _{1}\right) /2}}, \end{equation*}$$

where δe(δθ) = −[cos (z₁ + δθ) − cos z₁] is the change in energy due to the change in θ₁ and w₀ is a constant setting the overall timescale. This interaction functions as a bath with temperature β₂. The last spin, θ_N, is attached to a bath at temperature β₃, by applying a similar protocol with β₂ replaced by β₃.

We now calculate the first two moments of the energy transfer from the bath to the system. Note that the steady-state energy and fluctuations will be given by the combined effect of the two baths. We now use the assumptions described above, namely that the energy exchanges with the bath are slow compared with the relaxation time of the system. Here this enters by assuming that α ≪ 1 at a fixed w₀. As a consequence, the system at a given time is approximately in a microcanonical equilibrium. This means that up to 1/N corrections z₁ is distributed with a Boltzmann weight: P(z₁)∝exp (β₁cos z₁), where $\beta _{1}=\partial _{E_{1}}S_{1}$ is the microcanonical temperature of system 1. The average change in the energy is given by

$$\begin{eqnarray*} A_{12}\left( E_{1}\right) & =\int _{-\alpha }^{\alpha }{\rm d}\left( \delta \theta _{1}\right) \int _{0}^{2\pi }{\rm d}z_{1}P\left( z_{1}\right) w\left( \delta \theta _{1}\right) \delta e\left( \delta \theta _{1}\right) \\ & =\frac{\alpha ^{3}}{3}\frac{I_{1}\left( \beta _{1}\right) }{I_{0}\left( \beta _{1}\right) }\frac{\beta _{1}-\beta _{2}}{\beta _{1}}+O\left( \alpha ^{5}\right) \\ B_{12}\left( E_{1}\right) & =\int _{-\alpha }^{\alpha }{\rm d}\left( \delta \theta _{1}\right) \int _{0}^{2\pi }{\rm d}z_{1}P\left( z_{1}\right) w\left( \delta \theta _{1}\right) \left[ \delta e\left( \delta \theta _{1}\right) \right] ^{2}\\ & =\frac{A_{12}}{\beta _{1}-\beta _{2}}+O( \alpha ^{5}), \end{eqnarray*}$$

where I_n(β) is the modified Bessel of the first kind. The dependence on energy enters through the microcanonical energy–temperature relation: $E_{1}=-N\frac{I_{1}\left( \beta _{1}\right) }{I_{0}\left( \beta _{1}\right) }$. A₁₂, B₁₂ satisfy 2A₁₂ = (β₁ − β₂)B₁₂, which is precisely equation (2) for this system (as A_F = 0 for two coupled systems without any drive). Similarly, if the spin θ_N is attached to a bath at temperature β₃, then 2A₁₃ = (β₁ − β₃)B₁₃. The last two relations are valid more generally, as discussed in previous sections, and encode all the information on the system that is required to derive equations (6) and (7). For example, to obtain equation (7), we attach the system to two baths at temperatures β₂, β₃ acting independently at both sides of the chain. Then (using the conventions of figure 2), A₁ = A₁₂ − A₁₃ and B₁ = B₁₂ + B₁₃. At the steady state one has (see section 3) σ²₁ = B₁τ/2, which together with the relations 2A₁₂ = (β₁ − β₂)B₁₂, and 2A₁₃ = (β₁ − β₃)B₁₃ gives equation (7).

In what follows, we first derive equation (3) of the main text, and then proceed to carefully analyze the conditions for its validity. To derive equation (3) we look at two times t, t + Δt at which the driving protocol has returned to its original state, i.e. where H(t) = H(t + Δt), where H is the Hamiltonian of the combined system. As discussed in section 1 'framework and assumptions', we assume that both subsystems are relaxed in their respective energy shells, with energies E₁ and E₂, at times t and t + Δt. We denote the changes in E₁, E₂ during the time interval Δt by ΔE₁, ΔE₂ respectively, and define ΔE_B and ΔE_F via

$$\begin{equation*} \Delta E_{1}=\Delta E_{F}-\Delta E_{B},\quad \Delta E_{2}=\Delta E_{B} . \end{equation*}$$

ΔE_B is the energy transferred from system 1 to system 2, and ΔE_F is the work done on system 1 by the external drive.

Under these conditions, Liouville's theorem or the unitarity of the dynamics, together with micro-reversibility of the dynamics, imply a Crooks relation for the combined isolated system [14–16, 7]²

$$\begin{eqnarray*} &&\fl P_{E_{1},E_{2}}\left( \Delta E_{1},\Delta E_{2}\right) {\rm e}^{S_{1}\left( E_{1}\right) +S_{2}\left( E_{2}\right) }=\tilde{P}_{E_{1}+\Delta E_{1},E_{2}+\Delta E_{2}}\left( -\Delta E_{1},-\Delta E_{2}\right) {\rm e}^{S_{1}\left( E_{1}+\Delta E_{1}\right) +S_{2}\left( E_{2}+\Delta E_{2}\right) }\ , \end{eqnarray*}$$

where $P_{E_{1},E_{2}}\left( \Delta E_{1},\Delta E_{2}\right)$ is the probability of a transition (E₁, E₂) → (E₁ + ΔE₁, E₂ + ΔE₂), and $\tilde{P}$ is defined similarly, only with respect to the reversed protocol, defined by the dynamics generated by the time-reversed Hamiltonian $\tilde{H}( t^{\prime }) =H( t+\Delta t-t^{\prime })$. Here, we have used the assumption of the additivity on the entropy, valid for weak interactions, S(E₁, E₂) = S₁(E₁) + S₂(E₂), see the main text.

We approximate $\tilde{P}_{E_{1}+\Delta E_{1},E_{2}+\Delta E_{2}}\simeq \tilde{P}_{E_{1},E_{2}}$; this introduces errors which are of order 1/N, as shown below. By the assumption of independence of the driving and interaction mechanisms, ΔE_B and ΔE_F are statistically independent, and we can write $P_{E_{1},E_{2}}\left( \Delta E_{1},\Delta E_{2}\right) =P_{E_{1},E_{2}}^{F}\left( \Delta E_{F},0\right) P_{E_{1},E_{2}}^{B}\left( -\Delta E_{B},\Delta E_{B}\right)$. The Crooks relation then reads

$$\begin{eqnarray*} &&\fl P_{E_{1},E_{2}}^{F}\left( \Delta E_{F},0\right) P_{E_{1},E_{2}}^{B}\left( -\Delta E_{B},\Delta E_{B}\right) {\rm e}^{-\beta _{1}\Delta E_{F}-\left( \beta _{2}-\beta _{1}\right) \Delta E_{B}}\\ &&=\tilde{P}_{E_{1},E_{2}}^{F}\left( \Delta E_{F},0\right) \tilde{P} _{E_{1},E_{2}}^{B}\left( -\Delta E_{B},\Delta E_{B}\right) . \end{eqnarray*}$$

Integrating over ΔE_F, ΔE_B gives a Jarzynski relation

$$\begin{equation} \langle {\rm e}^{-\beta _{1}\Delta E_{F}}\rangle \langle {\rm e}^{{-}( \beta _{2}-\beta _{1}) \Delta E_{B}}\rangle =1. \end{equation} \tag{ A.1 }$$

Rearranging the equation, noting that $\tilde{P}^{B}\left( -\Delta E_{B},\Delta E_{B}\right) =P^{B}\left( \Delta E_{B},-\Delta E_{B}\right)$ and integrating gives

$$\begin{equation} \langle {\rm e}^{-\beta _{1}\Delta E_{F}}\rangle =\langle {\rm e}^{{-}( \beta _{2}-\beta _{1}) \Delta E_{B}}\rangle . \end{equation} \tag{ A.2 }$$

Together the two relations yield

$$\begin{equation*} \langle {\rm e}^{-\beta _{1}\Delta E_{F}}\rangle =\langle {\rm e}^{{-}( \beta _{2}-\beta _{1}) \Delta E_{B}}\rangle =1 . \end{equation*}$$

Taking the logarithm of these relations and expanding to second order gives relations (1) and (2) of the main text.

We now discuss the conditions for the validity of the above derivation. In it we have made the following approximations. First, we have assumed that $\tilde{P}_{E_{1}+\Delta E_{1},E_{2}+\Delta E_{2}}\simeq \tilde{P}_{E_{1} ,E_{2}}$. To find the next correction, we take $\tilde{P}_{E_{1}+\Delta E_{1},E_{2}+\Delta E_{2}}\simeq \tilde{P}_{E_{1},E_{2}}-\Delta E_{1} \partial _{E_{1}}\tilde{P}_{E_{1},E_{2}}-\Delta E_{2}\partial _{E_{2}}\tilde{P}_{E_{1},E_{2}}$. Following the above derivation, one obtains equation (1) together with the first correction

$$\begin{equation*} 2\langle \Delta E_{F}\rangle =(\beta _{1}+\partial _{E_{1} }) \langle \Delta E_{F}^{2}\rangle _{c}, \end{equation*}$$

where 〈ΔE²_F〉_c = 〈ΔE²_F〉 − 〈ΔE_F〉² is the second cumulant of ΔE_F. As β₁ is an intensive and E₁ extensive, the term $\partial _{E_{1}}\left\langle \Delta E_{F}^{2} \right\rangle _{c}$ introduces corrections of order 1/N which we neglect. Similarly, equation (2) becomes

$$\begin{equation*} 2\langle \Delta E_{B}\rangle =( \beta _{2}-\beta _{1} +\partial _{E_{2}}-\partial _{E_{1}}) \langle \Delta E_{B} ^{2}\rangle _{c}\ , \end{equation*}$$

and the corrections to equation (2) are of order 1/N for the same reason.

The second approximation involves taking S_i(E_i + ΔE_i) − S_i(E_i) ≃ β_iΔE_i. Again, we incorporate the next correction, $S_{i}\left( E_{i}+\Delta E_{i}\right) -S_{i}\left( E_{i}\right) \simeq \beta _{i}\Delta E_{i}+\frac{1}{2} S_{i}^{\prime \prime }\left( E_{i}\right) \left( \Delta E_{i}\right) ^{2}$. To this order equations (A.1) and (A.2) read

$$\begin{eqnarray*} &&\big\langle {\rm e}^{-\beta _{1}\Delta E_{F}+\frac{1}{2}S_{1}^{\prime \prime }\Delta E_{F}^{2}}\big\rangle \big\langle {\rm e}^{-( \beta _{2}-\beta _{1}) \Delta E_{B}+\frac{1}{2}( S_{1}^{\prime \prime }+S_{2}^{\prime \prime }) \Delta E_{B}^{2}}\big\rangle =1+S_{1}^{\prime \prime }\langle \Delta E_{F}\rangle \langle \Delta E_{B}\rangle \ ,\\ &&\big\langle {\rm e}^{-\beta _{1}\Delta E_{F}+\frac{1}{2}S_{1}^{\prime \prime }\Delta E_{F}^{2}}\big\rangle ( 1+S_{1}^{\prime \prime }\langle \Delta E_{F}\rangle \langle \Delta E_{B}\rangle ) =\langle {\rm e}^{-( \beta _{2}-\beta _{1}) \Delta E_{B} }\rangle , \end{eqnarray*}$$

respectively, which reduce to the above relations if S''₁〈ΔE_F〉〈ΔE_B〉 ≪ 1. Since $S_{i}^{\prime \prime }\left( E_{i}\right) =-\beta _{i}^{2}\left( E_{i}\right) /C_{v_{i}}\left( E_{i}\right)$ (here the Boltzmann constant k_B = 1, making $C_{v_{i}}$ unitless) this condition reads

$$\begin{equation*} \beta _{1}^{2}\langle \Delta E_{F}\rangle \langle \Delta E_{B}\rangle \ll C_{v_{1}} . \end{equation*}$$

Now we can expand $\ln \big\langle {\rm e}^{-\beta _{1}\Delta E_{F}+\frac{1}{2} S_{1}^{\prime \prime }\Delta E_{F}^{2}}\big\rangle =\ln \big\langle {\rm e}^{-( \beta _{2}-\beta _{1}) \Delta E_{B}+\frac{1}{2}( S_{1}^{\prime \prime }+S_{2}^{\prime \prime }) \Delta E_{B}^{2} }\big\rangle =0$ to next order and find

$$\begin{equation*} {-}\beta _{1}\left\langle \Delta E_{F}\right\rangle +\frac{1}{2}S_{1} ^{\prime \prime }\left\langle \Delta E_{F}^{2}\right\rangle +\frac{1}{2} \beta _{1}^{2}\left\langle \Delta E_{F}^{2}\right\rangle _{c}+\frac{1}{3!} \beta _{1}^{3}\left\langle \Delta E_{F}^{3}\right\rangle _{c}=0 . \end{equation*}$$

As 〈ΔE²_F〉 = 〈ΔE_F〉² + 〈ΔE²_F〉_c this becomes

$$\begin{equation*} -\beta _{1}\left\langle \Delta E_{F}\right\rangle +\frac{1}{2}S_{1} ^{\prime \prime }\left\langle \Delta E_{F}\right\rangle ^{2}+\frac{1}{2}\left( \beta _{1}^{2}+S_{1}^{\prime \prime }\right) \left\langle \Delta E_{F} ^{2}\right\rangle _{c}+\frac{1}{3!}\beta _{1}^{3}\left\langle \Delta E_{F} ^{3}\right\rangle _{c}=0 . \end{equation*}$$

The term β²₁ + S₁'' = β²₁[1 + O(1/N)], since $C_{v_{i}}=O\left( N\right)$. Comparing the first and last terms gives the condition $\frac{1}{2}\beta _{1}^{2}/C_{v_{1} }\left\langle \Delta E_{F}\right\rangle ^{2}\ll \beta _{1}\left\langle \Delta E_{F}\right\rangle$ or

$$\begin{equation*} \beta _{1}\left\langle \Delta E_{F}\right\rangle \ll C_{v_{1}} . \end{equation*}$$

In addition, the third cumulant must be small, β³₁〈ΔE_F³〉_c ≪ β₁〈ΔE_F〉 or

$$\begin{equation*} \beta _{1}^{2}\left\langle \Delta E_{F}^{3}\right\rangle _{c}\ll \left\langle \Delta E_{F}\right\rangle . \end{equation*}$$

A similar argument for $\ln \big\langle {\rm e}^{-( \beta _{2}-\beta _{1}) \Delta E_{B}+\frac{1}{2}( S_{1}^{\prime \prime } +S_{2}^{\prime \prime }) \Delta E_{B}^{2}}\big\rangle =0$ yields the requirement

$$\begin{equation*} \left\langle \Delta E_{B}\right\rangle \ll \frac{\left|\beta _{2}-\beta _{1}\right|}{\frac{\beta _{1}^{2}}{C_{v_{1}}}+\frac{\beta _{2}^{2} }{C_{v_{2}}}}=\frac{C_{v_{1}}C_{v_{2}}\left|\beta _{2}-\beta _{1}\right|}{\beta _{1}^{2}C_{v_{2}}+\beta _{2}^{2}C_{v_{1}}} . \end{equation*}$$

For β₂ − β₁ of order β₁ and β₂, this is satisfied if $\beta _{i}\left\langle \Delta E_{B}\right\rangle \ll C_{v_{j}}$ for i, j = 1, 2. Close to β₁ = β₂, both 〈ΔE_B〉 and |β₂ − β₁| go to zero, and this condition may be hard to test. However, 〈ΔE²_B〉_c remains regular (see discussion in section 2.2 in the main text), and using 2〈ΔE_B〉 = (β₂ − β₁)〈ΔE²_B〉_c, we obtain an equivalent relation

$$\begin{equation*} \left\langle \Delta E_{B}^{2}\right\rangle _{c}\ll \frac{C_{v_{1}}C_{v_{2}} }{\beta _{1}^{2}C_{v_{2}}+\beta _{2}^{2}C_{v_{1}}}\ ; \end{equation*}$$

as $\frac{C_{v_{1}}C_{v_{2}}}{\beta _{1}^{2}C_{v_{2}}+\beta _{2}^{2}C_{v_{1}} }>\min \left( \beta _{2}^{-2}C_{v_{1}},\beta _{1}^{-2}C_{v_{2}}\right)$, we obtain the condition

$$\begin{equation*} \left\langle \Delta E_{B}^{2}\right\rangle _{c}\ll \min \left( \beta _{2} ^{-2}C_{v_{1}},\beta _{1}^{-2}C_{v_{2}}\right) \end{equation*}$$

cited in the main text.