Self-assembly and rheology of dipolar colloids in simple shear - studied by multi-particle collision dynamics

TL;DR

This study uses advanced simulations to explore how dipolar colloids self-assemble under shear and how their microstructure influences rheological behavior, revealing universal structural dynamics across different systems.

Contribution

It introduces a hybrid simulation approach combining molecular dynamics and multi-particle collision dynamics to quantify the microstructure-rheology relationship in dipolar colloids.

Findings

Viscosity varies over several orders of magnitude.

Universal structural behavior governed by bonding and erosive stresses.

Quantitative agreement with literature data.

Abstract

Magnetic nanoparticles in a colloidal solution self-assemble in various aligned structures, which has a profound influence on the flow behavior. However, the precise role of the microstructure in the development of the rheological response has not been reliably quantified. We investigate the self-assembly of dipolar colloids in simple shear using hybrid molecular dynamics and multi-particle collision dynamics simulations with explicit coarse-grained hydrodynamics; conduct simulated rheometric studies and apply micromechanical models to produce master curves, showing evidence of the universality of the structural behavior governed by the competition of the bonding (dipolar) and erosive (thermal and/or hydrodynamic) stresses. The simulations display viscosity changes across several orders of magnitude in fair quantitative agreement with various literature sources, substantiating the universality of the approach, which seems to apply generally across vastly different length scales and a broad range of physical systems.

Figures

Tables

Table 1. Table 1: A summary of some of the previous rheological investigations of magnetic colloids (ferrofluids) from various literature sources: reference, designation, type of magnetic material and protective shell, estimated diameter of the magnetic core d m subscript 𝑑 𝑚 d_{m} of the interacting fraction, estimated particle hydrodynamic diameter D h subscript 𝐷 ℎ D_{h} (distance of closest approach σ 𝜎 \sigma ) including the stabilizing layer, type of solvent, estimated viscosity η 0 subscript 𝜂 0 \eta_{0} of the pure solvent, volume fraction of the dispersed magnetic material ϕ m subscript italic-ϕ 𝑚 \phi_{m} (total and estimating the amount of the interacting particles λ ≥ 1 𝜆 1 \lambda\geq 1 ), investigated range of shear rates γ ˙ ˙ 𝛾 \dot{\gamma} , highest investigated field B 𝐵 B , estimated dipole interaction parameter λ 𝜆 \lambda

		Dispersed phase				Medium		${\varphi }_m$, vol%
Ref.	Sample	Material	Shell	$d_m$, nm	$D_h$, nm	Solvent	$\eta$, mPa	Total	$\lambda \ge 1$	$\dot{\gamma }$, s-1	B, mT	$\lambda$
102	APG513A	Fe3O4	unknown	$14$	$20$	ester	$120$	$6.7$	$0.6$	$1-5$	$10$	$1.3$
109	APG513A	Fe3O4	unknown	$14$	$20$	ester	$120$	$7.2$	$0.6$	$15-37$	$40$	$1.3$
96	APG513A	Fe3O4	unknown	$16$	$20$	ester	$120$	$7.2$	$0.7$	$0.05-2$	$30$	$2.7$
110, 111	APG513A	Fe3O4	unknown	$17$	$21$	ester	$120$	$7.0$	$0.5$	$0.1-1$	$40$	$3.2$
111		Cobalt	Al2O3+KorantinSH	$10$	$16$	L9 oil	$65$	$0.4$	$0.4$	$0.05-1$	$40$	$3.5$
112	APG513A	Fe3O4	unknown	$16$	$20$	ester	$120$	$7.2$	$0.7$	$0.05-2$	$50$	$2.7$
112	APG513A	Fe3O4	unknown	$16$	$20$	ester	$120$	$7.2$	$0.7$	$0.05-2$	$50$	$3.2$
113	APG513A	Fe3O4	unknown	$16$	$20$	ester	$120$	$7.2$	$0.7$	$0.05-2$	$150$	$2.7$
113	APG513A	Fe3O4	unknown	$16$	$20$	ester	$120$	$7.2$	$0.7$	$0.05-2$	$150$	$2.5$
113, 114	Co87	Cobalt	Al2O3+KorantinSH	$10$	$16$	kerosene	$2.5$	$0.35$	$0.35$	$1-400$	$130$	$3.4$
103	Co87	Cobalt	Al2O3+KorantinSH	$13$	$19$	kerosene	$2.5$	$0.21$	$0.16$	$1-8$	$30$	$11$
104	S2	Cobalt	unknown	$12$	$18$	kerosene	$2.5$	$0.13$	$0.11$	$0.5-8$	$40$	$8.4$
115	NS FF	CoNi	L-$\alpha$-PC	$24$	$30$	mineral oil	$39.6$	$5.1$	$5.1$	${10}^{-5}-175$	$440$	$29$
116	FF $\varphi$=0.038	Fe3O4	oleic acid	$19$	$23$	kerosene	$2.5$	$3.8$	$2.2$	${10}^{-6}-{10}^3$	$200$	$6.7$
107	FF-A	unknown	unknown	$16$	$20$	unknown	$9$	$8$	$4.2$	${10}^2-{10}^3$	$280$	$1.7$
107	FF-B	unknown	unknown	$16$	$20$	unknown	$6$	$7$	$7$	$1-{10}^3$	$280$	$2.5$
108	FF	Fe3O4	oleic acid	$19$	$23$	kerosene	$2.5$	$8.8$	$4.1$	${10}^{-4}-{10}^3$	$340$	$6.4$
100	A0	CoFe2O4	citric acid	$35$	$40$	water	$1$	$0.23$	$0.23$	$3-27$	$40$	$28$
100	A0.75	CoFe2O4	citric acid	$28$	$33$	water	$1$	$0.91$	$0.39$	$3-27$	$40$	$13$
105		Fe3O4	unknown	$17$	$21$	ester	$34$	$6.3$	$3.1$	$0.5-10$	$50$	$3.2$
106	Co87	Cobalt	Al2O3+KorantinSH	$16$	$22$	kerosene	$2$	$0.72$	$0.66$	$0.5-10$	$60$	$22$
117	MagSilica	Fe3O4/Fe2O3	SiO2+C8TMS	$18$	$26$	AP201 oil	$68$	$1.3$	$1.3$	$0.1-100$	$40$	$2.5$
118	Sample 5	Fe3O4	OA/SA/ASA	$19$	$23$	PES-5	$250$	$3.6$	$1.6$	$1-100$	$130$	$5.1$

Equations

$$\lambda =\frac{1}{T^{\ast }}=\frac{{\mu }_0}{4\pi }\frac{{\mu }^2}{{\sigma }^{3}kT}$$

$$Mn=\frac{\pi {\eta }_0{\sigma }^3}{2\lambda kT}\dot{\gamma }$$

$$\beta U_{dd}\left({\mathbf{\mu }}_{\mathbf{i}},{\mathbf{\mu }}_{\mathbf{j}},{\mathbf{r}}_{\mathbf{ij}}\right)=-\lambda \left[\frac{3\left({\mathbf{\mu }}_{\mathbf{i}}\cdot {\mathbf{r}}_{\mathbf{ij}}\right)\left({\mathbf{\mu }}_{\mathbf{j}}\cdot {\mathbf{r}}_{\mathbf{ij}}\right)}{r_{ij}^5}-\frac{{\mathbf{\mu }}_{\mathbf{i}}\cdot {\mathbf{\mu }}_{\mathbf{j}}}{r_{ij}^3}\right]$$

$$\beta U_{WCA}\left(r_{ij}\right)=4\left[{\left(r_{ij}\right)}^{-2n}-{\left(r_{ij}\right)}^{-n}+\frac{1}{4}\right],r_{ij}\le \sqrt[n]{2}$$

$${\mathbf{r}}_i\left(t+\Delta t_{mpcd}\right)={\mathbf{r}}_i\left(t\right)+\Delta t_{mpcd}{\mathbf{v}}_i\left(t\right)$$

$${\mathbf{v}}_i\to {⟨\mathbf{v}⟩}_J+{\kappa }_J{\mathbf{\Omega }}_J\left(\alpha \right)\cdot \left({\mathbf{v}}_i-{⟨\mathbf{v}⟩}_J\right)$$

$$\sum _{n=1}^{\infty }{\int }_V{P}_n\left(\mathbf{r}\right)dV=1$$

$$g_n=\frac{\varphi }{n{v}_p}{\int }_V{P}_n\left(\mathbf{r}\right)dV$$

$$F=kT\sum _{n=1}^{\infty }{g}_n\ln \frac{g_n{v}_p}{e}+g_n{f}_n$$

$$g_n=\frac{X^n}{v_p}{e}^{-{ϵ}_{\ast }}$$

$$\sum _{n=1}^{\infty }n{g}_n=\frac{\varphi }{v_p}$$

$$X=\frac{1+2\varphi \exp {ϵ}_{\ast }-\sqrt{1+4\varphi \exp {ϵ}_{\ast }}}{2\varphi \exp {ϵ}_{\ast }}$$

$$⟨n⟩=\frac{{\sum }_{n=1}^{\infty }n{g}_n}{{\sum }_{n=1}^{\infty }{g}_n}=\frac{1}{2}+\sqrt{\frac{1}{4}+\varphi Z_2}$$

$$Z_2\left(\infty \right)\approx \frac{e^{2\lambda }}{3{\lambda }^2},Z_2\left(0\right)\approx \frac{1}{\lambda }{Z}_2\left(\infty \right)$$

$${\alpha }_{\lambda }=-\sqrt{\varphi }\ln X,{\beta }_{\lambda }=e^{{ϵ}_{\ast }}$$

$${\mathbf{\Gamma }}_{\mathbf{n}}^{\mathbf{m}}+{\mathbf{\Gamma }}_{\mathbf{n}}^{\mathbf{h}}=0$$

$$\tan {\theta }_n=\frac{1}{12}\left(n^2+n\right)Mn$$

$$n_{crit}\approx \sqrt{\frac{6\sqrt{2}}{Mn}}$$

$$G_{ij}^n=\frac{\varphi }{n{v}_p{g}_n}{\int }_V{r}_i{r}_j{P}_n\left(\mathbf{r}\right)dV$$

$$\tan \left(2{\theta }_n\right)=\frac{2G_{xy}^n}{G_{yy}^n-G_{xx}^n}$$

$$⟨\theta ⟩=\frac{{\sum }_{n=2}^{\infty }{\theta }_n{g}_n}{{\sum }_{n=2}^{\infty }{g}_n}$$

$$⟨\theta ⟩=sMn$$

$$s_0={\frac{\partial ⟨\theta ⟩}{\partial Mn}|}_{Mn\to 0}\approx \frac{1}{6}\varphi e^{2\lambda }\frac{180^{\circ }}{\pi }$$

$${\sigma }_{ij}=2{\eta }_0{\gamma }_{ij}+{\sigma }_{ij}^s+{\sigma }_{ij}^a$$

$${\sigma }_{xy}^a=\frac{1}{2}\sum _{n=2}^{n_{crit}}\left({\mathbf{\Gamma }}_{\mathbf{n}}^{\mathbf{m}}\cdot {\mathbf{e}}_{\mathbf{z}}\right)g_n={\eta }_0\dot{\gamma }\sum _{n=2}^{n_{crit}}\frac{9}{Mn}\frac{\left(n-1\right)\tan {\theta }_n}{1+{\tan }^2{\theta }_n}{v}_p{g}_n$$

$${\sigma }_{xy}^s={\eta }_0\dot{\gamma }$$

$$\displaystyle+\frac{1}{2}\beta_{n}\left(\left\langle e_{y}^{2}\right\rangle_{n}-\left\langle e_{x}^{2}\right\rangle_{n}\right)+\left(\chi_{n}-2\lambda_{n}\beta_{n}\right)\left\langle e_{x}^{2}e_{y}^{2}\right\rangle_{n}\Big{]}$$

$${\sigma }_{ij}=-\frac{1}{V}\sum _k\left[m^{\left(k\right)}{\tilde{u}}_i^{\left(k\right)}{\tilde{u}}_j^{\left(k\right)}-0.5\text{sideset}{}^{\ast }{\sum }_l\left(r_i^{\left(l\right)}-r_i^{\left(k\right)}\right)f_j^{kl}\right]+{\sigma }_{ij}^{{}^{\prime }}$$

$$\left[\eta \right]={\varphi }^{-1}\left(\frac{{\sigma }_{xy}}{{\eta }_0\dot{\gamma }}-1\right)$$

$${\sigma }_{ii}=\frac{1}{2V}\sum _k\sum _l\left(r_i^{\left(l\right)}-r_i^{\left(k\right)}\right)f_i^{\left(kl\right)}$$

$$N_1={\sigma }_{xx}-{\sigma }_{yy}=\frac{6\bar{f}}{Mn}{\eta }_0\dot{\gamma }\sum _{n=2}^{n_{max}}\left(n-1\right)\cos (2{\theta }_n)g_n{v}_p$$

Full text

\makeFNbottom

[TABLE]

††footnotetext: a* ETH Zurich, Computational Physics for Engineering Materials, Institute for Building Materials, Stefano-Franscini-Platz 3, CH-8093 Zurich, Switzerland.††footnotetext: b University of Latvia, Raina bulvaris 19, LV-1586 Riga, Latvia. E-mail: dmitrijs.zablockis@lu.lv*

1 Introduction

Self-assembly is a pathway to novel materials1 with structural hierarchy, multi-functionality and programmed response to mechanical stress or external field2. In magnetic colloids - ferrofluids - the magnetic nanoparticles are suspended in a non-magnetic liquid carrier and contain a spontaneously polarized monodomain 3. The ability of the ferrofluids to actively interact with a magnetic field, while retaining fluidity is the basis of their practical utility in various technological devices4 and novel applications in magnetic chemistry5 and medicine6. The suspensions of magnetic particles7 exhibit enhanced viscoelasticity as so far as to display a rapid, tractable and reversible field-induced rheological transition between a mostly Newtonian liquid and a plastic solid, enabling magnetomechanical devices and technologies, which take advantage of the tunable flow properties and can be continuously adjusted for improved performance in automotive8, mechatronic or biomechatronic/prosthetic applications 9 or for protection of civil constructions, machine tools and industrial installations against vibrations or seismic shock 10, 11. The electrical analogs, i.e. ferroelectric suspensions of highly polarizable particles in a nonpolarizable solvent, also display an appreciable enhancement of flow resistance 12 under electric field. Magnetic colloids are adaptable platforms for composite soft magnetic materials with customized physical properties 13, 14, 15 also exhibiting viscoelastic behavior16, 17. The new generation of magnetic colloids based on nanoparticles with accurately tailored size and/or oxide-free magnetic materials feature greatly enhanced interaction strength 18, 19, 20, 21, 22. The anisotropic bonding of the polarizable particles or spontaneously polarized nanoparticles and their directional assembly 2 under exposure to the external electromagnetic fields induces the thickening of the dipolar materials subjected to a macroscopic flow22, 12. The associated nonlinear rheology is sensitive to the flow-induced conformational changes of the dissolved microstructure. The safe and effective application of these systems in diverse areas, i.a. related to human health, relies on understanding and quantifying their microstructural and rheological properties. The characterization of the mechanisms providing the rheological specificity facilitates the development of more efficient materials and technological devices and processes, where their design extensively relies on phenomenology or structural rheology23. One of the difficulties is the absence of a rigorous micromechanical model for the associated micro-macro structure-property relationships. The involved self-assembly scenarios are difficult to quantify and despite a large amount of experimental data there is only limited reliable theoretical understanding beyond simpler models confined to the vicinity of the equilibrium state. Direct methods for correlating and interpreting experimental measurements across vastly different systems would deliver substantial advance to the field. In turn, for suspensions of permanent dipoles, e.g. ferromagnetic (or ferroelectric) colloids, where thermal fluctuations play a significant role, comparatively little is known about the structural transitions and conformational changes introduced by the competition of the orientational influence of the shear flow and an external electromagnetic field. The usual dimensional scaling is not applicable in this regime and the interpretation of experimental data is particularly challenging. In this paper we consider the microstructural transformations and the associated rheological response of sheared suspensions of thermalized spherical particles with embedded point dipoles as a model for dipolar colloids in an external field. We verify several developed micromechanical and micro-macro (rheological) structure-property relationships using the controlled environment of particle-level simulations to resolve the interplay between the anisotropic pair interactions, thermal fluctuations, solvent-mediated hydrodynamics and external fields. We employ a theoretical analysis and perform a validation against available experimental data to access a wide range of relevant scales and physical parameters for dipolar systems of substantially different physical nature to demonstrate the universality of the approach.

The paper is organized as follows: In Section 2 we describe our model and the simulation approach. We exploit a hybrid MD-MPCD method implementing coarse-grained hydrodynamics. In Section 3 we conduct simultaneous measurements of the micro (structural) and macro (rheological) characteristics: we start by revisiting some features of the equilibrium state (Section 3.1.1). In Section 3.1.2 we consider the conformational transformations and the erosion of microstructure in a simple shear flow. The viscometric properties are reported in Section 3.2 and the comparison with experiments is discussed in Section 4 in terms of the master curves. Finally, in Section 5 we consider the shear-induced mesophase transformations. To conclude, Section 6 provides the summary of this work.

2 Model

Non-equilibrium molecular dynamics simulations provide a direct quantitative link between the microstructure and the macroscopic properties. Free-draining Brownian dynamics has been used extensively in the past 24, 25, 26, 27, however, the idealized oversimplified hydrodynamics yields poor quantitative agreement with the experiments, underestimating the rheological properties by at least an order of magnitude 24, 28, 29, 30. Stokesian dynamics 31 includes a more rigorous treatment of the particle-solvent coupling and many body hydrodynamics, but is restricted to very small systems and/or reduced dimensionality: monolayer simulations32 (25 particles) of the ER effect reported a fair quantitative agreement with measured rheometric data and confirmed the relative importance of the hydrodynamic contributions to the total suspension stress. A lattice Boltzmann study 33 has shown that the role of the hydrodynamic interactions remains modest ($\phi=0.03-0.20$ and $\lambda=4,8$) in unsheared ferrofluids. In this work we use multi-particle collision dynamics34 employing a stochastic bath of point-particles for coarse-graining the solvent and thermalization. This method also resolves the hydrodynamic modes on the desired length scale and is computationally efficient.

2.1 Dimensionless parameters

The structural and rheological complexity of sheared dipolar colloids is sufficiently described in terms of the volume fraction $\phi$, the dipole-dipole interaction parameter $\lambda$ and the Mason number $Mn$. The interaction parameter $\lambda$ characterizes the ability of the colloid to form an equilibrium structure and expresses the ratio of the characteristic energy of the dipole-dipole interaction between the particles to the average energy of their thermal fluctuations:

[TABLE]

where $\mu_{0}$ is the magnetic constant, $\mu$ - magnitude of the particle dipole moment, $\sigma$ - distance of closest proximity, i.e. the effective diameter of the particles, including the thickness of the solvation shell. The Mason number35 is a ratio of the erosive hydrodynamic forces to the bonding polarization forces encompassing the aggregate stress and carries the meaning of a reduced shear rate

[TABLE]

where $\eta_{0}$ is the viscosity of the solvent and $\dot{\gamma}$ is the characteristic shear rate. The Mason number and $\lambda$ are related via the Peclet number $Mn=\frac{2}{3\lambda}Pe$. The temperature only enters the $\lambda$ parameter - for nanometric colloids and finite $\lambda$ the thermal effects remain an important factor.

2.2 Simulation method

The dipolar colloids are modeled as a monodisperse solution of spherical grains with a diameter $\sigma$ and an embedded point dipole $\bm{\mu}$. The hybrid approach 36 combines the molecular dynamics (MD) simulation of the interacting colloids immersed in a coarse-grained fluid governed by the multi-particle collision dynamics (MPCD) 34, 37. The MD-MPCD38, 39, introduced by Malevanets and Kapral, is an established approach to study particle suspensions 40, 41, instabilities 42, 43, colloids with isotropic attractions 36, 44, anisometric colloids 45, polymer-like models 46, 47, 48, 49; polyelectrolytes50, membranes 51, active matter 52, 53, 54 and self-assembly scenarios 55. For the first time this method is applied here to a system with anisotropic (dipolar56) interactions

[TABLE]

The excluded volume interactions between the colloids, incl. the steric hindrance of the solvation shell, are modeled by the repulsive Weeks-Chandler-Andersen (WCA) potential ($n=6$), which is a combination of soft potentials

[TABLE]

where $\bm{r_{ij}}$ is the center-to-center distance between particles $i$ and $j$, $\beta=\left(kT\right)^{-1}$. All lengths are given in units of particle diameter $\sigma$.

In the MPCD approach the coarse-grained solvent is implemented by an ensemble of point particles with mass $m_{f}$. Their dynamics proceeds in alternating streaming and collision steps. During the streaming the particles propagate ballistically during $\Delta t_{mpcd}$

[TABLE]

For the subsequent collision step the simulated volume is equipartitioned and the particles are assigned to cubic cells (bins) $a\times a\times a$, which define the effective collision environment. The multiparticle collisions are modeled by collective stochastic rotations, whereby the peculiar (thermal) velocities are transformed on a bin-wise basis

[TABLE]

where $\bm{v}_{i}$ is the velocity of ith MPCD-particle in the Jth bin containing $n_{J}$ particles, i.e. $i\in\left\{1,\dots,n_{J}\right\}$; $\left\langle\bm{v}\right\rangle_{J}=n_{J}^{-1}\sum_{i=1}^{n_{J}}{\bm{v}_{i}}$ is the velocity bias and $\bm{\Omega}_{J}\left(\alpha\right)$ is the matrix for a rotation by an angle $\alpha=\frac{\pi}{2}$ around one of the six orthogonal axes aligned with the coordinate axes of the simulation box34, 57, 36 randomly selected for every bin and collision step. The local scaling factor $\kappa_{J}$ is used to maintain a constant temperature $T$ in non-equilibrium simulations. The bin-wise conservation of the streaming momentum by the coarse-grained collision scheme leads to correct hydrodynamics on the length scales $r>a$. The stochastic multi-particle exchanges induce the viscosity of the MPCD fluid58. The colloid-solvent coupling is achieved at the streaming step 36 using the stochastic BCs 59, 36 to enforce the no-slip condition. The momentum exchanges are calculated explicitly36, imparting force and torque on the dissolved solutes. Hence, the colloids are correctly thermalized and the translational and rotational Brownian motion is recovered, because momentum conservation is ensured at each collision. The choice of the simulation parameters preserves the hierarchy of the physical and simulation timescales as described previously 36, 60 to ensure consistency with typical colloidal systems: the ratio of the mass of the colloid $M_{p}$ to the mass of the MPCD particle $m_{f}$ is $79$, the nominal number of the MPCD particles per bin is $\left\langle n_{J}\right\rangle=6$. The extent of the collision environment $a$ defines the spatial accuracy of the hydrodynamic flow60 and is set to $a=0.5\sigma$ to properly resolve the flow on the scale of the aggregates $n\geq 2$, whose contribution to the viscosity scales as $\propto n^{3}$. The mean free path of the MPCD particles is $l_{mpcd}=\Delta t_{mpcd}\left(\sqrt{\beta m_{f}}\right)^{-1}=0.3\sigma$ ensuring Galilean invariance58, 36. The time between the multi-particle collisions is $\Delta t_{mpcd}=0.034\tilde{t}$, where $\tilde{t}=\sigma\sqrt{\beta M_{p}}$. The simulations are performed with LAMMPS61, 62 in a rheometric slit channel of extent $L_{x}$, $L_{y}$, $L_{z}$ (volume $V=L_{x}L_{y}L_{z}$). For zero-field simulations $L_{x}=L_{y}=L_{z}=24\sigma$, whereas, for simulations with an applied field (in $y$ direction) the size of the box is increased for larger interaction parameters $\lambda$: $L_{x}=L_{z}=16\sigma$ ($\lambda<4$) or $24\sigma$ ($\lambda\geq 4$) and $L_{y}=32\sigma$ ($\lambda<4$), $48\sigma$ ($\lambda=4$), $72\sigma$ ($\lambda\geq 5$) to accommodate the growing microstructure. Periodic BCs are used in the $x$ (flow) and $z$ (vorticity) directions, in the $y$ (gradient) direction the channel is bounded by the free-slip walls realizing reflections $v_{n}\rightarrow-v_{n}$ of the incoming MPCD particles and a short-range repulsive potential of the form (4) is applied to the approaching colloids to prevent percolation of the microstructure at larger $\lambda$. In the flow simulations the shearing is induced by adjusting the streaming velocities $\left\langle\bm{v}\right\rangle_{J}$ of the MPCD particles in the bins nearest to the channel walls. Bin-wise sampling was used to measure the flow profile and determine the actual shear rate $\dot{\gamma}$. Every simulation was repeated between $4$ and $32$ times to produce statistically independent trajectories and reduce scatter in the measured observables. The error is reported in terms of the standard deviation.

In the following we consider the two limiting situations:

(a) strong external fields $B\rightarrow\infty$ and

(b) zero external field $B=0$ with permanent dipoles

. The strength of the dipolar interactions is varied by adjusting the $\lambda$ parameter.

A proximity-based clustering criterion is used to identify and quantify the self-assembled microstructure: the particles are considered linked if their center-to-center distance does not exceed $r_{c}=1.1\sigma$. Other definitions63, 64 can be employed to target the specific structural elements, e.g. linear chains. The strong anisotropy of the dipolar interaction implies that any two particles can only approach in a head-to-tail configuration forming linear structures (chains, rings or branched clusters), hence, $r_{c}$ corresponds to the attractive energy $-1.5\lambda kT$, i.e. $\frac{3}{4}$ of the equilibrium bond energy, consistent with the previous studies 64.

3 Results

3.1 Conformational properties

To characterize the conformational and statistical properties of the self-assembled microstructure the function $P_{n}\left(\bm{r}\right)$ is introduced, which is a positional probability density distribution function (PDF) due to the fact that the quantity $P_{n}\left(\bm{r}\right)dV$ expresses the probability of finding a colloidal particle (monomer) in the vicinity of point $\bm{r}$ relative to the center of mass of an identified aggregate ($n$-mer). It obeys a normalization condition

[TABLE]

Figure 1 shows the time-averaged aggregate PDF $P\left(\bm{r}\right)=\sum_{n=2}^{\infty}P_{n}\left(\bm{r}\right)$ in equilibrium at zero applied field (Note: single monomer contributions are substracted here for better visualization) representing compact globular clusters, which initially comprise relatively narrow shells of probability corresponding to the orbits of the monomers around the center of mass of the aggregates; and become rather diffuse with outward distance due to the flexibility of the fluctuating chains. The $n$-mer size (chain length) distribution $g_{n}$ representing the number of $n$-mers per unit volume of the colloid under the specific conditions of colloid concentration $\phi$, dipolar interaction strength $\lambda$ or external field is obtainable from the PDF:

[TABLE]

where $v_{p}=\frac{1}{6}\pi\sigma^{3}$ is the volume of a colloidal particle. The distribution $g_{n}$ determines the statistical and rheological properties of the ensemble of dipolar particles.

3.1.1 Equilibrium

The onset of chaining can be estimated from the second $\sim\phi^{2}$ virial coefficient $b_{2}$65 of the dipolar colloids in zero $-\frac{e^{2\lambda}}{3\lambda^{3}}$ and strong field $-\frac{e^{2\lambda}}{3\lambda^{2}}$. From $|b_{2}|\phi\approx 1$ for moderate volume fraction $\phi=5\%$: (a) in zero field $\lambda_{0}\approx 4$, (b) in strong field $\lambda_{\infty}\approx 3$. Hence, when the dipole-dipole interactions are strong enough to compete with the thermal fluctuations, the particles spontaneously aggregate, forming chains of various length (Figure 1).

On further cooling in zero field the colloid undergoes a hierarchical cascade of structural transformations mediated by sequentially emerging series of topological defects66, which will result in the percolation of the system of branched clusters (Figure S1, Supplemental Materials) 67.

Under a strong applied field the clusters are aligned and unbranched (Figure 2). The corresponding monomer density distribution displays a strongly anisotropic needle-like structure composed of narrow maximums in the monomer positional probability, which are spaced approx. $0.5\sigma$ apart due to the contributions from the chains with odd and even number of particles. The maximum probability is concentrated in the immediate vicinity of the center of mass of the aggregate, which corresponds on average to relatively short chains containing a few particles. The diffuse tail of the PDF persists significantly beyond the few central maxima.

At small lateral separation the interaction potential $\sim e^{-2\pi r}$ between straight dipolar chains of infinite length decays exponentially with $r$ - the interchain distance 68, 69. The parallel chains align with an offset by half a diameter $\sigma$, by which they can produce a short-range attraction confined to the near field ($r<2\sigma$) 68, 70. Finite size corrections 71, 69 amount to repulsive monopoles bound at the chain ends. The flexibility 71 or an offset by half the chain length 69 and polydispersity of the chain ensemble minimize the repulsion of the chain ends and enable the chain-chain approach into the near field (Figure 3). The far field lateral attraction $\propto nr^{-5}$ between the two field-aligned $n$-particle chains is originated by the dipole density fluctuations 68, 69 resulting in the condensation of chains into tapered bundled structures (Figure 3c). The zero-temperature $\lambda\rightarrow\infty$ ground-state for the ensemble of fully aligned dipolar particles should be a body-centered tetragonal (bct) solid 72 comprised by chains offset in the direction of the field 73, 74, 75. At ($\lambda>4-7$) the phase diagram of aligned dipolar particles shows a chain and bct phase coexistence 76.

3.1.1.1 Chain length distribution.

The chain ($n$-mer) length distribution has been calculated using chemical kinetics 77 or the density functional approach 78. According to the later method 79, 80, 81 the chain length distribution $g_{n}$ should provide the minimum of the free energy density functional of the polydisperse ensemble of chains:

[TABLE]

where the first term corresponds to the entropy of an ideal mixture of $n$-mers and the second absorbs intra-chain correlations and interaction with the external field 79, 80. The dimensionless internal free energy is $f_{n}=-\ln Z_{n}$, where $Z_{n}$ is the partition function of a self-assembled chain (S1) (Supplemental Materials). The minimization of eq. (9) is provided by79, 80

[TABLE]

where $\epsilon_{*}=\ln Z_{2}$ and $-\epsilon_{*}kT$ is the average bond energy 80 accounting for the dispersive thermal fluctuations ($\epsilon_{*}=2\lambda$ in the zero-temperature limit82) The coefficient $X$ is determined from the condition of particle conservation

[TABLE]

The average number of grains in a chain is 82

[TABLE]

The pair partition function $Z_{2}$ in zero and infinite field has been calculated 81 (S2)-(S3) (Supplemental Materials). The asymptotes for $\lambda\gg 1$ 79, 82

[TABLE]

coincide with the second virial coefficients 65. The equations (14) provide thermal corrections to the zero-temperature behavior $\epsilon_{*}=2\lambda$ by renormalizing the energy of the bond to account for the thermal fluctuations, which effectively inflate the chain. The predictions of the chain model at low to moderate $\lambda$ significantly depend on the nature of the dependence $\epsilon_{*}\left(\lambda\right)$. To determine the influence of the thermal fluctuations, we exploit the effective parameter $\lambda_{*}$: $\epsilon_{*}=2\lambda_{*}$, where $\lambda_{*}$ is used to fit the theoretical model to the simulation data and characterizes the effective structuring ability of the dipolar bond.

Under a strong field the particle chains are aligned and their length is increased by approximately a factor of $\propto\sqrt{\lambda}$ (13). The calculated chain length distribution (8) shows that very long chains become accessible at larger $\lambda$ at the cost of the reduced number of free floating monomers and shorter chains. The $n$-mer length distribution smoothly persists across many orders of magnitude and has an exponential character as predicted by eq. (10). There is a transition between two exponentially distributed regimes with different exponents, which is more apparent for larger $\lambda$. The deviation from the single exponential decay for longer chains has been observed experimentally 21. The average length of a chain calculated using $g_{n}$ (13) shows good agreement with the predictions of eqs. (10)-(14) (Figure 4a, right inset).

The appropriate scaling for $g_{n}$ follows from equation (10) and a master curve can be obtained by plotting in rescaled coordinates $\left(\alpha_{\lambda}n,\beta_{\lambda}g_{n}\right)$ (Fig. 4b), where

[TABLE]

For $\lambda>4$ only the first exponential regime collapses satisfactorily, whereas the tail of the $g_{n}$ progressively deviates, reflecting the lateral condensation of the chains, forming thick bundles (Fig. 3). For $\epsilon_{*}\gg 1$ to leading order $\alpha_{\lambda}\approx e^{-0.5\epsilon_{*}}$. The scaling coefficients estimated from collapsing the simulation data are plotted in the right inset of Fig. 4b. The asymptotes are provided by the zero temperature limit $\epsilon_{*}=2\lambda$, the thermal corrections (14) are indicated by the dashed lines. The scaling coefficients gradually approach the predicted limiting behavior as $\lambda$ is increased, but the influence of the thermal fluctuations is substantial in the whole investigated range.

3.1.2 Non-equilibrium

Simple shear $\bm{v}=\left(\dot{\gamma}y,0,0\right)$ is the main rheometric flow governed by a single viscosity coefficient. The majority of the rheological data is reported in terms of its shear rate $\dot{\gamma}$ dependence. Figure S2 (Supplemental Materials) presents a simulation snapshot of the field structured colloid in simple shear and the corresponding PDF is shown in a projection onto the shearing plane $P\left(x,y\right)=\int{P\left(\bm{r}\right)dz}$ (Figure S3). Qualitatively, the hydrodynamic stress erodes the microstructure: the number of longer chains is reduced in comparison to the equilibrium state, which is also accompanied by the visible collapse of the tail of the PDF under shear and a slight decrease of the average chain length. The deviation of the chains from their alignment in the direction of the field is reflected in the deformation of the PDF.

The microstructural model of chain-forming dipolar colloids is constructed based on the balance of the micromechanical stress79, 83. In the chain model83 the aggregates of dipolar particles are approximated by a straight rod-like chain (Figure 5), comprising particles in close contact, spaced equidistantly along the chain and the axis (director $\bm{e}_{r}$) of the linear $n$-mer forms an angle $\theta_{n}$ with the direction of the field. Here we assume a complete alignment of the individual dipole moments in the direction of the applied field. The total torque balance determines the orientational response of the rigid $n$-particle rod to the action of the hydrodynamic shear ($\bm{\Gamma^{h}_{n}}$) and polarization ($\bm{\Gamma^{m}_{n}}$)

[TABLE]

Tallying the dipolar ($\bm{F^{m}}$) and hydrodynamic ($\bm{F^{h}}$) force couples starting from the center ($i=1$) of the chain84: $\bm{\Gamma^{m,h}}=2\sum_{i=1}^{\nu_{n}}i\sigma\bm{e_{r}}\times\bm{F_{i}^{m,h}}$ where $\bm{e_{r}}$ is the director vector of the chain and $\nu_{n}=\frac{n-1}{2}$.

The energy of the magnetic bond $U_{dd}\left(\bm{e_{y}},\bm{e_{y}},\bm{e_{r}}\right)=-\lambda kT\left(3\cos^{2}\theta_{n}-1\right)$ produces a magnetic force couple: $\bm{F_{\theta}^{m}}=-\frac{1}{\sigma}\frac{\partial U_{dd}}{\partial\theta}\bm{e_{\theta}}=-6\lambda kT\sigma^{-1}\cos\theta_{n}\sin\theta_{n}\bm{e_{\theta}}$ and $-\bm{F_{\theta}^{m}}$ acting on the ith and (i-1)th particles. Summing the individual contributions along the chain, a realigning torque is provided by the end particles, where the magnetic force couples are not compensated $\bm{\Gamma^{m}}=-6\left(n-1\right)\lambda kT\cos\theta_{n}\sin\theta_{n}\bm{e_{r}}\times\bm{e_{\theta}}$.

In turn, assuming the free draining limit, the hydrodynamic force acting on the ith particle of the $n$-particle chain in simple shear $\bm{F^{h}_{i}}=\beta i\dot{\gamma}\sigma\cos\theta_{n}\bm{e_{x}}$ and the deviating torque $\bm{\Gamma^{h}}=\frac{1}{12}n\left(n^{2}-1\right)\beta\dot{\gamma}\sigma^{2}\cos^{2}\theta_{n}\bm{e_{r}}\times\bm{e_{\theta}}$.

The torque equilibrium yields the deviation angle of the $n$-mer

[TABLE]

The typical aggregation/disintegration scenario is tentatively governed by eq. (17), whereby the chains gradually grow in length by attaching monomers and other free-floating segments and deviate quadratically with $n$ from the direction of the field. The fragmentation is induced by the failure of the bond vicinal to the center of the chain, where the highest tensile stress is applied 85, and the fragments realign.

3.1.2.1 Chain erosion.

In shear flow the longer chains are destroyed by the hydrodynamically applied tensile stress. The simulations show that increasing the shear rate $\dot{\gamma}$ the chain length distribution is gradually eroded at $n\gg 1$, however, retaining the double-exponential form. This behavior suggests that a similar scaling may pertain and a universal curve can be obtained in a similar way as previously by collapsing the non-equilibrium chain-length distribution $g_{n}$ onto the equilibrium distributions at vanishing shear rates $Mn\rightarrow 0$ for each $\lambda$, which are then collapsed onto a single master curve for all $\lambda$ using the previously determined equilibrium scaling factors $\alpha_{\lambda}$ and $\beta_{\lambda}$ (15) (Figure 4).

In the result, Figure 6 shows the chain-length distributions $g_{n}$ calculated from the simulation data for various $\lambda=2-5$ and $Mn$ and plotted in the coordinates $\left(\alpha_{\lambda}\alpha_{Mn}n,\beta_{\lambda}\alpha_{Mn}^{-1}\beta_{Mn}g_{n}\right)$, where $\alpha_{Mn}$ and $\beta_{Mn}$ are the corresponding non-equilibrium horizontal and vertical scaling factors (Note: the vertical scaling is achieved by the factor $\alpha_{Mn}^{-1}\beta_{Mn}$). At $\lambda=4$ we observe an anomalous behavior for some shear rates from $Mn\approx 0.03$ to $Mn\approx 0.15$ (Fig. S5, Supplemental Materials), where instead of being eroded, the tail of the distribution $g_{n}$ recovers with increasing shear rate. This evidences that a magnetohydrodynamic instability is in effect and the system enters a lamellar phase (see also Section 5). In this regime the chains of nanoparticles rearrange themselves to form mesoscopic layers in the field-velocity plane. The associated shear-banding reduces the hydrodynamic stresses within the layered phase and the chains gradually extend. The Mason number where this happens is on the very edge of the experimental window for the magnetic colloids. In turn, for higher interaction strength $\lambda=5$ the phase boundary shifts towards lower shear rates and the delamination occurs before the shear flow can significantly erode the chain length distribution, hence, only a slight dispersion of the collapsed data is observed in Figure 6. Therefore, the same scaling as in equilibrium is also useful in the lamellar phase.

The scaling factors $\alpha_{Mn}$ and $\beta_{Mn}$ are plotted in Figure 7 as a function of the Mason number $Mn$ for various investigated interaction strengths $\lambda=2-5$. We noticed that $\alpha_{Mn}\cdot\beta_{Mn}\approx 1$ and $\alpha_{Mn}-1\approx 1-\beta_{Mn}$ for lower shear rates $Mn<10^{-2}$. The chain model, however, does not support the construction of the chain length master curve in non-equilibrium, so it is difficult to interpret this behavior.

Considering the radial force balance of a rigid chain (Figure 5): the $n$-particle chain is broken, when the strength of the bond $F_{r}^{m}=-\frac{\partial U_{dd}}{\partial r}=-3\lambda kT\sigma^{-1}\left(3\cos^{2}\theta_{n}-1\right)$ is overcome by the overall shear force extending the chain $F^{h}_{r}=\sum_{i=1}^{\nu_{n}}\bm{F_{i}^{h}\cdot\bm{e_{r}}}=\frac{1}{8}\left(n^{2}-1\right)\beta\dot{\gamma}\sigma\cos\theta_{n}\sin\theta_{n}$. The equation for the maximum length of a stable chain is produced from the radial mechanical balance of the bonding force versus the hydrodynamic erosive force $F_{r}^{m}+F_{r}^{h}=0$ accounting for eq. (17): $4n^{4}+5n^{3}-2n^{2}-3n-\frac{288}{Mn^{2}}=0$. To a leading order in $n$ the critical length is

[TABLE]

Hence, any $n$-mer with $n>n_{crit}$ is broken up by the shear flow.

Figure 8 shows the maximum number of particles in an identified aggregate that has been observed in a simulation trajectory at various interaction strengths $\lambda$ and shear rates $Mn$. This quantification is not exact, because it will depend on the duration of the production run, but considering the exponential decay of the chain length distribution it provides a good estimate of the actual critical length. The proportionality to $Mn^{-0.5}$ of the critical length $n_{crit}$ predicted by eq. (18) is not observed. In fact, we observe a non-monotonous dependence of $n_{crit}$ on the shear rate $Mn$ at some $\lambda\geq 4$, which is consistent with the transition into the shear-banded state and the associated decrease of the erosive hydrodynamic stresses within the layered phase. At $\lambda<4$ we indeed observe a power law decay $n_{crit}\propto Mn^{-p}$ but with a lower exponent than the predicted $0.5$ (18).

The concept of a critical chain 85, 86 enables to calculate the non-equilibrium chain length distribution in an empirical way83. The distribution $g_{n}$ can be estimated in simple shear by carrying out minimization (10) numerically with $\epsilon_{*}=-\beta U_{dd}\left(\theta_{n}\right)=\lambda\left(3\cos^{2}{\theta_{n}}-1\right)$ under the conservation condition (11), introducing the critical length (18) as an upper bound in eq. (11). The use of the global extremal principle is justified in the vicinity of the equilibrium, which for magnetic colloids corresponds to the whole range of the practically accessible shear rates $Mn\ll 1$.

Figure 8 shows the average number of particles in a chain (13) under shear and in strong applied field. For all investigated $\lambda=2-5$ the plateau, corresponding to the low shear rate regime $Mn\rightarrow 0$, is followed at some critical $Mn>10^{-2}$ by a near-power law decrease of the average chain length $\left\langle n\right\rangle$, which evidences a rapid breakdown of the microstructure. The theoretical predictions of the non-equilibrium chain model were obtained by fitting the theory to the simulation data in the plateau region, corresponding to the average chain length in equiilibrium, using a fitting parameter $\lambda_{*}^{\left\langle n\right\rangle}$, i.e. $\epsilon_{*}=2\lambda_{*}^{\left\langle n\right\rangle}$. The result is shown in comparison to the simulation data in Figure 8 exhibiting a good agreement. The theoretical model is not strictly applicable for $\lambda\geq 5$ because of the observed side to side aggregation of the chains forming thick bundles. Still, even in this regime it fairly well reproduces the behavior observed in the simulations and can be used with the appropriate selection of the effective interaction parameter $\lambda_{*}^{\left\langle n\right\rangle}$.

3.1.2.2 Chain deviation.

The field- and shear-induced conformational anisotropy of the microstructure bonded by the dipolar interactions is contained in the gyration tensor of the $n$-mer $G_{ij}^{n}=\frac{1}{n}\sum_{k=1}^{n}\left\langle r^{\left(k\right)}_{i}r^{\left(k\right)}_{j}\right\rangle$, where $r_{i}^{\left(k\right)}$ is the ith component of the position vector $\bm{r^{\left(k\right)}}$ of the kth particle relative to the center of mass of the $n$-mer. It is related to the positional probability distribution function $P_{n}\left(\bm{r}\right)$ through

[TABLE]

The deviation of the chain from the vertical ($y$) direction $\theta_{n}$ is obtainable from the gyration tensor

[TABLE]

In zero field a polymer-like behavior 87 is observed: under shear the spherical symmetry of the PDF is reduced (Figure S4, Supplemental Materials), the monomer probability distribution is compressed in the velocity gradient ($y$) direction and extended in the flow ($x$) direction. Thereby the averaged conformations resemble ellipsoids and the average steady alignment is induced by the flow. The conformations (birefringence) ellipsoid is oriented at an average angle $\left\langle\theta\right\rangle$ (measured relative to the direction of the velocity gradient). We observe signs of the limiting behavior $\left\langle\theta\right\rangle\rightarrow 45^{\circ}$ (extinction angle) for vanishing shear rates $Mn\rightarrow 0$88. For finite shear rates the chains are gradually aligned with the flow ($\left\langle\theta\right\rangle\rightarrow 90^{\circ}$).

In the case of a strong vertically applied field, according to eq. (17), the chain-length distribution $g_{n}$ produces the distribution of the deviation angle $\theta_{n}$ due to a significant quantity of subcritical chains with $n<n_{crit}$. The angle distributions obtained from the simulation data at different shear rates (shown in Figure 9 for $\lambda=5$ and is similar for all calculated $\lambda$) can be partially collapsed onto a single master curve when plotted in the coordinates $\left(n,\beta_{\theta}\theta_{n}\right)$, where $\beta_{\theta}$ is a coefficient, which depends on $Mn$, as predicted by eq. (17).

Initially, $\theta_{n}$ scales as $\propto\sqrt{n}$ in contrast to eq. (17), where a larger exponent was expected. For greater $n$ gradual deviations from the power law are observed, leading eventually to a plateau where the deviation angle is independent of the chain length but depends on the Mason number $Mn$. The coefficients $\beta_{\theta}$ used to obtain the partial collapse are plotted in Figure 10. The expected near-linear scaling of the deviation angle with respect to the Mason number $\theta_{n}\propto Mn$ predicted by eq. (17) is observed only for the smallest $\lambda$, whereas for larger $\lambda$ we observe a power law $\theta_{n}\propto Mn^{\kappa}$ with the exponent $\kappa<1$ (plotted in the inset of Figure 10), which decreases with increasing $\lambda$. The predicted maximum deviation angle of the critical chain is $\theta_{max}\approx 35^{\circ}$ (from (17) and (18)). This is not the case in our simulations. The largest angle we observe is $\approx 10^{\circ}$ for the largest shear rates. The observed discrepancies are most likely due to the flexibility of the chains, which is not accounted for in the chain model. Using the chain length distribution $g_{n}$ the statistical averaging of the inclination angle can be produced similar to eq. (13)

[TABLE]

The average angle $\left\langle\theta\right\rangle$ is presented in Figure 11 as a function of the Mason number $Mn$. At low shear rates ($Mn\rightarrow 0$) there is a linear regime, where the alignment angle increases proportionally to Mn: $\left\langle\theta\right\rangle\propto Mn$.

In contrast to the usual polymer theory where the angular response is solely due to the orientational influence of the shear flow, the orientational resistance87 has no meaning for the field-structured dipolar colloids, where the angle is determined by the simultaneous competition of the shear and external field. Instead, we define the angular compliance $s$:

[TABLE]

Figure 12 shows the orientational compliance calculated according to eq. (22) and plotted in the coordinates $\left(e^{\alpha_{s}}Mn,s_{0}^{-1}s\right)$. For vanishing shear rates and large dipolar strength $\lambda\gg 1$ it can be shown without difficulty that

[TABLE]

In practice, this relation holds starting from $\lambda>4$ as shown in the lower inset of Fig. 12. The simulations show that the horizontal shift factor $\alpha_{s}$ is approximately proportional to $\lambda$: $\alpha_{s}\approx 2.545\lambda$. Using these scaling factors a single master curve can be produced for the orientational compliance for all investigated dipolar strengths $\lambda=2-5$ and shear rates $Mn$.

Figure 11 shows the effective interaction strengths $\lambda^{\left\langle n\right\rangle}_{*}$ and $\lambda^{\left\langle\theta\right\rangle}_{*}$, which were used as fitting parameters to connect the predictions of the chain model with the structural parameters measured from the simulations, i.e. the average chain length $\left\langle n\right\rangle$ and the average deviation angle $\left\langle\theta\right\rangle$ calculated at various nominal $\lambda$. A limiting behavior $\lambda_{*}=\lambda$ corresponds to the rigid chain model79, 78 disregarding any influence of the thermal fluctuations, which is indicated by a solid line in Figure 11. The deviations from the linear dependence can be explained as the influence of the Brownian motion. The second line shows the thermal corrections (14) for a strong applied field. At $\lambda<2$ the asymptotic behavior (14) diverges and should not be taken into account. The simulations indicate that increasing the nominal interaction parameter $\lambda$ the influence of the Brownian motion on the inclination angle reduces faster than the influence upon the length of the chain, for which $\lambda^{\left\langle n\right\rangle}_{*}$ more or less follows the predicted curve. In fact, the rotational diffusivity $D_{rot}$ of a linear $n$-mer diminishes with larger $n$, i.e. $D_{rot}\propto n^{-1}$. The increase of the interaction strength $\lambda$ leads to the increase of the average chain length. Hence, the particles in the chain will still fluctuate at their positions, whereas the chains themselves are non-Brownian units for $\lambda>5$. Overall, the values of the effective interaction parameters, which were used to calibrate the chain model to the simulation data, are completely within the expected range, i.e. bounded by the zero-temperature limit from above and by the thermal corrections (14) from below. To add, the chain-chain interactions, which are unaccounted for in the current model, may also be acting and at least partially responsible for the lower values of the effective parameter $\lambda^{\left\langle n\right\rangle}_{*}<\lambda^{\left\langle\theta\right\rangle}_{*}$. The estimates show that the interchain interactions lead to a decease of the average length of the chains89.

3.2 Rheological properties

3.2.1 Stress tensor

The deviatoric stress tensor of the field-structured dipolar colloids comprises the usual symmetric and an additional antisymmetric component, which appears due to the breaking of the rotational symmetry by the imposed field and the associated violation of the angular momentum conservation79

[TABLE]

The antisymmetric stress increment $\sigma_{a}$ is generated by the system of couples induced by the noncentral dipolar interactions and is proportional to the number density $g_{n}$ of the torque-producing units in the polydisperse ensemble of aligned chains8583

[TABLE]

From eq. (17) $\tan\theta_{n}\approx\theta_{n}\propto s_{n}Mn$, hence, the antisymmetric stress increases linearly with the shear rate generating the low shear pseudo-Newtonian plateau. In turn, the symmetric stress increment $\sigma_{xy}^{s}$ is of hydrodynamic origin and is similar to a suspension stress, which originates from the perturbation of the shear flow introduced by the presence of the self-assembled microstructure. The explicit calculation of the flow field is only possible for simpler anisometric shapes90, 91 embedded in a fluid, e.g. a rigid ellipsoid circumscribing the aggregate92. Zubarev and Iskakova79, 83 approximate a linear $n$-particle chain by an ellipsoid of revolution with axes $\left(\sigma,n\sigma\right)$ and provide expressions for the symmetric stress increment of the ensemble of chains with a length distribution $g_{n}$:

[TABLE]

The total volume of the aggregate is conserved in this approach $V_{ellipsoid}=nv_{p}$. The coefficients $\alpha_{n}$, $\beta_{n}$, $\lambda_{n}$, $\xi_{n}$, $\chi_{n}$ are the geometrical parameters of a chain90, 78 and are given by eqs. (S4)-(S5) (Supplemental Materials). The moments $<\dots>_{n}$ are calculated with respect to the non-equilibrium orientational distribution function $\phi_{n}\left(\bm{e}\right)$, which is the solution of the corresponding Fokker-Planck equation. The general solution is impossible and the perturbative solution79, 78 has not yet been presented for strong fields. Assuming that $\phi_{n}\left(\bm{e}\right)$ is peaked around the deviation angle $\theta_{n}$ (17), justified by the considerations of mechanical stability, the moments can be roughly estimated $\left\langle e_{x}^{2}\right\rangle_{n}=\sin^{2}\theta_{n}$, $\left\langle e_{y}^{2}\right\rangle_{n}=\cos^{2}\theta_{n}$ and $\left\langle e_{x}^{2}e_{y}^{2}\right\rangle_{n}=\sin^{2}\theta_{n}\cos^{2}\theta_{n}$. From eq. (25) one can show that the contribution of the $n$-mer to the antisymmetric stress scales as the length cubed $\propto n^{3}$, whereas for the symmetric stress (26) the calculation is slightly more involved and yields $\propto\frac{n^{3}}{\ln{n}}$.

We have conducted simulated rheometric tests of the colloid in strong applied field and measured all components of the virial stress tensor in simple shear

[TABLE]

where $\bm{\tilde{u}^{\left(k\right)}}=\bm{u^{\left(k\right)}}-\bm{\bar{u}}$ is the thermal velocity of a colloidal/MPCD particle $k$ relative to the measured average velocity field $\bm{\bar{u}}$ and $\sigma^{{}^{\prime}}_{ij}$ is the component of the stress tensor peculiar to the MPCD solvent93, which is associated with the redistribution of momentum in multi-particle collisions $\sigma^{{}^{\prime}}_{ij}=-\frac{1}{V\Delta t_{mpcd}}\sideset{}{{}^{{}^{\prime}}}{\sum}_{k}{r_{i}^{\left(k\right)}\Delta p^{\left(k\right)}_{j}}$, where $\Delta p^{\left(k\right)}=m\left(\bm{u^{\left(k\right)}}-\bm{\hat{u}^{\left(k\right)}}\right)$ and $\bm{\hat{u}^{\left(k\right)}}$ is the velocity before the collision ($\Sigma^{*}$ is taken just over the colloids and $\Sigma^{{}^{\prime}}$- over the MPCD particles; $\Sigma=\Sigma^{{}^{\prime}}+\Sigma^{*}$).

The rheometric curves $\sigma_{xy}\left(\dot{\gamma}\right)$ and the viscosity $\eta=\frac{\sigma_{xy}}{\dot{\gamma}}$ in simple shear are the properties that are most frequently measured in a laboratory experiment. In turn, the intrinsic viscosity is defined as

[TABLE]

which is plotted in Figure 13.

At high interaction parameters $\lambda\gg 1$ the intrinsic viscosity is a universal function of the Mason number $\left[\eta\right]\propto Mn^{-\Delta}$, where the exponent $\Delta$ becomes independent of $\lambda$. Here, the Mason number produces a gradual collapse of the shear thinning region 94 as $\lambda$ increases. In turn, for intermediate $\lambda$ the Mason number alone is not sufficient to collapse the rheological response onto a single curve as it does not include the thermal effects and the exponent $\Delta\left(\lambda\right)$ becomes a function of $\lambda$, which leads to the branching of the viscometric curves. For low to moderate bonding energies $\lambda$ the calculated viscometric curves show at the lowest shear rates $Mn$ a pseudo-Newtonian plateau with a shear rate independent viscosity. At the higher shear rates a rapid shear thinning is observed, which is governed by a near power-law decrease of the intrinsic viscosity. For higher bond energies the power-law regions of the viscometric curves begin to converge, whereas the low shear plateau is elevated at approximately an exponential rate at the same time gradually moving to lower shear rates. At the largest $\lambda$ the plateau is outside the range of the Mason numbers accessible to the simulation and can no longer be resolved. The estimated shear thinning exponent $\Delta$ increases from about $\Delta\approx 0.3$ (at $\lambda=2$) to $\Delta\approx 0.7$ at $\lambda=7$. The limiting value $\Delta=1$ 85 is not reached here. The model of Zubarev and Iskakova 83 explains the deviations from the reciprocal scaling with $Mn$, assuming that the chains are not identical and taking into account the polydispersity of the chain ensemble. The shear thinning exponents $\Delta<1$ are in accord with the experimental observations 23 and Stokesian dynamics simulations in a monolayer 32.

3.3 Normal stresses

Considering the WCA (4) and dipole-dipole (3) interaction potential as a dissipative oscillator, the oscillatory motion of the particles is well damped by the viscous drag and the particles are mostly located at the bottom of the potential minimum, where the attraction force is compensated by steric repulsion. The thermal fluctuations introduce some tensile stresses even in equilibrium. A simple order of magnitude estimate shows that the kinetic contributions of the colloids to the virial stress are negligible. From eq. (27)

[TABLE]

Assuming that a net repulsive force $F_{r}=\frac{2\lambda kT}{\sigma}\bar{f}$ exists between each pair of particles in the chain aligned by a strong external field, the first normal stress difference can be straightforwardly estimated from eq. (29) for a rigid chain as

[TABLE]

where $\bar{f}$ is defined by the interaction potential. Assuming $\cos 2\theta_{n}\approx 1$

[TABLE]

which evidences that the first normal stress difference is, in fact, positive and a property of the bond and, hence, should be more or less independent of the shear rate, except for very large shear rates, when most of the larger bonded aggregates are eroded and the number of free monomers quickly increases $\left\langle n\right\rangle\rightarrow 1$.

The repulsive force in question is originated by the thermal fluctuations27. The added Brownian stress develops from the fluctuation induced strain of the chain structure. A trimer (i.e. three-particle chain) coupled to a Langevin environment is a minimal representative chain model if considering only the nearest neighbor interactions. The positional probability density distribution $P_{3}$ of the central particle in the trimer fluctuating along the director and bound by two neighbors is described by a stationary solution of the corresponding Fokker-Planck equation. In the vicinity of the center of the chain:

[TABLE]

where $\beta U\left(r\right)$ is the potential energy of the central particle in a trimer. The inset of Fig. 14 shows the probability density distribution of a trimer in equilibrium for $\lambda=4$ plotted along the axis connecting the particles. It can be seen that the reference particle is localized in the potential well and the probability practically becomes zero below the chosen clustering cutoff criterion $r_{c}=1.1$. We have also checked the positional PDF of an end particle in a dimer relative to its pair and, despite slight asymmetry, found it similar to $P_{3}$ in shape and amplitude.

The particles in a chain oscillate between repulsive and attractive force. While both contributions on average compensate in a symmetric potential, they do not in an anharmonic potential, e.g. the dipole-dipole interaction coupled with a steric repulsion

[TABLE]

The results (Figure 14) show that the estimate (31) over-predicts the normal stresses calculated from the simulation data, which is not surprising, considering the strength of the approximations. The second normal stress difference follows naturally $N_{2}=-N_{1}$. This is observed in the simulations within several percent (not shown here). Estimates show that the fluctuation induced stress survives the hard sphere limit of eq. (4) $n\rightarrow\infty$ with $\bar{f}\propto n$.

4 Experiments

4.1 Ferrofluids

4.1.1 Dimensional analysis

Magnetic colloids are inherently polydisperse systems, hence, the size variation translates ($\propto\sigma^{3}$) into a variation of the interaction parameter $\lambda$. The fraction of particles that are not large enough to form any structure $\lambda<1$ and exist predominantly in a non-aggregated state can be substantial. We assume a bidisperse structure of magnetic colloids, where the chains are primarily composed of larger particles $\lambda>1$ determining the rheological specificity 95, 96 and the smaller particles $\lambda<1$ are disregarded. Subsequently, we analyze the particle size distribution (PSD) of the investigated materials and renormalize it, discarding the influence of the smallest fraction. The estimated fraction of the large interacting particles is shown in column 10 of Table 1. The average diameter is determined by averaging ($\propto\phi$) with respect to the remaining large fractions. If the normalization of the PSD is not performed to display the influence of the aggregating fraction, the precise position of the experimental curve on the $\left(Mn,\left[\eta\right]\right)$ plot cannot be obtained. Still, this is a rather crude approximation, since the small particles can participate in the interactions with the bigger particles and be included in the chains 97, 98, thereby influencing their length 99. This may significantly impact the rheological response due to the lower tensile strength of the induced defects or otherwise 100.

An additional difficulty comes from the fact that only some of the published investigations report all the relevant material parameters necessary for modeling and quantitative comparison. The magnetometric measurements101 are useful to characterize the magnetic size $\sigma_{m}$ distribution and magnetic moment $\mu=\frac{1}{6}\sigma_{m}^{3}M_{s}$ (here $M_{s}$ is the spontaneous magnetization of the particle material) distribution of the spontaneously polarized particles 102, 103, 104, 105, 106 and estimate the dipole-dipole interaction parameter $\lambda$. One must take into account that the size of the magnetic core can be significantly smaller than the hydrodynamic diameter. The difference is associated with the existence of a stabilizing shell (e.g. a layer of surfactant) covering the monodomain as well as the magnetically "dead" layer on its surface. The total size of the solid grain can be estimated by transmission electron microscopy (TEM) 100, 107. In turn, the dynamic light scattering (DLS)107, 108 can be used to determine the hydrodynamic diameter directly. An exhaustive characterization of the ferrofluid using different methods facilitates the determination of $\sigma_{m}$ and $\sigma$. The hydrodynamic volume fraction $\phi$ contained in (28) is an awkward parameter, which is almost never reported. It is assumed that $\phi=\frac{\sigma^{3}}{\sigma_{m}^{3}}\phi_{m}$, where $\phi_{m}$ is the volume fraction of the dispersed magnetic material, which is usually 2-3 times smaller than $\phi$. The summary of the previously conducted rheological investigations of ferrofluids is shown in Table 1, which is more or less self-explanatory, and the results are presented in Figure 15 in a $Mn-\left[\eta\right]$ plot. In the analysis we use only the measurements at the highest reported magnitudes of the applied field, thereby assuming that the investigated ferrofluids are more or less at saturation, which may not always be completely true.

The experimentally measured viscosities of the various ferrofluids span 8 orders of magnitude across 8 orders of magnitude of the shear rate at quite moderate applied fields. Although some trends are already visible - e.g. increase of the viscosity with a larger dipolar interaction parameter - the interpretation of this chart is not straightforward due to the fact that the measurements were performed at different effective temperatures, i.e. dispersion in $\lambda$. The Mason number alone is not a reliable scaling parameter for ferrofluids and it is not sufficient to collapse the viscometric data. However, the collapse can be achieved with additional processing.

4.1.2 Viscosity master curve

The time-temperature superposition119 is performed to produce the viscosity master curve. The construction is achieved using the shift factors $\alpha_{\eta}$ and $\beta_{\eta}$ and is based on the observed analogy between thermal and hydrodynamic erosive stresses, i.e. the increase of the reduced temperature $T^{*}$ (or decrease of the $\lambda$ parameter) bears a similar effect on the viscoelactic properties of magnetic colloids as the increase of the shear rate expressed in terms of the Mason number Mn. Hence, the isothermal viscometric curves obtained at different temperatures (or $\lambda$) over a range of shear rates can be superposed by a scaling relation.

From eq. (25) and accounting for eqs. (10) and (12) it follows that for vanishing shear rates $Mn\rightarrow 0$: $\frac{\sigma^{a}}{\eta_{0}\dot{\gamma}}=0.75e^{-\epsilon_{*}}\left[Li_{-3}\left(X\right)-Li_{-1}\left(X\right)\right]$, where $Li_{n}\left(x\right)$ is a polylogarithm. For larger dipolar strengths $\lambda\gg 1$:

[TABLE]

In the chain model the low shear rate plateau is dominated by the antisymmetric stress, so equation (34) is a good approximation for eq. (28) at $Mn\rightarrow 0$. It also shows that the intrinsic viscosity is not independent of concentration, which is appropriate for self-assembled systems, and scales exponentially with the binding energy. Equation (34) is formally in accord with the Arrhenius model where $\epsilon_{*}$ carries the meaning of the activation energy. The vertical shift factor is defined using the low shear viscosity as a scaling parameter $\beta_{\eta}=\left[\eta\right]_{0}^{-1}$.

In turn, the horizontal shift factor $\alpha_{\eta}=Mn_{*}^{-1}$ is defined in terms of the critical Mason number $Mn_{*}$, i.e. the point (horizontal intercept), where the straight lines corresponding to the power-law segment $\propto Mn^{-\Delta}$ and the plateaus of the viscometric curves intersect on a log-log plot, which approximately marks a transition from the pseudo-Newtonian plateau, dominated by dipolar interactions, into the shear-thinning regime. Both shift factors scale exponentially with the dipole-dipole interaction parameter at $\lambda\gg 1$, whereas the small deviations at low to moderate values of $\lambda$ are the effect of the Brownian motion and the differences between the nominal and effective values of the interaction parameter discussed previously. It is assumed that for $\lambda\geq 7$ the value of the effective dipole-dipole interaction parameter is the same as the nominal $\lambda_{*}=\lambda$.

In Figure 16 we plot the simulation data using the time-temperature superposition in the normalized coordinates $\left(\alpha_{\eta}Mn,\beta_{\eta}\left[\eta\right]\right)$. As the interaction parameter is increased $\lambda\gg 1$, we see evidence of a limiting behavior, where both the data and the theoretical predictions by the chain model converge onto a single universal master curve. In the past a similar construction was tried by Baxter-Drayton et al.32 based on purely empirical considerations to interpret their monolayer simulations of electrorheological fluids.

Next, we re-plot the experimental data of Figure 15 in the normalized coordinates $\left(\alpha_{\eta}Mn,\beta_{\eta}\left[\eta\right]\right)$, where $\alpha_{\eta}=0.0825e^{2.57\lambda}$ and $\beta_{\eta}=3.58e^{-2\lambda}$ are determined from the time-temperature superposition of simulation data and theoretical predictions at our standard volume fraction $\phi=5\%$. Inspecting Figure 17 we find robust experimental evidence of a universal behavior predicted by the simulations and the chain model. By the use of this procedure the dispersion in $\lambda$ associated with the varying material properties (Table 1) is collapsed to a significant degree. The benefit of this scaling is that it shifts the experimental window, broadening the range of the accessible effective shear rates in terms of the Mason number $Mn$, which allows to examine greater portions of the master curve using different materials, which are not accessible with a single ferrofluid. There are several remaining sources of dispersion in the mapping of the experimental results:

(a) the uncontrolled contribution to the interaction potential comes from the residual short-range isotropic attraction, e.g. leaking through the steric barrier or due to the defects in the protective shell. The Van der Waals interaction is of the order $\propto kT$ and is likely significant in some weakly-associating ferrofluids $\lambda\approx 1-3$ based on synthetic magnetite (Fe3O4) 120

(b) difficulty of obtaining saturation with respect to the external field in some measuring arrangements, e.g. the group of data110, 103, 104, 100 directly underneath the master curve in Figure 17 was measured at $30-40$ mT, which is apparently not at saturation of the magnetic viscosity effect, hence, they display the same dependence, but do not quite reach the predicted magnitude

(c) incomplete compensation of the effect of polydispersity and the dispersion of the interaction parameter, especially at low values of $\lambda$. Similarly, incomplete characterization of the ferrofluids in general, requiring assumptions.

(d) the independence of the intrinsic viscosity on concentration is not completely established (Figure S7, Supplemental Materials) in the vicinity of the plateau, whereas the concentration of the chain-forming particles varies across at least one order of magnitude (Table 1) in the investigated samples.

.

There are many other uncontrolled parameters, which can significantly influence the viscoelastic properties of the magnetic colloid. For one, the amount of free surfactant or primary agglomerates 121 remaining in the solution from the preparation stage is difficult to quantify, but it can significantly influence the measured viscosity.

4.2 Other dipolar systems

The interpretation of the experimental data for ferrofluids is challenging because of the lack of a well-defined model system and the sensitivity of their structural and rheological properties to various parameters, which cannot be controlled or accurately measured. In turn, the dipolar systems, where the magnitude of the interactions can be varied, have proven to be excellent experimental benchmarks to test the theoretical models of dipolar materials23, 122, 123. The "inverse ferrofluids" (i.e. magnetic "holes") are obtained by dispersing inert particles in a ferrofluid matrix124, whereby a tunable permeability mismatch induced by the external field produces a virtual dipole associated with each "hole". The magnetorheology of inverse ferrofluids based on nonmagnetic silica nanoparticles (diameter $\sigma\approx 100-700$nm) dispersed at a volume fraction $\phi\approx 1-26\%$ has been previously reported as a function of the shear rate and magnetic field strength in terms of the Mason number125, 126, 127, 123. The corresponding intrinsic viscosity is reproduced in Figure 18 along with the fits to our simulation data using the chain model showing overall good agreement.

Unlike the magnetic "holes", the superparamagnetic colloids128 are magnetically inert microparticles impregnated with single-domain magnetic grains. Hence, they display a superparamagnetic response to the applied field. The viscoelastic behavior of the composite polystyrene spheres (diameter $\sigma\approx 1\mu m$) laden with maghemite (Fe2O3) and dispersed in an aqueous solution at a volume fraction $\phi=5.2\%$ was previously reported129 as a function of the Mason number $Mn$ and is also plotted in Figure 18. The dipole interaction parameter of the superparamagnetic colloids is rather large $\lambda\gg 1$ owing to the quantity of contained magnetic material ($>67$wt.%) and their response approaches the limiting behavior of the chain model.

5 Phase diagram

The microstructural transformations determine the complicated viscoelastic response of dipolar colloids. In Figure 19 a tentative phase diagram is presented showing various mesoscopic phases that we encountered in the simulations of the field-structured dipolar colloids under shear.

At small $\lambda$ (Brownian phase) the particles are free to diffuse through all spontaneously forming aggregates and the lamellar phase does not form. In turn, in systems dominated by dipolar interactions $\lambda\geq 5$ the chains approach laterally, aggregating into thick chain-like bundles (Figure 3), which can already be identified as persistent individual entities (rather than brief structural fluctuations). The bundles of chains form a hexagonal order (Figure 3a,b) in the velocity-vorticity plane and the lateral mobility of the particles is restricted. Under simple shear the stability of the hexagonal phase is diminished and the perturbation of the dissolved microstructure exposes new pathways in the form of a phase of flowing layers of chains interleaved by slip planes of pure carrier liquid. The competition of the thermal, polarization and shearing stresses directs the associated mesophase transformations. Above the Brownian transition ($\lambda\geq 5$) at weak shear ($Mn<10^{-5}-10^{-4}$) the effect of the flow is limited to affine deformation of the chain bundles (Figure 19a). Increasing the shear rate the hexagonal pattern transforms into a tilted lamellar phase (Fig. 19d and Fig. 20, left), where the individual layers are oriented at approx. $\pm 20^{\circ}$ to the direction of the shear flow. Upon further increase of the shear rate the lamellas abruptly reorient in the direction of the shearing velocity at $Mn\approx 10^{-3}$ (Fig. 19e and Fig. 20, right). The period of the shear banded structure is $L=8\sigma$ for $\lambda=4-7$ and is independent of the shearing rate. This conventional lamellar phase can be directly accessed from the Brownian phase by increasing the strength of the dipolar interactions $\lambda$ at high shear rates $Mn>10^{-3}$. In the vicinity of the critical point ($\lambda_{c}\approx 5$, $Mn_{c}\approx 10^{-4}$), where the Brownian, lamellar and tilted lamellar phases meet, we observe unstable lamellar phases oscillating in the shear flow between $+20^{\circ}$ and $-20^{\circ}$, where the chain bundles are slipping relative to each other in a complicated hexagonal motion. This behavior is masked by the intense Brownian fluctuations of the particles and is only visible on the averaged monomer density snapshots of the vorticity-velocity plane (Figure 19c). We note that any delamination transition is marked by an increase of the critical chain length $n_{crit}$ (Figure 8), evidencing that the transversal erosive stresses are significantly reduced. The previously calculated phase diagram 130 in $Pe-\lambda$ space has only limited qualitative agreement with Figure 19, most likely due to the idealized hydrodynamics used in the previous works. In particular, the existence of a tilted lamellar phase is a new development. Delaminated patterns were reported in an ER fluid under steady 131, 132, 133, 134 and oscillatory shear 135. In turn, suspensions of $0.5\mu m$ magnetic polystyrene particles showed the onset of a regular mesoscopic striped structure above a critical shear rate 136. The phase diagram (Figure 19) shows that the delamination transition is completely within the range of characteristic parameters of ferrofluids (Table 1).

6 Conclusions

We have used hybrid simulations with coarse-grained hydrodynamics to observe the self-assembly of dipolar colloids with significant interparticle interactions $\lambda>1$ and determine their rheological behavior in simple shear across six orders of magnitude of the Mason number $Mn$ for volume fractions $\phi=1-9\%$. The emergence of various aligned structures is observed across different length scales in equilibrium and non-equilibrium scenarios. We show by simulation and theoretically that a master curve exists with respect to the chain-length distribution at moderate values of the interaction parameter $\lambda\leq 4$. Non-equilibrium simulations with imposed shear flow revealed conformational changes in the monomer PDF and a structure defined by the simultaneous competition of field induced alignment, hydrodynamic torque, erosive stresses and dipolar bond strength. We provide evidence of the universality of the structural behavior governed by the competition of bonding (dipolar) and dissociating (thermal and/or hydrodynamic) stresses. The dynamical properties reflect the conformational transformations: for bond strengths $2\leq\lambda\leq 5$, the dipolar colloids display typical polymeric rheology, i.e. showing a pseudo-Newtonian plateau with constant viscosity for vanishing shear rates $Mn\rightarrow 0$ followed by a near-power law shear-thinning behavior $\left[\eta\right]\propto Mn^{-\Delta}$ at higher shear rates. For $\lambda>5$ the plateau is outside the simulation window resulting in an apparent yield stress and pseudo-plastic Bingham-like rheology. The simulations show that a meaningful description of the observed flow behavior can be obtained within a simple model of rigid rod-like chains with statistically distributed number of bonded particles. Quite surprisingly, despite several strong assumptions concerning the flexibility of the chains and disregarding chain-chain interactions (e.g. hydrodynamic), it provides a fair account of both the pseudo-Newtonian plateau and the shear thinning region of the viscometric curve. To ensure this, the model is calibrated to the simulation data using the dipolar interaction parameter $\lambda$ as a semi-empirical coefficient. We have discussed, which approximations should be improved to achieve better quantitative performance of the chain model.

The actual self-assembled microstructure persists across different length scales and is more complex than that of rigid single particle chains: above a critical interaction parameter $\lambda\geq 5$ there is an abrupt shear-induced delamination transition leading to the formation of densified layered (lamellar) phases. Depending on the Mason number $Mn$ three distinctly different flow regimes are observed: at low shear, the hexagonally arranged self-assembled bundles of chains support each other and resist the imposed stress, experiencing only affine deformations. Increasing the shear rate, the bundles start slipping, eventually leading to shear banding and forming of a tilted lamellar phase with layers of chains flowing at an angle $\pm 20^{\circ}$ to the direction of shear. At higher shearing rates the phases become fully aligned with the flow. The emergence of the layered phases affects the stress distribution in the system, which ultimately manifests itself in the lengthening of the chains.

In view of the used idealizations we obtain fair quantitative agreement between NEMD simulations and the viscoelastic response of magnetic colloids reported in various literature sources. This is achieved via the time-temperature superposition supported by the theoretical insight gained from the chain model and without using any fitting parameters. For the first time we have shown here evidence of a universal behavior of the rheology of ferrofluids obeying a single master curve. The observed collapse of the experimental viscometric curves substantiates the usefulness of the modeling approach. In fact, we provide evidence for the universality of the dipolar framework, which can be applied generally across very different length scales and physical systems.

Hopefully, this work will turn out to be a useful advancement towards developing a more general quantitative model of the rheology of dipolar colloids with a robust predictive strength for more complex geometries (incl. rheometric channels) and flows.

7 Acknowledgments

The financial support of Sciex-NMSch (Project 13.124) and the European Research Council (ERC) Advanced Grant 319968 FlowCCS is gratefully acknowledged. E.B. acknowledges the support of the LSR Program IMIS2.