Brownian systems with spatially inhomogeneous activity

TL;DR

This paper extends the Green-Kubo method to analyze inhomogeneous activity in Brownian systems, accurately predicting spatial variations in orientation and density, validated by simulations.

Contribution

It introduces a generalized Green-Kubo approach for inhomogeneous active particles and combines it with a Gaussian approximation and density functional theory.

Findings

Accurate analytical expression for average orientation.

Good agreement between theory and Brownian dynamics simulations.

Effective prediction of spatial density distribution.

Abstract

We generalize the Green-Kubo approach, previously applied to bulk systems of spherically symmetric active particles [J. Chem. Phys. 145, 161101 (2016)], to include spatially inhomogeneous activity. The method is applied to predict the spatial dependence of the average orientation per particle and the density. The average orientation is given by an integral over the self-part of the van Hove function and a simple Gaussian approximation to this quantity yields an accurate analytical expression. Taking this analytical result as input to a dynamic density functional theory approximates the spatial dependence of the density in good agreement with simulation data. All theoretical predictions are validated using Brownian dynamics simulations.

Equations

$${\dot{\mathbf{r}}}_i=v_0({\mathbf{r}}_i,t){\mathbf{p}}_i+{\gamma }^{-1}{\mathbf{F}}_i+{\mathbf{\xi }}_i,{\dot{\mathbf{p}}}_i={\mathbf{\eta }}_i\times {\mathbf{p}}_i,$$

$$\frac{\partial P(t)}{\partial t}={\Omega }_{\mathrm{a}}(t)P(t),$$

$$\displaystyle\Omega_{\rm eq}=\sum_{i=1}^{N}\boldsymbol{\nabla}_{i}\!\cdot\!\big{[}D_{\text{t}}\!\left(\boldsymbol{\nabla}_{i}\!-\beta\boldsymbol{F}_{i}\right)\big{]}\!+\!D_{\text{r}}\boldsymbol{R}_{i}^{2},$$

$$P(t)=P_{\mathrm{eq}}+{\int }_{-\infty }^{t}d{t}^{\prime }{e}_{+}^{{\int }_{t^{\prime }}^{t}ds{\Omega }_{\mathrm{a}}(s)}\delta {\Omega }_{\mathrm{a}}(t)P_{\mathrm{eq}},$$

$$K(t)$$

$$V(t)$$

$$⟨f⟩(t)=⟨f{⟩}_{\mathrm{eq}}-{\int }_{-\infty }^{t}dt⟨G(t^{\prime })e_{-}^{{\int }_{t^{\prime }}^{t}ds{\Omega }_{\mathrm{a}}^{†}(s)}f{⟩}_{\mathrm{eq}},$$

$$\mathbf{p}(\mathbf{r})=\frac{⟨{\sum }_i\delta (\mathbf{r}-{\mathbf{r}}_i){\mathbf{p}}_i⟩}{\rho (\mathbf{r})},$$

$$\mathbf{p}(\mathbf{r})=-{\int }_0^{\infty }dt\frac{e^{-2D_{r}t}}{3{\rho }_b}{⟨\sum _i{\nabla }_i{v}_0({\mathbf{r}}_i)e^{{\Omega }_{\mathrm{eq},\mathrm{tr}}^{†}t}\delta (\mathbf{r}-{\mathbf{r}}_i)⟩}_{\mathrm{eq},\mathrm{s}}$$

$${⟨\sum _i{v}_0({\mathbf{r}}_i){\mathbf{F}}_i\cdot e^{{\Omega }_{\mathrm{eq},\mathrm{tr}}^{†}t}\delta (\mathbf{r}-{\mathbf{r}}_i)⟩}_{\mathrm{eq},\mathrm{s}}=0.$$

$$\mathbf{p}(\mathbf{r})=-{\int }_0^{\infty }dt\frac{e^{-2D_{r}t}}{3}\int d{\mathbf{r}}^{\prime }\nabla v_0({\mathbf{r}}^{\prime })G_{\mathrm{vH}}^{\mathrm{s}}(|\mathbf{r}-{\mathbf{r}}^{\prime }|,t).$$

$$\mathbf{p}(\mathbf{r})$$

$$\chi (|\mathbf{r}-{\mathbf{r}}^{\prime }|,t)=\frac{e^{-2D_{r}t}}{3}\nabla G_{\mathrm{vH}}^{\mathrm{s}}(|\mathbf{r}-{\mathbf{r}}^{\prime }|,t).$$

$$G_{\mathrm{vH}}^{\mathrm{s}}(\mathbf{r},t)=\frac{1}{(4\pi D_{t}t{)}^{3/2}}{e}^{-\frac{r^2}{4D_{t}t}}.$$

$$\frac{\partial \rho (\mathbf{r},t)}{\partial t}$$

$$+D_t\nabla \cdot \left(\rho (\mathbf{r},t)\nabla \frac{\delta \beta \mathcal{F}[\rho ]}{\delta \rho (\mathbf{r},t)}\right),$$

$$\nabla \frac{\delta \beta \mathcal{F}[\rho ]}{\delta \rho (\mathbf{r},t)}=\frac{\nabla \rho (\mathbf{r},t)}{\rho (\mathbf{r},t)}-\nabla \int d{\mathbf{r}}^{\prime }{c}^{(2)}(|\mathbf{r}-{\mathbf{r}}^{\prime }|)\tilde{\rho }({\mathbf{r}}^{\prime },t),$$

$$\mathbf{p}(z)=-\frac{\nabla v_0(z)}{3(2D_r+D_t{\omega }^2)}=-{\hat{\mathbf{e}}}_{\mathbf{z}}\frac{v_a\omega \cos (\omega z)}{3(2D_r+D_t{\omega }^2)}.$$

Full text

active colloids, phase separation, wetting

Active Brownian particles (ABPs) are intrinsically nonequilibrium systems, constantly driven out of equilibrium by consuming energy from the local environment. Nonequilibrium statistical mechanics aims to calculate from the microscopic dynamics the relevant average quantities. However, this presents a difficult theoretical problem, even for simple model systems. One way to make systematic progress is to restrict attention to the linear response regime, where the problem becomes tractable. Within linear response, formally exact results for system averages can be obtained by integrating the time-correlation functions of the corresponding passive (equilibrium) system Green (1952); Kubo (1957).

We have recently applied the linear response (Green-Kubo) approach to a homogeneous system of ABPs Sharma and Brader (2016), focusing our attention on calculation of the average swim speed Cates and Tailleur (2015); Krinninger et al. (2016), a central quantity appearing in coarse-grained theories of active matter. We have demonstrated that the average swim speed, which describes how the motion of each particle is obstructed by its neighbours, can be obtained from a history integral over the equilibrium autocorrelation of tagged-particle force fluctuations Sharma and Brader (2016). The theory was tested using active Brownian dynamics simulations and provides a solid basis for the development of first-principles theoretical approaches.

Bulk systems have been the focus of much attention, due largely to the phenomenon of motility-induced phase separation, and several experimental studies of bulk ABPs have been performed Howse et al. (2007); Jiang et al. (2010); Palacci et al. (2013); Gomez-Solano et al. (2016); Palacci et al. (2010). However, there exist synthetic Saha et al. (2014); Hong et al. (2007); Magiera and Brendel (2015); Stenhammar et al. (2016) and living systems Berg et al. (1972); Jékely et al. (2008); Hoff et al. (2009); Jikeli et al. (2015) for which the propulsion strength is not a global constant, but dependent on the spatial location of the particles. For example, in a recent experimental study of synthetic microswimmers Lozano et al. (2016) position-dependent motility was implemented using an inhomogenous laser field, resulting in phototaxis. Position-dependent activity also features in the energy depot model Schweitzer et al. (1998).

Motivated by these considerations, we consider in this paper active systems for which the particle propulsion-speed varies in space. Our aim is to extend the Green-Kubo approach, previously applied in bulk Sharma and Brader (2016), to treat systems with inhomogeneous activity, thus providing a first-principles theoretical route to addressing more realistic situations. The physical observables to be considered are the average orientation, which features prominently in studies of inhomogenous systems Enculescu and Stark (2011); Lozano et al. (2016), and the density. We will show that both of these quantities become inhomogeneous in the presence of spatially varying activity and we calculate these explicitly for a simple test case. Our formalism also makes clear that an external activity field affects the system in a qualitatively different way than an external potential field; the former generates a linear response for the average orientation, whereas the latter does not.

In our Green-Kubo approach we employ a variation of the integration-through-transients approach, originally developed for treating interacting Brownian particles subject to external flow Morriss and Evans (2007); Fuchs and Cates (2002, 2005); Brader et al. (2012). We find that the average orientation is proportional to the local gradient of the activity field, even in the absence of a one-body torque. The relevant autocorrelation function is the well-known self part of the van Hove function Hansen and McDonald (1990), which can be very well approximated by a Gaussian. Within the Gaussian approximation a simple and accurate analytical expression can be obtained for the average orientation. Taking this as input to a Dynamic Density Functional theory we then proceed to develop a closed theory for the inhomogenous density. Our predictions are tested against data from Brownian dynamics simulations.

We consider a three dimensional system of $N$ active, interacting, spherical Brownian particles with coordinate $\boldsymbol{r}_{i}$ and orientation specified by an embedded unit vector $\boldsymbol{p}_{i}$. A space- and time-dependent self-propulsion of speed $v_{0}(\boldsymbol{r}_{i},t)$ acts in the direction of orientation and $i$ labels the particle. Omitting hydrodynamic interactions the motion can be modelled by the Langevin equations

[TABLE]

where $\gamma$ is the friction coefficient and the force on particle $i$ is generated from the total interparticle interaction energy according to $\boldsymbol{F}_{i}\!=\!-\nabla_{i}U_{N}$. For clarity of presentation we do not include an external potential field. The stochastic vectors $\boldsymbol{\xi}_{i}(t)$ and $\boldsymbol{\eta}_{i}(t)$ are Gaussian distributed with zero mean and have time correlations $\langle\boldsymbol{\xi}_{i}(t)\boldsymbol{\xi}_{j}(t^{\prime})\rangle=2D_{t}\boldsymbol{1}\delta_{ij}\delta(t-t^{\prime})$ and $\langle\boldsymbol{\eta}_{i}(t)\boldsymbol{\eta}_{j}(t^{\prime})\rangle=2D_{r}\boldsymbol{1}\delta_{ij}\delta(t-t^{\prime})$. The translational and rotational diffusion coefficients, $D_{t}$ and $D_{r}$, are treated as independent parameters. Note that in the second of Eqs. (1) the orientation vector does not couple to the activity field and hence no direct torque acts on the particle due to the position-dependent activity.

It follows exactly from (1) that the joint probability distribution, $\!P(\boldsymbol{r}^{N}\!\!,\boldsymbol{p}^{N}\!\!,t)$, evolves according to Gardiner et al. (1985)

[TABLE]

where $\Omega_{\rm a}$ is the time-evolution operator. We have used $\!P(\boldsymbol{r}^{N}\!\!,\boldsymbol{p}^{N}\!\!,t)\equiv P(t)$, and $\Omega_{\rm a}(\boldsymbol{r}^{N}\!\!,\boldsymbol{p}^{N}\!\!,t)\equiv\Omega_{\rm a}(t)$ to keep the notation compact. The time-evolution operator can be split into a sum of two terms, $\Omega_{\rm a}(t)=\Omega_{\rm eq}+\delta\Omega_{\rm a}(t)$, where the equilibrium contribution is given by

[TABLE]

with rotation operator $\boldsymbol{R}\!=\!\boldsymbol{p}\times\!\nabla_{\!\boldsymbol{p}}$ Morse et al. (1953) and $\beta=1/(k_{\rm B}T)$. Using $\Omega_{\rm eq}P_{\rm eq}=0$, we obtain a formal solution for the nonequilibrium distribution

[TABLE]

where $e_{+}$ is a positively ordered exponential function Brader et al. (2012). The active part of the dynamics is described by the operator $\delta\Omega_{\rm a}=-\sum_{i}\boldsymbol{\nabla}_{i}\!\cdot(v_{0}(\boldsymbol{r}_{i},t)\boldsymbol{p}_{i})$. The action of this operator on $P_{\rm eq}$ yields $\delta\Omega_{\rm a}(t)P_{\rm eq}=-P_{\rm eq}\left(K(t)+V(t)\right)$, where we have defined the quantities $K(t)$ and $V(t)$ as

[TABLE]

We obtain from Eqs. (4), (5) and (6) an exact expression for the nonequilibrium average of a test function $f\equiv f(\boldsymbol{r}^{N}\!\!,\boldsymbol{p}^{N}\!)$ as

[TABLE]

where we have defined $G(t)=K(t)+V(t)$ and the adjoint operator is given by ${\Omega}^{\dagger}_{\rm a}(t)={\Omega}^{\dagger}_{\rm eq}-\delta\Omega_{\rm a}(t)$, where ${\Omega}^{\dagger}_{\rm eq}=\sum_{i}D_{\text{t}}\!\left(\boldsymbol{\nabla}_{i}\!+\beta\boldsymbol{F}_{i}\right)\!\cdot\!\boldsymbol{\nabla}_{i}\!+\!D_{\text{r}}\boldsymbol{R}_{i}^{2}$. The integrand appearing in Eq. (7) involves the equilibrium correlation between $G$ at time $t^{\prime}$ and the observable $f$ which evolves from $t^{\prime}$ to $t$ according to the full dynamics. From here onwards we will consider only the linear response, obtained by replacing the full time-evolution operator in (7) by the time-independent equilibrium adjoint operator.

We will focus first on calculating to linear order in activity the average orientation per particle, defined as

[TABLE]

where $\rho(\boldsymbol{r})=\left\langle\sum_{i}\delta(\boldsymbol{r}-\boldsymbol{r}_{i})\right\rangle$ is the one-body density. We will henceforth assume that the activity does not vary in time, $v_{0}(\boldsymbol{r}_{i},t)=v_{0}(\boldsymbol{r}_{i})$; generalization to time-dependent situations is straightforward.

Our first remark is that any observable independent of the particle orientation vectors, $\boldsymbol{p}_{i}$, does not admit a linear response in $v_{0}(\boldsymbol{r})$. As the function $G(t)$ appearing in Eq. (7) is linear in $\boldsymbol{p}_{i}$ the angular integrals in the equilibrium average will yield zero by symmetry. Only odd functions of $\boldsymbol{p}_{i}$ will generate a nonzero linear-response. As the density is independent of $\boldsymbol{p}_{i}$ we can thus replace $\rho(\boldsymbol{r})$ in Eq. (8) by the bulk number density $\rho_{b}$ when working to first order in $v_{0}$.

The average orientation per particle is obtained by using Eq. (7) to evaluate the numerator of Eq. (8) to linear order in $v_{0}$. The equilibrium time evolution operator can be split into translational and rotational contributions, ${\Omega}^{\dagger}_{\rm eq}={\Omega}^{\dagger}_{\rm eq,r}+\,{\Omega}^{\dagger}_{\rm eq,tr}$. These operators commute (the rotational part acts only upon the particle orientation), which allows the angular integrals to be evaluated explicitly. This yields

[TABLE]

where $\langle..\rangle_{\rm eq,s}$ denotes an average over the spatial degrees of freedom. In deriving Eq. (9) we have employed two results: Firstly, the orientation decorrelates according to $\langle\boldsymbol{p}_{i}e^{{\Omega}^{\dagger}_{\rm eq,r}t}\boldsymbol{p}_{j}\rangle_{\rm eq,r}\!\!=\delta_{ij}\boldsymbol{1}e^{-2D_{r}t}/3$, where the equilibrium average is over rotational degrees of freedom (see the supplementary material). Secondly, in homogeneous equilibrium the interaction force on a particle is not correlated with its position

[TABLE]

Using $\boldsymbol{\nabla}_{i}v_{0}(\boldsymbol{r}_{i})=\int d\boldsymbol{r}\,\delta(\boldsymbol{r}-\boldsymbol{r}_{i})\boldsymbol{\nabla}v_{0}(\boldsymbol{r})$ and introducing the self-part of the equilibrium van Hove function Hansen and McDonald (1990), $G_{\rm vH}^{\rm s}(|\boldsymbol{r}-\boldsymbol{r}^{\prime}|,t)=\langle\delta(\boldsymbol{r}^{\prime}-\boldsymbol{r}_{1})e^{{\Omega}^{\dagger}_{\rm eq,tr}t}\delta(\boldsymbol{r}-\boldsymbol{r}_{1}\rangle_{\rm eq,s}$, enables us to express the average orientation in the compact form

[TABLE]

The average orientation of particles is antiparallel to the gradient of the activity field and has a magnitude determined by the equilibrium self van Hove function. This result can be rewritten in the alternative form

[TABLE]

where the space-time response function, $\chi(|\boldsymbol{r}-\boldsymbol{r}^{\prime}|,t)$, is given by

[TABLE]

The right hand side of Eq. (13) is simply the functional derivative of $\boldsymbol{p}(\boldsymbol{r})$ with respect to $v_{0}(\boldsymbol{r}^{\prime})$, evaluated at zero activity. Eqs. (12) and (13) are the key linear response results of this paper.

Recent experiments have shown that active particles tend to orient in an inhomogenous activity field Bickel et al. (2014); Lozano et al. (2016). In these studies, it was identified that the orientation is a consequence of an aligning torque acting on the particles. Within a first-principles theoretical approach, such a one-body torque would have to be explicitly incorporated into the orientational member of the Langevin equations (1). While it is clear that deterministic one-body torques can generate an average orientation, it is much less obvious that this will emerge from the present, torque-free Langevin equations as a purely statistical phenomenon. It is this aspect which is of primary interest in the present study. From Eq. (11) it is evident that ABPs do indeed tend to orient in an inhomogenous activity field, even in the absence of deterministic aligning torques. Moreover, our numerical results will demonstrate that this can be a large effect within certain parameter ranges. A notable feature of the present orientation mechanism is that it is independent of the particle diameter, in contrast to the torque-based mechanism identified in Refs. Bickel et al. (2014); Lozano et al. (2016).

We can use the linear response theory to distinguish the orientational response of ABPs to inhomogeneous activity from the response to an external potential field (with spatially constant activity), e.g. the sedimentation of ABPs under gravity Enculescu and Stark (2011). In the former situation the leading order contribution to $\boldsymbol{p}(\boldsymbol{r})$ is linear in $v_{0}$, as is evident from Eq. (12), whereas in the latter situation the leading order is quadratic. For a system with constant $v_{0}$ and external potential $v_{\rm ext}(\boldsymbol{r})$ the function $G(t)$ appearing in Eq. (7) is replaced by the time-independent function $G=K+D_{t}\beta\sum_{i}(\beta{\bf F}_{i}\cdot\nabla_{i}v_{\rm ext}(\boldsymbol{r}_{i})+\nabla^{2}_{i}v_{\rm ext}(\boldsymbol{r}))$. Use of Eq. (7) to calculate $\boldsymbol{p}(\boldsymbol{r})$ to linear order in $v_{0}$ thus yields three terms, all of which are zero; the first vanishes due to Eq. (10) and the others due to the symmetry of the angular integrands.

Given the exact linear response result Eq. (12), it is desirable to obtain from this a closed theory by approximating the self part of the van Hove function. A commonly employed approximation is the Gaussian Hansen and McDonald (1990)

[TABLE]

This approximation is known to be accurate, at least for hard-spheres, up to dimensionless densities as high as $\rho_{b}\!=\!0.6$ Hopkins et al. (2010). We will show below that this approximation proves very reliable for calculating the average orientation profiles of repulsive ABPs and, for certain choices of $v_{0}(\boldsymbol{r})$, enables evaluation of the integrals in Eq. (11) to yield an explicit analytical expression.

As discussed above, for reasons of symmetry the density remains unaltered from that in bulk to linear order; the first modification is quadratic in $v_{0}$. It is in principle straightforward to expand the exact expression in Eq. (7) to second order in activity (noting that $\Omega_{\rm a}^{\dagger}$ is a function of $v_{0}$), however, the formal expression thus generated for the density response is difficult to evaluate and yields little insight. We thus follow an alternative, simpler route to obtain the density approximately.

A coarse-grained expression for the density can be obtained by integrating Eq. (2) over all orientational and all but one translational degrees of freedom. However, this generates a term involving the two-body density. The dynamical density functional theory (DDFT) approximates this unknown term using an equilibrium free energy functional. In the case of passive particles this is sufficient to yield a closed theory Hopkins et al. (2010). Applying this procedure to Eq. (2) yields the following DDFT

[TABLE]

where $\mathcal{F}[\rho]$ is the Helmholtz free energy functional. We observe that Eq. (Brownian systems with spatially inhomogeneous activity) involves two unknown functions, $\rho(\boldsymbol{r},t)$ and $\boldsymbol{p}(\boldsymbol{r},t)$, and is thus not closed. Our strategy is to use as input to the DDFT the linear response result Eq. (12) with the Gaussian approximation Eq. (14). The resulting equation can then be solved self-consistently for the steady-state density. Note that the activity appears quadratically in the DDFT equation, such that the density is independent of the sign of $v_{0}$. It now only remains to specify the free energy functional.

The Helmholtz free energy can be split into two contributions, $\mathcal{F}=\mathcal{F}_{\rm id}+\mathcal{F}_{\rm ex}$. The ideal part is given exactly by $\mathcal{F}_{\rm id}[\rho]=k_{B}T\!\int d\boldsymbol{r}\rho(\boldsymbol{r},t)\left(\log(\Lambda^{3}\rho(\boldsymbol{r},t))-1\right)$, where $\Lambda$ is the thermal wavelength. The excess part, $\mathcal{F}_{\rm ex}$, encodes the interparticle interactions and is, in general, unknown. A commonly employed approximation, sufficient for our present purposes, is to use a functional Taylor expansion, truncated at second order in the density. Within this approximation the gradient of the functional derivative entering Eq. (Brownian systems with spatially inhomogeneous activity) is given by

[TABLE]

where $\tilde{\rho}(\boldsymbol{r},t)=\rho(\boldsymbol{r},t)-\rho_{b}$ and $c^{(2)}(r)$ is the bulk direct correlation function, easily obtained from standard liquid state integral equation theory. We employ here the Percus-Yevick integral equation Hansen and McDonald (1990).

Equations (12), (13) and (14) enable calculation of the average orientation per particle. Equations (Brownian systems with spatially inhomogeneous activity)-(16) take this as input and yield the density. We will next present numerical results and compare these with data from active Brownian dynamics simulation. The simulations are performed on a three-dimensional system of $N=500$ particles interacting via the pair-potential $\beta u(r)=4\varepsilon((\sigma/r)^{12}-(\sigma/r)^{6})$, where $\sigma$ sets the length scale and we set $\varepsilon=1$. The potential is truncated at its minimum, $r=2^{1/6}\sigma$ to yield a softly repulsive interaction. The system size $L$ is determined as $L\!=\!(N/\rho_{\rm b})^{1/3}$ in order to obtain the desired density. Periodic boundary conditions are applied in all three directions. The integration time step is fixed to $dt\!=\!10^{-5}\tau_{B}$ where $\tau_{B}\!=\!d^{2}/D_{t}$ is the time-scale of translational diffusion. The equation for time evolution of the orientation vector (Eq. (1)) is evaluated as an Ito Integral. Measurements are made after a time $20\tau_{B}$ to ensure equilibration. We choose the ratio of diffusion coefficients as $D_{r}/D_{t}\!=\!25$.

As a test case we impose the activity $v_{0}(z)=v_{a}\sin(\omega z)$ varying only in the $z$-direction, where $v_{a}$ and $\omega$ are parameters. Inserting this choice and the Gaussian approximation into Eq. (11) yields a simple theoretical prediction for the average orientation per particle

[TABLE]

In Fig. 1, we plot the average orientation per particle for two different bulk densities. For low density ($\rho_{b}=0.2$ in Fig. 1(a) and (b)), the analytical prediction of Eq. (17) is in good agreement with the numerics. At high density (Fig. 1(c) and (d)), the theory provides a slight overestimation. As expected, the average orientation increases with both the magnitude $v_{a}$ and angular frequency $\omega$ of the activity field. As can be seen in Fig. 1, the average orientation $\boldsymbol{p}(z)$ can attain significant values ($\approx 0.2$). This is despite the fact that the rotational diffusion occurs on a time scale much smaller than the translational diffusion ($D_{r}/D_{t}=25$). Since any activity field can be decomposed into a Fourier series, one can obtain the average orientation for a generic activity field. The only obvious limitation of the approach is that the theory predicts the linear order response and hence Eq. (17) is expected to be valid only for small activities.

We calculate the density by using Eqs. (16) and (17) as inputs to Eq. (Brownian systems with spatially inhomogeneous activity). In Fig. 2 we plot the theoretical prediction together with the numerically measured relative change in density for two different bulk densities. There are three noteworthy features: (1) the density shows peaks at the nodes of the activity field, (2) the change in density becomes increasingly asymmetric with activity and (3) the relative change in density decreases with increasing initial bulk density. The theoretical prediction is consistent with these features. However, the theory predicts a more symmetric change in density than observed in simulations. Nevertheless, it is remarkable that a fundamentally second-order effect, i.e., the change in density, can be calculated theoretically to a good accuracy using a simple adiabatic DDFT approach. We note that in the absence of interparticle interactions Eq. (Brownian systems with spatially inhomogeneous activity) can be easily solved in steady-state to yield a universal curve for the relative density change. Given that the density change is as large as $30\%$ in Fig. 2(a) for $v_{a}=20$, it is expected that nonlinear deviations appear in $\boldsymbol{p}(z)$. This is apparent in Fig. 1(a) which corresponds to the same set of parameters as Fig. 2(a). In this paper, we have not performed a systematic study of the range of validity of the linear response. Selective cases, with large activity, for which strong nonlinear deviations are observed in the system’s reponse are presented in the supplementary material.

The underlying mechanism of orientation and density change can be understood qualitatively as follows: A particle with an orientation antiparallel to the activity gradient experiences slowing down as it moves in the direction of the gradient whereas a particle with orientation parallel to the gradient speeds up. This asymmetric influence of the activity, therefore, leads to accumulation of particles in proportion to the magnitude of the local gradient of activity.

To summarize our main findings: we have derived a formally exact expression (Eq. (7)) for calculating averages in a system of interacting Brownian particles, subject to a position-dependent activity $v_{0}(\boldsymbol{r})$. From this we obtain the linear-response of the average orientation in Eq. (12) and identify the relevant time-correlation function as the self-part of the van Hove function. We find that linear response provides an accurate account of the orientiation $\boldsymbol{p}(\boldsymbol{r})$ over a significant parameter range. Taking the analytical prediction for $\boldsymbol{p}(\boldsymbol{r})$ as an input to the dynamic density functional theory, we can obtain the spatial dependence of the density in good agreement with simulation data. Our approach is perfectly suited to obtain the time-dependent response of system subjected to a space-time dependent activity. It will be interesting to study the response $\boldsymbol{p}(\boldsymbol{r},t)$ and the corresponding time-evolution of density for a time- and space-dependent activity.