How to derive a predictive field theory for active Brownian particles: a step-by-step tutorial

The study of active soft matter has developed into one of the most rapidly growing areas of physics. Field theories, which can be developed either via phenomenological considerations or by coarse-graining of a microscopic model, are a very useful tool for understanding active systems. Here, we provide a detailed review of a particular coarse-graining procedure, the interaction-expansion method (IEM). The IEM allows for the systematic microscopic derivation of predictive field theories for systems of interacting active particles. We explain in detail how it can be used for a microscopic derivation of active model B+, which is a widely used scalar active matter model. Extensions and possible future applications are also discussed.


Introduction
The physics of active soft matter [1,2] is among the most rapidly growing subdisciplines of condensed matter physics and statistical mechanics. While the name 'active matter' is applied to a quite diverse range of physical systems, we will, in this review article, use 'active particle' to denote a particle that converts energy into directed motion and 'active matter' for matter consisting of active particles. Examples for active particles include biological organisms like flying birds * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. or swimming bacteria, but also artificially produced objects like Janus particles driven by chemical reactions [3], lightpropelled particles [4], or ultrasound-driven particles [5]. A crucial motivation for the interest of researchers in active matter is, besides its significant potential for applications in medicine and materials science, the remarkable collective phenomena observed in active systems.
Field theories [6] are an extremely useful and frequently applied tool in active matter physics. They provide a coarsegrained description of active systems. Rather than calculating the coordinates of every single particle, one models the system using a small number of order-parameter fields, in many cases only the local particle number density ρ. Active field theories exist in very different ways, ranging from simple ones such as active model B (AMB) [7] and its extension active model B+ (AMB+) [8], which are primarily used for studying active phase separation, over active phase field crystal (PFC) models [9][10][11][12][13] that can describe crystallization, to complex nonlocal dynamical density functional theory (DDFT) [14][15][16]. All these models are overdamped, although extensions to systems with inertia are possible [13,[17][18][19].
Broadly speaking, active field theories can be obtained in two ways [6]. First, they can be derived using phenomenological considerations (top-down approach), for example as modifications of existing theories or using symmetry arguments. Second, one can obtain them by coarse-graining a microscopic particle-based model (bottom-up approach), which, in soft matter physics, is often given by Langevin equations [20]. This approach has the advantage that it provides microscopic expressions for the coefficients that occur in the derived model ('predictive theory'). Coarsegraining is facilitated by following an existing systematic procedure. Such procedures have the advantage that they provide a guidance for which steps to take, and that they allow the resulting theory to be compared to other ones obtained using the same procedure. Such a systematic procedure is given by the interaction-expansion method (IEM). 'IEM' is the name introduced in [21] for the coarse-graining procedure used in [22], which builds up on earlier work in [23][24][25]. The IEM was later applied also to three-dimensional systems [21], circle swimmers [26], active particles in external force fields [27], active particles with an orientation-dependent propulsion speed [28], and active particles with inertia [29]. Although the name 'IEM' is used only in [5,21,[26][27][28][29][30][31] 1 , closely related methods are used also in other derivations in the literature [35][36][37][38]. Therefore, a clear and accessible introduction to the way one can obtain a predictive field theory for active particles using the IEM is of broad relevance for the active matter community.
In this article, we present in detail the derivation of AMB+ from the microscopic Langevin equations governing the dynamics of the individual particles using the IEM. For pedagogical purposes, we show more intermediate steps than in usual presentations. Thereby, this article functions as a 'tutorial' in that it demonstrates all steps required to derive an active field theory from the microscopic dynamics. While the present review focuses on the IEM, these steps are also relevant more generally since the methods employed here, such as gradient expansions or Cartesian expansions, are used also in many other derivation methods. Moreover, we discuss several generalizations of this derivation that have been proposed in the literature. It is easily possible to further generalize the derivation presented here, such that our review is intended to be a starting point for future research employing methods of this type.

Active model B+
Field-theoretical models have a long tradition in active matter physics. Early works include the celebrated Toner-Tu model [39] describing flocks of birds and hydrodynamic models for bacterial turbulence [40]. While these models focus on orientational ordering phenomena, we here focus on theoretical descriptions of the state behavior of spherical active particles that use the particle number density field (rather than an orientational ordering field) as the central order parameter. An important model of this type is AMB, which was introduced in [7] based on earlier work in [23,41]. AMB+, introduced in [8], is an extension of AMB. Further extensions, such as active model H [42] for systems in a momentum-conserving solvent and active model I+ [29] for particles with inertia were developed later. Usually, AMB and AMB+ are introduced as stochastic field theories (although one frequently ignores the noise term [7]). We here present them in a deterministic form since the IEM gives rise to deterministic theories (see section 3.11). Whether deterministic or stochastic models are more useful depends on what one intends to use the model for.
A passive dynamics of type B 2 (overdamped conserved dynamics) is given byρ with a free energy F and a positive mobility M. The form (1), known as 'gradient dynamics' [44], ensures that the system evolves towards the minimum of F. This form is broken in AMB+, which is given bẏ with the free energy functional 3 Here, a, b, κ, λ, and ξ are constant parameters and d is the number of spatial dimensions. Compared to equation (1), AMB+ contains two terms λ ⃗ ∇ 2 ( ⃗ ∇ρ) 2 and ξ ⃗ ∇ · (( ⃗ ∇ρ)( ⃗ ∇ 2 ρ)) that cannot be expressed via the functional derivative of a free energy. The presence of these terms, which break timereversal symmetry, is what makes AMB+ an active field theory. For ξ = 0, AMB+ reduces to AMB. In contrast to AMB+, AMB does not allow for rotational currents since it can be written in the formρ = ⃗ ∇ 2 µ with a generalized chemical potential µ.
One can derive AMB+ phenomenologically, without any microscopic considerations, by noting that the most general theory for the dynamics of ρ that satisfies mass conservation (i.e. −ρ is the divergence of a current), is of second order in ρ and of fourth order in ⃗ ∇, and satisfies translational and rotational invariance is given bẏ where κ 0 , α, and λ 0 are constant coefficients. Adding , defining κ(ρ) = κ 0 + αρ and λ = λ 0 + α/2, and assuming a constant κ gives equation (2).
This sort of argument, however, does not give us microscopic expressions for the model parameters, these simply have to be used as adjustable constants. Moreover, we do not have a clear idea of the approximations required for deriving equation (4) and, therefore, of the range of validity of this model. These problems can be solved by a microscopic derivation, which is possible using the IEM.

The interaction-expansion method
The IEM is a method that combines three typical features of the microscopic derivation of predictive field theories that are applicable on macroscopic scales. The first of these steps is a projection onto order-parameter fields, e.g. the concentration (or density) and the polarization vector. Secondly, the IEM handles the convolution integral emerging in the derivation that stems from the interactions between the particles. Thirdly, it does so in a controlled way that gives analytical expressions for the coefficients which occur in the derived field theory and takes the specifics of the system, like particle length or angular correlations, into account.
All this is achieved via multiple expansions into convenient bases of eigenfunctions that are a good approximate representation of the collective dynamics and that are suitable for the physical system under consideration. Therefore, the IEM is, at its core, quite similar to perturbation theory in, e.g. quantum mechanics. A complicated and typically not analytically solvable problem is reduced to a superposition of known solutions with an adequate choice of prefactors that systematically tries to minimize the difference between the real solution and the approximated one.
In the following, we perform an example derivation stepby-step and explain the mathematics associated with the individual steps. Figure 1 shows the general procedure of the IEM as a flow chart. A microscopic description of an active matter system is typically given in terms of Langevin equations [20,45] for the motion of the individual particles, which can be easily constructed and have obvious physical interpretations. From them, the statistically equivalent Smoluchowski equation can be derived via the Fokker-Planck framework [46]. Integrating out the degrees of freedom of all particles except for one in the Smoluchowski equation gives a dynamic equation for the orientation-dependent one-particle density ϱ(⃗ r,û, t), where the vectors ⃗ r andû are the spatial and angular degrees of freedom and t is the time. This equation is closed by providing (from simulation results or from an analytical theory) an expression for the pair-distribution function [31] appearing in the interaction term of the dynamic equation. In general, this pair-distribution function is only known approximately. Handling the interaction term, which, due to the presence of spatial integrals (nonlocality 4 ) and angular integrals, is relatively complicated, is the primary purpose of the IEM. Several expansions are performed for this purpose-a gradient expansion [47] for the spatial nonlocality and a combined Fourier and Cartesian expansion [48] for the angular integrals. These expansions are exact as long as they are carried out to infinite order 5 . However, as soon as they are truncated at a finite order (that can be chosen depending on how precise one wants the theory to be), the resulting equations of motion are approximate. Truncating the gradient expansion gives rise to local dynamic equations that do not involve integrals. The Cartesian expansion replaces the angle-dependent density ϱ(⃗ r,û, t) by several order-parameter fields ρ(⃗ r, t) (density) and ⃗ P(⃗ r, t) (polarization) that depend only on position and time. Projecting onto ρ and ⃗ P by multiplyingρ with an orientation-dependent function and then integrating over the angular degrees of freedom and then eliminating also ⃗ P via a quasi-stationary approximation (QSA) then gives a closed equation of motion for ρ and a constitutive equation for ⃗ P.

Langevin equations
The starting point of a coarse-graining procedure is a microscopic description of the system of interest. For active matter, several popular options exist [24,49]: Active Brownian particles (ABPs), which move with a fixed speed in the direction of an orientation vectorû whose direction changes via a continuous diffusion process, run-and-tumble particles (RTPs), which undergo directed motion ('run') with constant orientation untilû is changed by a sudden random event ('tumble'), and active Ornstein-Uhlenbeck particles (AOUPs), where modulus and orientation ofû undergo Gaussian fluctuations. The relation of these models and possible unifications are discussed in [24,49]. Here, we focus on ABPs.
The spatial and angular Langevin equations of an ABP are an extension of the standard Langevin equations for passive particles [20]. For simplicity, an isotropic and overdamped particle in two spatial dimensions is considered, resulting in [2]  Here, ⃗ r i = (x 1,i , x 2,i ) T is the center-of-mass position, i is the particle index, a superscript T denotes the transpose, an overdot denotes a partial derivative with respect to t, v A is the bare active propulsion velocity,û(ϕ i ) = (cos(ϕ i ), sin(ϕ i )) T is the normalized orientation vector of a particle with orientation angle ϕ i , Figure 2. Visualization of the microscopic setup considered in the derivation and of the reduced description that exploits translational and rotational invariance. The unit vector in x-direction is denoted asêx. Reproduced from [22]. © IOP Publishing Ltd. All rights reserved.
is the velocity contribution that originates from the particleparticle interaction potential U 2 (∥⃗ r i −⃗ r j ∥) that solely depends on the absolute distance between the ith and the jth centerof-mass positions, ⃗ ∇ ⃗ ri = (∂/∂x 1,i , ∂/∂x 2,i ) T is the del operator with respect to⃗ r i , β is the thermodynamic beta, D T and D R are the scalar translational and rotational diffusion constants of a particle, respectively, and N is the total number of particles. Furthermore, ⃗ Λ T,i (t) and Λ R,i (t) are unit-variance Gaussian white noises with mean zero, i.e.
where ⟨·⟩ is the stochastic average, ⊗ the dyadic product, 1 the unit matrix, δ ij the Kronecker delta, and δ(t − t ′ ) the delta distribution. A visualization of the microscopic setup can be found in figure 2. The Langevin equations (5) are, at their core, force-balance equations: the friction force is directly proportional to the velocity, i.e.⃗ r, and must equal the superposition of all the additional forces. This includes the active propulsion force (∝ v A ), the interaction force (∝ ⃗ v int ), the Brownian force, which originates from the random-like collisions of the particle with the solvent molecules and is given by the terms 6 ∝ Λ, and possible additional forces that can be straightforwardly added to the dynamics. An external force ⃗ F E would result in an extra term +βD T ⃗ F E (⃗ r i , ϕ i , t) in equation (5) and an external torque T E in an extra term +βD R T E (⃗ r i , ϕ i , t) in equation (6). There exist, however, two caveats: • If the diffusivity of an ABP depends on the degrees of freedom of the particle itself, another term ∝ ⃗ ∇D, with the respective diffusivity D, has to be added to the Langevin equations so that they generalize correctly in the Fokker-Planck framework.
• The special case of isotropic two-dimensional particles reduces the angular Langevin equation (6) to scalar ones. More generally they would also be vector equations, which contain generalized noises Λ R,j (t ′ ) → ⃗ Λ R,j (t ′ ) that have a modified self-correlation equivalent to the translational noises.

Smoluchowski equation
The idea of the Fokker-Planck equation [46] is to describe the dynamics of a system by means of the probability density P( ⃗ X, t; ⃗ X 0 , t 0 ), which gives the probability that the system is in the state ⃗ X at time t given that it was in the state ⃗ X 0 at the initial time t 0 . For brevity, the initial sate is omitted in the notation from now on. If the probability dynamics is memoryless, i.e. it changes with only the current state in mind, the underlying stochastic process is a so-called continuous Markov process and can be expressed via the master equation [46] P( ⃗ X, t) =ˆdX ′ (W( ⃗ X, ⃗ X ′ , t) − W( ⃗ X ′ , ⃗ X, t))P( ⃗ X ′ , t). (12) Equation (12), which can be derived from the Chapman-Kolmogorov equation [50], simply states that the change of probability density is the difference of influx and outflux of probability. Consequently, W( ⃗ X, ⃗ X ′ , t) denotes the transition rate from ⃗ X ′ to ⃗ X at time t. By means of the Kramers-Moyal expansion [51,52], equation (12) can be rewritten in the form [46]Ṗ This is the general form of the Fokker-Planck equation with probability drift and probability diffusion Both ⃗ M 1 and M 2 can be calculated by using Langevin equations.
There are again two points to be noted: • The Pawula theorem [53] states that the Kramers-Moyal expansion contains either only two non-vanishing terms (in this case, equation (13) is equivalent to equation (12)) or an infinite number of non-vanishing terms (in this case, equation (13) is an approximation to equation (12)). For the Langevin model with Gaussian white noise considered here, all terms of higher order in the Kramers-Moyal expansion can be neglected. This implies that we can proceed with equation (13). • The state vector ⃗ X is a general object and the Fokker-Planck equation (13) holds generally, i.e. the dimensionality and shape of the particle only enters the exact form of the state vector but not the form of equation (13).
By making use of equations (5), (6) and (13)-(15), we finḋ where P = P( ⃗ X, t) is a short-hand notation and the state vector is given by Equation (16), which is a special case of the Fokker-Planck equation (namely the one describing overdamped systems), is known as the Smoluchowski equation [46].

The first equation of the BBGKY hierarchy
Equation (16) contains all microscopic information about the system-implying that it is typically impossible to solve in practice for a many-particle system, both because it is too complex and because the initial condition (positions and orientations of all particles) is not known. Therefore, one requires a coarse-graining procedure, where a complex microscopic dynamics such as (16) is replaced by a more tractable approximated one. For this purpose, we first define the orientation-dependent one-particle density ϱ(⃗ r, ϕ, t) (for two spatial dimensions) as We have integrated here over the degrees of freedom of all particles except for one. Without loss of generality, this particle can be chosen to be the one with index i = 1. Furthermore, the index is dropped for this specific case, i.e. ⃗ r 1 =⃗ r and ϕ 1 = ϕ. Note that the integration domains depend on the dimensionality. For, e.g. three spatial dimensions, the domain for the translational degrees of freedom changes to R 3 and the angular domain is the unit sphere S 2 that is characterized by the polar and azimuthal angle ϕ ∈ [0, 2π] and θ ∈ [0, π), respectively. One obviously has to take the Jacobian of the transformation for the angular representation into account, which is 1 for two spatial dimensions but sin(θ) for three spatial dimensions and its spherical coordinate representation. As a generalization of equation (18), the n-particle density can be defined as To find the dynamics of ϱ, we integrate equation (16) over the coordinates of all particles except for one. Using equation (18), this giveṡ with the interaction term [22] where U ′ (r) = dU/dr. Equation (20) is a dynamic equation for the one-body density ϱ that depends on the unknown twobody density ϱ (2) . Calculating a dynamic equation for ϱ (2) in the same way (integrating equation (16) over the degrees of freedom of all particles except for two) would, similarly, give an interaction term that involves the three-body density ϱ (3) . This generates a set of coupled differential equations known as the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy [54]. Equation (20) is the first equation of the BBGKY hierarchy. Further equations of this hierarchy have not been used so far in the IEM. In other contexts (such as DDFT [55]), further equations of the BBGKY hierarchy have been employed to obtain dynamic equations that are more accurate than the ones that are commonly used.

Interaction integral
The derivation up to now is rather general, approaches of this form are restricted neither to the IEM nor to active matter physics in general. For example, the derivation of DDFT by Archer and Evans [56] proceeds in exactly this way (although they do not consider orientational dependencies): derive a Smoluchowski equation from the Langevin equations and integrate this equation out to obtain a dynamic equation for the one-body density that is not closed due to the presence of an interaction term. The real difficulty in deriving field theories for systems of interacting particles lies in finding approximations that allow for a closed and solvable model. Different derivation methods employ different strategies to find a closure for the interaction term. Such a strategy allows to replace the expression (21), which involves ϱ (2) , by one in which only ϱ appears. DDFT, for example, relies on the adiabatic approximation [14]. Here, it is assumed that the relations between the two-and the one-body density which hold in equilibrium continue to hold in the dynamical case. These relations then allow for expressing the interaction term using ϱ. The adiabatic approximation is not exactly correct, but often relatively good. However, it corresponds to a close-to-equilibrium assumption, which is not what we want in the context of active matter (which is intrinsically far from equilibrium).
Here, we follow a different procedure. First, we define the pair-distribution function g as [57] Inserting equation (22) into equation (21) gives a closed dynamic equation for ϱ if we know what the pair-distribution function looks like. Fortunately, g is a well-studied object.
To approximate it, one usually assumes that the system is translationally and rotationally invariant. Moreover, the timedependence of g can be neglected for a stationary state. In this case, g can be written as g(r, θ 1 , θ 2 ) with the angles and the parameterization Translational invariance implies that g can only depend on the difference vector⃗ r ′ −⃗ r. This vector is determined by its length r and by the angle ψ R between it and the x-axis. Consequently, g generally depends on three angles (ψ R , ϕ, and ϕ ′ ). Rotational invariance implies that shifting all angles in the same way cannot have an effect on the behavior of the system. Therefore, we can choose to subtract ϕ from all angles. The angle ϕ itself, thereby, becomes irrelevant (zero), such that only two rather than three angles have to be taken into account. These two angles are defined in equations (23) and (24). In earlier applications of the IEM [25], the dependence on θ 2 was neglected, whereas later studies [22] included it. The physical meaning of the symmetry-based reduced variables is illustrated in figure 2. An analytical expression for the pair-distribution function of interacting Brownian spheres that respects translational and rotational invariance was obtained from Brownian dynamics simulations in [31] (see [58] for a Python implementation and [59] for an extension to three spatial dimensions). Inserting this expression into equation (22) gives a relation between ϱ and ϱ (2) , which can be used to close equation (20). The resulting equation of motion can then be solved or (this will be done in practice) approximated further. For the interaction term (21), we now have

Gradient expansion
While the dynamic equation for ϱ can be closed by providing an explicit expression for g, the result is not easy to work with. This has to do, among other things, with the fact that the dynamic equation is non-local, i.e. that the time derivative of the field at position ⃗ r depends on the values of the field at other positions. This non-locality can be removed via a gradient expansion [47], where one replaces the non-local expression (21) by one that depends only on the values of ϱ and its spatial derivatives at position ⃗ r. Gradient expansions are used very frequently in microscopic derivations, not only in active matter physics, but, for example, also in the derivation of PFC models [12,60,61]. We assume that U ′ 2 rapidly goes to zero (short-range interactions). Making the substitution ⃗ r ′ →⃗ r +⃗ r ′ (such that now gives, in the integral in equation (26), a density ϱ(⃗ r +⃗ r ′ ) (we suppress the dependence on angles and time for the moment), which can be Taylor expanded as Inserting equation (27) into equation (26) gives

Fourier expansion of the pair-distribution function
Having removed (or better: evaluated) the integrals over r in the interaction term via the gradient expansion, what remains are the angular integrals over θ 1 and θ 2 (or ψ R and ϕ ′ ). (This is why, in figure 1, we refer to the interaction term (28) as spatially local rather than just as local.) Ideally, what we want are differential equations (rather than integro-differential equations) for order-parameter fields depending solely on ⃗ r and t. Thus, we want to get rid of the integrals by evaluating them explicitly. This can be achieved by expanding the angular dependencies of both g and ϱ. We start with g, for which we use a Fourier expansion. Exploiting that g is real, this expansion reads with the r-dependent expansion coefficients [22] Inserting equation (29) into equation (28) gives

Cartesian order-parameter expansion
Ideally, we want to explicitly evaluate the integrals over ψ R and ϕ ′ in equation (31), since only in this way we can ensure that we do not have to deal with them in the final model. This is not immediately possible as long as ϱ still depends on ϕ ′ . To solve this problem, we expand also the angular dependence of ϱ. For this purpose, we use the Cartesian order-parameter expansion [48,62,63] with the fields that are symmetric traceless tensors for s ⩾ 2 and thus have two independent coefficients at every order s > 0. Here, u i is the ith component ofû and are the tensor Chebyshev polynomials of the first kind. The first three are given by [48] T 0 (û) = 1, Usually, the expansion (32) is truncated at order s = 2. In this case, it reads with the double tensor contraction:, the spatial density 7 the polarization and the nematic tensor The Cartesian order parameter expansion is orderwise equivalent to a Fourier expansion [48,62], i.e. the sum of all terms of a certain order in a Fourier expansion is equal to the sum of all terms of a certain order in a Cartesian expansion. In particular, this implies that the same orthogonality properties hold for both of them, and that one expansion can be rewritten into the other one (explicit conversion tables are provided in [48]). This leads to the question why we use a Fourier expansion for g and a Cartesian expansion for ϱ. While this has, to a certain extent, 'historical' reasons, a possible motivation are the relative advantages and disadvantages of these two types of expansions. The Fourier expansion is mathematically useful since it involves no redundancies (in contrast to the Cartesian expansion, which in two dimensions leads, e.g. to four coefficients at second order of which only two are independent) and since it is relatively easy to do calculations with it. Moreover, the Fourier expansion is much more familiar to most physicists. This makes it reasonable to use the Fourier expansion for g, which is an object that (at least in the present context) primarily serves a calculational purpose. The Cartesian expansion, on the other hand, has the advantage of being more easily interpretable as the expansion coefficients have a more intuitive physical meaning. If, for example, the polarization vector at a certain position points in a certain direction, then this is the direction that the particles located at this position are self-propelling towards. Consequently, it is useful to express results for ϱ in terms of a Cartesian rather than a Fourier expansion.

Approximate equations of motion for the order-parameter fields
To obtain from equation (20) equations of motion for the order-parameter field ϱ i1,...,is , we multiply equation (20) by T i1,...,is (û) and integrate the result over ϕ. This can be done at arbitrary orders [22]. Since our goal here is to give an introduction to this approach (and not to simply reproduce the general results from [22]), we here restrict ourselves to the standard order-parameter fields ρ and ⃗ P. In other words, we truncate the expansion (32) at first order, giving Including the nematic tensor Q, as done in [22], would make the result quantitatively more accurate, but does not change the general idea in any way.
In the Fourier expansion of g (equation (29)), we include only the contributions with n 1 , n 2 ∈ {−1, 0, 1} (as is reasonable since we have truncated the Cartesian expansion (32) at first order) and the gradient expansion (27) at order l = 3 (since we wish to derive AMB+, a theory of fourth order in derivatives-note that one derivative is already present for l = 0). Equation (20) can then be projected onto ρ by simply integrating it over ϕ and then dividing by 2π. The result iṡ ) .
Equation (43) may seem relatively complicated as it still involves several angular integrals. However, this is unproblematic since the integrals over ψ R and ϕ ′ can be easily calculated using a computer algebra system 8 (or by evaluating them by hand if one really likes solving integrals). This gives (dropping dependencies on ⃗ r and t and defining with the coefficients dr rU ′ 2 (r)(g 10 (r) + g −10 (r) 8 For evaluating the integrals in equation (43) Similarly, we can derive a dynamic equation for ⃗ P by multiplying equation (20) withû(ϕ)/π and then integrating over ϕ. The result iṡ with the coefficients dr r 3 U ′ 2 (r)(g 10 (r) + g −10 (r)).
Note that, in equation (52), we have only included terms up to order l = 2 (rather than of order l = 3 as in equation (43)). Even at order l = 2, we have already dropped all terms involving ⃗ P. The reason will be explained in section 3.10.

Quasi-stationary approximation
We have thus successfully derived a field theory for interacting ABPs that involves two fields, the density ρ and the polarization ⃗ P. What we want to derive, however, is AMB+, a theory involving only ρ. One might conclude that it was wrong to truncate the expansion (32) at order s = 1 rather than at order s = 0 and to thereby include a field (the polarization ⃗ P) that we do not wish to appear in our final result. However, if we had done this, then equation (44) would have reaḋ Equation (59) is a rather uninteresting and, in particular, a passive theory. The coefficients A 3 and A 7 are (as can be seen from their definitions in equations (47) and (51)) the ones that arise from the Fourier coefficient g 00 , i.e. the ones that are still nonzero if the pair-distribution function g has no angular dependence (as is the case in a passive system [23]). Physically, the observation that simply ignoring ⃗ P results in a passive theory can be understood based on the fact that, at least for ABPs, the activity is closely related to the orientationdependence. We thus require a strategy for incorporating the effects of ⃗ P into a theory in which ⃗ P does not explicitly appear. We can achieve this using a quasi-stationary approximation (QSA). The QSA is based on the observation that not all observables evolve on the same timescale. Conserved quantities such as ρ typically evolve slower than non-conserved quantities such as ⃗ P. The reason is that the local amount of, e.g. particle density can only change if the density spreads through the system [65]. Suppose now that, after a sufficiently long time, ⃗ P has reached a stationary state. This stationary state will be a solution of equation (52) for⃗ P = ⃗ 0. Solving equation (52) with⃗ P = ⃗ 0 gives ⃗ P as a function of ρ (and its spatial derivatives). Since ⃗ P evolves faster than ρ, this state will be reached on a timescale on which ρ is still changing. If, however, ρ changes, then the stationary solution for ⃗ P will also change. Therefore, ⃗ P relaxes to a new state corresponding to the new ρ. Thus, if we are interested in the time evolution of ρ only on timescales longer than the time that ⃗ P requires to relax, then we may assume that the ⃗ P appearing on the righthand side of equation (44) at time t is the ⃗ P that is the stationary solution of equation (52) for the density ρ that we have at time t. Inserting this solution for ⃗ P into equation (44) gives a closed dynamic equation for ρ only, which is precisely what we were aiming for.
If we, motivated by these considerations, make the QSȦ equation (52) gives Equation (61) may not seem especially useful since it gives ⃗ P in terms of a complicated expression that involves also ⃗ P and its derivatives. However, we can make use of the fact that we can ignore terms of higher than third order in derivatives (see below for an explanation). We simply insert equation (61) into itself (i.e. we replace every occurrence of ⃗ P on the righthand side of equation (61) by the entire right-hand side) and then drop terms of higher than third order in derivatives. The result will still contain terms involving ⃗ P. We thus repeat this procedure one more time and find the constitutive equation We have found an expression for ⃗ P in terms of ρ and its derivatives, valid in the long-time limit. To further simplify the result, we have also dropped terms of higher than second order in ρ. This is a so-called low-density approximation. If we assume that the density is low, then terms of higher order in densities are negligible. Re-arranging terms and applying the vector identities gives A dynamic equation for ρ can then be obtained by inserting equation (66) into equation (44) and dropping terms of higher than fourth order in ⃗ ∇ and second order in ρ. We finḋ This last step allows us to understand why we were able to ignore a number of terms (all of order l = 3 and the ones of order l = 2 that involve ⃗ P) in equation (52), and why we were able to drop all terms of higher than third order in ⃗ ∇ in equation (61). The sole purpose of equation (52) is to produce a constitutive equation for ⃗ P that we can plug into equation (44) in order to derive a dynamic equation of fourth order in gradients. The term with the lowest order in ⃗ ∇ in equation (44) is −(v A /2) ⃗ ∇ · ⃗ P. This term already involves one ⃗ ∇, implying that any terms of fourth order in ⃗ ∇ in equation (52) produce terms of fifth order in equation (44). These terms are irrelevant (we want a theory for ρ of fourth order), such that we can directly ignore terms of fourth order in equation (52). Moreover, equation (61) shows that ⃗ P is given by −(v A /D R ) ⃗ ∇ρ+ 'terms involving ⃗ P and/or more ⃗ ∇s'. If we then start recursively inserting equation (61) into itself, each ⃗ P in a term will lead to an additional ⃗ ∇. Consequently, terms of third order in ⃗ ∇ in equation (52) that involve ⃗ P would effectively be of (at least) fourth order in ⃗ ∇ and therefore irrelevant. This is an important observation since including terms of higher order in gradients and polarizations can very quickly make the derivation significantly more complicated. If we plan our derivation carefully before actually performing it, we can save a lot of time by avoiding unnecessary calculations 9 . 9 In principle, we even could have dropped the term proportional to A 6 in equation (44) right from the beginning by anticipating that ⃗ P is of first order Collecting and simplifying terms in equation (67) using the identities gives equation (4) with the coefficients If one, as is often done [22,25], includes also the nematic tensor Q as a relevant variable, the general procedure remains the same. In this case, one derives a dynamic equation also for Q, derives a constitutive equation for Q, inserts this equation into the dynamic equation for ⃗ P (giving⃗ P as a function of ρ, ⃗ P, and their derivatives), and then proceeds as shown here. Moreover, it is notable that we have combined here several truncated expansions (orientational expansions, gradient expansion, and low-density approximation). In this work, we have combined gradient expansions and low-density approximations simply by applying each of them separately, i.e. we have counted the number of ρs and ⃗ ∇s in a term and discarded it if there were more than two ρs or more than four ⃗ ∇s. An in gradients, such that the entire term will be of at least fifth order and therefore irrelevant. We have kept this term here for completeness since it does not really affect the difficulty of the calculations.
alternative, more sophisticated approach would be based on assigning dynamic weights to, e.g. ρ and ⃗ ∇ and to then demand that the sum of all dynamic weights of all ρs and ⃗ ∇s (the total dynamic weight) is not larger than a certain value [66]. For example, if one assigns the dynamic weight 1 to both ρ and ⃗ ∇ and discards terms with a total dynamic weight larger than three, then the terms ρ 3 or ⃗ ∇ 2 ρ would both be allowed (dynamic weight three), whereas the term ⃗ ∇ 2 (ρ 3 ) (dynamic weight five) would be neglected. In contrast, if we just separately count densities and gradients, then keeping ρ 3 and ⃗ ∇ 2 ρ would require us to keep also ⃗ ∇ 2 (ρ 3 ) (unless there are other reasons to drop this term).
To summarize: the microscopic derivation of the field theory (4) from the Langevin equations (5) and (6) is based on the following approximations: 1. Use of an approximate pair-distribution function g.

Some comments on noise
It is interesting to note that the dynamic equation (4) for the density ρ derived here does not contain noise terms. This might seem like a problematic omission since usual expositions of AMB+ and related models [6] do generally include and discuss such noise terms, and since other derivations [8] give rise to them. This problem has received very limited attention in active matter physics, but was discussed in more detail in the literature on DDFT [14,67,68]. DDFT exists in stochastic [69,70] and deterministic [56,71] variants (these variants differ in whether or not the dynamic equation forρ contains a noise term) and this difference initially gave rise to discussions concerned with which variant is the correct one. Eventually, it was found [67] that these variants differ in the physical interpretation of the density ρ. In deterministic variants, ρ is the ensemble-averaged density, obtained by averaging over various realizations of the Brownian noise [71]. The dynamics for this density does not contain noise terms since one has averaged over them. In contrast, the density ρ in stochastic DDFT is, depending on which variant is used, either the microscopic density operator defined asρ = N i =1 δ(⃗ r −⃗ r i ) [69] or a spatially smoothed field [70], but in any case the empirically observable density of an actual physical system.
In the present work, we have defined the density via an integral over the probability distribution P. The dynamics of P, given by equation (16), is deterministic (while P encodes the noise, it is not itself a noisy quantity), and therefore our density ρ is an ensemble-averaged quantity obeying a deterministic dynamics. Other microscopic derivations of active field theories [8] often rely on the Dean equation, a stochastic differential equation derived by Dean [69] which describes the dynamics of the density operatorρ. (The Dean equation is often referred to as 'Dean-Kawasaki equation', a name that confuses Dean's result with that of Kawasaki [70], who derived a similar equation for the spatially smoothed density [14].) As long as such a derivation does not contain an ensemble average (this is how deterministic DDFT was initially derived from the Dean equation [71]), the resulting dynamic equation will be stochastic and therefore (although its form may be very similar) have a different physical interpretation than the deterministic result obtained here.

Comparison to other derivation methods
We conclude our presentation of the IEM by comparing it to some selected other methods by which active field theories can be derived: • Derivation from the Dean equation: This method (reviewed in [14]), which was discussed already in section 3.11, is relatively popular in active matter physics. Examples where it is used to derive an active field theory include [72][73][74][75][76]. The Dean equation [69] provides an exact stochastic theory for the dynamics ofρ, which can then be coarse-grained. A major difference of this approach to the IEM is that the IEM gives a deterministic and the Dean equation a stochastic theory.

• Derivation via deterministic DDFT: In (deterministic)
DDFT [14], one applies a relation from equilibrium statistical mechanics ('adiabatic approximation'), which allows to calculate correlation functions based on an equilibrium free energy functional, to close the interaction term. This idea was applied to active systems in, e.g. [15,16,77], after having been introduced for passive systems in [56,71]. The advantage of DDFT is that one thereby does not require any further input, such as a pair-distribution function g obtained from simulations as used in the IEM. However, DDFT has the disadvantage that it is restricted to closeto-equilibrium systems since the adiabatic approximation breaks down far from equilibrium. Therefore, extensions of DDFT in which 'superadiabatic forces' (forces not captured in the adiabatic approximation) are included, known as 'power functional theory' [78], have also been applied to active matter [79,80]. • Derivation via PFC models: PFC models [60,81] can be derived as a limiting case of DDFT [61,82]. This also holds for the active PFC model [9,10,12], which has become a popular theoretical description of active matter systems [11,13,17,83,84]. While active DDFT is typically formulated as a dynamical theory for ϱ, the active PFC model uses ρ and ⃗ P as dynamical variables and is, like the IEM derivation, based on an orientational and a gradient expansion (although a QSA is usually not employed in PFC derivations). PFC models are based on the same close-to-equilibrium assumptions as DDFT and are therefore also restricted to the case of low activity. Since they use a local free energy, they are less accurate, but also easier to handle than DDFT. • Derivation via the Mori-Zwanzig formalism: The Mori-Zwanzig projection operator formalism [85][86][87][88][89], reviewed recently in [90][91][92], allows to derive mesoscopic and macroscopic equations of motion systematically from the microscopic dynamics by projection onto a set of (in principle arbitrary) relevant variables. The Mori-Zwanzig formalism is a standard tool in the microscopic derivation of field theories [88,93,94] and has more recently found some applications in active matter physics [95][96][97]. Compared to the IEM, it is much more general-one can use it to derive both deterministic and stochastic models [14], and it can be applied also in other fields of physics [98,99]-but also considerably more complicated. An explicit discussion of the relation between the Mori-Zwanzig formalism and the IEM was provided in [29], where it was shown that the Mori-Zwanzig formalism allows to justify certain approximations made in a derivation using the IEM.

Density-dependent swimming speed
In a system of many active particles, the swimming speed of a single particle can depend on the density of particles in its environment. This effect, which is crucial for the collective dynamics of active particles (see section 4.2), can have two physical reasons [100]. First, it is possible that the particles directly adapt their swimming speed to the particle density via biochemical effects. This is the case for quorum-sensing bacteria [101]. Second, the swimming speed can depend on the density in an effective way since particles collide with each other in a more dense region and are thereby slowed down. This is what happens in systems of ABPs. Within our field-theoretical approach, we can capture this effect by noting that equation (20) contains a term ⃗ ∇ · (v Aû (ϕ)ϱ) on the right-hand side, which accounts for selfpropulsion with a swimming speed v A . Suppose now that we could find a contribution Such a term would describe propulsion with an effective, density-dependent swimming speed v A,eff (ρ). Therefore, we can calculate the density-dependent swimming speed in the IEM framework by looking for a term of the form (78) in equation (20) [22]. Clearly, contributions to such a term can have two originsthe self-propulsion term ⃗ ∇ · (v Aû (ϕ)ϱ) and the interaction term. If we take a look at equation (31), giving an expression for the interaction term, we can further note that the interaction term has the form such that we need to find conditions under which 'something' is equal to f(ρ)û with f being some function of ρ. Looking further at equation (31), we observe that this can only be achieved if we truncate the sum over l at order l = 0 (otherwise there would be too many gradients to reproduce the form (78)) and if we further replace the function ϱ(⃗ r, ϕ ′ , t) by ρ (⃗ r, t), i.e. if we perform a Cartesian expansion truncated at zeroth order (otherwise we would not get a factor ρ). Doing this and then evaluating the integrals over ψ R and ϕ ′ , we find with a constant ζ = −2π 2 βD Tˆ∞ 0 dr rU ′ 2 (r)(g 10 (r) + g −10 (r)).

Motility-induced phase separation
Motility-induced phase separation (MIPS) [100], sometimes also referred to as 'active phase separation' [7], is one of the most widely studied phenomena in active matter physics. It corresponds to liquid-gas phase separation, i.e. to the separation of particles in regions of high and low density, in a system of purely repulsive identical particles. Such phase separation would be impossible in passive systems, but is very commonly observed in active ones. Physically, MIPS arises from a combination of two effects [100]. First, active particles tend to accumulate in regions where their speed is slow. Second, as discussed in section 4.1, the speed of active particles depends on the density. This leads to a positive feedback loop: if the density is higher in a certain region, then the particles are slower there, implying that the particles accumulate in this region and the density gets even higher, implying that the particles become even slower and so on.
Interesting insights into MIPS can be obtained already using a model of second order in derivatives. Dropping terms of fourth order in derivatives, equation (4) readṡ with the density-dependent diffusion coefficient A linear stability analysis of equation (82) gives the spinodal condition [22,23] Intuitively, the condition (84) can be understood as follows: For D(ρ) > 0, equation (82) essentially looks like a diffusion equation (although one with a density-dependent diffusion coefficient). In general, the diffusion equation predicts that there is a net motion of particles towards regions with lower density, such that density differences between two regions will vanish after some time. Mathematically, this can be seen from the fact that the current, in equation (82) given by −D(ρ) ⃗ ∇ρ, points in the opposite direction to ⃗ ∇ρ. This changes for D(ρ) < 0. In this case, the particle current changes its direction compared to the normal diffusion equation. Therefore, particles tend to move towards denser regions, such that initial density differences between two regions increase. The consequence is phase separation.
To explicitly calculate the spinodal resulting from equation (84), we require the numerical values of a and b. For this purpose, we use as a reference the simulation results from [31]. Here, the Weeks-Chandler-Andersen potential [102] with the effective particle diameter σ and the interaction strength ϵ was used. For simulations of ABPs interacting via such a potential, a convenient system of units is given by the Lennard-Jones units where σ is the unit of length, τ = σ 2 /(ϵβD T ) (Lennard-Jones time scale) is the unit of time, and ϵ is the unit of energy. As dimensionless parameters, the Péclet number and the packing density with the spatially averaged number densityρ can then be used.
One typically calculates the spinodal as a function of Pe and Φ, such that we should express a and b via these variables. The simulations from [31] set v A = 24 σ τ (88) and varied Pe by changing D T . Using also the relation (derivable from the Stokes-Einstein-Debye relation) gives Calculating b, given by equation (73), is more complicated as it depends on the parameters A 1 , A 3 , and A 10 , which are given by equations (45), (47) and (55) and thus depend on the pair-distribution function g. For the interaction potential (85), an analytical representation for g was obtained via a fit to particle-based simulations in [31]. From these results, it can be found that A 1 , A 3 , and A 10 are approximately given by 10 βD T A 10 = 147.8.
(A software implementation of the analytical representation for g can be found in [58].) What is notable here is that these expressions depend on the packing density Φ, which in 10 Note that [22,31] use a different notation convention. State diagram comparing simulation data for ABPs from [31] and the theoretical prediction for the spinodal and the critical point from [22] ('predicted spinodal' and 'predicted critical point') with earlier theoretical predictions for the spinodal (from [23,25]) and the critical point (from [104,105]). Reproduced from [22]. © IOP Publishing Ltd. All rights reserved.
the homogeneous phase-which is where the pair-distribution function was obtained [31]-is simply proportional to the density ρ. This implies that the assumption made in the derivation of AMB+ that coefficients like b are simply densityindependent constants is a strong approximation. When solving equation (84) to get the spinodal, it should be taken into account that b is a function of Φ.
Combining equations (73), (86), (88), (89) and (91) The fact that it is possible to calculate the parameters a and b from microscopic properties of the system is what makes theories such as the one discussed here predictive (since they are able to predict parameter values). In figure 3, the theoretical prediction for the spinodal for MIPS obtained in [22] (the derivation presented here is a simplified version of the one in [22]) is compared with spinodals obtained in earlier studies [23,25] and with simulation data for ABPs from [31]. Figure 3 shows the spinodal as a function of Pe and Φ. The color code indicates the characteristic length L c in units of σ as calculated from the simulations. L c measures the length scales on which density inhomogeneities occur [59,103]. A high value of L c implies density inhomogeneities on large scales and therefore phase separation. As can be seen, the agreement of the prediction (84) with the simulation results is good. What is also shown and compared to earlier predictions [104,105] is the predicted position of the critical point. (Although 'critical point' is a thermodynamic notion, one can give a generalized definition applicable also to active systems [106]. It is also possible to calculate the critical point of active systems directly from simulation data [28,105].) Also for the critical point the agreement of the theoretical prediction with the simulation data is good.
The application to MIPS also allows to understand why the pair-distribution function g is such an essential quantity in active matter physics. As discussed in section 3.10, assuming that g depends only on r gives (for our system of spherical ABPs) rise to a passive field theory, whereas the angular dependence is related to activity. In [23], it was found that the anisotropy of g is related to MIPS. The pair-distribution function measures the probability of finding two particles in a certain configuration [59]. In the case of MIPS, it is, for a given active particle, more likely that there is another particle in front of it than behind it. (Intuitively, particles swimming towards each other get stuck [107].) This can only be described by a function g that depends also on the particle orientations.

Extensions
We have presented here in detail a derivation for the simplest possible case-spherical overdamped ABPs in two spatial dimensions. The IEM can, however, be applied also in much more general setups. Therefore, we now present some extensions of the simple IEM shown above. Rather than re-doing the entire derivation, we explicitly discuss in which way these extensions differ from the 'standard' one. This is possible since the overall strategy remains the same in all these cases.

Three spatial dimensions
Given that we live in a three-dimensional world, an important extension of the derivation shown in section 3 is a generalization towards three spatial dimensions. Such a derivation was performed in [21]. The corresponding pair-distribution function was obtained in [59]. Here, we briefly summarize how the derivation from [21] differs from the one discussed in section 3.
The main difference to the two-dimensional case is the treatment of the orientational degrees of freedom. In two dimensions, one can always specify the orientation of a (hard) particle using one angle ϕ. In three dimensions, the required number of angles depends on the shape of the particle. We first consider (as is done in almost all theoretical investigations of both passive and active matter) only particles with an axis of continuous rotational symmetry. Examples of such particles are rods and active spheres. For such particles, known as uniaxial particles, two angles ϕ ∈ [0, 2π) and θ ∈ [0, π), the usual angles of spherical polar coordinates, are sufficient. In dynamic equations such as equations (16) and (20), a partial derivative ∂ ϕ with respect to ϕ is replaced by the angular momentum operator ⃗ R =û × ⃗ ∇û, where ⃗ ∇û is the del operator containing partial derivatives with respect to the elements ofû.
We thus have to extend both the angular expansion (29) and the Cartesian expansion (32) to the case of a dependence on two angles. The Fourier expansion (29) then generalizes to an expansion in spherical harmonics Y m l (û) [108] (here with Condon-Shortley phase convention), given by [21] with the expansion coefficients the Legendre polynomials P l (û), the Clebsch-Gordan coefficients C(ll ′ l r , mm ′ m r ), and the complex conjugation * . For the expansion of ϱ, we can use the three-dimensional Cartesian expansion [48,108] with the fields ϱ i1,...,is (⃗ r) = 2s + 1 4πˆ2 π 0 dϕˆπ 0 dθ ϱ(⃗ r,û)T i1,...,is (û) (98) and the tensor Legendre polynomials Also in the three-dimensional case, the angular and the Cartesian expansion are orderwise equivalent. A complication compared to the two-dimensional case is, however, that the number of independent expansion coefficients becomes larger at higher orders. In two dimensions, we have two independent coefficients at every order s > 0, whereas in three dimensions, we have 2s + 1 independent coefficients at every order.

Particles with arbitrary shape
For particles with a general shape (biaxial particles), three angles (Euler angles) are required. One can, for this purpose, further extend the expansion (97) to the case of a dependence on three angles. This procedure is explained in detail in [48].
In [66], the IEM was used to derive a field theory for active particles with arbitrary shapes.

Circle swimmers
In [26], the IEM was applied to systems of circle swimmers. These are characterized by an additional constant term ω (which is the angular velocity of a free particle) appearing on the right-hand side of equation (6). From the modified Langevin equations, one can then carry out the derivation in the way presented in section 3. The authors of [26] derived a model of second and fourth order in spatial derivatives. While the second-order model can be mapped onto the standard second-order model (82) by introducing an effective temperature and effective rotational diffusion coefficient, a direct mapping is no longer possible in the fourth-order model, where the additional contributions occurring in circle swimmer systems give rise to a chiral current involving a mixing of spatial derivatives.

Mixtures
Another extension, considered in [25] (with stronger approximations for the angular dependencies of g than in later studies [22]), is to study mixtures. A typical example is the investigation of mixtures of active and passive particles [12], but also mixtures of different active species can be considered [25]. Essentially, the main difference is that one has to consider separate fields ϱ ς , ρ ς , and ⃗ P ς for every species (ς is then an index for the species). Also, the pair-distribution function is then a function g ςν with two indices.

External fields
The dynamics of ABPs in external fields was studied using the IEM in [27]. Here, one adds an external force ⃗ F ext on the right-hand side of equation (5). Essentially, the derivation again goes through as presented in section 3. At second order in derivatives, one gets the expected resulṫ with the mobility Γ, i.e. the external field simply gives an additional contribution ρ ⃗ F ext in the density current. This is also the way an external field would manifest itself in, e.g. DDFT. At fourth order in derivatives, one finds, however, that the density current gets additional contributions in which ⃗ F ext is multiplied with nonlinear functions of ρ and ⃗ ∇ρ-a rather counterintuitive result that is in strong contrast to what one would expect based on theories such as DDFT. The origin of these contributions is the QSA discussed in section 3.10. Since the transport equation for ⃗ P will also involve ⃗ F ext , ⃗ F ext will appear in the constitutive equation for ⃗ P obtained from setting⃗ P = ⃗ 0, and by inserting the constitutive equation for ⃗ P into the contributions in the interaction term where ⃗ P is multiplied with nonlinear functions of ρ and ⃗ ∇ρ, one gets terms in which the same happens to ⃗ F ext .

Orientation-dependent propulsion speed
The derivation in section 3 has assumed that v A is simply a constant parameter that does not depend on ⃗ r or ϕ. While this assumption is made in most derivations of active field theories, a number of recent studies have investigated, both theoretically and experimentally, systems where v A depends on position [109][110][111] and orientation [28,112]. A dependence of v A on ϕ can arise, e.g. in light-driven [4] or ultrasound-driven [113] particles. Such dependencies provide a way of controlling the collective dynamics of active particles.
With this motivation, a field theory for the collective dynamics of active particles with orientation-dependent propulsion was derived using the IEM in [28]. This derivation is essentially parallel to the one shown above, with one crucial difference: If v A in equation (5) depends on ϕ, then, in order to remove all orientational dependencies from the final theory for ρ, we need to perform an orientational expansion also for v A . This expansion takes the (Cartesian) form v A (ϕ)û(ϕ) = ⃗ µ (1) with the orientation-averaged propulsion velocity and the symmetric velocity tensor This allows to derive the field theorẏ with the diffusion tensor D(ρ) = (D T + c 1 ρ + c 2 ρ 2 )1 + c 3 ρµ (2) + µ (2) · µ (2) 2D R (105) and constant coefficients c 1 , c 2 , and c 3 . Note that in deriving equation (104), only terms of up to second order in ⃗ ∇ have been taken into account.

Inertia
Active particles studied in experiments, such as microswimmers, are typically subject to very large friction forces and therefore well described by the overdamped Langevin equations (5) and (6). This, however, is not true for all active particles. Important counterexamples are vibrated granulates and flying insects [114]. In recent years, inertial active particles have attracted an increasing amount of attention (see [107] for a review). An extension of the IEM towards inertial dynamics was proposed in [29]. Here, we explain the main ideas of this extension.
If the translational (but not the rotational) motion is underdamped, the momentum ⃗ p of the particles has to be used as a state variable in addition to the position ⃗ r and orientation ϕ. Consequently, the many-body distribution P depends not only on the positions and orientations, but also on the momenta of all particles. Integrating the dynamic equation for P (the inertial extension of equation (16)) over the coordinates of all particles except for one gives a dynamic equation for the orientation-and momentum-dependent one-particle density P 1 (⃗ r,⃗ p, ϕ, t). One then defines the orientation-dependent oneparticle density as ϱ(⃗ r, ϕ, t) =ˆd 2 p P 1 (⃗ r,⃗ p, ϕ, t) (106) and makes the generalized local equilibrium approximation with the particle mass M, the Boltzmann constant k B , the temperature T, and the generalized velocity field ⃗ v. Equation (107) implies From the dynamic equation for P 1 (which constitutes the first equation in the BBGKY hierarchy), one then obtains dynamic equations for ϱ and ⃗ v, which play the role of equation (20) for the subsequent derivation-not because these equations are the first equations of the BBKGY hierarchy, but because they constitute unclosed governing equations for order-parameter fields depending on ⃗ r and ϕ that can be treated in a similar way as equation (20). The dynamic equation for ⃗ v contains an interaction term that is approximated via Fourier and gradient expansions in the standard way. Since we have two ϕdependent order-parameter fields, we require, in addition to the Cartesian expansion (32) of ϱ, a similar expansion for ⃗ v. This expansion reads with the local velocity ⃗ v(⃗ r) = 1 2πˆ2 π 0 dϕ⃗ v(⃗ r, ϕ) (110) and the local velocity polarization v ⃗ P (⃗ r) = 1 πˆ2 π 0 dϕû ⊗⃗ v(⃗ r, ϕ).
One thereby obtains coupled equations of motion for the four order-parameter fields ρ, ⃗ P, ⃗ v, and v ⃗ P . Due to the larger number of order-parameter fields, a choice needs to be made for the QSA regarding which fields are considered the slow ones. This depends on the properties of the physical system in question. In the usual case of strongly damped systems, such as microswimmers at low Reynolds number, one would assume the velocity field ⃗ v to relax quickly compared to the polarization field ⃗ P, whereas te Vrugt et al [29], in contrast, considered the case of weak damping (and large D R ) in which ⃗ v can be assumed to be slow compared to ⃗ P. This is a specific choice made in [29], in general the IEM allows to consider also other limits.

Conclusions
In this article, we have provided an introduction to the microscopic derivation of predictive field theories for active matter using the IEM. As an example, we have discussed in detail all steps required for a derivation of AMB+ using this method. Thereby, we have also covered a number of theoretical techniques that are important also beyond this specific application, such as orientational expansions or gradient expansions. Moreover, we have covered several extensions of the simple derivation for spherical ABPs in two dimensions from the literature.

Outlook
The IEM-and field-theoretical modeling of active matter in general-has significant potential for future research. As can be seen from the list of extensions, the method is quite flexible, allowing it to be applied also in contexts for which it was not originally developed. Possible further applications include field theories for particles with non-reciprocal interactions [38,115], less approximate models for mixtures [25], or the incorporation of hydrodynamic interactions [116]. Moreover, one could investigate in more detail the relation of the IEM to other approaches such as DDFT or the Mori-Zwanzig formalism in order to understand in more detail which approach is the best one in a certain given context.

Data availability statement
No new data were created or analyzed in this study.