Rheological investigation of gels formed by competing interactions: A numerical study

TL;DR

This numerical study explores how competing short-range depletion and long-range electrostatic interactions influence the rheological behavior and yielding of colloidal gels, revealing a systematic dependence on inter-particle potentials and flow conditions.

Contribution

It introduces a detailed numerical analysis of equilibrium and shear-induced states in gels with competing interactions, highlighting the impact of inter-particle potentials on rheological transitions.

Findings

Single yielding event at fixed strain in stress-strain curves

Stress overshoot follows a power-law with Péclet number, with a larger exponent than depletion gels

Steady states with the same Péclet number exhibit identical structural and rheological properties

Abstract

A transition from solid-like to liquid-like behavior occurs when colloidal gels are subjected to a prolonged exposure to a steady shear. This phenomenon, which is characterized by a yielding point, is found to be strongly dependent on the packing fraction. However, it is not yet known how the effective inter-particle potential affects this transition. To this aim, we present a numerical investigation of the rheology of equilibrium gels in which a short-range depletion is complemented by a long-range electrostatic interaction. We observe a single yielding event in the stress-strain curve, occurring at a fixed strain. The stress overshoot is found to follow a power-law dependence on the P\'eclet number, with an exponent larger than that found in depletion gels, suggesting that its value may depend systematically on the underlying colloid-colloid interactions. We also establish a mapping between equilibrium states and steady states under shear, which allows us to identify the structural modifications induced by the presence of the shear. Remarkably, we find that steady states corresponding to the same P\'eclet number, obtained by different combinations of shear rate and solvent viscosity, show identical structural and rheological properties. Our results highlight the importance of understanding the coupling between colloidal interactions, solvent effects, and flow to be able to describe the microscopic organization of colloidal particles under shear.

Figures

Equations

$$V(r)=4ϵ\left[{\left(\frac{\sigma }{r}\right)}^{2\alpha }-{\left(\frac{\sigma }{r}\right)}^{\alpha }\right]+A\frac{e^{-\kappa r}}{r/\sigma }.$$

$$Pe=\frac{F_{visc}}{F_{bond}}$$

$$F_{visc}=\xi v_s$$

$$F_{bond}=\frac{\Delta V\left(r\right)}{\Delta r}=\frac{V\left(r_{bond}\right)-V\left(r_{min}\right)}{r_{bond}-r_{min}}.$$

$${\Sigma }_{xz}=\frac{1}{V}⟨\sum _i[m_i{v}_{i,x}^{\prime }{v}_{i,z}+\sum _{j>i}{r}_{ij,x}{F}_{ij,z}]⟩$$

$$S_{eq}(q)=\frac{1}{N}⟨\sum _{ij}{e}^{-i\mathbf{q}\left({\mathbf{r}}_i-{\mathbf{r}}_j\right)}⟩,$$

$$S_{shear}\left(q_x,q_y,0\right)=\frac{1}{N}⟨\sum _{ij}{e}^{-i[{\mathbf{q}}_{\mathbf{x}}\left(x_i-x_j\right)+{\mathbf{q}}_{\mathbf{y}}\left(y_i-y_j\right)]}⟩.$$

$$g\left(\mathbf{r}\right)=g_s\left(r\right)+\sum _{l=1}^{\infty }\sum _{m=-l}^l{g}_l^m\left(r\right)Y_l^m\left(\theta \varphi \right)$$

$$Im{g}_2^{-2}\left(r\right)=\sqrt{\frac{15}{8\pi }}\frac{L^3}{N^2}\times$$

$$⟨\sum _i^N\sum _{j\ne i}^N\delta \left(\left|{\mathbf{r}}_i-{\mathbf{r}}_j\right|-r\right)\frac{\left(x_i-x_j\right)\left(z_i-z_j\right)}{r^4}⟩.$$

$${\overline{q}}_{lm}\left(i\right)=\frac{1}{{\tilde{N}}_b\left(i\right)}\sum _{j=0}^{{\tilde{N}}_b\left(i\right)}{q}_{lm}\left(j\right)$$

$${\overline{q}}_l\left(i\right)=\sqrt{\frac{4\pi }{2l+1}\sum _{m=-l}^l{\left|{\bar{q}}_{lm}\left(i\right)\right|}^2}$$

$$\begin{split}\overline{w}_{l}\left(i\right)=\frac{\underset{m_{1}+m_{2}+m_{3}=0}{\sum}\left(\begin{array}[]{ccc}l&l&l\\ m_{1}&m_{2}&m_{3}\end{array}\right)\overline{q}_{lm_{1}}\left(i\right)\overline{q}_{lm_{2}}\left(i\right)\overline{q}_{lm_{3}}\left(i\right)}{\left(\sum_{m=-l}^{l}\left|\overline{q}_{lm}\left(i\right)\right|^{2}\right)^{3/2}}\end{split}$$

Full text

colloidal gels, rheology

††preprint: AIP/123-QED

I Introduction

Colloidal particles often form disordered arrested states, such as glasses and gels. Depending on the colloid-colloid interactions, different kinds of glasses can be found, including attractive, repulsive or Wigner glasses Pham et al. (2002); Zaccarelli and Poon (2009); Bonn et al. (1999); Ruzicka and Zaccarelli (2011). Similarly colloidal gels can be formed by different routes Gleim, Kob, and Binder (1998). An important distinction can be made to distinguish out-of-equilibrium gels formed via arrested spinodal decomposition from cases where gelation is obtained in equilibrium Sciortino and Zaccarelli (2017). Non-equilibrium gels are found when colloids interact via hard-sphere-like excluded volume complemented by an isotropic short-range attraction, which is typical of depletion effects induced by non-adsorbing polymer chains. At high enough depletion strength, a colloid rich-colloid poor phase separation takes place, in which the dense phase undergoes dynamical arrest into a gel state which interrupts the spinodal decomposition process Lu et al. (2008). On the other hand, gels can be obtained in equilibrium from a homogeneous fluid state, when attraction is not isotropic (e.g. patchy or limited valenceSciortino and Zaccarelli (2017)) or when this is counter-acted by an additional long-range repulsion due to charge effects Groenewold and Kegel (2001); Stradner et al. (2004); Sciortino et al. (2004). In this case, the competition between short-range attraction and long-range repulsion is able to avoid phase separation Campbell et al. (2005); Sciortino, Tartaglia, and Zaccarelli (2005). Since colloidal gels are widely used for a variety of applications, including biomedical purposes Khetan and Burdick (2011); Guvendiren, Lu, and Burdick (2012) food processing Mezzenga et al. (2005), optical sensing, thermoelectrics or catalysis Gaponik, Herrmann, and Eychmuller (2011), it is fundamental to control the gelation process and to be able to discriminate among the wide zoology of colloidal gels.

To study the behavior of colloidal arrested states, one important experimental tool is rheology. In general, the application of a shear flow causes the occurrence of a solid-like to liquid-like transition that is preceded by a stress overshoot $\Sigma_{yield}$ at the yielding point. This indicates the maximum stress that the system can accumulate Liddel and Boger (1996); Mason, Bibette, and Weitz (1996). After yielding, the system is able to approach a steady state with liquid-like behavior. The way in which this steady state is reached of course depends on the studied system and on the shear conditions. However, different types of arrested states generally respond in a different manner. For the widely studied hard-sphere glasses, for which excluded volume interactions are responsible for the kinetic arrest, a single yield mechanism is observed Derec et al. (2003); Petekidis, Vlassopoulos, and Pusey (2004). This is normally attributed to cage breaking only or to particle exchange with their nearest neighbours. On the other hand, for attractive glasses, which can be induced by depletion interactions at relatively high packing fractions $\phi$, two different yielding points have been reported Pham et al. (2004, 2008). A first one ($\Sigma_{yield1}$) is associated to bond breaking at local level with the system retaining a solid-like character, while a second one ($\Sigma_{yield2}$) is related to a structural rearrangement (cage-breaking) after which the system is able to flow. Between the two yielding events, there is a local bond reorganization which, however, does not significantly alter the system. A double yielding mechanism has also been observed for depletion-induced gel-like samples at $\phi\sim 0.40$, Koumakis and Petekidis (2011) an effect that could possibly be associated to the relatively large packing fraction. Indeed, the caging effect disappears for $\phi\lesssim 0.2$,Koumakis and Petekidis (2011) leaving only one yielding point associated to bond-breaking for low density gels. In addition, it is interesting to note that for very large attraction strengths $\left(\sim 100k_{B}T\right)$, experiments on very dilute depletion gels also reported the occurrence of two yielding mechanisms Chan and Mohraz (2012): in this case the bond breaking yielding point observed for weaker attractions is preceded by another yielding point associated to the onset of bond rotation. Furthermore, it is important to consider that the presence of an imposed shear flow naturally induces anisotropy or heterogeneity in the structure of the system Varga and Swan (2018), which affects the rheological properties Vermant and Solomon (2005). This aspect has been tackled both in experiments Rajaram and Mohraz (2010); Min Kim et al. (2014); Eberle et al. (2014) and simulations Colombo and Del Gado (2014); Park, Ahn, and Wagner (2017); Moghimi, Jacob, and Petekidis (2017); Boromand, Jamali, and Maia (2017); Jamali, McKinley, and Armstrong (2017); Johnson, Landrum, and Zia (2018).

To our knowledge, no rheological studies —either numerical or experimental— have been performed to date on gels obtained via the competition of short-range attraction and long-range repulsion. The aim of the present work is to fill this gap by investigating this type of gels via numerical simulations performing start-up shear experiments for a wide variety of steady shear conditions. Mainly, we focus on Langevin Dynamics (LD) simulations, which is appropriate to treat colloids in an implicit solvent. With this approach, we provide evidence that, while some rheological features of depletion-induced gels are also found in our case, others reveal that the inter-particle potential plays an important role. Indeed the long-range repulsion makes our gels more resistant to shear flow, with a modified dependence of the stress overshoot with respect to shear (quantified by Péclet number). Similarly, we study the implication of the presence of shear flow at the microscopic level, focusing on the evolution of network bonds and the anisotropy. Since the competition between Brownian motion and shear flow modifies the gel response, we show that a different balance of the two effects allows to produce different anisotropies in the sheared system. Interestingly, we find that the resulting anisotropic patterns take a characteristic form when the gel undergoes crystallization under shear. As a next step we focus on the role of solvent on mechanical response and microscopic gel restructuring. To this end, we also implement Molecular Dynamics (MD) simulations, where the solvent is neglected. This has the aim to clarify the influence of the microscopic dynamics on the sheared systems, a practice that has been carried out for example in glassy systems, where it was found that the long-time dynamics does not depend on the presence of a solventGleim, Kob, and Binder (1998). Connecting equilibrium and steady states via a mapping at equal potential energy allows us to highlight the effect of the shear in counter-acting the long-range repulsion, while not significantly altering the local (average) structure as compared to gels without shear. We find that steady states obtained at equal Péclet number are identical to each other, independently of the effective solvent viscosity used in the simulations. However, the comparison with MD simulations reveals that a different microscopic dynamics inevitably alters the rheological response of the gel, but does not affect its average thermodynamic properties. To validate the robustness of our results, we also repeat some simulations for larger system sizes finding that no size effects are observed for the steady state properties that we have considered. Some differences however arise in the transient immediately after the start-up of the shear.

The paper is organized as follows. In Section II we describe the numerical simulations for the different methods that we use and we also define all the observables that are calculated in this study. In Sec. III, we present our results discussing the microscopic organization of the gel in the presence of shear (Sec. III.1), the anisotropy generated by the competition between Brownian dynamics and shear flow (Sec. III.2), and the shear-induced crystalline structures formed by the gel (Sec. III.3). We also show the results of the mapping between equilibrium states and stable states under shear (Secs. III.4-III.5) and we conclude the section by addressing the problem of system size effects (Sec. III.6). Finally, in Sec. IV, we discuss our results and present the conclusions of our study.

II Methods: Simulations and theory

II.1 Simulation details

We perform simulations of monodisperse colloids of diameter $\sigma$ and mass $m$ in presence of steady shear. Most of the simulations are carried out with $N=2000$ colloidal particles, but in order to assess size effects onto the shear results, we also repeat some simulations for $N=5000$ and $N=10000$ colloids. Particles interact via a potential $V(r)$ which is the sum of a short-range attraction and a long-range repulsion, Sciortino et al. (2004); Sciortino, Tartaglia, and Zaccarelli (2005) as

[TABLE]

Here the short-range attraction (mimicking depletion interactions) is modeled as a generalised Lennard-Jones potential Vliegenthart, Lodge, and Lekkerkerker (1999) with the potential depth $\epsilon$ and the particle diameter $\sigma$ being, respectively, the units of energy and length, while the long-range repulsion (representing a screened electrostatic contribution) is described by a Yukawa potential where $\kappa$ is the inverse of the Debye screening length and $A$ is the repulsion amplitude. Time is measured in units of $\sqrt{(m\sigma^{2})/\epsilon}$. Following Ref. Sciortino, Tartaglia, and Zaccarelli (2005), we fix $\alpha=18$, $A=4\epsilon$ and $\kappa=2\sigma^{-1}$. The resulting interaction potential is illustrated in Fig. 1. In our simulations we fix $k_{B}=1$ and use a cutoff of the interactions at $r_{c}=4\sigma$. Particles interacting with this potential are able to form an equilibrium gel at low/intermediate packing fractions $\phi$ and sufficiently low temperatures Sciortino, Tartaglia, and Zaccarelli (2005); Zaccarelli (2007). In this work we consider a packing fraction $\phi=\frac{\pi}{6}\sigma^{3}\frac{N}{V}$=0.16 and $T=0.1$, where $V$ is the volume of the cubic simulation box. To reduce statistical noise, data are always averaged over three independent realizations.

In our study we first prepare a gel state in equilibrium and then we perform a start-up shear test by applying a steady shear flow onto the gel imposing the so-called Lees-Edwards boundary conditions Lees and Edwards (1972). We consider the gradient velocity to be in the $\hat{z}$ direction while the shear velocity is in the $\hat{x}$ direction, so that the shear rate is defined as $\dot{\gamma}\equiv v_{x}/z$ and is measured in inverse time units. We employ a Langevin thermostat which acts on the so-called peculiar velocityShang, Kroger, and Leimkuhler (2017); Ruiz-Franco, Rovigatti, and Zaccarelli (2018), which is defined for particle $i$ as $v^{\prime}_{i,x}=v_{i,x}-u_{x}(z)$ where $v_{i,x}$ is the $x$-component of the particle thermal velocity and $u_{x}(z)=\dot{\gamma}z$ is the stream velocity. Simulations are performed in the canonical ensemble $NVT$ and the properties of the solvent are controlled by its viscosity $\eta$, which in turn determines the friction coefficient $\xi=3\pi\eta\sigma$.

To monitor how the microscopic dynamics affects the shear response of the gel, we also perform simulations in the absence of an implicit solvent. To this aim, we use a so-called Gaussian thermostat, for which the temperature is controlled by imposing a constant kinetic energy (iso-kinetic ensemble) Morriss and Dettmann (1998). The equations of motion are then solved by means of the SLLOD integrator Zhang et al. (1999).

II.2 Definition of a generalized Péclet number

The strength of the applied shear is usually quantified by the Péclet number $Pe$, which controls the balance between Brownian motion and shear effects: the behaviour of the system is essentially governed by Brownian motion for $Pe<$ 1, while it is dominated by the shear flow for $Pe>$1. However, in its standard definitionCloitre (2010), the Péclet number is taken to be equal to $\dot{\gamma}\tau_{B}$, where $\tau_{B}=\frac{\sigma^{2}}{4D}$ is the the Brownian time and $D=\frac{k_{B}T}{\xi}$ is the diffusion coefficient of the particle at infinite dilution. This definition is not appropriate to quantify the strength of the shear on the gel state because the system is far from dilute conditions and cannot be simply generalized because the diffusion coefficient in the gel tends to zero. To provide a meaningful definition of $Pe$, we adopt the modified definition of Ref.Koumakis and Petekidis (2011)

[TABLE]

which quantifies the resistance of the bonds between two colloids to the shear. Here $F_{visc}$ is the drag force that is able to displace two particles up to a distance larger than the attractive range of the potential:

[TABLE]

where $v_{s}=\dot{\gamma}r_{min}$ and $r_{min}\approx 1.05\sigma$ is the global minimum of potential, i.e. the equilibrium distance of two particles (see Fig. 1). Instead, $F_{bond}$ is the bonding force which is responsible for maintaining the bond between the two colloids. By increasing the interparticle distance $F_{bond}$ will increase until the two particles become more far than a maximum distance $r_{bond}$, which is the maximum of the potential, after which the effective force becomes repulsive. In the present case, $r_{bond}=1.28\sigma$ as shown in Fig. 1. We thus define the bonding force as the variation of energy in the attractive range of the effective potential $\Delta r=r_{bond}-r_{min}$, as

[TABLE]

Using the definition of eq.(2) we find that when the two forces are balanced, i.e. $Pe=1$, the original Péclet number calculated under dilute conditions would be much higher, i.e. $\approx 320$.

We notice that the Péclet number can be defined only for simulations in the presence of a solvent such as LD ones. For MD simulations we thus quantify the strength of the shear by varying the shear rate only.

II.3 Calculated observables

During application of the steady shear, we calculate the internal stress tensor $\Sigma_{xz}$ using the Irving-Kirkwood expression Irving and Kirkwood (1950):

[TABLE]

where $r_{ij}$ and $F_{ij}$ are, respectively, the distance and the force between particles $i,j$ and the brackets $\langle\ldots\rangle$ represent the ensemble average.

To provide a microscopic understanding of the stress tensor behavior, we monitor the time evolution of the bond organization between the particles, calculating: (i) the fraction of bonds $f_{b}$ defined as the number of bonds in a given configuration under shear divided by the number of bonds that were present in the system prior to switching on the shear; (ii) the fraction of unbroken bonds $f_{u}$ that were also present at zero-shear; (iii) the bond potential energy $E_{b}$ Landrum, Russel, and Zia (2016), defined as the absolute value of the average potential energy between all pairs of bonded particles, normalized to its initial value, again prior to switching on the shear flow.

Next, we examine the changes in the structure by calculating the static structure factor: in equilibrium this is defined as

[TABLE]

while when shear is applied we evaluate it in the velocity-vorticity plane at $q_{z}=0$, as

[TABLE]

For simplicity, we refer to both types of structure factors as $S(q)$ in the main text, since they depend only on the modulus of the wavevector, being it calculated in 3D or in 2D along the velocity-vorticity plane.

A complementary picture of the structure, which allows to identify the anisotropy induced by shear, is provided by using a suitable expansion of the pair correlation function $g\left(\mathbf{r}\right)$. In particular, we consider the expansion

[TABLE]

into spherical harmonics $Y_{lm}\left(\theta,\phi\right)$ Hess and Hanley (1983); Hanley, Rainwater, and Hess (1987), where $g_{s}\left(r\right)$ is the usual (isotropic) radial distribution function Hansen and MacDonald (2006) and the expansion coefficients are $g_{lm}\left(r\right)=\int g\left(\mathbf{r}\right)Y_{lm}^{*}d\Omega$ with $\Omega$ the solid angle and $d\Omega=sin\theta d\theta d\phi$. Our colloidal gel is made up of identical particles, implying that $g\left(\mathbf{r}\right)=g\left(-\mathbf{r}\right)$. This ensures that only even values of $l\geq 2$ have to be accounted for Hess and Hanley (1983); Hanley, Rainwater, and Hess (1987). However, coefficients with $l>2$ are of small amplitude, so that, in general, it is sufficient to consider only the term with $l=2$ Zausch (2009). In addition, due to the geometry of the planar Couette flow, only coefficients with $m=\pm 2$ are in the shear flow plane and the only non-zero contribution comes from their imaginary part. We thus focus on the imaginary part with $m=-2$, i.e. $Im\>g_{2}^{-2}\left(r\right)$, which is calculated as Zausch and Horbach (2009):

[TABLE]

This function is characterized by the presence of two peaks: the first one is a minimum signalling the accumulation of the particles along the compression axis, while the second one is a maximum which corresponds to the depletion of the particles along the extension axis Park, Ahn, and Wagner (2017); Khabaz et al. (2017).

To identify solid-like particles we use the local bond-order analysis introduced by Steinhardt et al. Steinhardt, Nelson, and Ronchetti (1983), where the complex vector $q_{lm}\left(i\right)$ of particle $i$ is defined as $q_{lm}\left(i\right)=\frac{1}{N_{b}\left(i\right)}\sum_{j=1}^{N_{b_{i}}}Y_{lm}\left(\hat{r}_{ij}\right)$, where $N_{b}\left(i\right)$ is the set of bonded neighbours of particle $i$ and $\hat{r}_{ij}$ is the unit vector specifying the orientation of the bond between particles $i$ and $j$. Using the complex vectors $q_{6m}$, we are able to assign a solid connection between particles $i$ and $j$ if $d_{6}\left(i,j\right)=\sum_{m=-6}^{6}q_{6,m}\left(i\right)\cdot q{}_{6,m}^{*}\left(j\right)\geq 0.7$. A particle is defined to be solid-like if it has 6 or more solid connections with its neighbours Pusey et al. (2009). The percentage of solid particles is thus defined as $X\left(t\right)=\frac{N_{X}}{N}$ with $N_{X}$ the number of solid-like particles. Following Ref. [60], we also calculate the rotationally invariant bond order parameters $\overline{q}_{l}\left(i\right)$ and $\overline{w}_{l}\left(i\right)$. To define these two parameters, it is necessary to compute the averaged local bond order parameters:

[TABLE]

where $\widetilde{N}_{b}\left(i\right)$ is the number of neighbours including the particle $i$ itself. In this way, the first invariant bond order parameter is defined as

[TABLE]

while the second one is defined as

[TABLE]

where the term in parentheses is the Wigner 3-j symbol. The integers $m_{1}$, $m_{2}$ and $m_{3}$ runs from $-l$ to $l$ but only the combination that meets the requirement $m_{1}+m_{2}+m_{3}=0$. Using $l=4$ and $l=6$ it is possible to establish a separation between BCC, FCC and HCP structures Russo and Tanaka (2012).

III Results

III.1 Microscopic organization of the gels under shear

We start by showing the behavior of the stress tensor under shear. To facilitate a comparison with experimental results on colloidal gels, we focus on the case of LD simulations, which take into account the presence of the solvent in an effective way, although neglecting hydrodynamic interactions. As discussed above, for several attractive arrested states Pham et al. (2008); Koumakis and Petekidis (2011); Chan and Mohraz (2012) two yielding mechanisms have been observed. It is now interesting to see what happens in the present case of equilibrium gels obtained by competing interactions. The behaviour of the normalized stress tensor $\Sigma_{xz}\sigma^{3}/8k_{B}T$ is reported as a function of strain $\gamma=\dot{\gamma}t$ for several values of $Pe$. We find that only one yielding mechanism takes place in our gel, confirming the results reported for depletion-induced gels at comparable packing fraction and attraction strength Koumakis and Petekidis (2011). Similarly to previous findings, the position of maximum of the stress tensor, defining the yielding point $\Sigma_{yield}$, is insensitive to shear rate Whittle and Dickinson (1997); Koumakis and Petekidis (2011). We also find a power law dependence of the yield stress on the shear rate, i.e., $\Sigma_{yield}\propto$ Peδ. However, for depletion-induced gels at intermediate $\phi$ the power law exponent has been found to be $\delta\sim 0.5$ both in simulations Whittle and Dickinson (1997) and experiments Koumakis and Petekidis (2011). On the other hand, numerical simulations for the Derjaquin-Landau-Verwey-Overbeek (DLVO) potential Park and Ahn (2013) have found $\delta\sim 0.56$. We find that for our model the yield stress follows a power law with an exponent $\delta\sim 0.62\pm 0.01$ (see inset in Fig. 2(a)). This result suggests that there is a systematic change of $\delta$ with the employed interaction potential between the colloids. Interestingly for dense colloidal glasses, an almost constant value of the stress overshoot with $Pe$ was foundLaurati et al. (2012). To provide an interpretation of these findings, we refer to early theoretical studies on polymer Groot and Agterof (1995) and colloidal gels Whittle and Dickinson (1997), which reported a link between $\delta$ and the fractal dimension $d_{f}$ of the system, i.e. $\delta=\frac{3-d_{f}}{2}$. Using this relationship, the fractal dimension would decrease from $\sim 2.0$ to $\sim 1.75$ as we add the long-range electrostatic repulsion to the short-range depletion one. This is in agreement with direct estimates of the fractal dimensions from assessment of the clusters in equilibrium Sciortino, Tartaglia, and Zaccarelli (2005), where a very low value of $d_{f}\sim 1.25$ was found. It is plausible that this value is slightly increased by the presence of the shear. Thus, from this type of result we can get an indirect estimate of the gel structure from rheological measurements. It would be interesting to test this relationship to systems with different fractal dimensions. We may also speculate that an increase of $\delta$ can be interpreted in terms of a stronger resistance to the flow of the gel under shear at structural level. This interpretation is confirmed by comparing our results with those presented in Ref. Johnson, Landrum, and Zia (2018), where numerical simulations of a gel were performed using a Morse potential with a minimum value of the energy $\sim 6k_{B}T$, similarly to our case, reporting an exponent $\delta\leq 0.5$. In our model, the presence of the long-range repulsion could possibly counteract the effect of the shear in breaking bonds, inducing an increase of $F_{bond}$ and thus, making our gel more resistant to shear flow, which is thereby reflected by an increase of $\delta$.

To deepen our microscopic understanding of the yielding point, we monitor the fraction of bonds $f_{b}$ and the bond energy $E_{b}$ in Fig. 2 (b) and (c) for the same $Pe$ values. We find that at yielding the total number of bonds does not change significantly, while the bond energy already starts to decrease. Interestingly, if we separately monitor only the fraction of unbroken bonds $f_{u}$ with respect to the zero-shear case, we find no broken bonds until yielding. Thus, the gel network remains essentially unaffected up to $\Sigma_{yield}$ in agreement with Refs. Koumakis and Petekidis (2011); Laurati, Egelhaaf, and Petekidis (2011); Johnson, Landrum, and Zia (2018). The decrease of the energy however might indicate that bonds become more and more stretched under the shear deformation, prior to eventual breakage. After the yielding point, when the bond energy has already dropped a significant amount, $f_{u}$ rapidly decreases until reaching a steady state where basically none of the initial bonds is left intact. However, new bonds appear and indeed, $f_{b}$ reaches a new plateau in the steady state. In all cases, the system is found to form a new percolating network, whose structure is very different from the initial gel state as we will show in section III.4.

III.2 Anisotropy

From the microscopic point of view, the anisotropy induced by the shear flow has been studied to understand how the balance between Brownian motion and shear flow can affect the microscopic structure for different systems. We have thus studied the evolution of the anisotropy of our gels by calculating Im$g_{2}^{-2}(r)$, for different Pe (changing both $\dot{\gamma}$ and $\xi$) in the LD simulations. We incidentally notice that for hard sphere systems studied within Brownian dynamics an inversion of the peak amplitudes in Im$g_{2}^{-2}(r)$[60] was observed upon increasing $Pe$, while for low shear contribution, the compressional peak is more pronounced than the extensional one, this situation was found to be reversed at high shear. This was interpreted as a consequence of the net contribution of the shear with respect to the Brownian dynamics. At high $Pe$, the convective motion of hard spheres increases and hence a larger anisotropy is found in the extensional axis. For the current system, anisotropy results are reported in Figs. 3 (a,b) for $\xi=10^{2}$ and in Figs. 3 (c,d) for $\xi=10$. For each solvent condition, we have monitored the evolution of anisotropy in the system and associate it to the analysis of the stress and bonds reported above, including the yielding manifestation up to the steady state.

For $\xi=10^{2}$ we find that the anisotropy starts to grow already well before the yielding point (occurring for $\gamma\sim 10^{-1}$). A characteristic two-peaked shape is observed, which is made of a negative peak followed, at larger distances, by a positive one. These two peaks indicates the increase of anisotropy along the compressional and extensional axis, respectively Park, Ahn, and Wagner (2017). Interestingly, the distance at which Im$g_{2}^{-2}\left(r\right)$ passes through zero corresponds to the minimum $r_{min}$ of the total interaction potential. In addition, we observe that the positive peak does not exceed the maximum bond distance $r_{bond}$, showing that this modification of the gel structure affects the bonds between the particles and not larger distances. At yielding, where still the network is intact but the maximum accumulation of stress occurs, the induced anisotropy is also maximum. It is found that anisotropy at this point is able to propagate at distances larger than the bond ones, affecting the whole structure, as shown by additional peaks arising for $r>r_{bond}$. When the bonds finally start to break, after yielding, the anisotropy also decreases and the peaks at larger distances disappear. Interestingly, after this happens, Im$g_{2}^{-2}\left(r\right)$ shows a long-range tail which is most evident in Fig. 3 (a,b). This tail could indicate that the dissipation of the anisotropy does not happen instantaneously but occurs within a finite time, in correspondence with the smooth decay of $f_{b}$ and $f_{u}$ after the yielding point. Indeed, at larger strains, the tail disappears. However, even in the steady state ($\gamma\gtrsim 10$), a significant amount of anisotropy still remains at local level. The reported behavior is enhanced for higher values of $Pe$, which in general induce a larger amount of anisotropy in the system Johnson, Landrum, and Zia (2018). The behavior described for hard-spheres systems above is not observed Foss and Brady (2000); Mewis and Wagner (2012), probably due to the fact that our short-range repulsion, although very steep, is not hard-sphere-like, thus allowing for a moderate compression of the particles themselves. This can be seen in the the compressional peak moving at slightly smaller distances with increasing $Pe$. Indeed, a similar feature was observed in the anisotropy of deformable particles Khabaz et al. (2017).

On the other hand, we find remarkable differences in the behavior of Im$g_{2}^{-2}(r)$ for $Pe\lesssim 3\cdot 10^{-2}$, as shown in Fig. 3(c,d). In this regime, the Brownian motion is much stronger than the shear flow. We monitor this behavior for the lower studied friction coefficient ($\xi=10$) and find that, after the yielding point, the extensional peak becomes much more pronounced while the compressional one tends to disappear. To clarify the effects due to the competition between Brownian dynamics and the shear flow, we compare the anisotropy distribution for two systems having the same shear rate $\dot{\gamma}=0.05$ but different friction coefficients: respectively $\xi=10^{2}$ in Fig. 3(b) and $\xi=10$ in Fig. 3(d). It is clear that, tuning the solvent properties at the same shear rate, we can manipulate the anisotropy induced in the system by the shear flow. While at high $Pe$ both extension and compression of the bonds take place, for low $Pe$ only the extensional axis grows. This feature, coupled to the lack of a clear yielding point (see Fig. 2), indicates that the system is undergoing a restructuring process. Under these conditions, the small perturbation acted by the shear flow onto the system is enough to allow the particles to reorganize towards a more ordered configuration. Indeed, a fluid-to-crystal transition is observed at higher strains, as shown in the next section. On the other hand, at large frictions (and hence higher $Pe$) the Brownian motion dominates the effect of the shear and crystallization does not occur, at least on the simulated time scales. Interestingly, a recent study of a jammed suspension under steady shear Khabaz et al. (2017) reported a higher accumulation of particles along the compression axes prior to crystallization, an opposite result to the present case. However, in that system, particles interact with a soft Hertzian repulsion and thus, the high deformability should be responsible for the observed behavior. Despite the differences, this study shows that such a large asymmetry of the two peaks of Im$g_{2}^{-2}\left(r\right)$ is a distinctive signature of an incipient crystallization. It would be interesting to confirm this feature in other works, either numerical or experimental ones.

III.3 Shear-induced crystallization

As anticipated in the previous section, we find that our gels undergo crystallization when we perform LD simulations at low enough friction coefficient ($\xi=10$). We have identified the important role of friction which, if too high, counteracts the effect of shear, and acts against ordering. To further strengthen this point, we also perform MD simulations, as described in the Methods, where the solvent is absent, which would thus mimic an atomic, rather than a colloidal, system.

Crystallization is found when a sudden drop in the potential energy occurs for a given trajectory. To quantify the transition, we monitor the fraction of solid particles $X$ and we calculate the bond local order parameters in order to discriminate between different crystal structures Lechner and Dellago (2008); Russo and Tanaka (2012) (see Methods for details). We report the strain evolution of $X$ for a few selected shear conditions in Fig. 4(a) for both LD and MD simulations. In order to quantify shear also in the case of MD, we refer to the value of the shear rate $\dot{\gamma}$, because in this case a $Pe$ cannot be defined. We find that crystallization occurs only in a narrow region of shear rates, which are not too large to be able to destroy the order and not too small in order to induce a significant rearrangement. Thus for the small $Pe$ difference considered in Fig. 4(a), we find that the system sheared with $Pe\approx 2\cdot 10^{-2}$ is able to crystallize, while the system with slightly larger $Pe\approx 3\cdot 10^{-2}$ is not able to crystallize within the simulated time window. However, an important point is that the final crystal state is the same for both LD and MD simulations. This is characterized by a predominant face-centered-cubic (fcc) structure, as shown in Fig. 4 (b1, b2) where the bond orientational parameter $\bar{w}_{6}$ and $\bar{w}_{4}$ are reported. A negative value of $\bar{w}_{6}$ can be used to discriminate a fcc arrangement from a body-centered-cubic (bcc) one. However, the fcc structure is quite similar to the hexagonal-close-packed (hcp) one in terms of $\bar{w}_{6}$. In order to discriminate between the latter two crystal structures, one needs to consider $\bar{w}_{4}$, which is predominantly negative for fcc. For even longer simulation times, we find that the system acquires a well-defined fcc crystal at $\gamma\geq 10^{4}$. We could only reach this long-time regime using MD simulations within our simulated time window.

III.4 Mapping between equilibrium and steady states under shear: structure

Up to this point we have investigated the microscopic structure of the system under shear, quantifying the anisotropy and detecting the onset of crystallization under specific shear conditions. After the yielding point, the system approaches a steady state whose microstructure will be different depending on $Pe$ or on the underlying microscopic dynamics. In particular, we have compared three cases — the absence of a solvent (MD) and two implicit solvents implemented through LD at respectively low ($\xi=10$) and high viscosity ($\xi=10^{2}$) — finding that the solvent can affect the kinetics of the deformation of the initial structure induced by the shear. It would thus be useful to have a way to compare these three cases when shear conditions are equivalent with respect to the underlying Brownian motion (or in the absence of it). To this aim, we can use the equilibrium states of the system as reference states and quantify the effect of the employed shear in each case with respect to them. In particular, we build a correspondence to equilibrium states using the potential energy as mapping observable. Thus, for any applied shear under MD and LD conditions, we consider the potential energy of the steady state and map it to the equilibrium state with the same potential energy. In this way, we establish a shear rate-temperature connection linking steady states obtained under different types of shear and equilibrium states. The obtained mapping is represented in Fig. 5(a), where we use $\dot{\gamma}$ as mapping variable in order to include also the MD simulations.

For each $T$ in equilibrium, a set of corresponding steady states arising from different simulation methods is defined by different values of $\dot{\gamma}$. The weaker is the effect of the solvent, the higher is the shear rate corresponding to the same equilibrium state. To see whether the mapping is meaningful we now consider two sets of states labeled in the Fig. 5(a) as Set I and Set II, corresponding to states in equilibrium with $T=0.24$ and $T=0.14$ respectively, and we compare the structure of these sets with their corresponding equilibrium states. In Fig. 6 we report the static structure factors, calculated as defined in Methods in equilibrium and under shear, for these two sets of corresponding states. Remarkably, we find that the $S(q)$ for all steady states obtained under different shear conditions and dynamics are superimposed onto each other within the statistical uncertainty of the numerical data. This is a confirmation of the efficiency of the mapping in connecting steady states at different shear rates among themselves, and with respect to equilibrium. Thus if one wants to compare different shear conditions, one needs to consider a different $\dot{\gamma}$ in order to arrive at a similar steady state structure.

In addition, Fig. 6 clearly shows the effect of shear on the microscopic structure of the system. While for equilibrium states a cluster peak is observed, extensively discussed in the literature as a generic feature of competing interactions, such peak disappears in the presence of shear in favour of a growing intensity of $S(q)$ for $q\rightarrow 0$. This allows us to deduce that, in the steady state, the shear flow acts essentially at large length scales, by screening the contribution of the long-range repulsion and enhancing the attractive interactions between the colloids as compared to equilibrium case. In this way the shear drives the system closer to phase separation and the presence of larger density fluctuations with respect to equilibrium also helps crystallization of the system. This is visible in the snapshots reported in Fig. 6(b), where the inhomogeneity of the structure is evident and confirms the findings of enhanced anisotropy after yielding reported in Fig. 3(d) prior to crystallization. It is important to notice that for steady states corresponding to lower $T$ (not shown), we find no growth of $S(q)$ at small wavevectors and consequently no crystallization. Finally, focusing on wavevectors larger than the nearest-neighbour peak, we see that the shear has a much weaker effect, leaving almost unaltered the local structure of the system. Hence, in our gels with competing interactions, shear essentially acts against the long-range repulsion and is able to strengthen the attractive interactions. It would be interesting to repeat this analysis for other types of gels in order to highlight the different effect of shear in those cases.

III.5 Comparing different steady states under shear through the established mapping: invariance on Péclet number, but dependence on microscopic dynamics

The established mapping not only allows us to identify the effect of the shear on the structure of the system, but also allows us to compare the rheological response of corresponding steady states. If the mapping is meaningful, this response should also be identical. We also notice that, from the definition of the Péclet number (Eq. 2), there are different possible combinations of $\xi$ and $\dot{\gamma}$ that allow the same $Pe$ value to be obtained. In Fig. 5(b) the potential energy per particle $E/N$ for LD simulations at the two different values of $\xi$ studied here is shown as a function of $Pe$ instead of shear rate. Thus it is evident that our mapping implies that state points with identical $Pe$ also have the same potential energy, confirming that it is only the balance between shear rate and Brownian motion that should determine the rheological response of the system.

We now compare steady states obtained for LD simulations with $Pe\approx 3\cdot 10^{-2}$ and varying $(\dot{\gamma},\xi)$ combinations, as highlighted in Fig. 5(a) with label Set II, where the steady states have the same potential energy. In Fig. 7(a) we show that steady states obtained under shear for the same $Pe$ not only possess an identical structure but also display identical anisotropy distribution at all strain values. We confirm that the same results also hold also for different sets of corresponding states showing asymmetric features of Im$g_{2}^{-2}\left(r\right)$ as in Fig. 3(d) (not shown). While this may seem like an obvious result, it offers the possibility to vary independently the two parameters in LD simulations (i.e. $\xi$ and $\dot{\gamma}$) in order to investigate different $Pe$ regimes. So, for example in our system, while using $\xi=10$ and too high values of shear rates, the temperature does not remain constant preventing us to explore high values of $Pe$ at this effective viscosity. However, building on the mapping, one can equivalently explore higher values of $Pe$ by increasing $\xi$ and using a smaller value of shear rate. Similarly, for example to study shear-induced crystallization may require very long simulation times at high solvent viscosities (and indeed we are not able to detect it within the duration of our simulations). However, a proper balance of the choice of $\dot{\gamma}$ and $\xi$ can be tailored for the specific needs of a given situation, allowing to explore the parameter space in a much more efficient way, without affecting the structure of the final state.

We also find that the mapping does not hold when we also consider MD simulations. Fig. 7(a) also shows that anisotropy is always higher in the absence of the solvent. This indicates that the presence of the Brownian motion mitigates the growth of the anisotropy and hence, different rheological response are observed when using MD and LD simulations.

Furthermore, in Fig. 7(b) we show the anisotropy persisting also in the steady state for $\xi=10^{2}$ for different values of Pe, corresponding to Set I and Set II state points in Fig. 6. For the larger $Pe$ value, the number of bonds in the steady state is small and the structure of the system is composed of many small groups. On the other hand, for the smaller Pe, where the system is closer to phase separation, the system is found in a new gel state, which is able to maintain a large amount of anisotropy, which propagates through the structure even well beyond the bond distance.

To understand the full rheological response of our gel, we also study the mechanical response versus strain for Set II and again we focus on the role played by the underlying microscopic dynamics. We plot in Fig. 8(a) the behavior of stress versus strain which clearly shows that the two LD simulations yield an identical behavior at large strains, while a dependence on the effective viscosity $\xi$ is observed for $\gamma\lesssim 10^{-2}$, where the effects of the Brownian dynamics balance the shear flow in a different way. On the other hand, the MD simulations display a much larger accumulated stress, with a yielding point that takes place at a different strain value, which is then followed by an oscillation before reaching a steady state with a stress larger than that found for LD simulations. To microscopically compare the three different simulation methods, we again consider the fraction of bonds and the fraction of unbroken bonds versus strain, shown in Fig. 8(b) and we find that these observables do not show a dependence on the shear conditions and on the solvent effects at all strains. However, the energy of the bonds reported in Fig. 8(c) does show significant differences between the MD and LD simulations. In particular, the MD data display an anticipated decrease of $E_{b}$, associated to their own yielding point, followed by an oscillation which reflects that observed in the stress. The long-time limit of $E_{b}$ in this case is different from that obtained in LD simulations, despite the potential energy being the same as imposed by our mapping and shown in Fig. 8(d). Thus, we conclude that, as expected in the absence of the solvent, the shear has a much stronger effect on the system: despite having the same potential energy and structure, the spatial configuration of the bonds is rather different as the stress tensor and $E_{b}$ show, indicating that the Brownian motion acts as an additional relaxation mechanism against shear flow.

III.6 Size effects

In this section, we investigate whether and how the size of the system affects the results reported so far. Stress versus strain curves are reported in Fig. 9 for three system sizes at two different values of $Pe$ for LD simulations with $\xi=10^{2}$. We find that the curves are superimposed onto each other at all investigated $N$, except for some small differences at short times. Likewise, the quantities $f_{b}$, $f_{u}$ and $E_{b}$ remain unchanged in the stationary state as shown in Fig. 9(b). Interestingly, in the transient regime, when the velocity profile has not yet matched the imposed one, we find that the bond energy slightly decreases on increasing $N$ (Fig. 9(c)) while the number of bonds remains constant. This occurs for $Pe>1$ in a small ${\gamma}$ window where shear flow acts on timescales comparable to those of particle diffusion, thus influencing the behavior of the system. In this regime bonds start to stretch as indicated by a decrease of $E_{b}$. However, on increasing the bond distance, particles feel the presence of a repulsive shoulder in their interaction thanks to which $E_{b}$ increases again as shown in the inset of Fig. 9(c) where it forms a sort of oscillation at small ${\gamma}$. Once the shear flow attains the imposed velocity profile, shear effects occur at a shorter time scale than the relaxation mechanisms and hence, the bonds stretch until they break. Such behaviour is more evident on increasing the system size due to the larger signal coming from bond stretching. The bond energy oscillations disappear for $Pe<1$, indicating that the shear flow do not affect considerably particle dynamics but can be considered as a perturbation of particle Brownian motion. The expression of different relaxation mechanisms in the transient regime has been recently reported in depletion gels Johnson, Landrum, and Zia (2018). However, in that case the oscillations in the bond energy were not reported. Finally we monitor the anisotropy distribution in (Fig. 9(d)) at the yielding point. We find that there is a small increase of anisotropy with $N$, but overall size effects are not very pronounced and do not qualitatively change the observed patterns. We also note that our results are also in qualitative agreement with simulations performed on much larger system sizes Moghimi, Jacob, and Petekidis (2017); Johnson, Landrum, and Zia (2018). Thus we can conclude that the results obtained for $N=2000$ particles are robust and qualitatively representative of larger system sizes for the considered properties, that can be considered “bulk” properties. Of course, we may expect some size dependence for the microscopic behavior and this will be addressed in future work.

IV Discussion and Conclusions

In this work we have investigated the rheological behavior of gels under steady shear with different numerical approaches. While several experimental and numerical studies have already addressed this problem for depletion-induced gels, which are formed out-of-equilibrium via an arrested spinodal decomposition, we have focused on equilibrium gels, obtained via the competition of depletion attraction and electrostatic repulsion. The main purpose of this work was thus to understand how the response of gels under shear is affected by the route by which the gel is obtained, and hence how it depends on the inter-particle potential between the particles. To reach this goal, we have performed three different types of simulations. Two sets of LD simulations were run for different friction coefficients, tuning in this way the effective viscosity of the implicit solvent. These were then compared to MD simulations where the presence of the solvent is neglected. In this way we could compare conditions which can describe Brownian colloidal motion with others which describe atomic dynamics. While the use of MD simulations in some cases is acceptable also for colloidal systems, for example when focusing on the slow dynamics onlyGleim, Kob, and Binder (1998), under shear the effect of the solvent becomes relevant. Although LD simulations do not take into account hydrodynamic interactions, they provide a more realistic approach than MD in order to compare with experiments. It is fair to say that our work represents one of the few examples providing a systematic comparison of the influence of microscopic dynamics under steady shear.

We have calculated the stress tensor for our gels with competing interactions, finding that they exhibit one yielding point prior to reaching a liquid-like steady state (Fig. 2). The strain at which the yielding point occurs is found to be independent of $Pe$. These results are in agreement with depletion-induced gels at comparable packing fractions and attraction strengths Koumakis and Petekidis (2011). Similarly to these studies, the stress overshoot displays a power-law dependence on Péclet number, but with an exponent that is higher than values obtained for depletion gels in both experiments and numerical simulations. We attribute this increment to a better resistance of the studied type gel to shear flow thanks to the long-range electrostatic contribution in the colloid-colloid potential. These findings suggest that the rheological response of a colloidal gel can be systematically varied by changing the effective interparticle interactions. Such a feature is very appealing for practical applications and for achieving a fine control of the rheological properties of a gel.

From the microscopic point of view, we have analysed the effect of strain on the gel structure for different $Pe$. In agreement with previous works on different types of gels Boromand, Jamali, and Maia (2017); Johnson, Landrum, and Zia (2018), we find that the initial bonds forming the gel are deformed after switching on the shear, and start to break only after the yielding point. However, they soon reorganize into a new network structure, whose characteristics depend on the Péclet number and also on the microscopic dynamics.

Also similarly to previous works Park, Ahn, and Wagner (2017); Johnson, Landrum, and Zia (2018); Jamali, McKinley, and Armstrong (2017), we find that the maximum anisotropy is reached at the yielding point. After this point, the initial gel structure is lost due to the bonds breakage, thereby decreasing the amount of stored stress and consequent anisotropy. At long times or large strains, the gel approaches a steady state which maintains some degrees of anisotropy at short length-scales Park, Ahn, and Wagner (2017); Jamali, McKinley, and Armstrong (2017). To quantify anisotropy, several works have reported the so-called fabric tensor Boromand, Jamali, and Maia (2017); Johnson, Landrum, and Zia (2018), even resolved along different directions Jamali, McKinley, and Armstrong (2017). Other works instead focus on the expansion of the pair correlation function Park, Ahn, and Wagner (2017); Moghimi, Jacob, and Petekidis (2017), including the present work. These different observables provide similar amount of information on the anisotropy, but the use of Im$g_{2}^{-2}\left(r\right)$ allows one also to obtain spatial resolution. In general, Im$g_{2}^{-2}\left(r\right)$ is found to have two roughly symmetric peaks, localized around the bond distance and of roughly maximum intensity at the yielding point. The balance between compressional (negative) and extensional (positive) peak depends on the employed $Pe$ as well as on the specific potential interaction. In our system, this holds for $Pe>10^{-1}$. For certain conditions, an asymmetric situation is found, where a large positive peak is accompanied by an almost absent negative peak. We find that this situation occurs for $Pe<10^{-2}$ and we suggest that this feature is a precursor of a fluid-to-crystal transition, which is obtained only in a narrow region of $Pe$ and at low enough solvent friction (at least within our simulation time window). Thus, the accumulation of particles in the extensional axis Im$g_{2}^{-2}\left(r\right)$ seems to be a pre-requisite for a transition from an amorphous to an ordered structure. Conversely, in those cases where the asymmetry is not found, crystallization is hampered because either the small crystal nuclei are not able to support the deformation induced by the shear flow or the used shear rate values are not able to dominate over the underlying Brownian dynamics. Interestingly, in Ref. Moghimi, Jacob, and Petekidis (2017), an increase of the extensional peak is also found for low $Pe$. Differently, from previous works on gels, our inter-particle potential includes a long-range repulsion, which is found to affect the anisotropy at the yielding point even at large length scales well beyond the bond distance.

Finally we also note that, while in general at low strains the anisotropy only acts at the direct level of bonds, at the yielding point it also affects large length scales well beyond the bond distance. The fact that the anisotropy is clearly observed also for $r>r_{bond}$ seems to be a distinctive feature of the present work, where the presence of an additional long-range Yukawa repulsion is responsible for the propagation of the anisotropy at larger scales, differently for depletion-like gelsPark, Ahn, and Wagner (2017); Moghimi, Jacob, and Petekidis (2017) but in agreement with numerical simulations of repulsive Yukawa glassesZausch and Horbach (2009).

To be able to connect results obtained with different types of simulations, we built on ideas borrowed from studies of the glass transition where states in equilibrium and in aging are connected through a time-temperature relation. In a similar fashion, we consider here a mapping between equilibrium and steady states under shear (so-called corresponding states) by considering state points at different $\dot{\gamma}$ and $T$ with the same potential energy (Fig. 5). This mapping confirms the intuitive expectation that, upon decreasing the effects of the solvent, higher values of shear rates are required in order to reach a similar steady state. We find that the static structure factors of the examined corresponding steady states are identical to each other (Fig. 6). In particular, the cluster peak observed in equilibrium, which is a distinctive feature of gels obtained by competing interactions Stradner et al. (2004), is destroyed by the presence of the shear in favour of a growth of $S\left(q\right)$ for $q\rightarrow 0$. This implies that the shear flow acts at large length scales, screening out the contribution of the long-range repulsion and hence, pushing the system closer to phase separation. Interestingly, these results are in agreement with those obtained by recent simulations of depletion-induced gels in the presence of hydrodynamic interactions Jamali, McKinley, and Armstrong (2017). In addition, the mapping puts forward the evidence that the response of the gel and the approach of the steady state is controlled solely by $Pe$ in LD simulations, providing identical results upon varying solvent conditions and shear rates. This confirms that the competition between Brownian motion and shear flow is the key parameter to control the behavior of the gel under shear. However, when we compare equal potential energy states under shear in the absence of the solvent, we find a stronger effect of the shear on the gel in terms of anisotropy and final steady state properties, signaling that the mapping does depend on the choice of the microscopic dynamics. These results can be useful in future simulation works because they allow an independent choice of solvent viscosity and shear rate respectively, maintaining the same $Pe$. This makes it possible to explore a wide range of shear conditions without affecting the resulting steady states. This is particularly valuable from the practical point to view to extend investigations at high $Pe$ for example ensuring temperature stability as well as to explore low $Pe$ where competing mechanisms inducing by shear may favour crystallization, as in the present case, or other underlying processes. It will be interesting to extend the mapping to different types of simulations where solvent is treated in a more accurate ways including hydrodynamics effects, e.g. using Dissipative Particle Dynamics Zausch et al. (2008); Ruiz-Franco et al. (2018) or Multi-Particle Collision DynamicsRipoll et al. (2005); Kapral (2008).

Finally, we have shown that our results are robust against the variation of the system size. Thus, the modification of the inter-particle interaction is able to provide a different rheological response with respect to the widely studied depletion gels. Our results call for experimental investigations probing the stress-strain behavior of gels with competing interactions, still missing so far. In addition, it would be interesting in the future to extend our study to different types of gels, particularly equilibrium gels resulting for patchy or limited-valence attractionsSciortino and Zaccarelli (2017) for which rheological investigations are scarce.

Acknowledgements.

JRF, EZ acknowledge support from ETN-COLLDENSE (H2020-MCSA-ITN-2014, Grant No. 642774). NG, EZ from ERC Consolidator Grant 681597 MIMIC. We thank Daniele Parisi for useful discussions.