{'fraction_non_alphanumeric': 0.035349420244437484, 'fraction_numerical': 0.016671889689752428, 'mean_word_length': 4.5374826469227205, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A nonlinear dynamical system approach for the yielding behaviour of a viscoplastic fluid

A nonlinear dynamical system model that approximates a microscopic Gibbs field model for the yielding of a viscoplastic material subjected to varying external stress recently reported in [1] is presented. The predictions of the model are in a fair agreement with the microscopic simulations and in a very good agreement with the microstructural semi-empirical model reported in [2]. With only two internal parameters, the nonlinear dynamical system model captures several key features of the solid-fluid transition observed in experiments: the effect of the interactions between microscopic constituents on the yield point, the abruptness of solid-fluid transition and the emergence of a hysteresis of the micro-structural states upon increasing/decreasing external forcing.The scaling behaviour of the magnitude of the hysteresis with the degree of the steadiness of the flow is consistent with previous experimental observations.


Introduction
A broad class of materials exhibit a dual response when subjected to an external tress. For low applied stresses they behave as solids (loosely speaking they may deform but they do not flow) but, if the stress exceeds a critical threshold generally referred to as the "yield stress", they behave as fluids (typically non-Newtonian) and a macroscopic flow is observed. This distinct class of materials has been termed as "yield stress fluids". During the past several decades one could notice an increasing level of interest of both theoreticians and experimentalists in yield stress. This has a twofold motivation. From a practical standpoint, such materials have found a significant number of applications in several industries (which include food, cosmetical, pharmaceutical, oil field engineering, etc.) and they are encountered in daily life in various forms such as food pastes, hair gels and emulsions, cement, mud etc.. More recently, hydrogels which exhibit a yield stress have found a number of future promising applications including targeted drug delivery Han et al. (1997); Qiu and Park (2001), contact lenses, noninvasive intervertebral disc repair Hou et al. (2004) and tissue engineering Beck et al. (2007).
From a fundamental standpoint, yield stress materials continue triggering intensive debates and posing difficult challenges to both theoreticians and experimentalists from various communities: soft matter physics, rheology, physical chemistry and applied mathematics. The progress in understanding the flow behaviour of yield stress materials made the object of several review papers Nguyen and Boger (1992); Coussot (2014); Balmforth et al. (2014); Bonn et al. (2015). The best known debate concerning the yield stress materials is undoubtedly that related to the very existence of a "true" yield stress behaviour Barnes (1999); Barnes and Walters (1985). During the past two decades, however, a number of innovative improvements in the rheometric equipment made possible measurements of torques as small as 0.1nNm and rates of deformation as small as 10 −7 s −1 ). Such accurate rheological measurements proved unequivocally the existence of a true yielding behaviour ; Putz and Burghelea (2009) ;Bonn and Denn (2009) ;Denn and Bonn (2011). The physics of the yielding process itself on the other hand remains elusive. The macroscopic response of yield stress fluids subjected to an external stress σ has been classically described by the Herschel-Bulkely model Herschel and Bulkley (1926a,b): Here σ y is the yield stress,γ is the rate of shear, i.e., the rate at which the material is being deformed, σ is the macroscopically applied stress (the external forcing parameter), K is a so-called consistency parameter that sets the viscosity scale in the flowing state and N is the power law index which characterises the degree of shear thinning of the viscosity beyond the yield point.
In spite of its wide use by rheologists, fluid dynamicists and engineers, the Herschel-Bulkley model (and its regularised variants, e.g. Papanastasiou Papanastasiou (1987)) is in fact applicable only for a limited number of yield stress materials, sufficiently far from the solid-fluid transition, i.e. when σ > σ y , and in the conditions of a steady state forcing, i.e. when a constant external stress σ is applied over a long period of time.
Thixotropy, which may be loosely understood 1 as a time dependence of the rheological parameters which results from a competition between destruction and rejuvenation of the soft material units subjected to stress, is considered to be a major reason for the departure from this simple yielding picture Möller et al. (2006). It has been recently suggested that a number of difficulties concerning the yielding behaviour of yield stress materials could be solved if a clear distinction between thixotropic and non-thixotropic yield stress fluids is made . However, it has been shown recently that a clear departure from the Herschel-Bulkley behaviour can be observed even for non-thixotropic yield stress fluids such as the Carbopol gels particularly during either controlled stressed rheological tests Putz and Burghelea (2009); Moyers-Gonzalez et al. (2011);Weber et al. (2012); Divoux et al. (2011Divoux et al. ( , 2013 or hydrodynamically "simple" flow problems such as a creeping motion of a spherical object Putz et al. (2008), the unsteady laminar pipe flow Poumaere et al. (2014) or the emergence of the Rayleigh-Bénard convection, Kebiche et al. (2014) in a yield stress fluid heated from below.
To overcome these difficulties, several phenomenological macroscopic models have been proposed Dullaert and Mewis (2006); Quemada (1998aQuemada ( ,b, 1999; Coussot et al. (2002a,b); Roussel et al. (2004); Putz and Burghelea (2009);de Souza Mendes (2009; Dimitriou et al. (2013); Dimitriou and McKinley (2014); Blackwell and Ewoldt (2014) which have a general form 2 : The particular feature of these models is that they describe the temporal evolution of a microstructural parameter a(t) as a function of the applied stress and a number of parameters C 1 , ..., C m . Part of these parameters describe the kinetics of the destruction/restructuration of the material and are difficult (or impossible!) to measure. The rest of the parameters are measurable via adequate macroscopic rheological tests (flow ramps, oscillatory measurements, creep/relaxation tests etc.). Thus, as the applied stress is increased, a(t) varies smoothly from 1 (the entire volume of material is in a solid state) to 0 (the entire volume of material is in a fluid state) and the combined solid and fluid rheological responses are accounted for accordingly into a constitutive relation. Finally, the problem is reduced to a system of coupled equation which can be solved either analytically if the structural evolution equation is sufficiently simple (see for example the λ model proposed by Coussot in Coussot et al. (2002b), Roussel et al. (2004)) or numerically.
Such approaches have several clear indisputable advantages and have contributed significantly to our current understanding of yield stress materials: 1. As opposed to the Bingham and the Herschel-Bulkley models which predict an abrupt and discontinuous solidfluid transition when a jumps from 1 to 0 at a well-defined value of the applied stress which coincides with the yield stress σ = σ y , such approaches which directly account for the evolution of a(t) are able to predict a smooth (gradual) solid-fluid transition which is often observed in experiments Putz and Burghelea (2009 Blackwell and Ewoldt (2014) and creep/relaxation flows. 3. They can flexibly model the irreversibility of the deformation states when the applied stresses are increased/decreased around the solid-fluid transition. Moreover, they can quantitatively describe the experimentally observed rheological hysteresis and the dependence of its magnitude on the degree of steadiness of the forcing Putz and Burghelea (2009) as reflected by the area of the hysteresis of the dependenceγ = σ (γ), on the degree of the steadiness of the external forcing (i.e. how fast is the external stress varied about the solid-fluid transition). 4. From a practical standpoint, they are relatively easy to implement. Though able to model sufficiently complex rheological data (ranging from controlled stress/strain unsteady flow ramps, creep tests and oscillatory tests in a wide range of frequencies and amplitudes), such phenomenological macroscopic models do have a number of limitations: 1. As the functional dependence F in Eq. 2 is generally chosen on an intuitive basis rather derived from first principles, these models can teach little about the microscopic scale physics of the yielding process. 2. They typically involve a rather large number of parameters some of which are not directly and easily measurable and can be obtained only by fitting the experimental data. 3. Second, such models are not inherently validated from a thermodynamical standpoint. The second law of thermodynamics is not guaranteed to be held and such a validation is not always straightforward as it requires the derivation of a thermodynamic potential Picard et al. (2002); Bautista et al. (2009);Hong et al. (2008).
Bearing in mind that the yield stress behaviour originates from the presence of a "soft" microstructure which can only sustain a finite local stress prior to its breakdown, an alternative way of assessing the dynamics of the yielding process is to focus on the evolution of the micro-structural soft material units as the external stress is gradually increased past the solid-fluid transition and next to assess the macroscopic scale behaviour from the perspective of statistical mechanics. A thermodynamic approach for the deformation of a physical gel has been recently proposed by An and coworkers An et al. (2010). By using a mean field approach, they construct a free energy functional and describe the microscopic scale dynamics of the gel network as a function of the applied stress in terms of the monomer volume fraction and an internal connectivity tensor characterising the gel network. Peshkov and his coworkers have employed the irreversible mechanics and thermodynamics of two-phase continua to describe the yielding process but a comparison with the experimental data seems to require future developments Peshkov et al. (2014). de Bruyn (2013 has modelled the restricted diffusion of small tracer particles in heterogeneous media by performing Monte Carlo simulations in a site-percolation model and his results partially agreed with the experimental observations Oppong et al. (2006); Oppong and de Bruyn (2007).
More recently, we have proposed in R. Sainudiin and Burghelea (2015) a microscopic picture of yielding inspired from the Ising model of magnetisation of a ferromagnet Ising (1925); Stanley (1987). The model was built on an analogy between the local agglomerative interactions in terms of assembly/disassembly of neighbouring microscopic constituents in a yield stress material subjected to an external stress and the local ferromagnetic interactions in terms of spin up (+1) / spin down (−1) of neighbouring particles in a microscopic ferromagnetic system subjected to an external magnetic field. First, our approach is fundamentally probabilistic and formalises Gibbs fields as time-homogeneous and time-inhomogeneous Markov chains over the state space of all microscopic configurations and thus it is thermodynamically validated. Second, the model has solely two parameters. Though able to capture several key physical features of the solid-fluid transition, from a practical perspective the applicability of this model to describe rheological measurements and real flow problems is somewhat limited. First, the microscopic constituents are assumed fixed over a lattice while the external stress is varied and thus there exists no direct of equivalent of the rate of deformationγ. Second, as it is a statistical model, its implementation is not trivial and its usage is time consuming. The aim of the current contribution is to derive an approximation of the model in the form of a classical structural approach similar to Eq. 2 which, together with the appropriate constitutive equation could ultimately describe rheological measurements.
The paper is organised as follows. In Sec. 2 we provide a brief description of the microscopic Gibbs field approach introduced in R. Sainudiin and Burghelea (2015) with a particular emphasis on its main predictions. A differential equation approximating the expected solid fraction of the material in the Gibbs model is derived and analysed in Sec. 3. The results of the simulations according to the microscopic model and the expected trajectories from the approximating differential equation are presented in Sec. 4. The paper concludes in Sec. 5 with a discussion of the main findings, their impact and their possible implications and extensions.

The microscopic Gibbs field model
First we present a summary of the microscopic Gibbs Field model by Sainudiin et. al. R. Sainudiin and Burghelea (2015). We model an idealised yield stress material or viscoplastic fluid as a network of particles in an appropriate solvent that are capable of assembling by "forming bonds" or disassembling by "breaking bonds" with their neighbours when an external stress σ is applied. As already mentioned in the introduction, this approach is inspired by the Ising model of ferromagnetism and its advantage is that, once formulated, it can fully benefit from the already developed tools of Statistical Physics Stanley (1987); Landau and Lisfshits (1980). We investigate the model when the network of particles is the regular graph given by the toroidal two-dimensional square lattice. Let x s ∈ Λ = {0, 1} denote the phase at site s. Phase 0 corresponds to being yielded or ungelled and phase 1 corresponds to being unyielded or gelled. The phase at a site directly affects its connectability with its neighbuoring sites. We assume that only two gelled sites can be connected with one another.
We consider the following Gibbs potential over the two types of cliques: and where {s} is the singleton clique, s, r is the doubleton clique with r ∈ N s (the set of four nearest neighbouring sites of a given site s), σ ≥ 0 is the external stress applied, α ≥ 0 is the site-specific threshold, and β ∈ (−∞, ∞) is the interaction constant between neighbouring sites. It is important to note that the "stress" σ has actually the dimensions of an energy transferred to the lattice which is an obvious consequence of the fact that the sites in the lattice are fixed and there is no equivalent of a deformation γ of the lattice.
Using the Gibbs potentials defined above one can write the associated energy of a site configuration as: As expected, the energy function above is very similar to the Ising Hamiltonian, Ising (1925); Stanley (1987). Here the external stress σ is the analogue of an external magnetic field and the interaction parameter parameter β plays the role of the coupling between neighbouring magnetic spins.
The probability distribution of interest on the site configuration space X n is then where Z kT is the normalizing constant or partition function

Local Specification
Let the number of neighbours of site s that are in phase 1 be x N s := r∈N s x r . Then, E s (x), the local energy at site s of configuration x, is obtained by summing the Gibbs potential V C (x) over all C s, i.e., over cliques C containing site s, as follows Let (λ, x(S \ s)) denote the configuration that is in phase λ at s and identical to x everywhere else. Then the local specification is where In this work we focus on the effect of varying external stress σ at a constant ambient temperature, and therefore without loss of generality, we take kT = 1 and work with π(x) = Z −1 1 exp(−E(x)).

Markov chain on configuration space
We can think of the microscopic Gibbs field model as an X n -valued Markov chain {X(m)} ∞ m=0 , where X(m) = (X s (m), s ∈ S n ) and X s (m) ∈ Λ, in discrete time m ∈ Z + := {0, 1, 2, . . .}. Let the initial condition, X(0) = x(0), be given by the initial distribution δ x(0) over X n that is entirely concentrated at state x(0). Then the conditional probability of the Markov chain at time-step m, given that it starts at time 0 in state x(0), is where, the |X n | × |X n | transition probability matrix P α,β,σ over any pair of configurations (x, x ) ∈ X n × X n is and θ = θ(s, α, β, σ), is indeed a function of the site s and the three parameters: α, β and σ. By ||x − x || = 1 we mean that the configurations x and x differ at exactly site s, i.e., x s x s . Similarly, by ||x − x || = 0 we mean that the two configurations are identical, i.e., x = x or x s = x s at every site s ∈ S n . We can think of our Markov chain evolving according to the following probabilistic rules based on (9) and (10): • given the current configuration x, we first choose one of the n 2 sites in S n uniformly at random with probability n −2 , • denote this chosen site by s and let the number of bondable neighbors of s be i = N s (x) ∈ {0, 1, 2, 3, 4}, and • finally change the phase at s to 1, i.e., set x s = 1 with probability and set x s = 0 with probability 1 − p i .
We emphasise the dependence of p i on the parameters α, β and σ by p i (α, β, σ). This is plotted in Fig. 1 for different parameter values. From the plots it is clear that α is a location parameter while β controls the scale of the relative difference between p i 's.

Sufficient Configuration Statistics for Energy
Two informative singleton clique statistics of a configuration x(m) at time m are the number and fraction of gelled sites, given respectively by: Similarly, two informative doubleton clique statistics of a configuration x are the number and fraction of connected pairs of neighboring sites, given respectively by: When the configuration is a function of time m and given by x(m), then the corresponding configuration statistics are also functions of time and are given by: The energy of a configuration x can be succinctly expressed in terms of a(x) and b(x) as and therefore where β ∈ (−∞, ∞) andσ = σ − α ≥ −α for a given α ≥ 0. Since the energy of a configuration x, given n, only depends on its a(x) and b(x), we can easily visualize any sample path ( x(0), . . . , x(m) ) ∈ X m+1 n in configuration space that is outputted by either Algorithm 1 or Algorithm 2 presented in R. Sainudiin and Burghelea (2015) as the following sequence of (m + 1) ordered pairs in the unit square: Finally, we reserve upper-case letters for random variables. Thus, A(X), A(X), B(X) and B(X) are the statistics of the random configuration X. And the notation naturally extends to A(m), A(m), B(m) and B(m) when X(m) is a random configuration at time m.

The main predictions of the Gibbs field model for the yielding of a yield stress material
The Gibbs field model has been tested by monitoring the evolution of the volume fraction of the solid microscopic constituentsā(t) during increasing/decreasing linearly stepped stress ramp that mimics a rheological flow ramp, see Putz and Burghelea (2009); Moyers-Gonzalez et al. (2011); Weber et al. (2012). In the context of the Gibbs field model, the closest equivalent of the characteristic forcing time t 0 or the time the stress is maintained constant during a controlled stress stepped ramped is the average number of hits per lattice site h in the Gibbs algorithm.
Corresponding to each step of the ramp the stress was kept constant during a time t 0 and the volume fraction of un-yielded lattice sites was obtained via the Gibbs algorithm detailed in the Appendix of Ref. R. Sainudiin and Burghelea (2015). By varying the characteristic forcing time t 0 we could test the dependence of the microscopic yielding dynamics on the degree of steadiness of the external forcing and attempt a qualitative comparison with the experimental results referred to in Sec. 1.
In spite of its very limited number of parameters, we have shown that this model can capture several key features of the solid-fluid transition: 1. In the limit of a steady state external forcing (t 0 very large), the solid-fluid transition is reversible upon increasing/decreasing applied stresses only if the interaction parameter β does not exceed a critical threshold β c , see Figs. 10 (a,b) in Ref. R. Sainudiin and Burghelea (2015). Beyond this threshold a "genuine" micro-structural hysteresis is observed even in the asymptotic limit of a steady state forcing, see Figs. 10 (c,d) in Ref. R. Sainudiin and Burghelea (2015). 2. During unsteady flow ramps (t 0 finite) a micro-structural hysteresis is observed even in the absence of interaction (β = 0). The dependence of the magnitude of the hysteresis on t 0 is non monotonic, see Figs. 11 in Ref. R. Sainudiin and Burghelea (2015). At large t 0 it scales as a power law t −ξ 0 with ξ decreasing as the interaction parameter β increases which is qualitatively similar the experimentally observed rheological hysteresis behaviour Putz and Burghelea (2009); Moyers-Gonzalez et al. (2011); Weber et al. (2012) For highly unsteady stress ramps (small t 0 ) the magnitude of the hysteresis follows a log-normal correlation which is equally consistent with the experimental scaling measured during rheological tests Divoux et al. (2013).
To conclude this part, we have been able to qualitatively describe several main features of the solid-fluid transition experimentally observed for yield stress materials subjected to an external stress using a Gibbs field statistical approach with only two internal parameters. This motivates us to derive in the following an approximate continuous version of this model which is similar in form to the classical micro-structural approaches generally described by Eq. 2 but remains thermodynamically validated.

An Approximating Nonlinear Dynamical System
Here we derive a nonlinear first-order differential equation to asymptotically approximate E(A(t)), the expected fraction of sites in the solid phase, in continuous time t that is measured in units of n 2 discrete time-steps as the number of sites n 2 → ∞, under a fixed externally applied stress σ and fixed rheological parameters α and β.
First consider the discrete-time Markov chain {X(m)} ∞ m=0 of (11) and (12) and recall that X(m) is the random site configuration of the chain at discrete time m and A(m) = s X s (m) is the number of sites that are in phase 1. We will derive the approximation first for the case when β = 0 in (12) and then for the general setting of β 0.
3.1. Non-interactive case with β = 0 If β = 0 then the probability of the phase in site s at the next time-step is independent of the current configuration, i.e., Therefore, the probability that the total number of sites in phase 1 increases by 1 in one time-step is obtained by adding the probability of a transition from phase 0 to phase 1 over every uniformly chosen site s as follows: Dividing both sides of the equality that defines the above event by n 2 we get Pr A(m + 1)/n 2 = a(m)/n 2 + 1/n 2 A(m)/n 2 = a(m)/n 2 By an analogous argument we can obtain the probabilities for the remaining two possibilities Finally by considering the instantaneous rate of change of the expected fraction of sites in phase 1 we get the limiting differential equation approximation as The simple relationship above is mathematically very similar to the so called "lambda-model" introduced in Coussot et al. (2002a,b) with the remark that we consider the stress σ as a forcing parameter rather than the rate of deformation.
Given the initial condition a(0) = a 0 , the analytic solution is with only one asymptotically stable fixed point Thus, a(t) in the above differential equation is the expected fraction of sites in phase 1 at time t in the limit of an infinite toroidal square lattice with |S n | = n 2 → ∞ and a realization of the continuous time Markov chain {A(t)} t≥0 is a(t).
Since β = 0, the probability of a site being in a given phase is independent of the phases of its neighboring sites. Thus, we can obtain b(t), the expected fraction of bonds, by simply multiplying a(t), the probability of finding a randomly chosen site in phase 1, by itself, i.e., 3.2. Interactive case with β 0 If β 0 then the probability of site s being in phase 1 at time m+1 depends on the configuration of the neighboring sites of s at time m through X N s (m) = r∈N s X r (m), the number of neighboring sites of s in phase 1 at time m.
Thus the probability that the phase changes from 0 to 1 in one time-step at site s given that a(m) is the total number of sites in phase 1 at time m is Since there are 4!/((4 − i)!i!) distinct neighborhood configurations with i of the four nearest neighbors of site s in phase 1, we can make the following binomial approximation for Pr X N s (m) = i X s (m) = 0, S = s, A(m) = a(m) in the above expression and obtain Therefore, the probability that the total number of sites in phase 1 increases by 1 in one time-step is obtained by adding the probability of a transition from phase 0 to phase 1 over every uniformly chosen site s as follows: Dividing both sides of the equality that defines the above event by n 2 we get Pr A(m + 1) = a(m) + 1/n 2 A(m) = a(m) By an analogous argument we can obtain the probability that A(m + 1) decreases by 1/n 2 as Using the same limiting approximation in the previous Section we can obtain the following differential equation This simplifies after factoring and extracting coefficients of a as follows: a(t) = p 0 − (4 p 0 − 4 p 1 + 1)a + 6 (p 0 − 2 p 1 + p 2 )a 2 −4 (p 0 − 3 p 1 + 3 p 2 − p 3 )a 3 +(p 0 − 4 p 1 + 6 p 2 − 4 p 3 + p 4 )a 4 .
We can understand (19) directly as a quartic polynomial in a whose coefficients are given by an alternating binomial series corresponding to the increase and decrease in a based on a combinatorial averaging over the transition diagram of site configurations at the four nearest neighbors of a given site. Next we characterize the qualitative asymptotic dynamics of the above nonlinear differential equation which reduces to the differential equation (16) if β = 0 and thereby p = p 0 = p 1 = p 2 = p 3 = p 4 .
If the externally applied stress σ is beyond α by 2β, i.e.
Therefore, alongσ = 2β our (19) is really just a cubic function of a as opposed to a quartic. The discriminant of this cubic is (19), and it takes negative values whenσ ∈ (0, 2.589145) (giving one real and two complex conjugate roots), takes positive values whenσ < 0 andσ > 2.589145 (giving three distinct real roots) and takes 0 whenσ ∈ {0, 2.589145} (giving three multiple real roots) as shown in Fig. 2.  In the left panel of Fig. 3, we present three different stability scenarios for the fixed points of equation (19) in the (σ, β) plane: (i) In the blue shaded region the right hand side of equation (19) has four real roots and only one of them is in [0, 1], this fixed point is stable. (ii) In the yellow region, starting at point (2.589145, 1.2945725), we have four distinct real roots with three of them in [0, 1]. Only one of the three distinct real roots is an unstable fixed point while the other two roots are stable fixed points. This naturally corresponds to a family of pitch-fork bifurcations and the associated hysteresis depending on where the system is initialised from. (iii) The unshaded region in the left panel of Fig. 3 corresponds to the parameter space where the quartic discriminant ∆ 4 is negative and thus implying the existence of two real roots (with one of them in [0, 1], stable fixed point) and two complex conjugate roots.
The real roots and their derivatives over each (σ, β) in a grid of parameter values from [−8, 12] × [−4, 4] were obtained through interval analytic methods using Hofschuster and Krämer (2003). Figure 4 shows the set of fixed points a * of the dynamical system as a function of (σ, β). The parameter space corresponding to the central shaded region of Fig. 3 containing the line β =σ/2 is evident in Fig. 4 with three fixed points in [0, 1]. The pitch-fork bifurcations along the planeσ = 2β or β =σ/2 determined by the non-negative sign of the cubic discriminant of Fig. 2 along the black line in Fig. 3 is displayed to highlight the dynamics with one unstable fixed point at 1/2 and two other stable fixed points that are equidistant on either side of 1/2.
We are interested in varying the externally applied stress σ for a given material characterized by fixed rheological parameters α and β. This amounts to varyingσ for a fixed β since the fixed α is absorbed intoσ = σ − α. The asymptotic dynamics when we apply a constant external stress for a long period of time are given by the fixed points a * in Fig. 4. We study such stress-dependent behavior from the Gibbs sampler and compare it with ODE approximation in the next Section. Note that the ODE model for β 0 is only in qualitative agreement with a(t), the expected volume fraction of the unyielded material at time t. This is because we are ignoring the dependent statistic b(t), the expected fraction of bonds or pairs of neighboring unyielded material at time t. Despite this simplification, as we will see in Sec. 4, there is qualitative agreement between the ODE and the Gibbs simulations. Furthermore, an admittedly ad hoc correction of the ODE through a translation of the vector field by (α 0 , β 0 ) even improves the quantitative approximation. We postpone a formal quantitative approximation of the ODE using perturbation theoretic methods to the future and focus here on obtaining insights from the Gibbs sampler that is in qualitative agreement with the ODE approximation. Figure 4: The fixed points a * as a set-valued function of the parametersσ = σ − α and β. The blue, black and azure points are the stable fixed points while the red and green points are the unstable fixed points of the system. There is a pitch-fork bifurcation alongσ = 2β that starts at (2.589145, 1.2945725) where the fixed point at 0.5 becomes unstable with two stable fixed points on either side.

Results
In this Section we mainly obtain various insights about the macroscopic behaviour of our model based on Monte Carlo simulations from Algorithms 1 and 2 in R. Sainudiin and Burghelea (2015), and make some comparisons with the approximating nonlinear ODE model of Sec. 3.

Comparison between Microscopic model and ODE approximation under varying stress
The energy of X(t), the random site configuration at time t, depends on two of its highly correlated statistics: A(t), the random fraction of gelled sites at time t, and B(t), the random fraction of connected sites at time t. One of our primary interests is to study A(t) and B(t) as X(t) is under the influence of time-varying externally applied stress σ(t).
Using Monte Carlo simulations from Algorithm 2 in R. Sainudiin and Burghelea (2015) of the time-inhomogeneous Markov chain {X(m)} Mh m=0 , under a time-dependent stress σ ramp, we can obtain multiple independent trajectories of A(σ), the fraction of gelled sites as a function of the external stress σ. This is to emulate conditions of an unsteady forcing during macroscopic rhelogical measurements. In the following, h is the average hits per site in the Gibbs sampler algorithm and we define it also as the characteristc forcing time t 0 for the stress ramp in our ODE simulations. We set h = 1000 in order to reach steady state for each value of σ. In Figure 5, the trajectories are shown as thin lines and the curves for the ODE approximation have the symbol on them. Note the reversibility of the response of the material when β ∈ {0, 1} (top row of Figure 5) upon increasing/decreasing applied stresses. The microscopic model and the ODE approximation quantitatively agree quite well when β < β c (β c ≈ 1.3), the threshold for three fixed points in [0, 1] for the ODE model. As we increase β beyond the aforementioned threshold β c we see that irreversible behaviour in the material appears and the comparison between the two models (discrete and continuous) is only qualitative in nature. This is due to the fact that our ODE approximation only models a, instead of modelling the dependent pair (a, b) that is sufficient for the energy, see Sec. 2.3. This effect can also bee seen if we compare the right panel of Fig. 3 with Fig. 8(c) in R. Sainudiin and Burghelea (2015). Clearly the light region of Fig. 8c (R. Sainudiin and Burghelea (2015)) corresponds to the yellow region where hysteresis is always present. The main discrepancy is the value of β c . In our ODE approximation, the calculated value is β c ≈ 1.3, on the other hand, from our Gibbs sampler simulations β GS c ≈ 1.5. As mentioned above this difference is due to the fact that in our approximation we disregards all bond interactions between neighbours.
As a qualitative remark one can note that even in the presence of strong interactions β > β c , both models predict an increase of the steepness of the solid fluid transition (defined as the slope of the dependenceā(t) on σ around the point whereā ≈ 1/2).

Comparison between model by Putz and Burghelea Putz and Burghelea (2009) and ODE approximation
In this section we will consider the model developed by Putz and coworkers in Putz and Burghelea (2009); Moyers-Gonzalez et al. (2011). This model is phenomenological in the sense that, unlike the Gibbs field model presented in Sec. 2 it is not derived from first principles. In this type of modelling one mimics the behaviour of the microstructure through the definition of a macroscopic structural variable with range in [0, 1], where 0 means completely unstructured or fluid and 1 means completely structured or solid. The structural variable a p satisfies a kinematic equation and usually depends explicitly on the stress and/or rate of strain. In the case of Putz and Burghelea (2009) where k r is the rate of recombination of micro-structural units, k d is the rate of destruction of the solid phase, σ y is the yield stress and w is a constant that controls how steep the change in the microstructure from solid to fluid and fluid to solid is.
In Figure 6 we present the simulations of equations (19) and (20) for three characteristic forcing times t 0 . As expected we have very good agreement between the models. This could be considered as a qualitative "proof " that the phenomenological models can actually approximate the behaviour of the microscopic models derived from first principles.

Determination of the yield point in the limit of a steady state forcing
A reliable estimation of the yield point is important to many practical applications involving yield stress materials.This is typically done by fitting steady state rheological measurements with models with various degrees of complexity ranging from the mathematically simple and classical Herschel-Bulkley correlation up to structural models. Thus, it appears natural to attempt in the following to obtain an estimate of the yield point for the case of a steady state forcing from the nonlinear dynamical system model presented herein.
To get an approximation for the yield point σ y during a steady state forcing process we will make the assumption (well supported by the results presented in Figs. 5, 6) that, corresponding to the yield point, the absolute value of the slope of the dependenceā * (σ) passes through a maximum: For simplicity, let us focus first on the non-interacting case, β = 0. From Eq. 17 on can readily show that the condition 21 reduces to σ y = α. Thus, in the non-interactive case, the yield point obtained during a steady state stressing practically coincides with the site specific threshold α of the Gibbs field model.
We now consider the interactive case β 0. To a leading order inā * and assuming that around the yield point a * ≈ 1/2 it can be shown using Eq. 19: The implicit dependence of the approximate yield stressσ y on the interaction parameter β may be obtained by solving numerically dā * ) dσ = 0. The result is presented in Fig. 7. For interactions weaker than the critical threshold β c , the apparent yield stress scales asσ y = σ − α = β (the dash-dotted line in Fig. 7). Beyond this threshold, the scaling becomes steeper,σ y = σ − α = 2β (the dashed line in Fig. 7). To conclude this part, the yield stress assessed via steady state controlled stress ramps is (according to our model) expected to depend linearly on both the site specific threshold α which may be intuitively understood as a measure of the strength of the microscopic constituents of the fluid and the strength β of their interaction and the slope of this behaviour switches when the strength of the interaction passes through the threshold β = β c .
In Fig. 8 we investigate the dependence of right hand side of eqn. 22 with respect toσ (left panel) and with respect to β (the right panel).
Regardless the value of the interaction parameter the stress dependence of the slope passes through a local maximum marked by a full symbol in Fig. 8(a). As previously explained, this may be considered as an indicator of the yield point. While β increases, the location of this maximum shifts towards larger stress values as already illustrated in Fig. 7. The value of this maximum slope increases monotonically with beta, the dashed line in Fig.8(a). As we approach β c the slope diverges Fig. 8(b). This is consistent with the fact that our steady solution becomes discontinuous as a function ofσ. Recall that we have a pitchfork bifurcation with stable fixed points {0, 1}, hence the value of σ y is not unique and depends on the initial condition.

Effect of the characteristic forcing time t 0 on the micro-structural hysteresis
In many practical situations the forcing time necessary to reach steady state in a given material cannot be reached. Perhaps one of the simplest such situation is that of a millimetre sized spherical object falling slowly (speeds of order of millimetres per second) in a Carbopol gel as investigated in Ref. Putz et al. (2008). The time scale associated to such motion is simply the ratio of the object's size to its speed and, in the case of the experiments reported inPutz et al. (2008) was of the order of a second or shorter. To understand such hydrodynamic problems steady state rheological measurements do not suffice. In Ref. Putz and Burghelea (2009) the rheological response of Carbopol gel was characterised during unsteady controlled stress linear ramps for various values of the holding time t 0 per stress values. A rheological hysteresis (which could, at least partially and qualitatively, explain the emergence of a fore-aft symmetry breaking of the flow pattern measured for the falling sphere experiment) was systematically found and its magnitude scaled as t −ξ 0 with ξ ≈ 1, see Fig. 11 in Putz and Burghelea (2009). The purpose of this section is to study the dependance of the area of the hysteresis of the micro-structural states observed upon increasing/decreasing forcing as a function of characteristic forcing time t 0 , compare the results with the predictions of the Gibbs field model R. Sainudiin and Burghelea (2015) and with the predictions of the predictions of the structural model by Putz and Burghelea Moyers-Gonzalez et al. (2011) as well as with the experiments Putz and Burghelea (2009);Divoux et al. (2013).
For this purpose, we have first run simulations using the nonlinear dynamical system described in Sec. 3 corresponding to several linear increasing/decreasing stress ramps for various values of the holding time per stress value t 0 and several values of the interaction parameter β. For each case we have calculated the area of the micro-structural hysteresis of the dependenceā(σ). The results obtained of the hysteresis area on the characteristic time t 0 obtained from the nonlinear dynamical system approach are represented in Fig. 9 as open symbols. For comparison, we calculated the same dependence by running the Gibbs field model for the same values of the interacting paramater and a number of hits per site h that matches t 0 . The results are represented in Fig. 9 as full symbols. The results obtained by the two approaches are in a good qualitative agreement: in both cases the magnitude of the microstructural hysteresis depends in a non-monotonic fashion on the degree of steadiness of the external forcing. In the steady state limit of large t 0 (h), a power law scaling in the form t −ξ 0 is observed. This finding is consistent with the experimental observations of a rheological hysteresis for a Carbopol gel subjected to an increasing/decreasing linear stress ramp

21
in a roughened plate-plate geometry Putz and Burghelea (2009). At a quantitative level, the agreement between the results obtained via the two approaches is only partial. The biggest differences in the scaling behaviour (see the inset in 9) are observed for the non-interactive case, β = 0 (the circles in Fig. 9). Eventhough that for the case β = 0 the approximation is exact, one should note that the Gibbs model is discrete. We expect that as n → ∞ (recall that n 2 is the number of sites in our lattice) the rate of decay in the hysteresis will converge to the value of the continuous model. As the interaction parameter β is gradually increased, the quantitative agreement between the two approaches improves (the circles, the squares and the rhombs). In all cases, however, the magnitude of the micro-structural hysteresis scales as a power law with the characteristic forcing time, t −ξ 0 , see the inset in Fig. 9. It is noteworthy that the values of the scaling exponent ξ are of the same order of magnitude as the ones measured experimentally in Putz and Burghelea (2009). It is equally interesting to note that both approaches predict a decrease of the scaling exponent with the interaction parameter. This finding is fully consistent with the experimental fact that, in the case of strongly interacting systems (e.g. laponite and bentonite suspensions) the rheological measurements exhibit a hysteresis (and, consequently, their reproducibility during subsequent tests is poor) even in the asymptotic case of a steady state forcing (very large waiting times t 0 ). More recently, a large hysteresis was systematically observed even for very large values of t 0 during controlled stress ramps performed with a suspension of a micro-alga with an electrically charged cellular membrane, Souliès et al. (2013).
In the case of a unsteady forcing i.e. small t 0 the simulations based on both approaches reveal a local maximum that shifts slightly towards larger values of t 0 as the interaction parameter β increases. This finding is consistent with the experimental results obtained by Divoux and his coworkers for three materials Divoux et al. (2013): a laponite suspension, for a carbon black suspension and mayonnaise. For a Carbopol gel, however, this local maximum was not observed neither in Putz and Burghelea (2009) (see Fig.11 therein) nor in Divoux et al. (2013) (see Fig. 3 (c) therein) because, most probably, it occurs at characteristic forcing time t 0 too small to be probed experimentally. An additional insight into the physical nature of this non monotonic behaviour was given in R. Sainudiin and Burghelea (2015) by showing that for values of t 0 below the maximum the lattice is only partially yielded corresponding to the maximal value of the stress reached during the ramp. Last, we compare the predictions of the nonlinear dynamical system model with the predictions of the micro-structural model by Putz and Burghelea (the dotted line in Fig. 9), Putz and Burghelea (2009). One can note a fair agreement with the result obtained for the non-interactive case β = 0.

Conclusions, outlook
We have presented a nonlinear dynamical system (ODE) approach for the solid fluid transition of a yield stress material subjected to an external stress that approximates the microscopic Gibbs field formulated from first principles introduced in Ref. R. Sainudiin and Burghelea (2015).
In spite of some quantitative differences mainly due to the fact that the ODE approximation does not properly account for the statistics of bonds between neighbouring microscopic constituents, both the ODE and the Gibbs field approach predict several key features of the solid-transition.
First, the transition is generally irreversible upon increasing/decreasing forcing an a micro-structural hysteresis is systematically observed, Fig. 5. A reversible transition may be observed solely in the non-interacting case β = 0 and in the limit of a steady state forcing (the top left panel in Fig. 5). By a systematic analysis of the stability of the fixed points of the nonlinear dynamical system we could show that a genuine hysteresis will be observed even in the asymptotic limit of steady forcing if the interaction parameter exceeds the threshold β c ≈ 1.3, (bottom panel in Fig. 5 ). The magnitude of the micro-structural hysteresis depends on both the level of interactions between the microscopic constituents and the degree of steadiness of the external forcing t 0 . In the limit of slow forcing the hysteresis decays as a power law t −ξ 0 , Fig. 9 and the power law exponents decreases with increasing strength of the interactions, the inset in Fig. 9. The monotonic decrease of the scaling exponent with β implies that, when strong interactions among the microscopic constituents are present, a strong irreversibility of the micro-structural states upon increasing/decreasing stresses even in the asymptotic limit of a steady state forcing.
Second, the abruptness of the solid fluid gradually increases with increasing interaction parameter β, Fig. 5. Third, we remark that the Gibbs Simulations as well as the approximating ODE are in qualitative agreement with the simple phenomenological model for the micro-structural hysteresis proposed in 6.
Fourth and equally important from a practical perspective, our model allows one to estimate the yield point and monitor its behaviour as a function of the interaction parameter, Fig. 7. A linear dependence is found and its slope changes corresponding to the critical point β = β c . When β > β c the value of the yield stress is no longer unique. This due to the fact that the process is no longer reversible. The steepness of the yielding transition diverges at the critical point β = β c , Fig. 8.
In closing, we believe there are several future directions worth pursuing. At a theoretical level, a more quantitative comparison between these models aimed at highlighting their differences may be useful. Ideally, perturbation theoretic methods should be used to improve the quantitative agreement between the nonlinear ODE model and the stochastic trajectories as opposed to the ad-hoc translations of the vector field done in this study. A more detailed model that simultaneously represents the fraction of gelled sites and the fraction of bonds in one dependent system would provide a better quantitative and qualitative approximation of the correlated site percolation model. Another interesting extension of our model could involve allowing for solvent effects through a model akin to correlated site-bond percolation of (Stauffer et al., 1982, Sec. D.II., p.136) but with our focus on external stress as opposed to temperature. In such a model we have an additional parameter that allows for a site to be occupied by a monomer with probability φ and by the solvent with probability 1 − φ.
From the more practical standpoint of the rheologist, it would be interesting to couple the nonlinear dynamical system approach to an appropriate elasto-viscoplastic constitutive relation and attempt either fitting experimental data or modelling industrially relevant non rheometric flow problems.