Asymptotic States of a Smoluchowski Equation

P. Constantin^{∗}, I. G. Kevrekidis^{†}, E. S. Titi^{‡}

∗ Department of Mathematics, The University of Chicago, Chicago, Il 60637. †Chemical Engineering, PACM and Mathematics, Princeton University, Princeton, NJ 08544. ‡Department of Mathematics, and Department of Mechanical and Aerospace Engineering,

University of California, Irvine, CA 92697-3875, Department of Computer Science and Applied Mathematics, Weizmann Institute of Science, Rehovot 76100, Israel.

We study the high concentration asymptotics of steady states of a Smoluchowski equation arising in the modeling of nematic liquid crystalline polymers.

1 Introduction

There are many levels of models describing the rheology of non-Newtonian complex fluids containing liquid crystalline polymers. Some descriptions combine macroscopic partial differential equations with microscopic stochastic differential equations ([10], [12], [14], [19], [16]).

A simple kinetic model of nematic liquid crystalline polymers - the rigid rod model - using the Maier-Saupe potential, gives rise to a Smoluchowski equation for the single particle distribution function on the unit sphere ([4]). In spite of its simplicity, this equation exhibits nontrivial nonlinear dynamical features, in contrast with classical Fokker-Planck equations for noninteracting particles ([11]). At high concentrations, the shape of the particles in suspension becomes important. The complex dynamical properties are then amplified considerably in the presence of symmetry breaking shear ( [5], [6], [7], [8], [13], [18]).

Our work addresses the transition to order first described by Onsager in his seminal paper ([15]) in which he developed a thermodynamic formalism for dilute colloidal solutions. Onsager calculated the free energy using cluster expansions and approximations of the forces between rod-like particles. He arrived at an expression for the free energy in terms of a configuration integral involving the distribution function ψ of particle orientations and an interparticle potential interaction kernel β. Onsager wrote the Euler-Lagrange equation for the variation of the configuration integral retaining the first nontrivial term in the cluster expansion. This nonlinear integral equation (see (2) below) is the same as the equation (10, 11) solving the time independent Smoluchowski equation (4). Different expressions used for the function β defining the interaction potential give rise to different models. The steady states of the Smoluchowski equation are obtained in the form

ψ = Z^{−1}e^{−V}

where Z is a constant used to enforce the normalization

∫_{S^{2}} ψ(ϕ)dσ = 1.

(dσ is the area element and ϕ are coordinates on the unit sphere S^{2} in R^{3}). The potential

V (ϕ) = −b ∫_{S^{2}} β(ϕ, ϕ′)ψ(ϕ′)dσ(ϕ′)  (1)

is given in terms of ψ. The function β embodies the interaction between the particles in suspension. The constant intensity b > 0 is expressed as the product b = cυ where c = N/V is the concentration (number of particles per volume) and υ is an excluded volume depending on the shape of the particles in the suspension. Onsager’s paper was concerned with the derivation of the function β and the study of the limit of high concentration. Taking the logarithm of the representation of the steady solutions of the Smoluchowski equation one arrives at

log(ψ(ϕ)) = log(Z^{−1}) + b ∫_{S^{2}} β(ϕ, ϕ′)ψ(ϕ′)dσ(ϕ′).  (2)

Equation (2) is precisely the equation studied by Onsager in his seminal paper [15] (equation (69), page 643). Onsager derived complicated empirical expressions for β but in the end resorted to a simple expression (equation (81) on page 647) which is proportional to − sin γ where γ ∈ [0, π] is the angle between the unit vectors x(ϕ), x(ϕ′). For the Maier-Saupe potential, which will be used in this work, the function β is

β(ϕ, ϕ′) = (x(ϕ) · x(ϕ′))^{2} − \frac{1}{3} = (cos γ)^{2} − \frac{1}{3}.  (3)

The important property shared by both the Maier-Saupe potential, the explicit example studied by Onsager and also by his empirically derived expressions is that −β is an increasing function of sin γ which has a minimum when the directions are parallel and a maximum when they are perpendicular ([15], pp. 644 and 646). The results in Onsager’s paper are based on an explicit ansatz for ψ (formula (80) on page 647) which “decreases rather too rapidly for large angles” but which “was, nevertheless, adopted as the best tractable function” ([15], p. 647). Using this ansatz Onsager was able to argue that in the limit of b → ∞ one has a transition from the isotropic uniform distribution to an ordered prolate distribution. His approach was variational, and, because he had to content himself with the ansatz in formula (80) of ([15]), the results were explicit, but not rigorously mathematically proved. In this paper we study the Smoluchowski equation on the unit sphere with the Maier-Saupe potential. This choice of the potential allows us to investigate rigorously the asymptotics of the steady state solutions for large values of the potential intensity, corresponding to large concentrations. We reduce the problem of finding steady state solutions of the Smoluchowski PDE with Maier-Saupe potential to the finite dimensional problem of finding the eigenvalues of a symmetric, traceless matrix. Linear combinations of these eigenvalues are critical points of a function associated to them. This description is key to the asymptotic analysis. We find multiple steady solutions which are clustered in three distinct groups. As the concentration is increased, these steady solutions converge to uniform, prolate and oblate states, confirming rigorously the transition discovered by Onsager. (On physical grounds one might suspect that the uniform state becomes unstable and that the oblate states are saddles, but we do not address this issue in this paper). Furthermore, the methods of study allow us to devise asymptotic expansions for the steady states, expansions that are valid a high but finite concentrations. These expansions are a first step towards a more comprehensive understanding of the long time dynamics of the Smoluchowski equation, and a preparation for the study of symmetry breaking perturbations.

2 Smoluchowski Equations

Consider a smooth compact connected Riemannian manifold without boundary (M^{n}, g) ([9]) and a real, symmetric smooth function

β : M × M → R,

β(m, p) = β(p, m). One can associate to β a linear operator

ψ ↦ V

given by

V (m) = −b ∫_{M} β(m, p)ψ(p)dσ(p)

where dσ is the Riemannian volume element. The Smoluchowski equation is, in local coordinates,

∂_{t}ψ = \frac{1}{\sqrt{g}} ∂_{i} (e^{−V} \sqrt{g}g^{ij}∂_{j}(e^{V} ψ)) .  (4)

We use the summation convention. The equation is a nonlinear Fokker-Planck equation (that is: it is a nonlinear equation, and it looks like a linear Fokker-Planck equation),

∂_{t}ψ = Δ_{g}ψ + div_{g}(ψ∇V)

where

Δ_{g} = \frac{1}{\sqrt{g}} ∂_{i} (\sqrt{g}g^{ij}∂_{j})

is the Laplace-Beltrami operator and the last term is

div_{g}(ψ∇V) = \frac{1}{\sqrt{g}} ∂_{i} (\sqrt{g}g^{ij}ψ∂_{j}V ) .

The Smoluchowski equation preserves mass,

∫_{M} ψ_{t}dσ = 0.

Smoluchowski equations have an energy functional

E = ∫_{M}(log ψ)ψdσ + \frac{1}{2} ∫_{M} V ψdσ

that is a non-increasing function of time when evaluated on solutions. Indeed, taking the derivative of E(ψ) when the time dependence comes from a smooth positive solution ψ = ψ(p, t) of the Smoluchowski equation, one obtains

\frac{d}{dt}E = − ∫_{M} |∇_{g}(V + log ψ)|^{2}ψdσ  (5)

In the expression above

|∇_{g}f|^{2} = g^{ij}∂_{i}f∂_{j}f.

One can also observe that log ψ + V is the formal density of the first variation \frac{δE}{δψ}.

\frac{δE}{δψ} χ = ∫_{M} (log ψ + V) χdσ

This follows if one assumes that the variations χ have vanishing integral (because the integral of ψ can be held constant). The fact that the map ψ → V is linear and symmetric is also needed for the above calculations. Therefore one may write, formally

\frac{d}{dt}E = − ∫_{M} |∇_{g} \frac{δE}{δψ} |^{2} ψdσ.

Many nonlinear equations share this dissipative structure, for instance lubrication approximations of Hele-Shaw problems ([2]), porous medium equations ([17]), and the Keller-Segel chemotaxis equation, ([1]). In fact, the latter is a Smoluchowski equation with non-smooth β.

It follows from the maximum principle that, if the initial datum f_{0} is a nonnegative (positive) function, then the solution of the Smoluchowski equation remains nonnegative (positive). The assumption of smoothness of β implies easily that:

Theorem 2.1 Let f_{0} be a nonnegative, continuous function on M. The solutions of (4) with initial data ψ(·, 0) = f_{0} exist for all positive time, are smooth (C^{∞}), nonnegative and normalized

∫_{M} ψ(m, t)dσ(m) = ∫_{M} f_{0}(m)dσ(m)  (6)

The proof can be done by successive approximations, and will be omitted. The smoothness of β is crucial: there are simple proofs of finite time blow up for Keller-Segel chemotaxis equations ([1]).

3 Steady States

We consider steady states of (4). From the decay of the energy functional (5) one deduces that any positive time independent solution of (4) must satisfy

ψ = Z^{−1}e^{−V}  (7)

for an appropriate positive constant Z.

From now on we are going to specialize to the problem of interest to us , in which M = S^{2} be the unit sphere in R^{3}. We will use without loss of generality the normalization

∫_{S^{2}} ψdσ = 1.  (8)

We consider local coordinates (ϕ = (θ, φ)). The Maier-Saupe potential is given by

V (x, t) = −bx_{i}x_{j}S^{ij}  (9)

where x_{i} are Cartesian coordinates in R^{3}, i, j = 1, 2, 3, and b is a positive constant. The matrix S is determined by

S^{ij}(t) = ∫_{S^{2}} x_{i}(ϕ)x_{j}(ϕ)ψ(ϕ, t)σ(dϕ) − \frac{1}{3}δ_{ij}  (10)

with σ(dϕ) = \sqrt{g}dϕ the surface area. Thus, V (x, t) is a homogeneous polynomial of second degree, restricted to the sphere. The coordinates on the two dimensional unit sphere are ϕ = (θ, φ), x_{1}(θ, φ) = sin θ cos φ, x_{2}(θ, φ) = sin θ sin φ, x_{3}(θ, φ) = cos θ. Recall also that

g^{11} = 1, g^{22} = (sin θ)^{−2}, g^{ij} = 0, i ≠ j

with ∂_{θ} = ∂_{1} and ∂_{φ} = ∂_{2} and that \sqrt{g} = sin θ.

In view of (7), (9) it follows that the steady states can be represented by

ψ(ϕ) = Z^{−1}e^{−V} = Z^{−1}e^{bS^{ij}x_{i}(ϕ)x_{j}(ϕ)}.  (11)

The matrix S is real, symmetric and traceless. This follows from the definition

S = ∫_{S^{2}} (x ⊗ x)ψdσ − \frac{1}{3}I.

The eigenvalues of S must lie between −1/3 and 2/3. Indeed, for any unit vector ξ one has from (10) that

Sξ · ξ = ∫_{S^{2}} (ξ · x(ϕ))^{2} ψ(ϕ)dσ(ϕ) − \frac{1}{3}  (12)

and because 0 ≤ (ξ · x(ϕ))^{2} ≤ 1 and the normalization (8) we deduce that the integral in the expression above has a value between 0 and 1.

The uniform distribution is the special solution for which the matrix S vanishes, Z is the area of S^{2} and ψ = Z^{−1}:

ψ_{0} = \frac{1}{4π} , S_{0} = 0.

In order to parametrize all steady solutions we consider ([3]) the real valued map

(S, b) ↦ Z(S, b)  (13)

defined for any real, symmetric, traceless matrix S and any positive b by the formula

Z(S, b) = ∫_{S^{2}} e^{bS^{ij}x_{i}(ϕ)x_{j}(ϕ)}dσ(ϕ).  (14)

We also consider the function

ψ_{S,b}(ϕ) = (Z(S, b))^{−1}e^{b(S^{ij}x_{i}(ϕ)x_{j}(ϕ))}  (15)

associated to any real, traceless, symmetric S and b > 0. Finally, for any real, traceless symmetric S and b > 0, denote

(Ŝ(S, b))^{ij} = ∫_{S^{2}} x_{i}(ϕ)x_{j}(ϕ)ψ_{S,b}(ϕ)dσ(ϕ).  (16)

Obviously Ŝ is a function of S and b. Actually, one can check that Z(S, b) depends only on the conjugacy class OSO^{−1}, O ∈ O(3). More specifically, if S_{1} = OSO^{−1} then the rotation invariance of the measure on the unit sphere implies that Z(S, b) = Z(S_{1}, b) and therefore ψ_{S,b}(ϕ) = ψ_{S_{1},b}(T ϕ) where T ϕ is the angle translation associated to the rotation O, Ox(ϕ) = x(T ϕ). The rotation invariance implies then that Ŝ(S_{1}, b) = O (Ŝ(S, b)) O^{−1}. Clearly, by construction

∫_{S^{2}} ψ_{S,b}(ϕ)dσ(ϕ) = 1.  (17)

In view of the considerations above we have:

Theorem 3.1 The positive, normalized steady solutions of (4) are in one-to-one correspondence with the solutions of the implicit trancendental matrix equation

Ŝ(S, b) = S + \frac{1}{3}I  (18)

where Ŝ(S, b) is associated to S and b by the formalism (14), (15), (16) above.

Because of rotation invariance, without loss of generality, we may restrict our attention to diagonal matrices

S^{ij} = λ_{i}δ_{ij}  (19)

with λ_{1}, λ_{2}, λ_{3} real eigenvalues that obey

λ_{1} + λ_{2} + λ_{3} = 0  (20)

and belong to the interval [−\frac{1}{3}, \frac{2}{3}]. The search for steady solutions then reduces to a search for the eigenvalues λ_{1}, λ_{2}, λ_{3}. It turns out that the eigenvalues solve a coupled system of equations that describe the critical points of a functional. In order to present the calculations it is convenient to change variables, from (λ_{1}, λ_{2}) to (v_{1}, v_{2}) defined by:

v_{1} = \frac{1}{2}(λ_{1} + λ_{2}), v_{2} = \frac{1}{2}(λ_{1} − λ_{2}).  (21)

We will also use the vector notation v = (v_{1}, v_{2}), and for convenience some calculations will be performed in the scaled variables u = (u_{1}, u_{2}) = bv

u_{1} = \frac{b}{2}(λ_{1} + λ_{2}), u_{2} = \frac{b}{2}(λ_{1} − λ_{2}).  (22)

We consider the convex compact

K = [−1, 1] × [0, 2π] = {(p, t); −1 ≤ p ≤ 1, 0 ≤ t ≤ 2π}  (23)

and we consider the pair of functions

y_{1}(p) = 1 − 3p^{2}  (24)

and

y_{2}(p, t) = (1 − p^{2}) cos t  (25)

defined for (p, t) ∈ K. We write

y = y(p, t) = (y_{1}(p), y_{2}(p, t)).

These functions and this compact are used to describe the transcendental equations obeyed by the eigenvalues of the matrices S corresponding to solutions of (18), steady states of (4) with Maier-Saupe potential.

Theorem 3.2 Consider, for any u = (u_{1}, u_{2}) the function

Z_{2}(u) = ∫_{K} e^{u·y(p,t)}dpdt  (26)

and associate to it the function

F(u) = log(Z_{2}(u)) − \frac{1}{b} (3u_{1}^{2} + u_{2}^{2}) .  (27)

Then the solutions of (18) coincide (via (19), (20), (22)) with the critical points u = (u_{1}, u_{2}) ∈ [−\frac{b}{3}, \frac{2b}{3}] × [0, \frac{b}{2}] of F.

The critical point equations

∇_{u}F = 0  (28)

can be written as

[y_{1}](u) = \frac{6u_{1}}{b} , [y_{2}](u) = \frac{2u_{2}}{b} .  (29)

where

[F] (u) = (Z_{2}(u))^{−1} ∫_{K} F (p, t)e^{u·y(p,t)}dpdt  (30)

(i). If 0 < b < 1/2 the function F is strictly concave and has a unique critical point at u = 0. The corresponding unique steady state of (4) is the uniform state ψ_{0}.

(ii). If b ≥ 8 then u = 0 is an isolated critical point. Consequently, no bifurcations from the uniform state ψ_{0} can occur in (4) for b ≥ 8.

(iii). If 0 ≤ b < 4 then on any line u_{1} = const there is at most one critical point. If b ≥ 4 then the number of critical points on each line u_{1} = const does not exceed 2 [\frac{b}{4}].

Proof. A simple computation using the definition (16) for diagonal matrices (19) shows that

Ŝ^{ij}(S, b) = 0 for i ≠ j.

The equations (18) reduce then to

∫_{S^{2}} x_{i}^{2}(ϕ)ψ_{S,b}(ϕ)dσ(ϕ) = λ_{i} + \frac{1}{3}

for i = 1, 2, 3. Because of the normalization (17), there are only two independent equations. The equations are

λ_{1} + \frac{1}{3} = Z^{−1} ∫_{0}^{2π} ∫_{0}^{π} cos^{2} φ sin^{3} θe^{bλ·X(θ,φ)}dθdφ  (31)

and

λ_{2} + \frac{1}{3} = Z^{−1} ∫_{0}^{2π} ∫_{0}^{π} sin^{2} φ sin^{3} θe^{bλ·X(θ,φ)}dθdφ  (32)

together with

Z(bλ) = ∫_{0}^{2π} ∫_{0}^{π} e^{bλ·X(θ,φ)} sin θdθdφ  (33)

where

λ = (λ_{1}, λ_{2}), X(θ, φ) = (X_{1}(θ, φ), X_{2}(θ, φ)) ,
λ · X(θ, φ) = λ_{1}X_{1}(θ, φ) + λ_{2}X_{2}(θ, φ)

and

X_{1}(θ, φ) = cos^{2} φ sin^{2} θ − cos^{2} θ  (34)

and

X_{2}(θ, φ) = sin^{2} φ sin^{2} θ − cos^{2} θ.  (35)

Notice that

X_{1} + X_{2} = 1 − 3 cos^{2} θ

can be used to express

cos^{2} φ sin^{2} θ = X_{1} + \frac{1}{3} (1 − X_{1} − X_{2})

and

sin^{2} φ sin^{2} θ = X_{2} + \frac{1}{3} (1 − X_{1} − X_{2}) .

We use the notation

[F] (bλ) = Z^{−1} ∫_{0}^{2π} ∫_{0}^{π} F (X(θ, φ))e^{bλ·X(θ,φ)} sin θdθdφ  (36)

Then, the system (31, 32) can be written as

λ_{1} = \frac{2}{3}[X_{1}] − \frac{1}{3}[X_{2}]  (37)

and

λ_{2} = \frac{2}{3}[X_{2}] − \frac{1}{3}[X_{1}].  (38)

Inverting this linear system, we have

[X_{1}] = 2λ_{1} + λ_{2}  (39)

and

[X_{2}] = λ_{1} + 2λ_{2}  (40)

Then we consider

Z_{2}(u) = ∫_{0}^{2π} ∫_{0}^{π} e^{u·Y (θ,φ)} sin θdθdφ  (41)

with Y (θ, φ) = (Y_{1}(θ, φ), Y_{2}(θ, φ)) defined by

Y_{1}(θ, φ) = sin^{2} θ − 2 cos^{2} θ  (42)

and

Y_{2}(θ, φ) = sin^{2} θ cos(2φ)  (43)

and with u = (u_{1}, u_{2}) ∈ [−\frac{b}{3}, \frac{2b}{3}] × [0, \frac{b}{2}]. The variables Y are related to X of (34, 35) via Y_{1} = X_{1} +X_{2}, Y_{2} = X_{1} −X_{2} The system (39), (40) which represents the steady solutions of (4), can be seen to be in one-to-one correspondence with the critical points of the function F(u) defined in (27) via (41, 42, 43). Indeed, the critical points satisfy the implicit equations

[Y_{1}] = \frac{6}{b}u_{1}  (44)

and

[Y_{2}] = \frac{2}{b}u_{2}  (45)

where

[F](u) = (Z_{2}(u))^{−1} ∫_{0}^{2π} ∫_{0}^{π} F (θ, φ)e^{u·Y (θ,φ)} sin θdθdφ.  (46)

The equations (44), (45) are equivalent to (39), (40). Changing variables p = cos θ and t = 2ϕ, he functions Y_{1} and Y_{2} become, in variables p, t, the functions y_{1}, y_{2} of (24), (25), the expected value (46) is the same as (30), and the equations (44) and (45) are the same as (29). When the parameters (u_{1}, u_{2}) are chosen to satisfy the implicit equations (29), then

[F](u) = ∫_{S^{2}} F (ϕ)ψ_{S,b}(ϕ)dσ(ϕ)  (47)

does represent the expected value of the function F at the corresponding steady state ψ_{S,b}.

In order to prove (i) we compute the Hessian of F. H(u) = (\frac{∂^{2}F}{∂u_{i}∂u_{j}} ) is given by

H(u) = ( \begin{matrix}[ξ^{2}] − \frac{6}{b} & [ξη] \\ [ξη] & [η^{2}] − \frac{2}{b} \end{matrix} )  (48)

where

ξ = Y_{1} − [Y_{1}],  (49)

and

η = Y_{2} − [Y_{2}].  (50)

Using the same notation, but changing variables to p, t, the Hessian of F(u) defined in (27), H(u) = (\frac{∂^{2}F}{∂u_{i}∂u_{j}} ) is given by

H(u) = ( \begin{matrix}[ξ_{1}^{2}] − \frac{6}{b} & [ξ_{1}ξ_{2}] \\ [ξ_{1}ξ_{2}] & [ξ_{2}^{2}] − \frac{2}{b} \end{matrix} )  (51)

where

ξ_{1} = y_{1} − [y_{1}],  (52)

and

ξ_{2} = y_{2} − [y_{2}].  (53)

The concavity for small b follows from

\begin{matrix}H(u) :: a ⊗ a = −\frac{2}{b}(3a_{1}^{2} + a_{2}^{2}) − (a · [y])^{2} \\ +a_{1}^{2}[y_{1}^{2}] + a_{2}^{2}[y_{2}^{2}] + 2a_{1}a_{2}[y_{1}y_{2}] \end{matrix}  (54)

Using the fact that the functions y_{1}, y_{2} have ranges included in the interval [−2, 1], and respectively, [−1, 1] and neglecting the non-positive (but unknown) contribution −(a · [y])^{2}, we arrive at an explicit sufficient condition b < \frac{1}{2} for the concavity. This concludes the proof of (i).

The proof of (ii) also uses the Hessian, but it requires a more careful analysis. Let us write

F(u) = F_{2}(u) − \frac{1}{b}(3u_{1}^{2} + u_{2}^{2}), with F_{2}(u) = log(Z_{2}(u))  (55)

and

H_{2}(u) = \frac{∂^{2}F_{2}(u)}{∂u_{i}∂u_{j}}.  (56)

Then we have, for arbitrary a = (a_{1}, a_{2})

H_{2}(u) :: a ⊗ a = [(a · ξ)^{2}] .  (57)

This shows that F_{2} is convex. In order to find explicit bounds we start by writing

(a · ξ)^{2} = (a · y)^{2} + (a · [y])^{2} − 2(a · y)(a · [y]).  (58)

We take a probability measure dπ on K = [−1, 1] × [0, 2π] and denote

⟨F ⟩ = ∫_{−1}^{1} ∫_{0}^{2π} F (p, t)dπ.  (59)

Integrating (58) with respect to dπ we obtain

⟨(a · ξ)^{2}⟩ = (a · ([y] − ⟨y⟩))^{2} + ⟨(a · (y − ⟨y⟩))^{2}⟩.  (60)

If we use normalized Lebesgue measure dπ = \frac{1}{4π} dtdp we note that

< y_{i} >= \frac{1}{4π} ∫_{−1}^{1} (∫_{0}^{2π} y_{i}(p, t)dt) dp = 0

for i = 1, 2. Because of this and (60) we deduce

\frac{1}{4π} ∫_{−1}^{1} (∫_{0}^{2π}(a · ξ)^{2}dt) dp = \frac{1}{4π} ∫_{−1}^{1} (∫_{0}^{2π}(a · y)^{2}dt) dp + (a · [y])^{2}

= \frac{4}{5}a_{1}^{2} + \frac{1}{3}a_{2}^{2} + (a · [y])^{2}.  (61)

Now it is easy to see, using the facts −2 ≤ y_{1} ≤ 1 and −1 ≤ y_{2} ≤ 1 that

e^{−4|u_{1}|−2|u_{2}|} \frac{1}{4π} ∫_{−1}^{1} (∫_{0}^{2π} F (p, t)dt) dp ≤ [F](u)  (62)

holds for any nonnegative function F, and in particular for F = (a · ξ)^{2}. We deduce the strict convexity inequality

H_{2}(u) :: a ⊗ a ≥ e^{−(4|u_{1}|+2|u_{2}|)} (\frac{4}{5}a_{1}^{2} + \frac{1}{3}a_{2}^{2}) .  (63)

Consequently

H(u) :: a ⊗ a ≥ c|a|^{2}  (64)

with c > 0 if

b ≥ 8e^{4|u_{1}|+2|u_{2}|}.  (65)

In particular, for b ≥ 8 the state u = 0 is an isolated critical point. This concludes the proof of (ii).

The proof of (iii) requires computing more explicitly Z_{2}(u).

Z_{2}(u) = ∫_{−1}^{1} e^{u_{1}(1−3p^{2})} (∫_{0}^{2π} e^{u_{2}(1−p^{2}) cos t}dt) dp.  (66)

The object in the inner paranthesis (encountered in the two dimensional study ([3]) is

Z_{1}(r) = ∫_{0}^{2π} e^{r cos t}dt.  (67)

This has an explicit expression

Z_{1}(r) = 2π ∑_{k=0}^{∞} r^{2k}2^{−2k} \frac{1}{(k!)^{2}} .

Substituting in the expression above we deduce

Z_{2}(u) = 2π ∑_{k=0}^{∞} C_{2k}(u_{1})u_{2}^{2k}2^{−2k} \frac{1}{(k!)^{2}}  (68)

where

C_{2k}(u_{1}) = ∫_{−1}^{1}(1 − p^{2})^{2k}e^{u_{1}(1−3p^{2})}dp.  (69)

Now we observe that (45) is equivalent to

\frac{∂Z_{2}}{∂u_{2}} − \frac{2u_{2}}{b} Z_{2} = 0.

The expression for this is

\begin{matrix}\frac{∂Z_{2}}{∂u_{2}} − \frac{2u_{2}}{b} Z_{2} = \\ −\frac{8π}{b} ∑_{k=1}^{∞} k [kC_{2(k−1)}(u_{1}) − \frac{b}{4}C_{2k}(u_{1})] (\frac{u_{2}}{2} )^{2k−1} \frac{1}{(k!)^{2}} . \end{matrix}  (70)

(This situation is very similar to the two dimensional situation, except that in that case the coefficients C_{2k}(u_{1}) were identically equal to 1). We observe that

0 ≤ C_{2k}(u_{1}) ≤ C_{2(k−1)}(u_{1})

holds. If b ≥ 4 we write

\frac{∂Z_{2}}{∂u_{2}} − \frac{2u_{2}}{b} Z_{2} = P (u) − Q(u)  (71)

where

P (u) = −\frac{8π}{b} \begin{matrix}[ & \frac{b}{4} & ] \\ ∑ \\ k=1 \end{matrix} k [kC_{2(k−1)}(u_{1}) − \frac{b}{4}C_{2k}(u_{1})] (\frac{u_{2}}{2} )^{2k−1} \frac{1}{(k!)^{2}}  (72)

and

Q(u) = \frac{8π}{b} ∑_{k=1+[\frac{b}{4}]}^{∞} k [kC_{2(k−1)}(u_{1}) − \frac{b}{4}C_{2k}(u_{1})] (\frac{u_{2}}{2} )^{2k−1} \frac{1}{(k!)^{2}} .  (73)

We observe therefore that

(\frac{∂}{∂u_{2}} )^{m} (\frac{∂Z_{2}}{∂u_{2}} − \frac{2u_{2}}{b} Z_{2}) < 0  (74)

holds for m ≥ 2 [\frac{b}{4}], which implies by Rolle’s theorem that there are at most 2 [\frac{b}{4}] distinct critical points on each line u_{1} = const. This concludes the proof of (iii) and of the theorem.

We will study now the asymptotics at large b, for fixed λ_{1}, λ_{2}, λ_{3}. We recall that

u_{1} = bv_{1}, u_{2} = bv_{2}.  (75)

We study thus the asymptotics as b → ∞, for fixed v = (v_{1}, v_{2}) ∈ [−\frac{1}{3}, \frac{2}{3}]×[0, \frac{1}{2}]. The system (29) determining the steady solutions of (4) with eigenvalues (21) is

[y_{1}](bv) = 6v_{1}, [y_{2}](bv) = 2v_{2}  (76)

Theorem 3.3 The steady solutions of (4), given by (11) can be parametrized by the intensity b > 0 and by the numbers v_{1} ∈ [−\frac{1}{3}, \frac{2}{3}],and v_{2} ∈ [0, 1] associated to the eigenvalues λ_{1}, λ_{2}, λ_{3} of the real, traceless, symmetric matrix S by the formulae

v_{1} = \frac{1}{2} (λ_{1} + λ_{2}) , v_{2} = \frac{1}{2} (λ_{1} − λ_{2}) .

If v = (v_{1}, v_{2}) is fixed and b is sufficiently large then the following cases are present. If v belongs to a compact subset of the region R_{1} = \{v = (v_{1}, v_{2}); −\frac{1}{3} ≤ v_{1} ≤ \frac{2}{3}, 0 < v_{2} ≤ \frac{1}{2}, 3v_{1} + v_{2} > 0\} then, for large enough b, all steady solutions in this region approach v = (\frac{1}{6}, \frac{1}{2}) as b → ∞ and consequently their eigenvalues converge to

λ_{1} = \frac{2}{3}, λ_{2} = −\frac{1}{3}, λ_{3} = −\frac{1}{3}.

The expected value of a function f(x),

[f] = ∫_{S^{2}} f(x(ϕ))ψ_{S,b}(ϕ)dσ(ϕ)

in the asymptotic steady state in this region is given by (77):

lim_{b→∞}[f] = f(e_{1}).

where e_{1} = (1, 0, 0) ∈ S^{2}. If (v_{1}, v_{2}) belongs to a compact subset of the region R_{2}∪R_{3}∪R_{4} where R_{2} = \{v = (v_{1}, v_{2}); −\frac{1}{3} ≤ v_{1} ≤ \frac{2}{3}, 0 < v_{2} ≤ \frac{1}{2}, 3v_{1}+v_{2} < 0\}, R_{3} = \{v = (v_{1}, v_{2}); −\frac{1}{3} ≤ v_{1} ≤ \frac{2}{3}, 0 < v_{2} ≤ \frac{1}{2}, 3v_{1} + v_{2} = 0\} and R_{4} = \{v = (v_{1}, 0); −\frac{1}{3} ≤ v_{1} < 0\}, then, for sufficiently large b, all steady solutions in this region approach v = (−\frac{1}{3}, 0) as b → ∞, and consequently their eigenvalues converge to

λ_{1} = −\frac{1}{3}, λ_{2} = −\frac{1}{3}, λ_{3} = \frac{2}{3}.

The expected value of a function f in the asymptotic steady state in this region is given by (85):

lim_{b→∞}[f] = f(e_{3})

where e_{3} = (0, 0, 1) ∈ S^{2}.

If v belongs to a compact subset of R_{5} = \{v = (v_{1}, 0); 0 < v_{1} ≤ \frac{2}{3}\} then, for sufficiently large b, all steady solutions in this region approach v = (\frac{1}{6}, 0) as b → ∞, and consequently their eigenvalues converge to

λ_{1} = \frac{1}{6}, λ_{2} = \frac{1}{6}, λ_{3} = −\frac{1}{3}.

The expected value of a function f in the asymptotic steady state in this region is given by (83)

lim_{b→∞}[f] = \frac{1}{2π} ∫_{0}^{2π} f(cos φ, sin φ, 0)dφ

Proof. In order to study the asymptotics for large b and fixed v we have to divide the v plane in five different regions.

Case I. If 3v_{1} + v_{2} > 0, v_{2} > 0 then, for any F ∈ C^{1}, 2π periodic function of t, we have

lim_{b→∞}[F](bv) = F (0, 0).  (77)

Indeed, if 3v_{1} + v_{2} > 0, v_{2} > 0 we multiply by e^{−b(v_{1}+v_{2})} both the numerator and the denominator of the ratio

[F](bv) = \frac{∫_{−1}^{1} ∫_{0}^{2π} F (p, t)e^{bv·y}dtdp}{∫_{−1}^{1} ∫_{0}^{2π} e^{bv·y}dtdp}.

Thus,

[F](bv) = \frac{∫_{−1}^{1} ∫_{0}^{2π} F (p, t)e^{bv·(y−(1,1))}dtdp}{∫_{−1}^{1} ∫_{0}^{2π} e^{bv·(y−(1,1))}dtdp}.

But

v · (y − (1, 1)) = −3p^{2}v_{1} − v_{2}(1 − (1 − p^{2}) cos t) =
−((3v_{1} + v_{2})p^{2} + v_{2}(1 − p^{2})(1 − cos t))

which is strictly negative, except when p = 0 and cos t = 1. For ϵ > 0, the contributions coming from regions |p| ≥ ϵ or |1 − cos t| ≥ ϵ are uniformly exponentially small. Choosing δ so that t ∈ [δ, 2π − δ] implies |1 − cos t| ≥ ϵ we have thus

[F](bv) = \frac{∫_{−ϵ}^{ϵ} ∫_{−δ}^{δ} F (p, t)e^{−b((3v_{1}+v_{2})p^{2}+v_{2}(1−p^{2})(1−cos t))}dtdp + O(e^{−cϵb})}{∫_{−ϵ}^{ϵ} ∫_{−δ}^{δ} e^{−b((3v_{1}+v_{2})p^{2}+v_{2}(1−p^{2})(1−cos t))}dtdp + O(e^{−cϵb})}.

Now we change variables in both integrals, x = \sqrt{b(3v_{1} + v_{2})}p, s = \sqrt{bv_{2}}t and obtain

∫_{−ϵ}^{ϵ} ∫_{−δ}^{δ} F (p, t)e^{−b((3v_{1}+v_{2})p^{2}+v_{2}(1−p^{2})(1−cos t))}dtdp = \frac{A}{b}F (0, 0) + O(\frac{1}{b^{\frac{3}{2}}} )

with A the same constant in both the numerator and the denominator ( A does depend on v_{2} and 3v_{1} + v_{2}). We pass to the limit b → ∞ and obtain (77). This calculation implies the asymptotic solution of (76):

v_{1} = \frac{1}{6}, v_{2} = \frac{1}{2}  (78)

if 3v_{1} + v_{2} > 0, v_{2} > 0. The asymptotic solution belongs to the region.

Case II. If v_{2} > 0 but 3v_{1} + v_{2} < 0 then we multiply both numerator and denominator by e^{2bv_{1}}. We have to study thus the limit of the ratio of the integral

∫_{−1}^{1} (∫_{0}^{2π} F (p, t)e^{b(1−p^{2})(3v_{1}+2v_{2} cos t)}dt) dp

and the same integral with F replaced by 1. Now fix ϵ > 0. The contributions

∫_{−1+ϵ}^{1−ϵ} (∫_{0}^{2π} F (p, t)e^{b(1−p^{2})(3v_{1}+2v_{2} cos t)}dt) dp

in both numerator and denominator are exponentially small, 0(e^{−cϵb}) with c uniform for (v_{1}, v_{2}) ∈ [−\frac{1}{3}, \frac{2}{3}] × [0, \frac{1}{2}]. In the region where 1 − p ∈ [0, ϵ] we change variables x = Ab(1−p) with A > 0, A = −2(3v_{1} +2v_{2}). The expressions there become

\frac{1}{Ab} ∫_{0}^{2π} ∫_{0}^{Abϵ} F (1 − \frac{x}{Ab}, t)e^{−(1+ \frac{4v_{2}(1−cos t)}{A})x(1− \frac{x}{2Ab})}dxdt.

Using a similar change of variables for p near −1 we obtain

[F] (bv_{1}, bv_{2}) = \frac{∫_{0}^{2π} \frac{(F (1,t)+F (−1,t))dt}{1+ \frac{4v_{2}}{A} (1−cos t)} + O(\frac{1}{b})}{∫_{0}^{2π} _{\frac{2dt}{1+ \frac{4v_{2}}{A} (1−cos t)}} + O(\frac{1}{b})}

and therefore, if 3v_{1} + v_{2} < 0 and v_{2} > 0 we get the nontrivial limit

lim_{b→∞}[F](bv_{1}, bv_{2}) = \begin{matrix}∫ & {2π \atop 0} & \frac{(F (1,t)+F (−1,t))dt}{1+ \frac{4v_{2}}{A} (1−cos t)} \\ ∫ & {2π \atop 0} & \frac{2dt}{1+ \frac{4v_{2}}{A} (1−cos t)} \end{matrix}  (79)

with A = 2 |3v_{1} + 2v_{2}|. Substituting F = y_{1} = 1 − 3p^{2} we obtain −2; substituting F = y_{2} = (1 − p^{2}) cos t we obtain zero. So, the asymptotic solution of (76) is

v_{1} = −\frac{1}{3}, v_{2} = 0  (80)

if 3v_{1} + v_{2} < 0 and v_{2} > 0. In this case the asymptotic solution does not belong to the region, but to its boundary.

Case III. If 3v_{1} + v_{2} = 0 but v_{2} > 0 then, after multiplying both numerator and denominator by e^{2bv_{1}} we arrive at the ratio of integrals of the form

∫_{−1}^{1} ∫_{0}^{2π} F (p, t)e^{−b(1−p^{2})(1−cos t)}dtdp.

The integrals are dominated by the behavior at (p, cos t) = (±1, 1), and we deduce the asymptotics

lim_{b→∞}[F](bv) = \frac{1}{2}(F (−1, 0) + F (1, 0)).  (81)

Substituting F = y_{1} and F = y_{2} we deduce that

v_{1} = −\frac{1}{3}, v_{2} = 0  (82)

is the asymptotic solution of (76) if the parameters v obey 3v_{1} + v_{2} = 0, v_{2} > 0. We note that the asymptotic solution does not belong to the region, and not even to its boundary.

Case IV. If v_{2} = 0 and v_{1} > 0, then the exponent is bv_{1}(1 − 3p^{2}), which, after amplification by e^{−bv_{1}} leads to ratios of integrals

∫_{−1}^{1} ∫_{0}^{2π} F (p, t)e^{−3bv_{1}p^{2}}dpdt.

The limit in this case is

lim_{b→∞}[F](bv) = \frac{1}{2π} ∫_{0}^{2π} F (0, t)dt  (83)

and substituting we obtain the asymptotic solution of (76)

v_{1} = \frac{1}{6}, v_{2} = 0  (84)

in the case v_{1} > 0, v_{2} = 0. The asymptotic solution belongs to the region. Case V. Finally, if v_{2} = 0 and v_{1} < 0 we amplify by e^{2bv_{1}} and deduce

lim_{b→∞}[F](bv) = \frac{1}{4π} ∫_{0}^{2π} (F (1, t) + F (−1, t)) dt.  (85)

Substituting F = y_{1} and F = y_{2} we obtain the asymptotic solution of (76)

v_{1} = −\frac{1}{3}, v_{2} = 0  (86)

if v_{1} < 0 and v_{2} = 0. The asymptotic solution belongs to the region. This concludes the proof of the theorem.

Remarks. In each of the regions R_{1}, R_{2}, R_{3}, R_{4}, R_{5} above, the limit

lim_{b→∞}[F](bv) = ∫ F (p, t)dμ

exists and is given by a probability measure dμ. It is easily seen that these limits are attained uniformly on compacts L, (v ∈ L) in each region. The probability measures dμ depend on the region but are the same for all v in the region, are concentrated on the boundary of the parameter set K and ∫ F (p, t)dμ are given by the right hand sides of (77, 79, 81, 83, 85). Correspondingly, in each compact \frac{1}{6}[y_{1}](bv) and \frac{1}{2}[y_{2}](bv) converge to the stated constant values when b → ∞. For instance, if v ∈ L a compact subset of R_{1} then \frac{1}{6}[y_{1}](bv) converges to \frac{1}{6} and \frac{1}{2}[y_{2}](bv) converges to \frac{1}{2} when b → ∞. If the compact L does not contain the point (\frac{1}{6}, \frac{1}{2}) then there are no solutions of the simultaneous equations (76) for b sufficiently large (how large depends on L). Note that there are only three qualitatively different behaviors: all eigenvalues equal to zero, two eigenvalues equal to −\frac{1}{3} with the third eigenvalue equaling \frac{2}{3} in this case, and two eigenvalues equal to \frac{1}{6} with the third eigenvalue equaling −\frac{1}{3}. For the case of eigenvalues (−\frac{1}{3}, −\frac{1}{3}, \frac{2}{3}) the corresponding asymptotic steady state is a delta function concentrated at a fixed direction on the unit sphere, the prolate nematic state. For the case of eigenvalues (\frac{1}{6}, \frac{1}{6}, −\frac{1}{3}), the asymptotic steady state is uniform measure concentrated on the equator, the oblate nematic state. The fact that the uniform state is the unique steady state for low concentrations can be strengthened to a dynamical stability statement. At large concentrations this uniform state is isolated, and other states are present: this suggests that the uniform state is nonlinearly dynamically unstable. We have not proved this fact in this paper.

One can get a more precise description by expanding the asymptotic analysis. For instance, if v ∈ L ⊂ R_{1} with L compact, then one can verify that

\frac{1}{6}[y_{1}] = \frac{1}{6} − \frac{3C_{1}}{6(3v_{1} + v_{2})b} + O_{2,1}(v, b)  (87)

where the constant C_{1} is independent of v, b and is given by C_{1} = \frac{∫_{−∞}^{∞} x^{2}e^{−x^{2}} dx}{∫_{−∞}^{∞} e^{−x^{2}}}.

The error term O_{2,1}(v, b) is small: there exists an absolute constant Γ_{2} > 0 so that

|O_{2,1}(v, b)| ≤ Γ_{2}b^{−2}  (88)

holds for all v ∈ L. For [y_{2}] we obtain

\frac{1}{2}[y_{2}] = \frac{1}{2} − \frac{C_{1}}{2(3v_{1} + v_{2})b} − \frac{C_{2}}{2v_{2}b} + O_{2,2}(v, b)  (89)

with the same constant C_{1} and with C_{2} independent of v and b. The remainder obeys:

|O_{2,2}(v, b)| ≤ Γ_{2}b^{−2}  (90)

for all v ∈ L. These relations are obtained using the Taylor expansions of the functions y_{1} = 1 − 3p^{2} and y_{2} = 1 − p^{2} − (1 − cos t) + p^{2}(1 − cos t) near p = 0, t = 0. Consequently, the asymptotic equations in R_{1} are

v_{1} = \frac{1}{6} − \frac{3C_{1}}{6(3v_{1} + v_{2})b}

and

v_{2} = \frac{1}{2} − \frac{C_{1}}{2(3v_{1} + v_{2})b} − \frac{C_{2}}{2v_{2}b}.

One can expand to higher order in these equations, resulting in higher order algebraic equations and smaller remainders, uniformly on compacts in R_{1}. The same can be done in any of the other regions, with similar kinds of asymptotic developments.

4 Conclusions

We have obtained asymptotic developments for the steady solutions ψ of the Smoluchowski equation (4). The steady solutions are parametrized by the intensity b of the interaction potential and by two real parameters describing the eigenvalues of a real, traceless symmetric 3 × 3 matrix S (14, 15). When the intensity b is small enough then the uniform solution ψ = \frac{1}{4π}, (S = 0) is the unique steady solution. At high intensities several steady solutions coexist. For very large b the eigenvalues of the matrices S are close to one of the three possibilities (0, 0, 0) (corresponding to the uniform state), (−\frac{1}{3}, −\frac{1}{3}, \frac{2}{3}) (corresponding to a state ψ concentrated on a single direction e ∈ S^{2}) and (\frac{1}{6}, \frac{1}{6}, −\frac{1}{3}) (corresponding to a state concentrated uniformly on a geodesic (big circle)).

Acknowledgments PC’s research partially supported by NSF DMS-0202531. IGK’s research partially supported by AFOSR (Dynamics and Control) and NSF (ITR). EST’s research partially supported by NSF DMS-0204794, the US CRDF under grant No. RM1-2343-MO-02 and by the S.M. Ulam Visiting Scholar Fellowship at the CNLS Los Alamos National Laboratory under D.O.E. contract W-7405-ENG-36. The comments of the referees are gratefully acknowledged.

References

[1] M. Brenner, P. Constantin, L. Kadanoff, A. Schenkel, S. Venkataramani, Diffusion, attraction and collapse, Nonlinearity 12 1071-1098, (1999).

[2] P. Constantin, T. Dupont, R. Goldstein, L. Kadanoff, M. Shelley, S. Zhou, Droplet breakup in a model of the Hele-Shaw cell, Phys. Rev. E 47 (6), 593-596 (1993).

[3] P. Constantin, I. Kevrekidis, E. S. Titi, Remarks on a Smoluchowski equation, Discrete and Continuous Dynamical Systems, to appear (2004).

[4] M. Doi, Molecular dynamics and rheological properties of concentrated solutions of rodlike polymers in isotropic and liquid crystalline phases, J. Polym. Sci., Polym. Phys. Ed. 19 229-243, (1981).

[5] V. Faraoni, M. Grosso, S. Crescitelli and P. L. Maffetone, The rigid rod model for nematic polymers: An analysis of the shear flow problem, J. Rheol. 43 (3) 829-843, (1999).

[6] G. Forest and Q. Wang, Monodomain response of finite-aspect-ratio macromolecules in shear and related linear flows, Rheologica Acta, 42 20-46, (2003).

[7] G. Forest, Q. Wang, R. Zhou, The weak shear kinetic phase diagram for nematic polymers, preprint May 2003.

[8] G. Forest, R. Zhou, Q. Wang, Scaling behavior of kinetic orientational distributions for dilute nematic polymers in weak shear, preprint July 17, 2003.

[9] S. Helgason, Differential Geometry, Lie Groups, and Symmetric Spaces Academic Press, London, 1978.

[10] R. M. Jendrejack, J. J. de Pablo and M. D. Graham, A method for multiscale simulation of flowing complex fluids, J. Non-Newtonian Fluid Mech. 108 123-142, (2002).

[11] R. Jordan, D. Kinderlehrer, F. Otto, The variational formulation of the Fokker-Planck equation, SIAM J. Math. Anal 29 (1998) 1-17.

[12] R. G. Larson, The Structure and Rheology of Complex Fluids, Oxford University Press, London, 1999.

[13] R. Larson and H. C. Öttinger, The effect of molecular elasticity on out-of-plane orientations in shearing flows of liquid crystalline polymers, Macromolecules 24 6270-6282, (1991).

[14] M. Laso and H. C. Öttinger, Calculation of viscoelastic flow using molecular models: the CONNFFESSIT approach., J. Non-Newtonian Fluid Mech. 47 1-20, (1993).

[15] L. Onsager, The effects of shape on the interaction of colloidal particles, Ann. N. Y. Acad. Sci 51, 627-659, (1949).

[16] H. C. Öttinger, Stochastic Phenomena in Polymeric Fluids, Tools and Examples for Developing Simulation Algorithms, Springer-Verlag, New York, 1996.

[17] F. Otto, The geometry of dissipative evolution equations: the porous medium equation, Comm. Partial Differential Equations, 26 101-174, (2001).

[18] C. Siettos, M. D. Graham and I. G. Kevrekidis, Coarse Brownian dynamics computations for nematic liquid crystals, J. Chem. Phys. 118 (22) 10149-10157, (2003).

[19] J. K. C. Suen, Y. L. Joo and R. C. Armstrong, Molecular orientation effects in viscoelasticity, Annu. Rev. Fluid Mech. 34 417-444, (2002).