Active Brownian particles in random and porous environments

Active matter concerns systems comprised of individual bodies undergoing motion via self-propulsion.¹ This description encompasses a plethora of systems over a wide range of lengthscales from bacteria,^2–4 biological microswimmers,⁵ schools of fish,⁶ bird flocks,⁷ to human crowds.⁸ 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 that, nevertheless, capture the key elements of the behavior of more complex systems is active colloids.^9–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 such as porous soils¹³ and organic tissues.¹⁴ 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;⁹ 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 toward 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 down the dynamics of such systems, providing access to rare states.^19–22 A special kind of localization has been explored in glass-forming systems^23,24 and may be realized using colloidal systems.^25–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 experiments or computer simulations.²⁸ 

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 their persistent motion, even in the absence of hydrodynamic interactions.^29,30 For many-body systems in two dimensions, the influence of disordered landscapes on the dynamics of active systems has been shown to manifest in clogging and localization transitions,^31–33 subdiffusion over long timescales,^34–36 destruction of flocking clusters,³⁷ suppression of Motility-Induced Phase Separation (MIPS), and prevention of uniform wetting at boundaries.³⁸ Furthermore, the manipulation of complex environments has been shown to provide a degree of control over the transport of active matter in the form of sorting³⁹ and over 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 behaviors, such as MIPS^42,43 and the formation of active crystals,^44,45 and fundamental properties—for example, pressure and the equation of state differ drastically from what might be expected from equilibrium systems.^46–48 

However, the question of how active Brownian spheres couple to a complex surrounding environment remains unanswered. Recently, there has been interest in experiments with mesoscale active matter in 3D. In one study, a random heterogeneous environment was found to impose strong inhibitions on the active transport of bacteria,¹⁵ 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.⁴⁹ This latter example is a well-controlled 3D colloidal model system that 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.⁵¹ Insights into the transport of active matter in 3D complex environments could provide a major step toward the control of such systems and aid in the progress toward 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 the 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 Sec. II, we outline the computational methods used to study these systems; in Sec. III, we discuss the results of the simulations; and finally, in Sec. IV, we will summarize and conclude our findings and discuss the future work in this area.

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

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 = 144 000 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 18 000 to 24 000. Analysis at a constant density is always conducted at ρ = 0.87. This state point is chosen such that it lies below the freezing line in the bulk.

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 and these bulk simulations use the approach outlined in Sec. II A.

In Sec. II B 1, 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 [Figs. 1(a) and 1(b)] is a heterogeneous system comprised of two distinct meso-phases that percolate through the entire simulation box comprised of a particle-rich phase and a particle-poor phase in which, for our parameters, no particles are found. The random environment [Figs. 1(c) and 1(d)] is comprised of randomly pinned particles that create a number of discrete obstacles dispersed throughout the system.

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], with the particles interacting via the Morse potential Eq. (4). This system is then evolved for 5 × 10⁷ 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)²/6D_T. An example of the resulting gel is shown in Figs. 1(a) and 1(b).

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 is 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 [Figs. 1(c) and 1(d)].

The fixing of the accessible free volume of the mobile particles enables the 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 $Vvoroi$⁠. The sum of these volumes provides the total volume accessible to the free particles $∑iNVvoroi$⁠, which are then averaged over ten independent simulation runs. This information is used to determine the fraction of particles to be pinned such that $∑iNVvoroi$ for the free particles in the pinned system matches that of the porous gel network at the same density.

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 the 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.

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 that one can use to understand the variation of timescales across different state points and environments.

This section is organized as follows: we first characterize the dynamical behavior 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. We characterize these with a number of measures—one-body measures, such as the Voronoi volume associated with each particle, and many-body properties, accessed with the topological cluster classification.⁶⁸ Finally, we probe the case of very high activity, where we observe a re-entrant mixing.

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 determine the regimes in which these systems become arrested. Figure 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 crystallize for state points that fall below the freezing line. This crystal regime was previously studied, from which the data for this phase boundary originate.⁴⁵ 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.⁶⁵ 

With the addition of a complex environment, there is an emergence of slow dynamics distinct from that of the crystal in the bulk system. Figure 2 (gel and random) shows 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 network and the random pinning systems exhibit structural relaxation times in excess of 10⁵τ_R at high densities, plotted as the darkest blue contour. These are state points that will not relax within the maximum run-time of the simulations and, 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.

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 χ₄ provides a measure of this, and 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 a shorter timescale than the complex environment systems. Finally, at a 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, χ₄ does not fully relax on the timescales we probe, as the environment enforces some degree of dynamical correlation over long timescales.

From the results presented in Sec. III A, 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 toward 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 insights into the relationship between self-propulsion and complex environments. Figure 4 displays the representative snapshots of the bulk, gel network, and random pinning systems in the steady-state at an activity of Pe = 100, alongside 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 (colored 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 a 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.

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 a 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. The active particles in this system at a 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 localized 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.

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 Sec. 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 time evolution of these data is also shown in the supplementary material.

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 a 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. In 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.

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 values of Pe. The probability density of the Voronoi volumes for each system is shown in Fig. 7 for various values of Pe. In the absence of activity, i.e., 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 toward smaller volumes and a broadening of the tail of the distribution toward 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 Sec. 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 localized particles N_D/N as a function of Pe.

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. Similar to the Voronoi volumes, we plot the probability densities of active particles in the three systems for various values of 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 about the behavior and interactions of the active particles.

In the bulk system [Fig. 9(a), bulk], the Δr distribution is a single peak 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 in 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 random), respectively. The distributions of both of these systems show the splitting of a single population into two or more populations with the progressive increase in Pe. The first of these is the emergent peak at very small displacements Δr/σ ∼ 10⁻¹. 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. An 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 Pe > 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 lengthscales, one at the interface and the other 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, and the result of this is plotted in Fig. 9(b). 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 creates 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 in the random pinning system. For the systems at high Pe values in Fig. 9(b), the particles that have V_voro < 0.8 and that 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 that 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.

So far, we have focused mainly on one-body structural properties via the Voronoi volumes. When studying amorphous systems, it is instructive to consider higher-order structural correlations as a precise means to probe smaller changes in the local structure. To this end, we use the topological cluster classification (TCC).⁶⁸ The TCC identifies local environments whose bond topology is identical to that of small minimum energy clusters of a suitable reference system, here is Lennard-Jones model. Some of us have shown that this is appropriate for passive WCA particles⁶⁶ and also have investigated the effect of activity on a similar (bulk) system.^45,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 immobilized particles in the bond network. Then, we 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 the activity is switched on, all cluster populations fall more abruptly compared to the bulk case. The first rapid decrease in the 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 the local structure population for Pe ≃ 140. As we will see in Sec. III D, this drop corresponds to a re-entrant mixing. Local structures are present in the lower population in active systems than in similar equilibrium systems as has been seen previously.⁶⁷ However, it is possible that phase separation or the influence of the environment could influence this.

In addition to the phase behavior discussed thus far, for our purposes, there is one more regime to be considered. At very high Pe values, ABPs will transition from a demixed state due to motility-induced phase separation to a mixed state, i.e., a homogeneous active fluid. A similar behavior has been observed previously in Refs. 43 and 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 values 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 values [Fig. 11(b)].

To study the re-entrant MIPS behavior in systems with the quenched disorder in Fig. 11(c), we plot the variance of the probability density of V_voro as a function of Pe: non-monotonic behavior in this quantity signals a re-entrant phase. Looking at these 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 the 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 values. For the gel system, we do not observe the re-entrant behavior within the considered range of Pe.

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 values, 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 values, the bulk system displays MIPS, where the systems form domains of dense/slow regions and low density/fast regions that nucleate and grow in a similar manner to the 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 form a complex structure where the active and inactive regions occupy different pores of the structure.

Finally, we considered how the local structure is perturbed by the activity and revealed that the re-entrant MIPS behavior is suppressed (or moved to higher Pe values) 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 the particle transport in the system.

We believe that these results could be important for better understanding the 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,⁴⁹ and indeed, combining this with some means to break the symmetry in the future might enable investigations of topotaxis.⁴¹ 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–48 which also might be measured experimentally in colloidal systems.⁷⁰ 

See the supplementary material, 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.

F.J.M. was supported by a studentship provided by the Bristol Centre for Functional Nanomaterials (EPSRC Grant No. EP/L016648/1). J.R. acknowledges the support from the European Research Council Grant No. DLV-759187. C.P.R. acknowledges the support from the European Research Council (ERC Consolidator Grant NANOPRS, Project No. 617266). C.P.R. and T.B.L. acknowledge support from Engineering and Physical Sciences Research Council, Award No. EP/T031077/1.

The authors have no conflicts to disclose.

Fergus Moore: Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Project administration (equal); Writing – original draft (equal); Writing – review & editing (equal). John Russo: Conceptualization (lead); Formal analysis (equal); Methodology (equal); Project administration (equal); Supervision (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Tanniemola B. Liverpool: Formal analysis (equal); Investigation (equal); Project administration (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). C. Patrick Royall: Conceptualization (supporting); Data curation (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

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 center 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 self-propulsion 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 the 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 and this is supported by the increase in the localized fraction N_loc/N with Pe in Fig. 8.

As 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, and this is also seen in the alpha-relaxation time in Fig. 2, where τ_α was measured to be of the order 10⁴ 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 and the active particles have a higher chance of becoming trapped. The random system converges to have a higher overlap value than the gel, and 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.