Competing Time Scales Lead to Oscillations in Shear-Thickening Suspensions

TL;DR

This paper demonstrates that the interaction between instrument inertia and particle contact formation causes shear-rate oscillations in shear-thickening suspensions, revealing new insights into flow behavior and contact dynamics.

Contribution

It introduces a model linking instrument inertia and contact formation time scales to oscillations, validated by experiments, and provides a method to measure contact relaxation times.

Findings

Oscillations depend on the interplay of two time scales.

Oscillation frequency measures contact relaxation time.

Flow curve access is extended beyond conventional rheometry.

Abstract

Competing time scales generate novelty. Here, we show that a coupling between the time scales imposed by instrument inertia and the formation of inter-particle frictional contacts in shear-thickening suspensions leads to highly asymmetric shear-rate oscillations. Experiments tuning the presence of oscillations by varying the two time-scales validate our model. The observed oscillations give access to a shear-jamming portion of the flow curve that is forbidden in conventional rheometry. Moreover, the oscillation frequency allows us to quantify an intrinsic relaxation time for particle contacts. The coupling of fast contact network dynamics to a slower system variable should be generic to many other areas of dense suspension flow, with instrument inertia providing a paradigmatic example.

Figures

Equations

$$\hat{f}(\sigma )=\exp \left(-({\sigma }^{\ast }/\sigma {)}^{\beta }\right),$$

$${\phi }_J(\hat{f})={\phi }_m\hat{f}+{\phi }_{\mathrm{rcp}}(1-\hat{f}).$$

$$\eta (\phi ,{\phi }_J)={\eta }_s{\left[1-\phi /{\phi }_J(\hat{f})\right]}^{-2},$$

$$\frac{\mathrm{d}f}{\mathrm{d}t}=-\frac{\dot{\gamma }}{{\gamma }_0}\left[f-\hat{f}(\sigma )\right],$$

$${\rho }_{A}h\frac{\mathrm{d}\dot{\gamma }}{\mathrm{d}t}={\sigma }_E-\eta (f)\dot{\gamma },$$

$$\frac{\mathrm{d}\dot{\Gamma }}{\mathrm{d}\tau }={\Sigma }_E-{\eta }_r(f)\dot{\Gamma }\equiv g_1(\dot{\Gamma },f),$$

$$\frac{\mathrm{d}f}{\mathrm{d}\tau }=-\frac{\dot{\Gamma }}{ϵ}\left[f-\hat{f}({\eta }_r(f)\dot{\Gamma })\right]\equiv g_2(\dot{\Gamma },f),$$

$$ϵ=\frac{t_c}{t_i}\equiv \frac{{\gamma }_0{\eta }_s^2}{{\rho }_{A}h{\sigma }^{\ast }}.$$

$$ϵ<{ϵ}_c=-\dot{\Gamma }(\mathrm{d}\dot{\Gamma }/\mathrm{d}{\Sigma }_{\mathrm{E}}),$$

$$ϵ\left[\frac{{\Sigma }_E}{\dot{\Gamma }}-{\eta }_r(f)\right]=-\left[f-\hat{f}({\eta }_r(f)\dot{\Gamma })\right]\frac{\mathrm{d}\dot{\Gamma }}{\mathrm{d}f}.$$

$$\frac{\mathrm{d}f}{\mathrm{d}t}=-\left(\frac{\dot{\gamma }}{{\gamma }_0}+\frac{1}{t_r}\right)\left[f-\hat{f}(\eta (f)\dot{\gamma })\right],$$

Full text

Concentrated suspensions of non-Brownian (or granular) particles in a Newtonian solvent occur widely in industry, e.g., concrete Van Damme (2018) mine tailings Boger (2009), and chocolate Afoakwa et al. (2007). Their viscosity, $\eta$, often increases with either shear rate, $\dot{\gamma}$, or stress, $\sigma$ Denn and Morris (2014). Such shear thickening is now understood as a transition from a low-viscosity state, with lubricated particle contacts, to a high-viscosity state, with frictional contacts, as the repulsive force between particles is overcome at a critical onset stress, $\sigma^{*}$ Mari et al. (2015a); Lin et al. (2015); Comtet et al. (2017); Clavaud et al. (2017).

A phenomenological model of this process by Wyart and Cates (WC) Wyart and Cates (2014) predicts three types of flow curve, $\sigma(\dot{\gamma})$. At low volume fraction, $\phi$, a smooth increase connects two constant-slope (= viscosity) branches in $\sigma(\dot{\gamma})$, giving rise to continuous shear thickening (CST). Above a critical $\phi_{\rm DST}$, $\sigma(\dot{\gamma})$ becomes S-shaped, with a backwards-bending (d$\sigma$/d$\dot{\gamma}<0$) region connecting the two branches, giving discontinuous shear thickening (DST). Finally, above some $\phi_{m}$, $\sigma(\dot{\gamma})$ has no flowing upper branch and it bends back to $\dot{\gamma}=0$: the system shear jams at high stresses.

In the CST regime, suspensions flow steadily and homogeneously, and the WC model fits data from nearly-monodisperse hard-sphere systems Guy et al. (2015); Royer et al. (2016). In the DST regime, there is a jump in $\sigma$ as the imposed $\dot{\gamma}$ is increased Bender and Wagner (1996), while under imposed $\sigma$, neither homogeneous nor shear-banded steady flow is possible Hermes et al. (2016). There is no general model for the system-specific flow in this regime. Recent experiments Saint-Michel et al. (2018); Rathee et al. (2017) and simulations Chacko et al. (2018) focus on banding: spatial variation with high-$\sigma$ and low-$\sigma$ regions. Many systems also show large temporal fluctuations Lootens et al. (2003), which sometimes begin as ‘relaxation oscillations’: $\dot{\gamma}$ periodically drops precipitously to a nearly-jammed state Larsen et al. (2014); Bossis et al. (2017, 2018), with a frequency that increases with applied stress Nagahiro et al. (2013); Hermes et al. (2016). We extend the WC model to account quantitatively for such oscillations.

The key physics is the competition between the dynamics of frictional contact formation and a ‘system variable’, here the acceleration of the rheometer geometry Larsen et al. (2014); Bossis et al. (2017). When the ratio of the time scale of the former to that of the latter is small, we predict homogeneous flow with relaxation-type $\dot{\gamma}$ oscillations. Fitting the observed $\sigma$-dependence of the oscillation frequency reveals and quantifies an additional time scale, that intrinsic to the relaxation of frictional contacts after formation. Thus, rheometer geometry inertia, often considered an artefact, can be used to probe suspensions near jamming.

WC introduced a stress-dependent steady-state fraction of frictional contacts, $\hat{f}$; simulations Mari et al. (2014); Chacko et al. (2018) find

[TABLE]

with $\beta\lesssim 1$. $\hat{f}$ controls the jamming point, $\varphi_{J}$, at which $\eta\to\infty$. Here we use weight fractions, $\varphi$, due to the porosity of our main model system (cornstarch) Han et al. (2017). At $\sigma\ll\sigma^{*}$, $\hat{f}=0$, and the system jams at random close packing, $\varphi_{J}=\varphi_{\rm rcp}$. When $\sigma\gg\sigma^{*}$, $\hat{f}\to 1$ and the system jams at some $\varphi_{m}<\varphi_{\rm rcp}$. The WC model linearly interpolates between these two limits:

[TABLE]

The distance to jamming then determines $\eta$ via

[TABLE]

with $\eta_{s}$ the solvent viscosity. At a given weight fraction, $\hat{f}$ increases with stress, lowering the jamming point, which in turn increases the viscosity, $\eta(\sigma)\equiv\eta\{\varphi,\varphi_{J}[\hat{f}(\sigma)]\}$. At $\varphi_{\rm DST}<\varphi<\varphi_{m}$, the flow curve, $\sigma(\dot{\gamma})$, becomes $\mathsf{S}$-shaped, Fig. 1(a) ($\cdot\cdot$), with a region where d$\sigma$/d$\dot{\gamma}<0$. At $\varphi>\varphi_{m}$, shear jamming (SJ) is predicted, with the flow curve doubling back to $\dot{\gamma}=0$ when $\varphi_{J}(\sigma)=\varphi$, Fig. 1(a) ( - -).

In non-steady flow, e.g. on reversal Lin et al. (2015); Peters et al. (2016), the contact network of the suspension takes finite time to adapt. Thus, the fraction of frictional contacts at any one instant, $f(t)$, may differ from its steady-state value, $\hat{f}$, given by Eq. 1, towards which $f(t)$ relaxes. Simulations show that $f$ evolves with the accumulated strain Mari et al. (2015a); so following previous work Mari et al. (2015a); Nakanishi et al. (2012), we write

[TABLE]

with a characteristic strain, $\gamma_{0}$ (and $\dot{\gamma}\geq 0$). We now use Eqs. 2 and 3 to relate $\eta$ to $f$, rather than just $\hat{f}$.

External stress, $\sigma_{E}$, is applied through the system boundaries. In a rheometer, this is the ‘geometry’, which has far higher mass than the suspension for a typical gap height, $h$, between the boundaries Läuger and Stettin (2016), Fig. 1(e). In the steady state, $\sigma_{E}=\eta(\hat{f})\dot{\gamma}$, the sample stress, $\sigma$. When $\mathrm{d}\dot{\gamma}/\mathrm{d}t\neq 0$, force balance between the geometry and the sample gives

[TABLE]

with $\rho_{A}$ the geometry’s areal density. Equations 4 and 5, being two-dimensional, cannot capture aperiodic flow, but can account for $\dot{\gamma}$-oscillations and elucidate the physics of unsteady flow in shear-thickening suspensions.

Measuring time in units of the geometry inertial time scale, $t_{i}=\rho_{A}h/\eta_{s}$, we rewrite Eqs. 5 & 4 as:

[TABLE]

where $\tau=t/t_{i}$. Other dimensionless variables are shear rate, $\dot{\Gamma}\equiv\mathrm{d}\Gamma/\mathrm{d}\tau=\dot{\gamma}\eta_{s}/\sigma^{*}$; applied stress, $\Sigma_{E}=\sigma_{E}/\sigma^{*}$; viscosity, $\eta_{r}=\eta(f)/\eta_{s}$; sample stress, $\eta_{r}(f)\dot{\Gamma}=\eta(f)\dot{\gamma}/\sigma^{*}$; and, strain, $\Gamma=\eta_{s}^{2}/(\rho_{A}h\sigma^{*})$.

The time scale for contact network formation, $t_{c}=\gamma_{0}\eta_{s}/\sigma^{*}$, competes with the inertial time, yielding our key dimensionless parameter,

[TABLE]

When $t_{i}\gg t_{c}$, i.e. $\epsilon\ll 1$, Eqs. 6 and 7 form a singular autonomous system Hinch (1991), which may undergo a Hopf bifurcation to show relaxation oscillations as the control parameter $\Sigma_{E}$ is varied Baer and Erneux (1986).

For a given $\varphi$ and $\Sigma_{E}$, a fixed point occurs where the nullclines $g_{1}=0$ and $g_{2}=0$ intersect, Fig. 2(a). Analysing the Jacobian Scheinerman (1996), $\begin{pmatrix}\partial g_{1}/\partial\dot{\Gamma}&\partial g_{1}/\partial f\\ \partial g_{2}/\partial\dot{\Gamma}&\partial g_{2}/\partial f\end{pmatrix}$, shows that this fixed point is unstable if

[TABLE]

which, since $\epsilon>0$, requires d$\dot{\Gamma}$/d$\Sigma_{E}<0$, i.e. a backwards-bending flow curve, see Supplemental Material for derivation SM . Thus, the DST-boundary (d$\dot{\Gamma}$/d$\Sigma_{E}$=0) forms the lower boundary of our region of potential instability, Fig. 2(b). The upper boundary of this occurs at shear jamming, $\varphi_{J}(\Sigma_{E})=\varphi$, where the flow curve touches the vertical axis so that $\dot{\Gamma}=0$. Above this boundary, no flow is possible. Between these two boundaries, $\epsilon_{c}(\varphi,\Sigma_{E})$ peaks at $\epsilon_{c}^{\max}=2\times 10^{-5}$: instability may occur between DST and shear jamming whenever $\epsilon<2\times 10^{-5}$ 111$\epsilon_{c}^{\max}$ is found from the maximum numerically calculated $\epsilon_{c}$ in Fig. 2(b) with $\delta\varphi\approx 10^{-5}$ and $\delta\log\Sigma_{E}=3\times 10^{-4}$.. Physically, at such small $\epsilon$ (i.e. $t_{i}\gg t_{c}$), the suspension thickens before the geometry slows, so the sample stress rises, driving $\hat{f}$ higher and causing further thickening in a vicious cycle, pushing the system away from the steady state.

We now describe our dynamical system by phase-plane trajectories that depend parametrically on $\tau$. Consider the regime $\varphi_{\rm DST}\leq\varphi<\varphi_{m}$ with S-shaped flow curves, Fig. 1(a) (- -). The $f$-nullcline, Fig. 2(a), reflects the shape of the steady-state flow curve 222As there is a one-to-one dependence of $\Sigma_{E}$ to $f$ on the $f$-nullcline, $g_{2}=0$, see Eq. 7.. Equations 6 and 7 show that trajectories point inwards everywhere on the rectangle defined by $\dot{\Gamma}=0$, $\dot{\Gamma}=\dot{\Gamma}^{\dagger}$ (where the $\dot{\Gamma}$-nullcline intersects the $\dot{\Gamma}$ axis), $f=0$ and $f=1$, Fig. 2(a). However, trajectories point outwards on any infinitesimally-small loop around the fixed point if it is unstable. The Poincaré-Bendixson Theorem Scheinerman (1996) then predicts a limit cycle in the region depicted in Fig. 2(a) if $\epsilon<\epsilon_{c}$.

A numerically-calculated limit cycle after the onset of DST is shown in Fig. 2(c). To understand this cycle, divide Eq. 7 by Eq. 6 to obtain

[TABLE]

If $\epsilon\rightarrow 0$, Eq. 10 requires d$\dot{\Gamma}/$d$f\to 0$ (vertical lines) or $f\rightarrow\hat{f}(\eta_{r}\dot{\Gamma})$ ($f$-nullcline). If $\epsilon\ll 1$, starting at $(0,0)$, the system follows the $f$-nullcline ($g_{2}=0,~{}g_{1}>0$), Fig. 2(c), at a rate controlled by $t_{i}$ (Eq. 6). At $B$, the system jumps vertically to join the ‘upper branch’ of the $f$-nullcline at $C$. It now follows the ‘upper branch’ of the $f$-nullcline ($g_{2}=0,~{}g_{1}<0$) until it reaches $D$, where it drops vertically to $A$, and the process repeats: we have a limit cycle. As a consistency check, the ‘jump’ $BC$ and hence the limit cycle relies on $\dot{\gamma}$ not changing ($t\ll t_{i}$) as a large number of frictional contacts form and the suspension shear thickens ($t>t_{c}$), i.e. $\epsilon\ll 1$, as assumed.

At $\varphi>\varphi_{m}$, Fig. 2(d), the ‘jump’ from $B$ takes the suspension towards jamming, $\eta\to\infty$ at $C$, whereupon $\dot{\gamma}$ abruptly goes to zero, giving a horizontal ‘jump’ to $D$, from where the system drops back to $A$ on the $f$-nullcline, again giving a limit cycle. Note that the $CD$ part of our limit cycle probes our system close to jamming. Unlike in conventional steady-state rheology Fall et al. (2010), our system should remain homogeneous: the time needed to traverse $BCD$ is simply too short to allow finite particle migration.

To validate our model, we first characterized a shear-thickening suspension known to show oscillations Hermes et al. (2016). Cornstarch (Sigma Aldrich, particle diameter $\approx$14\text{,}\mathrm{\SIUnitSymbolMicro m} and polydispersity $\approx$ 40% from static light scattering Hermes *et al.* ([2016](#bib.bib13))) was dispersed into 50 wt.% glycerol-water ($\eta_{s}=$6\text{\,}\mathrm{m}\mathrm{P}\mathrm{a}\,\mathrm{s}). We used a TA Instruments DHR-2 with roughened parallel plates (radius $R=$20\text{,}\mathrm{m}\mathrm{m} and $h=$1.0\text{\,}\mathrm{m}\mathrm{m} for flow curves, $1.5\text{\,}\mathrm{m}\mathrm{m}$ for time dependence), Fig. 1(f). Rim shear rates, $\dot{\gamma}=\Omega R/h$, and apparent stresses, $\sigma_{E}=2\mathcal{T}_{E}/\pi R^{3}$, come from the applied torque, $\mathcal{T}_{E}$, and measured angular velocity, $\Omega$. Cornstarch particles are porous Han et al. (2017); so we quote weight fractions, $\varphi$, using freshly-prepared samples and monitoring reproducibility.

The WC model captures credibly the time-averaged flow curves of this system for $\varphi<\varphi_{\rm DST}$, Fig. 1(a) (- -), with $\varphi_{m}$, $\varphi_{\rm rcp}$, $\sigma^{*}$ and $\beta$ determined from fitting the steady-state rheology (see SM SM ), averaging 3 upsweeps at 10pts./decade from $0.1\text{\,}\mathrm{Pa}$ to fracture with $10\text{\,}\mathrm{s}$ average and $5\text{\,}\mathrm{s}$ delay, separately ensuring reversibility. At $\varphi>\varphi_{m}$, the WC model works until the predicted flow curve bends backwards, Fig. 1(a) ( - -). Up to this point, the flow is steady: $\dot{\gamma}$ is constant in time, Fig. 1(b). At higher stress, the flow starts to oscillate, Fig. 1(c), before becoming aperiodic, Fig. 1(d) Hermes et al. (2016).

The measured geometry moment of inertia, $I$, gave $\rho_{A}\equiv 2I/\pi R^{4}=$175\text{,}\mathrm{kg}\text{,}{\mathrm{m}}^{-2} [SM ](#bib.bib30). Imposed-rate experiments gave $\gamma_{0}=\mathcal{O}(10^{-1})$ [SM ](#bib.bib30). Thus, $t_{c}=$1.1\text{\times}{10}^{-4}\text{\,}\mathrm{s}, $t_{i}=$44\text{,}\mathrm{s}$$, and $\epsilon=2.7\times 10^{-6}$, far below the $\epsilon_{c}^{\max}=2\times 10^{-5}$ for observing instability when $\varphi_{\rm DST}<\varphi<\varphi_{\rm rcp}$. Solving Eqs. 6 and 7 numerically at $\varphi=0.47$ and $\Sigma_{E}=0.93$, we find relaxation oscillations quantitatively matching experiments with no free parameters, Fig. 3(a).

Next, we varied $\epsilon\propto\eta_{s}^{2}$ by increasing the solvent glycerol proportion, see SM for time-averaged rheology SM . For $\eta_{s}=$15\text{,}\mathrm{mPa}\text{,}\mathrm{s}, $\epsilon\sim 2\times 10^{-5}\gtrsim\epsilon_{c}^{\max}$, only damped oscillations in a narrow stress range were observed, Fig. [3](#S0.F3)(c)-(e). For $\eta_{s}=$75\text{\,}\mathrm{mPa}\text{\,}\mathrm{s}, $\epsilon\sim 3\times 10^{-4}\gg\epsilon_{c}^{\max}$, no shear-rate oscillations are seen at stresses and weight fractions in the DST-regime, Fig. 3(f). Oscillations could also be eliminated by only reducing $h$ (increasing $\epsilon\gtrsim\epsilon_{c}^{\max}$, see SM SM ), however large variation of $\rho_{A}$ is restricted by rheometer design. We also studied shear-thickening silica suspensions (diameter $4\text{\,}\mathrm{\SIUnitSymbolMicro m}$) Royer et al. (2016), in which oscillations have not been reported before. Experiments were performed using an Anton-Paar MCR302 in a parallel-plate geometry ($R=$20\text{,}\mathrm{mm}$,~{}h=$1.5\text{,}\mathrm{mm}) with $\rho_{A}=$400\text{\,}\mathrm{k}\mathrm{g}\mathrm{m}^{-2}, see SM for details SM . In 87 wt.% glycerol-water with $\eta_{s}=$151\text{,}\mathrm{mPa}\text{,}\mathrm{s}$$ and $\epsilon=1\times 10^{-4}$, no oscillations were seen, Fig. 3(g). Reducing $\eta_{s}$ to $3.4\text{\,}\mathrm{mPa}\text{\,}\mathrm{s}$ using a dimethyl sulfoxide-water mixture, giving $\epsilon=2\times 10^{-7}$, we found relaxation oscillations, Figs. 3(h) and 3(i). All our available data are consistent with the predicted $\epsilon_{c}^{\max}=2\times 10^{-5}$ for instability.

Figure 3(a) pertains to $\sigma_{E}$ at the onset of DST. As $\sigma_{E}$ increases beyond this point, the oscillation frequency, $\nu$, increases Hermes et al. (2016), and the agreement between model and experiment worsens, Fig. 4. As the system comes ever closer to jamming at each precipitous drop in $\dot{\Gamma}$, the strain-dependent ansatz for $f$-relaxation, Eq. 4, becomes increasingly ineffective. The predicted time taken to traverse $DA$ in the limit cycle, Fig. 2(d), is lengthened compared to reality (cf. slow onset in Fig. 4 inset).

We therefore infer the existence of an additional intrinsic, strain-independent, mechanism for relaxing $f$ towards its steady-state value Maharjan and Brown (2017) and modify Eq. 4 to read

[TABLE]

with a new relaxation time $t_{r}$. There are now two contact relaxation mechanisms, dependent on strain ($\propto\dot{\gamma}/\gamma_{0}$) or time ($\propto 1/t_{r}$). The latter dominates as $\dot{\gamma}\to 0$, near jamming, so the time taken for $DA$ shortens, decreasing the period of the limit cycle, as observed.

Fitting the $\nu(\sigma_{E})$ data with this new model, Fig. 4 333The initial frequency spike predicted by the model, associated with small-amplitude limit cycles not seen experimentally SM , is neglected., gives $t_{r}\approx$ $0.24(5)\text{\,}\mathrm{s}$ SM . Since $t_{c}/t_{r}=5\times 10^{-4}\ll 1$, the strain-dependent mechanism dominates away from jamming SM . Interestingly, $t_{r}\approx$ $0.24\text{\,}\mathrm{s}$ is comparable to the relaxation time for cornstarch grains pushed into adhesive contact, $\sim$ $0.5\text{\,}\mathrm{s}$, so that surface chemistry matters Oyarte Gálvez et al. (2017).

The mechanism we propose for relaxation oscillations in shear-thickening suspensions, depending on flow-curve shape and geometry inertia, appears generic. It is therefore perhaps a puzzle why such oscillations have not been more widely reported. One reason is the use of high-viscosity solvents, thus giving $\epsilon\gg\epsilon_{c}^{\max}$. More prevalent could be the breakdown of simple shear flow where surface tension no longer confines the particles as sample stress peaks Brown and Jaeger (2012); Guy et al. (2015) at $C$ in the limit cycle, Fig. 2(d), causing fracture Strivens (1976). With only two dynamical variables, lacking spatial variation, our model cannot capture such inhomogeneous flow. It nevertheless well captures the development of relaxation oscillations en route to aperiodic unsteady flows, which are widely seen Saint-Michel et al. (2018); Rathee et al. (2017); Lootens et al. (2003); Pan et al. (2015).

Our model generalized to Eq. 11 has allowed us to extract an intrinsic contact-relaxation time scale, $t_{r}$, which is difficult to access using other methods such as shear reversal Lin et al. (2015) or cessation Maharjan and Brown (2017) due to instrument artifacts. Instead, our method of accessing $t_{r}$ relies on modelling the coupling with one such artifact, viz., geometry inertia. $t_{r}$ becomes important in modelling the flow properties whenever the suspension comes close to jamming and the shear rate drops. With our protocol for extracting this relaxation time, future work should be able to clarify the underlying physical mechanism, which may include particle softness Ness and Sun (2016), surface chemistry James et al. (2018) or long-range repulsion Mari et al. (2015b).

Finally, coupling between fast frictional-contact-network dynamics and a slower ‘system variable’, and hence the resulting types of behavior, should be found in many types of dense suspension flow. Thus, for example, in vorticity banding, particle migration is slow Chacko et al. (2018); in micro-channel oscillations, rearrangement due to fluid permeation is slow Isa et al. (2009); Kanehl and Stark (2017); in the settling of a ball in a suspension, the ball’s inertial dynamics are slow von Kann et al. (2011). Interestingly, relaxation-type oscillations, with periodic bursts of brief near-jamming episodes, have been observed in the pipe flow of polymethylmethacralate particles Isa et al. (2009), the settling velocity of a ball in cornstarch von Kann et al. (2011) and the shear rheology of polystyrene particles Larsen et al. (2014) (compare especially data presented in the latter two cases with, e.g., our Fig. 2(c)). It is therefore possible, perhaps likely, that the kind of physics we have modelled may be relevant far beyond the data sets presented here.

The data plotted in this work are available from Edinburgh DataShare Dat (2019).