Transient fractality as a mechanism for emergent irreversibility in chaotic Hamiltonian dynamics

Understanding irreversibility in macrophysics from reversible microphysics has been the holy grail in statistical physics ever since the mid-19th century. Here the central question concerns the arrow of time, which boils down to deriving macroscopic emergent irreversibility from microscopic reversible equations of motion. As suggested by Boltzmann, this irreversibility amounts to improbability (rather than impossibility) of the second-law-violating events. Later studies suggest that this improbability arises from a fractal attractor which is dynamically generated in phase space in reversible dissipative systems. However, the same mechanism seems inapplicable to reversible conservative systems, since a zero-volume fractal attractor is incompatible with the nonzero phase-space volume, which is a constant of motion due to the Liouville theorem. Here we demonstrate that in a Hamiltonian system the fractal scaling emerges transiently over an intermediate length scale. Notably, this transient fractality is unveiled by invoking the Loschmidt demon with an imperfect accuracy. Moreover, we show that irreversibility from the fractality can be evaluated by means of information theory and the fluctuation theorem. The fractality provides a unified understanding of emergent irreversibility over an intermediate time scale regardless of whether the underlying reversible dynamics is dissipative or conservative.

Understanding irreversibility in macrophysics from reversible microphysics has been the holy grail in statistical physics ever since the mid-19th century 1-3 . Here the central question concerns the arrow of time 4 , which boils down to deriving macroscopic emergent irreversibility from microscopic reversible equations of motion. As suggested by Boltzmann 5 , this irreversibility amounts to improbability (rather than impossibility) of the second-law-violating events. Later studies 6,7 suggest that this improbability arises from a fractal attractor which is dynamically generated in phase space in reversible dissipative systems. However, the same mechanism seems inapplicable to reversible conservative systems, since a zero-volume fractal attractor is incompatible with the nonzero phase-space volume, which is a constant of motion due to the Liouville theorem 3 . Here we demonstrate that in a Hamiltonian system the fractal scaling emerges transiently over an intermediate length scale. Notably, this transient fractality is unveiled by invoking the Loschmidt demon 8 with an imperfect accuracy. Moreover, we show that irreversibility from the fractality can be evaluated by means of information theory 9 and the fluctuation theorem 10 . The fractality provides a unified understanding of emergent irreversibility over an intermediate time scale regardless of whether the underlying reversible dynamics is dissipative or conservative.
In 1872, Boltzmann derived the H-theorem, which claims that the entropy is a nondecreasing function of time, and concluded that the state monotonically approaches equilibrium regardless of the initial state 1 . However, his proof is based on the assumption of molecular chaos, where collisions of molecules are assumed to be memoryless. Loschmidt refuted Boltzmann's idea by posing the irreversibility paradox that bears his name 8 : under the laws of reversible mechanics, any phase-space trajectory can be traced back if we start from the final state with the reversed velocity. Therefore, if we have a trajectory with positive entropy production, we should have its time-reversed trajectory with negative entropy production. Hence, monotonicity of entropy is incompatible with microscopic reversible dynamics. In response to Loschmidt's observation, Boltzmann argued that those states that result in positive entropy production appear much more frequently than those that lead to negative entropy production 5 . This Boltzmann's insight was far ahead of the timeit took a century before it was numerically verified. By studying a dissipative system obeying the time-reversible Nosé-Hoover equations of motion 11,12 , Holian et al. 6 observed that states in phase space collapse onto a fractal attractor and argued that one should sample a state precisely on the zero-volume fractal to cause a permanent violation of the second law. This implies that trajectories with negative entropy production cannot be observed, not because of the violation of the equations of motion but because of the vanishing probability to sample the trajectories 6,13 . It is illuminating to reconsider their discussion in view of the steady-state fluctuation theorem 7,14,15 . According to this theorem, the probability of negative entropy production P (−σ) is suppressed compared with that of the sign-reversed counterpart P (σ) by an exponential factor: P (−σ)/P (σ) e −σ . Therefore, the probability of second-law-violating events vanishes in the long run because of the time-extensive nature of entropy production. It is noteworthy that fractal structures again play key roles in the derivation of the steady-state fluctuation theorem 15,16 .
Thus, the Loschmidt paradox in reversible dissipative systems has been settled on the basis of the fractality in phase space. Here, the fractality arises because the equations of motion explicitly violate the Liouville theorem and therefore the phase-space volume contracts to zero. Thus, the very original paradox posed in a Hamiltonian system remains unsolved.
In the present paper, we reveal that a hitherto unnoticed type of fractality arises transiently in a chaotic Hamiltonian system. The Liouville theorem precludes infinitesimally structured fractal with zero volume but allows finite but exponentially small (e.g. sub-atomicscale) fractal structures to emerge. It is this type of transient fractality that yields emergent irreversibility in chaotic Hamiltonian systems.
To probe the fractality, we employ the standard fractal measure of the box-counting dimension d F . Specifically, we observe how the phase-space object of interest, such as the probability distribution function ρ, is vulnerable under small perturbations with typical length scale l. Let d E be the dimension of the embedding space, namely, phase space, and C l [ρ] denote the convolution of ρ with the d E -dimensional isotropic Gaussian with the standard deviation l. We can show that the probability that the state remains in the original support after the Gaussian convolution satisfies the scaling Figure 1 | Time reversal test by an imperfect Loschmidt demon in a Bunimovich billiard. a, Geometry of the Bunimovich stadium. A particle ballistically moves and elastically bounces on the boundary. b, Foliation in phase space. The phase-space structure evolves from a nonequilibrium state (a uniform distribution in a small cube) into a foliated structure, a typical feature of fractals. c, Schematic illustration of the time reversal test. The initial state ρ 0 undergoes time evolution (TE) over time T to reach the final state ρ T . Loschmidt's demon conducts the time reversal (TR) operation with accuracy l and prepares the stateρ l T , which then evolves over time T and ends up with the stateρ l 0 . Finally, we evaluate the Rényi-0 divergence D 0 (ρ 0 ||Tρ l 0 ) = D 0 (ρ T ||C l [ρ T ]), which reflects fractal vulnerability of ρ T against small perturbations. d, Time reversal test by the Loschmidt demon. When the evolution time T is 20, Tρ l 0 is almost identical to the initial state ρ 0 . However, as T increases further, Tρ l 0 deviates from ρ 0 more significantly, which indicates an increase of the Rényi-0 divergence between them. Although ρ T and Tρ l T look identical, they actually differ for large T , which can be measured by the difference between ρ 0 and Tρ l 0 . This points to the fact that the direct evaluation of D 0 (ρ T ||Tρ l T ) is extremely challenging without the time reversal test. See the Supplementary Information for details of numerical simulations. law as ρ(γ)>0 C l [ρ](γ)dγ ∼ (l 0 /l) dE−dF , where l 0 is the minimal length scale of ρ. Interestingly, the left-hand side can be rewritten by an information-theoretic measure called the Rényi-0 divergence 9 , which is defined as D 0 (ρ||C l [ρ]) := − ln ρ(γ)>0 C l [ρ](γ)dγ, and the scaling law is translated as where we define the fractal codimension by d C := d E −d F . Consequently, we can extract the fractal dimension of ρ by observing the behaviors of the Rényi-0 divergence against the variations of the length scale l. Although the scaling law (1) mathematically applies to the limit of l → 0, we utilize it for sufficiently small l( = 0) to probe fractality in a physical sense. See the Supplementary Information for details.
As a chaotic Hamiltonian system, we consider an elastic billiard in the Bunimovich stadium 17 (see Fig. 1a). The phase space of the Bunimovich billiard can be represented as a three-dimensional space (d E = 3) due to the conservation of energy; the coordinate of the phase space consists of two positions (x and y) and the direction of the velocity (θ). We use the Monte-Carlo method to study this system. First of all, we sample a phasespace point γ 0 from the probability distribution function at the initial time ρ 0 , which is prepared as a uniform distribution in a small cube. Then, we let it evolve over time T to obtain γ T , which collectively constitutes the final state ρ T . To evaluate the fractality, we calculate the Rényi-0 divergence between ρ T and the perturbed state C l [ρ T ]. To this aim, we add to γ T a Gaussian noise with the standard deviation l and judge whether the resulting point γ T belongs to the support of ρ T , contributing to the integral in D 0 . Unfortunately, this judgement is impractical because ρ T has a complicated support with an exponentially elongated boundary (see Fig. 1b). The imperfect Loschmidt demon plays a key role in evaluateing the Rényi-0 divergence (see Fig. 1c). The demon invokes the protocol called the time-reversal test 18,19 . For each point γ T in C l [ρ T ], the demon inverts the velocity to obtain the corresponding pointγ T = T γ T , which constitutes an ensembleρ l T = T C l [ρ T ], where T represents the velocity reversal. Then, the demon lets the point evolve over time T to reach the pointγ 0 and the corresponding ensembleρ l 0 . The fraction ofρ l 0 in the original support of ρ 0 contributes to the integral in D 0 . We note that our protocol is similar to the Loschmidt echo 20 , which is used to characterize decoherence in classically chaotic quantum systems. The difference is that the backward Hamiltonian (rather than the state ρ T ) is perturbed in the Loschmidt echo.
When l is small, the difference between ρ T and Tρ l T cannot be discerned (see Fig. 1d). However, after the reversed time evolution, the difference between ρ 0 and Tρ l 0 is manifested, which is quantitatively characterized by the Rényi-0 divergence D 0 . For T = 20, the state Tρ l 0 returns to the initial state almost completely. However, as T further increases, Tρ l 0 deviates from the initial state ρ 0 and expands over the entire phase space. This difference results in a larger value of D 0 , implying increasing vulnerability of the fractal structure of ρ T . Figure 2a plots the Rényi-0 divergence D 0 by changing the length scale l and time T . For intermediate values of l and T , we can observe D 0 is on a plane. This behavior indicates that D 0 satisfies the scaling law of a fractal in equation (1). In fact, by substituting l 0 we −ΛT into equation (1), we obtain To see the l-dependence more clearly, Fig. 2b gives the derivative of D 0 with respect to ln l. For a fixed T , a plateau can be observed in the intermediate regime of the length scale l, which is the evidence of the linear behavior of D 0 with respect to ln l. Moreover, since the value of the plateau corresponds to the fractal codimension d C , we find d C = 1.0. Therefore, we conclude that the structure with the fractal dimension d F = 2.0 is dynamically generated in the intermediate time scale for a fixed l. As a result, the proportionality factor for T amounts to Λ, which coincides with the Kolmogorov-Sinai (KS) entropy h KS of this system according to the Pesin formula 21 . We note that a similar linear growth, i.e., h KS T + const., in an intermediate time was found for a coarse-grained Shannon entropy in chaotic conservative maps 22 . However, no mention about the fractality is made in ref. 22.
In the Supplementary Information, we give a mathematical argument in support of this emergent fractality and a discussion on numerical accuracy. We stress that this transient fractality is compatible with the Liouville theorem. The Liouville theorem states that the phase-space volume is conserved. Meanwhile, the mathematical fractal with infinitely fine structures has zero volume. Therefore, the initial state with nonzero volume cannot evolve into an ideal mathematical fractal. Consequently, fractal structures seem impossible in the Hamiltonian system. However, we here consider the fractal scaling in a fixed length scale l. By fixing l, the phase-space object with nonzero volume can satisfy the fractal scaling law as long as it has a finer structure than l. Moreover, with l fixed, the state generates finer and finer structures over time until it finally looks uniformly nonzero. This is why fractality is transient in the Hamiltonian system with the Liouville theorem. We note that this transient fractal is different from the so-called fat fractal 23 , which has infinitely fine structures and positive volume at the same time.
We consider irreversibility based on a variation of the finite-time fluctuation theorem [24][25][26] . The detailed fluctuation theorem compares the path probability in the original physical dynamics P [Γ] with that in its appropriately chosen dual dynamicsP [Γ] 10,27 . The ratio of these probabilities corresponds to the entropy production:P [Γ]/P [Γ] = e −σ [Γ] . We here take the dual dynamics to be the imperfect time-reversed processρ l T →ρ l 0 (see Fig. 1c) and define an empirical entropy production by σ = − ln(P [Γ]/P [Γ]). Since the dynamics is deterministic in our case, the path probability is determined only by the initial point and therefore we obtain σ = − ln[Tρ l 0 (γ 0 )/ρ 0 (γ 0 )]. Note that the ensemble average σ is not the thermodynamic entropy production but an information-theoretic measure a b c +++++++ + + + + + + + + + + + + + + + + ++ +++ +++++++ +   (2), which emerges after the exponential growth. e, Saturation and finite-size effects. By reducing the initial phase-space volume (from below to top), we can increase the saturated value due to the finite-size effect. See the Supplementary Information for details of numerical simulations.
Notably, the empirical entropy production does not satisfy the conventional Jarzynski-type equality 24 : e −σ = 1. This is because σ is divergent when γ 0 is outside the support of ρ 0 as σ = − ln[Tρ l 0 (γ 0 )/0]. This phenomenon with divergent entropy production should be distinguished from ordinary irreversibility and referred to as absolute irreversibility 28 . In the presence of absolute irreversibility, the Jarzynski-type equality should be modified 28 as where λ is the probability of the singular part ofρ l 0 with respect to ρ 0 in the measure-theoretic sense (see e.g. refs. 29 and 30). In this case, since λ = ρ(γ)=0 Tρ l 0 (γ)dγ, we have 1−λ = e −D0 . By applying the Jensen inequality to equation (3), we obtain (4) We note that the magnitude relation between D 0 and D 1 is well-known in information theory 9 . In Fig. 3a, we numerically calculate the probability distribution function of the empirical entropy production σ, and verify the fluctuation theorem with absolute irreversibility in Fig. 3b.   Figures 3c and 3d compare the empirical entropy production σ with the Rényi-0 divergence D 0 = − ln(1 − λ) in logarithmic and linear scales, respectively. It can be seen that σ is bounded from below by D 0 .
Moreover, we can see that three different time regimes exist. Initially, D 0 exponentially increases as e ΛT , where Λ is the positive Lyapunov exponent (see Fig. 3c). Then, in the following intermediate time regime, D 0 grows linearly in time with the proportionality factor of h KS as shown in Fig. 3d (see equation (2)). Finally, D 0 gets saturated to be the value determined by the ratio of the initial-state volume to the entire phase-space volume. We note that the behaviors in these three time scales are analogous to the ones of the coarse-grained Shannon entropy 22 . We can postpone the occurrence of this finitesize effect to a later time by reducing the initial volume and enlarge the time domain for the transient fractal as shown in Fig. 3e. The linear growth in the transient time scale is reminiscent of the constant entropy production in a steady state of a dissipative system.
As discussed in the Supplementary Information, all of these arguments together with mathematical discussions lead us to conjecture that a generic chaotic Hamiltonian system with the spatial dimension d has transient fractality with d F = d and exhibits a transient linear growth of the empirical entropy with the rate identical to the KS entropy as long as the finite-size effect is negligible. Thus, we suggest that the transient fractality constitutes a universal mechanism of emergent irreversible phenomena in generic reversible dynamics.
We have demonstrated that transient fractality emerges in the Bunimovich billiard even though it obeys the Liouville theorem. An imperfect Loschmidt demon plays a crucial role in evaluating the Rényi-0 divergence characterizing fractality. By regarding the process generated by the demon as the dual process for the fluctuation theorem, we have evaluated the empirical entropy production σ and verified that it satisfies the fluctuation theorem with absolute irreversibility. The average of σ gives how the state ρ T is different from the state Tρ l T = C l [ρ T ] that the demon or an observer perceives by probing ρ T with nonzero accuracy l. Therefore, the state Tρ l 0 that we retrodict to be the initial state from our observation at time T becomes more distinct from the initial state ρ 0 . Thus, the averaged empirical entropy production σ captures this irreversible informational loss caused by the vulnerability of the transient fractal. We expect that this information-theoretic irreversibility is related to the thermodynamic irreversibility. In particular, it is worth investigating the relation between the empirical and thermodynamic entropy productions. An extension of our classical study to a quantum regime should merit further study because the natural length scale determined by the Planck constant should provide a natural cutoff in addition to l and there should be an interplay between them. Moreover, our fractal picture of irreversibility might give a unified framework to describe information loss in equilibriation regardless of whether a system is isolated or open, and classical or quantum.
Photon Science (ALPS) of JSPS. Author Contributions Y.M. designed the numerical simulations and developed the physical aspect of the theory. N.K. provided mathematical supports. Y.M. and N.K. wrote the manuscript of the main text and the Supplementary Information, respectively, with feedback from the others. All work was done under the supervision of M.U.
Supplementary Information for "Transient fractality as a mechanism for emergent irreversibility in chaotic Hamiltonian dynamics"

I. FRACTAL DIMENSION AND THE RÉNYI-0 DIVERGENCE
In our main text, the fractal dimension d F is determined from the scaling law of the Rényi-0 divergence: where d E is the dimension of the embedding space and l 0 is the smallest length scale of the fractal. From a viewpoint of rigorous mathematics, fractals should have an infinitesimal structure (l 0 → 0). However, from a viewpoint of physics, any fractal generated by a finite number of operations should have a nonzero (albeit extremely small) l 0 . We refer to fractals with an infinitesimal structure as mathematical fractals, and those with a nonzero minimal length scale as physical fractals. In this section, we verify equation (S1) for physical fractals on the basis of some natural assumptions. First, we review how the fractal dimensions are mathematically defined. Let S be a subset in the d E -dimensional Euclidean space. For a fixed a > 0, we consider a hypercubic lattice with spacing a, which partitions the space into d E -dimensional hypercubes (or "boxes") of volume a dE . We denote by N (a) the number of boxes that have nonvanishing overlap with S. The box-counting dimension d box is then defined as the limit of Meanwhile, given the probability distribution π supported on S, it is possible to obtain generalized fractal dimensions 1 including the correlation dimension. For that purpose, we cover the fractal S with N (a) boxes, which are labeled with integers from 1 to N (a). Let π i (a) denote the probability with which a sample taken from π falls on the i-th box. The correlation dimension d cor is defined as where C(a) is Although d box = d cor does not generally hold, this equality is satisfied in many practical cases. For instance, we may assume π to be uniform in the sense that the relation holds for all a > 0. In this case, we obtain C(a) = 1/N (a) + o 1/N (a) and the two fractal dimensions coincide with each other. As we discuss later, the assumption of uniformity is consistent with the time reversal test in the main text. Henceforth, we assume that π has a unique fractal dimension d F = d box = d cor .
It is also known that the function C(a) in equation (S3) can be replaced by the integral: Both C(a) and C (a) roughly estimate the probability with two samples γ, γ independently chosen from π fall within the distance of a. Now, we consider a physical fractal that reduces to a mathematical fractal S as long as we consider a length scale larger than l 0 . We define S to be the union of N (l 0 ) boxes covering S. The length scale l 0 describes the smallest structure in the physical fractal, with which S can be regarded to be the same as the mathematical fractal S. Furthermore, we consider a uniform probability distribution over S . We note that the uniformity of the distribution also holds for the final state ρ T in the time reversal test, since the initial state is uniform and the dynamics is volumepreserving. The probability density function ρ can be written as where 1 S is the indicator function of S , which returns 1 for γ ∈ S and 0 for γ / ∈ S . The Rényi-0 divergence D 0 ρ C l [ρ] can be computed as We further use (S7) to obtain The integral in (S11) coincides with C (l) if we replace the measure ρ(γ)dγ by π(dγ). This replacement is justified when l is sufficiently larger than l 0 , because the difference between ρ(γ)dγ and π(dγ) comes into play only in a length scale below l 0 . Therefore, for 1 l l 0 the Rényi-0 divergence is evaluated as Here the approximations ln N (l 0 ) ≈ −d box ln l 0 and ln C(l 0 ) ≈ d cor ln l 0 follow from the definition of the fractal dimensions (S2) and (S3). We obtain the desired relation (S1) by substituting d F = d box = d cor in the last equation.

II. CONJECTURE ABOUT THE FRACTAL DIMENSION IN A HAMILTONIAN SYSTEM WITH A GENERIC DIMENSION
In this section, we discuss the transient fractality in a generic isolated time-reversal-invariant Hamiltonian system. By a chaos-theoretic approach with differential geometry, we claim that a transient fractal of dimension d F = d should emerge, where d is the spatial dimension of the Hamiltonian system.
We denote by M the phase space of the Hamiltonian system. Under the constraint of energy conservation, M has the dimension of d E = 2d − 1 and can therefore be represented by some local coordinates (x 1 , . . . , x 2d−1 ). In addition, we assume that the support S 0 of the initial state ρ 0 be sufficiently small. More specifically, the spread of S 0 must be comparable to w ( 1) in all directions. We exclude the case in which the extent of S 0 in one direction is by far larger than that in the other directions. Cubes and spheres satisfy this condition.
Following the procedure of the time reversal test, we sample a phase-space point γ 0 from the initial state ρ 0 , which evolves into γ T = Φ T (γ 0 ), according to the time evolution Φ T : M → M over time T . We may use different local coordinates to represent γ 0 and γ T as d E -vectors; we describe them as γ 0 = (x 1 , . . . , x 2d−1 ) and γ T = (x 1 , . . . , x 2d−1 ).
As T goes to infinity, the singular values M j = M j (γ 0 , T ) are governed by the Lyapunov exponents Λ 1 ≤ · · · ≤ Λ 2d−1 as (S17) According to Oseledets' theorem 2 , in a volume-preserving dynamical system, the limit exists for almost all γ 0 . The eigenvalues of L = L(γ 0 ) are independent of the initial state γ 0 , the logarithms of which give the Lyapunov exponents. We may choose the coordinate system (x 1 , . . . , x 2d−1 ) such that L becomes diagonal, with which the linear approximation (S16) holds for M j ∼ e Λj T . Although the coordinate system (x 1 , . . . , x 2d−1 ) is determined locally on γ 0 , these local coordinates can be integrated to form one coordinate system 3 . In this coordinate system, the approximation in (S16) is as good as the convergence of the limit (S18).
We recall that the spread of S 0 , the support of the initial state, is as large as w in any direction. To estimate the Rényi-0 divergence, we roughly approximate S 0 by a hyper-rectangular region as S 0 ∼ {(y 1 , . . . , y 2d−1 ) | a j ≤ y j ≤ b j }. (S20) Here the interval [a j , b j ] contains x j and its width b j − a j is comparable to w. With the linear approximation in (S16), the support S T of the evolved state ρ T assumes the form of Therefore, the probability for γ T to lie in S T can be estimated as The behavior of integrals in (S22) depends on the ratio of M j w to l. When M j w l, this integral becomes of the order of unity, since the interval [M j a j , M j b j ] contains a significant part of the Gaussian peak. When M j w l, the integral contains only a small fraction of the Gaussian peak, and therefore remains of the order of M j w/l. From these observations, we obtain Therefore, the conditional probability in (S21) is estimated to be where we denote by α the greatest integer that satisfies the scale separation M α w/l 1. Although this criterion for the scale separation can be ambiguous for finite T , the ambiguity disappears as T goes to infinity. Now we employ the Lyapunov estimates in (S17) to obtain With sufficiently large T , the condition for the scale separation for α is independent of the initial state γ 0 and can be written as Therefore, we may remove the condition on γ 0 from the probability and finally obtain Hence we have a linear relation between D 0 and ln l (and also T ) in the regime of equation (S26). In particular, when T tends to infinity, the middle term of equation (S26) tends to zero and therefore α is the number of negative Lyapunov exponents. When l < w, this number is generally d − 1 in a time-reversal-invariant system with Lyapunov exponents (S19), unless some exponents other than Λ d happen to be zero. The linear relation (S27) collapses for extremely large T , where the Poincaré recurrence occurs. In the argument of this section, the long-time evolution affects the mapping from the phase space to the coordinate system (x 1 , . . . , x 2d−1 ). Although the perturbed state γ T and the unperturbed state γ T look different from each other in the coordinate representation, it is possible that they are close in the phase space, which leads to wrong conclusions about whether γ T resides in the support of ρ T (See Figure S1). This type of errors can emerge when the spread of S T is at least larger than the system size, but numerical results show that equation (S27) holds beyond this time scale. A reasonable estimate for T when equation (S27) ceases to be valid is yet to be found.

III. DETAILS AND PRECISION OF NUMERICAL SIMULATIONS
A. Details of numerical simulations Throughout our numerical simulations, the geometry of the stadium is set to R = 1 and L = 1. The velocity of the billiard is set to be unity. The initial state ρ 0 is prepared to be a uniform distribution in the cube with the side w and the least-valued vertex (x 0 , y 0 , θ 0 ) = (0.5, 0.5, π/4).
In Fig. 1 in the main text, we set w = 0.01 and l = e −15 . In Fig. 2a in the main text, we set w = 0.01 and calculate the Rényi-0 divergence D 0 by the Monte-Carlo method as described in the main text, while increasing the length scale l from e −20 to 1 by the multiplication of e 0.2 and varying T from 0 to 70 by the increment of 2. For each pair of l and T , we sample 10 9 events from a low-discrepancy sequence generated by the additive recurrence. Although the numerical precision decreases as the value of D 0 increases, the relative statistical error can be estimated to be at most two percent. In Fig. 2b, we calculate dD 0 /d ln l by locally applying the least-square fitting to D 0 with respect to ln l.
In Fig. 3a in the main text, we calculate the probability distribution of the empirical entropy production by setting w = 0.1 and l = e −10 . To evaluate the empirical entropy, we should calculate the probability distributionρ l 0 in the support of ρ 0 . To this aim, we divide the support of ρ 0 into 8 3 sub-cubes. In each sub-cube, we generate 2 24 ∼ 2 · 10 7 samples as the initial state and perform the time reversal test. The total sampling number is therefore 2 33 ∼ 9 · 10 9 . By accumulating the events that terminate in each sub-cube, we numerically evaluate the probability distributionρ l 0 . We increase T from 0 to 70 by an increment of 2. As T increases, the statistical error increases becauseρ l 0 diffuses and the number of events that return to each sub-cube decreases. Nevertheless, the relative statistical errors are at most five percent. From the probability distribution of the empirical entropy production in Fig. 3a, the data in Fig. 3b, c, and d are calculated. In Fig. 3e, we set l = e −3 w and decrease w from e −1 to e −5 by the multiplication of e −1 . For each w and T , the Rényi-0 divergence D 0 is evaluated from 10 9 samples.

B. Chaos and numerical precision
In this section, we argue that the double-precision arithmetics is enough to simulate our chaotic system under the present parameters. To numerically examine the effect of a finite precision, we conduct the time reversal test with no noise. Let 10 −p be the absolute precision of the floating number that we use in our simulations. The largest error originates from the error in the stable direction at the final state. When the velocity is reversed, the direction becomes unstable and the error grows exponentially as e ΛT . Therefore, the typical error is expected to be 10 −p e ΛT 10 0.2T −p , where we use Λ = 0.46 for our Bunimovich billiard. We numerically confirm this estimate by varying the precision of floating numbers as shown in Fig. SS2. Hence, the reliability condition for our numerical simulation is w 10 0.2T −p . Therefore, for w = 0.01 and the double precision (p = 16), the reliability is guaranteed for T 70, which covers the range of Fig. 2 in the main text.