A violation of universality in anomalous Fourier’s law

Since the discovery of long-time tails, it has been clear that Fourier’s law in low dimensions is typically anomalous, with a size-dependent heat conductivity, though the nature of the anomaly remains puzzling. The conventional wisdom, supported by renormalization-group arguments and mode-coupling approximations within fluctuating hydrodynamics, is that the anomaly is universal in 1d momentum-conserving systems and belongs in the Lévy/Kardar-Parisi-Zhang universality class. Here we challenge this picture by using a novel scaling method to show unambiguously that universality breaks down in the paradigmatic 1d diatomic hard-point fluid. Hydrodynamic profiles for a broad set of gradients, densities and sizes all collapse onto an universal master curve, showing that (anomalous) Fourier’s law holds even deep into the nonlinear regime. This allows to solve the macroscopic transport problem for this model, a solution which compares flawlessly with data and, interestingly, implies the existence of a bound on the heat current in terms of pressure. These results question the renormalization-group and mode-coupling universality predictions for anomalous Fourier’s law in 1d, offering a new perspective on transport in low dimensions.

with P the fluid's pressure.
(iii) Heat conductivity scaling: Due to the homogeneity of the interaction potential, the heat conductivity of the 1d diatomic hard-point gas exhibits a well-known density temperature separability [2]. Moreover, standard dimensional analysis arguments show that κ ∝ T /m [2], and the known dimensional anomaly for transport implies in turn that κ ∝ L α at leading order. We now raise these arguments to a formal scaling ansatz κ L (ρ, T ) = L α T /mk(ρ) , with k(ρ) a function solely dependent on density. Note that this ansatz discards possible subleading corrections in L.
where G (ρ) ≡ k(ρ)ρ −5/2 and denotes derivative with respect to the argument. This equation, together with the boundary conditions for the density field, ρ(x = 0, L) = ρ 0,L , which can be inferred from the constraints completely define the macroscopic problem in terms of ρ(x). Note that the externally controlled parameters in the problem are the temperatures of the boundary reservoirs, T 0,L , and the global number density η. The pressure and the heat current can be now obtained as P = T 0 ρ 0 and J = P 3/2 [G(ρ L ) − G(ρ 0 )]/(L 1−α √ m). A simple yet striking consequence of hypotheses (i)-(iii) can be now directly inferred from Eq. (S4). In fact, as both J and P are state-dependent constants, this immediately implies that G[ρ(x)] = ψL −α x + ζ, i.e. G[ρ(x)] is a linear function of position x ∈ [0, L], with ψ = J √ m/P 3/2 the reduced current and ζ = G(ρ 0 ) a constant. Equivalently, where we have assumed that the function G(ρ) has a well-defined inverse F (u) ≡ G −1 (u). This assumption seems reasonable as steady density profiles are typically well behaved and readily measurable in simulations and experiments, see e.g. Fig. 1 in the main text. Therefore, according to Eq. (S7), there exists a single universal master curve F (u) from which any steady state density profile follows after a linear spatial scaling x = L α (u − ζ)/ψ. This scaling behavior is automatically transferred to temperature profiles via the local EoS P = ρ(x)T (x), so 2 FIG. S1. Theoretical predictions. Density (a) and temperature (b) profiles as a function of T0 for η = 1, obtained from the full solution of the macroscopic heat transport problem for the 1d diatomic hard-point gas, see Eq. (S10) and the associated discussion. Also shown are the η-and T0-dependence of (c) the scaled reduced current, L 1−α ψ, and (d) the nonequilibrium fluid's pressure P , see Eqs. (S12)-(S13). All these curves are for a mass ratio µ = 3, for which α = 0.297(6) and a = 1.1633(9), see Table S1 below where all measured anomaly exponents for different µ's can be found.
These scaling laws are independent of the global density η or the nonequilibrium driving defined by the baths temperatures T 0 and T L , depending exclusively on the function k(ρ) controlling the fluid's heat conductivity. Alternatively, Eq. (S7) implies that any measured steady density profile can be collapsed onto the universal master curve F (u) by scaling space by the scaled reduced current L −α J √ m/P 3/2 measured in each case and shifting the resulting profile an arbitrary constant ζ (similarly for temperature profiles). This suggests a simple scaling method to obtain the universal master curves in simulations or experiments, a procedure that we implement in the main text. Note that these results are not limited to the 1d diatomic hard-point gas; equivalent results hold for general d-dimensional (non-critical) fluids driven arbitrarily far from equilibrium, see Ref. [2] for a proof.
The combination of our scaling ansatz for the heat conductivity and well-known dynamical invariances of 1d hard point fluids under scaling of different magnitudes (as e.g. temperature, velocities, mass, space, etc.) results in a well-defined density dependence for the heat conductivity, see main text, namely k(ρ) = aρ α , with a a constant of O(1). Such power-law dependence, which reflects the transport anomaly, is fully confirmed in local measurements of the density dependence of κ L , see Fig. 5 in the main text, from which we obtain precise estimates of the amplitude a(µ), see Table S1. Such clear-cut observation, together with the scaling formalism described above, allows now for a complete solution of the macroscopic transport problem for this model, written in terms of the external control parameters, namely T 0 , T L , η and L, together with α and a. In fact, recalling that G (ρ) = k(ρ)ρ −5/2 we obtain that G(ρ) = ν * (1 − ρ α−3/2 ), with ν * ≡ a/( 3 2 − α) and where we have chosen an arbitrary constant such that F (0) = 1 = G −1 (0). The universal master curve hence reads This prediction is compared with the measured master curves in Fig. 3 (right panel) of the main text, and the agreement is excellent for all mass ratios µ. Eq. (S9) implies in turn that density profiles can be written as while temperature profiles simply follow from T (x) = P/ρ(x), namely (S11)  The calculation is completed by expressing the heat current J and the pressure P in terms of the external parameters by using Eqs. (S5)-(S6) above. This yields The last equation for the current can be rewritten as These predictions are fully confirmed by simulations data, see Fig. 6 in main text and Section II below. As a self-consistent check, note that in the equilibrium limit T 0 → T L both the pressure and the heat current converge to their expected value, namely P = ηT L and J = 0. Fig.  S1 shows the density and temperature profiles predicted for a macroscopic diatomic hard-point fluid as a function of T 0 for η = 1, as well as the pressure and the scaled reduced current L −α ψ as a function of T 0 and η. These plots are obtained for a particular mass ratio µ = 3, for which α = 0.297(6) and a = 1.1633(9), see Table S1, and yield an excellent comparison with simulation data, see Fig. 1 in the main text and Fig. S3 in Section II. Interestingly, the master curve F (u) obtained above exhibits a vertical asymptote at u = ν * , see Eq. (S9), and this implies in turn the existence of a maximal scaled reduced current ψ * . Indeed, for the associated density profile to exist in its whole domain x ∈ [0, L], see Eq. (S10), the following condition must hold with P expressed as in Eq. (S12). This defines a maximal scaled reduced current ψ * , such that the scaled current defining an upper bound on the heat current in terms of the nonequilibrium pressure. The maximal scaled reduced current increases monotonously with T 0 , saturating to an asymptotic value in the T 0 → ∞ limit, namely (S16) Note however that both L 1−α J and P diverge as T 0 → ∞, though ψ * remains finite.
To end this section, we remark that Eqs. (S9)-(S13) constitute the solution of the macroscopic transport problem for this model. A comparison of the predicted density and temperature profiles, see Eqs. (S10)-(S11), with the finite-size data of Fig. 1 in the main text allow us to investigate in the main text the bulk-boundary decoupling phenomenon in detail by quantifying the jump between the effective boundary conditions imposed by the boundary layers on the bulk fluid and the empirical bath temperatures.

II. ADDITIONAL RESULTS
In this Section we provide additional data, obtained from our extensive simulations of the 1d diatomic hardpoint fluid model, which support our conclusions in the main text.

A. Macroscopic local equilibrium, pressure and reduced current
Our first goal is to test the macroscopic local equilibrium (MLE) property directly from our data. As described above, MLE conjectures that local thermodynamic equilibrium holds at the macroscopic level, in the sense that the stationary density and temperature fields are locally related by the equilibrium equation of state (EoS) [1], which for this model is simply the ideal gas EoS, In order to test MLE, we hence take the density and temperature profiles of Fig. 1 of main text measured for µ = 3 and different T 0 , N, η, and plot in Fig. S2 the local reduced temperature, T (x)/P , with P the finite-size pressure measured in each simulation (see below), as a function of the associated local density ρ(x). All data, comprising 2400 points from 80 different simulations for widely different systems sizes, temperature gradients and global densities, collapse onto a single curve which follows with high precision the expected ideal-gas behavior 1/ρ, see line in Fig. S2 and inset therein. Note that, interestingly, the excellent data collapse is maintained also for points within the boundary layers near the thermal walls. Moreover, similar results hold for all mass FIG. S2. Macroscopic local equilibrium. Measured local reduced temperature, T (x)/P , plotted as a function of the associated local density ρ(x) for µ = 3 and ∀ T0, η, N , corresponding to all profiles displayed in Fig. 1 of the paper and summing up to 2400 data points from 80 different simulations. An excellent data collapse is obtained which follows with high precision the expected ideal-gas behavior 1/ρ, plotted as a thin line. Inset: Scaling plot of ρ(x)T (x)/P vs ρ(x) for the same conditions. These data show that macroscopic local equilibrium is a very robust property, even in the presence of strong finite-size corrections on the hydrodynamic profiles.
ratios µ studied in this paper. In this way, the observed high-precision data collapses confirm the robustness of the MLE property far from equilibrium [1], even in the presence of important finite size effects, validating in an independent manner one of the hypotheses underlying the scaling picture of Section I.
We next focus on the nonequilibrium fluid's pressure P and the heat current J flowing through the system, that we measure both in the bulk and at the thermal walls. These observables are necessary in order to scale the spatial coordinate of the hydrodynamic profiles using the measured reduced current ψ = J √ m/P 3/2 in each case. Fig. S3 shows the measured P (a) and ψ (b) as a function of T 0 and η for µ = 3 and different system sizes N . These data refer to wall observables, though the associated bulk observables yield completely equivalent results (as otherwise expected). The comparison of these data with our predictions in Section I is excellent, see Fig. S1 above.

B. Slow relaxation in 1d transport
In order to test the robustness of the measured anomaly exponents against order-of-magnitude changes in the system size, we have made a considerable computational effort to characterize the steady-state heat transport in the diatomic hard-point fluid for very large N 's, namely N = 31623 and N = 10 5 + 1 (see main text). Of course it would be desirable to go even beyond N = 10 5 + 1. However, it is important to note that for such very large values of N obtaining reliable results from simulations of 1d heat transport is exceedingly difficult. The underlying reason is not only that more particles need more computer time to simulate, but also that relaxation (and correlation) times increase greatly with N . This is due to the appearance of current (and momentum) waves in 1d which bounce back and forth between the thermal walls at a constant and well-defined speed while they are slowly damped away, a result of the dimensional constraint which strongly suppress local fluctuations. As far as we know, this remarkably slow relaxation mechanism has not been described yet in the literature on 1d transport, so we add details about the relaxation process and timescales next.
In particular, we have performed a large number of relaxation experiments for different values of N ∈ [10 2 + 1, 10 4 + 1], for µ = 3, T 0 = 20 and η = 1. Initial states for these experiments are randomly drawn from a local equilibrium measure corresponding to the macroscopic density and temperature profiles (obtained from extensive simulations for intermediate system sizes after a long empirical relaxation time). These initial states hence display on average the stationary density and temperature stationary profiles, but lack the weak but long-range correlations which characterize any nonequilibrium steady state (and which are in fact responsible for heat trans- port). To characterize the relaxation process to the true nonequilibrium steady state, we measured the average instantaneous energy current as a function of time, with the average taken over many different realizations of the relaxation process starting from the random initial states defined previously. Fig. S4 shows the relaxation of J(t) for different N as a function of time, measured in units of t 0 = M/(2T L η 2 ), the mean free time of a heavy particle in a cold environment. Note the logarithmic scale in the time axis. From this figure it is clear that the building of the faint but long-range correlations associated to the nonequilibrium steady state proceeds via the formation of a current wave which bounces back and forth between the thermal walls at a constant velocity, while being slowly damped in the nonequilibrium fluid. This damped current wave is accompanied by a similar momentum wave. The period of these waves scales linearly with the system size, and hence the relaxation time to reach the steady state also grows linearly with L. In fact, by scaling time by L −1 and the excess current by L α , with α(µ = 3) = 0.297, a good collapse is obtained (at least for large N ), see inset to Fig. S4, which suggest that with G(τ ) some scaling function.
The most relevant point here is that relaxation (and correlation) times grow linearly with the system size, making very difficult to obtain good statistics of transport for large enough N . In fact, the computer time needed to perform one collision per particle on average in an optimized event-driven molecular dynamics simulation of N particles scales as (N log N ) × τ (due to the cost of re-ordering the collision time queue in a heap data structure), with τ ∼ 10 −5 s the typical timescale of an elementary event. As the fluid relaxation and correlation times scale linearly with N , the computer time needed to obtain reliable data averages hence scales as t sim ∼ n exp × (N 2 log N ) × τ , with n exp the number of measurements to obtain good statistics. For N = 10 5 and n exp in the few hundreds (at least), we thus have the mind-blowing timescale t sim ∼ 10 8 s, right at the edge of modern-day multiprocessor computer power. For these reasons going beyond N ∼ 10 5 does not seem feasible at this moment.

C. Finite-size corrections for the effective anomaly exponent measured with standard methods
One of the most popular methods to measure the heat conductivity of a 1d fluid and characterize along the way the associated anomaly exponent consists in setting the model fluid with N particles and density η under a small temperature gradient, with fixed wall temperatures T 0,L , and then increasing N at constant density. For large enough N the overall temperature gradient is small enough so one can approximate Fourier's law as with ∆T = T 0 − T L , J the measured current and L = N/η. In this way, the estimated heat conductivity follows asκ ≈ JL/∆T , which is expected to diverge as N γ for large enough values of N (though there is no way of knowing a priori which value of N is large enough).
What it's typically found however in actual, cutting-edge simulations is an effective heat conductivity diverging as κ ∼ Nγ (N ) , with an effective anomaly exponent which exhibits itself persistent finite-size corrections [3]. Indeed, this approximate method completely neglects the nonlinear density and temperature dependence of the heat conductivity, and by construction it can only yield meaningful results in the limit N → ∞. It is therefore no surprise that the effective anomaly exponent derived within this approach for finite N varies slowly with the system size, as in fact the very definition ofκ above (and hence γ) is correct only asymptotically. This weakness in the above definition is reinforced by the fact that, for a given N , the heat conductivity measured in this way differs from estimations of κ using alternative approaches, as e.g. the also popular Green-Kubo equilibrium method [3] (which, by the way, is again exact only in the N → ∞ limit). Moreover, not only the estimated value of κ for a given N differs among different approaches, but also its scaling with N , and hence the estimation of the anomaly exponent.
The overall situation is therefore rather unsatisfactory, making extremely difficult to characterize reliably and accurately anomalous Fourier's law in 1d with these standard linear response methods. In fact, when measuring γ with standard methods one needs to reach huge system sizes, as large as N ∼ 10 5 , to appreciate certain convergence, and even in this case the asymptotic behavior is not yet clearly defined, see e.g. Ref. [3].
In contrast with the standard methods described above, our scaling approach takes full advantage of the nonlinear character of the heat conduction problem and provides a fully consistent and highly accurate description of all measured data in a broad range of parameters, including three orders of magnitude in N ∈ [10 2 + 1, 10 5 + 1], but also a wide range of temperature gradients (from the linear response regime to the fully nonlinear domain) and densities. The new anomaly exponent α(µ) that we conjecture equal to γ and is measured from the striking collapse of large amounts of data is welldefined, exceptionally robust and does not change with the system size in the broad range explored (see analysis in the main text). This contrasts with the running, effective exponent obtained from the standard linear response methods described above, which varies widely (and far beyond our precision limits) in the same N -range, i.e. γ(N ) ∈ [0.25, 0.5], see e.g. Fig. 3.b in Ref. [3].
A natural question now is: why do standard methods measure an effective anomaly exponent which converges so slowly with N ? The answer to this question lies at the bulk-boundary decoupling phenomenon reported in this work (see main text), i.e. the fact that the bulk of the finite-size nonequilibrium fluid behaves as a macroscopic, infinite system subject to some effective boundary temperatures, T (ef) 0,L (N ). Indeed, we measured and characterized the effective T (ef) 0,L ∀µ, N, η and T 0 by comparing the measured temperature profiles with the theoretical prediction based on our scaling theory, see left panel of Fig. 6 in the main text and the associated discussion. The effective boundary temperatures so obtained turn out to be N -dependent, slightly differing from the wall temperatures T 0,L but converging to these values as N increases. Most surprisingly, however, this convergence is exceedingly slow. In fact, Fig. S5 shows the measured relative temperature gap at the hot thermal wall, defined as (T 0 − T with Λ a small amplitude. Similar results were obtained for other mass ratios µ and at the cold boundary. This demonstrates that the effective boundary temperatures the bulk fluid feels (induced by the boundary layers) approach the wall temperatures as N increases at a extremely slow, ∼ 1/ √ log N rate, and this explains the persistent deviations found in the effective anomaly exponent in literature. More in detail, as explained above standard methods lead to an effective heat conductivitỹ κ ≈ JL/∆T ∼ Nγ for large enough N . Using now Eq. (S18) describing the slow decay of boundary finite-size corrections, we have that T L . In this way, for the effective heat conductivitỹ Noting now that ∆T (ef) is the real temperature gradient felt by the bulk fluid, we thus expect JL/∆T (ef) ∼ N γ , so inserting this in the previous equation and taking logarithms we arrive for large N at S6. A metric to quantify data collapse. Sketch explaining the metric used to quantify data collapse, see Eq.
(S19). This metric estimates the distance between a curve k (2) and the reference curvek ( ) by measuring the average distance between each point in k and the interpolated point ink with the same y-coordinate (gray, shaded squares). Note that we restrict to points in k overlapping with the reference curvek, see filled squares. The distance corresponds in this example to the average length of the dashed segments.
i.e. an effective anomaly exponentγ(N ) which converges at a exceedingly slow rate toward the correct, asymptotic anomaly exponent γ, in a way that closely resembles actual measurements, see e.g. Ref. [3]. This confirms that the slowly-decaying boundary finite-size corrections associated to the boundary layers are responsible of the strong, persistent finite-size effects affecting the effective anomaly exponent measured with the standard linear response method. Moreover, as our scaling method is independent of the boundary temperatures driving the system out of equilibrium, this explains why our results for the new anomaly exponent α, that we conjecture is equal to the true asymptotic exponent γ, are free of these persistent finite-size corrections.

III. A METRIC TO QUANTIFY DATA COLLAPSE
In this section we briefly explain the standard metric used in this work to quantify data collapse. This metric is based on the collapse distance first proposed in Ref. [4] and widely used in physics literature, in particular in order to obtain scaling exponents via a distance minimization procedure.
We hence consider a set of K curves, each one containing M points, and we denote this set as {{(x The idea is now to choose an arbitrary curvek ∈ [1, K] as reference curve, and proceed to measure the distance of all other curves k =k to this reference curve along the x-direction. For that we measure the distance between each point in k and the interpolated point ink with the same y-coordinate. In order to do so, we have to restrict to points in k overlapping with the reference curvek. Note also that we choose to measure distances only along the x-direction because the scaling approach developed in this paper only affects the x-coordinates of the measured curves, see Section I and Figs. 2 and 3 (right panel) in the main text. Moreover, since the chosen reference curvek is completely arbitrary, we repeat this procedure for all curves as reference curve, and average the resulting distances. In this way, our collapse metric is defined as [4] D ≡ 1 i ) of curve k on curvek along the xaxis. The innermost sum over i in Eq. (S19) is restricted to points in curve k which overlap with curvek along the y-direction, i.e. those points in k whose y-coordinate is between the minimum and maximum y-coordinate of curvek. In order to obtain now the projectionx (k,k) i in Eq. (S19) any interpolation scheme can be used, though for our purposes the simplest linear interpolation works well. In particular, we choosē with A (k,k) i and B (k,k) i the slope and the y-intercept of the interpolating function, The points (x (k) i± , y (k) i± ) correspond to the points in thekcurve bracketing point i of k-curve along the y-direction, see sketch in Fig. S6. To normalize the distance metric, we divide the resulting sums by the total number of overlapping points, N overl . Moreover, because the L −α scaling in the x-coordinate of the measured density and temperature profiles may affect strongly the total span of the collapsed curves depending on the anomaly exponent α used, the collapse metric is also normalized by the total span in the x-direction of the curve cloud, max ≡ (x max − x min ) with x max = max k,i [{x In order to obtain the exponent α characterizing anomalous Fourier's law in our 1d fluid, we minimize the metric (S19) for varying mass ratios µ. In fact, the collapse metric D(α, µ) exhibits a deep and narrow minimum as a function of α for each µ, see inset to Fig. 3 in the main text, offering a precise measurement of the anomaly exponent. Moreover, an estimate of the exponent error can be obtained from the width and depth of this minimum [4]. By expanding ln D(α, µ) around the minimum at α = α 0 , the width can be estimated as [4] ∆α = α 0 2 ln D(α 0 ± α 0 , µ) D(α 0 , µ) , (S21) for a given level . Here we choose = 0.01, so the estimate for the anomaly exponent is α 0 ± ∆α with an errorbar reflecting the width of the minimum at the 1% level [4].