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

We numerically investigate the correlation function, the response and the breakdown of the Fluctuation-Dissipation Theorem (FDT) in active particles close to the motility-induced critical point. We find a strong FDT violation in the short time and wavelength regime, where the response function has a larger amplitude than the fluctuation spectrum. Conversely, at larger spatiotemporal scales, the FDT is restored and the critical slowing-down is compatible with the Ising universality class. Building on these results, we develop a novel field-theoretical description employing a space-time correlated noise which qualitatively captures the numerical results already at the Gaussian level. By performing a one-loop renormalization group analysis we show that the correlated noise does not change the critical exponents with respect to the equilibrium. Our results demonstrate that a correlated noise field is a fundamental ingredient to capture the features of critical active matter at the coarse-grained level.


INTRODUCTION
Active 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, should 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 non-equilibrium 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 spatio-temporal 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].
The aim of the present work is to collect and use the information from the critical correlation and the response of the order parameter field, in a data-driven approach, 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. To carry out this program, we perform numerical investigations based on an efficient field-free method allowing us to compute correlation and response functions at criticality probing different length and time scales [40]. 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. In particular we show that our results can be qualitatively explained by a scalar theory with conserved order parameter (Model B) and driven by a random field that is correlated both in space and in time. Interestingly, we find that already at the Gaussian level (i.e. ignoring the effect of fluctuations), our theory predicts a scale-dependent violation of FDT that is restored on large length-and time-scales in full agreement with the simulations. By further investigating this field-theoretical framework, we observe 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. Additionally, by taking into account the non-linear coupling, we show that the critical exponents do not change up to first-order in Wilson's RG, suggesting that the theory is consistent with the expected asymptotic equilibrium-like criticality, also close to the upper critical dimension. With respect to the active field theories considered so far, our model captures the non-equilibrium features of critical active particles already by considering only linear terms and without the addition of non-integrable (and non-linear) 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 in the field-theoretical framework that is needed to capture the non-universal features of near-critical active fluids.

Microscopic model
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 for its remarkable analytic properties [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 self-propelling force. F i = j =i f ij represents the force acting on the particle generated by two body interactions, i.e., f ij = −∇ ri φ(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 Perturbed and spontaneous dynamics of active particles. The field-free method allows to extract the correlation function of the spontaneous fluctuations (blue curve) and the integrated response function (orange curve). The latter coincides with the response obtained by directly applying the field (green curve).
white noise source, i.e. ξ α i (t) = 0 and ξ α i (t)ξ β j (s) = 2(D/µ 2 ) δ ij δ(t − s), where the greek indices indicate the Cartesian components. Here D is the diffusivity of the non-interacting 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]. Being the density 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 ∈ Z). The quantity which statistically characterizes the spontaneous fluctuations of the field is the auto-correlation function C(q, t) = 2N −1 ρ q (t + s) ρ −q (s) where brackets indicate the average over configurations. The associated linear response function is where h is the amplitude of an external force field f ext i which perturbs the i-th particle dynamics of Eq. (1) aṡ r i = µ (ψ i + F i + f ext 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 q-values requires a different simulation for each q. Secondly, to ensure that the field amplitude is small enough to avoid non-linear 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 Ref.s [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 Appendix 1. Fig. 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.

Correlation and response functions
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 T R (q, ω), where the double prime indicates the imaginary part. Fig. 2(a)-(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 q-values. Fig. 2(a) shows the evolution of ω C(q, ω) for different q-values. We observe that the correlator grows in amplitude and its peak shifts to lower frequencies upon lowering q, signalling a considerable slowing-down at large length-scales.

FIG. 2.
Correlation and response in the critical active system. (a) Fluctuation spectrum as a function of frequency ω at various q-s (colored curves). Simulating larger systems gives access to lower wave-vectors (see legend). (b) Frequencyresolved 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. (c) Scaled fluctuation spectra 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). (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)). Fig. 2(b) 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. 2(b)), 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. To appreciate the different frequency-dependence between correlation and response, we plot these functions scaled by their respective maximum in Fig. 2(c) and (d). Interestingly 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. 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. 2(c) and (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. 2(c)) 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. 3(a). 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. 3(a)) 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) in good agreement with the critical exponent z = 3.75 of the equilibrium Ising model with conserved magnetization in d = 2 [54]. 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. 3(a)), 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 [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, i.e. as done in Ref. [44] by employing the formula , since this quantity diverges upon performing the frequency integral. Interestingly, 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 . 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 favour of a theory characterized by colored noise.

Effective temperature
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 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. 3(b) 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 which can be interpreted as if, on small length scales, the systems appears to fluctuate with a temperature lower than the bath temperature. 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. 3(c). 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. 2(c) and (d)). Differently, for ω smaller than ω max , T eff (q, ω) remains close to unity and progressively approaches one as q decreases.
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 . 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 legend as in (a)). 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 . (d) data in (c) are scaled according to the theory. Collapsed data are well fitted by the Lorentzian predicted by the theory (full line).
quantities allow us to investigate the microscopic model on a coarse-grained level in terms of the density field having its own non-equilibrium fluctuating dynamics. The scenario found for these density relaxation functions can thus be modeled by considering 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.

Colored-noise driven field 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: 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 (|x − y|, |t − s|) 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 colored-noise driven field-theories. Note that by taking the limit T → 0 and ζ → 0 the model defined by Eq.s (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 Appendix 2.

Gaussian theory
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 Eq.s (4)-(7) and find correlator and response (details are given in Appendix 3) 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 Eq.s (9) and (10) and how it compares with the numerical results of the previous paragraph, we plot ω C 0 (q, ω) and R 0 (q, ω) in Figs. 4(a) 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 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 (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. 4(a),(b) we plot ω C 0 (q, ω) and R 0 (q, ω) at the same q, ω-values of the numerical simulations. The theoretical results for ω C 0 (q, ω) and R 0 (q, ω) show a striking similarity to the numerical results of Figs. 2(a),(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 R 0 (q, ω) also at low q if ω is high enough (see dashed lines in Fig. 4(b)). Furthermore, the shape of the scaled ω C 0 (q, ω), shown in Fig. 4(c), displays an evolution that is remarkably similar to the numerical one of Fig. 2(c). 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. 2(c). Also, for the scaled theoretical response of Fig. 4(d), the agreement with the numerical results ( Fig. 2(d)) 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. 4(e) (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. 4(e) (circles) at large q-values in analogy with Fig. 3(a). Note, however, that the field-theory 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. 3(a) 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 coarse-graining. By using (9) and (10) in the definition (3) we find the effective temperature of the Gaussian model (see also Appendix 4): 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. 3(b) we report T eff evaluated at the ω max of the correlator in Fig. 4(f). 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. 3(b) to estimate the parameter ζ. We find that the characteristic length-scale of the noise field ζ = 5.897(0.073). Interestingly this value is close to the characteristic length of the velocity correlation studied in [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. 4(g).
Similarly to the numerical results in Fig. 3(c), 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. 3(c). 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. 4(h)). Applying the same procedure to the data in Fig. 3(c) we find a good data collapse as shown in Fig. 3(d).
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 colored-noise 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 in the next paragraph.

Results from perturbation theory
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 Appendix 4). 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]. Then the perturbation theory to one loop yields where the analytical expression of the integral I is provided in Appendix 4. The corresponding correlation function reads 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 r = r b 2 , u = u b 4−d , ζ = ζ b −1 , and T = T b −z with z = 4 (see Appendix 5). While it is evident from the scaling analysis that d c = 4, it is also clear that T and ζ are irrelevant operators in the RG sense at the Gaussian fixed-point. This also suggests that the non-Gaussian fixed point is the usual Wilson-Fisher fixed point. In order to check this, we employ Wilson's RG up to one loop (further details are provided in Appendix 6). With = 4 − d, the calculation, that does not require any small T and ζ approximation, brings to the non-Gaussian fixed point r F.P. = − C(Λ, ζ , T )/6 and u F.P. = 16π 2 B(Λ, ζ , T ) −1 /3. On the critical surface ζ F.P. = T 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].

CONCLUSIONS
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 nonequilibrium 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 re-equilibration 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. 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 theory -already 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 [20] since the spatial correlations of the noise are important on a small scale and may affect the static critical exponents.
Finally, by using a dynamic RG approach, we have found that the colored-noise driven field-theory 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 it would be interesting 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. Finally, it would be worth trying to derive the correlated-noise driven 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 [63] for equilibrium systems.  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 (1), Eq.s (2) and (3) 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 Eq.s (2) and (3) is the same as the one in the particle's dynamics (Eq.s (1) and (2)) and no additional random numbers have to be generated. Secondly Eq.s (2) and (3) can be efficiently evaluated in parallel for all particles and summed to compute Eq. (1) only when needed.

Critical dynamics with correlated noise
In order to consider both cases of non conserved ad conserved scalar order parameter, we start with the following relaxation non-equilibrium dynamics with a = 0, 2 for Model A and B, respectively. The noise field η(x, t) is a fast degree of freedom but still correlated over a (non-zero) length-scale ζ a time-scale T . Its dynamics can be represented through the following Ornstein-Uhlenbeck evolution [58] Tη( with ξ(x, t) = 0 , After standard manipulations [57,[64][65][66][67][68], we arrive to the Janssen-De Dominicis dynamical action A[φ, ϕ] = A 0 + A int that is withφ(x, t) indicating the response field. We introduce the Fourier transform of the fields The dynamical action in the Fourier space is where we are adopting the notation κ = (q, ω) with q the wave vector (q ∈ R d ) and ω ∈ R the frequency, dκ a the ultraviolet cutoff. In the following we consider the behavior of the system close to the critical point, we assume that r ≥ 0 and vanishes at the critical point and we set the strength of the noise to T = 1.
In order to compute the response function we couple consider also the case in which the ϕ is coupled to an external field h and thus the Hamiltonian becomes The response function R(x − x , t − t ) is defined as follows [57] R(x − x , t − t ) = δ ϕ(x, t) δh(x , t ) h=0 = γ ϕ(x, t)(i∇) aφ (x , t ) .
In Fourier space, because of the translational and rotational invariance, we get R(q, ω) = γq a G(q, ω) with G(q, ω) the response propagator of the theory.

Gaussian Model
Neglecting the non-linear interaction we obtain the Gaussian model whose generating functional can be written as where we have introduced a doublet of scalar field Φ ≡ (φ(κ), ϕ)(κ)) and a doublet of external sources J ≡ (Ĵ(κ), J(κ)).
From the scaling analysis we get As usual we haveχ + χ + d = 0 and thus z = 2 + a meaning that z = 2 for model A, and z = 4 in model B. The couplings r and u scale in the same way in both models Above d = 4 the only relevant parameter is r. Around and below d = 4, the coupling of the non-linear interaction becomes relevant.

Non-Gaussian fixed point
In order to provide a quantitative analysis of the stability of the equilibrium non-Gaussian fixed point with respect to fluctuations due to correlated noise, we employ the usual RG machinery. This program can be accomplished using recursion relations near d = 4 as in the case of static critical phenomena.
As usual, we consider a scale parameter b > 1 and the corresponding scale transformation R s b that is a change of scale of factor b. The second operation R i b that has to be applied to the diagrammatic expansions is the integration over internal wave vectors k in the domain b −1 Λ < |k| < Λ and internal frequency ω ranging from −∞ < ω < ∞. The two commutative operations define the RG transformation R = R s b R i b . The renormalized coupling constants r R and u R we are looking for are then related with vertex functions Γ 1,1 (q, ω) and Γ 1,3 (q, ω) in the following way [57] r R = lim and the corresponding RG equations shall be obtained by performing the scaling r R → r and u R → u . For computing the vertex functions Γ(q, ω) (a,b) we employ diagrammatic perturbation theory up to one loop. The corresponding expressions of Γ (1,1) (q, 0) and Γ (1,3) (q, 0) are Γ (1,1) (q, 0) = γq 2 r + q 2 + u 2 I B (r) where we have defined the following quantities