Structure formation in suspensions under uniform electric or magnetic field

The structure formation of particles with induced dipoles dispersed in a viscous fluid, under a spatially and temporarily uniform external electric or magnetic field, is investigated by means of Brownian Dynamics simulations. Dipole–dipole interactions forces, excluded volume forces and thermal fluctuations are accounted for. The resulting structures are characterized in terms of average orientation of their inter-particle vectors (second Legendre polynomial), network structure, size of particle clusters, anisotropy of the gyration tensor of every cluster and existence of (cluster) percolation. The magnitude of the strength of the external field and the volume fraction of particles are varied and the structural evolution of the system is followed in time. The results show that the characteristic timescale calculated from the interaction of only two dipoles is also valid for the collective dynamics of many-particle simulations. In addition, the magnitude of the strength of the external field in the range of values we investigate influences only the magnitude of the deviations around the average behavior. The main characteristics (number density of branch-points and thickness of branches) of the structure are mainly affected by the volume fraction. The possibility of 3D printing these systems is explored. While the paper provides the details about the case of an electric field, all results presented here can be translated directly into the case of a magnetic field and paramagnetic particles.


Introduction
When dielectric or conductive particles are exposed to an external electric or magnetic field, dipoles are induced to the particles due to the difference in the dielectric permittivity or magnetic susceptibility between them, and the medium (Jones 1995). Due to these dipoles, the particles interact (Jones 1995) and the result of this interaction is the formation of structures of particles (Klingenberg et al. 1989). Depending on the field conditions, the particles can arrange into strings (Martin et al. 1998), that evolve, and aggregate in time (Martin et al. 1998), planes (Martin et al. 2000), or networks (Martin and Snezhko 2013). In addition, the short-time aggregation dynamics has been studied (Dominguez-Garcia et al. 2007;Promislow et al. 1995), as well the manipulation of a single particle (Lee et al. 2020). The enhancement of the properties of the material in the direction along which the structures are created originates from the arrangement of the particles (Martin et al. 2000;Martin and Gulley 2009). The transport properties are greatly affected by the structures formed, due to enhancement by percolation paths (Martin and Snezhko 2013). The suspensions that contain filler particles responsive to an external field, are called electro-/magnetorheological fluids and have been a subject of research due to their unique characteristics. Their rheological response occurs in O(1s) timescales and is exploited in force-feedback sensors (Zhang et al. 2018) and in robotics (Karasawa and Goddard 1989). The rapid viscosity increase in the direction perpendicular to the applied field has also been widely investigated in the past (Klingenberg et al. 1989;Brady 1990, 1992;Klingenberg et al. 1991a, b;Mohebi et al. 1996). The response of these systems under shear has been a topic of extensive research (Ginder and Davis 1994;Satoh et al. 1998). The effect of the Mason number (Melle et al. 2000), which is the ratio of viscous to magnetic forces, has been studied for the investigation of mixing on micrometer scales (Gao et al. 2012;Melle et al. 2003;Calhoun et al. 2006) with a rotating field. The behavior of such systems under the influence of an oscillating magnetic field has been investigated extensively over the past decade (Sánchez and Rinaldi 2010;Liu et al. 2019;Soto-Aquino and Rinaldi 2015;Ido et al. 2016;Jonasson et al. 2019;Ruta et al. 2015). The structures created in a spatially uniform external field have also been studied at large timescales both athermally (no Brownian motion) (Martin et al. 1998), and thermally (including Brownian motion) (Martin 2001).
Additive manufacturing is nowadays a subject of extended research (Ngo et al. 2018), although it has been on the surface for a long time (Yan and Gu 1996). The materials used include metals (Visser et al. 2015), polymers (Ligon et al. 2017), ceramics (Warnke et al. 2010) and combinations of them (Ngo et al. 2018). Printing of polymer-matrix composite materials is used for various purposes such as rapid prototyping (temporary substitutes for parts of higher mechanical performance) (Czyzewski et al. 2009), or topology optimization (introduction of shape memory in a passive matrix) (Maute et al. 2015). They have a wide range of applications: aerospace applications (cabin interiors) (Raja et al. 2010), medicine (mimics of living tissue) (Villar et al. 2013), anthropology (reconstruction of medieval skulls) (Massimiliano and de Crescenzio Francesca 2008) and design (spatially dependent elasticity) (Oxman et al. 2012). Another reason for using composite systems is the incorporation of dielectric, magnetic or conductive functionality to the matrix material (Castles et al. 2016;Kokkinis et al. 2015;Czyzewski et al. 2009). Additionally, if the microstructure of the composite can be controlled during printing by the use of an external electric or magnetic field, one can achieve specific structures of particles inside a complex geometry. This can be an important application of this work, especially if photo-reactive resins (Anastasio et al. 2018) filled with electrically of magnetically active particles are concerned. Resins are widely used as media for particles (Kim et al. 2003;Kim and Shkel 2004), so that the structure of the particles is captured with the solidification of the resin. The low viscosity of the resins before curing makes the structure formation possible within the usual timescales (O(1s)) of 3D printing. The targeted technique is stereolithography (SLA) (Bártolo 2011). Potential uses of these systems in applications include personalized hearing aids (Dodziuk 2016), flat lenses with a gradient of concentration of particles that have the functionality of their curved counterparts (Kurochkin et al. 2018;Viskadourakis et al. 2018), piezoelectric or Hall effect sensors (Quanlu 2002), and direction-specific thermally or electrically conductive composites (Martin and Snezhko 2013;Subramanian et al. 2019). In the following, we focus on the physical behavior of these systems, in terms of structural evolution.
The motion of particles serving as dipoles inside a fluid can be simulated with a variety of methods depending on the effects (hydrodynamic interactions) and scales (length or time) that one desires to resolve. The most accurate technique is molecular dynamics, as it simulates the motion of the atoms (Frenkel and Smit 2002;Vogiatzis and Theodorou 2014) by integrating the (frictionless) Newton's equations of motion. However, the system we are interested in, exhibits disparity in length-and timescales, as it consists of monomers of the liquid matrix whose average size is on the order of nm and particles whose size is on the order of μm. Molecular Dynamics can resolve the nano-scales in time and space. The effect of the matrix on the particles is governed by the collision of the monomers on the surface of the particle, and can be effectively represented by thermal fluctuations on the μm scale. Stokesian Dynamics (Bonnecaze and Brady 1990) is another option for simulating these systems and especially resolving the hydrodynamic interactions at close distances. One of the best known techniques is the Brownian Dynamics (BD) (Soto-Aquino and Rinaldi 2010), where the equations of Newton are solved in the overdamped limit, so no acceleration is present (i.e., correlations in velocity vanish within a single timestep of the integration) (Dhont 1996). However, the great advantage with respect to the system in discussion is the incorporation of thermal fluctuations. BD can resolve spatial effects on the order of the size of the particle and times of up to O(10s). The Finite Element method (Gao et al. 2012;Kang et al. 2008) is also suitable for solving the equations describing the problem, although with this method the system is treated either as a macroscopic medium, or at a much finer spatial resolution than the size of the particles. The macroscopic approach is totally neglecting the microstructure, i.e., solving the equations of Maxwell for resolving the electromagnetic response of the material and/or the Navier-Stokes equations for the rheological response. The finer resolution approach solves the same equations implementing the mass and momentum balance equations, but also resolves the local effects of non-uniformities of the electromagnetic field and/or the flow field, respectively. The downside of this method is the "limited" number of particles that can be studied due to the computational cost. To study the structures in short/intermediate timescales, the BD technique is used. The motion of the particles is mainly affected by the dipole-dipole interactions (Jones 1995;Klingenberg et al. 1989;Parthasarathy and Klingenberg 1996).
Our method results in a structure of particles, and our characterization depends on the nature of these structures. As already mentioned, one can get various structures depending on the field conditions, like strings (Martin et al. 1998), planes (Martin et al. 2000), or networks (Martin and Snezhko 2013). The goal of characterizing such structures is to extract their main features and to describe these features in a quantitative manner. For example, for strings, one should know the size, orientation, thickness and correlations between these measures, for planes the size and orientation, and for networks the number density of branch-points (BP), the degree of BP, the thickness of the branches, the number of clusters present, and the existence of percolation. There are a lot of techniques used for the characterization of an ensemble of points/particles, like studying their gyration tensor (eigenvalues and eigenvectors)  and Voronoi tesselation (Voronoi 1908). The Voronoi tessellation is widely used in the characterisation of disordered systems (Montoro and Abascal 1993), polymers (Greenfield and Theodorou 1993;Starr et al. 2002;Damasceno et al. 2012;Vogiatzis and Theodorou 2014), colloidal gels (Varadan and Solomon 2003), and granular materials (Li and Li 2009). Other techniques like the (radial or cylindrical) pair-correlation function, bond-angle distribution (Hütter 2000), and quermass integrals (Hütter 2003) have also been developed for the characterization of the structure of networks. Here, the gyration tensor of clusters, the Voronoi polyhedra, and a previously developed skeletonization method (Kerschnitzki et al. 2013;Kollmannsberger et al. 2017;Manikas et al. 2020) are employed for the quantitative characterization of the structures.
In this paper, we use BD with dipole-dipole interactions to simulate the motion of particles (dipoles) within a fluid under a spatially uniform field. The characterization of the structures is performed by grouping particles in clusters and by reducing the structure to its skeleton (skeletonization). Our goal is to identify the most relevant physical parameters, and relate the corresponding physical input with the structure and eventually the transport properties, specifically the thermal conductivity of the material. The paper is organized as follows. Section 2 provides the basic methodological tools that were used for the production of the structures. Moreover, we present a dimensionless analysis with all the relevant dimensionless groups. In Sect. 3, we introduce the tools we use to characterize the obtained structures. In Sect. 4, we present the structural evolution in time in dimensionless terms, and the characteristic features of our formed structures. Finally, the paper is concluded with a discussion in Sect. 5.

Methodology
In this section, our methodology is discussed, including the simulation details and dimensionless analysis. Our system consists of particles with particle radius R p and volume fraction φ, and the monomeric fluid that serves as medium with constant viscosity η in a finite box of edge-length L, where periodic boundary conditions are applied in all three spatial directions. The system is influenced by an external field that is applied to the system. In the current paper, the system consists of particles that have dipoles induced by a uniform external electric or magnetic field. The interactions governing the dynamics are the dipole-dipole interactions (Jones 1995;Klingenberg et al. 1989). One can define the potential energy of N particles with induced dipoles under an external field (Jones 1995;Gao et al. 2012), where p i is the dipole moment of particle i, defined in Eq.
(2), E e (r i ) is the field intensity at the position r i , m is the dielectric permittivity of the medium, C(r) = ∇ r ∇ r 1 r = 1 r 3 ( 3 r 2 rr − I) is the second derivative of the position with r = |r|, and r i j = r i − r j is the inter-particle separation vector. The appearance of dipole moments is a result of the polarization induced by the external field (Jones 1995) where K = p − m p +2 m is the Clausius-Mossotti constant (Böttcher 1973), with p the dielectric constant of the particles.
The gradient of the potential energy with respect to the position r i , results in the expression for the forces In Eq.
(3), one can observe that dipole i is different from dipole j, however, in our study, we only considered dipoles that are induced by the field, so all the dipoles are identical (p i = p j = p) since no local polarization is considered, only a spatially uniform field is studied, and only particles are used that are identical in size, shape, and dielectric permittivity (or magnetic permeability). For the rest of this paper, we are using different indices for different dipoles, as it corresponds to the general case.
The presented form of the energy, Eq. (1), is the pointdipole approximation (Calhoun et al. 2006;Keaveny and Maxey 2008). The validity of this approximation depends on the dielectric mismatch p / m , and the volume fraction (Yang et al. 2006); the lower these two values, the better is the approximation. In the current paper, we neglect the multipolar interactions (Keaveny and Maxey 2008;Jones 1995) that become significant at close distances (r i j < 4R p ). At these distances, the dipole approximation is not adequate, as it neglects the spatial variation of the field around the particles, which at close distances is significant and the importance of higher order multipoles increases. The dipole approximation is also ignorant of the polarization arising from the neighboring dipoles; the dipoles can influence the field intensity locally, causing an alternation of the polarization of the neighboring dipoles. Mathematically, this is considered in the mutual dipole approximation (Keaveny and Maxey 2008), where the dipole moments depend on their relative positions of the dipoles within the system. Despite the limitations of Eq. (1), it is still used it in this paper to simulate the structure formation of particles under a spatially uniform field for two reasons. First, the systems considered in this paper have a low dielectric mismatch, which will be absorbed in a dimensionless group further below. And second, this approximate interaction is appropriate at the low volume fractions studied here, while at the higher values of volume fraction, it nevertheless can give a qualitative picture of the overall behavior of the many-particle system. Also, it is mentioned in the literature that the calculation of the multipole interactions is computationally inefficient with respect to the accuracy obtained in comparison to the simple dipole approximation (Keaveny and Maxey 2008).
In Eq.
(1), we present the potential energy in terms of electric dipoles. However, if one uses magnetic dipoles, m instead of p, the equivalent expression of the field intensity, H e instead of E e , and adjusts the prefactor of the second term of the right-hand side of Eq.(1), m replaced by 1/μ m , one directly obtains the equivalent case for magnetic dipoles in the presence of a magnetic field (Gao et al. 2012). The nature of the interactions is similar when the particles used exhibit magnetization that remains parallel to and linear in terms of the applied field, e.g., paramagnetic particles, and therefore, in terms of dimensionless parameters the difference between the two cases becomes obvious. However, one should be careful with ferromagnetic behavior (Morimoto and Maekawa 2000), as it is not included in the potential energy presented in Eq. (1).
Due to the finite size of the simulation box, every interaction is truncated at a certain cut-off radius, r cut < L/2, with L being the edge-length of the simulation box. The problem of finite truncation of the interactions is dealt with by calculating the forces neglected outside a sphere of radius r cut , and correcting the final result (tail corrections). If one deals with electrostatics or pairwise potentials that have n < 4 in 1/r n in terms of energy, c.f. Eq. (1), then the neglected energy is infinite and one should account for long-range corrections with other methods; the same translates to the forces as well.
In the case of the presence of partial charges, this problem is known for almost a century (Ewald 1921), and has been dealt with the method of Ewald summation. The Ewald summation resolves the issue by turning part of the forces (and/or potential) in the reciprocal space with a Fourier transform, taking advantage of the periodicity of the system. No estimations are made, the solution results from an analytical transformation. Other techniques have been developed like the reaction field method or the Wolf summation (Allen and Tildesley 1987;Fukuda and Nakamura 2012), however, these techniques treat the system as a medium making approximations in the calculation of the quantities, which exhibit systematic errors (Allen and Tildesley 1987). In reciprocal space, the sum converges much faster with the use of wave vectors. Studies about the accuracy of the method have been performed both for charges (de Leeuw et al. 1980;Perram et al. 1988;Karasawa and Goddard 1989), and dipoles (Wang and Holm 2001). In this paper, the Ewald summation method is used to account for the long-range corrections of the forces, and we use the parameter α = 7.5 L −1 that defines the split between real and reciprocal space, where L is the box length, and k = 10, where k is the amount of wave vectors used.
The potential that was chosen is purely attractive and there is need for an extra force that will prevent the particles from overlapping. This problem was encountered in the past (Klingenberg et al. 1989). The choice is usually made between a power law ( 1/r 12 ), and an exponential ( e −κ(r i j /(2R p )−1) ). Our choice is the exponential, Eq. (4) and has to do with the particle configurations observed experimentally. The power law alters the potential in close distances preventing particles from formatting chains that are thicker than one time the particle diameter (Klingenberg et al. 1989). The form of the force for preventing particles from overlapping is where κ = 30, is a constant that determines the interaction range of the excluded volume, while the pre-exponential factor defines the overall strength of the interaction (Gao et al. 2012), andr is the unity vector of r. One can observe that the softness of the particles (hard particles are difficult to handle numerically) depends on the magnitude of the dipole moments, which is rather unexpected due to the different physical origin of the excluded volume and dipole-dipole interactions. This forcefield has been used extensively in other studies in the literature (Martin et al. 1998;Mohebi et al. 1996;Gao et al. 2012;Melle et al. 2003). As the excludedvolume force (4) depends on the dipole moments, it vanishes for infinitely weak dipole moments; however, the expression (4) is suitable for the cases in this paper, where dipole moments always have finite values, because there is always an imposed external field that induces the dipoles. Additionally, the magnitude of the excluded volume force is always larger than the other physical phenomena, i.e., Brownian and dipole-dipole interaction forces; therefore, it produces results very similar as if an alternative excluded-volume force was used. The motion of a particle inside a fluid is influenced by friction forces, originating from the collisions of the molecules of the surrounding fluid with the particle. These forces can be incorporated by a drag force, where ζ = 6πηR p is the friction coefficient, and v i is the velocity of particle i. In our case, we use the Stokes drag force (Larson 1999), while many-particle hydrodynamic interactions are neglected. The reason for neglecting the latter is that otherwise the computational cost would be prohibitive, in particular, since the main goal of this paper is to conduct an extensive study of the effect of dipole-dipole interactions and thermal fluctuations on the structure evolution and to resolve the physical behavior of the system. We expect that correctly accounting for hydrodynamics will influence primarily the time scale, however, there could also be an effect of hydrodynamic interactions on the structure formed. Lubrication forces are much larger than those predicted by Stokes' law, but are only important at small particle-separations where the particles are already locked into their local structure by the dipole-dipole forces. Hence, lubrication effects are not expected to significantly influence the structures formed (Klingenberg et al. 1989). The Brownian force is given by where k B is the Boltzmann constant, T is the temperature, and W is a Wiener process and dW its increment (Öttinger 1996). This has the properties dW(t) = 0 and dW(t)dW(t ) = δ(t − t )dtdt 1. So its increment can be expressed as dW i = √ dtξ i , where ξ i is a vector whose components are drawn from a random-number distribution with mean value of 0 and standard deviation of 1. The realization of the thermal noise through a Wiener process is such that the Brownian force has a mean value of zero, and it is uncorrelated in time.
For all the above, we have chosen to use Brownian dynamics (BD) simulations (Dhont 1996). This kind of simulations is an approach to the mathematical modeling of molecular systems by the use of stochastic differential equations (Öttinger 1996). The effect of inertia is considered negligible for systems where the relaxation time for the particle velocity, τ p = m/ζ , with m being the mass of the particle and ζ the friction coefficient, is much smaller than the other timescales.
To reduce the set of physical parameters that influence the behavior of the system, and to introduce dimensionless quantities with physical meaning, we re-write Eq. (7) in dimensionless form. The procedure that is followed begins by introducing scaling quantities for the principal variables of the problem, Here, n d is the number density, and the index "sp" refers to the spatial average over the simulation box. For the lengthscale r c , our choice is based on the system properties. We would like the lengthscale to be independent of the system size (box length), and the size of the particles (radius or diameter), so we select the third root of the average volume per particle in the system. This quantity is a measure of the inter-particle distances in the system, and independent of extensive variables of the system. For the scaling of the dipole moments, we choose the dipole moment that corresponds to an induced dipole with the characteristic field intensity value. For the characteristic value of E, we choose the spatially ("sp") averaged field intensity in the simulation box. Finally, the timescale is chosen such so that the electromagnetic effects occur on a dimensionless timescale of order unity in Eq. (7), so the prefactor of this term is set to 1 (Dantzig and Tucker 2001), see Eq. (12). It has to be noted that the excluded volume forces have the same scaling with the dipole-dipole interactions (square of the magnitude of the dipole moment). The characteristic timescale is thus given by As a result, Eq. (7) turns into (Öttinger 1996;Gardiner 2004) with If there was an externally imposed flow, e.g. imposed simple-shear deformation, another characteristic timescale associated with that flow would enter the analysis. However, imposed flow is not considered in this study. Furthermore, we neglect hydrodynamic interactions, i.e., the effect of the flow created by the movement of one particle on the neighboring particles is neglected (see also discussion after Eq. (5)).
The BD algorithm employed in our work has already been used in several studies before (Hütter 1999;Zakhari et al. 2017Zakhari et al. , 2018a, without the dipole interactions. Here, the dipole-dipole interactions are tested by considering some limiting cases. The alignment of two particles in the absence of Brownian forces is tested and compared with the analytical solution. In Fig. 1, one can see the interparticle distance versus time for an attractive (interparticle vector aligned with external field) and a repulsive case (interparticle vector perpendicular to the external field) for a simulation of two particles with no Brownian motion and no excluded volume interactions. Good agreement between the BD simulations and the analytical solutions is obtained. In "Appendix B", the simulation of many particles is shown in the absence of dipole interaction, but with Brownian motion and excluded volume interactions included. The obtained structures agree with the expectation.
For solving Eq. (13), one should use time discretization, with a finite value for the time increment dt. The timestep is set, so that important information of the fastest physical process is preserved. In the current study, the fastest process is governed either by the dipole-dipole interactions or the excluded volume interactions. At large distances, the forces decay, so one should consider only the case with the largest forces, which is the case of particles at small distance. In this case, we are going to consider r * i j = 1.9R p , in view of thermal fluctuations. For determining which of the two interactions is dominant, we check Eqs. (14), (15). If one calculates the char- If one estimates these timescales for φ = 30%, and the chosen value of r * i j , then one gets t * c,dip = 3.87607 × 10 −1 and t * c,exv = 5.30915 × 10 −2 . The timescale concerning the excluded volume interaction is faster, so one should consider the excluded volume effect for the choice of timestep. In the dimensionless form of the model used here, the resulting timescale corresponding to the excluded volume depends on the average characteristic distance between the particles, so one should estimate the timescale value for different volume fractions. Taking this fact into account, we select a timestep of dt * = 2.4 × 10 −4 , which is sufficient for all volume fractions considered in this study. The details of this choice are presented in "Appendix A". It has to be mentioned, that the tools we apply only need the position vectors (Cartesian coordinates) of the particles.
Most measures concern two particles, see Sect. 3.1. Although most of the times two-particles measures are not enough to resolve the whole structure characteristics, so one can look either in the neighbourhood (nearest neighbours), see Sect. 3.2, or clusters/structures of particles, see Sect. 3.3. If one has in mind particles organised in strings, then one would think in terms of touching particles (bonded) creating an anisotropic structure. This characteristic of the structure led to the distinction of the particles to clusters. The condition for characterising two particles as bonded is r i j ≤ 2.2R p . This choice was based on the ratio U exv /k B T , where at distance of 2.2 R p this ratio is lower than 10%. This means that the effect of the excluded volume has ceased and the displacement due to Brownian motion is dominant at this limit.

Measures for pairs of particles
There exist a lot of measures concerning two-particle properties, the most well known being the pair-correlation function (Allen and Tildesley 1987). However, due to the nature of our structures (anisotropic) in this paper, we discuss a measure of the orientation of the inter-particle direction, S 2 , with respect to an external direction, that being the external field, E This is the second Legendre polynomial for every pair of particles in the system independent of the distance between them, however, periodic boundary conditions are applied and only the primary image is considered. Eq. (17) is sensitive to the orientation of the inter-particle distances with respect to the direction of the field. S 2 measure has a value of 1 if the average orientation is parallel to the external direction, a value of 0 if there is a random orientation with respect to the direction of the field, and a value of −1/2 if the orientation is perpendicular to the external direction. If two, parallel with the field, one-particle-thick strings are present, then the value depends on the distance between them. However, despite the fact that they are ideally oriented the value of the measure is lower than 1, as the inter-particle vectors between particles belonging to different strings are also considered.

Measures for nearest neighbours of a particle
If a neighbourhood of particles is to be considered, then the average volume per particle comes up naturally. The Voronoi polyhedron of a particle is uniquely defined by the positions of its neighbours. This quantity is the essence of the calculation of the Voronoi polyhedra (Voronoi 1908).
The Voronoi tessellation is widely used in the literature Vogiatzis and Theodorou 2014;Damasceno et al. 2012;Starr et al. 2002;Montoro and Abascal 1993). The result of the tessellation is polyhedra. The volume, shape (number of faces), and orientation of these polyhedra can give us insight on characteristic features of the underlying particle structure. The characteristics of the polyhedra depend on the relative positions of their first neighbours. We choose to use Mitrich's method (Mirtich 1996), to calculate the gyration tensor of the polyhedra, and extract information as the main direction of the polyhedron and the shape of it. The gyration tensor where the brackets indicate the arithmetic average over the vertices k of the polyhedron of particle i, and r cm,i is the position of the center of mass of the Voronoi polyhedron of particle i. The gyration tensor defines an ellipsoid . The relative magnitude of the eigenvalues of the gyration tensor give the relative difference in size of the three axes of the ellipsoid.
The measures that will be used for characterising the polyhedra out of the gyration tensor and conclusively the structures, are the volume of the polyhedron, V i , and the relative shape anisotropy Damasceno et al. 2012) where b = (λ 1 − 1 2 (λ 2 + λ 3 )) is the asphericity of the ellipsoid defined by the gyration tensor, c = (λ 2 − λ 3 ) is the acylindricity, and the eigenvalues of the gyration tensor S are in descending order, λ 1 ≥ λ 2 ≥ λ 3 . The value of κ 2 for a sphere is 0, for an infinitely long rod is 1, and for a flat disk is 1/4.

Measures for clusters/networks of particles
Here, we discuss morphology measures concerning clusters/structures of particles. First, we introduce the gyration tensor for clusters of particles where the brackets indicate the arithmetic average over the particles i, the index "I " indicates the cluster number, and r cm,I is the position of the center of mass of the particles belonging to the cluster I . The gyration tensor defines an ellipsoid , and its geometrical features can give insight about the structure of the clusters.
The first measure is the same as in Eq. (17), the only difference being that the average is limited to the particles belonging to the same cluster where the index I denotes the cluster index, so for cluster I the average is defined as the average of the particles that participate in I . The second measure quantifies the anisotropy of each cluster. For that we calculate the gyration tensor . The gyration tensor is calculated per cluster, so each cluster can be described as an ellipsoid , where the eigenvalues of the gyration tensor give the relative difference in size of the three axes of the ellipsoid. Therefore, we introduce a measure that quantifies anisotropy λ * : where λ max is the largest eigenvalue of the gyration tensor, and det(S I ) is the determinant of the gyration tensor of cluster I . This measure has a value of 1 if the anisotropy is high, and a value of 0 if the system is isotropic (sphere). Note that this measure does not depend on the direction of the anisotropy. The third measure that we pick, has to do with the size of the cluster. This measure is, where N I is the number of particles that belong to the cluster I , N L is the minimum amount of particles that stacked together can span the edge-length of the box in one direction. This measure has a value of N * < 1 if the cluster size is smaller than the N L , and a value of N * ≥ 1 if the particlestring has enough particles to span the box, even if the box is not actually spanned in that specific configuration (flocculates). A schematic representation of the described measures can be found in Fig. 2. The number of clusters can be retrieved from N * , if one knows the volume fraction, φ, and the number of particles, N : Based on φ and N , one can calculate the ratio L/(2R p ), which is equal to N L , and identify the average number of particles per cluster. With this, in turn, the number of clusters is obtained if one divides the number of particles N by the average number of particles per cluster. In this paper, the dimensionless parameters used were chosen with regard to the characterisation of the resulting structures, e.g. oriented, anisotropic and percolating clusters. Focus is intentionally given on single cluster properties in normalized terms from which other quantities can be inferred. The averages of the measures S 2,I , λ * I , and N * I are reported in the following, and these averages will be denoted by the same symbol as their per-cluster counterparts, but without the cluster index. The size of the cluster serves also as a weighting factor for these averages, so that every particle of the structure has equal effect on the measured quantity. However, it is much different from calculating these quantities for the entire system.
In previous work of ours (Manikas et al. 2020), we introduced a method to characterise large structures of particles, including networks. That method uses a 3D-binary image of the initial structure, and thins out this structure to an infinitely thin skeleton. This procedure is called skeletonization (Kollmannsberger et al. 2017), and it results to an ensemble of voxels (skeleton); there are still complex connections in the skeleton representing the abstract connectivity inherent to the actual structure provided as input. This issue is tackled by several post-processing steps that result in a topologically equivalent skeleton with essentially the same connectivity (Manikas et al. 2020). The whole scheme we described, results to a simplified, and easier to characterise skeleton structure. The simplified skeleton is characterised in terms of: -branch-points (BP) (number density, degree) -branches (thickness) -existence of percolation Details concerning these points follows.
The BPs are identified by the number of their connections to neighboring voxels in the simplified skeleton. Every voxel bearing three or more connections is accounted for as a BP. Therefore, one can identify the number of BPs. A BP reduction step is required to prevent unphysical results of BPs lying in neighboring voxels. This step concerns BPs located inside the volume of a single primary particle, when BPs like the ones described are identified they are grouped together as one BP. Then, the number density of branch points per volume is given where N BP is the absolute number of BPs, and N p is the amount of primary particles. Except the number density, also the average degree of BPs can be studied, d BP . The degree is defined by counting the bonds of all the BPs in the simplified skeleton, and then taking their average. The average thickness of the branches is calculated by cutting a slice of the initial binary image at the position of each voxel and taking as the normal vector the orientation of the skeleton at this point (tangent of the skeleton) (Manikas et al. 2020). The percolation in each direction of the box is also checked. One tries to establish in-box pairing between any couple of particles each of them being close to two opposite faces of the simulation box ("close to" implying they have a bond crossing that respective face of the box).

Results
The results of particle-based simulations using the tools mentioned in Sect. 3, will be reported. At this point, we would like to note that the error-bars in this section will always refer to different random numbers sequences. Our system consists of 1000 particles, and our protocol indicates that we use three simulations with different random numbers sequences to produce the error-bars. In the rest of the section, we discuss the effect of the variation of random initial configurations, see Sect. 4.1, the B * variation, see Sect. 4.2, and the φ variation, see Sect. 4.3.

Variation of initial configuration
In this section, we study the effect that different random initial configurations have on the structure with the measures introduced before. We obtain the initial configurations by setting the particles to a simple cubic lattice, and letting them re-arrange into the box. We use a Brownian Dynamics scheme and the deterministic part for the equilibration loop consists of the repulsive part of a Lennard-Jones potential (Jones and Chapman 1924;Smit 1992) with σ = 2R p /r c , and ε = 1. This randomness of the structure has been investigated by examining the number density of particles in a tessellated box and the pair correlation function. The box contains 1000 particles, it is tessellated in 125 cubic boxes, and the number density is defined as the average number of particles that appear in the cubic box. The average number density over the cells is constant by definition, so we look at the standard deviation of this measure over the cells and Fig. 3 A random configuration of particles (red) is presented with their Voronoi cells (blue) produced using the Voro++ library (Rycroft 2009) through (dimensionless) time of 90t c . The samples are equilibrated, until the standard deviation between the cells shows a plateau with a value smaller than 0.2, so that a homogeneous distribution over the cells is reached. At the same point, the pair-correlation function is calculated and has the same characteristics as a Lennard-Jones fluid (Allen and Tildesley 1987). The details can be found in "Appendix B".
Once these random configurations are obtained, we study the impact of varying them in the structure evolution in terms of S 2 , λ * , and N * . We choose these measures mainly, because they are applied to the primary structure and no transformation of the structure has to occur (skeletonization). We used four different random initial configurations, each of them evolved by three simulation runs with different random number sequences. Since the initial structures have the same characteristics and the same average inter-particle distance, it is expected to see an impact on the same order of magnitude as the thermal noise (error-bars). Our expectation was not verified, as the variations observed were larger than the ones produced by different realizations of the thermal noise. This variation is not trivial, so for the rest of this paper we are using three different initial configurations in combination with the different random number sequences, for every set of parameters. The results can be found in the "Appendix B", as they are not essential part of the essence of this paper the physical behavior of our systems. A random configuration of particles with their Voronoi cells (Voronoi 1908) is presented in Fig. 3.

Variation of B *
The variation of the single physical dimensionless parameter, B * , will be discussed in this section. The variation of the parameter B * , has to do with varying physical parame-ters like the field intensity or temperature and many-particle simulations are needed to investigate these effects.
Different B * = 0.01, 0.0316, 0.1, 0.316 and φ = 1, 2, 5, 10, 20, and 30 values were used to unveil the effect of the only dimensionless parameter to the structure formation. The response of the system in terms of our measures (S 2 , λ * , and N * ), is expected to depend on the value of the parameter B * , as B * defines the magnitude of the thermal noise, see Eq. (13). However, B * is a parameter introduced based on two-particle interactions and we expect to observe some deviations in the generalisation to many-particle systems. In Fig. 4, one can see the response over dimensionless time for φ = 10 %, is overlapping to a single curve independent of the B * . This means that the characteristic timescale governing the twoparticle dynamics is of relevance also for systems with many particles, and no adaptation of the definition of B * is needed to study the structure formation. In the rest of this paper, we are going to use the value of B * = 0.01. Only one value of φ is presented here, although we used the same B * values for φ = 1, 2, 5, 10, 20, and 30, we display some characteristic values of low and high φ in "Appendix C".
If one would increase the value of B * further, e.g., to values as high as order unity, one should expect to see "deviations" from the pattern of Fig. 4, the reason for these "deviations" being the following. In order for the scaling with t c to collapse the simulation results obtained under different conditions, it is assumed that the dipole-dipole interactions are dominant over the thermal fluctuations; if this does not hold, the structure formation will be changed qualitatively. No hard limits can be defined, however, we believe that for B * > 0.5 the thermal fluctuations dominate, and a straightforward collapse of the results by scaling can not be achieved anymore.
It should be noted that the simulations presented in Fig. 4, were completed in dimensionfull units. The difference in dimensionless time comes from the fact that equal overall dimensionfull times were used. One should expect equal simulation times in t * for the different B * -values, however, our simulations results were independent of the B * -value, as they overlap onto a single curve. To this end, we concluded that conducting additional simulations would not be a significant contribution to the essence of this paper.

Variation of volume fraction
In this section, we would like to focus on the φ variation and its impact on the structure. Our goal is to determine the structures created in different systems, and explore the applicability of our dimensionless analysis, and morphology tools. The goal of this study is the determination of structureformation timescales and structure that can be achieved.

Characterisation of clusters
The volume fraction, φ, was varied (1, 2, 5, 10, 20, and 30) for the minimum B * value (0.01). We expect that structure formation will occur at the same dimensionless time as φ increases. This occurs due to the lower inter-particle distances met in higher φ-values, the timescale is decreasing, see Eq. (11). However, we also expect significant differences among different φ-values as our scaling was performed for two-particle interactions and in this section, we study manyparticle systems with different inter-particle distances. In Fig. 5, one can observe that all the φ values used show initial response in our measures at t * ∼ 0.05. That means that our expectation was correct and the structure formation for higher φ in terms of time is resolved from the scaling analysis. The difference observed for different φ implies that the two-particle scaling is not enough to resolve the structural evolution of many-body systems when the inter-particle distances are modified. That means that more complex physics take place. We also observe a significant drop in our measures as time increases especially for φ > 10%. This can be justified by the fact that our measures are applied to whole clusters, and all the particles create a single cluster in higher φ, see Fig. 5c. For this purpose, further investigation is required, to this end we are going to explore the neighbourhoods of particles.
At this point, it should be mentioned that the structure formation occurs in O(1s) for all φ-values, Eq. (12). This timeframe is suited for applications like 3D printing, and more specifically stereolithography. In recent studies (Anastasio 2019), it has been observed that the gel-point of acrylate systems is reached at times within the same order of magnitude.

Characterisation of neighbourhood
We study the neighbourhoods of particles through the Voronoi tessellation technique for the same configurations we investigated before. Except for the Voronoi volume distributions, the relative shape anisotropy of the cells, defined in Eq. (19), is also explored.
The main issue of our tools is the insufficiency for identifying the main characteristics of a structure that is formed by one cluster (probable network). We expect that the Voronoi cells of the particles are going to distinguish the core of the structure (V cor < V vor and isotropic shape). We speculate that particles belonging to the external layer of the structure will exhibit V ext > V vor , and pretty anisotropic shape (distant neighbours). Although, conclusions can be drawn about the structure, no identification of the core was made. In Fig. 6, one can see that the evolution of the structure in time, and the distributions are far from random (broader), no quantitative identification of the characteristics (BP, branches) is possible. A well-known (Voronoi 1908;Montoro and Abascal 1993) technique in many different fields (Greenfield and Theodorou 1993;Starr et al. 2002;Damasceno et al. 2012;Vogiatzis and Theodorou 2014;Varadan and Solomon 2003;Li and Li 2009;Duyckaerts and Godefroy 2000) has proven inefficient for the quantitative characterisation of our systems.

Characterisation of skeleton
Since, our attempt at identifying the main structure characteristics concerning networks or late stages of the structure formation did not yield the desirable results, we will focus at transforming our structure to its bare skeleton with our previously introduced (Manikas et al. 2020) skeletonization methodology. This technique gives us the possibility to identify the main characteristics of the structure and create a skeleton from which we can extract all the important information (BPs, thickness of branches, existence of percolation).   . 6 a The distribution of scaled Voronoi volumes is plotted as a function of dimensionless time for a system of φ = 10% and B * = 0.01. The scaling is performed with the average Voronoi volume of the configuration.The distribution of scaled Voronoi volumes (b), and relative shape anisotropy (c) are plotted for dimensionless time t * = 10t c and φ = 1, 2, 5, 10, 20, and 30. The distributions of a random configuration of particles, and of an isotropic network (isonet) are also plotted for comparison Here, we discuss the results obtained form the skeletonization of our structures. This technique uses a 3D-binary image as an input, in this image the cavities are filled in for avoiding complications The results concerning the filling of the cavities (Manikas et al. 2020), are presented in "Appendix D" for all φ, and B * = 0.01. The number density, and degree of BPs are presented in Fig. 7. The average thickness of branches for all φ used are presented in Fig. 8.
We expect to find an initial increase of BPs due to the initial structure formation and then a decreasing amount of BPs for increasing t * . The particles create a structure at early stages, and after they rearrange microstructurally to the optimal configuration with respect to the potential (more compact, less short branches). One can observe in Fig. 7, that for φ = 10% our expectation is verified, as after an initial increase until t * = 2, n BP decreases and eventually reaches a plateau. We also expect the plateau value of n BP to increase with increasing φ, as the smaller the inter-particle distances the more network-like structures created. Our expectation is confirmed in the final structure (t * = 10), the n BP is increasing with increasing φ.
As the structure is formed we expect more complex pathways being present initially, and the degree of BPs to be reduced as thicker components and simpler pathways are formed over time. In Fig. 7b, for φ = 10%, we observe an initial increase following the same trend as n BP , and after that, a plateau around its minimum value (three) is reached within errorbars. We also expect the degree to be unaffected for varying φ, as the structure at advanced times is evolved and the complex pathways have been replaced by thicker components of the structure. In Fig. 7b, no significant influence of φ is observed for the degree at the final structure (t * = 10) as expected.
The evolution of the structure is expected to be followed by increasing thickness in time. In Fig. 8, for φ = 10%, we observe an initial increase until the structure is formed (t * = 1), and then deviates around a plateau value. An exemption is observed for φ = 30%, the thickness is increasing until dimensionless time of approximately 4, and then decreases to approximately a stable value. This decrease is related with the filling of the cavities shown in "Appendix D", as the volume of cavities filled decreases an order of magnitude at t * > 4 the average distance of a skeleton voxel to an empty voxel is decreased as the percentage of empty voxels is larger in the box. Therefore, thinner branches are present. We also expect an increasing thickness for increasing φ, as the higher the φ the smaller the inter-particle distances and the thicker the structures formed. In Fig. 8, for t * = 10, one can observe that indeed the thickness is increasing proportionally to φ.
The existence of percolation is also checked. It is expected to achieve percolation along a single direction at the same time, as N * reaches a value of one. Percolation occurs for lower (1%, 2%, 5%, 10%) φ at time of 1 in the direction of the field z. Percolation in not reached in the other directions (x, y), in the limit we investigate in this paper. On the contrary, the higher φ (20%, 30%), reach percolation in z at t = 1 and in x, y at t = 2. This is explained by the fact that the unoccupied space in larger φ is limited and results to the creation of percolation paths also on the x y plane, which is perpendicular to the main direction of anisotropy.    Fig. 8 The average thickness of branches is plotted against time for φ = 1, 2, 5, 10, 20, and 30

Thermal conductivity
At this point, we would like to present Fig. 9, where the enhancement factor of the thermal conductivity of the medium for the principal directions for the φ used in this article is shown. We did not develop this theoretical analysis (Martin and Gulley 2009), but we use it to demonstrate the enhancement that one can achieve by aligning the particles  Fig. 9 Thermal conductivity enhancement factor (Martin and Snezhko 2013) against φ for final configuration in a medium. This analysis of the thermal conductivity of composites produced under an external field (magnetic or electric) shows that the particle structuring could improve thermal transport, however, other transport properties have isomorphic behavior (e.g. electrical conductivity, dielectric permittivity) (Martin and Snezhko 2013). This theory is an adaptation and extension of the Maxwell theory that takes into account spatial correlations between particles (Eucken 1932), it originates from the equivalent theory concerning the magnetic susceptibility (Martin et al. 2000). In Fig. 9, one can see how much the thermal conductivity coefficient K eff can be enhanced by the alignment of particles compared to random addition of particles and no addition of particles, K poly . The enhancement ratio is calculated as (Martin and Gulley 2009) where the subscript z refers to the direction of the external field, w is an index that indicates any of the three principal axes, x, y, or z, of the sample, δ x = δ y = 1, and δ z = 2. This expression differs from the Maxwell theory only through the structural order parameter: where P 2 (x) = 3 2 x 2 − 1 2 is the second Legendre polynomial. We use Eq. (25), developed in Martin and Gulley (2009), to relate the actual structure with the thermal conductivity, as an example of the enhancement of transport properties due to the structure formation. A more direct way of calculating the thermal conductivity would be performing full finite-element method (FEM) calculations, solving the energy balance (Guan et al. 2014). Another way of obtaining the thermal conductivity from simulations would involve the Fluctuation Dissipation Theorem which relates thermodynamic response functions to the appropriate timeautocorrelation functions. The thermal conductivity is related to the time-autocorrelation function of the heat-flux operator. The effective contribution of particles towards the overall thermal conductivity of the system is calculated using the appropriate Green-Kubo relation in relation to the heat-flux operator. Once this contribution to the thermal conductivity is obtained, the parallel model of thermal conductivity is used to calculate the effective conductivity as described in Bhattacharya et al. (2004).

Comparison to literature
In this paper, a relation between the physical input and the final structure of suspensions under a uniform external field is made. Our numerical results can be compared to experimental findings as follows. The formation of strings in the direction of the external field is observed in our simulations, and it has also been observed experimentally both for a few particles (Halsey 1993;Loudet and Poulin 2001) and many particles (Sprecher et al. 1987). Beyond such a general agreement, we would also like to compare experiments with our numerical results in the regime of volume fractions for which our model is most appropriate, φ ≤ 10%. In this regime, it is observed in experiments (at φ ∼ 2%, see Fermigier and Gast (1992)) that chains form, which then might percolate depending on the magnitude of the field and the time the suspension is exposed to this field. This behavior is analogous to the regime of structure formation t * ≤ 1 of our suspensions in Fig. 5. The thickening of the clusters observed in our simulations has also been observed experimentally (Mohapatra et al. 2020), where the chains aggregate over time and form thicker chains. Additionally, taking the experimental parameters studied in Mohapatra et al. (2020) and using them in the calculation of the characteristic timescale according to our procedure, one obtains indeed the times that correspond to the structure formation observed experimentally in Mohapatra et al. (2020). This supports our hypothesis that two-particle scaling is appropriate for describing many-particle dynamics.

Summary and discussion
The structure formation and evolution of electrically/magnetically (EM) active particles under a spatially uniform external electric of magnetic field was studied. EM active particles create dipoles when exposed to an external field. When two dipoles are sufficiently close, their interaction results in the formation of structures, e.g. chains. The physical parameters of relevance for this formation were organized in dimensionless groups, and the effect of these dimensionless numbers on the structure evolution was investigated. The possibility of 3D printing these systems with stereolithography (Stansbury and Idacavage 2016) was explored in terms of the relevant timescales. In this section, the results concerning the physical behavior of the system presented in Sect. 4, are discussed.
We used BD simulations to track and characterize structures of particles under a spatially uniform external field. We simulate the electro-magnetic interactions with the simple dipole approximation, the Stokes drag and Brownian force are used for the interaction with the surrounding medium. A method was presented to capture the structural dynamics of these systems, with special attention to the effect of the dipole-dipole interactions on the structure formation and evolution. We introduced a combination of measures for quantifying the evolution of the structure in time.
In Sect. 4.1, the initial configuration of the system was varied, and exhibited a larger impact on the results than the effect of the variation in random numbers sequence. The variation of the initial structure is required for the observation of the structural evolution of such systems. One should be careful with that observation, as an initial configuration with specific structure (e.g., FCC crystal structure) could have a large impact on a study of this kind.
The observations about B * -variation were the following. Our scaling concerned only two-particle interactions, so our parameters like B * and t c are based on that assumption. This scaling should be obeyed for simulations of two-particle systems. However, a generalization to many-particle systems cannot be performed easily, as it is not clear if the same scaling is obeyed. To this end, we performed simulations of 1000 particles with various values of B * , where we observed that the scaling is obeyed for many-particle systems as well, and the magnitude of the thermal fluctuations can be considered negligible in the regime we investigated, the only effect being that the fluctuations around the average increase.
Additionally, one can identify three main regimes, in Fig. 4. During the increase of S 2 in Fig. 4a, we have the formation of the chains. After that and until dimensionless time of 10, the anisotropy stays at the maximum level, Fig. 4b. The orientation measure S 2 reaches a maximum and then decreases to a plateau, as it is sensitive to thickening of the strings that are being formed. This conclusion is based on Fig. 4c, were at the same time the increase begins t * ∼ 0.1, the chains span the box in length, N * ∼ 1, after that the chains keep growing and the S 2 drops due to the thickening. At larger times, t * > 8, it seems that the system reaches an equilibrium state, where S 2 exhibits a plateau.
The φ-variation was performed in steps. The first step involved the same three measures we used before: S 2 , λ * , and N * . There, we identified that different φ-values exhibit different behavior, as our two-particle scaling is not enough to resolve this variation. However, one can see that with increasing φ, the S 2 and λ * decrease to really low values indicative of isotropy, at similar values of dimensionless time. From N * , we observe that at high volume fractions of particles, φ > 20%, all particles of the system belong to the same cluster. These big aggregates are not characterized well by the above measures. The reason is that these measures were developed for string-like structure and at this point we do not have distinct strings. The structures look more like a network. To characterize them, we focused on the investigation of neighbourhoods of particles using the Voronoi tessellation. Although, some qualitative observations were made, one could not distinguish particles that belong to the core of the structure by the characteristics of their Voronoi cells. The main reason is that too many particles belong to the layer surrounding the core layer, resulting to large fraction of particles with the characteristics we are seeking for.
At high volume fractions, attention was paid on transforming the structure into a form that is easier to characterize, namely a skeleton with simplified connectivity. In this way, we were able to calculate the branch points (BPs), the thickness of branches, and the existence of percolation in all three directions and study their evolution in time. After the formation of the structure at early times, the local rearrangement of particles results in more compact structures as time advances. The BPs and thickness of branches seem to be mainly affected by the volume fraction φ, the higher φ the higher the number density of BPs and the thickness of branches. The existence of percolation and the anisotropy of the structure were investigated. For φ > 10%, the structures formed percolate in all three directions (x, y, z). A lower values of φ, the structure percolates only in the direction of the external field, z.
The thermal conductivity is calculated by way of a theoretical analysis, and it appears to be effectively increased by the structuring of the particles (see Fig. 9). The relative enhancement ratio increases monotonously with the volume fraction for the final configurations of our structures. In addition to that, the thermal conductivity in all directions (x, y, z) is enhanced in relation to the random arrangement of particles. It has to be noted that this enhancement is significantly larger if polydispersity is introduced to the system, as the interparticle distances are decreasing and the effective thermal conductivity coefficient is growing, Eq. (25). In this paper, the thermal conductivity is presented, however, other effective properties such as the magnetic susceptibility, dielectric permittivity, and electrical conductivity exhibit isomorphic behaviour (Martin and Snezhko 2013).
Acknowledgements This research forms part of the research program of the Brightlands Materials Center (BMC, www.brightlandsmaterials center.com).
Funding KM has received funding by the Brightlands Materials Center (BMC, www.brightlandsmaterialscenter.com).

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.
The failure of the simulations occurs, because of really large forces produced by large overlaps between the particles. This forces cannot be resolved with longer timesteps than the one calculated by t * c,exv .

B Initial configurations
The protocol followed for the creation of equilibrated initial configuration will be presented here. Initially, we work with a system of 1000 particles in φ = 40%, so if the configuration is equilibrated for this system it will be also for every system with the same positions of the particles and lower volume fraction. We initially place the particles in a simple cubic lattice. Then, we let them equilibrate with a repulsive part of a Lennard-Jones potential. On top of that, we use a higher temperature T = 493K, so that the diffusion timescale reduces. The equilibration is checked with the standard deviation value of the number density of particles per cubic cell. We observe the standard deviation, which drops with time, Fig. 11a. Our criterion for terminating the equilibration loop is the standard deviation to reach a value lower than 0.2. For further checking this assumption we calculate the radial dis-tribution function, Fig. 11b, and one can observe that it looks like a Lennard-Jones fluid with excluded volume interaction (Frenkel and Smit 2002). As we obtained these random initial configurations, we studied the effect on the structure formation. The volume fraction of the selected system for study is φ = 10%. In  Fig. 12, we present the measure introduced in the previous section for four different initial configurations, indicated by different colours. One could observe that in one dimensionless unit of time the difference in the measured quantities is low compared to the values of the measures, but visible. If one compares with the variation of the random numbers (error-bars), the variation of the initial configurations has a larger impact on the measured quantities.

C Variation of B * and volume fraction
In this part, we present the results for different B * , and φ. These results correspond to the same values used in Sects. 4.2 and 4.3 for B * = 0.01, 0.0316, 0.1, 0.316. One can see that the results for different φ's presented collapse to a single master curve as in Fig. 4. The three measures S 2 , λ * , and N * are presented for φ = 1, 30% in Figs. 13 and 14 respectively. In Table 1, one can find all the values of B * and φ used in this paper and their relevant location in the paper.

D Filling of cavities
When an image is processed to its skeleton all closed cavities (cavities that do not connect to the side of the simulation box with empty voxels) are represented as surfaces of voxels, to avoid this artefact of the technique it is common to fill the cavities before processing (Kerschnitzki et al. 2013;Kollmannsberger et al. 2017). These cavities are characterized in terms of number density and volume fraction with respect to the simulation box, as they are being part of the structure. In Fig 15, one can see that the φ dominates both the number density (normalized with number density of primary structure), and volume fraction of cavities. The number density   . 15 a Number density of cavities for different volume fractions, φ = 1, 2, 5, 10, 20, and 30 is plotted against time. b Volume fraction of cavities is plotted against time for the same φ shows an initial increase, which is interpreted with the structure formation happening. The general trend is the decrease of both number density, and volume fraction of cavities are decreasing in a consistent way, until they reach a plateau and fluctuate around this plateau.