# Gel and glass transition in fragile colloidal clays

R. Angelini${}^{\textsf{\footnotesize{\mbox{?},\penalty 1000\mbox{?}}}}$, G.
Ruocco${}^{\textsf{\footnotesize{\mbox{?},\penalty 1000\mbox{?}}}}$, B.
Ruzicka${}^{\textsf{\footnotesize{\mbox{?},\penalty 1000\mbox{?}}}}$

(Received July 9, 2019)

###### Abstract

Dynamic Light Scattering (DLS) measurements were performed on colloidal
suspensions of Laponite® at different concentrations in the range
$C_{\text{w}}=(1.5{\div}3.0)$%. The slowing down of the dynamics induced by
aging was monitored by following the temporal evolution of autocorrelation
functions at different concentrations towards the gel and glass transition.
Exploiting analogies with supercooled liquids approaching their glass
transitions, an Angell plot for the structural relaxation times was drawn.
Finally, the fragility of Laponite® suspensions at different concentrations,
in different solvents, at two salt concentrations and with the addition of a
polymer was reported and discussed.

Key words: glass transition, colloids, Laponite®, gels, fragility

PACS: 64.70.Pf, 82.70.Dd

###### Abstract

Âèìðþâàííÿ äèíàìíîãî ðîçñþâàííÿ ñâòëà áóëè ïðîâåäåí íà êîëî¿äíèõ ñóñïåíçÿõ
Ëàïîíòó® ïðè ðçíèõ êîíöåíòðàöÿõ â îáëàñò $C_{\text{w}}=(1.5{\div}3.0)$%.
Ñïîâëüíåííÿ äèíàìêè, íäóêîâàíå ñòàðííÿì, ìîíòîðèëîñü åðåç àñîâó åâîëþöþ
àâòîêîðåëÿöéíèõ ôóíêöé ïðè ðçíèõ êîíöåíòðàöÿõ â íàïðÿìêó ïåðåõîäó â ñòàí ãåëþ
òà ñêëà. Âèêîðèñòîâóþè àíàëîã¿ ç ïåðåîõîëîäæåíèìè ðäèíàìè ïðè íàáëèæåíí äî ¿õ
ïåðåõîäó â ñòàí ñêëà, áóâ íàìàëüîâàíèé ãðàôê Àíäæåëà äëÿ àñó ñòðóêòóðíî¿
ðåëàêñàö¿. Ïîâäîìëÿòüñÿ òà îáãîâîðþòüñÿ êðèõêñòü ñóñïåíçé Ëàïîíòó® ïðè ðçíèõ
êîíöåíòðàöÿõ, â ðçíèõ ðîçèííèêàõ, ç äâîìà êîíöåíòðàöÿìè ñîë òà ç äîäàâàííÿì
ïîëìåðó.

Ключов слова: ïåðåõä â ñòàí ñêëà, êîëî¿äè, Ëàïîíò®, ãåë, êðèõêñòü

##  1 Introduction

In the last decades the glass transition has been largely investigated in
complex systems such as molecular liquids and colloids [1, 2]. These studies
have shown how the behaviour of these systems and their dynamical arrest are
characterized by several analogies. The glassy state can be obtained driving
them out of equilibrium by decreasing temperature $T$ in the case of molecular
liquids and by increasing volume fraction $\phi$ or waiting time
$t_{\text{w}}$ for colloidal systems. Changing these control parameters the
dynamics slows down enormously and a dramatic increase of the characteristic
relaxation time ($\tau(T)$, $\tau(\phi)$, $\tau(t_{\text{w}})$) is observed up
to the glassy state.

The growth of several orders of magnitude of both structural relaxation time
and viscosity with decreasing temperature when the glassy state is approached,
has been deeply investigated in supercooled molecular liquids [3], and the
rapidity with which these quantities increase has been quantified by the
fragility index $m$ defined as:

$m=\left[\frac{\partial\log\tau}{\partial(T_{\text{g}}/T)}\right]_{T=T_{\text{g}}},$
(1.1)

where $\tau$ is the relaxation time and $T_{\text{g}}$ is the glass transition
temperature. The “universal” Angell plot that classifies supercooled molecular
liquids in two classes of materials has been derived. In “fragile” liquids,
typically consisting of molecules interacting through non-directional,
noncovalent interactions, the relaxation time and viscosity are highly
sensitive to changes in $T$ with a super-Arrhenius (Vogel-Fulcher-Tammann)
behaviour:

$\tau=\tau_{0}\mathrm{e}^{\bar{D}/(T-T_{0})},$ (1.2)

where $T_{0}$ is the temperature at which $\tau$ diverges.

In “strong” liquids, typically characterized by three dimensional network
structures of covalent bonds, the above relation is simplified considering
$T_{0}=0$ and the systems have a much lower $T$ sensitivity with an Arrhenius
dependence:

$\tau=\tau_{0}\mathrm{e}^{D/T}$ (1.3)

with the constant $D$ equal to $E/k_{\text{B}}$ where $E$ is the activation
energy and $k_{\text{B}}$ is the Boltzmann constant [4].

This unifying concept describes and classifies the behaviour of supercooled
molecular liquids [3, 5, 6]. Moreover, a large variation in fragility as that
observed in the $T$ dependence of molecular liquids has been also found in the
case of the $\phi$ dependence of colloidal system [7]. In this case, the
fragility of the system can be externally tuned through the softness of the
particles [8, 9].

Therefore, the possibility that the classification along a ‘strong’ to
‘fragile’ scale could be applied not only to supercooled molecular liquids but
also to any glass forming systems should be explored more in detail. For this
reason, in this work we revisit our previous data on Laponite® aqueous
dispersions by performing a new data analysis to evidence how this system
behaves with respect to this universal glassy system description.

The investigated system is a synthetic clay that dispersed in water originates
a charged colloidal system characterized by a strong anisotropy due to the
shape of the nano-disc particles (with a 25:1 ratio) and due to the
inhomogeneity of the charges (negative and positive respectively distributed
on the faces and rims [10]) that determines a microscopic competition between
directional attractive and repulsive interactions. Laponite® peculiar phase
diagram has been largely investigated both experimentally [11, 12, 13, 14, 15,
16, 17, 10, 18] and theoretically [19, 20, 21, 22, 23, 24, 25, 26]. The system
spontaneously evolves towards different arrested states depending on the clay
concentration and ionic strength. In particular, in salt free water
conditions, an equilibrium gel [27] and a Wigner glass [28, 29] can be
recognized at weight clay concentrations $1.0\%\leqslant C_{\text{w}}<2.0\%$
and $2.0\%\leqslant C_{\text{w}}\leqslant 3.0\%$, respectively. Moreover, at
an increasing waiting time, a spontaneous glass-glass transition from a Wigner
to a Disconnected House of Cards (DHOC) glass stabilized by attractive
interactions has been found [30].

The aging evolution of dynamical properties of Laponite® has been broadly
investigated through different scattering techniques [31, 32, 33, 34, 35, 36,
30, 37]. In particular, Dynamic Light Scattering (DLS) has been widely
exploited to probe the evolution of fast and slow relaxations as a function of
waiting time ($t_{\text{w}}$) [38, 39, 40, 33, 14, 41, 42, 43] in aqueous
Laponite® suspensions. Moreover, DLS studies at different ionic salts [14,
17], adding polymers [44, 45, 46] and in different solvents such as D2O [47,
48, 49] have been performed to investigate their effects on the systems.

##  2 Materials and methods

###  2.1 Sample preparation

Laponite® is a synthetic layered silicate with chemical structure of the
unitary cell constituted by six octahedral magnesium ions sandwiched between
two layers of four tetrahedral silicon atoms groups, obeying the empirical
molecular formula
$\text{Na}^{+0.7}[(\text{Si}_{8}\text{Mg}_{5.5}\text{Li}_{0.3})\text{O}_{20}(\text{OH})_{4}]^{-0.7}$
[50]. The unitary cell is repeated around 1500 times in two dimensions to form
each Laponite® disc with a diameter of 25 nm and a thickness of 1 nm [10]. The
substitution of magnesium ions by lithium ions in the octahedral sheet
originates a negative charge of $-0.7e$ in a unit cell. These elementary
charges are uniformly distributed over the surface of Laponite® discs and,
therefore, the dissociation of OH- ions from the rims raises the pH of the
solution and leads to their positive charge [51]. Thus, the suspension is
composed by nanosized discs with inhomogeneous charge distribution: negative
on the surface, around several hundred unit charge $e$, and positive on the
rims, ten times lower for a salt-free system, as in our case [17, 31, 51].

The samples at different weight concentrations in the range 1.5% to 3.0% were
prepared, following the detailed protocol described in [10]: The oven-dried
Laponite® RD manufactured by Laporte Ltd was dispersed under stirring for 30
minutes in ultra pure deionized water ($C_{\text{s}}\approx 1\times 10^{-4}$
M). Soon after stirring, the system was filtered through a 0.45 µm pore size
Millipore filter. The origin of the waiting time ($t_{\text{w}}=0$) is the
time at which the suspension is filtered and sealed in glass tubes with 10 mm
of diameter for DLS measurements. All the procedure has been carried out in a
glovebox under N2 flux to prevent CO2 contamination.

###  2.2 DLS measurements

DLS measurements were performed with a five angles setup. A solid state laser
with wavelength $\lambda=642$ nm and power of 100 mW is focused in the center
of a cylindrical cuvette. Single mode fibers collect the scattered light at
five different scattering angles at around
$\theta=30^{\circ},50^{\circ},70^{\circ},90^{\circ}$ and $110^{\circ}$, which
correspond to five different scattering vectors, according to the relation
$Q=(4\piup n/\lambda)\sin(\theta/2)$, where $n$ is the refractive index of the
solvent. Time autocorrelation functions are then simultaneously computed at
different $Q$ by calculating the intensity autocorrelation function in the
time range $10^{-6}{\div}1$ s as:

$g_{2}(Q,t)=\frac{\left\langle I(Q,0)I(Q,t)\right\rangle_{t}}{\left\langle
I(Q,0)\right\rangle^{2}_{t}}\,,$ (2.1)

where $t$ is the delay time and $\left\langle\cdots\right\rangle_{t}$ denotes
the time average. This technique permits to explore both fast and slow
relaxations in Laponite® suspensions.

##  3 Results and discussions

Figure 1: (Colour online) Normalized intensity autocorrelation functions vs
delay time $t$, for Laponite® aqueous suspensions at weight concentration
$C_{\text{w}}=3.0$% at $T=298$ K, at increasing waiting time from
$t_{\text{w}}=11$ min (black symbols) to $t_{\text{w}}=325$ min (violet
symbols) and at a scattering vector $Q=1.8\times 10^{-3}$ Å-1 corresponding to
a scattering angle $\theta=90^{\circ}$. Full lines superimposed to symbols are
fits obtained through equation (3.1).

The slowing down of the dynamics induced by aging in Laponite® approaching the
gel and glass transition has been monitored by following the evolution of
autocorrelation functions with waiting time at different concentrations in the
range $C_{\text{w}}=(1.5{\div}3.0)$%. In figure 1, the normalized intensity
autocorrelation functions as obtained through DLS are reported, as an example,
for a high concentration sample at $C_{\text{w}}=3.0$% and $T=298$ K. By
increasing $t_{\text{w}}$, correlation curves slow down and become a more and
more stretched, typical signature of aging. For the longest investigated
waiting time ($t_{\text{w}}=325$ min), there is a qualitative change in the
correlation function with an evident crossover from a complete to an
incomplete decay to zero that indicates a strong ergodicity breaking due to
the formation of a glassy state. Figure 1 also shows that the correlation
functions present a two step decay. This behaviour has been attributed to the
presence of two different relaxation processes: a fast or microscopic
relaxation related to the interactions between an atom and the “cage” of its
nearest and a slow or structural relaxation associated with the structural
rearrangements of the particles. It has been observed in several systems
ranging from simple monoatomic liquids [52] to hydrogen bonds liquids [55, 53,
54], to glass-former systems [56].

Figure 2: (Colour online) Mean relaxation time $\tau_{\text{m}}$ defined
through equation (3.2) for aqueous Laponite® suspensions at different clay
concentrations in the glassy state $C_{\text{w}}=2.0$%, 2.5%, 2.8% and 3.0%,
and in the gel state $C_{\text{w}}=1.5$% and 1.8%, at $T=298$ K and scattering
vector $Q=1.8\times 10^{-3}$ Å-1. Solid lines represent the best fits through
equation (3.4).

In order to describe this two step decay, we used a fitting expression that is
the squared sum of an exponential and a stretched exponential function [40]:

$g_{2}(q,t)-1=b\left[a\mathrm{e}^{-t/\tau_{1}}+(1-a)\mathrm{e}^{-(t/\tau_{2})^{\beta}}\right]^{2},$
(3.1)

where $b$ represents the coherence factor, $\tau_{1}$ is the microscopic
relaxation time, $\tau_{2}$ is the structural relaxation time and $\beta$ is
the shape parameter that together with $\tau_{2}$ describes the slow part of
the autocorrelation function. The fits are shown as full lines in figure 1.

Figure 3: (Colour online) Angell plot for aqueous Laponite® suspensions in the
glassy state for $C_{\text{w}}=2.0$%, 2.5%, 2.8%, 3.0% and in the gel state
$C_{\text{w}}=1.5$% and 1.8%. Symbols and full line, obtained through the fits
with equation (3.4), represent the typical behaviour of fragile glass formers
while the dashed line represents the behaviour of strong systems. Figure 4:
(Colour online) Fragility index ($m$) obtained from the Angell plot as a
function of Laponite® suspensions concentration ranging from the gel to the
glassy state. Fragility values have been reported for Laponite® in salt free
water (squares) and deuterated water (circles), at two different salt
concentrations $C_{\text{S}}=10^{-3}$ M (up triangles) and
$C_{\text{S}}=2\times 10^{-3}$ M (down triangles) and in presence of Poly
ethylene oxide (PEO) at $C_{\text{PEO}}=0.37$% and $C_{\text{w}}=2.0$%
(exagons).

To characterize the slowing down of the dynamics towards the gel and glass
transition, we focus on the behaviour of the parameters related to structural
relaxation, namely $\tau_{2}$ and $\beta$ which permits to define the “mean”
relaxation time $\tau_{\text{m}}$ according to the relation:

$\tau_{\text{m}}=\tau_{2}\frac{1}{\beta}\Gamma\left(\frac{1}{\beta}\right),$
(3.2)

where $\Gamma$ is the usual Euler gamma function. The dependence of
$\tau_{\text{m}}$ on the waiting time is reported in figure 2 for aqueous
Laponite® suspension samples with increasing concentration. The sharp growth
observed clearly evokes the typical behaviour of supercooled liquids rapidly
quenched below their glass transition temperatures $T_{\text{g}}$. For these
systems, the structural relaxation time $\tau_{2}$ shows a Vogel-Fulcher-
Tammann (VFT) behaviour:

$\tau_{2}=\tau_{0}\mathrm{e}^{m/(T-T_{0})},$ (3.3)

where $T_{0}$ sets the apparent divergence of the structural relaxation and m
is the Fragility index. The evident analogy between colloidal and supercooled
systems suggests that the waiting time plays the same role of inverse
temperature. In this context, the fitting function for the “mean” relaxation
time $\tau_{\text{m}}$ can be written as:

$\tau_{\text{m}}=\tau_{0}\mathrm{e}^{mt_{\text{w}}/(t_{\text{w}\infty}-t_{\text{w}})},$
(3.4)

where $\tau_{0}$, $m$, and $t_{\text{w}\infty}$ are fitting parameters, in
particular, $t_{\text{w}\infty}$, as $T_{0}$ for supercooled liquids, sets the
divergence of $\tau_{\text{m}}$ associated with the formation of an arrested
state: at the beginning, the sample is liquid and then, with time, the
viscosity strongly increases until an arrested state is reached and the sample
does not flow anymore if turned upside down. Figure 2 clearly shows that the
mean relaxation time is strongly dependent on the clay concentration. In this
frame, similarly to a glass transition temperatures $T_{\text{g}}$ and
following the definition proposed by Angell for supercooled liquids, one can
define the glass transition time $t_{\text{g}}$ as the time at which
$\tau_{\text{m}}=100$ s for each Laponite® concentration. Properly normalizing
the relaxation times of figure 2, we obtained an Angell plot for Laponite®
shown in figure 3. Our data collapse into two master curves with two different
growths: one corresponding to samples that undergo a transition toward the
glassy state ($C_{\text{w}}=2.0\%{\div}3.0$%) and the other ones to samples
undergoing a transition toward the gel state ($C_{\text{w}}=1.5$% and 1.8%).
Both curves show the same behaviour expected for fragile supercooled liquids
while the dashed line corresponds to strong glass formers behaviour.

From the fits through equation (3.4) we obtained the fragility index $m$
reported in figure 4, it shows a step behaviour assuming values almost
constant below and above $C_{\text{w}}=2.0$%. This highlights that $m$ is not
dependent on the specific Laponite® concentration but only depends on the
sample final state, gel or glass. Our data are in good agreement with the
results reported in [43] in the high concentration range $C_{\text{w}}>2.0$%.
Moreover, in figure 4 fragility for aqueous Laponite® in salt free water
(squares) is compared with fragility values found in deuterated water
(circles) from [48], in water at two different salt concentrations
$C_{\text{S}}=10^{-3}$ M (up triangles) and $C_{\text{S}}=2\times 10^{-3}$ M
(down triangles) from [14] and in presence of Poly Ethylene Oxide (PEO) at
$C_{\text{PEO}}=0.37$% in the case of a sample at Laponite® concentration
$C_{\text{w}}=2.0$% (exagons) from [46]. One can observe that, despite the
change of solvent, the addition of salt or the presence of a polymer like PEO
that can fasten or slow down the aging dynamics, what really matters to
identify the fragility of the system is the final state reached by the sample.

##  4 Conclusions

In this work we revisited our previous data on Laponite® aqueous dispersions
by performing a new analysis to evidence analogies with supercooled liquids
approaching their glass transitions. In particular, the slowing down of the
dynamics induced by aging was monitored by following the temporal evolution of
autocorrelation functions measured through DLS. Colloidal suspensions of
Laponite® at different concentrations in the range
$C_{\text{w}}=(1.5{\div}3.0)$% were investigated and the structural relaxation
times reported as a function of the waiting time. Making a parallelism with
the glass transitions in supercooled liquids, a glass transition time
$t_{\text{g}}$ was defined and an Angell plot drawn. Finally, if changing
concentration, solvent, salt content or adding a polymer could fasten or slow
down the aging dynamics of Laponite® suspensions, what really matters to
identify the fragility of the system is the final state reached by the sample.

## References

  * [1] Debenedetti P.G., Stillinger F.H., Nature, 2001, 410, 259, doi:10.1038/35065704.
  * [2] Zaccarelli E., J. Phys.: Condens. Matter, 2007, 19, 323101, doi:10.1088/0953-8984/19/32/323101.
  * [3] Angell C., Proc. Natl. Acad. Sci. U.S.A., 1995, 92, 6675–6682, doi:10.1073/pnas.92.15.6675.
  * [4] Ediger M., Angell C., Nagel S.R., J. Phys. Chem., 1996, 100, 13200–13212, doi:10.1021/jp953538d.
  * [5] Angell C., Ngai K.L., McKenna G.B., McMillan P.F., Martin S.W., J. Appl. Phys., 2000, 88, 3113–3157, doi:10.1063/1.1286035.
  * [6] Böhmer R., Ngai K.L., Angell C., Plazek D., J. Chem. Phys., 1993, 99, No. 5, 4201, doi:10.1063/1.466117.
  * [7] Mattsson J., Wyss H.M., Fernandez-Nieves A., Miyazaki K., Hu Z., Reichman D., Weitz D.A., Nature, 2009, 462, No. 5, 83–86, doi:10.1038/nature08457.
  * [8] Nigro V., Angelini R., Bertoldo M., Bruni F., Ricci M., Ruzicka B., Soft Matter, 2017, 13, 5185–5193, doi:10.1039/C7SM00739F.
  * [9] Nigro V., Angelini R., Rosi B., Bertoldo M., Buratti E., Casciardi S., Sennato S., Ruzicka B., J. Colloid Interface Sci., 2019, 545, 210–219, doi:10.1016/j.jcis.2019.03.004.
  * [10] Ruzicka B., Zaccarelli E., Soft Matter, 2011, 7, 1268, doi:10.1039/c0sm00590h.
  * [11] Mourchid A., Lécolier E., Van Damme H., Levitz P., Langmuir, 1998, 14, 4718–4723, doi:10.1021/la980117p.
  * [12] Tanaka H., Meunier J., Bonn D., Phys. Rev. E, 2004, 69, 031404, doi:10.1103/PhysRevE.69.031404.
  * [13] Mongondry P., Tassin J.F., Nicolai T., J. Colloid Interface Sci., 2005, 283, 397, doi:10.1016/j.jcis.2004.09.043.
  * [14] Ruzicka B., Zulian L., Ruocco G., Langmuir, 2006, 22, 1106, doi:10.1021/la0524418.
  * [15] Cummins H.Z., J. Non-Cryst. Solids, 2007, 353, 3891, doi:10.1016/j.jnoncrysol.2007.02.066.
  * [16] Ruzicka B., Zulian L., Ruocco G., Philos. Mag., 2007, 87, 449, doi:10.1080/14786430600962732.
  * [17] Jabbari-Farouji S., Tanaka H., Wegdam G.H., Bonn D., Phys. Rev. E, 2008, 78, 061405,
doi:10.1103/PhysRevE.78.061405.

  * [18] Suman K., Joshi Y.M., Langmuir, 2018, 34, 13079–13103, doi:10.1021/acs.langmuir.8b01830.
  * [19] Dijkstra M., Hansen J.P., Madden P.A., Phys. Rev. Lett., 1995, 75, 2236–2239,
doi:10.1103/PhysRevLett.75.2236.

  * [20] Kutter S., Hansen J., Sprik M., Boek E., J. Chem. Phys., 2000, 112, 311–322, doi:10.1063/1.480582.
  * [21] Trizac E., Bocquet L., Agra R., Weis J., Aubouy M., J. Phys.: Condens. Matter, 2002, 14, 9339–9352, doi:10.1088/0953-8984/14/40/322.
  * [22] Odriozola G., Romero-Bastida M., de J. Guevara-Rodríguez F., Phys. Rev. E, 2004, 70, No. 2, 021405, doi:10.1103/PhysRevE.70.021405.
  * [23] Mossa S., De Michele C., Sciortino F., J. Chem. Phys., 2007, 126, 014905, doi:10.1063/1.2408418.
  * [24] Jönsson B., Labbez C., Cabane B., Langmuir, 2008, 24, 11406–11413, doi:10.1021/la801118v.
  * [25] Delhorme M., Jönsson B., Labbez C., Soft Matter, 2012, 8, 9691, doi:10.1039/c2sm25731a.
  * [26] Jabbari-Farouji S., Weis J.J., Davidson P., Levitz P., Trizac E., Sci. Rep., 2013, 3, 3559–3565,
doi:10.1038/srep03559.

  * [27] Ruzicka B., Zaccarelli E., Zulian L., Angelini R., Sztucki M., Moussaïd A., Narayanan T., Sciortino F., Nat. Mater., 2011, 10, 56–60, doi:10.1038/nmat2921.
  * [28] Bonn D., Tanaka H., Wegdam G., Kellay H., Meunier J., Europhys. Lett., 1999, 45, 52–57,
doi:10.1209/epl/i1999-00130-3.

  * [29] Ruzicka B., Zulian L., Zaccarelli E., Angelini R., Sztucki M., Moussaïd A., Ruocco G., Phys. Rev. Lett., 2010, 104, 085701, doi:10.1103/PhysRevLett.104.085701.
  * [30] Angelini R., Zaccarelli E., de Melo Marques F.A., Sztucki M., Fluerasu A., Ruocco G., Ruzicka B., Nat. Commun., 2014, 5, 4049, doi:10.1038/ncomms5049.
  * [31] Martin C., Pignon F., Piau J.M., Magnin A., Lindner P., Cabane B., Phys. Rev. E, 2002, 66, 021401, doi:10.1103/PhysRevE.66.021401.
  * [32] Bandyopadhyay R., Liang D., Yardimci H., Sessoms D.A., Borthwick M.A., Mochrie S.G.J., Harden J.L., Leheny R.L., Phys. Rev. Lett., 2004, 93, 228302, doi:10.1103/PhysRevLett.93.228302.
  * [33] Schosseler F., Kaloun S., Skouri M., Munch J.P., Phys. Rev. E, 2006, 73, 021401,
doi:10.1103/PhysRevE.73.021401.

  * [34] Ruzicka B., Zulian L., Angelini R., Sztucki M., Moussaïd A., Ruocco G., Phys. Rev. E, 2008, 77, 020402, doi:10.1103/PhysRevE.77.020402.
  * [35] Angelini R., Zulian L., Fluerasu A., Madsen A., Ruocco G., Ruzicka B., Soft Matter, 2013, 9, 10955, doi:10.1039/c3sm52173g.
  * [36] Hansen E.L., Jabbari-Farouji S., Mauroy H., Plivelic T.S., Bonn D., Fossum J.O., Soft Matter, 2013, 9, 9999, doi:10.1039/c3sm52103f.
  * [37] De Melo Marques F.A., Angelini R., Zaccarelli E., Farago B., Ruta B., Ruocco G., Ruzicka B., Soft Matter, 2015, 11, 466–471, doi:10.1039/C4SM02010C.
  * [38] Nicolai T., Cocard S., J. Colloid Interface Sci., 2001, 244, 51–57, doi:10.1006/jcis.2001.7930.
  * [39] Bellour M., Knaebel A., Harden J.L., Lequeux F., Munch J.P., Phys. Rev. E, 2003, 67, 031405,
doi:10.1103/PhysRevE.67.031405.

  * [40] Ruzicka B., Zulian L., Ruocco G., Phys. Rev. Lett., 2004, 93, 258301, doi:10.1103/PhysRevLett.93.258301.
  * [41] Jabbari-Farouji S., Wegdam G.H., Bonn D., Phys. Rev. Lett., 2007, 99, 065701,
doi:10.1103/PhysRevLett.99.065701.

  * [42] Pujala R.K., Bohidar H.B., Soft Matter, 2012, 8, 6120, doi:10.1039/c2sm25337b.
  * [43] Saha D., Joshi Y.M., Bandyopadhyay R., Soft Matter, 2014, 10, 3292–3300, doi:10.1039/c4sm00187g.
  * [44] Mongondry P., Nicolai T., Tassin J.F., J. Colloid Interface Sci., 2004, 275, 191–196,
doi:10.1016/j.jcis.2004.01.037.

  * [45] Zulian L., Ruzicka B., Ruocco G., Philos. Mag., 2008, 88, 4213–4221, doi:10.1080/14786430802581876.
  * [46] Zulian L., de Melo Marques F.A., Emilitri E., Ruocco G., Ruzicka B., Soft Matter, 2014, 10, 4513–4521, doi:10.1039/c4sm00172a.
  * [47] Tudisca V., Ricci M., Angelini R., Ruzicka B., RSC Adv., 2012, 2, 11111–11116, doi:10.1039/c2ra21486e.
  * [48] De Melo Marques F.A., Angelini R., Ruocco G., Ruzicka B., J. Phys. Chem. B, 2017, 121, 4576–4582, doi:10.1021/acs.jpcb.6b12596.
  * [49] Tudisca V., Bruni F., Scoppola E., Angelini R., Ruzicka B., Zulian L., Soper A.K., Ricci M.A., Phys. Rev. E, 2014, 90, No. 3, 032301, doi:10.1103/PhysRevE.90.032301.
  * [50] Lal J., Auvray L., J. Appl. Crystallogr., 2000, 33, 673–676, doi:10.1107/S0021889899013308.
  * [51] Tawari S.L., Koch D.L., Cohen C., J. Colloid Interface Sci., 2001, 240, 54–66, doi:10.1006/jcis.2001.7646.
  * [52] Balucani U., Zoppi M., Dynamics of the Liquid State, Clarendon Press, Oxford, 1994\.
  * [53] Giura P., Angelini R., Ruocco G., Sette F., J. Chem. Phys., 2007, 127, 084508, doi:10.1063/1.2753161.
  * [54] Angelini R., Giura P., Fioretto D., Monaco G., Ruocco G., Sette F., Phys. Rev. B, 2004, 70, 224302, doi:10.1103/PhysRevB.70.224302.
  * [55] Angelini R., Giura P., Monaco G., Ruocco G., Sette F., Verbeni R., Phys. Rev. Lett., 2002, 88, 255503, doi:10.1103/PhysRevLett.88.255503.
  * [56] Götze W., Mayr M., Phys. Rev. E, 2000, 61, 587, doi:10.1103/PhysRevE.61.587.

Ukrainian 0

Ïåðåõä â ñòàí ãåëþ òà ñêëà â êðèõêèõ êîëî¿äíèõ ñóñïåíçÿõ Ð.
Àíäæåëí${}^{\textsf{\footnotesize{\mbox{?},\penalty 1000\mbox{?}}}}$, Äæ.
Ðóîêêî${}^{\textsf{\footnotesize{\mbox{?},\penalty 1000\mbox{?}}}}$, Á.
Ðóçöêà${}^{\textsf{\footnotesize{\mbox{?},\penalty 1000\mbox{?}}}}$

${}^{\textsf{\footnotesize 1}}$ ISC-CNR, I-00185 Ðèì, òàëÿ
${}^{\textsf{\footnotesize 2}}$ Ôçèíèé ôàêóëüòåò, Óíâåðñèòåò Ðèìó
‘‘Sapienza’’, I-00185, òàëÿ ${}^{\textsf{\footnotesize 3}}$ Öåíòð íàíîíàóê ïðî
æèòòÿ, IIT@Sapienza, òàëéñüêèé íñòèòóò òåõíîëîãé, 00161, Ðèì, òàëÿ


# Axisymmetric pseudoplastic thin films in planar and spherical geometries

Chris Reese
Department of Physics
Lewis and Clark Community College
Godfrey, IL 62035
<EMAIL_ADDRESS>
###### Abstract

A simplified, pseudoplastic rheology characterized by constant viscosity
plateaus above and below a transition strain rate is applied to axisymmetric,
gravitationally driven spreading of a thin fluid film with constant volume
flux source in planar and spherical geometries. The model admits analytical
solutions for flow velocity and volume flux. Shear thinning influence on layer
evolution is investigated via numerical simulation. Isoviscous, asymptotic
behaviors are recovered in small and large transition stress limits. The
effect of viscosity ratio on layer extent agrees with scaling arguments. For
intermediate transition stress, a flow behavior adjustment is observed
consistent with heuristic arguments. Planar and spherical geometry solutions
are in agreement for sufficiently small polar angle.

##  1 Introduction

Gravitational spreading of one fluid into another of differing density, where
flow is predominantly unidirectional (e.g., horizontal or down-slope)
constitutes a gravity current. Thin film, viscous gravity currents have been a
subject of interest since the seminal work of Huppert [1982a, b] resulting in
a significant analytical, numerical, and experimental literature addressing
many aspects of layer dynamics, see e.g., [Simpson, 1999; Huppert, 2006;
Ancey, 2007; Ungarish, 2009].

While first developed for isoviscous fluids, thin film theory has been
extended to include many rheological behaviors. For example, temperature
dependent viscosity is of particular importance in geophysical contexts
[Stasuik et al., 1993; Bercovici, 1994; Balmforth et al., 2004; Algwauish,
2019]. Also, strain rate dependent viscosities (i.e., generalized newtonian
fluids) are ubiquitous in natural and engineering applications. Viscoplastic
rheologies such as Bingham and Herschel-Bulkley models find application to
lava dome evolution [Blake, 1990; Balmforth and Craster, 1999; Balmforth et
al., 2000, 2006]. Another behavior of practical interest is pseudoplasticity.
A pseudoplastic, or shear thinning, fluid is characterized by an effective
viscosity which decreases with increasing strain rate. This behavior is
observed in polymeric [Bird et al., 1987] and geophysical [Lavallée et al.,
2007; Vasseur et al., 2023] fluids. Shear thinning, power law constitutive
models admit various analytical approaches in the thin layer limit [Gratton et
al., 1999; Perazzo and Gratton, 2003; Myers, 2005; Neossi Nguetchue and
Momoniat, 2007]. Other shear thinning rheologies include Sisko, Cross, and
Carreau models. Often, approximations to these models are adopted (e.g.,
[Wrobel, 2020; James et al., 2021]) in the interest of analytical simplicity
(however, cf. Pritchard et al. [2015]).

Geometries addressed in the literature are primarily two-dimensional and
axisymmetric flows on a plane. However, thin film theory has also been
extended to curved substrates [Oron et al., 1997; Roy et al., 2002; Kang et
al., 2016, 2017; Taranets, 2019]. Takagi and Huppert [2010] considered flow
and fingering instability of thin films on a cylinder and sphere with vertical
gravitational acceleration. Fluid spreading on a sphere with radial
gravitational acceleration was considered in a geodynamical context [Reese et
al., 2010, 2011]. This geometry may find application to planetary scale (i.e.,
low spherical harmonic degree) superplume head evolution (e.g., Bercovici and
Lin [1996]; Kerr and Mériaux [2004]) and/or global scale thermochemical diapir
dynamics [Watters et al., 2009] on terrestrial planets and icy satellites.

In this work, a simplified pseudoplastic rheology [James et al., 2021] is
adopted to investigate axisymmetric gravity currents in planar and spherical
geometries. The generalized newtonian viscosity (Sec 3.1) is a three
parameter, piecewise linear function characterized by small and large strain
rate viscosity plateaus and a rheological transition stress. Section 2 is a
heuristic, analysis of flow dynamics. The model is reviewed in section 3.
Section 4 presents numerical results addressing effects of pseudoplastic
transition stress and viscosity ratio variation. A summary and concluding
comments are provided in section 5. The appendix (section 6) outlines
benchmarks of the numerical methods employed in the study.

##  2 Scaling analysis

Vertical and horizontal length scalings can be derived from dimensional
analysis (e.g., Griffiths and Fink [1993]) based on the fundamental, thin film
caveat that layer thickness $H$ is small with respect to lateral extent $L$.
To first order in the aspect ratio $H/L$, scaling analysis implies: flow is
primarily tangential to the surface, normal stress is negligible resulting in
a hydrostatic pressure distribution, and the tangential pressure gradient is
balanced by the shear stress gradient perpendicular to the surface.

###  2.1 Small strain rate

In the isoviscous, small strain rate limit,
$\dot{\varepsilon}\ll\dot{\varepsilon_{c}}$, spreading is controlled by
viscosity $\eta$ (see Sec. 3.1). For lateral velocity scale $V$, continuity
implies that the vertical velocity is first order small in the aspect ratio
$\dfrac{H}{L}\,V$. The pressure scale $p\sim\rho gH$. The pressure gradient is
balanced by vertical shear stress gradient,

$\frac{\rho gH}{L}\sim\frac{\tau}{H}\;,$ (2.1)

where the shear stress scale $\tau\sim\eta\dfrac{V}{H}$ implying

$\frac{\rho gH}{L}\sim\eta\frac{V}{H^{2}}\;.$ (2.2)

Kinematics require that the lateral velocity scale

$V\sim L/t\;,$ (2.3)

where $t$ is time. For a constant volume flux $Q$, mass conservation requires

$Qt\sim L^{2}H\;.$ (2.4)

Eliminating $H$ and $V$ between Eqs. (2.2,2.3,2.4) yields the lateral length
scaling as a function of time

$L\sim\left(\frac{Q^{3}\rho g}{\eta}\right)^{1/8}\;t^{1/2}\,.$ (2.5)

It follows that the layer height scales as

$H\sim\left(\frac{\eta\,Q}{\rho g}\right)^{1/4}\;.$ (2.6)

###  2.2 Large strain rate

For large strain rate, $\varepsilon\gg\varepsilon_{c}$, spreading is
controlled by the high strain rate viscosity $\mu$ (see Sec. 3.1). In this
limit, the lateral and height scales are expected to be

$L_{\mu}\sim\chi^{1/8}\;L\,,\qquad\qquad H_{\mu}\sim\chi^{-1/4}\;H\,,$ (2.7)

where the viscosity ratio $\chi=\eta/\mu$. For intermediate strain rate, i.e.
$\tau\sim\tau_{c}$, there is no asymptotic scaling for $L$ and $H$ throughout
layer evolution as illustrated by numerical results described below.

##  3 Model

Consider a fluid of density $\rho$ and strain rate dependent generalized
newtonian viscosity $\eta_{\rm eff}(\dot{\varepsilon})$ spreading on a rigid
surface. For characteristic velocity and length scales $V$ and $L$
respectively, a Reynolds number can be defined as the ratio of the momentum
diffusion time to advection time,

${\rm Re}=\frac{L^{2}/\nu}{L/V}=\frac{\rho VL}{\eta_{\rm eff}}\;,$ (3.1)

where $\nu=\eta_{\rm eff}/\rho$ is the kinematic viscosity. Sufficiently small
Re guarantees non-inertial flow reducing the Cauchy momentum equation

$-\nabla p+\nabla\cdot\boldsymbol{\tau}+\rho\,\boldsymbol{g}=0\;,$ (3.2)

where $p$ is pressure and $\boldsymbol{g}$ is gravitational acceleration.

For an incompressible fluid, the continuity equation

$\nabla\cdot\boldsymbol{u}=0\;.$ (3.3)

Free surface boundary conditions are zero pressure $p=0$, zero traction
$\boldsymbol{\tau}\cdot\hat{\boldsymbol{n}}=0$, and the kinematic condition
specifying that a fluid parcel on the boundary remain on the boundary. Basal
boundary conditions are no-slip (i.e., zero velocity component tangential to
the basal surface), and matching of fluid velocity perpendicular to the layer
base with source velocity $w_{s}$.

###  3.1 Pseudoplastic constitutive relation

For a generalized newtonian rheology, the effective viscosity is a function of
strain rate. The deviatoric stress tensor

$\boldsymbol{\tau}=2\,\eta_{\rm
eff}(\dot{\varepsilon})\,\boldsymbol{\dot{\varepsilon}}\;,$ (3.4)

where

$\boldsymbol{\dot{\varepsilon}}=\frac{1}{2}\left[\nabla\boldsymbol{u}+\nabla\boldsymbol{u}^{\rm
T}\right]\;,$ (3.5)

is the strain rate tensor and

$\dot{\varepsilon}=\left[\frac{1}{2}\;\operatorname{Tr}\left(\boldsymbol{\dot{\varepsilon}}\boldsymbol{\dot{\varepsilon}}^{\rm
T}\right)\right]^{1/2}\,,\qquad\tau=\left[\frac{1}{2}\;\operatorname{Tr}\left(\boldsymbol{\tau}\boldsymbol{\tau}^{\rm
T}\right)\right]^{1/2}\;,$ (3.6)

are the strain rate and stress invariants, respectively.

The piecewise linear approximation to a pseudoplastic rheology adopted in this
study is

$\eta_{\rm
eff}(\dot{\varepsilon})=\begin{cases}\eta&\dot{\varepsilon}\leq\dot{\varepsilon}_{\rm
c}\\\
\dfrac{\tau_{c}}{2\dot{\varepsilon}}\left(1-\dfrac{\mu}{\eta}\right)+\mu&\dot{\varepsilon}>\dot{\varepsilon}_{\rm
c}\\\ \end{cases}$ (3.7)

where $\dot{\varepsilon}_{\rm c}$ is the critical strain rate for onset of
shear thinning, $\tau_{c}=2\,\eta\,\dot{\varepsilon}_{\rm c}$ is the critical
stress invariant, $\eta$ is the low strain rate viscosity, and $\mu$ is the
high strain rate viscosity. This rheological model is compared to the
isoviscous case and a Cross fluid in Fig. 1.

Figure 1: A qualitative comparison of various rheologies. (left) The stress
invariant versus strain invariant for isoviscous (blue), simplified
pseudoplastic (red), and Cross fluids (green). All fluids have the same small
strain rate viscosity. The simplified pseduoplastic and Cross fluids approach
the same asymptotic large strain rate viscosity. The parameter $\chi$ is the
ratio of small strain rate to large strain rate viscosities, $\chi=\eta/\mu$.
Solid and dashed curves correspond to $\xi$ = 3 and 10, respectively. (right)
Effective viscosities of the rheological models.

###  3.2 Planar geometry

In this section, the governing equations are non-dimensionalized and expanded,
to leading order, in the thin layer limit for the rheological model described
above. This constitutive relation admits determination of the radial velocity
profile and radial volume flux. Vertical integration of the continuity
equation yields the layer height evolution equation. For cylindrical
coordinates ($r$, $z$, $\phi$), the axisymmetric velocity

$\boldsymbol{u}=u(r,z,t)\,\hat{\mathbf{r}}+w(r,z,t)\,\hat{\mathbf{z}}\;.$
(3.8)

The vertical locations of the flow base and free surface are $z=0$ and
$z=h(r,t)$, respectively. The gravitational acceleration
$\boldsymbol{g}=-\,g\hat{\boldsymbol{z}}$.

####  3.2.1 Non-dimensionalization

Motivated by heuristic arguments (Sec. 2), the governing equations are non-
dimensionalized and expanded in the small layer aspect ratio limit. For a
rigorous formulation of the expansion in planar axisymmertry the reader is
referred to Balmforth et al. [2000]. As in Section 2, let $H$ be the
characteristic layer thickness and $L$ be the radial length scale. The
vertical coordinate $z$ and layer thickness $h$ are scaled by $H$ while
pressure is scaled by $\rho gH$. The radial coordinate $r$ is scaled by $L$.
Continuity implies that the vertical velocity scale is first order small with
respect to radial velocity. To leading order, the radial velocity scale
$V=\rho gH^{3}/\eta L$ representing a balance between lateral pressure
gradient and vertical shear stress gradient. Time is scaled by $L/V$. In the
thin layer limit, the aspect ratio $H/L$ is the small parameter in which the
governing equations are expanded.

####  3.2.2 Pseudoplastic transition height

To leading order, normal stress is negligible and gravitational body force is
balanced by hydrostatic pressure subject to the free surface boundary
condition, $\left.p\right|_{z=h}=0$,

$p=h\,\left(1-\frac{z}{h}\right)\;.$ (3.9)

The radial pressure gradient due to layer height variation is balanced by
vertical shear of the radial flow

$\frac{\partial p}{\partial r}=\frac{\partial\tau_{rz}}{\partial z}\;.$ (3.10)

Substituting for $p$, integrating, and applying the boundary condition
$\left.\tau_{rz}\right|_{z=h}=0$ yields

$\tau_{rz}=-hh^{\prime}\,\left(1-\frac{z}{h}\right)\;,$ (3.11)

where $h^{\prime}=\dfrac{\partial h}{\partial r}$. The stress invariant
$\tau=h|h^{\prime}|\,\left(1-\dfrac{z}{h}\right)$ increases as $z$ decreases
from $z=h$ and there is a vertical level $Y$ where $\tau=\tau_{c}$,

$Y=h\left[1-\frac{\tau_{c}}{h|h^{\prime}|}\right]=h\left[1-\frac{2\,\dot{\varepsilon}_{c}}{h|h^{\prime}|}\right]\;.$
(3.12)

####  3.2.3 Radial velocity

The constitutive equation implies

$\tau_{rz}=2\,\eta_{\rm
eff}(\dot{\varepsilon})\,\dot{\varepsilon}_{rz}=\eta_{\rm
eff}\,(\dot{\varepsilon})\dfrac{\partial u}{\partial z}\;,$ (3.13)

where the strain rate invariant
$\dot{\varepsilon}=\dfrac{1}{2}\,\left|\dfrac{\partial u}{\partial z}\right|$.
Let $u_{-}(z)$ be the radial velocity for $z<Y$. Substituting for $\eta_{\rm
eff}$, with radial velocity increasing monotically with height
$\left(\dfrac{\partial u_{-}}{\partial z}>0\right)$ yields

$\frac{\partial u_{-}}{\partial z}=-\chi
hh^{\prime}\left(1-\frac{z}{h}\right)-2\,(\chi-1)\,\dot{\varepsilon}_{c}\;,$
(3.14)

where $\chi=\dfrac{\eta}{\mu}$ is the low strain rate to high strain rate
viscosity ratio. Integrating with respect to $z$, and applying the no slip
boundary condition $\left.u_{-}\right|_{z=0}=0$,

$u_{-}(z)=-\chi
h^{2}h^{\prime}f\left(\frac{z}{h}\right)-2\,(\chi-1)\,h\,\dot{\varepsilon}_{c}\frac{z}{h}\;,\qquad
0\leq z\leq Y\;,$ (3.15)

with $f(x)=x-x^{2}/2$.

Above $z=Y$, $\eta_{\rm eff}=\eta$, and the momentum balance reduces to the
isoviscous case. Letting $u_{+}(z)$ be the velocity for $z>Y$,

$\frac{\partial u_{+}}{\partial z}=-hh^{\prime}\left(1-\frac{z}{h}\right)\;.$
(3.16)

Integrating and matching velocities $u_{-}(Y)=u_{+}(Y)$ gives,

$u_{+}(z)=-h^{2}h^{\prime}f\left(\frac{z}{h}\right)-(\chi-1)\,h^{2}h^{\prime}f\left(\frac{Y}{h}\right)-2\,(\chi-1)\,h\,\dot{\varepsilon}_{c}\frac{Y}{h}\;,\qquad
Y\leq z\leq h\;.$ (3.17)

Eliminating the transition strain rate using Eq. (3.12),

$u(z)=\begin{cases}-h^{2}h^{\prime}\left[\chi
f\left(\dfrac{z}{h}\right)-(\chi-1)\left(1-\dfrac{Y}{h}\right)\dfrac{z}{h}\right]&0\leq
z\leq Y\\\\[10.00002pt]
-h^{2}h^{\prime}\left[f\left(\dfrac{z}{h}\right)+(\chi-1)\dfrac{Y^{2}}{2h^{2}}\right]&Y\leq
z\leq h\end{cases}$ (3.18)

This velocity distribution is shown in Fig. 2. In the limit $Y\rightarrow 0$,
the velocity $u(z)=-h^{2}h^{\prime}f(z/h)$ which corresponds to the $\eta_{\rm
eff}\rightarrow\eta$ isoviscous limit. Likewise, when $Y\rightarrow h$,
$u(z)=-\chi h^{2}h^{\prime}f(z/h)$ which is the $\eta_{\rm eff}\rightarrow\mu$
isoviscous limit.

Figure 2: Radial velocity profiles for the cases $\chi=3$ (left), $\chi=10$
(right), and different pseduoplastic transition levels (left to right) $Y/h$ =
(0, 0.3, 0.5, 0.9, 1). The black curves correspond to the low strain rate and
high strain rate isoviscous limits (see discussion in text).

The height evolution equation (see next section) depends on the vertically
integrated radial velocity, i.e., the radial volume flux per unit azimuthal
length

$U=\int_{0}^{Y}u_{-}(z)\,dz+\int_{Y}^{h}u_{+}(z)\,dz\;.$ (3.19)

Evaluating,

$U=-h^{3}h^{\prime}\left[\frac{1}{3}+\frac{1}{2}\,(\chi-1)\left(\frac{Y}{h}\right)^{2}\left(1-\frac{Y}{3h}\right)\right]\;.$
(3.20)

In the low strain rate $Y\rightarrow 0$ limit, $U$ reduces to the isoviscous
$\eta_{\rm eff}\rightarrow\eta$ case, i.e., $U=-h^{3}h^{\prime}/3$. Likewise
for $Y\rightarrow h$, $U=-\chi h^{3}h^{\prime}/3$ which is the $\eta_{\rm
eff}\rightarrow\mu$ isoviscous limit.

####  3.2.4 Evolution equation

The free surface kinematic condition

$\frac{\partial h}{\partial t}+u(r,h,t)\frac{\partial h}{\partial
r}-w(r,h,t)=0\;.$ (3.21)

Vertically integrating the continuity equation over the layer height subject
to basal and free surface boundary conditions and identifying
$w(r,0,t)=w_{s}(r,t)$ yields the layer height evolution equation representing
mass conservation,

$\frac{\partial h}{\partial t}+\frac{1}{r}\frac{\partial}{\partial
r}\left(rU\right)=w_{s}(r,t)\;,$ (3.22)

where $U$ is given by Eq. (3.20).

###  3.3 Spherical geometry

The following section is a brief outline the thin layer expansion in spherical
geometry. It is shown that, to leading order, local flow on a sufficiently
large sphere is insensitive to substrate curvature resulting in polar velocity
and volume flux identical to the planar case. In spherical coordinates
($\theta$, $r$, $\phi$), the axisymmetric velocity field

$\boldsymbol{u}=u(\theta,r,t)\,\hat{\boldsymbol{\theta}}+w(\theta,r,t)\,\hat{\boldsymbol{r}}\;.$
(3.23)

The radial locations of the flow base and free surface are $r=R$ and
$r=R+h(\theta,t)$, respectively. The gravitational acceleration
$\boldsymbol{g}=-\,g\,\hat{\boldsymbol{r}}$. The layer is considered to extend
laterally along an arclength $R\,\theta_{f}(t)$. The thin layer approximation
$h\ll R\,\theta_{f}$ is assumed to hold throughout spreading given
sufficiently large substrate curvature $R$. Because $\theta_{f}\sim 1$, it
follows that $h\ll R$.

The radius of curvature $R$ provides an intrinsic lengthscale for non-
dimensionlization. Lengths are scaled by $R$ and pressure by $\rho gR$. The
polar velocity scale $V=\rho gR^{2}/\eta$ and timescale is $R/V$. Upon non-
dimensionalizing, the quantity in which the governing equations are expanded
is $h\ll 1$. A change of variables is introduced

$r=1+\xi\;,$ (3.24)

where $0\leq\xi\leq h$ is the non-dimensional radial coordinate measured from
the spherical surface.

The radial pressure gradient balances the gravitational body force subject to
the free surface boundary condition $\left.p\right|_{\xi=h}=0$,

$p=h\,\left(1-\frac{\xi}{h}\right)\;.$ (3.25)

The polar pressure gradient due to layer height variation is balanced by
radial shear of the polar flow. In terms of $\xi$,

$\frac{1}{(1+\xi)}\frac{\partial
p}{\partial\theta}=\frac{1}{(1+\xi)^{2}}\frac{\partial}{\partial\xi}\left((1+\xi)^{2}\tau_{\xi\theta}\right))\;.$
(3.26)

Substituting for $p$, applying the boundary condition
$\left.\tau_{\xi\theta}\right|_{\xi=0}=0$, and dropping terms
${\cal{O}}(h^{2})$ and higher,

$\tau_{\xi\theta}=-hh^{\prime}\left(1-\frac{\xi}{h}\right)\;,$ (3.27)

where $h^{\prime}=\dfrac{\partial h}{\partial\theta}$. This linear stress
distribution is identical to the planar case Eq. (3.11). Thus, the radial
level $\xi=Y$ below which the stress invariant exceeds $\tau_{c}$ is given by
Eq. (3.12).

To leading order, the constitutive equation

$\tau_{\xi\theta}=2\,\eta_{\rm
eff}(\dot{\varepsilon})\,\dot{\varepsilon}_{\xi\theta}=\eta_{\rm
eff}\,(\dot{\varepsilon})\dfrac{\partial u}{\partial\theta}\;,$ (3.28)

where the strain rate invariant
$\dot{\varepsilon}=\dfrac{1}{2}\,\left|\dfrac{\partial
u}{\partial\theta}\right|$. Analysis proceeds as in the planar case. The polar
velocities $u_{\pm}(\xi)$ above and below $Y$ satisfy Eqs. (3.16) and (3.14),
respectively. The polar velocity profile $u(\xi)$ and polar volume flux per
unit azimuthal length $U$ are given by Eqs. (3.18) and (3.20), respectively,
with $h^{\prime}\rightarrow\dfrac{\partial h}{\partial\theta}$ and
$z\rightarrow\xi$.

To leading order, in terms of $(\theta,\xi,t)$, the free surface kinematic
condition,

$\frac{\partial h}{\partial t}+u(\theta,h,t)\frac{\partial
h}{\partial\theta}-w(\theta,h,t)=0\;.$ (3.29)

Radially integrating the continuity equation subject to basal and free surface
boundary conditions, identifying $w(\theta,0,t)=w_{s}(\theta,t)$ and dropping
terms $\mathcal{O}(h^{2})$ and higher yields

$\frac{\partial h}{\partial
t}+\left[\frac{1}{\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\,U\right)\right]=w_{s}(\theta,t)\;,$
(3.30)

where $U$ is given by Eq. (3.20)

##  4 Numerical results

###  4.1 Planar geometry

Variation of rheological transition stress and viscosity ratio are
investigated numerically. The transition stress range is chosen to include
anticipated end-member behaviors. The non-dimensional source function

$w_{s}(r)=w_{0}\,(r_{0}^{2}-r^{2})\,H(r_{0}-r)\;,$ (4.1)

where $r_{0}$ = 0.15, $w_{0}$ = 0.1, and $H$ is the unit step function.

####  4.1.1 Height field

The effect of varying the pseudoplastic transition stress is considered for
viscosity ratio $\chi$ = 10. For large transition stress ($\tau_{c}=300\times
10^{-4}$), the location of the rheological transition surface $Y$ is
indistinguishable from zero except for small regions near the flow front (Fig.
3). Also, the solution converges to the similarity solution in the similarity
variable $\zeta=r/t^{1/2}$. away from the source. Thus, the fluid behaves
isoviscously with $\eta_{\rm eff}\approx\eta$.

Figure 3: Evolution for $\tau_{c}=300\times 10^{-4}$ and $\chi=10$. (top)
Layer height field $h(r,t)$ (blue) together with the transition surface
$Y(r,t)$ (red, dashed) for non-dimensional times (1, 2, 3, 4, 5) $\times$ 103.
The transition surface location $Y=0$ except for small regions around the flow
front. (bottom) Layer height field versus the similarity variable
$\zeta=r/t^{1/2}$. Away from the source, solution convergence onto the
similarity form indicates that the fluid behaves isoviscously with $\eta_{\rm
eff}\approx\eta$.

For intermediate transition stress ($\tau_{c}=30\times 10^{-4}$), the location
of the transition surface is initially approximately half of the layer height
but decreases as evolution proceeds (Fig. 4). In this case, it is not possible
to collapse the solution onto a similarity form.

Figure 4: As Fig. 3 with $\tau_{c}=30\times 10^{-4}$. (top) Height field and
pseudoplastic transition surface located at approximately half the layer
height. (bottom) Height field as a function of the similarity variable
$\zeta$. For intermediate transition stress, a similarity solution of the
isoviscous form is not admitted.

For small transition stress ($\tau_{c}=3\times 10^{-4}$), the transition
surface location is approximately equal to layer height (Fig. 5). Scaling
layer height and radius by the appropriate viscosity ratio factors (Sec. 2)
indicates that the solution is converging to an isoviscous similarity solution
with high strain rate viscosity $\eta_{\rm eff}\approx\mu$.

Figure 5: As Fig. 3 with $\tau_{c}=3\times 10^{-4}$. (top) Height field and
pseudoplastic transition surface with location approximately equal to the
layer height. (bottom) Height field scaled by the viscosity ratio factor
$\chi^{1/4}$ as a function of the similarity variable $\zeta$ scaled by the
viscosity ratio factor $\chi^{-1/8}$. The approximate convergence of the
solution onto the similarity form suggests the flow is effectively isoviscous
with $\eta_{\rm eff}\rightarrow\mu$.

####  4.1.2 Flow front

In Fig. 6, flow front location as a function of time is shown for viscosity
ratio $\chi$ = 10 and three values of transition stress. The flow front is
defined by $h(r_{f},t)$ = 0.01. For sufficiently large $\tau_{c}$, flow is
controlled by the low strain rate viscosity, i.e., $\eta_{\rm eff}\approx\eta$
and flow radius exhibits the asymptotic, isoviscous, time scaling. Also, for
small transition stress, the flow behaves isoviscously with $\eta_{\rm
eff}\approx\mu$ throughout the evolution time considered. In the case of
intermediate $\tau_{c}$, flow radius behavior exhibits a transition.
Initially, layer evolution is consistent with the $\eta_{\rm eff}\approx\mu$
regime. As the layer spreads, the stress invariant decreases, and the behavior
approaches that for the small strain rate limit.

Scaling analysis (Sec. 2) suggests that, for spreading initially controlled by
$\mu$, shear stress decreases with time like $t^{-1/2}$. That is, the
dimensional stress invariant,

$\tau\sim\left(\rho^{3}g^{3}\mu^{5}Q\right)^{1/8}\;t^{-1/2}\;.$ (4.2)

A flow transition time scale $t_{c}$ can be defined as the time when
$\tau\sim\tau_{c}$

$t_{c}\sim\frac{\left(\rho^{3}g^{3}\mu^{5}Q\right)^{1/4}}{\tau_{c}^{2}}\;.$
(4.3)

Non-dimensionalizing (Sec. 3.2.1),

$t_{c}\sim\frac{Q^{1/4}}{\chi^{5/4}\tau_{c}^{2}}\;.$ (4.4)

For non-dimensional source function Eq. (4.1),
$Q=\dfrac{\pi}{2}w_{0}r_{0}^{4}$

$t_{c}\sim 10^{3}\,\left(\frac{10}{\chi}\right)^{5/4}\,\left(\frac{30\times
10^{-4}}{\tau_{c}}\right)^{2}\;.$ (4.5)

Thus, the intermediate transition stress case is expected to undergo this flow
transition during the non-dimensional evolution time of the calculation
$T=5\times 10^{3}$. The low transition stress case considered would not
exhibit such behavior until a non-dimensional time $t\sim 10^{2}\,T$.

For fixed transition stress and evolution time, increasing viscosity ratio
$\chi$ implies increasing flow radius (Eq. 2.7). For $\tau_{c}=3\times
10^{-4}$ and $T=5\times 10^{3}$, flow radius was calculated as a function of
$\chi$. Results agree with the scaling analysis (Fig. 6).

Figure 6: (left) Variation in flow front radius with time for $\chi$=10 and
pseudoplastic transition stress $\tau_{c}=(300,30,3)\times 10^{-4}$ (blue,
greed, red). The black dashed line shows the characteristic, isoviscous
scaling. (right) Flow front location after fixed evolution time $T=5\times
10^{3}$ for $\tau_{c}=3\times 10^{-4}$ as a function of viscosity ratio
$\chi$. Results are in good agreement with the expected asympotic scaling
(green dashed line).

###  4.2 Spherical geometry

In spherical geometry, the non-dimensional source function

$w_{s}(\theta)=w_{0}\,(\theta_{0}^{2}-\theta^{2})\,H(\theta_{0}-\theta)$ (4.6)

with $\theta_{0}$ = 0.15 and $w_{0}$ = 0.1. This source corresponds to non-
dimensional volume flux $Q$ identical to the planar geometry volume flux to
$\mathcal{O}(\theta_{0}^{5})$. In the following sections, the effects of
varying rheological transition stress and viscosity ratio are investigated
numerically.

####  4.2.1 Height field

The range of rheological transition stresses investigated in spherical
geometry is the same as that for the planar geometry case. Likewise, the
viscosity ratio $\chi$ = 10. Results are summarized in Fig. (7). For large
transition stress ($\tau_{c}=300\times 10^{-4}$), the location of the
rheological transition surface $Y\sim 0$. In this limit, the fluid behaves
isoviscously with $\eta_{\rm eff}\approx\eta$. For intermediate transition
stress ($\tau_{c}=30\times 10^{-4}$), the low viscosity part of the flow is
initially approximately half of the layer height. The low viscosity layer
height fraction decreases as evolution proceeds. For small transition stress
($\tau_{c}=3\times 10^{-4}$), $Y\sim h$.

Figure 7: Evolution in spherical geometry for $\chi=10$ and, from top to
bottom, $\tau_{c}=(300,30,3)\times 10^{-4}$, respectively. Layer height fields
$h(\theta,t)$ (blue) together with the transition surface $Y(\theta,t)$ (red,
dashed) for non-dimensional times (2, 4, 6, 8, 10) $\times$ 103. (top) The
transition surface location is indistinguishable from $Y=0$. (middle) The low
viscosity region constitutes a decreasing fraction of the layer height
throughout evolution. (bottom) The high viscosity part of the flow is only a
small fraction near the free surface.

####  4.2.2 Flow front

The flow front location is defined by $h(\theta_{f},t)$ = 0.01. As in the
planar case, the flow front location is controlled by high (low) viscosities
for large (small) rheological transition stress (Fig. LABEL:eq:sphr_front). In
the intermediate $\tau_{c}$ case, flow dynamics transition during evolution
due to decreasing low viscosity layer fraction. In spherical geometry,
isolation of the viscosity ratio effect is complicated by transition to a
converging flow front when spreading proceeds past the equatorial polar angle.
To isolate rheological influence on layer dynamics, the time for spreading to
$\theta_{f}\sim\pi/2$ is calculated as a function of viscosity ratio $\chi$.
Scaling analysis (Sec. 2) suggests that the time for spreading to polar angle
$\theta_{f}$ for $\eta_{\rm eff}\sim\eta$,

$t^{\ast}\sim\left(\frac{Q^{3}\rho g}{\eta}\right)\theta_{f}^{2}\;.$ (4.7)

In the asymptotic limit, $\eta_{\rm eff}\sim\mu$,

$t_{\mu}^{\ast}\sim\chi^{-1/4}\;t^{\ast}\;.$ (4.8)

Numerical results (Fig. 8) are in good agreement with the anticipated scaling.

Figure 8: (left) Variation in flow front polar angle with time for $\chi$=10
and pseudoplastic transition stress $\tau_{c}=(300,30,3)\times 10^{-4}$ (blue,
greed, red). The black dashed line shows the characteristic, isoviscous
scaling. (right) Time for spreading to equatorial polar angle as a function of
viscosity ratio $\chi$. Flow behavior is in agreement with asymptotic scaling
analysis (green dashed line).

Figure 9: (top,left) Height fields $h(\theta,t)$ and $h(r,t)$ in spherical
(blue) and planar (green) geometries for the isoviscous (or large transition
stress) and non-dimensional times (0.2, 0.4, 0.6, 0.8, 1, 2, 3, 4, 5) $\times$
103. (top,right) Flow front location as a function of time. For polar angle
$\theta\lesssim 0.5$, the solutions are in good agreement. (bottom) As for the
top figures for the case $\tau_{c}=3\times 10^{-4}$ and $\chi$ = 10.

###  4.3 Small polar angle approximation

In the small polar angle limit $\theta\ll 1$, the non-dimensional evolution
equation for layer height in spherical geometry reduces to that for planar
geometry. Planar and spherical solutions for identical source functions should
be equal in the small polar angle approximation. Layer height fields and flow
front location for the two geometries are compared in Fig. (9) for the
isoviscous and small transition stress cases. Solutions are in agreement for
polar angle $\theta\lesssim 0.5$. For example, the relative difference in flow
front location at this angle is $\sim$ 0.3 %

##  5 Conclusions

A pseudoplastic rheology was applied to axisymmetric thin film evolution in
planar and spherical geometries with constant volume flux source. Closed form
expressions for velocity profile and volume flux were derived. The numerical
approach utilized in the study was benchmarked against previous analytical and
numerical solutions. Influence of rheological transition strain rate and
viscosity ratio on layer evolution was explored numerically. In the limits of
large and small transition strain rate, approximately isoviscous evolution was
observed. For intermediate transition strain rate, control of layer spreading
can undergo an adjustment from high to low strain rate viscosity. Planar and
spherical geometry solutions agree for sufficiently small polar angle (see
e.g., Takagi and Huppert [2010]).

While admittedly simplified, the rheological model captures bulk features of
shear thinning and allows for efficient exploration of parameter space. It may
find applicability in geodynamical contexts including, but not limited to,
plume head evolution Bercovici and Lin [1996], isostatic adjustment to
thermochemical diapirs [Watters et al., 2009], and glacial dynamics [Schoof et
al., 2010]. The rheology is also adaptable to channelized flow [Sochi, 2015]
and, as silicic magmas exhibit shear thinning behavior [Jones et al.,, 2020;
Vasseur et al., 2023], could be implemented in models of dike emplacement
and/or magma flow in volcanic conduits [Gonnermann and Manga, 2007]. The
pseudoplastic model is readily extendable to spreading on substrates with
constant radius of curvature and vertical gravitational acceleration [Takagi
and Huppert, 2010] relevant to industrial coating applications. Finally,
modification of the model to more complex, non-planar surfaces [Lin et al.,
2012, 2021] is also a possibility.

##  6 Appendix: Numerical method

In this appendix, the numerical method is benchmarked against analytical
solutions and other numerical approaches. Layer height evolution equations
(Eqs. 3.22,3.30) are of the form

$\frac{\partial h(r,t)}{\partial t}=\mathcal{D}\left[h(r,t)\right]$ (6.1)

where $\mathcal{D}$ is a non-linear differential operator. The Python package
py-pde [Zwicker, 2020] provides methods for solving partial differential
equations of this form. The py-pde, method-of-lines scheme utilizes implicit,
Adams backward differentiation [Hindmarsh, 1983; Petzold, 1983]. To benchmark
the numerical approaches, the isoviscous, constant volume flux source case is
considered and compared to analytical and numerical solutions.

###  6.1 Planar geometry

The source function used for pseudoplasticity (Eq. 4.1) is adopted for the
isoviscous benchmark. Likewise, the initial and boundary conditions for
$h(r,t)$ are the same as those for the pseudoplastic rheology. In the limit of
a constant, point source flux $w_{s}$, the isoviscous case admits an
analytical, similarity solution [Huppert, 1982a] in the similarity variable
$\zeta=r/t^{1/2}$. Sufficiently far from the source, numerical solutions
converge to the similarity form (Fig. 10). The flow front location scales with
$t^{1/2}$ as expected from the similarity solution and asymptotic scaling
analysis [Griffiths and Fink, 1993].

Figure 10: (left,top) Isoviscous height field evolution in planar geometry.
Snapshots of the flow height profile every 103 non-dimensional times units are
shown. (left,bottom) Height field evolution scaled by $t^{1/2}$ as a function
of the similarity variable $\zeta=r/t^{1/2}$. Away from the source, solutions
converge to the similarity form. (right) Flow front location $r_{f}$ defined
as the non-dimensional height where $h(r_{f},t)=0.01$. After initial transient
behavior [Ball and Huppert, 2019], numerical results agree with the expected
asymptotic scaling.

Figure 11: Comparison of the numerical method for spherical geometry with
previous results. (left) Isoviscous height field evolution. Results are shown
for a case with $\theta_{w}=0.3$ and non-dimensional source amplitude
$w_{0}^{\prime}=7.9\times 10^{-13}$. Snapshots of the flow height profile for
non-dimensional times (0.614, 2.46, 4.77) $\times$ 1010. (green) yin-yang
grid, (red) adaptive grid, (blue) Python py-pde. (right) Flow front location
$\theta_{f}(t)$ defined as the angular location where non-dimensional height
$h(r_{f},t)=0.01\;h(0,t)$. Results are in good agreement throughout spreading.
Line colors correspond to the top figure.

###  6.2 Spherical geometry

In the absence of a similarity solution in spherical geometry, the numerical
method adopted in this work is benchmarked against previous numerical results
[Reese et al., 2010, 2011]. One previous approach uses an adaptive grid scheme
designed for nonlinear parabolic equations [Blom and Zegeling, 1994] adopted
for axisymmetric spreading on a spherical surface. Another method utilizes a
composite, overlapping [Chesshire and Henshaw, 1990] “yin-yang” grid [Kageyama
and Sato, 2004; Kageyama, 2005]. This method decomposes the sphere into two
component grids, one being a low to mid-latitude portion of a standard
spherical-polar grid and the other a rotation of the first and is explicitly
two-dimensional in ($\theta$, $\phi$). That is, in the axisymmetric,
isoviscous case (or small strain rate limit for pseudoplastic rheology), the
layer height evolution equation Eq. (3.30) reduces to

$\frac{\partial h}{\partial
t}=\frac{1}{12}\nabla^{2}_{\theta}h^{4}+w_{s}(\theta,t)\;,$ (6.2)

where $\nabla_{\theta}$ is the polar angle part of the Laplacian on the unit
sphere. The “yin-yang” grid method integrates the evolution equation for flow
thickness in time using an explicit Euler scheme and standard, centered,
second-order approximations for the full tangential Laplacian operator
$\nabla^{2}_{\theta}+\nabla^{2}_{\phi}$. The solution at the component grid
boundaries are determined by bilinear interpolation from neighboring points.
Source axisymmetry results in axisymmetric layer spreading.

To accommodate the benchmark, the constant volume flux source function is
modified to that used in Reese et al. [2010]. In that study, the source
function is a truncated gaussian,

$w_{s}(\theta)=w_{0}\;\exp(-\theta^{2}/\theta_{w}^{2})\;H(\theta_{0}-\theta)\;.$

As for the pseudoplastic cases, the initial condition $h(\theta,t=0)$ is a
small, finite, smoothly varying function. The gradient of $h$ at the edges of
the domain is set to zero.

Good agreement between methods is found for both flow front location as a
function of time and flow profile as a function of $\theta$. Also shown is the
characteristic scaling (dashed line) for the planar geometry, constant volume
flux case [Huppert, 1982a]. The flow front eventually converges on the source
antipode, $\log\theta_{f}\approx$ 0.5. For sufficiently small distance from
the anitpodal axis of symmetry, a regime transition occurs [Gratton and
Minotti, 1990; Diez et al., 1992] which appears to be resolved in the
solutions (Fig. 11).

## References

  * Algwauish [2019] Algwauish, G.M.S.A. (2019), The non-Newtonian and non-isothermal spreading of liquid domes using mathematical and numerical methods, Doctoral dissertation, Keele University.
  * Ancey [2007] Ancey, C. (2007), Plasticity and geophysical flows: A review, JṄon-Newton. Fluid Mech., 142, 4-35.
  * Ball and Huppert [2019] Ball, T.V. and Huppert, H.E. (2019), Similarity solutions and viscous gravity current adjustment times, J. Fluid Mech., 874, 285-298.
  * Balmforth and Craster [1999] Balmforth, N.J. and Craster, R.V. (1999), A consistent thin-layer theory for Bingham plastics, J. Non-Newton. Fluid Mech., 84, 65-81.
  * Balmforth et al. [2000] Balmforth, N.J., Burbridge, A.S., Craster, R.V., Salzig, J. and Shen, A. (2000), Visco-plastic models of isothermal lava domes, J. Fluid Mech., 403, 37-65.
  * Balmforth et al. [2004] Balmforth, N.J., Craster, R.V., and Sassi, R. (2004), Dynamics of cooling viscoplastic domes, J. Fluid Mech., 499, 149-182.
  * Balmforth et al. [2006] Balmforth, N.J., Craster, R.V., Rust, A.C. and Sassi, R. (2006), Viscoplastic flow over an inclined surface, J. Non-Newton. Fluid Mech., 139, 103-127.
  * Bercovici [1994] Bercovici, D. (1994), A theoretical model of cooling viscous gravity currents with temperature-dependent viscosity, Geophys. Res. Lett., 21, 1177-1180.
  * Bercovici and Lin [1996] Bercovici, D. and Lin, J. (1996), A gravity current model of cooling mantle plume heads with temperature dependent buoyancy and viscosity, J. Geophys. Res., 101, 3291-3309.
  * Bird et al. [1987] Bird, R., Armstrong, R., and Hassager, O. (1987), Dynamics of Polymeric Liquids. Wiley.
  * Blake [1990] Blake, S. (1990), Viscoplastic models of lava domes, Lava Flows and Domes: Emplacement Mechanisms and Hazard Implications (ed. J.H. Fink), IAVCEI Proceedings in Volcanology, 2, 88-128.
  * Blom and Zegeling [1994] Blom, J.G., and Zegeling, P.A. (1994), Algorithm 731: A moving grid interface for systems of one dimensional time dependent partial differential equations, ACM Trans. Math. Software, 20, 194-214.
  * Chesshire and Henshaw [1990] Chesshire, G., and Henshaw, W. (1990), Composite overlapping meshes for the solution of partial differential equations, J. Comp. Phys., 90, 1-64.
  * Diez et al. [1992] Diez, J.A., Gratton, G., and Gratton, J. (1992), Self-similar solution of the second kind for a convergent viscous gravity current, Phys. Fluids A, 4, 1148-1155.
  * Gonnermann and Manga [2007] Gonnermann, H.M. and Manga, M. (2007), The fluid mechanics inside a volcano, Ann. Rev. Fluid Mech., 39, 321-356.
  * Gratton and Minotti [1990] Gratton, J., and Minotti, F. (1990), Self-similar viscous gravity currents: Phase plane formalism, J. Fluid Mech., 210, 155-182.
  * Gratton et al. [1999] Gratton, J., Minotti, F., and Mahajan, S.M. (1999), Theory of creeping gravity currents of a non-Newtonian liquid, Phys. Rev. E, 60, 6960-6967.
  * Griffiths and Fink [1993] Griffiths, R.W. and Fink J.H. (1993), Effects of surface cooling on the spreading of lava flows and domes, J. Fluid Mech., 252, 667-702.
  * Hindmarsh [1983] Hindmarsh, A.C. (1983), ODEPACK, A systemized collection of ODE solvers, Scientific Computing (ed. R.S. Stepleman et al.), IMACS Trans. Sci. Comp. 1, 55-64.
  * Huppert [1982a] Huppert, H.E. (1982a), The propagation of two-dimensional and axisymmetric viscous gravity currents over a rigid horizontal surface, J. Fluid Mech., 121, 43-58.
  * Huppert [1982b] Huppert, H.E. (1982b), Flow and instability of a viscous current down a slope, Nature, 300, 427-229.
  * Huppert [2006] Huppert, H.E. (2006), Gravity currents: a personal perspective, J. Fluid Mech., 554, 299-322.
  * James et al. [2021] James, F., M’ Baye, M., Msheik, K. and Nguyen, D. (2021), A lubrication equation for a simplified model of shear thinning fluid, ESAIM: Proceedings and Surveys, 70, pp.147-165.
  * Jones et al., [2020] Jones, T.J., Llewellini, E.W., and Mader, H.M., (2020), The use of a shear-thinning polymer as a bubbly magma analogue for scaled laboratory experiments, J. Volcanol. Geotherm. Res., 392, 106768.
  * Kageyama and Sato [2004] Kageyama, A. and Sato, T. (2004), “Yin-Yang grid”: An overset grid in spherical geometry, Geochem. Geophys. Geosyst., 5, Q09005.
  * Kageyama [2005] Kageyama, A. (2005), Dissection of a sphere and yin-yang grids, J. of the Earth Sim., 3, 20-28.
  * Kang et al. [2016] Kang, D., Nadim, A., and Chugnova, M. (2016), Dynamics and equilibria of thin viscous coating films on a rotating sphere, J. Fluid Mech., 791, 495-518.
  * Kang et al. [2017] Kang, D., Nadim, A., and Chugnova, M. (2017), Marangoni effects on a thin liquid film coating a sphere with axial or radial thermal gradients, Phys. Fluids, 29, 1-15.
  * Kerr and Mériaux [2004] Kerr, R.C. and Mériaux, C. (2004), Structure and dynamics of sheared mantle plumes, Geochem. Geophys. Geosyst., 5, Q12009.
  * Lavallée et al. [2007] Lavallée, Y., Hess, K.U., Cordonnier, B., and Dingwill, D.B. (2007), Non-Newtonian rheological law for highly crystalline dome lavas, Geology, 35, 843-846.
  * Lin et al. [2012] Lin, T-S, Kondic, L., and Filipov, A. (2012), Thin films flowing down inverted substrates: Three dimensional flow, Phys. Fluids, 24, 022105.
  * Lin et al. [2021] LIn, T-S, Dijksman, J.A., and Kondic, L. (2021), Thin liquid films in a funnel, J. Fluid Mech., 924, A26.
  * Myers [2005] Myers, T.G., (2005), Application of non-Newtonian models to thin film flow, Phys. Rev. E, 72, 066302.
  * Neossi Nguetchue and Momoniat [2007] Neossi Nguetchue, S.N. and Momoniat, E. (2007), Axisymmetric spreading of a thin power law fluid under gravity on a horizontal plane, Nonlinear Dyn., 52, 361-366.
  * Oron et al. [1997] Oron, A., Davis, S.H., and Bankoff, G. (1997), Long-scale evolution of thin liquid films, Rev. Modern Phys., 69, 931-980.
  * Perazzo and Gratton [2003] Perazzo, C.A. and Gratton, J. (2003), Thin film of non-Newtonian fluid on an incline, Phys. Rev. E, 67, 016307.
  * Petzold [1983] Petzold, L. (1983), Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations, SIAM J. Sci. Stat. Comput., 4, 136-148.
  * Pritchard et al. [2015] Pritchard, D., Duffy, B.R., and Wilson, S.K. (2015), Shallow flows of generalised Newtonian fluids on an inclined plane, J. Eng. Math., 94, 115-133.
  * Reese et al. [2010] Reese, C.C., Orth, C.P., and Solomatov, V.S. (2010), Impact origin for the Martian crustal dichotomy: Half emptied or half filled?, J. Geophys. Res., 115.
  * Reese et al. [2011] Reese, C.C., Orth, C.P., and Solomatov, V.S. (2011), Impact megadomes and the origin of the martian crustal dichotomy, Icarus, 213, 433-442.
  * Roy et al. [2002] Roy, R.V., Roberts, A.J., and Simpson, M.E. (2002), A lubrication model of coating flows over a curved substrate in space, J. Fluid Mech., 454, 235-261.
  * Schoof et al. [2010] Schoof, C., and Hindmarsh, C.A. (2010), Thin-film flows with wall slip: An asymptotic analysis of higher order glacier flow models, Q. J. Mech. Appl. Math., 63, 73-114.
  * Simpson [1999] Simpson, J.E. (1999), Gravity Currents: In the Environment and the Laboratory, 2nd ed. Cambridge University Press.
  * Sochi [2015] Sochi, T. (2015), Analytical solutions for the flow of Carreau and Cross fluids in circular pipes and thin slits, Rheol. Acta, 54, 745-756.
  * Stasuik et al. [1993] Stasuik, M.V., Jaupart, C., and Sparks, R.S.J. (1993), Influence of cooling on lava flow dynamics, Geology, 21, 335-338.
  * Takagi and Huppert [2010] Takagi, D. and Huppert, H.E. (2010), Flow and instability of thin films on a cylinder and sphere, J. Fluid Mech., 647, 221-238.
  * Taranets [2019] Taranets, R.M. (2019), Thin-film flows on a solid sphere, J. Math. Sci, 243, 493-504.
  * Ungarish [2009] Ungarish, M. (2009), An introduction to Gravity Currents and Intrusions. CRC Press.
  * Vasseur et al. [2023] Vasseur, J., Wadsworth, F.B., and Dingwell, D.B. (2023), Shear thinning and brittle failure in crystal-bearing magmas arise from local non-Newtonian effects in the melt, Earth Planet. Sci., 603, 117988.
  * Watters et al. [2009] Watters, W.A., Zuber, M.T., and Hager, B.H. (2009), Thermal perturbations caused by large impacts and consequences for mantle convection, J. Geophys. Res., 114, E02001.
  * Wrobel [2020] Wrobel, M. (2020), An efficient algorithm of solution for the flow of generalized Newtonian fluid in channels of simple geometries, Rheol. Acta, 59, 651-663.
  * Zwicker [2020] Zwicker, D. (2020), py-pde: A Python package for solving partial differential equations, J. Open Source Softw., 5(48), 2158.


# Decay properties of the $d^{*}(2380)$ hexaquark multiplet

M. Bashkanov<EMAIL_ADDRESS>Department of Physics, University
of York, Heslington, York, Y010 5DD, UK  G. Clash  Department of Physics,
University of York, Heslington, York, Y010 5DD, UK  M. Mocanu  Department of
Physics, University of York, Heslington, York, Y010 5DD, UK  M. Nicol
Department of Physics, University of York, Heslington, York, Y010 5DD, UK
D.P. Watts  Department of Physics, University of York, Heslington, York, Y010
5DD, UK

###### Abstract

The recently discovered $d^{*}(2380)$ hexaquark is expected to be the lightest
member of an extended SU(3) antidecuplet of hexaquark states. The experimental
search for the other heavier and strange partners of the $d^{*}(2380)$ in the
antidecuplet is a challenging task. Evaluating the most appropriate
methodologies necessitates some understanding of the underlying properties of
the decuplet states such as mass, width and decay branches. In this paper we
provide estimates of these key properties for all decuplet states,
extrapolating from information and insights garnered for the $d^{*}(2380)$.
The predictions form a basis for the design of future discovery experiments.

##  I  Introduction

The properties of the six-quark containing $d^{*}(2380)$ (hexaquark) have been
established quite rigorously in recent years following its observation in
proton-neutron scattering and pionic fusion reactions [1, 2, 3, 4, 5, 6, 7, 8,
9]. It has a mass of $M_{d^{*}}=2380$ MeV, vacuum width $\Gamma=70$ MeV and
quantum numbers $I(J^{P})=0(3^{+})$. All of the strong decay branches have
been identified and measured in experiment [9]. The electromagnetic properties
of the $d^{*}(2380)$ were also investigated recently from measurements of its
photoexcitation from deuteron targets [10, 11, 12], with further programmes
planned [13]. The particle has also recently been shown to have potential
impact in astrophysics[14, 15, 16].

The search for a strange SU(3) partners of the $d^{*}(2380)$ is a natural next
step for hexaquark studies. Characterisation of additional members of the
antidecuplet would provide valuable insights into the underlying physics of
hexaquark systems. However, as there is an abundance of states and associated
decay channels from which to focus the initial searches, some prioritisation
needs to take place. In this paper we develop a simple model to identify the
most appropriate final state reaction channels for all members of the
antidecuplet, including estimates of the partial decay widths.

##  II $d^{*}(2380)$ multiplet

From the unitary group theory of the strong interaction for the light quark
(u,d,s) sector any strongly interacting particle, such as the $d^{*}(2380)$,
should be part of an SU(3) multiplet. For the case of the $d^{*}(2380)$ it
would be expected to be a member of an antidecuplet, Fig. 1. The spectroscopic
study of the other multiplet members would provide important new constraints
on the $d^{*}$ internal structure, complimentary to that achievable in Form-
Factor studies. In a simple molecular picture, the $d^{*}$ SU(3) multiplet
would derive from the coupling of two baryon decuplet members bound by long-
range pion exchange: corresponding to $\Delta\Delta$ for the $d^{*}(2380)$ to
$\Delta\Omega$ for the $d_{sss}$. However, since the pion does not have a
coupling to strange quarks, there is an expectation in such molecular systems
that the binding energy should decrease with increasing strangeness content of
the state. Conversely in a compact hexaquark (non-molecular) picture the
binding energy should increase with increasing strangeness content as the
presence of heavier s-quarks would imply stronger binding.

Figure 1: $d^{*}(2380)$ multiplet

Regardless of the inferred structure the $10^{*}$ antidecuplet decay
properties are largely driven by the SU(3) symmetry and available phase space
for the decay products. Being a $10^{*}$ SU(3) antidecuplet, the decay
couplings are limited to either octet+octet, $8\oplus 8$ or decuplet+decuplet,
$10\oplus 10$ baryons. The possible decay branches, deriving from the
established $8\oplus 8$ and $10\oplus 10$ baryonic members members are
summarised in table 1.

Table 1: Expected decay branches of the $d^{*}(2380)$ SU(3) multiplet Particle | $8\oplus 8$ | $10\oplus 10$
---|---|---
$d^{*}$ | $pn$ | $\Delta\Delta$
$d_{s}$ | $\frac{1}{\sqrt{2}}(N\Lambda-N\Sigma)$ | $\Delta\Sigma^{*}$
$d_{ss}$ | $\frac{1}{\sqrt{6}}(\sqrt{2}N\Xi-\Sigma\Sigma+\sqrt{3}\Sigma\Lambda)$ | $\frac{1}{\sqrt{3}}(\sqrt{2}\Delta\Xi^{*}+\Sigma^{*}\Sigma^{*})$
$d_{sss}$ | $\Sigma\Xi$ | $\frac{1}{\sqrt{2}}(\Delta\Omega-\Sigma^{*}\Xi^{*})$

##  III Formalism

In our calculations we assume that the Breit-Wigner width of all the
resonances is energy dependent and that the coupling constants for the decays
into $8\oplus 8$ and $10\oplus 10$ are independent of the decaying particle’s
hypercharge/strangeness. This enables the total width of the decaying state to
be expressed as

$\Gamma_{tot}=\Gamma_{8}+\Gamma_{10},$ (1)

with $\Gamma_{8}$ corresponding to the partial width for hexaquark decay into
$8\oplus 8$ and $\Gamma_{10}$ for $10\oplus 10$. In case of $8\oplus 8$ decay
the final state particles are stable against strong decay so the $\Gamma_{8}$
decay width can be expressed as

$\Gamma_{8}=g_{8}^{2}p^{2L+1}F_{8}(p)$ (2)

with

$F_{8}(p)=\frac{R^{2L}}{1+R^{2L}p^{2L}}$ (3)

where $g_{8}$ is the coupling constant of the hexaquark decay into $8\oplus
8$, $p$ is the momentum of the ejectile in the hexaquark rest frame, $L$ is
the angular momentum in the system and $F(p)$ is a Form-Factor, usually
introduced to account for potential barriers, similar to Ref. [17, 18]. For a
$J^{p}=3^{+}$ particle decaying into two $J^{p}=1/2^{+}$ there are two
possible scenarios - a ${}^{3}D_{3}$ partial wave with $L=2$ and ${}^{3}G_{3}$
partial wave with $L=4$. From partial wave analysis of the $d^{*}(2380)\to pn$
reaction we know that the majority ($\sim 90\%$) of the $8\oplus 8$ decay
proceeds through the ${}^{3}D_{3}$ [7, 8, 19]. In our calculation we will
assume that all $8\oplus 8$ proceed with $L=2$. This assumption leads to minor
corrections, but allows a significant reduction in the required number of
coupling constants. Since all the $3^{+}$-hexaquarks lie far above the
$8\oplus 8$ threshold the $\Gamma_{8}$ width is expected to be nearly constant
throughout the resonance.

The $10\oplus 10$ channel is a lot more challenging since the baryon decuplet
contains resonant states with rather sizeable width an associated strong
energy dependence. The $\Gamma_{10}$ width is expected to vary significantly
between different resonances. We have calculated the $\Gamma_{10}$ as follows

$\displaystyle\Gamma_{10}=\gamma_{10}\int
dm_{1}^{2}dm_{2}^{2}F(q_{10})|D_{D_{1}}(m_{1}^{2})D_{D_{2}}(m_{2}^{2})|^{2}$
(4) $\displaystyle
F(q_{10})=\frac{\Lambda^{2}}{\Lambda^{2}+q^{2}_{10}/4},\;\;\;\Lambda=0.16~{}GeV/c$
(5) $\displaystyle
D_{D}=\frac{\sqrt{m_{D}\Gamma_{D}(q_{M})/q_{M}}}{M_{BM}^{2}-m_{D}+im_{D}\Gamma_{D}^{tot}(q_{M})}$
(6)
$\displaystyle\Gamma_{D}=\gamma(q_{M})^{3}\frac{R^{2}}{1+R^{2}(q_{M})^{2}}$
(7)

here $F(q_{10})$ is a Form-Factor which depends on the relative momentum
($q_{10}$) between the two decuplet baryons in the hexaquark decay. We follow
the prescription of Ref. [17], and parametrize this dependence in a monopole
form with a cut-off parameter $\Lambda$ [20]. $D_{D_{1}}/D_{D_{2}}$ are the
propagators for the decuplet of baryons, with $m_{1}/m_{2}$ the baryon masses
(or correspondingly the invariant mass of the meson-baryon system $M_{BM}$
from the decuplet baryon decay into the octet of mesons and the octet of
baryons) and $m_{D}$ being their nominal Breit-Wigner masses. The energy
dependent width of the baryon decuplet decay, $\Gamma_{D}$, is parametrised in
a standard form taking a $P$-wave decay resonance with Blatt-Weisskopf barrier
factors of $R=6.3$ GeV/c and $\gamma=0.74,0.4,0.14$ for the
$\Delta,\Sigma^{*}$ and $\Xi^{*}$ respectively [21]. $\Gamma_{D}^{tot}$ is the
sum of all partial widths. The width of the $\Omega$ is considered to be zero,
since it is stable with respect to strong decays. The Form-Factor parameters
are fixed based on the $d^{*}\to\Delta\Delta\to d\pi\pi$ invariant mass
distributions of Ref. [2]. The $g_{8}$ constant is fixed to the value
extracted from $\Gamma(d^{*}\to pn)=8$ MeV. The $\gamma_{10}$ is taken as a
normalisation factor and is fixed to reproduce
$\Gamma_{10}=\Gamma_{D_{1}}+\Gamma_{D_{2}}$ at zero binding energy, e.g
$\Gamma(d^{*}\to\Delta\Delta)=2\cdot\Gamma_{\Delta}$ for $M_{d^{*}}=2\cdot
M_{\Delta}$.

##  IV Results

We firstly explored the validity of the adopted cut off parameter
$\Lambda=0.16$ GeV/c in the model. The Form-Factor of the form Eq. 5. with a
cut-off parameter $\Lambda=0.16$ GeV/c in the $d^{*}\to\Delta\Delta\to
d\pi\pi$ was first introduced in a Ref. [2] to explain the so-called ABC-
effect, an enhancement in $M_{\pi\pi}$ close to threshold. Indeed, for the
case where the nucleons in the deuteron have a very small relative momentum,
there is a correspondence in the relative momentum between the $\Delta$’s and
the relative momentum between the pions (and hence with the pion invariant
mass). It was later speculated by A. Gal, that reduction of the $\Delta\Delta$
system size within the $d^{*}$ can lead to a further reduction of the $d^{*}$
width [21]. To clarify the situation we have studied the predicted
$d^{*}\to\Delta\Delta$ width dependence as a function of the cut of parameter
$\Lambda$, Fig. 2. The results indicate that the adopted $\Lambda=0.16$ GeV/c
can not only reproduce the ABC effect, but also gives agreement with the
measured $\Gamma(d^{*}\to\Delta\Delta)=62$ MeV. As evident in Fig 2 we also
reproduce the trend in which a higher value adopted for the cut-off parameter
leads to a smaller width.

Figure 2: $d^{*}(2380)$ width as a function of the Form Factor cut-off
parameter $\Lambda$ (left) and partial width $\Gamma_{d^{*}\to\Delta\Delta}$
as a function of binding energy (right) for various $\Lambda$ values,
$\Lambda=0.1,0.16,0.2,0.5,1$ GeV/c in a rainbow order from violet to red.

We therefore adopt $\Lambda=0.16$ GeV/c for subsequent calculations. The
predicted width for all decuplet members is shown on Fig. 3. In Fig. 4 we show
the predicted branching ratios as a function of binding energy. The results
are summarised in Table. 2

Figure 3: $d^{*}$ multiplet total width as a function of binding energy
(relative to the Decuplet-Decuplet pole) for the $d^{*},d_{s},d_{ss},d_{sss}$
from left to right split into major decay branches (note the log scale). The
vertical dashed line (common to all figures) shows the nominal expected mass,
obtained under the assumption of the same binding for all multiplet members.

Figure 4: $d^{*}$ multiplet branching ratio as a function of binding energy (relative to Decuplet-Decuplet pole) for the $d^{*},d_{s},d_{ss},d_{sss}$ from left to right. The vertical dashed line (common to all figures) shows the nominal expected mass, obtained under the assumption of the same binding for all multiplet members. Table 2: $d^{*}$ multiplet width results for the $B=-84$ MeV.  $d^{*}$ | $d_{s}$ | $d_{ss}$ | $d_{sss}$
---|---|---|---
decay |  $\Gamma$, [MeV] | BR [%] | decay |  $\Gamma$, [MeV] | BR [%] | decay |  $\Gamma$, [MeV] | BR [%] | decay |  $\Gamma$, [MeV] | BR [%]
$\Delta\Delta$ | 63 | 89 | $\Delta\Sigma^{*}(\Sigma^{*}\to\Lambda\pi)$ | 21 | 69 | $\Sigma^{*}\Sigma^{*}(\Lambda\Lambda\pi\pi)$ | 0.6 | 3.8 | $\Xi^{*}\Sigma^{*}(\Sigma^{*}\to\Lambda\pi)$ | 0.34 | 3
$pn$ | 8 | 11 | $\Delta\Sigma^{*}(\Sigma^{*}\to\Sigma\pi)$ | 1.7 | 6 | $\Sigma^{*}\Sigma^{*}(\Lambda\Sigma\pi\pi)$ | 0.05 | 0.3 | $\Xi^{*}\Sigma^{*}(\Sigma^{*}\to\Sigma\pi)$ | 0.0007 | 0.05
|  |  | $N\Lambda$ | 4.1 | 13 | $\Sigma^{*}\Sigma^{*}(\Sigma\Sigma\pi\pi)$ | 4.06e-6 | 0 | $\Omega\Delta$ | 5.2 | 41
|  |  | $N\Sigma$ | 3.7 | 12 | $\Delta\Xi^{*}$ | 7.8 | 48.9 | $\Xi\Sigma$ | 7.1 | 56
|  |  |  |  |  | $N\Xi$ | 2.6 | 16.4 |  |  |
|  |  |  |  |  | $\Sigma\Sigma$ | 1.1 | 7.1 |  |  |
|  |  |  |  |  | $\Sigma\Lambda$ | 3.8 | 23.5 |  |  |
total | 71 |  | total | 30.5 |  | total | 15.95 |  | total | 12.6 |

One can clearly see that the $8\oplus 8$ decays are predicted to be
increasingly important for the higher strangeness states, while both partial
and total widths reduces substantially. For all hexaquark members only
$10\oplus 10$ with $\Delta$ in a final state are important. For the $d_{ss}$
non-$\Delta$ channels ($\Sigma^{*}\Sigma^{*}\to$ anything) cover 4% only. For
the $d_{sss}$ the non-$\Delta$ ($\Xi^{*}\Sigma^{*}$) decay has 3% branching
ratio. It is interesting to note that the semi-electromagnetic decay branch
$d_{sss}\to\Omega\Delta\to\Omega N\gamma$ is predicted to have a higher
probability than the purely hadronic
$d_{sss}\to\Xi^{*}\Sigma^{*}\to\Xi\Sigma\pi\pi$ decay.

Hexaquarks are expected to be produced copiously in heavy ion collisions,
however our estimations indicate the width for all $d^{*}$ multiplet states is
rather large. Unfortunately this indicates their clean identification in the
tough background conditions in typical heavy ion collisions may be
challenging. The most feasible channel for such studies, having both a large
partial width and a convenient isospin $3/2$ separation, is the $\Xi\Sigma$
branch in $d_{sss}$ decay. This could potentially be tested using $\Xi-\Sigma$
correlation functions.

For all channels the most prominent final states have high particle
multiplicity, indicating hermetic detector apparati and exclusivity conditions
would need to be established. The $d_{sss}$ is the only member with a dominant
octet-octet decay, however the necessity of associated kaon production to
conserve strangeness in the production mechanism also provides challenges.The
requirements may be more easily achieved using beams containing intrinsic
strangeness (e.g. Kaon beams as proposed in Ref. [22]).

##  V Summary

We have developed a theoretical model, employing experimentally constrained
parameters, to predict the possible decay branches and partial widths for all
members of the $d^{*}$ hexaquark antidecuplet. For all strange-quark
containing members of the antidecuplet the predicted widths are rather large,
with the most promising decay channels including the broad $\Delta$ resonance
in the final state. We demonstrated that a $d^{*}$ Form-Factor, which was
first introduced to explain peculiarities in the $d^{*}\to d\pi\pi$ decay, can
also explain the smallness of the $d^{*}$ width, in agreement with the
qualitative dimensional arguments of A. Gal. The results of the paper will be
an important guide for the ongoing search for the $d^{*}$ anti-decuplet
members employing photon-, pion-, and kaon-induced reactions, as well as in
high energy collider experiments.

###### Acknowledgements.

This work has been supported by the U.K. STFC ST/V002570/1, ST/P004008/1 and
Strong2020 grants.

## References

  * [1] M. Bashkanov et al., Phys. Rev. Lett. 102, 052301 (2009).
  * [2] P. Adlarson et al., Phys. Rev. Lett. 106, 242302 (2011).
  * [3] P. Adlarson et al., Phys. Lett. B 721, 229
  * [4] P. Adlarson et al., Phys. Rev. C 88, 055208 (2013).
  * [5] P. Adlarson et al., Phys. Lett. B 743, 325 (2015).
  * [6] P. Adlarson et al., Eur. Phys. J. A 52, 147 (2016).
  * [7] P. Adlarson et al., Phys. Rev. Lett. 112, 202301 (2014).
  * [8] P. Adlarson et al., Phys. Rev. C 90, 035204 (2014).
  * [9] M. Bashkanov, H. Clement, T. Skorodko, Eur. Phys. J. A 51 (2015) 7, 87.
  * [10] M. Bashkanov et al., Phys. Lett. B789, 7, (2019).
  * [11] M. Bashkanov et al., Phys. Rev. Lett. 124, 132001, (2020)
  * [12] M. Bashkanov, D.P.Watts, A.Pastore, Phys. Rev. C 100, 012201 (2020).
  * [13] D.P.Watts, M.Bashkanov, M. Ostrick,”Photodsintegration of transversely polarised deuterons”, A2 proposal, (2020)
  * [14] I. Vidaña, M. Bashkanov, D.P. Watts, A. Pastore, Phys. Lett. B 781, 112-116, (2018).
  * [15] A. Mantziris, A. Pastore, I. Vidaña, D.P. Watts, M. Bashkanov, A. Romero, Astron.Astrophys. 638, A40, (2020).
  * [16] M. Bashkanov, D.P. Watts, J. Phys. G 47, no.3, 03LT01, (2020).
  * [17] A. Pricking, M. Bashkanov, H. Clement arxiv:1310.5532
  * [18] V. Kukulin et al., Eur. Phys. J. A 56 (2020) no.9, 229
  * [19] P. Adlarson et al., Phys. Rev. C 102, 015204, (2020).
  * [20] The appropriateness of the value for this cutoff parameter is explored in Ref [2], where it was shown to provided agreement with measured invariant mass distributions in the region of the so-called ABC-effect.
  * [21] A. Gal, Phys. Lett. B769, 436, (2017).
  * [22] M. Amaryan et al., arXiv:2008.08215.