Scale invariance of phase separation dynamics of colloidal suspensions

Here we first consider some characteristic features of colloidal phase separation, through its comparison with phase separation in other materials. For ordinary binary liquid mixtures (mixtures of two liquids with similar dynamics, i.e., dynamically symmetric mixtures), the physics of phase separation is rather well understood on the basis of the concept of dynamic scaling,36 and thus we can predict the process in a universal framework. This scaling is a direct consequence of the presence of unique characteristic length and time scales in the problem, i.e., the correlation length ξ and lifetime τ ξ of critical composition fluctuations, respectively, which are related through τ ξ ~ 6πηξ3/(k B T), where η is the viscosity, k B T is the thermal energy. For dynamically asymmetric mixtures (mixtures of two components with very different dynamics), on the other hand, it is not so obvious whether we can extract such unique length and time scales controlling dynamic phenomena. For example, in polymer solutions there is no scale-invariant nature for dynamics,37,38,39 although the static phase diagram can be universally described after a proper scaling;36,40 more specifically, it is known that for polymer solutions the effective dynamic critical exponent and the pattern evolution during phase separation are significantly influenced by the characteristics of polymers, i.e., the degree of polymerization N.37,38,39 This is a consequence of a non-universal relation between the characteristic lifetime of composition fluctuations τ ξ and the structural (or rheological) relaxation time of the system τ α : the ratio τ ξ /τ α crucially depends on the size of polymer (N),37,38,39 since the internal degrees of polymers affect differently the two types of relaxations.40

For colloidal suspensions, on the other hand, there are no internal degrees of freedom in a colloidal particle. This implies that any characteristic time scales of collective motions of colloids should be proportional to that of the diffusive motion of an individual particle; in other words, τ ξ /τ α should be independent of the colloidal size. Therefore, once we can precisely match the intercolloid potential as well as the control parameters that determine the thermodynamic states (the volume fraction ϕ and the inverse temperature βε in our case), the phase separation dynamics of colloidal suspensions is expected to be scaled using a characteristic time of colloid diffusion (i.e., the Brownian time τ B ) as a unit time. This scale-invariant nature of colloidal dynamics suggests a possibility of direct numerical prediction of dynamical phenomena in real colloidal suspensions without any adjustable parameters.

Numerical simulation based on direct computation of the Navier–Stokes equation

To realize this, the key problem is how to describe the dynamics: in order to precisely predict the kinetics of phase demixing in colloidal suspensions, we need a model with which we can perform simulations within realistic computational time while retaining the essential physical laws governing the phenomena.

To reduce a numerical cost stemming from a huge separation in length and time scales between colloids and solvent molecules, a coarse-grained description of the solvent is essential. Brownian dynamics (BD) simulations may be the most widely used coarse-grained method. In this method, the dynamic couplings of solvent molecules to colloids are simplified as random forces and local friction acting on colloids, and the motions of colloidal particles are described by a set of Langevin equations. For this simplicity, BD simulations are one of the most popular methods. However, this simple version of BD simulations (see Methods) completely ignores HI. This is a crucial deficiency for studying structural formation in colloidal phase separation, since HI may play a crucial role as we mentioned in the Introduction.

This problem can be overcome by introducing HI tensors in the above Langevin-type equations, whose exact expression is known in terms of multipole expansion. This approach is widely known as the Stokesian dynamics (SD) methods.10,11 Since we need not deal with the degrees of freedom of the solvent in an explicit manner in this scheme, SD methods realize high computational efficiency. For example, efficient algorithms for SD simulations whose simulation cost scales as \({\cal O}(N\,{\mathrm{log}}\,N)\)41 for a periodic boundary condition and \({\cal O}(N)\)42 for a free boundary condition (N being the number of colloids) have been reported. Furthermore, thermal noise can be straightforwardly incorporated as random force acting on colloids such that HI tensors introduced are to be consistent with the fluctuation–dissipation theorem. The lowest rank HI tensor is known as the Oseen tensor, which decays slowly as r−1 (r being intercolloid distance). Thus, this becomes the leading-order term only for large r; in other words, if we use only lower rank HI tensors, we can reproduce HI at long distance but cannot capture HI at short distance, or lubrication effect. In principle, inclusion of higher rank tensors allows us to capture HI at shorter distance;42 however, it significantly decreases numerical efficiency, since the number of the tensor elements become explosively larger with an increase in the rank of a tensor. Thus, the near-field HI is often approximated by two-body interaction assuming the additivity.11,14,24,41

Another popular approach is simulations based on molecular dynamics (MD) method such as the hybrid-MD method known as multi-particle collision dynamics (MPCD) or stochastic rotational dynamics (SRD)16,18,19,43 and dissipative particle dynamics (DPD).12 In this type of approach, the dynamics is described by solving Newton’s equations of motion for solvent molecules and colloids. Thus, momentum conservation is satisfied on a local level, and furthermore thermal noise can be naturally introduced. In this case, however, we inevitably need to deal with at least three time scales: the time required for sound wave in solvent to propagate, τ s , the time required for solvent diffusion, τ ν , and the translational Brownian diffusion time, τ B . Since these time scales differ by many orders of magnitude (i.e., τ s ≪ τ ν ≪ τ B ) in real systems (e.g., 1 μm colloids in water at room temperature), we need to make a special care for the separation of these time scales.43

The effect of sound wave can be removed by describing the solvent molecules as an incompressible fluid which obeys the Navier–Stokes equation. In this approach, however, we suffer from a difficulty originating from the (nonslip) solid–liquid boundary condition, which has to be satisfied on the surfaces of all colloid particles in motion. If we directly treat such moving boundary conditions in the Navier–Stokes equation, we need to generate complicated adaptive mesh for every time steps,44,45 which is numerically costly. We have overcome this difficulty by treating a solid colloidal particle as an undeformable fluid particle with high viscosity. We call this method “Fluid Particle Dynamics (FPD) method” (see Methods).20 This idea of embedding the solid–fluid boundary conditions into the Navier–Stokes equation by smoothly connecting the solid and fluid with finite-thickness interface20 have now been adopted to various direct numerical simulations (see, e.g., refs 46,47) By combining this approach with fluctuating hydrodynamics,48,49 we have shown that FPD simulation can almost perfectly reproduce the motion of colloids under thermal noise at low volume fractions while retaining statistical mechanical consistency between colloids and a solvent.22,50

Difference in the Schmidt number S c between experiments and simulations

However, there still remains a difficulty coming from the presence of the inertia term in such direct simulation approaches: we need to deal with two time scales associated with the momentum diffusion of the solvent and the diffusion of colloidal particles at the same time.

The Schmidt number S c = ν/D is defined as the ratio of the two time scales, where ν is the kinetic viscosity of the solvent and D is the diffusion coefficient of an isolated colloidal particle. In typical colloidal systems, these two diffusion coefficients differ by many orders of magnitude (for example, S c ~ 106 for μm-size colloids suspended in water at room temperature), and this condition (S c ≫ 1) is widely assumed in the laws regarding dynamics of colloids such as the Stokes–Einstein relation.49 In other words, the equilibration of the momentum degrees of freedom is a necessary condition to have the scale invariance of the phenomena in the diffusive regime. For such a large Schmidt number, however, numerically accessing the diffusive time regime, in which structural ordering of colloids takes place, starting from the ballistic regime, is almost impossible with current computational power. Indeed, in most of simulation studies, the Schmidt number has been set to modest values.49 Thus, it might look almost hopeless to realize numerical simulations with a true predictive power, or to retain the scale invariance, while incorporating full hydrodynamic interactions and thermal noise in our scheme (see Fig. 1a).

Fig. 1 Difference in the Schmidt number S c between experiments and simulations. a The scaling functions of the mean square displacement of a free particle, h (see SM A.2.1, Eq. (S7), on the functional form) for experiments (EXP) (S c ~ 107: brown curve) and fluid particle dynamics simulations (FPD) (S c = 8.0: blue curve). The brown and blue arrows represent the time scales covered by EXP and FPD respectively. The gray shaded region is the time span over which we compare the coarsening behavior between EXP and FPD. In this time region, needless to say, h in EXP shows the diffusive behavior (h = 1). On the other hand, although we may claim that h in FPD is close to 1, it is not easy to make a rigorous justification that our simulation completely reaches the diffusive regime. b An estimation of the characteristic length l c over which the momentum diffusion of a solvent takes place within τ B . For S c ~ 10, which is a typical value in FPD method, we obtain l c /σ ~ 1 Full size image

However, we show here that this is not necessarily the case. Let us consider a characteristic time required for the momentum of the solvent to diffuse over a characteristic distance l c from a particle, which is given by \(\tau _

u (l_c) \cong l_c^2/

u\). Denoting the colloid diameter as σ and the Brownian time as \(\tau _{\mathrm{B}} = \left( {\frac{\sigma }{2}} \right)^2/6D\), then \(\tau _

u /\tau _{\mathrm{B}} \cong 24(l_c/\sigma )^2S_{\mathrm{c}}^{ - 1}\). This means that how far the momentum diffuses away within time scale of τ B depends on a setting of S c . For S c ~ 10 and τ ν /τ B ~ 1, for example, we obtain l c /σ ~ 1, implying that a flow field around a particle is established by momentum diffusion at least up to a distance of about the particle diameter from the particle (see Fig. 1b). This indicates that, for S c ~ 10, our method can properly describe near-field many-body hydrodynamic interactions, which is expected to play a dominant role in the aggregation process of colloids undergoing Brownian motion. We show that the effect of the inertia term is negligibly small for S c ≳ 10 in Supplementary Material (SM) A.2.3 (Fig. S5). Furthermore, the short interaction range of our potential, 0.13σ (see below), is much smaller than l c ~ σ at S c ~ 10. These support the validity of our simulations with S c ≳ 10 for the problem of colloidal phase separation in our system (see also SM A.2.3 and Fig. S4 on the details). Hereafter, we show that a precise numerical prediction of colloidal phase separation is indeed possible with this approach.

Detailed comparison between experiments and simulations

To verify numerical predictability, it is essential to have experimental data that allow us to precisely compare with the corresponding simulation results on a quantitative level. To this end, we experimentally study the kinetic processes of phase demixing of colloids interacting with reversible short-range attractions at semi-dilute volume fractions (ϕ = 2–10%) in real time by three-dimensional (3D) confocal microscopy observation (see SM B on the details). We prepare two types of colloids with different particle diameters (σ = 1.9 μm (EXP1) and 2.9 μm (EXP2)) but with almost the same interaction range with respect to the colloid size (Δ ~0.13). Matching of the intercolloid potential allows us to map the thermodynamic states in the two systems onto a single phase diagram (see Fig. S11b for phase diagram). By comparing the kinetics of colloidal phase separation between EXP1 and EXP2, we experimentally examine the scalability of the phenomena through the particle-size dependence (see SM B on the details of the samples and a special experimental protocol to initiate phase separation). Then, we further compare the experimental results with those of two types of simulations, BD and FPD20,50 (see Methods). We choose the system size L as L/σ = 34.6 and confirm that it is large enough to avoid finite-size effects in our study (see SM A.1). Since FPD can deal with many-body hydrodynamic interactions whereas our BD simulation completely neglects them, the comparison allows us to examine the impact of hydrodynamic interactions on phase separation.

We carefully match the intercolloid potential and the inverse temperature between the experiments and simulations using two fitting parameters. The detail is described in SM C. We will show below that by choosing the Brownian time, τ B as a unit time, the two sets of experimental results and the FPD simulation results show excellent agreement on phase demixing kinetics even on a quantitative level. We emphasize that we do not use any adjustable parameters in this comparison of the dynamic processes after matching the thermodynamic behavior between experiments and simulations.

Cluster formation in a dilute suspension

Now we turn our attention to the dynamics of phase separation taking place at two state points, point A (ϕ ≅ 2%, βε ≅ 7) (see Fig. 2) and point B (ϕ ≅ 10%, βε ≅ 6) (see Fig. 4), at which clusters and percolated network structures are formed at the end of observations respectively. In the following, we take an average over three independent phase-demixing processes in both experiments (EXP) and simulations (SIM) for gaining high statistical accuracy. Here we stress that the reproducibility was extremely high for both EXP and SIM (see SM B.3, Fig. S8 for EXP).

Fig. 2 Cluster-forming phase separation at state point A (ϕ ≅ 2%, εβ ≅ 7). We compare the time development of three-dimensional (3D) structures obtained from Brownian dynamics (BD) and fluid particle dynamics (FPD) simulations and confocal microscopy observation in EXP2 at t w /τ B ≃ 0, 50.4, 100.8 and 201.6 (see also Supplementary Movie 1). The particle color represents the number of particles in the cluster (see the color bars on the top right side); for example, clusters involving more than 10 particles appear red. The box sizes are commonly set as L/σ = 34.3 Full size image

First we discuss cluster-forming phase separation in a dilute suspension at point A (ϕ ≅ 2%, βε ≅ 7). To quantitatively characterize the coarsening dynamics during phase separation, we calculate the temporal change of the characteristic wave number of the structure factor S(q, t w ) as its first moment: \(\langle q(t_{\mathrm{w}})\rangle = {\int} dq\,qS(q,t_{\mathrm{w}})/{\int} dq\,S(q,t_{\mathrm{w}})\) (see SM B.4 on the details). The results are shown in Fig. 3a. First we can see that the two sets of experimental results with different sizes of colloids (EXP1 and EXP2) exhibit almost identical coarsening behavior after scaling the length and time by σ and τ B respectively. This implies that colloidal phase separation with a short-range attraction is scale invariant for the same set of (ϕ, βε, Δ), for which the equilibrium phase diagrams can also be scaled.51 This result is a direct experimental support for the scale invariance of the phenomena.

Fig. 3 Analysis of cluster-forming phase separation. a Temporal change of the characteristic wave number estimated from the structure factor, 〈q〉. b Temporal change of the fraction of clusters of size i, n i /n p . Gray curves are the theoretical predictions of Eq. (1).27 Five curves and different symbols correspond to i = 1–5 from top to bottom. In both (a, b), we compare experiments and simulations: EXP1 (black cross), EXP2 (black square), Brownian dynamics (BD; blue curve) and fluid particle dynamics (FPD; brown curve). c Dependence of the viscous drag coefficient for two particles approaching with each other, normalized by that of a free isolated particle, C 2 , on the interparticle distance r scaled by σ. Brown cross symbols are the results calculated by FPD method. The gray curve is the theoretical prediction.52 Here r is chosen as r = r s + σ H , where r s is the distance between the surfaces of the two colloids and σ H is the hydrodynamic radius of the colloids (σ H = 1.1σ). In the BD simulation, C 2 = 1 (blue line) since hydrodynamic interactions are ignored Full size image

Next we can see clearly that FPD almost perfectly reproduces the experimental data, but BD fails even qualitatively (see Supplementary Movies 1). This strongly indicates not only the numerical predictability of our FPD simulation but also the crucial role of the hydrodynamic degrees of freedom of the solvent, or hydrodynamic interactions, in the kinetic process of phase demixing of a colloidal suspension and the resulting colloidal aggregation. We can also see that 〈q〉 decreases with time t w as \(\langle q\rangle \propto t_{\mathrm{w}}^{ - \alpha }\) (α: the domain growth exponent) with α ≅ 0.33 for ϕ ≅ 2% for both EXP and FPD, but the exponent is considerably smaller for BD. The growth exponent α ~ 1/3 seen at ϕ ≅ 2% is characteristic of Smoluchowski’s Brownian coagulation mechanism26,27 (see the gray line in Fig. 3a), indicating that the phase demixing observed belongs to typical spinodal decomposition with droplet morphology.

In the above, we see the significance of hydrodynamic interactions in the macro-scale coarsening process of colloidal phase separation. Next, we focus on the microscopic process of aggregation to make a more rigorous comparison between experiments, simulations and theory. We regard particles as connected and belonging to the same cluster if the interparticle distance is less than (1 + Δ)σ. According to Smoluchowski’s Brownian coagulation theory,26,27 the number of clusters with size i at time t per unit volume, n i (t), can be given by

$$n_i(t)/n_{\mathrm{p}} = (n_{\mathrm{p}}Kt)^{i - 1}(1 + n_{\mathrm{p}}Kt)^{ - i - 1},$$ (1)

where n p (=n 1 (0)) is the total number of particles per unit volume and K is the so-called coagulation rate. Here we employ K as a fitting parameter to the experimental data. In Fig. 3b, we plot n i (t)/n p for different i(=1∼5) for EXP and SIM. We can see that FPD again well reproduces the experimental data of each type of cluster. By fitting Eq. (1) to the experimental data, we obtain K/K s = 0.50, 0.93, 1.2, 1.4 and 1.7, respectively for i = 1, 2, 3, 4 and 5 (gray lines in Fig. 3b). Here we use K s = 4k B T/3η, which is the collision rate derived by Smoluchowski, as a value to normalize K. Note that n i (t)/n p ∼ (n p Kt)−2 for a long-time limit and the height of the asymptotic behavior is proportional to K−2. The rate of collision between two equal-size spherical particles, which interact with the potential U(r), in a viscous liquid is given by27

$$\frac{K}{{K_{\mathrm{s}}}} = \left\{{{\int}_{\!\sigma} ^\infty d r\frac{\sigma }{{r^2}}C_2(r)e^{\beta U(r)}} \right\}^{ - 1},$$ (2)

where C 2 (r) = ζ 2 (r)/ζ 2 (∞) is the ratio between the viscous drag coefficient for two particles at the center-of-mass distance r approaching each other (ζ 2 (r)) to the one for a single spherical particle (ζ 2 (∞)). Estimating the collision rate with relation (2), we obtain K/K s = 0.56, which is reasonably close to the experimental result for i = 1 (K/K s = 0.50 (see above)). Collision rates for i ≥ 2 are higher than this value, likely reflecting the aspherical shape of clusters of i ≥ 2. Such details of cluster morphology and the resulting modification of hydrodynamic interactions make theoretical predictions of the collision rate extremely challenging.

Figure 3c also shows that there is a huge mismatch between BD and EXP and the aggregation speed is much faster for BD than for EXP. In BD, the collision rate for i = 1 is K/K s ~ 1, which is almost equivalent to the one for the case where hydrodynamic interactions are completely neglected (if we assume C 2 (r) = 1 in Eq. (2), we get K/K s = 1.1). These results clearly indicate that collisions between colloids are strongly hindered by hydrodynamic interactions. In Fig. 3c, we plot the functional form of C 2 (r) obtained by FPD together with the theoretical prediction,52 which shows an excellent agreement. We can see that there is a steep increase in the friction coefficient between the two particles (or, the “strength” of hydrodynamic interactions between them) for \(r/\sigma \lesssim 2\). This tells us that hydrodynamic effects at such a “short-range” distance, or “lubrication effects”, are responsible for the significant suppression of the collision rate. We note that the results of BD can never be mapped on those of EXP (FPD) by rescaling time by other time units, which is obvious from the difference in the growth exponent between BD and EXP (FPD).

Network-forming phase separation in a semi-dilute suspension

As the volume fraction of colloids increases, it is expected that the many-body nature of hydrodynamic interactions plays a more important role. Now we focus on network-forming phase separation in a semi-dilute suspension (see Fig. 4). To check how precisely FPD can capture such effects, we make detailed comparison between EXP and SIM at the volume fraction ϕ ≅ 10% (point B), where space-spanning network structures were formed. This formation of space-spanning network structures of the minority colloid-rich phase in phase separation is characteristic of viscoelastic phase separation.38,39,53 It is worth stressing here that all the quantities we measured (see Fig. 5a–d) are almost the same between EXP1 and EXP2, which is a direct experimental support of the scale invariance of the phenomena.

Fig. 4 Network-forming phase separation at state point B (ϕ ≅ 10%, εβ ≅ 6). We compare the time development of three-dimensional (3D) structures obtained from Brownian dynamics (BD) and fluid particle dynamics (FPD) simulations and confocal microscopy observation in EXP2 at t w /τ B ≃ 0, 30, 60 and 120 (see also Supplementary Movie 2). A clear network structure is already formed at t w /τ B = 30 for BD, but not for EXP2 and FPD. Particles are colored to distinguish front particles from back ones Full size image

Fig. 5 Structural analysis of network-forming phase separation. We compare the time evolution of structures between experiments and simulations, EXP1 (black cross), EXP2 (black square), Brownian dynamics (BD; blue curve) and fluid particle dynamics (FPD; brown curve): the characteristic wave number a, the number of bonds per particle b, the total length of the skeleton l c and the Genus number of the skeleton d. We made network analyses for the box whose side length is L/σ = 34.6 in common for experiments and simulations Full size image

Considering the hierarchical nature of network structures, we perform structural analysis at various length scales. First we examine the temporal change in the characteristic wave number 〈q〉 during phase separation. The results are shown in Fig. 5a. We can see that FPD almost perfectly reproduces the experimental data, but BD fails even qualitatively (see Supplementary Movies 2). This strongly indicates again the numerical predictability of our FPD method and the crucial role of the hydrodynamic degrees of freedom of the solvent and the resulting many-body hydrodynamic interactions in the kinetic process of phase demixing of a colloidal suspension. We emphasize that the coarsening behavior obtained by FPD simulations is reliable and does not suffer from a smallness of the Schmidt number S c . We have checked the S c dependence of the simulation results and confirmed that for S c ≥8.0, there is no dependence of the phase separation kinetics on S c (see SM A.2.3, Figs. S4 and S5).

We discuss the coarsening dynamics observed in the wave number space in more detail. In SM B.4, Fig. S9a, we analyze the early stage of phase separation in an experiment and confirm that in the initial stage the scattering intensity grows exponentially with time for every wave number, which are characteristic of the Cahn’s linear regime of spinodal decomposition in a thermodynamically unstable state.36 In SM B.4, Fig. S9b, we also show the time evolution of S(q, t w ) in the experiment, where we can see that the scattering peak appears in the very early stage, then its position continuously shifts towards a lower wave number with time. Consistently, we can see in Fig. 5a that, for both EXP and FPD, 〈q〉 is almost constant with time in the initial stage but then starts to decrease by obeying a power law \(\langle q\rangle \propto t_{\mathrm{W}}^{ - \alpha }\) (α ≅ 0.46). This domain growth exponent α is different from the well-known exponents for ordinary spinodal decomposition in a binary critical mixture, α = 1/3 for the Brownian coagulation mechanism and α = 1 for Siggia’s hydrodynamic coarsening mechanism.36 The physical meaning of α ~0.5 is not clear at this moment, but we note that this exponent has often been observed for systems with strong dynamic asymmetry, e.g., gas–liquid phase separation in molecular systems54,55,56 and viscoelastic phase separation of colloidal suspensions,22,57,58,59 protein solutions58 and lyotropic liquid crystals.60 Here we emphasize that network formation at ϕ ≅ 10% clearly indicates that phase separation of colloidal suspensions does not belong the dynamical universality class of binary mixtures of simple liquids (model H36) which predicts droplet formation at this ϕ, but belong to viscoelastic phase separation.38,39,53 We also note that 〈q〉 deviates from the power law in the late stage. This behavior may reflect the dynamic arrest due to gelation, but we do not focus on this trend here since it is out of scope of this work.

Next we study the aggregation process at the particle scale. Figure 5b shows the temporal change of the number of interparticle bonds per particle, n b , where the bond definition is the same as the one for clusters. The time development of n b is almost identical between EXP and FPD, but again shows much faster increase, or faster contraction, for BD (see Fig. 4 and Supplementary Movie 2). We can also see that in the entire time range of the observation n b < 7 (note that n b ~12 for a compact structure in 3D), indicating that most particles are located in the domain interface. This lack of separation between the domain size and the interface thickness makes a coarse-grained theoretical description of phase separation quite difficult.

Now we turn our attention to the topological characteristics of network structures in real space. For further characterization of the network structures, we coarse-grain a phase-separated structure by operating Gaussian blurring of each particle, extract the network domain by binarization and then extract the skeleton of the network structure by applying a skeletonization procedure (see SM D on the details). For skeleton structures obtained by this procedure, we first calculate the total length of the skeleton, l. Figure 5c shows the temporal change of l normalized by the size of the box L. To capture the topological feature of a network, we also calculated the Genus number G (Fig. 5d), which is equivalent to the total number of holes (or pores) existing in the network structure, by analyzing the connectivity of the skeleton (see SM D on the details). All these results show excellent agreements between EXP and FPD. On the other hand, the results of BD exhibit much faster changes than EXP and FPD, again reflecting the lack of hydrodynamic lubrication effects as in the dilute system (at point A). In the late stage of the inset of Fig. 5d, we can see that the G value of BD is approximately half of EXP and FPD. This clearly indicates that the presence or absence of hydrodynamic interactions has a significant influence on the kinetic pathway and the resulting topological characteristics of the final network structures. We stress that it is these topological features that control the physical properties of network materials, such as elastic, transport and surface properties.

To summarize, we systematically study the phase separation process of colloidal suspensions by comparing single-particle-level time-resolved 3D confocal microscopy observation, BD simulations without hydrodynamic degrees of freedom and FPD simulations with full many-body hydrodynamic interactions. We find that the coarsening dynamics of cluster-forming and network-forming phase separation in the experiments using colloids with two different sizes can be almost perfectly reproduced by FPD simulations without any adjustable parameters after careful matching of the intercolloid potential and the temperature (or, thermal noise), but cannot at all by BD simulations. This finding demonstrates the fundamental importance of many-body hydrodynamic interactions in the dynamical structural formation of colloidal suspensions. More importantly, it indicates that simulation methods based on direct computation of the Navier–Stokes equation including FPD simulation has a high predictive power for nonequilibrium processes in colloidal suspensions, which may significantly contribute to not only the basic physical understanding of these phenomena but also the computer-aided design of colloidal materials.