# On the Generalization of the Hébraud -Lequeux Model to Multidimensional Flows

Julien Olivier, Michael Renardy

► To cite this version:

Julien Olivier, Michael Renardy. On the Generalization of the Hébraud -Lequeux Model to Multidimensional Flows. Archive for Rational Mechanics and Analysis, 2013, 208 (2), pp.569-601. 10.1007/s00205-012-0603-7. hal-01263132

HAL Id: hal-01263132

https://hal.science/hal-01263132

Submitted on 27 Jan 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the Generalization of the Hébraud-Lequeux Model to Multidimensional Flows*

DRAFT

OLIVIER Julien† RENARDY Michael

December 19, 2012

## Abstract

In this article we build a model for multidimensional flows based on the idea of Hébraud and Lequeux for soft glassy materials. The construction of the model is based on the ideas of HÉBRAUD and LEQUEUX but care is taken to build a frame indifferent multi-dimensional model. The main goal of this article is to prove that the methodology we have developed to study the well-posedness and the glass transition for the original Hébraud-Lequeux model can be successfully generalized. Thus this work may be used as a starting point for more sophisticated studies in the modeling of general flows of glassy materials.

## Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>1.1</td><td>Soft glasses</td><td>2</td></tr><tr><td>1.2</td><td>Scalar and tensorial model</td><td>3</td></tr><tr><td>1.3</td><td>The model</td><td>3</td></tr><tr><td>2</td><td>Main results</td><td>6</td></tr><tr><td>2.1</td><td>Hypotheses</td><td>6</td></tr><tr><td>2.2</td><td>Well-posedness and regularity</td><td>7</td></tr><tr><td>2.3</td><td>Glass transition</td><td>8</td></tr><tr><td>3</td><td>Construction of the model</td><td>9</td></tr><tr><td>3.1</td><td>The simple shear Hébraud-Lequeux model</td><td>9</td></tr><tr><td>3.2</td><td>Frame invariance</td><td>10</td></tr><tr><td>3.3</td><td>Full tensorial modeling</td><td>11</td></tr></table>

<footnote>*partially supported by ANR project Maniphyc ANR-08-SYSC-010 and partially funded by the European Union's Seventh Framework Programme (FP7-REGPOT-2009-1) under grant agreement n 245749. Supported by the National Science Foundation under Grant DMS-1008426</footnote>

<footnote>†The author would like to thank D. BRESCH and T. COLIN for fruitful discussions on this paper.</footnote>
<table>
  <tbody>
    <tr>
      <td>4</td>
      <td>Well-posedness of the model for steady homogeneous flow</td>
      <td>12</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Proof of Propositions 1 and 2</td>
      <td>13</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Study of the function F when a or b become zero</td>
      <td>15</td>
    </tr>
    <tr>
      <td>4.2.1</td>
      <td>Study of the function b ↦ F(0, b)</td>
      <td>15</td>
    </tr>
    <tr>
      <td>4.2.2</td>
      <td>Study of the function a ↦ F(a, 0)</td>
      <td>16</td>
    </tr>
    <tr>
      <td>4.2.3</td>
      <td>Study of the case μ = μ<sub>c</sub></td>
      <td>18</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Verification of Hypothesis H2</td>
      <td>25</td>
    </tr>
    <tr>
      <td>5</td>
      <td>Proof of Corollary 1</td>
      <td>27</td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>The case μ &gt; μ<sub>c</sub></td>
      <td>27</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>The case μ &lt; μ<sub>c</sub></td>
      <td>27</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>The case μ = μ<sub>c</sub></td>
      <td>27</td>
    </tr>
    <tr>
      <td>6</td>
      <td>The form of A and A<sub>0</sub> and examples</td>
      <td>28</td>
    </tr>
    <tr>
      <td>6.1</td>
      <td>The form of the drift term</td>
      <td>28</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Couette flow</td>
      <td>28</td>
    </tr>
    <tr>
      <td>6.3</td>
      <td>Elongational flows</td>
      <td>29</td>
    </tr>
  </tbody>
</table>

# 1 Introduction

## 1.1 Soft glasses

In this article we are interested in the modeling of soft glassy materials. Soft glassy materials are named after glasses because they exhibit similar behavior: very complicated free energy landscapes and a transitional behavior akin to phase transition. Soft glassy materials can be loosely described as “particles” dispersed inside a Newtonian fluid which may be liquid or gas. When the particles are solids we obtain what are called suspensions, when the particles are liquid bubbles we obtain emulsions. Foams and granular flows may also be included in this category.

Numerous models exist to describe specific materials: one can cite the model by JOP, FORTERRE and POULIQUEN [11] for dense granular flows or the one by BÉNITO, BRUNEAU, COLIN, GAY and MOLINO [5] for foams and emulsions. We want to call these two particular models *tensorial* or *multidimensional* models because they give constitutive equations that are not subject to any geometrical constraint on the flow type: the constitutive law links the full stress tensor to kinematic quantities such as the deformation tensor, the mass density or the temperature.

On the other hand we have models for generic soft glassy materials derived from their analogy to real molecular glasses. The first to have been conceived was by P. SOLLICH, LEQUEUX, HÉBRAUD and CATES and simply called SGR (for Soft Glassy Rheology) [16] and was originally a model restricted to simple shear. This model was later generalized to a tensorial version by M. E. CATES and P. SOLLICH [8]. However, this model suffers from one drawback: the constitutive relation depends on an *ad hoc* phenomenological parameter (an effective temperature) which is hard to interpret. This is why HÉBRAUD and LEQUEUX came up with a modified version of the SGR model [9] whose parameters have a more direct physical interpretation. However, this model is designed for simple shear and this article is an attempt to give and analyze a tensorial version of this model. To avoid misunderstandings, we emphasize that, although our tensorial
model is constructed in a fashion analogous to the one-dimensional HL model, it does not reduce to this model even if the flow is parallel shear.

## 1.2 Scalar and tensorial model

Let us make precise what we mean by a tensorial or, respectively, simple shear model. It is well known that in a continuous medium there is a tensor $T$, that is a field of $d \times d$ matrices ($d$ is the space dimension) which describes the internal force that holds the medium together. This tensor is called the stress tensor. With $\rho$ denoting the mass density and $u$ the Eulerian velocity field, the conservation of momentum for an incompressible medium reads:

$$ \rho(\partial_t u + u \cdot \nabla u) = \operatorname{div} T - \nabla p, $$

where the divergence of a matrix tensor is the vector of the divergences of its columns. Here $p$ denotes the pressure, which is implicitly determined by the incompressibility constraint

$$ \operatorname{div} u = 0. $$

A constitutive law for a continuous material is a relation between $T$ and kinematic quantities. For instance Newton's law of simple fluids reads

$$ T = 2\mu D(\nabla u), \qquad (1) $$

where $\mu$ is the viscosity, and

$$ D(\nabla u) = \frac{1}{2}(\nabla u + (\nabla u)^T) $$

is the deformation rate tensor.

On the other hand, rheological measurements often emphasize simple shear flow situations, and constitutive laws are often formulated and verified against experiments specifically for such flows. If $x$ is the shear direction and $y$ is its orthogonal, then the velocity field has the special form $u = (\dot{\gamma}(t)y, 0, 0)$. In this configuration, $\dot{\gamma}(t)$ is the shear rate and has the dimension of an the inverse time. We then suppose that no quantity depends on $x$ (except for the pressure which may be a linear function of $x$) and the conservation of momentum reads:

$$ \rho \partial_t u^x = \partial_y T_{12} - \partial_x p. \qquad (2) $$

Thus only one component of the stress tensor is relevant in this situation.

In general, the rheology of a fluid in shear tells us little about its behavior in other geometries. Thus, there are many ways a model might be generalized. In our attempt, we are guided by the heuristics underlying the formulation of Hébraud and Lequeux, and we attempt to preserve the reasoning and formal structure which led to their model.

## 1.3 The model

To end this introduction, we give the model we have designed in the original, dimensionless, physical variables and then transform it by various changes of dependent and independent variables to put it in a form suited for mathematical analysis. We assume a three-dimensional flow. In general, our model is a
Fokker-Planck like equation for a probability density $p$ over the space $S_3$ of real symmetric matrices whose generic point is denoted by $\Sigma$ and represents a stress. We represent a symmetric matrix in the form

$$ \begin{pmatrix} \Sigma_1 \\ \vdots \\ \Sigma_6 \end{pmatrix} \mapsto \begin{pmatrix} \Sigma_1 & \Sigma_4/\sqrt{2} & \Sigma_5/\sqrt{2} \\ \Sigma_4/\sqrt{2} & \Sigma_2 & \Sigma_6/\sqrt{2} \\ \Sigma_5/\sqrt{2} & \Sigma_6/\sqrt{2} & \Sigma_3 \end{pmatrix}. \quad (3) $$

The Euclidean metric in this six dimensional space turns out to be invariant under transformations which correspond to a rotation in physical space. All notations like integration, Laplacian etc. will be relative to this metric. We shall denote by $\Sigma$ the vector $(\Sigma_1, \Sigma_2, ..., \Sigma_6)$, and by $\mathcal{M}(\Sigma)$ the matrix associated with it.

Our model has the following form

$$ \begin{aligned} \partial_t p(t, x, \Sigma) + (u \cdot \nabla_x)p + G(\Sigma, \nabla u) : \nabla_{\Sigma} p \\ = -\mathbf{1}_{|\Sigma|>1}(\Sigma)p(t, x, \Sigma) + \Gamma(p)(t, x)\rho(\Sigma) + \mu\Gamma(p)(t, x)\Delta_{\Sigma}p(t, x, \Sigma), \end{aligned} \quad (4) $$

where $\Gamma$ is called the fluidity and is defined by

$$ \Gamma(p)(t, x) = \int_{|\Sigma|>1} p(t, x, \Sigma)d\Sigma. \quad (5) $$

Here $\rho$ is a rotationally invariant measure supported inside the unit ball and having unit integral; in the original HL model $\rho$ is a delta function located at the origin, and we shall be particularly interested in this case. The term $G(\Sigma, \nabla u)$ represents a drift term. That is, in the absence of stochastic effects, $\Sigma$ would satisfy the equation

$$ \partial_t \Sigma + (u \cdot \nabla_x)\Sigma = G(\Sigma, \nabla u). \quad (6) $$

Below, we shall use $h(|\Sigma|)$ to denote the characteristic function of $|\Sigma| > 1$. For an explanation of the terms and their physical meaning we refer to Section 3.

All the mathematical analysis we shall do in this paper is concerned with a situation which is stationary in time and homogeneous in space. We set

$$ \nabla_x u = \varepsilon M, \quad (7) $$

where $\varepsilon$ is the deformation rate and $M$ the deformation type. Some examples of deformation types are given in Section 6. Our mathematical analysis will only be concerned with the situation where $\varepsilon$ is small. Moreover, the form of $G$ we shall use is linear with respect to $\Sigma$ and $\varepsilon$, i.e.

$$ G(\Sigma, \varepsilon M) = \varepsilon(\mathcal{A}\Sigma + \lambda\mathcal{A}^0), \quad (8) $$

where $\mathcal{A}$ is a $6 \times 6$ matrix and $\mathcal{A}^0$ a vector of $\mathbb{R}^6$ which depend on $M$. Again see Section 6 for the constructions of these elements for several deformation types. As a consequence of the incompressibility condition, the matrix $\mathcal{A}$ will always be traceless.

Consequently the physical system now reads

$$ \left\{ \begin{array}{l} -\mu\Gamma\Delta p + \varepsilon(\mathcal{A}\Sigma + \lambda\mathcal{A}^0) \cdot \nabla p + h(|\Sigma|)p = \Gamma\rho, \\ \int_{\Sigma \in \mathbb{R}^6} p(\Sigma)d\Sigma = 1. \end{array} \right. \quad (9) $$
Since $h$ is discontinuous, it is advantageous to split the problem into an equation inside the unit ball $B$ and another one outside. This leads to the system

$$ egin{cases} -\mu\Gamma\Delta q + \varepsilon(\mathcal{A}\Sigma + \lambda\mathcal{A}_0) \cdot \nabla q = \Gamma\rho & \text{in } B = \{|\Sigma| < 1\}, \\ -\mu\Gamma\Delta r + \varepsilon(\mathcal{A}\Sigma + \lambda\mathcal{A}_0) \cdot \nabla r + r = 0 & \text{in } |\Sigma| > 1, \\ q = r & \text{on } |\Sigma| = 1, \\ \partial_n q = \partial_n r & \text{on } |\Sigma| = 1, \end{cases} \quad (10) $$

and the integral constraint

$$ \int_B q(\Sigma) d\Sigma + \int_{\mathbb{R}^6 \setminus B} r(\Sigma) d\Sigma = \int_B q(\Sigma) d\Sigma + \Gamma = 1. \quad (11) $$

Thus $p$ is a solution to (9) if, and only if $(p|_B, p|_{\mathbb{R}^6 \setminus B})$ is a solution to (10)-(11). In this form the stress tensor attached to $(q, r)$ (or $p$) is the vector obtained by integration:

$$ T^{\varepsilon} = \int_B q(\Sigma) \mathcal{M}(\Sigma) d\Sigma + \int_{\mathbb{R}^6 \setminus B} r(\Sigma) \mathcal{M}(\Sigma) d\Sigma. \quad (12) $$

We know from our analysis in [15] that the limit $\varepsilon \to 0$ we are studying is mathematically singular but can be regularized by a change of variables. This is why we introduce the new parameters

$$ \mathfrak{a} = \frac{\varepsilon}{\mu\Gamma} \qquad \mathfrak{b} = \sqrt{\mu\Gamma} \quad (13) $$

Then we define for $\mathfrak{b} > 0$ the diffeomorphisms $\Psi_{\mathfrak{b}}$ as

$$ \Psi_{\mathfrak{b}}: \begin{cases} ]0, +\infty[ \times \partial B & \rightarrow \mathbb{R}^6 \setminus \bar{B} \\ (\theta, \eta) & \mapsto (1 + \mathfrak{b}\theta)\eta \end{cases} \quad (14) $$

Here $E$ is the function defined in the unit ball $B$, solution to the system

$$ \left\{ \begin{array}{ll} -\Delta E = \rho & \text{in } |\Sigma| < 1, \\ E = 0 & \text{on } |\Sigma| = 1. \end{array} \right. \quad (15) $$

Finally, we look at the functions $\tilde{q}$ and $\tilde{r}$ given by

$$ \tilde{q}(\Sigma) = \mu q - E, \quad (16) $$

$$ \tilde{r}(\theta, \eta) = \mu(r \circ \Psi_b)(\theta, \eta) \quad (17) $$

We shall use the notation

$$ L = (\mathcal{A}\Sigma + \lambda\mathcal{A}^0) \cdot \nabla. \quad (18) $$

Moreover, for $|\Sigma| > 1$, we shall use the formulae

$$ \Delta = \frac{1}{\mathfrak{b}^2} \partial_{\theta}^2 + \frac{5}{\mathfrak{b}(\mathfrak{b}\theta + 1)} \partial_{\theta} + \frac{1}{(\mathfrak{b}\theta + 1)^2} \Delta_{\tan}, $$

$$ \nabla_{\Sigma} = \frac{1}{1 + \mathfrak{b}\theta} \nabla_{\tan} + \frac{1}{\mathfrak{b}} \eta \partial_{\theta}. \quad (19) $$
Here $\Delta_{\text{tan}}$ and $\nabla_{\text{tan}}$ are the tangential parts of Laplacian and gradient on the unit sphere.

System (10) is now transformed into

$$ egin{cases} 
-\Delta \tilde{q} + aL\tilde{q} = -aLE, \\
-\partial_{\theta}^2 \tilde{r} - \frac{b^2}{(1+b\theta)^2} \Delta_{\text{tan}} \tilde{r} + \frac{ab^2}{1+b\theta}((1+b\theta)A\eta + \lambda A_0) \cdot \nabla_{\text{tan}} \tilde{r} \\
+ \left(-\frac{5b}{1+b\theta} + ab((1+b\theta)A\eta + \lambda A_0) \cdot \eta \right) \partial_{\theta} \tilde{r} + \tilde{r} = 0, \\
\tilde{q}(\eta) = \tilde{r}(0, \eta), \\
b\partial_n \tilde{q}(\eta) - \partial_{\theta} \tilde{r}(0, \eta) = -b\partial_n E(\eta). 
\end{cases} 
\quad (20) $$

while Eq. (11) is formally transformed into $F = \mu$ where

$$ F(a, b) = \mu_c + \int_B \tilde{q}^{a,b}(\Sigma) d\Sigma + b^2. \quad (21) $$

The outline of the paper is the following: after presenting the results in Section 2 and the construction of the model (4) in Section 3, the system of (20) and (21) is analyzed in Section 4. From it we deduce the well-posedness of (10)–(11) in Section 5. The glass transition in our model follows from the asymptotics of the solution. The formal procedure to compute coefficients in the asymptotic expansions is discussed at length in our previous paper on the one-dimensional case [15] and will not be repeated here. In Section 6 we show that the hypotheses required on the deformation types are valid for standard deformation types.

# 2 Main results

## 2.1 Hypotheses

In our analysis, we shall use the following assumptions:

**H1** We assume that $\rho$ is a nonnegative Radon measure on the unit ball which belongs to $W^{-1,\beta}(B)$ for some $\beta$ with $1 < \beta \le 2$ and has unit integral. Also, $\rho$ is rotationally invariant, and the singular support of the distribution $\rho$ is bounded away from the unit sphere.

**H1** implies that $E$ as defined by (15) is in $W_0^{1,\beta}(B)$ but is also $C^\infty$ near the unit sphere. Note also that from (15), we also have $E$ rotation invariant. We are particularly interested in the case where $\rho$ is the delta function.

**H2** Let us define $\tilde{q}_0^a$ as the $W_0^{1,\beta}(B)$ solution of

$$ egin{cases} 
-\Delta \tilde{q}_0^a + aL\tilde{q}_0^a = -aLE & \text{in } |\Sigma| \le 1, \\
\tilde{q}_0^a = 0 & \text{on } |\Sigma| = 1. 
\end{cases} 
\quad (22) $$

Then we assume that the function

$$ a \mapsto \int_{|\Sigma| \le 1} \tilde{q}_0^a(\Sigma) d\Sigma \quad (23) $$
has a negative derivative over $]0, +\infty[$. This hypothesis is in general not easy to check. We shall later verify it numerically for the case when $\rho = \delta_0$ and $\lambda$ is large.

**H3** The trajectory of the vector field $\Sigma \mapsto \mathcal{A}\Sigma + \lambda\mathcal{A}_0$ originating from every point of the open unit ball leaves the unit ball in the both directions. This is evidently true if $\mathcal{A}_0$ is nonzero and $\lambda$ is sufficiently large.

## 2.2 Well-posedness and regularity

We first study the well-posedness of the differential system (20) with the theory of elliptic equations and systems. This allows us to define the function $F(\mathbf{a}, \mathbf{b})$ of (21). Then we study the limiting cases $\mathbf{a} \to 0$ and $\mathbf{b} \to 0$ and set up an argument based on the implicit function theorem to treat the case when $\mathbf{a}$ or $\mathbf{b}$ is small. Finally, another application of the implicit function theorem allows us to analyze the behavior of (10)–(11) for small nonzero $\varepsilon$.

**Proposition 1.** For every positive $\mathbf{a}$ and $\mathbf{b}$, the system (20) has a unique solution such that $\tilde{q} \in W^{2,\beta}(B)$, $\tilde{r} \in H^1(\mathbb{R} \times \partial B)$. Moreover, this solution is of class $C^\infty$ outside the singular support of $\rho$. Finally, for any integer $n$, there exists a $C_n$ such that $\theta^n \tilde{r} \to 0$ as $\theta \to \infty$ if $ab^2 < C_n$. An analogous statement holds for derivatives of $\tilde{r}$. All bounds of $\tilde{q}$ and $\tilde{r}$ remain uniformly valid as $\mathbf{a}$ and/or $\mathbf{b}$ tends to zero. Moreover, the solution depends analytically on $\mathbf{a}$ and $\mathbf{b}$.

We are particularly interested in the limit where $\mathbf{a}$ or $\mathbf{b}$ is zero. In this case, we still have

**Proposition 2.** The solutions given by the previous proposition remain $C^\infty$ functions of $\mathbf{a}$ and $\mathbf{b}$ in the limit where $\mathbf{a}$ or $\mathbf{b}$ or both tend to zero.

We can now define the function $F$ given by

$$F(\mathbf{a}, \mathbf{b}) = \mu_c + \int_B \tilde{q}^{\mathbf{a}, \mathbf{b}}(\Sigma) d\Sigma + b^2.$$

Our goal is now to study $F(\mathbf{a}, \mathbf{b}) = \mu$, which is related to (11): the couple $(\tilde{q}^{\mathbf{a}, \mathbf{b}}, \tilde{r}^{\mathbf{a}, \mathbf{b}})$ for which the parameters satisfy $F(\mathbf{a}, \mathbf{b}) = \mu$ would exactly be the couple obtained from a solution of (10)–(11) by the change of variables given by (13), (16) and (17).

**Theorem 1.** The function $F$ has the following properties:

* $F$ is analytic for $\mathbf{a} > 0, \mathbf{b} > 0$.

* $F$ can be continued with $C^\infty$ regularity up to the boundaries $\mathbf{a} = 0$ and $\mathbf{b} = 0$.

* The function

$$b \mapsto F(0, b)$$

is monotonically increasing from $[0, +\infty[$ to $[\mu_c, +\infty[$. Consequently, for any $\mu \ge \mu_c$, there is a unique $b_0 \ge 0$ (with $b_0 = 0$ when $\mu = \mu_c$) such that $F(0, b_0) = \mu$ and when $b_0 \ge 0$, we also have $\partial_b F(0, b_0) > 0$.
• Because of hypothesis **H2** and **H3**, the function

$$a \mapsto F(a, 0)$$

is monotonically decreasing from $[0, +\infty[ \text{ to } ]0, \mu_c]$. Consequently, for any $\mu \le \mu_c$, there is a unique $a_0 \ge 0$ (with $a_0 = 0$ when $\mu = \mu_c$) such that $F(a_0, 0) = \mu$ and when $a_0 > 0$, we also have $\partial_a F(a_0, 0) < 0$.

• We have $F(0, 0) = \mu_c$. Moreover, $\partial_a F(0, 0) = 0$ and $\partial_{aa}^2 F(0, 0) < 0$.

Now we note that when our original parameter $\varepsilon$ is small, then either **a** or **b** must be small. From the above properties and the implicit function theorem, we conclude that near $a = 0$, we can solve for **b** as a function of **a** and $\mu$, while near $b = 0$, and for $a > 0$, we can solve for **a** as a function of **b** and $\mu$. From this, we can deduce the following corollary concerning the original parameters $(\varepsilon, \Gamma)$:

**Corollary 1.** • We have the following asymptotic expansions near $\varepsilon = 0$:

- if $\mu > \mu_c$,

$$\Gamma \sim \sum_{k=0}^{+\infty} \tilde{c}_k \varepsilon^k,$$

and we have $\tilde{c}_0 > 0$;

- if $\mu < \mu_c$,

$$\Gamma \sim \sum_{k=1}^{\infty} \tilde{c}_k \varepsilon^{k/2},$$

and we have $\tilde{c}_1 > 0$;

- if $\mu = \mu_c$,

$$\Gamma = \sum_{k=4}^{\infty} \tilde{c}_k \varepsilon^{k/5},$$

and we have $\tilde{c}_4 > 0$.

• For all $\mu > 0$ and $\varepsilon > 0$ small enough, there exist a unique solution to the system (10)-(11).

## 2.3 Glass transition

From the results of the previous subsection, we can deduce the asymptotic form of the solution. We shall state these results in terms of the original functions $q$ and $r$.

**Proposition 3.** Let us denote by $(q^\varepsilon, r^\varepsilon)$ the unique solution of (10)-(11). Then under hypotheses **H1**–**H3** we have asymptotic expansions of the following form:

• If $\mu > \mu_c$,
$$ \begin{cases} q(\Sigma) \sim \bar{Q}^0 + \varepsilon \bar{Q}^1 + \varepsilon^2 \bar{Q}^2 + \dots, \\ r(\Sigma) \sim \bar{R}^0 + \varepsilon \bar{R}^1 + \varepsilon^2 \bar{R}^2 + \dots, \end{cases} \quad (24) $$

where the profiles $\bar{Q}^k$ and $\bar{R}^k$ are functions of $\Sigma$.
• If $\mu < \mu_c$,
$$ \begin{cases} q(\Sigma) \sim \bar{Q}^0 + \varepsilon^{1/2}\bar{Q}^1 + \varepsilon\bar{Q}^2 + \dots, \\ r(\Sigma) \sim \varepsilon^{1/2}R^1 + \varepsilon R^2 + \dots, \end{cases} \qquad (25) $$

where the profiles $\bar{Q}^k$ are functions of $\Sigma$ and the $R^k$ are functions of $(\Sigma/|\Sigma|, (|\Sigma| - 1)/\varepsilon^{1/2})$.

• If $\mu = \mu_c$,
$$ \begin{cases} q(\Sigma) \sim \bar{Q}^0 + \varepsilon^{1/5}\bar{Q}^1 + \varepsilon^{2/5}\bar{Q}^2 + \dots, \\ r(\Sigma) \sim \varepsilon^{2/5}R^2 + \varepsilon^{3/5}R^3 + \dots, \end{cases} \qquad (26) $$

where the profiles $\bar{Q}^k$ are functions of $\Sigma$ and the $R^k$ are functions of $(\Sigma/|\Sigma|, (|\Sigma| - 1)/\varepsilon^{2/5})$.

The result of physical interest of this paper is

**Theorem 2.** Let $T^\varepsilon$ be the macroscopic stress vector computed from $(q^\varepsilon, r^\varepsilon)$ through (12) By integration we can deduce from Proposition 3 the following results:

• If $\mu > \mu_c$,
$$ T^\varepsilon \sim \lambda\kappa(\mu, \rho)\varepsilon M(\mathcal{A}^0) + \varepsilon^2 T^2 + \dots \qquad (27) $$

• If $\mu < \mu_c$,
$$ T^\varepsilon \sim T^0 + \varepsilon^{1/2}T^1 + \dots \qquad (28) $$

• If $\mu = \mu_c$,
$$ T^\varepsilon \sim \varepsilon^{1/5}T^1 + \varepsilon^{2/5}T^2 + \dots \qquad (29) $$

*Remark 1.* Note that when $\mu > \mu_c$, we have Newtonian behavior, with a viscosity $\kappa$ which depends on $\mu$ and $\rho$. For $\mu < \mu_c$, we have a yield stress, while at $\mu = \mu_c$, we have a shear thinning power law fluid. This is a generalization of a known result on the original HL model [9, 14, 15].

# 3 Construction of the model

This section is dedicated to the modeling aspects of our work. We first review the main features of the Hébraud-Lequeux (HL) model of [9] which we try to generalize. Then we explain the minimal symmetry constraint that our equations must satisfy for the required *frame indifference* property. Finally we give a possible generalization of the HL model.

## 3.1 The simple shear Hébraud-Lequeux model

We will now review the (HL) model which is the foundation of our model. The HL model comes from the mean field theory. It describes the state of the material at a mesoscopic level and then takes some average of this mesoscopic description to give the properties of the material at the macroscopic scale. From a mathematical point of view the model is given by a Fokker-Planck type equation:
$$
\left\{
\begin{array}{l}
\displaystyle \partial_t p(t, \sigma) = -G_0 \dot{\gamma}(t) \partial_{\sigma} p - \frac{1}{T_0} H(|\sigma| - \sigma_c) p(t, \sigma) \\
\displaystyle \qquad + \Gamma(p(t)) \delta(\sigma) + \alpha \Gamma(p(t)) \partial_{\sigma}^2 p(t, \sigma), \\[2ex]
\displaystyle \Gamma(p) = \frac{1}{T_0} \int H(|\sigma| - \sigma_c) p(t, \sigma) d\sigma, \\[2ex]
\displaystyle \int_{\sigma \in \mathbb{R}} p(t, \sigma) d\sigma = 1.
\end{array}
\right.
\qquad (30)
$$

Let us now describe the rationale behind this equation. First note that $p$ is a function of the time $t$ and of a variable $\sigma$ which is a “mesoscopic stress”. Then $\sigma_c$ and $\alpha$ are constants depending on the fluid: $\sigma_c$ is a stress threshold above which the fluid relaxes to a zero stress state with typical time $T_0$ (H designates the Heaviside function), and $\alpha$ acts as a control parameter: the higher $\alpha$, the less the fluid can structure itself and the more it behaves like a Newtonian fluid.

The idea behind this model is to imagine that in the material is composed of mesoscopic particles which undergo the shear rate $\dot{\gamma}(t)$. Each of these particle has its own stress. Then $p(t, \sigma)d\sigma$ is the number of particles whose stress is in an interval $d\sigma$ around $\sigma$. Now when submitted to a shear rate $\dot{\gamma}(t)$ the stress of a particle evolves with the following rule: if it is smaller than the stress threshold $\sigma_c$, then it evolves linearly in time. This behaviour is called an elastic behaviour since it mimics the behaviour of macroscopic elastic Hookean solids. If the stress grows beyond $\sigma_c$, then the particle will enter into a relaxation phase and its stress may instantaneously drop to 0 in a random decay process with a certain relaxation time. The last mechanism is the following: when a particle relaxes to 0 it will induce a random modification of the stress to all the other particles. This is modeled as a diffusion in stress space..

For a given probability distribution $p$, we recover the macroscopic stress $\tau$ as

$$ \int_{\sigma \in \mathbb{R}} \sigma p(t, \sigma) d\sigma. $$

For more details on the underlying physical ideas, we refer to [9] and [15].

## 3.2 Frame invariance

An extension to a multidimensional setting needs to be consistent with frame invariance. Frame invariance is the statement that when working in different frames one should observe the same behaviour of the material. Specifically, a rotation of the medium given by a matrix $Q$ should not change the stress except for rotating its principal axes with the medium, i.e. the stress should change to

$$ T^* = Q T Q^T. \qquad (31) $$

For this reason, we cannot simply replace Hooke's linear law with

$$ \partial_t T + (u \cdot \nabla) T = G_0 \nabla u. \qquad (32) $$

However, J. OLDROYD showed in [13] that this can be rectified by adding non-linear terms. More precisely let us define

$$ g_a(T, M) = T W(M) - W(M)T - a(D(M)\tau + \tau D(M)), \qquad (33) $$
where $D(M)$ and $W(M)$ are the symmetric and skew-symmetric part of the
$3 \times 3$ matrix $M$ and $a$ is a parameter usually taken in $[-1, 1]$ for stability issues.
We can now define

$$
\mathcal{D}_a T = \partial_t T + u \cdot \nabla T + g_a(\nabla u, T) \quad (34)
$$

and prove that an admissible law is given by

$$
\mathcal{D}_a T = 2G_0 D(\nabla u). \tag{35}
$$

We shall adopt this law for the deterministic part of stress evolution.

3.3 Full tensorial modeling

The modeling of the stochastic terms is less straightforward. We are guided
by a desire to keep the equations as simple as possible. Thus, we shall replace
the yield criterion by a critical level of a stress magnitude $|T|$, although other
possibilities exist. Also, we make no attempt, for instance, to build positive
definiteness constraints into our admissible stress space, even though our deter-
ministic model of stress evolution is naturally associated with the condition that
$G_0 I + aT$ is positive definite. We note that, for instance, traditional dumbbell
models such as described by R.B. BIRD, O. HASSAGER, R.C. ARMSTRONG
and C.F. CURTISS [6], are given in terms of a probability distribution for an
orientation vector $R$ and then the stress is given in terms of the dyadic product
$RR^T$, which naturally enforces positive definiteness. On the other hand, the
materials for which the Hébraud-Lequeux model is intended are not naturally
modeled by dumbbells. Moreover, stress tensors violating positive definiteness
will occur with low probability if our dimensionless parameter $\lambda$ is large, which,
in any case, is the only situation for which we have definitive results.

To make sense of integrals, diffusion operators etc., we need to define a
Euclidean metric in stress space. We note that the quantity

$$
\sum_{i,j} T_{ij}^2
\quad (36)
$$

is frame invariant, and if we set

$$
T = \begin{pmatrix}
\Sigma_1 & \Sigma_4/\sqrt{2} & \Sigma_5/\sqrt{2} \\
\Sigma_4/\sqrt{2} & \Sigma_2 & \Sigma_6/\sqrt{2} \\
\Sigma_5/\sqrt{2} & \Sigma_6/\sqrt{2} & \Sigma_3
\end{pmatrix},
\quad (37)
$$

then this frame invariant measure is simply the usual Euclidean norm in Σ-space.
The deterministic evolution transforms to a law of the form

$$
\partial_t \Sigma = G(\Sigma, \nabla u). \tag{38}
$$

We thus obtain a frame indifferent model of the same general structure as
HL by the equation

$$
\begin{equation}
\begin{split}
\partial_t p(t, x, \Sigma) + (u \cdot \nabla_x) p + G(\Sigma, \nabla u) : \nabla_{\Sigma} p = \\
\qquad - \frac{\mathbf{1}_{|\Sigma| > \sigma_c}(\Sigma)}{T_0} p(t, x, \Sigma) + \Gamma(p)(t, x) \rho(\Sigma) + \alpha \Gamma(p)(t, x) \Delta_{\Sigma} p(t, x, \Sigma).
\end{split}
\tag{39}
\end{equation}
$$
Here $\rho$ will usually taken to be the delta function at the origin, but it could be a more general radially symmetric function of unit integral.

We shall throughout assume that the flow is incompressible, i.e. $u$ has zero divergence. It can then be shown that the normalization

$$ \int p(t, x, \Sigma) \, d\Sigma = 1 \qquad (40) $$

is preserved if we set

$$ \Gamma = \frac{1}{T_0} \int_{|\Sigma| > \Sigma_c} p \, d\Sigma. \qquad (41) $$

To rewrite the system in non dimensional variables, we set

$$ \Sigma' = \frac{\Sigma}{\sigma_c}, \qquad p' = \sigma_c^6 p, \qquad \Gamma' = \int_{|\Sigma'|>1} p'(\Sigma') d\Sigma', \qquad (42) $$

for mesoscopic variables and

$$ x' = \frac{x}{L}, \qquad t' = \frac{t}{T_0}, \qquad u' = \frac{T_0}{L} u \qquad (43) $$

for macroscopic variables. Here $L$ is a length scale which actually disappears in the final equation. Finally, we set

$$ G'(\Sigma', \nabla' u') = \frac{T_0}{\sigma_c} G(\Sigma, \nabla u). \qquad (44) $$

Now we introduce the dimensionless number

$$ \mu = \frac{\alpha}{\sigma_c^2}. \qquad (45) $$

We obtain the system (4)

$$ \partial_t p(t, x, \Sigma) + (u \cdot \nabla_x)p + G(\Sigma, \nabla u) : \nabla_{\Sigma} p = - \mathbf{1}_{|\Sigma|>1}(\Sigma)p(t, x, \Sigma) + \Gamma(p)(t, x)\rho(\Sigma) + \mu\Gamma(p)(t, x)\Delta_{\Sigma} p(t, x, \Sigma), $$

completed with the constraint

$$ \int_{\mathbb{R}^6} p \, d\Sigma = 1, \qquad (46) $$

the definition of the fluidity

$$ \Gamma(p)(t, x) = \int_{|\Sigma|>1} p(t, x, \Sigma) d\Sigma, $$

and some initial condition.

# 4 Well-posedness of the model for steady homogeneous flow

This section is devoted to the proof of Propositions 1 and 2 and Theorem 1.
## 4.1 Proof of Propositions 1 and 2

We first address the solvability of the problem when a and b are positive and ignore questions about the limit when a or b tends to zero. In this case, it is most convenient to use the original form of the equation

$$ -\mu\Gamma\Delta p + \varepsilon Lp + h(|\Sigma|)p = \Gamma\rho. \quad (47) $$

Let $\hat{q}$ be the solution of the equation

$$ -\mu\Gamma\Delta\hat{q} + \varepsilon L\hat{q} = \Gamma\rho \quad (48) $$

inside B, with Dirichlet conditions on $\partial B$. The function $\hat{q}$ is in $W^{1,\beta}(B)$ and is $C^\infty$ outside the singular support of $\rho$. We can thus smoothly extend $\hat{q}$ outside B such that the extended function has compact support. Let this extended function be denoted by $\hat{p}$ and let $P = p - \hat{p}$. Then P satisfies the equation

$$ \mathcal{L}P := -\mu\Gamma\Delta P + \varepsilon LP + h(|\Sigma|)P = h(|\Sigma|)(\mu\Gamma\Delta\hat{p} - \varepsilon L\hat{p} - \hat{p}) := \mathcal{R}. \quad (49) $$

By construction, the right hand side $\mathcal{R}$ of this equation is $C^\infty$ outside B with a jump across $\partial B$; in particular it is in $H^{-1}(\mathbb{R}^6)$.

If we formally multiply (49) by P and integrate by parts, we obtain

$$ \mu\Gamma \int_{\mathbb{R}^6} |\nabla P|^2 d\Sigma + \int_{\mathbb{R}^6 \setminus B} P^2 d\Sigma = \int_{\mathbb{R}^6} \mathcal{R}P d\Sigma, \quad (50) $$

which yields the estimate

$$ \|P\|_{H^1(\mathbb{R}^6)} \le C \|\mathcal{R}\|_{H^{-1}(\mathbb{R}^6)}. \quad (51) $$

The only nonstandard feature in the justification of this estimate is that the operator L has unbounded coefficients. To get around this difficulty, we temporarily replace L by an operator of the form $V(\Sigma) \cdot \nabla_\Sigma$, where $V(\Sigma)$ is divergence free and has compact support. We note that once we know (51), we can immediately sharpen it to

$$ \|P\|_{H^1} + \|LP\|_{H^{-1}} \le C \|\mathcal{R}\|_{H^{-1}}. \quad (52) $$

Next, we consider the behavior at infinity. From (49), we find

$$ \mathcal{L}(P|\Sigma|^n) = \mathcal{R}|\Sigma|^n + \tilde{\mathcal{R}}, \quad (53) $$

where $\tilde{\mathcal{R}}$ has a bound of the form

$$ P(n^2|\Sigma|^{n-2} + n\varepsilon|V(\Sigma)||\Sigma|^{n-1}) + n|\nabla P||\Sigma|^{n-1}. \quad (54) $$

Note that $\mathcal{R}$ has compact support. As long as V has a bound of the form $|V(\Sigma)| \le C|\Sigma|$ and $\varepsilon$ is small relative to $1/n$, we can repeat the estimate above and obtain a bound on the $H^1$ norm of $P|\Sigma|^n$. We now take a sequence $V_k$ which converges to $A\Sigma + \lambda A_0$. All the estimates are uniform in $k$ and hence persist in the limit. In particular, if $\varepsilon$ is small enough, the solutions will be uniformly integrable, which means we still have

$$ \int_{\mathbb{R}^6 \setminus B} p d\Sigma = \Gamma \quad (55) $$
in the limit.

Elliptic regularity (see S. AGMON, A. DOUGLIS and L. NIRENBERG [2, 3]) yields the statement that the solutions are of class $C^\infty$ on both sides of $\partial B$ and outside the singular support of $\rho$. In particular, regularity near the interface $\partial B$ can be obtained by mapping both sides to a half plane and then regarding the equations on both sides with the interface conditions as an elliptic system [3]. Moreover, the analytic dependence on $\varepsilon$ and $\Gamma$ is a straightforward consequence of (52) and standard perturbation theory.

We further note that all the estimates discussed so far remain valid uniformly as $\varepsilon \to 0$ while $\Gamma$ remains positive, with the only exception of (52). It is the failure of (51) in the limit which leads to a potential loss of analyticity for $\varepsilon = 0$. However, we can take derivatives of (47) with respect to $\varepsilon$ and repeat the arguments above. If we take one derivative, for instance, we find

$$ \mathcal{L}\partial_{\varepsilon}p = -Lp + \partial_{\varepsilon}\mathcal{R}. \qquad (56) $$

This problem is of the same form as the equation for $p$ itself, and we can repeat the same estimates. We can repeat the procedure for higher derivatives. As a consequence, we find bounds on all derivatives of $p$ which remain valid as $\varepsilon \to 0$.

On the other hand, ellipticity is lost for $\Gamma \to 0$, and we need a different argument. This is why we transformed our variables to get the system (20), which remains nondegenerate in the radial direction in the limit. We now show how (20) can be used to obtain estimates which are uniform as $b \to 0$.

In an analogous fashion as above, we subtract from $\tilde{q}$ a reference function which satisfies $-\Delta\hat{q} + aL\hat{q} = -aLE$, with Dirichlet boundary conditions, and from $\tilde{r}$ a smooth reference function which has compact support and satisfies the interface conditions. We then end up with a problem of the form

$$ \left\{ \begin{array}{l} -\Delta Q + aLQ = 0, \\ -\partial_{\theta}^2 R - \frac{b^2}{(1+b\theta)^2} \Delta_{\tan} R + ab^2((A\eta) \cdot \nabla_{\tan})R \\ \qquad + ab(1+b\theta)(A\eta \cdot \eta)\partial_{\theta}R + R + SR = f, \\ Q(\eta) = R(0, \eta), \\ b\partial_n Q(\eta) - \partial_{\theta}R(0, \eta) = 0. \end{array} \right. \qquad (57) $$

Here $f$ is a smooth function of compact support which depends smoothly on $a$ and $b$, and $S$ denotes a linear operator which is easily dealt with as a perturbation in the estimates.

We multiply the first equation in (57) by $bQ$, the second equation by $R$, and integrate. After an integration by parts, this yields an estimate of the form

$$ \sqrt{b}\|Q\|_{H^1(B)} + \|R\|_{L^2(\mathbb{R}\times\partial B)} + \|\partial_{\theta}R\|_{L^2(\mathbb{R}\times\partial B)} \leq C. \qquad (58) $$

Here and in the following $C$ is a generic constant which does not depend on $a$ and $b$. Now take $\chi$ to be any smooth radial function on $B$ which is 1 in a neighborhood $U$ of $\partial B$ and 0 in a neighborhood of the origin. Moreover, let $\partial_t$ be any derivative tangent to $\partial B$. We can then apply $\partial_t$ to all the equations of (57), multiply the first equation by $\chi b \partial_t Q$, the second equation by $\partial_t R$ and repeat the same steps that led to (58) to obtain

$$ \sqrt{b}\|\partial_t Q\|_{H^1(U)} + \|\partial_t R\|_{L^2(\mathbb{R}\times\partial B)} + \|\partial_t \partial_{\theta} R\|_{L^2(\mathbb{R}\times\partial B)} \leq C. \qquad (59) $$
This yields a bound for $||R||_{H^1(\partial B)}$. Using the first equation of (57) and the first of the interface conditions, we get a bound for $||Q||_{H^{3/2}(B)}$. It is now easy to see how to bootstrap this argument to get estimates for higher derivatives.

For the behavior at infinity, multiply the second equation of (57) by $\theta^n R$ and integrate by parts. This yields a bound of the form

$$ \|\theta^{n/2} R\|_{L^2(\mathbb{R}\times\partial B)} \leq C. \qquad (60) $$

Here $n$ can be taken arbitrarily large as $ab^2 \to 0$.

We now have estimates for $Q$ and $R$ which hold uniformly in $a$ and $b$. Now we make take derivatives of (57) with respect to $a$ and $b$ and observe that the derivatives of $Q$ and $R$ satisfy elliptic systems of the same form, with inhomogeneous terms that depend on $Q$ and $R$. We can thus repeat the same estimates to get bounds for the derivatives. The same procedure can be repeated for higher order derivatives.

## 4.2 Study of the function F when a or b become zero

### 4.2.1 Study of the function $b \mapsto F(0, b)$

When we set $a = 0$ in (20) we get the system

$$ \left\{ \begin{array}{l} -\Delta \tilde{q} = 0 \\ -\partial_{\theta}^2 \tilde{r} - \frac{b^2}{(1+b\theta)^2} \Delta_{\tan} \tilde{r} - \frac{5b}{1+b\theta} \partial_{\theta} \tilde{r} + \tilde{r} = 0 \\ \tilde{q}(\eta) = \tilde{r}(0, \eta) \\ b\partial_n \tilde{q}(\eta) - \partial_{\theta} \tilde{r}(0, \eta) = -b\partial_n E(\eta) \end{array} \right. \qquad (61) $$

which has the property of being “rotation invariant” in the sense that the solution does not depend on $\eta$. Indeed $\rho$ being itself invariant by rotation, $E$ is a radial function and $\partial_n E(\eta)$ is actually a constant, denoted by $E'(1)$ in the sequel. Precisely we can compute

$$ \tilde{q}^{0, b} = \frac{bE'(1)}{B'(1/b)} B(1/b), \qquad (62) $$

$$ \tilde{r}^{0, b} = \frac{bE'(1)}{B'(1/b)} B\left(\frac{1}{b} + \theta\right). \qquad (63) $$

where $B(x) = K_2(x)/x^2$, and $K_2$ is the modified Bessel function of the second kind and order 2. Then plugging these expressions into (21) yields,

$$ F(0, b) = \mu_c - \frac{1}{6} \frac{b}{B'(1/b)} B(1/b) + b^2. \qquad (64) $$

We note that

$$ -\frac{B(\omega)}{B'(\omega)} \xrightarrow{\omega \to +\infty} 1, \qquad (65) $$

by classical estimates on the modified Bessel function of the second kind (see M. ABRAMOWITZ and I.A. STEGUN [1]), and hence

$$ \lim_{b \to 0^+} F(0, b) = \mu_c. \qquad (66) $$
We shall prove that the function $\mathfrak{b} \mapsto -\mathfrak{b}\mathcal{B}(1/\mathfrak{b})/\mathcal{B}'(1/\mathfrak{b})$ is monotone increasing.
We immediately conclude that

$$F(0, \mathfrak{b}) \geq \mu_c + \mathfrak{b}^2, \qquad (67)$$

and thus

$$\lim_{\mathfrak{b} \to +\infty} F(0, \mathfrak{b}) = +\infty. \qquad (68)$$

It remains to show the monotonicity of the function $\mathfrak{b} \mapsto -\mathfrak{b}\mathcal{B}(1/\mathfrak{b})/\mathcal{B}'(1/\mathfrak{b})$.
Obviously, this equivalent to showing that the function

$$\phi(x) = - \frac{x\mathcal{B}'(x)}{\mathcal{B}(x)} = \frac{xK_3(x)}{K_2(x)} \qquad (69)$$

is monotone increasing. We find

$$\phi'(x) = \frac{x^2(K_1(x)^2 + K_3(x)^2) - (16 + 2x^2)K_2(x)^2}{2xK_2(x)^2}.$$

The positivity of this function follows from the recurrence relation (again see [1])

$$K_3(x) = K_1(x) + \frac{4}{x}K_2(x),$$

and the Turan inequality (see M.E.H ISMAIL and M.E. MULDOON [10])

$$K_2(x)^2 < K_1(x)K_3(x) = K_1(x)^2 + \frac{4}{x}K_1(x)K_2(x).$$

### 4.2.2 Study of the function $\mathfrak{a} \mapsto F(\mathfrak{a}, 0)$

We compute $F(\mathfrak{a}, 0)$ via (21) to obtain:

$$F(\mathfrak{a}, 0) = \mu_c + \int_{|\Sigma|>1} \tilde{q}^{\mathfrak{a},0}(\Sigma)d\Sigma \qquad (70)$$

and Hypothesis **H2** exactly states that this function has a strictly negative derivative with respect to $\mathfrak{a}$. When $\mathfrak{a}$ vanishes we obtain

$$\begin{cases} -\Delta \tilde{q}^{0,0} = 0, \\ \tilde{q}^{0,0}(\eta) = 0, \end{cases} \qquad (71)$$

thus $\tilde{q}^{0,0} = 0$ and

$$F(0, 0) = \mu_c. \qquad (72)$$

On the other hand we will show that

$$\lim_{\mathfrak{a} \to +\infty} \tilde{q}^{\mathfrak{a},0} = -E, \qquad (73)$$

which proves that

$$\lim_{\mathfrak{a} \to +\infty} F(\mathfrak{a}, 0) = 0. \qquad (74)$$

What we are studying is actually a singular limit of a diffusion-transport equation to a stationary transport equation for vanishing viscosity. We require the following lemma:
**Lemma 1.** There is a function $\zeta$ of class $C^2$, such that we have the following estimates:

$$m_1 \le \zeta \le m_2, \quad -(\mathcal{A}\Sigma + \lambda\mathcal{A}^0) \cdot \nabla\zeta \ge m_3, \quad |\nabla\zeta| \le m_4, \quad |\Delta\zeta| \le m_5,$$

where $m_i$ are positive constants.

The proof of Lemma 1 requires that all trajectories of the vectorfield $\Sigma \mapsto \mathcal{A}\Sigma + \mathcal{A}_0$ which start inside the unit ball leave the ball in both directions. In this case, we can simply initialize $m$ on the part of the unit ball where these trajectories exit. Let $\zeta$ equal any smooth positive function at these points. Inside the ball, we then set $m(\zeta) = m(\zeta^+) + l(\zeta)$. Here $\zeta^+$ is the point where the trajectory leaves the ball and $l(\zeta)$ is the length of the trajectory between $\zeta$ and $\zeta^+$. See D. BRESCH and J. SIMON [7].

Now, recall that $\mathcal{A}$ depends on $a$ and $M$. Recall also that $\mathcal{A}^0$ depends only on $M$. The following proposition gives a general result:

**Proposition 4.** For every deformation type $M$ with $M + M^T \ne 0$, and for every $a$, there exists $\lambda_m \ge 0$, so that, for $\lambda \ge \lambda_m$, the characteristic curves of the vector field $\Sigma \mapsto \mathcal{A}\Sigma + \lambda\mathcal{A}^0$ cross the unit sphere.

This proposition is essentially noting that the vector field depends very simply on $\lambda$: when $\lambda$ is large, the vector field is in the limit a constant vector field and characteristic curves are parallel lines.

Now we apply the lemma to prove (73). We consider a sequence of approximate problems

$$-\frac{1}{a}\Delta H_a^j + (\mathcal{A}\Sigma + \lambda\mathcal{A}_0) \cdot \nabla H_a^j = \frac{1}{a}\rho^j, \quad (75)$$

with Dirichlet boundary conditions, and nonnegative smooth $\rho^j$. By the maximum principle, $H_a^j$ is also nonnegative. Let $\zeta$ be a smooth function verifying the estimates of Lemma 1. Multiply the equation by $\zeta$ and integrate. We find, after integrating by parts,

$$-\frac{1}{a} \int_B (\Delta \zeta) H_a^j - \frac{1}{a} \int_{\partial B} \zeta \partial_n H_a^j - \int_B ((\mathcal{A}\Sigma + \lambda \mathcal{A}_0) \cdot \nabla \zeta) H_a^j = \frac{1}{a} \int_B \rho^j \zeta. \quad (76)$$

On the left side, the first term is bounded from below by $-m_5/a \int_B H_a^j$, the second term is nonnegative by the maximum principle, and the third term is bounded below by $m_3 \int_B H_a^j$. The right hand side is bounded by $m_2/a \int_B \rho^j$. We conclude that

$$\int_B H_a^j \le \frac{m_2}{m_3 a - m_5} \int_B \rho^j. \quad (77)$$

This is true for any nonnegative function $\rho^j$. Now let $\rho^j \to \rho$ when $j \to +\infty$. What we need is the convergence of both sides of (77) as $j \to +\infty$. The right-hand side can be written

$$\int_B \rho^j = \langle \rho^j, \mathbf{1} \rangle_{\mathcal{M}, \mathcal{C}^0(\bar{B})}, \quad (78)$$

and as such converges to $\langle \rho, \mathbf{1} \rangle_{\mathcal{M}, \mathcal{C}^0(\bar{B})}$. On the other hand we actually want to know if $\int_B H_a^j \to \int_B (\tilde{q}^{a,0} + E)$. We note first that since $H_a^j$ is positive, (77) gives us a uniform bound in $L^1(B)$ when $a$ is fixed. Consequently any
j-subsequence of ($H_a^j$) weakly converges to a positive measure $\chi_a$. Moreover we have that

$$ \chi_a(B) \le \frac{m_2}{m_3 a - C_1}, \qquad (79) $$

so that $\chi_a(B)$ tends to 0 as $a \to +\infty$.

We now use the fact that when $a$ is fixed, the linear operator $(-\Delta + aL)^{-1}$ is continuous over $W^{-1,\beta}$ with value in $W_0^{1,\beta}$. Consequently $H_a^j$ converges to $\tilde{q}^{a,0} + E$ at least strongly in $L^\beta$ and thus in $L^1(B)$ when $j \to +\infty$. Which means that for any j-subsequence $\int_B H_a^j$ converges to $\int_B (\tilde{q}^{a,0} + E)$ and thus

$$ \int_B (\tilde{q}^{a,0} + E) = \chi_a(B) \to_{a \to +\infty} 0. \qquad (80) $$

This is what we wanted to prove.

### 4.2.3 Study of the case $\mu = \mu_c$

From (72) we get that

$$ F(0,0) = \mu_c. \qquad (81) $$

We are left with the task of proving

$$ \partial_b F(0,0) > 0, \qquad (82) $$

$$ \partial_a F(0,0) = 0, \qquad (83) $$

$$ \partial_{aa}^2 F(0,0) < 0. \qquad (84) $$

**Computation of $\partial_b F(0,0)$**. We use (64) to write:

$$ \frac{F(0, b) - F(0, 0)}{b} = -\frac{1}{6} \frac{1}{\mathcal{B}'(1/b)} \mathcal{B}(1/b) + b, \qquad (86) $$

and by the limit (65) we get

$$ \partial_b F(0,0) = \frac{1}{6}. \qquad (87) $$

**Computation of $\partial_a F(0,0)$**. We differentiate $(\tilde{q}^{a,0}, \tilde{r}^{a,0})$ with respect to $a$ and set $a = 0$. We obtain the relations (recall that $\tilde{q}^{0,0} = 0$):

$$ \left\{ \begin{array}{ll} -\Delta \partial_a \tilde{q}^{0,0} = -LE & \text{for } |\Sigma| < 1, \\ \partial_a \tilde{q}^{0,0}(\eta) = 0 & \text{for } |\eta| = 1, \\ \partial_a \tilde{r}^{a,0}(\theta, \eta) = 0 & \text{for } \theta > 0 \text{ and } |\eta| = 1, \end{array} \right. \qquad (88) $$

and thus by (21),

$$ \partial_a F(0,0) = \int_B \partial_a \tilde{q}^{0,0}(\Sigma) d\Sigma. \qquad (89) $$

To compute this integral let us introduce the test function

$$ \psi_0 = \frac{1}{12}(1 - |\Sigma|^2), \qquad (90) $$
This function has the property that

$$ \begin{cases} -\Delta\psi_0 = 1 & \text{for } |\Sigma| \le 1, \\ \psi_0(\eta) = 0 & \text{for } |\eta| = 1. \end{cases} \tag{91} $$

Consequently multiplying the equation of $\partial_a \tilde{q}^{0,0}$ from (88) by $\psi_0$ and integrating by parts yields:

$$ \begin{aligned} \int_B \partial_a \tilde{q}^{0,0}(\Sigma) d\Sigma &= - \int_B L E(\Sigma) \psi_0(\Sigma) d\Sigma \\ &= - \int_{r=0}^1 \frac{1-r^2}{12} E'(r) r^5 \left( \int_{|\eta|=1} (A\eta r + \lambda A^0) \cdot \eta d\eta \right) dr, \end{aligned} \tag{92} $$

using that by Hypothesis **H1**, $E$ is radial. We now use the fact that $A$ is traceless so that

$$ \int_{|\eta|=1} (A\eta r + \lambda A^0) \cdot \eta d\eta = 0, \tag{93} $$

which ends the proof.

**Computation of** $\partial_{aa}^2 F(0, 0)$. We differentiate twice $(\tilde{q}^{a,0}, \tilde{r}^{a,0})$ with respect to $a$ and set $a=0$. We obtain the relations (recall that $\tilde{q}^{0,0} = 0$):

$$ \left\{ \begin{aligned} -\Delta \partial_{aa}^2 \tilde{q}^{0,0} &= -L \partial_a \tilde{q}^{0,0} & \text{for } |\Sigma| < 1, \\ \partial_{aa}^2 \tilde{q}^{0,0}(\eta) &= 0 & \text{for } |\eta| = 1, \\ \partial_{aa}^2 \tilde{r}^{a,0}(\theta, \eta) &= 0 & \text{for } \theta > 0 \text{ and } |\eta| = 1. \end{aligned} \right. \tag{94} $$

Note that from (21) and taking into account that $\partial_{aa}^2 \tilde{r} = 0$ we have

$$ \partial_{aa}^2 F(0, 0) = \int_B \partial_{aa}^2 \tilde{q}^{0,0}(\Sigma) d\Sigma, \tag{95} $$

so we have to prove that this integral is negative. Again we use the test function $\psi_0$ defined by (90) to write

$$ \begin{aligned} \int_B \partial_{aa}^2 \tilde{q}^{0,0}(\Sigma) d\Sigma &= - \int_B L \partial_a \tilde{q}^{0,0}(\Sigma) \psi_0(\Sigma) d\Sigma \\ &= - \int_B (A\Sigma + \lambda A^0) \cdot \nabla \partial_a \tilde{q}^{0,0}(\Sigma) \psi_0(\Sigma) d\Sigma \end{aligned} \tag{96} $$

We decompose $A = A^s + A^a$ where $A^s$ is the symmetric part of $A$ and $A^a$ its skew-symmetric part. This decomposition is somewhat reminiscent of the decomposition of $\nabla u$ into its symmetric and skew-symmetric part in the study of the well posedness of some stationary Fokker-Planck equation arising in the study of polymeric fluid that is done by JOURDAIN, LE BRIS, LELIÈVRE and OTTO [12] or ARNOLD, CARILLO and MANZINI [4]. We can then write

$$ A^s \Sigma + \lambda A^0 = \nabla \left( \frac{1}{2} A^s \Sigma \cdot \Sigma + \lambda A^0 \cdot \Sigma \right), \tag{97} $$
and integrate by parts to find:

$$
\begin{align}
& - \int_B (\mathcal{A}\Sigma + \lambda\mathcal{A}^0) \cdot \nabla\partial_a \tilde{q}^{0,0}(\Sigma)\psi_0(\Sigma)d\Sigma \nonumber \\
&= - \int_B \nabla\left(\frac{1}{2}\mathcal{A}^s\Sigma \cdot \Sigma + \lambda\mathcal{A}^0 \cdot \Sigma\right) \cdot \nabla\partial_a \tilde{q}^{0,0}(\Sigma)\psi_0(\Sigma)d\Sigma \nonumber \\
&\quad - \int (\mathcal{A}^a\Sigma) \cdot \nabla\partial_a \tilde{q}^{0,0}(\Sigma)\psi_0(\Sigma)d\Sigma \nonumber \\
&= \int_B \left(\frac{1}{2}\mathcal{A}^s\Sigma \cdot \Sigma + \lambda\mathcal{A}^0 \cdot \Sigma\right) \Delta\partial_a \tilde{q}^{0,0}\psi_0(\Sigma)d\Sigma \nonumber \\
&\quad + \int_B \left(\frac{1}{2}\mathcal{A}^s\Sigma \cdot \Sigma + \lambda\mathcal{A}^0 \cdot \Sigma\right) \nabla\partial_a \tilde{q}^{0,0} \cdot \nabla\psi_0(\Sigma)d\Sigma \tag{98} \\
&\quad - \int_B (\mathcal{A}^a\Sigma) \cdot \nabla\partial_a \tilde{q}^{0,0}(\Sigma)\psi_0(\Sigma)d\Sigma \nonumber \\
&= \int_B \left(\frac{1}{2}\mathcal{A}^s\Sigma \cdot \Sigma + \lambda\mathcal{A}^0 \cdot \Sigma\right) \Delta\partial_a \tilde{q}^{0,0}\psi_0(\Sigma)d\Sigma \nonumber \\
&\quad - \int_B \partial_a \tilde{q}^{0,0} \operatorname{div} \left(\left(\frac{1}{2}\mathcal{A}^s\Sigma \cdot \Sigma + \lambda\mathcal{A}^0 \cdot \Sigma\right) \nabla\psi_0(\Sigma)\right) d\Sigma \nonumber \\
&\quad - \int_B (\mathcal{A}^a\Sigma) \cdot \nabla\partial_a \tilde{q}^{0,0}(\Sigma)\psi_0(\Sigma)d\Sigma \nonumber \\
&= \mathcal{I}_1 - \mathcal{I}_2 - \mathcal{I}_3 \nonumber
\end{align}
$$

We now compute each of the $\mathcal{I}_k$ separately.

Computation of $\mathcal{I}_3$. We first note that

$$
\begin{align}
\mathcal{I}_3 &= - \int_B \psi^0(\Sigma) \Sigma \cdot \mathcal{A}^a \nabla \partial_a \tilde{q}^{0,0}(\Sigma) d\Sigma \nonumber \\
&= - \int_B \nabla \left( \int_B^{|\Sigma|} r \psi^0(r) dr \right) \cdot \mathcal{A}^a \nabla \partial_a \tilde{q}^{0,0}(\Sigma) d\Sigma \tag{99} \\
&= \int_B \left( \int_B^{|\Sigma|} r \psi^0(r) dr \right) \operatorname{div} (\mathcal{A}^a \nabla \partial_a \tilde{q}^{0,0}(\Sigma)) d\Sigma \nonumber \\
&= 0, \nonumber
\end{align}
$$

because $\operatorname{div}(\mathcal{A}^a \nabla \partial_a \tilde{q}^{0,0}) = \mathcal{A}^a : d^2 \partial_a \tilde{q}^{0,0}$ and this contracted product is 0 since $d^2 \partial_a \tilde{q}^{0,0}$ is a symmetric matrix and $\mathcal{A}^a$ a skew-symmetric one. We also note that the boundary term in the last integration by parts vanishes, since $\nabla \partial_a \tilde{q}^{0,0}$ is normal to the boundary, and then $\mathcal{A}^a$ maps this vector to a vector which is tangent to the boundary.

Computation of $\mathcal{I}_1$. We use the equation obeyed by $\partial_a \tilde{q}^{0,0}$ from (88) to write

$$
\begin{align}
\mathcal{I}_1 &= \int_B \left( \frac{1}{2} \mathcal{A}^s \Sigma \cdot \Sigma + \lambda \mathcal{A}^0 \cdot \Sigma \right) (\mathcal{A} \Sigma + \lambda \mathcal{A}^0) \cdot \nabla E(\Sigma) \psi_0(\Sigma) d\Sigma \tag{100} \\
&= \int_B \left( \frac{1}{2} \mathcal{A}^s \Sigma \cdot \Sigma + \mathcal{A}^0 \cdot \Sigma \right) (\mathcal{A} \Sigma + \lambda \mathcal{A}^0) \cdot \left( E'(|\Sigma|) \frac{\Sigma}{|\Sigma|} \right) \psi_0(\Sigma) d\Sigma \nonumber
\end{align}
$$
We can simplify this last integral by using the change of variable $\Sigma \to -\Sigma$ to remove the odd part of the function. This leads to:

$$ \begin{aligned} & \int_B \left( \frac{1}{2} \mathcal{A}^s \Sigma \cdot \Sigma + \lambda \mathcal{A}^0 \cdot \Sigma \right) (\mathcal{A}\Sigma + \lambda \mathcal{A}^0) \cdot \left( E'(|\Sigma|) \frac{\Sigma}{|\Sigma|} \right) \psi_0(\Sigma) d\Sigma \\ = & \frac{1}{2} \int_B \frac{\psi_0(\Sigma) E'(|\Sigma|)}{|\Sigma|} (\mathcal{A}^s \Sigma \cdot \Sigma)^2 d\Sigma \\ & + \lambda^2 \int_B \frac{\psi_0(\Sigma) E'(|\Sigma|)}{|\Sigma|} (\mathcal{A}^0 \cdot \Sigma)^2 d\Sigma \end{aligned} \quad (101) $$

We note that since $\mathcal{A}^s$ is symmetric we can use an orthogonal matrix to change variables and write:

$$ \begin{aligned} & \frac{1}{2} \int_B \frac{\psi_0(\Sigma) E'(|\Sigma|)}{|\Sigma|} (\mathcal{A}^s \Sigma \cdot \Sigma)^2 d\Sigma \\ = & \frac{1}{2} \int_B \frac{\psi_0(\Sigma) E'(|\Sigma|)}{|\Sigma|} \left( \sum_i \lambda_i \Sigma_i^2 \right)^2 d\Sigma \quad (102) \\ = & \frac{1}{2} \sum_{i,j} \lambda_i \lambda_j \int_B \frac{\psi_0(\Sigma) E'(|\Sigma|)}{|\Sigma|} \Sigma_i^2 \Sigma_j^2 d\Sigma, \end{aligned} $$

where the $\lambda_i$ are the eigenvalues of $\mathcal{A}^s$, the symmetric part of $\mathcal{A}$. But it is clear that

$$ \int_B \frac{\psi_0(\Sigma) E'(|\Sigma|)}{|\Sigma|} \Sigma_i^2 \Sigma_j^2 d\Sigma $$

does not depend on $(i, j)$ when $i \neq j$ and does not depend on $i$ when $i = j$. Moreover, because $\mathcal{A}^s$ is traceless when $\mathcal{A}$ is, we have $\sum_{i \neq j} \lambda_i \lambda_j = (\sum_i \lambda_i)^2 - \sum_i \lambda_i^2 = -\sum_i \lambda_i^2$, and thus

$$ \begin{aligned} & \sum_{i,j} \lambda_i \lambda_j \frac{1}{2} \int_B \frac{\psi_0(\Sigma) E'(|\Sigma|)}{|\Sigma|} \Sigma_i^2 \Sigma_j^2 d\Sigma \\ = & \frac{1}{4} \left( \sum_i \lambda_i^2 \right) \int_B \frac{\psi_0(\Sigma) E'(|\Sigma|)}{|\Sigma|} (\Sigma_1^4 + \Sigma_2^4 - 2\Sigma_1^2\Sigma_2^2) d\Sigma \quad (103) \end{aligned} $$

Note that by Maximum principle, $E'$ is negative. Consequently, we get

$$ \frac{1}{2} \int_B \frac{\psi_0(\Sigma) E'(|\Sigma|)}{|\Sigma|} (\mathcal{A}^s \Sigma \cdot \Sigma)^2 d\Sigma = -C_1 \left( \sum_i \lambda_i^2 \right). \quad (104) $$

where

$$ \begin{aligned} C_1 = & -\frac{1}{4} \int_B \frac{\psi_0(\Sigma) E'(|\Sigma|)}{|\Sigma|} (\Sigma_1^4 + \Sigma_2^4 - 2\Sigma_1^2\Sigma_2^2) d\Sigma \\ = & -\frac{1}{4} \int_B \frac{\psi_0(\Sigma) E'(|\Sigma|)}{|\Sigma|} (\Sigma_1^2 - \Sigma_2^2)^2 d\Sigma \end{aligned} \quad (105) $$

is a positive constant depending only on $\rho$. In the same way, we can change variables so that $\mathcal{A}^0 = |\mathcal{A}_0|(1, 0, \dots, 0)$. We would obtain that there is a positive constant $C_2$ with expression
$$C_2 = \int_B \frac{\psi_0(\Sigma) |E'(|\Sigma|)|}{|\Sigma|} \Sigma_1^2 d\Sigma. \qquad (106)$$

such that

$$\mathcal{I}_1 = - \left( C_1 \sum_i \lambda_i^2 + C_2 \lambda^2 |\mathcal{A}^0|^2 \right). \qquad (107)$$

Computation of $\mathcal{I}_2$. Recall that we have

$$\mathcal{I}_2 = \int_B \partial_a \tilde{q}^{0,0} \mathrm{div} \left( \left( \frac{1}{2} \mathcal{A}^s \Sigma \cdot \Sigma + \lambda \mathcal{A}^0 \cdot \Sigma \right) \nabla \psi_0(\Sigma) \right) d\Sigma. \qquad (108)$$

With some elementary differential calculus we find

$$\mathrm{div} \left( \left( \frac{1}{2} \mathcal{A}^s \Sigma \cdot \Sigma + \lambda \mathcal{A}^0 \cdot \Sigma \right) \nabla \psi_0(\Sigma) \right) = - \left( \frac{2}{3} \mathcal{A}^s \Sigma \cdot \Sigma + \frac{7}{6} \lambda \mathcal{A}^0 \cdot \Sigma \right) \qquad (109)$$

and thus:

$$\mathcal{I}_2 = - \int_B \left( \frac{2}{3} \mathcal{A}^s \Sigma \cdot \Sigma + \frac{7}{6} \mathcal{A}^0 \cdot \Sigma \right) \partial_a \tilde{q}^{0,0} d\Sigma. \qquad (110)$$

We need to decompose into even and odd parts (with respect to $\Sigma$) which means we have to identify the even and odd parts of $\partial_a \tilde{q}^{0,0}$.

First note that we have from (88) that

$$\begin{aligned} -\Delta \partial_a \tilde{q}^{0,0} &= -LE \\ &= -(\mathcal{A}\Sigma + \lambda\mathcal{A}^0) \cdot \frac{E'(|\Sigma|)\Sigma}{|\Sigma|} \qquad (111) \\ &= -\frac{E'(|\Sigma|)}{|\Sigma|}(\mathcal{A}^s\Sigma \cdot \Sigma) - \lambda \frac{E'(|\Sigma|)}{|\Sigma|}\mathcal{A}^0 \cdot \Sigma. \end{aligned}$$

The even part of $\partial_a \tilde{q}^{0,0}$ is thus given by $G_1$ and the odd part by $G_2$, where

$$\left\{ \begin{array}{ll} -\Delta G_1 = -E'(|\Sigma|) \frac{A^s \Sigma \cdot \Sigma}{|\Sigma|} & \text{in } |\Sigma| < 1, \\ G_1 = 0 & \text{on } |\Sigma| = 1, \end{array} \right. \qquad (112)$$

and

$$\left\{ \begin{array}{ll} -\Delta G_2 = -\lambda \frac{E'(|\Sigma|)}{|\Sigma|} \mathcal{A}^0 \cdot \Sigma & \text{in } |\Sigma| < 1, \\ G_2 = 0 & \text{on } |\Sigma| = 1, \end{array} \right. \qquad (113)$$

With $G_1$ and $G_2$ we can split $\mathcal{I}_2$ in two terms:

$$\begin{aligned} \mathcal{I}_2 &= - \int_B \left( \frac{2}{3} \mathcal{A}^s \Sigma \cdot \Sigma \right) G_1(\Sigma) d\Sigma - \frac{7}{6} \int_B (\mathcal{A}^0 \cdot \Sigma) G_2(\Sigma) d\Sigma \qquad (114) \\ &= \mathcal{J}_1 + \mathcal{J}_2. \end{aligned}$$

Computation of $\mathcal{J}_1$. We now introduce an orthogonal matrix $Q$ such that $Q^T A^s Q$ is a diagonal matrix. We change variables by setting $\Sigma = QS$ and
consider the function $G_1^*$ such that $G_1^*(S) = G_1(\Sigma)$. The function $G_1^*$ thus satisfies the following equation:

$$ -\Delta G_1^* = - \sum_i \lambda_i \frac{E'(|S|)}{|S|} S_i^2. \qquad (115) $$

Now we introduce $\overline{G}_1$ which satisfies:

$$ \left\{ \begin{array}{l} -\Delta \overline{G}_1 = - \frac{E'(\sqrt{y^2 + |y'|^2})}{\sqrt{y^2 + |y'|^2}} y^2, \\ \overline{G}_1 = 0, \end{array} \right. \qquad (116) $$

where $(y, y')$ belongs in the unit ball of $\mathbb{R}^6$. $\overline{G}_1$ can be seen as a function of $y$ and $|y'|$ or a vector function over the unit ball of $\mathbb{R} \times \mathbb{R}^5$. We can now write, by linearity,

$$ G_1^*(S_1, \dots, S_6) = \sum_i \lambda_i \overline{G}_1(S_i, |\hat{S}_i|), $$

where $\hat{S}_i$ denotes any vector of $\mathbb{R}^5$ with components $S_j$ with $j \neq i$. Using now the matrix $Q$ to change variable in the integral we can write:

$$ \begin{align} & \int_B \left( \frac{2}{3} A^s \Sigma \cdot \Sigma \right) G_1(\Sigma) d\Sigma = \frac{2}{3} \int_B \left( \sum_i \lambda_i S_i^2 \right) G_1^*(S) dS \nonumber \\ &= \frac{2}{3} \int_B \left( \sum_i \lambda_i S_i^2 \right) \left( \sum_j \lambda_j \overline{G}_1(S_j, |\hat{S}_j|) \right) dS \tag{117} \\ &= \frac{2}{3} \sum_{i,j} \lambda_i \lambda_j \int_B S_i^2 \overline{G}_1(S_j, |\hat{S}_j|) dS \nonumber \end{align} $$

Once again we see that $\int_{|S| \le 1} S_i^2 \overline{G}_1(S_j, |\hat{S}_j|) dS$ does not depend on $(i, j)$ if $i \neq j$ and on $i$ for the terms $i = j$. This proves that there is a numerical constant

$$ C_3 = \frac{2}{3} \int_{|S| \le 1} (S_1^2 - S_2^2) \overline{G}_1(S_1, \hat{S}_1) dS. \qquad (118) $$

so that

$$ \mathcal{J}_1 = -C_3 \sum_i \lambda_i^2. \qquad (119) $$

To evaluate this integral we introduce the following test function:

$$ \begin{align} \psi_1(y, y') &= -\frac{3}{64}y^4 - \frac{7}{160}y^2|y'|^2 + \frac{1}{24}y^2 + \frac{1}{320}|y'|^4 - \frac{1}{120}|y'|^2 + \frac{1}{192} \nonumber \\ &= (1 - y^2 - |y'|^2) \left( \frac{3}{64}y^2 - \frac{1}{320}|y'|^2 + \frac{1}{192} \right), \tag{120} \end{align} $$

where $y \in \mathbb{R}$ and $y' \in \mathbb{R}^5$. One can check that this function satisfies $-\Delta \psi_1 = y^2$ and obviously, $\psi_1$ vanishes on the unit sphere. Now we have the following relation:

$$ \begin{align} & \int_{|S| \le 1} (S_1^2 - S_2^2) \overline{G}_1(S_1, \hat{S}_1) dS \tag{121} \\ &= - \int_{|S| \le 1} \frac{E'(|S|) S_1^2}{\sqrt{S_1^2 + |\hat{S}_1|^2}} (\psi_1(S_1, \hat{S}_1) - \psi_1(S_2, \hat{S}_2)) dS. \tag{122} \end{align} $$
One can check using the definition of $\psi_1$ that

$$ \psi_1(S_1, \hat{S}_1) - \psi_1(S_2, \hat{S}_2) = \frac{1}{20}(S_1^2 - S_2^2)(1 - |S|^2), \quad (123) $$

so that what is left to compute is

$$ \int_{|S| \le 1} (S_1^2 - S_2^2) \overline{G}_1(S_1, \hat{S}_1) dS = - \int_{|S| \le 1} \frac{E'(|S|) S_1^2 (S_1^2 - S_2^2)(1 - |S|^2)}{20|S|} dS. \quad (124) $$

Now using the change of variable which exchanges the variables $S_1$ and $S_2$ and writing the average we have that the integral is also

$$ \int_{|S| \le 1} (S_1^2 - S_2^2) \overline{G}_1(S_1, \hat{S}_1) dS = - \int_{|S| \le 1} \frac{E'(|S|)(S_1^2 - S_2^2)^2 (1 - |S|^2)}{40|S|} dS. \quad (125) $$

Therefore

$$ C_3 = \frac{1}{5} \int_B \frac{\psi_0(\Sigma) |E'(|S|)|}{|S|} (S_1^2 - S_2^2)^2 dS. \quad (126) $$

By comparing this expression with the expression of $C_1$ from (105) we have $C_3/C_1 = \frac{4}{5} < 1$. This implies that $C_5 = C_1 - C_3 > 0$.

Computation of $\mathcal{J}_2$. The integral $\mathcal{J}_2$ is treated using the same method. We introduce an orthogonal matrix $Q^0$ such that $(Q^0)^T A^0$ is the vector $(|A^0|, 0, \dots, 0)$ and we write $S = Q^0\Sigma$ and $G_2^*(S) = G_2(\Sigma)$. We thus have by a change of variable that:

$$ \mathcal{J}_2 = -\frac{7}{6} \int_{|\Sigma| \le 1} (\lambda A^0 \cdot \Sigma) G_2(\Sigma) d\Sigma = -\frac{7\lambda|A^0|}{6} \int_{|S| \le 1} S_1 G_2^*(S) dS, \quad (127) $$

with

$$ \begin{cases} -\Delta G_2^* = -\lambda |A^0| \frac{E'(|S|)}{|S|} S_1 & \text{in } |\Sigma| < 1, \\ G_2^* = 0 & \text{on } |S| = 1. \end{cases} \quad (128) $$

Consequently, there is a numerical constant $C_4$ such that:

$$ \mathcal{J}_2 = -C_4 \lambda^2 |A^0|^2, \quad (129) $$

and this constant has the following expression:

$$ C_4 = \frac{7}{6} \int_{|S| \le 1} S_1 \overline{G}_2(S_1, \hat{S}_1) dS, \quad (130) $$

where

$$ \begin{cases} -\Delta \overline{G}_2 = - \frac{E'(\sqrt{y^2 + |y'|^2})y}{\sqrt{y^2 + |y'|^2}}, \\ \overline{G}_2 = 0. \end{cases} \quad (131) $$

We now introduce the test function $\psi_2(y, y') = \frac{1}{16}y(1 - (y^2 + |y'|^2))$ which satisfies $-\Delta\psi_2 = y$ and $\psi_2$ vanishes on the sphere. Thus we have:

$$ \begin{align} C_4 &= \frac{7}{6} \int_{|S| \le 1} S_1 \overline{G}_2(S_1, \hat{S}_1) dS \nonumber \\ &= -\frac{7}{6} \int_{|S| \le 1} \frac{E'(|S|)S_1}{|S|} \frac{S_1}{16} (1 - |S|^2) dS \tag{132} \\ &= \frac{7}{8} \int_{|S| \le 1} \frac{\psi_0(S) |E'(|S|)|}{|S|} S_1^2 dS. \nonumber \end{align} $$
Comparing $C_4$ to $C_2$ from (106) we can see that $C_4/C_2 = \frac{7}{8} < 1$ and thus $C_6 = C_2 - C_4 > 0$.

Finally we have proven,

$$ \mathcal{I}_2 = - \left( C_3 \sum_i \lambda_i^2 + C_4 \lambda^2 |\mathcal{A}^0|^2 \right). \qquad (133) $$

In conclusion, one gets

$$ \begin{aligned} \partial_{\mathbf{a}\mathbf{a}}^2 F(0,0) &= \int_B \partial_{\mathbf{a}\mathbf{a}}^2 \tilde{q}^{0,0}(\Sigma) d\Sigma \\ &= (\mathcal{I}_1 - \mathcal{I}_2 - \mathcal{I}_3) \\ &= -(C_5 \sum_i \lambda_i^2 + C_6 \lambda^2 |\mathcal{A}_0|^2) < 0. \end{aligned} \qquad (134) $$

## 4.3 Verification of Hypothesis H2

We shall assume that $\lambda$ is large, so we can neglect $\mathcal{A}\Sigma$ relative to $\lambda\mathcal{A}_0$. Without loss of generality, we may assume that $\mathcal{A}_0 = (1,0,0,0,0,0)$. The equation satisfied by $q_0^a$ is then

$$ -\Delta q_0^a + \mathbf{a}\lambda \frac{\partial}{\partial \Sigma_1} q_0^a = -\mathbf{a}\lambda \frac{\partial E}{\partial \Sigma_1}, \qquad (135) $$

with Dirichlet boundary condition. We set $q_0^a = q - E$ to find

$$ -\Delta q + \mathbf{a}\lambda \frac{\partial q}{\partial \Sigma_1} = -\Delta E = \rho. \qquad (136) $$

We may further set $q = p \exp(\mathbf{a}\lambda\Sigma_1/2)$. This leads to the axisymmetric problem

$$ -\Delta p + \frac{\mathbf{a}^2 \lambda^2}{4} p = \rho \exp(-\mathbf{a}\lambda\Sigma_1/2). \qquad (137) $$

To simplify, we set $\mathbf{a}\lambda/2 = a$. With $\rho = \delta$, and Dirichlet conditions for $p$, this has the solution

$$ p(r) = \frac{a^2 K_2(ar)I_2(a) - I_2(ar)K_2(a)}{8\omega_5 I_2(a)r^2}. \qquad (138) $$

Here $\omega_5$ is the surface area of the 5-dimensional sphere. To verify the hypothesis, we need to evaluate

$$ \int_B q d\Sigma = \int_B p(r) \exp(a\Sigma_1) d\Sigma. \qquad (139) $$

We now set

$$ \rho = \sqrt{r^2 - \Sigma_1^2}, \quad \Sigma_1 = r \cos(\theta), \quad \rho = r \sin \theta. \qquad (140) $$

The integral then becomes

$$ \int_0^1 \int_0^\pi p(r) \exp(ar \cos \theta) \omega_4 r^5 (\sin \theta)^4 d\theta dr. \qquad (141) $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Plot of the function $S(a)$</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Plot of the integrand</figcaption></figure>

We can carry out the $\theta$ integration explicitly and finally end up with

$$ \frac{1}{a^2} \int_0^1 3\pi\omega_4 r^3 p(r) I_2(ar) dr. \qquad (142) $$

We omit irrelevant positive factors which do not depend on $a$, and we finally have to investigate the function

$$ S(a) = \int_{0}^{1} r \frac{K_{2}(ar)I_{2}(a) - I_{2}(ar)K_{2}(a)}{I_{2}(a)} I_{2}(ar) dr. \qquad (143) $$

Figure 1 shows a plot of this function, which is clearly decreasing with $a$. Indeed, it appears that the integrand in (143), for any fixed $r$, is a decreasing function of $a$. Figure 2 shows a plot of the function

$$ s(a, r) = \frac{K_{2}(ar)I_{2}(a) - I_{2}(ar)K_{2}(a)}{I_{2}(a)} I_{2}(ar). \qquad (144) $$
# 5 Proof of Corollary 1

This proof is the same as in the 1d case which can be found in [15]. We reproduce it here for the sake of completeness.

## 5.1 The case $\mu > \mu_c$

If $\mu > \mu_c$, let be $\mathfrak{b}_0$ be the unique positive number such that $F(0, \mathfrak{b}_0) = \mu$. Since we have $\partial_\mathfrak{b}F(0, \mathfrak{b}_0) > 0$ we can apply the implicit function theorem which defines a function $g$ such that $F(\mathfrak{a}, \mathfrak{b}) = \mu$ is locally equivalent to $\mathfrak{b} = g(\mu, \mathfrak{a})$. In terms of the parameters $\varepsilon$ and $\Gamma$ (see (13)) one finds the equation

$$ \sqrt{\mu\Gamma} = g\left(\mu, \frac{\varepsilon}{\mu\Gamma}\right). \qquad (145) $$

Fix $\mu$ and define $\tilde{g}(\varepsilon, \Gamma) = \mu\Gamma - g(\mu, \varepsilon/(\mu\Gamma))^2$ and you find that

* $\tilde{g}(0, \mathfrak{b}_0^2/\mu) = 0$,

* $\partial_\Gamma \tilde{g}(0, \mathfrak{b}_0^2/\mu) = \mu$,

so that $\Gamma$ can be expressed as an analytic function of $\varepsilon$ for small $\varepsilon$. Moreover at $\varepsilon = 0$ one finds $\Gamma = \mathfrak{b}_0^2/\mu$.

## 5.2 The case $\mu < \mu_c$

If $\mu < \mu_c$, let $\mathfrak{a}_0$ be the unique positive number such that $F(\mathfrak{a}_0, 0) = \mu$. Since we have $\partial_\mathfrak{a}F(\mathfrak{a}_0, 0) < 0$ we can apply the implicit function theorem which defines a function $g$ such that $F(\mathfrak{a}, \mathfrak{b}) = \mu$ is locally equivalent to $\mathfrak{a} = g(\mu, \mathfrak{b})$. In terms of the parameters $\varepsilon$ and $\Gamma$ (see (13)) one finds the equation

$$ \varepsilon = \Gamma\mu g\left(\mu, \sqrt{\Gamma\mu}\right). \qquad (146) $$

Since we have

$$ \partial_{\sqrt{\Gamma}} \left( \sqrt{\Gamma\mu} \sqrt{g(\mu, \sqrt{\Gamma\mu})} \right) \Big|_{\sqrt{\Gamma}=0} = \sqrt{\mu g(\mu, 0)} = \sqrt{\mu \mathfrak{a}_0} > 0, \qquad (147) $$

we can use the local inversion theorem to write $\sqrt{\Gamma}$ as a function of $\sqrt{\varepsilon}$ for small $\varepsilon$. Consequently, $\Gamma$ itself is a $C^\infty$ function of $\sqrt{\varepsilon}$ near $\varepsilon = 0$ with leading order $\Gamma \sim \varepsilon/(\mu\mathfrak{a}_0)$.

## 5.3 The case $\mu = \mu_c$

If $\mu = \mu_c$ we have $F(0, 0) = \mu_c$. Moreover,

$$ \partial_\mathfrak{b} F(0, 0) > 0, \qquad (148) $$

so that, by applying the implicit function theorem, we obtain that $F(\mathfrak{a}, \mathfrak{b}) = \mu_c$ is locally equivalent to $\mathfrak{b} = h(\mathfrak{a})$, and hence

$$ \sqrt{\mu\Gamma} = h\left(\frac{\varepsilon}{\mu\Gamma}\right). \qquad (149) $$
We set $\sqrt{\tilde{\Gamma}} = \tilde{\Gamma}\varepsilon^{2/5}$. In terms of $\varepsilon$ and $\tilde{\Gamma}$ we have the relation:

$$ \tilde{\Gamma}^5 = \frac{\tilde{\Gamma}^4}{\sqrt{\mu\varepsilon^{2/5}}} h\left(\frac{\varepsilon^{1/5}}{\mu\tilde{\Gamma}^2}\right) := \tilde{h}\left(\frac{\varepsilon^{1/5}}{\tilde{\Gamma}^2}\right), \quad (150) $$

where $\tilde{h}(a) = h(a/\mu)/(\sqrt{\mu a^2})$. We have $h(0) = 0$ and we calculate

$$ h'(0) = -\frac{\partial_a F(0, 0)}{\partial_b F(0, 0)} = 0 \quad (151) $$

$$ h''(0) = -\frac{\partial_{aa}^2 F(0, 0)}{\partial_b F(0, 0) \partial_b g(0, 0)} > 0. \quad (152) $$

Consequently $\tilde{h}$ is a $C^\infty$ function with a positive limit $\tilde{h}(0) = c > 0$. Finally, we apply the implicit function theorem one last time to express $\tilde{\Gamma}$ as a function of $\varepsilon^{1/5}$ which is possible since

$$ \partial_{\tilde{\Gamma}} \left( \tilde{\Gamma}^5 - \tilde{h} \left( \frac{\varepsilon^{1/5}}{\tilde{\Gamma}^2} \right) \right) \Big|_{\varepsilon^{1/5}=0, \tilde{\Gamma}=c^{1/5}} = 5c^{4/5} > 0. \quad (153) $$

We finally get that

$$ \Gamma = \varepsilon^{4/5} \tilde{\Gamma} (\varepsilon^{1/5})^2, \quad (154) $$

with $\tilde{\Gamma}$ a $C^\infty$ function converging to $c^{1/5}$ as $\varepsilon^{1/5}$ goes to 0.

# 6 The form of $\mathcal{A}$ and $\mathcal{A}_0$ and examples

## 6.1 The form of the drift term

The deterministic stress evolution, in dimensionless variables, has the form

$$ D_t T = \partial_t T + (u \cdot \nabla_x) T \\ = \varepsilon \left( \frac{1+a}{2}(MT + TM^T) - \frac{1-a}{2}(TM + M^T T) + \lambda(M + M^T) \right). \quad (155) $$

Here $\lambda = G_0/\sigma_c$, where $G_0$ is a stress modulus. In components, this reads

$$ D_t T_{ij} = \varepsilon \left( \sum_k \left[ \frac{1+a}{2}(M_{ik}T_{kj} + M_{jk}T_{ki}) - \frac{1-a}{2}(T_{ik}M_{kj} + T_{jk}M_{ki}) \right] + \lambda(M_{ij} + M_{ji}) \right). \quad (156) $$

We then need to transform this to the components of $\Sigma$. Instead of doing this in generality, we shall discuss the specifics for three classical rheological flows: the Couette flow and two types of elongational flows.

## 6.2 Couette flow

We are interested in comparing the results of our model with the original HL model. In stationary Couette flow, we have

$$ M = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}. \quad (157) $$
This leads to the equations

$$
\begin{align}
D_t T_{11} &= (1+a)T_{12}, \\
D_t T_{22} &= -(1-a)T_{12}, \\
D_t T_{33} &= 0, \\
D_t T_{12} &= \frac{1+a}{2} T_{22} - \frac{1-a}{2} T_{11} + \lambda, \\
D_t T_{13} &= \frac{1+a}{2} T_{23}, \\
D_t T_{23} &= -\frac{1-a}{2} T_{13}. \tag{158}
\end{align}
$$

This leads to a matrix $\mathcal{A}$ which is

$$
\mathcal{A} = \begin{pmatrix}
0 & 0 & 0 & \frac{\sqrt{2}}{2}(1+a) & 0 & 0 \\
0 & 0 & 0 & \frac{\sqrt{2}}{2}(-1+a) & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 \\
\frac{\sqrt{2}}{2}(-1+a) & \frac{\sqrt{2}}{2}(1+a) & 0 & 0 & 0 & \frac{1+a}{2} \\
0 & 0 & 0 & 0 & 0 & \frac{-1+a}{2} \\
0 & 0 & 0 & 0 & \frac{-1+a}{2} & 0
\end{pmatrix}, \quad (159)
$$

and a vector

$$
\mathcal{A}^0 = (0, 0, 0, \sqrt{2}, 0, 0)^T. \tag{160}
$$

6.3 Elongational flows

Elongational flow are characterized by a diagonal deformation rate tensor. We
distinguish two main types of elongational flows. The first kind is axisymmetric
flow in which

$$
M = \begin{pmatrix} 2 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & -1 \end{pmatrix}. \tag{161}
$$

In this case the matrix $\mathcal{A}$ is given by

$$
\mathcal{A} = \begin{pmatrix}
4a & 0 & 0 & 0 & 0 & 0 \\
0 & -2a & 0 & 0 & 0 & 0 \\
0 & 0 & -2a & 0 & 0 & 0 \\
0 & 0 & 0 & a & 0 & 0 \\
0 & 0 & 0 & 0 & a & 0 \\
0 & 0 & 0 & 0 & 0 & -2a
\end{pmatrix},
\quad (162)
$$

The vector $\mathcal{A}^0$ is given by

$$
A^0 = (4, -2, -2, 0, 0, 0)^T. \tag{163}
$$

The second kind of elongational flows are planar elongational flows for which

$$
M = \begin{pmatrix} 1 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & 0 \end{pmatrix}. \qquad (164)
$$
This leads to

$$A = \begin{pmatrix} 2a & 0 & 0 & 0 & 0 & 0 \\ 0 & -2a & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & a & 0 \\ 0 & 0 & 0 & 0 & 0 & -a \end{pmatrix}, \qquad (165)$$

and a vector

$$A^0 = (2, -2, 0, 0, 0, 0)^T. \qquad (166)$$

## References

[1] M. Abramowitz and I.A. Stegun. *Handbook of Mathematical Functions*. Dover, New York, 1965.

[2] S. Agmon, A. Douglis and L. Nirenberg. Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions. I. *Comm. Pure Appl. Math.*, 12:623-727, 1959.

[3] S. Agmon, A. Douglis and L. Nirenberg. Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions. II. *Comm. Pure Appl. Math.*, 17:35-92, 1964.

[4] A. Arnold, J.A. Carrillo, and C. Manzini. Refined long-time asymptotics for some polymeric fluid flow models. *Commun. Math. Sci.*, 8(3):763-782, 2010.

[5] S. Bénito, C.-H. Bruneau, T. Colin, C. Gay and F. Molino. An elastovisco-plastic model for immortal foams or emulsions. *Eur. Phys. J. E*, 25:225-251, 2008.

[6] R.B. Bird, O. Hassager, R.C. Armstrong and C.F. Curtiss. *Dynamics of Polymeric Liquids*, Vol. 2. Wiley, New York, 2nd ed. 1987.

[7] D. Bresch and J. Simon. Western boundary currents versus vanishing depth. *Discrete Contin. Dyn. Syst. Ser. B*, 3(3):469-477, 2003.

[8] M.E. Cates and P. Sollich. Tensorial constitutive models for disordered foams, dense emulsions and other soft non-ergodic materials. *J. Rheol*, 48:193-207, 2004.

[9] P. Hébraud and F. Lequeux. Mode-coupling theory for the pasty rheology of soft glassy materials. *Phys. Rev. Lett.*, 81(14):2934-2937, 1998.

[10] M.E.H. Ismail and M.E. Muldoon. Monotonicity of the zeros of a cross-product of Bessel functions. *SIAM J. Math. Anal.*, 9:759-767, 1978.

[11] P. Jop, Y. Forterre and O. Pouliquen. A constitutive law for dense granular flows. *Nature*, 441:727-730, 2006.
[12] B. Jourdain, C. Le Bris, T. Lelièvre and F. Otto. Long-time asymptotics of a multiscale model for polymeric fluid flows. *Arch. Ration. Mech. Anal.*, 181(1):97-148, 2006.

[13] J.G. Oldroyd. On the formulation of rheological equations of state. *Proc. Roy. Soc. London. Ser. A.*, 200:523-541, 1950.

[14] J. Olivier. Asymptotic analysis in flow curves for a model of soft glassy rheology. *Z. Angew. Math. Phys*, 61(3):445-466, 2010.

[15] J. Olivier and M. Renardy. Glass transition seen through asymptotic expansions. *SIAM J. Appl. Math.*, 71(4):1144-1167, 2011.

[16] P. Sollich, F. Lequeux, P. Hébraud, and M.E. Cates. Rheology of soft glassy materials. *Phys. Rev. Lett.*, 78(10):2020-2023, 1997.