Confined active Brownian particles: theoretical description of propulsion-induced accumulation

The stationary-state distribution function of confined active Brownian particles (ABPs) is analyzed by computer simulations and analytical calculations. We consider a radial harmonic as well as an anharmonic confinement potential. In the simulations, the ABP is propelled with a prescribed velocity along a body-fixed direction, which is changing in a diffusive manner. For the analytical approach, the Cartesian components of the propulsion velocity are assumed to change independently; active Ornstein–Uhlenbeck particle (AOUP). This results in very different velocity distribution functions. The analytical solution of the Fokker–Planck equation for an AOUP in a harmonic potential is presented and a conditional distribution function is provided for the radial particle distribution at a given magnitude of the propulsion velocity. This conditional probability distribution facilitates the description of the coupling of the spatial coordinate and propulsion, which yields activity-induced accumulation of particles. For the anharmonic potential, a probability distribution function is derived within the unified colored noise approximation. The comparison of the simulation results with theoretical predictions yields good agreement for large rotational diffusion coefficients, e.g. due to tumbling, even for large propulsion velocities (Péclet numbers). However, we find significant deviations already for moderate Péclet number, when the rotational diffusion coefficient is on the order of the thermal one.


Introduction
The distinct feature of active systems is their autonomous motion powered by an internal energy source or by utilizing energy from their environment [1][2][3][4][5][6][7][8][9]. Prototypes of such self-propelled objects are omnipresent in biology and range from the macroscopic scale of flocks of birds and mammalian herds [3] down to moving bacteria [2,6,10] on the micrometer scale. Various synthetic microswimmers have been designed [7,[11][12][13][14] in order to mimic the behavior of biological microorganisms or to serve as micro-and nano-machines with promising applications in technology and healthcare [7,11]. The rational design of synthetic microswimmers for such applications greatly benefits from a thorough understanding of the relevant physical mechanisms. This comprises the propulsion mechanism of individual microswimmers, their collective behavior, and their properties in external fields and confinement [6,7]. In order to tackle, in particular, the last two aspects, simplified and generic models of microswimmers are applied such as the active Brownian particle (ABP) [4,6,7,[15][16][17][18][19][20][21][22][23]. This model has been exploited to characterize the phase behavior of active systems [17][18][19][20][21]. More importantly, ABPs serve as model systems to characterize the nonequilibrium properties of active systems. Here, the existence of an equation of state has been discussed intensively [8,[24][25][26][27][28][29][30]]. An even more interesting aspect is the existence of a stationary-state distribution function of such systems. Due to their nonequilibrium character, active systems violate detailed balance [2,8,21,[31][32][33], which complicates the calculation of a stationary-state distribution function or might even render it impossible.
Various theoretical approaches have been adopted to describe the nonequilibrium properties of ABPs, specifically their phase behavior [8,21,[34][35][36][37][38], but also effective potentials have been determined [39][40][41][42]. In order to arrive at analytical results, approximations are typically necessary to overcome the complexity of the underlying equations. Thereby, somewhat simpler theoretical models are considered, which focus on particular aspects of the original ABP typically treated in computer simulations.
In this article, we perform Brownian dynamics simulations and analytical calculations of active particles in three dimension, which are confined spatially by a radially symmetric potential, and determine their spatial distribution function. In simulations, we apply the well-known ABP description of an active particle, namely a spherical colloid with given propulsion of fixed magnitude, which changes direction in a diffusive manner [6-8, 17, 18, 20]. A somewhat simpler theoretical model is studied analytically, with independent stochastic processes along the Cartesian coordinates of the active velocity, which corresponds to a Gaussian, but non-Markovian, colored-noise process for the activity [6,8,17,22,33,[43][44][45]. This model is denoted as active Ornstein-Uhlenbeck particle (AOUP) [33], a notation we also adopt here. The AOUP description has been applied in wide range of simulation and theoretical studies of confined and interacting active particles, such as motility-induced phase transitions [33,40] and alterations of the glass behavior [46], or the accumulation phenomena at walls [39], to name just a few.
Here, we particularly consider ABPs and AOUPs confined in a harmonic potential [4,[47][48][49][50][51][52][53]. As for passive systems, harmonically bound active particles are an excellent model system to extract the generic features of confined nonequilibrium systems. Specifically, the equations of motion of AOUPs can be solved analytically, and thus provide valuable insight into the out-of-equilibrium dynamics of active systems. Correspondingly, several theoretical studies of AOUPs have been performed, predominantly in one dimension [49,51,53,54]. As an extension, we provide the full solution of the dynamical equations in terms of position and velocity for arbitrary dimensions, taking into account thermal as well as self-propulsion contributions. The spatial distribution, obtained by integration of the general distribution function over velocity, of an AOUP is a Gaussian centered at the origin of the potential. However, the off-center accumulation observed for ABPs can be captured analytically by considering the spatial conditional probability distribution function of AOUPs at a given propulsion velocity. We like to emphasize that only the conditional probability distribution function for AOUPs is suitable for a comparison with the spatial distribution for ABPs. As a result, the analytical approach qualitatively reproduces the features obtained in ABP simulations.
For anharmonic potentials, we apply the unified colored noise approximation (UCNA) to find an approximate stationary-state solution [55,56]. The comparison between the simulation and analytical results reveals the quantitative suitability of the applied approximation schemes. Our studies exhibit good agreement between theory and simulation results for large rotational diffusion coefficients D R of the active process for a wide range of propulsion velocities. In contrast, for small D R , agreement is obtained for small propulsion velocities only. In general, our studies are useful in the attempt to find a suitable analytical description of the nonequilibrium properties of systems of ABPs.

Active Brownian particle
An ABP is a spherical colloidal particle of diameter σ, which is propelled with constant velocity v 0 along its unit orientation vector e in d dimensions [4, 17-20, 22, 27, 44]. Its translational motion is described by the Langevin equation Here, r andṙ are the particle position and velocity, respectively, = - F U is the force by the external potential U, and g T is the translational friction coefficient. Thermal fluctuations are captured by the Gaussian and Markovian random force G( ) t with zero mean and the second moments where k B is the Boltzmann constant and T the temperature. Strictly speaking, the stochastic process G( ) t does not have to be of thermal origin only, and k B T not necessarily the thermal energy. It accounts for white-noise stochastic processes affecting the translational motion of the ABP and may comprise active contributions. Independent of that we will refer to it as thermal motion in the following.
The orientation e changes independently as where ĥ is a Gaussian and Markovian stochastic process with zero mean and the second moment in 3D; σ is the diameter of the colloid 2 . In general, D R can be larger than the thermal value due to activity-induced rotation. As a particular result, the correlation function of the propulsion vector is given by [27,57] We solve the translational equations of motion (1) via the Ermak-McCammon algorithm [58]. The equation of motion (3) for the orientation vector is solved as follows. By a unnormalized estimation of the vector e for the time interval Dt is obtained, where cos sin , sin sin , cos 7 T in spherical coordinates [27]. The x e (x f q Î { } , ) are unit vectors, which follow from equation (7) by differentiation and normalization, and the h D x are again Gaussian and Markovian stochastic processes with zero mean and the second moments i.e.,

Active Ornstein-Uhlenbeck particle
The propulsion velocity e v 0 , with the equation of motion (3), obeys the strict condition = | ( )| e t 1. For our analytical considerations, we lift this restriction and consider the frequently used equation for the propulsion velocity [4,6,22,33,44,51]; here, the Cartesian velocity components are independent. The damping factor g R is related to the rotational diffusion coefficient according to g = -( ) d D 1 R R , and h is a Gaussian and Markovian stochastic process with zero mean and the second moments , , . The Langevin equations (9) with noise amplitudes (v D 0 2 R ), which are independent of the stochastic variables (v), describe a process denoted as Ornstein-Uhlenbeck process [59]. Hence, we will refer to the active particle obeying equations (1) and (9) as AOUP in the following, adopting the notation of [33].
With the Gaussian character of h and the correlation function (11), equation (1) is the Langevin equation of a particle in the presence of white noise (G) and, in addition, non-Markovian colored noise. The non-Markovian character prohibits a general solution of the Langevin equations (1), and no closed Fokker-Planck equation can be derived in general [55]. Analytical solutions exist for a constant force and a harmonically bound particle (linear force), typically for one dimension [49,51,[53][54][55]. For nonlinear forces approximations have to be applied. 2 For strictly 2D systems, the hydrodynamic diffusion coefficient D T diverges in the limit of system size  ¥ L ( s ( ) D L ln T ). However, in simulations and experiments-for the latter, there is no strictly 2D system-typically slit geometries, i.e. quasi 2D, or 'semi-infinite' systems with self-propelled particles absorbed at a wall are considered. Our simulations of a colloid show that the rotational diffusion coefficient in a quasi-2D system is approximately equal to the bulk value. However, the translational diffusion coefficient is reduced [66,67]. Hence, a ratio s can be defined also for such systems, however, Δ is smaller than the value D = 1 3 of a 3D system.

Equations of motion and Fokker-Planck equation
The radially symmetric harmonic potential leads to independent equations of motion for the Cartesian components of the position vector, i.e.
Complementary, the AOUP dynamics is described by the distribution function , , refers to the Cartesian coordinate directions, which obeys the Fokker-Planck equation [59] y k y y g y g y g y

Solution of the equations of motion
Due to the Gaussian nature of the stochastic processes and the fact that a Gaussian is determined by its first and second moment, the full time-dependent solution of the Fokker-Planck equation can be obtained [59].
Explicitly, the conditional distribution function is given by T and the initial values r 0 and v 0 at the time t=0. Note that we will use = a v v d 0 0 in the following. The symmetric matrix L is the covariance matrix rr rv rv vv 2 2 L -1 its inverse, and L | |the determinant. The appearing moments conveniently follow from the linear Langevin equations (9) and (13), which yields  Since the dynamics of the velocity v is independent of the position of the particle, the relaxation behavior of equations (18) and (21) is naturally determined by g R only. Two independent relaxation mechanisms are present for the moments involving positions. For k g  1 R , the long-time behavior is dictated by the relaxation time g 1 R of the propulsion velocity. In the opposite limit k g  1 R , the relaxation behavior of L rv is still governed by g R . Only for á ñ a r and L rr relaxation due to confinement is important, whereby L rr changes with time as similar to a passive Brownian particle [59], but with an activity-dependent amplitude. Two relaxation processes have already been found in [49] for the particle's positional autocorrelation function.
In the stationary state, the distribution function turns into The distribution function (23) reduces to the expression derived in [49,53] for d=1. It is important to note that equation (23) is different from the stationary-state distribution function of a passive Brownian particle [4,59] even in the limit  v 0 0 , where the latter yields independent distributions for the potential and the kinetic energy, since the term does not vanish. The result as such is not surprising, because of the differences in the underlying models. For AOUPs and ABPs, the overdamped dynamics of the particle position is considered and the velocity (orientation) is introduced as an independent degree of freedom. As a consequence, the forces on the AOUP/ABP determine its velocity. In contrast, in case of passive Brownian particles, the forces rule the acceleration.
In general, the coupling of position and velocity prevents the definition of a suitable effective temperature based on the stationary-state distribution function, aside from inconsistencies with dynamical quantities [49]. However, an effective temperature T e can be introduced in the limit g k  1 T 0 2 R as in [49]. Aside from k B T and the dimension d, T e agrees with the definition in [49]. The other Λs in equation ( R B e . In the limit g  ¥ R , Lrv 1 vanishes and independent distribution functions are obtained for r and v, identical with those of a passive Brownian particle. This is not surprising, since in this limit the stochastic process (11) turns into a Markovian process with a finite amplitude. Hence, v can be neglected in equations (1) and (13), which leads to a decoupling of the equations of motion for the position and velocity. However, as already pointed out in [49], the interesting case for active particles is the limit of strong propulsion g k  1 R .
We would like to mention that the stationary-state second moments following from equations (20) and (21) for AOUPs are identical with those of ABPs, since they involve only up to second moments of the velocity v (respectively e), i.e. the correlation functions (5) and (11). However, it has to be kept in mind that the distribution function for an ABP is non-Gaussian in general. An alternative derivation of the second moment for the particle position has been presented in [54].
Integration of equation (23) over the positional coordinates yields the velocity distribution function Equations (29) and (30) also simply follow from the Fokker-Planck equation of the Langevin equation (9) for the velocity.
Integration of equation (23) over the velocity coordinates yields the positional distribution function , 0 with (see equation (20)) Naturally, we obtain a simple Gaussian distribution function for y ( ) r of the AOUP, with a maximum at the minimum r=0 of the potential. Computer simulations of the ABP model obeying equations (1) and (3) yield a different distribution with, depending on activity, an off-center maximum [47,51].
Strictly speaking, the distribution function (32) should not be compared with that of an ABP. In equation (32), all velocity contributions have been integrated out, whereas for an ABP, the magnitude of the velocity is fixed. The latter corresponds to a distribution function y ( ) r at fixed v 0 . The corresponding distribution function of an AOUP is the conditional probability distribution function y ( | ) r v 0 . Hence, this is the suitable distribution function of an AOUP to be compared with the spatial distribution function of an ABP.

Conditional probability distribution function
The conditional probability density y ( | ) r v of the position r for a fixed velocity v is defined by the relation [59] Using equations (23) and (30), we find (35) over the orientation of the velocity v. Evaluation of the integrals yields: is the modified Bessel function of the first kind.
• Three dimensional (3D) system We introduce the Péclet number Pe, the ratio Δ between the translational and rotational diffusion coefficient, and the strengthk of the harmonic potential as dimensionless numbers to characterize the system, where s s Alternatively, the Péclet number can be introduced. However, we prefer equation (38), since we want to resolve the effect of the two independent quantities v 0 (Pe) and D R (Δ) separately.
In both expressions, the termk 1 corresponds to the thermal (white-noise) contribution. Activity dominates when the Péclet number exceeds the critical value The following asymptotic dependencies are obtained for equations (36) and (37): , and the Gaussian function in equations (36) and (37) turns into the Boltzmann distribution - of a passive particle. The relevance of the Bessel or the hyperbolic sine function depends on the values ofk and Δ, respectively. However, as long as Pe 1 c , the positional distribution is solely determined by the Boltzmann distribution of a thermal passive particle.
Interestingly, Pe c increases then linearly with k, and can assume large values. Hence, for strongly confined AOUPs and moderate Δ,  Pe Pe 1 c and the termk 1 in equations (40) and (41) dominates even for Péclet numbers significantly larger than unity. . Thus, the argument of the Bessel and hyperbolic sine function is proportional to Pe 1 , i.e. is small for large Pe, and these functions will only weakly contribute to y ( ⃒ ) r v 0 . Consequently, we obtain a Boltzmann distribution with the effective temperature of equation (28) , because  Pe 1. We like to emphasize that Δ is typically smaller than unity. In case of thermal motion, D = 1 3 for a colloid in a fluid. Activity can cause larger rotational diffusion coefficients and, hence, smaller Δ, which has been confirmed experimentally for E. coli bacteria [6,[60][61][62] and Chlamydomonas reinhardtii cells [60].
ABPs confined in a 2D harmonic potential have been considered in [48]. The analytical approximations provided in that article also follow from the distribution function (36). Inserting equations (38) in equation (36), we find in the limit This expression coincides with that derived in [48]. In the opposite limit g  ¥ , Taylor expansion of the Bessel function to lowest order yields N c is the normalization coefficient. In terms of the dependence on the potential U h (r) and the Péclet number, this result agrees with the expression derived in [48]. However, we find a different dependence on g = D R R and the strength of the potential. In [48], the term s ( ) U h is missing and Δ appears only linearly. This might be related to the alternative approach of [48]. The advantage of our result is that we not only obtain limiting expressions, but can also describe the crossover from small to large g R .
Finally, the conditional probability density y ( | ) v r (y y y ) is given by

Comparison of ABP and AOUP in harmonic potential
Here, we will compare the theoretical results derived in sections 4.2 and 4.3 for an AOUP in a harmonic potential with respective simulation results for an ABP.

One-dimensional system
Typically one dimensional (1D) harmonically confined AOUPs are considered [49,51,53,54]. In contrast, our approach presented so far, with the definitions (4) and (10), is only meaningful for > d 1. However, equations (9) and (10) can simply be adapted to a pure 1D motion by setting The general solution presented in section 4 still applies, we just have to set d=1 and replace - Integration of the (1D) conditional probability (35) over the two possible (plus/minus) directions of motion yields the conditional probability distribution function R . An estimation of x m , respectively r m of higher dimensional systems, follows from the (general) condition of vanishing (average) velocity = r 0 (equation (1)), i.e. the balance of the active and confining force, which yields T and, hence, This is identical with the expression following from equation (51) for strong confinement D  k 1 1 . An activity-induced off-center-peaked distribution function for harmonically confined run-and-tumble particles (RTPs) [63] has been obtained in [47], and a power-law distribution function is derived, with an exponent depending on the tumbling rate (see also [51]). This power-law distribution exhibits a singularity at the position x m for strong confinement. In contrast, the Gaussian function is typically finite at its maximum, and | | x values larger than | | x m are accessible. Thus, the distribution functions for ABPs and RTPs are qualitatively different. However, common to both results is the presence of off-center peaks of the probability distribution function at the same location for strong confinement. Figure 1 displays the radial conditional probability distribution function, y ( | )

Three-dimensional system
r v 0 , of AOUPs and the radial distribution function, y ( ) r , of ABPs obtained from simulations for various force constants and Péclet numbers. The distribution functions of the ABPs exhibit the expected activity-induced off-center accumulation with a maximum shifting to larger radial distances with increasing Péclet number as already obtained in [51] for a 2D system. A qualitatively similar behavior is obtained for AOUPs, most pronounced for k g D »  k 1 R ( figure 1(b)). The critical Péclet numbers (42) for = k 1, 10 1 , and 10 2 are Pe c = 5.1, 11.2, and 71.6. As discussed in section 4.  (37) is still close to a Boltzmann distribution. Yet, the argument of the hyperbolic sine function is not necessarily small anymore ( D  k 1), which leads to a pronounced off-center maximum.
In the opposite limit,  Pe Pe c , the distribution function y ( | ) r v 0 qualitatively describes the behavior of ABPs. The quantitative agreement depends on the stiffness of the harmonic potential. For = k 1, reasonable agreement is obtained for all considered Péclet numbers < Pe 10 2 (not shown). This is not surprising, since this corresponds to the limit g k > 1 R , and a Boltzmann distribution is obtained for y ( ⃒ ) r v 0 , with the effective temperature of equation (28), as discussed in section 4.3 (ii). For stronger confinement, significant deviations appear between y ( | ) r v 0 of AOUPs and the distribution function of ABPs, specifically for = k 10, the more the larger Pe. In contrast, the differences are smaller again for = k 100, where Pe c is larger. Estimations of the radial positions with vanishing velocity according to equation (53)-we obtain the same (radial) value in any dimension-for the various Péclet numbers are indicated by vertical short, black bars in figure 1, and the respective maxima of the distribution functions from simulations by colored bars. For small Péclet numbers, the exact position of the maximum of the simulation data is difficult to identify due to statistical inaccuracies. As soon as a clear maximum appears, its position agrees well with the theoretical prediction (53) at large potential stiffness ( = k 100). Here, the difference between the position of the maximum and r m decreases from approximately 10% to 5% with increasing  Pe 40. The deviations for = k 10 are in the range of -35% 20% for  Pe 20. Hence, equation (53) provides an upper estimate of the maximum of the distribution function, whose accuracy improves with increasing potential stiffness and Péclet number. The maxima of the analytical approach are shifted toward smaller radial distances, as is evident from figure 1, but the linear dependence on Pe is maintained. Figure 2(a) illustrates the dependence of the radial conditional distribution function, y ( | ) r v 0 , of AOUPs on Δ at a fixed Pe. Thereby, the inset of figure 2(a) shows that the maximum of y ( | ) r v 0 shifts to large radial distances with increasing Δ, or decreasing g = D 2 R R . The hyperbolic sine term in equation (37) suggests that y ( | ) r v 0 increases exponentially for appropriate r as long as the Gaussian function is changing only slowly. Such an exponential dependence has been presented in [51]. To highlight this asymptotic behavior, we present the distribution function where x 0 is given by The latter reduces to x = D - Interestingly, y x ( | ) v 0 is independent of the Péclet number for  Pe Pe c . Note that this is always the case when the thermal contribution by G is neglected. Thus, we obtain a unique distribution, y x ( | ) v 0 , independent of Pe for a given Δ. The dependence on the Péclet number is adsorbed in the variable ξ.
Similarly, the radial probability distribution function of an ABP, with y ( ) r the probability distribution function from figure 1(a), is displayed in figure 2(b). Again, we obtain close agreement of the various curves for large Péclet numbers, with the same scaling variable ξ as for AOUPs. The range of exponential growth is rather short, but would increase for D  1.
However, D < 1 is more relevant for self-propelled particles, as pointed out before. Here, the range of radial distances over which the distribution functions increase exponentially is rather shorter than longer. The presence of an exponential regime depends on D k . For D  0, x  0 0 2 and the Gaussian function in equation (54) decays fast with increasing ξ, suppressing an exponential regime. In contrast, for very strong potentials with D  k 1, we observe a pronounced, Pe-independent exponential regime. The differences between the ABP and AOUP results of figure 1 reflect the variance in the equations of motion (3) and (9). Since the second moment of both models is identical, we characterize their differences by the kurtosis K. Note that both, y ( | ) r v 0 of an AOUP and the distribution function of an ABP are non-Gaussian in general. With the conditional distribution function (37) for a AOUP, we obtain in 3D  Figure 3 shows kurtoses for the parameters of figure 1. K decreases with increasing Péclet number, corresponding to a platykurtic kurtosis; this means that the distribution functions are flatter than a Gaussian and decay faster at the tails. Thereby, the kurtoses of the ABPs are smaller than those of the equivalent AOUPs. As discussed before, at smallk, the distribution functions are essentially Gaussian with K close to 5/3.

Unified colored noise approximation
An analytical solution of equations (1) and (9) is the exception rather than the rule. Naturally, a Fokker-Planck equation can be derived for the coupled Langevin equations due to their Gaussian and Markovian character [59]. The stochastic process for r only, however, is no longer Markovian, which is reflected by the colored noise correlation function equation (11). This implies the lack of a Fokker-Planck-type equation for the (reduced) probability density function y ( ) r t , [55]. At the moment, at best an approximate Fokker-Planck-type equation can be derived [55,64,65]. One of the approximation schemes is denoted as UCNA [55]. This approximation has been applied to ABPs in [30,39,41,52,56]. Here, our goal is to elucidate the applicability and validity of the UCNA for the ABP model of section 2 confined in an anharmonic potential.

UCNA distribution function of radial potential
A derivation of the multidimensional stationary-state probability distribution function within the UCNA is presented in the appendix (equation (A.8); see also [55,56]). In case of a radially symmetric potential This distribution function agrees with that presented in [56] except of the factor in front of the exponential function-equation (6) of [56] contains an additional term with a first derivative of the potential. Note that the normalization factor Z is In the limit g  ¥ R , which corresponds to the limit of an overdamped dynamics (see equation (A.2)), the potential term in the exponent of equations (59) and (A.8) dominates and the distribution function reduces to with the respective normalization factor Z 0 . Hence, a Boltzmann distribution with the potential U and the effective temperature (28) is obtained. Note that the white-noise process G has been neglected in UCNA. This result is consistent with the general expectation for the Langevin equation (1) in the limit g  ¥ R discussed in section 4.2 following equation (28).
The conditional distribution function for the velocity y ( | ) u r (A.12) is y g g p g g g g g g  In addition, the mean square velocity (A.13) is given by

Harmonic potential
For the harmonic potential U h of equation (12), the distribution function (59) is equal to y ( ) r of equation (

Comparison of simulation and UCNA results for anharmonic potential
We consider the potential where s = r r , as an example for an anharmonic potential. Here, we scale the potential energy with the thermal energy k B T, which provides a consistent definition of dimensionless parameters with the previous sections.
The radial distribution function (59) now reads in terms of the Péclet number and Δ. Hence, y ( ) r is determined by the two independent parameters, Pe and Δ. In the following, we elucidate the extent to which the UCNA provides a solution of the dynamics of AOUPs and ABPs, respectively. Figure 4 displays probability distribution functions of an AOUP without thermal noise for various Péclet numbers. Obviously, the simulation data agree very well with the UCNA for  Pe 10, and they are still qualitatively reproduced for Pe = 100. For  Pe 1, the distribution function monotonously decreases with increasing radial distance r, whereas a peak appears at larger Pe. The peak for Pe = 10 is caused by the factor with the second derivative in front of the exponential function of equation (65). This term increases strongly whenr approaches the singularity at = r 5. For  Pe 10, the exponential factor decays fast with increasingr , before the term D -

Comparison of AOUP simulations with UCNA
Here, the distribution function can well be described by  (1), without thermal noise, and equation (9). The dashed lines are analytical results within UCNA (equation (65)).
equation (61). With increasing Pe, the exponential factor decays at larger radii and the pre-factor gains importance, resulting in the peak. Deviations between simulations and theory appear for > Pe 10. These differences are even more pronounced for D > 0.1. At D = 1, good agreement is still obtained for  Pe 1, but larger Pe lead to marked discrepancies between theory and simulation curves, with peak heights of the simulation data exceeding those of the theoretical prediction substantially, similar to those displayed in figure 5(b) for ABPs. Here, theory predicts the behavior only rather qualitatively. Considering smaller Δ, we find good agreement for even larger Péclet numbers.
In the limit D  1, equation (65) reduces to Interestingly, this function is independent of Δ for any Pe. A comparison shows that equation (66) indeed describes the simulation data very well for  Pe 1. Good agreement is already obtained for » Pe 1 in the vicinity of the peak of the distribution function. Only for smallr deviations are found due to the fast decay of ¶ ¶Ũ r 2 2 with decreasingr . Figure 5 shows simulation results for the radial positional probability distribution function of ABPs. Qualitatively, the shape of the curves is very similar to the AOUP results. However, the appearing peaks for ABPs are sharper and higher. Quantitatively, we obtain close agreement with the theoretical prediction over a range of Péclet numbers, where the range depends on Δ. For D = 0.01 (not shown), we observe deviations for  Pe 100. At D = 0.1 (Figure 5(a)), the curves for Pe = 10 already digress from each other, and for D = 1 ( figure 5(b)), deviations are present for  Pe 1. The inset of figure 5(a) displays results for a 2D system. The qualitative behavior is the same as in 3D. However, the peaks are somewhat higher and narrower. The agreement with the theoretical approach is equally good.

Comparison of ABP simulations and UCNA
Hence, the simulation and theoretical results are consistent with the assumption of UCNA, and an adequate theoretical representation is obtained for g  ¥ R , or D  0. Alternatively, the UCNA applies for  Pe 1 in the limit D  ¥. However, the latter limit is of little practical importance, since typically D < 1 for biological and synthetic microswimmers. The limitations of the UCNA for microswimmers are evident for D » 0.3, where thermal fluctuations determine the rotational diffusion. Here, we already see deviations between simulation results and the theoretical prediction for  Pe 1 (see figure 5(b)). We neglect thermal fluctuations in the current study. For Péclet numbers on the order of unity, thermal fluctuations may matter and determine the properties of the active system. Hence, the distribution function of the radial position would be given by the Maxwell-Boltzmann distribution y b -(˜) (˜) r U exp rather than the UCNA. Considering the condition (52) for the balance of active and confining forces, we find very good agreement of the peak position from simulations with the theoretical prediction for all Péclet numbers.  (1), without thermal noise, and equation (3). The dashed lines are analytical results within the UCNA (equation (65)). The inset in (a) shows results for a 2D system for the same parameters.

Mean square velocity
The radial mean square velocity á ñ = á ñ u r u u 2 2 (63) for the potential (64) reads in terms of Δ. Evidently, this expression is independent of the Péclet number. Figure 6 shows simulation and theoretical results for the radial dependence of á ñ u v u 2 0 2 . Evidently, the mean square ABP velocity is constant over most of the potential range. Only at large radial distances, á ñ u u 2 drops rapidly. Consistent with the radial probability distribution y ( ) r , the simulation data drop to zero at s = r r m m determined by the condition (52) for the various Péclet numbers. Interestingly, the simulation results are independent of Δ. In contrast, the analytical expression (67) is independent of Pe, but depends on Δ. Hence, the approximate conditional velocity distribution function (62) (see also equation (A.12)), qualitatively reproduces the simulation result, but misses important dependencies on Pe and Δ.

Summary and conclusions
We have investigated the properties of ABPs confined in a radially symmetric potential. Two potentials have been considered, a harmonic and an anharmonic potential with a singularity at a finite radius. In the analytical approach, the AOUP model has been adopted, and the general time-dependent nonequilibrium distribution function for the harmonic potential has been determined. In case of the anharmonic potential, the stationarystate multidimensional distribution function within the UCNA has been considered [56]. In computer simulations, stationary-state distribution functions have been determined for ABPs in three dimensions.
A primary goal of our study is to clarify the extent to which the simulation results can be reproduced by the theoretical approaches. The study of ABPs and AOUPs confined in a harmonic potential clearly reveals the strong effect of the condition = = | | v v const. 0 on the stationary-state properties. Thereby, it is essential to compare appropriate distribution functions, which is the conditional probability distribution function y ( | ) r v 0 for AOUPs. We find good agreement between AOUP and ABP results for weak harmonic potentials  k 1 up to rather large Péclet numbers. For larger stiffness (  k 10), deviations appear when Pe exceeds the critical value Pe c of equation (42), where activity dominates over thermal fluctuations. However, for large potential stiffness very good agreement is achieved over a wide range of Péclet numbers, which is partially related to the increasing critical Péclet number, Pe c , with increasingk. In particular, for AOUPs, we find a distribution function y x ( | ) v 0 , where x~Pe r 2 , which is independent of the Péclet number in the limit  Pe Pe c . The dependence on the propulsion velocity is absorbed in the variable ξ. Only a dependence on D k remains, i.e., the potential strength and the rotational diffusion coefficient. Interestingly, a similar scaling is obtained for the radial distribution function of ABPs.
The good, at least qualitative, agreement between AOUP and ABP results reveals that the AOUP approach is a valuable approximation for an ABP system. The easier analytical tractability of the AOUP equations of motion renders the model advantages for an analytical description of the nonequilibrium statistical properties of active matter. Our studies of an anharmonic system within the UCNA yield good agreement with simulations for an AOUP. The probability distribution functions agree well for  D 0.1 and  Pe 100. Thereby, the agreement improves with decreasing Δ, i.e. increasing rotational diffusion coefficient. The distribution functions obtained from simulations of ABPs exhibit larger differences to the UCNA results; visible disparities appear already for Pe = 10 at D = 0.1, which increase with increasing Pe and Δ.
In general, the colored-noise correlation function applied for the UCNA is g á ñ= g - t R R with independent parameters g R and the (translational) diffusion coefficientD [55,56]. The relation to our studies is established by setting g =v D 3 0 2 R (see equation (11)) with the independent parameters g R and v 0 . The two different approaches lead to distinctly different regimes of applicability of the UCNA. As discussed in [55], the UCNA is suitable, i.e., the second derivative in equation (A.1) can be neglected, in the limits g  ¥ R and g  0 R . In contrast, for our independent parameters, only the limit g  ¥ R applies. For smaller g R , or larger Δ, the propulsion velocity v 0 , or the Péclet number, has to be small (see equation (5.37) of [55]). This is consistent with our presented results.
Our studies shed light onto the applicability of UCNA for ABPs. The good agreement with simulations for D  1 suggest that UCNA can very well be applied to describe the properties of ABPs with a large rotational diffusion coefficient, corresponding typically to run-and-tumble motion rather than thermal rotational diffusion. With the help of the stationary-state distribution function, other properties can be calculated, such as an effective free energy [39]. By the latter, the thermodynamic framework of passive systems can be applied and equations of states be derived [30].

Acknowledgments
Financial support by the Deutsche Forschungsgemeinschaft (DFG) within the priority program SPP 1726 'Microswimmers-from Single Particle Motion to Collective Behavior' is gratefully acknowledged.
Appendix. UCNA-distribution functions A.1. Multidimensional stationary-state positional distribution function Neglecting the random force G, the time derivative of equation (1) and insertion of equation (9)