Monodisperse suspensions of Brownian colloidal spheres crystallize at high densities, and ordering under shear has been observed at densities below the crystallization threshold. We perform large-scale simulations of a model suspension containing over $105$ particles to quantitatively study the ordering under shear and to investigate its link to the rheological properties of the suspension. We find that at high rates, for $Pe>1$, the shear flow induces an ordering transition that significantly decreases the measured viscosity. This ordering is analyzed in terms of the development of layering and planar order, and we determine that particles are packed into hexagonal crystal layers (with numerous defects) that slide past each other. By computing local $\psi 6$ and $\psi 4$ order parameters, we determine that the defects correspond to chains of particles in a squarelike lattice. We compute the individual particle contributions to the stress tensor and discover that the largest contributors to the shear stress are primarily located in these lower density, defect regions. The defect structure enables the formation of compressed chains of particles to resist the shear, but these chains are transient and short-lived. The inclusion of a contact friction force allows the stress-bearing structures to grow into a system-spanning network, thereby disrupting the order and drastically increasing the suspension viscosity.

## I. INTRODUCTION

Suspensions are ubiquitous in nature and play an important role in a wide variety of environmental and technical processes. Examples include concrete, pharmaceuticals, blood, and more. Dilute suspensions typically exhibit liquidlike behavior, but when the volume fraction of suspended particles is high, complex rheological properties arise [1–4]. Understanding and controlling this behavior is of academic interest and crucial for many industrial applications.

Colloidal suspensions are known to order under shear, giving rise to shear thinning properties [5–12]. Dense hard sphere suspensions will also order at rest due to thermal motion, given enough time, forming face-centered cubic (fcc) or hexagonal close-packed (hcp) lattices [13–15]. The dynamics of this crystallization process are complex, depending strongly on the volume fraction of particles, particle size distribution, boundary conditions, and the applied shear rate [12,15–18].

In addition, at high densities, many suspensions exhibit shear thickening properties [1,4,19–21]. This is generally associated with cluster formation and contact forces between particles. At higher shear rates, more frequent contact between particles gives rise to large microscopic forces and a high macroscopic viscosity. This rise in viscosity is called continuous shear thickening (CST) or discontinuous shear thickening (DST), depending on how rapidly the viscosity rises with shear rate. There is substantial evidence that frictional forces, due to particle roughness, play a key role in the shear thickening phenomenon, with DST being attributed to a sharp rise in the number of frictional contacts [1,2,4,22–24].

Using a model monodisperse suspension of hard spheres, we investigated the ordering transition for a system in the liquid-crystal coexistence region of the equilibrium phase diagram. We analyzed the types of ordered states that arise under shear, as well as how the ordered domains nucleate and grow. By computing the distribution of stresses within the suspension, we discovered that the largest contributions to the shear stress come from the defects in the ordered state. The defect structure enables the formation of compressed chains resisting the flow. While these chains are short-lived in the frictionless system, the addition of contact friction between particles disrupts the order and allows the stress-carrying structures to grow into a system-spanning network.

This paper is organized as follows. First, we give a detailed introduction to the model and describe the simulation and analysis performed in Sec. II. Next, we discuss the ordering transition the sheared suspension undergoes and its effect on viscosity in Sec. III A. In Sec. III B, we quantitatively analyze the ordering and discuss how it develops over time, as the suspension is sheared. The formation of defects and analysis of their structure are explained in Sec. III C. The microstructure analysis is then related to the mechanics by analyzing the local particle stresses in Sec. III D. Finally, we comment on the role of friction in the disruption of ordering and redistribution of stresses in Sec. III E.

## II. METHODS

### A. Simulation details

We performed discrete element method (DEM) simulations of a suspension of monodisperse, spherical particles, as shown in Fig. 1. The suspension contained 108 000 particles in a periodic simulation box with the particle volume fraction ranging from $\varphi =54%$ to $58%$, with most of the analysis presented being focused on the $\varphi =54%$ case. This is a range of $\varphi $ where we observe substantial ordering and which spans CST and DST regimes as observed in simulations with friction [25,26], while below the static jamming volume fraction, the suspensions can shear jam at high shear rates when $\varphi \u227356%$ [25,26]. The solvent was treated implicitly, with a viscosity of $\eta 0=1$ in reduced units of $m\epsilon /d$. The particle mass $m$ and diameter $d$, as well as the simulation energy scale $\epsilon $, were also set to 1. Note that the particles have finite mass because we are solving Newton’s equations of motion, which include inertia. However, we have ensured we are in the overdamped limit by running simulations with suspension viscosity increased up to $\eta 0=10$ and finding that its value does not affect our results.

All simulations were performed with LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) [27]. While generally known for molecular dynamics simulations, the LAMMPS code has specialized modules that allow for the modeling of hard sphere suspensions of Brownian particles. For this application, LAMMPS utilizes a simple physics-based DEM model that ignores the detailed flow behavior of the suspension solvent in exchange for computational efficiency. Hydrodynamic interactions between spheres are largely controlled by lubrication forces, as discussed in Sec. II B. Shear flow was imposed along the $x$ direction (with gradient along $z$) based on the Lees–Edwards boundary condition, with an additional Stokesian drag force that causes particles to follow the imposed shear profile over time. The robustness of this approach improves with increasing volume fraction of spheres as a less detailed knowledge of flow of the background fluid is needed. Indeed, it has been found that at volume fractions of approximately $40%$ and higher, the flows produced are reasonably consistent with fully detailed simulations [28]. This is due to the fact that, at higher volume fractions, the surfaces of the solid inclusions are close enough such that the lubrication forces dominate over the long-range hydrodynamics of the background fluid [29].

In addition to the lubrication forces, a short-ranged Yukawa potential was used to give the particles an exponentially decaying steric repulsion following $UYukawa=(A/\kappa )e\u2212\kappa h$. Here, $h$ is the separation between particle surfaces, while $A=1000\epsilon /d$ and $\kappa =100d\u22121$ are set so that the energy is $10\epsilon $ at contact and falls below $kBT$ at a surface separation of $h=0.02d$. Thermal fluctuations were introduced using a random Brownian force to maintain a finite temperature such that $kBT=\epsilon $. This sets the time scale in the simulations as $\tau =md2/kBT$. We tested the model parameters for a slightly more dilute suspension ($\varphi =45%$) and found results in excellent agreement with the earlier work of Foss and Brady where they performed simulations of hard spheres with hydrodynamic interactions [5].

All shear simulations were run from a disordered initial state. This state was prepared by first placing the particles on a cubic lattice and then running a high temperature DEM simulation without any shear. Once the crystal had melted, the system was cooled to $kBT=1\epsilon $, and the pair correlation was checked to ensure the configuration was amorphous. For all rates considered, a copy of this initial state was sheared until the system reached steady state, as evidenced by steady viscosity, pressure, and energy. To ensure that we remained in the viscous regime rather than inertial (see [30] for more on this transition), we performed additional simulations with $\eta 0=10$ and found that the suspension relative viscosity was unchanged.

### B. Lubrication forces

Hydrodynamic lubrication forces arise due to the confinement of a solvent between moving particles [31]. The interactions due to the relative motion of the spheres have been described by Ball and Melrose [29], who computed forces ($F\u2192$) and torques ($\tau \u2192$) acting on spherical particles as

where the forces and torques denoted by $sq$, $sh$, and $pu$ are due to squeezing, shearing, and pumping, respectively, of the fluid between the spheres (see Fig. 2 for a visual depiction of these modes). The velocity of the spheres relative to each other is split into normal and tangential components as $vn\u2192$ and $vt\u2192$, where the normal is defined as the direction along the center-centerline between the spheres. $\omega t\u2192$ is the angular velocity difference in directions tangent to the normal direction. The coefficients for Eqs. (1)–(4), as computed by Ball and Melrose [29], are

where $r$ is the radius of the spheres, $h$ is the surface-surface separation between spheres, and $\eta 0$ is the viscosity of the solvent.

As one can see from the coefficients, these forces are singular at contact ($h=0$), and this singularity is typically handled in simulations by setting an inner cutoff $hin$. This sets a finite minimum separation such that the lubrication forces use $hin$ in lieu of $h$ when $h<hin$. We tested the dependence of the stress response on this threshold and found that for values of $hin>10\u22124d$, the hydrodynamic stresses were underestimated at high rates. The exact threshold is, of course, dependent on the presence of other particle interactions, and a different cutoff can be appropriate if there are other forces, such as electrostatic repulsion, between particles. For all data presented, we used $hin=10\u22124d$.

### C. Microstructural analysis

Several analytic tools were used to quantitatively study the ordering of particles in the suspension. We computed the density profiles $\rho (z)$ and found evidence of particle layering. Steinhardt’s bond orientational order parameters were used to characterize the 3D ordering [32], but no 3D crystal structure was found (see the supplementary material [33]). Due to the nature of the order, we focus on 2D correlation functions and order parameters that are computed in layers of particles.

Generally, all of the results that are presented as functions of $z$, such as the density profiles $\rho (z)$ in Fig. 4, are computed by averaging the plotted property over all particles in a slice of the system along the gradient direction. With the density profile as reference, we use a slice width of $2d$ when averaging the correlation functions and order parameters (i.e., the average is over all particles $i$ such that $\Delta z=|zi\u2212z|\u2264d$). Particles in the same layer are identified by a threshold separation of $\Delta z\u22640.4d$, which corresponds to the width of the density peaks in Fig. 4. To clarify this, let us consider the case of the correlation function computed as

Here, the Heaviside step function $H$ is used to count particles $j$ within a certain distance $rxy$ of particle $i$, with discrete bin size $rbin=0.05$. The quantity is averaged over all $Nz$ particles within the slice of width $2d$ for better statistics, and the calculation is confined to $ij$ pairs within the same layer.

A similar procedure is followed with the order parameters $\psi 6$ and $\psi 4$. The individual particle values are computed, following [34], as

The sum is limited to the $Nn$ nearest neighbors of particle $a$ in the same layer, and $\theta ab$ measures the angle of the bond between $a$ and $b$. To obtain $\psi 6(z)$, the individual particle $\psi 6a$ values are averaged for all particles in a slice of width $2d$ centered at $z$.

In addition to analysis of the microstructure, we calculate the stress in the suspension to connect the ordering to the mechanics. The stress is calculated from the particle contributions to the stress tensor, following [35], as

where $i$ and $j$ denote particles, $\alpha $ and $\beta $ denote Cartesian directions, $V$ is the system volume, $Fij$ is the force between two particles, $rij$ is the distance between two particles, and $v$ is the particle velocity. The suspension viscosity is computed from the shear stress $\sigma xz$ and shear rate $\gamma \u02d9$ as $\eta =\sigma xz/\gamma \u02d9$. In Sec. III D, we discuss how the total stress can be divided into particle stresses to give local information.

### D. Friction

When including frictional forces, we follow the contact model of Mari *et al*. [25,26]. The steric repulsion between particles is modeled as a Hookean force, with a normal force of $FN=kh$. Unlike the Yukawa-like potential used previously, this allows for particle contact (which is the criterion for activating frictional forces). However, to maintain a hard contact, we use a high spring constant of $k=104\gamma \u02d9mkBT$. With these parameters, we observe particle overlap $<6%$, comparable to the criteria used in other recent simulations [24]. The dependence on shear rate arises from the fact that a higher spring constant is needed to limit overlap at higher shear rates, due to collisions happening more frequently and at greater velocities. Alternatively, one could use a fixed, high value of $k$, corresponding to what is required for the highest shear rate, at all rates. However, this requires a smaller time step to adequately resolve the collisions, and the variable spring constant approach was found to reduce computational cost while producing results that matched the results from the constant $k$ approach [24].

The frictional forces act tangentially to the particle contact and also follow a simple Hookean model. Any tangential displacement of particles after making contact, $\Delta rt$, is acted upon by a restoring frictional force $Ft=k\Delta rt$, with the constraint that $Ft\u2264FN$ (Coulomb’s friction law with a friction coefficient of $1$). The contact/friction model we use is relatively simple compared to some newer studies [36,37]. However, all these models exhibit qualitatively similar behavior, and the simpler model also compares favorably with experiments [38]. In this study, we focus mainly on the differences between simulations with/without friction, rather than the quantitative details related to the exact frictional model chosen.

## III. RESULTS

All results are presented in reduced units, meaning distance is in units of particle diameter $d$, mass in units of particle mass $m$, and energy in units of $\epsilon =kBT$. The solvent is treated as a Newtonian fluid with viscosity $\eta 0=1$ in reduced units of $m\epsilon /d$, and all results on viscosity refer to the relative viscosity of the suspension compared to the solvent. Shear flow is applied in the $x$ direction, with gradient along $z$, and the results on 2D correlations or order parameters are in the flow-vorticity ($xy$) plane.

### A. Shear-induced ordering

We investigate the steady shear response of the model suspension of monodisperse colloidal spheres just described at volume fraction $\varphi =54%$. At this density, monodisperse hard sphere suspensions are in the liquid/crystal coexistence region of the phase diagram [13,15]. However, our initial configuration is obtained by quenching the suspension to a metastable glassy state that does not exhibit any crystallization. We subject replicas of this configuration to steady shear at seven different shear rates (ranging from $10\u22122\tau \u22121$ to $10\tau \u22121$). The Peclet number $Pe$, defined as the ratio of motion due to shear compared to motion due to Brownian forces, is $Pe=6\pi \eta 0r3\gamma \u02d9/(kBT)=(3/4)\pi \gamma \u02d9$ for the parameters chosen in our simulations.

As the simulation time is too short for crystallization to occur in the stationary suspension, the low $Pe$ flow exhibits a microstructure and viscosity that are independent of the elapsed time (and applied strain). The shear stress, and thus the apparent viscosity, fluctuate due to the Brownian motion of the suspended particles, but the steady state is reached almost immediately. As shear rate is increased, we observe shear thinning behavior and the viscosity drops significantly as seen in Fig. 3(a). In addition, for $Pe>1$, there is a substantial difference in the initial viscosity and the steady state viscosity.

For $\gamma \u02d9=1\tau \u22121$ and $\gamma \u02d9=10\tau \u22121$, the viscosity at the start of shear is about $\eta /\eta 0=45$, but as the applied strain increases, the viscosity drops toward a final steady state value of $\eta /\eta 0\u224310$ [Fig. 3(b)]. System snapshots are displayed in Fig. 3(c) to show that this drop in $\eta /\eta 0$ is coupled to microstructural changes in the suspension. The pictured differences in local density are coupled to the formation of an ordered phase: the applied shear, at $Pe>1$, promotes a switch from the metastable glassy state to an ordered state. Signs of a similar transition have been seen in a previous study even with smaller applied strains [39]. Surprisingly, the velocity profile remains linear, despite the apparent banding in the local volume fraction.

### B. Analysis of ordering

We note that the ordering observed here is different from the equilibrium crystallization of similar suspensions at rest. As the system is sheared, the particles assemble into a layered configuration, with layers in the flow-vorticity ($xy$) plane spaced along the gradient direction ($z$). In Fig. 4, we show a plot of the local density $\rho (z)$. The density profile is flat at early times (low strains) but develops an oscillatory pattern at large $\gamma $. This shows the formation of particle layers, corresponding to the peaks in density, with center-to-center distance between the layers of $\u22430.9d$. The shear flow imposes the requirement that these layers slide past each other, which disrupts the typical fcc or hcp crystal structure in the equilibrium crystals [13,15]. Thus, the ordered state we observe consists of domains of stacked, 2D crystals coexisting with domains of lower density. This is similar to what is seen in earlier experiments and simulations with oscillatory shear when going to high amplitudes [6,9,40,41]. To better characterize this microstructure, we focus our quantitative analysis on the structure within these layers and compute 2D correlation functions and order parameters.

We first calculate a 2D version of the pair correlation function $g(rxy)$ in the flow-vorticity plane. This is computed within the previously described layers of particles, at various $z$ positions along the gradient direction, using the $xy$ distance between those particles. For the purposes of this calculation, particles $i$ and $j$ are considered to be in the same layer if $|zi\u2212zj|<0.45d$, which is half the peak-peak distance in the oscillatory density profile seen in Fig. 4. The results of this calculation are shown in Fig. 5, with plots of $g(rxy)$ for various values of $\gamma $ and $\gamma \u02d9$.

In Fig. 5(a), we show the liquidlike ordering observed at lower rates (which is independent of $\gamma $ for the strains considered). There is a clear first neighbor peak, but the correlations rapidly decay beyond that. This also matches what we obtain from our initial state corresponding to a glassy, amorphous state. Note also that all the curves for different $z$ values lie on top of each other—the lack of ordering is consistent for all the layers along the gradient direction.

Next, we consider the scenario at a higher rate ($\gamma \u02d9=1\tau \u22121$). Here, the system transitions to a much lower viscosity state when sheared, accompanied by an increased layering of particles and the formation of domains with different local densities. In Fig. 5(b), we see that this substantially changes the 2D ordering of the particles. In particular, particle positions within layers are more strongly correlated and remain correlated over much larger distances. Moreover, the curves are no longer $z$-independent: the local structure varies across the system. At the highest shear rate of $\gamma \u02d9=10\tau \u22121$, the peaks of $g(rxy)$ remain in similar positions and heights, but the variation with $z$ increases further [see Fig. 5(c)].

To better understand the changes in the local ordering, both in space and over time, we compute the 2D order parameter $\psi 6$ in the flow-vorticity plane. Figure 6 shows a snapshot of particles in the dense region of the sheared suspension at $\gamma \u02d9=10\tau \u22121$ arranged in a hexagonal lattice, so the $\psi 6$ order parameter (which quantifies the degree of hexagonal order and goes to 1 in a perfect hexagonal lattice) is a natural choice to measure the local order [34]. Like for the 2D pair correlations, we compute this parameter in the $xy$ plane for particles within a layer. $\psi 6i$ is computed following Eq. (9).

The particle order parameter $\psi 6i$ can be averaged over the layer to obtain $\psi 6(z)$. This is the quantity we plot in Fig. 6 for $\gamma \u02d9=10\tau \u22121$ and different values of $\gamma $. By doing so, we see that the order parameter is uniformly $0$ in the disordered state at early times. As the strain accumulates, parts of the system begin to exhibit a sharp increase in $\psi 6$, indicating the nucleation and growth of ordered domains. The position of these ordered domains matches the higher local density regions seen in the snapshots of Fig. 3, suggesting the ordering allows for more efficient packing of the particles. Also, as we saw for the density in those snapshots, $\psi 6$ is highly heterogeneous when the structure first begins to change (and the viscosity begins to drop), but it becomes more uniform in the final steady state. The high amount of hexagonal order we observe in the steady state is in strong agreement with previous measurements made via scattering experiments [6,9,41], which reveal the same type of order in colloidal silica suspensions. A threshold of $\psi 6>0.6$ was used to determine the fraction of ordered particles in Fig. 3(a), and it is clear that this is much higher for $Pe>1$. Interestingly, as the shear rate is increased beyond that, there seems to be a small decrease in the fraction of particles exhibiting order, consistent with the disruption of order at high enough rates as seen by Ackerson [6].

### C. Formation of defects

In steady state, with increased ordering and decreased viscosity, the relatively low average value of $\psi 6$ obscures the fact that the vast majority of particles have $\psi 6\u22431$. In Fig. 7, we show a snapshot for $\gamma \u02d9=1\tau \u22121$ where each particle is colored by its $\psi 6$ value. Here, it is clear that there are two populations: the majority of particles exhibit nearly perfect hexagonal order ($\psi 6\u22431$), while there still exists a non-negligible minority that form domains of low $\psi 6$. This second population has a tendency to form chainlike structures where the hexagonal lattice is disrupted, and the particles seem to arrange in a more squarelike lattice.

Evidence for chain structures, combined with hexagonal structures, has been seen before in the Stokesian dynamics simulations of Xu *et al*. [42]. However, given their limited system size (<200 particles), it was unfeasible to monitor the formation of these large-scale defects. Interestingly, similar defects have also been seen experimentally in particles assembled on a curved surface [43]. In that case, line defects, or “pleats,” formed as a way to handle the curvature-induced geometric frustration of the hexagonal packing. It is possible that we are observing a similar phenomenon but with stresses due to shear flow instead of curvature. If the squarelike clusters are as common as this snapshot suggests, it should be possible to identify them with a different order parameter, $\psi 4$.

$\psi 4$, defined in Eq. (10), selects for fourfold symmetry in the $xy$ plane, with nearest neighbors at 90$\xb0$ angles to each other. To suppress small fluctuations in the order parameter and highlight only larger groups of particles, we now average the particle order parameters over the nearest neighbor shells. This coarse-grained order parameter is shown in Fig. 8. Here, the two snapshots correspond to the exact same configuration, colored by $\psi 6$ (left) and $\psi 4$ (right). By comparing these two, it becomes clear that the defects in the hexagonal crystal correspond to elongated domains of enhanced $\psi 4$.

Crucially, these domains do not encompass entire layers of particles, which makes them difficult to detect via the quantitative layer analysis we performed earlier. However, the variations in the ordering are coupled to spatial fluctuations in the density. In the third snapshot of Fig. 8, we show the same configuration characterized by the local volume fraction $\varphi =NlVp/Vl$, where $Nl$ is the number of particles within five particle diameters, $Vp$ is the volume of a particle, and $Vl=(4/3)\pi 53$. The larger regions of low $\psi 6$/high $\psi 4$ coincide with lower local density. Hexagonal order allows for denser packing of the spherical particles, but due to the shear flow a fully ordered state is not formed. Instead, we observe chainlike defects characterized by the $\psi 4$ order parameter and lower local density. As seen in Fig. 3(a), ordering is correlated with a substantial decrease in viscosity at $Pe>1$, but further increasing the shear rate decreases both the number of ordered particles and the viscosity. This implies that these elongated defects actually help us to sustain steady, unimpeded flow.

### D. Impact on stresses

Recognizing that all these microstructural changes (fluctuations in density, enhanced layering/ordering) are accompanied by a substantial drop in the viscosity of the suspension, we calculate the particle contributions to the stress tensor to quantify the link between the order/density and stresses [35]. The virial formulation of the stress tensor is broken up into particle contributions as

where $\alpha $ and $\beta $ can be $x$, $y$, or $z$ to generate the components of the stress tensor, $F\u2192ij$ and $r\u2192ij$ are the force and position vectors between particles $i$ and $j$, and $v\u2192i=v\u2192totali\u2212\gamma \u02d9zx^$ is the deviation of the particle velocity from the flow profile set by the shear rate. The total stress tensor is obtained by summing $\sigma i$ for all particles as in Eq. (11). The distribution of particle contributions to the shear stress is plotted in the inset in Fig. 9(a), and it shows a peak at the system stress of $\sigma xz\u2243100kBT/d3$ and long tails with a small fraction of particles contributing stresses that are very large in magnitude. These large values come from particles that instantaneously experience a strong force. Due to the computation of these stresses at a particle level, and a large system size, we observe these extremes that are averaged out when considering the macroscopic stress. This raises the question of which particles are responsible for these large stresses.

Computing the local volume fractions in the neighborhood of each particle as before, we now analyze their distribution in Fig. 9(a). The blue curve shows a bimodal distribution, with the larger peak around $\varphi \u22430.58$ and a smaller peak at $\varphi \u22430.49$. If we instead restrict our analysis to the high stress particles, selected by using a stress threshold of $\sigma xzi>\sigma t=15000kBT/d3$ corresponding to $1%$ of the particles, we see a single peak at the lower density. In other words, the high stress particles are almost exclusively in low density regions, i.e., the defects. This is represented visually in Fig. 9(b), where only the high stress particles are rendered. Notably, the high stress particles are almost all blue, corresponding to lower local volume fraction. In addition, they tend to form short chains in the compressive direction, which have been observed in previous simulations [44] and are highly reminiscent of the force chains in sheared granular materials [45]. In our simulations, these chains form and break up rapidly as particles rearrange due to Brownian fluctuations and shear motion. Despite their transience, the presence of these structures shows a possible link between dense suspensions and granular systems even in the absence of friction.

At first, it seems somewhat counterintuitive that higher stresses are associated with lower local density, but this reveals the tendency of the locally stressed structures to form in regions of defects. Yet, at a macroscopic level, we see that the presence of defects is driven by the applied shear rate [see the fraction of ordered particles decreasing at $\gamma \u02d9>1\tau \u22121$ in Fig. 3(a)]. Does stress cause defects or do defects cause stress? A possible explanation draws on previous work in 2D systems. The experiments of Irvine *et al*., for example, showed that hexagonally packed particles on curved surfaces formed chainlike defects as an optimal way to relax the curvature-induced stresses [43]. Additionally, the simulations of Schwenke *et al*. showed that nanoparticles adsorbed to an interface formed a densely packed ordered phase with defects where adsorption could continue [46]. Similarly, we observe growth of chainlike defects as we increase the macroscopic stress by increasing $\gamma \u02d9$. These defects may enable the larger ordered domains to slide past each other with minimal resistance. A similar effect may be taking place in suspensions under oscillatory shear which undergo a change in crystal orientation when the strain amplitude is increased [40]. While they also allow for the formation of the compressed force chains seen in Fig. 9(b), those chains are transient as the contacts are lubricated and frictionless. In this way, the defects are an efficient response to the deformation-induced stresses, and their growth allows for a net decrease in the effective viscosity of the suspension.

### E. Frictional forces disrupt order

Thus far, all our results have been for entirely frictionless particles, but we have also investigated how the presence of friction might affect our findings. We use a frictional model based on the work of Mari *et al*. [25,26], where a springlike force acts against the rotation of particles in contact. A key parameter of this model is the friction coefficient, which we set equal to 1 so our simulations are in the regime where shear thickening has been observed [25,26].

To test the effect of frictional forces on the ordering, we perform simulations with friction starting from the sheared ordered states discussed in Secs. III A–III D. In Fig. 10, we show results from simulations at four shear rates ($\gamma \u02d9=0.01\tau \u22121$, $0.1\tau \u22121$, $1\tau \u22121$, and $10\tau \u22121$) where the initial state has been sheared without frictional forces, meaning that for the high rate cases ($\gamma \u02d9\u22651\tau \u22121$), the system exhibits the layering and $\psi 6$ ordering discussed previously.

In Fig. 10(a), we plot the shear stress as a function of applied shear strain. The first half shows the decrease in stress due to the ordering transition for frictionless spheres. After the dashed line, frictional forces are activated and there is a rapid increase in stress to a new plateau. This rise is correlated with the disruption of the ordered states. This disruption is very different from the defect regions discussed in Sec. III C, which were still found to form layers. Figure 10(b) shows a snapshot from simulations at $\gamma \u02d9=10\tau \u22121$, and the layered structure seen in the frictionless states is not preserved at all.

In terms of the rheology, the addition of frictional forces dramatically increases the viscosity of the suspension and causes it to steadily shear thicken, with viscosity increasing by a factor of 2 over a decade increase in shear rate [see Fig. 10(c)]. With respect to the frictionless case, therefore, the rise in viscosity is associated with two factors: the extra stress from frictional forces and the disruption of the low-viscosity ordering. This behavior is in stark contrast with the thinning we observed without friction (see Fig. 3) and in agreement with the results of previous studies, which indicate that friction is the key to shear thickening [1,4,22–26].

Analysis of the particle contributions to stresses in our system shows that there is a fundamental change in which particles contribute to the increased shear stress. Mirroring the analysis in Sec. III D, we show the particle histogram of local $\varphi $ in Fig. 11(a). The presence of a single peak indicates that we no longer have two coexisting states (order/defect) as in the frictionless simulations. Additionally, the particle contributions to the shear stress are distributed normally, with a mean of $\sigma xz\u22432400kBT/d3$, and using the same criteria for high stress particles, $\sigma xzi>\sigma t=15000kBT/d3$, now selects $\u224335%$ of the particles instead of just $1%$. Interestingly, the histogram for the local volume fraction of these particles matches that for all particles, indicating that the link between local volume fraction and local stress is much less obvious here.

Visualizing the contacts between particles as bonds [Fig. 11(b)], we find that the high stress particles form a percolating network that spans both lower (blue) and higher (red) $\varphi $ regions. The ability of stresses to accumulate in a large-scale structure, instead of only in highly localized regions of defects, explains the substantial increase in the suspension viscosity. The existence of a percolated stress/contact network matches previous simulation studies of frictional systems [24–26].

We find that the inclusion of friction has not only changed the microstructure of the suspension, it has also changed how stress is carried within that structure. Instead of the majority of the stress response arising from small regions of defects, stress is distributed in a space-spanning structure, akin to the percolated frictional contact networks that are necessary for DST [47]. The formation of this larger-scale, percolating network structure raises some interesting new lines of inquiry. For one, studying the heterogeneity of the network, in terms of structure and local stresses, could provide new insight into the shear thickening behavior and the stress fluctuations seen in experiments [19,48]. Similarly, the dynamic evolution of the network under shear is another avenue to explore in future work. A better understanding of the dynamics would assist in efforts to tune the thickening, which might be achieved through the application of vibrational forces as done by Sehgal *et al*. [49].

## IV. CONCLUSIONS

We performed simulations of a dense colloidal suspension undergoing shear. The simulation approach included hydrodynamic lubrication forces and short-ranged repulsive interactions. The use of LAMMPS allowed us to perform large-scale simulations with $105$ particles in each sample. At high shear rates, $Pe>1$, we find that the shear flow induces an ordering transition in the suspension that led to a significant reduction in the overall viscosity. Shear-induced crystallization has also been seen previously in experiments [6,9,12] and simulations [8,10,11] of colloidal suspensions, as well as in polymer solutions [50]. By computing 2D pair correlations and the $\psi 6$ order parameter, we have determined that the ordered state consisted of hexagonally packed layers that slide past each other as the system is sheared in the steady state. Due to our use of a large system size, we were able to observe the nucleation and growth of large ordered domains, which coexist with lower density defect regions. The defects typically formed chainlike structures along the flow direction and can be identified by the $\psi 4$ order parameter. Surprisingly, these lower density regions actually contribute the most to the overall stress response due to the formation of short, transient chains of particles resisting compression. These “stressed chains” are dynamic, rapidly forming, and breaking apart under shear but are reminiscent of the force chains observed in jammed granular systems. The fact that these are mainly restricted to the defected regions could explain why the overall shear stress decreases when the ordered phase forms.

Upon introducing frictional forces between particles after reaching the ordered state, we observe that the ordering is rapidly disrupted. Friction changes the rheological properties of the suspension by drastically increasing viscosity and causing the suspension to shear thicken instead of shear thin. Our analysis of the particle contributions to the stress shows that these phenomena are accompanied by a major redistribution of the stresses within the microstructure. Instead of being localized in less dense regions as in the frictionless case, the particles with large stress contributions form a percolating network that spans low and high density regions. It is possible that the additional constraints on particle motion due to friction may serve to stabilize and grow the stress-carrying structures that were formed in the frictionless case [51], thus, disrupting the order and leading to the observed network. These results show how the nature of particle contact, for otherwise identical suspensions, can control the transition from steady state order to disorder, or thinning to thickening. Additionally, they shed new light on how the link between particle interactions and microstructure controls both the stress distribution and the macroscopic flow properties.

## ACKNOWLEDGMENTS

The authors acknowledge the NIST PREP Gaithersburg Program (No. 70NANB18H151), the National Science Foundation (No. NSF DMR-2026842), and Georgetown University for support.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.