From a microscopic inertial active matter model to the Schrödinger equation

Active field theories, such as the paradigmatic model known as ‘active model B+’, are simple yet very powerful tools for describing phenomena such as motility-induced phase separation. No comparable theory has been derived yet for the underdamped case. In this work, we introduce active model I+, an extension of active model B+ to particles with inertia. The governing equations of active model I+ are systematically derived from the microscopic Langevin equations. We show that, for underdamped active particles, thermodynamic and mechanical definitions of the velocity field no longer coincide and that the density-dependent swimming speed plays the role of an effective viscosity. Moreover, active model I+ contains an analog of the Schrödinger equation in Madelung form as a limiting case, allowing one to find analoga of the quantum-mechanical tunnel effect and of fuzzy dark matter in active fluids. We investigate the active tunnel effect analytically and via numerical continuation.


I. INTRODUCTION
The study of active particles has become one of the fastest-growing fields of research in soft matter physics and statistical mechanics due the enormous number of novel and interesting effects that active matter can exhibit. Among these effects are a plethora of analogies between active matter and quantum mechanics. This includes Bose-Einstein condensation [1][2][3], Hall viscosities [4,5], orientational order in systems of fully symmetric particles [6,7], Schrödinger-type dynamics in polar liquids [8], spin-orbit coupling [9], time crystals [10], and topological effects [11]. A very simple yet extremely powerful description for active matter is given by scalar active field theories such as active model B (AMB) [12] and the more general active model B+ (AMB+) [13]. These provide a minimal description for many remarkable effects such as active phase separation and have led to crucial insights into the thermodynamics of active matter [14][15][16][17][18].
While such field theories have also been coupled to the momentum-conserving dynamics of the solvent [16,17,19,20], the inertia of the active particles themselves has been ignored in this context. However, recent experiments [21][22][23] have found that the inertia of active particles is important in a variety of contexts. Moreover, theoretical and experimental studies have found a number of remarkable effects associated with inertial active matter [24], ranging from self-sustained temperature gradients [25] through restored equilibrium crystallization [26] to damping-dependent phase boundaries [27]. * u.thiele@uni-muenster.de; http://www.uwethiele.de † Corresponding author: raphael.wittkowski@uni-muenster.de Consequently, there has been a strongly increasing recent interest in inertial active matter [28][29][30][31][32]. Field theories for inertial active matter have been derived in Refs. [27,33,34] as extensions of the active phase field crystal (PFC) model [35][36][37][38]. Active PFC models can be derived as an approximation of dynamical density functional theory (DDFT) [39], and have two disadvantages compared to AMB+. First, they rely on the closeto-equilibrium ("adiabatic") approximation that DDFT is based on, and second, they require two order parameter fields (density ρ and polarization P ) rather than just one, making them more complex. In contrast to PFC models, to the best of our knowledge, no extension of AMB+ to the underdamped case has been derived yet. A second gap is that up to now the relation between inertial active matter and quantum mechanics remains unexplored.
In this work, we close both of these gaps. First, we obtain, via a microscopic derivation using the wellestablished interaction expansion method [40][41][42][43], an extension of AMB+ to particles with inertia that we refer to as active model I+. The derivation systematically generalizes previous work on this subject and also reveals some remarkable new physics. In particular, it is found that thermodynamic and mechanical definitions of the velocity field lead to different results in the active case, and that the density-dependent swimming speed of active particles gives rise to an effective viscosity of the active fluid. Second, we show that active model I (the underdamped analogon of AMB) contains, as a limiting case, the Madelung equations [44,45], which constitute a hydrodynamic representation of the Schrödinger equation [46,47]. This allows us to find analoga of the quantummechanical tunnel effect and of fuzzy dark matter in an active fluid. A numerical investigation of the active tunnel effect using continuation methods shows that it also occurs when the approximations required for the active-quantum mapping are not exactly satisfied. This implies its robustness as physical phenomenon.

A. Active model I+
Our starting point is AMB+ [13], which is given bẏ (1) with a local particle number density ρ( r, t), a mobility M , an (overdamped) free-energy density f o (ρ) typically assumed to be a fourth-order polynomial, the notation f o = ∂ ρ f o , and constants κ, λ, and ξ. An overdot denotes a partial derivative with respect to time t. The model (1) is overdamped. Typically, one makes the simplifying assumption of a constant mobility M ρ ≈ M 0 , which is valid only in uniform states or close to a critical point, but qualitatively reasonable also in other cases. The purpose of this approximation is to get a noise that is additive rather than multiplicative [18]. Here, we do not have a noise term since our microscopic derivation interprets ρ as an ensemble-averaged density [48]. By setting ξ = 0 in Eq. (1), one obtains AMB [12]. AMB, in turn, can be thought of as a minimal extension of the Cahn-Hilliard equation [49] to the active case.
A remarkable feature of AMB and AMB+, which distinguishes them from passive field theories, is that the right-hand side of Eq. (1) cannot be written as a gradient dynamics, i.e., in terms of the functional derivative of a free energy [12,50]. In addition, AMB+ allows (unlike AMB) for circulating currents in steady state [13]. One can derive AMB+ either phenomenologically by writing down a general theory of a certain order in gradients and fields (top-down approach) or microscopically by explicit coarse-graining of the microscopic equations of motion of active particles (bottom-up approach) [18]. Here, the bottom-up approach has the advantage of providing explicit expressions for the coefficients appearing in the model (predictive theory) [42,43] and giving a clearer insight into the origin of the various terms and the approximations required to get them.
Since AMB+ is overdamped, it does not take the inertia of the active particles into account. In this work, we obtain via a microscopic derivation an extension of AMB+ to the underdamped case, which we will refer to as active model I+ (AMI+), with "I" standing for "inertial". It is given bẏ with the velocity field v, the rotational diffusion coefficient D R , the free energy density f , its derivative f = ∂ ρ f , the particle mass m, the friction coefficient γ = 1/(mM ), and the local density-dependent swimming speed Here, v 0 is the propulsion speed of a free particle and A 1 is a constant (see Eq. (37) in Section II B). We have also included an external potential U 1 for generality. The form (4) agrees with the expression derived in Ref. [43]. In overdamped fluids, the existence of a densitydependent swimming speed -that can arise, e.g., from particle collisions in the case of active Brownian particles (ABPs) considered here or from quorum-sensing in the case of bacteria -is essential for the phenomenon of motility-induced phase separation (MIPS), where repulsively interacting particles phase-separate (which would not be possible in a passive system) [51]. From AMI+, we can see that v ld (ρ) plays two interesting roles in the underdamped case: 1. It leads to a second term in the continuity equation (2) in addition to the well-known passive term ∇ · (ρ v). This second term is related to the selfpropulsion term known from the active PFC model (see Section IV A).
2. It gives rise to an effective viscosity 1 v ld (ρ) 2 /γ. This implies that a system of underdamped active particles should behave more like a viscous fluid for larger activity (larger v ld ) and more like an ideal fluid for larger density (smaller v ld ).
AMI+ contains AMB+ as a limiting case. Showing this requires two approximations: 1. We assume the system to be overdamped (large γ), i.e., we set the material derivative˙ v + ( v · ∇) v in Eq. (3) to zero, solve the resulting equation for v and insert the result into Eq. (2). (This is analogous to the procedure required for deriving overdamped from underdamped DDFT [39].) 2. Using Eq. (4), we write in Eq. (2) with the effective free energy density where Λ is the (irrelevant) thermal de Broglie wavelength, and then define f o = f + f e and f e = ∂ ρ f e . Equation (6) shows that we can interpret v 2 0 /(2M D R ) as a shift of the temperature [53], since the first term on the right-hand side has the form of an ideal gas free energy.
Interestingly, as shown in Section IV A, the microscopic derivation reveals another form of "effective temperature" that is a feature of inertial active matter. The free energy f appearing in Eq. (3) and microscopically given by Eq. (106) has, as a prefactor in the ideal gas term, a factor k B T + mv 2 0 /2 rather than k B T as in the passive case. This shows that the active kinetic energy mv 2 0 /2 plays the role of a thermal energy in inertial active matter.
By taking the curl of Eq. (3) and defining the vorticity ω = ∇ × v, we can obtain the active vorticity equatioṅ Starting from AMI+, we can again make two approximations: 1. We assume that v ld (ρ) is small. From Eq. (4), we can see that this assumption is justified if v 0 and A 1 are both small (i.e., in the case of weak activity) or, for larger activities, if v 0 ≈ A 1 ρ/(γm).
2. We drop the term proportional to ξ, such that the material derivative of v is given by the sum of the gradient of a generalized chemical potential and a damping term. Setting ξ = 0 is the usual approximation by which one gets from AMB+ to AMB.
We then obtain the simpler active model I (AMI), which is given byρ Equation (9) can be written aṡ with a generalized chemical potential It is straightforward to obtain AMB from AMI by taking the overdamped limit.

B. Microscopic derivation of active model I+
Microscopically, a two-dimensional system of N underdamped ABPs is described by the Langevin equations [33]˙ where r i (t), p i (t), and ϕ i (t) are position, momentum, and orientation (direction of self-propulsion force) of the i-th particle,û(ϕ) = (cos(ϕ), sin(ϕ)) T is its orientation vector, m is its mass, γ and γ R are the translational and rotational friction coefficients, D and D R are the translational and rotational diffusion coefficients, U = U 2 +U 1 the potential consisting of interaction potential U 2 and external potential U 1 , η i (t) and ξ i (t) are translational and rotational Gaussian white noises with zero mean and unit variance, and v 0 is the particles' self-propulsion velocity. The corresponding Fokker-Planck equation is given by [33] where P N is the N -body probability distribution and is the Liouvillian. The dependence on t is not written explicitly to simplify the notation. By integrating Eq. (15) over the coordinates of all except for one particle, we find [33] ∂ t P 1 ( r, p,û) +ˆd 2 r 2ˆd 2 p 2ˆd ϕ 2 ( ∇U 2 ) · ∇ p P 2 ( r, r 2 , p, p 2 , ϕ, ϕ 2 ) with ∇ = ∇ r and the one-and two-particle distribution functions P 1 and P 2 defined as [54] P n = N ! (N − n)!ˆd 2 r 1 · · ·ˆd 2 r N −nˆd 2 p 1 · · ·ˆd 2 p N −n where n ∈ {1, . . . , N }. The index 1 is dropped for the coordinates. We define the particle density ( r,û) =ˆd 2 p P 1 ( r, p,û) (19) and apply the local equilibrium approximation [54] with the Boltzmann constant k B , the temperature T , and a generalized velocity field v( r,û). The expression (20) may also be seen as the zeroth-order term of an expansion of a general distribution P 1 around local equilibrium; the higher-order terms then give viscous corrections [54]. Equation (20) implies Although Eq. (20) is a common approximation, the fact that we allow v to depend onû (i.e., that we use a generalized velocity) distinguishes our approach from previous passive [54] and active [33] inertial field theories. We now drop arguments of the fields unless unclear. Integrating Eq. (17) over p and using Eqs. (19) and (21) yields Similarly, we can multiply Eq. (17) by p/m, integrate over p, and use Eqs. (19)- (22) to get with the interaction term I =ˆd 2 r 2ˆd ϕ 2 2 ( r, r 2 , ϕ, ϕ 2 ) ∇U 2 (r), where 2 =´d 2 p´d 2 p 2 P 2 is the two-body distribution function and r = r− r 2 with the Euclidean norm · is a distance. The derivation of Eq. (23) generally follows the standard procedure of deriving hydrodynamic equations from microscopic dynamics [54]. Our result differs from the standard form of velocity transport equations by the presence of the term −D R v(∂ 2 ϕ )/ , which arises from the term D R ∂ 2 ϕ in Eq. (22). We can define the pair-distribution function g as [40,43] g( r, r 2 , ϕ, ϕ 2 ) = 2 ( r, r 2 , ϕ, ϕ 2 ) ( r, ϕ) ( r 2 , ϕ 2 ) .
Following the treatment in Ref. [43], we assume the pair-distribution function to be translationally and rotationally invariant, implying that it can be written as g(r, θ 1 , θ 2 ) with the angles θ 1 = ϕ R − ϕ and θ 2 = ϕ 2 − ϕ and the parametrization r − r 2 = rû(ϕ R ). Then, we can perform a Fourier and a gradient expansion [43,55] of g and find with the r-dependent expansion coefficients [43] We now carry out the Cartesian orientational expansions [6] ( r,û) = ρ( r) +û · P ( r), with the non-orientational particle density the local polarization and the local velocity polarization with the dyadic product ⊗. Here, our treatment differs in an important way from standard treatments of active overdamped [42,43], passive underdamped [54] and even active underdamped [33] particles. Since we have a velocity field v that also depends onû, we have to perform the orientational expansion not only for the density, but also for the velocity. We now insert Eq. (28) into ln(Λ 2 ) and Taylor expand around P = 0. This gives where we have replaced ρ by a spatially and temporally constant reference density 0 in the last step. Similarly, we insert Eqs. (28) and (29) into v(∂ 2 ϕ )/ and Taylor expand around P = 0 to find Finally, an orientational expansion of the interaction term gives with the coefficients These coefficients can be time-dependent by inheriting a time-dependence of g [43], but we will assume them to be constant. From Eqs. (22), (23), (28), (29), and (34)-(36), we obtain the general local field theory for underdamped ABPṡ with the unit matrix 1.
As is standard in passive model B [12], we make the simplifying assumption that κ is constant. Also, we defineλ + δ/2 = λ. Then, Eq. (45) reduces to Eq. (3) of AMI+. In the passive limit, we have v 0 = 0, A 1 = 0, and A 3 = A 4 = 0 since the pair-distribution function is isotropic in the passive case [57]. In this case, most terms on the right-hand side of Eq. (3) vanish, such that it can be written as − 1 m ∇ δF δρ with a free energy F =´d 2 r (f (ρ) + ρU 1 ). (Of course, this does not mean that there can be no gradient terms in the passive case. These can be obtained by a more sophisticated treatment of interaction terms.)

C. Density-dependent swimming speed
The existence of a density-dependent swimming speed is crucial for the phenomenon of MIPS in active matter. In the case of ABPs, the particles are slowed down by collisions with other particles in a dense phase, which reduces the effective swimming speed in such a phase and creates a feedback that leads to phase separation [51].
One can calculate the density-dependent swimming speed v D from the interaction-expansion method by looking for a contribution of the form ∇ · (v D [ρ]û ) in the dynamic equation for [42]. Inserting Eq. (29) into Eq. (22) gives˙ This shows that the role of the density-dependent swimming speed is, in our extended theory, played by the tensorial quantity v P . The part of this tensor that is proportional to 1 then directly gives us v D . As shown in Section IV A, this tensor is given by Equation (48) provides a microscopic expression for the density-dependent swimming speed in the active fluid. Interestingly, the "density-dependent" swimming speed depends not only on the density ρ, but also on the velocity v. This suggests that v also has to be taken into account when describing the emergence of MIPS in underdamped active fluids.

D. Mechanical vs thermodynamic velocity
It is worth briefly discussing here the definition of the velocity field v. In the theory of classical passive fluids, this can be done in two ways: 1. Mechanically: The density ρ obeys a continuity equationρ = − ∇ · g with the momentum density g =´d 2 p pP 1 ( r, p). We can then define the velocity field as v = g/ρ [58], implying thaṫ 2. Thermodynamically: We assume that the one-body distribution function P 1 takes the local equilibrium form and then define the velocity field to be the field v appearing in Eq. (50) [59]. Multiplying the dynamic equation for P 1 by p, inserting Eq. (50), and integrating over p also leads to Eq. (49).
One might question here whether the thermodynamic definition is sufficiently general as it appears to depend on the fluid being in a local equilibrium state. However, this second definition is actually very general, as can be seen from its application in the Mori-Zwanzig formalism [60][61][62][63][64]. When deriving the equations of hydrodynamics in the Mori-Zwanzig formalism (as done, e.g., in Ref. [59]), one introduces a relevant distributionP N that has the local equilibrium form [65] with the Hamiltonian H, the thermodynamic conjugates a i of the mean values a i of the relevant variables A i , and the number of relevant variables w. The Hamiltonian typically has the form H with the total interaction potential U . In fluid mechanics, we use as a relevant variable the total momentum density g = N i=1 p i δ( r − r i ) with the Dirac delta distribution δ. Inserting g into Eq. (51), writing v for g , and integrating over the phase-space coordinates of all except for one particle gives which is proportional to the result (50). Thus, the velocity field v is simply the thermodynamic conjugate for the momentum density [59]. If the relevant (local equilibrium) distribution (52) always equals the actual distribution, we have an ideal fluid, otherwise we have dissipation.
Since both definitions give the same result (49) in the passive case, their difference is usually not even mentioned. However, they are not equivalent for the active fluid considered here. While the mechanical route leads to Eq. (49) also in the active case, the thermodynamic route considered here gives the different result (2). From the Mori-Zwanzig point of view, our local equilibrium distribution (20) differs from the standard one (50) due to the presence of additional relevant degrees of freedom. An essential parameter for active phase separation in overdamped [51] and underdamped [25] active fluids is the average ofû · p (which corresponds to the average of the projection of the particle momentum onto the direction of self-propulsion). Motivated by this observation, we use the "momentum density polarization" g P = N i=1û i ⊗ p i δ( r − r i ) as a relevant variable in addition to g. Using the same steps as before, Eq. (51) then gives where : denotes a double tensor contraction. Hence, v P is the thermodynamic conjugate for g P . The form (53) is proportional to our local equilibrium form (20), as can be seen by inserting the expansion (29) into Eq. (20). When turning to a reduced description in terms of ρ and v only (namely AMI+), the deviation of (53) from (52) gives rise to a "viscosity term" in Eq. (3) with a viscosity proportional to v ld (ρ), consistent with the fact that deviations from the distribution (50) give rise to the standard viscosity in a passive fluid [54]. We are using the thermodynamic route in this work because it makes the microscopic derivation significantly easier as it is more natural from a statistical mechanics point of view. Whichever route we choose, it is a remarkable observation that they do not agree. This result is reminiscent of the one obtained for pressure in the case of AMB (where it is found that the mechanical and the thermodynamic definition no longer coincide for active fluids [12]).
The limit γ → 0 has to be taken in a very careful way here to avoid a blow-up resulting from the fact that the system is then active (implying energy influx) but not damped. This problem can be addressed by letting damping and activity go to zero in such a way that the damping term −γ v vanishes, while the active terms remain finite. In the passive limit, the parameters A 1 , A 3 , and A 4 will vanish since g has no angular dependence in this case. From the microscopic definitions given by Eqs. (106)-(110), it can be seen that we can avoid a blow-up by letting γ, v 0 , A 1 , A 3 , and A 4 go to zero at the same rate. In this case, we have More generally, we can -as will be discussed later -still derive an interesting quantum-mechanical model (resembling models used in Refs. [66,67]) for γ = 0. Thus, here, the assumption γ = 0 is primarily made to simplify the derivation and will be dropped later. Moreover, the static mapping used in Section II G can be made already for the overdamped AMB, such that our later analysis of the tunnel effect does not hang on the assumption γ → 0.
We define ρ q = 2ρ (ρ q will later be interpreted as the quantum-mechanical density) and assume that ρ q has only small deviations from a spatially and temporally constant reference density ρ 0 . Noting that adding a constant to ρ q has no influence on the dynamics once we have set f = 0, we can then approximately write (55) As Eq. (8) is linear in ρ, it is left unchanged by the replacement ρ → ρ q /2, i.e., it holds for ρ q in exactly the same way as for ρ. In Eq. (54), we replace ρ by ρ 0 ln(ρ q /ρ 0 )/2 (motivated by Eq. (55)) and assume The last step uses the identity [46] where ξ is a function, and the fact that ln(ξ) = 2 ln( √ ξ). Moreover, we set 2 2m = κρ 0 with the reduced Planck constant . We then arrive at the Madelung equations [45] Next, we assume that v is a potential flow such that we with a phase S, and that v satisfies the condition (see Ref. [68]) with a closed loop L and n ∈ Z. We can then substitute where ψ is (an analogue of) the wave function and i is the imaginary unit. Combining Eqs. (59)- (61) and (63) then finally yields the Schrödinger equation It is also interesting to see what happens if we do not set f (ρ q ) = 0. In this case, Eq.
Applying the substitutions (58), (61), and (63) then gives the nonlinear Schrödinger equation [69] i For example, if we set f = aρ q , we find [67] i which is the Gross-Pitaevskii equation. This equation has a wide range of applications, such as modeling Bose-Einstein condensates (see Ref. [67] and references therein).
In the considered active matter system, an even more realistic case would be to have 2 f = k B T ln(ρ q /2) (as suggested by the microscopic derivation in Section II B) and γ = 0. In this case, Eq. (56) (68) The Madelung transformation then leads to [66] i which is the Schrödinger-Langevin equation [70].
From a mathematical perspective, the fact that AMI allows to derive the Madelung equations is essentially a consequence of the fact that AMI is a compressible Euler equation with the pressure being given by the most general expression of a certain order in gradients and densities. Since the Madelung equations are of this order in gradients and densities (if we can approximate densities by their logarithm), they must be contained in AMI. The reason why we use AMI and not the more general AMI+ is that AMI+ contains (some form of) viscosity and leads to velocity fields with non-vanishing rotation, which does not make sense if we want to think of the velocity as the gradient of a phase. It has been shown in the overdamped case that coarse-graining the dynamics of quorum-sensing particles (which have a densitydepending swimming speed) leads to active model B rather than B+, i.e., the rotational term is not relevant for these particles [18]. This indicates that if one wants to study the "active quantum" effects discussed in the present paper experimentally or via particle-based simulations, quorum-sensing bacteria could be a more promising system than ABPs. What is particularly interesting is that we have found a classical system where the "quantum potential" from the Madelung equations gives, in a certain limit, an appropriate expression for the classical (nonequilibrium) chemical potential, and that this classical system is active. On the other hand, the derivation shows that AMI provides an approximate description of a quantum system if f = 0 and γ = 0 and if κ = −λρ 0 is chosen according to Eq. (58) (ρ 0 should then be the average probability density of the quantum-mechanical particle).

F. Analogy to dark matter
An important field of application for the Madelung equations is the study of dark matter [71]. Recently, there has been an increase of interest in so-called fuzzy dark matter (FDM), which consists of ultralight scalar 3 The expressions ln(ρq/2) and ln(ρq) differ only by an additive constant ln (2), which vanishes after applying ∇. Therefore, we can simply write ln(ρq) here.
particles. Research on FDM is motivated by the lack of evidence for other dark matter candidates, by the fact that such ultralight particles are predicted by various models from particle physics (such as string theory), and by the interesting pattern formation effects that it leads to. The FDM particles are mostly in the ground state and can thus be described by a single macroscopic wavefunction as in a Bose-Einstein condensate [72]. This wavefunction can (in the nonrelativistic limit without cosmic expansion) be described by the Madelung equations (59) and (60) with the gravitational constant G. Equation (70) determines U 1 , which is here the gravitational potential, via the density ρ q . Let us now consider the dynamics of a system of ABPs with an electric charge q. By the reasoning from Section II B, its dynamics would be given by AMI in the limit where v ld and ξ are small. The (electrostatic) potential U 1 could be calculated from the charge distribution ρ via the Poisson equation with the permittivity . As shown in Section II E, AMI contains the Madelung equations as a limiting case. Therefore, in the "quantum limit" considered in Section II E, an underdamped charged active matter system would be described by equations of the same form as a fuzzy dark matter system, suggesting an interesting parallel between active and astrophysical systems. Interestingly, also generalized Madelung equations with f = 0 and γ = 0 have been used in the context of dark matter physics [67], implying that the analogy between dark and active matter persists also in this more general and more realistic case. The analogy between dark and active matter is further supported by the fact that (as mentioned above) fuzzy dark matter models are based on Bose-Einstein condensates, which have been found also in active matter [1][2][3].
Note, however, that there is also an important difference, namely the fact that the density appears with a different sign in the gravitational Poisson equation (70) and the electrostatic Poisson equation (71). This is a consequence of the fact that gravity is a purely attractive force, whereas electrostatic forces are repulsive for particles of the same charge. Therefore, the patterns observed in fuzzy dark matter and in charged active systems might be quite different.

G. Tunnel effect
In this section, we restrict ourselves to one-dimensional systems.

Tunnel effect in quantum mechanics
Time-independent problems in quantum mechanics can be described by the stationary Schrödinger equation with the energy E. One of the most remarkable phenomena of quantum mechanics is the tunnel effect, where a particle travels through a potential barrier that it could not pass through classically. It can be described theoretically by solving Eq. (72) for the potential where V 0 is the height and 2L the width of the potential barrier. As is well known, the solution of Eq. (72) with the potential (73) is given by with the wavenumbers the transmission coefficients T 2 and T 3 , and the reflection coefficients R 1 and R 2 . (Explicit expressions for these coefficients are given in Ref. [74].) The physical interpretation of the solution (74) is that it describes the wavefunction of a particle with energy E that approaches a rectangular potential barrier of height V 0 > E from the left, decays within the barrier, and continues to travel as a wave on the right of the barrier. Since the particle has thereby passed through a barrier that it could not have passed through classically, this phenomenon is known as the "tunnel effect". Due to the linearity of Eq. (72), another solution is given by which is simply the superposition of the solution given by Eq. (74) and the same solution mirror reflected at x = 0 (corresponding to a particle coming from the right). Such a symmetric tunneling solution has advantages in a numerical treatment (as it allows to use periodic boundary conditions) and captures the same physics. The quantum-mechanical density ρ q = |ψ| 2 for the solution (77) is given by where we have written R 1 + T 3 = Re iα with the modulus R and the phase α of the complex number R 1 + T 3 .
Using the Madelung transformation, Eq. (72) can be rewritten as [74] The tunnel effect was studied in a Madelung framework in Ref. [74], where it was shown that, from a hydrodynamic point of view, it can be seen as arising from a pressure jump at the boundary of the barrier that balances the jump in the potential.

Tunnel effect in active matter
We now show how an analogon of the tunnel effect can be found in active matter. For simplicity, we set v = 0, such that AMI reduces to µ = const. with µ given by Eq. (11). Solutions of Eq. (11) with µ = const. are also stationary solutions of AMB, such that all of the following considerations apply to both AMI and AMB. Defining ρ = ρ 0ρ , µ = µ 0μ , f = µ 0f ,Ũ 1 = µ 0 U 1 , and x = x 0x (the tildes denote dimensionless quantities and ρ 0 , µ 0 , and x 0 are constants) gives which is simply Eq. (11) in dimensionless form. We could have eliminated the parameters κ and λ by an appropriate choice of the constants ρ 0 , µ 0 , and x 0 , but we assume here that these constants have already been used to eliminate other parameters, e.g., in the free energy. We now consider the special case with κ = −λ and f = 0, in which Eq. (82) reduces to A solution of Eq. (83) for the potential (73) is given by with the wavenumbers the phase shift and two constants A and B that satisfy Equations (87) and (88) ensure that ρ and ∂ x ρ are continuous at the boundaries of the potential barrier. We have thus found an analytical solution of Eq. (83) for an active system at a potential barrier, namely Eq. (84). At the boundary of the potential barrier, a discontinuity in κ∂ 2 x ρ balances the discontinuity in U and thereby ensures that µ is constant (as required for a stationary solution). Following the analysis in Ref. [74], we can define I = λ(∂ x ρ) 2 (which is here the active term) and a pressure Π = −ρκ∂ 2 x ρ. Π is the pressure one would get from the thermodynamic expression µρ − f for λ = 0 and f = 0 [12]. (Here, we have λ = 0, so Π is in general not equal to the thermodynamic or mechanical pressure in the active system.) With these definitions, Eq. (82) gives for f = 0 At the boundaries of the potential barrier, ρ and I are continuous. The "tunneling" is thus a consequence of a pressure jump ∆Π = ρV 0 at the boundaries.
For v = 0, the stationary Madelung equation (80) coincides with the stationary form of AMI or AMB given by Eq. (83) if we identify ρ = 1 2 ln(ρ q ) (cf. Eq. (55)), κ = 2 /(2m) (cf. Eq. (58)), and µ = E. Therefore, we do not even need to employ the approximation (55) from the dynamical case, we can just straightforwardly map the quantum onto the classical problem. 4 Taking the logarithm of the quantum solution (78) does indeed give us something that (apart from phases and prefactors) looks like Eq. (84), indicating that similar physical mechanisms act here. In particular, the change in the potential leads to a shift in the wavenumber from k to κ that gives a density decay within the barrier for µ < V 0 (or E < 0) both in the quantum and in the active case (compare Eqs. (75) and (76) to Eqs. (85) and (86)).
The strong mathematical analogy between AMI and the Madelung equations (or between AMB and the stationary Schrödinger equation) holds only for the rather special case f (ρ) = 0 and κ = −λ. In a real experiment, these equalities will, of course, be realized at most approximately. Therefore, it is interesting to investigate how robust the analogy between active matter and quantum mechanics is if these equalities are modified. For this purpose, we consider the more general model where U 1 is still given by Eq. (73). For a = 0 and κ = −λ, the analytical solution (84) is known. Starting from these parameter values and this solution, we can find solutions for Eq. (90) for other parameter values via numerical continuation (see Section IV B). We wish to ensure that the density ρ is always positive and that ρ and ∂ x ρ take identical values on both boundaries of the domain, allowing us to use periodic boundary conditions. This determines the one-dimensional domain Ω = [−α/k − L, L + α/k]. Furthermore, we set κ = 1, V 0 = 5, µ = 1, and α = π/4 and use the analytical solution (84) as starting solution for the continuation. Using Eqs. (85)-(87) we obtain L = arctanh(1/2)/2. Moreover, we set B = 0.5 (an arbitrary positive constant can be chosen here); A is then determined by Eq. (88). Note that with this also the averaged rescaled particle density is determined asρ =´Ωdx ρ(x)/|Ω| ≈ 0.7945 (where |Ω| is the domain length). It can be chosen arbitrarily using different values of B. Hence, the following result does not depend on the particle number. The starting state is now continued, changing various parameters while keepingρ fixed. This in turn determines µ as corresponding Lagrange multiplier. Alternatively, one could keep µ fix, in which caseρ would change during the continuation. However, this is not pursued here. Figure 1 shows bifurcation diagrams and solution profiles that illustrate the tunnel effect that can be observed in model (90). In the left and center columns of Fig. 1, we see how the L 2 -norm ´Ω dx ρ 2 (x)/|Ω| and the generalized chemical potential µ, respectively, depend on the parameters (top) a, (center) κ, and (bottom) λ. Finally, the right column of Fig. 1 presents solution profiles for the states indicated by orange and blue circles in the corresponding bifurcation diagrams to their left. The dashed black curve in each solution plot indicates the analytical solution given by Eq. (84) for comparison.
In general, the plots in the right column of Fig. 1 show that the general form of the solution does not change significantly if the parameters' values are not exactly those used for the analytical mapping. This indicates that the tunnel effect, and the general active-quantum analogy presented here, are not an artifact of picking the parameter values in such a way that it works, but rather a robust phenomenon that can be investigated also in microscopic simulations and experiments. Furthermore, according to a linear stability analysis that is performed during the continuation (see Section IV B), the solution is linearly stable with respect to perturbations compatible with mass conservation for all considered parameter values. Despite this limitation and the fact that we consider a small domain, the stability of all solutions emphasizes the relevance for experiments.
We can also get a more detailed idea of the effect that changing the various parameters has on the solution (84). In general, a steep decrease of ρ towards x = 0 indicates that the field cannot penetrate far into the potential barrier, whereas a more flat curve is a sign of a strong tunnel effect. Changing λ has no strong effect on the form of the solution (Fig. 1(c3)). The tunnel effect becomes more pronounced for positive values of a, whereas it is suppressed by negative ones (Fig. 1(c1)). Since positive values of a are more plausible on physical grounds (one would typically expand f around a local minimum rather than around a local maximum), we can expect this "tunneling" to be even more significant in real systems. Note that for sufficiently large values of a, we get µ > V 0 , such that strictly speaking we do not have tunneling anymore (since tunneling requires E < V 0 and µ corresponds to E). For µ > V 0 , κ becomes imaginary (see Eq. (86)) such that ρ has the form ln(cos(x)) also within the barrier. The strongest effect can be found by varying κ (Fig. 1(c2)). If it is small (close to zero), we observe a sharp decrease and thus very weak tunneling. For larger κ, on the other hand, the field can pass through the barrier much more easily. This result is plausible since, as indicated above, it is the discontinuity in the κ term that balances the discontinuity of the potential. Also, larger values of κ imply that gradients, which are smaller if the fluid passes through the barrier (i.e., if the tunnel effect is present), are associated with an energetic cost, implying that "tunneling" is more likely to occur for larger κ.

III. CONCLUSIONS
In this work, we have systematically derived an extension of common scalar active matter models to the underdamped case which we refer to as active model I+. This model and its derivation reveal some interesting and novel features of inertial active matter, such as the fact that mechanical and thermodynamic definitions of the velocity give different results and that the particles' density-dependent swimming speed acts as an effective viscosity. Moreover, we have shown that AMI+ contains (a nonlinear extension of) the Madelung equations and therefore the (nonlinear) Schrödinger equation as a special case, such that the Schrödinger equation can be seen as an active field theory. This allows to study quantum effects in active-matter systems, as has been demonstrated for the tunnel effect and for fuzzy dark matter. A numerical investigation of the active tunnel effect shows that this active-quantum analogy has no sensitive depen-dence on the assumptions that have been made to derive it, indicating that it is of broader relevance for both theory and experiment.

A. Quasi-stationary approximation
Here, we explain in more detail the microscopic derivation of Eq. (45). Using the quasi-stationary approximationv Eq. (44) gives where v ld is defined in Eq. (4). By inserting Eq. (92) recursively into itself and neglecting terms of higher than second order in gradients, of higher than first order in velocities, or that involve products of polarizations with velocities, we find The motivation behind these approximations is that we wish to derive a theory of third order in gradients and of second order in velocities and that we assume both polarizations and velocities to be small. (By "velocity", we mean v, whereas v P is always referred to as "velocity polarization".) The velocity polarization v P appears in Eq. (43) only in the term (v P · ∇) · v P /2 (quadratic in v P and of first order in gradients) and in the term D R P · v P /(2ρ) (product with the small polarization). If we insert Eq. (93) into Eq. (42) and drop again terms containing products of velocities with polarizations (in particular the advection term), we finḋ where we used the vector identity If we have v ld (ρ) ≈ v 0 , the first term on the right-hand side of Eq. (94) reduces to the self-propulsion term known from the active PFC model [35]. We now make the further quasi-stationary approximation˙ P = 0 and find Inserting Eq. (96) into itself and neglecting terms of higher than third order in gradients and second order in densities gives where we have used Eq. (95) again. This agrees with the result from Ref. [57] if we neglect terms of higher order in gradients and the velocity term in Eq. (97). We can insert Eq. (97) into Eq. (93) and neglect terms of higher than third order in gradients to get Eqs. (47) and (48). Inserting Eqs. (47), (48), and (97) into Eq. (41) and neglecting terms of higher than second order in gradients gives Eq. (2). The reason that terms of second order in gradients are sufficient is that all third-order terms would include also v, which (as is evident from Eq. (3)) is of at least first order in gradients.
Deriving Eq. (3) is slightly more involved. First, we deal with the term D R P · v P /(2ρ) appearing in Eq. (43). Inserting Eqs. (47), (48), and (97), dropping terms of higher than third order in gradients, terms quadratic in v that are of higher than first order in gradients (since v is of first order in gradients), terms of higher than second order in fields, and products of density gradients and velocities (these approximations will be referred to as "standard approximations" from here on) gives where we have used ( ∇ ⊗ ∇ρ) · ∇ρ = 1 2 ∇( ∇ρ) 2 .
We have not expanded the expression v ld (ρ) 2 in the last term of Eq. (98) to simplify the notation even though this term thereby contains terms up to third order in fields. The first term on the right-hand side of Eq. (98) can be rewritten using ( ∇ρ)/ρ = ∇ ln(ρ). In the fourth-fromlast, third-from-last, and penultimate terms, we replace ρ by 0 in the denominator such that these terms are of second order in ρ as required. This yields Next, we consider the term (v P · ∇) · v P /2. Inserting Eqs. (4), (47), and (48) gives with the standard approximations where we used Eq. (99) and ρ ∇ ∇ 2 ρ = ∇(ρ ∇ 2 ρ) − ( ∇ρ)( ∇ 2 ρ).
Finally, using Eqs. (97) and (102) and the standard approximations, we find We have not dropped the higher-order contributions in ρ for the logarithmic term, which is consistent with the fact that we do not make a constant-mobility approximation [75]. It is interesting that, instead of the thermal energy k B T we would have in the passive case, the ideal gas contribution f is proportional to k B T + mv 2 0 /2, implying that the active contribution to the kinetic energy effectively shifts the temperature by mv 2 0 /(2k B ). This is a different sort of "effective temperature" than the one reported for active systems in Refs. [41,53]. The additional term mv 2 0 ln(ρ)/2 originates from Eq. (100). Equation (105) has the form of Eq. (45), such that we can read off the expressions f (ρ) = k B T + mv 2 0 2 ln(ρ)

B. Numerical path continuation
In Section II G, we use numerical path continuation via the Matlab package pde2path [76]. Starting from the analytical solution (Eq. (84)) of model (90), pde2path subsequently applies tangent predictors and Newton correctors to track a branch of steady states through parameter space. A numerical linear stability analysis during the continuation yields the stability of the corresponding solution and enables the detection of bifurcations. pde2path uses the finite element method and the model is implemented in a weak formulation. We have used a primary control parameter (a, κ, or λ) and the chemical potential µ as a secondary one which is adapted freely during the continuation to ensure mass conservation.

DATA AVAILABILITY
The raw data corresponding to the figure shown in this article are available as Supplementary Material [77].

CONFLICTS OF INTEREST
There are no conflicts of interest to declare.