Active particles under confinement and effective force generation among surfaces

TL;DR

This paper investigates how geometric confinement affects the behavior of a one-dimensional active suspension with thermal noise, revealing the formation of polar order and a Casimir-like force due to active forces and thermal fluctuations.

Contribution

It combines numerical and analytical methods to analyze active suspension behavior under confinement, highlighting the emergence of polar order and active force-induced Casimir-like forces.

Findings

Active suspension forms a polar ordered layer near walls.

Thermal noise and active forces influence system structure.

Active forces generate a Casimir-like attractive force.

Abstract

We consider the effect of geometric confinement on the steady-state properties of a one-dimensional active suspension subject to thermal noise. The random active force is modeled by an Ornstein-Uhlenbeck process and the system is studied both numerically, by integrating the Langevin governing equations, and analytically by solving the associated Fokker-Planck equation under suitable approximations. The comparison between the two approaches displays a fairly good agreement and in particular, we show that the Fokker-Planck approach can predict the structure of the system both in the wall region and in the bulk-like region where the surface forces are negligible. The simultaneous presence of thermal noise and active forces determines the formation of a layer, extending from the walls towards the bulk, where the system exhibits polar order. We relate the presence of such ordering to the mechanical pressure exerted on the container's walls and show how it depends on the separation of the boundaries and determines a Casimir-like attractive force mediated by the active suspension.

Equations

$$\gamma \dot{x}=-{\varphi }^{\prime }(x)+\gamma u+\gamma \sqrt{2D_t}\xi ,$$

$$\tau \dot{u}=-u+\sqrt{2D_a}\eta ,$$

$$\dot{x}=v+\sqrt{2D_t}\xi$$

$$\displaystyle\dot{v}=-\frac{1}{\tau}\Bigl{[}v+\frac{\phi^{\prime}(x)}{\gamma}-\sqrt{2D_{a}}\eta\Bigr{]}-\frac{\phi^{\prime\prime}(x)}{\gamma}v-\frac{\sqrt{2D_{t}}}{\gamma}\phi^{\prime\prime}(x)\xi.$$

$$v\frac{\partial f(x,v)}{\partial x}-D_t\frac{{\partial }^2}{\partial x^2}f(x,v)$$

$$-\frac{1}{\tau }\frac{\partial }{\partial v}\left(\frac{D_a}{\tau }\frac{\partial }{\partial v}+\left(1+\left(\theta (-x)+\theta (x-2L)\right)\frac{\tau k}{\gamma }\right)v\right)f(x,v)=$$

$$\left(\theta (-x)\frac{kx}{\tau \gamma }+\theta (x-2L)\frac{k(x-2L)}{\tau \gamma }\right)\frac{\partial f(x,v)}{\partial v}$$

$$+\left(\theta (-x)+\theta (x-2L)\right)\frac{D_t}{\gamma }k\left(\frac{k}{\gamma }\frac{{\partial }^2}{\partial v^2}f(x,v)-2\frac{{\partial }^2}{\partial v\partial x}f(x,v)\right),$$

$$p_h(x,u)=\mathcal{N}{e}^{-\left(\frac{k{x}^2}{2D_a\gamma }\frac{1}{\frac{1}{\Gamma }+\Delta }\right)}{e}^{-\left(\frac{\tau \Gamma }{2D_a}{\left(u-\frac{k}{\gamma }\frac{x}{1+\Delta }\Gamma \right)}^2\right)},$$

$$p_l(x,u)\approx p_h(x,u)\theta (-u)\theta (-x)$$

$$n_0^{left}(x)=n_w{e}^{-\left(\frac{k{x}^2}{2D_a\gamma }\frac{1}{\frac{1}{\Gamma }+\Delta }\right)}\text{erfc}\left(\sqrt{\frac{\tau \Gamma }{2D_a}}\frac{1}{(1+\Gamma \Delta )}\frac{k}{\gamma }x\right)\theta (-x),$$

$$U_{eff}(x)$$

$$\approx \frac{1}{(\frac{1}{\Gamma }+\Delta )}\left(\frac{k}{2D_a\gamma }{x}^2-\frac{2}{\sqrt{\pi }}\sqrt{\frac{\tau }{2\Gamma D_a}}\frac{k}{\gamma }|x|\right).$$

$$x_p=-\sqrt{\frac{2}{\pi \Gamma }{D}_a\tau },$$

$$f_c(x,v)=(\frac{\tau }{2\pi D_a}{)}^{1/2}\sum _{\nu \ge 0}{n}_{\nu }(x)h_{\nu }(v)\exp (-\frac{\tau }{2D_a}{v}^2),$$

$$h_{\nu }(v)=(-1{)}^{\nu }(\frac{D_a}{\tau }{)}^{\nu /2}\exp (\frac{\tau }{2D_a}{v}^2)\frac{d^{\nu }}{d{v}^{\nu }}\exp (-\frac{\tau }{2D_a}{v}^2).$$

$$\frac{\partial n_{\nu -2}(x)}{\partial x}+\nu \frac{D_a}{\tau }\frac{\partial n_{\nu }(x)}{\partial x}=-\frac{\nu -1}{\tau }{n}_{\nu -1}(x)+D_t\frac{{\partial }^2}{\partial x^2}{n}_{\nu -1}(x),$$

$$m(x)=\int dvv{f}_c(x,v)=\frac{D_a}{\tau }{n}_1(x).$$

$$\frac{\partial m(x)}{\partial x}=D_t\frac{{\partial }^2}{\partial x^2}{n}_0(x),$$

$$\frac{1}{\tau }(1+\frac{D_a}{D_t})m(x)-D_t\frac{{\partial }^2}{\partial x^2}m(x)+2\frac{D_a^2}{{\tau }^2}\frac{\partial n_2(x)}{\partial x}=0,$$

$$\frac{\tau }{D_a}\frac{\partial m(x)}{\partial x}+\frac{2}{\tau }{n}_2(x)-D_t\frac{{\partial }^2}{\partial x^2}{n}_2(x)+3\frac{D_a}{\tau }\frac{\partial n_3(x)}{\partial x}=0$$

$$n_0(x)=(n_w-n_M)\frac{\cosh (\frac{x-L}{\lambda })}{\cosh (\frac{L}{\lambda })}+n_M,$$

$$m(x)=m_w\frac{\sinh (\frac{x-L}{\lambda })}{\cosh (\frac{L}{\lambda })}=\frac{D_t}{\lambda }(n_w-n_M)\frac{\sinh (\frac{x-L}{\lambda })}{\cosh (\frac{L}{\lambda })}$$

$${\lambda }^2=D_t\tau \frac{D_t}{D_a+D_t}$$

$$\frac{\partial n_0(x)}{\partial x}+2\frac{D_a}{\tau }\frac{\partial n_2(x)}{\partial x}=0.$$

$$⟨v^5{⟩}_x\sim 10⟨v^2{⟩}_x⟨v^3{⟩}_x,$$

$${\mathcal{P}}_{wall}=-N{\int }_{-\infty }^{0}dx{n}_0(x){\varphi }^{\prime }(x),$$

$${\mathcal{P}}_{wall}=N{n}_w{D}_a\gamma (\frac{1}{\Gamma }+\Delta )s(\Gamma ,\Delta ),$$

$$s(\Gamma ,\Delta )=\frac{1}{1+\frac{\Gamma -1}{1+\Gamma \Delta }\left[1+\sqrt{1+\frac{1+\Gamma \Delta }{\Gamma -1}}\right]}\le 1,$$

$$\frac{D_a}{\tau }{n}_1(x)=D_t\frac{\partial }{\partial x}{n}_0(x),$$

$$\frac{\partial }{\partial x}\left(\left(\frac{D_a}{\tau }+\frac{D_t}{\tau }\right)n_0(x)+2\frac{D_a^2}{{\tau }^2}{n}_2(x)-D_t^2\frac{{\partial }^2}{\partial x^2}{n}_0(x)-\frac{D_t}{\gamma }{\varphi }^{\prime \prime }(x)n_0(x)\right)$$

$$=-\frac{{\varphi }^{\prime }(x)}{\tau \gamma }{n}_0(x).$$

$$\displaystyle\tau\gamma\,\Bigl{(}\frac{(D_{a}+D_{t})}{\tau}n_{0}(x)+2\frac{D_{a}^{2}}{\tau^{2}}n_{2}(x)-D_{t}^{2}\,\frac{\partial^{2}}{\partial x^{2}}n_{0}(x)-\frac{D_{t}k}{\gamma}\theta(-x)n_{0}(x)\Bigr{)}_{-\infty}^{\bar{x}}$$

$$=\frac{{\mathcal{P}}_{wall}}{N},$$

Full text

I Introduction

Self-propelled particles, motile organisms such as bacteria, and artificial micro-swimmers display a characteristic tendency to aggregate, a phenomenon which is currently the object of vivid interest among physicists and biologists Ramaswamy (2010); Romanczuk et al. (2012); Marchetti et al. (2013); Saintillan and Shelley (2015). Unlike standard molecular systems, where aggregation is induced by attractive forces and/or entropic interactions Likos (2001) due to volume exclusion, active particles may spontaneously produce regions of higher density because their dynamical properties change if they interact with other particles. These phenomena have been investigated experimentally, by numerical simulation and theoretically Bechinger et al. (2016); Stenhammar et al. (2014) and led in the case of self-propelled particles endowed with only repulsive inter-particle interactions to the concept of motility induced phase separation (MIPS) Cates and Tailleur (2015) analogous to the liquid-gas coexistence in standard liquids. On the other hand, the accumulation phenomenon in the proximity of a purely repulsive confining wall, i.e. the aggregation with an external object, occurs even when active particles are not subject to mutual interactions Lee (2013). Such a behavior is of great practical importance since experiments are often conducted on systems where the size of the experimental apparatus could be of an order of magnitude comparable to the persistence length of the active particles, which is the typical distance over which particle’s orientation persists. The explanation of the underlying mechanism attributes the accumulation to the reduction of the particles’ mobility in the presence of the walls and is captured by some existing theories Ezhilan et al. (2015); Maggi et al. (2015); Nikola et al. (2016); Fily et al. (2014); Smallenburg and Löwen (2015); Marconi et al. (2017).

Confining surfaces besides triggering particle accumulation in a very thin adjacent region, may also create a diffuse layer where neither the density is constant nor polar order field vanishes as in bulk systems. In a series of recent articles Brady and coworkers Yan and Brady (2015, 2018) have thoroughly investigated such an inhomogeneous layer by means of a mesoscopic approach where these effects were captured by a simple system of differential equations and the action of the walls was taken into account by prescribing the appropriate boundary conditions.

In the present article, we consider the effect of a confining potential, $\phi(x)$, varying along the single $x$-direction, on the steady-state behavior of an assembly of non-interacting self-propelled particles described by means of the so-called active Ornstein-Uhlenbeck particle (AOUP) model Fily and Marchetti (2012); Szamel (2014) The AOUP is driven by an active force of variable intensity and direction assimilated to a Gaussian colored noise process Hanggi and Jung (1995) sharing the same exponential two-time correlator as the active force of the ABP model.

The characteristic time, $\tau$, of the process represents the average persistence of the trajectory along a given spatial direction, i.e. the crossover time from a ballistic to a diffusive behavior. In fact, in both models the mean square displacement evolves ballistically at short time and diffusively at long times with an effective diffusion coefficient given by the sum of an active contribution, $D_{a}$ plus a thermal contribution, $D_{t}$, stemming from the microscopic collisions with the solvent molecules. Regarding the difference between the two models, in ABP the absolute value of the active speed is constant, whereas in AOUP each component independently fluctuates according to a Gaussian distribution. Our choice to use the AOUP, instead of other popular alternatives such as the ABP and the Run and Tumble Tailleur and Cates (2008) models is motivated by the possibility of applying straightaway methods and concepts similar to those employed in the study of the Kramers equation Titulaer (1978). In the last few years, approximate treatments of the AOUP, such as the Fox method of Ref. Farage et al. (2015); Wittmann and et al (2017) and the so-called unified colored noise approximation (UCNA) Marconi and Maggi (2015), have been developed: by using an adiabatic approximation, which is tantamount to impose a detailed balance condition Cates (2012), the UCNA allowed making reasonably accurate predictions about the steady-state properties of a rather general class of active systems Marconi et al. (2016); Marini Bettolo Marconi et al. (2017). Nevertheless, the UCNA results regarding the structure of active suspensions in the proximity of a confining surface disagree with those obtained by mesoscopic treatments of the ABP model, where the container wall is treated as an infinitely sharp interface and a set of boundary conditions on the density and polar fields are imposed on it. Such a discrepancy is more severe when the finite diffusivity of the solvent, $D_{t}$, is not negligible and a polar order appears close to the surface. In the present theory, we go beyond the UCNA and do not impose the detailed balance condition in deriving the form of the steady-state solution. In contrast with mesoscopic approaches  Ezhilan et al. (2015); Yan and Brady (2015); Duzgun and Selinger (2018), we treat the wall and bulk regions on equal footing and instead of considering the wall as a sharp boundary Fily et al. (2017) we study the distribution function in each region thus providing a microscopic description of why and how particles accumulate at the boundaries and form a diffuse layer near it.

At variance with the ABP model which is well defined only for two or more dimensions, the AOUP model can be implemented also in one-dimension. In a system with a simple geometry, such as infinite parallel plates - a situation which can be realized assuming periodic boundary conditions - the coordinates parallel to the plate play just a minor role, as a constant factor in the definition of the control parameters. In practice, provided we restrict to a region far enough from the edges of the plates the present treatment applies also to the case of finite plates and the one-dimensional description is valid.

Besides clarifying the mechanism causing the enhancement of the density and polar order near the walls in self-propelled systems, we discuss the role of the activity in determining the forces that the particles exert on inclusions, a topic of current interest. In fact, several groups by numerical simulation of RnT Ray et al. (2014), ABP Ni et al. (2015) and swimmer suspensions Parra-Rojas and Soto (2014), have recently reported evidence about effective interactions arising between two plates placed in a bath of active particles, a phenomenon akin to the Casimir-like attractive force Casimir (1948) observed in the presence of non-equilibrium diffusive dynamics Aminov et al. (2015); Brito et al. (2007).

The paper is organized as follows: in Section II, we introduce the model of confined active particles and consider a truncated parabolic confining potential Solon et al. (2015); Sandford et al. (2017). In Section III, we illustrate the numerical method and study the model by numerical simulations, in Section IV we present our theory which goes beyond the UCNA concerning an important aspect: for a system of non-interacting particles the UCNA predicts that the distribution function has a local dependence on the potential, thus if the potential and its derivatives vanish in some region of space the distribution is uniform. By using a hierarchy of equations for the velocity-moments of the phase-space distribution we are able to describe non-local effects and obtain predictions which are in better agreement with the numerical simulation results. In section V, we discuss the pressure exerted on the walls by the active suspension using the results of Section IV. As an application we derive a new expression for the effective force between two plates induced by the activity when the molecular diffusion is finite. Finally in Sec. VI we present the conclusions.

II Model

The model consists of a system of $N$ non-interacting active particles suspended in a fluid, driven by an active force $\gamma u$, where $\gamma$ is a Stokes friction constant, and subjected to an external potential $\phi(x)$ and to a random force representing the effect of the collisions with the molecules of the fluid. The self-propulsion force, originating from an internal mechanism and fluctuates both in intensity and direction and is modeled by a colored noise term, $u(t)$, evolving according to an Ornstein-Uhlenbeck process of correlation time, $\tau$. The resulting governing equations read:

[TABLE]

where $\xi$ and $\eta$ are two independent white noises with unitary variance and zero average and $(D_{t},\gamma)$ and $(D_{a},\tau)$ refer to the interactions with the solvent and with the active bath, respectively. The term $\gamma u$ represents the self-propulsion mechanism, an internal degree of freedom converting energy into motion and has the following self-correlation $\langle u(t)u(t^{\prime})\rangle=\frac{D_{a}}{\tau}\exp(-|t-t^{\prime}|/\tau)$, with variance $D_{a}/\tau$ identified with the active power. It is well-known that this system is out of equilibrium whenever $\tau>0$ Marconi et al. (2017) and that in the limit $\dot{u}\approx 0$ Eq. (2) reduces to $u\approx\sqrt{2D_{a}}\eta$, a Wiener process so that Eq. (1) describes a Brownian passive particle where the term $\gamma u$ merely produces an extra contribution to the diffusion.

We, now, generalize to the case $D_{t}\neq 0$ the change of variable of ref Marconi et al. (2016) which allows a hydrodynamic study of the model: we define the new variable $v=\dot{x}-\sqrt{2D_{t}}\xi$ and replace Eqs. (1) and (2) by the following set of equations:

[TABLE]

Thus the AOUP dynamics Eqs. (1)-Eq. (2) has been mapped onto the underdamped dynamics of a fictitious Brownian particle of position $x$ and velocity $v$ and effective mass $\mu=\tau\gamma$ evolving with a space dependent Stokes force and experiencing a delta-correlated thermal noise acting additively on the $x$ and multiplicatively on the velocity Caprini et al. (2018). Given the presence of multiplicative noise terms we use the Stratonovich interpretation of the stochastic differential equation Hanggi and Jung (1995).

For mathematical convenience, we shall restrict our study to the case where the effect of the confining walls is represented by two repulsive truncated harmonic potentials, $\theta(-x)kx^{2}/2$ and $\theta(x-2L)k(x-2L)^{2}/2$ ($\theta$ being the Heaviside function) acting only in the regions $(-\infty,0)$ and $(2L,\infty)$, whereas the central region ($0,2L)$ is a force-free region. The harmonic force is proportional to $k$, modeling the penetrability of the wall: since the range and strength of the force are both finite it could describe an elastic membrane of stiffness $k$ allowing the particles to explore the regions $x\leq 0$ and $x\geq 2L$. On the other hand, if the spatial resolution of the experimental device is low or the penetrability of the wall is small the use of a sharp interface model, obtained by imposing no-flux boundary conditions to prevent particle crossings as in ref. Yan and Brady (2015), is well justified.

Finally, with the aim of developing the theoretical methods of Sec. IV we introduce the stationary Fokker-Planck equation (FPE) Risken (1984) for the phase-space distribution $f(x,v)$ providing an equivalent statistical description of the system (3)-(4). We substitute $\phi^{\prime}(x)=k[x\theta(-x)+(x-2L)\theta(x-2L)]$ and $\phi^{\prime\prime}(x)=k[\theta(-x)+\theta(x-2L)]$ and obtain the following equation:

[TABLE]

where $\gamma/k$ is the characteristic time of the $x$-process.

III Numerical methods and results

In our numerical simulations equations (1)-(2) have been integrated by using the Euler-Maruyama algorithm neglecting terms of order $dt^{5/2}$, where $dt$ is the time-step size of the numerical integration Mannella and Palleschi (1989). Except where noted, each simulation has been run until time $2\cdot 10^{3}\tau$, with $dt\sim O(10^{-4}-10^{-6})$, depending on the $\tau$ values. The observables, such as the probability distribution and its momenta, have been computed by using both time and ensemble averages: we usually perform simulations with $N=10^{4}$ particles, waiting for a transient time $10^{3}\tau$, in such a way the system is fully thermalized.

In Fig. 1 we display the density profile obtained by numerical simulation in the case of two walls separated by a distance $2L=8$ keeping constant the ratio $D_{a}=1$ and $\tau=1$ and varying the intensity of the thermal noise $D_{t}$ as shown in the legend. One can observe that the density profile, $n_{0}(x)$, is continuous for all values of $D_{t}$ including the value $D_{t}=0$, at variance with the UCNA which predicts a finite jump at $x=0$ and $x=2L$ Fily et al. (2017). The effect of increasing the thermal diffusion, $D_{t}$, is to broaden the distribution with respect to the case $D_{t}=0$ and is best appreciated in the inset which shows that the profiles corresponding to the larger values of $D_{t}$ have slower decay. Such a scenario is similar to the one observed in the ABP model Yan and Brady (2015), where the presence of thermal noise has two consequences: a) it reduces the accumulation near the walls and b) it determines an exponential decay of the density profile in the force-free region and an associated screening length, $\lambda$, roughly dependent on the ratio $\Delta=D_{t}/D_{a}$. It is also interesting to see that the accumulation phenomenon near a repulsive wall, a typical non-equilibrium effect, survives upon the addition of thermal noise and disappears only in the limit $D_{t}\gg D_{a}$ when the particles behave as passive ones. We may conclude that $D_{t}$ has a double role in the potential region: on one hand, reduces the accumulation, decreasing the height of the peak; on the other hand, it favors the dispersion for $x>0$ as shown in the inset of Fig. 1.

In Fig. 2 we analyze the system with $D_{t}=0$ discussing for the sake of simplicity just the left wall: the particles accumulate approximately in the region close to $x=0$ and their profile does not have a Brownian counterpart: in fact, in the Brownian limit $\tau\to 0$, they would be uniformly distributed between [math] and $2L$ and depleted within the repulsive regions according to the Boltzmann weight $\propto\exp(-\phi(x)/T)$ at a uniform temperature $T=(D_{a}+D_{t})\gamma$. If $\tau>0$, the accumulation can be understood by considering the Eq.(5) (with $D_{t}=0$): the Stokes force is discontinuous, being $\gamma\Gamma v$ for $x<0$, with $\Gamma=1+\tau k/\gamma$, and $\gamma v$ for $x>0$. Hence, on one hand particles slow down in the regions $[-\infty,0]$ and on the other hand the repulsive wall pushes the particles towards the edge $x=0$. The interplay between the slow-down and the repulsion determines the observed accumulation. In the left panel of Fig. 2 each line corresponds to a system with $\tau=10$ and $D_{a}$ ranging between $1$ and $100$. Notice that the peak broadens and shifts towards more negative values of the x-coordinate with increasing active power $D_{a}/\tau$. The right panel of Fig. 2 shows that the location of the peak does not change if $D_{a}/\tau$ remains constant but its height increases when $\tau$ increases: indeed, a larger $\tau$ corresponds to a longer time spent by the particle in front of the wall and has no influence on the peak dispersion.

IV Theoretical treatment

Now, we present a theoretical analysis of the confined active system using a velocity-moment expansion to derive approximate solutions of the FPE (5) and compare the theoretical predictions with the numerical solutions of Eqs. (1)-(2). To this purpose it is mathematically convenient to study separately the two boundary regions characterized by a finite value of the external field from the central force-free region and write the stationary distribution associated with the FPE (5) as $f(x,v)=\theta(-x)f_{l}(x,v)+\theta(x)\theta(L-x)f_{c}(x,v)+\theta(x-L)f_{r}(x,v)$, where $f_{l}$ and $f_{r}$ are the distribution functions in the left and right regions, respectively, while $f_{c}$ is the distribution in the central region.

IV.1 Density profile in the wall region

In order to derive a theoretical expression for the probability density $n_{0}(x)$ in agreement with the numerical results above illustrated, we first consider the region $x<0$ and leave the treatment of $n_{0}(x)$ in the central region, $0\leq x\leq L$, until Sec. IV.2. Neglecting the truncated shape at $x=0$ of the potential $\frac{kx^{2}}{2}\theta(-x)$ and using the $(x,u)$ representation, we begin by approximating the probability distribution $p_{l}(x,u)$ by the stationary distribution, $p_{h}(x,u)$ of an AOUP confined to a symmetric harmonic trap $(\frac{kx^{2}}{2})$ given by  Das et al. (2018):

[TABLE]

${\cal N}$ being a normalization factor and $\Delta=D_{t}/D_{a}$. A similar expression for $f_{h}(x,v)$ is reported in Eq. (34) if we employ the $(x,v)$ representation. The reduced probability distribution $n_{0}(x)=\int dup_{h}(x,u)$ computed using formula (6) shows a poor agreement with the numerical results of Fig. 2, since it is a Gaussian centered at $x=0$, in contrast with the numerical evidence illustrated in Figs. 1 and 2 where such a peak is shifted towards negative x-values. The domain of validity of the harmonic approximation (6) can also be tested by comparing the simulation results for the conditional probability distribution function, $p(x|u)=p(x,u)/p(u)$ (where $p(u)=\int dxp(x,u)$) with the corresponding quantity obtained from the theoretical expression Eq. (6). In the left panel of Fig. 3, one observes significant deviations between the two curves when $x>0$.

Instead, in the right panel of Fig 3, we display the comparison between the theoretical (computed from the Gaussian formula (6)) and numerical x-variance, $\int dxx^{2}p(x|u)$: one can see that the first decreases as $u$ increases, whereas the latter remains nearly constant. The departure from the Gaussian prediction, becomes more and more relevant when $u$ increases, while it is negligible for negative values of $u$.

The comparisons shown in Fig. 3 , indicate that the Gaussian formula (6) agrees quite reasonably with the numerical results for the conditional probabilities $p(x|u)$ and $p(u|x)$, when $u<0$ and $x<0$, respectively, but the same formula does not provide an adequate prediction for the density profile $n_{0}(x)$ for $x<0$. Based on these evidences we improve the Gaussian approximation (6) by modifying the left wall $(x,u)$ distribution in the following way:

[TABLE]

and an analogous expression for right wall distribution, i.e. $p_{r}(x,u)\approx p_{h}(x,u)\theta(u)\theta(x-L)$. The rationale for such an assumption is the following remark: the wall can be regarded as a very massive body having zero speed and the active particle as a moving object with self-propulsion force $\gamma u$. A collision between the left wall and the active particle occurs only if $u$ assumes negative values and $x<0$. These two conditions are encapsulated in formula (7) and we, now, use it to derive an expression for $n_{0}(x)$ in the wall region $x\leq 0$ by integrating with respect to $u$. The resulting density, $n_{0}(x)$, in the region $x\leq 0$, reads:

[TABLE]

where $n_{w}$ is the density at $x=0$ and $x=2L$. Formula (8) shows a fairly good agreement with the simulation data both for $D_{t}=0$ and for $D_{t}>0$, as the left panel of Fig. 4 reveals.

Eq. (8), which describes the space-density in the potential regions, generalizes the result of Yan and Brady (2015) to a soft-wall, modeled as an external truncated potential. Moreover, we overcome the unphysical results of the UCNA approximation at $D_{t}=0$ and for hard walls, i.e. a discontinuous space density discussed in ref.Fily et al. (2017).

In order to assess the ansatz $p_{l}(x,u)=p_{h}(x,u)\theta(-u)\theta(-x)$ , we analyzed the numerical joint probability distribution $p(x|u)$, at fixed $u$ and estimated how important is the neglected contribution due to the population with $u\geq 0$. The right panel of Fig. 4 shows that in the region $x\leq 0$ the population characterized by a positive sign of the active force $p(x,u>0)$ represents only a small contribution to the density $n_{0}(x)$ for $x\leq 0$, thus roughly validating the approximation leading to Formula (8). One can observe that the $p(x,u\geq 0)$ decreases much faster than $p(x,u\leq 0)$ as $x$ becomes more negative. Interestingly, a similar scenario was reported by Widder and Titulaer for a related model Widder and Titulaer (1989): these authors studied the distribution function in the presence of a partially absorbing wall with specular reflection and found that $p(x,u)$ at the boundary $x=0$ was peaked at negative values of $u$ and rapidly decreasing towards zero for positive $u$.

Let us remark that the argument of the complementary error function in Eq. (8) is proportional to the ratio between the wall force $-kx$ and the average absolute value of the active force, $\gamma\sqrt{D_{a}/\tau}$. On the other hand, if $\tau=0$ there is no shift and as we shall see below the accumulation phenomenon is completely suppressed and on the contrary, one observes a depletion of the density controlled by the standard Maxwell-Boltzmann weight. It is possible to define an effective potential as $U_{eff}(x)=-\ln n_{0}^{left}(x)$ and obtain:

[TABLE]

For small $|x|$ and not too stiff walls ($k\gg 1$ and/or $D_{a}/\tau\ll 1$) we find that the effective force vanishes when

[TABLE]

whereas for strong walls $x_{p}\approx 0$. The position of the peak of the distribution does not depend on the ratio $\Delta$ of the two diffusion coefficients, but its width does. The position, $x_{p}$, of the peak gives a measure of the stiffness of the wall and we expect that for $k$ large enough the wall is quite impenetrable and $x_{p}\approx 0$, while for smaller values of $k$ we have $x_{p}<0$, a situation describing a soft wall or a floppy membrane Rodenburg et al. (2017).

IV.2 Central region

As illustrated in the previous section, the predictions for the density profile relative to the potential regions are in good agreement with the numerical data. Nevertheless, in the force-free regions (both for $D_{t}=0$ and $D_{t}\neq 0$) one observes a phenomenology which cannot be captured by simply setting to zero the external potential in Eq. (6), with the result of producing a constant $n_{0}(x)$. In particular, the numerical $n_{0}(x)$ displays a smooth decay from the wall value towards the value at midpoint $x=L$. Hereafter, we develop a hydrodynamic approach in order to find an approximation scheme for $n_{0}(x)$ and for this reason we consider appropriate to switch again to the $(x,v)$ representation of the distribution function. Since for a uniform system $f(x,v)$ is a Maxwell-Boltzmann distribution, in the force-free region we may expect to find a good approximation by expanding $f(x,v)$ in Hermite functions of the velocity and taking into account only the first few terms. To this purpose, we seek for an approximate solution of the FPE (5) in the central region by employing the following Hermite expansion:

[TABLE]

where the Hermite polynomials are

[TABLE]

By substituting the expansion in the FPE  (5) when $k=0$ we obtain the following recursion relation for the amplitudes $n_{\nu}$:

[TABLE]

with the condition $n_{-1}=0$ and $n_{0}(x)=\int dvf_{c}(x,v)$. We can, now, define the steady-state average polarization $m(x)$ as the first velocity moment of the distribution function:

[TABLE]

Under stationary conditions and $D_{t}=0$, $m(x)$ is proportional to the local average of the active force (being $u=v$ in this region) and thus vanishes in the absence of external fields in virtue of Eq. (1); this is seen by considering Eq. (11) for $\nu=1$ together with (12) :

[TABLE]

and assuming a zero current condition. On the other hand, if $D_{t}>0$, the fact that the external force is zero in some region of space does not necessarily imply the corresponding vanishing of $m(x)$. To show that, let us use again Eq. (13) and consider a finite density gradient term, $D_{t}\frac{\partial n_{0}(x)}{\partial x}$: a non uniform density profile is now sufficient to induce a polarization even where the external force acting on the AOUP is locally zero. As we shall see in section V, such a coupling between standard diffusion and polar order, represented by the $n_{0}(x)$ and $m(x)$ terms, respectively, determines an effective force between inclusions immersed in active suspensions.

Since Eq. (11) represents an infinite hierarchy of coupled differential equations for the coefficients $n_{\nu}$, we need to introduce a suitable truncation able to capture the phenomenology discussed in the previous Sections. Our truncation procedure comes easily in the Hermite-basis and consists in setting $n_{\nu}(x)=0$ for all $\nu\geq n_{max}$, i.e. in assuming corrections around a Gaussian-like approximation. The simplest possibility is to set $\nu_{max}=2$ (case A), which leads to the so called screening-approximation. Instead, by considering $\nu_{max}=4$ (case B), the resulting approximation is equivalent to the hydrodynamic treatment based on the first three moments of the velocity distribution Huang (2009) together with the idea that the term $n_{3}(x)$ (analogous to the heat flux) can be eliminated in favor of the spatial gradient of $n_{2}(x)$, i.e. to the gradient of a kinetic temperature.

Thus we write the following equations:

[TABLE]

[TABLE]

where we used eq. (13) to eliminate $n_{0}(x)$. In case A the closure is $n_{2}(x)=0$, while in case B, in analogy with the phenomenological procedure followed in hydrodynamic treatments, one assumes $n_{3}(x)=-\frac{\tau}{3}\frac{\partial n_{2}(x)}{\partial x}$. Both approximations predict exponential solutions and for $D_{a}/D_{t}$ small enough a typical length over which the moments vary scaling proportionally to $\sqrt{D_{t}\tau\,\frac{D_{t}}{D_{a}+D_{t}}}$. In the following, for the sake of simplicity, we shall report only results concerning the so called screening approximation (case A) first employed in reference Yan and Brady (2015) in the framework of the ABP model. The solution reads

[TABLE]

where $n_{w}=n_{0}(0)=n_{0}(2L)$ and $n_{M}$ depends on the geometry of the problem. The "polarization field" is given by:

[TABLE]

where

[TABLE]

and $m_{w}$ is the polarization field at the wall.

Notice that if $D_{t}=0$ the polarization field vanishes. The comparison in Fig.5 between the numerical $n_{0}(x)$ and the analytic prediction displays a fair agreement if $D_{t}$ is not too small with respect to $D_{a}$.

This result is consistent with the one obtained in ref. Yan and Brady (2015) for the ABP model, confirming that AOUP is a useful approximation of ABP which captures all the physical aspect of the accumulation near the walls also in the presence of thermal noise. In addition, our study sheds some light on the closure employed in the hydrodynamic-like approach, by considering the description in terms of the particle velocity.

As a consequence of the simultaneous presence of two baths (active and thermal), $\langle v\rangle_{x}$ is non-zero as shown in Fig. 6 (c) for different values of $D_{t}$: in a region close to the wall, $\langle v\rangle_{x}\neq 0$ and decays monotonically to 0. The decay length decreases as $D_{t}$ decreases until it disappears when $D_{t}\ll D_{a}$ and the thermal noise becomes negligible. The role of thermal noise is not trivial producing a non-monotonic behavior: in a thin space region, close to the wall, $\langle v^{2}\rangle_{x}$ decreases reaching a minimum and then increases until it reaches the limit value $D_{a}/\tau$, as shown in the Fig.6(d). Such an effect can be accounted for theoretically by going beyond the Gaussian closure $n_{2}=0$, i.e. by truncating the coupled system of equations (11) at a higher level, but for space reasons we do not include this possibility in the present analysis.

IV.2.1 Beyond the Gaussian approximation for systems without thermal noise.

Indeed, in the absence of translational noise the screening approximation Eq. (16) cannot be used because the limit $D_{t}\to 0$ is singular. We also find numerically that in this limit $\langle v\rangle_{x}=0$ as expected. In Fig. 7(a) we observe that for $\tau=1$ and $D_{a}=0.1$ and $1$ the second velocity moment $\langle v^{2}\rangle_{x}$ grows until reaches the constant value $D_{a}/\tau$, the quadratic velocity moment of a uniform system. We can distinguish two different behaviors in the interval $[0,2L]$: a persistent region where $\langle v^{2}\rangle_{x}\neq D_{a}/\tau$ where the influence of the wall remains important, and a far region characterized by $\langle v^{2}\rangle_{x}\approx D_{a}/\tau$ which basically is bulk-like. Only in the case $D_{a}=10$ the curve does not saturate and there is no separation between the two regions, since the persistence length $\sqrt{D_{a}\tau}$ is comparable with $L$.

Let us go back to the theory and see that, when $D_{t}=0$ and $\nu=2$, Eq. (11) predicts that, being $n_{1}(x)=n_{0}(x)\langle v\rangle_{x}/(D_{a}/\tau)=0$, the profile $n_{0}(x)$ is simply related to $n_{2}(x)$ by

[TABLE]

Hence, the Gaussian approximation $n_{2}(x)=0$, which was employed to derive to formula (16), fails when $D_{t}=0$ because it would predict a constant profile $n_{0}$ in the force-free region: this is consistent with the UCNA approximation, but not with the numerical result. The breakdown of the Gaussian approximation can be further ascertained by checking the numerical data against the following relations which hold for a Gaussian distribution of velocities: $\langle v^{2n+1}\rangle_{x}=0$ and $\langle v^{2n}\rangle_{x}=C_{2n,2}\langle v^{2}\rangle_{x}$, being $C_{2n,2}$ the binomial coefficient. The numerical study of the third and fourth moments of the velocity distribution, for a system with $D_{t}=0$, displays evident deviations from the Gaussian predictions as illustrated in Fig. 7(b) where a non zero third moment of the velocity is reported and Fig. 7(c) where the distribution has a non-vanishing kurtosis. Having established that the Gaussian closure, $n_{2}(x)=0$, is unfit to capture the observed behavior, a possible remedy to such a deficiency could be a higher order truncation of the series (9), taking into account terms $n_{\nu}(x)$ with $\nu=4$ in the Hermite expansion. Such a procedure leads to a solution for the space-density resembling Eq. (16) with a screening length, $\lambda=\sqrt{5D_{a}\tau/6}$, proportional to the persistence lenght and capturing the phenomenology of the case $D_{t}=0$. This possibility is briefly discussed in Appendix B, for space reasons, whereas here we report the following approximate factorization of the fifth moment of the velocity in terms of lower moments that we found empirically from our numerica data:

[TABLE]

where the factor $10$ is a combinatorial factor which takes into account the number of ways of factorizing the average, is correct. As shown in Fig. 7(d) the comparison between the numerical estimates of $\langle v^{5}\rangle_{x}$ and $10\langle v^{2}\rangle_{x}\langle v^{3}\rangle_{x}$. corroborates the validity of the hypothesis expressed by Eq. (19). In appendix B, we present an argument supporting this factorization of the average.

V Forces on the confining walls

We turn, now, to consider the mechanical properties of the confined active system and derive a formula for the mechanical pressure, $\mathcal{P}_{wall}$, exerted on a harmonic wall by the active gas. To achieve that, we use the following equation expressing the mechanical balance condition between the pressure exerted by the particles on the wall and the one exerted by the wall on the $N$ particles:

[TABLE]

where the upper limit of the integral takes into account the fact that the left wall potential vanishes for any $x>0$. Now, we compute $\mathcal{P}_{wall}$ at the left wall by substituting $\phi^{\prime}(x)=kx$ and obtain:

[TABLE]

where $Nn_{w}$ represents the numerical density at the wall $x=0,2L$, the factor $D_{a}\gamma(\frac{1}{\Gamma}+\Delta)$ has the dimensions of a temperature and is the effective temperature of an active particle confined in a harmonic trap Szamel (2014); Das et al. (2018). The last factor contained in the pressure formula (21)

[TABLE]

is the result of two effects which can be observed when $\tau$ decreases: a) the shift of the peak of the density distribution towards more negative values of $x$ and b) its broadening. For $\tau=0$, being $\Gamma=s(\Gamma,\Delta)=1$, Eq. (21) reduces to $Nn_{w}(D_{a}\gamma+D_{t}\gamma)$ which is the pressure of a suspension of Brownian particles against a wall, so that (21) can be seen as the generalization to the active case of the ideal gas formula.

V.1 The wall boundary conditions

In the pressure formula (21), the constant $n_{w}$ is yet undetermined, and as we shall see below it can be fixed by specifying the system set-up, i.e. its geometrical and physical properties. By multiplying the FPE (5) by $v$ and integrating with respect to $v$ and using the no particle flux condition

[TABLE]

we obtain the following expression relating a total derivative to the wall force:

[TABLE]

We now integrate with respect to $x$ Eq. (23) between $-\infty$ and $\bar{x}$ where $\bar{x}>0$:

[TABLE]

where we have taken into account Eq. (20).

V.1.1 Semi-infinite system.

Now, we restrict our analysis to the case $D_{t}$ not too small with respect to $D_{a}$, thus excluding the singular limit $D_{t}=0$. In order to simplify the analysis, we consider the limit $\bar{x}\gg\lambda$, which allows us to assume that the term $n_{2}(\bar{x})\approx 0$ and we evaluate the left-hand side of Eq. (24) using Eq. (16) with the result:

[TABLE]

Using the explicit representation of $\mathcal{P}_{wall}$, Eq.(21), we express the probability density at the wall, $n_{w}$, in terms of the probability density $n_{M}$ and of the parameters of the model:

[TABLE]

Let us remark that we always have $n_{w}>n_{M}$ because both factors $\frac{1+\Delta}{\frac{1}{\Gamma}+\Delta}$ and $\frac{1}{s(\Gamma,\Delta)}$ are larger than $1$ if $\tau>0$, so that the wall density is higher than the density at midpoint and we may argue that there is a positive surface excess. Only in the Brownian limit $\tau=0$ we obtain $n_{w}=n_{M}$ for all values of $k$. If now we take the limit of a semi-infinite system $L\to\infty$, we have $n_{0}(x)=(n_{w}-n_{M})e^{-x/\lambda}+n_{M}$ and we can make the identifications $\rho=Nn_{M}$ and $\rho_{w}=Nn_{w}$, where $\rho$ and $\rho_{w}$ are the bulk and wall numerical densities, respectively. Finally, using the condition that $\mathcal{P}_{bulk}=\mathcal{P}_{wall}$, necessary in order to have mechanical equilibrium, we identify the r.h.s. (25) with the bulk pressure, $\mathcal{P}_{bulk}$, of a uniform system at density $\rho$, i.e. \mathcal{P}_{bulk}=D_{a}\gamma\,\Bigl{(}1+\Delta\Bigr{)}\rho.

V.1.2 Wall Pressure in a slit system.

We turn, now, to study the pressure in a slit-like geometry with the purpose of understanding how the force acting between two parallel plates immersed in a solution of active particles, depends on the wall separation. To find the density at the wall, $Nn_{w}$ we must compute $n_{0}(x)$ in the whole space, match the expression (8) with (16) at each wall and finally normalize the profile. The probability density profile $n_{0}(x)$ can be written as:

[TABLE]

where $n^{left(right)}_{0}(x)$, the probability density distribution with walls located at $0\,(2L)$, is given by Eq. (8) and $n_{0}^{f}(x)$ the density in the free region $[0,2L]$ is given by Eq. (16). From the normalization of the probability distribution we have the condition:

[TABLE]

where the factor $2$ in the last term takes into account the symmetry of the two walls. After performing the integrals and eliminating $n_{M}$ we obtain the relation:

[TABLE]

We point out that the first term takes into account the finite width of the peak of $n_{0}(x)$ in the regions $x\leq 0$ and $x\geq 2L$ due to the softness of the walls: it vanishes for all values of the remaining parameters when $k\rightarrow\infty$. As a check we consider the equilibrium limit $\tau\to 0$

[TABLE]

as expected, since it corresponds to the situation of an overdamped passive system in contact with two independent white noise sources. In the hard wall limit $k\to\infty$ and $\tau=0$ we find $n_{w}=1/2L$ and $\mathcal{P}_{wall}=\frac{N}{2L}D_{a}\gamma\left(1+\Delta\right)$, the ideal gas equation of state corresponding to a system subjected to two white noise baths.

V.2 Effective force between plates and Casimir effect

Under the same hard-wall limit $k\to\infty$, but with $\tau>0$ one can see that there is an accumulation effect at the wall: in fact, according to Eq.(29) we have $n_{w}\geq 1/2L$. Finally, the wall pressure can be computed inserting Eq. (29) in Eq. (21). For the sake of simplicity, we write its expression in the limit $k\to\infty$ but $\tau$ finite:

[TABLE]

In order to ascertain the activity induced force acting on parallel plates, we must compare the pressures of two systems having different sizes, $L_{1}$ and $L_{2}$, with $L_{1}<L_{2}$ but the same average density of each system, $\bar{\rho}=\frac{N_{i}}{2L_{i}}$, is identical. Clearly, the system with the smaller size according to Eq. (31) will exert the smaller pressure on the walls. We may conclude that two parallel plates surrounded by a sea of active particles and separated by a distance $2L$ will experience an effective attraction according to Eq. (31). Precisely, if $L_{2}\to\infty$, but $\bar{\rho}$ is fixed and $L_{1}/\lambda\gg 1$ we have

[TABLE]

Thus, in the low-density regime we consider we find that the force increases linearly the active power $D_{a}$ and depends monotonically on plate separation in agreement with the simulation results by Ni et al. Ni et al. (2015) and the theoretical prediction of Vella et al. Vella and et al (2017) of a decay $\propto 1/L$.

Such an effect can also be illustrated by the following gedanken-experiment: let us consider a finite system and insert a third hard wall C, identical to the first two, at an arbitrary position $y\in(0,2L)$, in such a way that the average numerical densities in the two resulting compartments are equal: $N_{l}/y=N_{r}/(2L-y)$ with $N_{l}+N_{r}=N$. According to the present theory the pressure difference between the left and right compartment is given by:

[TABLE]

Of course, for a Brownian system, if the average densities in the left and in the right compartments are set to be equal for any choice of the wall position $y$, i.e. if the condition $N_{l}=Ny/2L=N_{r}=N(2L-y)/2L$ is satisfied, the pressure difference, $\Delta\mathcal{P}$ vanishes.

The physical reason of such a phenomenon is strictly related to the accumulation of active particles in front of a wall and clearly emerges in the Eq. (29): increasing $L$ the constant $n_{w}$ grows, meaning that more particles push on the wall and exert a larger pressure.

On the other hand, our prediction suggests a completely different situation in the active case which we show in the Fig. 8. In particular $\Delta\mathcal{P}>0$ if $y>L$ and $\Delta\mathcal{P}<0$ if $y<L$: in fact, the particles in the small compartment exert a smaller pressure on the wall C than the one exerted by those in the larger compartment, in spite of the fact that the numerical densities are equal. For small separations $2L$ the approximation $n_{2}(x)\approx 0$ ceases to be correct and we do not expect that the force obeys anymore the scaling $L^{-1}$, however there is some room for improvement, for instance, by employing higher order closure approximations, such as including terms $n_{3},n_{4}$ etc, but in the present study we do not pursue such a possibility.

VI Conclusions

To conclude, let us remark that the physical effect of activity is twofold: i) it determines a non-uniform density profile because active particles accumulate near the wall and, in the case of deformable boundaries, penetrate inside them; on the contrary, a system of non-interacting Brownian particles in the hard wall limit would not develop any density gradient; ii) the pressure exerted on the walls of a slit of width $2L$ by a system of average density $\bar{\rho}$ depends on the wall separation. The second phenomenon is relevant when the confinement length becomes comparable with the persistence length. Hereafter, we summarize the main achievements of the present work.

1. We have extended the study of the AOUP to the case of one one-dimensional non-interacting particles under confinement, and, going beyond the UCNA approximation, we do not integrate out the active noise but retain it as a "velocity" variable. By considering a simple parametrization of the bounding potential, we have been able under reasonable approximations to derive simple expressions for the density profile and polar field and eliminate some negative features of the UCNA solution, such as its jumps in correspondence of discontinuities of the potential, and its failure to account for polar order near a boundary. Our treatment introduces a healing length which produces smoother density profiles which have been found in good agreement with the results of numerical simulations.

2. The moment method employed to approximate the non equilibrium distribution function $f(x,v)$ predicts an effective force between two plates immersed in an active suspension, similar to the force in the classical Casimir effect: indeed, the attraction between the plates is due to the active particles, which could represent active bacteria, while low-density Brownian particles, i.e. colloidal particles, do not exert any appreciable force on the plates. The possibility of generating and controlling the force between immersed objects, for instance by tuning the illumination of an active suspension or modifying its concentration and/or temperature, is quite interesting and offers an alternative to other techniques which instead require the chemical modification of the surfaces. Finally, the effective temperature of such suspensions, which determines the intensity of such a force, can be higher than the solvent temperature Rohwer et al. (2017).

Future work will concern the extension of the theory to higher dimensions in order to treat active solution-mediated interactions between inclusions of more general shape Vishen et al. (2018) or moving pistons Caprini et al. (2017). Including the interaction among the particles is also a challenge and we may expect that with increasing density the excluded volume effects could lead to an effective repulsion for some values of the plate separation.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

We thank Andrea Puglisi and Angelo Vulpiani for illuminating discussions.

Appendix A Stationary distribution for the harmonic oscillator in the $(x,v)$ variables

We know the full stationary solution of (5) for the harmonic oscillator. It can be written as a double Gaussian with a velocity distribution whose peak changes with $x$. This peak corresponds to an x-dependent velocity, $\langle v\rangle_{x},$ which is also the mean velocity at fixed $x$.

[TABLE]

with

[TABLE]

Let us remark that $f_{h}(x,v)$ due to the presence of a non-vanishing average velocity $\langle v\rangle_{x}=$ has the form of the distribution function of a system in local but not global equilibrium, in contrast with the case $\Delta=0$.

Appendix B How to rationalize the non-Gaussian closure

As we can see from the structure of the solution in the potential free-region the screening length vanishes when $D_{t}\tau\to 0$ even though $D_{a}$ remains fixed. In order to remove such a nonphysical feature, we must consider carefully the limit $D_{t}=0$. In this case, the study hierarchy (11) becomes relatively simple and allows to predict a non-vanishing decay length and sheds some light on the form of the closure (19). We begin by writing explicitly the hierarchy assuming $n_{1}(x)=0$ in this limit:

[TABLE]

An option is to break the hierarchy by setting $n_{5}(x)=0$ and after eliminating $n_{4}(x)$ we write

[TABLE]

We obtain a closed set of linear differential equations that can be solved by combinations of exponentials of the form $e^{\pm\mu x}$, with $\mu$ determined by a simple algebraic equation. We find $\mu^{2}=\frac{6}{5}\,\frac{1}{D_{a}\tau}$ and the profile is given by $n_{0}(x)=A\cosh(\mu x)+C$. We now try to verify the working hypothesis (19). Using the Hermite expansion (9) when $n_{1}=0$ and $v_{a}^{2}=\frac{D_{a}}{\tau}$ we obtain the relations:

[TABLE]

Since we have assumed $n_{5}=$ we obtain the equality:

[TABLE]

Such a relation is compatible with Eq.(19) only in the regime when $\langle v^{2}\rangle_{x}$ can be replaced the constant factor $v_{a}^{2}$ in Eq (46). This is possible if the space dependent average $\langle v^{2}\rangle_{x}-v_{a}^{2}=2v_{a}^{4}\frac{n_{2}(x)}{n_{0}(x)}\approx 0$, i.e. in a regime of small $\tau$.

The empirical relation (19) instead is consistent with the choice $n_{5}(x)=\frac{n_{3}(x)n_{2}(x)}{n_{0}(x)}$. However, the substitution of such a relation into Eqs. (45) leads to a closed set of non-linear equations which cannot be solved by simple analytic methods.