Nonequilibrium Green-Kubo relations for hydrodynamic transport from an equilibrium-like fluctuation-response equality

Near equilibrium, Green-Kubo relations provide microscopic expressions for macroscopic transport coefficients in terms of equilibrium correlation functions. At their core, they are based on the intimate relationship between response and fluctuations as embodied by the equilibrium fluctuation-dissipation theorem, a connection generically broken far from equilibrium. In this work, we identify a class of perturbations whose response around far-from-equilibrium steady states is linked to steady-state correlation functions via an equilibrium-like fluctuation-response equality. We then utilize this prediction to substantiate linearized hydrodynamic transport equations that describe how spatial inhomogeneities in macroscopic nonequilibrium systems relax. As a consequence, we derive nonequilibrium Green-Kubo relations for the transport coefficients of two types of hydrodynamic variables: local conserved densities and broken-symmetry modes. A byproduct of this work is to provide a theoretical foundation for the validity of Onsager's regression hypothesis around nonequilibrium steady states. Our predictions are analytically and numerically corroborated for two model systems: density diffusion in a fluid of soft, spherical active Brownian particles and phase diffusion in the noisy Kuramoto model on a square lattice.


I. INTRODUCTION
The fluctuation-dissipation theorem (FDT) is a cornerstone of statistical mechanics [1][2][3][4]. It quantifies a fundamental correspondence between two experimental procedures for interrogating an equilibrium system: measuring a system's response to a weak perturbation and passively observing its fluctuations both contain the same information. A major consequence has been in refining our understanding of the material coefficients that determine how spatial inhomogeneities in near-equilibrium macroscopic systems relax via hydrodynamic transport. The resulting predictions, known as Green-Kubo relations, equate these macroscopic transport coefficients D to the microscopic equilibrium correlation functions of local current observables j r (t), whose fluctuations depend on space r and time t in a volume V [4][5][6][7][8][9][10]: Here, β is the inverse temperature and χ is the static susceptibility (or thermodynamic derivative). Even for far-from-equilibrium steady states, response can still be related to a nonequilibrium correlation function by a modification of the FDT's original derivation [11,12]. The resulting correlation function turns out to be rather formal, often requiring detailed microscopic knowledge of the steady state or its dynamics [13][14][15][16][17]. Even still, this nonequilibrium modification of the FDT allows one to link integrals of nonequilibrium correlation functions to microscopic currents [18][19][20], in much * jmhorow@umich.edu the same spirit as the macroscopic Green-Kubo relation in (1).
Nonequilibrium Green-Kubo relations that relate macroscopic transport coefficients to local current observables have appeared in the literature in particular situations. These studies can be categorized according to the method. For a two-dimensional nonequilibrium viscous fluid [21,22], Green-Kubo relations were deduced by assuming that Onsager's regression hypothesis [23] remains valid around nonequilibrium steady states. An alternative approach utilizes the projection operator method [24,25] adapted for non-Hamiltonian dynamics [26][27][28][29]. The resulting Green-Kubo relations incorporate a time-reversed dynamics, apparently obscuring the interpretation of the resulting correlation functions. This obstacle has been overcome for at least one specific model of a nonequilibrium active fluid [30].
In this work, we demonstrate that generally Green-Kubo relations for macroscopic transport coefficients maintain their equilibrium form arbitrarily far from equilibrium. The key step in this derivation is the elucidation of a class of perturbations with explicit conjugate variables whose response is given by simple nonequilibrium correlation functions, akin to the equilibrium FDT. This observation generalizes our previous work on static response to time-dependent perturbations [31]. We then exploit this equilibrium-like fluctuation-response equality to provide a theoretical foundation for linearized hydrodynamic equations governing transport in homogenous nonequilibrium fluids. The technique we use was originally developed by Kubo to analyze conductivity [6], but elaborated for simple equilibrium fluids by Oppenheim and collaborators [32][33][34][35]. We consider two kinds of slow hydrodynamic modes, local densities of conserved variables and the Nambu-Goldstone modes that emerge from Active Brownian Particles Dynamical generator: Perturbation: FIG. 1. Perturbing the dynamics of a fluid of active Brownian particles: Particles (red) are self-propelled with a velocity v0 whose direction θ diffuses. They interact pair-wise through a repulsive force F and experience translational noise. The surrounding environment induces diffusive behavior of the particles characterized by diffusion coefficients Dt and Dr. The stochastic equations of motion (21) correspond to the Fokker-Planck operator L. The type of perturbation in (3) corresponds to multiplying all system parameters by 1 − λQ.
a broken continuous symmetry [4]. Our theory is illustrated by numerical simulations of two examples, a fluid of soft active Brownian particles and a noisy Kuramoto model.

II. FLUCTUATION-RESPONSE EQUALITY
Consider a classical system with microscopic state z, whose dynamics, either stochastic or deterministic, are well-modeled as Markovian. As such, the system's probability density P (z, t) evolves according to the Master equation [36] ∂ t P (z, t) = LP (z, t), with generator L. We will further assume that left undisturbed the system will relax to a unique steady-state distribution P ss (z) given by the solution of LP ss = 0. This model incorporates deterministic Hamiltonian dynamics, where LP = {P, H} is the Poisson bracket with a Hamiltonian H, as well as non-Hamiltonian deterministic and stochastic dynamics. A stochastic nonequilibrium example that will be a useful illustration is the fluid of active Brownian particles (ABPs) pictured in Fig. 1, where L is the Fokker-Planck operator [36]. Our focus is on how averages of observables O(t) = O(z)P (z, t)dz change in response to weak perturbations. While this can be analyzed quite generally [11,12,14,15], we have identified a specific class of dynamical perturbations with interesting and useful properties of the form with a small time-dependent parameter λ(t) and a timeindependent conjugate coordinate Q(z). For an equilibrium system evolving under Hamiltonian dynamics, the perturbation takes the form λL(QP ) = λ{QP, H}, which at first glance appears to differ from the classical linear response perturbation λ{P, Q} generated by modifying the Hamiltonian H = H − λQ. However, at equilibrium, the distribution is solely a function of the Hamiltonian, like for example the canonical distribution P eq ∝ e −H (with β = 1). In this case, both perturbations have the same effect on an equilibrium distribution, which can be seen by noting that {QP eq , H} = {Q, H}P eq = {P eq , Q}, since {P eq , H} = 0. Thus, this approach encompasses equilibrium linear response theory. Clearly, multiple perturbations can then lead to the same response, and in Appendix A, we delineate in detail this equivalence class. For nonequilibrium dynamics, implementing the perturbation embodied in (3) can require a rather complicated coordinated change in the system's parameters, as can be seen in Fig. 1. Now, imagine starting the system in its nonequilibrium steady state at t = −∞ and then turning on the perturbation. A standard first-order perturbation analysis [11,12,37] detailed in Appendix A reveals that the average of an observable O(t) at time t over a statistical ensemble initially at steady state will deviate from the steady-state average O ss = O(z)P ss (z)dz by Here, we have introduced the time-translationallyinvariant steady-state correlation function, defined for any two observables, O 1 (z) and O 2 (z), by in terms of the transition probability P (z, t|z , s) obtained as the solution of (2) with a delta-function initial condition δ(z − z ).
Noteworthy is that the response to these perturbations is given by the steady-state correlation between the observable O and the conjugate observable Q that explicitly appears in the perturbation (3), just like the equilibrium FDT. By contrast, previous studies that addressed generic perturbations, arrived at correlation functions where instead ofQ, the conjugate observable that appeared required knowledge of the entire steady-state distribution (via ∂ λ ln P ss ) or of the generator of the dynamics [11,12,14,15].
While executing such a perturbation in the lab will generically be challenging, the theoretical utility is born out in the simplicity of the response function and the applications that entails. Firstly, the fluctuation-response equality in (4) implies we can measure any correlation function we choose by observing the appropriate response (3). We can exploit this in a computational experiment, where there is no issue applying any perturbation we choose. Such an approach would be the far-fromequilibrium equivalent of the "nonequilibrium perturbation method" utilized for measuring equilibrium correlation functions in simulation [38]. The second major application, and the one we focus on for the rest of this paper, is to derive Green-Kubo relations. Some of the earliest derivations of Green-Kubo relations exploited the equilibrium FDT [8,9,32]. It was very quickly realized that the external perturbations played only a formal intermediary role in the derivation and therefore need not be realizable in the lab [32]. A concrete example of this point is the derivation of the Green-Kubo relation for the shear viscosity of a simple fluid due to Jackson and Mazur where a non-gradient force was introduced into a Hamiltonian dynamics to mimic a shear [39]. With this observation, we now employ this previously developed program to study transport in macroscopic systems, adapting it to nonequilibrium steady states using our fluctuation-response equality (4).

III. HYDRODYNAMIC TRANSPORT
We now focus our attention on an N -particle macroscopic system in a volume V . We will assume that its macroscopic dynamics on long length-and time-scales can be accurately captured by a few hydrodynamic variables, which are many-bodied observables with a relaxation time that diverges with system size. Typically they are local densities of conserved quantities or Nambu-Goldstone modes due to broken continuous symmetries [4]. For example, the near-equilibrium dynamics of a simple fluid are completely captured by the conserved number, momentum, and energy densities [4]. Away from equilibrium, hydrodynamics has been quite successful in modeling active matter [40]-systems composed of constituents that individually consume energy. Theories of flocking [41], active gels [42], and experiments on cellular spindle dynamics [43] further reinforce the need to describe both conserved densities as well as brokensymmetry modes. For clarity of presentation, we will focus on systems that possess a single hydrodynamic variable A r (t); the generalization to multiple hydrodynamic variables is presented in Appendix B.
In the homogeneous steady state, our hydrodynamic variable will be spatially uniform and constant in time with valueĀ. Upon the emergence of a small spatial inhomogeneity, either due to an external stimulus or internal fluctuation, the deviation from the steady-state δA r (t) = A r (t) −Ā will by assumption relax (or regress) via a linear hydrodynamic equation with drift v and transport coefficient matrix D, ∂ t δA r (t) = −∇ · (v δA r (t)) + ∇ · D · ∇δA r (t), (6) or in Fourier space (δA k (t) = V δA r (t)e ik·r dr), Near equilibrium when there is a single hydrodynamic variable v = 0, due to time-reversal symmetry, as there can be no preferred direction of motion. To have a nontrivial drift v = 0 near equilibrium, multiple hydrodynamics are required, in which case v is called the Eulerian term or the reactive coupling [44]. Nonequilibrium steady states, by contrast, can support spatial flow of a single variable and therefore we allow for a nonzero v here. Equation (6) (or equivalently (7)) serves as the definition of the parameters v and D. As such, at the macroscopic hydrodynamic level, they can only be measured by first setting up a spatially inhomogeneous profile and fitting the subsequent relaxation to Eq. (6) (or (7)).
We can now view (7) as the beginning of a long wavelength (small k = |k|) expansion of a generalized linear transport equation modeling the exponential regression of the hydrodynamic variable: where M k is a generalized transport coefficient that for small-k must behave as M k ≈ −ik · v + k · D · k + · · · . In the following section, we will use our equilibriumlike fluctuation-response relation to provide a statisticalmechanical basis for this expansion and as a result extract microscopic expressions for v and D, known as Green-Kubo relations. But first we need to connect the macroscopic to the microscopic. To this end, we specify our hydrodynamic system's microscopic state. For each of the i = 1, . . . , N particles, we denote its position as r i and any other of its degrees of freedom, like momentum or polarity, as s i , such that z i = (r i , s i ) and z = (z 1 , · · · , z N ). Microscopic expressions for hydrodynamic variables are then typically formed as densities of single-particle observables a j (z j ) through a sum of the form [4] A r (z) = N j=1 a j (z j )δ(r − r j ).
For example, the choice a j = 1 defines the particle number density. Notice that the {Â r (z)} r∈V are a family of state-space observables parameterized by the spatial location r. The link to the macroscopic hydrodynamic variable then emerges on average, Â r (t) = A r (t).
The two types of hydrodynamic variables we will consider-conserved densities and Nambu-Goldstone modes-are identified by their microscopic dynamics. For the conserved density, its evolution along a single dynamical trajectory is distinguished by a continuity equation ∂ tÂr + ∇ r · j r = 0, with corresponding local current j r . Re-expressing this continuity equation in Fourier space, ∂ tÂk = ik ·j k , is particularly illuminating, since it makes explicit that the small-k modes vary slowly in time with a relaxation time that diverges as k → 0. By contrast, Nambu-Goldstone modes need not be conserved. Instead, they are distinguished by diverging static fluctuations.

IV. GREEN-KUBO RELATIONS
The generalized transport equation in (8) describes the relaxation from an initially inhomogeneous profile. To connect this relaxation to the microscopic dynamics, we follow the procedure laid out in [32][33][34]. Underpinning this procedure is the conceptual insight that if we start to slowly turn on an external perturbation conjugate to A r from t = −∞, by t = 0 we have generated an inhomogeneous profile that is minimally disturbed from the steady state. At t = 0, we switch off the perturbation and track the subsequent evolution. If this microscopic experiment is to be consistent with the macroscopic transport equation, the two must both give the same evolution on long length-and time-scales. This consistency requirement leads to nonequilibrium Green-Kubo relations connecting the relaxation dynamics to fluctuations.
Our first step is to specify the coordinate in our perturbation (3). The required form is where we eventually take → 0 + , and the Heaviside step function Θ(−t) turns off the perturbation at t = 0. The choice of f r is immaterial as long as all integrals converge, as it will shortly drop out of the calculation. Using this choice of conjugate coordinate in the fluctuation-response equality (4) leads to an expression for the evolution ofÂ for t ≥ 0, which after an integration by parts and sending → 0, becomes We now remove the dependence on f r in favor of the nonuniform profile generated by this perturbation at t = 0. This is facilitated by first Fourier transforming (11) and exploiting the assumed translational invariance of the homogenous steady state to arrive at in the large system-size limit [32]. This equality is true at every positive time, including at t = 0, which allows us to solve for f k = V δÂ k (0) / Â k (0)Â −k (0) ss in terms of the initial value δÂ k (0) and substitute back in to find Equality (13) demonstrates that the time evolution of the relaxation of the average is identical to that of the correlation function, at least for the particular ensemble produced at t = 0 by this perturbation. Importantly, equality (13) depends only on the behavior of the sole slow hydrodynamic variable. (Here, there is just one variable, but in general we need to include all slow degrees of freedom as described in Appendix B.) Thus, at long enough length-and time-scales, any vagaries of the initial preparation will die out exponentially fast, and (13) should continue to remain true when the generalized transport equation in (6) accurately describes the exponential relaxation of A r (t) = Â r (t) . From this we can conclude that the steady-state correlation function of the hydrodynamic mode is governed by the same macroscopic linear equation, at least for long times and small k. This prediction is a key contribution of this study as it provides a statisticalmechanical rationale for an Onsager's regression hypothesis around nonequilibrium steady states, as was conjectured in [21,22]. Consistency between the macroscopic and microscopic descriptions embodied by (14) leads to the conclusion that M k can be inferred from the dynamics of steadystate correlation functions [32][33][34]. Thus, we can exploit (14) to solve for the generalized transport coefficient systematically for small k, when we expect hydrodynamics to be an accurate description, and thus deduce the expansion parameters M k ≈ −ik · v + k · D · k + · · · . The steps are detailed in Appendix B, but follow closely [32][33][34]. Here, we report the results.
The precise form of the resulting Green-Kubo relations depends on the specifics of how the hydrodynamic variable's microscopic dynamics behave for small wave numbers. Let us first address conserved local densities whose microscopic relaxation rate diverges with small k, ∂ tÂk = ik · j k , but whose static correlation function remains finite,χ = lim k→0 Â kÂ−k ss . The static correlation functionχ is related to the static structure factor by dividing by an extensive quantity such as V or N . In equilibrium,χ/V can be related to a thermodynamic susceptibility via the equilibrium FDT, which can then be deduced solely from the equilibrium equation of state [4,8]. Equality (12) at t = 0 shows thatχ can still be connected to the static response induced by (3). However, such an interpretation loses much of its utility away from equilibrium due to the absence of any macroscopic nonequilibrium thermodynamics. Even still, recent work suggests there might be a useful thermodynamic structure, at least for spherical ABPs [45].
Carrying out the expansion of the generalized transport coefficient, we find for the drift which captures the conservative (or Eulerian) part of the dynamics. The transport coefficient then captures the fluctuations around this steady drift, which motivates the introduction of the dissipative current as the relative flow: I r = j r − vÂ r . The dissipative current I r represents the rapidly fluctuating part of j r , in light of the orthogonality lim k→0 I kÂ−k ss = 0, suggesting that its fluctuations decay on microscopic timescales [34]. The resulting Green-Kubo relation is most simply expressed using the half-Fourier transform in time I kω = ∞ 0 I k (t)e iωt dt: where the first term is the dissipative-current correlation function and the second is the matrix of first-order corrections to the drift The dissipative-current correlation function is assumed to decay fast enough so that the time-integral converges to a finite value. As such, we are assuming that the plateau value problem, a long-standing problem of Green-Kubo relations [24,[46][47][48][49], does not occur in systems of interest.
Comparison with the typical equilibrium expression in (1) brings to light two apparent differences. The first is the use of the dissipative current I r instead of the current j r , and the second is the appearance of E. Both terms had already appeared in early works on near-equilibrium transport [33,34,39], and were rediscovered for hardcore and Langevin dynamics [26,27,50], an observation we recount in more detail in Sec. VI. Historically E was set to zero since early studies focused on nearequilibrium systems evolving with deterministic Hamiltonian dynamics where time-reversal symmetry forces j kÂ−k ss = 0 [6,39]. However, if the dynamics are stochastic it is possible for E = 0 even in time-reversalsymmetric equilibrium systems, as we will see for an equilibrium fluid of Brownian particles.
Nambu-Goldstone modes, by contrast, are identified by the divergence of their static correlation function Â kÂ−k ss as k → 0. Near equilibrium, rigorous arguments based on Goldstone's theorem and the Bogoliubovinequality lead to a minimum divergence of 1/k 2 [4]. However, for stochastic nonequilibrium dynamics, the microscopic foundation of these modes is shakier. The connection between symmetries and conservation laws arising from Noether's theorem is not robust enough for stochastic dynamics to encompass typical nonequilibrium matter [51]. Without that connection there are no general grounds for a nonequilibrium Goldstone theorem. Progress has been made on specific models using a field-theoretic approach [52], whose general applicability is unclear to us. Even still our approach allows us to demonstrate that diverging static fluctuations will generically lead to hydrodynamic modes in a nonequilibrium system. To be as agnostic as possible, we take a generic divergence of the form Â kÂ−k ss ∼ k −q with q arbitrary. To have a consistent small k expansion of (14) when the static correlation function diverges, the correlation functions Ȧ kωȦ−k ss and Ȧ kÂ−k ss must scale in such a way that the limits defining the drift and the leading order transport coefficient are finite: These expressions are actually the most general forms for v and D, and specialize to (15) and (16) when the hydrodynamic mode is conserved.

V. ILLUSTRATIONS
The Green-Kubo relations in (15), (16), (19), and (20) predict that two distinct experimental procedures must give the same result: (i) measure the transport parameters v and D directly via their defining equation (6) by perturbing the system and then measuring the rate of exponential relaxation, or (ii) extract the same information by passively observing the steady-state fluctuations. In this section, we corroborate and illustrate this equivalence. Our first model is a fluid of ABPs with a single conserved density, and the second is a noisy Kuramoto model, where the breaking of the continuous rotational symmetry upon synchronization leads to a Nambu-Goldstone mode.

A. Active brownian particles
We first consider a fluid of N spherical ABPs in a twodimensional box of size L × L, with periodic boundary conditions. Interactions are modeled through a repulsive, short-ranged, pair potential φ(|r i − r j |). Activity enters by each particle being self-propelled with a velocity v 0 e(θ i ) = v 0 (cos θ i , sin θ i ), whose orientation θ i diffuses with diffusion coefficient D r . Including translational noise with diffusion coefficient D t leads to an evolution governed by the pair of overdamped Langevin equations [53,54] where ξ i and η i are independent Gaussian white noises, µ is the bare mobility, and is the the total force acting on the i-th particle due to pair interactions.
The only conserved variable is the total number of particles. Momentum and energy are not conserved due to the self-propulsion and noise. Accordingly, the only hydrodynamic variable is the local particle density ρ r = i δ(r − r i ), with steady-state averageρ = N/L 2 . Since the particles do not prefer any particular direction, the density transport exhibits an unbiased isotropic diffusion with v = 0 and transport coefficient proportional to the identity matrix D = DI: with a diverging relaxation time τ (k) = 1/(k 2 D).
We first validate the Green-Kubo relation (16) through direct numerical simulation. This is accomplished by first measuring the transport coefficient D as the rate of exponential relaxation from an inhomogeneous initial condition per its definition and then comparing that to the prediction of the Green-Kubo relation. As detailed in Appendix C, we use a harmonic interaction potential φ(r) = (K/2)(r − a) 2 Θ(a − r), with interaction strength K and length a, and we set D t = 0 to enhance the nonequilibrium effects and facilitate the numerical analysis. To measure D directly, we first observed the evolution of the Fourier modes ρ k from a nonuniform initial condition with all the particles localized to the band L/4 ≤ x i ≤ 3L/4, as depicted in Fig. 2(a). Figure  2(b) displays the expected exponential relaxation of a hydrodynamic mode whose slope is proportional to k 2 D. Repeating this experiment 10 times and then averaging gives us our estimate of D and its error (Fig. 2(c)). Next, using a single long steady-state simulation, we estimated D from the steady-state correlation functions appearing in the Green-Kubo relation. To gain further insight, we complement the numerical validation with an analytic analysis. The microscopic dynamics of the density ρ r is given by a nonlinear stochastic differential equation with multiplicative noise, known as the Dean equation [55,56]. It is further coupled to the self-propulsion orientation through the polarization density p r = i e(θ i )δ(r − r i ) as well as higher harmonics. To make analytic progress we take the limit of a dense, weakly-interacting fluid. Following [56], we then linearize these equations about the steady state (ρ andp = 0), and truncate the orientation harmonics at p r , since higher harmonics do not contribute in the limit k → 0. The resulting linear Gaussian dynamics derived in Appendix C can be solved analytically allowing us to determine each term in the Green-Kubo relation, which are all diagonal (C = CI, E = EI) with elements where φ 0 = lim k→0 φ k is the zero wave-vector limit of the pair potential. Combining, we arrive at a prediction for the transport coefficient Weak-interactions and activity enhance the diffusion. We further see for an equilibrium Brownian fluid where v 0 = 0, the current-current correlation function vanishes (C = 0), and the transport coefficient is determined solely by E. Thus in stochastic models even in equilibrium, E is required for an accurate prediction of the transport coefficient.
The predictions from this linearized theory (24) are compared to the simulation in Fig. 2(c) where the dashed line is for the weaker interaction (K = 0.5) and the dotted is the stronger interaction (K = 1.0). For the weaker interaction (K = 0.5), the linear approximation agrees well with the simulations when v 0 is large. When v 0 is small, the predictions of the linear theory for both the strong and weak interactions fall outside the error bars, overestimating the transport coefficient. Even still, the coincidence of the two numerical measurements of D suggest the macroscopic Green-Kubo relation remains valid. This illustrates how even when the microscopic dynamics are nonlinear, the emergent macroscopic dynamics can still be linear, as emphasized by Van Kampen [57].

B. Stochastic Kuramoto model
Our second illustration is a variant of the Kuramoto model, which is a canonical model for synchronization among a large collection of phase oscillators [58]. In the original Kuramoto model, all oscillators interacted with each other and evolved deterministically [59], though later variations included a variety of modifications such as allowing for dynamical noise [60]. Here, we consider a noisy version with N = L x L y oscillators evenly spaced on a two-dimensional square lattice of size L x × L y with nearest-neighbor interactions and periodic boundary conditions. The time evolution of the i-th oscillator's phase θ i is described by the stochastic equatioṅ where the sum extends over neighbors of i, and K is the interaction strength. There are additionally two sources of noise. The oscillator's intrinsic frequency Ω i is sampled from a Gaussian distribution of meanΩ and variance σ 2 . Each oscillator also experiences white Gaussian dynamical noise ξ i with a strength characterized by T . When the strength of noise, σ and T , is weak and the interaction K is strong, a finite system synchronizes and all the oscillators rotate coherently at a common frequencyΩ. The effect of the noise on this synchronization transition has recently been investigated in [61].
In the synchronized state, the rotational symmetry of the equations of motion (25) is spontaneously broken. As a result, we expect that localized perturbations in phase will relax diffusively via a Nambu-Goldstone mode. We can see this explicitly if we consider a small perturbation to the synchronized state θ i (t) = θ 0 +Ωt+δθ i (t), where θ 0 is an arbitrary offset that emerges when the symmetry is broken. Linearizing the interaction in (25) by assuming neighboring phases are close (|δθ i −δθ j | 1), and Fourier transforming (θ k = j θ j e ik·rj ), shows that each mode oscillates independentlẏ for small k. Upon averaging over the noises, we find that the average phase relaxes as In other words, the phase regresses diffusively with relaxation time τ (k) = 1/(k 2 K) characterized by a transport coefficient D = K that does not depend on the strength of the noises. The transport coefficient D can also be deduced via the Green-Kubo relation (20). The exact calculation based on the linear equation (26) gives for small k, Consequently, the Green-Kubo relation (20) yields the consistent prediction D = K. See Appendix D for more details on the linearized theory.
To corroborate this analysis, we estimated the transport coefficient by two independent numerical simulations. We first setup an initial nonuniform profile of phase by rotating the leftmost tenth of oscillators ahead by π/2 and observed the subsequent regression of the first few smallest k modes θ k in the x-direction, as illustrated in Fig. 3(a). Figure 3(b) displays the expected linear regression of θ k , averaged over 100 realizations, whose slope we use to extract D = 1/(k 2 τ (k)) (inset). Steady-state simulations are then used to deduce D from the correlation functions in the Green-Kubo relation (20). Figure 3(c) shows the coincidence of the two different measurements of D in comparison to the prediction from the linearized theory (dotted line). As long as the strength of the noises is weak enough for the system to reach the synchronized state, the transport coefficient D is only determined by K. See Appendix D for thorough information about the simulation methods and the error analysis.

VI. COMPARISON WITH EARLIER RESULTS
The near-equilibrium FDT and Green-Kubo relations are deep, long-standing results, and consequently have been extensively analyzed and extended in various directions. In this section, we compare and contrast our present developments with pertinent earlier literature.

A. Fluctuations and response
In Sec. II, we identified a class of dynamical perturbations whose response can be identified as a simple correlation function, akin to the equilibrium FDT. Graham has also identified a group of perturbations that result in an equilibrium-like fluctuation-response equality [62]. Here, we elucidate the connection between these two approaches. To make this connection, let us take a look at the steady-state distribution of our perturbed dynamics in Eq. (3) with external control parameter fixed λ(t) =λ. The resulting steady-state distribution, given as the solution of L P ss (z) = 0, is then approximately accurate up to first order inλ. We see that the steadystate of the perturbed dynamics is approximately an exponential re-weighting of the unperturbed steady-state P ss (z) by the conjugate coordinate Q.
For nonequilibrium systems modeled using a Fokker-Planck equation, Graham previously demonstrated that a perturbation that changes the steady-state via an exponential re-weighting as in Eq. (29) when applied dynamically will satisfy an equilibrium-like FDT, equivalent to Eq. (4). Graham's analysis, however, was performed by first deriving a covariant form of the Fokker-Planck equation and then identifying a decomposition of the dynamics where it is more natural to work with exponential shifts in the steady-state distribution. As noted in Ref. [63], it is perhaps the rather formal mathematical language of Ref. [62] that has obscured this earlier prediction.
Our analysis complements this previous work in two important ways. Our approach is valid for any Markovian dynamics and is not restricted to continuoustime continuous-space diffusion processes modeled via a Fokker-Planck equation As such, we identify the pertinent perturbations using the original governing Master equation, without recourse to a covariant formulation. This allows us to identify the changes in system parameters directly, as in Fig. 1.

B. Green-Kubo expressions
Green-Kubo relations for non-Hamiltonian dynamics -such as hard-core classical fluids or systems described by stochastic Langevin equations -have been deduced using the Mori-Zwanzig projection operator method independently by Ernst and Brito [26,27] as well as by Español [28,29]. While those authors explicitly only considered systems near equilibrium, their arguments trivially extend to homogenous nonequilibrium fluids like those considered here.
The work of Ernst and Brito most closely parallels the present analysis. They projected directly onto the linear Langevin equation for the hydordynamic variables. Español, by contrast, employed the projection operator method to obtain the nonlinear Fokker-Planck equation describing the evolution of the full distribution of the hydrodynamic fluctuations. From this more general approach one can in principle obtain an equation for the linear regression of the hydrodynamic variables through a suitable expansion, as detailed for example by Zwanzig [25,64]. Thus, both analyses lead to similar conclusions. In this section, we demonstrate how our linear response theory derivation of the Green-Kubo relations agrees with the relations obtained using the projection operator method. For simplicity and clarity of presentation, we will assume in this section there is a single conserved hydrodynamic variable and j k A −k = 0, so that v = E = 0.
To state the prediction obtained from the projection operator method, we first need to introduce a couple of concepts. The first object we will need is the generator of the time-reversed [26][27][28][29] or dual dynamics [65,66] where L † is the adjoint operator of L defined by f (z)L † g(z)dz = g(z)Lf (z)dz. We have also assumed that there are only even variables under time-reversal. If there are odd variables, one must additionally reverse their sign by including the time-reversal operator (see Ref. [67] for an exposition in the quantum context). In the dynamics generated byL, every trajectory of the original dynamics appears with the same probability except run in reverse [37]. WhenL = L, the dynamics are said to be detailed balance, and every trajectory occurs with the same probability as its time-reverse [65]. Apart from this interpretation, we will also find useful the following property of the dual generator, which allows us to "commute" the generator L past the steady-state distribution at the expense of introducing the adjoint of the dual generator. With the generator and its dual in hand, we can use them to define timedependent hydrodynamic variables in the Heisenberg picture going forwards and backwards in time as Here, the superscript H stands for Heisenberg picture.
We have also taken this opportunity to introduce local currents for both the forward j H k and time-reversedj H k dynamics. This should be contrasted with the conservation equations introduced earlier, where the derivative was taken along a single dynamical trajectory. The projection operator method then leads to the following Green-Kubo relation for the transport coefficient [26][27][28][29] where the subscript on the average O Pss = O(z)P ss (z)dz emphasizes that this is a static average over the steady-state distribution as opposed to the average over steady-state dynamical trajectories employed in the rest of the paper, denoted by · ss . Much has been made over the fact that the correlation function in Eq. (34) contains a mixture of the forward and timereversed dynamics [26][27][28][29]. Indeed, it leads one to believe that to measure this correlation function one needs access to both the forward and the time-reversed dynamics in the same experimental setup. Admittedly, this appears challenging in an experiment, though the effect of timereversal can be readily be extracted computationally, as was done in Ref. [50].
The equality in Eq. (34) appears identical to our expression for the Green-Kubo relation (16) except under the integral is a two-time average of Heisenberg operators going in two time directions as opposed to the forwardin-time current-correlation function in Eq. (16). In fact, these two expressions are mathematically identical, and Eq. (34) is simply the Heisenberg representation of a cur-rent correlation function. To see this, observe that where in the last line we used Eq. (31). Next, we can recognize the resulting operator expression as the time derivative of the steady-state correlation function upon comparison with Eq. (5), to find Thus, the Green-Kubo relations obtained using the projection operator method contains the same currentcorrelation functions as obtained here, but expressed in terms of the Heisenberg picture.
When the dynamics are stochastic or not time-reversal symmetric, it is not well appreciated that when standard correlation functions are expressed in the Heisenberg picture they include the time-reversed or dual dynamics. This point may have interfered with the widespread understanding that the projection operator method when applied to the study of nonequilibrium fluids leads to predictions equivalent in form to equilibrium.

VII. CONCLUSION
We have identified a class of perturbations whose response verifies an equilibrium-like fluctuation-response equality. This allowed us to rationalize Onsager's regression hypothesis around nonequilibrium steady states and served as the foundation of a systematic method to extract linearized hydrodynamic transport equations around homogenous nonequilibrium steady-states. The key resulting predictions were Green-Kubo relations that link hydrodynamic transport coefficients to steady-state fluctuations, which were verified in two models both analytically and numerically. As for the traditional equilibrium Green-Kubo relations, our approach here complements results for non-Hamiltonian and stochastic systems based on the projection operator method [26][27][28][29]. In the projection operator method, one must conjecture a criterion to single out a projected state that contains all relevant macroscopic dynamical information. On the other hand, we make no assumptions about the microscopic distribution corresponding to the inhomogeneous state and thus avoid the choice of projected state; however, we do assume that our perturbation generates a response that captures the long-time correlations. As a result, we naturally find Green-Kubo expressions in terms of simple steady-state correlation-functions.
Our derivation of Green-Kubo relations is valid for any system as long as the microscopic dynamics is Markovian and the steady-state is statistically translationally invariant. Relaxing the Markovian assumption may be possible, since non-Markovian dynamics can be made Markovian by introducing auxiliary variables [36]. Another important direction is to inhomogenous boundary-driven steady-states, when the environmental interactions can be modeled via Markovian stochastic processes.
Finally, near-equilibrium time-reversal symmetry implies that cross-transport coefficients are equal, a prediction known as Onsager reciprocity [23]. The microscopic expressions for transport coefficients valid farfrom-equilibrium derived here open the door to studying the violation of Onsager reciprocity as well as its connection to time-reversal-symmetry breaking and dissipation.

Appendix A: Equilibrium-like fluctuation-response relation
In this section, we derive and analyze the equilibriumlike fluctuation-response equality (4).
To proceed, we will solve to linear order the Master equation of the perturbed dynamics ∂ t P (z, t) = LP (z, t) − λ(t)L(Q(z)P (z, t)). (A1) We assume the unperturbed system (λ = 0) relaxes to a unique steady-state distribution P ss (z) given as the solution of LP ss (z) = 0. Taking the initial condition at t = −∞ as the steady-state distribution P (z, −∞) = P ss (z), the formal solution of (A1) is (A3) Here, we have identified the time-translationallyinvariant steady-state correlation function, which is de-fined for any two observables, O 1 (z) and O 2 (z), by in terms of the transition probability P (z, t|z , 0) obtained as the solution of (A1) with a delta-function initial condition δ(z − z ).
The defining characteristic of this perturbation is that the conjugate coordinate Q(z) appearing in the perturbation appears naturally in the response correlated with the observable. Since this attribute is central to our analysis, we next deduce the class of perturbations with this property.
To this end, consider the dynamics in the presence of a generic linear perturbation for an arbitrary linear operator M suitably well behaved. Following the identical linear perturbation theory analysis as above, we find that the response of an observable to this perturbation is MP ss (z) = L(Q(z)P ss (z)). (A8) In Sec. II, we demonstrated that our dynamical perturbation (3) is equivalent to perturbing the energy of a deterministic Hamiltonian dynamics. We can now use our equivalence condition (A8) to demonstrate that perturbing the potential energy of an equilibrium system generates the same response as our perturbation (3) for a broader class of nondeterministic dynamics.
To this end, we consider a general underdamped Langevin dynamics whose generator is given by (arguments are suppressed for clarity) where q and p are collective notations for the positions and momenta of all constituent particles, H = H(q, p) is a Hamiltonian-like energy function, f nc is a nonconservative force, G is a matrix describing friction, and T is a matrix describing thermal fluctuations. The underdamped Langevin dynamics (A9) encompasses (i) Hamiltonian dynamics as a special case where f nc = G = T = 0 and (ii) overdamped Langevin dynamics as a limiting case where friction and thermal fluctuation are strong. We define a potential energy perturbation as a change of the energy function H(z) → H(z) − λ(t)U (q), whose generator is given by L − λ(t)M with MP (z, t) = (∇ q U (q)) · (∇ p P (z, t)). (A10) Let us compare this with our perturbation (3) with the choice Q(z) = U (q): L(U (q)P ss (z, t)) = −P ss (z)(∇ p H(z)) · (∇ q U (q)), (A11) The equivalence condition (A8) then implies the potential energy perturbation is equivalent to L − λ(t)LU (q) if or equivalently, ln P ss (z) = −H(z) + h(q) for any arbitrary position-dependent function h(q). The condition (A12) is satisfied at equilibrium where the steadystate distribution is given by the Boltzmann distribution P ss (z) = e −H(z) /( dze −H(z) ) in units of β = 1. Therefore, the potential energy perturbation H(z) → H(z) − λ(t)U (q) is equivalent to L − λ(t)LU (q) at equilibrium.

Appendix B: Derivation of Green-Kubo relations
In this section, we elaborate the connection between the equilibrium-like fluctuation-response equality and Green-Kubo relations in the presence of multiple hydrodynamic variables.
Consider a system with n hydrodynamic variables, {Â α r (z)} r∈V with α = 1, · · · , n, locally defined at each location r. In a homogeneous steady state, the hydrodynamic variables have uniform mean valuesĀ α = A α r ss . Hydrodynamic transport coefficients describe the rate at which an inhomogeneity relaxes away as the system approaches the homogeneous steady state. In order to exploit (4) to derive Green-Kubo relations for hydrodynamic transport coefficients, we choose the form of the perturbation functions to be where is a positive constant which we will eventually take to zero so that that perturbation is turned on slowly and Θ(t) is the Heaviside step function forcing the perturbation off at t = 0. The conjugate fields f α r are arbitrary as long as all integrals converge. The Einstein summation convention is adapted for repeated Greek indices (1 ≤ α, β, γ ≤ n) throughout this section. Under these choices, equality (4) for t ≥ 0. We eliminate the time integral by first integrating by parts and then taking the limit → 0 (adiabatically turned-on perturbation), Using the translation invariance of a homogeneous steady state Â α r (t)Â β r (0) ss = Â α r−r (t)Â β 0 (0) ss and an exchange of integration range V dr V dr ≈ V d(r − r ) V dr for a large system size V twice, we obtain a compact symmetric expression for the Fourier modeŝ The choice of f β k is immaterial since it can be eliminated by the initial condition δÂ α k (0) = Â α k (0)Â β −k (0) ss f β k /V . Although equality (B4) is based on the particular ensemble generated at t = 0, the long-time behavior of the hydrodynamic modes is insensitive to this initial ensemble. At the same time, the long length-and time-scale relaxation behavior of hydrodynamic variables δA α k = δÂ α k towards the homogeneous steady state will be well captured by the linear regression equation with a generalized transport coefficient M k . The connection between microscopic and macroscopic dynamics can be made by noticing that particularities of the initial ensemble at t = 0 will die out fast and (B4) should become universal at long length-and time-scales. Therefore, comparing (B4) and (B5), we conclude that the equality holds at long times for small k = |k|.
In the following, we make use of the procedure to drive Green-Kubo relations for equilibrium transport established by Oppenheim and collaborators [33,34]. We begin by manipulating (B6) into a form that makes as many of the time-derivates explicit as possible, allowing us to do an expansion for long times and small k. To this end, we note that the long-time-scale dynamical information is contained in the small frequency modes of the half-Fourier transform δÂ kω = ∞ 0 dte iωtÂ k (t). Taking the half-Fourier transform in time, we can express (B6) in two ways. First, by directly taking the half-Fourier transform or, second, by using the transform to replace a time derivative with a frequency as For brevity, from now on, we will suppress the time argument t = 0 unless especially necessary for clarity, thuŝ from which we can derive Green-Kubo relations for hydrodynamic transport coefficients. It is useful to define a dissipative rateḂ α . From equality (B7), the correlation Ḃ α kωÂ γ −k ss vanishes in the hydrodynamic limit, suggesting that the fluctuations ofḂ α k decay on microscopic times. Thus, the dissipative rateḂ α k represents the rapidly fluctuating part ofȦ k . Since More generally, time-retardation effects in the macroscopic dynamics by allowing the generalized transport coefficient to frequency dependent M kω [4]. Repeating the same algebraic procedure in this case leads to the same end result (B11) except for the replacement of M αβ k by M αβ kω . For simplicity, we adhere to the frequencyindependent case.
The last step in the derivation of Green-Kubo relations is to connect hydrodynamic transport coefficients to the microscopic correlation functions using equality (B11). To this end, we assume that the generalized transport coefficient can be expanded as a series in ik in the hydrodynamic limit, which is equivalent to a gradient expansion in real space. The first two (multi-component) coefficients v and D of the expansion M αβ k ≈ −ik · v αβ + k · D αβ · k + · · · describe the drift and diffusive behavior of the macroscopic dynamics, respectively. Inserting this expansion of M k into (B11) and comparing the terms order by order in k, lead to the Green-Kubo relations. Although this procedure can be done generally without any assumption on the hydrodynamic variablesÂ α k , it can be simplified for local densities of conserved variables thanks to corresponding continuity equations. Thus, we derive Green-Kubo relations for local densities and Nambu-Goldstone modes separately.

Local densities
Local densities of conserved variables satisfy continuity equations ∂ tÂ α r = −∇ · j α r , which define the corresponding local currents j α r . The continuity equations in Fourier space,Ȧ α k = ik · j α k , make explicit that the rates of change of local densities are at least linear in k. The dissipative ratesḂ α k =Ȧ α r + M αβ kÂ β k ≈ ik · (j α k − v αβÂ β k ) + · · · are also proportional to k in the small-k limit, which suggests the definition of dissipative currents I α k = j α k − v αβÂ β k [33,34]. Plugging the definitions of currents into (B11), leads to Defining the static correlation function asχ βγ = lim k→0 Â β kÂ γ −k ss and comparing linear terms in k, we The dissipative-current correlation term k · I α kω I γ −k ss · k does not appear, since it is at least quadratic in k. By introducing a first-order correction matrix E αγ via we can also derive a Green-Kubo relation for D αβ by comparing second order terms,

Nambu-Goldstone modes
Nambu-Goldstone modes, which emerge when a continuous symmetry is broken, are not necessarily conserved. Thus the corresponding local currents are not defined in general due to the absence of continuity equations, and thus the k-dependence of the rateȦ α k is not simple. Moreover, steady states with a broken continuous symmetry are characterized by a long-range correlation, which is indicated by the divergence of the limit lim k→0 Â β kÂ γ −k ss . Due to the absence of local currents and finite static correlation functions, which made the derivations for local densities simpler, we need to consider the small-k expansion of (B11) more generally. Thus, we take a divergence of the form Â β kÂ γ −k ss ∼ k −q with an arbitrary exponent q. We recall (B10) here for conve-nience: To be consistent with macroscopic regression, the leading orders must balance in the small-k limit. Therefore, comparing the leading order terms of k after taking the following considerations, we obtain a Green-Kubo relation for v αβ : First, the term M αβ −k ss is clearly a higher-order correction. Second, in contrast to local densities for conserved quantities, there is no guarantee that Ȧ α kωȦ γ −k ss is higherorder than Ȧ α kωÂ γ −k ss . Lastly, we define inverse matrix withk = k/k. The next order correction leads to a Green-Kubo relation for D αβ , It is noticeable that (B17) and (B18) are reduced to (B14) and (B15) when the mode is conserved. In fact, (B17) and (B18) are generally valid for any hydrodynamic variables, since they are derived without any prior knowledge aboutÂ α k .
Appendix C: Active Brownian Particles

Linearized Dean equation
In this section, we derive a coupled linearized Dean equation for a system of active Brownian particles. For the sake of readability, time arguments are suppressed in this section unless strictly necessary. The resulting linear equations predict a Green-Kubo relation for the diffusive transport coefficient D for the case of soft interactions at a high density. We recall the equations of motion for active Brownian particles (21): where F i = −∇ ri j =i φ(|r i − r j |) with interaction potential φ(r); simulations are performed with the specific choice φ(r) = (K/2)(r − a) 2 Θ(a − r). The time evolution equation of the local particle number density ρ r = N j=1 δ(r − r j ), known as the Dean equation [55], is not closed. Instead, it is determined by an infinite hierarchy involving the local harmonic densities whose nth-orders are defined by and The particle density corresponds to the 0th order harmonic density ρ r = c r ) = N j=1 e(θ j )δ(r − r j ). However, it was shown out that the infinite hierarchy of equations can be truncated at the first order in the hydrodynamics limit [68]. In what follows, we shall ignore the higherorder harmonics.
The truncated coupled Dean equations are then given by where the noise fields are constructed as The nonlinear terms in (C4) and (C5) can be linearized by assuming the fluctuations around the homogeneous state are small. Since the homogeneous state does not exhibit any spatial discrete symmetry, the homogenous densities areρ = N/V andp = 0. The fluctuations around this homogenous state are then denoted by ρ r =ρ + δρ r and p r = δp r . Neglecting the second-order terms in the fluctuations δρ r δρ r and δp r δρ r , we can approximate the interaction as for X ∈ {ρ, p}. Applying (C9) to (C4) and (C5), and taking the Fourier transform in space, we obtain the following linearized equations: . This type of linearization of the Dean equation is known to be accurate for soft interactions at a high density [56].
Since all the steady-state covariances of the noise fields in (C6-C8) are constant, we can replace multiplicative noises in (C10) and (C11) with additive noises in the homogeneous steady-state. Consequently, in the homogeneous steady state, (C10) and (C11) become coupled linear stochastic equations with additive noises: The noise fields Λ Equation (C12) belongs to the class of the multivariate Ornstein-Uhlenbeck processes, whose covariance matrix is analytically solvable. We define the steady-state covariance matrix of δρ k and ik · p k as [36]. The two-time correlation function of δρ k is given by δρ k (t)δρ −k (0) ss = [e −tL · Σ] 11 . Thus we havẽ In conclusion, the linearized Dean equation predicts that the transport coefficient is given by Noteworthy is that δρ k is the same as ρ k as long as k is nonzero since the offset V dr e ik·rρ = N δ(k) only contributes at k = 0.

Numerical simulations
We measure the transport coefficient D in two independent numerical simulations. ABP simulations were performed using bespoke computer code implementing the Euler algorithm in Python v2.7.16. The length scale was set by a = 1, the force scale by setting µ = 1, and the time scale by setting D r = 10 and D t = 0. All simulations were performed at a density ofρ = 1. The integration time step are chosen to be ∆t = 0.001 for relaxation simulations and ∆t = 0.005 for steady-state simulations. The simulations are performed for 10 combinations of two interaction strengths K ∈ {0.5, 1} and five self-propulsion speeds v 0 ∈ {1, 2, 3, 4, 5}.

a. relaxation simulation
We first measure D directly by observing the regression of ρ k . Due to the noise in the data, we only use the smallest-k mode to measure D. Initially, we place N = 1600 particles uniformly in the region between x = L/4 and x = 3L/4 and y from 0 to L, with L = 40, and then simulate the particle dynamics for 10 5 time steps. The time series of the Fourier modes of the local particle number density is calculated from the position data. We only observe the x-directional mode, i.e., k = (2π/L, 0). For a given parameter set (K, v 0 ), we perform 10 independent simulations to get 10 time series of ρ k . Each mode ρ k decays exponentially with an exponent −1/τ . We first measure the decay exponent 1/τ for 10 samples, and then estimate D from τ /(2π/L) 2 . Taking the average and standard error on the mean over 10 independently measured D, we estimate the transport coefficient and its error.
b. steady-state simulation Second, we deduce the transport coefficient from the Green-Kubo relation Dχ = C + E wherẽ Since the particles do not prefer any particular direction, the drift coefficient vanishes v = 0. To this end, we measure the steady-state correlation functions ρ k ρ −k ss , j k ρ −k ss , and j k (t)j −k (0) ss numerically. The local current j k is calculated by taking the x-component of (C20) The system size is chosen to be L = 20 for the steadystate simulations. To estimate the limiting value of the correlation functions for k → 0, we use the smallest-k mode as an approximation instead of taking the extrapolation from several modes. This is because the extrapolation is not reliable due to the noise in data and the small system size. To calculate the steady-state correlation functions, we assume that the system is ergodic and replace the ensemble average with a time average, i.e., We run a single long simulation for 4 × 10 6 time steps and drop the first 8 × 10 5 steps to discard transient dynamics.
The details of estimating the three quantities in (C19) are as follows. First, we estimateχ as the correlation function ρ k ρ −k ss of the smallest-k mode. Second, we estimate E by the slope of Im{ j k ρ −k } as a linear function of k in the small-k regime. To measure the slope, we use the first two smallest-k modes. Lastly, assuming exponential decay of j k (t)j −k (0) ss = j k j −k ss e −t/τ , we estimate C from τ j 2 0 ss . The two-time correlation j 0 (t)j 0 (0) ss is calculated using the Wiener-Khinchin theorem [36], and τ is estimated from a least-square linear fitting on ln j 0 (t)j 0 (0) ss as a function of t. To analyze the error, we divide the time series of ρ k and j k into 10 blocks and repeat the estimation of D for each block of data. The error of D is estimated by the standard error on the mean over the 10 estimations.
Appendix D: Noisy Kuramoto model

Linear approximation
In this section, we derive a linear approximation to the noisy Kuramoto model near the synchronization state. Thanks to the linearity, correlation functions can be obtained analytically, and thus the transport coefficient of the Nambu-Goldstone mode can be deduced from a Green-Kubo relation.
We recall the stochastic equation that governs the time evolution of the i-th oscillator (25): where the sum extends over neighbors of i. The oscillators are evenly spaced on a two-dimensional square lattice and the position of the i-th particle is denoted by r i . The intrinsic frequency Ω i is sampled from a Gaussian distribution with meanΩ and variance σ 2 . In the synchronized state, the oscillators rotate coherently at a common frequencyΩ, so that each oscillator deviates slightly θ i (t) = θ 0 +Ωt + δθ i (t), where θ 0 is an arbitrary offset. Assuming the phase differences between neighboring oscillators are small, we approximate the interaction in (D1) as Then near the synchronized state, the time evolution of the Fourier modes δθ k = j δθ j e ik·rj is governed by the equatioṅ For long-wavelength modes, we can linearize (D3) in k to obtain the approximated linear equatioṅ with solution Utilizing this solution, we can obtain the steady-state equal-time correlation function from the long time limit which depends on the disorder. Since the intrinsic frequency of each oscillator is independent of the others, the frequency correlation function is given by leading to Other correlation functions involving time-derivatives can be obtained from the two-time correaltion function δθ k (t)δθ −k (0) ss = δθ k δθ −k ss e −tk 2 K = σ 2 N k 4 K 2 1 + KT σ 2 k 2 e −tk 2 K . (D9) Explicitly, we have δ θ k δθ −k ss = lim t→0 ∂ t δθ k (t)δθ −k (0) ss

Numerical simulations
We measure the transport coefficient D in two independent numerical simulations. Simulations of the noisy Kuramoto model were performed using bespoke computer code implementing the Heun algorithm [69] in Python v3.7.4. The integration time step are chosen to be ∆t = 0.001 for the relaxation simulations and ∆t = 0.01 for the steady-state simulations. The phases of oscillators are defined in the range of [−π,π], and the average phase is set to zero. This is equivalent to choosing a co-moving frame of reference. The simulations are performed for 8 combinations of four interaction strengths K ∈ {1, 2, 3, 4} and two noise strengths (T, σ) ∈ {(0.002, 0.005), (0.01, 0.01)}.

a. relaxation simulation
We first measure D directly by observing the regression of | θ k | for the first five modes with the smallest k values. To obtain the average θ k , 100 independent simulations are performed for a given parameter set (K, T, σ). The intrinsic frequencies Ω i and the noise trajectories ξ i (t) vary by sample. The initial phases are assigned to be π/2 for a leftmost tenth of the oscillators and 0 for the others. The system size is L x = 1000 and L y = 20 and simulations run for 10 5 steps. The time series of θ k is calculated from the phase data. Only the x-direction modes are observed so as to see the slowest regressions. That is, the wave vectors of the observed modes are k = (2π/L x )n with n ∈ {(1, 0), (2, 0), (3, 0), (4, 0), (5, 0)}. The averaged modes | θ k | decay exponentially with exponent −1/τ (k). We first measure 1/τ (k) for the five different k by least-square linear fittings of ln | θ k |. Then a quadratic fitting of 1/τ (k) provides an estimate for D. The error is estimated from the standard deviation of the quadratic fitting.
The relaxation from the inhomogeneous initial state to the synchronized homogeneous steady state is also observed from two other quantities. The order parameter  . 4(a) shows the order parameter of the system increases with time, indicating a relaxation to the synchronized steady-state. Also, the phase distribution P (θ, t) in FIG. 4(b) shows the unsynchronized phases (θ = π/2) gradually absorb into around the synchronized phases (θ = 0). b. steady-state simulation Second, we deduce the transport coefficient from the Green-Kubo relation (20). Since the linear theory predicts θ kωθ−k ss is much smaller than θ k θ −k ss by the factor k 2 , we estimate the transport coefficient by the ratio To this end, we measure the steady-state correlation functions θ k θ −k ss and θ k θ −k ss numerically. The system size is chosen by L x = L y = L = 20 for the steady-state simulations. The first 18 modes with the smallest k values are observed. That is, the wave vectors of the observed modes are k = (2π/L)n with |n| ∈ {1, √ 2, 2, √ 5, 3, √ 8, √ 10}. The degeneracy of the modes are {2, 2, 2, 4, 2, 2, 4}, respectively. To calculate the steady-state correlation functions, we assume that the system is ergodic and replace the ensemble average with a time average, i.e., O 1 (t)O 2 (s) ss = lim T →∞ 1 T T 0 dt O 1 (t + t )O 2 (s + t ) for any observables O 1 and O 2 . Since the linear theory predicts the relaxation time is τ (k) = 1/(k 2 K), we run a single long simulation of (5 × 10 6 + τ (k = 2π/L x )/∆t) time steps and drop the first τ (k = 2π/L x )/∆t steps to discard transient dynamics. We perform a weighted quadratic fitting on θ k θ −k ss /(k 2 θ k θ −k ss ) with respect to k so as to get its limiting value for k → 0. The weighting factor is the inverse square of the error. We estimate D and its error by the fitting value a of the fitting function a + bk 2 and its standard deviation, respectively.
To estimate errors of correlation functions, we found the effective independent number of data points n eff = t run /(2t corr (k)) by measuring the correlation times of θ k (t)θ −k (t) andθ k (t)θ −k (t) directly, where t run = 5×10 3 is the simulation time and t corr (k) is numerically measured correlation time. The errors are estimated by the standard deviations of θ k (t)θ −k (t) andθ k (t)θ −k (t) divided by the factor √ n eff [38].