Anomalous Collective Dynamics of Auto-Chemotactic Populations

While the role of local interactions in nonequilibrium phase transitions is well studied, a fundamental understanding of the effects of long-range interactions is lacking. We study the critical dynamics of reproducing agents subject to autochemotactic interactions and limited resources. A renormalization group analysis reveals distinct scaling regimes for fast (attractive or repulsive) interactions; for slow signal transduction, the dynamics is dominated by a diffusive fixed point. Furthermore, we present a correction to the Keller-Segel nonlinearity emerging close to the extinction threshold and a novel nonlinear mechanism that stabilizes the continuous transition against the emergence of a characteristic length scale due to a chemotactic collapse.

While the role of local interactions in nonequilibrium phase transitions is well studied, a fundamental understanding of the effects of long-range interactions is lacking. We study the critical dynamics of reproducing agents subject to auto-chemotactic interactions and limited resources. A renormalization group analysis reveals distinct scaling regimes for fast (attractive or repulsive) interactions; for slow signal transduction the dynamics is dominated by a diffusive fixed point. Further, we present a correction to the Keller-Segel nonlinearity emerging close to the extinction threshold and a novel nonlinear mechanism that stabilizes the continuous transition against the emergence of a characteristic length scale due to a chemotactic collapse.
Nonequilibrium phase transitions encompass a broad class of systems, including absorbing-state phase transitions [1,2], roughening transitions [3,4], and ordering transitions in active matter [5,6]. Most theoretical studies of these paradigmatic model systems focus on the role of local interactions. However, in addition to shortranged interactions, several biological and synthetic systems exhibit many-body long-range interactions between agents [7]. For example, the social amoeba Dictyostelium discoideum uses chemical signaling and chemotaxis to control aggregation under harsh conditions [8], signaling molecules mediate intercellular communication in microbial populations [9], and microrobots and robotic fish use infrared, electrical, and acoustic signals to communicate [10].
Studying long-ranged interactions has a longstanding history in the context of equilibrium continuous phase transitions [11][12][13]. Their nonequilibrium counterparts are, however, less well explored. Most attention has been paid to systems where the long-rangedness results from Lévy-flight-like motion, nonlocal effects due to an underlying network architecture or spatially-dependent reaction rates [14][15][16][17]. There, the additional interactions may lead to a new universality class [14,16] or change the nature of the phase transition [17]. Here, we are interested in the role of long-range chemical signaling on classical models of population dynamics.
For this purpose, we consider agents emitting a signal in the form of a chemical substance which spreads by diffusion and can be sensed by other agents that respond by adapting their direction of motion, a process known as chemotaxis. The dynamics of such populations has been analyzed in terms of drift-diffusion models for the agent density coupled to a chemical field, termed Keller-Segel (KS) models [18][19][20]. These studies have identified a plethora of different phenomena -ranging from aggregation [21,22] to the formation of complex patterns [19,[23][24][25]. While the role of thermal fluctuations [26,27] and fluctuations around a constant back-ground density [28,29] have been investigated, the role of large-scale demographic noise -which is particularly important close to the extinction threshold [1,2,30] remains largely unexplored. In this letter, we investigate how long-ranged chemical signaling affects the collective behaviour of a population consisting of a single type of reproducing agents close to extinction.
We consider a generic model of a population of diffusing cells (agents) and chemicals in terms of two fluctuating density fields ρ(x, t) and c(x, t). The population dynamics is assumed to follow logistic growth, i.e., cells proliferate at a rate µ, die at a rate λ, and resource availability limits population growth to a finite carrying capacity. In addition, we consider the effect of an auto-chemotactic interaction, where each cell is capable of responding to a chemical signal, while simultaneously sourcing it with strength α [31,32]. We are interested in an effective, hydrodynamic description of this system, valid on macroscopic scales and in the presence of demographic noise. The corresponding Langevin equations are dρ dt = (D ρ ∇ 2 + θ) ρ − γ ρ 2 + 2Λρ ξ + I[ρ, ∇c] , (1) where θ = µ − λ is the net growth rate, θ/γ the carrying capacity, λ c the degradation rate of the signaling molecules, and D ρ,c are the diffusion constants. The macroscopically relevant noise is multiplicative with amplitude 2Λρ(x, t) and Gaussian white noise ξ(x, t).
Higher order nonlinearities and other noise terms are irrelevant close to the absorbing state [33]. Without the additional interaction I[ρ, ∇c], Eq. (1) corresponds to the noisy Fisher-Kolmogorov equation [34,35], whose universal properties fall into the universality class of directed percolation (DP) [1,2]. The interaction term I[ρ, ∇c] -which we assume to only depend on gradients in c [33] -accounts for the directed motion of cells along chemical gradients. Its form not only depends on cellular details but also on the level of coarse graining. In particular, the absence of global mass-conservation in the population dynamics allows for a nonconservative effective interaction.
At mean-field level, the dynamics exhibits two length scales, a diffusion length of the agents ξ ρ = D ρ /|θ| and of the chemicals ξ c = D c /λ c . The latter is linked to the interaction range, since λ c inhibits signal transduction over long distances. For long-ranged chemotactic interactions (ξ c → ∞, see Appendix A) [33], the only relevant scale is ξ ρ . Below this scale, demographic processes only play a minor role and the chemotactic interaction can be formulated in terms of a conserved current I[ρ, ∇c] = ∇J , where J = χ[ρ, ∇c]ρ∇c and the sensitivity function χ[ρ, ∇c] encodes details of the sensing process [19,36,37].
However, close to the extinction threshold the system is dominated by a divergent correlation length ξ > ξ ρ and strongly enhanced fluctuations. Further, coarse graining to large scales inevitably 'mixes' the effects of chemotaxis and birth-death processes. Whereas the net production by the linear birth-death term is independent of the density distribution, the net degradation due to the growth limiting term is enhanced by density fluctuations. Thus, the evolution of the total mass is coupled to the chemotactic interaction by the interplay of resource limitation and chemotactic drift, which alters the dynamics of density fluctuations. Therefore, an explicit coarse graining procedure is needed to determine all the relevant contributions. This is achieved by a renormalization group (RG) analysis (see [33]), which reveals that close to the extinction threshold the effective chemotactic interaction, correctly accounting for birth-death processes, is given by It consists of a conservative interaction -the classical Keller-Segel (KS) [36] nonlinearity -and an additional nonconservative term. A dimensional analysis shows that all other contributions are irrelevant at the pertinent length scales [33]. Importantly, Eq. (3) does not imply that the chemotactic interaction explicitly breaks particle number conservation. Rather it accounts for the fact that close to the extinction threshold the interplay between strong density fluctuations, chemotactic drift and population dynamics require an effective description of the form (3). Conversely, if fluctuation corrections are weak, i.e., far away from the extinction threshold, a conserved current yields the proper description.
To analyze Eqs. (1) and (2), we first neglect the noise term and study the resulting mean-field equations. They yield two homogeneous stationary solutions: the absorbing state ρ 0 = c 0 = 0 corresponding to the inactive phase and a state corresponding to the active phase with the agent density equal to the carrying capacity: ρ 1 = θ/γ and c 1 = αρ 1 /λ c . From a linear stability analysis of these homogeneous states one infers that there are three distinct phases (Fig. 1). For θ < 0, only the absorbing state is stable. In contrast, the homogeneous active state is stable for θ > 0 and In the case of θ > 0 and χ 2 below this threshold, however, both homogeneous solutions are unstable against spatial perturbations. This Turing-type [38] instability indicates the onset of pattern formation [24,25], as explicitly confirmed by numerical simulations shown in Fig. 1. At θ = 0 one finds a transcritical bifurcation, indicating a continuous, absorbing-state phase transition with θ acting as the control parameter. Close to the extinction threshold (θ → 0) and for long-ranged interactions (ξ c → ∞) the system becomes intrinsically scale invariant. In particular, the correlation length of density fluctuations should diverge as ξ ∝ θ −ν , and for a cell cluster emerging from a single seed, its mean-squared radius and survival probability at criticality should scale as ⟨R 2 ⟩(t) ∝ t 2/z and P (t) ∝ t −δ , respectively [1,2]. The mean-field critical exponents are given by ν = 0.5, z = 2 and δ = d/4. By dimensional analysis one identifies the following effective parameters: In addition to the DP coupling u (representing resource limitation) two new chemotactic couplings g 1 and g 2 emerge. The parameter w measures the time delay in the chemotactic interaction due to the finite diffusion speed of the signaling molecules. Employing field theoretical RG and a systematic perturbation expansion around the upper critical dimension d c = 4, we derive the flow equations [33] µ The flow functions f 1 -f 4 contain all information about the dependence of the theory on the arbitrary momentum scale µ in d = 4 − ε dimensions. Scale invariance is implied by the existence of IR-stable (µ → 0 stable) fixed points [30].
In contrast to previous studies [28], all calculations are performed by approaching the phase transition from the inactive phase, the full dynamics of the chemical concentration field are taken into account, and the limiting case of DP is correctly recovered.
Inspecting Eq. (6c), one observes that w = 1 is an invariant manifold of the RG flow. Moreover, systems where w ≲ 1 only slowly evolve away from this hyperplane. Therefore, we first focus on this quasi-static limit of infinitely fast diffusing chemicals [33].
We begin by investigating the case of a classical KS interaction. This implies starting the RG coarse graining at a scale where the chemotactic nonlinearities are equal, i.e., g 1 = g 2 = g 0 (gray plane in Fig. 2). In addition to the anticipated Gaussian and DP fixed points, the RG flow exhibits a stable fixed point (CA) and a stable fixed line (CR) (Fig. 2). They represent two different types of scale-invariant dynamics, corresponding to chemo-attractive (CA) and chemo-repellent (CR) systems. Only if g 0 = 0 the flow reaches the DP fixed point, which is unstable under the inclusion of chemotaxis, highlighting the importance of long-ranged interaction for the agents' critical behaviour. Further, irrespective of the sign of the interaction, the flow leaves the plane of KS interactions and terminates in either the stable subdiffusive CA fixed point (z = 2 + ϵ/18) for chemoattraction (g 0 < 0) or the stable superdiffusive CR fixed line (z = 2 − ϵ/2) for chemo-repulsion (g 0 > 0).
We conclude that accounting for long-range chemotactic interactions quantitatively changes the nature of the phase transition compared to DP, leading to two new universality classes of absorbing-state phase transitions. The values of the associated dynamical exponents z (Tab. I) match the expectation that chemorepellent agents accelerate and chemo-attractant agents decelerate colony dispersal compared to DP.
Further, the fact that all flow lines leave the g 1 = g 2 plane confirms that a KS interaction is not sufficient to model the universal dynamics near criticality. Fluctuation-generated terms are a generic phenomenon close to critical points [39,51]. Similarly, in our case the nonconservative part of Eq. (3) is 'generated' even if not included from the beginning and the effective chemotactic interaction can in general not be given in terms of a conserved current. Consequently, close to criticality g 1 ̸ = g 2 is of great physical interest. In particular, the question arises how the RG analysis relates to the meanfield analysis, which identified a band of linearly unstable modes for g 2 < −u (in the long ranged limit).
Indeed, the RG flow equations can be rewritten as a set of only two equations forū = u + g 2 and g 1 : In the quasi-static limit, the solution of the resulting Poisson equation (see Appendix B) allows to eliminate the chemical field, leading (among other terms) to an effective growth-limiting term with the shifted coupling constantū = u + g 2 [33].
We find that the domain of attraction of the CA and CR fixed points are separated by an invariant manifold at g 1 = 0 ( Fig. 3(a)), leading to two different types of dynamical scaling behaviors for g 1 < 0 and g 1 > 0, respectively. This further stresses the difference between the two chemotactic couplings: While the term ∼ g 2 ρ∇ 2 c can be absorbed into an effective growth-limiting term, only the nonlinearity ∼ g 1 ∇ρ∇c qualitatively changes the RG flow. In addition to the separatrix at g 1 = 0, the RG flow is organized by the critical manifolds containing the CA and CR fixed points, given (to one-loop order) by the linesū = g 1 /6 andū = g 1 , respectively ( Fig. 3(a)). These lines are also the boundaries of the domains of attraction of the CA (orange) and CR (blue) fixed points, separating them from runaway flow.
Given that, in the long-ranged limit, the instability condition (4) simplifies toū < 0, one might have antici-  pated runaway flow in this entire region. Strikingly, the RG analysis predicts scaling for g 1 <ū < 0, which seems contradictory at first. However, the linear stability analysis does not allow any conclusions about the steady state of the dynamics. Crucially, g 1 does not affect the linear dynamics, but only contributes to nonlinear effects. In particular, it enters the following exact relation for the time evolution of the average mass (see Appendix C) whereρ indicates a spatial and ⟨·⟩ an ensemble average with respect to the noise ξ. Equation (7) implies that, depending on the sign of g 1 −ū, fluctuations drive the system either toward or away from the absorbing state. It applies to the dynamics both above and below the absorbing-state phase transition, and especially when approaching the phase transition at θ,ρ → 0. This rationalizes why forū − g 1 > 0 (including all KS models) nonlinear effects combined with demographic noise lead to a continuous absorbing-state phase transition, despite the band of linear unstable modes forū < 0. In contrast, for g 1 > 0, the system is attracted by the CR fixed point forū > g 1 /6, and exhibits runaway flow whenū < g 1 /6 ( Fig. 3(a)). The region 0 <ū < g 1 is particularly interesting: Eq. (7) implies that the linear stability of the spatially uniform, active state is counteracted by a nonlinear term (∼ g 1 −ū) disfavoring a homogeneous state. Our RG analysis indicates that the antagonism between these two effects leads to flow towards the CR fixed point in the regime g 1 >ū > g 1 /6 but to runaway flow forū < g 1 /6. Since the nonlinear instability is dominant in the latter regime, the observed runaway flow is possibly indicative of a fluctuation-driven first-order transition.
The agents' active motion can result in an effective diffusion constant D ρ of similar magnitude as D c [31,40,41]. Therefore, it is crucial to study the case w ̸ ≈ 1. In this case, the full flow equations (6a)-(6c) exhibit an additional fixed point of mixed stability we call critical fixed point (CP) at (u, g 1 , g 2 , w) = (0.08ε, −0.45ε, −0.16ε, 0.64) and a second DP fixed point at w = 0. Our RG analysis shows that the CP fixed point has a dynamic critical exponent z = 2 to all loop orders [33], implying purely diffusive dynamics, akin to the critical fixed point characterizing the roughening transition of the Kardar-Parisi-Zhang equation [4,42,43]. As before, we first consider the case of g 1 = g 2 = g 0 ; the resulting basins of attraction for the various fixed points are depicted in Fig. 3(b). All points located on the invariant manifold g 0 = 0 flow towards the second DP fixed point at w = 0. Since the CP fixed point is located at w < 1 and unstable in w-direction, it separates the parameter space g 0 < 0 into two parts. Points above this separatrix flow to CA, whereas below it the system exhibits a new type of runaway flow (striped dark gray). In contrast to CA, the basin of attraction of CR does not extend to w < 1. As pointed out above, a chemo-repellent implies superdiffusive motion (z < 2), equivalent to ∂ µ w < 0 near the fixed line (CR). This renders the fixed line unstable in w-direction. However, in the emerging runaway region (striped blue) the projection of the fixed line to w < 1 is still a strong attractor which separates it from other regions of runaway flow ( Fig. 3(b) and (c)). The typical shape of the phase diagram for general g 1 and g 2 , at fixed values of u and w, is shown in Fig. 3(c). It features all four, possibly different, kinds of runaway flow and bears a strong resemblance to Fig. 3(a).
Altogether, the analyzed model reveals a correction to the well known Keller-Segel nonlinearity in the presence of large fluctuations and exhibits a rich phase diagram with two new absorbing-state phase transitions and various types of runaway regions. The emergence of fixed points associated to either a chemoattractant or -repellent, demonstrates the relevance of auto-chemotactic interactions for the collective behavior of cells at their extinction threshold. In particular, they highlight the impact of chemotactic signaling for the survival probability and spreading velocity of single colonies (Tab. I). For w = 1 we have presented a possible mechanism by which the runaway flow found in Fig. 3(a) can be related to a fluctuation-induced first-order transition (cf. (7)).
Further, the emergence of the CP fixed point not only gives rise to an unexpected type of purely diffusive scaling behavior, it also highlights the importance of the time delay introduced by the finite diffusion speed of the signaling substance. The reminiscence of the CP fixed point to the critical fixed point describing the roughening transition of the KPZ equation suggests the intriguing scenario of a strong coupling fixed point below the separatrix.
Naturally, the multitude of theoretical predictions presented calls for a numerical study. Additionally, we hope that our work will stimulate nonperturbative approaches [44,45] that help to unravel the observed anomalous dynamics. From a broader perspective, our results suggest that by combining known universality classes of nonequilibrium population dynamics [2,30] with various types of auto-chemotactic feedbacks, a broad class of novel scale-invariant dynamics could be discovered. Since scale invariance can only be observed if no length scale is introduced by the chemotactic interaction, the long-ranged limit λ c → 0 is of particular interest. However, simply inserting λ c = 0 into Eqs. (1) and (2) leads to a divergent chemical density and a steady state condition ρ 1 = −D c ∇ 2 c 1 (x)/α. Thus, there would no longer be a homogeneous steady state for the chemical density, which leads to an unphysical shift to the homogeneous steady state density ρ 1 = θ/(γ + αχ 2 /D c ) of the agents. This deviates from the actual carrying capacity θ/γ and shows that one needs to take into account the 'charge- where we subtracted the homogeneous, albeit time dependent average production of the signaling molecule withρ(t) denoting the spatial average of ρ at time t. Importantly, this homogeneous shift does not alter the dynamics of ρ. However, the evolution of the chargeneutral chemical density is now given by with no overall net production, i.e., More details on this limit are provided in the supplemental material [33].

APPENDIX B: QUASI-STATIC LIMIT
Another important limit is the so-called quasi-static limit, where D c /D ρ → ∞ and the chemical field thus instantly adjusts to changes in the density field ρ. Assuming that α/D c remains finite [21], Eq. (9) leads to the Poisson equation For more details we refer to the supplemental material [33].

APPENDIX C: MASS EVOLUTION
One way to analyze the impact of different interactions is to study their effect on the time evolution of the average density ⟨ρ⟩, where ⟨·⟩ signifies an ensemble average with respect to the noise ξ. To derive this evolution we first note that Eqs. (1) and (2) are Itô Langevin equations and thus Further we split the agents' density into ρ(x, t) =ρ(t) + ρ(x, t), integrate Eq. (1) over space and insert Eq. (11).
For the deterministic terms, this yelds where we used that Vρ = 0 and used the definitions (5) of the effective couplings, as well asū = u + g 2 . Taking the ensemble average of Eq. (13) leads to the exact result of Eq. (7). This result highlights the difference between the linear growth term ∝ θ and the nonlinearity ∝ γ modelling resource limitation. While the former contributes a distribution independent term to Eq. (13), the latter leads to a mass evolution which is dependent on the density profile. Thus, it is the resource limitation which makes the mass evolution susceptible to the influence of chemotaxis.

DERIVATION OF LANGEVIN EQUATIONS
The model analyzed in the main text consists of two parts: Diffusive particles A that obey logistic growth dynamics and a chemical which is secreted by A-particles and whose gradients influence the motion of A-particles. In order to derive a set of effective equations, we first treat the dynamics of demographic and chemotactic processes separately. To this end one may think of the following set of microscopic reactions A coarse-grained stochastic description in terms of the continuous density ρ(t) can be derived by a Kramers-Moyal expansion [46] of the corresponding master equation. This yields where θ = µ − λ denotes the effective growth rate, Λ = µ + λ − γ the noise amplitude and ξ Gaussian white noise. Equivalently, one may apply operator based approaches [47,48] and a subsequent Cole-Hopf transformation of the resulting field theory. Note that while different approaches strictly speaking correspond to different realizations of the stochastic process, all rely on Itô calculus and the underlying master equation. Beyond the continuous limit, the only approximation involved in all of these approaches is the truncation at second order in fluctuations which enables the description in terms of a Langevin/Fokker-Planck equation. However, higher orders can be shown to be irrelevant close to the absorbing state (see section 'The Response Functional'). While the above equation implies a well-mixed system, any description of chemotaxis requires a spatially extended description. Since chemotaxis implies that agents (A) adjust their motion to their surroundings, some form of active swimming is required. Even though the microscopic details of the biological processes leading to chemotaxis may vary significantly in different settings -which are of no particular interest to the present study -one may derive an effective descriptions for the chemotactic interaction. One way such an effective interaction can be formulated is to assume a chemotactic drift whose local velocity is given by the product of the local gradient in the chemical density c(x, t) and the sensitivity function χ[ρ, ∇c]. This results in a generalized form of the stochastic Keller-Segel (KS) model [26] (which in its classical form assumes a constant sensitivity) with an effective diffusion constant D ρ and Gaussian white noise η 1 . A more rigorous derivation of Eq. (16) can be given in terms of the Dean-Kawasaki approach [49,50] or a lattice gas with modified hopping rates. Further, the dynamics for the chemical density c(x, t) can straightforwardly be deduced from the microscopic reactions -secretion by A and decay, with rates α and λ c , respectively -and is given by with diffusion constant D c and η 2,3 again representing Gaussian white noises. To continue, one has to combine equations (15)- (17) and, even though the coarse graining procedures used are different and incorporate distinct effects, the first guess is to simply combine all the appearing terms; this yields It is important to emphasize that this is a highly coarse-grained description: While some parameters -like the growth rate θ -have a clear interpretation in terms of a microscopic model, others -especially the sensitivity function χ[ρ, ∇c] -have no such interpretation and are purely phenomenological. Therefore, Eqs. (18) and (19) can at most be valid at a finite range of scales. By inspecting Eqs. (18) and (19) one can already identify two important length scales: the diffusion length of the agents ξ ρ = D ρ /|θ| and the chemicals ξ c = D c /λ c . On length scales below ξ ρ , the population dynamics of particles A only play a minor role and can be neglected. However, since the goal of this manuscript is to extract the critical behavior of the presented model at the largest length (and time) scales and to evolve the dynamics to these scales, one has to be aware that the effective equations of motion might change as one continuously changes the scale; especially since fluctuations play an important role close to the phase transition where the correlation length diverges. Possible ways how the equations of motion can change upon coarse graining, as long as no global symmetries are violated, are: 1. Parameters tend to zero under coarse graining and are irrelevant for the critical dynamics. The dimensional analysis detailed in 'The Response Functional' shows that this is the cases for the noise terms η 1,2,3 and that the only relevant contribution of the sensitivity function is a constant sensitivity χ 0 .
2. Microscopic and mesoscopic relations between parameters might change. For example the microscopic relation between the reaction rates and the noise amplitude, i.e., Λ = µ + λ − γ no longer holds for the effective reaction rates and the effective noise amplitude on larger scales.
3. New types of effective interactions may arise at larger scales due to the presence of strong fluctuations close to criticality [39,51]. A well-known example for this is seen during real space RG schemes for the two dimensional Ising system where effective next-to-nearest neighbor interactions arise [52].
One example relevant to our model of points 2 and 3 is the effective chemotactic interaction, which, after reducing the sensitivity χ[ρ, ∇c] to χ 0 , is given by χ 0 ∇(ρ∇c). Even though this looks like a single term, one can, tentatively, separate this interaction into two terms: χ 0 on all scales since both terms combined are supposed to model a particle number conserving chemotactic current. However, this only holds in absence of any processes that explicitly break particle number conservation. In the combined model particle numbers are not conserved; hence, there is a priori no reason to expect the relation χ 0 , which might hold at some scales, to also be valid at macroscopic scales. The RG calculations below, indeed, show that in our model χ (1) 0 is in general not equal to χ (2) 0 , giving rise to a contribution (χ 0 )ρ∇ 2 c to the effective chemotactic interaction (also see 'Confirmation of Non-Conservative Interaction' and Fig. 2 in the main text). Thus, it is prudent to keep both terms separate with coupling parameters χ 1 and χ 2 . Combining all of the above finally results in the Langevin equations we rely on for the mean-field and RG analysis Equations (20) and (21) conclude the derivation of the correct macroscopic Langevin equations close to the continuous phase transition. It is important to emphasize that this is result is not specific to RG. RG is a systematic way of analyzing how -close to a continuous phase transition -effective interactions change with larger scales; therefore, the observation that the population dynamics and the chemotactic interaction mix upon coarse graining to give rise to an effective non-particle-number-conserving chemotactic interaction strongly suggests that this effect should be accounted for in any coarse-grained description of our model near the critical point. However, we make no prediction about the strength of this effect away from criticality or at small scales. In such cases the chemotactic interaction can most likely be formulated in terms of a conserved current and other noise terms (η 1,2,3 ) as well as higher orders of the sensitivity function χ[ρ, ∇c] may be important.

MEAN-FIELD ANALYSIS
As in the main text, we here analyze a model of chemotactic cells (ρ) reacting to gradients in a chemical signalling substance (c) given by the following Langevin equations: Here, D ρ represents the diffusion constant of the cells, θ is the effective linear growth rate, γ models a competition for resources, χ 1 and χ 2 are the chemotactic response parameters, and Λ the noise amplitude. The dynamics of the chemical is characterized by its diffusion constant D c , its degradation rate λ c , and its production by the cells at rate α. The noise term ξ has zero mean and is delta-correlated in time and space, i.e. ⟨ξ(x, t)⟩ = 0 and ⟨ξ(x, t) ξ(y, t ′ )⟩ = 2 δ(x − y) δ(t − t ′ ).
In this section, we first study the mean-field behavior of the model by ignoring the impact of the demographic noise. In this case, one finds two qualitatively different homogeneous steady states. The trivial solution ρ 0 = c 0 = 0 is coined the absorbing phase and the non-trivial, active phase solution is given by ρ 1 = θ/γ, c 1 = αρ 1 /γλ c . To determine their stability against small density fluctuations δρ(x, t) and δc(x, t), we expand around the homogeneous states such that Transforming to momentum space one finds an equation of the form ∂ t ϕ =Â · ϕ, with ϕ = (δρ, δc) and Note that only χ 2 but not χ 1 enters the equation at linear order in the perturbations. For the homogeneous solutions to be stable, we require ℜ(σ ± (k)) < 0 ∀k for the eigenvalues σ ± (k) ofÂ. Inserting ρ 0 = c 0 = 0, one findŝ from which we can read off the eigenvalues σ + (k) = −D ρ k 2 + θ and σ − (k) = −D c k 2 − λ c . Trivially, σ ± < 0 if and only if θ < 0 -i.e. the homogeneous absorbing state is always linearly stable if death dominates birth. In the case of the active solution, one getsÂ By calculating the eigenvalues of the zero mode, i.e. ofÂ 1 (k = 0), σ + (k = 0) = −θ and σ − (k = 0) = −λ c , one finds the active state to be unstable if θ < 0. For general k the eigenvalues read where the auxiliary functions b and c are given by For ρ 1 to be stable in a certain parameter regime, one requires ℜ(σ + (k)) < 0 ∀k, which is equivalent to c(k) > 0 ∀k. Thus, for θ > 0 the active state is unstable if and only if c(k) has a real root. This is the case only if B < 0 and B 2 > 4AC, which is equivalent to When this condition is satisfied, there exists a k for which σ + is positive and where the homogeneously active phase is thus unstable.
As we already argued for in the main part, θ ≈ λ c ≈ 0 is a necessary condition for scale invariant dynamics. However, it is not intuitively clear how to treat the fraction λ c /θ in the instability condition (31) in this limit. Accounting for fluctuation-induced shifts to the transition values θ * , λ * c , we write where we introduced the relative control parameters τ = θ − θ * and λ = λ c − λ * c Right at the transition this Eq. (32) reduces to θ * /λ * . As we show in the perturbative analysis, there is no diagram contributing to the renormalization of λ c . Thus, λ * c = 0, whereas ρ experiences a finite shift in its critical temperature. In this sense λ c /θ → 0 holds as one approaches the transition and equation (31) reduces to αχ 2 < −γD c .

The Long-Ranged and Quasi-Static Limit
For the chemotactic interaction to be intrinsically long-ranged -and thus for the system's dynamics to exhibit scale invariance -one needs to be in the limit where the decay rate of the signalling molecules becomes vanishingly small, i.e. λ c → 0. However, this limit is more subtle than one might expect. From equation (23), a formal solution for the time evolution of the chemical density c can be given in terms of the spatially Fourier-transformed density If we now take λ c → 0, we obtain for the homogeneous mode which is in general divergent for t → ∞. Note that this is the same divergence one encounters for the stationary solution c 1 when taking λ c → 0. Fortunately, the dynamics of ρ do not depend on the homogeneous mode of the chemical. This allows us to use the reduced quantityc which obeys which differs from the dynamics of c only by a homogeneous, albeit time dependent, term. Note that this homogeneous shift changes the steady state density toc 1 = 0 but does not alter the above linear stability analysis. By setting λ c = 0 in the instability condition (31), one obtains the simpler condition In addition to the long-ranged limit, we are also interested in the quasi-static limit D c /D ρ → ∞. Assuming D c ≫ D ρ , the approximate solution of Eq. (23) reads withρ = ρ −ρ and the average densityρ. In the second line we replaced the argument of ρ k with t -since the exponential, for D c /D ρ → ∞, only contributes to the integral at s = t -and then performed the integral. Given that α/D c is finite,c k is a solution to the Poisson equation This implies thatc(x, t) is quasi-stationary as it instantly adjusts toρ(x, t). Further, we notice that in order to make sense of this limit, one has to simultaneously assume that α ∼ D c (as also pointed out in Ref. [21]). To see that it is indeed necessary to handle the divergence in Eq. (34) with care, one can impose that the quasistatic limit must not alter the active steady state density. However, taking the limit without shifting toc would result in ∇ 2 c = −αρ/D c , thereby shifting ρ 1 to θ/(γ + αχ 2 /D c ).

THE RESPONSE FUNCTIONAL
To investigate the equations of motion (Eqs. (22) and (23)) beyond their mean field behavior -i.e. including the demographic noise term -we employ the dynamical renormalization group in form of the response functional formalism [53][54][55][56]. Following the approach displayed in [30], all moments of the fields can be written in the form of a path integral where we introduced the response fieldsρ (not to be confused with the average density) andc. The statistical weight is given by the action S = S 0 + S int which consists of the Gaussian part and the nonlinear interaction term where all terms after the integral signs are integrated over. Note that we included the linear term αcρ into S int rather than treating it as part of S 0 . This is not necessary, but simplifies identifying the effective couplings and calculating the Feynman diagrams. Our goal is to show that the terms included in Eq. (41) are the only relevant interactions for the renormalization procedure.
To determine whether an interaction is relevant or irrelevant, we introduce the momentum scale µ and calculate the naive scale dependence of the coupling parameters. From Eq. (40) it follows that [ρρ] = [cc] = µ d . By rescaling ρ → µ ξ ρ andρ → µ −ξρ some freedom in choosing the individual dimensions of the fields is left. Indeed, this freedom implies that it is still possible to choose the naive scale dependence of the different interaction terms; a seeming contradiction to the fact that the relevance and irrelevance of couplings cannot be arbitrary. However, this is only at first glance contradictory since it is a priori not obvious how the different vertices -and therefore the different field dimensionshave to be combined to create valid Feynman diagrams. For the action (41) we find that a γ-vertex can only appear together with a Λ-vertex and χ 1,2 -vertices only in combination with an αand a Λ-vertex. Consequently, only specific combinations of vertices have a defined relevance/irrelevance under the RG procedure. To correctly analyze which interactions are relevant for the RG analysis, it is therefore prudent to rescale the action in such a way that it contains as few dimensionfull parameters as possible. We choose where s 1 = Λγ −1 and s 2 = αD −1 c . Additionally we introduce the effective couplings which should, respectively, be interpreted as the standard coupling from directed percolation, the strength of the two different chemotactic interactions and a measure of the interaction time delay introduced by the finite diffusion speed of the chemical. Then the action reads Note that by rescaling time as t → tD −1 ρ and writing D c D −1 ρ = w(1 − w) −1 the action can be expressed as a function of u, g 1/2 and w only. We choose not to do so here, since this simplifies introducing all required renormalization factors. Note, however, that the perturbation series will only depend on w and not on the individual diffusion constants. Due to the specific form of the rescaling chosen in Eq. (42) the field dimensions can now be determined as  [30], the action S contains a third term of the form whose appearance can be traced back to diffusive noise (h 2 ) and a higher order contribution in the demographic noise (h 1 ). Dimensional analysis yields [h 1 ] = µ 2−d and [h 2 ] = µ −d/2 , which renders both irrelevant in d = 4 − ε dimensions. Thus, our initial choice of neglecting these terms in Eq. (22) is justified. The same holds true for any possible noise term in Eq. (23). Taking into account all the reactions associated to c and applying a Kramers-Moyal expansion, one can again follow the derivation of noise terms in [30] to obtain the additional contribution The first two terms result from demographic noise in c, the last from diffusional noise. is an analytic function (which is not necessarily the case [20]), the general contribution to the action reads The dimension of the coupling is then given by [h n1,n2,n3 ] = µ f , with Since, at d = 4, f (n 1 , n 2 , n 3 ) < 0 whenever n 1 or n 2 are greater than zero, all interactions resulting from such choices of n 1,2 are irrelevant for the RG calculations. For n 1 = n 2 = 0, on the other hand, all couplings are marginal at d = 4 dimensions as f (0, 0, n 3 ) = 0, independent of n 3 . This means that in order to fully renormalize the theory to infinite loop order, one, in principle, has to include infinitely many interactions and therefore infinitely many counter terms (see below). In that sense the action would no longer be renormalizable and one had to employ nonperturbative methods [44,45] to properly treat this set of relevant interactions. To avoid this subtlety, one has to make the stronger assumption χ[ρ, c] = χ[ρ, ∇c]. In this case, interactions with higher orders in c always come with at least one derivative acting on each c. It follows that n 3 ≤ n 1 and f < 0 at d = 4 unless n 1 = n 2 = n 3 = 0, implying that for all renormalizable theories no relevant contributions besides χ[ρ, ∇c] = χ 0 exist. Even though this assumption is unlikely to hold in a strict biological sense, chemotaxis is known to be robust over orders of magnitude of chemical concentration [57] and thus χ[ρ, c] = χ[ρ, ∇c] may be a reasonable approximation. Note that the above analysis completely relies on the naive scaling dimension of the involved coupling parameters. For instance, care has to be taken for the predictions in d = 2, since the vertex associated to h 1 is relevant for d ≤ 2. Altogether, we have that the actions displayed in (40), (41) and (44) indeed contain all relevant vertices (displayed in Fig. 4) and yield a minimal model for bacterial chemotaxis.

RENORMALIZATION
Having identified the relevant minimal action, one can apply a standard graphical perturbation expansion in terms of Feynman diagrams to calculate the flow equations to first non-trivial order. Details of such a calculation can be found in [30,58] or any textbook on quantum-field theory. Note that we perform all calculations in the absorbing state (θ < 0) close to the transition.
The central element of any perturbation expansion are the bare two-point Green's functions or propagators associated to the Gaussian part of the action (40), as well as the vertices displayed in Fig. 4. However, Green's functions calculated from this action turn out to be divergent. To renormalize the above theory, we rely on a multiplicative renormalization scheme, where we introduce the Z-factors for the quadratic part of the action as The nonlinear terms on the other hand, are renormalized via the choice This leaves us with a total of seven independent renormalization factors. The subscript R indicates dimensionless renormalized parameters (except for the temperatures θ and λ c , which still have a dimension). Note that from now on we will drop this subscript R as we exclusively work with the renormalized couplings. Following the usual steps, introducing the Z-factors gives rises to counter-terms which take the same form as the vertices shown in Fig. 4.
In principle the Z-factors can be determined iteratively to arbitrary loop orders by requiring every relevant vertex function to be finite. Yet, there is some freedom as one can always add finite terms to impose additional renormalization constraints. We choose not to do so and employ the so-called minimal subtraction (MS) scheme. Determining the Z-factors correctly to one loop order requires the calculation of all possible one loop Feynman diagrams. To keep the calculations as concise as possible, we use a list of frequently used standard integrals, identities and abbreviations listed in the section Standard Integrals and Identities.

Propagator Renormalization
In the following, we give a summary of the calculations of all one loop Feynman diagrams. More precisely, we determine their contributions to the renormalization of the corresponding Green's functions of the renormalized parameters. Consequently, we only give the divergent parts of the integrals since all finite contributions are irrelevant for the renormalization. We start with the renormalization of the propagator, which requires the calculation of the two diagrams shown in Fig. 5. To calculate the first diagram I P 1 we only have to identify the combinatorial factor associated with the diagram, the involved coupling constants and the propagators. Then we perform the ω-integration  and use Tab. III to give the final result. Thus, I P 1 can be evaluated to where we used the abbreviation ∆ = k 2 4 − iΩ 2Dρ − θ Dρ and employed the result for I 0,1 (∆) in the last line. The calculation of I P 2 needs more work and is given here step by step. The first few steps, albeit more tedious, are the same as before Since two propagators are left after the frequency integration one needs to employ the Feynman parameter identity Eq. 92. This was done in the last line with x = 1 − x and the use of the abbreviation Here we have defined the auxiliary functions δ(x) = x 2 + (1 − w)x k and Hence, the whole integral can be written as In the second line we got rid of all anti-symmetric parts of the integral, as well as the non divergent contributions. In the last line we made use of the fact that the divergent part of I 0,2 (∆) does not depend on ∆ and thus can be pulled out of the x integral. To get to the final result the following intermediate integrals need to be calculated: Additionally inserting I 0,2 (∆) = 1

Renormalization of DP couplings
For the renormalization of theρρρ-vertex we need to consider six diagrams (cf. Fig. 6). Their respective contributions are calculated in the following. Note that we set all external momenta and frequencies to zero since the DP vertices are momentum and frequency independent. From a dimensional analysis it can be further deduced that none of the divergences depend on the terms of O(q) and lower in the propagators (which would contribute to the momentum shift δ and the argument of the standard integrals ∆ in previous calculations). These terms are abbreviated as ·.
Also note that none of the Feynman parameter integrals depend on the Feynman parameters and the additional constants introduced by Eq. (92) and (93) always cancel each other. Having this in mind, the following calculations are performed in four steps: First all vertex factors, combinatorical prefactors and propagators are introduced; then the frequency integral is performed by identifying the poles of the propagators; subsequently all diffusion constants are extracted and the appropriate standard integral from Tab. III introduced; finally the result for the integral is inserted and the prefactors used to formulate the result in terms of the effective coupling constants. Following these steps one can give the results for the contributions of theρρρ-vertex renormalization as: To calculate the results for theρρρ-renormalization, we point out that there are three contributing diagrams which can be obtained by replacing theρρρ-vertex with aρρρ-vertex in I u 1 , I u 2 and I u 3 . Hence, the analytical results for the diagrams can be retrieved by just adding a minus sign in the respective calculations of theρρρ renormalization.

Renormalization of Chemotactic Couplings
The only remaining diagrams are the ones required for the renormalization of the chemotactic couplings (Fig. 7). Since these diagrams depend on external momenta, one cannot set them to zero and has to keep track of their contributions to the terms proportional to the loop momentum q in the propagators. The terms of order O(1) in q are again denoted by · and can be neglected. To calculate the contributions of I g 1 we, write down all the coupling constants and propagators, perform the frequency integration and introduce the Feynman parameters x and x = 1 − x (92): Where, in the fourth line, we defined δ(x) = − 1 2 (k + xp) and shifted q → q + δ. We now turn our attention to I g 2 . Simply inserting all the coupling constants and propagators yields Here we introduced A(q, k, p) as the product of the vertex factors for the twoρρc-vertices in the diagram: Note that no terms in A(q, k, p) are of order O(q 4 ). Additionally, there are four propagators in the diagram and after performing the ω integral three will be left. Thus, all diverging parts are proportional to I 2,3 (∆) orĨ 2,3 (∆) and, therefore, independent of ∆. Hence, all the parts contributing to ∆ are only denoted by · and the ∆ dependence of I 2,3 andĨ 2,3 dropped in the following calculation. Then, performing the frequency integral results in We continue by pulling out the diffusion constants and using the Feynman parameter trick: In the last line we shifted the loop momentum in both integrals by δ 1 and δ 2 , respectively, where Now, we need to calculate the shifted numerator A(q + δ): A(q + δ) = (g 2 1 − g 1 g 2 ) (δp) q 2 + 2(δq)(pq) − g 2 1 (pq)(kq) + (g 1 g 2 − g 2 1 )(kp) q 2 − g 2 1 (pq) 2 + (g 2 2 − g 1 g 2 )p 2 q 2 Inserting δ 2 = 2(1 − w) δ 1 into the integral, one can anticipate the appearance of the following expressions: With these results and x,y,z 1 = 1/2, one can separate the δ 1 dependent terms and get x,y,z 2Ĩ 2,3 (p, δ 1 ) + I 2,3 (pδ 1 ) .
Inserting the results for I 2,3 ,Ĩ 2,3 and x,y,z δ 1 = 1 12 (2k + p) yields: Collecting and grouping all the contributions yields the final result: Finally, we need to calculate where A(q, k, p) again denotes the product of the chemotactic vertices.

Flow Equations II
All fixed points, together with the associated eigenvalues and scaling exponents are shown in Tab. II. Notably, the exponents do not vary along the fixed line (CR). However, this was to be expected, since the fixed line collapses to a fixed point upon a change of variables, as shown in the main text. As the flow equations describe a four dimensional space, it is in general not possible to give an exhaustive visualization of the flow of the system. Only in specific cases does the flow remain in a lower dimensional hyperplane (such as the g 1 -ū-plane at w = 1 as described in the main text). In the rest of the cases we are obliged to ignore the exact flow behavior and focus on where a flow line starting at a point with initial coordinates (u, g 1 , g 2 , w) i ends up. The collection of points that flow towards a certain fixed point makes up its basin of attraction, and we can visualize two-dimensional slices of this space to ascertain which regions in parameter space are controlled by which fixed point. In the main text we extensively treat several cases. In Fig.8 we plot the basins of attraction and the topology of the four dimension flow in different g 1 -ū-planes at various fixed values of u and w. Notably, the topology of the flow diagram is independent of u for w = 1. This relates to the fact that in this limit the proper effective variable is given byū = g 2 + u as argued for in the main part of this letter. One observes that once the (w = 1)-plane is left, a new region of runaway flow appears at g 1 < 0. Interestingly, this region grows in a winding fashion, increasing in size as w decreases, at the same time causing the basin of attraction of the CA fixed point to shrink. The dependence of the size of this 'wedge' of runaway flow on the parameter w can be studied by defining an angle ψ between the line defined by g 1 = 0 and the boundary between the runaway flow and the basin of attraction of the CA fixed point. In Fig. 9 we observe that for a relatively large range of w this angle is very small indicating a negligible region of runaway flow and a phase diagram that is not very different from that at w = 1. The CR fixed line becomes unstable for w < 1, but as this instability is relatively weak, one can expect the large scale behavior to be similar to that of w = 1 on both sides of the g 1 = 0 invariant manifold. For smaller w the influence of u starts to become more pronounced, distorting the boundaries between the different runaway regions.
Note that apart from the basins of attraction of the CA and CR fixed points we define four different types of runaway flow. The dark gray region is defined as the runaway that lies below the basin of attraction of CA for g 1 < 0. The light gray region and the striped blue region lie in the g 1 > 0 plane, and are divided on the basis of flow behavior. Flow in the striped blue region is affected by the attractive nature of the projection of the CR fixed point below w < 1, whereas flow originating in the light gray region is not and is, therefore, associated to the runaway flow already present at w = 1. The gray striped region corresponds to runaway linked to the CP fixed point, and is given by the runaway that lies above the basin of attraction of the CA fixed point for g 1 < 0.
From Fig. 8 it is obvious that the CR fixed line becomes unstable for w < 1. By noting that we can relate the flow equation of w to the one of D ρ via dw d ln(µ) we can easily explain why this was to be expected. In the vicinity of a super-diffusive fixed point (β * D < 0), Eq. (86) implies that w decreases as µ → 0. By the same argument, the opposite is to be expected at a sub-diffusive fixed point. Therefore, the super-diffusive fixed line CR has to be unstable in w-direction. On the other hand, the subdiffusive fixed point CA is expected to be stable in w-direction, explaining why its basin of attraction extends in the w-direction. Moreover, it is apparent that z = 2, i.e. β D = 0 has to hold for any fixed point at w ̸ = [0, 1]. Thus it is clear that the CP fixed points obeys z = 2 to all loop orders.
Apart from the fixed points, another important feature of the flow equations (74)-(77) are its invariant manifolds. Inspecting Eq. (76) it is clear that ∂g 1 /∂lnµ| g1=0 = 0; thus, g 1 = 0 is such an invariant manifold. Moreover, one can show that this result is true to any loop order which can be understood by inspecting (62): The part of the result contributing to the g 1 renormalization is proportional to g 1 , whereas the contribution to the g 2 renormalization is not proportional to g 2 . Consequently, dividing by g 1 and multiplying with g 1 (which essentially leads to the contribution of this diagram to Eq. (76)) results in something proportional to g 1 . Repeating this procedure for g 2 , one realizes that this results in a contribution to the g 2 -flow that is not proportional to g 2 , thus allowing the flow to cross the g 2 = 0 hyperplane. Therefore it is sufficient to show that to all loop orders all divergences proportional to kp are also proportional to g 1 . This is the case, because one had to take the g 2 term of the vertex factor of every chemotactic vertex for the contrary to be possible. In particular this includes the g 2 p 2 part from the vertex where the incoming c-field connects with the rest of the diagram. However, since this is already proportional to p 2 , it can no longer renormalize g 1 , proving that at least one factor of g 1 is included in every kp-divergence. Hence, the g 1 = 0 hyperplane can never be crossed. The behavior of the angle ψ between the line defined by g1 = 0 and the boundary between the runaway flow given by the gray striped region and the CA basin of attraction in Fig. 8 as w is varied. We note that the ψ interpolates between ψ = 3/4π, implying no orange region, at small w and ψ = 0, implying no runaway, at large w. In must be remarked that ψ is defined in the case of u = 1 where the boundary can be approximated by a straight line and the angle can thus be taken as a good representation of the size of the striped runaway region. However, as in Fig. 8 it is clear that the growth of the region is very similar for all u, we believe the behavior to hold qualitatively in general.

CONFIRMATION OF NONCONSERVATIVE INTERACTION
In section 'Derivation of Langevin Equations' the different impacts of the RG flow on the effective equations of motion are explained and a non-particle-number-conserving effective chemotactic interaction is proposed. In this section we demonstrate how the generalized chemotactic interaction χ 1 ∇(ρ∇c) + (χ 2 − χ 1 )ρ∇ 2 c arises already at one-loop level. To this end, we analyze how a conserved chemotactic interaction is modified during the RG step. Starting from the classical Keller-Segel nonlinearity χ 0 ∇(ρ∇c), one can derive the RG flow functions by inserting g 1 = g 2 = g 0 into Eqs. (74)-(77). If the resulting flow equations for g 1 and g 2 are identical, i.e., the g 1 = g 2 hyperplane constitutes an invariant manifold, it is possible to renormalize the theory with a single effective coupling constant. There are three Feynman diagrams that contribute to these flow functions and their respective values for g 1 = g 2 = g 0 are I g 1 g0 = √ 32π 2 u −1 D ρ µ ε/2 ug 0 ε 2(kp) + 3p 2 (87) (−3 + 2w) kp + p 2 (88) Here, all the terms proportional to kp and p 2 contribute to the renormalization of g 1 and g 2 , respectively. Importantly, one recognizes that while I g 2 and I g 3 contribute equally to both flow equations -being consistent with the Keller-Segel nonlinearity -only the diagram I g 1 which couples the chemotactic vertex with the resource limiting nonlinearity breaks this relation. This is crucial, since it shows that performing a single RG step in the presence of resource limitation generates a nonconservative contribution to the chemotactic interaction also if it is not included from the beginning. Thus, a consistent coarse graining of the theory with a conservative effective chemotactic interaction close to criticality is not possible and an additional term needs to be included.