Structure evolution of suspensions under time-dependent electric or magnetic field

The effect of time-dependent external fields on the structures formed by particles with induced dipoles dispersed in a viscous fluid is investigated by means of Brownian Dynamics simulations. The physical effects accounted for are thermal fluctuations, dipole-dipole and excluded volume interactions. The emerging structures are characterised in terms of particle clusters (orientation, size, anisotropy and percolation) and network structure. The strength of the external field is increased in one direction and then kept constant for a certain amount of time, with the structure formation being influenced by the slope of the field-strength increase. This effect can be partially rationalized by inhomogeneous time re-scaling with respect to the field strength, however, the presence of thermal fluctuations makes the scaling at low field strength inappropriate. After the re-scaling, one can observe that the lower the slope of the field increase, the more network-like and the thicker the structure is. In the second part of the study the field is also rotated instantaneously by a certain angle, and the effect of this transition on the structure is studied. For small rotation angles (θ≤20∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta \le 20^{{\circ }}$$\end{document}) the clusters rotate but stay largely intact, while for large rotation angles (θ≥80∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta \ge 80^{{\circ }}$$\end{document}) the structure disintegrates and then reforms, due to the nature of the interactions (parallel dipoles with perpendicular inter-particle vector repel each other). For intermediate angles (20<θ<80∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$20<\theta <80^{{\circ }}$$\end{document}), it seems that, during rotation, the structure is altered towards a more network-like state, as a result of cluster fusion (larger clusters). The details provided in this paper concern an electric field, however, all results can be projected into the case of a magnetic field and paramagnetic particles.

permittivity or magnetic susceptibility, respectively, between them and the medium (Jones 1995). Due to the induced dipoles, the particles interact (Jones 1995) and, as a result a structure of particles emerges (Klingenberg et al. 1989). The particles, depending on the field conditions, can form strings (Martin et al. 1998) that develop by either thickening or aggregating (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). The properties of the entire composite material are enhanced especially in the direction along which the structures are created (Martin et al. 2000;Martin and Gulley 2009). Particularly, the formed structures greatly affect the transport properties, due to the existence of percolation paths (Martin and Snezhko 2013). Electro-/magneto-rheological fluids are suspensions that contain filler particles that are responsive to an external field; these suspensions have become a subject of research due to their unique characteristics. The characteris-tic time-scale of their rheological response is O(1s) and is exploited in robotics and more specifically in force-feedback sensors (Karasawa and Goddard 1989). A significant amount of studies have investigated the rapid viscosity increase in the direction perpendicular to the applied field (Klingenberg et al. 1989Brady 1990, 1992;Mohebi et al. 1996;Ginder and Davis 1994;Satoh et al. 1998). The structures created and their behavior in a external spatially uniform field have been studied on large time-scales both in athermal (without Brownian motion) (Martin et al. 1998), and thermal (with Brownian motion) (Martin 2001;Manikas et al. 2021) settings.
Additive manufacturing has been a subject of scientific interest for quite some time (Yan and Gu 1996), and it has received increased attention more recently (Ngo et al. 2018). Different kinds of materials are used including metals (Visser et al. 2015), ceramics (Warnke et al. 2010), polymers (Ligon et al. 2017) and combinations of them (Ngo et al. 2018). Polymer-matrix composite materials are printed and used for various purposes such as topology optimization (Maute et al. 2015), and rapid prototyping (Czyzewski et al. 2009). The range of applications includes aerospace applications (cabin interiors) (Raja et al. 2010), medicine (mimics of living tissue) (Villar et al. 2013), anthropology (reconstruction of medieval skulls) (Massimiliano et al. 2008), and design (spatially-dependent elasticity) (Oxman et al. 2012). Another advantage of composite materials is the functionalization of the matrix material by introducing dielectric, magnetic or conductive functionality by adding filler particles (Czyzewski et al. 2009;Castles et al. 2016;Kokkinis et al. 2015). Additionally, specific structures of particles could possibly be achieved inside a complex geometry by controlling the microstructure of the composite while printing, through the use of an external electric or magnetic field. This would be possible if photo-reactive resins (Anastasio et al. 2018) filled with electrically of magnetically active particles are used: Particles are dispersed in the resin (Kim et al. 2003;Kim and Shkel 2004), so that by solidification of the resin one could fixate the desired structure of particles; before curing, resins exhibit low viscosity allowing the structure formation within the time-scales (O(1s)) of stereolithography (SLA) (Bártolo 2011). Potential applications of these systems include flat lenses with a gradient of concentration of particles that have the functionality of their curved counterparts (Kurochkin et al. 2018), personalised hearing aids (Dodziuk 2016), piezoelectric or Hall effect sensors (Quanlu 2002) and direction-specific thermally or electrically conductive composites (Martin and Snezhko 2013). In the following, the formation of particle structures in a time-dependent imposed field is examined numerically.
A combination of electromagnetic suspensions and curable resins with application to SLA would be possible with sufficient knowledge of the structure formation under time-dependent fields. Different kinds of structures can be obtained by time-dependent fields, more network-like structures being obtained if the field is increasing gradually (Martin et al. 2000). The effect of this slope with respect to the characteristic time-scales of structure formation (Manikas et al. 2021) is investigated in this paper. Multidirectional external fields are interesting due to the transitions of the fields either inside a layer or from layer to layer during 3Dprinting, if one wants to create a path of particles inside a layer or a 3D structure of particles inside the printed object. In the literature, rotating fields has been studied before. The effect of the Mason number (Melle 2000), which is the ratio of viscous to magnetic forces (Gao et al. 2012;Melle et al. 2003;Calhoun et al. 2006), is studied in a system with a rotating field. In this paper, an already formed structure is subjected to a rotated, with respect to the initial field orientation, field and the structure evolution is followed.
In this paper, we use Brownian dynamics (BD) simulations to simulate the motion of particles (dipoles) in a fluid under a spatially uniform time-dependent field. This results in the formation of structures of particles which are analyzed in detail. As already mentioned, strings (Martin et al. 1998), planes (Martin et al. 2000), or networks (Martin and Snezhko 2013) of particles can be obtained depending on the field conditions. There are a lot of techniques used for the characterisation of an assembly of points/particles (Theodorou and Suter 1985;Voronoi 1908;Li and Li 2009;Varadan and Solomon 2003;Montoro and Abascal 1993;Greenfield and Theodorou 1993;Starr et al. 2002;Damasceno et al. 2012;Vogiatzis and Theodorou 2014;Hütter 2003Hütter , 2000. In previous studies, we achieved quantitative characterisation of structures of particles by using a combination of techniques (Manikas et al. 2021(Manikas et al. , 2020, namely skeletonization and cluster analysis. Here, these techniques are used for the quantitative characterisation of the structures formed under a time-dependent imposed field. The goal of this paper is twofold. In the first part, we investigate the effect of unidirectional time-dependent fields on the structure formation, and compare it with the case of a constant field (Manikas et al. 2021). In the second part, we investigate the effect of rotating the field on the structure evolution. The paper is organised as follows. Section 2 describes the tools that were used for the production and characterisation of the structures. In Sect. 3, results are presented for the evolution of the structure in the course of (dimensionless) time, and the characteristic features of the formed structures for both uniand multi-directional fields are discussed. Finally, the paper is concluded with a discussion in Sect. 4. system consists of N p spherical particles with radius R p at volume fraction φ, dispersed in a fluid with viscosity η, and contained in a cubic box of edge-length L; periodic boundary conditions are applied in all three directions. An external field is applied to the system. In the particles, dipoles are induced due to the mismatch in the dielectric permittivity of the matrix ( m ) and the particles ( p ), under a time-dependent but spatially homogeneous external field.

Simulation details
In our model, we include dipole-dipole interactions (Jones 1995;Klingenberg et al. 1989), excluded volume interactions (Klingenberg et al. 1989;Gao et al. 2012), single particle Stokes drag between the particle and the surrounding liquid, and thermal fluctuations. For reducing the physical and numerical parameters of our model to the minimum, we introduce characteristic scales for length, dipole moment, and time. The positions of the particles are scaled with the characteristic length-scale r c = 1/ 3 √ n p , where n p is the number density of particles (note: φ = (4π/3)R 3 p n p ), and the dipole moments with the characteristic dipole moment magnitude p c = 4π m K R 3 p E c , with E c the maximum value of the electric field and K the Clausius-Mossotti constant (Jones 1995). The time is scaled with the characteristic time-scale introduced in Manikas et al. (2021), Equation 1 is the result of writing the particle dynamics in dimensionless form with respect to the dipole-dipole interactions. The dipole-dipole interactions are considered the most dominant and interesting interactions with respect to the focus of this study, i.e., the effect of time-dependent fields. For values of r c that correspond to φ = 1−30%, m = 4 0 , and R p = 1¯m, one finds that the timescale t c is on the order of O(0.1−10s). In contrast, the Brownian timescale for the present system is given by t B = 4R 2 p /(6D 0 ), where D 0 = k B T /(6πηR p ); using R p = 1¯m, η = 0.5 Pa s, and T = 293 K, one finds t B = 1.56 × 10 −3 s, which is orders of magnitude smaller than the characteristic timescale t c .
Note that the particles under consideration are not equipped with rotational degrees of freedom for the following reasons: First, the particles are spherical, and second the dipole moment that dictates the particle-particle interaction is determined solely by the imposed field, i.e., the dipole moment of a particle does not change if the imposed field is constant even though the particle may rotate. For these reasons, in the present study rotational degrees of freedom have been omitted.
The corresponding dimensionless time is denoted by (2) To study the system numerically, we use BD (Öttinger 1996;Gardiner 2004) simulations of the following evolution equations, in dimensionless form (see Manikas et al. 2021 for the derivation), with While ξ i are statistically independent vectors with statistically independent components that are drawn independently from a Gaussian distribution with zero mean and a standard deviation of unity, F em * i is the force due to the dipole-dipole interactions, F exv * i is the force preventing the particles from overlapping, and p * i is the scaled dipole moment of particle i, with scaled electric field E * = E/E c ; p * i in general depends on t (1) , however, this dependence is omitted from Eq. (4) for convenience. The notation r * = r/r c is used to denote the position vectors r i and difference vectors r i j = r i − r j in rescaled form, andr is the unit vector,r = r/r , where r = |r|; this notation is used for all quantities in the rest of the paper for scaled vectors, unit vectors and magnitude of vectors. The range of excluded volume is chosen as κ = 30; further details can be found in Gao et al. (2012) and Manikas et al. (2021). The Boltzmann constant is denoted by k B , the absolute temperature by T , and the friction coefficient is ζ = 6πηR p . The long-range corrections of the dipoledipole interactions are resolved using the Ewald-summation technique (Wang and Holm 2001).
The only parameter that is varied in this study is the dimensionless quantity B * , which denotes the ratio of the thermal effects to the dipole-dipole interaction. It summarizes the two major physical effects affecting the structure formation in one parameter. Its value determines whether and when structure formation occurs, and it scales as B * ∼ √ T /E c . More details concerning the justification of the choices made for the model, physical effects and the validity of the interaction forces and algorithm can be found in previous work of the authors (Manikas et al. 2021).

Characterisation of morphology
In this section, we shortly introduce the measures that we are going to use for the rest of this paper, for more details on these measures the reader is referred to Manikas et al. (2021). Morphology measures are separated into cluster measures and structure/network measures. Figure 1 gives a schematic representation of certain quantities in relation to the measures introduced below. The first measure concerns the orientation of the inter-particle vectors in a cluster with respect to an external direction, that being the directionÊ of the imposed electric field, where the index I indicates that the average is taken over all particles in cluster I . Furthermore, the anisotropy of each cluster is quantified by with the gyration tensor (Theodorou and Suter 1985) (superscript "T" denotes the transpose, and r cm,I is the position of the center of mass of cluster I ) and λ max,I the largest eigenvalue of S I ; λ * I = 0 for isotropic and 0 < λ * I ≤ 1 for anisotropic clusters. The number of particles in a cluster, N I , is reported in the form where N L is the minimum amount of particles that are needed to span the box-length L in one direction. The number of particles per cluster is related to the number of clusters: Since the number of primary particles in the system is constant, an increase in the average number of particles per cluster (i.e., increasing N * ) implies that the number of clusters decreases. Throughout this paper, the averages of the above quantities over all clusters will be reported, weighted by N I , Fig. 1 Schematic representation of the minimum number of particles needed to span the box, N L , the misalignment of the particle connector r i j with the imposed field E, and the ellipsoid of the gyration tensor S with A the quantity of interest and N cl the number of clusters. Furthermore, the spatial arrangement of particles is characterized with a view on the emerging network of all particles, without subdividing the set of all particles into clusters. To that end, we have recently developed a method based on the image analysis technique called skeletonization (Kollmannsberger et al. 2017) in combination with graph analysis tools (Manikas et al. 2020). The simplified skeleton is characterised in terms of (1) the number density of branch-points (BP), see also Manikas et al. (2021Manikas et al. ( , 2020, with N BP the number of BP in the simulation box, (2) the degree of the BP averaged over all BP, d BP , which is defined as the average amount of bonds of a BP, (3) the average thickness of a branch d B (for details, see Manikas et al. 2020), and (4) the existence of percolation.

Results
In this section, we present simulation results concerning suspensions under uni-directional time-dependent fields, and multi-directional time-dependent fields. The system simulated consists of N p = 1000 particles, averaging over three simulations with different random numbers sequences, and three different random initial configurations to generate the error-bars (Manikas et al. 2021). In all simulations, the ratio between the thermal fluctuations and the (final) strength of the electric field, B * , and the characteristic time-scale t c are kept constant; only the time-dependence of the field-strength of the electric field is varied (i.e., we vary the increase of the field strength until the final (constant) value E c is reached) through the simulations. This means that the magnitude of the thermal fluctuations relative to the final field strength

Uni-directional time-dependent fields
For uni-directional time-dependent fields, we use a spatially constant field with variation of its strength over time.
In Fig. 3, three snapshots of the configuration of particles during early (a), intermediate (b) and late stages of structure formation in the uni-directional field are provided for illustration purposes. More quantitatively, in Fig. 4, one can observe the measures S 2 , λ * , and N * for t (1) s = 0.0, 0.1, 1.0, and 10.0. As one can observe, the initial response of the measures presented depends strongly on t (1) s , so the constant t (1) -scaling is not optimal for this case. We would like to scale the time corresponding to the field value at each point in time. In this way, one can compare the structure evolution in terms of equal strength of the external field and identify the differences induced by the different magnitude of the thermal fluctuations in time. In view of first contribution on the right-hand side of the dynamics (3) with scaled dipole moment (7), this can be achieved by introducing where for a field of the type presented in Fig. 2, one finds, Departing from Eq. (3), one finds with where the dipole momentsp i are unit vectors. The deterministic term in Eq. (19) is of order unity and does not show the explicit time-dependence of the electric field anymore, however, the stochastic term now depends on the actual value of the field strength E * , see Eq. (22). That means that if there were no thermal fluctuations, the cases for different t (1) s were indistinguishable in the t (2) -representation, the latter representing the case of a time-independent field. However, this is not the case in this paper, as thermal fluctuations are included in our model. Thermal fluctuations are relevant particularly at early times, where the field strength is small and B * * is large (see Eq. (22)). The t (2) -scaling is only accurate if the dipole interactions are dominant over the thermal fluctuations (B * * close to 0). B * * ranges from zero, where the thermal fluctuations are absent and our t (2) -scaling is sufficient, to a large value, where the thermal fluctuations dominate over the dipole-dipole interactions. In our case, we start with a large value for B * * , which then decreases gradually to a constant value. One needs to define a transition value, below which the B * * -value is low enough for the field to be dominant over the thermal fluctuations, so that the formation of clusters occurs at the same dimensionless time. To that end, we shift the curves such that we "neglect" the time before the strength of the field reaches E * negl = 0.4933. If one substitutes this value in Eq. (22), then a value of B * * = 0.0203 is obtained, meaning that the prefactor of the dipole-dipole interactions is almost 50 times larger than the one of the thermal fluctuations. The shift is the third time-transformation we apply and present it as t (3) , where c negl = E * negl 3 /3 = 0.04, and originates from the value of t (2) when E * = E * negl . For the sake of simplicity, we will refer to the different conditions of electric field in terms of the value of t (1) s (i.e., in scaling level t (1) rather than in scaling level t (3) ).
In Fig. 5, one can see the results presented in Fig. 4 after the transformation of time to t (3) . The initial increase of the cluster measures occurs at the same time, meaning that our scaling works as expected. This means that the evolution of the structure due to the dipole-dipole interactions is actually suppressed until the field is strong enough to dominate over the thermal fluctuations. However, our t (3) -scaling is still not ideal, as the transition between the two regimes described (large versus zero B * * -value) occurs gradually with a finite t The skeletonization measures are also examined for the same set of simulations. We expect that with higher t (1) s the structure becomes more network-like rather than string-like, because if the thermal fluctuations are stronger relative to the dipole-dipole interactions, the particles are arranging in a more isotropic way. Our expectation is verified by the higher number density of BPs achieved with higher t (1) s , see Fig. 6a. The degree of BPs is approximately the same for all t (1) s , i.e., no significant difference in the complexity of the network structure is observed, Fig. 6b. By complexity of the network we refer to a structure becoming more complex if the number and/or degree of the branchpoints increase. However, the standard deviation of the degree is clearly larger for the timedependent cases, except the case of t (1) s = 1.0. The thickness is also influenced by t (1) s : the higher t (1) s , the higher the thickness, see Fig. 6c. It is also observed that larger variations in thickness occur, in comparison to the constant-field case (Manikas et al. 2021). The observed overshoot and decrease in thickness for t  1, and 10. b Degree of branch-points versus time for the same φ and slopes. c The average thickness of branches versus time the course of time, these cavities disappear (defects heal) due to relaxation and local energy-minimization, which in turn implies that the branches densify and thereby their thickness decreases. In fact, the volume of cavities drops an order of magnitude (for t (1) s = 0, Manikas et al. 2021) due to the rearrangements of particles after the structure formation occurs (t (3) > 2). In contrast, for higher values of t (1) s , the gradual increase of the strength of the external field allows for a better rearrangement of the particles already during the early stages of structure formation, and therefore we do not see the overshoot and decrease in thickness at the higher values of t (1) s .
In Fig. 7, results are shown for t (1) s = 1.0, and different values of the volume fraction φ (1%, 2%, 5%, 10%, and 20%). The structure formation is expected at the same dimensionless time, since a major effect of φ has already been taken into account when re-scaling towards t (1) . In Fig. 7, one can observe that independent of φ an initial response in the structure measures occurs at t (3) ∼ 0.1, however, different behavior between the φ-values is observed thereafter. The discrepancy observed for varying φ implies that the twoparticle scaling, on which Eq. (1) and (2) are based, is not sufficient for collapsing the structural evolution of manybody systems for different φ (Manikas et al. 2021). That means that more complex processes take place than what could be anticipated based on a two-particle picture. A significant drop in S 2 and λ * is observed as time increases, especially for φ > 10%, which implies the formation of isotropic structures. The measures presented here are calculated as cluster averages, see Eq. (12). In the presence of clusters comparable to the total amount of particles of the system for φ > 10%, no orientation/anisotropy can be detected as the measures account for the whole system at once.
With respect to the network analysis, the effect of the volume fraction φ can be described as follows. The number of BPs is expected to increase initially due to the emerging structure formation followed by a decrease of the number of BPs, as t (3) grows. A structure is created early in the process and the evolution occurs through microstructural rearrangements towards the optimal configuration prescribed by the interaction (more compact, less short branches). In Fig. 8a, for φ ≤ 20% after an initial increase, n BP decreases and settles to a plateau, meaning our expectation is verified. For φ = 30%, the number of BPs rises gradually over time; thicker and less mobile branches are formed due to the lower interparticle distances, resulting in forming connections between different branches in larger time. The final value of n BP is expected to raise with raising φ, since more network-like structures are created when smaller inter-particle distances are present. The final structure exhibits higher n BP for higher φ, as expected.
Initially, more complex structures are expected to appear as the structure is formed. As the structure evolves, thicker components and simpler structures are formed and the degree of BPs d BP is expected to decrease again. In Fig. 8b, for φ = 5% the same trend as n BP is observed, meaning that an initial increase followed by a plateau around a value of 4 is observed. No significant influence of the volume fraction φ on the final value of the degree d BP (at t (3) ∼ 9) is observed.
The average branch thickness d B is expected to increase during the initial formation of the structure, and it is also expected that the final value for the thickness would increase monotonically with volume fraction φ, due to the lower interparticle distances. In Fig. 8c, for all φ an initial increase is observed until the structure is formed, followed by fluctuations around a plateau value, as expected. It seems that for φ = 1%, 2%, 5% the time needed for reaching the plateau is shorter than for φ = 10%, 20%, 30%. The thickening is related to microstructural rearrangements leading to thicker components due to the elimination of small gaps between the particles, which is more difficult to occur when the particles are essentially immobilized. In Fig. 8c, for t (3) ∼ 9 one can observe that indeed the final value for the thickness d B increases monotonically with increasing volume fraction φ.

Transformation back to real time
In Sect. 3.1, the results are presented in terms of transformed time, t (3) . However, one could translate the results to a "real" (dimensionless) time, t (1) . The relation between t (1) and t (3) is given by which are Eqs. (23) and (24) inverted. Fig. 9 shows examples of this relation. The main difference between the t (1) -and t (3) -representations is the amount of time the system spends at B * * -values higher than B * * = 0.01. Our transformation reduces this discrepancy by scaling each moment in time with the time scale corresponding to the instantaneous value of the field strength. In this sense, the time it takes for reaching the maximum fieldstrength is basically compressed, see Figs. 4a and 5a. If one transforms back to t (1) , the higher the t (1) s the more time it takes the structure to form.

Multi-directional time-dependent fields
For multi-directional time-dependent fields, we depart from a configuration obtained by applying a field that is constant in time and space for a duration of t (1) = 10.0, then rotate the field instantaneously to a different direction (see Fig. 10) and follow the structure evolution for a duration of t (1) = 5.0. The reason that we use t (1) = 5.0 is that we are interested in the main transition, rather than the evolution after that; a time of approximately 5.0 should suffice for providing this information, as the characteristic time for two particles to touch is 1.0 and the structure formation occurs at the same time (Manikas et al. 2021). We use θ to denote the angle of rotation of the field around the x-axis; the initial field orientation is in the z-direction. The values of θ examined here range from θ = 10 • to θ = 90 • in steps of 10 • . In this part, the variation of the volume fraction is limited to φ = 1%, 2%, 5% and 10%. The choice for lower volume fraction is made based on the expectation that rotation at higher volume fraction will be significantly less effective due to the lack of space. Note that the S 2 reported in this section always uses as reference direction the field-orientation before the rotation, which is the z-direction.
The expectation about the physical evolution of the system is based on the interaction potential. The radial component of the interaction force changes its sign at an angle of θ c = 54.7 • between the external field and the interparticle vector. Concretely, starting from Eq. (4) and using that the dipoles are always oriented in the direction of the imposed field, one can show thatr i j · F em * i j = 1 − 3(cos θ) 2 /r * i j 4 , which changes sign at θ = θ c . This implies that two particles initially aligned with the external field would rotate trying to align their interparticle vector with the new orientation of the external field if θ < θ c , and would repel each other and reform their contact if θ > θ c . Therefore, we expect to have rotation of structural entities until the described angle, while disintegration and reformation in the new orientation are expected for larger angles.
The effect of the angle of rotation of the field orientation on the structure evolution is investigated for a system of φ = 5%, see Fig. 11. Before proceeding with the interpretation of the results for the orientation S 2 , a thought experiment is useful: If the structure before the reorientation consisted of perfect chains in the initial direction of the imposed field and if all these chains followed perfectly the rotation of the field, see Fig. 12, the value for S 2 of the reoriented chains with respect to the original direction of the field would be given byS 2 (θ ) = (3(cos θ) 2 − 1)/2, in analogy to Eq. (8). In Fig. 11a, for each angle of field rotation, the value on the curve which is closest to the corresponding value ofS 2 is included (triangle). For low angles of rotation, θ < 20 • , we observe a constant plateau initially for S 2 in Fig. 11a, followed by a decrease indicating thickening of the clusters. The anisotropy (Fig. 11b) remains at high levels throughout the process, a fact that indicates the presence of anisotropic clusters (chains). In Fig. 11c, one can observe (for θ < 20 • ) that the size of the clusters does not decrease, indicative of the clusters being rotated. For large angles of rotation, θ > 80 • , we observe that the orientation reaches the value of S 2 (90 • ) and then reaches a plateau. We interpret this as disintegration of the clusters/strings to smaller string-shaped clusters that rotate and then merge to form larger string clusters in the new field-orientation. Our explanation is supported by the results for the anisotropy, see Fig. 11b, where a drop is observed at intermediate times, indicative of smaller and less anisotropic clusters (λ * ∼ 0.6). For intermediate angles of rotation, 20 • ≤ θ ≤ 80 • , we observe scenarios intermediate to the ones just discussed, whereS 2 (θ ) is reached, but not maintained as the clusters are increasing in size. The anisotropy, λ * , remains in high levels (λ > 0.8) for these angles, thus the anisotropic shape is maintained. For the same range of θ , we also observe that the size decreases initially and then increases. The decrease depends on θ in a monotonic way: higher θ means larger decrease in size. For θ > 60 • substantial decomposition occurs, see Fig. 11c. After the decrease, the size increases again, however, the structure is formed in another direction, see Fig. 11a. The increase in size is greatly Fig. 12 Perfect rotation of two chains following the rotation of the external field for a (small) angle of rotation θ influenced by the rotation as merging of rotating clusters occurs for 20 • ≤ θ ≤ 80 • , see Fig. 11. The effect of the angle of rotation θ , for different volume fractions φ, on the structure after the transient regime is shown in terms of percolation in Table 1. For φ ≤ 10%, we expect percolation to occur only in the direction of the field before the switch, and then eventually to occur instead along the final direction of the field. For φ = 1%, 2%, 5%, 10%, our expectation is met for the 10 • , 20 • , 80 • and 90 • of rotation, as we have percolation in z direction for 10 • and 20 • , and percolation in y direction for 80 o and 90 • . No percolation in any direction is observed for all other angles for φ = 1%, 2%, as the clusters formed are not large enough to span the box in any direction. For φ = 5%, percolation is observed in both y and z directions for intermediate angles of rotation 30 • , 40 • , 50 • , 60 • and 70 • in the final state. For φ = 10%, percolation occurs in all three directions (x, y, z) for θ = 30 • , 60 • , which is related to the large clusters formed at these angles. For θ = 40 • , 50 • , 70 • , percolation is exhibited in the z and y directions, but not in the x direction. It is noted that the behavior for φ = 10% is non-monotonous; this is a signature of the fact that percolation is a binary measure that strongly depends on local rearrangements/paths.
The effect of the rotation of the field orientation on the structure formation is studied for low values of φ; the final values (after t (1) = 5.0) of the morphology measures for different φ are plotted versus the field orientation θ in Fig. 13. We expect that the higher φ, the less space is available for the clusters to move, and therefore cluster merging will result in more network-like structures, while for low φ ≤ 2% we rather expect rotation of the clusters. The orientation S 2 of φ = 1% and φ = 2% exhibit almost perfect alignment with the rotated field (i.e., S 2 being close toS 2 ), see Fig. 13. The anisotropy λ * maintains a high value and also the size N * remains almost constant, meaning that our expectation of rotating chains is met for low values of φ. The results for 5% and 10% show similar behavior at low and high angles. However, at intermediate angles (20 • ≤ θ ≤ 80 • ), S 2 seems independent of the rotation angle, showing a plateau close to random orientation, S 2 ∼ 0.25, see Fig. 13a. This effect is related to the fusion of clusters, since larger clusters result in isotropic S 2 -values (Manikas et al. 2021). The anisotropy remains high for all φ, except at certain angles for φ = 10%. The data for φ = 10% show two dips at 30 o and 60 • , that are related to fusion of clusters resulting in larger clusters. Correspondingly, one observes in N * two maxima at the same angles, Fig 13c. In general, N * increases monotonously with φ.
As far as the skeletonization measures are concerned, we expect simpler structures when rotation or disintegration and reformation occurs, and more complex ones when partial disintegration and cluster fusion occurs. In Fig. 14a, one can observe that the number of BPs depends only weakly on θ and, with the exception of φ = 10%, increases with increasing φ. For φ = 2% and θ = 30 • , there is a higher value, which may be due to the number of BPs strongly depending on local fusion of clusters. The lower BP-values of φ = 10% in comparison to φ = 5% can be explained on the basis of the resulting structure as follows (see Fig. 15): after rotation, the sample at φ = 5% forms interconnected chains, while the sample at φ = 10% forms sheet-like structures (Martin et al. 2000), resulting in less BPs. The complexity of the BPs, see Fig. 14b, does not exhibit any significant systematic trend over θ and φ, which indicates that the structures created are highly dependent on incidental local cluster fusion; the error bars are relatively large. In Fig. 14c, one can observe the average thickness of the branches over time, where the thickness seems almost independent of θ , within error bars. The thickness increases monotonously with increasing φ, which is expected since less space is available at higher φ. It is noted that the standard deviations of the thickness are large as compared to the average.

Summary and outlook
The structure formation and evolution of electrically/ magnetically (EM) active particles under a unidirectional or a rotated time-dependent external electric or magnetic field  Final configurations for φ = 5% (a) and φ = 10% (b), after rotating the field by θ = 20 • was investigated. Dipoles are induced by an external field to the EM-active particles. The interaction between the dipoles leads to rearrangements and eventually results in the formation of structures, e.g. chains. The main physical parameters that govern this formation were re-organized in dimensionless groups with respect to the external field, and the effect of these dimensionless numbers on the structure evolution was investigated (see Manikas et al. 2021 for the case of a constant external field).
Brownian dynamics (BD) simulations have been used to track and characterize structures of particles under a unidirectional time-dependent or a rotated external field. The dipole-dipole interactions are calculated using the simple dipole approximation, and the Stokes drag and Brownian force are used for the interaction with the surrounding medium (Manikas et al. 2021). The structural dynamics of these systems were captured, and special attention was given to the effect of the dipole-dipole interactions on the formation and evolution of the structure. A combination of measures introduced earlier (Manikas et al. 2021) (cluster properties, network analysis) was applied for assessing in a quantitative way the evolution of the structure in time under different conditions.
The effect of using a unidirectional time-dependent field on the structure formation was investigated. A spatially constant field with increasing field strength is used. The rate of increase in the field strength is determined by the time t (1) s it takes to reach the maximum field strength, after which the field is kept constant at its maximum value for t (1) = 10.0. The values for t (1) s examined are 0.1, 1.0, and 10.0. The dynamics are re-mapped nonlinearly in time in such a way as to scale out the time-dependence on the fieldstrength, however, this transformation is not ideal (lack of complete overlap of results for different conditions) due to the presence of thermal fluctuations. Our mapping would work perfectly in the absence of thermal fluctuations, or, more generally, if the dipole interaction dominates over the thermal fluctuations. To make this more apparent in the transformed time, the transformed time is also shifted in order to match the instances when the field becomes dominant over the thermal fluctuations, B * * < 0.0203 (dipole-dipole interactions 50 times larger than thermal fluctuations). It is noted that this additional shift would not be necessary if the increase in field strength from the initial to the final value occurred instantaneously, however, in this paper a gradual increase in field strength is of interest.
The results after the time-transformation exhibit a delay in structure formation for larger t (1) s . For larger t (1) s , more isotropic configurations of particles are created, namely network-like structures with thicker branches, and one also observes a larger standard deviation of the local structure in comparison to t (1) s = 0 (Manikas et al. 2021). The variation of the volume fraction φ revealed that after the formation of the structure at early times, the rearrangement of particles locally results in more compact structures as time advances, especially for higher values of φ.
An already formed configuration is used and the field is rotated by different angles θ , keeping the field constant thereafter for a time of 5.0. This imposed rotation results in either rotation of the formed substructures, e.g. chains at low angles (θ < 20 • ), or decomposition and reformation at large angles (θ > 80 • ). At intermediate angles, we observe cluster merging and partial decomposition and reformation of the formed strings. The effect of varying the volume fraction of particles φ was also studied. The φ-study showed that for φ < 5% rotation of the structures as a whole occurs, while for φ ≥ 5% less space is available for the clusters to move, therefore cluster merging is observed which results in more network-like structures.
The effect of a time-dependent multi-directional field was examined in this paper. The structure evolution strongly depends on the angle of rotation, however, it is possible to take advantage of an already formed structure for creating a new structure in a different direction. This possibility may benefit the potential 3D-printing practitioner aiming for different orientation from layer to layer or inside a layer. Structure characteristics that were not observed in unidirectional field (Manikas et al. 2021) can be generated, as rotating an already formed structure could result in more network-like structures for low volume fractions of φ ≤ 10%, or percolation in multiple directions can be achieved. Additionally, the structure can be rotated preserving the desired characteristics if the rotation angle is small enough.
In this study, the effect of time-dependent fields on the structure formation/evolution was assessed with a simplified and computationally efficient model. However, one could wonder about the influence of additional effects, e.g., many-particle hydrodynamic interactions, multipoles, or different geometries (i.e. wall barriers). All these effects would improve the model to simulate the physical processes more accurately and realistically; however, it is not expected that this would result in qualitative differences as compared to the results presented in this paper. Another direction for possible extension in the future could be to consider other experimental protocols, e.g., looking at other time-dependencies of the field, incl. field-strength decrease and double rotation. Such studies could enrich the toolbox for the development of novel processes in 3D printing, with control of the microstructure.
Acknowledgements KM gratefully acknowledges fruitful discussions and technical support by Georgios Vogiatzis. This research forms part of the research program of the Brightlands Materials Center (BMC, www.brightlandsmaterialscenter.com).
Funding KM has received funding by the Brightlands Materials Center (BMC, www.brightlandsmaterialscenter.com).

Declarations
Compliance with ethical standards Yes.

Conflict of interest
On behalf of all authors, the corresponding author states that there is 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/.