Feedback-Induced Phase Transitions in Active Heterogeneous Conductors

An active conducting medium is one where the resistance (conductance) of the medium is modified by the current (flow) and in turn modifies the flow, so that the classical linear laws relating current and resistance, e

An active conducting medium is one where the resistance (conductance) of the medium is modified by the current (flow) and in turn modifies the flow, so that the classical linear laws relating current and resistance, e.g., Ohm's law or Darcy's law, are modified over time as the system itself evolves. We consider a minimal model for this feedback coupling in terms of two parameters that characterize the way in which addition or removal of matter follows a simple local (or nonlocal) feedback rule corresponding to either flow-seeking or flow-avoiding behavior. Using numerical simulations and a continuum mean field theory, we show that flow-avoiding feedback causes an initially uniform system to become strongly heterogeneous via a tunneling (channel-building) phase separation; flow-seeking feedback leads to an immuring (wallbuilding) phase separation. Our results provide a qualitative explanation for the patterning of active conducting media in natural systems, while suggesting ways to realize complex architectures using simple rules in engineered systems. DOI Flow through heterogeneous conductors is important in many problems in physics, biology, geology, and engineering. While most studies are limited to transport through a static medium, transport and flow can often modify the medium itself, which then modifies transport. Thus, these active heterogeneous conductors coevolve with the transport through them. For example, in geophysics, branching patterns are formed through the interplay of erosion, transport, and deposition [1,2], while networks of electrical fuses form complex patterns of failure due to the interplay in changes in current and resistance [3,4]. In biological systems, gradients and physical flows often arrange matter through feedback mechanisms to control transport at the cellular, organismal, and societal level. Specific examples of active heterogeneous conducting media in biology include network formation of slime molds [5][6][7][8][9], formation and remodeling of vascular networks [10][11][12][13][13][14][15], and the nest architectures of social insects [16][17][18][19][20][21].
Elements universal to these active conducting media are conservation of current, feedback that changes resistance, and stochasticity. A minimal distillation of the common elements in these different systems corresponds to a stochastically evolving network driven uniformly by fluxes and forces at the boundary due to pressure, voltage, or concentration gradients. Since problems involving steady state diffusion of heat, concentration gradients, or flow through a porous material are mathematically analogous to current flow through electrical circuits, we will use the language of circuit theory from now on.
We focus on a translationally symmetric case with periodic boundary conditions and a uniform driving voltage in the verticalẑ direction. The vertices are arranged in a square network, with the current I between neighboring vertices i; j given by where V i is the voltage at vertex i, g is the driving force,ẑ is the vertical direction, andr ij is the relative position of i; j. Each vertex i either contains a particle (ρ i ¼ 1) or is empty (ρ i ¼ 0), and the resistance between two vertices Ω ij increases by ΔΩ when full; Particles are removed from their vertices at a flow-dependent rate proportional to rðv i Þ, where v i ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi q is the current through the cell [22]. They are then added to an empty vertex with probability proportional to aðv j Þ, leading to a simple algorithm for the evolution of the medium [23]: 1) Remove a particle from filled vertex i randomly selected with probability proportional to rðv i Þ.
2) Solve for the new current through the network [24] [25]. 3) Add the particle to an empty vertex j randomly selected with probability proportional to aðv j Þ. 4) Solve for the new current through the network. We emphasize that these four steps conserve particle number, but the movement of particles is nonlocal in that the distance between i; j may be arbitrarily large [26]. Furthermore, we note that the these dynamics lack detailed balance [ Fig. 1(d)].
Since the removal and addition rates can either increase or decrease with local current in the active systems described earlier, we explore this range of possibilities with two parameters α R , α A . For the removal process, we choose rðv i Þ ¼ v i −α R ; at positive α R , particles in high current will rarely be removed (flow-seeking), and vice versa. For the addition process, we choose aðv j Þ ¼ v j α A ; at positive α A , empty vertices with high current will likely be filled (flow-seeking) and vice versa.
For simplicity, we have chosen our functional forms such that in any circuit, rðv i Þ ∝ g −α R , aðv j Þ ∝ g α A . Because only the ratios of currents are important, we may set this system of equations to be dimensionless by making the substitutions: Ω 0 → 1, ΔΩ → ΔΩ=Ω 0 , g → 1. This gives four dimensionless parameters: ΔΩ; α R ; α A , andρ, the mean particle density. We want the difference between filled and unfilled vertices to be large; here we choose ΔΩ ¼ 19 such that the resistance between filled vertices is 20 times higher than empty vertices.
Simulations.-For each α R ; α A , we start the system at a uniform densityρ and evolve it so that each particle or hole moves an average of one thousand times. This procedure is repeated twenty times for the parameters α R ; α A ¼ ð−6; −3; …; 3; 6Þ,ρ ¼ 1=4, computing the current through the entire network at every step. The results in Fig. 1 show that the system spontaneously forms channels (with high conductivity) at sufficiently negative α R ; α A , consistent with natural systems sharing these qualitative biases. In drainage networks, erosion increases with current (α R < 0), while deposition decreases (α A < 0), while in biology, ants have been observed to remove corpses from high-wind areas (α R < 0) and place them in low-wind areas (α A < 0); both of these systems form channels. Likewise, the system forms walls at sufficiently positive α R ; α A , consistent with fuse networks (α A > 0). This kind of phase transition is also similar to those seen in driven lattice gas models [27,28].
When the system channelizes due to negative α R , we find thin channels; negative α A gives thick channels. To understand these transitions, we consider their robustness to perturbations. When the system has formed a set of parallel channels, occasionally a channel gets blocked [ Fig. 1(c)]. When the channel is thin, current must go through the clog blocking the channel, and therefore total current through the channel is reduced while the clogging particle has much of this current forced through it. On the other hand, when a channel is thick, current will go around the clog, and so is barely impeded. Therefore, at negative α A , the thin channel has reduced current, and this clogging will cause the thin channel to fill, while a large channel is much more robust and will not be filled. On the other hand, for negative α R , the clog in the wide channel has little current through it, and lingers, allowing the wide channel to eventually be filled; the clog in the thin channel has high current forced through it and is quickly removed. A similar argument shows that positive α A leads to thin walls, while positive α R leads to thick walls. The system exhibits large hysteretic effects which are especially strong at very negative α A ; α R , when the system is "frozen" and fluctuations are suppressed.  vertices are covered by gray squares. Current is driven in the upwardsẑ direction; direction and magnitude of current between neighboring vertices is indicated by red arrows. (b) Phase diagram forρ ¼ 1=4. Each individual box represents a single system that has equilibrated for a particular (α A ; α R ), where α A ; α R have values (−6, 3, 0, 3, 6). At positive α R , the system forms thick walls; negative α R gives thin channels. Positive α A leads to a series of thin walls; negative α A gives thick channels. At positive α R , negative α A , a phase separation occurs at both orientations, giving a set of clumps.
(c) Robustness of thick and thin channels. (d) Lack of detailed balance: α A < 0, and K 2 ≈ K 3 , as the current through the lower right vertex has only weak dependance on the occupation of the upper left vertex. However K 1 ≉ K 4 , as the current through the upper left vertex strongly depends on the occupation of the lower right vertex; therefore, Interestingly, when α R is positive and α A is negative, both channelization and wall-building phase separations occur, as the system phase separates into thick clumps, which can be explained through a continuum model.
To characterize the wall-building phase separation, we use the row density ρ z ¼ P x ρ xz =L x as an order parameter, consistent with the fact that when the distribution of row densities becomes bimodal [ Fig. 2(d) example, wall building has occurred] [29]. We characterize the channelization phase separation in terms of the column density ρ x ¼ P z ρ xz =L z ; when this is bimodal, channelization has occurred (Fig. 1). We characterize a clumping phase separation through local densityρ r ¼ ð1=AÞ where Θ is the Heaviside function [30]; when this is bimodal and neither column or row density are bimodal, clumping has occurred.
Characterizing the distribution of filled vertices with a mean density field ρ, the continuum version of Eq. (1) is where κ, the conductivity, is a function of density. Similarly, in the continuum limit, the discrete addition and removal activity are replaced with a stochastic equation for density evolution, where Rðρ; J; α R Þ is the mean removal rate from a region of with density ρ, current J, and a bias of α R . A is the mean addition rate, and J is itself a functional of ρ obtained by solving Eq. (2). N ¼ ∬ R=∬ A acts as a sort of chemical potential, which is set to conserve total particle number. η is the stochastic noise term, whose form will be discussed later [31]. The predictions of the continuum model depend strongly on the functions κ; A; R. To determine them, we use a hybrid approach, sampling via numerical experiment using randomly placed particles and then varying the density to approximate the entire functions, leaving us with no fitting parameters [32]. While the discrete and continuum model both lack detailed balance and thus we cannot write down a free energy function associated with their dynamics, we can use equilibrium considerations when certain limits or symmetries are assumed. Nearly disordered limit.-In the limit where α R ; α A → 0, only the linear response is important. We write Eq. (3) as where T ðρ; J;αÞ ¼ −Rðρ; J; α R Þ þ N Aðρ; J; α A Þ is the time derivative functional andα ¼ ðα R ; α A Þ. In Fourier space, we separate the effects of J and ρ on T : When k ∝x, density fluctuations are horizontal (channels) and cause fluctuations in current; when k ∝ẑ, fluctuations are vertical (walls) and current is uniform. The full relation can be shown [33] to be ½∂J z ðkÞ=∂ρðkÞ ¼ ð∂κ=∂ρÞk 2 x =jkj 2 . To first order _ ρðkÞ ¼ −½dT ðkÞ=dρðkÞρðkÞ þ η, where η is independent of k. This allows us to predict the mean amplitude of fluctuations: We note that this is independent of the magnitude of k, and is a function of its direction alone [ Fig. 2(a)], because  of the dipolelike interactions between particles which inhibit current upstream and downstream, while increasing it laterally. Similar interactions have been noted in fuse networks [34][35][36].
x translation symmetry.-To predict wall building, we assume translation symmetry in thex direction, so that Jðx; zÞ is constant throughout the system. In the discrete model the current through any vertex is proportional to J, so that we may make the simplification Rðρ; J; α R Þ → Rðρ; α R ÞJ −α R ; Aðρ; J; α A Þ → Aðρ; α A ÞJ α A [37].
We can view a horizontal slice containing A vertices as having uniform density, obeying the dynamics shown in Fig. 2, with a mean addition rate Aðρ; α A ÞJ α A N , and a mean removal rate of Rðρ; α R ÞJ −α R . The first criteria for a phase separation to occur is mass balance, i.e., no net particle transfer, between two horizontal slices of densities ρ 1 , ρ 2 : where we have definedÑ in order to separate the dependence of J and ρ. Assuming each slice is large (A → ∞), we may find a recursion relation for the equilibrium distribution of densities: giving conditions for free energy balance The continuum model predicts the system to form walls whenρ falls between ρ 1 ; ρ 2 which satisfy both mass and free-energy balance. Note that the wall phase separation is independent of J. z translation symmetry.-To predict channelization, we assume translation symmetry in theẑ direction, so Jðx; zÞ ¼ κ(ρðxÞ)ẑ. As before, each vertical slice will obey the dynamics in Fig. 2, except the mean removal rate is now Rðρ; α R ÞκðρÞ −α R ; the mean addition rate becomes N Aðρ; α A ÞκðρÞ α A . Following the same procedure, the criteria for mass balance becomes and the condition for free energy balance becomes Z ρ 2 We note that when α A ¼ −α R , the criteria for wall building and channelization become identical; if a phase separation occurs, it will occur in both orientations, giving rise to a clumping phase separation. This is consistent with our observations in simulations Figs. 1(b) and 1(c). Our model for active heterogeneous conductors relies on a small number of very simple elements: a conserved current in a medium where flow and resistance are coupled to each other, in the presence of noise. This allows us to provide a unified coarse-grained approach that links a number of different physical and biological systems with different underlying fine-grained mechanisms that are often considered disparately. Despite its minimal complexity, our numerical simulations of the model produce the same channeling and immuring phase separations through the variation of only two parameters corresponding to biases in material addition and removal, consistent with the biases of the systems it is inspired by; flow-avoiding drainage networks and ant corpse piles lead to channels, and flow-seeking fuse breaks that lead to walls. A complementary continuum model corroborates our numerical simulations and also leads to the formation of walls, channels, and clumps. However, the functions characterizing conductivity and activity used in this continuum model come from numerical experiments which neglect microscopic correlations, resulting in a prediction of a continuous phase transition as opposed to the observed discontinuous one. In addition, because there is no inherent length scale to the continuum model, it cannot explain the transition between thin and thick structures. A continuum model considering the formation of the thinnest structures also predicts channelization and walling [38], but a theory combining both elements has no additional predictive power. Additionally the model predicts scale-free dipolelike correlations observed in the discrete model, which ought to exist in all nearly disordered systems where conductivity and current are coupled. While in these systems, addition and removal biases complement each other, our model shows that only one type of bias is needed for a phase separation. Opposing removal and addition biases, if they could be engineered, offer new possibilities for patterning that explore the entire phase diagram.
In the first section, we derive the strength of fluctuations for the low-⃗ α limit claimed in the paper. In the second section, we explain how the bimodality of a distribution was determined. All other sections are included for the sake of completeness. We start off from Darcy's Law and a Langevin equation,ρ Where T is the time derivative functional, ⃗ α is short for (α R , α A ) = and J is itself a functional of ρ obtained by solving (1). We decompose this: We note that, due to symmetry, the third term is zero and may be removed. Moving into fourier space, where ρ = ρ 0 + ∫∫ ρ(k)e ik·r , we we have: We now must find: To do so, we start off with a uniform density ρ 0 and then apply a sinusoidal perturbation ∆ρe ik·r . The conductivity is, to within first order Giving us a mean current We set ∆V to conserve current up to first order: Giving us our change in current, Plugging 3 into 2 yields:ρ where ⟨η(k, t)η * (k, t ′ )⟩ = 2δ(t − t ′ )D. The Einstein relation then predicts the strength of fluctuations to within first order:

II. DETERMINING IF A DISTRIBUTION IS BIMODAL
Each simulation i at a particular parameter value gives us a P i (N ), the probability of measuring a N particles in a row, column, or clump in simulation. Averaging many individual simulations allows us to calculate an average P (N ), as well as σ(N ), the estimated standard deviation of this average measurement. We use a 3-sigma threshold of statistical significance-we are significantly more likely to measure N particles than N ′ iff If a local maximum P (N ) can reach a larger P (N ′′ ) without moving through a valley where (4) holds, it is considered to be a false peak. If not, it is considered to be a true peak. A histogram with at least two true peaks is considered to be bimodal.

One true peak
Two true peaks FIG. 1 The distribution on the left has many false peaks, but is not considered to be bimodal. The right distribution is.

III. FOURIER TRANSFORM, TWO-POINT CORRELATIONS, AND CONDUCTIVITY
We can view the two-point correlations in real space and Fourier space, as well as the effect of bias on the conductivity of the medium.

IV. COMPARISON TO LOCAL DYNAMICS
For local dynamics, we we use a slightly modified time integration step. 1. Remove a particle from filled vertex i randomly selected with probability proportional to r(v i ). 2. Solve for the new current through the network.
3. Add the particle to an empty vertex j randomly selected with probability proportional to a(v j ) e − (r j −r i ) 2 2σ 2 . 4. Solve for the new current through the network. This prohibits a removed particle from traveling non-locally. Non Local σ = 1

FIG. 3 Comparison of Nonlocal and Local Dynamics
The phase diagram generated is very similar (Fig. 3).

V. ANALYTICAL MEAN FIELD THEORY
An alternate mean field theory produces some of the same qualitative behavior. Instead of relying on numerics to find the values of A(ρ, α A ), R(ρ, α R ), κ(ρ), we rely on a very simple model which gives analytical results.

Channels Walls
FIG. 4 Schematic of approximations made to obtain analytic forms for A, R, κ for channelization and wall-building Channels: For predicting a channelization phase transition, we assume that current is unable to travel in thex direction (Fig. 4). Therefore, the equation for conductivity is: Because current cannot flow laterally, an equal current of J is pushed through the filled and empty vertices, and so We note that is a function of α A + α R , and has no individual dependence on α A , α R .
Walls: For predicting a wall phase transition, we assume that, between rows, current can freely flow in thex direction without any resistance (Fig. 4). Therefore, the equation for conductivity is The total driving across a wall is J κ wall , and thus the current across an empty vertex is J κ wall . The total current across a filled vertex is J κ wall (1+∆Ω) . Therefore: We note that A wall R wall is also function of α A + α R , and has no individual dependence on α A , α R . Whenρ = .25, a channelization phase transition occurs when α A + α R ≲ −1.55. A wall-building phase transition occurs when α A + α R ≳ 2.25.

VI. SHORT LENGTH SCALE CONTINUUM MODEL
We can also create a continuum model on a short length scale. To do so, we select a periodic structure with two regions labeled 1 and 2. Region 1 comprises a fraction V 1 of the squares, while region 2 comprises a fraction V 2 of squares.
The density will originally be uniform, s.t.ρ 1 =ρ 2 =ρ, and a mean density of s (Fig. 5) can transfer between squares such that At a particular imbalance s, the probability of an particle moving from region 2 to region 1 divided by the probability of the opposite process gives us a "fugacity": N(ρ, s, α A , α R ) = R 2 (ρ 2 , κ(ρ, s), α R ) · A 1 (ρ 1 , κ(ρ, s), α A ) R 1 (ρ 1 , κ(ρ, s), α R ) · A 2 (ρ 2 , κ(ρ, s), α A ) Therefore, the free energy of a state with an imbalance s is − ∫ s 0 Ln[N(ρ, s ′ , α A , α R )]ds ′ s will be set to minimize free energy, and when the optimal s is nonzero, the continuum model predicts the system to spontaneously "crystallize" into a form where regions 1 and 2 have different density. For thin channels, region 1 is set by δ x mod d , where d some integer which sets the spacing between channels or pillars. This walls are the same except that region 1 is now set by δ y mod d . The short-length scale continuum model gives similar behavior to the uniform continuum model, although the change in free energy is nearly always lower.

VII. METHODS
The code used was written in a combination of C++, MATLAB, and Objective-C, and can be downloaded at http://web.mit.edu/socko/Public/PublishedCode/ActivePorousMediaCode.zip. Currents were solved using the Eigen