Topological properties of a self-assembled electrical network via ab initio calculation

Interacting electrical conductors self-assemble to form tree like networks in the presence of applied voltages or currents. Experiments have shown that the degree distribution of the steady state networks are identical over a wide range of network sizes. In this work we develop a new model of the self-assembly process starting from the underlying physical interaction between conductors. In agreement with experimental results we find that for steady state networks, our model predicts that the fraction of endpoints is a constant of 0.252, and the fraction of branch points is 0.237. We find that our model predicts that these scaling properties also hold for the network during the approach to the steady state as well. In addition, we also reproduce the experimental distribution of nodes with a given Strahler number for all steady state networks studied.

This model correctly reproduces the experimentally measured results, and also predicts the topological structure of the emerging network during the formation process. Surprisingly, we find that the observed steady state degree distribution relations are also obeyed during the approach to steady state.

Results
Physical model. We consider the case of i = 1, 2, … , N electrically conducting spheres of radius R. Since the current has been experimentally measured to be quite low 21 , the conductor interactions can be modeled using electrostatics alone. In the empty regions between conductors, the electric potential → V r t ( , ) is determined by the Laplace equation with appropriate boundary conditions.
The total charge on conductor i is related to the potential via Here the integral is taken over the surface S i of conductor i and n i is the surface normal vector of the surface element da i . The permittivity of the medium is ε.
The conductor charges are taken to be the boundary conditions for Eq. 1 with the exception of an additional conductor i = 0, the ground electrode, which held at a fixed position with zero electric potential. When two or more conductors are in electrical contact, they form a new conductor with a charge equal to the sum of all the conductors in contact. If a conductor is in contact with the ground electrode, it will also have zero electric potential.
The solution of Eq. 1 with these boundary conditions allows the calculation of the electric force → F i on the i th conductor as another surface integral of the form This is enough to write down the dynamics of the system: where m i , γ i is the mass and drag coefficient of conductor i, respectively. A charge source term J s has been added to account for the charge supplied to the system by external means. F(r i − r j ) is the contact force between conductor i and j. In this work we take the contact force between two spherical objects to be for r i − r j < 2 R, and F(r i − r j ) = 0 otherwise 22 .

Network analysis.
We computed the positions of the N conductors in the system as a function of time for nine values of N between N = 100 and N = 324. These numbers were chosen to be comparable to previous experiments 11,12 . A comparison between the steady state produced by the experimental system and the state produced by the numerical solution of the model after t = 120 s can be seen in Fig. 1. This time was chosen such that the system reached an approximately steady state in all cases. From the set of conductor positions obtained from the numerical solutions of the equations of motion, Eq. 4, an N node graph t ( )  was constructed that represents the electrical network at a given time t in which each node corresponds to one conductor. For the analysis, two conductors i, j are considered to be electrically connected at This is the connection criterion used in the experimental work 12 . The weight of the connection between nodes i, j is w ij = 1 if the conductors are electrically connected and w i,j = 0 otherwise. We then define the degree d i of each node i to be It is also possible to define an anti-arborescence  t ( ) rooted at the ground electrode. To do this, each conductor i is assigned a direct successor D i with the interpretation that the flow of charge in the network is from i to D i . The successors are then computed iteratively by defining the successor of all conductors connected to the ground electrode to be the ground electrode. In each subsequent iteration, the successor D i of the i th conductor is defined to be the nearest conductor that is connected to i and has a successor provided that i does not already have a successor. This process is iterated until no new successors can be assigned. Depending on the connectivity of  t ( ), it is possible that not all of the N conductors will be in t ( )  and so we will use M to denote the number of nodes in  t ( ). The total number of nodes which are directly or indirectly connected to the ground electrode and have degree We define a branch node to be a node i that has a degree d i ≥ 3. The total number of nodes B(t) which are branch nodes at time t can be calculated as with the successors D i defined, each node i can then be assigned a Strahler number s i . These are determined using the standard procedure. First, we assign s i = 1 for all  ∈ i t ( ) with d i = 1. Then for each node that doesn't have a Strahler number defined we set s i equal to the maximum Strahler number of the nodes which have node i as their direct successor. If more than one node with i as its direct successor have the maximum Strahler number, s i is incremented by s i → s i + 1. The number of nodes of with Strahler number j is σ j :

Comparison to experimental results. Experimental results suggest a linear relationship between
) 0 252 1 was found in steady state networks 11 . Figure 2 shows a plot of Δ 1 (t) vs. N computed from the model after t = 120 s of run time compared to experimentally observed steady state relationship.
In addition, the relation ) 0 237 was also observed in steady state networks 11 . Figure 3 shows a plot of B(t) vs. N computed from the model after t = 120 s of run time compared to experimentally observed steady state relationship.
A similar result was observed to hold for the distribution of nodes of a given Strahler number 12 . Experimental results show that σ 1 , σ 2 , σ 3 are each linearly related to N. The relations are σ 1 = 0.455 N, σ 2 = 0.275 N and σ 3 = 0.169 N. Figure 4 shows a comparison of the computed distribution of Strahler numbers as compared to the experimentally observed relations.

Discussion
It can be seen from the model that any stationary state of the system must correspond to a connected graph. This is because any conductor that lacks a connection (either directly or through contact with other conductors) to the ground electrode will eventually experience a force directed towards the ground electrode or another conductor in contact with the ground electrode due to the accumulation of charge from the source term J s . Only in the event that all conductors have a connection to the ground electrode do the forces on the conductors vanish.
Experimentally, it has been noted that the stationary networks rarely have closed loops 11,12 . This may be due to the form of Eq. 3, which shows that the force on any surface element of a conductor is normal to the surface and directed outwards. A closed loop can be thought of a single conductor with electric field → = E 0 inside the loop. Any such loop may experience a force that acts to expand the loop, and thus separate the conductors that comprise it. This force can only be zero in the case that → = E 0 everywhere on the outside surface of the loop as well. Therefore, closed loops in the conductor network are at best neutrally stable, and unstable in the presence of any external electric field.
The stationary states of the system are thus expected to be connected acyclic graphs, or trees. For any M node tree  , we can define the fraction f i (t) of nodes in  that have degree i as Second, there is a constraint related to the number of branches and endpoints. This can be interpreted as a statement that every branch in the tree ultimately has an end associated with it, and the number of additional branches created by a node of degree i is i − 2. This can be expressed in terms of the fractions f i as For spherical particles in 2D, the maximum degree any node can have i max = 4, as any degree larger than this would be a closed loop if connections are defined using 5. Thus there are only two independent numbers that specify the degree distribution of the network. The degree distribution can be expressed in terms of the fraction of endpoints f 1 and the fraction of branch nodes b = f 3 + f 4 as follows.
Using the experimental values of f 1 = 0.252 and b = 0.237 we obtain the full degree distribution = .
= . We note the apparent similarity of this degree distribution with that minimum spanning trees for random sets of nodes 13 which have f 1 = 0.253, f 2 = 0.527, f 3 = 0.204, and f 4 = 0.016. This may suggest that the network evolves in accordance with some optimization principle, such as the resistance minimizing principle 16 , or the maximum entropy production principle 17 .
The model is able to reproduce the experimental degree distribution for self assembled electrical networks starting from the physical interaction between conductors as well as correctly explain the rarity of closed loop structures. In addition, the model predicts that the degree distribution of the network remains constant during the approach to the steady state. We are also able to reproduce the experimentally observed distribution of Strahler numbers in the network. Methods Calculation of electric potential. Solutions to Eq. 1 were generated in parallel with the red-black Gauss Seidel method 23 . The region of interest was a 2D square L = 5.12 cm on each side which was discretized into a square grid with 512 × 512 sites. Neumann boundary conditions were taken at the edges of the square region, and all sites within conductor i = 0 were held at zero potential. In addition, a constraint was imposed such that the total charge on each conductor with i > 0 remained fixed at it's known value from Eq. 4. This constraint corresponds to finding a set of potentials {φ i } on the remaining conducting surfaces which were used as Dirichlet boundary conditions. This is accomplished by first defining the function i is the solution of Eq. 1 subject to the Dirichlet boundary conditions given by the set of φ i 's. Equation 16 gives the charge on the i th conductor if all the conductors were held at the known fixed voltages {φ i }. With this, we construct the functions Along with the total error associated with these initial conditions The problem of solving for the conductor voltages is now a root finding problem which can be solved by gradient descent. Gradient descent for this problem would involve making the update for some positive γ i . However, this update requires N sums over N terms, and also requires knowledge of the coefficients c ij . Instead, we make a local update by considering only one term of the sum Here we have made the redefinition γ i c ii → γ since c ii is always positive as well 24,25 . This update does not require the sum over N terms, and does not explicitly require knowledge of c ii . It only requires that γ be chosen to ensure convergence. At each iteration, the conductor potentials are changed by an amount In general this can be either positive or negative. The total error E then changes on the (n + 1) th iteration by Since c ij is a positive matrix 24,25 , the quantity ∑ = q c q i j N i n ij j n , 1 ( ) ( ) is always positive. Therefore, the total error is always decreasing provided γ is chosen such that Thus the update rule in Eq. 20 converges provided γ is chosen to satisfy Eq. 24. In this work a value of γ = 7.95 nV/C was found to be sufficient for convergence. This method was iterated until a solution with E/N < 2.3 pC was obtained.
Integration of the equations of motion. From the solution to Eq. 1 on the discretization grid, a bilinear interpolation was used to compute the potential between grid points. This is given by Scientific RepoRts | 7:41621 | DOI: 10.1038/srep41621 int for x, ∈ y (0,511). Here x is in units of grid spacing (px). This interpolation allows the computation of the electric potential in the area between grid points by approximating the behavior in the local region as changing linearly in x and y between the known values on the grid points. The surface integral in Eq. 3 was calculated from this interpolation using the trapezoid method of integration with N s = 500 evenly spaced points around each conductor. Here the charge is computed as the scalar sum of the charge on each of the N s discrete surface elements.
The equations of motion in Eq. 4 were numerically integrated using Euler's method with a timestep of Δ t m = 0.01 ms for the mechanical degrees of freedom and Δ t e = 0.001 s for the electrical degrees of freedom.
The mass of all conductors was set to m i = 16 g. The fluid drag was assumed to be Stokes drag, and so γ i = 6πμR with μ = 650 cP for castor oil 26 . The permittivity of the medium was set to be ε = 4.7ε 0 , where ε 0 is the permittivity of free space. This is near the dielectric strength of castor oil 27 . For the contact dynamics we used i j i j 6 1/2 3/2 . Conductors were considered to be in contact if |r i − r j | < = 2 R + 2 px. The charge source was J s = 3.8 pC/s. The conductors were initially laid out in a square grid centered in the region of interest with a separation between their centers of 3 R and the radius of all conductors was set to be R = 0.8 mm. The ground electrode position was fixed at (x 0 , y 0 ) = (496, 257). The position of the ground electrode is shifted one grid site off of the center line in the y direction to explicitly break symmetry.