# Dynamics of soft interacting particles on a comb

Davide Venturelli, Pierre Illien, Aurélien Grabsch, Olivier Bénichou

► To cite this version:

Davide Venturelli, Pierre Illien, Aurélien Grabsch, Olivier Bénichou. Dynamics of soft interacting particles on a comb. Journal of Physics A: Mathematical and Theoretical, 2025, 58 (21), pp.215001. 10.1088/1751-8121/add590 . hal-05075884v2

HAL Id: hal-05075884

https://cnrs.hal.science/hal-05075884v2

Submitted on 25 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
**PAPER • OPEN ACCESS**

You may also like

# Dynamics of soft interacting particles on a comb

- Exponential speed up in Monte Carlo
sampling through radial updates
Johann Ostmeyer

To cite this article: Davide Venturelli et al 2025 J. Phys. A: Math. Theor. 58 215001

- On the diagonals of rational functions: the
minimal number of variables
S Hassani, J M Maillard and N Zenine

View the [article online](#) for updates and enhancements.

- The gravitational index and allowable
complex metrics
Pietro Benetti Genolini and Sameer Murthy
# Dynamics of soft interacting particles on a comb

Davide Venturelli¹,²,* ID, Pierre Illien³ ID, Aurélien Grabsch¹ ID and Olivier Bénichou¹

¹ Laboratoire de Physique Théorique de la Matière Condensée, CNRS/Sorbonne Université, 4 Place Jussieu, 75005 Paris, France

² Laboratoire Jean Perrin, CNRS/Sorbonne Université, 4 Place Jussieu, 75005 Paris, France

³ Laboratoire PHENIX (Physico-Chimie des Electrolytes et Nanosystèmes Interfaciaux), CNRS/Sorbonne Université, 4 Place Jussieu, 75005 Paris, France

E-mail: davide.venturelli@cnrs.fr

Received 24 February 2025; revised 23 April 2025

Accepted for publication 7 May 2025

Published 20 May 2025

<figure><img src="image_2.png" /></figure>

CrossMark

## Abstract

We study the dynamics of overdamped Brownian particles interacting through soft pairwise potentials on a comb-like structure. Within the linearized Dean–Kawasaki framework, we characterize the particle density fluctuations by computing their one- and two-point correlation functions. For a tracer particle constrained to move along the comb backbone, we determine the spatial correlation profile between its position and the density of surrounding bath particles. Furthermore, we derive the correction to the diffusion coefficient of the tracer due to interactions with other particles, validating our results through numerical simulations.

Keywords: soft interacting particles, Brownian dynamics, comb-like structures, inhomogeneous media

## Contents

<table><tr><td>1. Introduction</td><td>2</td></tr><tr><td>2. The model</td><td>3</td></tr><tr><td>3. Joint dynamics of the tracer and bath particles density</td><td>5</td></tr><tr><td>3.1. Linearized Dean–Kawasaki equation</td><td>6</td></tr><tr><td>3.2. On the comb</td><td>7</td></tr><tr><td>4. Dynamics of the density fluctuations</td><td>8</td></tr><tr><td>4.1. Propagator</td><td>8</td></tr></table>

<figure><img src="image_3.png" /></figure>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
<table><tbody><tr><td>4.2. Two-point function: general evolution equation</td><td>9</td></tr><tr><td>4.3. Correlator with Dirichlet initial conditions</td><td>10</td></tr><tr><td>4.4. Equilibrium correlator</td><td>10</td></tr><tr><td>4.5. Proof that the stationary and equilibrium correlators coincide</td><td>13</td></tr><tr><td>5. Tracer particle dynamics</td><td>13</td></tr><tr><td>5.1. Generalized correlation profiles</td><td>14</td></tr><tr><td>5.2. Effective diffusion coefficient</td><td>16</td></tr><tr><td>5.3. Validity of the linearized theory</td><td>19</td></tr><tr><td>6. Conclusions</td><td>19</td></tr><tr><td>Data availability statement</td><td>20</td></tr><tr><td>Acknowledgments</td><td>20</td></tr><tr><td>Appendix A. Details on the dynamics of the density fluctuations</td><td>20</td></tr><tr><td>A.1. Derivation of the propagator</td><td>20</td></tr><tr><td>A.2. Derivation of the response function</td><td>22</td></tr><tr><td>A.3. Derivation of the correlator for Dirichlet initial conditions</td><td>23</td></tr><tr><td>A.4. Derivation of the equilibrium correlator</td><td>25</td></tr><tr><td>A.5. Proof of the uniqueness of the stationary solution</td><td>27</td></tr><tr><td>A.6. Double Laplace transform of causal/stationary functions</td><td>29</td></tr><tr><td>A.7. Particles in homogeneous space: analysis of a fluctuating Gaussian field in the Laplace domain</td><td>29</td></tr><tr><td>Appendix B. Structure factor and ISF</td><td>31</td></tr><tr><td>Appendix C. Details on the generalized correlation profiles</td><td>33</td></tr><tr><td>C.1. Derivation of the generating function</td><td>33</td></tr><tr><td>C.2. Derivation and large-distance behavior of correlation profiles</td><td>35</td></tr><tr><td>C.3. On higher-order cumulants</td><td>35</td></tr><tr><td>Appendix D. Brownian dynamics simulations</td><td>36</td></tr><tr><td>References</td><td>37</td></tr></tbody></table>

# 1. Introduction

Particles diffusing in complex environments can undergo anomalous diffusion, characterized by a nonlinear growth with time of their mean squared displacement $\langle x^2(t) \rangle$. The ubiquity of anomalous diffusion in soft and biological systems [1] has, over the past decades, stimulated significant interest in developing minimal theoretical models to capture this behavior [2, 3]. In this context, diffusion of particles in systems with geometrical constraints, such as fractal or disordered lattices, has been the subject of an intense theoretical scrutiny [4, 5], also due to its relevance to the description of transport in porous media, polymer mixtures, and living cells.

Among the simplest inhomogeneous structures is the *comb*: this can be visualized as a line (called the *backbone*) spanning the system from one end to the other, and connected to infinite structures (the teeth), as depicted in figure 1. Comb structures have been originally introduced to represent diffusion in critical percolation clusters [6–8], with the backbone and teeth of the comb mimicking the quasi-linear structure and dead ends of the clusters, respectively. In general, one expects that a particle will spend a long time exploring a tooth, resulting in a sub-diffusive motion along the backbone—e.g. $\langle x^2(t) \rangle \sim t^{1/2}$ in two spatial dimensions [9]. Since then, comb-like models have been used to describe real systems such as cancer proliferation [10] and dendronized polymers [11], transport in spiny dendrites [12], and the diffusion of cold atoms [13] or in crowded media [1].
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Schematic representation of interacting Brownian particles evolving on a comb structure in $d=2$. Bath particles (blue) can move horizontally along the backbone, or vertically along the teeth. Note that, in the continuum model considered here, the teeth are continuously distributed along the horizontal direction (in the sketch they are instead represented as spatially separated, only for graphical purposes). The tracer particle (red) is constrained to diffuse only along the backbone.</figcaption></figure>

From a theoretical perspective, the single-particle comb model is sufficiently simple to allow for the derivation of several exact analytical results, both in the continuum [14–19], and on the lattice [20–25]—see also [26] for an overview. By contrast, the many-body problem, consisting of interacting particles evolving on a comb, has received much less attention, and so far restricted to the case of hard-core lattice gases [27–30]. The natural step forward, which we aim to address here, is to explore the dynamics of interacting particles in continuum space, subject to the comb constraint.

To this end, in this work we consider a system of overdamped Brownian particles interacting via soft pairwise potentials, and constrained to move on a comb, as described in section 2. Using the Dean–Kawasaki formalism [31, 32], we first derive in section 3 the exact equations that describe the dynamics of the particle density field $\rho(\mathbf{x},t)$. Expanding the latter around a constant background density, in section 4 we then derive the one- and two-point functions that completely characterize the density fluctuations within the Gaussian approximation. To the best of our knowledge, this represents the first application of the Dean–Kawasaki theory to non-homogeneous media [33, 34]. Next, in section 5, we single out a tagged tracer from the bath of interacting particles, assuming that it is constrained to move only along the comb backbone. In this setting, we first derive the spatial correlation profiles that describe the cross-correlations between the tracer position and the surrounding bath density, and then use them to estimate the effective diffusion coefficient of the tracer, which we finally test using Brownian dynamics simulations.

## 2. The model

We consider a system of $(N+1)$ Brownian particles at positions $\mathbf{r}_i(t)$, which evolve according to the overdamped Langevin equations

$$ \dot{\mathbf{r}}_i(t) = -\hat{\mu}^{(i)}(\mathbf{r}_i(t)) \sum_{j \neq i} \nabla_{\mathbf{r}_i} U_{ij}(\mathbf{r}_i(t) - \mathbf{r}_j(t)) + \eta_i(\mathbf{r}_i(t), t). \quad (1) $$

Here, the Gaussian noises $\eta_i$ satisfy $\langle \eta_i(\mathbf{x}, t) = 0 \rangle$ and have variance

$$ \langle \eta_i(\mathbf{x}, t) \eta_j^\mathrm{T}(\mathbf{x}', t') \rangle = 2k_B T \hat{\mu}^{(i)}(\mathbf{x}) \delta_{ij} \delta(\mathbf{x} - \mathbf{x}') \delta(t - t'), \quad (2) $$
so that the fluctuation-dissipation theorem is satisfied (henceforth we will set the Boltzmann constant $k_B$ to unity, for simplicity). For future convenience, we single out the particle with $i=0$, which we denote here and henceforth as the tracer. We allow for the particles' mobility to generically depend on their position in space, and we assume that the tracer can have a different mobility with respect to that of the other bath particles; hence, the mobility matrices in equation (1) read

$$ \hat{\mu}^{(i)}(\mathbf{x}) = \begin{cases} \hat{\mu}^{(0)}(\mathbf{x}), & i=0, \\ \hat{\mu}(\mathbf{x}), & i \ge 1. \end{cases} \qquad (3) $$

Finally, the inter-particle interaction potentials in equation (1) read

$$ U_{ij}(\mathbf{x}) = \begin{cases} U_0(\mathbf{x}), & i=0 \text{ or } j=0, \\ U(\mathbf{x}), & i,j \ge 1, \end{cases} \qquad (4) $$

between tracer-bath or bath-bath particles, respectively.

So far, the features of the underlying space in which the particles are evolving have not yet been specified. To enforce the comb geometry, which is the subject of this work (see figure 1), one can introduce in equation (1) the anisotropic mobility matrices<fnref n="4" />

$$ \hat{\mu}(\mathbf{x}) = \mu \begin{pmatrix} \delta(x_2) \dots \delta(x_d) & 0 \\ 0 & 1_{d-1} \end{pmatrix} = \mu \begin{pmatrix} \delta(\mathbf{x}_{\perp}) & 0 \\ 0 & 1_{d-1} \end{pmatrix}, \quad \text{and} \quad \hat{\mu}^{(0)}(\mathbf{x}) = \mu_0 \begin{pmatrix} 1 & 0 \\ 0 & 0_{d-1} \end{pmatrix}. \qquad (5) $$

This constrains the tracer particle to move only on the backbone (which we chose without loss of generality to be oriented along the Cartesian direction $\hat{e}_1$), while all other bath particles can move along the backbone only if $\mathbf{x}_{\perp} = \mathbf{0}$. In $d=2$, equation (5) reduces to

$$ \hat{\mu}(\mathbf{x}) = \mu \begin{pmatrix} \delta(y) & 0 \\ 0 & 1 \end{pmatrix}, \quad \text{and} \quad \hat{\mu}^{(0)}(\mathbf{x}) = \mu_0 \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix}, \qquad (6) $$

corresponding to the comb schematically represented in figure 1. For $d \ge 3$, the resulting structure may be rather visualized as a backbone sliced by orthogonal hyperplanes.

Note that, for nontrivial choices of $\hat{\mu}^{(i)}(\mathbf{x})$ in equation (2), the Gaussian noises $\eta_i$ are white and multiplicative—to fix their physical meaning, we adopt here the Itô prescription. Choosing $\hat{\mu}(\mathbf{x})$ as in equation (6), it is then simple to check that the Langevin equation (1) reduces, in the single-particle case, to

$$ \partial_t P(\mathbf{r}, t) = \mu T [a\delta(y) \partial_x^2 + \partial_y^2] P(\mathbf{r}, t), \qquad (7) $$

which is the Fokker-Planck equation for a Brownian particle at position $\mathbf{r} = (x, y)$ on a $2d$ comb, as commonly adopted in the literature [9, 26] (and where we have replaced $\delta(y) \mapsto a\delta(y)$, where $a$ is a length scale—see the previous footnote). From the well-known solution [9] of equation (7) (recalled in appendix A, see cf equation (A.7)), one finds at long times $t$ that

<footnote n="4">⁴ We stress that, by dimensional consistency, the delta functions that appear in the expressions above should be replaced by $\delta(x_i) \mapsto a\delta(x_i)$, where $a$ is a length scale. In the following, we will set this length scale to unity for simplicity; indeed, note that $a$ can in any case be reabsorbed via a suitable rescaling of the corresponding Fourier variable $q_1$ (see cf section 4), e.g. $q_1^2 \mapsto q_1^2 a^{d-1}$.</footnote>
$\langle y^2(t) \rangle \approx 2\mu Tt$, while $\langle x^2(t) \rangle \approx 2a\sqrt{\mu Tt/\pi}$. Thus, a free Brownian particle on a comb dif-fuses along the teeth, but subdiffuses along the backbone.

Crucially, in the presence of interactions or external potentials [19], inserting $\delta(y)$ only in the diffusion coefficient (as done in equation (7) for the free-particle case) would not prevent the particles from leaving the comb. This is why the term $\delta(y)$ has been included also in the particle's mobility (5).

In the following section, we set out to compute, starting from the set of Langevin equations (1), the evolution equation satisfied by the fluctuating particle density.

### 3. Joint dynamics of the tracer and bath particles density

The dynamics of the microscopic particle density can be derived by standard methods within the Dean-Kawasaki formalism, also known as stochastic density functional theory (SDFT) [31, 32, 34]. To this end, here we generalize the derivation of [35] to the case in which the particles' mobility encodes geometrical constraints, as in equation (5). Following [31, 35], we first introduce the fluctuating density of bath particles

$$ \rho(\mathbf{x}, t) = \sum_{i=1}^{N} \delta(\mathbf{x} - \mathbf{r}_i(t)), \qquad (8) $$

in terms of which the evolution equation of the tracer (i.e. $i=0$) becomes

$$ \dot{\mathbf{r}}_0(t) = -\hat{\mu}^{(0)}(\mathbf{r}_0(t)) \nabla_{\mathbf{r}_0} \int dy U_0(\mathbf{r}_0(t) - \mathbf{y}) \rho(\mathbf{y}, t) + \eta_0(\mathbf{r}_0(t), t), \qquad (9) $$

while the bath density can be shown to follow the Dean equation [31, 32]

$$ \partial_t \rho(\mathbf{x}, t) = \nabla \cdot \left[ \hat{\mu}(\mathbf{x}) \rho(\mathbf{x}, t) \nabla \frac{\delta \mathcal{F}}{\delta \rho(\mathbf{x}, t)} \right] + \nabla \cdot \left[ \rho^{\frac{1}{2}}(\mathbf{x}, t) \boldsymbol{\xi}(\mathbf{x}, t) \right], \qquad (10) $$

where $\boldsymbol{\xi}$ is a Gaussian noise field with zero mean and correlations

$$ \langle \boldsymbol{\xi}(\mathbf{x}, t) \boldsymbol{\xi}^{\mathrm{T}}(\mathbf{x}', t') \rangle = 2T\hat{\mu}(\mathbf{x}) \delta(\mathbf{x} - \mathbf{x}') \delta(t - t'), \qquad (11) $$

and where $\mathcal{F}$ takes the form of a pseudo free energy,

$$ \mathcal{F}[\rho, \mathbf{r}_0] = T \int d\mathbf{x} \rho(\mathbf{x}) \log \left( \frac{\rho(\mathbf{x})}{\rho_0} \right) + \frac{1}{2} \int d\mathbf{x} d\mathbf{y} \rho(\mathbf{x}) U(\mathbf{x} - \mathbf{y}) \rho(\mathbf{y}) \\ + \int d\mathbf{y} \rho(\mathbf{y}) U_0(\mathbf{y} - \mathbf{r}_0), \qquad (12) $$

with $\rho_0$ being a uniform background density.

By choosing $U_0(\mathbf{x} - \mathbf{y}) = U_0(\mathbf{y} - \mathbf{x})$—i.e. by assuming that the interaction between the tracer and the other bath particles is reciprocal—then we may also rewrite

$$ \dot{\mathbf{r}}_0(t) = -\hat{\mu}^{(0)}(\mathbf{r}_0(t)) \nabla_{\mathbf{r}_0} \mathcal{F}[\rho, \mathbf{r}_0] + \eta_0(\mathbf{r}_0(t), t), \qquad (13) $$

so that the system formally admits the joint stationary distribution

$$ \mathcal{P}[\rho, \mathbf{r}_0] \propto e^{-\frac{1}{T} \mathcal{F}[\rho, \mathbf{r}_0]}. \qquad (14) $$
Crucially, however, reaching this solution dynamically requires ergodicity, and this is not generically granted in the presence of geometrical constraints—such as those encoded in $\hat{\mu}(\mathbf{x})$ for the comb (see section 3.2). Indeed, a simple counter-example can be produced by choosing $\hat{\mu}(\mathbf{x}) = \text{diag}(1, 0)$ in $d=2$: clearly, the system cannot relax along the $\hat{y}$ direction, because there are no forces or noise along that direction. This point will be further discussed in section 4.5.

### 3.1. Linearized Dean-Kawasaki equation

The derivation is so far exact; however, solving equation (10) is challenging due to the presence of nonlinearities in both the deterministic and the stochastic term. To make progress, here we linearize it assuming small bath density fluctuations around a uniform background $\rho_0$ [35]. This strategy is by now quite standard in the physics literature [34–43], and has allowed (in spite of the highly singular nature of the Dean-Kawasaki equation [44]) to derive approximate but accurate predictions for several distinct systems of particles interacting in uniform geometries, under suitable assumptions (discussed below and in section 5.3). To this end, we consider

$$ \rho(\mathbf{x}, t) = \rho_0 + \rho_0^{\frac{1}{2}} \phi(\mathbf{x}, t), \qquad (15) $$

and plug it back into equation (10); we then discard small terms<fnref n="5" /> according to the approximation $\rho_0^{-1/2} \phi \ll 1$. The result reads

$$ \dot{\mathbf{r}}_0(t) = -\rho_0^{-\frac{1}{2}} \hat{\mu}^{(0)}(\mathbf{r}_0(t)) \nabla_{\mathbf{r}_0} \int dy v(\mathbf{r}_0(t) - \mathbf{y}) \phi(\mathbf{y}, t) + \eta_0(\mathbf{r}_0(t), t), \qquad (16) $$

$$ \partial_t \phi(\mathbf{x}, t) = \nabla \cdot \left\{ \hat{\mu}(\mathbf{x}) \nabla \left[ T\phi(\mathbf{x}, t) + \int dy u(\mathbf{x}-\mathbf{y}) \phi(\mathbf{y}, t) + \rho_0^{-\frac{1}{2}} v(\mathbf{x}-\mathbf{r}_0) \right] + \xi(\mathbf{x}, t) \right\}, \qquad (17) $$

where we rescaled the interaction potentials as

$$ v(\mathbf{x}) = \rho_0 U_0(\mathbf{x}), \quad u(\mathbf{x}) = \rho_0 U(\mathbf{x}). \qquad (18) $$

Although this approximation becomes formally exact in the dense limit (i.e. upon sending $\rho_0 \to \infty$ by keeping the product $\rho_0 V$ fixed, with $V$ being the volume of the system), it has proven to be accurate even away from the dense limit provided that the interaction potentials $U(\mathbf{x})$ and $U_0(\mathbf{x})$ are sufficiently soft, so that particles can overlap completely at a finite energy cost due to thermal fluctuations [35, 42]. This point will be thoroughly discussed in section 5.3. Note that we can still introduce a pseudo free energy for the linearized system,

$$ \mathcal{F}_2[\phi, \mathbf{r}_0] = \frac{T}{2} \int d\mathbf{x} \phi^2(\mathbf{x}) + \frac{1}{2} \int d\mathbf{x} dy \phi(\mathbf{x}) u(\mathbf{x}-\mathbf{y}) \phi(\mathbf{y}) + \rho_0^{-1/2} \int dy \phi(\mathbf{y}) v(\mathbf{y}-\mathbf{r}_0), \qquad (19) $$

so that equation (17) can be rewritten as

$$ \partial_t \phi(\mathbf{x}, t) = \nabla \cdot \left[ \hat{\mu}(\mathbf{x}) \nabla \frac{\delta \mathcal{F}_2}{\delta \phi(\mathbf{x}, t)} + \xi(\mathbf{x}, t) \right]. \qquad (20) $$

<footnote n="5">⁵ Note that in fact we are not simply linearizing the equation. For instance, the term $\mu \rho_0^{-1} \nabla \cdot \phi(\mathbf{x}, t) \nabla v(\mathbf{x}-\mathbf{r}_0)$ turns out to be suppressed by this approximation, although it is linear in $\phi$. However, we adhere here to this commonly adopted terminology.</footnote>
As we did for equation (13), if $v(\mathbf{x}) = v(-\mathbf{x})$ we can also rewrite equation (16) as

$$ \dot{\mathbf{r}}_0(t) = -\rho_0^{-1/2} \hat{\mu}^{(0)}(\mathbf{r}_0(t)) \nabla_{\mathbf{r}_0} \mathcal{F}_2[\phi, \mathbf{r}_0] + \eta_0(t), \quad (21) $$

and similar considerations apply.

Our derivation is so far quite general, in that the spatial geometry is encoded in the mobility matrices $\hat{\mu}$ and $\hat{\mu}^{(0)}$ which have not yet been specified. In the following section (and throughout the rest of this paper), we specialize our calculation to the comb geometry.

## 3.2. On the comb

Choosing the mobilities as in equation (5), the linearized Dean–Kawasaki equation (16) for the tracer particle specializes to

$$ \dot{r}_0(t) = -\rho_0^{-1/2} \mu_0 \partial_{r_0} \int dy v(\mathbf{r}_0(t) - \mathbf{y}) \phi(\mathbf{y}, t) + \eta_0(t), \quad (22) $$

where $\eta_0$ is a scalar Gaussian noise with zero mean and variance

$$ \langle \eta_0(t) \eta_0(t') \rangle = 2T\mu_0\delta(t-t'). \quad (23) $$

Indeed, since the tracer can only move along the backbone, its position can be specified by using a single scalar coordinate $r_0(t)$, i.e. $\mathbf{r}_0(t) = r_0(t)\hat{\mathbf{e}}_1$. In the following sections, we will focus first on the evolution of the density $\phi$ alone—namely, we will temporarily switch off the interaction with the tracer in equation (17) by setting $v(\mathbf{x}) = 0$, to obtain

$$ \partial_t \phi(\mathbf{x}, t) = \nabla \cdot \left\{ \hat{\mu}(\mathbf{x}) \nabla \left[ T\phi(\mathbf{x}, t) + \int dy u(\mathbf{x}-\mathbf{y}) \phi(\mathbf{y}, t) \right] + \xi(\mathbf{x}, t) \right\}, \quad (24) $$

where the variance of the noise $\xi$ was given in equation (11). This will give access to the dynamical space-time correlations of the density field of interacting bath particles, independ- ently of our initial choice of singling out a tracer particle. As usual, note that we can rewrite equation (24) as

$$ \partial_t \phi(\mathbf{x}, t) = \nabla \cdot \left[ \hat{\mu}(\mathbf{x}) \nabla \frac{\delta \mathcal{H}}{\delta \phi(\mathbf{x}, t)} + \xi(\mathbf{x}, t) \right], \quad (25) $$

with

$$ \mathcal{H}[\phi] \equiv \frac{T}{2} \int d\mathbf{x} \phi^2(\mathbf{x}) + \frac{1}{2} \int d\mathbf{x} dy \phi(\mathbf{x}) u(\mathbf{x}-\mathbf{y}) \phi(\mathbf{y}). \quad (26) $$

In turn, the dynamical propagator and correlator of $\phi$ generally serve as building blocks of the perturbation theory in the presence of a tracer [35, 45–49], whose resulting dynamics will be the subject of section 5. Note that the theory in equations (25) and (26) is quadratic in the density field $\phi$, yet its analysis is made difficult by the presence of the $\delta(\mathbf{x}_{\perp})$ term in equation (5), which manifestly breaks its spatial translational invariance. In the next section, we will show how this difficulty can be overcome by resorting to a self-consistent approach.
## 4. Dynamics of the density fluctuations

The evolution equation (24) for the density field reads in Fourier space

$$ \partial_t \phi_q(t) = -\mu q_{\perp}^2 (T + u_q) \phi_q(t) - \mu q_1^2 \int \frac{d^{d-1}p_{\perp}}{(2\pi)^{d-1}} (T + u_{q_1, -\mathbf{p}_{\perp}}) \phi_{q_1, \mathbf{p}_{\perp}}(t) + \nu_q(t), \quad (27) $$

where we introduced the Fourier transform $\nu_q(t)$ of the scalar Gaussian noise

$$ \nu(\mathbf{x}, t) \equiv \nabla \cdot \xi(\mathbf{x}, t), \quad (28) $$

with correlations (see equation (11))

$$ \langle \nu_q(t) \nu_p(t') \rangle = 2\mu T \delta(t-t') [q_1^2 \delta(q_1+p_1) + q_{\perp}^2 \delta^d(\mathbf{q}+\mathbf{p})]. \quad (29) $$

Here we have used the Fourier convention $f_q(t) = \int d^d x e^{-i\mathbf{q} \cdot \mathbf{x}} f(\mathbf{x}, t)$, and normalized the Dirac delta in Fourier space as $\int \frac{d^d q}{(2\pi)^d} \delta^d(\mathbf{q}) = 1$. Below, we will frequently need to adopt mixed representations involving the Fourier or Laplace transform (often with respect to only a few among the spatial and temporal variables on which the observables depend). With a slight abuse of notation, whenever not ambiguous, we will adopt the same symbol (e.g. $f$) for all these transforms, to avoid the proliferation of hat symbols (such as $\hat{f}, \tilde{f}$ and so on). As a general convention, we will adopt $\mathbf{x}, \mathbf{y}, \mathbf{z}$ for physical space, $\mathbf{q}, \mathbf{p}, \mathbf{k}$ for Fourier momenta, $t, t'$ for time, and $s, s'$ for Laplace variables.

To characterize the fluctuations of the particle density $\phi(\mathbf{x}, t)$, we now set out to compute its one- and two-point functions. This extends to the case of the comb geometry previous analyses conducted on Brownian particles evolving in homogeneous space [35, 42, 45]—see also appendix A.7, where some of these expressions are reproduced for convenience.

### 4.1. Propagator

We start by addressing the propagator, i.e. the solution of

$$ \partial_t \phi(\mathbf{x}, t) = \nabla \cdot \hat{\mu}(\mathbf{x}) \nabla \frac{\delta \mathcal{H}}{\delta \phi(\mathbf{x}, t)} + \delta(\mathbf{x}, t), \quad (30) $$

or equivalently in Fourier space

$$ \partial_t G_q(t) = -\mu q_{\perp}^2 (T + u_q) G_q(t) - \mu q_1^2 \int \frac{d^{d-1}p_{\perp}}{(2\pi)^{d-1}} (T + u_{q_1, -\mathbf{p}_{\perp}}) G_{q_1, \mathbf{p}_{\perp}}(t) + \delta(t). \quad (31) $$

We solve this equation in appendix A.1; the key step is to realize that, by applying the operator $\int \frac{d^{d-1}q_{\perp}}{(2\pi)^{d-1}} (T + u_{q_1, -q_{\perp}})$ to both members of equation (31), we can generate a closed self-consistency relation, which then gives access to the solution. It proves convenient to write the result in the Fourier–Laplace domain, i.e. $G_q(s) = \int_0^{\infty} dt e^{-st} G_q(t)$, where it reads

$$ G_q(s) = g_{\perp}(\mathbf{q}, s) g_1(q_1, s), \quad (32) $$

with

$$ g_{\perp}^{-1}(\mathbf{q}, s) \equiv s + \mu q_{\perp}^2 (T + u_q), \quad (33) $$
$$g_1^{-1}(q_1, s) = 1 + \mu q_1^2 \int \frac{d^{d-1}q_{\perp}}{(2\pi)^{d-1}} g_{\perp}(\mathbf{q}, s) (T + u_{q_1, -\mathbf{q}_{\perp}}). \quad (34)$$

As expected, by setting $u_q = 0$ we correctly recover the moment generating function of a single Brownian particle on the comb, starting at the origin at time $t = 0$—see appendix A.1.

Note that, for translationally invariant differential equations, the propagator coincides with the Green's function (i.e. the solution of equation (30) in which $\delta(\mathbf{x})$ is replaced by $\delta(\mathbf{x} - \mathbf{x}')$), being the latter a function of $\mathbf{x} - \mathbf{x}'$ only [30]. By contrast, here the system is only translationally invariant along the backbone $\mathbf{x}_{\perp} = \mathbf{0}$; accordingly, the solution of

$$\partial_t \phi(\mathbf{x}, t) = \nabla \cdot \hat{\mu}(\mathbf{x}) \nabla \frac{\delta \mathcal{H}}{\delta \phi(\mathbf{x}, t)} + f(\mathbf{x}, t), \quad (35)$$

in the presence of a generic source term $f(\mathbf{x}, t)$, turns out to read

$$\phi_q(s) = g_{\perp}(\mathbf{q}, s)f_q(s) - \mu q_1^2 G_q(s) \int \frac{d^{d-1}q_{\perp}}{(2\pi)^{d-1}} (T + u_{-\mathbf{q}}) g_{\perp}(\mathbf{q}, s)f_q(s), \quad (36)$$

which only reduces to $\phi_q(s) = G_q(s)f_q(s)$ if $f(\mathbf{x}, t) \propto \delta(\mathbf{x}_{\perp})$.

## 4.2. Two-point function: general evolution equation

Computing the two-point function

$$C_{q,p}(t, t') = \langle \phi_q(t) \phi_p(t') \rangle \quad (37)$$

proves more challenging. To this end, let us first introduce the response function

$$R_{q,p}(t, t') = \langle \phi_q(t) \nu_p(t') \rangle. \quad (38)$$

Upon multiplying the evolution equation (27) by $\phi_p(t')$ or $\nu_p(t')$ and taking their average, we obtain the coupled Schwinger-Dyson equations [50, 51]

$$\partial_t C_{q,p}(t, t') = -\mu q_{\perp}^2 (T + u_q) C_{q,p}(t, t') - \mu q_1^2 \int \frac{d^{d-1}k_{\perp}}{(2\pi)^{d-1}} (T + u_{q_1, -k_{\perp}}) C_{(q_1, k_{\perp}), p}(t, t') \\ + R_{p,q}(t', t), \quad (39)$$

$$\partial_t R_{q,p}(t, t') = -\mu q_{\perp}^2 (T + u_q) R_{q,p}(t, t') - \mu q_1^2 \int \frac{d^{d-1}k_{\perp}}{(2\pi)^{d-1}} (T + u_{q_1, -k_{\perp}}) R_{(q_1, k_{\perp}), p}(t, t') \\ + 2\mu T \delta(t - t') [q_1^2 \delta(q_1 + p_1) + q_{\perp}^2 \delta^d(q + p)]. \quad (40)$$

We note that, while $C_{q,p}(t, t')$ is enslaved to $R_{q,p}(t, t')$, by contrast equation (40) provides a closed equation for $R_{q,p}(t, t')$. We thus solve first equation (40) self-consistently in the Laplace domain, as detailed in appendix A.2, to obtain

$$R_{q,p}(s, s') = \frac{2\mu T}{s+s'} g_{\perp}(\mathbf{p}, s) [q_{\perp}^2 \delta^d(\mathbf{q} + \mathbf{p}) + q_1^2 \delta(q_1 + p_1) s G_q(s)]. \quad (41)$$

Here and henceforth, we have assumed for simplicity the interaction potential $u(\mathbf{x})$ to be rotationally invariant. The response function in equation (41) is symmetric under $(\mathbf{p} \leftrightarrow \mathbf{q})$, which was not obvious a priori, while the lack of symmetry ($s \leftrightarrow s'$) is expected on causality grounds.
Indeed, invoking the properties of the double Laplace transform of causal functions summarized in appendix A.6, we recognize in equation (41) the structure

$$R_{q,p}(s, s') = \frac{\tilde{R}_{q,p}(s)}{s + s'}, \qquad (42)$$

corresponding in the time domain (upon taking the inverse Laplace transform with respect to both s and s') to

$$R_{q,p}(t, t') = \tilde{R}_{q,p}(t - t') \Theta(t - t'). \qquad (43)$$

### 4.3. Correlator with Dirichlet initial conditions

The expression of $R_{q,p}(s, s')$ found in equation (41) can now be inserted into equation (39) to obtain a closed equation for the correlator, which we solve self-consistently in appendix A.3 for the case of Dirichlet initial conditions $\phi_q(t=0) = \phi_q^{(0)}$. Recalling that $\phi(x, t)$ represents the fluctuation with respect to a uniform background density $\rho_0$ (see equation (15)), this amounts to imposing quenched initial conditions $\rho(x, t=0) = \rho_0 + \rho_0^{1/2} \phi^{(0)}(x)$ on the particle density (not necessarily flat). The corresponding correlator reads

$$C_{q,p}(s, s') = \frac{2\mu T}{s+s'} g_{\perp}(\mathbf{q}, s) g_{\perp}(\mathbf{p}, s') \left\{ \begin{aligned} & p_{\perp}^2 \delta^d(\mathbf{q}+\mathbf{p}) + q_1^2 \delta(q_1+p_1) \left[ s' G_{\mathbf{q}}(s') + s G_{\mathbf{p}}(s) \right. \\ & \left. \qquad + \frac{s' g_1(q_1, s') - s g_1(p_1, s)}{s-s'} \right] \\ & + \phi_{\mathbf{p}}^{(0)} G_{\mathbf{p}}(s') g_{\perp}(\mathbf{q}, s) \left[ \phi_{\mathbf{q}}^{(0)} - \mu q_1^2 g_1(q_1, s) M^{(0)}(q_1, s) \right] \end{aligned} \right\}, \qquad (44)$$

with

$$M^{(0)}(q_1, s) = \int \frac{d^{d-1}q_{\perp}}{(2\pi)^{d-1}} (T + u_{q_1, -\mathbf{q}_{\perp}}) g_{\perp}(\mathbf{q}, s) \phi_{\mathbf{q}}^{(0)}. \qquad (45)$$

In the case of flat initial conditions, i.e. $\phi_q^{(0)} \equiv \phi^{(0)}$, it can be shown that the last term in equation (44) simply reduces to (see appendix A.3)

$$[\phi^{(0)}]^2 G_{\mathbf{q}}(s) G_{\mathbf{p}}(s'), \qquad (46)$$

while the first term represents the connected part of the correlation function, i.e. $\langle \phi_q(s) \phi_p(s') \rangle_c$. Note that the so-obtained correlator is symmetric under the exchange of $(\mathbf{q}, s) \leftrightarrow (\mathbf{p}, s')$, as expected from its definition in equation (37), but it is not diagonal in the momenta $\mathbf{q}$ and $\mathbf{p}$ (i.e. it is not proportional to $\delta^d(\mathbf{q} + \mathbf{p})$), because of the spatial anisotropy of the system. However, a common dependence on $\delta(q_1 + p_1)$ remains, as a consequence of the translational invariance along the spatial direction parallel to the backbone.

### 4.4. Equilibrium correlator

At long times, we expect the two-time correlation function to attain a stationary form such that

$$C_{q,p}(t, t') = C_{q,p}(|t - t'|). \qquad (47)$$
Using the properties of the double Laplace transform summarized in appendix A.6, this corresponds to

$$C_{q,p}(s, s') = \frac{C_{q,p}(s) + C_{q,p}(s')}{s + s'}, \qquad (48)$$

for some function $C_{q,p}(s)$ that we now set out to compute. (By contrast, the Dirichlet correlator given in equation (44) does not have a time-translational invariant structure, as expected). The simplest way to proceed is to invoke the fluctuation-dissipation theorem that relates the correlation function at equilibrium to the linear susceptibility [52]

$$\chi(\mathbf{x}, \mathbf{y}, t, t') \equiv \left. \frac{\delta \langle \phi(\mathbf{x}, t) \rangle}{\delta j(\mathbf{y}, t')} \right|_{j=0}. \qquad (49)$$

The latter is obtained by adding to the Hamiltonian $\mathcal{H}[\phi]$ in equation (26) a source field $j(\mathbf{x}, t)$ linearly coupled to $\phi(\mathbf{x}, t)$ as

$$\mathcal{H}_j[\phi] \equiv \mathcal{H}[\phi] - \int d\mathbf{x} j(\mathbf{x}) \phi(\mathbf{x}), \qquad (50)$$

and by modifying the equation of motion (24) and its solution $\langle \phi(\mathbf{x}, t) \rangle$ accordingly.

As detailed in appendix A.4, one can show that the susceptibility $\chi$ follows an evolution equation very similar to equation (40) satisfied by the response function $R$, and that the two are in fact simply proportional to one another:

$$R_{q,p}(t, t') = 2T\chi_{q,p}(t, t'). \qquad (51)$$

This is expected, since the noise $\nu_q(t)$ enters linearly in the equation of motion (27). Like the response function, the linear susceptibility is also causal and time-translational invariant, that is

$$\chi_{q,p}(t, t') = \tilde{\chi}_{q,p}(t-t') \Theta(t-t'). \qquad (52)$$

By analogy with equations (41) and (42), we can thus write in the Laplace domain

$$\tilde{\chi}_{q,p}(s) = \mu g_{\perp}(p_{\perp}, s) [q_{\perp}^2 \delta^d(\mathbf{q} + \mathbf{p}) + q_1^2 \delta(q_1 + p_1) s G_q(s)]. \qquad (53)$$

Next, we recall the fluctuation-dissipation theorem [52]

$$\frac{d}{d\tau} C_{q,p}(\tau) = -T \tilde{\chi}_{q,p}(\tau), \quad \text{for } \tau \equiv t - t' > 0, \qquad (54)$$

or equivalently, in the Laplace domain,

$$C_{q,p}(s) = \frac{1}{s} [C_{q,p}(\tau = 0^+) - T\chi_{q,p}(s)]. \qquad (55)$$

The stationary fluctuations $C_{q,p}(\tau = 0^+)$ of the field can be accessed by noting that its evolution equation (25) admits the equilibrium distribution (see the discussion in section 3)

$$\mathcal{P}[\phi] \propto e^{-\frac{1}{2}\mathcal{H}[\phi]}, \qquad (56)$$
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Scaling form assumed at long times $t$ by the intermediate scattering function $\tilde{\mathcal{F}}(\mathbf{q}, t)$—see section 4.4 and appendix B for the details. The amplitude of the ISF is zero in correspondence of the dark-blue regions, and maximum in correspondence of the red regions. The plot is obtained for a system of size $V = 100$ in $d=2$, where the combination $\mu(T + u_{q=0})$ and the maximum amplitude of the ISF were normalized to unity for illustrative purposes.</figcaption></figure>

whence (see appendix A.4, and compare e.g. with equation (43) in [45])

$$C_{q,p}(\tau = 0^+) = \langle \phi_q(t) \phi_p(t) \rangle_{t \to \infty} = \frac{T}{T + u_q} \delta^d(\mathbf{q} + \mathbf{p}). \qquad (57)$$

Using equation (55) then gives, after some straightforward algebra$^6$,

$$C_{q,p}(s) = T_{g_{\perp}}(\mathbf{p}, s) \left[ \frac{\delta^d(\mathbf{q} + \mathbf{p})}{T + u_q} - \mu q_1^2 \delta(q_1 + p_1) G_q(s) \right]. \qquad (58)$$

This expression admits a closed-form analogous in the time domain, for simple choices of $u_q$—see appendix A.4.

The stationary correlator in equation (58), obtained within the linearized SDFT, is one of the key results of this work: it encodes all dynamical two-point correlations of the fluctuating particle density. Various dynamical observables, such as structure factors and first passage times, can be accessed starting from equation (58); as an application, here we consider the intermediate scattering function (ISF)

$$\tilde{\mathcal{F}}(\mathbf{q}, t) = \frac{1}{N} \langle \rho_{\mathbf{q}}(t) \rho_{-\mathbf{q}}(0) \rangle - \rho_0 \delta^d(\mathbf{q}) = \frac{\rho_0}{N} \langle \phi_{\mathbf{q}}(t) \phi_{-\mathbf{q}}(0) \rangle = \frac{\rho_0}{N} C_{\mathbf{q},-\mathbf{q}}(t). \qquad (59)$$

As we demonstrate in appendix B, for sufficiently long times $t$ the ISF assumes a scaling form $\tilde{\mathcal{F}}(q_1 t^{\frac{1}{4}}, q_{\perp} t^{\frac{1}{2}})$ that is independent of the specifics of the interaction potential $u(\mathbf{x})$, and which we compute explicitly for the case $d=2$ (see equation (B.9)). The resulting scaling function is plotted in figure 2 for selected choices of the various parameters, and with its amplitude rescaled to unity for graphical convenience. Such ISF is manifestly reminiscent of the anisotropy of the comb structure.

<footnote>$^6$ As a consistency check, one can integrate equation (58) along the backbone by setting $q_1 = 0$. As expected, one recovers this way the known expression of the two-point function for particles diffusing in homogeneous space—see the discussion in appendix A.7.</footnote>
### 4.5. Proof that the stationary and equilibrium correlators coincide

In the previous section we computed the equilibrium correlator, see equation (58). However, as we stressed under equation (14), there is no reason to believe a priori that the stationary state of the system is also an equilibrium one. To this end, here we prove that the solution obtained in section 4.4 is unique—namely, that no other function exists for this system with the property of time-translational invariance required in equation (47).

We start by noting that an alternative and more direct way to impose the time-translational invariance $C_{q,p}(t,t') = C_{q,p}(t-t')$ of the correlator (without invoking the fluctuation-dissipation theorem, which holds only at equilibrium) is to require

$$ (\partial_t + \partial_{t'}) C_{q,p}(t,t') = 0. \qquad (60) $$

This condition can be used to construct the stationary correlator as follows. First, we change indices in the field equation (27) and write an evolution equation for $\partial_{t'} \phi_p(t')$. Next, we multiply it by $\phi_q(t)$ and take its average over the noise. This gives another Schwinger-Dyson equation,

$$ \partial_{t'} C_{q,p}(t,t') = -\mu p_{\perp}^2 (T+u_p) C_{q,p}(t,t') - \mu p_1^2 \int \frac{d^{d-1}k_{\perp}}{(2\pi)^{d-1}} (T+u_{p_1,-k_{\perp}}) C_{q,(p_1,k_{\perp})}(t,t') \\ + R_{q,p}(t,t'), \qquad (61) $$

akin to equation (39), which can then be summed to the latter. Setting the result equal to zero as per equation (60), and upon introducing

$$ \mathcal{G}_{q,p}^{-1} \equiv \mu [q_{\perp}^2 (T+u_q) + p_{\perp}^2 (T+u_p)], \qquad (62) $$

$$ S_{q,p}(t,t') \equiv R_{q,p}(t,t') + R_{p,q}(t',t), \qquad (63) $$

we find that the stationary correlator $C_{q,p}(t,t')$ must satisfy

$$ C_{q,p}(t,t') = \mathcal{G}_{q,p} S_{q,p}(t,t') - \mu \mathcal{G}_{q,p} \left[ q_1^2 \int \frac{d^{d-1}k_{\perp}}{(2\pi)^{d-1}} (T+u_{q_1,-k_{\perp}}) C_{(q_1,k_{\perp}),p}(t,t') \right. \\ \left. + p_1^2 \int \frac{d^{d-1}k_{\perp}}{(2\pi)^{d-1}} (T+u_{p_1,-k_{\perp}}) C_{q,(p_1,k_{\perp})}(t,t') \right]. \qquad (64) $$

This integral equation is analyzed in appendix A.5; there, we demonstrate that in fact equation (64) admits a unique solution, corresponding to the equilibrium correlator found in equation (58). This concludes our analysis of the one- and two-point functions within the linearized Dean-Kawasaki theory on the comb, which completely characterize (within this Gaussian approximation) the particle density fluctuations.

## 5. Tracer particle dynamics

In section 4 we have analyzed the fluctuations of the bath density, assuming that all bath particles are equal. In this section, instead, we go back to equations (16) and (17), where a tagged tracer (singled out from the particles bath) evolves in contact with the fluctuating density of surrounding particles. The resulting coupled equations for the tracer and the bath density
read, in Fourier space,

$$ \dot{r}_0(t) = -h\mu_0 \int \frac{d^d q}{(2\pi)^d} i q_1 v_q \phi_q(t) e^{i q_1 r_0(t)} + \eta_0(t), \quad (65) $$

with $\langle \eta_0(t) \eta_0(t') \rangle = 2\mu_0 T \delta(t - t')$, and

$$ \begin{aligned} \partial_t \phi_q(t) = & -\mu q_{\perp}^2 (T + u_q) \phi_q(t) - \mu q_1^2 \int \frac{d^{d-1}p_{\perp}}{(2\pi)^{d-1}} (T + u_{q_1, -\mathbf{p}_{\perp}}) \phi_{q_1, \mathbf{p}_{\perp}}(t) + \nu_q(t) \\ & - h\mu [q_{\perp}^2 v_q + q_1^2 v_{q_1}(\mathbf{x}_{\perp} = \mathbf{0})] e^{-i q_1 r_0(t)}, \end{aligned} \quad (66) $$

with the correlations of $v_q(t)$ given in equation (29). Above we denoted for brevity $h \equiv \rho_0^{-1/2}$.

With this in hands, in the next two sections we tackle two distinct (but related) problems using perturbation theory for small $h$:

(i) In section 5.1 we analyze the generalized correlation profiles between the tracer position and the density of surrounding bath particles [42]. As explained below, these quantities encode the response of the bath to the presence of a moving inclusion.

(ii) In section 5.2 we use such correlation profiles to access the effective tracer particle dynamics, and in particular its effective diffusion coefficient [35].

Note that, being constrained to move only along the backbone, we expect that the tracer particle will always exhibit diffusive behavior in spite of the comb geometry (whereas the other bath particles sub-diffuse along the backbone, as recalled in section 2).

## 5.1. Generalized correlation profiles

The simplest measure of the cross correlations between the tracer position and the density of surrounding bath particles is the average density profile in the reference frame of the tracer:

$$ \psi(\mathbf{x}, t) \equiv \langle \phi(\mathbf{x} + \mathbf{r}_0(t), t) \rangle \quad \rightarrow \quad \psi_q(t) = \langle \phi_q(t) e^{i\mathbf{q} \cdot \mathbf{r}_0(t)} \rangle, \quad (67) $$

where $\mathbf{r}_0(t) = r_0(t) \hat{e}_1$. This quantity has been first analyzed in [35] for the case of interacting Brownian particles evolving in uniform space. In particular, in the presence of an external bias applied to the tracer (which drives the system out of equilibrium), the average stationary density $\psi(\mathbf{x})$ has been proven in [35] to exhibit an algebraic tail in the direction opposite to that of the bias. This long-range effect is in fact a signature of the many-body interacting nature of the system. In the absence of a bias (i.e. at equilibrium), $\psi(\mathbf{x})$ was instead shown to be a rapidly decaying function.

For the case of particles interacting on a comb (with the tagged tracer constrained along the backbone), the average profile $\psi(\mathbf{x}, t)$ can be evaluated by using Stratonovich calculus and the Novikov theorem, as delineated in appendix C. In the stationary state attained by the system at long times, the result reads in Fourier space

$$ \psi_q = -h \frac{v_q}{T + u_q} + \mathcal{O}(h^2). \quad (68) $$

Crucially, this prediction is exactly the same as that derived in [35] for particles evolving in homogeneous $d$-dimensional space, i.e. without the comb constraint [42]. Although this may
sound surprising, in fact a non-perturbative calculation<fnref n="7" /> reveals that equation (68) follows directly from the equilibrium distribution in equation (14), showing once again that the latter correctly represents the stationary state of the system.

However, signatures of the comb structure are expected to show up by considering higher-order correlation profiles such as

$$g(\mathbf{x}, t) \equiv \langle r_0(t) \phi(\mathbf{x} + \mathbf{r}_0(t), t) \rangle \quad \rightarrow \quad g_{\mathbf{q}}(t) = \langle r_0(t) \phi_{\mathbf{q}}(t) e^{i\mathbf{q} \cdot \mathbf{r}_0(t)} \rangle. \qquad (69)$$

The information encoded in $g(\mathbf{x}, t)$ is the following: how does the bath density at distance $\mathbf{x}$ from the tracer fluctuate, in response to a fluctuation of the tracer position $\mathbf{r}_0(t)$? Similarly, all correlations profiles of the form $\langle [r_0(t)]^n \phi(\mathbf{x} + \mathbf{r}_0(t), t) \rangle$ (for integer $n$) are compactly encoded in the generating function [53]

$$w(\mathbf{x}, \lambda, t) \equiv \frac{\langle \phi(\mathbf{x} + \mathbf{r}_0(t), t) e^{\lambda r_0(t)} \rangle}{\langle e^{\lambda r_0(t)} \rangle} \quad \rightarrow \quad w_{\mathbf{q}}(\lambda, t) = \frac{\langle \phi_{\mathbf{q}}(t) e^{(\lambda+i\mathbf{q}_1)r_0(t)} \rangle}{\langle e^{\lambda r_0(t)} \rangle}, \qquad (70)$$

from which they can be retrieved by taking derivatives with respect to $\lambda$, before setting $\lambda = 0$.

An analysis of these profiles for particles interacting in homogeneous space has been recently presented in [42], with a special focus on the stationary profile $g(\mathbf{x})$. The main finding was that $g(\mathbf{x})$ exhibits a nontrivial large-distance behavior also in the absence of a bias: in fact, it displays a long-range power-law decay $g(x) \sim x^{1-d}$, whose exponent does not depend on the details of the interaction potential $u(\mathbf{x})$, but actually only on the spatial dimensionality $d > 1$ of the system. This surprising result has been checked to hold for a variety of systems (including hard-core lattice gases and Lennard-Jones suspensions), using both analytical calculations and numerical simulations.

On the comb, these generalized correlation profiles can again be evaluated by using Stratonovich calculus and the Novikov theorem, as shown in appendix C. In particular, their generating function (70) reads in the stationary limit<fnref n="8" />

$$w_{\mathbf{q}}(\lambda) = hG_0(\mathbf{q}, \lambda) \left\{ \frac{v_{\mathbf{q}}}{T + u_{\mathbf{q}}} \left[ i\mu_0 T q_1 (\lambda + iq_1) - \mu q_{\perp}^2 (T + u_{\mathbf{q}}) \right] - \mu q_1^2 \frac{v_{q_1} (\mathbf{x}_{\perp} = \mathbf{0}) + M^{(4)}(q_1, \lambda)}{1 + \mu q_1^2 M^{(3)}(q_1, \lambda)} \right\} + \mathcal{O}(h^2), \qquad (71)$$

where we introduced

$$G_0^{-1}(\mathbf{q}, \lambda) \equiv \mu q_{\perp}^2 (T + u_{\mathbf{q}}) - i\mu_0 T q_1 (2\lambda + iq_1), \qquad (72)$$

$$M^{(3)}(q_1, \lambda) \equiv \int \frac{d^{d-1}q_{\perp}}{(2\pi)^{d-1}} (T + u_{-\mathbf{q}}) G_0(\mathbf{q}, \lambda), \qquad (73)$$

$$M^{(4)}(q_1, \lambda) \equiv \int \frac{d^{d-1}q_{\perp}}{(2\pi)^{d-1}} v_{\mathbf{q}} G_0(\mathbf{q}, \lambda) [i\mu_0 q_1 T(\lambda + iq_1) - \mu q_{\perp}^2 (T + u_{\mathbf{q}})]. \qquad (74)$$

Similarly, the correlation profile $g_{\mathbf{q}}(t)$ introduced in equation (69) turns out to reduce in the stationary state to (see appendix C)

$$g_{\mathbf{q}} = h \frac{i\mu_0 T q_1}{\mu q_{\perp}^2 (T + u_{\mathbf{q}}) + \mu_0 T q_1^2} \left[ \frac{q_1^2 M^{(2)}(q_1)}{1 + q_1^2 M^{(1)}(q_1)} - \frac{v_{\mathbf{q}}}{T + u_{\mathbf{q}}} \right] + \mathcal{O}(h^2), \qquad (75)$$

<footnote n="7">⁷ See section II.D.4 in the supplemental material of [42].</footnote>

<footnote n="8">⁸ Its counterpart in homogeneous space is reported for comparison in equation (C.14). Similarly, the profile $g_{\mathbf{q}}$ in equation (75) should be compared with equation (12) in [42], valid in homogeneous space.</footnote>
where we called

$$M^{(1)}(q_1) = \int \frac{d^{d-1}q_{\perp}}{(2\pi)^{d-1}} \frac{\mu(T+u_{-\mathbf{q}})}{\mu q_{\perp}^2(T+u_{\mathbf{q}}) + \mu_0 T q_1^2}, \qquad (76)$$

$$M^{(2)}(q_1) = \int \frac{d^{d-1}q_{\perp}}{(2\pi)^{d-1}} \frac{\mu v_{\mathbf{q}}}{\mu q_{\perp}^2(T+u_{\mathbf{q}}) + \mu_0 T q_1^2}. \qquad (77)$$

By symmetry in the momenta $\mathbf{q}$, we deduce that $g(-\mathbf{x}) = -g(\mathbf{x})$, as expected from its definition (69)—see also the inset of figure 3(a). Moreover, if $\mu_0 = \mu$, then the stationary correlation profile $g_{\mathbf{q}}$ becomes independent of the particles mobility. In $d=2$ and setting $v_{\mathbf{q}} = u_{\mathbf{q}}$, it is then simple to derive its large-distance behavior (see appendix C)

$$g(x_1) \sim \frac{h\sqrt{T}u_{\mathbf{q}=0}}{2\pi (T+u_{\mathbf{q}=0})^{3/2}x_1^{-1}}. \qquad (78)$$

This asymptotic behavior is checked in figure 3(a). Remarkably, it is consistent with the behavior $g(x_1) \sim x_1^{1-d}$ observed in uniform $d$-dimensional space (i.e. without the comb) in [42]. In this sense, equation (78) further validates the claim of universality of the large-distance behavior of the correlation profile $g(\mathbf{x})$, whose decay exponent was found in [42] to depend only on the spatial dimensionality $d$ of the system.

## 5.2. Effective diffusion coefficient

In this section we focus on the effective dynamics of the tracer particle, and we derive the first perturbative correction to its diffusion coefficient, due to the presence of the other bath particles. First, note that the leading-order correction to the generalized profiles computed in section 5.1 turned out to be of $O(h)$. Conversely, any correlation function involving $\mathbf{r}$ but not $\phi$ must exhibit corrections at least of $O(h^2)$—the simplest way to prove this fact is to note that the system of equations (65) and (66) is invariant under the transformation $\{h \mapsto -h, \phi \mapsto -\phi\}$ [46, 54]. In general, the perturbative calculation of the effective diffusion coefficient in homogeneous space can be compactly addressed within the path-integral formalism developed in [45] (and later adopted in e.g. [35, 55]). However, it is not evident how to extend such path-integral techniques to the case of non-homogeneous space analyzed here. By contrast, we showed in [42] that the tracer-bath correlation profiles—that we already derived in the previous section—actually encode the tracer statistics, and thus offer a straightforward alternative way to access its diffusion coefficient.

To show this, let us focus on the moment generating function

$$\varphi(q_1, t) \equiv \langle e^{-iq_1 r_0(t)} \rangle, \qquad (79)$$

that we aim to compute in the stationary regime. To this end, one can first use Stratonovich calculus to derive, starting from equations (65) and (66), the relation<fnref n="9" />

$$\partial_t \langle e^{\lambda r_0(t)} \rangle = \lambda \langle \eta(t) e^{\lambda r_0(t)} \rangle - h\mu_0\lambda \cdot \int \frac{d^d p}{(2\pi)^d} i p_1 v_{\mathbf{p}} \langle e^{(\lambda+ip_1)r_0(t)} \phi_{\mathbf{p}}(t) \rangle. \qquad (80)$$

<footnote n="9">⁹ See section II.G in the supplemental material of [42].</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** (a) Tracer–bath spatial correlation profile $g(x_1)$, whose Fourier transform is predicted by equation (75) in the long-time stationary state. Its asymptotic scale-free behavior is given in equation (78), while the inset shows its symmetry $g(-x) = -g(x)$. In the plot we chose a Gaussian interaction potential $u_q = \exp(-q^2/2)$, and we set $h$ and $T$ to unity for simplicity. (b) Comparison between the effective diffusion coefficient $D_{eff}$ predicted by equation (89), and the results of numerical simulations of a system of particles interacting via Gaussian potentials of strength $\epsilon$. In the plot, we vary $\epsilon$ (and hence the dimensionless parameter $\theta_1 \equiv \epsilon\rho_0\sigma^d/T$, see section 5.2), with all other parameters kept fixed. Details of the numerical simulation are reported in appendix D.</figcaption></figure>

The second expectation value on the r.h.s. of equation (80) is reminiscent of a quantity we have already met, i.e. the generating function $w_p(\lambda, t)$ of tracer–bath correlation profiles introduced in equation (70). Since the evolution equation of $w_p(\lambda, t)$ (derived perturbatively in appendix C.1) admits as a stationary solution $w_p(\lambda)$ given in equation (70), here it is reasonable to assume that

$$ \langle e^{(\lambda+ip_1)r_0(t)} \phi_p(t) \rangle = w_p(\lambda, t) \langle e^{\lambda r_0(t)} \rangle \approx w_p(\lambda) \langle e^{\lambda r_0(t)} \rangle, \quad (81) $$

in the stationary regime. Plugging this into equation (80) then gives

$$ \partial_t \langle e^{\lambda r_0(t)} \rangle = \lambda \langle \eta(t) e^{\lambda r_0(t)} \rangle - h\mu_0 \langle e^{\lambda r_0(t)} \rangle \lambda_1 \int \frac{d^d p}{(2\pi)^d} i p_1 v_p w_p(\lambda). \quad (82) $$

On the other hand, the first expectation value on the r.h.s. of equations (80) and (82) can be computed by standard methods using the Novikov theorem [56], and the result reads

$$ \langle \eta(t) e^{\lambda r_0(t)} \rangle = \mu_0 T \lambda \langle e^{\lambda r_0(t)} \rangle, \quad (83) $$

which is the same as for a non-interacting tracer (in spite of the interaction with the bath density, see [42] for further details). By setting $\lambda = -iq_1$ and using equation (79), we thus obtain

$$ \partial_t \varphi(q_1, t) = -q_1^2 \mu_0 T \varphi(q_1, t) - h \mu_0 \varphi(q_1, t) \int \frac{d^d p}{(2\pi)^d} (p_1 \cdot q_1) v_p w_p(-iq_1), \quad (84) $$
whose solution<fnref n="10" />

$$ \ln \varphi(q_1, t) = -\mu_0 t \left[ q_1^2 T + h \int \frac{d^d p}{(2\pi)^d} (p_1 \cdot q_1) v_p w_p (-iq_1) \right] \quad (85) $$

gives the leading-order correction to the cumulant generating function of the tracer position. Thus, as anticipated above, equation (85) links the statistical properties of the tracer position $r_0(t)$ to the knowledge of the tracer-bath correlation profiles encoded in $w_p(\lambda)$. In fact, this expression can be made explicit by using the leading-order estimate of $w_p(\lambda)$ found in equation (71). Several comments are now in order. First, we note that the tracer position $r_0(t)$ has non-Gaussian statistics (indeed, the r.h.s. of equation (85) is not simply quadratic in $q_1$). In particular, one can check that $\varphi(q_1 = 0, t) = 1$ and that the average tracer position vanishes,

$$ \langle r_0(t) \rangle = i \frac{\partial \varphi(q_1, t)}{\partial q_1} \bigg|_{q_1=0} = 0, \quad (86) $$

as expected by the symmetry of the problem under $x_1 \mapsto -x_1$ (indeed, we have assumed that the tracer starts its motion at the origin, i.e. $r_0(t=0) = 0$). Conversely, its variance reads

$$ \begin{aligned} \langle r_0^2(t) \rangle &= -\frac{\partial^2}{\partial q_1^2} \varphi(q_1, t) \bigg|_{q_1=0} \\ &= 2\mu_0 T t \left\{ 1 - h^2 \int \frac{d p_1}{2\pi} p_1^2 \left[ M^{(5)}(p_1) - \frac{\mu_0}{\mu} \frac{[p_1 M^{(2)}(p_1)]^2}{1 + p_1^2 M^{(1)}(p_1)} \right] \right\} + \mathcal{O}(h^4), \end{aligned} \quad (87) $$

where

$$ M^{(5)}(p_1) \equiv \int \frac{d^{d-1}p_{\perp}}{(2\pi)^{d-1}} \frac{\mu_0 v_{\perp}^2}{(T+u_{\perp}) [\mu q_{\perp}^2 (T+u_{\perp}) + \mu_0 T q_1^2]}, \quad (88) $$

while the functions $M^{(1,2)}(p_1)$ were given in equations (76) and (77), respectively. Calling $D_0 = \mu T$ the bare diffusion coefficient of the non-interacting tracer, from equation (87) we thus obtain $\langle r_0^2(t) \rangle = 2D_{\text{eff}}t$ with the effective diffusion coefficient

$$ \frac{D_{\text{eff}}}{D_0} = 1 - h^2 \int \frac{d p_1}{2\pi} p_1^2 \left[ M^{(5)}(p_1) - \frac{\mu_0}{\mu} \frac{[p_1 M^{(2)}(p_1)]^2}{1 + p_1^2 M^{(1)}(p_1)} \right] + \mathcal{O}(h^4). \quad (89) $$

This prediction represents another key result of this work, and is the spatially heterogeneous counterpart to the diffusion coefficient first obtained in [35] for homogeneous media, using path-integral methods (see equation (52) therein).

Again, we note that for $\mu_0 = \mu$ the second term in equation (89) actually becomes $\mu$-independent. To be concrete, we now set $v_q = u_q$, and assume for the interaction potential in real space $U(\mathbf{x})$ the general rotationally invariant form

$$ U(\mathbf{x}) = \epsilon f(x/\sigma), \quad (90) $$

<footnote n="10">¹⁰ Note that equation (85) apparently implies that not only the variance, but actually all other higher-order cumulants of $r_0(t)$ also grow linearly with $t$. However, we stress that equation (85) has been derived under the assumption that $w_q(\lambda, t)$ becomes stationary, and must thus be taken with caution—we comment further on this point in appendix C.3.</footnote>
where $x = |\mathbf{x}|$, $\epsilon$ is an energy scale, $\sigma$ is a microscopic length scale, and the function $f$ is dimensionless. For example, within the Gaussian core model [57] one has $U(x) = \epsilon \exp[-x^2/(2\sigma^2)]$. In Fourier space, the corresponding rescaled interaction potential $u_q = \rho_0 U_q$ then reads

$$u_q = \rho_0 \epsilon \sigma^d f(q\sigma). \qquad (91)$$

By inspecting equation (89) with $\mu_0 = \mu$, and upon reinstating the microscopic length scale $a$ by replacing $q_1^2 \mapsto a^{d-1}q_1^2$ (see the discussion in footnote 4), it is then simple to see that $D_{\text{eff}}/D_0$ only depends on the dimensionless combinations $\theta_1 \equiv \epsilon\rho_0\sigma^d/T$, $\theta_2 = \rho_0\sigma^d$, and $\theta_3 = R/a$, rather than on these parameters taken separately. In figure 3(b) we thus compare our prediction in equation (89) to the results of numerical simulations obtained by varying the interaction strength $\epsilon$ with all other parameters fixed (see appendix D for the details), which amounts to varying the effective parameter $\theta_1$. The prediction is shown to be accurate for small interaction strengths $\epsilon \ll k_B T$, as usually expected within the linearized Dean–Kawasaki theory [35, 42]—this point is further discussed in the next section.

### 5.3. Validity of the linearized theory

We conclude our discussion by briefly recalling the assumptions under which our calculation is expected to accurately reproduce the results of numerical simulations. These assumptions are not specific to the comb geometry, but actually underlie similar treatments (even in uniform space) based on the linearized Dean–Kawasaki theory [35]: (i) the system should admit a well-defined average background density $\rho_0$ (see equation (15)); note that, for very dilute systems with attractive potentials, clustering phenomena may take place and undermine this assumption; (ii) the interaction potentials should admit a Fourier transform, which excludes strong short-distance repulsion (such as the one provided by Lennard–Jones-type potentials); (iii) the system should not be too dilute, our results being expressed in a power series in the small parameter $h = \rho_0^{-1/2}$. More precisely, the latter parameter not being dimensionless, an estimate of the limits of validity of the perturbative expansion can be rather formulated as [42]

$$\frac{\epsilon}{k_B T} \ll \sigma^{-d} \left[ \rho_0 \int \frac{d^d q}{(2\pi)^d} f^2(q\sigma) \right]^{-1/2}, \qquad (92)$$

depending on the specific details of the interaction potential (90). This effectively conveys that the linearized Dean–Kawasaki theory is an expansion around the weakly interacting case, $\epsilon$ being the typical strength of these interactions.

## 6. Conclusions

In this work we considered a system of overdamped Brownian particles interacting via soft pairwise potentials, and constrained to move on a comb-like structure (see figure 1). Within the Dean–Kawasaki formalism, we first derived in section 3 the exact equations that describe the dynamics of the particle density field $\rho(\mathbf{x}, t)$. The latter can be expanded around a uniform background density, which allowed us to characterize the Gaussian fluctuations of the density field. In particular, in section 4 we computed the density–density correlator both starting from fixed initial conditions (see equation (44)), and in the stationary limit assumed by the system at long times (see equation (58)). The latter can be used to construct for instance the ISF shown in figure 2. To the best of our knowledge, this represents the first application of the Dean–Kawasaki theory to non-homogeneous media [33, 34].
Furthermore, in section 5 we singled out a tagged tracer from the bath of interacting particles, under the assumption that it is constrained to move only along the backbone. In this setting, we derived the spatial correlation profile between the tracer position and the surrounding bath density (see equation (75) and figure 3(a)), and the effective diffusion coefficient of the tracer, which we tested using Brownian dynamics simulations (see equation (89) and figure 3(b)).

Further extensions of this work could address the unconstrained problem in which the tracer particle is allowed to leave the backbone (which is technically more challenging than the constrained case considered in this work, even within the perturbation theory delineated in section 5). This would ideally provide access to the effective subdiffusion coefficient of the tracer along the backbone, i.e. to its correction due to the presence of the other bath particles. Moreover, the framework delineated in section 3 can in principle be extended to other non-homogeneous geometries (such as fractal networks and random percolation structures), which can be encoded with a suitable choice of the particle's mobility matrix $\hat{\mu}(x)$ in equation (1).

## Data availability statement

All data that support the findings of this study are included within the article (and any supplementary files).

## Acknowledgments

DV thanks Leticia F Cugliandolo, Pietro Luigi Muzzeddu, and Gilles Tarjus for illuminating discussions.

## Appendix A. Details on the dynamics of the density fluctuations

In this appendix we completely characterize the density fluctuations on the comb by deriving the results reported in section 4. In particular, in appendices A.1–A.3 we provide the derivations of the propagator, response function, and Dirichlet correlator, respectively. In appendix A.4 we then derive the stationary correlator, proving in appendix A.5 that it coincides with the equilibrium correlator (in spite of the geometrical constraints of the system, which could in principle preclude ergodicity). In appendix A.6 we recall a few causality properties of the double Laplace transform, and finally in appendix A.7 we analyze in the Laplace domain the free (Gaussian) field [52], as a test bench for the techniques employed on the comb.

### A.1. Derivation of the propagator

We start from equation (31), which rules the evolution of the Fourier transform $G_q(t)$ of the propagator $G(x, t)$. We then Laplace-transform the last equation: using that $G(x, t = 0^+) = 0$ by construction, we find

$$G_q(s) = g_{\perp}(q_{\perp}, s) [1 - \mu q_1^2 F(q_1, s)], \qquad (A.1)$$

where $g_{\perp}$ was defined in equation (33), while we introduced the auxiliary function

$$F(q_1, s) = \int \frac{d^{d-1}q_{\perp}}{(2\pi)^{d-1}} (T + u_{q_1, -q_{\perp}}) G_q(s). \qquad (A.2)$$
Crucially, the terms enclosed in square brackets in equation (A.1) do not depend on $\mathbf{q}_{\perp}$. Upon applying $\int \frac{d^{d-1}\mathbf{q}_{\perp}}{(2\pi)^{d-1}}(T+u_{q_1,-\mathbf{q}_{\perp}})$ to equation (A.1), we thus obtain the self-consistent equation

$$F(q_1, s) = M(q_1, s) [1 - \mu q_1^2 F(q_1, s)], \qquad (A.3)$$

where we introduced the function

$$M(q_1, s) = \int \frac{d^{d-1}q_{\perp}}{(2\pi)^{d-1}} (T + u_{q_1, -q_{\perp}}) g_{\perp}(\mathbf{q}_{\perp}, s). \qquad (A.4)$$

We thus solve equation (A.3) to obtain $F$, which can then be replaced into equation (A.1) to give

$$G_{\mathbf{q}}(s) = \frac{g_{\perp}(\mathbf{q}_{\perp}, s)}{1 + \mu q_1^2 M(q_1, s)}. \qquad (A.5)$$

Upon calling (compare with equation (34))

$$g_1^{-1}(q_1, s) \equiv 1 + \mu q_1^2 M(q_1, s) \qquad (A.6)$$

the denominator in equation (A.5), we obtain the expression of $G_q(s)$ reported in equation (32).

Note that equation (31) reduces, for $u_q = 0$ and in $d=2$, to the Fourier transform of the Fokker-Planck equation (7) for a single (non-interacting) Brownian particle on the comb, placed initially at the origin at time $t=0$. The corresponding solution $G_q(s)$, i.e. its moment generating function written in the Laplace domain, then reduces to

$$G_{\mathbf{q}}(s) = (s + \mu T q_2^2)^{-1} \left[ 1 + \frac{a \sqrt{\mu T q_1^2}}{2\sqrt{s}} \right]^{-1}, \qquad (A.7)$$

where we computed explicitly the integral in the definition (34) of $g_1(q_1, s)$, and we reinstated the length scale $q_1^2 \mapsto aq_1^2$ (see the discussion in footnote 4). In particular, one can use the Tauberian theorems [58] to deduce the long-$t$ behavior of the inverse Laplace transform of the particle variance along the two orthogonal directions:

$$-\left. \frac{\partial^2 G_q(s)}{\partial q_2^2} \right|_{q=0} = 2\mu Ts^{-2} \quad \leftrightarrow \quad \langle y^2(t) \rangle \approx 2\mu Tt, \qquad (A.8)$$

and

$$-\left. \frac{\partial^2 G_q(s)}{\partial q_1^2} \right|_{q=0} = a \sqrt{\mu T} s^{-\frac{3}{2}} \quad \leftrightarrow \quad \langle x^2(t) \rangle \approx 2a \sqrt{\frac{\mu Tt}{\pi}}. \qquad (A.9)$$

Alternatively, one can invert explicitly equation (A.7) to the time domain, finding

$$G_q(t) = \frac{4q_2^2 e^{-q_2^2 \mu Tt} - 4aq_1^2 \sqrt{\frac{q_2^2}{\pi}} \text{DawsonF}\left(\sqrt{q_2^2 \mu Tt}\right) + a^2 e^{\frac{q_1^4}{4} a^2 \mu Tt} q_1^4 \text{Erfc}\left(\frac{q_1^2}{2} a \sqrt{\mu Tt}\right)}{a^2 q_1^4 + 4q_2^2}, \qquad (A.10)$$

which leads to the same asymptotics in the long-$t$ limit.
## A.2. Derivation of the response function

We start by Laplace transforming the corresponding Schwinger–Dyson equation (40) with respect to both $t \to s$ and $t' \to s'$, finding

$$s R_{q,p}(s, s') = -\mu q_{\perp}^2 (T + u_q) R_{q,p}(s, s') - \mu q_1^2 \int \frac{d^{d-1}k_{\perp}}{(2\pi)^{d-1}} (T + u_{q_1, -k_{\perp}}) R_{(q_1, k_{\perp}), p}(s, s') \\ + \frac{2\mu T}{s+s'} [q_1^2 \delta(q_1+p_1) + q_{\perp}^2 \delta^d(q+p)] . \qquad (A.11)$$

Note that the first Laplace transform, $t \to s$, acts on the first term in equation (40) as

$$\partial_t R_{q,p}(t, t') \to s R_{q,p}(s, t') - R_{q,p}(t=0^+, t') \Theta(-t'), \qquad (A.12)$$

where the $\Theta(-t')$ stems from the causality of the response function—see its definition in equation (38). The second Laplace transform, $t' \to s'$, then suppresses the term proportional to $\Theta(-t')$.

Equation (A.11) can be rewritten as

$$R_{q,p}(s, s') = g_{\perp}(q, s) \left\{ \frac{2\mu T}{s+s'} [q_1^2 \delta(q_1+p_1) + q_{\perp}^2 \delta^d(q+p)] - \mu q_1^2 F^{(1)}(q_1; s, s') \right\}, \qquad (A.13)$$

where we introduced the auxiliary function

$$F^{(1)}(q_1; s, s') \equiv \int \frac{d^{d-1}q_{\perp}}{(2\pi)^{d-1}} (T + u_{q_1, -q_{\perp}}) R_{q,p}(s, s'). \qquad (A.14)$$

Applying the operator $\int \frac{d^{d-1}q_{\perp}}{(2\pi)^{d-1}} (T + u_{q_1, -q_{\perp}})$ to equation (A.13) we thus obtain the self-consistency relation

$$F^{(1)}(q_1; s, s') = M(q_1, s) \left[ \frac{2\mu T}{s+s'} q_{\perp}^2 \delta^d(q+p) - \mu q_1^2 F^{(1)}(q_1; s, s') \right] \\ + \frac{2\mu T}{s+s'} q_1^2 \delta(q_1+p_1) (T+u_p) g_{\perp}(-p, s). \qquad (A.15)$$

From now on, we will assume for simplicity $u(x)$ to be rotationally invariant, so that $u_q$ depends only on $|q| = q_1^2 + q_{\perp}^2$. This implies in particular

$$g_{\perp}(q, s) = g_{\perp}(-q, s), \qquad M(q_1, s) = M(-q_1, s), \qquad (A.16)$$

as can easily be evinced from their definitions in equations (33) and (A.4). Solving for $F^{(1)}$ in equation (A.15) and replacing the result into equation (A.13) then renders the response function given in equation (41), upon noting that

$$1 - \mu p_{\perp}^2 (T + u_p) g_{\perp}(p, s) = s g_{\perp}(p, s), \qquad (A.17)$$

and by using the propagator $G_q(s)$ introduced in equation (32).
Using the causality structure of the response function (see cf appendix A.6), we can finally invert equation (41) with $u(\mathbf{x}) = A\delta(\mathbf{x})$ and $d = 2$, which gives

$$R_{q,p}(t, t') = \frac{2\mu T \Theta(\tau) q_1^2 \delta(q_1 + p_1)}{(q_2^2 - p_2^2)(4p_2^2 + p_1^4)(4q_2^2 + q_1^4)} \left\{ \begin{aligned} & 2q_2^3 (4p_2^2 + p_1^4) e^{-q_2^2 \mu T \tau} \left[ 2q_2 \right. \\ & \left. -q_1^2 \text{Erfi}(q_2 \sqrt{\mu T \tau}) \right] - (\mathbf{p} \leftrightarrow \mathbf{q}) + e^{\frac{q_1^4}{4} \mu T \tau} (q_2^2 - p_2^2) q_1^8 \text{Erfc}\left(\frac{q_1^2}{2} \sqrt{\mu T \tau}\right) \end{aligned} \right\} \\ & + 2\mu T \Theta(\tau) p_2^2 \delta(\mathbf{q} + \mathbf{p}) e^{-q_2^2 \mu T \tau}, \quad (A.18)$$

with $\tau = t - t' > 0$ and $T \equiv T + A$.

### A.3. Derivation of the correlator for Dirichlet initial conditions

As in the previous sections, we start by Laplace-transforming the corresponding Schwinger-Dyson equation (39) with respect to $t \to s$, finding

$$C_{q,p}(s, t') = g_{\perp}(\mathbf{q}, s) \left[ C_{q,p}(t=0^{+}, t') - \mu q_1^2 \int \frac{d^{d-1}k_{\perp}}{(2\pi)^{d-1}} (T+u_{q_1, -\mathbf{k}_{\perp}}) C_{(q_1, \mathbf{k}_{\perp}), \mathbf{p}}(s, t') \right. \\ \left. + R_{p,q}(t', s) \right]. \quad (A.19)$$

Note that this time the quantity $C_{q,p}(t=0^{+}, t') \neq 0$. On the contrary, this term can be used to enforce any particular initial condition at $t_0 = 0$, because by construction

$$C_{q,p}(t = t_0, t') = \phi_q(t_0) \phi_p(t_0) G_p(t' - t_0), \quad (A.20)$$

where $G_q(t)$ is the propagator in equation (32). We remark that the seemingly arbitrary choice $t_0 = 0$ is necessary to allow the use of the Laplace transform, but in fact it entails no loss of generality. Indeed, once a particular solution $C_{q,p}(t, t')$ is found, one can simply introduce the new coordinates

$$\tau = t + t_0, \quad \tau' = t' + t_0, \quad (A.21)$$

which practically amounts to shifting the initial conditions to the time $\tau = t_0$. In particular, this may allow to find the correlator in the stationary state by eventually taking the limit $t_0 \to -\infty$.
Upon further Laplace-transforming equation (A.19) with respect to $t' \to s'$, we thus obtain

$$ C_{q,p}(s,s') = g_{\perp}(\mathbf{q},s) \left[ \phi_{\mathbf{q}}^{(0)} \phi_{\mathbf{p}}^{(0)} G_{\mathbf{p}}(s') - \mu q_1^2 \int \frac{d^{d-1}k_{\perp}}{(2\pi)^{d-1}} (T+u_{q_1,-\mathbf{k}_{\perp}}) C_{(q_1,\mathbf{k}_{\perp}),\mathbf{p}}(s,s') \right. \\ \left. + R_{\mathbf{p},\mathbf{q}}(s',s) \right], \qquad (A.22) $$

where we called for brevity $\phi_{\mathbf{q}}^{(0)} \equiv \phi_{\mathbf{q}}(t=0)$, and where the response function $R_{\mathbf{q},\mathbf{p}}(s,s')$ was given in equation (41) (note, however, the exchange of its arguments). Again, our strategy consists in integrating equation (A.22) with respect to $\mathbf{q}_{\perp}$ to obtain a self-consistency equation. To this end, we first collect the various term in equation (A.22) according to their dependence on $\mathbf{q}_{\perp}$, and use the property

$$ g_{\perp}(\mathbf{q},s)g_{\perp}(\mathbf{q},s') = \frac{g_{\perp}(\mathbf{q},s) - g_{\perp}(\mathbf{q},s')}{s' - s}, \qquad (A.23) $$

which is readily proven by inspection of the definition (33) of $g_{\perp}(\mathbf{q},s)$. We can thus rewrite

$$ C_{\mathbf{q},\mathbf{p}}(s,s') = g_{\perp}(\mathbf{q},s) \phi_{\mathbf{q}}^{(0)} \phi_{\mathbf{p}}^{(0)} G_{\mathbf{p}}(s') + g_{\perp}(\mathbf{q},s) g_{\perp}(\mathbf{q},s') \delta^d(\mathbf{q}+\mathbf{p}) \frac{2\mu T}{s+s'} p_{\perp}^2 \\ + g_{\perp}(\mathbf{q},s) \left[ \frac{2\mu T}{s+s'} \frac{p_1^2 \delta(q_1+p_1) s' G_{\mathbf{p}}(s')}{s'-s} - \mu q_1^2 F^{(2)}(q_1;s,s') \right] \\ - g_{\perp}(\mathbf{q},s') \frac{2\mu T}{s+s'} \frac{p_1^2 \delta(q_1+p_1) s' G_{\mathbf{p}}(s')}{s'-s}, \qquad (A.24) $$

where we introduced the auxiliary function

$$ F^{(2)}(q_1; s, s') = \int \frac{d^{d-1}q_{\perp}}{(2\pi)^{d-1}} (T + u_{q_1, -\mathbf{q}_{\perp}}) C_{\mathbf{q},\mathbf{p}}(s, s'). \qquad (A.25) $$

Applying $\int \frac{d^{d-1}q_{\perp}}{(2\pi)^{d-1}} (T+u_{q_1,-\mathbf{q}_{\perp}})$ to equation (A.22) renders

$$ F^{(2)}(q_1; s, s') = \phi_{\mathbf{p}}^{(0)} G_{\mathbf{p}}(s') M^{(0)}(q_1, s) + M(q_1, s) \left[ \mathcal{A} - \mu q_1^2 F^{(2)}(q_1; s, s') \right] \\ - \mathcal{A} M(q_1, s) + g_{\perp}(\mathbf{p}, s) g_{\perp}(\mathbf{p}, s') \delta(q_1 + p_1) \frac{2\mu T}{s+s'} p_{\perp}^2 (T+u_{\mathbf{p}}), \qquad (A.26) $$

where the functions $M^{(0)}$, $M$ were introduced in equations (45) and (A.4), respectively, and where we called for brevity

$$ \mathcal{A} \equiv \frac{2\mu T}{s+s'} \frac{p_1^2 \delta(q_1+p_1) s' G_{\mathbf{p}}(s')}{s'-s}. \qquad (A.27) $$

Solving for $F^{(2)}$ in equation (A.26) gives

$$ \frac{F^{(2)}(q_1; s, s')}{g_1(q_1, s)} = \phi_{\mathbf{p}}^{(0)} G_{\mathbf{p}}(s') M^{(0)}(q_1, s) + \mathcal{A}[M(q_1, s) - M(q_1, s')] \\ + g_{\perp}(\mathbf{p}, s) g_{\perp}(\mathbf{p}, s') \delta(q_1 + p_1) \frac{2\mu T}{s+s'} p_{\perp}^2 (T + u_{\mathbf{p}}), \qquad (A.28) $$
which can be further simplified by rewriting

$$ \mu p_{\perp}^2 (T + u_{\mathbf{p}}) = g_{\perp}^{-1}(\mathbf{p}, s) - s, \qquad (A.29) $$

$$ \mu q_1^2 [M(q_1, s) - M(q_1, s')] = g_1^{-1}(q_1, s) - g_1^{-1}(q_1, s'), \qquad (A.30) $$

as one can easily check starting from the definitions of $g_{\perp}$, $g_1$, and $M$ given in equations (33), (34) and (A.4), respectively. Replacing the expression of $F^{(2)}$ found in equation (A.28) back into equation (A.24) renders, after some algebra, the correlation function $C_{q,p}(s, s')$ reported in equation (44).

Note that specializing the latter to flat initial conditions $\phi_q^{(0)} \mapsto \phi^{(0)}$ yields

$$ \phi_{\mathbf{p}}^{(0)} G_{\mathbf{p}}(s') g_{\perp}(\mathbf{q}, s) [\phi_{\mathbf{q}}^{(0)} - \mu q_1^2 g_1(q_1, s) M^{(0)}(q_1, s)] \\ \mapsto [\phi^{(0)}]^2 G_{\mathbf{p}}(s') g_{\perp}(\mathbf{q}, s) [1 - \mu q_1^2 g_1(q_1, s) M(q_1, s)] = [\phi^{(0)}]^2 G_{\mathbf{q}}(s) G_{\mathbf{p}}(s'), \qquad (A.31) $$

where the second step follows from the definitions of $M^{(0)}$, $M$ given in equations (45) and (A.4), while in the third step we recognized

$$ 1 - \mu q_1^2 g_1(q_1, s) M(q_1, s) = g_1(q_1, s), \qquad (A.32) $$

and we used the definition (32) of $G_q(s)$.

### A.4. Derivation of the equilibrium correlator

First, we note that simpler field theories generally allow to find the stationary correlator starting from the Dirichlet correlator computed in correspondence of fixed initial conditions at time $t = t_0$. Indeed, if a time-domain expression of $C_{q,p}(t, t')$ is available, one can formally send $t_0 \to -\infty$, as we remarked at the beginning of appendix A.3. We exemplify this route in appendix A.7 for the simple case of the free (Gaussian) field [52]. However, directly inverting the Laplace expression of $C_{q,p}(s, s')$ found in equation (44) proved challenging in the present context.

The strategy we adopt in this section hinges instead on the assumption that the system reaches at long times the equilibrium distribution given in equation (56). Here the fluctuation-dissipation theorem (in the form of equation (54)) links the stationary correlator $C_{q,p}(\tau)$ to the linear susceptibility

$$ \chi_{q,p}(t, t') = \frac{\delta \langle \phi_q(t) \rangle}{\delta j_{-p}(t')} \Big|_{j=0}, \qquad (A.33) $$

computed in the presence of a field $j(\mathbf{x}, t)$ linearly coupled to the Hamiltonian as in equation (50). The corresponding equation of motion follows from equation (25) as

$$ \partial_t \phi(\mathbf{x}, t) = \nabla \cdot \left\{ \hat{\mu}(\mathbf{x}) \nabla \left[ T\phi(\mathbf{x}, t) + \int dy u(\mathbf{x} - \mathbf{y}) \phi(\mathbf{y}) - j(\mathbf{x}, t) \right] + \xi(\mathbf{x}, t) \right\}. \qquad (A.34) $$

Taking its average, its spatial Fourier transform, and finally its functional derivative as in equation (A.33), leads to the following evolution equation for the susceptibility:
$$ \partial_t \chi_{\mathbf{q},\mathbf{p}}(t,t') = -\mu q_{\perp}^2 (T+u_{\mathbf{q}}) \chi_{\mathbf{q},\mathbf{p}}(t,t') - \mu q_1^2 \int \frac{d^{d-1}k_{\perp}}{(2\pi)^{d-1}} (T+u_{q_1,-\mathbf{k}_{\perp}}) \chi_{(q_1,\mathbf{k}_{\perp}),\mathbf{p}}(t,t') \\ + \mu \delta(t-t') [q_1^2 \delta^d(\mathbf{q}+\mathbf{p}) + q_{\perp}^2 \delta^d(\mathbf{q}+\mathbf{p})]. \quad (A.35) $$

The latter is formally analogous (up to a multiplicative factor $2T$) to the evolution equation (40) satisfied by the response function $R_{\mathbf{q},\mathbf{p}}(t,t')$, whose solution has been computed in appendix A.2, hence we can write immediately its solution

$$ \chi_{\mathbf{q},\mathbf{p}}(s) = \frac{\mu}{s+s'} g_{\perp}(p_{\perp},s) [q_{\perp}^2 \delta^d(\mathbf{q}+\mathbf{p}) + q_1^2 \delta^d(q_1+p_1) s G_{\mathbf{q}}(s)] = \frac{\tilde{\chi}_{\mathbf{q},\mathbf{p}}(s)}{s+s'}. \quad (A.36) $$

Finally, by integrating the fluctuation-dissipation relation (54) one can obtain the stationary correlator $C(s)$, provided that one first can fix the initial condition

$$ C_{\mathbf{q},\mathbf{p}}(\tau=0^+) = \langle \phi_{\mathbf{q}}(t) \phi_{\mathbf{p}}(t) \rangle_{t \to \infty}. \quad (A.37) $$

The latter encodes the equal-time fluctuations of the field in the equilibrium state described by the distribution $P[\phi]$ in equation (56). Since $P[\phi]$ features the Hamiltonian $H[\phi]$ in equation (50), which is Gaussian, the equilibrium fluctuations of $\phi$ can be simply obtained by first constructing the generating functional [59]

$$ \mathcal{Z}[j] = \int \mathcal{D}\phi \exp\left\{-\frac{1}{T}\mathcal{H}_j[\phi]\right\} = \exp\left(\frac{1}{2T}\int \frac{d^d q}{(2\pi)^d} \frac{j_q j_{-q}}{T+u_q}\right), \quad (A.38) $$

where $\mathcal{H}_j[\phi]$ was introduced in equation (50), and where the functional integral is assumed to be normalized so that $\mathcal{Z}[j=0]=1$. We can then compute

$$ \langle \phi_{\mathbf{q}} \phi_{\mathbf{p}} \rangle = T^2 \frac{\delta^2 \mathcal{Z}[j]}{\delta j_{-\mathbf{q}} \delta j_{-\mathbf{p}}} \Big|_{j=0} = \frac{T}{T+u_{\mathbf{q}}} \delta^d(\mathbf{q}+\mathbf{p}), \quad (A.39) $$

which is the result reported in equation (57). Using the fluctuation-dissipation theorem written in the Laplace domain as in equation (55), we then find the stationary correlator given in equation (58).

The latter can eventually be inverted to the time domain for selected choices of $u_{\mathbf{q}}$. Note that the so-obtained $C_{\mathbf{q},\mathbf{p}}(\tau)$ is only expected to reproduce the branch with $\tau > 0$, while we physically expect it to be symmetric for $\tau < 0$, i.e. $C_{\mathbf{q},\mathbf{p}}(\tau) = C_{\mathbf{q},\mathbf{p}}(|\tau|)$. The result depends in general on the spatial dimensionality, hence we find for instance in $d=2$, with $u(\mathbf{x}) = A\delta(\mathbf{x})$,

$$ C_{\mathbf{q},\mathbf{p}}(\tau) = \frac{T}{\mathcal{T}} \delta(\mathbf{q}+\mathbf{p}) e^{-q_2^2 \mu \mathcal{T} \tau} - \frac{T}{\mathcal{T}} q_1^2 \delta(q_1+p_1) \left\{ \frac{2q_2 e^{-q_2^2 \mu \mathcal{T} \tau} [2q_2 - q_1^2 \text{Erfi}(q_2 \sqrt{\mu \mathcal{T} \tau})]}{(p_2^2 - q_2^2)(4q_2^2 + q_1^4)} \right. \\ \left. + (\mathbf{q} \leftrightarrow \mathbf{p}) + \frac{4q_1^4 e^{\frac{q_1^4}{4} \mu \mathcal{T} \tau} \text{Erfc}(\frac{q_1^2}{2} \sqrt{\mu \mathcal{T} \tau})}{(4p_2^2 + p_1^4)(4q_2^2 + q_1^4)} \right\}, \quad (A.40) $$

where we called again $\mathcal{T} \equiv T+A$.

As we stressed, the derivation proposed in this section relied on the assumption that the system reaches thermal equilibrium described by the Boltzmann distribution. In the next section we thus prove that the equilibrium solution represents the only admissible time-translational invariant solution for this system.
### A.5. Proof of the uniqueness of the stationary solution

In section 4.5, the sole requirement of time-translational invariance led to the condition (64) for the correlator $C_{q,p}(t, t')$. For simpler field theories, this condition is sufficient to construct explicitly the stationary correlator—we exemplify this for the free field in appendix A.7. However, in the case of the comb the integral equation (64) is hard to solve as it stands, and thus we proceed differently.

First, we use the causality structure of the response function and of the stationary correlator to reduce the time dependence in equation (64) to a single variable. To this end, we step to the Laplace domain, where the response functions in equation (63) take the form

$$S_{q,p}(s, s') = \frac{\tilde{R}_{q,p}(s) + \tilde{R}_{p,q}(s')}{s + s'} = \frac{\tilde{R}_{q,p}(s) + \tilde{R}_{q,p}(s')}{s + s'}. \quad (A.41)$$

In the first step we used the causal form of the response function $R_{q,p}(s, s')$ reported in equation (42), while in the second step we used the symmetry under $q \leftrightarrow p$ of the particular solution found in equation (41). Similarly, the correlator $C_{q,p}(s, s')$ can be written in terms of $C_{q,p}(s)$ as in equation (48). By plugging this decomposition into equation (64) we obtain two distinct (and equivalent) relations, depending either on $s$ or $s'$, the first of which reads

$$C_{q,p}(s) = G_{q,p} \left\{ \tilde{R}_{q,p}(s) - \mu \int \frac{d^{d-1}k_{\perp}}{(2\pi)^{d-1}} (T + u_{q_1, -k_{\perp}}) [q_1^2 C_{(q_1, k_{\perp}), p}(s) + p_1^2 C_{q, (p_1, k_{\perp})}(s)] \right\}. \quad (A.42)$$

Next, we note that equation (A.42) fixes completely the diagonal part of the correlator $C_{q,p}(s)$ in momentum space. Indeed, let us decompose (without loss of generality)

$$C_{q,p}(s) = c_{q}^{(1)}(s) \delta^{d}(\mathbf{q} + \mathbf{p}) + c_{q,p}^{(2)}(s) \delta(q_1 + p_1). \quad (A.43)$$

We have additionally assumed that $C_{q,p}(s) \propto \delta(q_1 + p_1)$, as expected from the translational invariance of the system along the backbone. Plugging this form into equation (A.42) renders

$$c_{q}^{(1)}(s) = \frac{T}{T + u_{q}} g_{\perp}(q, s), \quad (A.44)$$

$$c_{q,p}^{(2)}(s) = \mu q_1^2 G_{q,p} \left\{ 2Tg_{\perp}(\mathbf{p}, s) s G_{q}(s) - \left[ (T + u_{q}) c_{q}^{(1)}(s) + (T + u_{p}) c_{p}^{(1)}(s) \right] \right. \\ \left. - \int \frac{d^d k}{(2\pi)^d} \delta(k_1 - q_1) (T + u_{-k}) \left[ c_{k,p}^{(2)}(s) + c_{q,k}^{(2)}(s) \right] \right\}. \quad (A.45)$$

Note that the diagonal dependence encoded in $c_{q}^{(1)}(s)$ is exactly the same as that of the equilibrium correlator reported in equation (58). This suggests to express the off-diagonal part $c_{q,p}^{(2)}(s)$ as a deviation with respect to the equilibrium case: we thus introduce

$$c_{q,p}^{(2)}(s) = -\mu q_1^2 T g_{\perp}(\mathbf{p}, s) G_{q}(s) + c_{q,p}^{(3)}(s). \quad (A.46)$$

We then plug this into equation (A.45), and use equation (A.44) to recognize

$$[(T + u_q) c_q^{(1)}(s) + (T + u_p) c_p^{(1)}(s)] = T[g_{\perp}(q, s) + g_{\perp}(p, s)]$$
$$= g_{\perp}(\mathbf{q}, s) g_{\perp}(\mathbf{p}, s) [2s + \mathcal{G}_{\mathbf{q}, \mathbf{p}}^{-1}], \quad (A.47)$$

where in the last step we used the definition of $g_{\perp}$ given in equation (33), and that of $\mathcal{G}_{\mathbf{q}, \mathbf{p}}$ in equation (62). This leads, after some algebra, to a simpler relation of the form

$$c_{\mathbf{q}, \mathbf{p}}^{(3)}(s) = -\mu q_1^2 \mathcal{G}_{\mathbf{q}, \mathbf{p}} \int \frac{d^d k}{(2\pi)^d} \delta(k_1 - q_1) (T + u_{-\mathbf{k}}) [c_{\mathbf{k}, \mathbf{p}}^{(3)}(s) + c_{\mathbf{q}, \mathbf{k}}^{(3)}(s)]. \quad (A.48)$$

The latter integral equation is homogeneous, and thus it is clearly solved by $c_{\mathbf{q}, \mathbf{p}}^{(3)}(s) \equiv 0$. In the following, we will prove that such solution is unique: as a consequence, the equilibrium correlator given in equation (58) is the unique solution of equation (A.42), which defines the stationary correlator. To this end, we introduce another transformation as

$$c_{\mathbf{q}, \mathbf{p}}^{(4)}(s) \equiv \mathcal{G}_{\mathbf{q}, \mathbf{p}}^{-1} c_{\mathbf{q}, \mathbf{p}}^{(3)}(s), \quad (A.49)$$

in terms of which equation (A.48) becomes

$$c_{\mathbf{q}, \mathbf{p}}^{(4)}(s) = -\mu q_1^2 \int \frac{d^d k}{(2\pi)^d} \delta(k_1 - q_1) (T + u_{-\mathbf{k}}) [\mathcal{G}_{\mathbf{k}, \mathbf{p}} c_{\mathbf{k}, \mathbf{p}}^{(4)}(s) + \mathcal{G}_{\mathbf{q}, \mathbf{k}} c_{\mathbf{q}, \mathbf{k}}^{(4)}(s)]. \quad (A.50)$$

This form makes clear that the dependence of $c_{\mathbf{q}, \mathbf{p}}^{(4)}(s)$ on $\mathbf{q}$ and $\mathbf{p}$ is additive, namely

$$c_{\mathbf{q}, \mathbf{p}}^{(4)}(s) = f_{\mathbf{q}}(s) + f_{\mathbf{p}}(s), \quad (A.51)$$

for some unknown function $f_{\mathbf{q}}(s)$. Plugging this decomposition into equation (A.50) and introducing the auxiliary function

$$N_{\mathbf{q}} \equiv \mu q_1^2 \int \frac{d^d k}{(2\pi)^d} \delta(k_1 - q_1) (T + u_{-\mathbf{k}}) \mathcal{G}_{\mathbf{k}, \mathbf{q}}, \quad (A.52)$$

we obtain (using the parity of $\mathcal{G}_{\mathbf{q}, \mathbf{p}}$)

$$\begin{aligned} & \left[ (1 + N_{\mathbf{q}}) f_{\mathbf{q}}(s) + \mu q_1^2 \int \frac{d^d k}{(2\pi)^d} \delta(k_1 - q_1) (T + u_{-\mathbf{k}}) \mathcal{G}_{\mathbf{k}, \mathbf{q}} f_{\mathbf{k}}(s) \right] \\ & = - \left[ (1 + N_{\mathbf{p}}) f_{\mathbf{p}}(s) + \mu q_1^2 \int \frac{d^d k}{(2\pi)^d} \delta(k_1 - q_1) (T + u_{-\mathbf{k}}) \mathcal{G}_{\mathbf{k}, \mathbf{p}} f_{\mathbf{k}}(s) \right]. \end{aligned} \quad (A.53)$$

The left-hand side of this equation depends on $\mathbf{q}_{\perp}$, whereas the right-hand side depends on $\mathbf{p}_{\perp}$—hence they must be equal to a quantity that is independent of both $\mathbf{q}_{\perp}$ and $\mathbf{p}_{\perp}$:

$$(1 + N_{\mathbf{q}}) f_{\mathbf{q}}(s) + \mu q_1^2 \int \frac{d^d k}{(2\pi)^d} \delta(k_1 - q_1) (T + u_{-\mathbf{k}}) \mathcal{G}_{\mathbf{k}, \mathbf{q}} f_{\mathbf{k}}(s) = a(q_1, s), \quad (A.54)$$

$$(1 + N_{\mathbf{p}}) f_{\mathbf{p}}(s) + \mu q_1^2 \int \frac{d^d k}{(2\pi)^d} \delta(k_1 - q_1) (T + u_{-\mathbf{k}}) \mathcal{G}_{\mathbf{k}, \mathbf{p}} f_{\mathbf{k}}(s) = -a(q_1, s). \quad (A.55)$$
By symmetry under the exchange $\mathbf{q}_{\perp} \leftrightarrow \mathbf{p}_{\perp}$, we argue that it can only be $a(q_1, s) = 0$. This way we have reduced equation (64) to the single-variable integral equation

$$ (1 + N_{\mathbf{q}}) f_{\mathbf{q}}(s) + \mu q_1^2 \int \frac{d^d k}{(2\pi)^d} \delta(k_1 - q_1) (T + u_{-\mathbf{k}}) G_{\mathbf{k}, \mathbf{q}} f_{\mathbf{k}}(s) = 0. \quad (A.56) $$

Taking another convolution against $(T + u_{-\mathbf{q}})$ finally renders

$$ \int \frac{d^d q}{(2\pi)^d} (T + u_{-\mathbf{q}}) (1 + 2N_{\mathbf{q}}) f_{\mathbf{q}}(s) = 0, \quad (A.57) $$

which is only satisfied if $f_{\mathbf{q}}(s) = 0$. This concludes the proof of the uniqueness.

### A.6. Double Laplace transform of causal/stationary functions

The causality structure of a function in the time domain carries signatures on its Laplace transform. Indeed, it is straightforward to prove that [60–62]

$$ F(t_1, t_2) = f(t_2 - t_1) \Theta(t_2 - t_1) \quad \rightarrow \quad \hat{F}(s_1, s_2) = \frac{\hat{f}(s_2)}{s_1 + s_2}, \quad (A.58) $$

and similarly

$$ \begin{aligned} G(t_1, t_2) &= g(|t_2 - t_1|) = g(t_2 - t_1) \Theta(t_2 - t_1) + g(t_1 - t_2) \Theta(t_1 - t_2) \\ &\rightarrow \hat{G}(s_1, s_2) = \frac{\hat{g}(s_1) + \hat{g}(s_2)}{s_1 + s_2}. \end{aligned} \quad (A.59) $$

This simplifies significantly the inverse Laplace transformation of the response function $R_{\mathbf{q},\mathbf{p}}(s, s')$ and the stationary correlator $C_{\mathbf{q},\mathbf{p}}(s, s')$ derived in appendices A.2 and A.4, respectively: we basically have to transform only with respect to $s$, and not with respect to both $s$ and $s'$. Unfortunately this is not the case for the correlator $C_{\mathbf{q},\mathbf{p}}(s, s')$ in equation (44) for Dirichlet initial conditions, because the corresponding $C_{\mathbf{q},\mathbf{p}}(t, t')$ is not expected to be time-translational invariant.

### A.7 Particles in homogeneous space: analysis of a fluctuating Gaussian field in the Laplace domain

Here we consider the relaxational dynamics of a Gaussian field [52],

$$ \dot{\phi}_{\mathbf{q}}(t) = -\alpha_{\mathbf{q}} \phi_{\mathbf{q}}(t) + \nu_{\mathbf{q}}(t), \quad (A.60) $$

$$ \langle \nu_{\mathbf{q}}(t) \nu_{\mathbf{q}'}(t') \rangle = \Omega_{\mathbf{q}} \delta^d(\mathbf{q} + \mathbf{q}') \delta(t - t'), \quad (A.61) $$

to which equation (27) for the fluctuating particle density essentially reduces in homogeneous space (i.e. for $\hat{\mu}(\mathbf{x}) \equiv \mu \mathbb{1}_d$), upon identifying $\alpha_{\mathbf{q}} = \mu q^2 (T + u_{\mathbf{q}})$ and $\Omega_{\mathbf{q}} = 2\mu T q^2$. Although this problem is easily solvable in the time domain, here we will characterize it instead in the Laplace domain. This serves as a benchmark for the methods employed above in the case of the comb, which in contrast cannot be analyzed directly in the time domain.

First, the propagator corresponding to equation (A.60) can be immediately checked to give

$$ G_{\mathbf{q}}(s) = \frac{1}{s + \alpha_{\mathbf{q}}} \quad \leftrightarrow \quad G_{\mathbf{q}}(t) = \Theta(t) e^{-\alpha_{\mathbf{q}} t}. \quad (A.62) $$
To find the correlation function, we first write the Schwinger–Dyson equations (see section 4.2)

$$ \partial_t C_{q,p}(t, t') = -\alpha_q C_{q,p}(t, t') + R_{p,q}(t', t), \qquad (A.63) $$

$$ \partial_t R_{q,p}(t, t') = -\alpha_q R_{q,p}(t, t') + \Omega_q \delta^d(\mathbf{q} + \mathbf{p}) \delta(t - t'). \qquad (A.64) $$

Solving the latter in the Laplace domain (as we did in appendix A.2 for the comb) yields

$$ R_{q,p}(s, s') = \frac{\Omega_q \delta^d(\mathbf{q} + \mathbf{p})}{(s + s')(s + \alpha_q)} = \frac{\Omega_q \delta^d(\mathbf{q} + \mathbf{p})}{(s + s')} G_q(s), \qquad (A.65) $$

which can be inverted to the time domain (with respect to both its Laplace variables) to give

$$ R_{q,p}(t, t') = \Omega_q \delta^d(\mathbf{q} + \mathbf{p}) \Theta(t') \Theta(t - t') e^{-\alpha_q (t - t')}. \qquad (A.66) $$

To find the correlator, we Laplace-transform equation (A.63) and plug in the solution (A.65) for the response function:

$$ C_{q,p}(s, s') = \frac{1}{s + \alpha_q} \left[ C_{q,p}(t = 0^+, s') + \frac{\Omega_q \delta^d(\mathbf{q} + \mathbf{p})}{(s + s')(s' + \alpha_q)} \right]. \qquad (A.67) $$

Transforming back to the time domain yields

$$ C_{q,p}(t, t') = C_{q,p}(t = 0^+, t') e^{-\alpha_q t'} + \frac{\Omega_q \delta^d(\mathbf{q} + \mathbf{p})}{2\alpha_q} \left[ e^{-\alpha_q |t-t'|} - e^{-\alpha_q (t+t')} \right], \qquad (A.68) $$

where the term $C_{q,p}(t = 0^+, t')$ must be fixed self-consistently by using the initial conditions

$$ C_{q,p}(t = t_0, t') = \phi_q(t_0) \phi_p(t_0) G_p(t' - t_0), \qquad (A.69) $$

specified in correspondence of $t_0 < t'$. Indeed, computing equation (A.68) in $t = t_0$ and solving for $C_{q,p}(t = 0^+, t')$ gives

$$ C_{q,p}(t = 0^+, t') = \phi_q(t_0) \phi_p(t_0) e^{-\alpha_q (t' - 2t_0)} + \frac{\Omega_q \delta^d(\mathbf{q} + \mathbf{p})}{2\alpha_q} e^{-\alpha_q t'} (1 - e^{2\alpha_q t_0}), \qquad (A.70) $$

which can then be replaced back into equation (A.68) to find the Dirichlet correlator [52]

$$ C_{q,p}(t, t') = \phi_q(t_0) \phi_p(t_0) e^{-\alpha_q (t+t'-2t_0)} + \frac{\Omega_q \delta^d(\mathbf{q} + \mathbf{p})}{2\alpha_q} \left[ e^{-\alpha_q |t-t'|} - e^{-\alpha_q (t+t'-2t_0)} \right]. \qquad (A.71) $$

Note that, starting from the flat initial condition $\phi_q(0) = 0$, the Dirichlet correlator reads in the Laplace domain

$$ C_{q,p}(s, s') = \frac{\Omega_q \delta^d(\mathbf{q} + \mathbf{p})}{(s + s')(s + \alpha_q)(s' + \alpha_q)}. \qquad (A.72) $$

Note also that, unfortunately, an explicit Laplace expression of the correlator starting from initial conditions set at a generic time $t_0$ is not available. Indeed, $t=0$ plays a special role for the Laplace transform, which in some sense breaks time-translational invariance. This is why in appendix A.3 we had to fix the initial conditions at time $t_0 = 0$.
Finally, to find the stationary correlator in the Laplace domain we start by writing another Schwinger–Dyson equation, akin to equation (A.63):

$$ \partial_{t'} C_{q,p}(t, t') = -\alpha_p C_{q,p}(t, t') + R_{q,p}(t, t'). \qquad (A.73) $$

We then impose the stationarity condition

$$ (\partial_t + \partial_{t'}) C_{q,p}(t, t') \equiv 0 \quad \rightarrow \quad C_{q,p}(t, t') = \frac{R_{q,p}(t, t') + R_{p,q}(t', t)}{\alpha_q + \alpha_p}. \qquad (A.74) $$

Plugging in the Laplace expression (A.65) of the response function yields

$$ C_{q,p}(s, s') = \frac{\Omega_q \delta^d(\mathbf{q} + \mathbf{p})}{2\alpha_q(s + s')} \left( \frac{1}{s + \alpha_q} + \frac{1}{s' + \alpha_q} \right), \qquad (A.75) $$

whose time-domain version

$$ C_{q,p}(t, t') = \frac{\Omega_q \delta^d(\mathbf{q} + \mathbf{p})}{2\alpha_q} e^{-\alpha_q |t - t'|} \qquad (A.76) $$

clearly coincides with the formal limit $t_0 \to -\infty$ taken in equation (A.71). It also coincides with the equilibrium correlator obtained by using the fluctuation–dissipation theorem [47, 52].

## Appendix B. Structure factor and ISF

In this appendix we recall some basic concepts from the physics of liquids [63], so as to make contact with the notation used in SDFT.

Let us first introduce the (static) structure factor

$$ S(\mathbf{q}) = \frac{1}{N} \sum_{i,j=1}^{N} \langle e^{-i\mathbf{q} \cdot (\mathbf{r}_i - \mathbf{r}_j)} \rangle = \frac{1}{N} \langle \rho_{\mathbf{q}} \rho_{-\mathbf{q}} \rangle, \qquad (B.1) $$

where in the second step we used the definition of the density field $\rho$ given in equation (8). First, note that one can always decompose equation (B.1) as

$$ S(\mathbf{q}) = 1 + \frac{1}{N} \sum_{i \neq j} \langle e^{-i\mathbf{q} \cdot (\mathbf{r}_i - \mathbf{r}_j)} \rangle, \qquad (B.2) $$

where the first term accounts for self-correlations. For large $\mathbf{q}$, the second term rapidly oscillates and averages out to zero, so that $S(\mathbf{q}) \sim 1$; moreover, from equation (B.1) it follows that $S(\mathbf{q} = \mathbf{0}) = N$. As such, by construction, $S(\mathbf{q})$ does not admit an inverse Fourier transform. Using equation (B.1) and the definition of the density fluctuation $\phi$ introduced in equation (15), we can rewrite

$$ \begin{aligned} S(\mathbf{q}) &= \rho_0 \left[ \delta^d(\mathbf{q}) + \frac{1}{N} \langle \phi_{\mathbf{q}} \phi_{-\mathbf{q}} \rangle \right] \\ &= 1 + \rho_0 \delta^d(\mathbf{q}) + \rho_0 h(\mathbf{q}). \end{aligned} \qquad (B.3) $$

Above, in the first line we used $[\delta^d(\mathbf{q})]^2 = V\delta^d(\mathbf{q})$, where $V$ denotes the volume of the system (i.e. $\rho_0 = N/V$), while the second line defines the pair correlation function $h(\mathbf{q})$—note that
the latter does admit an inverse Fourier transform [35, 63]. Note also that the term $\rho_0\delta^d(\mathbf{q})$ in equation (B.3) is only due to the presence of a uniform nonzero background density $\rho_0$; since this only entails a difference in $\mathbf{q} = \mathbf{0}$, one typically rather plots

$$ \mathcal{G}(\mathbf{q}) = S(\mathbf{q}) - \rho_0 \delta^d(\mathbf{q}) = \frac{\rho_0}{N} \langle \phi_{\mathbf{q}} \phi_{-\mathbf{q}} \rangle. \qquad (B.4) $$

Using equation (57), we deduce that in our case

$$ \mathcal{G}(\mathbf{q}) = \frac{\rho_0}{N} C_{\mathbf{q},-\mathbf{q}}(\tau=0^+) = \frac{T}{T+u_{\mathbf{q}}}, \qquad (B.5) $$

where again we used $\delta^d(\mathbf{q} = \mathbf{0}) = \mathcal{V}$. As noted in earlier works [35, 64], this result—stemming from the linearization of the Dean–Kawasaki equation that we performed in section 3— corresponds to the random phase (or Debye–Hückel) approximation in the physics of liquids or electrolytes [63]. Under such approximation, we note that the property $\mathcal{G}(\mathbf{0}) = 0$ is lost, whereas $\mathcal{G}(\mathbf{q} \to \infty) = 1$, as expected.

Similarly, the temporal evolution of density fluctuations can be described using the intermediate scattering function [63]

$$ F(\mathbf{q}, t) = \frac{1}{N} \sum_{i,j=1}^{N} \langle e^{-i\mathbf{q} \cdot (\mathbf{r}_i(t) - \mathbf{r}_j(0))} \rangle = \frac{1}{N} \langle \rho_{\mathbf{q}}(t) \rho_{-\mathbf{q}}(0) \rangle. \qquad (B.6) $$

Similar considerations to those spelled out above for the structure factor apply also here, so that in practice it is useful to subtract the constant background and plot instead

$$ \mathcal{F}(\mathbf{q}, t) = F(\mathbf{q}, t) - \rho_0 \delta^d(\mathbf{q}) = \frac{\rho_0}{N} \langle \phi_{\mathbf{q}}(t) \phi_{-\mathbf{q}}(0) \rangle = \frac{\rho_0}{N} C_{\mathbf{q},-\mathbf{q}}(t), \qquad (B.7) $$

where in the last step we recognized the two-point function introduced in equations (37) and (47). For the comb geometry, the latter is only known explicitly in the Laplace domain, i.e. $C_{\mathbf{q},\mathbf{p}}(s)$ given in equation (58). However, its form suggests to rescale the momenta as $z_1^2 = q_1^2/\sqrt{s}$ and $z_\perp^2 = q_\perp^2/s$; in particular, the Fourier transform of the interaction potential $u(\mathbf{x})$ this way reads

$$ u_{\mathbf{q}} = u_{(z_1 s^{1/4}, z_\perp s^{1/2})}. \qquad (B.8) $$

Upon formally taking the inverse Laplace transform, it then becomes evident that $C_{\mathbf{q},-\mathbf{q}}(t)$ col lapses, for large times $t$, on a scaling function that does not depend on the details of the interac tion potential $u_{\mathbf{q}}$, but merely on $u_{\mathbf{q}=0}$—see equation (B.8) with $s \sim 0$. As noted in appendix A.4, this is precisely the limit in which the inverse Laplace transform can be computed explicitly: this way we find, for $d=2$,

$$ \begin{aligned} \mathcal{F}(\mathbf{q}, t) & \sim \mathcal{F}(y_1 = q_1 t^{\frac{1}{4}}, y_2 = q_2 t^{\frac{1}{2}}) \\ &= \frac{T}{\mathcal{T}} e^{-\mu \mathcal{T} y_2^2} \left\{ 1 + \frac{4y_1^2}{L(y_1^4 + 4y_2^2)} \left[ \frac{y_1^4}{(y_1^4 + 4y_2^2)} \left( 1 - e^{\mu \mathcal{T} (y_2^2 + y_1^4/4)} \text{Erfc} \left( \frac{y_1^2}{2\sqrt{\mu \mathcal{T}}} \right) \right) \right. \right. \\ & \qquad \left. \left. - y_2^2 \mu \mathcal{T} - \frac{e^{\mu \mathcal{T} y_2^2} y_1^2 \sqrt{\mu \mathcal{T}}}{2\sqrt{\pi}} + \frac{y_1^2 (4y_2^2 - y_1^4 + 2y_2^2 (y_1^4 + 4y_2^2) \mu \mathcal{T})}{y_2 (y_1^4 + 4y_2^2)} \text{Erfi} \left( 4y_2 \sqrt{\mu \mathcal{T}} \right) \right] \right\}, \end{aligned} \qquad (B.9) $$
where we called $\mathcal{T} = T+u_{q=0}$. This function is plotted in figure 2. In particular, note that the term $\propto \delta(q_1+p_1)$ in equation (58) gives rise to a term $\propto L$ when computed for $\mathbf{p}=\mathbf{q}$ (as prescribed by equation (B.7)), where $L$ is the linear size of the system, whereas the term $\propto \delta^d(\mathbf{q}+\mathbf{p})$ gives rise to a term $\propto \mathcal{V} = L^d$. Intuitively, this is because the structure factor and the ISF defined above entail averages over all particles in the system, among which only a subextensive fraction resides on the backbone. As a consequence, the correction to the ISF due to the particles subdiffusing along the backbone appears as a subleading contribution in equation (B.9), which is expected to vanish in the thermodynamic limit $\mathcal{V} \to \infty$ (while it remains relevant in a finite system, see figure 2).

## Appendix C. Details on the generalized correlation profiles

Here we detail the derivation of the generalized correlation profiles between the tracer position and the density of all other bath particles, discussed in section 5.1.

Several steps in the following derivation resemble the ones carried out in section II.D of the supplemental material of [42], to which we refer the reader. In particular, the profiles $\psi_q$ and $g_q$ introduced in equations (67) and (69) can be computed directly by using these very same methods, starting from the coupled equations of motion (65) and (66) for the tracer and the bath density. In the following, we limit ourselves for simplicity to the derivation of the generating function $w_q(\lambda, t)$ introduced in equation (70), from which the previous two profiles can actually be generated.

### C.1. Derivation of the generating function

We start from the definition (70) of the generating function $w_q(\lambda, t)$, and note that

$$w_q(\lambda, t) = \frac{\langle \phi_q(t) e^{(\lambda+iq_1)r_0(t)} \rangle}{\langle e^{\lambda r_0(t)} \rangle} = \frac{\langle \phi_q(t) e^{(\lambda+iq_1)r_0(t)} \rangle}{\langle e^{\lambda r_0(t)} \rangle_0} + \mathcal{O}(h^2), \quad (C.1)$$

meaning that the leading-order contribution to $w_q(\lambda, t)$ can be captured by computing the denominator only up to $\mathcal{O}(h^0)$. Since equation (65) reduces to 1d Brownian motion if $h=0$, such denominator simply reads [52]

$$\langle e^{\lambda r_0(t)} \rangle_0 = e^{\lambda^2 \mu_0 T t}. \quad (C.2)$$

At long times, we expect $w_q(\lambda, t)$ to reach a stationary state satisfying

$$0 \equiv \partial_t w_q(\lambda) = \frac{(\partial_t - \lambda^2 \mu_0 T) \langle \phi_q(t) e^{(\lambda+iq_1)r_0(t)} \rangle}{\langle e^{\lambda r_0(t)} \rangle_0} + \mathcal{O}(h^2), \quad (C.3)$$

where we used equations (C.1) and (C.2). (Note that the latter is an assumption, to be verified a posteriori.) We then focus on

$$\partial_t \langle \phi_q(t) e^{(\lambda+iq_1)r_0(t)} \rangle = \langle \dot{\phi}_q(t) e^{(\lambda+iq_1)r_0(t)} \rangle + (\lambda+iq_1) \langle \dot{r}_0(t) \phi_q(t) e^{(\lambda+iq_1)r_0(t)} \rangle, \quad (C.4)$$

which we obtained by using Stratonovich calculus. Indeed, note that both noise terms in equations (65) and (66) are additive, hence the result of the calculation does not depend on the choice of the stochastic calculus convention. The two terms on the r.h.s. of equation (C.4)
can then be computed by using the coupled equations of motion; in particular, the first term gives

$$ \frac{\langle \dot{\phi}_{\mathbf{q}}(t) e^{(\lambda+iq_1)r_0(t)} \rangle}{\langle e^{\lambda r_0(t)} \rangle_0} = -\mu q_{\perp}^2 (T+u_{\mathbf{q}}) w_{\mathbf{q}}(\lambda) - \mu q_1^2 \int \frac{d^{d-1}q_{\perp}}{(2\pi)^{d-1}} (T+u_{\mathbf{q}}) w_{\mathbf{q}}(\lambda) \\ - h\mu [q_{\perp}^2 v_{\mathbf{q}} + q_1^2 v_{q_1}(\mathbf{x}_{\perp}=\mathbf{0})], \qquad (C.5) $$

because by using the Novikov theorem [56, 65] one can prove that [42]

$$ \langle \nu_{\mathbf{q}}(t) e^{(\lambda+iq_1)r_0(t)} \rangle = 0. \qquad (C.6) $$

Similarly, the second term yields

$$ \frac{\langle \dot{r}_0(t) \phi_{\mathbf{q}}(t) e^{(\lambda+iq_1)r_0(t)} \rangle}{\langle e^{\lambda r_0(t)} \rangle_0} \\ = \frac{\langle \eta_0(t) \phi_{\mathbf{q}}(t) e^{(\lambda+iq_1)r_0(t)} \rangle}{\langle e^{\lambda r_0(t)} \rangle_0} - h\mu_0 \int \frac{d^d p}{(2\pi)^d} i p_1 v_{\mathbf{p}} \frac{\langle \phi_{\mathbf{p}}(t) \phi_{\mathbf{q}}(t) e^{(\lambda+iq_1+ip_1)r_0(t)} \rangle}{\langle e^{\lambda r_0(t)} \rangle_0} \\ = (iq_1 + \lambda) \mu_0 T w_{\mathbf{q}}(\lambda) + iq_1 h \mu_0 T \frac{v_{-\mathbf{q}}}{T+u_{\mathbf{q}}}. \qquad (C.7) $$

Above we noted that

$$ \langle \phi_{\mathbf{p}}(t) \phi_{\mathbf{q}}(t) e^{(\lambda+iq_1+ip_1)r_0(t)} \rangle = \langle \phi_{\mathbf{p}}(t) \phi_{\mathbf{q}}(t) \rangle_0 \langle e^{(\lambda+iq_1+ip_1)r_0(t)} \rangle_0 + \mathcal{O}(h) \\ = \frac{T}{T+u_{\mathbf{q}}} \delta^d(\mathbf{q}+\mathbf{p}) \langle e^{\lambda r_0(t)} \rangle_0 + \mathcal{O}(h), \qquad (C.8) $$

where in the first step we used that the stochastic processes $r_0(t)$ and $\phi_{\mathbf{q}}(t)$ are independent at $\mathcal{O}(h^0)$ and thus they factorize, while in the second step we inserted the stationary two-point function of the field found in equation (57). Moreover, the expectation value involving the noise $\eta_0(t)$ can be derived as in [42] by using the Novikov theorem, yielding

$$ \langle \eta_0(t) \phi_{\mathbf{q}}(t) e^{(\lambda+iq_1)r_0(t)} \rangle = (iq_1 + \lambda) \mu_0 T \langle \phi_{\mathbf{q}}(t) e^{(\lambda+iq_1)r_0(t)} \rangle. \qquad (C.9) $$

Plugging equations (C.5) and (C.7) back into equation (C.3) yields the defining equation of the stationary $w_{\mathbf{q}}(\lambda)$, namely

$$ G_0^{-1}(\mathbf{q}, \lambda) w_{\mathbf{q}}(\lambda) = -\mu q_1^2 \int \frac{d^{d-1}q_{\perp}}{(2\pi)^{d-1}} (T+u_{\mathbf{q}}) w_{\mathbf{q}}(\lambda) - h\mu [q_{\perp}^2 v_{\mathbf{q}} + q_1^2 v_{q_1}(\mathbf{x}_{\perp}=\mathbf{0})] \\ + iq_1 h \mu_0 T \frac{v_{-\mathbf{q}}}{T+u_{\mathbf{q}}}, \qquad (C.10) $$

where the function $G_0(\mathbf{q}, \lambda)$ was introduced in equation (72). Equation (C.10) can as usual be solved self-consistently by applying $\int \frac{d^{d-1}q_{\perp}}{(2\pi)^{d-1}} (T+u_{\mathbf{q}})$ to both its sides; a tedious but straightforward calculation then renders the expression of $w_{\mathbf{q}}(\lambda)$ reported in equation (71).
## C.2. Derivation and large-distance behavior of correlation profiles

Other generalized correlation profiles can be promptly generated starting from $w_q(\lambda)$ given in equation (71). For instance, the profile $\psi_q$ introduced in equation (67) can be simply computed as $\psi_q = w_q(\lambda = 0)$. Similarly, since $g_q = \frac{dw_q(\lambda)}{d\lambda}|_{\lambda=0}$ (see equation (69)), by deriving implicitly equation (C.10) we obtain

$$G_0^{-1}(\mathbf{q}, \lambda = 0) g_{\mathbf{q}} = -\mu q_1^2 \int \frac{d^{d-1}q_{\perp}}{(2\pi)^{d-1}} (T+u_{\mathbf{q}}) g_{\mathbf{q}} + iq_1\mu_0 T \left[ \frac{hv_{-\mathbf{q}}}{T+u_{\mathbf{q}}} + 2\psi_{\mathbf{q}} \right]. \quad (C.11)$$

The self-consistent solution of equation (C.11) renders the stationary profile $g_q$, which is reported in equation (75). Its large-distance behavior can then be inspected e.g. by using the method presented in section I.C of the supplementary material in [42], or more simply by formally replacing the interaction potential by $u(\mathbf{x}) \mapsto u_{\mathbf{q}=0}\delta(\mathbf{x})$ (indeed, note that we have assumed everywhere the interaction potential $u(\mathbf{x})$ to decay rapidly, so that it admits a Fourier transform). Either way leads to the asymptotic result reported in equation (78).

## C.3. On higher-order cumulants

In section 5.2 we noted that our prediction for the cumulant generating function (85) of the tracer position seems to imply that all cumulants of $r_0(t)$ grow linearly with time. Crucially, however, equation (85) has been derived under the assumption that $w_q(\lambda, t)$ becomes stationary. To see why this assumption may be delicate, in this appendix we consider the simpler case of interacting particles in uniform space (i.e. without the comb constraint), such as the ones considered in e.g. [35, 42]—and corresponding to the model described in section 2, but with $\hat{\mu}(\mathbf{x}) = \hat{\mu}_0(\mathbf{x}) = \mu \mathbf{1}_d$, and $U_0(\mathbf{x}) = U(\mathbf{x})$. First, using Stratonovich calculus one can derive the exact relation [42]

$$\partial_t \Psi(\lambda, t) = \lambda^2 \mu T - h \mu \lambda \cdot \int \frac{d^d q}{(2\pi)^d} i \mathbf{q} u_{\mathbf{q}} w_{\mathbf{q}}(\lambda, t), \quad (C.12)$$

which is analogous to equation (84) but for $\Psi(\lambda, t) = \ln \langle e^{\lambda \cdot r_0(t)} \rangle$, and with the generating function

$$w(\mathbf{x}, \lambda, t) \equiv \frac{\langle \phi(\mathbf{x} + \mathbf{r}_0(t), t) e^{\lambda \cdot \mathbf{r}_0(t)} \rangle}{\langle e^{\lambda \cdot \mathbf{r}_0(t)} \rangle}. \quad (C.13)$$

Following analogous steps to the ones we presented in appendix C.1, one finds in the stationary state

$$w_q(\lambda) = \frac{-h u_q}{T + u_q} \left[ 1 + \frac{T i q \cdot \lambda}{q^2 (2T + u_q) - 2 T i q \cdot \lambda} \right] + \mathcal{O}(h^2), \quad (C.14)$$

where $q = |\mathbf{q}|$. From equation (C.12), this would imply

$$\Psi(\lambda, t) = \mu t \left[ \lambda^2 T - h \lambda \cdot \int \frac{d^d q}{(2\pi)^d} i \mathbf{q} u_{\mathbf{q}} w_{\mathbf{q}}(\lambda) \right]. \quad (C.15)$$
By inserting the prediction for $w_q(\lambda)$ found in equation (C.14), we can obtain estimates of the cumulants of $\mathbf{r}_0(t)$ up to and including $O(h^2)$, upon taking derivatives with respect to $\lambda$. For instance, the fourth cumulant would read

$$ \frac{\langle r_0^A(t) \rangle_c}{t} = 96h^2\mu T^3 \int \frac{d^d q}{(2\pi)^d} \frac{u_q^2 q_1^4}{(T+u_q)(2T+u_q)^3 q^6} + O(h^4), \quad (C.16) $$

where we denoted $r_0 = \hat{e}_1 \cdot \mathbf{r}_0$, while $\langle \cdot \rangle_c$ is the connected correlation function. However, we note that in $d \le 2$ this integral requires regularization via the introduction of a lower (infrared) cutoff $q \ge 2\pi/L$, implying that the fourth cumulant depends explicitly on the size $L$ of the system—in particular, the integral in equation (C.16) becomes divergent for $L \to \infty$, signaling the breakdown of our perturbative estimate.

In fact, our prediction in equation (C.15) has to be compared with equation (33) in [66], where a cumulant generating function analogous to $\Psi(\lambda, t)$ has been computed perturbatively for an analogous system, and for a generic time $t$ (i.e. not necessarily in the stationary state). In particular, the fourth cumulant reported in equation (52) of [66] coincides with the one in equation (C.16) only for $d > 2$, whereas for $d \le 2$ the scaling of $\langle r_0^A(t) \rangle_c$ with $t$ turns out to be faster than linear.

We conclude that, in general, higher-order cumulants of $r_0(t)$ only scale as $\propto t$ provided that the integral that defines the prefactor is convergent [45], which must (somewhat inconveniently) be checked a posteriori. By contrast, the integral defining the diffusion coefficient turns out to be always finite, no matter the spatial dimension $d$, both in uniform space [35] and for the case of the comb analyzed here (see section 5.2).

## Appendix D. Brownian dynamics simulations

In this appendix we describe the Brownian dynamics simulation of the system of interacting particles introduced in section 2. To this end, we focus on dimension $d=2$, and call $\mathbf{r}_i(t) = (X_i(t), Y_i(t))$. When trying to simulate the system on a delta-like comb such as the one introduced in section 2, the first obvious problem is that one never has $Y_i(t) = 0$ exactly, hence the particles would never actually diffuse horizontally. However, there are well-established strategies in the literature to circumvent this issue: for instance, in [17] the backbone is represented by a strip of spatial width $\epsilon_c$, while in [19, 67, 68] the $\delta(y)$ function that appears in the mobility matrix (6) is smoothed out and replaced by a Gaussian. The latter strategy is the one we chose to adopt here: in practice, we simulate

$$ \dot{X}_i(t) = -\mu \mathcal{A}(Y_i(t)) \sum_{j \neq i} \frac{\partial}{\partial X_i} U(\mathbf{r}_i(t) - \mathbf{r}_j(t)) + \sqrt{2\mu T \mathcal{A}(Y_i(t))} \eta_i^x(t), \quad (D.1) $$

$$ \dot{Y}_i(t) = -\mu \sum_{j \neq i} \frac{\partial}{\partial Y_i} U(\mathbf{r}_i(t) - \mathbf{r}_j(t)) + \sqrt{2\mu T \eta_i^y(t)}, \quad (D.2) $$

where the $\eta_i$ are uncorrelated random variables with

$$ \langle \eta_i^x(t) \eta_j^x(t') \rangle = \delta_{ij} \delta(t-t'), \quad (D.3) $$

and similarly for $\eta_i^y(t)$. The function $\mathcal{A}(y)$ can in principle be any smooth representation of a nascent delta function, such as

$$ \mathcal{A}(y) = \frac{1}{\sqrt{2\pi\epsilon_c}} \exp\left(-\frac{y^2}{2\epsilon_c^2}\right). \quad (D.4) $$
Choosing an integration time step $\Delta t$ such that $\sqrt{2\mu T\Delta t} \lesssim 2\epsilon_c$ is in general sufficient to make the measurements independent of $\epsilon_c$ [19]. (We verified this by running several simulations with decreasing values of $\epsilon_c$, until we found that the $\epsilon_c$-dependence had actually been lost). Note that equations (D.1) and (D.2) are multivariate stochastic differential equations with multiplicative noise, that we chose to integrate using a simple Euler-Maruyama algorithm (indeed, higher-order methods would not significantly speed up the numerical integration in such cases [69, 70]).

In particular, choosing Gaussian interparticle potentials $U(r) = \epsilon e^{-\frac{r^2}{2R^2}}$, the force terms in equations (D.1) and (D.2) reduce to

$$ \frac{\partial U(\|\mathbf{r}_i - \mathbf{r}_j\|)}{\partial X_i} = - \frac{X_i - X_j}{R^2} U(\|\mathbf{r}_i - \mathbf{r}_j\|), \qquad (D.5) $$

$$ \frac{\partial U(\|\mathbf{r}_i - \mathbf{r}_j\|)}{\partial Y_i} = - \frac{Y_i - Y_j}{R^2} U(\|\mathbf{r}_i - \mathbf{r}_j\|). \qquad (D.6) $$

We start the simulation by preparing a system with 200 particles at density $\rho_0 = 1$, and fix all other parameters to unity apart from the interaction strength $\epsilon$, which we vary within the interval $\epsilon \in [0, 0.5]$. After a thermalization time $T_{\text{term}} = 10^5$ with a time step $\Delta t = 10^{-2}$, we record the trajectory of the tracer $r_0(t)$ from time $t=0$ up to time $T_{\text{max}} = 4 \cdot 10^5$, and repeat the whole process $n=10$ times. This way we end up with a collection of trajectories $\{r_0^{(i)}(t)\}_{i=1}^n$. We then use each of these trajectories to compute the time-averaged mean squared displacement [43, 71]

$$ \overline{\delta_i^2(\tau, T_{\text{max}})} = \frac{1}{T_{\text{max}} - \tau} \int_0^{T_{\text{max}} - \tau} dt |r_0^{(i)}(t + \tau) - r_0^{(i)}(t)|^2. \qquad (D.7) $$

Since our system is ergodic, we can then ensemble-average over the $n$ independent trajectories to obtain the estimate for the MSD

$$ \langle |r_0(t) - r_0(0)|^2 \rangle = \lim_{T_{\text{max}} \to \infty} \frac{1}{n} \sum_{i=1}^{n} \overline{\delta_i^2(\tau = t, T_{\text{max}})}, \qquad (D.8) $$

whose linear fit finally leads to the estimates of the diffusion coefficient $\langle |r_0(t) - r_0(0)|^2 \rangle \approx 2D_{\text{eff}}t$ reported in figure 3(b).

## ORCID iDs

Davide Venturelli https://orcid.org/0000-0001-6952-9568

Pierre Illien https://orcid.org/0000-0002-0769-2068

Aurélien Grabsch https://orcid.org/0000-0003-4316-5190

## References

[1] Höfling F and Franosch T 2013 *Rep. Prog. Phys.* **76** 046602

[2] Metzler R and Klafter J 2000 *Phys. Rep.* **339** 1–77

[3] Condamin S, Tejedor V, Voituriez R, Bénichou O and Klafter J 2008 *Proc. Natl Acad. Sci. USA* **105** 5675–80

[4] ben-Avraham D and Havlin S 2005 *Diffusion and Reactions in Fractals and Disordered Systems* (Cambridge University Press)
[5] Bouchaud J and Georges A 1990 *Phys. Rep.* **195** 127-293

[6] White S R and Barma M 1984 *J. Phys. A: Math. Gen.* **17** 2995-3008

[7] Weiss G H and Havlin S 1986 *Physica A* **134** 474-82

[8] Havlin S, Kiefer J E and Weiss G H 1987 *Phys. Rev. A* **36** 1403-8

[9] Arkhincreev V E and Baskin E 1991 *Sov. Phys. - JETP* **73** 161-5 (available at: http://www.jetp.ras.ru/cgi-bin/dn/e_073_01_0161.pdf)

[10] Iomin A 2006 *Phys. Rev. E* **73** 061918

[11] Frauenrath H 2005 *Prog. Polym. Sci.* **30** 325-84

[12] Méndez V and Iomin A 2013 *Chaos Solitons Fractals* **53** 46-51

[13] Sagi Y, Brook M, Almog I and Davidson N 2012 *Phys. Rev. Lett.* **108** 093002

[14] Arkhincreev V 2002 *Physica A* **307** 131-41

[15] Baskin E and Iomin A 2004 *Phys. Rev. Lett.* **93** 120603

[16] Lenzi E K, da Silva L R, Tateishi A A, Lenzi M K and Ribeiro H V 2013 *Phys. Rev. E* **87** 012121

[17] Ribeiro H V, Tateishi A A, Alves L G A, Zola R S and Lenzi E K 2014 *New J. Phys.* **16** 093050

[18] Iomin A 2011 *Phys. Rev. E* **83** 052106

[19] Trajanovski P, Jolakoski P, Zelenkovski K, Iomin A, Kocarev L and Sandev T 2023 *Phys. Rev. E* **107** 054129

[20] Burioni R and Cassi D 2005 *J. Phys. A: Math. Gen.* **38** R45-R78

[21] Villamaina D, Sarracino A, Gradenigo G, Puglisi A and Vulpiani A 2011 *J. Stat. Mech.* L01002

[22] Rebenshtok A and Barkai E 2013 *Phys. Rev. E* **88** 052126

[23] Agliari E, Blumen A and Cassi D 2014 *Phys. Rev. E* **89** 052147

[24] Agliari E, Sartori F, Cattivelli L and Cassi D 2015 *Phys. Rev. E* **91** 052132

[25] Illien P and Bénichou O 2016 *J. Phys. A: Math. Theor.* **49** 265001

[26] Iomin A, Méndez V and Horsthemke W 2018 *Fractional Dynamics in Comb-Like Structures* (World Scientific) (http://doi.org/10.1142/11076)

[27] Bénichou O, Illien P, Oshanin G, Sarracino A and Voituriez R 2015 *Phys. Rev. Lett.* **115** 220601

[28] Poncet A, Grabsch A, Bénichou O and Illien P 2022 *Phys. Rev. E* **105** 054139

[29] Grabsch A, Berlioz T, Rizkallah P, Illien P and Bénichou O 2024 *Phys. Rev. Lett.* **132** 037102

[30] Berlioz T, Venturelli D, Grabsch A and Bénichou O 2024 *J. Stat. Mech.* 113208

[31] Dean D S 1996 *J. Phys. A: Math. Gen.* **29** L613-7

[32] Kawasaki K 1998 *J. Stat. Phys.* **93** 527-46

[33] te Vrugt M, Löwen H and Wittkowski R 2020 *Adv. Phys.* **69** 121-247

[34] Illien P 2024 The Dean-Kawasaki equation and stochastic density functional theory (arXiv:2411.13467)

[35] Démery V, Bénichou O and Jacquin H 2014 *New J. Phys.* **16** 053032

[36] Dean D S and Podgornik R 2014 *Phys. Rev. E* **89** 032117

[37] Démery V and Dean D S 2016 *J. Stat. Mech.* 023106

[38] Poncet A, Bénichou O, Démery V and Oshanin G 2017 *Phys. Rev. Lett.* **118** 118002

[39] Krüger M and Dean D S 2017 *J. Chem. Phys.* **146** 134507

[40] Krüger M, Solon A, Démery V, Rohwer C M and Dean D S 2018 *J. Chem. Phys.* **148** 084503

[41] Mahdisoltani S and Golestanian R 2021 *Phys. Rev. Lett.* **126** 158002

[42] Venturelli D, Illien P, Grabsch A and Bénichou O 2024 Universal scale-free decay of tracer-bath correlations in d-dimensional interacting particle systems (arXiv:2411.09326)

[43] Muzzeddu P L, Kalz E, Gambassi A, Sharma A and Metzler R 2024 Self-diffusion anomalies of an odd tracer in soft-core media (arXiv:2411.15552)

[44] Cornalba F and Fischer J 2023 *Arch. Ration. Mech. Anal.* **247** 76

[45] Démery V and Dean D S 2011 *Phys. Rev. E* **84** 011148

[46] Venturelli D, Ferraro F and Gambassi A 2022 *Phys. Rev. E* **105** 054125

[47] Venturelli D and Gambassi A 2022 *Phys. Rev. E* **106** 044112

[48] Venturelli D and Gambassi A 2023 *New J. Phys.* **25** 093025

[49] Venturelli D, Loos S A M, Walter B, Roldán E and Gambassi A 2024 *Europhys. Lett.* **146** 27001

[50] Itzykson C and Zuber J B 1980 *Quantum Field Theory (International Series In Pure and Applied Physics)* (McGraw-Hill)

[51] Cugliandolo L F and Kurchan J 1993 *Phys. Rev. Lett.* **71** 173-6

[52] Täuber U C 2014 *Critical Dynamics: A Field Theory Approach to Equilibrium and Non-Equilibrium Scaling Behavior* (Cambridge University Press)

[53] Grabsch A, Rizkallah P, Poncet A, Illien P and Bénichou O 2023 *Phys. Rev. E* **107** 044131

[54] Basu U, Démery V and Gambassi A 2022 *SciPost Phys.* **13** 078
[55] Benois A, Jardat M, Dahirel V, Démery V, Agudo-Canalejo J, Golestanian R and Illien P 2023
*Phys. Rev. E* **108** 054606

[56] Łuczka J 2005 *Chaos* **15** 026107

[57] Lang A, Likos C N, Watzlawek M and Löwen H 2000 *J. Phys.: Condens. Matter* **12** 5087–108

[58] Hull T E and Froese C 1955 *Can. J. Math.* **7** 116–25

[59] Le Bellac M 1991 *Quantum and Statistical Field Theory* (Clarendon)

[60] Debnath L 2015 *Int. J. Appl. Comput. Math.* **2** 223–41

[61] Ditkin V, Prudnikov A and Wishart D 2017 *Operational Calculus in Two Variables and Its Applications (Dover Books on Mathematics)* (Dover Publications) (available at: https://books.google.fr/books?id=tAnDwAAQBAJ)

[62] Wald S, Henkel M and Gambassi A 2021 *J. Stat. Mech.* 103105

[63] Hansen J P and McDonald I R 2013 *Theory of Simple Liquids* (Elsevier) (http://doi.org/10.1016/C2010-0-66723-X)

[64] Dean D S 2018 Stochastic density functional theory (Lecture notes from the Bangalore School on Statistical Physics 2018) (available at: www.icts.res.in/sites/default/files/bssp2018-David-Dean-notes.pdf)

[65] Novikov E A 1965 *Sov. Phys. - JETP* **20** 1290–4 (available at: http://www.jetp.ras.ru/cgi-bin/dn/e_020_05_1290.pdf)

[66] Démery V and Gambassi A 2023 *Phys. Rev. E* **108** 044604

[67] Domazetoski V, Masó-Puigdellosas A, Sandev T, Méndez V, Iomin A and Kocarev L 2020 *Phys. Rev. Res.* **2** 033027

[68] Sandev T, Domazetoski V, Iomin A and Kocarev L 2021 *Mathematics* **9** 221

[69] Greiner A, Strittmatter W and Honerkamp J 1988 *J. Stat. Phys.* **51** 95–108

[70] Klauder J R and Petersen W P 1985 *SIAM J. Numer. Anal.* **22** 1153–66

[71] He Y, Burov S, Metzler R and Barkai E 2008 *Phys. Rev. Lett.* **101** 058101