Field theory of active Brownian particles in potentials

The active Brownian particle (ABP) model exemplifies a wide class of active matter particles. In this work, we demonstrate how this model can be cast into a field theory in both two and three dimensions. Our aim is manifold: we wish both to extract useful features of the system, as well as to build a framework which can be used to study more complex systems involving ABPs, such as those involving interaction. Using the two-dimensional model as a template, we calculate the mean squared displacement exactly, and the one-point density in an external potential perturbatively. We show how the effective diffusion constant appears in the barometric density formula to leading order, and determine the corrections to it. We repeat the calculation in three dimensions, clearly a more challenging setup. Comparing different ways to capture the self-propulsion, we find that its perturbative treatment results in more tractable derivations without loss of exactness, where this is accessible.


I. INTRODUCTION
Active Brownian particles (ABPs) [1] are one of the paradigmatic models of active matter [2]. Such particles move with a constant self-propulsion speed w 0 , in a time-varying direction, represented by a vector subject to rotational diffusion with constant D r . At the same time, the particles are subject to thermal diffusion with constant D t [3,4]. As opposed to Run-and-Tumble (RnT) motion [3], both degrees of freedom, the vectorial velocity w(t) as well as the position r(t), are continuous functions of time t. ABPs approximate so-called Janus particles particularly well [1,2]. Similarly, monotrichous bacteria vary their orientation more smoothly and are much more reliant on rotational diffusion [5]. In contrast, the motion of the much-studied E. coli is closer to RnT.
The motion of ABPs is most directly described by a pair of Langevin equations for their position r(t) as a function of time t, ∂ t r(t) = w(t) + ξ(t) with ξ(t) · ξ(t ′ ) = 2dD t δ(t ′ − t) (1a) where the angular brackets • indicate an ensemble average. Both ξ(t) and ζ ⊥w(t) (t) represent white noise with vanishing mean and correlators as stated. The translational noise ξ(t) is the commonly used Gaussian displacement in d dimensions. The rotational noise ζ ⊥w(t) (t) is confined to the (d − 1)-dimensional tangential plane orthogonal to w(t) at any time [6], so that w(t) stays on the surface of a sphere of radius w 0 at any time. This is guaranteed by w(t) · ζ ⊥w(t) (t) = 0, which requires an Ito interpretation of Eq. (1b). In Section I A we cast the above in Fokker-Planck form, avoiding any such issues of stochastic calculus. The motion of ABPs in two dimensions has been analysed in various ways using classical methods [2][3][4]7], while ABPs in three dimensions have received much less attention. This might be related to the complications arising from the curvature of the surface the rotational diffusion takes place on. The aim of the present work is to develop the mathematical toolkit that will allow us to study ABPs using field theory. Some of us have studied ABPs in two dimensions before [8]. We will now study ABPs in two dimensions from a different perspective, primarily to provide us with a template for the treatment in three dimensions.
Studying ABPs in two and three dimensions through field theory will not only give direct insight about ABPsdeveloping and verifying this framework will also open the door to add interaction. We take a first step in this work by considering external potentials.
There are, in principle, always two ways of describing single-particle dynamics using field theory. Either, one can use a suitable eigensystem in which to expand the fields, or proceed by considering the motion as a perturbation about pure diffusion. The former approach is, in theory, more powerful, as the particle dynamics is fully contained in a set of eigenfunctions, whose eigenvalues determine the "decay" rates. However, a combination of the two approaches is sometimes called for, particularly when there are multiple intertwined degrees of freedom, not all of which can be included in an appropriate eigensystem. The motion of Run-and-Tumble particles in a harmonic potential is an example of such a system [9]. In the present case of ABPs, the choice seems obvious in two dimensions, as, after suitable simplifications, we encounter the Mathieu equation. However, in three dimensions, no such simplifications are available. This necessitates exploring new methods, in particular determining the eigenfunctions themselves perturbatively, or treating self-propulsion as a perturbation about pure diffusion. This latter method is first tried out in the easier case of two dimensions.
The complexity of the methods increases further when we allow for an external potential, which must be treated perturbatively. As a test bed for the methods, we return to ABPs in two dimensions [8], before moving on to three dimensions. In both cases, we determine the mean squared displacement (MSD) in closed form, as well as the stationary density in the presence of an external potential. These two observables are related: the MSD naturally gives rise to the notion of an effective diffusion constant, which may or may not feature as the effective temperature in a Boltzmann-like distribution of particles in the stationary state.
In the following, we introduce the ABP model, first through a Fokker-Planck description, and subsequently in the form of an action. In Section II, we proceed to characterise the MSD and stationary distribution in two dimensions. We do this first using Mathieu functions, and then perturbatively. These methods and results serve as the template for the treatment of ABPs in three dimensions in Section III.

A. Active Brownian Particles
Active Brownian Particles (ABP) are particles that move by self-propulsion by a constant speed w 0 as as well as due to thermal noise characterised by translational diffusion constant D t . The direction of the self-propulsion is itself subject to diffusion. In the following, we will generally denote this direction by Ω, and more specifically in the plane the polar angle by ϕ ∈ [0, 2π) and in spherical coordinates the azimuthal angle by ϕ ∈ [0, 2π) and the polar angle by θ ∈ [0, π). The direction of the self-propulsion is then in two and three dimensions respectively. The rotational diffusion is parameterised by the constant D r , specifically by a Fokker-Planck equation of the form ∂ t P (Ω, t) = D r ∇ 2 Ω P (Ω, t) with in two and three dimensions respectively. The Fokker-Planck equation of active Brownian particles is ∂ t P (r, Ω, t) = LP (r, Ω, t) with L = D t ∇ 2 r − w(Ω) · ∇ r + ∇ r · ∇ r Υ(r) + D r ∇ 2 Ω where P (r, Ω, t) is the probability to find a particle after some initialization at time t at position r with the director pointing to Ω and ∇ r the usual spatial gradient as opposed to ∇ Ω acting on the director, Eq. (3). The self-propulsion is implemented by the term w · ∇ r , which is the only term mixing the degrees of freedom Ω and r. The effect of an additional external potential Υ(r) which acts equally on particles without orientation is implemented by ∇ r · (∇ r Υ), which is structurally identical to that of a drift w and reduces to that term when ∇ r Υ is constant in space. We write Υ ′ for ∇ r Υ in the following. When the operator ∇ r · Υ ′ acts on P in Eq. (4), only the leftmost spatial gradient ∇ r acts on the product of Υ ′ and P , whereas ∇ r Υ results solely in Υ ′ . For ease of notation below, we further introduce From the single-particle Fokker-Planck Eq. (4) the multiple, non-interacting particle action A[χ, χ] = d d r dΩ dtχ(r, Ω, t) − ∂ t + L − r χ(r, Ω, t) follows immediately [10], with an additional mass r ↓ 0 as a regulator. At this stage, we leave the dimension d ∈ {2, 3} and the integral dΩ unspecified. Expectations of the annihilator field χ(r, Ω, t) and the creator fieldχ(r, Ω, t) can be calculated in a path integral [10] • = Dχ Dχ e A[χ,χ] • = •e −AP [χ,χ] which allows for a perturbative treatment by splitting the action into a harmonic part A 0 , whose path-integral is readily taken, and a perturbative part A P , whose exponential is expanded, where we have introduced • 0 as an average on the basis of the harmonic part of the action only. Which part of the action Eq. (6) is treated perturbatively depends crucially on the eigensystem the fields χ andχ are written in, as it determines which part of the operator becomes diagonal. Because the external potential is arbitrary, it cannot be diagonalised in its present generality. We may thus consider two schemes: In the first, the fields are written in eigenfunctions of L 1 , Eq. (5), with the harmonic and perturbative part of the action A 01 = d d r dΩ dtχ(r, Ω, t) − ∂ t + L 1 − r χ(r, Ω, t) (9a) A P 1 = d d r dΩ dtχ(r, Ω, t)∇ r · ∇ r Υ(r) χ(r, Ω, t) , which is more efficiently dealt with using Gauss' theorem in A P 1 to make ∇ r act onχ rather than Υ ′ χ. The benefit of including the self-propulsion in the harmonic part of the action A 01 , Eq. (9a), via L 1 , Eq. (5a), at the expense of having to find suitable eigenfunctions, is to have only one perturbative vertex to worry about, namely that representing the external potential. The major disadvantage is that the diagonalising effect of the eigenfunctions on L 1 relies on momentum conservation, which is notably broken in the presence of an external potential, so that the external potential becomes significantly more difficult to handle. A similar problem would be encountered if we were to allow for pair interaction, in which case the action contains terms that are not bilinear, which result in difficult integrals involving products of more than two eigenfunctions.
The eigenfunctions of the angular part of L 1 in two spatial dimensions, after some manipulation, are the π-periodic Mathieu functions [11], as discussed further below. Compared to other orthonormal systems, they are far less well studied [12]. In three dimensions, we shall call them "three-dimensional Mathieu functions", but as far as we are aware, they have not been introduced in the literature of orthonormal functions in their own right. We have characterised them to some extent in App. A using similar tools as outlined in [11]. Of course, the complete framework of Sturm-Liouville theory provides broad theoretical backing. Once their existence is established, the eigenfunctions need to be determined only to the extent that the observables require it.
In the second scheme, the fields are written in terms of eigenfunctions of the simpler L 2 in Eq. (5), with the harmonic and perturbative part of the action This results in simpler expressions overall, but suffers from the significant disadvantage of an additional vertex mediating the self-propulsion. One may expect a significantly enlarged number of diagrams to be considered in any interesting observable, but this is not the case for the observables considered here. We proceed by determining the MSD and the stationary density for ABPs in two dimensions. In the course, we will discuss the details of the action, the diagrammatics and the eigenfunctions. The starting point of the present derivation is the operator L 1 in Eq. (5a) and the action Eq. (9a), restated here with the explicit parameterisation of the director for convenience. The key difficulty here is to find eigenfunctions of the operator −w · ∇ r + ∇ 2 Ω . This is greatly helped by measuring the azimuthal angle relative to the reciprocal vector k, as we will discuss in the following. We first introduce where u ℓ (γ, q(k)) and u ℓ (γ, q(k)) are suitably normalised π-periodic Mathieu functions [11,12], depending on the dimensionless, purely imaginary parameter which is a function of the (absolute) magnitude k = |k| = | − k| of the k vector. The function σ(k), hitherto undetermined, allows a crucial simplification of the action further below. For the time being it suffices to state orthogonality even in the presence of the shift by σ, as the δ¯(k + k ′ ) forces this shift to be identical, σ(k) = σ(−k ′ ), which can therefore be absorbed into the dummy variable and subsequently into a change of the integration limits, from where it disappears, because of the periodicity of the integrand. Using Eq. (12) in the action Eq. (11) produces where α is the polar angle of k = k(cos α, sin α) T , so that k · w(ϕ) = kw(cos α cos ϕ + sin α sin ϕ) = kw cos(ϕ − α). Choosing σ(k) = α allows a change of variables and a change of integration limits, so that the last integral in Eq. (15) reads where in the last equality we have used that the Mathieu functions obey and the orthonormality Eq. (14). With this simplification the action Eq. (11) becomes diagonal and the (bare) propagator can immediately be read off where we have introduced G ℓ (k, ω) as a shorthand, as well as the diagrammatic representation to be used in Section II A 2.

Mean squared displacement
The mean-squared displacement (MSD) r 2 (t) of a two-dimensional ABP (2DABP) can be calculated on the basis of a few algebraic properties of the u ℓ (γ, q). The probability of finding a 2DABP with director ϕ at position r at time t, after being placed at r 0 with director ϕ 0 at time t 0 is χ(r, ϕ, t)χ(r 0 , ϕ 0 , t 0 ) , which is easily determined from Eq. (19) as up to an inverse of Fourier-transform of k. Making use of the Fourier transform, the MSD of a 2DABP irrespective of the final director is Having integrated over ϕ, it is obvious that the MSD cannot depend on the initial orientation of the director, ϕ 0 , as changing it, amounts to a rotation of the isotropic plane. This is indeed trivial to show when w 0 = 0, in which case q(k) = 0 and the Mathieu functions degenerate into trigonometric functions. However, demonstrating this when w 0 = 0 is far from trivial. We proceed by rewriting the MSD as an average dϕ 0 (2π) −1 of Eq. (21) over the initial angle ϕ 0 , so that both Mathieu functions, u ℓ andũ ℓ reduce to a single coefficient via which vanishes in fact for all odd indices ℓ. We thus arrive at The u ℓ (γ, q) are in fact the even π-periodic Mathieu functions for even ℓ and odd π-periodic Mathieu functions for odd ℓ [11]. The u ℓ (γ, q) differ from the u ℓ (γ, q) only by a factor 1/π [8], which is a convenient way to meet the orthogonality relation Eq. (14) while maintaining that the non-tilde u ℓ (γ, q) are standard Mathieu functions. At this stage, all we need to know about the Mathieu functions is [11] A 0,0 (q) = 1 for j ∈ N 0 , which means that only two of the A ℓ,0 in the sum Eq. (23) contribute, namely ℓ = 0 and ℓ = 2, as all other A ℓ,0 vanish at q(k = 0) = 0 even when differentiated twice. With Eq. (25) the MSD becomes [8] identical to that of Run-and-Tumble particles, Eq. (49) of [8], if the tumble rate α is D r . In asymptotically large t, the MSD of a 2DABP is equivalent to that of conventional diffusion with effective diffusion constant and RnT motion with tumbling rate α = D r [3,4,8].
We will use the above as a template for the derivation of the MSD of 3DABP. There, we will make use again of the isotropy of the initial condition ϕ 0 and further of the isotropy of the displacement in any of the three spatial directions. In fact, picking a "preferred direction" can greatly reduce the difficulty of the present calculation.

External Potential
Next, we place the 2DABP in a (two-dimensional) potential Υ(r). Without self-propulsion, the density ρ 0 (r) at stationarity is given by the barometric formula ρ 0 (r) ∝ exp (−Υ(r)/D t ). As a sanity check for the present fieldtheoretic framework, we derive it in App. C. In the present section we want to determine the effect of the self-propulsion, which we expect is partially covered by D t in the Boltzmann distribution being replaced by D 2D eff , Eq. (27). As the external potential is treated perturbatively, stationarity even at vanishing potential requires a finite volume, which we assume is a periodic square, i.e. a torus, with linear extent L. This requires a small adjustment of the Fourier representation introduced in Eq. (12), which now reads where L −2 n∈Z 2 has replaced the integration over k in Eq. (12), with the couple n = (n x , n y ) T running over the entire Z 2 and k n = nL/(2π). The bare propagator Eq. (19) is adjusted simply by replacing δ¯(k+k 0 ) by, say, L 2 δ n+p,0 . With the introduction of Fourier sums, we implement periodic boundary conditions, effectively implementing periodic observables and, below, a periodic potential.
As far as the Mathieu functions are concerned this adjustment to a finite domain has no significant consequences. Due to the Mathieu functions, however, the perturbative part of the action A P 1 Eq. (9b) gets significantly more difficult, as the functions are orthogonal only if the k-modes match, where each n 1 , n 2 and n 3 runs over all Z 2 in the sum and the Kronecker δ-function δ n1+n2+n3,0 enforces n 1 +n 2 +n 3 = 0. Similarly, ℓ 1 and ℓ 3 run over all non-negative integers indexing the Mathieu functions. The potential Υ(r) now enters via its modes Υ n , more explictly rendering the potential effectively periodic. Of course, a non-trivial Υ(r) spoils translational invariance and correspondingly Υ n provides a source of momentum in Eq. (29b), so that k n3 is generally not equal to −k n1 , with the sole exception of n 2 = 0. As a result, generally σ(k n3 ) = σ(−k n1 ) and q(|k n3 |) = q(|k n1 |), Eq. (13). The integral however, is generally diagonal in ℓ 1 , ℓ 3 only when k n1 = −k n3 . It is this projection, Eq. (31), that is at the heart of the complications to come. It is notably absent in the case of purely diffusive particles, App. C 4. With the help of Eq. (31), the perturbative action Eq. (29b) can be written as The perturbative part of the action, Eq. (29b), may diagrammatically be written as where a dash across any line serves as a reminder that the vertex carries a factor of k carried by the leg. Unless there are further singularities in k, any such line vanishes at k = 0. The two factors of k are multiplied in an inner product. The dangling bauble in Eq. (33) represents the external potential, which is a source of momentum. We will determine in the following the stationary density at position r, starting with a single particle placed at r 0 with orientation ϕ 0 at time t 0 . We will determine this density to first order in the potential and second order in q = 2iw 0 k/D r , Eq. (13), starting from the diagrammatic representation in inverse space, where the first diagram represents the bare propagator as introduced in Eq. (19) and the second diagram the contribution to first order due to the external potential Eq. (33). The limit t 0 → −∞ needs to be taken after an inverse Fourier transform in ω and ω 0 . The structure of the propagator Eq. (19), essentially of the form (−iω + λ + r) −1 with some non-negative λ ≥ 0 dependent on the spatial and angular modes, means that every such Fourier transform results in an expression proportional to exp (−(λ + r)(t − t 0 )). Removing the mass by r ↓ 0 and taking t 0 → −∞ results in the expression vanishing, except when ℜ(λ) ≤ 0 [13], which is the case, by inspection of Eq. (19), only when k = 0 and ℜλ ℓ (q(k)) ≤ 0, which requires ℓ = 0, as Eq. (36) is easily verified as the periodic solutions of the Mathieu Eq. (17) are trigonometric functions. Any propagator carrying a factor of k therefore vanishes under this operation. The only pole of the propagator to be considered is the one at ω = −ir. Because only the rightmost, incoming leg in any of the diagrams is undashed, taking the inverse Fourier transform followed by the limit t 0 → −∞ thus amounts to replacing the incoming legs by δ n,0 and every occurance of ω in the diagram by 0. This is equivalent to an amputation of the incoming leg [13]. At this point, the inverse Fourier transforms from k to r and the transforms of azimuthal degree of freedom from ℓ to ϕ are still to be taken. We may summarise by writing Using Eqs. (19) and (32), the last diagram evaluates to The projection ∆ ℓ,ℓ0 (k n , k p ) needs to be known only for ℓ 0 = 0 and k p = 0, and because u ℓ=0 (γ, q = 0) = 1/ √ 2, Eqs. (22) and (31) give To carry out the integral over ϕ in Eq. (34), we have to transform Eq. (37) from ℓ, ℓ 0 to ϕ, ϕ 0 using Eq. (12) first, In the first term, k n = k p = 0, so that the product of the two Mathieu functions degenerates to 1/(2π), cancelling with the integration over ϕ still to be performed. In the second term, integration over ϕ can again be done with the help of Eq. (22), while the final Mathieu function u 0 (41) Since A ℓ,0 (q), Eq. (25), vanish for odd ℓ and otherwise behave to leading order like q ℓ/2 , Eq. (25), expanding these coefficients to second order means to retain only ℓ = 0, 2 in the sum, so that with explicit G ℓ (k n , 0), Eq. (19), we arrive at where q 2 = −4w 2 0 k 2 n D −2 r and k n = |k n |. The limit r ↓ 0 is yet to be taken and . . . on the right of Eq. (42) refers to the higher orders in the external potential. Expanding the eigenvalues λ ℓ beyond Eq. (25), suggests that k n · k n+p δ p,0 can be cancelled with k 2 n in the denominator when ℓ = 0. However, for n = 0 the denominator does not vanish as r > 0. Performing the inverse Fourier transform to return real space gives after r ↓ 0, where n = 0 is explicitly omitted from the sum and the first term in the curly brackets shows the expected D 2D eff = D t + w 2 0 /(2D r ), Eq. (27), amended, however, by a correction of order q 4 /k 2 n . To identify the lowest order correction in k n and weak potentials, we expand further than Eq. (43), in the first (of three terms, in the curly bracket) correction term in the sum in Eq. (44), which gives by dropping the last term in Eq. (44), as it is O k 4 n , and not expanding the second term any further, as it is O k 2 n as it stands. The expression above it to be compared to the first order term in the barometric formula, Eq. (C20), which produces only −Υ n /D t as the first-order summand when setting w 0 = 0. The corrections beyond that are generically due to the stochastic self-propulsion and amount to more than replacing D t by D 2D eff . Eq. (46) completes the derivation in the present section. In the next section, we will re-derive MSD r 2 (t) and density ρ 0 (r) of active Brownian particles in two dimensions at stationarity, in an attempt to recover Eqs. (26) and (44) or (46) in a perturbation theory in small w 0 and small potential.

B. Perturbation in small w0
The starting point of the present derivation is the operator L 2 in Eq. (5b) and the harmonic part of the action Eq. (10a), which becomes diagonal with Compared to Eq. (12), in Eq. (48) the Mathieu functions have been replaced by exponentials, which are orthonormal irrespective of k. The bare propagator of Eq. (47) can immediately be read off, The perturbative part of the action Eq. (10b) includes the self-propulsion, but its precise form depends on whether the system is finite or infinite and thus left to be stated in the sections ahead.

Mean squared displacement
It is instructive to derive the MSD with harmonic action Eq. (47), orthogonal system Eq. (48) and perturbative action Eq. (10b), rewritten as in the absence of an external potential. To ease notation, we may introduce the tumbling vertex The expressionik x cos ϕ +ik y sin ϕ is due to w · ∇ r of the self-propulsion hitting exp (ikr), Eq. (2). For the diagrammatics we introduce the self-propulsion vertex so that the full propagator reads using the notation of Eq. (50). A marginalisation over ϕ amounts to evaluating Eq. (54) at ℓ = 0, Eq. (48). The MSD can be determined by double differentiation with respect to k as done in Eqs. (21) and (26). This calculation is simplified by calculating the mean squared displacement in x only, i.e. taking ∂ 2 kx rather than ∂ 2 kx + ∂ 2 ky . This "trick", however also requires uniform initialisation in ϕ 0 to avoid any bias. In summary where the integral over k 0 is the inverse Fourier transform with x 0 = 0 and is trivially taken by using δ¯(k + k 0 ) due to translational invariance. The two integrals over ϕ and ϕ 0 in Eq. (55) amount to setting ℓ = 0 and ℓ 0 = 0, with all factors of 2π cancelling. Of the three terms in Eq. (54), the second one, which is the first order correction in the self-propulsion carrying a single factor of the tumbling vertex W ℓ,ℓ0 (k), Eq. (52a), is odd in k x . Differentiating twice with respect to k x and evaluating at k x = k y = 0 will thus make it vanish. The third term is quadratic in k x in the numerator and not singular in k x in the denominator. The only way for the third term not to vanish at k x = k y = 0 is when W ℓ,ℓ1 (k)W ℓ1,ℓ0 (k) is differentiated twice with respect to k x , No higher order terms in w 0 can contribute, as they all carry higher powers of k x as a pre-factor. We thus arrive at where the term 2/(−iω + D r + r) originates from H ℓ1 (k, ω) at ℓ 1 = −1, 1 according to Eq. (56). Taking the inverse Fourier transform from ω to t− t 0 then indeed reproduces Eq. (26) in the limit r ↓ 0 with the repeated poles producing the desired factors of t − t 0 . This confirms that even when the self-propulsion is dealt with perturbatively, observables can be derived in closed form. This is not a triviality, as w 2 0 t/D t is a dimensionless quantity and thus might enter to all orders. We proceed in the present framework by allowing for an external potential.

External Potential
In this section, we aim to re-derive the density ρ 0 (r) of active Brownian Particles in an external potential, as initially determined in Section II A 2. As the self-propulsion is a perturbation, the full perturbative action Eqs. (10b), (51) reads using the notation k n1 = (k n1x , k n1y ) T . The action Eq. (58) is based on discrete Fourier modes similar to Eq. (28) as the system's volume is finite. It deviates from Eq. (28) but using Fourier modes for the director as in Eq. (48) Compared to Eq. (29b), Eq. (58) has the great benefit of the diagonality of the angular modes in the potential term.
To calculate the particle density ρ 0 (r) in the stationary state, Eq. (34), we expand the full propagator diagrammatically, As in Section II A 2, taking the limit t 0 → −∞ amounts to an amputation upon which any diagram vanishes that carries a dashed propagator before (i.e. to the right of) a source of momentum. As a result, only the five leftmost diagrams in Eq. (60) contribute at stationarity. In the following, we will restrict ourselves to diagrams with at most one bauble, i.e. to first order in the potential. Further, we will incorporate the self-propulsion to second order, leaving us with the four leftmost diagrams in Eqs. (60a)-(60d). As seen in Section II A 2, Eq. (60a) contributes to the density merely with 1/L 2 . The other three diagrams require more detailed attention. Firstly, Eq. (60b) gives using Eq. (50) for H ℓ (k n , ω). Upon taking the inverse Fourier transform in ω and ω 0 , the limit t 0 → −∞ can be taken, replacing the incoming leg H ℓ0 (k p , ω 0 ) effectively by δ p,0 δ ℓ0,0 . After transforming from ℓ, ℓ 0 to ϕ, ϕ 0 and taking the integral over ϕ all factors of 2π have cancelled, leaving of Eq. (61) as H 0 (k n , 0) = (D t k 2 n + r) −1 , whose factor k −2 n cancels with k n · k n unless k n = 0, in which case the expression vanishes. Eq. (62) is the lowest order correction in both Υ and w 0 to the density.
Of the remaining two diagrams, the leftmost in Eqs. (60c) and (60d), only the latter contributes once t 0 → −∞ forces ℓ 0 = 0 and dϕ forces ℓ = 0, because the single tumbling vertex, Eq. (52a), requires the angular modes to be offset by 1. The only diagram left to consider is thus After taking t 0 → −∞ and dϕ , the only two contributions from the sum are ℓ 1 = −1, 1, so that from Eq. (52a), so that W 0,1 (k)W 1,0 (k) = W 0,−1 (k)W −1,0 (k) = −w 2 0 k 2 /4. After taking the limit r ↓ 0 in Eq. (64) the factor k 4 n cancels with the same factor in (D t k 2 n + r) −2 whenever n = 0, so that finally to be compared to the purely drift-diffusive Eq. (C20) and to Eqs. (44) or (46), derived using Mathieu functions, rather than the perturbative approach to self-propulsion here. By inspection of the two expressions, Eqs. (44) and (66b), the two are equivalent provided which can indeed be shown by expanding all terms on the right hand side to order w 2 0 only, absorbing all other terms into O w 4 0 , for example This concludes the discussion of ABPs in two dimensions. In the following sections, we will calculate MSD and density in an external potential for ABPs in three dimensions, using the above as a template.

III. ACTIVE BROWNIAN PARTICLES IN 3 DIMENSIONS
In three spatial dimensions, the diffusion of the director takes place on a curved manifold, which renders the Laplacian rather unwieldy, Eq. (3). The spherical harmonics provide a suitable eigensystem for this operator, but the self-propulsion term spoils the diagonalisation. As in the two-dimensional case, Section II, there are two possible ways ahead: Either find suitable special functions or deal with the self-propulsion term perturbatively. In the following subsection, we focus on the former.

A. "Three-dimensional Mathieu functions"
The harmonic part of the action is Eq. (9) with Eqs. (2) and (3) in Eq. (5a), more explicitly using the notation r = (x, y, z) T . The eigensystem that we will use now is modelled along the two-dimensional Mathieu functions, Eq. (12). Specifically, it is orthonormal, similar to Eq. (14) and q(|k|) of the two-dimensional case, Eq. (13), replaced by and two instead of one arbitrary functions σ(k) and τ (k) cf. Eq. (14). The benefit of those, however, is rather limited compared to the two-dimensional case, where k · w was rewritten as kw 0 cos(ϕ − α). We therefore choose σ ≡ τ ≡ 0 and demand These eigenfunctions are characterised in some detail in App. A. For now, it suffices to know that they exist and that the eigenvalues λ m ℓ (q) are indeed discrete and indexed in both ℓ and m. Having relegated the details of the eigensystem to the appendix, the bare propagator can now be determined without much ado, and we proceed with the MSD.

Mean squared displacement
The MSD derives from the propagator Eq. (75a) via a double derivative. As in the two-dimensional case, the derivation simplifies considerably with uniform initialisation. Further, because the eigenequation (74) is particularly simple whenever k x = k y = 0, we may write the MSD in terms of derivatives with respect to k z only, as derived in more detail in Eq. (A27). These coefficients we need to determine at most up to order q 2 in order to express r 2 (t) in closed form. Using the notation q =iw 0 k z /D r in q(k z e z ) = qe z , Eq. (76) produces with Eq. (77) from the expansion of A 0,0 0,0 (qe z ), Eq. (A22a), and A 0,0 1,0 (qe z ), Eq. (A23a), with A 0,0 ℓ,0 (qe z ) ∈ O q 2 for ℓ ≥ 2, Eq. (A24). The two eigenvalues λ 0 0 (qe z ) = −q 2 /6 + . . . and λ 0 1 (qe z ) = 2 + . . ., Eqs. (A18) and (A20), need to be known only to leading order in q. Performing the derivative and re-arranging terms then produces the final result structurally similar to the result in two dimensions, Eq. (26), and indeed identical to the MSD of Run-and-Tumble particles in three dimensions, Eq. (49) of [8], if the tumble rate α is replaced by 2D r [3,4]. The effective diffusion constant in three dimensions is to be compared to Eq. (27).

External Potential
Allowing for an external potential and confining the particles to a finite space means that the fields are now written as similar to Eqs. (28) and (70) and any expectation needs to be taken with the perturbative action similar to Eq. (31). The resulting diagrammatics of the full propagator to leading order in the potential is as in Eq. (35) but with the index ℓ replaced by the couple ℓ, m. The arguments that follow to extract the stationary density are similar to those in Section II A 2. Yet, to be able to draw on the results in App. A, all k-vectors featuring in the three-dimensional Mathieu functions need to be parallel to e z . Because amputated diagrams draw all kdependence from spatial variation of the potential, we thus demand Υ(r) = Υ(z), i.e. that the potential depends only on the z-component of r = (x, y, z) T . Its Fourier transform can therefore be written in the form where n = (n x , n y , n z ) T . The notation of Eq. (85) is slightly ambigious, as Υ nz on the right is a Fourier transform in the z-direction only, while Υ n on the left is a Fourier transform in all space, yet the alternatives seem to obfuscate even more. Just like in Section II A 2, the stationary limit enforces k p = 0, as well as ℓ 0 = 0, as otherwise λ m ℓ (q(0)) = 0, Eq. (A5b), in the propagator Eq. (75a).

B. Perturbation in small w0
Along the lines of Section II B, we return now to calculate the MSD and the stationary density while treating the self-propulsion as a perturbation over diffusion, so that we can use the eigensystem denote the well-known spherical harmonics with P m ℓ denotes the standard associated Legendre polynomials which carry a factor (−1) m in their definition, so that, for example Y 1 1 (θ, ϕ) = − 1 and obey the eigenvalue equation The harmonic part of the action Eq. (10a) now reads so that the bare propagators are differing from Eq. (50) only by the additional factor δ m,m0 and by the eigenvalue of the spherical Laplacian being ℓ(ℓ + 1) rather than ℓ 2 . In the absence of an external potential the perturbative part of the action reads similar to Eq. (51) with tumbling vertex In principle, any product like sin θ cos ϕY m ℓ (θ, ϕ) can be written in terms of a sum of spherical harmonics (see App. B), where c L,M ℓ ′ ,m ′ ,ℓ,m are the Clebsch-Gordan coefficients, further detailed [15] in App. B. The coupling matrix W m ′ ,m ℓ ′ ,ℓ (k) can therefore be written as the simplest among the three terms is with the other two stated in App. B. With this framework in place, we proceed to calculate the MSD and the density at stationarity.

Mean squared displacement
We follow Section II B 1 to calculate the MSD of active Brownian particles in three dimensions, which is corresponding to Eq. (55). With the help of the diagrammatics Eq. (54) and the reasoning before Eq. (56), this can essentially be done by inspection. The derivative of the bare propagator will contribute −3(−2D t t) to the MSD. The only non-trivial contribution is due to , which has a very simple structure because the integration over θ and ϕ, forces ℓ = 0 and m = 0, given that Y m ℓ are orthogonal and Y 0 0 = 1/ √ 4π and so integrals over Y m ℓ vanish unless ℓ = m = 0. The same reasoning applies to ℓ 0 = m 0 = 0. The summation that appears in the last term of Eq. (54) therefore becomes because only ℓ 1 = 1 contributes, resulting in a factor ℓ 1 (ℓ 1 + 1) = 2 in front of D r in the propagator. Differentiating twice with respect to k z and taking the inverse Fourier transform of the expression with the bare propagators attached, finally reproduces Eq. (79) exactly.

External potential
Calculating the stationary density ρ 0 (r) in the presence of an external potential perturbatively in the self-propulsion and the potential follows the pattern outlined in Section II B 2. Because the volume is finite, the representation we choose for the fields is rather than Eq. (91), with the sums over n ∈ Z 3 replacing the integrals over k.
Given the orthogonality of the spherical harmonics, the perturbative part of the action now reads corresponding to Eq. (58) in two dimensions, but with the tumbling vertex W m ′ ,m ℓ ′ ,ℓ (k) given by Eq. (98). In principle, Υ n is arbitrary, but in order to stay with the simple form of W m,m ′ ℓ,ℓ ′ (k z e z ), Eq. (101b), we restrict Υ(r) to vary only with z, making again use of the notation Eq. (85). The sole contribution that needs our attention is due to the single term Eq. (63), which draws on Eq. (103).

IV. DISCUSSION AND CONCLUSION
There are two main perspectives on the present results, first about the concrete observables calculated and second about the formalism. a. Observables: In this work we calculated the mean squared displacement and the stationary density of ABPs in two and three dimensions using two distinct methods. The MSD calculations were performed primarily for methodological purposes, although the full time-dependence may be a new result. [2,3,7,8].
In two dimensions, we determined the exact MSD of ABPs using Mathieu functions with full time dependence in Eq. (26) as well as using the perturbative approach leading to Eq. (57). In three dimensions, similarly to the two dimensional case, we used "three-dimensional Mathieu functions" for computation of the exact full time-dependent MSD in Eq. (79) and spherical harmonics in Eq. (104). The effective diffusion constants for both, the two and three dimensional case in Eq. (27) and Eq. (80) respectively, are identical to the Run-and-Tumble motion case using appropriate conversion [3,8].
Generalizing to d > 3 dimensions using the perturbative approach and so-called hyperspherical harmonics [16,17], the expression for the MSD of ABPs is [17] with the details to be found in [18]. The effective diffusion constant is consistent with [3].
The MSD was calculated irrespective of the final orientation of the ABP and integrated about the initial angle without loss of generality. This leads to a significant reduction in the number of terms in the final MSD expression as discussed after Eq. (21) in two dimensions. The same method was used in the three dimensional case with "three-dimensional Mathieu functions" in Section III A 1. For the self-propulsion perturbative case, it is a matter of convenience to determine the mean squared displacement in only one spatial direction, which however requires uniform initialisation in order to avoid bias, as discussed in two dimensions after Eq. (59). The same approach was taken in three dimensions in Section III B 1.
The calculation of the stationary density in the presence of the external potential is, to our knowledge, a new result with the present generality. The stationary density in two dimensions was computed perturbatively in presence of small potentials in Eqs. (44), (46) and (66) to first order in the potential. These expressions differ only in that they are expanded in different parameters, i.e. the parameter q of the Mathieu function in (44), the spatial momentum k n in (46) and k, w 0 in (66). We further assumed a finite volume L d , to guarantee a stationary state in the absence of a potential, similar to drift-diffusive case discussed in Section C. As a result of using Fourier sums, the potentials and the resulting observables are periodic. The two-dimensional MSD is in agreement with experimental results [1].
In three dimensions, the stationary density calculation with the "three dimensional Mathieu-functions" in the external potential was simplified by allowing the potential to vary in only one direction, Eq. (85), App. A. In two dimensions, this assumption was avoided by the re-orientation of the coordinate system in Section II A. In the perturbative approach, Section III B, the above-mentioned assumption helped avoid having to compute general Clebsch-Gordan coefficients in Eq. (98), but determined nevertheless in App. B.
The stationary density in three dimensions in Eq. (89) is calculated to first order in the potential and expanded to second order in q. The perturbative calculation in the self-propulsion results in Eq. (107). The restriction on the potential in three dimensions is more easily lifted in the perturbation theory using spherical harmonics, since the Clebsch-Gordan coefficients can be derived without much ado. On the other hand, characterising the "threedimensional Mathieu functions" requires at the very least a careful analysis of the completeness of the eigensystem, more easily done for q e z .
Despite the restrictions on the potential, the perturbation theory covers a wide class of potentials. As is analysed in detail for drift-diffusive systems in Section C, the resummed pertubation theory recovers known closed form expressions, including a confining potential. The most interesting case of "active sedimentation" in a gravitational potential in a finite vessel can be realised within the restrictions mentioned above by long-stretched, periodic ratchet in the z-direction with period L, with particle mass m, gravitational acceleration g and mobility µ. With r confined to [0, L) d by periodicity, r·e z = −ǫL with ǫ ∈ [0, 1) is mapped to Υ((1 − ǫ)Le z ) by periodicity in Eq. (30), so that the potential "energy" Eq. (111) jumps to the maximum right after reaching the minimum for decreasing t, implementing a sharp potential wall at z = 0, the bottom of the vessel. The "barometric formulas" in Eqs. (46) and (90) show the corrections to the density that are generically due to the activity. Those cannot possibly affect the 0-mode, which is explicitly excluded from the summation. To leading order the density of ABPs behaves like that of an equilibrium system with no self-propulsion, w 0 = 0 and the translational diffusion constant D t , Eq. (C20), replaced by the effective diffusion constant D eff , Eq. (27) and (80). Terms beyond that are generically due to the non-equilibrium nature of active matter. The lowest order correction is quadratic in k n , namely w 2 0 k 2 n /(2D 2 r ) in two dimensions and w 2 0 k 2 n /(12D 2 r ) in three dimensions. It might be interesting to investigate how this can be measured in slowly varying potentials where contributions in large n are suppressed because Υ n vanishes there.
b. Methodology We have introduced all ingredients to set up a field theory of ABPs, in particular: the action, the propagators, the projections of the special functions such as ∆ ℓ1,ℓ3 (−k n1 , −k n3 ), Eqs. (31) and (83)), and the tumble vertices (such as W ℓ ′ ,ℓ (k), Eqs. (52a) and (98). As well as providing the concrete mathematical framework, our derivations above highlight the advantages and disadvantages of using special functions rather than a perturbation theory. We will briefly summarise these here.
Using special functions, i.e. the Mathieu functions, Section II A and Section III A, has the benefit that the bare propagator already contains the self-propulsion. As long as space r and its boundary conditions can be implemented through a Fourier-transform, this is the most efficient way to deal with free problems, wrapping the self-propulsion in the properties of certain special functions. The resulting bare propagator is the exact and complete characterisation of the free stochastic particle movement. When the special functions are well-known, literature provides the relevant properties of functions and eigenvalues. When they are not well-known, they can be determined perturbatively in (hyper-)spherical harmonics, effectively performing the perturbation theory that would otherwise be done explicitly at the level of diagrams. In return, the wealth of literature on Sturm-Liouville problems is readily applicable, as touched on in Eq. (A8).
We were able to avoid having to determine the "three-dimensional Mathieu functions" in greater detail by using a preferred direction, which generally requires integration of final and/or initial orientation. Thanks to the arbitrary function σ(k) in Eq. (12), this was avoided two dimensions, where we could use the well-known Mathieu functions. This may be related to the metric of a d-dimensional sphere to be curved only in d > 2. In three dimensions, despite having two arbitrary functions σ(k) and τ (k) in Eq. (70), an elegant simplification of the eigenproblem does not seem to be available.
Using a perturbation theory, i.e. expanding diagrams in w 0 , Section II B and Section III B, means that even the free particle has to be treated explicitly perturbatively. Although the full propagator now contains infinitely many terms, this does necessarily translate to infinitely many terms for every observable, as seen, for example, in the fully time-dependent MSD calculated exactly using perturbation theory in two and three dimensions, Section II B 1 and Section III B 1 respectively. The big advantage of using a perturbation theory is that orthogonality is not spoiled by a lack of momentum conservation, as is caused by an external potential or interaction. Such a lack of momentum conservation is what necessitates the introduction of the projection ∆ ℓ1,ℓ3 (−k n1 , −k n3 ) in the case of special functions, Eq. (32) and similarly (82), but this does not happen when using a perturbation theory, Eqs. (58) and (106).
It is striking that the perturbation theory in weak potentials recovers the barometric formula for a binding potential (App. C). The calculation of such stationary features is facilitated by the amputation mechanism as discussed in Section II A 2 and similarly used in [8].
c. Outlook: We are currently working on a field theory of interacting ABPs displaying features of the Vicsek Model [19,20]. As the upper critical dimension there is d c = 4, we aim to generalise the spherical harmonics to higher dimensions. This would also allow the use of dimensional regularisation in a renormalised field theory. While some characteristics of such hyperspherical functions follow a systematic pattern, allowing the spatial dimension to vary continuously like d = 4 − ǫ represents a significant challenge.
In the following, we want to characterise the "three-dimensional Mathieu functions", which we call that only because the eigenfunctions of the equations below in two dimensions are the well-known Mathieu-functions [11,12].
For q = 0 vanishing, Eq. (A1) is the eigen-equation of the spherical harmonics [12]. We therefore choose to express the three-dimensional Mathieu functions in spherical harmonics Y k j (θ, ϕ), with the factor of √ 4π solely to ease notation and P m ℓ (z) denoting the associated Legendre polynomials. For q = 0, where we assume that the coefficients A m,k ℓ,j (q) and eigenvalues λ m ℓ (q) can be expanded in small q. Explicitly, for with the asterisk denoting the complex conjugate and indeed We focus henceforth on those q that have no projection in the x and y directions, q = (0, 0, q) T = qe z , which simplifies Eq. (A1) dramatically, rendering it similar to the expression in two dimensions, Eq. (17). Only the θdependence is affected by the presence of q e z , so that A m,k ℓ,j (q) ∝ δ m,k is diagonal in m, k for such q, resulting in a simple structure and indexing of the azimuthal dependence, i.e. on ϕ, Because we will integrate u m ℓ (θ, ϕ, qe z ) over ϕ, we can further narrow down the scope of this section to m = 0, which simplifies the square root in Eq. (A8) to √ 2j + 1. Because the Legendre polynomials obey the eigen-equation Eq. (A1) for m = 0 and q = qe z with Eq. (A4) is simply Projecting out the Legendre polynomials one by one, the coefficients A 0,0 ℓ,j (qe z ) obey where we have used with some normalisation N 0 ℓ to be used below, and λ ℓ = λ 0 ℓ (qe z ) to ease notation.
Eq. (A11) can of course be written in matrix form  By demanding that this equation is to be solved to a certain order in q and by assuming α ℓ,j ∈ O q |ℓ−j| , only a finite number of matrix elements and vector elements need to be considered at a time, while also constraining the order of the eigenvalue. Using Eq. (A5) as the starting point, the first order equation for ℓ = 0 is where we have labelled the unknown term of α 0,1 by u 1 q and the unknown term of λ 0 0 (qe z ) by vq. The 0 in the eigenvalue on the right hand side of Eq. (A14) is that of λ 0 ℓ at ℓ = 0, Eq. (A5b). Focussing solely on terms of order q produces two equations, 0 = vq and q/ √ 3 + 2u 1 q = 0, determining both v = 0 and u 1 = −1/(2 √ 3). With those in place, the next order equation is  determining u 1 , u 2 and v, now labelling terms of order q 2 . The only additional difficulty for ℓ > 0 is that the number of rows and columns of the matrix to be considered going from one order to the next may increase by 2 for the first ℓ orders. For example the first non-trivial order for ℓ = 1 is  With this process, we have determined  which can be improved further by demanding that the term of order q 6 in λ 0 0 is consistent with α 0,2 = q 2 /(18 For ℓ = 1, we found  i.e. λ 0 1 (qe z ) = 2 + Matrix expressions such as Eqs. (A17) and (A19) are easily verified with a computer algebra system, simply by determining whether the remainder has the expected order.

(A27)
Appendix C: Field Theory of the Drift Diffusion in an Arbitrary Potential In the following we derive the field theory for a drift-diffusion particle in a finite volume with potential Υ(r) and determine the stationary particle density ρ 0 (r). Our motivation is two-fold: Firstly, we want to reproduce the barometric formula ρ 0 (r) ∝ exp (−Υ(r)/D t ) in the absence of drift, which is the equilibrium situation of pure diffusion, where the diffusion constant D t plays the role of the temperature. The barometric formula is Boltzmann's factor, which we want to see reproduced by the field theory as a matter of consistency, Section C 4. The derivation will provide us with a template for the more complicated case of a drift. Secondly, we want to determine precisely those corrections due to a constant drift with velocity w 0 . We will compare our field-theoretic, perturbative result with the classical calculation in one dimension in Section C 5.
We study this phenomenon in a finite volume with linear extent L, so that the system develops into a stationary state even in the absence of an external potential, which will be treated perturbatively. In a finite volume, the Fourier series representation of any reasonable (finite, continuous, more generally fulfilling the Dirichlet conditions) potential exists. In an infinite volume, we would have to demand that the potential vanishes at large distances, which may result in the stationary density no longer being normalisable even in the presence of the potential.
To simplify the calculations, we further demand the domain to be periodic. The annihilation field χ(r, t) and the Doi-shifted creation fieldχ(r, t) at position r and time t can then be written as with the notation in keeping with the rest of the present work. The finite volume we have chosen to be a cube with linear extension L, easily generalisable to a cuboid or similar. The Fourier modes are k n = 2πn/L with n = (n 1 , n 2 , . . . , n d ) T a d-dimensional integer in Z d . We further use the notation d = d/(2π) and later, correspondingly, δ¯(ω) = 2πδ(ω) to cancel that factor of 2π. The time-independent potential Υ(r) follows the convention Eq. (30), restated here for convenience in arbitrary dimensions

Action
The Fokker-Planck equation of the fully time-dependent density ρ(r, t) of particles at position r and time t, subject to diffusion with constant D t and drift with velocity w, in a potential Υ(r) can be written as removing all ambiguity about the action of the gradient operators by stating that they act on all objects to their right. The action of particles with diffusion constant D t and moving ballistically with velocity w 0 is [8] with k n = |k n | and r ↓ 0 a mass to regularise the infrared and establish causality. An additional external potential results in an additional space-dependent drift term and thus in the action e.g. Eqs. (9b) or (29a) to be treated perturbatively, Eq. (8).
The bare propagator follows immediately from Eq. (C4), and the perturbative part of the action provides a single vertex with a "bauble" to be summed over.

Propagator
With the help of Eqs. (C6) and (C7) the full propagator can be written perturbatively, More systematically, the full propagator may be expanded as where we changed the sign of k p in the creator field compared to Eq. (C8) to improve notation below. The expansion terms E j are defined as follows with the general recurrence amounting to a convolution to be analysed further below.

Stationary state
The density ρ 0 (r) in the stationary or steady state is found by taking of the propagator the limit t 0 → −∞ of the time of the initial deposition, expected to be independent of r 0 and t, provided the system is ergodic and possesses a (unique) stationary state. Provided all internal legs are dashed, i.e. they carry a factor k, the diagrammatic mechanics of this limit amounts to an amputation of the incoming leg in any composite diagram. This is most easily understood by taking the inverse Fourier transform from ω, ω 0 to t, t 0 , while leaving the spatial dependence in k, starting from the propagators that have the form Eq. (C6). The inverse Fourier transform will result in exponentials of the form exp (−it 0 s), where s is the location of the pole of the propagators, here s =i(D t k 2 n +iw 0 · k n + r) with ℑ(s) ≥ 0. Only if the imaginary part of s vanishes, the exponential will not vanish as t 0 → −∞. This happens only when n = 0 and in the limit of r ↓ 0. The resulting residue is thus determined by the remainder of the diagram evaluated at the pole ω = 0 and with the incoming leg replaced by δ p,0 , so that an inverse Fourier transform in k p to r 0 results in a factor 1/L d . As for the individual propagator Eq. (C11a) in Eq. (C10), it contributes indeed only this background, 1/L d . Writing therefore E i with only one arguments as the desired limit, i baubles allows the stationary density to be written in the form via Eqs. (C10) and (C13) and gives from Eq. (C11a) and using Eq. (C6) and where accounts for the subtlety that the right hand side of Eq. (C11b) strictly vanishes for n = 0, because of the k n prefactor and all other terms being finite as r = 0. However, in Eq. (C14) the limit r ↓ 0 has to be taken so that t 0 → −∞ produces a finite result, resulting in using the convention 0∞ = 0. The indicator function Eq. (C18) needs to be in place for every prefactor k that also features in the denominator of a propagator G (k, 0) so that the limits r ↓ 0 and t → ∞ can safely be taken. For easier comparison, we state the expansion of the density to first order in the potential for drift-diffusion particles in a finite volume determined explicitly so far, ei kn·r k 2 n Υ n D t k 2 n +iw 0 · k n + . . .
using Eqs. (C15), (C16) and (C17). With the indicator function Eq. (C18) suppressing unwanted singularities, the recurrence Eq. (C12) remains essentially unchanged at staionarity, as only the outgoing propagator needs to be adjusted to ω = 0, using Eq. (C6). In one dimension we can further simplify, From Eq. (C21) it can be gleaned that in real-space E j (r) obeys the recurrence where we have refrained from introducing a separate symbol for the inverse Fourier transform E j (r) of E j (k). It is instructive to demonstrate that this implies that Eq. (C15) therefore solves the Fokker-Planck Eq. (C3) at stationarity, i.e. we want to show that using Eq. (C23) in Eq. (C15) provides a solution of Eq. (C3) with ∂ t ρ = 0, To see this, we sum both sides of Eq. (C23) over j = 0 . . ., which on the left-hand side produces ρ 0 of Eq. (C15), including the vanishing term as E 0 (r) = L −d , Eq. (C16), so that Eq. (C23) implies and thus Eq. (C24), with the E j in Eq. (C15) determined by Eq. (C21). Eq. (C21) is the central result of the present section. By means of Eq. (C15), it allows the calculation of the stationary density starting from E 0 , Eq. (C16). Without drift, the density should produce the Boltzmann factor, which we will recover in the next section. While this result is easy to obtain from a differential equation approach in real-space r, it is very instructive to perform this calculation in k-space. With drift, on the other hand, we can compare the density ρ 0 resulting from Eqs. (C21) to the known stationary density in one dimension [22], which is done in Section C 5.

Barometric formula
In the absence of a drift, w = 0, the density will evolve into the Boltzmann factor which demands that the probability current is spatially uniform [23]. Provided the potential is bounded from above and from below and with boundary conditions given, the solution is unique by Eq. (C27), including the normalisation constant given by Eq. (C28).
In the following, we show the equivalence of Eq. (C27) and the field-theoretic result Eqs. (C15), (C16) and (C21), by expanding the density ρ 0 from Eq. (C27) in terms of some B j , to be characterised next.

a. Properties of Bj(kn)
We expand the exponential in Eq. (C27) order by order in the potential and introduce Fourier-modes While generally E j (k n ) = B j (k n )/N , we will show in the following that taking the Fourier transform on both sides of Eqs. (C15) and (C32), and using Eq. (C33). Equivalently, we show that B j (k n ) is the sum of the terms of jth order in Υ of N j E j (k n ), We will thus show, beyond Eq. (C26), that the E j are correctly normalised, providing us with an expansion of the density ρ 0 order by order in the potential. Showing this, will further provide us with a useful, non-trivial algebraic identity.
To characterise B j more generally, we introduce the recursion formula by elementary calculations on the basis of Eq. (C31) and similar to Eq. (C21). Starting from Eq. (C36a), it can be used to produce for j ≥ 2 This expression can be simplified further by substituting dummy variables according to p i = n i − n i−1 for i = j − 1, j − 2, . . . , 1 with the convention n 0 = 0. As for this can be allowed to run in the sum by enforcing Eq. (C42) via a Kronecker δ-function, so that eventually for all i ∈ {1, 2, . . . , j}, which labels the index of k pi on the right-hand side as the p i are all equivalent dummy variables. Summing this sum for i = 1, 2, . . . , j therefore gives j times the same result, while the sum of the k pi gives k p1 + . . . + k pj = k n inside the sum, according to the Kronecker δ-function. In conclusion for j ≥ 1 j pj,...,p1 k pj Υ pj Υ pj−1 . . . Υ p1 δ p1+...pj ,n = k n pj ,...,p1 and thus pj ,...,p1 where the indicator function on the left avoids the divergence of k n · k pj /k 2 n . Using Eq. (C46) in Eq. (C43) gives where the second equality is a matter of labelling the indices again in n and n 1...j−1 using Eq. (C42) in reverse, and thus reverting from the indexing in Eq. (C43) to that in Eq. (C41), which does not carry a Kronecker δ-function.
In this form the summation can again be written recursively, similarly to Eq. (C40), which is done after the third equality, by absorbing the summation over n j−2 , . . . , n 1 into B j−1 (k nj−1 ). Eq. (C47c) is the key outcome of the present section.
b. Relationship between Ej(kn) and Bj (kn) We proceed with confirming Eq. (C35) by induction, using Eq. (C47c) in the course. To access this identity, we firstly rewrite Eq. (C35) using Eq. (C39), so that for j ≥ 0 Of the two terms on the right hand side, the first can be rewritten using Eq. (C35), which is assumed to hold for j. This allows us to replace B j (k nj ) by a sum. In the second term we replace δ n,0 by E 0 (k n ) according to Eq. (C16), resulting in The summation over n j ∈ Z d can now be performed on E j−i using Eq. (C21) at w = 0, completing our demonstration that Eq. (C35) holds for j + 1 if it holds for j. With the induction basis Eqs. (C37) or (C38) this concludes our demonstration that the diagrams Eq. (C14) or equivalently (C21) produce a normalised expansion of the density ρ 0 (C15) identical to the Boltzmann factor, order by order in the external potential. The field-theoretic approach thus determines the stationary density without the need of an a posteriori normalisation. And while the field theory is an expansion necessarily in a sufficiently weak potential and necessarily in a finite volume, no such restriction applies to the resulting Boltzmann factor. In the next section we will similarly demonstrate the correspondence between classical and field-theoretic stationary density in the presence of a drift w.

Density in one dimension with drift
In this section, we will demonstrate that the expansion Eq. (C15) in terms of diagrams Eq. (C14) produces the correct, normalised, stationary result in the presence of a potential and a drift in one dimension. The restriction to one dimension is because we are not aware of a closed-form expression beyond.
Following the pattern in Section C 4, we want to show that the diagrammatic expansion Eq. (C15) order by order in the potential reproduces a known expression once this is written out order by order in potential. The difficulty here is similar to that in Section C 4, namely that it is fairly easy to expand the known expression for the stationary density without the normalisation. Working in Fourier-space, we will once again write both the expansion and the normalisation in modes.
The first step is to characterise the expansion terms F and to bring them in a form that lends itself more naturally to a comparison to the field theory. This characterisation involves some rather involved manipulations. In a second step, we will then demonstrate that the resulting expansion agrees with the diagrammatic one.
The classical stationary solution of the Fokker-Planck Eq. (C3) of a drift-diffusion particle in one spatial dimension on a periodic domain and subject to a periodic potential defined so that Υ(x + L) = Υ(x) for x ∈ [0, L) is [22,23], similar to Eq. (C27) defined with an explicit normalisation, The density can be written systematically in orders of the potential up to the normalisation to be determined below. If F j (k n ) is the term of order j in Υ, expanding Z ± (x) = exp (∓w 0 x/D t ) n (±Υ(x)/D t ) n /n! in Eq. (C51), immediately gives (C54) An easy route to perform the Fourier transformation on the right is to split it into two and perform separate Fourier transformations of (−Υ(x)/D t ) i and of the second term, x+L x dy . . . tying both together again by a convolution afterwards. While the first term produces essentially Eq. (C31), the second term involving the integral is more complicated. Writing however which is well-defined as −ik n + w 0 /D t = 0 for all n by virtue of k n and w 0 being real and w 0 = 0, allows the second integral to be removed by an integration by parts. Surface terms then drop out as Υ(x) is periodic. Following at first the pattern of Eq. (C41) the resulting Fourier series reads This expression has some shortcomings that are mostly down to notation, which can be seen in the first sum by setting i = 0 or in the second sum by setting i = j. In both cases, the sums are ill-defined, with the first sum running over indeces n 1 , . . . , n −1 for i = 0 and similarly the last sum for i = j. Choosing i = 1 or i = j − 1 causes similar problems. However, Eq. (C54) allows for all of these choices. The problem is resolved by changing the indeces similar to Eq. (C42), first in the individual Fourier-transforms and then after tying them together in a convolution f j,i (k n ) = 1 (LD t ) j p1,...,pj Υ p1 . . . Υ pj Lδ p1+...+pj ,n 1 − exp (−w 0 L/D t ) −ik p1+...+pj−i + w 0 /D t (C57) with the convention of k p1+...+pj−i = 0 for i = j, as can be verified by direct evaluation of the corresponding Fourier transform. Similarly, we demand δ p1+...+pj ,n = δ 0,n for j = 0, so that f 0,0 (k n ) = Lδ 0,n (1 − exp (−w 0 L/D t ))D t /w 0 , confirmed by direct evaluation of Eq. (C54). We state in passing F j for j = 0, 1 by direct evaluation of Eqs. (C54) and (C57), We proceed by deriving a recurrence formula for F j .
where we have used the equivalence of the i terms in Eq. (C64) to write Eq. (C65a) which carries a factor of i in the numerator in the square bracket, taken the summation over p j and other terms outside the rest of the expression in Eq. (C65b), shifted i by unity in Eq. (C65c), to finally re-express the sum as a convolution of a pre-factor with F j−1 (k n ) using Eq. (C59). Eq. (C65d) is the central result of thise section, providing a simple recurrence formula for F j (k n ) for all j ≥ 1. We will use Eq. (C65d) in the following section to demonstrate the equivalence of the classical result Eq. (C53) and the field theoretic result Eq. (C15) with Eq. (C14) in the presence of drift in one dimension.
b. Relationship between Ej(kn) and Fj(kn) Similar to Section C 4 b, we will now demonstrate that the density ρ 0 expanded in F j (k n ) according to Eq. (C53) with recurrence Eq. (C65d) and with normalisation Eq. (C52), is equal to that expanded in E j (k n ) according to Eqs. (C15) and, in one dimension, (C22). We will do so by demonstrating that the jth order in Υ of N ρ 0 with ρ 0 according to Eq. (C15) equals F j (k n ). Expressing N in terms of F j (k n ) from Eqs. (C52) and (C53), this amounts to the task of showing j i=0 F i (0)E j−i (k n ) = F j (k n ) .
(C66) similar to Eq. (C35), but now allowing for w 0 = 0 and restricting ourselves to one dimension. Eq. (C66) clearly holds for j = 0, as using Eqs. (C16) and (C58a). It also holds for the more complicated j = 1, Eq. (C58b), using Eqs. (C17) and (C58b) at k n = 0, the former in one dimension and the latter producing F 1 (0) = 0. In the last equality of Eq. (C68) we have dropped the indicator function as k n I(n = 0) = k n , as the denominator has no zero at k n = 0 given that w 0 = 0. As in Section C 4 b, we assume Eq. (C66) to hold for j and use Eqs. (C22) and (C65d) to show in the following that Eq. (C66) holds for j + 1. Writing therefore for j ≥ 0 F j+1 (k n ) = (I(n = 0) + δ n,0 )F j+1 (k n ) (C69a) = I(n = 0) LD t qi k n−q (−ik n + w 0 /D t ) Υ n−q F j (k q ) + E 0 (k n )F j+1 (0) (C69b) as in Eq. (C48) using Eq. (C39) in the first line and in the second line using Eq. (C65d) with p j replaced by n − q to rewrite F j+1 (k n ) and δ n,0 replaced by E 0 (k n ) according to Eq. (C16). The F j (k q ) can now be replaced by Eq. (C66), F j+1 (k n ) = I(n = 0) LD t qi k n−q (−ik n + w 0 /D t ) Υ n−q j i=0 F i (0)E j−i (k q ) + E 0 (k n )F j+1 (0) (C70) and the convolution performed on E j (k n ), Eq. (C22), which gives confirming Eq. (C66) for j + 1, provided it holds for j. With Eq. (C71) as the induction step and Eqs. (C67) or (C68) as the base case, this confirms Eq. (C15) with (C22) in the case of drift diffusion in one dimension.

Summary
In the present appendix, we have introduced a general formula to determine the stationary density ρ 0 (x) Eq. (C15) digrammatically, recursively and perturbatively in the potential, Eq. (C21). We have validated this expression for the case of pure diffusion producing the barometric formula, Eq. (C27), in Section C 4 and for the case of drift-diffusion in one dimension, Eq. (C51), in Section C 5. In both cases, we have used induction to show the equivalence of the density derived from theory and the density derived from classical considerations. The key challenge is to compare the field theoretic terms to a certain order in the potential to the corresponding terms in the classical result, which is straight-forwardly expanded in powers of the potential only up to a normalisation, which needs to be expanded itself, Eqs. (C34), (C35) and (C66).
In order to treat the potential perturbatively, we had to allow for vanishing potential, which allows for a stationary state only if the system is finite. However, the resulting perturbation theory reproducing the Boltzmann factor in the case of no drift means that it applies equally for confining potentials. This applies similarly to the case with drift, as the thermodynamic limit L → ∞ may be taken after resumming the perturbative result into the integral expression Eq. (C51).