Fluctuation-Dissipation Relations in Active Matter Systems

We investigate the non-equilibrium character of self-propelled particles through the study of the linear response of the active Ornstein-Uhlenbeck particle (AOUP) model. We express the linear response in terms of correlations computed in the absence of perturbations, proposing a particularly compact and readable fluctuation-dissipation relation (FDR): such an expression explicitly separates equilibrium and non-equilibrium contributions due to self-propulsion. As a case study, we consider non-interacting AOUP confined in single-well and double-well potentials. In the former case, we also unveil the effect of dimensionality, studying one, two, and three-dimensional dynamics. We show that information about the distance from equilibrium can be deduced from the FDR, putting in evidence the roles of position and velocity variables in the non-equilibrium relaxation.


I. INTRODUCTION
The fluctuation-dissipation relations (FDR) played a pivotal role in the development of non-equilibrium statistical mechanics. Starting from the regime of weak departure from equilibrium, with the pioneering work of Einstein on mobility-diffusivity relation, through the central work of Onsager on the reciprocal relations [1], a unification of further theoretical results was obtained in the Kubo's linear response theory [2], which includes also Green-Kubo relations useful for transport coefficients in spatially extended systems. In Onsager's and Kubo's theories, all relations remain "simple", provided that some kind of time-reversal symmetry holds in the unperturbed state. As soon as this symmetry is lost, the simplicity of linear response and transport coefficients is no more guaranteed [3].
Time-reversal symmetry is broken in widespread natural phenomena: a constraint against equilibration is usually given by non-equilibrium boundary conditions, or in the presence of very slow relaxation dynamics, as in glassy systems. Hindrance to equilibration can also be caused by internal non-conservative forces, such as those in granular systems and in self-propelled particles. In all those cases response theory must be generalized, paying a price in simplicity. Several generalised relations have been proposed in the last decades: they still connect the perturbed system to the unperturbed one, but often require some detailed microscopic knowledge of the system. Some generalised relations may take a simple form in particular situations, e.g. an effective temperature may replace the thermostat temperature in a range of well-separated timescales for aging glassy systems, but this scenario is far from being general [4,5]. More frequently, one faces a situation where the equilibrium FDR is modified by additional contributions of a more complex nature.
The several approaches to FDR present in the literature belong to two main classes. Class A requires the knowledge of the stationary distribution, while class B only requires the knowledge of the microscopic dynamical rules (transition rates or Langevin equations). First examples of FDR of class A were obtained for chaotic systems and Brownian dynamics with non-conservative forces [6,7]. In those cases, the knowledge of the steady-state probability distribution (at least, in some approximations) still maintains a fundamental importance: in the absence of this information the generalized FDR remains, somehow, an implicit relation which can be useful to test ansatz on the steady-state properties [8]. More recent FDR of this kind express the non-equilibrium contributions in terms of stochastic entropy production [9,10]. Examples of FDR of class B are derived from the Malliavin weight sampling [11], or the Novikov theorem [12], which have been intensively employed in the context of glassy physics and represent a powerful tool to calculate the susceptibility and, thus, the so-called effective temperature of a system [4,13]. However, in this approach, the FDR involves correlations with the noise and the physical meaning of these terms remain often difficult to catch. Other schemes in this class include the description in terms of "frenetic" contributions [14,15], focussing on the role of time-symmetric fluctuations out of equilibrium (see also [16] for derivations of analogous FDR for discrete spin variables).
The generalization of the FDR to active matter represents a challenging issue. Systems of active particles are usually far from equilibrium and many of them move in the solvent through complex mechanisms involving chemical reactions or mechanical agents, such as cilia or flagella [17][18][19][20]. In the spirit of minimal modelling, these systems could be described through simple stochastic dynamics that resembles that of passive colloids except for the addition of a coarse-grained time-dependent force, often called self-propulsion or simply active force [21,22]. This force replaces the microscopic details of the system that, in general, are related to complicate internal mechanisms of energy transduction and involve an intrinsic source of strong deviation from thermodynamic equilibrium. Except for special cases, the steady-state properties of these non-equilibrium models are not known analytically or without approximations and, thus, FDR of class A remain in implicit forms [23,24]. Explicit expressions can be worked out in the limit of small persistence time τ , the first example obtained in [25]. On the other hand, the Malliavin weight sampling procedure has been applied to the case of active particle dynamics [26] and, in particular, employed to numerically calculate susceptibility and effective temperature of suspensions of active particles [27][28][29][30][31], even, in phase-separated configurations [32]. We remark that usually the concept of effective temperature is thermodynamically meaningful only in systems with well-separated time-scales [33]. This approach has been also employed to calculate the transport coefficients, such as the mobility, in combination with an approximate method valid at low-density values or small persistence regimes [34,35]. An attempt to generalise FDR to active systems has been recently presented in [36] and in [15,37]: in the latter case a relation involving a double-time derivative of correlators appears, making the meaning of the formula and its numerical (and experimental) application not immediate.
In this paper, we provide a clear example in the framework of active matter for which a generalized FDR of class B can be explicitly obtained in terms of simple steady-state correlations, involving functions of the observables (such as positions and velocities), which do not require a knowledge of the steady-state properties and do not involve any approximation (e.g. for small persistence regimes, etc.). Moreover, we propose a new quantity to measure the departure from equilibrium, which is defined in terms of the generalized FDR and takes into account the features of the system under time-reversal. We study the behavior of such a quantity for non-interacting self-propelled particles confined via different potentials, ranging from non-harmonic single-well to double-well potentials. In the former case, we also investigate the effect of the dimensionality on the non-equilibrium features.
The article is structured as follows: in Sec. II, we introduce model and notations, while in Sec. III we report the FDR relation obtained in our approach. Sec. IV is dedicated to the numerical study of the linear response, exploring both convex and non-convex potentials, in one, two, and three dimensional systems. Finally, in Sec. VI, we compare our results with two popular approximated results showing their failure, while the final section is dedicated to discussions and conclusions.

II. SELF-PROPELLED PARTICLES
We consider a well-known scheme to describe the behavior of self-propelled particles, the Active Ornstein-Uhlenbeck (AOUP) model [38][39][40][41][42][43][44][45] (also known as Gaussian Colored Noise (GCN)). The AOUP has been employed to reproduce the phenomenology of passive colloids immersed in a bath of active particles [46][47][48] but also -perhaps at a more approximate level -the dynamics of self-propelled particles themselves. Its connection with other popular models for self-propelled particles has been addressed by some authors [49,50]. For instance, the AOUP model can reproduce the accumulation near the boundaries of channels and obstacles [40,51], the non-equilibrium clustering or phaseseparation [25,52,53] typical of active matter and the spatial velocity correlation spontaneously observed in dense active systems [54,55]. According to the AOUP scheme, the position in d dimensions x of the self-propelled particle evolves with the following stochastic equation: where γ is the drag coefficient. In Eq. (1), we have neglected the thermal noise due to the solvent since the thermal diffusion is usually smaller than the effective diffusion due to the active force, in several experimental active systems [18]. The term F(t) = −∇U [x(t)] models the force due to an external confining potential, while the term f a (t) represents the self-propulsion of the particle i and is described by an Ornstein-Uhlenbeck process: where ξ is a white noise vector with zero average and unit variance. The parameter τ is the persistence time while D a is the diffusion coefficient due to the self-propulsion. Finally, the term h(t) is an external small perturbation used to probe the linear response properties of the system. In what follows, the perturbed dynamics will be denoted by the superscript h to be distinguished by the unperturbed one. The response function can be computed as the observed (normalised) variation of some observable when the perturbation is an impulse at time 0, i.e. h(t) = γδx(0)δ(t) (with δx(0) a vector with infinitesimal displacements representing the impulsive perturbation). In general, with the exception of singular measures, it is possible to derive a compact formula for the response function of one of the degrees of freedom x i of the system to the perturbation of a degree of freedom x j , R xixj , defined -in the steady-state -as where t > 0 and, in the last equality, we have introduced the functional derivative with respect to the perturbative force h j (0). In the following, for the one-dimensional case, we will replace R xx (t) with R x (t) for notational convenience.

III. GENERALIZED FLUCTUATION-DISSIPATION RELATION FOR SELF-PROPELLED PARTICLES
The model introduced so far reaches a non-equilibrium stationary state even in the absence of perturbation displaying a non-vanishing entropy production [25,[56][57][58], with a few special exceptions: in the absence of external potential, or in the presence of linear forces, the detailed balance is restored and the system recovers an apparent equilibrium state with vanishing entropy production [59]. This is an imperfection of the AOUP model [60]. Except for these special cases, the probability distribution, p(x, f a ), is unknown. An approximation can be obtained, valid in the proximity of equilibrium (τ = 0) [25,42,59,61]. Therefore, in general, this model of self-propelled particles cannot be treated within the FDR of class A, because the FDR remains implicit: where the average · is obtained considering the unperturbed system and p(x, f a ) is the stationary probability distribution in the absence of perturbations that depends on the whole set of particle positions and velocities. Eq. (4) can be used only perturbatively in powers of the persistence time since the expression of p(x, f a ) is only known perturbatively in powers of τ [23]. For small τ , the system is near the equilibrium, and, thus, the relation (4) is restricted to the same regimes. A similar approach has been employed in [25]. We propose an alternative approach, inspired by FDR of class B, to calculate the response function in terms of suitable steady-state correlations. Below, we report the final outcome of our technique applied to an AOUP confined by an external potential, while the detailed calculations are accurately described in Appendix A: where all the correlations have been evaluated in the steady-state and we have introduced the particle velocity v j =ẋ j . According to our notation, repeated indices are summed, U (s) = U (x(s)) and ∇ i = ∂ xi . Let us comment in detail on Eq. (5): first, we observe that it is symmetric under time reversal (i.e. when times t and 0 are swapped) and involves both position and velocity correlations. The first line of Eq. (5) has the same form of the equilibrium FDR holding for passive particles. Indeed, when the detailed balance condition holds (for τ → 0), the first and the second terms on the right-hand side of Eq. (5) coincide, while the second line disappears, in such a way that D a γR xixj = x i (t)∇ j U (0) . The second line provides two additional terms that are exquisitely non-equilibrium contributions to the response function. Remarkably, these non-equilibrium correlations involve the particle velocity, without containing direct simple dependence on the perturbed observable, namely the particle position (but of course position is involved through the potential). At variance with the well-known equilibrium scenario, in active systems, R xixj (t) is not only determined by a time correlation involving the position but is strongly affected by the correlations between the other variables involved in the dynamics, the velocity in this case. The harmonic case -which is a special case where the AOUP model satisfies the detailed balance even when τ > 0 -is treated in Appendix B.
Understanding what are the leading terms in Eq. (5) is a challenging issue that cannot be performed in general but requires a numerical analysis. In addition, since the detailed balance does not hold in systems of active particles, in general, the reversibility is broken by the presence of non-zero steady-state currents that produce non-vanishing entropy. The structure of Eq. (5) suggests a way to measure how much the system is far from equilibrium by quantifying the unbalance between the couples of terms, to see the effect of irreversibility on the response function. Therefore, we define: where all the correlations are calculated in the steady-state. D x ij (t) and D v ij (t) vanish at the initial time, t = 0 and for any t in any equilibrium configurations where the detailed balance holds. For notational convenience, we replace D x ij and D v ij with D x and D v in the one-dimensional case.

IV. NUMERICAL RESULTS
In this section, we numerically check the FDR (5), in two typical cases that cannot be analytically solved: A) the quartic potential, B) the double-well potential. On the one hand, A) represents the simpler convex case of study, except for the harmonic confinement where correlation functions and responses can be analytically computed. In the harmonic case, we remark that the expression (5) is consistent with previous results [23,62] that, in particular, predict an exponential time-decay of the response function with a typical time γ/k (see the Appendix B). In this simple case, the response function does not show any dependence on the parameters of the active force. On the other hand, B) is the simpler non-convex case able to reveal the interplay between the self-propulsion and the non-convex curvature of the potential. This feature has already manifested many dynamical anomalies leading to the emergence of regions with effective negative mobility in the large persistence regime [63] showing a bifurcation-like behavior with non-local stationary probability distribution [64]. In the quartic case, we also unveil the role of the dimensions comparing the results for one, two and three dimensional systems.
The numerical study is performed keeping fixed D a and the potential strength, to investigate the role of the persistence time, τ . Time is measured in units of a typical time t * , ruling the relaxation of the response function for the passive system (obtained for τ → 0), which is chosen as a reference case. In particular, t * = γ/k for the harmonic potential and t * = 1/ √ D a k for the quartic case. In all the cases, the system is perturbed via a small force along the x-axis, with amplitude δx = 10 −3 . This choice guarantees the linear regime of the response calculated using Eq. (3).
A. Quartic potential Fig. 1 shows the response function, R x , in the case of a quartic potential, U (x) = kx 4 /4 for several values of τ exploring both the small and the large persistence regime, i.e. the near and far equilibrium regimes, respectively. The figure reveals the agreement between R x (t) calculated from its definition, given by Eq. (3), and using the correlation functions appearing in Eq. (5) adapted to the one-dimensional case. This numerically confirms the extension of the FDR to AOUP active particles, even far from equilibrium. Specifically, in the small τ regime, the different curves of R x (t) collapse onto the same curve that corresponds to the response function of a passive Brownian particle at temperature D a γ. Instead, in the large τ regime, the larger is τ the slower is the relaxation of the response function, R x (t), that roughly displays two decay regimes. Indeed, it is known that AOUP in a single-well non-harmonic potential accumulate far from the potential minimum roughly near the points where the active force is balanced by the external one, in such a way that the spatial density shows two symmetric peaks [65]. Somehow, the system behaves as if an effective double-well confinement control the particle dynamics, an effect captured by many approximated results [65]. In the quartic case, the τ -dependence is a non-equilibrium consequence caused by the interplay between the non-linearity of the potential and the persistence of the active force. Indeed, in the harmonic case, the response is τ -independent and uniquely determined by the relaxation time of the potential [62], namely t * = k/γ.
In Fig. 1, panels (b)-(g), the terms of the FDR (Eq. (5)) are separately studied and compared with R x (t), for different values of τ = 10 −1 , 1, 10, respectively. As emerges from Fig. 1 , only weakly affect the first time decay of R x (t) also vanishing approximatively at t ∼ τ . When τ is increased, the contribution of the velocity correlation starts growing even if its decay remains roughly determined by τ , while the potential correlation decreases slowly. Interestingly, the velocity correlation reaches negative values and increases very slowly towards zero almost balancing the potential correlation, as shown in Fig. 1 , becomes dominant in the large persistence regime, as shown in Fig. 1 Finally, to understand the role of the detailed balance violation in the temporal decay of the response, we plot the terms D x (t) (solid red lines) and D v (t) (solid blue lines), in panels Fig. 1 (e)-(g). If these two terms are close to zero, then the detailed balance holds. As expected, this occurs for small τ while the detailed balance is strongly violated as far as τ is increased. Moreover, the irreversibility manifests in different ways: D x (t) is positive and reaches  a peak for t ∼ τ , while D v (t) shows an oscillation from positive to negative values. As a final remark, the violation of the detailed balance strongly suggests that any equilibrium-like approaches that assume this condition, such as the Unified Colored Noise approximation (that will be explicitly evaluated in Sec. V), cannot work in the large persistence regime.

Higher dimensional systems
We explore the role of the system dimensions, both on response functions and FDR. In the case of harmonic confinement (in equilibrium both for dimensions d = 1, 2, 3), where the detailed balance holds, the dimensions of the system do not affect the time-decay of the response function.
In Fig. 2 (a), we compare R x (t) for one-, two-and three-dimensional self-propelled particles confined through the quartic potentials, U (x) = k|x| 4 /4. We explore both the small and the large persistence regime reporting two reference cases, τ = 10 −3 , 10, respectively. In the small persistence regime (τ = 10 −3 ) that roughly coincides with the passive case, the larger is the dimension the faster is the decay of R x (t), an effect that is entirely due to the non-linearity of the potential. The increase of the system dimensions increases the effective trapping of the particle because of the confinement, an effect not related to the active force. Instead, in the large persistence regime (τ = 10), the correlations show a first temporal decay which does not depend on the dimension of the system (that coincides for d = 1, 2, 3), while, in the second stage of the decay, the higher dimensional systems decrease faster than their corresponding in lower dimensions. Somehow, the active force can suppress the dimensional dependence of the response function, at least for a small time-window. Fig. 2 (b) and (c) show the functions D x (t) and D v (t) for d = 1, 2, 3 at τ = 1, for simplicity. We highlight that results are similar for other values of τ (not shown). The qualitative functional forms of D x (t) and D v (t) are unchanged with d, resembling the shape described in the one-dimensional system: a single peak around t ∼ τ for D x (t) and an oscillation from positive to negative values for D v (t). The higher is the dimension of the system, the smaller is the maximal amplitude of both D x (t) and D v (t). This means that the increase of d, somehow, reduces the breaking of the detailed balance, producing smaller currents. This is also consistent with the increasing difficulties to observe collective phenomena when the dimensions are increased, that in three dimensions typically require larger values of the active forces compared to two-dimensional systems [66].

B. Double-well potential
The response function, in the double-well potential case, U (x) = k(x 4 /4−x 2 /2), is reported in Fig. 3 (a) for different values of τ , exploring both the small and the large persistence regimes. In the small persistence regime, R x (t) are collapsed onto the same curve for a broad range of τ (for τ ≤ 10 −2 ). As also occurs in the quartic potential case, the activity is the faster degree of freedom and can be approximated by a white noise, f a ≈ √ 2D a ξ, so that τ does not play any role. In this case, the relaxation towards zero is mainly determined by two time-regimes, as usual in the case of passive Brownian particles. The first is determined by the relaxation towards the minimum of one of the two wells, the second accounts for the jump from a well to another. When τ is increased, the first time-regime starts decreasing faster while the second time-regime decays more slowly approaching zero for very larger times t/t * τ . Moreover, the increase of τ also reduces the role played by the second time-regime in the temporal decay of R x (t) (the larger τ , the smaller the value of R x (t) when the second time regime starts occurring). The second time-regime is even suppressed when τ is sufficiently large (τ ≥ 5, with the potential setting employed in Fig. 3 (a)). This suppression is consistent with the result reported by Fily [67] in the infinite persistence regime where the steady-state density distribution exactly vanishes in the regions where the potential is concave (meaning that no jumps occur). More generally, the statistical relevance of the second time-regime decreases with the increase of the persistence time since the jumps from a minimum to the other occurs rarely when τ is increased both in the small τ regime [68][69][70] and in the large τ regime [64].
As occurs in the case of the quartic potential, the time-decay of R x (t) is dominated by the positional correlation or the velocity correlations, in the small and large persistence regimes, respectively (not shown). In Fig. 3 (b), (c), (d) and (e), we report D x (t) and D v (t) for different values of τ , revealing an interesting scenario. In analogy with the quartic potential case, the D x (t) assumes positive values via a single peak that shifts for larger t/t * when τ is increased. The D v (t) displays an oscillation around zero similarly to the quartic potential case. In a first range of τ , the amplitudes of both D x (t) and D v (t) increase with τ (as shown in panels (b) and (d)). A further increase of τ (panels (d) and (e)) produces the amplitude decrease of both D x (t) and D v (t), until both the functions become flat, approximatively for τ ≥ 10. This non-monotonic behavior with the persistence time implies that the system shows an optimal value of τ that maximizes the departure from the equilibrium via the breaking of the detailed balance. We justify this non-monotonic behavior through the phenomenology reported in [63]. When the particle moves close to one of the two minima, the particle explores a near-equilibrium regime where the local detailed balance almost holds and the local currents are almost zero both in the small and large persistence regimes. This explains why, in this space region, effective equilibrium approaches (that assume the detailed balance condition) give good predictions for the local steady-state properties. Instead, when the particle overcomes the inflected point of the potential (for which ∂ 2 x U = 0), the particle enters into a non-equilibrium region that, in the large persistence regime, displays an effective negative mobility. There, the particle velocity increases exponentially in time (with very small fluctuations) until the second minimum is reached through an accelerated motion. Thus, most of the non-equilibrium currents are generated when the particle jumps. As we already discussed, the number of jumps from a minimum to the other decreases with τ (in particular, until to be suppressed for τ = 10 and our setting of the potential). Our measure of D x and D v confirms that the detailed balance is mostly broken when the jumps occur and, thus, it is reasonable to assume that the detailed balance is almost restored for τ large enough where no jumps occur.

C. Measuring the non-equilibrium in active systems
The study of the correlations and, in particular, the time behavior of D x (t) and D v (t), suggests introducing a measure to quantify the breaking of the detailed balance and, thus, how much the system is far from equilibrium. Since D v (t) (and, in principle, also D x (t)) could assume negative values, we need to consider |D v (t)| and |D x (t)|. We propose a measure based on the time integral of these observables and introduce the quantity: where where t * is the typical time that rules the relaxation of the response function for the passive case (see above Sec. IV).
The observable µ has the following meaning: starting from an initial configuration (which is averaged), µ measures: i) how much the system is far from the equilibrium during the whole trajectory history; ii) the global impact of the detailed balance breaking on the time-decay of the response function. Indeed, when the detailed balance holds, both µ x and µ v vanish since D x (s) = D v (s) = 0. In general, the larger is µ, the greater is the aumount and how long the detailed balance has been broken starting from an initial configuration. These observables are shown in Fig. 4 (a) and (b) for the quartic and double-well potentials, respectively, for different values of τ /t * . In the quartic potential case, µ is an increasing function of τ meaning that the persistence of the activity increases the departure from the equilibrium. While µ ≈ µ x in the small persistence regime, this is no longer true in the large persistence regime where µ ≈ µ v . The double-well potential case confirms the non-monotonic behavior already observed in Fig. 3 (b)-(e): µ shows a peak and then starts decreasing for larger values of τ until to become almost zero. Moreover, µ is mainly determined by µ x for the whole set of τ values. Indeed, even if D v (t) assumes also larger value than D x (t), it is different from zero just for a narrow initial time-window while D x (t) remains larger for a longer time.

V. FAILURE OF THE APROXIMATED APPROACHES
To provide a qualitative explanation of several phenomena typical of active matter, such as particle accumulation near boundaries or the dynamics in confining potentials, different approximation schemes have been successfully introduced. Among the others, the Unified Colored Noise approximation (UCNA) has been crucial for theoretical purposes, providing the first analytical results for the probability distribution function of both interacting and noninteracting confined systems [71]. This approach provides a scheme to replace the self-propulsion with effective interactions and is quite similar to the so-called Fox approximation [72]. Moreover, it is derived assuming vanishing currents and, thus, the detailed balance, providing the best equilibrium-like predictions to describe the AOUP nonequilibrium dynamics.
Despite these approaches are very useful to understand the static properties of self-propelled particles, at least when confined through convex potentials or interacting through convex interactions, it has been shown that they fail to describe the active time-dependent properties, such as time-correlations and response functions. This concept has been already stressed in [23] where the authors show that the UCNA approximation is able to reproduce the response function just in small persistence regimes, where the system is near the equilibrium and a generalized FDR could be obtained perturbatively in τ using Eq. (4). In the large persistence regime, the numerical study for non-linear potentials shows the failure of the UCNA approach.
The aim of this Section is to enforce this idea, through additional analytical arguments, showing that the breakdown of the detailed balance plays an important role in correlations and response functions and, as a consequence, also for susceptibility and effective temperature.
A. Assuming the detailed balance A possible approximated approach could consist in assuming the detailed balance condition in Eq. (5). This choice means that D x (t) = 0 and D v (t) = 0 and should lead to results similar to those obtained in the UCNA approximation. Thus, the following relations hold: in such a way that the response can be approximated by where the superscript D stands for detailed balance. Using the property v applying the derivative and using the equation of motion, we obtain where we have also used the reversibility condition holding because of the detailed balance assumption. The results of this approximation are shown in Fig. 5 for the quartic potential and the double-well potential. As expected, in the case of a quartic potential the approximation holds for small τ = 10 −1 , i.e. in near equilibrium regimes, while is less accurate for large τ = 10. In the case of a double-well potential, the approximation holds for small and large values of τ , while is not accurate in the intermediate regime. We remark that the failure of R D x (t) to reproduce the decay of R x (t) is much evident in correspondence of the ranges of τ such that the measure of the non-equilibrium proposed in this work, Eq. (9), assumes large values.
We also stress that R D ij (t) is expressed in the following form: where C j (0) is defined by Eq.(14) that apparently looks like similar to Eq.(4). Moreover, we also point out that it is not possible to express C j as ∝ d/dx j ν(x, v) if we require that ν is a normalized probability as in Eq.(4). This observation alone implies that any approximation for the probability distribution could not agree with the detailed balance assumption.

B. Unified Colored Noise approximation
According to the UCNA scheme [65,71,73], we can explicitly derive an approximate solution for the steadystate probability distribution holding for general potentials that has been used to reproduce the accumulation near boundaries [71] (modeled as soft truncated potentials) and the particle accumulation far from the potential minimum of a single-well potential [65]. Moreover, as shown in [23], this procedure obtained neglecting the particle velocity leads to wrong results for time-correlations and response functions even in the harmonic case (that can be analytically checked). Recently, the UCNA has been successively extended to include the particle velocity [65,74], that displays a Gaussian-like shape with a space-dependent velocity variance. The extended version of the UCNA has the following form: where the effective Hamiltonian H and Γ(x) are given by: Here, the matrix Γ(x) plays the role of a space-dependent Stokes force that affects both the effective potential and the variance of the particle velocity. We remind that these interpretations should be limited to the convex potential (or non-convex potential in the small persistence regime) since the approximation cannot work if Γ(x) assumes nonpositive values in some regions of space. Thus, there is no hope of applying the method in the case of double-well potentials.
Using the UCNA, it is possible to derive an expression for the linear response function, R x (t). Indeed, plugging the log-derivative of Eq. (15) into Eq. (4), we obtain: where the superscript U stands for UCNA and K ij = v j (s)v k (s) is the generalized kinetic energy. The derivation of this result is reported in Appendix C. The results of the UCNA approximation are shown in Fig. 5(a) for the quartic potential. R U x (t) is a good approximation of R x (t) just in the small persistence regime (here, shown for τ = 10 −1 ), while is not accurate for large τ = 10 since R U x (t) decays slower than R x (t). We stress that R U x (t) in general does not coincide with R D x (t) because the last terms in the second lines of Eq. (18) and Eq. (14) are different. However, according to the UCNA prediction, τ Γ kj (x) K kj (x) ≈ D a , where K ij (x) is the space-dependent kinetic temperature introduced in Ref. [59,65] for the scalar case (and the last average is realized at fixed x). Eq. (18) coincides with Eq. (14) only by using this reasonable identification. We also stress that the UCNA dramatically fails to describe active particles in non-convex potentials because the matrix Γ −1 is ill-defined.

VI. CONCLUSION
Our study provides an explicit extension of the generalized FDR to non-equilibrium systems of active matter, that are numerically checked in many cases of interest, in particular, the case of the quartic potential explored for one, two and three dimensional systems and the case of a double-well potential. At variance with previous approaches, we derive simple expressions for the correlators involved in the FDR, which could be measured in experimental systems. Our results can be extended to more general dynamics, within the framework of non-equilibrium systems.
Through our expression, we are also able to extrapolate a measure to determine how much the system is far from equilibrium evaluating its influence on the time-decay of the response function. We quantify how much the breaking of the detailed balance affects the response decay looking directly at the correlations of the generalized FDR. This analysis agrees with previous qualitative observations showing a monotonic increase of the departure from the equilibrium with the increase of the persistence time in the quartic potential case but a reentrant (non-monotonic) behavior in the case of a double-well potential where the equilibrium is almost restored for large enough persistence times. Finally, as the intuition suggests, increasing the dimensionality reduces the departure from equilibrium.
Our expression of the FDR shows that the breaking of the detailed balance plays a crucial role to determine the time-decay of the response function and, thus, that any approximated descriptions that assume the detailed balance (or similar hypothesis to neglect the dynamical properties of active particles) to calculate responses, susceptibilities or mobilities are expected to yield poor results, unless in regimes of small activity, or when interactions are weak and confinement is absent or harmonic (all cases where detailed balance is weakly or not violated by this model).

VII. ACKNOWLEDGEMENTS
This research was funded by MIUR PRIN 2017 grant number 201798CZLJ and by Regione Lazio through the Grant "Progetti Gruppi di Ricerca" N. 85-2017-15257. The authors thank A. Vulpiani for useful discussions.
Appendix A: Derivation of the FDR for active particles As a first step, we manipulate the original dynamics (i.e. Eqs. (1) and (2)) to eliminate f a in favor ofẋ = v, through an exact change of variables. Applying the time-derivative to Eq. (1), using Eq.(2) to eliminateḟ a and, finally, replacing f a with v (through Eq. (1)), we get: where χ = √ 2D a ξ is a noise vector such that The term Γ(x) is a two-dimensional matrix acting as a space-dependent friction whose components read: where the space-dependence is provided by the potential curvature. The original dynamics, perturbed on the particle's position, is mapped onto an underdamped dynamics with an effective perturbation that is not simply h t as in Eqs. (1). The dynamics of v i is affected by an effective perturbation H t that reads: Despite the strangeness of the transformed dynamics, it is straightforward to apply the Novikon theorem [12] to get an exact expression for R xixj (t − s) in terms of noise correlations. Indeed, since the functional derivative with respect to H t /γ is equivalent to the functional derivative with respect to the noise χ: where we have used the chain rule in the last equality. The inversion of this relation reads: Thus, starting from the definition of R xixj (t − s), with t > s, we get: a function of time and the steady-state distribution is known. Since the system is linear, each component of the dynamics evolves independently with the others and, thus, studying the one-dimensional system is enough. For this reason, we drop the subscripts in what follows. Choosing U (x) = k/2x 2 , the steady-state distribution is a multivariate Gaussian, namely: where v =ẋ corresponds to the particle's velocity, such that: and the coefficient Γ is given by and represents the effective viscosity of the dynamics which is spatial independent in the harmonic case. We start evaluating Eq. (4), that explicitly leads to the following expression for the response function in terms of correlations: where, we d/dx is the total derivative (not just the partial derivative) and the function needs to evaluated as a function of x and f a because of Eq. (B2). We observe that R x reduces to the well-known relation for passive overdamped systems in the limit τ → 0, since in this limit Γ → 1. In this special case where the detailed balance holds, the correlations can be calculated as a function of time, in such a way that, in the steady-state, we get: Combining these results, we obtain which does not depend on τ . As a consequence, the response function in the presence of a harmonic force is not influenced by the value of τ . In this special case, also the generalized FDR, given by Eq. (5), can be evaluated explicitly and reads: Using the time-dependent expressions for the correlation functions, we obtain Eq. (B6), confirming that the equivalence between Eq. (4) and Eq. (5) in the harmonic case.
We remind that the spatial derivative comparing in Eq.(4) and, thus, in Eq.(C1), is a total derivative [23]. As already shown in [23], one needs to replace the velocity v with its whole expression as a function of the position before taking the derivative, i.e. γv = −∇U + f a , otherwise one gets inconsistent results. Following these procedures, we get: The spatial log-derivative of the Γ determinant can be evaluated in components and reads: Finally, collecting the results all together in Eq.(4), we obtain Eq. (18).