Critical active dynamics is captured by a colored-noise driven field theory

Active matter may sometimes behave almost indistinguishably from equilibrium matter. This is particularly evident for some particle-based models and active field-theories close to a critical point which falls in the Ising universality class. Here we show however that, even when critical, active particles strongly violate the equilibrium fluctuation-dissipation in the high-wave-vector and high-frequency regime. Conversely, at larger spatiotemporal scales the theorem is progressively restored and the critical dynamics is in effective equilibrium. We develop a field-theoretical description of this scenario employing a space-time correlated noise field finding that the theory qualitatively captures the numerical results already at the Gaussian level. Moreover a dynamic renormalization group analysis shows that the correlated noise does not change the equilibrium critical exponents. Our results demonstrate that a correlated noise field is a fundamental ingredient to describe critical active matter at the coarse-grained level. Active particles can display a critical behaviour indistinguishable from an equilibrium one. The authors numerically study an active particles system close to the motility-induced critical point, and demonstrate that a nonequilibrium coloured noise field captures the coarse-grained behaviour of the system.

A ctive matter may display striking non-equilibrium phenomena such as the unidirectional propulsion of ratchet motors driven by active particles 1,2 or the spontaneous accumulation of self-propelled bacteria or colloids interacting with asymmetric obstacles 3 . However, there are situations where the non-equilibrium features of an active system are not immediately evident as, for instance, when self-propelled particles exhibit a collective motion or self-organization on large scales 4 similar to what observed in equilibrium systems. The archetypal example of this kind of behavior is the Motility-Induced Phase Separation (MIPS) 5 . This phenomenon can be microscopically understood in terms of active particles that move slower in denser regions, thus triggering an effective attraction that brings the system close to a spinodal decomposition [6][7][8][9][10][11][12][13][14] . Although the MIPS can be observed even in active particles without attractive interaction forces, it shares many similarities with the gas-liquid coexistence in equilibrium systems. Such analogies can be captured through effective equilibrium approaches which allow to reduce the non-equilibrium fluctuating forces to an effective interaction potential 5,[15][16][17][18][19] . The mapping to the gas-liquid phase transition suggests that MIPS ends in a critical point. Although previous numerical results have suggested that, in some active systems, such a critical point displays non-Ising exponents 20 , recent studies have shown that some other active on-lattice and off-lattice models have a critical point belonging to the Ising universality class 21,22 . As a consequence in the latter case, equilibrium-like φ 4 field theories 23 should provide the correct asymptotic description of MIPS around its critical point, at least at sufficiently large spatial dimensions. However it has been pointed out that additional terms, which break the time reversal symmetry, could be included in the standard φ 4 theory to capture the non-equilibrium and non-universal features near the MIPS critical point. Although these terms turn out to be irrelevant from the point of view of Renormalization Group (RG) transformations, they might yield, for example, a non-zero entropy production rate 24 . In this context it is thus crucial to understand if and how the critical dynamics of an active system becomes effectively identical to the one of an equilibrium critical system.
One of the most direct and natural way to unveil the nonequilibrium nature of a system is to look at the response and the correlation function of the observable of interest. These dynamical functions are related in equilibrium by the Fluctuation Dissipation Theorem (FDT). Its violation can be thus used to quantify how far the system is from equilibrium at various spatiotemporal scales. This approach has been widely employed to study the properties of several off-equilibrium systems such as glasses, gels and granulars [25][26][27][28][29][30][31][32][33][34] . Simultaneous measurements of response and correlation functions have also been used to reveal non-equilibrium fluctuations in active particle simulations [35][36][37][38][39][40] and in active experimental systems, such as living red-blood cell membranes 41 and suspensions of swimming bacteria probed by passive tracers 42,43 . Despite the deep interest on this topic, to our knowledge, no previous work has investigated systematically (over several spatio-temporal scales) the FDT violation in a microscopic model of active particles close to the critical point.
In the present work we collect and use the information from the critical correlation and the response of the order parameter to build a field-theoretical model that is able to faithfully reproduce the non-universal features of active particles close to the MIPS critical point. We find that the FDT is strongly violated at high frequencies and large wave-vectors. In this regime we find that the frequency-resolved correlator is much lower in amplitude than the corresponding response. Upon lowering the frequency and the wave-vector, we find that the response and correlator tend to coincide, thus satisfying FDT and validating the effective equilibrium picture. To rationalize this numerical evidence, we put forward a colored-noise-driven dynamical field theory, that is able to explain the coarse-grained behavior of the active system. We find that already at the Gaussian level, our theory predicts a scale-dependent violation of FDT that is restored on large lengthand time-scales in full agreement with the simulations. We further find that the additional parameters describing the noise field -its spatial correlation length and correlation time-are irrelevant under RG transformation and the system belongs to the Ising universality class. With respect to the active field theories considered so far, our model captures the non-equilibrium features of critical active particles without the addition of nonintegrable contributions to the field dynamics 24,44,45 . The overall picture, stemming from simulations and theory, unveils that a colored noise-field is a fundamental ingredient to describe near-critical active fluids.

Results
Simulations. We consider a system composed by N self-propelled active Ornstein-Uhlenbeck particles 16,46 (AOUPs) in d = 2. This model has been intensively studied analytically 38,[47][48][49][50][51] and, as discussed in the following, it is the ideal model for characterizing the perturbed and spontaneous critical dynamics. The equations of motion of AOUPs read where r i indicates the i-th particle's position and ψ i is the selfpropelling force. F i = ∑ j≠i f ij represents the force acting on the particle generated by two body interactions, i.e., f ij ¼ À∇ r i ϕðr ij Þ, with r ij = |r i − r j |. The two-body potential is modeled as a steep inverse power-law ϕ(r) = (r/σ) −12 /12 with a cut-off at r = 2.5 σ.
Here σ represents the diameter of the particle and is set to 1. In Eq.
(2)τ is the persistence time of the active force and ξ is a standard white noise source, i.e. hξ α i ðtÞi ¼ 0 and hξ α i ðtÞξ β j ðsÞi ¼ 2ðD=μ 2 Þ δ αβ δ ij δðt À sÞ, where the greek indices indicate the Cartesian components. Here D is the diffusivity of the noninteracting particles and μ is the particles mobility (set to 1 in simulations). Note that, by taking the limit τ → 0 at constant D, AOUPs reduce to passive Brownian particles in equilibrium at temperature T = D/μ (k B = 1 units are used throughout the present work).
We have recently located the MIPS critical point for this model and demonstrated that the resulting static critical exponents are in agreement with the Ising universality class 22 . In the present work we focus our study on the dynamics of the near-critical state-point with density ρ = 0.95 and active-force parameters τ = D = 16.5 and by employing various system sizes. Particles are enclosed in a rectangular box of sides (L x , L y = L x /3) and periodic boundary conditions apply, as in ref. 22 . Since the local density is the order parameter field of our critical system, we study the real part of its Fourier transform at time t, i.e. ρ q ðtÞ ¼ ∑ i cosðq Á r i ðtÞÞ, where the sum runs from i = 1 to i = N and q = 2π(n x /L x , n y /L y ) represents the wave-vector (here n x;y 2 Z). The quantity which statistically characterizes the spontaneous fluctuations of the field is the auto-correlation function where brackets indicate the average over configurations. Note that we assume that translational and rotational invariance hold in our system so that C(q, t) is identical to the correlator computed from the complex density field, i.e. hρ q ðt þ sÞρ Àq ðsÞi withρ q ðtÞ ¼ ∑ j e iqÁr j ðtÞ .
The linear response function, associated with C(q, t), is where h is the amplitude of an external force field f ext i which perturbs the i-th particle dynamics of Eq. (1) as i Þ and has the form f ext i ðtÞ ¼ À2 h ΘðtÞ q sinðq Á r i ðtÞÞ with Θ(t) being the Heaviside step function.
The definition of χ(q, t) immediately poses two main difficulties. The first is that computing χ(q, t) at different qvalues requires a different simulation for each q. Secondly, to ensure that the field amplitude is small enough to avoid nonlinear effects, one should repeat the simulations for various values of h. Fortunately, applying the field is not required for measuring χ(q, t) in AOUPs. Indeed, among the "family" of active models with exponentially correlated noise 52 , AOUPs have the unique feature that the exact linear response of any dynamic variable can be computed in unperturbed simulations as demonstrated in refs. 39,40 . In the present work we employ the method developed by Szamel 40 , which generalizes the Malliavin weights method 53 to systems driven by persistent Gaussian noise and efficiently combines it with parallel (GPU-based) simulations. The key idea of the method is to derive the appropriate fluctuating variables, encoding the linear response, by taking the limit of vanishing external field before taking the statistical average. Details on the method implementation can be found in Methods. Figure 1 shows the functions of interest C(q, t) and χ(q, t) at q = 2π(6/L x , 2/L y ) for a system of N = 3750 at (ρ = 0.95, τ = D = 16.5). The response function χ(q, t) is evaluated both with a small external field (h = 0.1) and with the field-free method, showing a good agreement between the two and validating the implementation of the latter technique.
We report here the results for correlation and response of the critical active system at different spatio-temporal scales. We focus on the frequency-resolved correlation and response functions, which we indicate as C(q, ω) and R(q, ω) respectively, that are obtained as the time-Fourier transforms C(q, ω) = ∫dt e iωt C(q, t) and R(q, ω) = i ω ∫dt e iωt χ(q, t). Note that using the factor i ω yields the frequency-resolved impulsive response function which, in equilibrium, is related to C(q, ω) by the FDT: ω Cðq; ωÞ ¼ 2 TR 00 ðq; ωÞ where the double prime indicates the imaginary part. Figure 2a, b reports the frequency-resolved correlation and response functions at the near-critical state-point (ρ = 0.95, τ = D = 16.5). To access a wide range of q-values we employ different system sizes ranging from N = 3750 to N = 60 × 10 3 , rather than simulating one single large system, in order to speed-up the simulations at high qvalues. Figure 2a shows the evolution of ω C(q, ω) for different qvalues. We observe that the correlator grows in amplitude and its peak shifts to lower frequencies upon lowering q, signaling a considerable slowing-down at large length-scales. Figure 2b shows that the response 2 T R″(q, ω) has a similar qualitative behavior. However, by reporting ω C(q, ω) on the same plot (gray dashed lines Fig. 2b), it is immediately evident that the response has a much higher amplitude than the correlator at large q-values, implying a dramatic violation of the FDT at these length-scales. Moreover, the high-q correlator decays faster than the response at high ω, further suggesting that the FDT-violation is also frequency-dependent. Differently, we observe that at low q-values ω C(q, ω) and 2 T R″(q, ω) become much closer and almost perfectly coincide for small frequencies.
The differences between ω C and R″ can be evidenced also by plotting these functions as color maps in the whole (q, ω)-range, as shown in the insets of Fig. 2a, b. From these plots is evident that ω C(q, ω) sensibly decreases in amplitude upon increasing q and ω while R″(q, ω) does not.
To appreciate the different frequency-dependence between correlation and response, we plot these functions scaled by their respective maximum in Fig. 2c, d. We find that, while the shape of ω C(q, ω) significantly changes upon changing q, R″(q, ω) does not vary and a good collapse is found for the scaled response at all q-s. This is also evident from the color maps shown in the insets of Fig. 2c, d. More specifically, in the low-ω regime, both correlation and response show a good collapse and they both grow as ω as highlighted by the dashed-dotted power law fits in Fig. 2c, d. Differently, in the high-ω region, the response decays slightly slower than ω −1 at all q-s, while the correlator behavior shows a progressive change from a single ω −3 decay at high q (purple curve in Fig. 2c) to a two-step decay at intermediate q (magenta curve), first showing a ω −1 dependence followed by a crossover to a faster decay. Finally, at very low q (orange curve), the decay is almost fully captured by a slow (approximately ω −1 ) power-law.
In addition, the peak positions of ω C(q, ω) and R″(q, ω), which can be identified as the system relaxation frequency, also show some differences. We denote the frequency where the peak is found as ω max and we plot it as a function of q for both functions in Fig. 3a. We observe that the correlator has an ω max lower than the one of the response in the high-q regime, implying that the response relaxes faster than the correlator at these length scales. However, at small q-s the two relaxation frequencies almost perfectly coincide, nicely following a power-law decay as ω max $ q z .
A direct power law fit of the low q-values (orange points in Fig. 3a corresponding to the largest system with N = 60 × 10 3 ) yields z = 3.78(0.13) for the response and z = 3.80(0.14) for the correlator (the fit error is reported in brackets), which are The field-free method allows to extract the correlation function of the spontaneous fluctuations C(q, t) (blue curve) and the integrated response function χ(q, t) (orange curve). The latter coincides with the response obtained by directly applying the field (green curve). The curves are evaluated at the near-critical state point τ = 16.5 and ρ = 0.95. compatible with the critical exponent z = 3.75 of the equilibrium Ising model with conserved magnetization in d = 2 54 . It is worth to stress that extending the fitting to larger q-s leads to a value of z that is lower than the Ising z, but also the quality of the power-law fit deteriorates since the data-points deviate from a straight line.
The deviation from the power-law behavior observed at large q-values signals that we are probing the microscopic relaxation and hence we are far from the scaling regime. To better understand this deviation, we have also simulated the critical equilibrium triangular lattice gas 22 , shown in the inset of Fig. 3a, which displays a qualitative similar deviation at large q-s. This suggests that, independently on the fact the the system is in equilibrium or not, at large enough q-values one is probing the dynamics at short scales where the microsocpic details play an important role and the relaxation frequency does not follow a scaling law.
To conclude this paragraph it is worth noting that, by using the active field theoretical framework proposed in refs. 24,44,45 , one could in principle evaluate the entropy production rate directly from the numerical correlation and response exploiting the Harada-Sasa (HS) relation 55 . However, our numerical results show that, while ω C(q, ω) decays relatively fast (ω C(q, ω)~ω −3 ) at high frequencies, the response R″(q, ω) decays much slower (R″(q, ω)~ω −1 ). This implies that the entropy production rate S cannot be computed using the HS relation for white noise, i.e. as done in ref. 44 by employing the formula S / R dq q dÀ1 R 1 À1 dω q À2 ω ½ω Cðq; ωÞ À 2 T R 00 ðq; ωÞ, because this quantity diverges upon performing the frequency integral (unless one introduces an upper frequency cut-off). Since this result is a consequence of assuming a white-noise-driven field theory, it suggests that our data can not be described by white noise only. Moreover, a similar divergence is obtained when studying one single harmonic oscillator driven by a colored (Ornstein-Uhlenbeck) noise, in which one also finds that R″(ω)~ω −1 and ω C(ω)~ω −3 . In the colored-noise case the HS relation has to be modified, as shown in ref. 36 , removing the divergence in S. These results suggest that (i) any field theory, including only white noise, is not sufficient to capture the behavior observed in numerical simulations and (ii) that our numerical results could be rationalized in terms of colored noise. We will return to this point in more detail at the end of the next paragraph after showing further evidenced in favor of a theory characterized by colored noise.
From the observations presented above it emerges that the FDT violation depends on both q and ω. It is therefore desirable to quantify the FDT violation by a single function, which Correlation and response in the critical active system. a Frequency-resolved correlation functions at various q-s (colored curves). Simulating larger systems gives access to lower wave-vectors (see legend). The inset shows ln½ωCðq; ωÞ as color map (value increasing from blue to yellow) at all q and ω-s. The gray area represents the highest and lowest ω measured at each q determined by the simulation duration and sampling frequency. The horizontal lines indicate the q-values selected for the main panel. b Frequency-resolved response function evaluated by the field-free method at the same q-s as in a (colored curves, same legend as in a). The correlation spectra for the highest and lowest wave-vectors in a are reproduced here as dashed lines to show that the deviation between response and correlation is more pronounced at high q. The inset is the same of the inset of a for ln½R 00 ðq; ωÞ. c Frequencyresolved correlation functions from a revealing an evident deviation from scaling in the high-q/high-ω regime. The asymptotic behaviors of these functions are highlighted by power-law fits (straight lines, see legend). The inset is the same of the inset of a for ln½ωC=ðωCÞ max . d Scaled response functions from b. Differently from c, a good data-collapse is observed for the response function. The asymptotic behaviors are evidenced by power law fits (straight lines, same legend as in c). The inset is the same of the inset of (a) for ln½R 00 =R we identify as the (normalized) frequency-dependent "effective temperature": The normalization ensures that, when T eff (q, ω) = 1, the system can be considered as being in effective equilibrium at the bath temperature T. We emphasize that here we do not assign any deeper meaning to the T eff than a convenient "violation factor". While it has been shown that, in some glassy systems, the separation of timescales of aging processes allows one to interpret T eff as a "thermodynamic" temperature 28,56 , here we rather use Eq. (3) only as a suitable measurement of the FDT violation.
To first understand the q-dependence of T eff , we show in Fig. 3b the quantity T eff ðq; ω max Þ that is the T eff evaluated at the characteristic frequency ω max of the correlator. We find that the large-q regime is dominated by a strong violation of the FDT with a T eff ðq; ω max Þ < 1. Upon lowering q, however, T eff ðq; ω max Þ asymptotically tends to unity suggesting that, over large scales, the system relaxes as if it were in equilibrium at temperature T. Moreover, to understand the ω-dependence of T eff , we plot T eff (q, ω) as a function of ω in Fig. 3c. We note that, even at low q, where the relaxation frequencies of response and correlator coincide, T eff (q, ω) still displays clear violation of the FDT at high ω, where the decay of T eff (q, ω) is well captured by a power law~ω −2 . This result is consistent with Eq. (3) considering that, at high frequencies, we have found the approximate decays ω C(q, ω)~ω −3 and R″(q, ω)~ω −1 (see Fig. 2c, d). Differently, for ω smaller than ω max , T eff (q, ω) remains close to unity and progressively approaches one as q decreases. The full (q, ω) dependence of T eff is depicted in the inset of Fig. 3c where we see T eff decreasing considerably at high (q, ω).
To summarize our numerical findings, we found that the spontaneous fluctuations (measured by the correlator) are weaker than those induced by the external field (measured by the response) at short time-and length-scales. These quantities allow us to investigate the microscopic model on a coarsegrained level in terms of the density field having its own nonequilibrium fluctuating dynamics. The scenario found for these density relaxation functions can thus be modeled by considering Fig. 3 Relaxation frequency and effective temperature of the critical active system. a Peak frequency ω max of the correlator (circles) and of the response (squares) as a function of q. Error-bars represent the fit error on the peak position. Different colors indicate different systems sizes (see legend). The straight full line is a power-law fit of the low-q correlator peak-frequencies (i.e. the yellow points for which N = 60 × 10 3 ), yielding a critical exponent z very close to the Ising one (z = 3.75), which is lower than the mean-field value (the best fit with z = 4 is indicated by the dashed line). The inset (same axes as the main panel) shows the peak frequency for the equilibrium triangular lattice gas. Deviations from the scaling ω max $ q 4 are also found at high q-values. b Effective temperature evaluated at the correlator peak-frequency T eff ðq; ωÞ max as a function of q (different colors indicate different system sizes with the same as in a. Error-bars represent the linear propagation of the error on the peak amplitude. The black curve is a fit of the low-q data by means of the expression found in the one-loop colored-noise-driven field theory (Eq. (11)). c Frequency resolved effective temperature T eff (q, ω) in the low-q regime (colored curves), different colors indicate different q-values (see legend). The dashed line represents the power-law ω −2 predicted by the theory for asymptotic behavior of T eff . The inset shows ln½T eff ðq; ωÞ as color map (value increasing from blue to orange) at all q and ω-s. The gray area represents the highest and lowest ω measured at each q determined by the simulation duration and sampling frequency. The dashed black lines indicate the of low q-values range selected for the main panel. d Data in c are scaled according to the theory. Collapsed data are well fitted by the Lorentzian predicted by the theory (full line). The inset is the same of the inset of c but for the scaled effective temperature appearing in the main panel.
COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-022-00830-5 ARTICLE COMMUNICATIONS PHYSICS | (2022) 5:55 | https://doi.org/10.1038/s42005-022-00830-5 | www.nature.com/commsphys a coarse-grained noise-field which does not "excite" enough the density modes at high q and ω and therefore is "colored". Differently the external field, inducing the response, can vary on an arbitrarily high frequency and wavelength equally exciting all available modes. As discussed in the following paragraph, this picture can be made rigorous by means of a field-theoretic approach and describes qualitatively well the simulation data.
Theory. We now show that the results of the previous paragraph can be rationalized using an appropriate non-equilibrium relaxation dynamics of a conserved scalar field φ(x, t) at position x = (x 1 ,.., x d ) in the d-dimensional space at time t. The dynamics of the field is governed by the equations: hηðx; tÞ ηðy; t 0 Þi ¼ À i∇ ð Þ 2 K ζ;T ðjx À yj; jt À t 0 jÞ ð5Þ where the parameter γ represents a macroscopic mobility (setting the time-scale of the relaxation dynamics) and H LG is the standard φ 4 Landau-Ginzburg Hamiltonian In Eq. (5) the kernel K ζ;T ðjx À yj; jt À sjÞ indicates that the noise η(x, t) is time and space translational invariant and correlated over a length-scale ζ and over a time-scale T . If we further assume that the noise relaxes exponentially both in space and in time then, in the Fourier domain, the noise kernel takes the Lorentzian form where κ compactly indicates 57 the Fourier-vector κ = (q, ω) and T is the noise strength. As usual we assume that r ≥ 0 (r vanishes at the Gaussian critical point) and u > 0. These type of models have been studied in the past (but only in the case of a non-conserved order parameter 58 ) as the simplest colorednoise-driven field-theories. Note that by taking the limit T ! 0 and ζ → 0 the model defined by Eqs. (4-7) reduces to the standard version of Model B which describes the relaxational dynamics of a conserved scalar field in equilibrium at a temperature T 57 . More details about the field theoretical framework in presence of the correlated noise field are given in Supplementary Note 1.
We start by showing that, already at the Gaussian level (u = 0), the model (4-7) predicts a scenario in qualitative agreement with the simulation data. To show this, we solve Eqs. (4-7) and find correlator and response (details are given in the Supplementary Note 2) as, where G 0 (q, ω) is the standard Gaussian response propagator and the subscript zero refers to the Gaussian theory. To understand the behavior of response and correlation functions given by Eqs. (9) and (10) and how it compares with the numerical results of the previous paragraph, we plot ω C 0 (q, ω) and 2T R 00 0 ðq; ωÞ in Fig. 4a and b, respectively. Since we are interested in the critical point we set r = 0. As we will see in the following the peak frequency of R 00 0 ðq; ωÞ and ω C 0 (q, ω) for low q is given by ω max ¼ γ q 4 , therefore we choose γ such that the ω max of the Gaussian colored field-theory matches the one found in particles simulations at the lowest q (this yields γ ≈ 10 2 ). The noise parameters are set to T ¼ 16:47 and ζ = 5.59 (this particular choice will be justified in the following).
In Fig. 4a, b we plot ω C 0 (q, ω) and 2T R 00 0 ðq; ωÞ at the same q, ω-values of the numerical simulations. The theoretical results for ω C 0 (q, ω) and 2T R 00 0 ðq; ωÞ show a striking similarity to the numerical results of Fig. 2a, b, i.e. the correlator and the response are considerably different at large q but they almost coincide at low q-values. As in simulations we observe some deviations between ω C 0 (q, ω) and 2T R 00 0 ðq; ωÞ also at low q if ω is high enough (see dashed lines in Fig. 4b). Furthermore, the shape of the scaled ω C 0 (q, ω), shown in Fig. 4c, displays an evolution that is very similar to the numerical one of Fig. 2c. The correlator shows a fast (~ω −3 ) decay at large q, while it shows a slower decay (~ω −1 ) at intermediate and low q in full analogy with the numerical results of Fig. 2c. Also, for the scaled theoretical response of Fig. 4d, the agreement with the numerical results ( Fig. 2d) is evident as a good data-collapse is obtained.
From Eq. (9) we readily find that ω max ¼ γ q 4 for the response function, as also shown in Fig. 4e (squares). Differently for ω C 0 (q, ω), the function ω max ðqÞ is more complicated but it can be approximated at low q as ω max % γ q 4 À 2γ 3 T 2 q 12 . This implies that relaxation frequency for the correlator is slightly lower than the one of the response as shown also in Fig. 4e (circles) at large q-values in analogy with Fig. 3a. Note, however, that the fieldtheory cannot capture the deviation from the perfect scaling ω max $ q z for the response function observed both in numerical simulations and also in the equilibrium lattice gas (see Fig. 3a and its inset). Indeed the response (9), for r = 0, represents a truly scale-free relaxation in which the microscopic details of the dynamics have been completely washed-out by the coarsegraining.
By using (9) and (10) in the definition (3) we find the effective temperature of the Gaussian model (see also the Supplementary Note 2): We first notice that the effective temperature (11) contains the essential information about the noise kernel of Eq. (7). This implies that, according to this model, the study of T eff allows us to characterize the noise-field. In analogy with the analysis of Fig. 3b we report T eff evaluated at the ω max of the correlator in Fig. 4f. We find that this has the same tendency, observed in simulations, to approach unity from below as q decreases. In the regime where ω max ¼ γq 4 , Eq. (11) implies that, at low q-values, we can approximate T eff ðq; ω max Þ % ð1 þ ζ 2 q 2 Þ À2 . We use this formula to fit the low-q data of Fig. 3b to estimate the parameter ζ. We find that the characteristic length-scale of the noise field ζ = 5.897(0.073). This value is close to the characteristic length of the velocity correlation studied in ref. 22 for the microscopic critical active system. This suggests that the coarse-grained noise field embodies the velocity correlation that develops at the particle level. We report the T eff dependence on ω at different q-s in Fig. 4g. Similarly to the numerical results in Fig. 3c, at low ω, T eff is approximately constant (its value approaches unity upon lowering q) while T eff decays as ω −2 at high ω. Given the strong analogies between the theory and the numerical model, we use Eq. (11) to fit the data Fig. 3c. In this way we estimate T ¼ 16:47ð0:25Þ and ζ = 5.592(0.028) (this justifies the choice of parameters mentioned at the beginning of this paragraph).
Although we do not attempt here a microscopic derivation of the parameter T , we note that this is very close to the relaxation time of the microscopic active force. Finally we note that Eq. (11) also implies that T eff data for different q and ω should collapse on the same curve ð1 þ T 2 ω 2 Þ À1 when we plot ð1 þ ζ 2 q 2 Þ 2 T eff ðq; ωÞ as a function of (1 + ζ 2 q 2 ) ω (as shown in Fig. 4h). Applying the same procedure to the data in Fig. 3c we find a reasonably good data collapse as shown in Fig. 3d in the low-q regime, while some deviation at very high q is evidenced by the color-map shown in the inset of Fig. 3d. It is worth to stress that this analysis sets some strict bounds on the noise-field. For example, the noise cannot have long-ranged correlations as these would likely cause deviations from the Ising exponents. Moreover, the noise must have also a finite correlation length (and not just a finite correlation time) to reproduce the results of the microscopic model. At this point, two questions naturally arise: (i) since the Gaussian theory is strictly valid only for d above the critical dimension d = d c (which turns out to be d c = 4 as in equilibrium) why does the T eff of Eq. (11) fit well the simulation data in d = 2 ? (ii) is the colorednoise-driven field theory consistent with the critical exponent of the Ising universality class in d < 4 ? To partially answer these questions we then study the field theory perturbatively in d = 4 − ϵ dimensions as detailed below.
To show that the T eff (q, ω) given by Eq. (11) also applies below d c , we calculate the correlator and the response to first order in perturbation theory (see Supplementary Note 2). By inserting these functions in the definition (3), we find the same exact result of Eq. (11). The key steps leading to this conclusion can be summarized in the following way. We start by writing the Dyson equation of the form Gðq; ωÞ ¼ G 0 ðq; ωÞ À1 þ Σðq; ωÞ Â Ã À1 with Σ(q, ω) being the self-energy 23  where the prime indicates the real part of R 0 The crucial observation here is that the response (12) depends on the correlated-noise parameters only via the term I , while the Lorentzian kernel of the noise factors out in C(q, ω) as in C 0 (q, ω). Therefore, T eff (q, ω) turns out to be the same as the one computed in the Gaussian theory. This suggests that the expression (11) can be used as a suitable approximation to model the effective temperature also below the upper critical dimension.
We now turn to the question regarding the universality of the colored-noise-driven field theory. A scaling analysis of the model readily shows that, under the Kadanoff transformation x → bx and t → b z t, we have for the couplings and u 0 F:P: ¼ ϵ 16π 2 BðΛ; ζ 0 ; T 0 Þ À1 =3. On the critical surface ζ 0 F:P: ¼ T 0 F:P: ¼ 0 and we obtain CðΛ; 0; 0Þ ¼ Λ 2 , B(Λ, 0, 0) = 1, i. e., the Wilson-Fisher fixed point 59 .
This study confirms the idea that the universality class of the colored-noise-driven field theory is the same as the Ising universality class which is compatible with the observation reported in the present work and in previous investigations 21,22 .

Discussion
In this work we have studied numerically and analytically the dynamical properties of an active system around its MIPS critical point. We have found that the FDT is strongly violated at short time and length scales and that effective equilibrium is progressively restored at large spatiotemporal scales.
It is interesting to qualitatively compare this scenario with the breakdown of the FDT found in other types of non-equilibrium systems. For example in glassy systems, the FDT violation is typically stronger at large spatiotemporal scales 60 , which can be interpreted as the rapid equilibration of the fast degrees of freedom at the bath temperature followed by a gradual reequilibration of the slow degrees of freedom (associated with the collective rearrangements) which "remember" another temperature. Contrarily, in our critical active system the stronger violation occurs at small length-and time-scales highlighting a very different type of non-equilibrium behavior.
We have found that the the frequency dependence of the correlator and the response in the particle simulations cannot be rationalized in terms of a field theory driven by white noise only. Differently we have shown that the scenario resulting from simulations is captured by a field-theory where the order-parameter dynamics is driven by a noise field correlated in time and space. The model directly derives from the numerical observation that the spontaneous fluctuations appear to be weaker (at high q and ω) than the ones induced by the external field. In this context the correlated noise field turns out to be a crucial ingredient to model the critical dynamics at the coarse-grained level, as the theoryalready in its Gaussian version-qualitatively reproduces most of the non-equilibrium features observed numerically. In particular the Gaussian model predicts a scale-dependent effective temperature T eff (q, ω) that tends to a constant for (q → 0, ω → 0) and therefore equilibrium is restored asymptotically. We remark that this is consistent with our previous numerical results 22 showing that the Ising critical exponents are observed for large system sizes. In addition our model also phenomenologically justifies the deviations from the Ising exponents found in ref. 20 since the spatial correlations of the noise are important on a small scale and may affect the static critical exponents. In summary the main advantage of our colored-noise model is that it captures both the FDT-violation at short scales and the Ising-like critical behavior at large scales observed in simulations of critical AOUPs.
Regarding the problem of computing the entropy production rate in critical active particles numerical results show that, for AOUPs, the standard HS relation cannot be applied. It would thus be interesting to try to generalize the HS formula to the case of colored-noise-driven fields, as done before with the HS-type equation proposed in ref. 36 for one single degree of freedom driven by colored-noise. From the theoretical side, by using a dynamic RG approach, we have found that the colored-noise-driven fieldtheory falls in the Ising universality class also below the upper critical dimension. This happens because the noise memory kernel introduces two operators that are RG irrelevant. It is worth noting that, in the present work, the analytical computation has been done using Wilson's RG scheme up to one loop. A detailed study of our field theoretical framework, at higher orders in perturbation theory, is in progress also to understand how the effective temperature may change its functional form beyond first order. As a further perspective one could try to measure correlation and response functions also in the critical active lattice models 21,61 and to check if the current scenario applies also in these cases. Although no field-free method has been yet developed for measuring the response in these systems, it could be that their numerical efficiency still allows one to obtain quickly the response function by directly applying the field. Moreover, it would be also interesting to explore the consequences of having a correlated noise field deep in the phase separation region of the MIPS and to understand, for example, what the role such a complex noise may have in the formation of interfaces at a coarse-grained level. It would also be important to investigate if the colored-noise-driven field theory can predict a micro-phase separation, for some parameter range, and then probe these regimes numerically to understand how the critical point may be destabilized as suggested in previous theoretical studies 45 .
Finally, it would be worth trying to derive the correlated-noisedriven field-theory directly from the microscopic dynamics by using, for example, the Mori-Zwanzig formalism 62 or by applying the Ito's lemma to the density field as done in ref. 63 for equilibrium systems.

Methods
Implementation of the Malliavin weights for AOUPs. To evaluate the response over unperturbed trajectories we follow ref. 40 according to which the response function of interest is given by where the Malliavin variables Q and P can be rewritten in terms of single-particle variables, i.e. Q = ∑ i Q i and P = ∑ i P i . The evolution of the Q i and P i is governed by the equations: where H i is the Hessian matrix H i ¼ À2 q q cosðq Á r i ðtÞÞ defined by the dyadic product q ⊗ q. To compute the response (14), Eqs. (15) and (16) must be integrated according to the Ito's rule from t ≥ 0 with initial conditions Q(t ≤ 0) = 0 and P(t ≤ 0) = 0. This method is particularly convenient first because the white noise ξ i in Eqs. (15) and (16) is the same as the one in the particle's dynamics (Eqs. (1) and (2)) and no additional random numbers have to be generated. Secondly Eqs. (15) and (16) can be efficiently evaluated in parallel for all particles and summed to compute Eq. (14) only when needed.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code availability
The simulation codes which have been use to produce the data of this study are available from the corresponding author upon reasonable request.
Received: 14 October 2021; Accepted: 7 February 2022; Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/.