The role of mobility in epidemics near criticality

Beatrice Nettuno^*, Davide Toffenetti^*, Christoph Metzl^*, Linus Weigand^*, Florian Raßhofer, Richard Swiderski, and Erwin Frey
Arnold Sommerfeld Center for Theoretical Physics and Center for NanoScience, Department of Physics, Ludwig-Maximilians-Universität München, Theresienstraße 37, D-80333 Munich, Germany

$^*$ These authors contributed equally. B Nettuno and D Toffenetti performed the analytical work, C Metzl and L Weigand did the numerical one.

August 27, 2024

The general epidemic process (GEP), also known as susceptible-infected-recovered model (SIR), describes how an epidemic spreads within a population of susceptible individuals who acquire permanent immunization upon recovery.
This model exhibits a second-order absorbing state phase transition, commonly studied assuming immobile healthy individuals.
We investigate the impact of mobility on disease spreading near the extinction threshold by introducing two generalizations of GEP, where the mobility of susceptible and recovered individuals is examined independently.
In both cases, including mobility violates GEP's rapidity reversal symmetry and alters the number of absorbing states.
The critical dynamics of the models are analyzed through a perturbative renormalization group approach and large-scale stochastic simulations using a Gillespie algorithm.
The renormalization group analysis predicts both models to belong to the same novel universality class describing the critical dynamics of epidemic spreading when the infected individuals interact with a diffusive species and gain immunization upon recovery.
At the associated renormalization group fixed point, the immobile species decouples from the dynamics of the infected species, dominated by the coupling with the diffusive species.
Numerical simulations in two dimensions affirm our renormalization group results by identifying the same set of critical exponents for both models.
Violation of the rapidity reversal symmetry is confirmed by breaking the associated hyperscaling relation.
Our study underscores the significance of mobility in shaping population spreading dynamics near the extinction threshold.

§ INTRODUCTION

The recent COVID-19 pandemic has sparked an increased interest in the mathematical modeling of disease and virus spreading [1, 2, 3].
Various quantitative models have been successfully applied to aid vaccination and containment strategies [4, 5, 6, 7].
This underscores the critical importance of understanding epidemic processes.
However, realistic modeling efforts remain challenging: Many factors such as human mobility [8] and social behavior [9, 10] as well as spatial heterogeneity [11, 12] are known to play pivotal roles in shaping the course of epidemics.

Naturally, this calls for exploring simpler spreading models of epidemics that may be applied to an even broader class of phenomena, including opinion dynamics [13, 14], forest fires [15, 16] and even catalytic surface reactions [17, 18].
While inherently less intricate, these conceptual models are crucial in identifying common themes and providing qualitative insights [19, 20, 21].
Typically, one considers lattice-gas models, where active agents ($B$), such as a burning tree or an infected individual, activate nearby susceptible agents ($A$) [22, 23].
In its most basic form, this may be accounted for by a local reaction

\begin{eqnarray}
A+B \overset{\sigma}{\rightarrow} B + B \, \nonumber
\end{eqnarray}

at constant rate $\sigma$. Over time, agents lose their activity, e.g. by recovery of infected individuals, with deactivation being either permanent or transitory.
If the activation process proves stronger, activity (e.g. disease or fire) may persist and spread throughout the system.
Conversely, if deactivation dominates and there is no spontaneous creation of active or susceptible agents, the system approaches an inactive absorbing state where active individuals are absent.
The transition between these two phases is coined an absorbing state phase transition [23, 24], which is a class of non-equilibrium phase transitions.

Like their equilibrium counterparts, the near-critical behavior of absorbing state phase transitions can be categorized into a limited number of universality classes, the most important ones being directed percolation (DP) [24, 25, 26] and dynamical percolation (DyP), here referred to as general epidemic process (GEP) [27, 28].
The contact process (CP) [29] gives one particular realization of the DP universality class, which acts as a minimal model of disease spreading when there is no immunization.
It is formulated in terms of local reactions

\begin{eqnarray}
A_{i}+B_{j} \overset{\sigma}{\rightarrow} B_{i} + B_{j} \, , \quad
B_{i}\overset{\zeta}{\rightarrow} A_{i} \, , \nonumber
\end{eqnarray}

with at most one agent per site and ($i,j$) denoting adjacent lattice sites.
Due to the absence of immunization, it is often referred to as the susceptible-infected-susceptible (SIS) model.
Whenever immunization is permanent, one reverts to the susceptible-infected-recovered (SIR) model [30] which falls into the GEP universality class [27, 28].
Its associated reaction scheme reads

\begin{eqnarray}
A_{i}+B_{j} \overset{\sigma}{\rightarrow} B_{i} + B_{j} \, , \quad
B_{i}\overset{\zeta}{\rightarrow} C_{i} \, , \nonumber
\end{eqnarray}

with immune (recovered) individuals being represented by an additional species $C$.
The crucial difference between these two reaction schemes is that for the CP/SIS dynamics, an epidemic can survive locally. In contrast, in the GEP/SIR model, activity may only persist as a traveling wave spreading through the population.
A generalized model that allows for reinfection at a varying rate contains both of these scenarios as limiting cases [31].

Hitherto, one crucial aspect has been omitted in these classical conceptual models:
In general, activity may not only propagate due to the infection dynamics but also strongly depend on the mobility of susceptible and infected agents.
In its most elementary form, this can be accounted for by each agent undergoing diffusive motion.
It is known that more long-ranged forms of mobility, such as Lévy-flights [32] or chemotactic motion [33] can affect the behavior close to a non-equilibrium phase transition.
In this work, we exclusively focus on diffusive motion, which alters the local environment of active agents.

Previous studies show that the inclusion of diffusion alone can significantly change the observed critical behavior:
While the ordinary pair contact process (PCP) [34] falls into the DP universality class [35, 36], the pair contact process with diffusion (PCPD) [37, 38] shows deviating critical behavior.
However, the nature of its universality class remains a subject of ongoing discussions [39, 40].
Moreover, numerical simulations indicated that SIR dynamics with symmetric exchange hopping culminates in a small albeit significant deviation from the known GEP universality class [41].

For mobile agents, it is further necessary to revisit the occupation number constraint:
If each site can be occupied by at most one agent, diffusion may only occur in the form of two-agent exchange hopping, as illustrated in reference [41]; it is impossible to account for the individual mobility of different species.
Further, in a coarse-grained description, nearest neighbor (n.n.) infections give rise to gradient terms akin to a classical diffusion law for the active species [42].
Thus, a clear distinction between the two modes of activity spreading (individual mobility and nearest-neighbor infections) cannot be made on a mesoscopic level.
This rationalizes why the contact process (SIS model) with exchange hopping of infected and susceptible agents and a constrained agent number remains in the DP universality class [42].

In contrast, relaxing the exclusion constraint and, thus, allowing for distinct mobilities of susceptible and infected agents gives rise to the diffusive epidemic process (DEP) [43].
Its emergent critical behavior is only partially understood and crucially depends on the relative mobilities of susceptible and recovered agents:
The limit of immobile susceptible agents was shown to fall into the universality class of quenched Edwards-Wilkinson growth by an exact mapping [44, 45].
Field-theoretical studies have identified two new classes of critical behavior [43, 46, 47] for infected agents diffusing faster or equally as fast as the susceptible ones.
When diffusion of susceptible agents is faster, neither perturbative [43] nor non-perturbative methods [48] could identify an entirely conclusive fixed point structure.
DEP was even speculated to undergo a fluctuation-induced first-order transition [49].
However, follow-up numerical studies have questioned this possibility, reporting a continuous transition for the latter regime [50, 51, 52].
The ongoing debate about the nature of the phase transition and the associated critical exponents was addressed in a recent publication for one-dimensional systems [53].
By unambiguously showing that the transition is continuous, this numerical study supports the hypothesis of a strong coupling fixed point and further reported sub-diffusive activity spreading.

The above examples illustrate how unraveling the role of mobility in epidemic spreading remains a field of ongoing research, with a conclusive answer yet to be found.
Therefore, in this study, we investigate how the mobility of susceptible and recovered agents affects critical behavior in a generalized version of the SIR/GEP model.
For that purpose, we employ a field-theoretical Renormalization Group (RG) analysis [22] as well as extensive numerical simulations.

Our modifications aim to explore different possibilities of coupling infected individuals with a diffusive species.
In particular, an additional caregiving reaction is introduced.
This enables us to consider two limiting cases in which the role of mobile susceptible and recovered agents can be studied independently.
Moreover, we exclusively consider systems without occupation number constraints, allowing for a clear distinction of different mobilities for different species and a consistent field-theoretical description.
Regardless of how the GEP model is generalized by introducing additional diffusivity of the different species, our results show a consistent critical behavior that deviates from the conventional GEP paradigm and thus establishes a new universality class.

Our theoretical prediction of a new universality class is further corroborated through exact stochastic simulations.
Not considering any exclusion constraint, our numerical approach complements and shows remarkable agreement with earlier numerical results obtained in [41].
This enables us to state beyond doubt that mobility is a critical perturbation to the important GEP universality class.
Moreover, the presented analysis offers a clear understanding of how mobile agents alter the universality class and thereby may shed light on analogous cases such as DEP or PCPD.

This paper is organized as follows:
In section <ref>, we introduce our modifications of the SIR/GEP reaction scheme and discuss their phenomenological implications.
Particular emphasis is put on limiting cases that we consider to address our questions. We further introduce the observables analyzed in this paper.
Section <ref> provides a mean-field analysis to gain intuition on our models' absorbing state phase transitions.
We elucidate how diffusing agents influence deterministic dynamics, specifically examining the emergence of diffusive currents close to the propagating front of infected agents in the active phase.
Section <ref> presents a perturbative renormalization group analysis, providing analytical evidence for a new universality class in the form of a new IR-stable fixed point.
In addition, we discuss under which circumstances the mobility of either susceptible or recovered individuals violates an intrinsic symmetry of GEP and alters the absorbing state in the model, ultimately modifying the RG flow functions.
Section <ref> concisely describes the algorithm used in the simulations and presents the key numerical results, including critical scaling exponents and their estimated error.
We compare our results to the literature where possible and conclude that the numerical analysis underpins all our analytical predictions of the new universality class.
The last section summarizes our findings and discusses their implications with possible future research avenues.
Further, we provide detailed appendices that include a derivation of the field-theoretical action, a complete summary of the perturbative RG calculations, a description of our numerical approach and error estimates, and additional data.

§ REACTION-DIFFUSION MODELS OF EPIDEMIC PROCESSES

The epidemic models considered in this paper are formulated as reaction-diffusion models on a bosonic $d$-dimensional lattice—wherein an arbitrary number of individuals can occupy each site. We consider a population consisting of three species: $A$ represents the individuals susceptible to the disease, $B$ the infected, and $C$ the recovered members immune to the disease.
We refer to susceptible and recovered individuals collectively as healthy individuals.
The general class of epidemic models  we consider is given by the stochastic reaction scheme:

\begin{eqnarray}
\fl \mbox{infection} \qquad &&A_{i}+B_{i}\overset{\sigma}{\rightarrow} 2 B_{i}, \label{eq:chemical.rates1} \\
\fl \mbox{recovery}\qquad &&B_{i}\overset{\zeta}{\rightarrow} C_{i}, \label{eq:chemical.rates2} \\
\fl \mbox{caregiving} \qquad &&C_{i}+B_{i}\overset{\gamma}{\rightarrow} 2 C_{i} \label{eq:chemical.rates3}
\end{eqnarray}

\begin{eqnarray}
\fl  \mbox{diffusion of }A \qquad &&A_{i} \overset{D_A}{\rightarrow}A_{j},
\label{eq:diffusion.rates1} \\
\fl \mbox{diffusion of }B \qquad  &&B_{i} \overset{D_B}{\rightarrow}   B_{j},
\label{eq:diffusion.rates2} \\
\fl \mbox{diffusion of }C \qquad  &&C_{i}\overset{D_C}{\rightarrow}C_{j}.
\label{eq:diffusion.rates3}
\end{eqnarray}

The lower indices $i,j$ indicate the lattice sites with $(i,j)$ being neighboring sites.
The first reaction eq:chemical.rates1 describes the infection of a healthy individual ($A$) by a sick individual ($B$) with rate $\sigma$.
The second reaction eq:chemical.rates2 represents the recovery from the disease with rate $\zeta$. After recovery, full immunization ($C$) is assumed, implying that each individual can be sick only once.
Compared to the standard GEP reaction scheme, we introduce an additional “caregiving" reaction eq:chemical.rates3, where immune individuals ($C$) help the infected ($B$) to recover with a rate $\gamma$.
Finally, equations eq:diffusion.rates1–eq:diffusion.rates3 represent the mobility of each of the three species $A$, $B$ and $C$, which can hop from a site i to a nearest-neighbor site j with rates $D_A$, $D_B$, and $D_C$, respectively.
While the reactions eq:chemical.rates1–eq:chemical.rates3 are on-site processes, the nearest-neighbor hopping eq:diffusion.rates1–eq:diffusion.rates3 introduces a non-trivial spatial dependence.
All the reactions described above are stochastic, meaning that the reaction rates represent the probability per unit of time for the corresponding reaction to occur.
This general model exhibits an absorbing state phase transition between active and inactive phases:
In the active phase, the disease has a non-zero probability of spreading through the whole system, and the mean number of infected individuals increases.
In contrast, in the inactive phase, the epidemic subsides, and the system reaches the absorbing state, where infected individuals are no longer present.
The spatiotemporal distribution of the three types of individuals in this system is dictated by the dynamics of the infected population, denoted as $B$.
Therefore, the density of $B$ functions akin to an order parameter in phase transitions, and we will adopt this terminology throughout this paper. However, in the active phase, the density of $B$ does not uniformly reach a non-zero value across space.
Instead, it is non-zero along an expanding annulus since the epidemic spreads through the system as a traveling wave, as depicted in fig:1.

This paper analyzes how accounting for the healthy species' diffusivity affects the order parameter's critical dynamics.
We study three limiting cases of the general reaction scheme eq:chemical.rates1–eq:diffusion.rates3, all exhibiting a second-order absorbing state phase transition: the general epidemic process (GEP), the diffusive general epidemic process (DGEP) and the general epidemic process with caregiving (CGEP).
The reaction scheme of the GEP consists of eq:chemical.rates1, eq:chemical.rates2 and eq:diffusion.rates2, i.e., it assumes immobile healthy individuals and disregards the caregiving reaction (${D_A=D_C=\gamma=0}$).
This model was first proposed in 1927 [54] but was formulated in terms of stochastic processes only much later in [27].
At first, the critical properties of the GEP were believed to be different from other known universality classes (see, e.g. [55]). Still, later, it was established that the GEP transition lies in the universality class of Dynamical Percolation (DyP) [56, 28], even though the properties of these two systems far away from criticality differ [57].

GEP is the paradigmatic model for disease spreading with immunization, but it disregards the mobility of healthy individuals.
To study the role of mobility of the healthy species $A/C$, we consider two distinct limiting cases of the general reaction scheme eq:chemical.rates1–eq:diffusion.rates3. In these, mobility of susceptible and recovered individuals is examined independently.
To account for the mobility of susceptible individuals, we introduce a model variant of GEP with a finite diffusion constant $D_A \neq 0$, consisting of the reactions eq:chemical.rates1, eq:chemical.rates2, eq:diffusion.rates1 and eq:diffusion.rates2;  we refer to it as DGEP. Note that in both GEP and DGEP, the recovered individuals $C$ can effectively be treated as dead individuals since they do not interact with other species; thus, the mobility of $C$ cannot influence the dynamics of $B$.

To account for the mobility of recovered individuals, we introduce the CGEP model (eq:chemical.rates1–eq:chemical.rates3, eq:diffusion.rates2 and eq:diffusion.rates3), where a “caregiving" reaction eq:chemical.rates3 is taken into account, coupling the infected population with the recovered one.
We consider the case where the recovered population is mobile (${D_C\neq 0}$), while the healthy population is immobile (${D_A= 0}$).
This results in infected individuals ($B$) interacting with both mobile individuals ($C$) through caregiving (see equation eq:chemical.rates3) and immobile individuals ($A$) through infection (see equation eq:chemical.rates1). Consequently, CGEP allows us to investigate whether the dynamics at criticality are predominantly influenced by coupling with a diffusing or non-diffusing field.

Models and phenomenology.
(a) Reactions in the general model: On-site reactions eq:chemical.rates1–eq:chemical.rates3 (left) and diffusion eq:diffusion.rates1–eq:diffusion.rates3(right) as a stochastic nearest-neighbor hopping process on a 2D lattice (right).
The reaction scheme of the studied sub-models can be obtained in the following limits: GEP ($D_C=D_A=\gamma=0$), DGEP ($D_C=\gamma=0$), CGEP ($D_A=0$). (b) Exemplary stochastic simulation of an epidemic with recovery in the active phase, starting from a single seed (left-most panel) and showing an annular spreading of infected individuals depicted as blue dots (other species not shown).
The inset shows the 2D square lattice on which the simulations were performed. In this simulation, ten infected individuals $B$ are placed initially at the center of the system with a homogeneous density of susceptible individuals $A$ (not shown) and no recovered individuals $C$.

§.§ Observables

Epidemic models intend to describe how a disease spreads through a population of healthy individuals.
To quantitatively characterize the spreading, it is sufficient to consider a finite set of observables defined as an average over many stochastic realizations. Here we introduce three observables, which are typically considered for epidemic spreading:
The average number of infected individuals $\langle N\rangle(t)$ at time $t$, the survival probability $P_{\rm s}(t)$ of the epidemic, which represents the probability of having at least one infected individual left at time $t$, and the mean square radius $\langle R^2 \rangle (t)$ of the cluster of infected individuals.
The qualitative behavior of these observables depends on whether the system is in the active or inactive phase.
In the active phase, where the disease has a non-zero probability of spreading through the whole system, $\langle N\rangle(t)$ and $\langle R^2 \rangle (t)$ grow indefinitely, while $P_{\rm s}(t)$ approaches a non-zero value for $t\rightarrow\infty$.
shows the typical behavior of an epidemic process with immunization in the active regime.
In contrast, infected individuals go extinct in the inactive phase, such that both $\langle N\rangle(t)$ and $P_{\rm s}(t)$ approach zero for long times, and the mean square radius is not defined anymore.
At the phase transition, the system exhibits scale invariance, with observables following a power-law scaling [21, 23].
Consequently, the three observables introduced previously take the following asymptotic form [22]:

\begin{eqnarray}
\langle N \rangle (t,\sigma)
&=& t^{\theta} \hat{N}\left( (\sigma-\sigma_c)^{\nu_{\parallel}} t \right),
\label{eq:observable1}\\
P_{\rm s}(t,\sigma)
&=& t^{-\delta} \hat{P}_{\rm s}\left((\sigma-\sigma_c)^{\nu_{\parallel}} t\right) \, ,
\label{eq:observable2}
\\
\langle R^2 \rangle (t,\sigma)
t^{\frac{2}{z}} \hat{R}^2\left((\sigma-\sigma_c)^{\nu_{\parallel}} t\right)
\, .
\label{eq:observable3}
\end{eqnarray}

where $\sigma_c$ is the critical infection rate at which the system exhibits the absorbing state phase transition.
The four exponents $\theta, \delta, \nu_{\parallel}$, and  $z$—called critical exponents—characterize a universality class and are common to all phase transitions that belong to it. The functions $\hat{N}( x )$, $\hat{P}_{\rm s}(x)$ and $\hat{R}^2(x)$—called scaling functions—contain both a universal as well as a model-specific part; they are, however, difficult to infer analytically. Therefore, they are not the object of study in this paper.

§ MEAN-FIELD ANALYSIS

To gain some intuition on the various epidemic models and the role of diffusion, we first discuss the dynamics in continuous space without considering stochastic effects.

Mean-field analysis. The spreading of the density of infected $B$ (a), susceptible $A$ (b) and recovered $C$ individuals (c). The data was obtained through a finite element method (FEM) simulation of equations eq:deterministicequations1–eq:deterministicequations3 in a parameter regime corresponding to the active phase and using seed initial conditions. The blue squares, yellow circles, and purple triangles in (a)–(d) correspond to phase space points where the wave of $B$ respectively has not yet passed, it is passing or has passed. (d) The well-mixed phase space of the densities $A(t)$ and $B(t)$ of equations eq:deterministicequations1 and eq:deterministicequations2, in the ${C=0}$ plane. The dynamics exhibit a fixed line ${B^*=0}$, changing stability at the threshold value ${A^{ \rm mf}_c=\zeta/\sigma}$ at the boundary between the pink and blue line. (e)–(f) The influence of mobility on the profile of the densities $A$ and $C$ as a function of the radius $r$. The value approached at a large radius by $A$ and $C$ is determined by the initial conditions and is therefore equal for the diffusive and non-diffusing case. In (e), profiles with and without diffusivity for the susceptible species are shown in orange and gray, respectively, while in (f), they are shown in purple and gray. The orange (e) and purple (f) arrows show the direction of diffusive currents of susceptible and recovered individuals.

We use a mean-field theory where the system's description is no longer based on the number of individuals at each lattice site,
but formulated in terms of the three scalar densities $A(\bi{r},t),B(\bi{r},t)$, and $C(\bi{r},t)$.
The corresponding mass-conserving reaction-diffusion equations [58] for the general reaction scheme eq:chemical.rates1–eq:diffusion.rates3 are given by:

\begin{eqnarray}
\partial_{t}A (\bi{r} ,t)
D_A \nabla^2 A -\sigma A B \, , \label{eq:deterministicequations1} \\
\partial_{t}B (\bi{r} ,t)
D_B \nabla^2 B +\sigma A B- \gamma B C-\zeta B \, ,\label{eq:deterministicequations2}\\
\partial_{t}C (\bi{r} ,t)
D_C \nabla^2 C+ \gamma B C+\zeta B \, . \label{eq:deterministicequations3}
\end{eqnarray}

The equations of motion for the different models introduced above are obtained by setting specific parameters to zero: for the GEP, ${D_A = D_C = \gamma = 0}$; for the DGEP, ${D_C = \gamma = 0}$; and the CGEP, ${D_A = 0}$.

§.§ Well-mixed system

To better understand the type of phase transition exhibited by these epidemic models, we first focus on the well-mixed case, disregarding any spatial dependence and gradients in equations eq:deterministicequations1–eq:deterministicequations3.
Since the total number of individuals ${A+B+C}$ is conserved, it is sufficient to consider the dynamics of densities $A(t)$ and $B(t)$ only.
We choose the initial conditions such that no recovered individuals are present, $C(0)=0$, imitating the initial conditions used later in stochastic simulations (see e.g. fig:1(b)).
fig:2(d) illustrates the dynamics of the well-mixed system in the phase space spanned by $A$ and $B$.
Independent of their initial values, the dynamics ultimately reach a stationary state where infected individuals are extinct (${B^* = 0}$).
This is because an epidemic with immunization cannot survive in loco, i.e., in a well-mixed analysis, but has to spread through space to survive. In particular, the line ${B=0}$ is unstable for ${A^*> \zeta / \sigma}$ (light blue line in fig:2(d)) and stable for $A^*< \zeta / \sigma$ (pink line in fig:2(d)).
Stated differently, an epidemic transiently emerges in a well-mixed system only when ${A^*> \zeta / \sigma}$, as illustrated by the flow lines shown in fig:2(d).
This scenario differs from an epidemic lacking permanent immunization, such as the diffusive epidemic process (DEP) defined by the reactions ${A+B \to 2 B}$ and ${B \to A}$ [53].
In that case, a well-mixed system exhibits a transcritical bifurcation marked by the appearance of a stable active fixed point (${B^{*} \neq 0}$), allowing the disease to persist in loco.

§.§ Spreading dynamics

Next, we performed
simulations of the spatially extended model  eq:deterministicequations1–eq:deterministicequations3 using a Finite Element Method (FEM) [59].
We are interested in the spreading of a disease through a population of susceptible individuals. Therefore, we chose initial conditions with a spatially uniform distribution of susceptible individuals ${A(\bi{r},0) = A_0}$, no recovered agents ${C(\bi{r},0)=0}$ and a small density of infected individuals—a seed—placed in the center of the system. In continuous space the seed is represented by a Gaussian distribution ${B(\bi{r},0)= B_0 \exp(- |\bi{r}|^2 / \rho^2)}$, where $B_0$ and $\rho$ are small compared to the other model's parameters.
We term this type of initial conditions seed initial conditions.

As suggested by the phase space flow of the well-mixed system (fig:2(d)), the spreading dynamics are qualitatively different depending on the initial homogeneous density $A_0$ and the values of the infection and recovery rates, $\sigma$ and $\zeta$, which are the tuning parameters of the deterministic phase transition.
Suppose the system is initialized with $A_0$ above the threshold density ${A_0 > A^{\rm mf}_c = \zeta/\sigma}$. In that case, the seed grows, and the infection spreads throughout the system as an annular traveling wave, as depicted in figure <ref>(a). This happens because the system, except for the location of the seed,  is initialized at uniform density values $(A_0,0,0)$ corresponding to an unstable fixed point in the well-mixed analysis (depicted as light blue squares in figures <ref>(a)–(d)); thus, enabling the perturbation to grow.
The susceptible individuals \(A\) at the leading edge of the wavefront are gradually converted by infected \(B\) individuals.
After the wave has passed, these infected individuals recover such that the density of the susceptible individuals drops below the threshold, ${A < A^{\rm mf}_c}$ (pink triangle in figures <ref>(a)–(d)), resulting in the extinction of the epidemic behind the trailing edge of the wavefront.
This spreading scenario corresponds to the active phase.
Conversely, if the initial concentration of
susceptible individuals is below the threshold, ${A_0 < A^{\rm mf}_c}$, the perturbation induced by the seed of sick individuals decays and the epidemic goes extinct:
The system is initialized at uniform density values $(A_0,0,0)$, corresponding to a stable fixed point in the well-mixed analysis so that any perturbation decays.
The model is the inactive phase in this scenario. The absorbing state phase transition is located on the active and inactive phase boundary at the deterministic threshold ${A_0 = A^{\rm mf}_c}$.

§.§ Diffusive currrents

Last, we discuss how the mobility of healthy individuals ($A$/$C$) influences the spreading dynamics in the generalized model.
fig:2(e) and (f) illustrate the density profiles of susceptible ($A$) and recovered ($C$) individuals in the spatially extended system.
The density of susceptible individuals is low where the wave of infected individuals has already passed and high where it has not since the susceptible individuals get infected by the $B$ individuals in the wave.
On the other hand, the spatial profile of the recovered individuals ($C$) is inverted compared to the spatial profile of susceptible individuals ($A$):
High-density regions in $A$ coincide with low-density areas of $C$ and vice versa.
This is because infected individuals are converted into recovered ones by the caregiving and recovering reactions, increasing their number after the wave of infected individuals has passed.
A comparison between fig:2(a),(b), and (c) shows that the boundary between the low/high-density regions of $A$ and $C$ is located at the wavefront of infected individuals $B$.
If $A$ or $ C$ are mobile, the difference in density induces a diffusive current from regions of high to low density (see figures <ref>(e) and (f)).
These currents have two notable effects:
(i) The diffusivity of $A$ and $C$ tends to minimize gradients by smearing out the interface between high and low-density phases.
(ii) The currents transport individuals from high-density regions to low-density ones, resulting in a higher density of $A$ and a reduced density of $C$ after the wave has passed as compared to the case where $A/C$ individuals are immobile (see fig:2(e) and (f)).

In the DGEP model, susceptible individuals are diffusing. Therefore, currents of $A$ individuals are present. Conversely, the CGEP model exhibits currents of $C$ individuals.
The currents of $A$ and $C$ individuals in DGEP and CGEP, respectively, have the same orientation but reversed direction due to the inverted configuration of the densities.
Likewise, the impact on infected individuals resulting from the coupling with $A$ and $C$ is also reversed: Susceptible individuals transform into infected ones (${A + B \rightarrow 2B}$), whereas recovered individuals aid in the recovery of infected ones (${C + B \rightarrow 2C}$).
The overall effect of the two currents on the spreading of the infection might, however, be similar as there is compensation for the opposite direction of currents and opposite reactive effects.
In the next section, we demonstrate through an RG analysis that the DGEP and CGEP systems lie in the same universality class, which could be linked to the similar effects of the respective diffusive fluxes.

§ FIELD-THEORETICAL ANALYSIS

In the previous section, we discussed the spreading of infected individuals within the mean-field approximation of the stochastic epidemic model. We observed that, above a threshold for the initial density of susceptible individuals, the seed of infected individuals spreads throughout the entire system. Below this threshold, infected individuals initially spread, but eventually, the epidemic goes extinct, reaching the absorbing state.
A similar phenomenology, typical of epidemic models featuring an absorbing state phase transition (see for example [53]), persists when stochasticity is considered.
Stochastic effects become increasingly important near criticality when the system is close to the extinction threshold, and demographic noise cannot be neglected.
To analyze the stochastic critical dynamics, we use a field-theoretical description, also termed Fock space formulation, in which the stochastic lattice gas model given by eq:chemical.rates1–eq:diffusion.rates3 is formulated in terms of a master equation for the probability distribution of individuals' densities and then converted into a path integral representation [22, 60].
This enables the use of standard perturbative renormalization group tools to determine the universality class and critical exponents of the GEP, DGEP, and CGEP models [21].

§.§ Action of the field theories

The actions of the three analyzed models containing only relevant terms are derived in appendix A
and can be expressed as limiting cases of the following general action:

\begin{eqnarray}
\fl \mathcal{S}= \int_{\bi{q}}
\tilde{B}(-\bi{q}) \,
\bigg[
\rmi\omega+D_{B}k^{2}+\tau
\bigg] \,
B(\bi{q}) - \frac{U}{2} \int_{\bi{q}_1,\bi{q}_2} \!\!\! B(\bi{q}_1+\bi{q}_2)\tilde{B}(-\bi{q}_1)\tilde{B}(-\bi{q}_2) \nonumber
\\
\Bigg[ \frac{V \, B(\bi{q}_1)}{\rmi\omega_1+D k_1^2} + \frac{W \, B(\bi{q}_1)}{\rmi\omega_1} \Bigg]
\tilde{B}(-\bi{q}_1-\bi{q}_2) \, ,
\label{eq:action}
\end{eqnarray}

where we denote frequency $\omega$ and momentum $\bi{k}$ (${|\bi{k}|=k}$) by a 4-vector ${\bi{q}=(\omega, \bi{k})}$ and write the integral as ${\int_{\bi{q}}:=\int_{\mathbb{R}^{d+1}}   \rmd \omega \, \rmd^{d}\bi{k}} /(2 \pi)^{d+1}$ .
Here, ${B(\bi{q})= B(\omega, \bi{k})}$ denotes the Fourier transform of the density of infected individuals. In contrast, ${\tilde{B} (\bi{q}) = \tilde{B} (\omega, \bi{k})}$ is its dual response field; for an interpretation on its physical meaning see [21, 60].

The first term contains the inverse propagator for the field $B$, where the critical tuning parameter $\tau$ plays the role of an effective infection rate, encoding the dynamical interplay of infection and recovery.
At criticality, one expects that these processes balance out such that $\tau= 0$ and the propagator contains no intrinsic length scale.
From the Fock space formulation, one finds that $\tau=\zeta-A_0 \sigma$, which can be tuned to zero by adjusting the values of $\zeta$, $\sigma$ and $A_0$ such that the deterministic threshold is reached $A^{\rm mf }_c=A_0=\zeta / \sigma$, as found in section <ref>.
Note, however, that the actual critical value of $\tau$ is shifted with respect to ${\tau = 0}$ when stochasticity is considered.
The interaction term
proportional to $U$ originates from the demographic noise in the model and has the same functional form as in the DP and GEP actions [25, 28].
The two terms proportional to $V$ and $W$ arise from the non-linear contributions of the infection and caregiving reactions mediated by $A$ and $C$ individuals.
Notably, the details on how the healthy species ($A, C$) interact with the infected one ($B$) do not enter these two terms, but only whether the interacting species is diffusing or immobile.
Indeed, the two terms contain an explicit momentum contribution—${(\rmi\omega_1+D k_1^2)^{-1}}$ for $V$ and $(\rmi\omega_1)^{-1}$ for $W$—corresponding to the propagator of a diffusing field and a non-diffusing field, respectively. Which of the fields $A$ or $C$ corresponds to which propagators depends on the model. Accordingly, $D$ corresponds to $D_A$ (DGEP) or $D_C$ (CGEP) depending on the considered model. Table <ref> summarizes the values of $U$, $V$, $W$, and $D$ expressed in terms of the microscopic parameters for the different models.

For the GEP model, where only a coupling with a non-diffusing field ($A$) is present, $V$ equals zero, and $W$ describes the infection reaction.
In the DGEP model, where only a coupling with a diffusing field ($A$) is present, $W$ equals zero, and the infection reaction is described by $V$.
Finally, in the CGEP model, the infection reaction is described by $W$ (coupling with a non-diffusing species $A$), and the caregiving one is described by $V$ (coupling with a diffusing species $C$).
Note that in the limit of ${D=0}$, the term proportional to $V$ becomes identical to the one proportional to $W$, and thus, it describes a coupling with a non-diffusing field.
Following this reasoning, we recover the GEP action from the DGEP and CGEP models in the limit ${D=0}$.

Interaction vertices in terms of microscopic parameters. Values of the interaction vertices in the action eq:action in terms of physical parameters for the three models. For a derivation of those values, please see appendix A. Note that the parameter $D$ is not defined in the GEP model, as no other diffusing species other than the infected one is present.
GEP      DGEP      CGEP

$U$     $2 \sigma A_0$     $2 \sigma A_0$     $2 \sigma A_0$

$V$     $0$     $\sigma^2 A_0$      $\gamma \zeta$

$W$     $\sigma^2 A_0$     $0$     $\sigma^2 A_0$

$D$     $-$     $D_A$     $D_C$

§.§ Rapidity reversal symmetry

The action of the GEP model, (equation eq:action with ${V=0}$), is invariant under the following transformation [28]:

\begin{eqnarray}
B( \bi{q})
\rightarrow
-\frac{U}{W} \, \rmi\omega \,
\tilde{B}(\bi{q})
\, ,
\qquad
{\rm and}
\qquad
\tilde{B}(\bi{q})
\rightarrow
\, .
\label{eq:rapidity.reversal.symmetry}
\end{eqnarray}

This symmetry is known as rapidity reversal symmetry.
Generally, symmetries constrain the degrees of freedom of physical systems. In the context of criticality, this happens by reducing the number of independent critical exponents through a hyperscaling relation [23]

\begin{eqnarray}
\theta=\frac{d}{z}-2\delta-1, \label{eq:hyperscaling.relationII}
\end{eqnarray}

where $\theta$, $\delta$, and $z$ were defined in equations eq:observable1–eq:observable3 and $d$ is the spatial dimension of the system.
The number of independent exponents in the GEP universality class is reduced to three by the above hyperscaling relation.

Because of the vertex proportional to $V$, encoding the coupling of the field $B$ with a diffusive field, action eq:action violates the rapidity reversal symmetry, similarly as in [41]. Consequently, the critical exponents of models featuring a coupling with a diffusing species are not linked by the hyperscaling relation eq:hyperscaling.relationII and the four critical exponents appearing in eq:observable1–eq:observable3 are all independent.
This suggests a scenario in which at least one of the critical exponents deviates from the corresponding GEP, making DGEP and CGEP belong to a separate universality class.

§.§ Number of absorbing states

Another significant difference between coupling infected individuals with a diffusing and with a non-diffusing species lies in the number and type of absorbing states the system exhibits. This is an essential property because it is generally believed that the nature of the absorbing states has a decisive influence on the critical dynamics and the universality class of a model (see, e.g. [26, 22]).
A famous example is the DP conjecture, according to which all absorbing state phase transitions described by a single order parameter decoupled from any slow variable and with a single absorbing state belong to the same DP universality class [37, 25].

In the GEP model, where the infected individuals are coupled to a non-diffusing species only, the entire dynamics of the system are encoded in the movement of the infected individuals.
Once these entities are extinct (i.e., ${B(\bi{r},t)=0}$), the dynamics come to a halt, resulting in a perpetual frozen state for the system.
Consequently, each arrangement of the non-diffusing coupled fields represents a distinct absorbing state, reached when no remaining $B$ individuals are left. In the particular case of GEP, the distribution of the frozen-in cluster of individuals $A$ at criticality is described by the universality class of isotropic percolation [28].
The situation changes if the infected individuals are coupled to a mobile species instead, as in the DGEP model. In this case, the susceptible species can still diffuse freely after the extinction of the infected individuals and approach a homogeneous distribution in the limit ${t \rightarrow \infty}$. Therefore, in the presence of mobility, the absorbing state only exists in a state with homogeneous density in $A$, reached by diffusive currents. Note that the distribution of recovered individuals doesn't influence the dynamics of the infected population. Therefore, its configuration in the absorbing state can be disregarded.
The classification is more subtle in the CGEP model since it features absorbing states where the recovered individuals are mobile and reach a homogeneous distribution while the susceptible individuals are frozen in time.
We perform a perturbative renormalization group analysis to examine if and how these differences affect the critical behavior of the different models.

§.§ Renormalization group analysis and flow equations

The renormalization group
theory is applied to study how the parameters of a theory change as the typical length scale on which the system is observed is increased by a coarse-graning procedure [61, 62]. This scale dependence is encoded in a set of differential equations termed flow functions [63].
Scale-invariant behavior in equilibrium systems is typically observed at the critical points of second-order phase transitions or in the low-temperature phase of systems featuring Goldstone modes [21]. In non-equilibrium systems, scale invariance is obtained in the presence of Goldstone modes as for the KPZ equation [64] and—fairly generically—also in reaction-diffusion systems with absorbing states [23, 22]. In all these cases, the critical behavior is described by a fixed point of the flow equations, such that the theory's parameters remain invariant at different length scales. This fixed point is reached at large length scales and is therefore denoted an IR-stable fixed point.

Here, we utilize a perturbative field-theoretical ($Z$-factor)
RG approach to calculate these flow functions to one-loop order. The expansion parameter is $\epsilon = d_c - d$, the deviation from the upper critical dimension of the theory, which in our case is $d_c = 6$.
The $Z$-factors are determined through dimensional regularization in a minimal subtraction scheme (for a review of these methods, see [65, 21]).
To perform a perturbative RG study, it is crucial to identify the effective couplings, which are combinations of the theory's coupling parameters.
In appendix B, we show that one possible choice of effective couplings is given by

\begin{eqnarray}
\fl
\frac{UW}{(4 \pi D_B)^3} \, \mu^{-\epsilon}=:g_1 \, ,
\qquad
\frac{UV}{(4 \pi D_B)^3} \,  \mu^{-\epsilon}=:g_2 \, ,
\qquad {\rm and}
\quad
\frac{D}{D_B}=:g_D \, . \label{eq:CGEP.eff.coupling}
\end{eqnarray}

The renormalized momentum scale $\mu^{-\epsilon}$ is inserted to ensure that the effective couplings $g_1$ and$g_2$ are dimensionless at ${d=6-\varepsilon}$ dimensions. They encode the interplay between
the demographic noise (vertex $U$) and the non-linear reactions, which couples the infected species with a diffusing (vertex $V$) or non-diffusing species (vertex $W$). As previously discussed, the information on how the healthy species $A$ and $C$ are coupled to the infected one does not enter the effective couplings. It only matters whether the coupled species is diffusing (effective coupling $g_2$) or immobile (effective coupling $g_1$). Finally, the effective coupling $g_D$ describes the mobility of diffusing healthy individuals
relative to the mobility of infected ones.
In appendix B, a detailed derivation of the RG flow functions for the action eq:action is given. We ensure those flow functions describe the system dynamics at criticality by restricting the analysis to the critical plane $\tau=0$.
We find:

\begin{eqnarray}
\fl  \beta_{g_{1}}(g_{1},g_{2},g_{D})
:=\mu \frac{d}{d\mu}g_{1}=-\epsilon g_{1}+g_{1}\left(\frac{7g_{1}}{2}+g_{2}\frac{5+5g_{D}+2g_{D}^2}{2(1+g_{D})^3}\right) \, ,
\label{eq:flow.function.CGEP1}
\\
\fl  \beta_{g_{2}}(g_{1},g_{2},g_{D})
:=\mu \frac{d}{d\mu}g_{2}= -\epsilon g_{2}+g_{2}\left(\frac{5g_{1}}{2}+g_{2}\frac{5+5g_{D}+2g_{D}^2}{2(1+g_{D})^2}\right) \, ,
\label{eq:flow.function.CGEP2}
\\
\fl  \beta_{g_{D}}\left(g_{1},g_{2},g_{D}\right) :=\mu \frac{d}{d\mu}g_{D}=-g_{D}\left(\frac{7g_{1}}{12}+g_{2}\frac{7+4g_{D}+g_{D}^2}{12(1+g_{D})^3}\right) \, .
\label{eq:flow.function.CGEP3}
\end{eqnarray}

fig:3 illustrates the RG flow resulting from eq:flow.function.CGEP1–eq:flow.function.CGEP3 for dimensions below ${d_c = 6}$ (${\epsilon>0}$).
It exhibits three different fixed points: the Gaussian fixed line, the well-known GEP fixed point, and a new fixed point we term the DGEP fixed point.

Renormalization group flow.
A sketch of the RG flow of the effective couplings $g_i$ (solid black lines with arrows) obtained from the RG flow equations eq:flow.function.CGEP1–eq:flow.function.CGEP3 where the direction of the arrows denotes the flow towards larger length scales. Note that we consider only the physical region of the parameter space where all effective couplings are positive.
The flow exhibits a Gaussian fixed line  (brown shaded line) at ($g_1=g_2=0$), the GEP fixed point (purple dot) at (${g_1^*=2\epsilon/7, g_D^*=g_2^*=0}$) and the new DGEP fixed point (green dot) at (${g_1^*=0=g_D^*}$, ${g_2^*=2\epsilon/5 }$).
In GEP, the infected
individuals ($B$) are exclusively coupled to immobile susceptible individuals via the infection reaction, such that the model's RG
flow is restricted to the $g_D$-$g_1$ plane (gray-shaded plane). It approaches the GEP fixed point. Here, $g_D$ represents the mobility of an uncoupled species, and
the flow is indicated by a curved line as the dynamics of $g_1$ is independent of $g_D$, but the one of $g_D$ depends on $g_1$ (see eq:flow.function.CGEP1–eq:flow.function.CGEP3).
On the contrary, in DGEP, the infected individuals are only coupled to diffusing susceptible individuals via the infection reaction, restricting the model's RG flow to the $g_2$-$g_D$ plane (green shaded plane). The flow reaches the DGEP fixed point along this plane.
In CGEP, the infected individuals are coupled to both immobile susceptible individuals via the infection reaction, and mobile recovered individuals via the caregiving reaction. The RG flow starts from ($g_1\neq 0, \ g_2 \neq 0, \ g_D \neq 0$) (coordinates of the point indicated by gray dashed lines) and approaches the DGEP fixed point.

The Gaussian fixed line (${g_1=g_2=0}$) corresponds to a model where the noise amplitude is effectively zero ${U=0}$, i.e., the mean-field limit, or equivalently where both the infection and caregiving reactions are absent (${V=W=0}$), effectively decoupling the dynamics of the infected individuals $B$ from that of the healthy ones. Note that at the Gaussian fixed line (brown line in fig:3), ${g_2}$ is equal to zero, indicating that the effective coupling $g_D$ no longer enters in action eq:action, since it represents the mobility of a species which does not interact with the infected one, like the recovered individuals \(C\) in the CGEP model.

In GEP, where infected individuals only interact with a non-diffusing species, the RG flow starts from ${g_2=0}$ and is constrained to the $g_1$–$g_D$ plane (the gray plane in fig:3). The flow reaches the GEP fixed point (${g_1^*=2\epsilon/7,g^*_2=g^*_D=0}$) as identified by [28]. Since ${g_2=0}$, $g_D$ represents the mobility of a species whose dynamics are not coupled to that of \(B\) and any $g_D$-dependence vanishes in $\beta_{g_1}$. Therefore, $g_1$ is the only effective coupling in the GEP model, and its value approaches the GEP fixed point under the RG flow given by $\beta_{g_1}(g_1,0,0)$, which is the same flow function as found in [28].
The GEP fixed point is stable in the $g_1$ direction but unstable along the $g_2$ direction.
For finite $g_2$, i.e., when the infected individuals are coupled to a diffusing field, the RG flow approaches the new DGEP fixed point (${g_1^*=0=g_D^*}$, ${g_2^*=2\epsilon/5 }$). This fixed point is stable along all three directions $g_1$, $g_2$ and $g_D$. Indeed, both the flow corresponding to DGEP—described by $\beta_{g_2}(0,g_{2},g_{D})$ and $\beta_{g_D}(0,g_{2},g_{D})$—and CGEP—which requires all three effective couplings—approach this new IR-stable fixed point.
Thus, both models belong to the same universality class, different from GEP, whose critical exponents are given in table <ref>.
This universality class describes the critical dynamics close to the extinction threshold of models where the infected individuals are coupled to a diffusing species, as in the DGEP and CGEP models.
The surprising fact is that, at criticality, the interactions with a non-diffusing species are entirely suppressed by the coupling with a diffusive field since at the DGEP fixed point ${g_1=0}$.

In summary, the universality class of the critical dynamics is altered if the infected population is coupled to at least one diffusing species, suppressing all other interactions with non-diffusing fields. Previous studies found similar results [66], suggesting that the same principle can be valid in other generic epidemic models where the coupling between a diffusive and a non-diffusive species is present.
Critical exponents. The value of the critical exponents eq:observable1–eq:observable3 of the GEP and DGEP universality classes calculated both with an $\epsilon$-expansion below the upper critical dimension $d_c=6$ (theory) and through stochastic simulations of the GEP (sim. GEP), DGEP (sim. DGEP) and CGEP (sim. CGEP) models in dimension $d=2$.  The obtained exponents are compared with the numerical [67] and theoretical [28] literature values of the GEP universality class.
3cGEP         3cDGEP

theory     lit. [67]     sim. GEP         theory     sim. DGEP     sim. CGEP

$z$     $2-\frac{1}{6}\epsilon$     $1.1295$     $1.13 \pm 0.004$       $2-\frac{7}{30}\epsilon$     $1.14 \pm 0.03$     $1.128 \pm 0.004$

$\delta$     $1-\frac{5}{28}\epsilon$     $0.092$     $0.093 \pm 0.001$        $1-\frac{1}{4}\epsilon$     $0.091 \pm 0.009$     $0.093 \pm 0.003$

$\theta$     $\frac{3}{28}\epsilon$     $0.586$      $0.582 \pm 0.003$        $\frac{3}{20}\epsilon$     $0.54 \pm 0.02$     $0.554 \pm 0.005$

$v_{\parallel}$     $1+\frac{1}{28}\epsilon$     $1.506$      $1.5 \pm 0.1$         $1+\frac{1}{20}\epsilon$     $1.4 \pm 0.1$     $1.4 \pm 0.1$

§.§ Vanishing diffusivity implies additional divergences

At first sight, it may appear counterintuitive that at the DGEP fixed point, the effective coupling $g_D$ is equal to zero since this is the value it takes when the coupled species is immobile.
To clarify this, we look at the limit ${D\rightarrow 0}$, which helps to elucidate how this limit leads to a different universality class and why $g_D$ is zero at both the GEP and DGEP fixed point.

In the minimal subtraction scheme used in this work, the divergences arising in the perturbative calculation of the $n$-point correlation functions completely determine the flow functions and, thus, the fixed point structure of the model. In appendix B, we show that
almost all divergences arising from the vertices $V$ are $W$ are symmetric, meaning that divergences
proportional to $g_1$ can be obtained from the ones proportional to $g_2$ in the limit ${D\rightarrow 0}$.
The only exception stems from a single diagram (diagram (A), in appendix B), which turns out to have a divergent contribution only proportional to $g_1$ and not $g_2$. This divergence is responsible for the asymmetry between the flow functions $\beta_{g_1}$ and $\beta_{g_2}$ eq:flow.function.CGEP1–eq:flow.function.CGEP2.
In appendix B we evaluate this diagram, which results in two different contribution: $g_2 V D^3_B  I(D)+g_1 W D^3_B I(D=0)$, where:

\begin{eqnarray}
\fl I(D):= \int
\frac{\rmd^d \bi{p}}{(2 \pi)^d}
\frac{1}{(D-D_B) p^2 +\tau} \, \frac{1}{\big[\rmi\omega_1 +(D_B+D)p^2 +\tau\big]^2} \, \frac{1}{\rmi\omega_1 +2 D p^2}
\, .  \label{eq:anonamlous.divergent .diagram}
\end{eqnarray}

The integral over the momentum $ \bi{p}$ is evaluated at dimension ${d=6-\epsilon}$, as dimension regularization requires. The key point is that only $I(D=0)$ is divergent, not $I(D\neq 0)$. The reason can be inferred by simple power counting. If ${D\neq0}$, the denominator scales for large momentum as $p^{8}$ and the integral converges. If ${D=0}$, the scaling changes to $p^{6}$ and the integral diverges.
Consequently, the only divergent contribution of diagram (A), calculated through dimensional regularization, is given by ${g_1 W D^3_B I(0)=g_1 W/(\rmi\omega_1 \epsilon)+\mathcal{O}(\epsilon^0)}$.
This divergence has the same structural form as the vertex $W$ in action eq:action; thus, it
has to be reabsorbed by the $Z$-factor of the vertex $W$ in the renormalization procedure. The implication is that the $Z$-factor of the vertex $W$
and $V$ differ, leading to different flow functions for the effective couplings $g_1$ and $g_2$. In particular, the additional divergence arising from diagram (A) shifts the $7g^2_1 /2$ in eq:flow.function.CGEP2 to $5g^2_2 /2$ in eq:flow.function.CGEP1 and creates an asymmetry. The asymmetry in the flow functions of $g_1$ and $g_2$ is responsible for changing the value of the GEP fixed point to
(${g_1^*=2\epsilon/7,g^*_2=g^*_D=0}$) as compared to the DGEP one  (${g_2^*=2\epsilon/5,g^*_1=g^*_D=0}$), leading to a different universality class.

We believe that a similar mechanism by which some Feynman diagrams acquire additional divergences in the limit of no diffusivity can equally explain why many models featuring absorbing state phase transitions change universality class depending on the diffusivity of the various species. One example is the DEP model with immobile healthy species, which lies in a different universality class than the diffusive DEP model [68, 69].

Having explained how setting ${D = 0}$ leads to an additional divergence, we can now resolve the seeming ambiguity that the effective coupling ${g_D=D / D_B}$ vanishes at the DGEP fixed point.
The important difference between the GEP and DGEP fixed points is how $g_D$ reaches its zero value.
When infected individuals are coupled to a diffusing species, like in DGEP, $g_D$ approaches zero by $D_B$ flowing to infinity under the RG flow. In contrast, $D$ does not flow and remains fixed at its non-zero initial value.
On the other hand, if $B$ is coupled to an immobile species, like in the GEP model, the effective coupling $g_D$ is zero due to the value of the diffusion constant ${D=0}$.
However, we showed before that the flow functions acquire an additional contribution only if ${D=0}$, independent of the value of $D_B$. Therefore, even if $g_D = 0$ at the DGEP fixed point, we have $D\neq0$, thus guaranteeing that the DGEP flow functions differ from the GEP ones.

§ STOCHASTIC SIMULATIONS OF EPIDEMIC MODELS

In the preceding section, we discerned a novel universality class through an RG analysis by identifying a relevant modification of the GEP model.
However, since we are using a perturbative method expanding around the upper critical dimension ${d_c=6}$, the quantitative predictions concerning the values of the critical exponents may not be accurate in physical dimensions ${d\leq 3}$.
Therefore, to determine the values of the critical exponents, we carried out stochastic simulations.
Specifically, we used a Gillespie algorithm based on an algorithm published in previous work [53] to study the critical dynamics in the vicinity of the absorbing state phase transition for two-dimensional lattice models as defined in section <ref>.

§.§ Gillespie algorithm for spatially extended systems

Most previous studies implement non-diffusive processes on a regular lattice with an exclusion constraint and simulate the dynamics by standard Monte Carlo updates [23, 27, 67, 70, 24, 71, 72], where two-individual reactions are implemented as nearest-neighbor interactions [73, 27].
The exclusion constraint considers steric interactions between the individuals and allows, at most, a single occupancy for each lattice site.
With this constraint, the mobility of individual agents can only be implemented by pairwise exchange processes [41]. This, however, implies a symmetry in the mobilities of different species and, thus, disallows exploring asymmetric diffusion constants. Especially, the interesting limiting cases with certain immobile species cannot be implemented without further effort. Beyond that, the above class of Monte Carlo algorithms does not bear an intrinsic notion of time.
Rather, time is associated with the number of updating steps.
These limitations can be overcome by resorting to a different class of stochastic updating schemes.

Here, we employ a modified version of the stochastic simulation algorithm (SSA), commonly referred to as  Gillespie algorithm [74, 75]. In this approach, the underlying lattice structure of the models is retained while the exclusion constraint is abandoned. Thus, multiple individuals can occupy a single site.
Two-individual reactions occur on a single site, and diffusion is implemented as nearest-neighbor hopping. This decouples the individual species' mobility and allows their hopping rates to differ. In particular, to account for the complexity of the stochastic many-body dynamics, we resort to the hybrid next reaction scheme [76, 53] which, in principle, allows to solve the problem in sub-polynomial time.
The benefit of this choice is that it generates exact stochastic trajectories with a clear notion of time and does not require any exclusion constraint. We are specifically building on an approach that has recently proven successful in resolving a longstanding debate regarding the critical dynamics of the DEP model [53]. While the hybrid next reaction scheme provides certain benefits, as outlined above, it also requires substantially more computational resources than standard Monte Carlo updates on a lattice with exclusion constraint.

§.§ Initial conditions and units.

As outlined in section <ref>, we study the spreading dynamics starting from a single seed of infected individuals. To that end, we initialize a square lattice with a homogeneous density of $A$ individuals and an infectious seed of $B$ individuals placed on the central site.
For all of our simulations, the initial density of $A$ individuals is one individual per lattice site, and the infectious seed consists of ten $B$ individuals. Additionally, finite-size effects can be largely eliminated with single seed initial conditions [23]. Finite-size effects arise mostly from the interaction of $B$ individuals with the boundary of a finite lattice. In all the simulations, we ensure that no such finite-size effects occur by choosing a large enough system size $L$, such that no $B$ individual ever reaches the boundary during the finite simulation time. Note that this is only made possible by the algorithm's sub-linear scaling with system size. Lastly, time is measured in units of $10$ times the recovery time, setting ${\zeta = 0.1}$.
Length scales are measured in multiples of the lattice spacing ${a = 1}$.

§.§ GEP model as a benchmark

Critical exponents of GEP model.
(a) Averaged local slope for the number of infected individuals for a set of values of the infection rate $\sigma$ as indicated in the legend. Based on these results, we estimate the critical point to lie at $\sigma_c = 1.3765$ (between the two lines at $\sigma=1.376$ and $\sigma=1.377$) and the value of $\theta=0.582 \pm 0.003$ as indicated by the circular marker and error bars.
The triangular marker indicates the current best literature estimate for $\theta = 0.585$ [23].
For each value of $\sigma$, results were obtained by averaging over 24,000 stochastic trajectories starting from a single seed at the origin.
For each realization, the dynamics were monitored for 50,000-time units in a 2D system with linear expansion of 10,000 lattice sites. Other parameters were chosen as ${D_A=D_C=0}$, ${D_B=0.1}$, ${\zeta=0.1}$, and ${\gamma=0}$.
(b)–(d) Double logarithmic plots for the key observables of the spreading dynamics for a range of infection rates indicated in the graph:
(b) survival probability $P_{\rm s} (t)$,
(c) mean number of infected individuals $\langle N \rangle (t)$, and (d) mean square radius of the infected cluster $\langle R^2 \rangle (t)$.

To validate the chosen algorithm and to establish a benchmark for the critical exponents and their errors, we first reinvestigate the GEP model to test how well our code reproduces the GEP critical exponents from the literature [67, 71, 72, 41].
Figures <ref>(b)–(d) show double logarithmic plots of the three key observables defined in equations eq:observable1–eq:observable3 of the spreading dynamics for a range of values of the infection rate $\sigma$.
Inferring the critical value of $\sigma$ from these plots amounts to finding the curves that do not bend upwards or downwards, i.e., are straight and thus represent power laws.
Given the proximity of the curves, achieving high accuracy proves difficult.
Therefore, a common approach is to use the averaged local slope [73, 70]. It is defined as follows for the mean number $\langle N\rangle(t)$ of infected individuals and a local approximation $\theta(t)$ of the corresponding critical exponent,

\begin{eqnarray}
\theta(t)
= \frac{\ln[\langle N \rangle(t) ] - \ln[\langle N \rangle (t/\Delta t) ]}{\ln[\Delta t]}
\,  ,
\label{eq:local.slope.theta}
\end{eqnarray}

and analogously for all other observables and exponents.
Here, the time interval $\Delta t$ must fulfill the condition ${\Delta t > 1}$. <ref> elaborates on our choice of ${\Delta t=10}$, which is comparable to previous studies [70, 18, 34, 77, 40].

The intercept with the ${1/t=0}$ axis—which signifies the measured critical exponent—is estimated by a linear extrapolation of $\theta(t)$ for our best estimate of the critical infection rate [23, 70]. It corresponds to the line of the infection rate, for which a linear extrapolation is possible and thus results in a finite intercept.
Figure <ref>(a) shows the averaged local slope for the mean number of infected individuals. In appendix E, we explain how we determined the error of the exponents.
Comparing our results for the critical exponents with literature values [67, 71, 72, 41] we find excellent agreement, cf. table <ref>.
The larger size of our error bars is due to using a different simulation algorithm [53] combined with our cautious and conservative approach to error estimation as outlined in appendix E.
Our simulation algorithm convincingly reproduces well-established results from prior studies [67, 71, 72, 41], providing a robust benchmark for evaluating the outcomes of the DGEP and CGEP models in the subsequent section.

§.§ Simulation results for the DGEP and CGEP models

Critical behavior of the DGEP and CGEP models.
(a) Averaged local slope for the number of infected individuals for the DGEP model for a set of values for the infection rate $\sigma$ indicated below the graph.
Based on these results, we estimate the critical point to lie at ${\sigma_c = 0.93}$ and the value of ${\theta=0.54 \pm 0.02}$ as indicated by the circular marker with error bars. The triangle shows our estimate for GEP for comparison, where we omitted the error bars because they are smaller than the marker.
For each value of $\sigma$, results were obtained by averaging over $20,000$ stochastic trajectories starting from a single seed at the origin. The dynamics were monitored for $10,000$ time units in a 2D system with linear expansion of $2,600$ lattice sites for each realization. Other parameters were chosen as ${D_C=0}$, ${D_A=D_B=0.1}$, ${\zeta=0.1}$, and ${\gamma=0}$.
(b) Averaged local slope for the number of infected individuals for the CGEP model for a set of values for the infection rate $\sigma$ indicated below the graph.
Based on these results, we estimate the critical point to lie at ${\sigma_c = 0.43537}$ and the value of ${\theta=0.554 \pm 0.005}$ as indicated by the circular marker. A triangular marker represents our estimate for GEP. Both results' errors are smaller than the respective marker.
For the critical value $\sigma_c$, results were obtained by averaging over $34,000$ stochastic trajectories starting from a single seed at the origin. The dynamics were monitored for $40,000$ time units in a 2D system with linear expansion of $8,800$ lattice sites for each realization.
For each other value of $\sigma$, results were obtained by averaging over 34,000 stochastic trajectories starting from a single seed at the origin. The dynamics were monitored for $20,000$ time units in a 2D system with linear expansion of $4,400$ lattice sites for each realization.
In both cases, the other parameters were chosen as ${D_A=0}$, ${D_B=D_C=0.1}$, ${\zeta=0.05}$, and the rate ${\gamma=0.1}$ for the new caregiving reaction.

Data collapse for the DGEP and CGEP models. (a) Rescaled number of infected individuals in the DGEP simulations. The data was collapsed using the critical ${\sigma_c=0.93}$, ${\theta=0.54}$, and ${\nu_\parallel=1.4}$. (b) Rescaled number of infected individuals in the CGEP simulations. The data was collapsed using the critical ${\sigma_c=0.43537}$, ${\theta=0.55}$, and ${\nu_\parallel=1.4}$. The underlying data and simulation parameters are identical to the data shown in figure <ref>.

After benchmarking our stochastic simulation algorithm using the GEP model, we next performed stochastic simulations of the spreading dynamics of the DGEP and CGEP models.
The corresponding results for the averaged local slope of the number of infected individuals are shown in figure <ref>. For the figures of the other observables, we refer the reader to appendix F.
Using the methods outlined above and in appendix E, we determined the critical exponents $\theta$, $\delta$ and $z$ of the three key observables defined in equations eq:observable1–eq:observable3.
They are summarized in table <ref> together with their error estimates.

There are two key observations:
Our numerical estimates of the exponents $z$ and $\delta$—the critical exponents related to the temporal scaling of the mean squared radius and the survival probability—match their respective values of our GEP simulations within error bars as well as those of prior work [67, 71, 72, 41].
However, the numerical estimate of the critical exponent $\theta$—related to the temporal scaling of the mean number of infected individuals—of the DGEP ($\theta = 0.54 \pm 0.02$)  and CGEP ($0.554 \pm 0.005$) models distinctly differ from the corresponding GEP ($0.582 \pm 0.003$) value. Importantly, they lie well outside of the error bars.
Secondly, the numerical estimates of the three critical exponents for DGEP and CGEP are identical within error bars.
Taken together, the numerical results affirm our analytical findings from the RG analysis that the critical dynamics of the DGEP and CGEP models belong to the same universality class, which differs from the GEP universality class.

Remarkably, we find very good agreement between our RG predictions of the exponents in dimension ${d=2}$ (i.e., ${\epsilon = 4}$) and our numerical estimates.
The analytic prediction for $\theta$ of ${\theta_{\rm theory} = 0.6}$ deviates less than $10 \%$ from the numerical estimates of ${\theta_{\rm num} = 0.54 \pm 0.02}$ (DGEP) and ${\theta_{\rm num} = 0.554 \pm 0.005}$ (CGEP).
This is also the case for the dynamical exponent $z$.
Here, the RG calculation predicts ${z_{\rm theory} = 1.0\overline{6}}$ compared to the numerical estimates of ${z_{\rm num} = 1.14 \pm 0.03}$ (DGEP) and ${z_{\rm num} = 1.128 \pm 0.004}$ (CGEP).
For $\delta$, the analytical prediction is ${\delta_{\rm theory} = 0}$, which deviates from the numerical estimates of ${\delta_{\rm num} = 0.091 \pm 0.009}$ (DGEP) and ${\delta_{\rm num} = 0.093 \pm 0.003}$ (CGEP). Nevertheless, the absolute deviation is still minimal considering the value of ${\epsilon =4}$.
However, our simulations can neither confirm different values for the $\delta$ and $z$ exponents of the DGEP universality class when compared to the GEP one nor the analytic prediction of ${\theta_{\rm DGEP} > \theta_{\rm GEP}}$.

Moreover, using our simulation results, we confirm that the hyperscaling relation eq:hyperscaling.relationII is indeed violated for the CGEP model:

\begin{eqnarray}
0.554 \pm 0.005 = \theta \neq \frac{d}{z}-2\delta-1 = 0.59 \pm 0.02 \, .
\label{eq:broken.hyperscaling.relation.CGEP}
\end{eqnarray}

Conversely, for the DGEP model, our carefully estimated errors do not permit a definite conclusion.
The absence of a rapidity reversal symmetry eq:rapidity.reversal.symmetry enforces the need to infer one additional independent critical exponent to fully characterize the universality class [23].
Figure <ref> illustrates the estimation of the independent critical exponent $\nu_\parallel$ via a data collapse [21, 24] for both models.
The specifics of the methodology and the error estimation are discussed in appendix E.
This analysis confirms the analytical result that the DGEP and CGEP models exhibit identical values for the $\nu_\parallel$ exponent as indicated in table <ref>.
The RG calculation in two dimensions predicts ${\nu_{\parallel, {\rm theory}} = 1.2}$, while the numerical estimate for both models is ${\nu_{\parallel, {\rm num}} = 1.4 \pm 0.1}$.
This is the largest absolute deviation between the two methods, though it is still smaller than $20\%$ in relative terms.

Moreover, our results show excellent agreement with the ones recently reported by Deng and Odor [41] for a version of the GEP model with symmetric exchange hopping and an exclusion constraint, which is a different microscopic implementation of the DGEP model.
Generally reporting smaller error bars, even a significant difference in the exponent $\delta$ could be identified. Our algorithm decouples the mobility of the individual species. With that, we could reinforce our analytical claim that the new DGEP universality class is also robust against keeping $C$ individuals immobile for the GEP reaction scheme with diffusing $A$ species.
Taking these results obtained on a system with exclusion constraint together with ours without such, one may faithfully claim the existence of a new universality class that deviates from the well-known GEP universality class. Surprisingly, the values of the exponents of the new DGEP universality class are very similar to those of the GEP universality class. This effect occurs both in the RG calculations as well as the stochastic simulations in two dimensions.

§ DISCUSSION

In this paper, we explored the impact of mobility on the critical dynamics in epidemics through a combined approach of renormalization group analysis and stochastic simulations.
Specifically, we focused on epidemics with immunization and considered two modifications—called DGEP and CGEP—of the general epidemic process (GEP): DGEP generalizes the GEP model by considering mobile susceptible individuals $A$, while in the CGEP model a new “caregiving” reaction (${B+C \rightarrow 2C}$) is introduced, where mobile recovered individuals ($C$) can aid the recovery of infected individuals ($B$).
In both models, the infected individuals are coupled to a diffusing species via a non-linear reaction. Diffusivity of the coupled species disrupts a symmetry of the GEP model, the so-called rapidity reversal symmetry.
The universality class of the absorbing state phase transition of the respective models was studied by a one-loop perturbative RG analysis in $d=6-\epsilon$ dimensions, just below the upper critical dimension. We found that the flow functions of both the CGEP and DGEP models reach the same new RG fixed point, indicating that they belong to a new universality class different from the GEP one.
Stochastic simulations affirmed these findings, which identified the same set of critical exponents for both systems.
The critical exponents' values do not fulfill the hyperscaling relation associated with the rapidity reversal symmetry, providing evidence that this symmetry is violated in the new universality class. This new universality class describes the spreading of an epidemic near the extinction threshold, where infected individuals interact with a diffusing species and gain immunization upon recovery.

Based on these analytical and numerical studies, we draw the following three main conclusions.
Firstly, the emergence of a new universality class shows that mobility plays a crucial role in the critical dynamics of epidemic models,
or, more generally, of non-equilibrium phase transitions.
This highlights the significance of exploring this phenomenon [40].
Diffusivity is important in other paradigmatic epidemic models for which a complete analytical understanding is still lacking [53].
Secondly, in the CGEP model—which couples infected individuals to both a diffusing and a non-diffusing species—the RG flow approaches a fixed point where the coupling to the non-diffusing species vanishes, i.e., the coupling to immobile species becomes increasingly irrelevant approaching the critical regime.
We hypothesize that this is one example of a more general effect: On large length and time scales and close to criticality, couplings to mobile species dominate the couplings to immobile ones.
Lastly, we showed that DGEP and CGEP exhibit the same critical behavior, indicating that, though sensitive to the introduction of mobility, universality classes might be surprisingly indifferent to how and with which diffusing species the order parameter is coupled.
This claim warrants a broader investigation of different ways mobility might be introduced to other non-equilibrium phase transitions.
In our models, one could test the robustness of this effect by adding new reactions, for instance, a vaccination reaction (${A+C \rightarrow 2C}$), coupling the susceptible species $A$ with the recovered $C$, and verify if the critical dynamics is again only determined by the mobility of the healthy species, or if the specific manner in which a healthy species is coupled to the infected one matters.

The existence of the new DGEP universality class was suggested by the violation of rapidity-reversal symmetry and by the different numbers of absorbing states between GEP and DGEP.
In particular, the GEP model features an infinite number of absorbing states, corresponding to arbitrary arrangements of frozen-in susceptible individuals after the disease has gone extinct.
In contrast, DGEP exhibits absorbing states of uniformly distributed susceptible individuals $A$, reached by diffusive currents.
The absorbing states in the CGEP model encompass both immobile and mobile species, making classification based on this criterion challenging.
However, having shown that the interaction with the immobile species becomes irrelevant at criticality, we can disregard the spatial distribution of the immobile species in the absorbing state and correctly classify the CGEP phase transition into the DGEP universality class.

To understand why adding mobility to the GEP model affects its critical behavior, we compare our result with previous works
[22], which studied the density decay of a two-species annihilation process (${A+B \rightarrow \emptyset}$).
They found that the corresponding decay rate $\lambda$ acquires non-trivial RG corrections when both $A$ and $B$ are mobile.
As the annihilation reaction is local, the RG corrections can be understood as an altered meeting probability for two diffusing individuals due to stochasticity.
We believe that similar considerations can be applied to our models. In particular, through a mean-field analysis (section <ref>), we showed that diffusion creates diffusive currents of healthy individuals.
When stochasticity is included, the effect of these currents can be understood analogously as a change in the probability of encounters between infected and healthy species, which leads to the new DGEP universal class.
We believe that the above explanation in terms of meeting probabilities applies generically, explaining why diffusion can alter the critical behavior and highlighting the robustness of the universality class featuring diffusive individuals.

Finally, our numerical methods, based on the hybrid next reaction scheme version of the SSA [76, 74, 75, 53], may open avenues for further numerical studies, enhancing our understanding of the role of diffusion in absorbing state phase transitions and non-equilibrium phase transitions in general.
The stochastic simulations of the DGEP model highlight the importance of decoupling the mobility of different species to allow for a more general analysis of the impact of mobility (see [53]).
The crucial step is to abandon the exclusion constraint in the lattice structure. This also aligns with the standard derivation of a field-theoretical
action in the path-integral formalism that strongly relies on bosonic density fields [21].
To conclude, we suggest that our stochastic simulation algorithm holds the potential for exploring the impact of diffusion in a broader range of problems beyond critical dynamics.
In particular, it could be used to analyze models with anomalous diffusion, for example, extracting anomalous diffusion exponents [78] and calculating higher order moments [79].

From a broader perspective, our study highlights the importance of mobility for the spreading dynamics of populations close to the extinction threshold. It suggests general mechanisms that may lead to new critical dynamics. Future studies may focus on different types of population dynamics [80, 81] and non-local interactions between species [33, 82, 83, 84].
We also see future potential in investigating the role of different types of individuals mobility like anomalous diffusion realized via Lévy flights [85, 86] and long-range dispersal [87, 88].

We thank Richard Dengler for inspiring discussions and Patrick Wilke and Moritz Striebel for their valuable support in getting started with the numerical code.
We acknowledge financial support by the German Research Foundation (DFG) through SFB 1032 (Project ID No. 201269156), the Excellence Cluster ORIGINS under Germany’s Excellence Strategy (EXC-2094-390783311), the European Union (ERC, CellGeom, project number 101097810), and the Chan-Zuckerberg Initiative (CZI).

§ DERIVATION OF THE ACTION

In this appendix, we derive a field-theoretical description for our models starting from the general reaction scheme eq:chemical.rates1–eq:chemical.rates3.
In section <ref>, we discussed the absorbing state phase transition within the mean-field approximation, where the effect of noise was disregarded.
To account for stochasticity, a field-theoretical analysis is imperative and is developed throughout this section.
We derive a field-theoretical action describing infected individuals' dynamics close to the extinction threshold.

§.§ The Coherent State Path
Integral method.

We start by reviewing how to map the reaction scheme eq:chemical.rates1–eq:diffusion.rates3 to a field-theoretical description for the well-mixed case and then generalize to a spatially extended system.
Given the probabilistic nature of the stochastic process, we aim to find the time evolution of the probability ${p(n_A,n_B,n_C,t)=: p(\bi{n},t)}$ that there are $n$ individuals for each given species.
A common way to find it is through the Coherent State Path Integral (CSPI) method, a specific type of Fock space formulation [60].
In this method, the probability distribution at time $t$ is marginalized over a Poissonian initial condition $p(\bi{n},t|\bi{x}_0, t_0)$, with initial mean ${\bi{x}_0=(x_{{0}_A},x_{{0}_B},x_{{0}_C})}$.
Probability distributions other than a Poissonian distribution could be used for marginalization, but this choice often is the most convenient in calculations [60].
The time evolution of the densities $\bi{n}(t)$, given by a backward master equation, can be mapped to the time evolution of the Poissonian means $\bi{x}(t)$ with initial conditions $\bi{x}_{0}$.
The different $x_i(t)$ capture the averaged statistical behavior of the density of each species $n_i(t)$. Following the standard procedure [60], we convert the backward master equation for $\bi{x}(t)$ into a path integral representation for $p(\bi{n},t|\bi{x}_0, t_0)$:

\begin{eqnarray}
p(\bi{n},t|\bi{x}_0,t_0) \propto \int \hspace{1mm}\mathcal{D}[\bi{x}(t'), \bi{q}(t')]
\exp\left(-\int_{t_0}^{t}\rm{d}t\mathcal{L}[\bi{x}(t'),\bi{q}(t')]\right) \label{eq:marginal.probability}
\end{eqnarray}

The marginalized probability distribution $p(\bi{n},t|\bi{x}_0,t_0)$ is expressed as a integral over all “paths" of $\bi{x}(t)$ and its response field $\bi{q}(t)$. Those paths are weighted with a stochastic Lagrangian density $\mathcal{L}$ of the form [21]

\begin{eqnarray}
\mathcal{L}[\bi{x}(t),\bi{q}(t)]= \rmi \bi{q}\partial_t \bi{x}(t)-\mathcal{Q}(\bi{x}(t),\rmi\bi{q}(t)+1) ,\label{eq:shifted:lagrangian}
\end{eqnarray}

where $\rmi \bi{q}\partial_t \bi{x}(t)$ captures the time evolution of the mean $\bi{x}(t)$ and the Liouville operator $\mathcal{Q}$ encodes the dependence on the specific reaction scheme. For an illustration, we look at a simple reaction involving three species $A$, $B$, and $C$:

\begin{eqnarray}
kB+mA  \overset{\lambda}{\rightarrow} lC.
\label{eq:reaction}
\end{eqnarray}

Calling $A(t)$, $B(t)$, $C(t)$ the time evolution of the Poissonian means for individuals $A$, $B$, $C$ and $\tilde{A}(t)$, $\tilde{B}(t)$, $\tilde{C}(t)$ their response fields, the Liouville operator for reaction eq:reaction reads

\begin{eqnarray}
\mathcal{Q}= \lambda B^{k}A^{m}(\tilde{C}(t)^{l}-\tilde{B}(t)^{k}\tilde{A}(t)^{m}) \label{eq:standardLiouville}
\end{eqnarray}

When studying a multi-reaction system, the contributions of the different reactions get summed up, as the Liouville operator is linear.

The above procedure can be generalized to a spatially extended system, where the dependence on space can be included by allowing individuals to hop on a lattice with lattice constant $a$.
To be able to use eq:standardLiouville, individuals' hopping processes between different sites must be formulated as reactions.
This is done by considering individuals on two different sites to belong to different species, coupled through the reaction $B_i {\rightarrow} B_j$ with the hopping rate specified by $D_B$.
For a sufficiently large system, a continuous description can be achieved by taking the limit $a\rightarrow 0$, resulting in the Poissonian means $\bi{x}(\bi{r},t)$  which are now additionally space-dependent.
Explicitly for our system, we rename $\bi{x}(\bi{r},t)$ in terms of the different species—e.g. $x_A(\bi{r},t)=A(\bi{r},t)$ and $q_A(\bi{r},t)=\tilde{A}(\bi{r},t)$. Using the standard form eq:standardLiouville for $\mathcal{Q}$ for the general reaction scheme eq:chemical.rates1–eq:diffusion.rates3, we get the following Lagrangian in continuous space:

\begin{eqnarray}
\fl \mathcal{L}= \tilde{A}(\bi{r},t) \big[\partial_{t}-D_{A}\nabla^{2}\big]
+ \tilde{B}
\big[\partial_{t}-D_{B}\nabla^{2}\big]
B(\bi{r},t) \label{eq:action.from.coherent} \\
+ \tilde{C}(\bi{r},t)
\big[ \partial_{t}-D_{C}\nabla^{2} \big]
C(\bi{r},t) -\zeta B(\bi{r},t)
\Big[ \tilde{C}(\bi{r},t)-\tilde{B}(\bi{r},t)\Big] \nonumber  \\
- \nonumber
\sigma A(\bi{r},t)B(\bi{r},t)\Big[ \big( \tilde{B}(\bi{r},t)+1 \big)^2-\big(\tilde{B}(\bi{r},t)+1 \big) \big(\tilde{A}(\bi{r},t)+1 \big) \Big] \\
-  \gamma \, C(\bi{r},t)B(\bi{r},t) \Big[ \big(\tilde{C}(\bi{r},t)+1\big)^2 -\big(\tilde{B}(\bi{r},t)+1\big)\big(\tilde{C}(\bi{r},t)+1\big)
\Big] .
\nonumber
\end{eqnarray}

The first three terms, containing the time and spatial derivatives, correspond to the propagators for the fields $A$, $B$, and $C$.
The Laplace operator originates from nearest-neighbor hopping.
The terms proportional to $\sigma$, $\zeta$, and $\gamma$ correspond to the on-site reactions of infection, recovery, and caregiving.
Those terms arise from the Liouville operators corresponding to each reaction.
In contrast to the general form for the Liouville operator eq:standardLiouville, all the response fields are shifted by 1, e.g. ${\tilde{A}\rightarrow \tilde{A}+1}$, in accordance to eq:shifted:lagrangian. This shift is known as the “Doi shift". The same kind of dynamic action functional can also be derived via other formalisms, for example, the Martin-Siggia-Rose-Bausch-Janssen-Wagner-de Dominicis response functional [89, 90, 91, 92, 60].
The strength of the CSPI approach is that it does not require any assumptions; all the information is contained in the specific form of the reaction scheme.
However, CSPI characterizes the evolution of fields that represent the statistically averaged behavior of a density (specifically, the Poissonian means of its distribution) rather than the density itself.
Consequently, some of the system's stochasticity is inherently blurred within this formalism and cannot be captured. An instance of this is “diffusional noise", also known as conserved noise.
In a Langevin description, this usually corresponds to a term of the type $\nabla( \sqrt{\rho} \eta)$, where $\rho$ is the density and $\eta$ a white Gaussian noise. To obtain this term in a Fock Space formulation, a description in terms of densities must be recovered by the so-called Cole-Hopf transformation, converting the CSPI action into a Kramers-Moyale type action [60].
Although diffusional noise cannot be derived within the CSPI formalism, its absence is inconsequential in our context as it is RG irrelevant at our system's upper critical dimension ${d_c=6}$; thus, it can be neglected.

§.§ Collecting  relevant terms.

Close to the extinction threshold, the dynamics predicted by the mean-field equations eq:deterministicequations1–eq:deterministicequations3 become invalid as fluctuations grow larger in the vicinity of an absorbing state phase transition.
We aim to calculate the corrections arising from stochasticity with a perturbative renormalization group method.
The first step is to expand the action around the deterministic critical values for the fields. In particular, we shift the fields
${A(\bi{r},t)\rightarrow A_0 + A(\bi{r},t)}$, and tune $A_0$ towards the critical $A_c^{\rm{mf}}$ (see section <ref> ).
This way, the shifted fields represent the deviations from the deterministic critical solution due to noise.
Via a dimensional analysis, one can determine the naive scaling dimension of the different interactions in eq:chemical.rates1–eq:chemical.rates3.
Thus, one can determine which interactions are relevant near criticality and which are irrelevant, i.e., their reaction strength goes to zero under the RG flow and can be neglected [21].
To ensure the renormalizability of the action, we have to retain only the relevant terms,  i.e., the ones with a positive momentum dimension.
The relevant part of action eq:action.from.coherent, where the fields have been already shifted, reads:

\begin{eqnarray}
\fl  \mathcal{S}= \int \rmd t \rmd^{d}\bi{r} \tilde{A}(\bi{r},t)\big[\partial_{t}-D_{A}\nabla^{2}\big]A(\bi{r},t) + \tilde{B}(\bi{r},t)\big[\partial_{t}-D_{B}\nabla^{2} \big]B(\bi{r},t) \nonumber  \\
+ \tilde{C}(\bi{r},t)[ \partial_{t}-D_{C}\nabla^{2} ]C(\bi{r},t) -\zeta B(\bi{r},t)\Big[ \tilde{C}(\bi{r},t)-\tilde{B}(\bi{r},t)\Big] \nonumber  \\ -
\sigma A(\bi{r},t)B(\bi{r},t)  \tilde{B}(\bi{r},t)-\sigma A_0 B(\bi{r},t) \Big[  \tilde{B}^2(\bi{r},t)+\tilde{B}(\bi{r},t)-\tilde{A}(\bi{r},t) \Big] \nonumber  \\ +   \gamma \, C(\bi{r},t)B(\bi{r},t) \tilde{B}(\bi{r},t) \label{eq:ren.action}
\end{eqnarray}

§.§ Integrating out the fields $A$ and $C$.

Action eq:ren.action can be simplified by integrating out the fields $A$ and $C$, as their corresponding response fields appear only linearly in the interaction terms since the non-linear terms, raising from stochastic noise, are irrelevant.
In particular, we show that the integral over the response fields results in a $\delta$-function constraining the evolution of $A(\bi{r},t)$ and $C(\bi{r},t)$ to obey the linearized part of the deterministic equations eq:deterministicequations1–eq:deterministicequations3. Note that this simplification is not necessary since once could equally renormalize the above action and would obtain the same flow functions eq:flow.function.CGEP1–eq:flow.function.CGEP3.
We proceed by collecting the terms proportional to $\tilde{A}$ and $\tilde{C}$

\begin{eqnarray}
\fl \int \mathcal{D}[B,\tilde{B},A ,C] \exp \left(- \mathcal{S}_{h} \right) \times \nonumber\\
\fl \int \mathcal{D}[ \tilde{A}] \exp  \left(- \int \rmd t \rmd^{d} \bi{r} \tilde{A}(\bi{r},t)\Big[ \partial_{t}A(\bi{r},t)-D_{A}\nabla^{2} A(\bi{r},t) +\sigma A_0 B(\bi{r},t) \Big] \right) \times \nonumber \\
\fl  \int \mathcal{D}[\tilde{C}]\exp   \left(- \int \rmd t \rmd^{d}\bi{r} \tilde{C}(\bi{r},t)\Big[ \partial_{t} C(\bi{r},t)-D_{C}\nabla^{2}C(\bi{r},t) -\zeta B(\bi{r},t)\Big] \right)= \nonumber \\
\fl  \int \mathcal{D}[B,\tilde{B},A,C] \exp  \left(- \mathcal{S}_{h} \right) \times \nonumber
\\
\delta_{\bi{r},t}\left(  \partial_{t}A(\bi{r},t)-D_{A}\nabla^{2} A(\bi{r},t) +\sigma A_0 B(\bi{r},t) \right) \times \nonumber\\
\delta_{\bi{r},t}\left(  \partial_{t} C(\bi{r},t)-D_{C}\nabla^{2}C(\bi{r},t) -\zeta B(\bi{r},t) \right), \nonumber
\end{eqnarray}

where $\delta_{\bi{r},t}$ denotes a delta-function imposing its argument to be zero at each point in space and time. To simplify notation, we introduced the action $\mathcal{S}_{h}$

\begin{eqnarray}
\fl \mathcal{S}_{h}= \int \rmd t \rmd^{d}\bi{r} \Big[ \tilde{B}(\bi{r},t)\left( \partial_{t}-D_{B}\nabla^{2}+ \zeta-A_0 \sigma \right) B(\bi{r},t)
-\sigma A_0 B(\bi{r},t)  \tilde{B}^2(\bi{r},t) \nonumber \\ - \sigma A(\bi{r},t)B(\bi{r},t)  \tilde{B}(\bi{r},t)
+\gamma \,  C(\bi{r},t)B(\bi{r},t) \tilde{B}(\bi{r},t)\Big].
\end{eqnarray}

As anticipated, integrating out the fields $A$ and $C$ results in a delta-function requiring that the two fields follow the mean-field equations eq:deterministicequations1–eq:deterministicequations3 linearized around the homogeneous solution $( A_0, 0, 0 )$. The non-linear part of the mean-field equations do not enter the action eq:ren.action as they are RG irrelevant at ${d_c=6}$.
The solution of the linearized mean-field equations can be easily obtained in Fourier space:

\begin{eqnarray}
A(\omega,\bi{k})= -\sigma A_0 \frac{ B(\omega,\bi{k})}{\rmi \omega+D_A k^2 } \qquad
C(\omega,\bi{k})=\zeta  \frac{B(\omega,\bi{k})}{\rmi \omega+D_C k^2 } \nonumber
\end{eqnarray}

By inserting the values of $A$ and $C$ in the Fourier transformed action $\mathcal{S}_{h}$, we get the final action

\begin{eqnarray}
\fl \mathcal{S}=\int_{\bi{q}}  \tilde{B}(-\bi{q}) \left(\rmi \omega +D_{B}k^{2}+\zeta-A_0 \sigma  \right)B(\bi{q}) -\int_{\bi{q}_1,\bi{q}_2 }\sigma A_0 B(-\bi{q}_1-\bi{q}_2)  \tilde{B}(\bi{q}_1)\tilde{B}(\bi{q}_2)   \nonumber \\
+ \int_{\bi{q}_1,\bi{q}_2} \Bigg[  \frac{\sigma^2 A_0 B(\bi{q}_1)}{\rmi \omega_1+D_A k_1^2 } +  \frac{\gamma \zeta B(\bi{q}_1)}{\rmi \omega_1+D_C k_1^2 } \Bigg] B(\bi{q}_2)  \tilde{B}(-\bi{q}_1-\bi{q}_2). \label{eq:final.action}
\end{eqnarray}

where we denote frequency $\omega$ and momentum $\bi{k}$ (${|\bi{k}|=k}$) by a 4-vector ${\bi{q}=(\omega, \bi{k})}$ and write the integral as ${\int_{\bi{q}}:=\int_{\mathbb{R}^{d+1}}   \rmd \omega \, \rmd^{d}\bi{k}/(2\pi)^{d+1}}$.
Here, ${B(\bi{q})= B(\omega, \bi{k})}$ denotes the density of the infected individuals in Fourier space, while ${\tilde{B} (\bi{q}) = \tilde{B} (\omega, \bi{k})}$ is its dual response field; for an interpretation on its physical meaning see [93, 21, 60].
To recover the right action for the three models discussed in the main text (GEP, DGEP, and CGEP) we set ${\tau= \zeta-A_0 \sigma}$ in the first term.
The couplings for the non-linear terms and the value of $D$ differ for each model.
Only the first non-linearity is the same in all cases ${U=2 \sigma A_0}$.
In GEP, the caregiving reaction is absent ($\gamma=0$) and $A$ immobile ($D_A=0$). In the general action eq:action this corresponds to ${V=0, W=\sigma^2 A_0}$.
Analogously, in DGEP, the caregiving reaction is not present (${\gamma=0}$), but ${D_A \neq 0}$.
In the general action eq:action this corresponds to ${V=\sigma^2 A_0, W=0, D=D_A}$.
In CGEP, all reactions are present, but the susceptible species $A$ is immobile ${D_A=0}$.
The action of CGEP can be obtained from eq:action setting ${V=\zeta\gamma, W=\sigma^2 A_0, D=D_C}$.

§ RENORMALIZATION

In this appendix, we present how to calculate the RG flow functions eq:flow.function.CGEP1–eq:flow.function.CGEP3 for the general epidemic model presented in the main text of this paper.
The flow functions are calculated using the $Z$-factor RG approach and by performing an $\epsilon$-expansion around the upper critical dimension ${d_c=6}$, below which the critical exponents acquire non-trivial values [94].
The idea of the $Z$-factor approach is to redefine couplings and fields such that the observables—which seemingly acquire infinite values for ${\epsilon \rightarrow 0}$—are well-defined. Those redefined parameters are termed renormalized, whereas the original parameters are referred to as bare. Throughout this section, we represent the bare parameters with an overscore ring (e.g. $\mathring{\tau}$ and $\mathring{D_B}$), while the renormalized parameters are denoted in a standard manner (e.g. $\tau$ and $D_B$).

Before renormalizing the action eq:action, it is convenient to rescale it in terms of effective couplings, allowing us to express the coefficients $(U,V,W,D)$ in terms of  $(g_2,g_1,g_D)$, otherwise, the resulting equations for the $Z$-factors are underdetermined. In the field-theoretical approach, the fields themselves are “integration parameters" integrated over each path (see e.g. eq:marginal.probability); thus, they can be rescaled by appropriately changing the integration measure.
To express the action in terms of effective couplings, we rescale the bare fields in the following way:

\begin{eqnarray}
\mathring{B}\rightarrow \left( \mathring{U} \mathring{D}_B \mathring{V}^{-1} \right)^{\frac{1}{2}} \mathring{B} \qquad \mathring{\tilde{B}} \rightarrow \left( \mathring{U} \mathring{D}_B \mathring{V}^{-1} \right)^{-\frac{1}{2}} \mathring{\tilde{B}}
\end{eqnarray}

Since the field and response field are inversely rescaled, the determinant of the Jacobian for this change of variable is unitarian, and the general expression for the probability of a density configuration is unchanged eq:marginal.probability.
The rescaled action reads:

\begin{eqnarray}
\fl \mathcal{S}= \int_\bi{q} \mathring{\tilde{B}} (-\bi{q})  \left( \rmi\omega+\mathring{D}_B k^{2}+\mathring{\tau} \right)
\mathring{B}(\bi{q})-\int_{\bi{q}_1,\bi{q}_2} \frac{(\mathring{g}_2)^{\frac{1}{2}} \mathring{D}_B}{2} \mathring{\tilde{B}}(-\bi{q}_1)\mathring{\tilde{B}}(-\bi{q}_2) \mathring{B}(\bi{q}_2+\bi{q}_1) \nonumber\\
\fl +\int_{\bi{q}_1,\bi{q}_2} \Bigg[ \frac{ \mathring{g}_1  \mathring{D}^2_B }{ (\mathring{g}_2)^{\frac{1}{2}}} \frac{\mathring{B}(\bi{q}_1)}{\rmi\omega_1}
+ (\mathring{g}_2)^{\frac{1}{2}} \mathring{D}^2_B\frac{\mathring{B}(\bi{q}_1)}{\rmi\omega_1+\mathring{g}_D \mathring{D}^2_B k_1^{2}}  \Bigg] \mathring{B}(\bi{q}_2)  \mathring{\tilde{B}} (\bi{-q}_1-\bi{q}_2) \, , \label{eq:action.reparam}
\end{eqnarray}

where we inserted the value of the bare effective couplings:

\begin{eqnarray}
\mathring{g}_{1}= \frac{\mathring{U} \mathring{W}}{(4 \mathring{D}_{B} \pi)^{3}},  \qquad
\mathring{g}_{D}:=\frac{\mathring{D}}{\mathring{D}_{B}},\qquad {\rm{and}} \nonumber \qquad
\mathring{g}_{2}:= \frac{\mathring{U} \mathring{V}}{(4 \mathring{D}_{B} \pi)^{3}}. \label{eq:effectiv.coupling}
\end{eqnarray}

In eq:action.reparam, we wrote the integral over a single momentum as ${\int_{\bi{q}}:=\int_{\mathbb{R}^{d+1}}   \rmd \omega \, \rmd^{d}\bi{k}} /(2 \pi)^{d+1}$, but we redefine the integral over the two momenta such that the factor $(4 \pi)^{3}$ appearing from the effective couplings is reabsorbed into the integral measure to lighter the notation: ${\int_{\bi{q}_1, \bi{q}_2 }:=~(4 \pi)^{3/2}\int_{\mathbb{R}^{2d+2}}(\rmd \omega_1 \rmd \omega_2 \rmd^{d}\bi{k}_1 \rmd^{d}\bi{k}_2)/(2\pi)^{2d+2}}$

The dimension of the effective couplings can be found imposing that the action is dimensionless. With respect to an arbitrary momentum scale $\mu$, this gives:

\begin{eqnarray}
[\mathring{g}_1]=[\mathring{g}_2]=\mu^{6-d} \qquad {\rm{and}} \qquad [\mathring{g}_D]=\mu^{0}
\end{eqnarray}

The effective couplings are irrelevant for dimensions greater than $d_c=6$ and therefore, flow to 0 under the RG flow. For dimensions smaller than $d_c=6$, instead,  they become relevant and get renormalized. Therefore, the upper critical dimension of action eq:action.reparam is ${d_c=6}$.

§.§ Reparametrization of the action and minimal subtraction scheme

We proceed to rewrite the action in terms of renormalized quantities. Those are defined using the $Z$-factors:

\begin{eqnarray}
\fl \mu^{\epsilon} Z_{g_1} g_1 := \mathring{ g_1}, \quad  Z^{\frac{1}{2}}_{B} B := \mathring{B}, \quad
\frac{Z_{\tau}}{\left( Z_{B}Z_{\tilde{B}}\right)^{\frac{1}{2}}}\tau  :=\mathring{\tau}, \quad \frac{\left( Z_{B}Z_{\tilde{B}}\right)^{\frac{1}{2}} }{Z_{D}} g_D := \mathring{ g}_D,   \nonumber \\
\mu^{\epsilon} Z_{g_2} g_2 := \mathring{ g_2},
\quad
Z^{{\frac{1}{2}}}_{\tilde{B}} \tilde{B} :=\mathring{\tilde{B}},    \qquad
\frac{Z_{D}}{\left( Z_{B}Z_{\tilde{B}}\right)^{\frac{1}{2}} } D_{B}:= \mathring{D}_{B}.
\label{eq:Z.factor.definition}
\end{eqnarray}

Note that the renormalized couplings $g_D$ and $D_B$ are defined through the same $Z$-factor, as the diffusion constant $D$ does not renormalize, i.e., ${\mathring{D}=D}$.
Inserting these values in the action eq:action.reparam, we can split it in two parts ${\mathcal{S}=\mathcal{S}_{R}+\mathcal{S}_{\rm{c}}}$.
The renormalized action $S_{R}$ has the same functional form as the bare one eq:action.reparam, containing renormalized couplings and fields. It reads:

\begin{eqnarray}
\fl \mathcal{S}_{R}=\int_\bi{q} \tilde{B} (-\bi{q})  \left( \rmi\omega+{D}_B k^{2}+\tau \right)
B(\bi{q})-\int_{\bi{q}_1,\bi{q}_2}  \frac{(g_2)^{\frac{1}{2}} D_B}{2} \tilde{B}(-\bi{q}_1)\tilde{B}(-\bi{q}_2) B(\bi{q}_2+\bi{q}_1) \nonumber\\
\fl +\int_{\bi{q}_1,\bi{q}_2} \Bigg[  \frac{g_1 D^2_B}{(g_2)^{\frac{1}{2}}} \frac{B(\bi{q}_1)}{\rmi\omega_1}
+ (g_2)^{\frac{1}{2}} D^2_B \frac{B(\bi{q}_1)}{\rmi\omega_1+g_D D_B k_1^{2}} \Bigg] B(\bi{q}_2)  \tilde{B} (\bi{-q}_1-\bi{q}_2). \label{eq:action.reparam.R}
\end{eqnarray}

The extra dependence on the momentum scale $\mu^{\epsilon}$ appearing in the effective couplings $g_1,g_2$ is reabsorbed in the integral measure of the interacting vertices to lighten the notation, i.e., $ {\int_{\bi{q}_1, \bi{q}_1 }:=~(4 \pi)^{3/2}\mu^{\epsilon/2}\int_{\mathbb{R}^{2d+2}}(\rmd \omega_1 \rmd\omega_2 \rmd^{d}\bi{k}_1 \rmd^{d}\bi{k}_2)/(2\pi)^{2d+2}}$.
This notation will be consistently assumed for integrals over two momenta in what follows.
The action $\mathcal{S}_{\rm{c}}$ is made of “counter-terms" containing the $Z$-factors:

\begin{eqnarray}
\fl \mathcal{S}_{\rm{c}}= \int_\bi{q} \tilde{B} (-\bi{q})\left( \big[(Z_{\tilde{B}}Z_{B})^{\frac{1}{2}}-1\big]\rmi\omega
+\big[Z_{D}-1\big]D_B k^{2} +\big[Z_{\tau}-1\big] \tau_{R} \right) B_{R}(\bi{q})
\nonumber\\  -\int_{\bi{q}_1,\bi{q}_2} \Big[Z^{\frac{1}{2}}_{g_2} Z_D Z^{\frac{1}{2}}_{\tilde{B}}-1\Big] \frac{(g_2)^{\frac{1}{2}} D_B}{2} \tilde{B}(-\bi{q}_1)\tilde{B}(-\bi{q}_2) B(\bi{q}_2+\bi{q}_1) \nonumber \\
+\int_{\bi{q}_1,\bi{q}_2} \Big[Z^{-\frac{1}{2}}_{g_2}  Z_{g_1} Z^2_D Z^{-\frac{1}{2}}_{\tilde{B}}-1 \Big] \frac{g_1 D^2_B}{(g_2)^{\frac{1}{2}}\rmi\omega_1} B(\bi{q}_1)  B(\bi{q}_2)  \tilde{B} (\bi{-q}_1-\bi{q}_2) \nonumber \\
+ \int_{\bi{q}_1,\bi{q}_2}  \Big[Z^{\frac{1}{2}}_{g_2} Z^2_D Z^{-\frac{1}{2}}_{\tilde{B}}-1 \Big]\frac{(g_2)^{\frac{1}{2}} D^2_B }{\rmi\omega_1+g_D D_B k_1^{2}} B(\bi{q}_1)
B(\bi{q}_2)  \tilde{B} (\bi{-q}_1-\bi{q}_2) \, . \nonumber
\end{eqnarray}

In the field-theoretical formalism, all physical observables can be expressed as expectation values of a product of fields [94].
These expectation values correspond to $n$-point correlation functions, or equivalently, vertex functions $\Gamma^{(m,n)}$ [21].
The $Z$-factors are defined such that the counter-terms of action $\mathcal{S}_c$ cancel the divergent contributions of the vertex functions $\Gamma_R$ calculated from the action $\mathcal{S}_R$.
This procedure ensures that observables are divergence-free.

Vertex functions are calculated through a diagrammatic expansion over one-particle irreducible (1PI) diagrams [65], organized by the number $n$ of internal loops of a diagram (or, equivalently, in terms of the deviation from the upper critical dimension $\epsilon^n$).
The $Z$-factors are perturbatively defined such that the (${n-1}$)-loop order expansion in $\mathcal{S}_{\rm{c}}$ cancels the $n$-loop divergences of the renormalized vertex functions $\Gamma_R$.
In our case, only $\Gamma_{R}^{(1,1)}$, $\Gamma_{R}^{(2,1)}$ and $\Gamma_{R}^{(1,2)}$ contain divergences. So, the $Z$-factors are determined through the following conditions:

\begin{eqnarray}
\label{eq:Z.factors.condition}
\fl \mathcal{O}(\epsilon^0) \overset{!}{=} \frac{(g_2)^{\frac{1}{2}} D_B( 1-(Z_{g_2}Z^{-1}_{\tilde{B}})^{\frac{1}{2}} Z^2_D  )}{\rmi\omega_1+g_D D_B k_1^2} +
\frac{(1- Z^2_D Z_{g_1}(Z_{\tilde{B}}Z_{g_2})^{-\frac{1}{2}} )}{\rmi\omega_1 } \frac{g_1 D_B}{(g_2)^{\frac{1}{2}}} -\Gamma_R^{(2,1)},\label{eq:Zfactorscondition1} \\
\fl \mathcal{O}(\epsilon^0) \overset{!}{=} \left( 1-Z_{\tau} \right) \tau_{R}-(Z_{D}-1)D_{B} k^{2}- \left( (Z_{\tilde{B}}Z_{B})^{{\frac{1}{2}}}-1\right)\rmi\omega-\Gamma_R^{(1,1)},  \label{eq:Zfactorscondition2}  \\
\fl \mathcal{O}(\epsilon^0) \overset{!}{=} \left( Z^{\frac{1}{2}}_{g_2} Z_D Z^{\frac{1}{2}}_{\tilde{B}}-1 \right) (g_2)^{\frac{1}{2}} D^2_B - \Gamma_R^{(1,2)} .\label{eq:Zfactorscondition3}
\end{eqnarray}

To lighten notation, we will consistently use vertex functions rescaled by a factor $\mu^{ \epsilon/2}(4 \pi)^{3/2}$ to compensate the same factor previously reabsorbed in the integral measure of $\int_{\bi{q}_1,\bi{q}_2}$.
Note that the factor $2$ in front of the vertex $ ((g_2)^{1/2} D_B)/2$ arises from combinatorial factors.

§.§ Calculation of the integral of diagram (A)

Next, we give an example of how to evaluate the divergent part of a Feynman diagram needed to obtain the $Z$-factors. We choose to evaluate diagram (A) since it is responsible for creating an asymmetry in the RG flow functions of $\beta_{g_1}$, and $\beta_{g_2}$ eq:flow.function.CGEP1–eq:flow.function.CGEP3, making the RG flow approach the new DGEP fixed point.

fig:figureB1 reports the Feynman rules for the renormalized action eq:action.reparam.R together with the diagram (A), which contributes to the vertex function $\Gamma_R^{(2,1)}$. In the diagram, the external frequencies are denoted with $(\omega,\bi{k})$ and the internal momentum over which to integrate by $(\omega_{p},\bi{p})$. The integral corresponding to diagram (A) reads:

\begin{eqnarray}
\fl (A)= (D_B)^4 \mu^{\frac{3 \epsilon}{2}} (4 \pi)^{\frac{9}{2}} \int \frac{\rmd^{d} \bi{p}}{(2 \pi)^{d}} \frac{\rmd w_{p}}{(2 \pi)}
\frac{1}{\rmi(w_{2}+w_{p})+D_{B}| \bi{k}_{1}+\bi{p}|^{2}+\tau} \frac{1}{\rmi w_{p}+D_B k^2+\tau} \times           \nonumber\\
\fl \frac{1}{\rmi(w_{1}-w_{p})+D_{B}|\bi{k}_{2}-\bi{p}|^{2}+\tau}
\left( \frac{g_1 }{\rmi(g_2)^{\frac{1}{2}}(w_{1}-w_{p})}+ \frac{(g_2)^{\frac{1}{2}}}{\rmi(w_{1}-w_{p})+D|\bi{k}_{1}-\bi{p}|^{2}} \right) \times \nonumber \\
\left(\frac{g_1 }{\rmi w_{p}}+ \frac{g_2 }{\rmi w_{p}+D p^2} \right).\label{eq:diagram.A}
\end{eqnarray}

We inserted ${D=g_D D_B}$ to lighten the notation.
Note that the factor $\mu^{3 \epsilon/2} (4 \pi)^{9/2}$ was previously reabsorbed in the measure for integrals over two momenta.
We will show that the only divergence arising from this integral is renormalizing the vertex proportional to $g_1D^2_B\left((g_2)^{1/2}\rmi\omega_1\right)^{-1}$ in eq:action.reparam.
Therefore, in the following calculation, $\omega_1$ is kept non-zero as it explicitly enters the vertex structure, which we aim to renormalize. At the same time, we can set the other external momenta to zero—i.e., ${\omega_2=\bi{k}_1=\bi{k}_2=0}$—since the divergences cannot depend on their values.

We first evaluate the integral over $\omega_p$ by using the residue theorem in the complex plane.
We close the contour in the lower half of the complex plane, where there are three poles $z_i$:
${z_{1} = \omega_1}$, ${z_{2}= \omega_1-\rmi D p^2}$ and ${z_{3}= \omega_1- \rmi (D_B p^2+\tau)}$. We use a forward (Itô) prescription where the pole $z_1$ is included in the lower half plane through the shift ${z_1= \lim_{\delta \rightarrow 0}\omega_1-\rmi \delta}$ [21].

The integral over momentum $ \bi{p}$ is evaluated just below the upper critical dimension in $d=6-\epsilon$ dimensions with a finite diffusion constant $D\neq 0$. The only divergent contribution of eq:diagram.A in the limit $\epsilon \rightarrow 0$ is the  residue of $z_1$ inserted in the term proportional to $g_1(g_2)^{1/2}$ of eq:diagram.A, which gives rise to the integral $I(D=0)$ with $I(D)$ defined as:

\begin{eqnarray}
\fl  I(D)=  \int \frac{\rmd^{d} \bi{p}}{(2\pi)^d}
\frac{1}{(D_B-D) \bi{p}^2 +\tau} \, \frac{1}{\big[ \rmi\omega_1 +(D_B+D)\bi{p}^2 +\tau \big]^2} \, \frac{1}{\rmi\omega_1 +2 D \bi{p}^2} \label{eq:niceint}
\end{eqnarray}

To better understand the limit $D \rightarrow 0$, we also consider the contribution of the residue of $z_2$ inserted in the term proportional to $g_2^{3/2} $, giving rise to the integral $I(D)$.
Summarizing, diagram (A) gives, among others, the two contributions:

\begin{eqnarray}
\fl  \mu^{\frac{3 \epsilon}{2}} (4 \pi)^{\frac{9}{2}}g_2^{\frac{3}{2}} D^4_{B}I(D) \quad {\rm{and}} \qquad \mu^{\frac{3 \epsilon}{2}} (4 \pi)^{\frac{9}{2}}g_1^2g_2^{-\frac{1}{2}} D^4_B I(D=0).
\end{eqnarray}

In the main text, we referred to these contributions as $V D^3_B g_2 I(D)$ and $W D^3_B g_1 I(0)$, as these are the value they would assume in the old parametrization of the action eq:action.
Note that the first term, if divergent, could renormalize the vertex proportional to $(g_2)^{1/2}D_B$ in the second line of eq:action.reparam.R (in which case setting $k_1=0$ would not be justified). In contrast, the second renormalizes the vertex proportional to  $(g_1)^2/(g_2)^{1/2}$. By simple power counting, we expect only $I(0)$ to be divergent so that the two vertices in the second line of eq:action.reparam.R renormalize differently, leading to asymmetric flow functions for $g_1$ and $g_2$. In the following, we evaluate $I(D)$ and $I(0)$ to prove that only the latter is relevant to the renormalization and to illustrate how such integrals are performed.
We achieve this by first collecting the diffusion constant appearing in the propagators and then using identity eq:identiti2, known as the Feynman trick, to bring the integrals in the following form:

\begin{eqnarray}
\fl   I(0)=\frac{\Gamma(3)}{\rmi\omega_1 D^3_B } \int_{0}^{1} \rmd x \, x \int \frac{\rmd^{6-\epsilon} \bi{p}}{(2 \pi)^{d}}  \left( \bi{p}^2+(1-x)\frac{\tau}{D_B}+x\frac{\rmi\omega_1+\tau}{D_B}\right)^{-3} \\
\fl I(D)=\frac{\Gamma(4)}{D(D_B-D)(D_B+D)^2} \times \nonumber\\ \int_0^1 \rmd x  \rmd y \,x \int\frac{\rmd^{6-\epsilon} \bi{p}}{(2 \pi)^{d}} \left( p^2+ \frac{(1-x-y)\tau}{(D_B-D)} + x\frac{\rmi\omega_1+\tau}{D_B+D}+y\frac{\rmi\omega_1}{D}\right)^{-4}
\end{eqnarray}

The integral over momentum $\bi{p}$ can now be expressed in function of the Euler-Gamma function $\Gamma$—the analytical continuation of the factorial function—using identity eq:identiti1:

\begin{eqnarray}
I(0)=&\frac{1}{\rmi\omega_1 D^3_B } \int_{0}^{1} \rmd x \, x   \frac{\Gamma(\frac{\epsilon}{2}) }{(4 \pi )^{3}} \left((1-x)\frac{\tau}{D_B}+x\frac{\rmi\omega_1+\tau}{D_B}\right)^{-\frac{\epsilon}{2}}
\label{eq:value.of.I1}\\
I(D)=&\frac{1}{D(D_B-D)(D_B+D)^2}\int_{0}^{1} \rmd x \rmd y \, x \frac{\Gamma(1+\frac{\epsilon}{2})}{(4\pi )^3} \times \nonumber \\ &\left(  (1-x-y)\frac{\tau}{(D_B-D)} +x \frac{\rmi\omega_1+\tau}{D_B+D}+y\frac{\rmi\omega_1}{D}\right)^{-1-\frac{\epsilon}{2}}
\end{eqnarray}

To obtain the leading order of the integrals $I(0)$ and $I(D)$, we expand in $\epsilon$. In particular, we use

\begin{eqnarray}
\Gamma(\epsilon)=\frac{1}{\epsilon}+\mathcal{O}(\epsilon^0)\qquad
{\rm{and}} \qquad \Gamma(1+\epsilon)=\mathcal{O}(\epsilon^0) \nonumber
\end{eqnarray}

such that the divergent part of $I(0)$ in the limit $\epsilon \rightarrow 0$ can be expressed as a $1/ \epsilon$ pole.

\begin{eqnarray}
I(0)=&\frac{1}{\rmi\omega_1 D^3_B } \int_{0}^{1} \rmd x \, x   \frac{2 }{\epsilon (4 \pi )^{3}}= \frac{1}{\rmi\omega_1 D^3_B  (4 \pi )^{3}} \frac{1}{\epsilon}
\label{eq:value.of.I}\\
I(D)=&\frac{1}{D(D_B-D)(D_B+D)^2} \int_{0}^{1} \rmd x \rmd y \frac{x}{(4\pi )^3} \times \nonumber \\
&\left(  (1-x-y)\frac{\tau}{(D_B-D)} +x \frac{\rmi\omega_1+\tau}{D_B+D}+y\frac{\rmi\omega_1}{D}\right)^{-1}
\end{eqnarray}

The integral over $x,y$ in $I(D)$ is finite as long as $D\neq0$, as was claimed initially. In the limit $D \rightarrow 0$, the integral over $x,y$ scales as $D \ln(D)$, implying that the integral $I(D \rightarrow 0) \propto \lim_{D\rightarrow 0}\ln(D)$ becomes divergent. On the other hand, the divergent part of $I(0)$ is parameterized in terms of a $1/\epsilon$ pole. This way, we showed that the $1/\epsilon $ pole corresponds to a logarithmic divergence $\lim_{D\rightarrow 0}\ln(D)$ if the limit $D\rightarrow 0$ is taken after the integral is performed. Overall, the only divergent contribution of (A) to the vertex functions in the limit $\epsilon \rightarrow0$ is reported below. To account for the right dimensionality, we multiply by the factor $\mu^{d-d_c}$ to fix the dimensionality of the integral $I(0)$ in eq:value.of.I.

\begin{eqnarray}
\mu^{3 \epsilon/2} (4 \pi)^{\frac{9}{2}} \frac{(g_1)^2}{(g_2)^{\frac{1}{2}}} D^4_B \mu^{-\epsilon}I(0)= \mu^{ \epsilon/2}(4 \pi)^{\frac{3}{2}} \frac{(g_1)^2 D_B }{(g_2)^{\frac{1}{2}}}\frac{1}{\rmi\omega_1 \epsilon}
\end{eqnarray}

§.§ $Z$-factors and flow functions

Feynman rules and diagrams. The Feynman rules and one-loop diagrams contributing to $\Gamma_R^{(2,1)}$ (diagrams (A),(B),(C) and (D)), $\Gamma_R^{(1,2)}$ (diagrams (E) and (F)) and $\Gamma_R^{(1,1)}$ (diagram (G)).

fig:figureB1 shows all Feynman diagrams contributing to the vertex functions  $\Gamma_{R}^{(1,1)}$, $\Gamma_{R}^{(2,1)}$ and $\Gamma_{R}^{(1,2)}$. We evaluate the divergent part of those following the same procedure as was done for diagram (A), i.e., the integrals are evaluated at dimension ${d=d_c-\epsilon}$ and solved with the help of the two identities eq:identiti1 and eq:identiti2.
In the limit ${\epsilon \rightarrow 0}$, the divergent contributions appear as poles proportional to $1/\epsilon$.
The divergent contribution of the Feynman diagrams shown in figure <ref> are reported below:

\begin{eqnarray}
(A)=& \mu^{ \frac{\epsilon}{2}}(4 \pi)^{\frac{3}{2}}\frac{g_1 D^2_B}{(g_2)^{\frac{1}{2}} \rmi\omega_1} \frac{g_1}{ \epsilon}\nonumber \\
(B)=0  \nonumber \\
(C)=&\mu^{ \frac{\epsilon}{2}}(4 \pi)^{\frac{3}{2}} \left( \frac{(g_2)^{\frac{1}{2}} D^2_B} {\rmi\omega_1+ D k_1^2} + \frac{g_1 D^2_B}{(g_2)^{\frac{1}{2}} \rmi\omega_1} \right) \left( \frac{g_1}{\epsilon} \frac{3}{4}+ \frac{g_2}{\epsilon} \frac{3+g_{D}}{4(1+g_{D})^2} \right)\nonumber
\\
(D)=&\mu^{ \frac{\epsilon}{2}}(4 \pi)^{\frac{3}{2}}\left( \frac{(g_2)^{\frac{1}{2}} D^2_B} {\rmi\omega_1+ D k_1^2} + \frac{g_1 D^2_B}{(g_2)^{\frac{1}{2}} \rmi\omega_1} \right)  \left( \frac{g_1}{\epsilon}\frac{1}{4}+\frac{g_2}{\epsilon}\frac{1}{4(1+g_{D})} \right) \nonumber
\\
(E)=&-\mu^{ \frac{\epsilon}{2}}(4 \pi)^{\frac{3}{2}}(g_2)^{\frac{1}{2}} D_B \left( \frac{g_1}{\epsilon}\frac{1}{4}+\frac{g_2}{\epsilon}\frac{1}{4(1+g_{D})}  \right) \nonumber
\\
(F)=&- \mu^{ \frac{\epsilon}{2}}(4 \pi)^{\frac{3}{2}}(g_2)^{\frac{1}{2}} D_B\left( \frac{g_1}{\epsilon}\frac{3}{4}+\frac{g_2}{\epsilon}\frac{3+g_{D}}{4(1+g_{D})^2}  \right) \nonumber \\
(G) =& \mu^{ \frac{\epsilon}{2}}(4 \pi)^{\frac{3}{2}} \Big[\frac{g_1}{\epsilon}\left(\frac{3}{4}\rmi\omega+\tau_R+D_B \frac{k^2}{6}\right)  + \nonumber \\ &\frac{g_2}{\epsilon}\left( \frac{3+g_{D}}{4(1+g_{D})^2}\rmi\omega+\frac{2+g_{D}}{2(1+g_{D})^2} \tau_R + \frac{1+4g_{D}+g_{D}^2}{6(1+g_D)^3}D_B k^2 \right) \Big]\nonumber
\end{eqnarray}

In the above equations, each divergent contribution, i.e., each $1/\epsilon$ pole, proportional to $g_2$ becomes equal to the one proportional to $g_1$ in the limit of ${g_D=0}$, except for the diagram ($A$), where only the term $I(0)$ proportional to $g_1$ is present.

To obtain the divergent part of the vertex functions, we sum up the corresponding contributions of the Feynman diagrams.
As we mentioned before, we rescale the vertex function by $\mu^{ \epsilon/2}(4 \pi)^{3/2}$ to be consistent with the choice of the integral measure $\int_{\bi{q}_1,\bi{q}_2}$.

\begin{eqnarray}
\fl \Gamma_{R}^{(2,1)} =&-(C)-(D)-(A)= -\frac{ g_1 D_B}{(g_2)^{\frac{1}{2}} \rmi\omega_1} \left(\frac{2 g_1}{\epsilon}+\frac{g_2}{\epsilon} \frac{g_{D}+2}{2(g_{D}+1)^2} \right) \nonumber \\
\fl  &-  \frac{ (g_2)^{\frac{1}{2}} D_B }{\rmi\omega_1+ g_D D_B k_1^2} \left( \frac{g_1} {\epsilon}+ \frac{g_2} {\epsilon} \frac{g_{D}+2}{2(g_{D}+1)^2} \right),
\\
\fl \Gamma_{R}^{(1,2)} =&-2(E)-2(F)=
(g_2)^{\frac{1}{2}} D^2_B \left( \frac{2g_1}{\epsilon} + \frac{g_2}{\epsilon}   \frac{g_{D}+2}{(1+g_{D})^2} \right),   \\
\fl \Gamma_{R}^{(1,1)} =&-(G)=
-\frac{g_1}{\epsilon}\left(\frac{3}{4}\rmi\omega+\tau_R+D_B \frac{k^2}{6}\right) \nonumber \\ \fl &+\frac{g_2}{\epsilon}\left( \frac{3+g_{D}}{4(1+g_{D})^2}\rmi\omega+\frac{2+g_{D}}{2(1+g_{D})^2} \tau_R + \frac{1+4g_{D}+g_{D}^2}{6(1+g_D)^3}D_B k^2) \right) \, ,
\end{eqnarray}

where the $2$ in front of (E) and (F) in $\Gamma_{R}^{(1,2)}$ is due to a combinatorial factor.
To find the $Z$-factors, we employ the conditions given by the minimal subtraction scheme eq:Zfactorscondition1–eq:Zfactorscondition3. This yields:

\begin{eqnarray}
\fl Z_{g_{1}}=1+\frac{1}{\epsilon}\left(\frac{7}{2}g_{1}+g_{2}\frac{5+5g_{D}+2g_{D}^2}{2(1+g_{D})^3}\right) \qquad
&Z_{D} =1+\frac{1}{\epsilon}\left(\frac{g_{1}}{6}+g_{2}\frac{1+4g_{D}+g_{D}^2}{6(1+g_{D})^3}\right)  \nonumber\\
\fl Z_{g_{2}}=1+\frac{1}{\epsilon}\left(\frac{5}{2}g_{1}+g_{2}\frac{5+5g_{D}+2g_{D}^2}{2(1+g_{D})^3}\right)
\qquad &Z_{\tau}=1+\frac{1}{\epsilon}\left( g_{1}+g_{2}\frac{2+g_{D}}{2(1+g_{D})^2} \right)
\nonumber\\
\fl Z_{\tilde{B}} =1+\frac{1}{\epsilon} \left( \frac{7}{6}g_1+ g_2{\epsilon}\frac{7+13g_{D}+4g_{D}^2}{6(1+g_{D})^3}\right)  \qquad
&Z_{B}=1+\frac{1}{\epsilon} \left(\frac{1}{3} g_1 +\frac{g_2}{\epsilon}\frac{2-g_{D}-g_{D}^2}{6(1+g_{D})^3}\right) \nonumber
\end{eqnarray}

Having calculated the value of the $Z$-factors, we are finally in the position to calculate the flow functions for the general action eq:action. To do so, we rewrite the effective couplings in terms of $Z$-factors and bare couplings and exploit the fact that the bare couplings are independent of the momentum scale $\mu$. Therefore, the derivatives of the bare parameters with respect to $\mu$ have to be zero, e.g.,  $\frac{d}{d \mu} \mathring{g}_2=0$.
The flow functions at one-loop order read:

\begin{eqnarray}
\fl \beta_{g_{1}}:=\mu \frac{d}{d\mu}g_{1}=\mu \frac{d}{d\mu} \left(  \mu^{-\epsilon} Z^{-1}_{g_1} \mathring{g}_{1}\right)=-\epsilon g_{1}+g_{1}\left(\frac{7}{2}g_{1}+g_{2}\frac{5+5g_{D}+2g_{D}^2}{2(1+g_{D})^3} \right)\label{eq:flowfunctions1C}\\
\fl \beta_{g_{2}}:=\mu \frac{d}{d\mu}g_{2}=\mu \frac{d}{d\mu} \left(  \mu^{-\epsilon}  Z^{-1}_{g_2} \mathring{g}_2\right)=-\epsilon g_{1}+g_{1}\left(\frac{5}{2}g_{1}+g_{2}\frac{5+5g_{D}+2g_{D}^2}{2(1+g_{D})^3}\right)\label{eq:flowfunctions2C}\\
\fl \beta_{D}:=\mu \frac{d}{d\mu}g_D=\mu \frac{d}{d\mu} \left(    \frac{\left( Z_{B}Z_{\tilde{B}}\right)^{\frac{1}{2}}}{Z_D} \frac{D}{ \mathring{D}_B } \right)=-g_{D}\left(\frac{7}{12}g_{1}+g_{2}\frac{7+4g_{D}+g_{D}^2}{12(1+g_{D})^3}\right) \label{eq:flowfunctions3C}
\end{eqnarray}

A fixed point in the flow function of the effective couplings $\beta_{g_1},\beta_{g_2},\beta_{g_D}$ represents a parameter regime where the system is scale-invariant, as expected from a second-order phase transition. The above flow function posses three different fixed points:  GEP fixed point (${g_1^*=2\epsilon/7,g^*_2=g^*_D=0}$), the  DGEP fixed point (${g_1^*=0,g^*_2=2 \epsilon /5,g^*_D=0}$) and the gaussian fixed line (${g_1^*=g^*_2=0,g^*_D}$) depicted in figure <ref> in the main text.

It appears unusual that, at the GEP fixed point, the vertex $g_1 / (g_2)^{1/2}$ diverges. This divergence is a consequence of the selected parametrization and does not relate to physical divergences. The problem can be overcome by changing the parametrization of the action so that no divergences arise at the GEP fixed point. One possibility is:

\begin{eqnarray}
\mathring{B}\rightarrow \left( \mathring{U} \mathring{D}_B \mathring{W}^{-1} \right)^{\frac{1}{2}} \mathring{B}, \quad \mathring{\tilde{B}} \rightarrow \left( \mathring{U} \mathring{D}_B \mathring{W}^{-1} \right)^{-\frac{1}{2}} \mathring{\tilde{B}} . \nonumber
\end{eqnarray}

§.§ Standard integrals

Lastly, we report two standard identities, which were used to evaluate the divergent part of integrals (A)–(G) through dimensional regularization

\begin{eqnarray}
\fl  \int_{-\infty}^{\infty} \frac{\rmd^{d} q}{(2 \pi)^{d}} \frac{q^{n}}{(q^2+\tau)^{a}}=\frac{1}{(4 \pi)^{\frac{d}{2}}\Gamma(\frac{d}{2})} \frac{\Gamma(\frac{d+n}{2})\Gamma(a-\frac{d+n}{2})}{\Gamma(a)} \tau^{\frac{d+n}{2}-a} \label{eq:identiti1}
\end{eqnarray}

\begin{eqnarray}
\fl  \frac{1}{P_1^{a_1}...P_n^{a_n}}=\frac{\Gamma(a_1+...+a_n)}{\Gamma(a_1)...\Gamma(a_n)} \int_{0}^{1} \rmd x_1...\rmd x_n \left( \frac{\delta(1-x_1...-x_n)x_1^{a_1-1}...x_n^{a_n-1}}{(x_1 P_1+...+x_n P_n)^{a_1+...+a_n}}\right) \label{eq:identiti2}
\end{eqnarray}

where $\Gamma$ stands for the Euler-Gamma function, an analytical continuation of the factorial function.

§ CRITICAL EXPONENTS

This section aims to derive the value of the critical exponents $z,\nu_{\parallel},\theta,\delta$ describing the scaling of the number of infected individuals, the survival probability, and the mean-squared radius:

\begin{eqnarray}
\langle N \rangle \left(t,\sigma \right)
&=& t^{\theta} \hat{N}\left( (\sigma-\sigma_c)^{\nu_{\parallel}} t \right),
\label{eq:observable1C}\\
P_{\rm s}(t,\sigma)
&=& t^{-\delta} \hat{P}_{\rm s}\left((\sigma-\sigma_c)^{\nu_{\parallel}} t\right) \, ,
\label{eq:observable2C}
\\
\langle R^2 \rangle (t,\sigma)
t^{{2}/{z}} \hat{R}^2\left((\sigma-\sigma_c)^{\nu_{\parallel}} t\right)
\, .
\label{eq:observable3C}
\end{eqnarray}

While the $Z$-factors of the effective couplings $(Z_{g_2}, Z_{g_1},Z_{g_D})$ determine the RG flow functions and fixed points, the other $Z$-factors $(Z_{B}, Z_{\tilde{B}}, Z_{\tau}, Z_{D})$ encode the scaling correction of the renormalized parameters as compared to the bare ones and determine the value of the critical exponents.

§.§ Anomalous dimension

By imposing that the action eq:action.reparam.R is dimensionless, the dimension of the bare parameters can be obtained:

\begin{eqnarray}
[\mathring{B}]=\mu^{(d+2)/2} \quad [\mathring{\tilde{B}}]=\mu^{(d-2)/2} \quad
[\mathring{\tau}]=\mu^{2} \quad [\mathring{D}_B]=\mu^{0}
\end{eqnarray}

where dimensionality is measured with respect to the momentum scale $\mu$.
Close to criticality, the scaling of physical parameters, i.e., the renormalized parameters, is not determined only by dimensional analysis; they acquire non-trivial correction termed anomalous dimension. To account for these corrections, we give a fictional dimensionality to the $Z$-factors—the anomalous dimensions $W_D, W_\tau, \sigma $ and $\tilde{\sigma}$—:

\begin{eqnarray}
\fl  Z_{D}:=\mu^{W_D}\qquad  Z_{\tau}:=\mu^{W_\tau} \qquad  Z_{B}:=\mu^{\sigma} \qquad  Z_{\tilde{B}}:=\mu^{\tilde{\sigma}} \label{eq:dimensionZ.factor}
\end{eqnarray}

This way, the scaling of any renormalized parameter is given in terms of dimensional analysis and the corresponding RG corrections, encoded in the scaling of the $Z$ factors eq:dimensionZ.factor. This is done using the definitions for the renormalized parameters eq:Z.factor.definition, such that the correct scaling for the physical parameters is recovered.
To obtain the anomalous dimensions $(W_D,W_{\tau},\sigma,\tilde{\sigma})$, we take the derivative of the $Z$-factors with respect to $\ln(\mu)$, ,i.e., $\mu \frac{ d }{d \mu} Z$. As the $Z$-factors depend only on effective couplings, the derivative $\mu \frac{d}{d \mu}$ can be evaluated using the flow functions eq:flowfunctions1C–eq:flowfunctions3C. At a fixed point, the effective couplings stop to flow, and the anomalous dimensions eq:dimensionZ.factor approach a finite constant value, encoding the RG correction to the parameters' scaling. Theories flowing to the same RG fixed point are characterized by parameters with the same scaling behavior, leading to the concept of universality class. Below are calculated the anomalous dimensions at the DGEP ($g^*_2=2\epsilon /5, g^*_1=g^*_D=0$) fixed point:

\begin{eqnarray}
W^*_D &= \mu \frac{d}{d \mu} \ln Z_D \bigg|_{FP}= g^*_{2}\frac{1+4g^*_{D}+g_{D}^{*2}}{6(1+g^*_{D})^3}=\frac{1}{15} \epsilon \nonumber\\
W^*_{\tau} &= \mu \frac{d}{d \mu}Z_{\tau} \bigg|_{FP} =-g^*_{R}\frac{2+D^*}{2(1+D^2)^2}=-\frac{5}{2} \epsilon \nonumber\\
\sigma^* &= \mu \frac{d}{d \mu } Z_B \bigg|_{FP} = -g^* \frac{2-D^*-(D^*)^2}{6(1+D^*)^3}=-\frac{2}{15} \epsilon \hspace{0,8cm} \nonumber \\
\tilde{\sigma}^* &= \mu \frac{d}{d \mu } Z_{\tilde{B}}\bigg|_{FP}=- g_2^* \frac{7+13g_{D}^*+4g_{D}^{*2}}{6(1+g_{D}^*)^3}= -\frac{7}{15}\epsilon \hspace{0,3cm} \nonumber
\end{eqnarray}

where the asterisk denotes that the anomalous dimensions
are evaluated at the DGEP fixed point.
From the knowledge of the scaling correction from the $Z$-factors, we calculate the critical scaling of the renormalized fields and parameters, inferred from eq:Z.factor.definition:

\begin{eqnarray}
\lbrack B\rbrack =[\mathring{B}]Z_{B}^{-\frac{1}{2}}=\mu^{\frac{d+2-\sigma^{*}}{2}} \label{appendixC:scaling1} \\
\lbrack \tilde{B} \rbrack=
=\mu^{\frac{d-2-\tilde{\sigma}^{*}}{2}} \label{appendixC:scaling2} \\
\lbrack\tau \rbrack=[\mathring{\tau}]
=\mu^{2-W^*_{\tau}+\frac{\sigma+\tilde{\sigma}}{2}} \label{appendixC:scaling3}\\
\lbrack D_B \rbrack=[\mathring{D}_B]Z_D^{-1}Z_B^{\frac{1}{2}}Z_{\tilde{B}}^{\frac{1}{2}}=\mu^{-W^{*}_D+\frac{\sigma+\tilde{\sigma}}{2}} \label{appendixC:scaling4}
\end{eqnarray}

§.§ Scaling of Green's functions

At criticality, the scaling of the renormalized parameters changes the scaling exponents of the observables eq:observable1C–eq:observable3C as compared with their scaling away from criticality, given by simple dimensional analysis. Using the field-theoretical formalism, any physical observable can be expressed as an expectation value of fields. Those expectation values are the connected Green functions $G^{(m,n)}_R$ of the theory:

\begin{eqnarray}
G^{(m,n)}_R(\tau,\{\bi{q}_i\},{g_i})=\Bigl \langle B(\bi{q}_1)...B(\bi{q}_m) \tilde{B}(\bi{q}_{m+1})...\tilde{B}(\bi{q}_{m+n}) \Bigr \rangle
\nonumber
\end{eqnarray}

where $\langle \cdot \rangle$ denotes a path integral expectation value weighted with the action $S_R$.
The scaling of the Green functions can be found using the method of characteristics and equations appendixC:scaling1–appendixC:scaling4
This procedure is called the Callan-Symanzik formalism [21] and allows the expression of the scaling of the renormalized Green functions at different momentum scales $\mu$.
Here we state the final expression for $\mu \rightarrow 0$, i.e., the IR critical regime:

\begin{eqnarray}
\ G^{(m,n)}_R(\tau,\bi{r},t,\{g_{i}\})= ( t)^{-n\frac{d-2+\tilde{\sigma}^*}{2z}-m\frac{d+2+\sigma}{2z}}\Phi^{(m,n)} \Bigl(\frac{\bi{x}}{\xi} , \frac{t}{\xi^{z}} \Bigr) \label{eq:Green.function.scaling}
\end{eqnarray}

where $\Phi^{(m,n)}$ is a generic non-linear scaling function just like $\hat{P_s}, \hat{N}, \langle \hat{R}^2 \rangle  $ in eq:observable1C–eq:observable3C. The correlation length $\xi$ and the critical exponent $z$ appear in its argument.
On the left of the expression, the power of time comes from the scaling of the renormalized fields $B_R$ and $\tilde{B}$, given by appendixC:scaling1 and appendixC:scaling2.
Here, we will present a heuristic argument of how the scaling of $\xi$ and the critical exponent $z$ can be found. For a formal derivation, please refer to [21].

Requiring that all terms in the propagator of the field $B$ have the same dimension, e.g. $[D_B k^2]=[\omega]$, and using appendixC:scaling1–appendixC:scaling3, one can infer

\begin{eqnarray}
[\omega]:= \mu^{z}\qquad \rm{with}\qquad  z=2-W^{*}_D+\frac{\sigma^*+\tilde{\sigma}^*}{2}.
\end{eqnarray}

A similar argument can be made for the scaling of the correlation length $\xi$. The correlation length of the mean-field system can be expressed in terms of bare parameter (by simple dimensional analysis): $\mathring{\xi}=(\mathring{\tau}/\mathring{D}_B)^{1/2}$.
Assuming that the scaling of $\xi$ is corrected by anomalous dimensions, while its functional form remains the same at criticality, we write $\xi=(\tau / D_B)^{\nu}$ where $\nu$ is the critical exponent. One can determine the value $\nu$ by requiring that the renormalized correlation length has the right dimension:

\begin{eqnarray}
[\tau/ D_B]^{\nu}=\mu^{-1} \rightarrow \nu=(2+W_D-W_{\tau})^{-1}.
\end{eqnarray}

§.§ Scaling of physical observables

To complete the derivation of the scaling form of the observables eq:observable1C–eq:observable3C, we have to express them in terms of $G^{(m,n)}_R$. Here, we only report the final results. For a derivation of the expression for the survival probability $P_s$ see [95], for the derivation of $\langle N \rangle$ and $\langle R^2 \rangle$ see [26]

\begin{eqnarray}
P_s (t,\tau) \propto \langle \tilde{B} (\bi{r}=0,-t) \rangle =G^{(0,1)}_R(\tau,\bi{r}=0,t) \label{appendixC:scaling.of.P} \\
\langle N  \rangle (t,\tau) = \int \rmd^d \bi{r} \langle \tilde{B}(0,0) B(\bi{r},t) \rangle =\int \rmd^d \bi{r} G^{(1,1)}_R(\tau,\bi{r},t) \label{appendixC:scaling.of.N} \\
\langle R^2 \rangle (t,\tau) = \frac{1}{\int \rmd^d \bi{r} \langle B(\bi{r},t) \rangle} \int \rmd^d \bi{r} \, r^2 \langle B(\bi{r},t) \rangle \nonumber \\
\textcolor{white}{\langle R^2 \rangle (t,\tau)} =\frac{1}{\int \rmd^d \bi{r} G^{(1,0)}_R( \tau,\bi{r},t) } \int \rmd^d \bi{r} \, r^2 G^{(1,0)}_R(\tau,\bi{r},t) \label{appendixC:scaling.of.R}
\end{eqnarray}