Fergus Moore

John Russo

Tanniemola B Liverpool

C Patrick Royall

Active Brownian Particles in Random and Porous Environments

The transport of active particles may occur in complex environments, in which it emerges from the interplay between the mobility of the active components and the quenched disorder of the environment. Here we explore structural and dynamical properties of Active Brownian Particles (ABPs) in random environments composed of fixed obstacles in three dimensions. We consider different arrangements of the obstacles. In particular, we consider two particular situations corresponding to experimentally realizable settings. Firstly, we model pinning particles in (non-overlapping) random positions and secondly in a percolating gel structure, and provide an extensive characterization of the structure and dynamics of ABPs in these complex environments. We find that the confinement increases the heterogeneity of the dynamics, with new populations of absorbed and localized particles appearing close to the obstacles. This heterogeneity has a profound impact on the motility induced phase separation (MIPS) exhibited by the particles at high activity, ranging from nucleation and growth in random disorder to a complex phase separation in porous environments.

I. INTRODUCTION

Active matter concerns systems comprised of individual bodies undergoing motion via self-propulsion 1 . This description encompasses a plethora of systems over a wide range of lengthscales from bacteria [2][3][4] , biological microswimmers 5 , schools of fish 6 , bird flocks 7 , to human crowds 8 . While these systems can all be classified as active matter, accurate models must be tailored to the specifics of each system and its environment. One class of much simpler model systems which nevertheless captures key elements of the behavior of more complex systems is active colloids [9][10][11][12] .

Yet if we are to apply such model systems in a biological context, it is essential that we study the dynamics of active particles in environments that are relevant to their real-world counterparts. For biological active matter on mesoscopic lengthscales this means environments like porous soils 13 and organic tissues 14 . Environments such as these have several qualities in common, they are often crowded, random and irregular. This of course has an impact on the transport or displacement of the active bodies inside these spaces 9 , for example in the (biological) case of bacteria, their run-and-tumble dynamics can be drastically altered 15,16 .

In equilibrium systems, the dynamics of fluids in dense and complex environments has long been an area of interest. For example, the inclusion of specific structures into dense fluids has proven significant for progress towards a) Electronic mail: john.russo@uniroma1.it b) Electronic mail: paddy.royall@espci.psl.eu understanding the glass transition and in liquids 17,18 . Furthermore, the addition of randomly pinned particles within a dense ensemble is known to greatly slow the dynamics of such systems, providing access to rare states [19][20][21][22] . A special kind of localisation has been explored in glass-forming systems 23,24 , and which may be realized using colloidal systems [25][26][27] . Here pinning, i.e. immobilizing a fraction of the particles, provides access to the so-called ideal glass, a putative amorphous state of very low configurational entropy whose diverging timescales render it otherwise inaccessible to experiment or computer simulation 28 .

For active systems, in the simplest case of a confining wall, self-propelling spheres will accumulate at a wall as a consequence of the timescale of its persistent motion, even in the absence of hydrodynamic interactions 29,30 . For many-body systems in two-dimensions, the influence of disordered landscapes on to the dynamics of active systems has been shown to manifest in clogging and localisation transitions [31][32][33] , subdiffusion over long timescales [34][35][36] , destruction of flocking clusters 37 , suppression of Motility-Induced Phase Separation (MIPS) and the prevention of uniform wetting at boundaries 38 . Furthermore, manipulation of complex environments has been shown to provide a degree of control over the transport of active matter in the form of sorting 39 , and the intriguing phenomenon of topotaxis (control over net flow directions by controlling the topology of the environment 40,41 ). Active Brownian particles (ABPs) exhibit rich phase behavior such as MIPS 42,43 , and the formation of active crystals [START_REF] Omar | Active Crystallization[END_REF][START_REF] Moore | [END_REF] and fundamental properties such as pressure and the equation of state differ drastically from what might be expected from equilibrium systems [46][47][48] .

However, the question of how active Brownian spheres couple to a complex surrounding environment remains unanswered. Recently, there is been interest in experiments with mesoscale active matter in 3d. In one study a random heterogeneous environment was found to impose strong inhibitions on active transport of bacteria 15 , and in another study the impact of dimensionality was made clear, with the dominance of three-dimensional structure in the presence of an anisotropic potential 49 . This latter example is a well-controlled 3d colloidal model system which provides the inspiration for our work because it is possible to confine such systems using pinning 26,50 or allowing a subset of particles to undergo gelation [START_REF] Varrato | Proc. Nat. Acad. Sci[END_REF] . Insight into the transport of active matter in 3d complex environments could provide a major step towards control of such systems and aid progress towards applications such as drug delivery.

In the present article, we perform three-dimensional molecular dynamics simulations of active Brownian particles in complex heterogeneous environments. To model such environments, we choose two example structures which, as noted above, may be realized in experiment: a random homogeneous array of pinned particles, providing an extension of disordered random obstacle studies to 3D; and a continuous, percolating porous network (a gel), simulating the complex environments typical of active matter under confinement. Furthermore, we will investigate the structural and dynamical properties of ABPs within these structures, with a focus on phase separation and how this varies from the MIPS observed in bulk suspensions.

This article is organized as follows: in section II we outline the computational methods used to study these systems, in section III we discuss the results of the simulations, and finally in section IV we will summarise and conclude our findings and discuss future work in this area.

II. MODEL AND METHODS

A. Active Particles

We model active colloids as active Brownian particles, which propel with a constant velocity V 0 along their individual direction vectors e, which in turn are subject to rotational diffusion. We implement this model through molecular dynamics simulations using a customized version of the open source LAMMPS package [START_REF] Moore | [END_REF][START_REF] Plimpton | [END_REF] , which integrates the following equations of motion:

ṙ = V 0 e + βD T F + 2D T η (1) ė = 2D R ξ × e (2)
Here ṙ is the particle velocity, V 0 is the magnitude of the constant active velocity, and F is the inter-particle force. The thermal fluctuations promoting translational diffusion are included in the Gaussian white-noise term η, where η = 0, and D T is the translational diffusion coefficient. Thermal noise driving rotational diffusion of the direction vector e is represented by an independent Gaussian noise term ξ, where ξ = 0, and D R is the rotational diffusion coefficient. The two diffusion coefficients are related via D T = D R σ 2 /3, where σ is the particle diameter. Time is scaled in units of the characteristic rotational diffusion time τ R = 1/(2D R ) 42 . The active particles are modelled as being similar to hard spheres and to achieve this we include a Weeks-Chander-Andersen (WCA) inter-particle potential, which takes the form:

βu(r ij ) =    4βε σ rij 12 -σ rij 6 + βε r ij ≤ 2 1 6 σ 0 r ij > 2 1 6 σ
(3) where ε = 5 is the interaction strength r ij is the interparticle distance, and β = 1/k B T the thermal energy. Since we use the WCA interaction, it is hard to define a volume fraction. Methods that determine an effective diameter such as Barker-Henderson 53 , may not hold outside of equilibrium systems. Therefore, as in ref. 54 , we use the total density ρ = N/V , where N is the number of particles and V is the volume of the system.

We use the Péclet number to refer to the relative strength of the activity in the system, which we define as: Pe = V 0 /σD R . Throughout this work we keep D R constant via D T = 1, and vary Pe by changing the propulsion velocity V 0 .

Simulations are performed with periodic boundary conditions. The majority of the work is carried out in a cubic box of dimension length L = 55σ, with a total number of N = 144000 particles. In some cases, there was a need to sample from many state points and for these a smaller system was used where L = 27.5σ and N ranges from 18000 to 24000. Analysis at constant density is always conducted at ρ = 0.87. This state point is chosen such that it lies below the freezing line in the bulk.

B. Complex Environments

The complex environments relevant to microscopic biological systems are often irregular and random in nature. To investigate the dynamical properties of ABPs under these conditions, one must prepare obstacle geometries that satisfy these requirements. Here we consider two primary structures; porous gel networks and randomly pinned particles. As noted above, these may be realized in experiments In addition to these, we include simulations studying the bulk dynamics of ABPs as a reference, these bulk simulations use the approach outlined in the previous section.

In the following section, we describe how the two complex environments are created and characterized. A schematic depicting these two environments is displayed in Fig. 1, featuring 3d renderings of each environment type along with a cross-sectional slice. The porous gel network (Fig. 1a-b), is a heterogeneous system, comprised of two distinct meso-phases which percolate through the entire simulation box comprised of the particle-rich phase and a particle-poor phase in which for our parameters no particles are found. The random environment (Fig. 1c-d), is comprised of randomly pinned particles that create number of discrete obstacles dispersed throughout the system.

Preparation of a porous network-In this work the porous network is modeled as a colloidal gel, specifically a colloid-polymer mixture 55,56 . We do this because we seek to connect our work to experiments where such a network could be realized 49,56 . Therefore, the preparation protocol in our simulations is as close to one which might be experimentally realizable as possible. For suitable parameters, colloids with such an attractive interaction begin to phase separate via spinodal decomposition which is then arrested, leaving a bicontinuous network, ie a gel 56,57 . The (polymer-induced) attraction between the colloidal particles that would be used in an experiment is here modeled with the Morse potential 58,59 ,

βu (r ij ) = βε exp [a 0 (σ -r ij )] (exp [a 0 (σ -r ij )] -2)
(4) where a 0 = 33 is a range parameter.

To create the gel structures in simulation, we begin with particles in a simple cubic crystal at the desired number density, and then evolve this system according to Brownian dynamics (Eq. ( 1), for V 0 = 0), with the particles interacting via the Morse potential Eq. (4). This system is then evolved for 5 × 10 7 integration steps, which is equivalent to 1200τ B after which the system is frozen and no further movement is allowed. Here τ B is the Brownian time τ B = (σ/2) 2 /6D T . An example of the resulting gel is shown in Fig. 1(a-b).

With the gel in place, the positions of the free particles are initialized via the Lubachevsky-Stillinger algorithm 60 . This method comprises the following steps: First, initial particle positions are randomly assigned. Then, the particles are slowly grown and displaced from an initial diameter σ in = 0.1, to the desired diameter σ = 1, such that these particles experience minimal overlaps with themselves and with the frozen gel particles. Additionally, to guard against the presence of any small but sufficient remaining particle overlaps, a pre-run simulation with a soft potential is performed:

u(r ij ) = A 1 + cos πr ij r c r ij < r c (5)
where r c = 2 1 6 is the potential cut-off, and the constant A is ramped from 0 to 100 over 1.2τ R without activity. Following this the system is equilibrated again without activity, for particles following the equations of motion outlined in Eq. (1) and Eq. ( 2), and the WCA interparticle potential Eq. (3). A low density example of this system is shown in Fig. 1(e-f).

Random pinning -For the random pinning case, an arbitrary configuration of particles is generated in the simulation box at the desired density. These particles then follow Brownian dynamics with the soft potential Eq. ( 5) to eliminate any significant overlaps before the system is equilibrated with the WCA potential. A fraction of the particles are then chosen at random and frozen. This creates the random obstacles. This fraction is chosen such that the volume accessible to the free particles is the same in both the gel network and the random pinning systems (Fig. 1(c-d)).

The fixing of the accessible free volume of the mobile particles enables comparison of observables in both environments. Fixing this is a necessary step as the difference in structure between the gel network and the random pins could result in the free particles operating at two different effective densities in the case of the same number of frozen particles.

For this, the density of the gel network is held constant and the number of pinned particles is varied as a function of the total density ρ. The total free volume available to the mobile particles is determined by taking the Voronoi tessellation of the instantaneous configuration and associating with each particle the volume of its Voronoi cell V i voro . The sum of these volumes provides the total volume accessible to the free particles

N i V i
voro , which are then averaged over ten independent simulation runs. This information is used to determine the fraction of particles to be pinned, such that N i V i voro for the free particles in the pinned system matches that of the porous gel network at the same density.

Lengthscales in the complex environments -The interplay between the lengthscale of the two environments and the persistent motion of the active particles will have a large impact on the dynamics. Therefore, to characterize the lengthscale the obstacles impose on the active particles we will use the pore chord length. The chord length is a measure of the distance between two interfaces in a homogeneous phase of a heterogeneous system. A chord is defined as the distance between two interfaces in a heterogeneous system, where the chord lies wholly within one phase. The chord length distribution p( ), defines the probability of finding a chord of length between and + d within one phase. We characterize the environments in this work by the mean pore chord length L c . In practice, this is determined by measuring chords of varying lengths along each axis of the three dimensional sample that lie wholly within the pore phase 61,62 . The mean chord length was measured and averaged for six independent configurations for each environment species. The gel networks have a mean pore chord length L c = 6.62σ, whereas for the random pinning system L c = 3.24σ, almost half that of the gel system. The difference in this measurement derives from arrangements of the particles comprising these structures: in the case of the gel, particles are arranged locally into dense branches and therefore the branches provide the relevant lengthscale in this system. Conversely, in the random system the particles are arranged in a non-overlapping random configuration, and the surrounding free space is then dependent on the shorter lengthscale of the average particle separation.

C. Dynamical analysis

The addition of obstacles into dense fluids greatly influences the dynamics, and in some cases a system may become arrested. The structural relaxation time τ α provides a useful metric which one can use to understand the variation of timescales across different state points and environments.

The relaxation time τ α is determined via the self part of intermediate scattering function:

F s (k, t) = 1 N N j=1 exp i k • ( r j (t) -r j (0)) ( 6
)
where k is the wavevector k = | k|, taken as 2π. We define

τ α as F s (k = 2π, τ α ) = e -1 .
Here the index j runs over all the particles. The persistent motion of active particles induces clustering and aggregation at boundaries. This will likely cause density fluctuations where some fraction of particles are in dense and crowded regions while others are in locally dilute regions. To quantify the degree to which these variations are taking place within different environments, we use the four-point dynamic susceptibility χ 4 . To calculate χ 4 , we follow the methodology of Lačević et al. 63 . For this, one must first define an overlap function w (|r j (0) -r i (t)|), where i, and j are particle indices. This measures the degree of spatial similarity between configurations of a system as a function of time. The overlap is unity inside a region |r j (0) -r i (t)| ≤ a and 0 otherwise, where a = 0.3σ. The fraction of overlapping regions in a system of particles at times 0 and t is given by:

Q(t) = 1 N N j=1 N i=1 w (|r j (0) -r i (t)|) . (7)
The fluctuation of Q(t) then defines χ 4 . This quantity is a susceptibility and measures dynamic heterogeneity [START_REF] Berthier | Dynamical heterogeneities in glasses, colloids, and granular media[END_REF] .

4 (t) = V N 2 k B T Q 2 (t) -Q(t) 2 . ( 8
)

III. RESULTS

Our results section is organized as follows. We first characterize the dynamical behaviour at low activity characterized by the Péclet number. We then move on to moderate activity, where significant structural changes are observed, related to motility induced phase separation. These we characterize with a number of measures, one-body measures such as the Voronoi volume associated with each particles, and many-body properties, accessed with the topological cluster classification 68 . Finally, we probe the case of very high activity, where we observe a re-entrant mixing.

A. Low Pe: crystal suppression and heterogeneous dynamics

Structural relaxation -As mentioned in the introduction, active systems can undergo a localization transition in the presence of quenched disorder. This phenomenon is often only present in the homogeneous phase i.e. for activity below that associated with the boundary of motility induced phase separation. To assess the extent to which complex environments impose localization on active Brownian particles, we measure the structural relaxation time (τ α ). This enables us to to determine the regimes in which these systems become arrested. 2 shows the variation of the structural-relaxation time in the bulk, gel network and random pinning systems as a function of Pe and total number density ρ.

In the bulk system (Fig. 2-bulk), the particles will crystallise for state points that fall below the freezing line. This crystal regime was studied previously from which the data for this phase boundary originates [START_REF] Moore | [END_REF] . Outside of the crystalline regime, the active particles exhibit a monotonic decrease in τ α as Pe rises for a given density, which is qualitatively similar to a temperature increase in passive systems [START_REF] Malins | [END_REF] .

With the addition of a complex environment, there is an emergence of slow dynamics distinct to that of the crystal in the bulk system. Figure 2-gel and Fig. 2random show the variation in τ α for the gel network and the random pins respectively. Comparing the two panels, we can see that the two systems operate on significantly different timescales, with the system with random obstacles exhibiting slower dynamics than the gel network for the majority of state points considered. Both the gel net-work and the random pinning systems exhibit structural relaxation times in excess of 10 5 τ R at high densities, plotted as the darkest blue contour. These are state points which will not relax within the maximum run time of the simulations and are thus for the purposes of this work, we define these as arrested. In this arrested regime, active particles are localized due to the constraints imposed by the gel or the pins, even at Péclet numbers of the order of 10 for dense ensembles in the random pins.

Dynamic Susceptibility -The interplay of the complex environment and the self propulsion will cause a degree of dynamical correlation in these systems, as clusters of active particles become aligned or absorbed at the boundaries. The dynamic susceptibility χ 4 provides a measure of this, it is plotted for the three systems in Fig. 3. The dynamic susceptibility manifests a peak at the timescale for which the particle dynamics are maximally correlated on the chosen lengthscale a (Eq. 8). This measure shares some similarities with the structural relaxation time. For example, in the passive systems (Pe = 0), the positions of the peaks indicate the same relationship as measured in τ α in Fig. 2, with the gel system being the fastest, followed by the bulk and then the random system coming significantly later. At low activity (Pe = 10), these peaks become more clearly defined and move to shorter times. Furthermore the relative positions of these peaks switch, with the dynamic correlations in the bulk system happening on shorter timescale than the complex environment systems. Finally, at higher Pe (the regime of motility-induced phase separation), the dynamic correlations are appreciably more prominent in the bulk system. Conversely, it is clear that for the gel and random system, χ 4 does not fully relax on the timescales we probe, as the environment enforces some degree of dynamical correlation over long timescales.

B. Intermediate Pe: MIPS in random environments

From the results presented in the previous section it is clear that the different confinement types cause significant and varying perturbations to the dynamics of active systems at low Pe. What remains unclear, is the mechanism by which active particles overcome localization due to the environment, and what behavior these particles exhibit when driven towards higher Péclet numbers, particularly in the regime of motility induced phase separation (MIPS). Compared to the lengthscale of MIPS, we know that the complex environment restricts the active particles to motion over shorter lengthscales; from the chord length measurements, the confining lengthscale is 6.62σ in the gel and 3.24σ in the case of random pinning.

We use the individual particle Voronoi volumes (V voro ) as a measure of the local density to provide insight into the relationship between self-propulsion and complex environments. Figure 4 displays representative snapshots of the bulk, gel network and random pinning systems in the steady-state at an activity of (Pe = 100), along- side a slice through each system. In the bulk system (top panel) the particles have undergone motility-induced phase separation. Within the large dense region are particles (coloured dark blue, which have V voro ≤ 0.8 in Fig. 4), clustering due to their persistent motion, which are surrounded by a gas of active particles at a lower density.

Figure 4 (middle panel) displays the behavior of active particles in the gel at high Pe. The structure of the gel is non-trivial featuring extreme variations in surface curvature and channel width. It is clear that this has a strong impact on the dynamics, with a distinct pattern emerging which is clearly dependent on the gel structure. In particular MIPS leads to a highly complex structure, with some pores being filled by the slow/dense phase and others by the fast/less dense phase. The location of the regions which acquire a low or high density under MIPS is fixed by the initial configuration, with independent runs always exhibiting the same demixing pattern.

The influence of the random pinning (bottom panel) on the active particle dynamics at high Pe is distinct from that of the bulk and gel systems. Here, the particles phase separate into a large cylindrical droplet somewhat akin to the bulk system. However, this droplet has random pins throughout its structure holding the dense phase in space. Furthermore, the presence of the random pins brings about an unexpected result: the arrangement of the random pins pre-defines the location in which the dense MIPS phase will form. For the same configuration of pinned particles, but distinct alternate initial positions and velocities; MIPS occurs in roughly the same region of the system, with only small fluctuations from run to run. Therefore, we believe that the location of the MIPS is somehow encoded in the initial positions of the randomly pinned particles.

! = 0$ ! ! = 1$ ! ! = 5$ ! ! = 10$ ! ! = 50$ ! ! = 100$ ! ! = 0$ ! ! = 1$ ! ! = 5$ ! ! = 10$ ! ! =
The three-dimensional snapshots in Fig. 4 provide a good indication of the phase behavior of these systems. However, they do not contain any information on the stability of the dense phase. Therefore, in Fig. 5 we plot the average Voronoi volume V voro as a function of space in both the gel network and random pinning systems. These plots give information regarding the way in which each system phase separates. For the gel, the inclusion of the network into the space does not allow for the formation of a single dense droplet as seen in bulk systems experiencing MIPS. Instead, the structure of the gel determines the locations in which the system will phase separate. Given the disordered nature of the gel, there are sections where the pores are more constricted, have tighter curvature or are less connected; these are the locations that will trap active particles. The spatial distribution of V voro in the random pinning system tells an alternate story. Active particles in this system at sufficient Pe will form a large droplet around a subset of pinned particles. The random pins in this droplet keep it stable over long time periods in a steady state where particles are exchanged between the droplet and the surrounding active gas.

The Voronoi volumes give information into the phase separation of these systems but not into the transport dynamics of the individual particles in this environment. Some insight into this aspect can be gained by looking at the average single particle displacements ∆r for the same system. The distributions of ∆r are plotted in Fig. 5 for displacements over the time period t = 6τ R , which is a time period of the order of τ α in the passive bulk system at this density. For the gel system, these displacements are largely uniform in the wider channels. Close to the surfaces of the gel particles the displacements are significantly less, indicating shorter movements along the surface. Furthermore, there are several locations where particles are localised with displacements less than the particle diameter. These are all located in regions of high surface curvature.

For displacements in the random system, there are four identifiable populations of particles, each with its distinct environment. The first is the group of localized particles. These are particles that have become trapped between the random pins and other particles in the dense phase and are recognizable as the dark spots dispersed through the dense phase. These are surrounded by the second group of particles that are not localized but remain trapped within the dense phase and are moving very slowly (∆r < σ). Beyond this are particles in the interface, which undertake mid-range displacements. Finally, outside the dense phase is the active gas where particles are completing relatively large and uniform displacements.

Time-evolution of MIPS in random pinning -Our protocol enables us to investigate the process by which the system undergoes MIPS. To this end, we show a time-sequence of snapshots of a pinned system undergoing MIPS. In Fig. 6 we show the formation of MIPS, starting from the passive WCA system as described in section II. Here we consider the same state point as in the previous figures (ρ = 0.87, Pe = 100) and plot the Voronoi volumes as in Fig. 5 (top row). The timeevolution of this data is also shown in the supplementary movie.

At time t = 0, we see that the system is largely uniform in density. Visual inspection indicates the following. Even at quite small times (t τ R ), there are fluctuations in density which are small, both in the change of local volume per particle and also in their spatial extent. Over time, the larger fluctuations (in terms of their spatial extent) appear to grow at the expense of smaller fluctuations and even at the quite short time of 10τ R , the pattern seems to be largely fixed. Longer times correspond to an increase in density difference between the particle-rich regions (blue) and the more dilute regions. During this time regime, the interface between these regions becomes better defined. While we have been able to observe the time-evolution of MIPS under random pinning, exactly what causes the spatial distribution of the MIPS phases and how this is encoded in the random pins seems to be a challenging problem for the future.

Density fluctuations as a function of Pe -So far we have seen examples of how these systems phase separate at Pe ≈ 100, now we will look at how the local density fluctuates in these systems for different Pe. The probability density of the Voronoi volumes for each system are shown in Fig. 7 for various Pe. In the absence of activity ie. Pe = 0, all systems feature an approximate Gaussian distribution. Common to all the environments there is the fact that the addition of activity produces a shift in the peak towards smaller volumes and a broadening of the tail of the distribution towards larger volumes. In the bulk system (Fig. 7-bulk) the distributions feature a non-monotonic trend in the spread. First increasing as a function of Pe up to a maximum at Pe = 60, before decreasing. This is the first sign of re-entrant MIPS mixing, which will be the focus of Section III D.

For the gel network (Fig. 7-gel) and the random pins (Fig. 7-random), this broadening is monotonic, with the gel network covering a wider range of volumes. However, for the random pins and the bulk systems we observe the emergence of twin-peaked distributions at higher activity as a result of phase separation. For both the gel and the random pins, the influence of the complex environment leads to a splitting of the active particles into more than one population, with a proportion of active particles aggregating or becoming localized because of interactions with the environment.

The persistent motion of active particles induces symmetry breaking that causes them to aggregate at surfaces and walls. In these systems the surfaces could be the edge of a MIPS dense phase, the surface of the gel network, or a dense cluster of pins. These surfaces collect active particles. To determine how the environment structure affects the collection of particles at surfaces, we count the number of particles located in locally dense regions N D , defined for particles where V voro < 0.8. In Fig. 8 we plot the fraction of localised particles N D /N as a function of Pe.

Active transport in complex environments -We have seen so far that with a progressive increase in Pe, all three systems undergo dramatic changes in terms of local density. We have also seen that at high Pe the average single particle displacements reveal information concerning the interaction the active particles have with their environments. Like with the Voronoi volumes, we plot the probability densities of active particles in the three systems for various Pe [Fig 9(a)]. At first glance, it is clear that for all systems the particles complete larger displacements as Pe is increased. Moreover, these distributions are highly featured and reveal a great deal of information to the behaviour and interactions of the active particles.

In the bulk system [Fig. 9(a-bulk)], the ∆r distribution is a single peaked at Pe=0. As Pe increases, this distribution shifts to higher displacements and we observe the growth of a second peak; indicating the presence of MIPS in the system. With a further increase of Pe there is a continuous transition between the relative heights of the peaks as the fraction of particles in the dense phase grows. This is corroborated by the growth of regions of locally dense particles over this range of Pe observed previously in Fig. 8.

The displacements for the gel network and the random pinned system are plotted in Fig. 9(a-gel) and Fig. 9(arandom) respectively The distributions of both of these systems show the splitting of a single population into two or more populations with the progressive increase of Pe. The first of these is the emergent peak at very small displacements ∆r/σ ∼ 10 -1 . These particles move only a small fraction of their diameter and have become localized as a result of the interplay between their activity and the environment.

Interestingly, both systems in complex environments [Fig. 9(a-gel/random)] feature a strong peak at ∆r/σ = 1, with some smaller features at subsequent integer displacements. Examination of the average displacements in Fig. 5, shows that they are located along the surfaces of the gel network and in small pockets of lower density in the dense phase of the pinning system. The location of these displacements makes it clear that they are a result of particle re-arrangements at the obstacle interface, and in dense particle clusters. For the gel network [Fig. 9(a-gel)] at P e > 0, the remaining particles are in a single large population, moving comparable distances to particles in the bulk system. However, in the case of the random pinned system [Fig. 9(a-random)], at Pe ≥ 40 there is a splitting of larger displacements across two length scales, one at the interface and the second in the active gas.

Thus far we have primarily considered two observables V voro and ∆r, both of which tell part of the story. These two observables can be correlated to complete this picture, the result of this is plotted in Fig. 9b. The result is a series of contours stacked logarithmically, each layer denotes the level at which a percentage of the data lies below. These plots provide some insight into the population splitting phenomena we have seen so far. We see that the combination of confinement and activity create a subpopulation of particles that are arrested and have very little free space, a feature not found in the bulk system. Interestingly, this arrested group is relatively larger FIG. 10. Higher-order structural analysis using the topological cluster classification. Here particles in specific local environments corresponding to clusters of 5-13 particles, and the hexagonal close-packed and face centered cubic crystals are considered. The number of particles in each environment Nc is then plotted as a function of Pe for the three geometries under consideration, bulk, gel and random systems at ρ = 0.87. Colors of the lines and data points correspond to the clusters depicted in the legend. The colors of the particles in the renderings in the key correspond to geometric properties of the clusters. In particular, the grey particles are in 3,4 or 5-membered rings and the yellow particles correspond to so-called "spindles" 68 . in the random pinning system. For the systems at high Pe in Fig. 9(b), the particles that have V voro < 0.8 and in some form belong to dense clusters, these particles still cover a wide range of displacements. This is likely a combination of two phenomena, first is that of the particles that have been localized over a longer time frame. These are particles that are not moving and have very little free space. The other case, is that of particles which have been mobile but have very little space. These will be particles that have moved from a position at a previous time and are now located in a dense cluster or at the object interface.

C. Local structure

So far, we have focussed mainly on one-body structural properties via the Voronoi volumes. It is instructive when studying amorphous systems to consider higherorder structural correlations, as a precise means to probe smaller changes in the local structure. To this end, we use the topological cluster classification (TCC) 68 . The TCC identifies local environments whose bond topology is identical to that of small minimum energy clusters of a suitable reference system, here the Lennard-Jones model. Some of us have shown that this is appropriate for passive WCA particles 66 and also have investigated the effect of activity on a similar (bulk) system [START_REF] Moore | [END_REF]67 . Here, we consider clusters of differing sizes (5-13 particles) as indicated in Fig. 10. The number of particles in each of these clusters N c , scaled by N is then plotted as a function of Pe for ρ = 0.87. We further consider the influence of the local environment where N c is the number of particles participating in a particular cluster. To identify the cluster population we need the bond network, which we identify with a modified Voronoi decomposition (specifically we set the Voronoi parameter of Ref. 68 to f c = 0.82). Now we are interested in the active particles, but in order to identify the clusters, we include the immoblized particles in the bond network. We then run the TCC, but when analyzing data, we only consider active particles.

Figure 10 shows a common trend between the bulk, gel and random environments: the equilibrium local structure population is disrupted by increasing Pe. Both of the confined environments show cluster populations that are lower compared to the bulk case. The random environment, where the position of the obstacles is taken from equilibrium configurations, at Pe = 0 has the same population numbers as the bulk case, but as soon as activity is switched on all cluster populations fall more abruptly compared to the bulk case. The first rapid decrease in local structure at Pe 30 coincides with the formation of MIPS. Unlike the gel and random pinning systems the bulk case experiences another large drop in local structure population for Pe 140. As we will see in the next section, this drop corresponds to a re-entrant mixing. Local structures are present in lower population in active systems than in similar equilibrium systems as has been seen previously 67 . However, it is possible that phase separation or the influence of the environment could influence this.

D. Motility-Induced Mixing

In addition to the phase behavior discussed thus far, for our purposes, there is one more regime to be considered. At very high Pe, ABPs will transition from a demixed state due to motlity-induced phase separation to a mixed state, ie a homogeneous active fluid. A similar behavior has been observed previously in Ref. 43,69 in bulk systems. We confirm these observations in the bulk case. The presence of the transition is apparent in the probability density of V voro in the bulk system [Fig. 11(a)]. These distributions show two populations at intermediate Pe, but a single population at Pe = 0 and at Pe > 120. Looking at the correlation of V voro with ∆r confirms the single re-entrant phase at high Pe [Fig. 11(b)].

To study the re-entrant MIPS behaviour in systems with quenched disorder in Fig. 11c we plot the variance of the probability density of V voro as a function of Pe: non monotonic behaviour in this quantity signals a re-entrant phase. Looking at this data we can see that for the bulk and random systems, as Pe increases, the variance also increases up to a maximum; beyond which it decays to an intermediate value. The peak corresponds to the state where there is a roughly equal fraction of particles in the dilute and dense phases N D /N ≈ 1/2. We plot the same for the systems with confinement in Fig. 11(c). Notably, the random pinning system also undergoes a transition to a re-entrant fluid. However this transition is delayed relative to the bulk, indicating that the presence of the pins works to stabilize MIPS at high Pe. For the gel system we do not observe the re-entrant behavior within the considered range of Pe.

IV. CONCLUSION

Suspensions of Active Brownian Particles show a rich range of dynamical behavior, and our goal was to explore the influence of complex confinement at high densities.

To do this, we prepared confining geometries with different static properties: randomly pinned particles from an equilibrium bulk configuration and from a porous gel structure. Both confining geometries are constructed to have the same free volume available to the mobile particles, thus allowing a direct comparison of the effects of the static lengthscale of the confinement on the behavior of active particles.

We first explored the phase behavior at low Pe revealing how pinning suppressed the crystallization of the fluid at high densities. The relaxation time of the particles is slowed down by the obstacles (more for the random case compared to the gel case) and also it becomes more dynamically heterogeneous, as confirmed by the study of the overlap function and four-point susceptibility.

At intermediate Pe the bulk system displays MIPS, where the systems forms domains of dense/slow regions and low density/fast regions that nucleate and grow in a similar manner to equilibrium phase separation of two disordered phases. Surprisingly, we find that the obstacles not only do not suppress MIPS formation, but actually act to stabilize it. Random obstacles display a very similar phase-separation pattern as the bulk case, but the domains do not appear homogeneously in the system, but are always formed from the same regions of the sample. The mechanism of MIPS nucleation from random obstacles is still unknown, and should be addressed in future studies. In particular, the effect of system size and changing the state point would be most important to explore. The system size that we have studied here fully demixes to two distinct phases, but the timescale for this for different system sizes, and indeed whether larger systems fully demix would benefit from a detailed finite size analysis. Furthermore, exactly how the pinned particle environment encodes the spatial distribution of the MIPS patterns remains an outstanding challenge. In the case of the gel, the MIPS domains change completely, and forms a complex structure where the active and inactive regions occupy different pores of the structure.

Finally we considered how local structure is perturbed by the activity, and revealed that re-entrant MIPS behavior is suppressed (or moved to higher Pe) by the random environments. In particular the gel environment seems to be the most effective in stabilizing the density and activity fluctuations. In this case, we attributed this behavior to the absorption at the rough walls, where the persistent motion of active particles creates localized and highly dense regions, which in turn frees space inside the pores that increase particle transport in the system.

We believe that these results could be important for better understanding transport in biological environments and for guiding future studies of active matter particles in random media. The confining environments that we consider are experimentally realizable 49 , and indeed combining this with some means to break the symmetry in the future might enable investigations of topotaxis 41 . Other possibilities include the possibility to explore the interplay of the complex environments such as those we have investigated and state functions such as pressure [46][47][48] which also might be measured experimentally in colloidal systems 70 .

SUPPLEMENTARY MATERIAL

Supplementary movie 1: This movie corresponds to the stills in Fig. 6. The state point is ρ = 0.87 and Pe = 100 for the random pinning system. The movie shows the emergence of dense and dilute regions.

ACKNOWLEDGMENTS

FJM was supported by a studentship provided by the Bristol Centre for Functional Nanomaterials (EP-SRC grant EP/L016648/1). JR acknowledges support from the European Research Council Grant DLV-759187. CPR acknowledges support from the European Research Council (ERC Consolidator Grant NANOPRS, project number 617266).

APPENDIX Overlap

The overlap function Q(t) compares a particle configuration with itself at a later time. An important feature of this measure of similarity is that it is not particle specific, it matters not whether it is the same particle occupying that space, only that there is a particle there. With this in mind, it is clear why the bulk system [Fig. 12(a)] behaves as it does. In this system, there are no fixed obstacles, and thus there are no points at which a cluster of particles could be anchored. Therefore, even in a system with MIPS present the centre of mass of a dense cluster remains diffusive and thus Q(t) decays to a small value. Q(t) does not decay to 0 in this case in the bulk due to the relatively high density, as there will always be some degree of overlap.

For active particles in the other two systems, there is clear evidence of the interplay between the selfpropulsion and the complex environment. Both systems converge to an overlap value, higher than that of the bulk system. This is a product of the fixed geometry of the obstacles. Although the obstacle particles are not considered, there will be regions in structure that will be more likely to trap active particles, and even though particles are mobile, there is a high probability that there will be some particles in these regions. Interestingly, the gel and the random pins approach their convergent overlaps from different directions as Pe increases. For the gel system [Fig. 12(b)], the overlap value Q(t) increases with Pe at longer times, while decreasing with an increased rate at shorter times. The increased rate on short timescales is explained by the increased propulsion velocity promoting a quicker re-configuration of the active particle populations. However, this increase in Pe has the added effect of increasing the likelihood of a particle becoming trapped against the walls of the gel network, this is supported by the increase in the localised fraction N loc /N with Pe in Fig. 8.

Like in the bulk and the gel network, the overlap of the random system also decays at an increasing rate as Pe increases at short times (Fig. 12 c). In the absence of activity, the random system particles have a very slow decay in Q(t). This is a result of the pinned particles dramatically slowing down the dynamics of passive particles, this also seen in the alpha-relaxation time in Fig. 2, where τ α was measured to be of the order 10 4 larger than that of the bulk or gel systems. Unlike in the gel system, Q(t) for the random pinning approaches its convergent value from above. This is likely due to the arrangement of the obstacles in the random pinning system, since the pinned particles are dispersed through the entire space, the active particles have a higher chance of becoming trapped. The random system converges to have higher overlap value than the gel, this is due to the larger fraction of particles becoming localized in the random system compared to the gel at the same activity and density.

FIG. 1 .

 1 FIG. 1. Preparation of complex environments for ABPs. (a) Percolating gel network at number density ρ = 0.38. (b) Crosssection through the gel network of depth 2σ. (c) A collection of randomly pinned obstacles at ρ = 0.31. (d) Cross-section through the random obstacles of depth 2σ. (e-f) Gel network (grey) filled with active particles (blue) to a total density ρ = 0.42. (g-h) Random obstacles filled with active particles to a total number density ρ = 0.42.

FIG. 2 .

 2 FIG. 2. Structural relaxation time (τα) as a function of number density ρ and Péclet number Pe. Here we consider the bulk, gel network and random pinning systems from left to right. Colours indicate the common logarithm of τα. Crystalline states are marked in white. The number density ρ = 0.87 is emphasized with a black dotted line. Due to the large amount of sampling required, this data was measured using the smaller simulation box size.

  χ

Figure

FIG. 4 .

 4 FIG. 4. (Left) 3d Snapshots of the three systems at Pe = 100, particles are coloured by their Voronoi volumes (Vvoro), obstacles are not rendered. (Right) A slice through each system (depth 4σ).

FIG. 5 .

 5 FIG. 5. Heat map of Voronoi volumes and single particle displacements in the gel network (left), and random pinning (right) systems respectively. Note that displacements are measured over ∆t = 6τR. All systems are at ρ = 0.87, and Pe = 100.

  FIG. 6. Time-evolution of the pinned system undergoing MIPS. The state point is that of Figs. 5 (right), ρ = 0.87, and Pe = 100. The color map represents the Voronoi volumes of the active particles and the pale yellow particles are the pinned particles. Each snapshot is taken at the time noted underneath from the start of active motion. Data are taken from a slice of depth σ.

FIG. 7 .FIG. 8 .

 78 FIG. 7. Probability density of the Voronoi volumes (Vvoro), for the bulk, gel, and random pinning systems respectively. All systems are at ρ = 0.87, and various Pe (see figure legend).

FIG. 9 .

 9 FIG. 9. (a) Probability density of the Voronoi volumes (Vvoro), for the bulk, porous network, and random pinning systems respectively. All systems are at ρ = 0.87, and various Pe (see figure legend). (b) Correlation of the Voronoi volumes and the single particle displacements of motile particles in the gel network, random pins and the bulk system. All systems are at a density of ρ = 0.87, and plotted for Pe = 0, 20, and 120. Black lines contain 90% of data. Colours show contour levels below which the indicated percentage of the data will lie beneath.

FIG. 11 .

 11 FIG. 11. (a) Probability density of the Voronoi volumes in the bulk system at various Pe. As Pe increases the distributions show first one, then two populations, and then a re-entrant single phase at high Pe. (b) Correlation of Vvoro and ∆r in the bulk system at Pe = 160 shows a single population. (c) Variance of the distribution of Voronoi volumes in the porous network (σ 2 voro ), random pinning and the bulk systems as a function of Pe, all at density ρ = 0.87.

FIG. 12 .

 12 FIG.12. The overlap Q(t) Overlap is shown for for particles in the bulk, gel, and random systems at ρ = 0.87 for various Pe (see legend).