Conditions under which a natural iterative method for calculating the orientation distribution of rodlike particles decreases the free energy at each step

Abstract

It is shown that under certain conditions the iterative method suggested by the calculus of variations for calculating a cylindrically symmetric orientation distribution function (ODF) for rodlike particles strictly decreases the free energy at each full step. This monotonic behavior has strong implications for convergence of the sequence, or a subsequence, of the calculated ODFs. The result is valid not only for the reference system of hard core particles with indistinguishable ends, but also for free energy functions of similar type. The effect of an applied field is also permitted. Since the behavior of the iteration is intrinsic to the general form of the free energy function it applies to rod mixtures and may extend to a wider class of free energy functions. Outside the conditions that guarantee a monotonically decreasing free energy, we find that the iteration can fail to converge.

---

1 Introduction

This paper concerns the spontaneous alignment of rodlike particles. Mean field studies in terms of excluded volume date back to Onsager,¹ and have been extended to include long-range interparticle interactions (as reviewed by Franco-Melgar, Haslam and Jackson²). In calculating the orientation distribution functions (ODFs), an iterative procedure proposed in 1984³ has proven widely effective on a case-by-case basis, e.g., ref. 4–11. The purpose of this paper is to identify the remarkably broad range of conditions under which this iteration is guaranteed to decrease the free energy at each step, which is a substantial indication for convergence. The presentation is divided into three major sections.

The Background section develops necessary concepts and notation, beginning with a review of the origins of the iterative procedure for finding the free energy minima and concluding with a review of the expansion of the particle interactions and orientation distributions in terms of Legendre polynomials. Whereas the discussion in ref. 3 and 12 focused on excluded volume effects, while also noting possible application to long range interactions, the present analysis takes a general point of view, including long-range interactions and applied field effects, in addition to volume exclusion.

The subsequent Results and discussion section shows that for physically reasonable potentials, the iteration for cylindrically symmetric ODFs strictly decreases the free energy at each full step (unless one is already at a converged solution, i.e., a fixed point of the iteration). While Herzfeld et al.³ proved that each iterative step gives a descent direction (meaning that the free energy will decrease for a sufficiently small move in that direction), they provided only empirical evidence of the convergence of full steps in particular instances, as have others who have since used the method, e.g., ref. 4–11. The main benefit of the present analysis is in specifying when overshoot of the full step can be excluded. It also identifies conditions under which the iteration will converge to a “nearby” local minimum of the free energy function.

2 Background

2.1 The model

We consider a model for the orientation distribution function f(Ω) of identical rodlike particles, where the orientations are given in terms of spherical coordinates Ω = (θ,ϕ) with 0 ≤ θ ≤ π and 0 ≤ ϕ ≤ 2π, and f is required to be a continuous function of (θ,ϕ) and 2π-periodic in ϕ.

An orientation distribution function (ODF) is, by definition, a continuous probability density function for the orientation of any given particle. Thus

f(Ω) ≥ 0 (1)

andFollowing the approaches of Cotter,^13–15 Gelbart and Barron¹⁶ and Gelbart and Gelbart,¹⁷ and allowing for the possibility of an applied field, the mean field free energy (per rod, in units of kT) depends on the orientation distribution according towhere A is a constant, V and W are continuous functions, the entropy contribution appears in the second term, the contribution from interaction with an applied field appears in the third term, and the contribution due to interparticle interactions appears in the last term (with the positive constant B dependent on the particle geometry and proportional to the particle volume fraction). The corresponding Euler–Lagrange equation that local minima of F will satisfy iswhere λ is the constant such that f(Ω) satisfies eqn (2).

Notice that the mean field interactions between rodlike particles depend only on the angle γ₁₂ in [0,π] between their two orientation vectors (i.e., between the vector from the origin to the point P₁ on the unit sphere given by Ω₁ = (θ₁,ϕ₁) and the vector from the origin to the point P₂ on the unit sphere given by Ω₂ = (θ₂,ϕ₂)). For future reference, we note that, according to the law of cosines on the triangle whose vertices are P₁, the origin, and P₂,

γ₁₂ = Arc cos[sin(θ₁)sin(θ₂)cos(ϕ₁ − ϕ₂) + cos(θ₁)cos(θ₂)] (5)

For the reference system of hard core particles with no distinctions between the two ends, V = 0 and W is given by:

W(γ₁₂) = sin(γ₁₂) (6)

Elsewhere, this W is sometimes referred to as the Onsager kernel. In this case, the isotropic ODF (f(Ω) ≡ 1) is always a solution to eqn (4). Note that sin(γ₁₂) is positive within (0,π) with a maximum at γ₁₂ = π/2 which corresponds to orthogonal particles. As the coefficient B in the interparticle interaction term gets larger, this term increasingly favors particle alignment where γ₁₂ occurs less frequently near π/2.

2.2 Properties of the solutions of the Euler–Lagrange equation

Other workers have shown that, even in the absence of an applied field (i.e., V = 0), eqn (4) has a wide variety of solutions. It is worthwhile to review those results in order to focus the rest of this paper on the solutions that are of primary physical interest, i.e., those that represent minima in the free energy.

For the hard core reference system, Kayser and Raveché¹² considered ODFs having cylindrical symmetry and found that there is a branch of anisotropic solutions of the Euler–Lagrange eqn (4) at B = 32/π (ρ = 4 in their normalization and notation) off of the isotropic solution (f(Ω) ≡ 1), and beyond that value of B the isotropic solution is no longer a local minimum of the free energy. Furthermore, ODFs from bifurcations off of the isotropic solution at larger values of B are also not local minima. Vollmer¹⁸ has done an extensive analysis of the general case where the kernel W(γ₁₂) in the particle interaction term is a function of cos(γ₁₂), while V = 0, and allowed for solutions of the Euler–Lagrange equation that are not necessarily cylindrically symmetric. In ref. 18 it is proven, in particular, that for W(γ₁₂) equal the Onsager kernel sin(γ₁₂), the solutions branching off the isotropic solution at B = 32/π are axisymmetric.

For W(γ₁₂) equal to the Maier–Saupe kernel , Liu, Zhang, and Zhang,¹⁹ and Fatkullin and Slastikov²⁰ proved that there is only one bifurcation from the isotropic solution, showed that the solutions of the Euler–Lagrange equation are all axisymmetric, and provided explicit formulas for the solutions. Fatkullin and Slastikov²⁰ also proved the same holds for the dipolar kernel W(γ₁₂) = −cos(γ₁₂).

For the coupled dipolar/Maier–Saupe kernel W(γ₁₂) = −σ cos(γ₁₂) − κ cos²(γ₁₂) (with positive constants σ and κ), Zhou, Wang, Wang and Forest²¹ proved that solutions of the Euler–Lagrange equation that are local minima of F are axisymmetric. They also showed there are solutions that are not axisymmetric (but these are not local minima).

For more general W, Ball²² proved existence of non-axisymmetric solutions under certain conditions, and noted that it is an open question whether there exist non-axisymmetric solutions of eqn (4) (for general W, and particularly for the Onsager kernel) that are local minima of F.

Yin, Zhang and Zhang²³ have conducted a wide ranging, and beautifully illustrated, numerical investigation of the solution landscape for the W kernels noted above, and found a rich variety of non-axisymmetric solutions of (4), but none were local minima of the free energy F. They point out that, as far as currently known, it is possible there could be a secondary bifurcation branch from an axisymmetric primary bifurcation that could lead to non-axisymmetric minima, however, all secondary bifurcation solutions that they observed were not minima.

2.3 Assumption of cylindrical symmetry

The upshot of the studies summarized in Section 2.2 is that, in the absence of an applied field, all known minima of F correspond to cylindrically symmetric ODFs. The same should be true in the presence of an applied field, as long as that field does not break cylindrical symmetry, i.e., as long as V(Ω) is independent of ϕ. Therefore, since the focus of this study is the behavior of the iterative method³ for finding solutions of the Euler–Lagrange eqn (4) that are minima of the free energy F, it will henceforth be assumed that V depends only on θ, and the orientation distribution functions f are cylindrically symmetric, i.e.,

f is independent of ϕ. (7)

In this case, the set S of orientation distribution functions will be all continuous functions on [0,π] satisfying

f(θ) ≥ 0 (8)

and the free energy function in eqn (3) will have the formwhere
Since cos(ϕ₁) is 2π periodic, the inner integral over ϕ₁ in [0,2π] is independent of ϕ₂, and K simplifies toNote that, by inspection, K is symmetric, i.e.,

K(θ₂,θ₁) = K(θ₁,θ₂) (12)

The assumption of cylindrical symmetry should be accompanied by recognition that there can be ODFs that are local minima of the free energy F relative to the set of axisymmetric ODFs, but are saddle points, i.e., have directions of decrease in F, when non-axisymmetric ODFs are included. Such cases will be discussed below.

It will be useful in certain instances to express an orientation distribution function f(θ) as a linear combination of Legendre polynomials P_m(cos θ), which is natural in the setting of spherical coordinates:

Since(see, for example, Chapter 11 of Weber and Arfken²⁴), it follows thatSimilarly, it will be useful to express the interparticle potential W(γ) aswhereA noteworthy feature of the Onsager kernel W(γ) = sin(γ) is that w_m ≤ 0 for all m > 0, and the sum of the w_m is absolutely convergent, i.e., is finite (see Appendix 1 in ESI†). This is also the case for the dipolar kernel W(γ) = −P₁(cos γ), the Maier–Saupe kernel W(γ) = −(2/3)P₂(cos γ), as well as the coupled dipolar/Maier–Saupe kernel. Generally, it will be seen below that w_m ≤ 0 for m > 0 plays an essential role in the proof that each iterative step strictly decreases the value of the free energy function (10) and (11) (unless one is already at a converged solution). The physical significance of this condition is illustrated in Fig. 1, which shows that each term w_mP_m(cos γ) for m > 0 with a negative coefficient w_m contributes to a minimum of W(γ) at γ = 0. This minimum favoring alignment counteracts the effect of entropy.

Fig. 1 Legendre polynomial plots. (a) −1 times the first three odd index Legendre polynomials. Each has a minimum at γ = 0 and a maximum at γ = π. Terms of this form in the energy potential W contribute to it being energetically favorable for particles to point in the same direction and unfavorable to point in opposite directions. (b) −1 times several even index Legendre polynomials. Each has a minimum at γ = 0 and also at γ = π. Terms of this form in the energy potential W contribute to it being energetically favorable for particles to point in the same direction or in opposite directions.

In view of these features of the −P_m(cos γ) functions, it will henceforth be assumed that W(γ) has an expansion in terms of Legendre polynomials that satisfies these specific conditions on the coefficients:

In the case of particles with no distinctions between the two ends, so that W(π − γ) = W(γ), we also have

K(π − θ₁,θ₂) = K(θ₁,π − θ₂) = K(θ₁,θ₂) (19)

because if Z is in [−1,1], then Arc cos(−Z) = π − Arc cos(Z), whence replacing ϕ by ϕ + π in the integrand of eqn (11) gives the result.

2.4 The iteration suggested by the calculus of variations

To set the stage for the proof in Results and discussion, and have a key equation in hand, we now formally derive the iterative method given in ref. 3 and 12, extended to include a term for the possibility of an applied field. For the rest of this section, assume the ODF f(θ) is a local minimum of the free energy function and positive on [0,π] (to avoid problems with ln f). Then for any continuous function g(θ) on [0,π] for which max |g(θ)| is sufficiently small andit will be the case that the function f(θ) + tg(θ) will be an ODF (i.e., satisfy eqn (8) and (9)) and a positive function on [0,π] for −1 ≤ t ≤ 1. Since it is being assumed that f is a local minimum of F, the derivative of F(f(θ) + tg(θ)) with respect to t must be 0 at t = 0. Since f(θ) + tg(θ) > 0, we can “differentiate under the integrals” in F, evaluate at t = 0, use (20), recall that K is symmetric, and find that, for all such functions g, f must satisfy:which suffices to show that the quantity in the square brackets must be a constant, i.e.,where the constant λ corresponds to the Lagrange multiplier for the normalization constraint .

From eqn (22) we have

This suggests (ref. 3 and 12) trying the well known successive substitution method: if fⁿ(θ) is a starting value satisfying eqn (8) and (9), or is the “current approximation” to a local minimum f of the free energy function, then the next iterate fⁿ⁺¹(θ) (the approximation succeeding fⁿ(θ)) is given by I(fⁿ)(θ) defined to be

fⁿ⁺¹(θ) = I(fⁿ)(θ) = N(fⁿ)(θ)/D(fⁿ) (24)

This is shorthand for eqn (16) in Herzfeld et al.³ for cylindrically symmetric ODFs.

Note that fⁿ⁺¹ and any function f that satisfies eqn (23) will be positive on [0,π].

Several special cases deserve particular attention: if V(π − θ) = V(θ) and the ODF f(θ) is a solution of eqn (23), then g(θ) defined to be f(π − θ) will also be a solution of eqn (23). This follows from the rotational invariance of γ₁₂ and can be verified directly by employing eqn (22):

and then using the change of variable η₂ = π − θ₂ in the integral over θ₂ in (25), and the definition of K(θ₁,θ₂) in eqn (11).

The same steps verify that if V(π − θ) = V(θ) and fⁿ(π − θ) = fⁿ(θ), then fⁿ⁺¹(π − θ) will equal fⁿ⁺¹(θ). Therefore, in cases where there may be a local minimum f of F that is not symmetric about θ = π/2, one should be sure to use initial ODFs f⁰ that are not symmetric about θ = π/2.

Similarly, if W(π − γ) = −W(γ) and there is no applied field (V(θ) = 0), when fⁿ(π − θ) = fⁿ(θ), will be 0, and fⁿ⁺¹ will be the isotropic solution. Therefore, in this case one should thus also use initial ODFs that are not symmetric about π/2.

In addition, when W(π − γ) = W(γ) and V(π − θ) = V(θ), by eqn (19), fⁿ⁺¹ and any function satisfying eqn (23) will satisfy f(π − θ) = f(θ).

2.5 The expansion of K(θ₁,θ₂) in terms of Legendre polynomials

This expansion was developed in Kayser and Raveché;¹² we summarize it here in terms of our notation and normalizations since it is a key part of the proof that the iteration decreases the free energy. As shown in Appendix 2 in ESI,† forand K as in eqn (11)Moreover, if [scr B, script letter B](f) is defined to be the last term (the particle interaction term) in eqn (10), andthen it is shown in Appendix 2 in ESI† thatand in Appendix 3 in ESI† it is shown that is finite. Note that for f to satisfy eqn (9), η₀ must be 1. Furthermore, if g is the difference between two ODFs, then η₀ for g is 0. Notice also that, if all the w_m are ≤0 for m > 0, then [scr B, script letter B](g) ≤ 0. This plays a key role in the proof that the iteration decreases the free energy.

The result (27) together with the fact that the Legendre polynomials are orthogonal (14) shows that, when V = 0, the isotropic ODF f = 1 is always a solution of eqn (22) and (23).

When V = 0, the values of B that are potential locations of bifurcations of anisotropic solutions of the Euler–Lagrange equation off of the isotropic solution are B_m = (2m + 1)/(−w_m) (see Appendix 5 in ESI†). The smallest of these B_m, denoted by B*, is the location at which one would generally expect to have a branch of anisotropic local minima of the free energy. If there is an odd number of the B_m equal to B*, then there will be such a branch at B* (Vollmer,¹⁸ Rabinowitz²⁵).

3 Results and discussion

3.1 Critical relationships

To make the presentation of new results somewhat self-contained, this section repeats a few primary equations.

Let S denote the set of orientation distribution functions (ODFs), i.e., functions f(θ) that are continuous on [0,π] and satisfy

f(θ) ≥ 0 (30)

Assume that the conditions in eqn (18) on W are satisfied and the free energy function F defined on functions f in S is given by eqn (10) and (11), where A is a constant and B is a positive constant.

As noted after eqn (11), K(θ₂,θ₁) = K(θ₁,θ₂). By eqn (19), when W(π − γ) = W(γ), K(π − θ₁,θ₂) = K(θ₁,π − θ₂) = K(θ₁,θ₂).

The iterative method (as in Kayser and Raveché¹² and Herzfeld et al.³) takes a function fⁿ(θ) in S and produces the next iterate fⁿ⁺¹(θ) = I(fⁿ)(θ) defined in terms of our normalization and notation by:

fⁿ⁺¹(θ) = I(fⁿ)(θ) = N(fⁿ)(θ)/D(fⁿ) (32)

where the operators N and D are given byThe alternative formulation of eqn (32)where the constant λ equals ln(D(fⁿ)), will be useful in the proof that the iteration decreases the free energy when the condition (18) holds.

Note that if f⁰(θ) is any function in S (actually any continuous (or L¹([0,π])) function), then, directly from its definition and the properties of K noted above, f¹(θ) = I(f⁰)(θ) will be in S and in fact will be positive on [0,π]. In addition, if W(π − γ) = W(γ) and V(π − θ) = V(θ), then f¹(π − θ) = f¹(θ).

Now assume f⁰(θ) is in S, and for each non-negative integer n, define fⁿ⁺¹(θ) to be I(fⁿ)(θ). For a given n, define

g(θ) = fⁿ⁺¹(θ) − fⁿ(θ) and fⁿ(t)(θ) = fⁿ(θ) + tg(θ) for t ∈ [0,1] (36)

Note for each positive integer n, fⁿ(t)(θ) will be in S and will be positive for each t in [0,1].

3.2 The main result

3.2.1 Statement of the main result. When the conditions in eqn (18) on W are satisfied, for each non-negative integer n either

fⁿ⁺¹(θ) = fⁿ(θ) (37)

or

F(fⁿ(t)) is a strictly decreasing function of t for t ∈ [0,1]; in particular, F(fⁿ⁺¹) < F(fⁿ). (38)

An immediate consequence of this result is that, if f(θ) in S is a local minimum (over functions in S) of F, then I(f) = f, and hence f(θ) > 0 on [0,π].

Furthermore, from the result that F(fⁿ(t)) is a strictly decreasing function of t for n = 0, 1, 2, 3,… (unless fⁿ⁺¹ = fⁿ for some n) we will prove below that a subsequence of {fⁿ(θ)} converges uniformly to a function f(θ) in S for which I(f) = f, and, with some further assumptions, that the entire sequence converges to f.

3.2.2 Proof of the main result. Let n be a non-negative integer. Then fⁿ and fⁿ⁺¹ are orientation distribution functions, i.e., in S, and fⁿ⁺¹ is positive for θ in [0,π]. Suppose fⁿ⁺¹ ≠ fⁿ (otherwise there is nothing left to prove). We will show that (see the definitions in eqn (36)) the derivative with respect to t of F(fⁿ(t)(θ)) is <0 for 0 < t < 1. Then the mean value theorem of calculus (if y(x) is a continuous function on the finite closed interval [a,b] and y has a finite derivative y′(x) for each x in the open interval (a,b), then there is an x̂ in (a,b) for which y(b) − y(a) = (b − a) × y′(x̂)) shows F(fⁿ(t)) is strictly decreasing for 0 ≤ t ≤ 1 completing the proof.

In what follows, we write down the derivative; use the fact that (since fⁿ and fⁿ⁺¹ both satisfy eqn (31)); utilize eqn (35) as was done in eqn (17) of ref. 3 where it was shown that the derivative of F(fⁿ(t)) is <0 at t = 0 (as long as fⁿ(θ) > 0 on [0,π]); and employ eqn (29) (with f replaced by g), and the sentence below it to show that the contribution to the derivative of F(fⁿ(t)) from the term involving [scr B, script letter B](g) is ≤0 so it needs no further consideration. This is the place where the signs of the {w_m} in the expansion of W_m come into play.

The first step is to write out explicitly what F(fⁿ(t)) is, and then take its derivative with respect to t. For this we can take the derivative under the integral signs, since fⁿ(t)(θ) is bounded away from 0 for any given t in (0,1] so x ln(x) not being differentiable at x = 0 does not cause any problems.

Recalling that g(θ) = fⁿ⁺¹(θ) − fⁿ(θ) and fⁿ(t) = fⁿ + tg, for t in [0,1] we have

Then for t in (0,1] the derivative is:Now recall that . Since K is symmetric (K(θ₂,θ₁) = K(θ₁,θ₂)) the third and fourth terms on the right side of eqn (40) are equal. The last term on the right side of eqn (40) equals 2t[scr B, script letter B](g), which from eqn (29) and the sentence below it, is ≤0, so it needs no further treatment. We thus haveNow from the key eqn (35), we haveUsing this in the last term in eqn (41) we see thatSince λ is a constant, the integral of λg sin(θ) is 0. Denote the right hand side of (43) by [scr R, script letter R](t). Recalling that g = fⁿ⁺¹ − fⁿ, we then haveSince fⁿ⁺¹(θ) > 0 on [0,π], [scr R, script letter R](t) is well defined for t in (0,1] (and also for t = 0 if fⁿ > 0 (which is always the case for n > 0), though we don't need [scr R, script letter R] defined at t = 0 for our proof).

Since ln x is a strictly increasing function, and here fⁿ⁺¹ ≠ fⁿ; as was seen in eqn (17) of Herzfeld et al.,³ [scr R, script letter R](0) < 0 (as long as fⁿ > 0). What really matters for our proof is that, by inspection, we have the striking result that [scr R, script letter R](1) = 0. Hence if we could show that [scr R, script letter R](t) was a strictly increasing function for t in (0,1], we would have shown that [scr R, script letter R](t) < 0 for t in (0,1), and thereby dF(fⁿ(t))/dt is <0 for t in (0,1) and (by the mean value theorem) the proof of the main result would be done. So let's look at the derivative of [scr R, script letter R] for t in (0,1):

since here fⁿ⁺¹ ≠ fⁿ and fⁿ⁺¹ > 0 and fⁿ ≥ 0 and t is in (0,1) and hence the denominator which equals tfⁿ⁺¹ + (1 − t)fⁿ is positive, and so indeed [scr R, script letter R](t) < 0 for t in (0,1), and the proof is complete.

This result, and the proof, seem to be intrinsic to the form of the free energy function, and the corresponding form of the iterative method suggested by the calculus of variations. As will be seen in the next subsection, F(fⁿ(t)) being strictly decreasing for t in [0,1] has significant implications for convergence of the iteration. Note that the convergence results below depend on the form of the iteration (22)–(24), not just that F(fⁿ(t)) is strictly decreasing (unless the iteration has converged). It may be the case that the approach used here to show that the iteration decreases the free energy is applicable in other settings.

3.2.3 A note on polydisperse systems. If one is considering a polydisperse system as in eqn (18) in ref. 3, and one carries out the iteration updating the ODF for one particle type at a time, then under the same conditions on each “diagonal” W^σσ as used here for W, the iteration at each step will decrease the free energy (unless the iteration left the ODF for that particle type fixed); the proof is essentially the same, but with more algebra.

3.3 Convergence consequences of F(fⁿ(t)) being strictly decreasing (unless fⁿ⁺¹ = fⁿ)

3.3.1 An infinite subsequence of the iterates {fⁿ} will have a uniformly convergent subsequence. If fⁿ⁺¹ = fⁿ for some n, then trivially, {fⁿ} converges to some f in the set S of ODFs, and I(f) = f, so we don't need to further address that case. We next show that any infinite subsequence of {fⁿ} is conditionally compact meaning it will have a subsequence that converges uniformly to some function f that by continuity will be in S. This follows in general for the type of integration involving K(θ₁,θ₂) in the definition of the iteration. Specifically, since K(θ₁,θ₂) and V(θ) are continuous (and hence uniformly continuous and bounded) for θ₁ and θ₂ and θ in [0,π], and fⁿ is non-negative on [0,π] and for each n, the magnitude of the argument of the exponential in eqn (33) is uniformly bounded by (B max(|K(θ₁,θ₂)|) + max(|V(θ)|)), so N(fⁿ)(θ) is uniformly bounded (i.e., by the same constant for all n), and, importantly, D(fⁿ) is uniformly bounded away from 0, so {fⁿ} is uniformly bounded. Since K and V are uniformly continuous, {fⁿ} is equicontinuous, meaning given ε > 0 there is a δ > 0 such that |fⁿ(θ₁) − fⁿ(θ₂)| < ε whenever |θ₁ − θ₂| < δ and this holds for every n and every θ₁, θ₂ in [0,π] for which |θ₁ − θ₂| < δ. The Arzelà–Ascoli theorem (e.g., page 158 in Rudin²⁶) then applies and every infinite subsequence of {fⁿ} is conditionally compact (i.e., has a subsequence that converges uniformly to a continuous function on [0,π] which will be in S since each fⁿ is).

The same reasoning shows that the set of functions {I(s)} for s ranging over the set of ODFs S is conditionally compact.

3.3.2 The first convergence result. Statement: let {f^n_i} be a uniformly convergent subsequence of {fⁿ} (at least 1 exists from the discussion above), converging to the function f(θ) as i → ∞. Then I(f) = f.

Proof: if not, then from eqn (38) in the main result, we have F(I(f)) < F(f) which will lead to a contradiction. Define f(t)(θ) as follows. Let t = n + r where n is a non-negative integer and 0 ≤ r ≤ 1. Then f(t)(θ) is defined as fⁿ(r)(θ), as in eqn (36). This continuously “patches together” the fⁿ(θ) functions, and by the main result, F(f(t)) is strictly decreasing for t in [0,∞) (unless for some n, fⁿ⁺¹ = fⁿ, but then there would be nothing left to prove). From the discussion in the first paragraph of Section 3.3.1, {fⁿ} is uniformly bounded, and so we have F(f(t)) is uniformly bounded over t, and so has a finite lower bound. Hence F(f(t)) converges monotonically “down to” a finite limit which we denote by β (which is the greatest lower bound for the values F(f(t))).

Now F is a continuous function on S (with ‖f₁ − f₂‖ = max(|f₁(θ) − f₂(θ)|)), so F(f^n_i) must converge to F(f), which from the previous paragraph equals β. Since D(fⁿ) is uniformly bounded away from 0, I(f^n_i) must converge to I(f) and then F(I(f^n_i)) must converge to F(I(f)). But by definition, I(f^n_i) = f^n_i+1 so F(I(f^n_i)) must be ≥β. Hence F(I(f)) ≥ β, contradicting F(I(f)) < F(f).

Discussion: note in general, a subsequence {f^n_i} from an iteration fⁿ⁺¹ = I(fⁿ) might possibly converge to a solution of the calculus of variations eqn (22), or equivalently, (23) at which F has one or more directions of decrease, rather than a local minimum. However, in practice, that would be “unlikely” since the iteration strictly decreases F (unless fⁿ⁺¹ = fⁿ for some n). Similarly, unless f⁰ is a local maximum, a subsequence {f^n_i} can not converge to a local maximum, since f(t) is strictly decreasing (unless fⁿ⁺¹ = fⁿ for some n).

We will see that with strong but natural assumptions, the full sequence {fⁿ} will converge to a solution of I(f) = f.

The result below suggests (but does not prove) that if one starts with f⁰ “close enough” to an isolated local minimum f of F, then the full sequence {fⁿ} will converge to f. Note Kayser and Raveché¹² and Scheurle²⁷ (which is ref. 6 in ref. 12) examine convergence using the spectral radius of the Fréchet derivative of the mapping under consideration which here is I(f), and bifurcation theory. Herzfeld, Berger and Wingate³ observed that computational results indicated that the mapping f → I(f) was Lipschitz continuous with a Lipschitz constant L below 1 near solutions of I(f) = f for the cases for which computations were carried out, with L (perforce) approaching 1 near bifurcation points.

What one would “like” (and might expect) to be the case (but is not proven here) is that one would often have the situation of an isolated local minimum f of F contained in an open set U on whose boundary ∂S, F is ≥ some constant k which is greater than F(f), and f is the only solution of I(f) = f in U. Note I have not yet specified what norm is to be used to define “open” and “boundary” and “local minimum”. Then if f⁰ is in U and F(f⁰) < k, one would have the full sequence {fⁿ} converging to f, because these assumptions and the main result would show that f(t) can't “escape” U (since F(f(t)) is decreasing but would have to increase to “cross” the boundary of U), and the entire sequence fⁿ would have to converge to f by the assumption that f was the only fixed point of I in U.

For this purpose, the appropriate (and natural) function space and norm to use is L²([0,π],dμ = sin θdθ), i.e., functions z for which the norm ‖z‖₂ defined to be is finite. The maximum norm (‖z‖ = maximum over θ in [0,π] of |z(θ)|) on the set C[0,π] of continuous functions on [0,π] is too strong in that one can use the same type of construction as in Appendix 4 in ESI† with hat functions to show that if f was a local minimum of F then for any ε > 0 and ν > 0 there is a function f_ν in S with |F(f_ν) − F(f)| < ε and max |f_ν(θ) − f(θ)| = ν. The idea is that with narrow hat functions, one can construct a “perturbation” of f satisfying these conditions that is still in S. Note that while max |f_ν(θ) − f(θ)| = ν; ‖f_ν − f‖₂ would be small if the “hat” functions were sufficiently narrow (had small support). This shows that the desired situation in the previous paragraph can not take place using the maximum norm to define “open” and “boundary”. Also note that if z is in C[0,π], then it will also be in L²([0,π],dμ = sin θdθ), and

3.3.3 The second convergence result. Definitions: U is an open set in S, if U is a subset of S and for each u in U, there is an ε > 0 (depending on u) such that if s is in S and ‖s − u‖₂ as defined in the previous paragraph is <ε then s is in U. The boundary of U, denoted by ∂U is the set of all functions b in S such that for any ε > 0 there are members u₁ and u₂ of S with ‖u₁ − b‖₂ < ε and ‖u₂ − b‖₂ < ε and u₁ is in U and u₂ is not in U (u₁ or u₂ can equal b) (this is the relative topology of L²([0,π],dμ = sin θdθ) on S).

A comment: let U be an open set in S and let β be the greatest lower bound for F(∂U). Note since closed bounded sets in L²([0,π],dμ = sin θdθ) are not necessarily compact, there might not be a function b in ∂U for which F(b) = β (e.g., the set of orthonormal Legendre polynomials [scr P, script letter P]_n(cos θ) defined in eqn (S.21) in Appendix 3 in ESI† has no convergent subsequence).

Statement: under the conditions and definitions given in the previous 2 paragraphs; if f⁰ is in U and F(f⁰) < β it will be the case that the limit f of any uniformly convergent subsequence {f^n_i} of fⁿ (starting with f⁰) is in U (from the results above, there is at least one subsequence converging uniformly to some function f̄ in S, i.e., ‖f^n_i − f̄‖ converges to 0).

Proof: using the notation from above, let f^n_i be a convergent subsequence of fⁿ converging to some function f̄ (which from the first convergence result satisfies I(f̄) = f̄). We need to show f̄ is in U. Consider f(t), and F(f(t)) which is decreasing (actually strictly decreasing unless some fⁿ⁺¹ = fⁿ). The idea is that the assumptions imply that U is “an energy well” from which f(t) can't “escape”, since to do so, f(t) would have to “cross” the boundary of U, where F ≥ β.

Suppose f̄ is not in U. Now since F(f(t)) is continuous and decreasing and we assumed F(f(0)) < β it can not be the case that f̄ is in ∂U where F ≥ β. Let T be the least upper bound of the t values for which f(t) is in U for t in [0,t]. Since U is open, by continuity and eqn (46)T > 0. We next need to show T is finite. Since we are working with the assumption that f̄ is not in U, and we just noted that f̄ can't be in ∂U, there must be an ε > 0 such that if u satisfies ‖u − f̄‖₂ < ε then u is not in U (otherwise we would have f̄ in the boundary of U). However, we have f(n_i) = f^n_i converging uniformly to f̄, so (by eqn (46)) T must be finite. By the definition of T, and the same type of reasoning used above, f(T) must be in the boundary of U, which is a contradiction to F(f(t)) < β (if f(T) is not in U and not in ∂U then T would not be the least upper bound; if f(T) is in U then T would not be an upper bound). Thus we have shown f̄ is in U.

3.3.4 Corollary of the second convergence result. Statement: under the assumptions of the second convergence result, if f⁰ is in U and F(f⁰) < β, and there is at most 1 solution of I(f) = f in U, then the full sequence fⁿ converges uniformly to a function f in U for which F(f) < F(u) for all u ≠ f in U.

Proof: from the second convergence result there is at least 1 subsequence converging uniformly to some function f and the limit of any uniformly convergent subsequence is in U (and from the first convergence result the limit satisfies I(f) = f). If the full sequence doesn't converge to f, then one could extract a second subsequence converging to another function f₂ in U satisfying I(f₂) = f₂, contradicting the assumption that there was only one solution of I(f) = f in U. Suppose there was some function u ≠ f in U for which F(u) ≤ F(f). Then we could set f⁰ = u, and by the main result, either I(u) = u; or F(I(u)) < F(u) and there would be a convergent subsequence of the iterates starting with u converging to some function in U at which the free energy was <F(f), and so either way we would get the contradiction of a second solution of I(f) = f in U.

3.4 The behavior of the iteration

3.4.1 Numerical implementation. In Herzfeld et al.,³W and ODFs were represented by their values on evenly spaced grids of θ and ϕ values, and integrals were approximated by trapezoidal numerical integration, which is well suited for use with periodic functions. This allowed for readily estimating the accuracy of the results by successively doubling the number of intervals and comparing resulting calculated values of interest. R code for calculating ODFs this way, extended to more general W and an applied field, is available in the GitHub repository,²⁸ which also contains code for sample runs, with corresponding results. This approach has been found to be useful in various settings, e.g., ref. 4–11.

Representative results are illustrated in the next section. Unless otherwise noted, the iterative method was applied using trapezoidal numerical integration with 512 evenly spaced θ intervals on [0,π] and 2048 evenly spaced ϕ intervals on [0,2π]. Initial ODFs used in calculations for the hard core reference system were either “axial” (peaks at θ = 0 and θ = π and a trough at θ = π/2), or “planar” (peak at θ = π/2 and troughs at θ = 0 and θ = π). (Further details are available in ref. 28.) As noted below eqn (25), when there may be local minima of F that are not symmetric about π/2, one should include initial ODFs for the iteration that are not symmetric about π/2.

Our computational procedure avoids the approximations made in previous work. Earlier approaches for obtaining minima of the free energy (in the hard particle system where W(γ) = sin(γ) and V = 0) have included replacing the ODF f(θ) either by a finite expansion in terms of Legendre polynomials P_m(cos θ) with even index m, and determining the coefficients that minimize the free energy (Lasher²⁹), or by a function with just one optimizable parameter (Onsager,¹ Baron and Gelbart³⁰). Another approach replaced the particle interaction term in the free energy with the expansion (29) and used the resulting form of the Euler–Lagrange eqn (22) in terms of the unknown {η_m} (note η₀ must be 1 for f to be an ODF):

where λ is the Lagrange multiplier determined by the constraint , together with the consistency requirementAlben³¹ and Cotter¹³ used these conditions with M = 2 to solve for the one unknown η₂ (η₀ = 1, and here η₁ = 0 by symmetry when W(γ) = sin(γ)). Kayser and Raveché¹² used (47) and (48) to iteratively determine the {η_m} in a fashion analogous to (24), with M taken large enough to obtain good accuracy (for W = sin(γ), all the odd index η_m are 0). Equivalent forms of (47) and (48) were used by Lekkerkerker et al.³² to obtain minima.

The iteration method studied here is complementary to the uniform sampling method of Yin, Zhang, Zhang.²³ Both procedures will accurately determine ODFs that minimize the free energy. By its nature, the iterative method will in general not find ODFs that are solutions of the Euler–Lagrange eqn (22) but are not local minima of the free energy (such as were obtained by Kayser and Raveché¹² using a different method). On the other hand, Yin, Zhang and Zhang²³ give an extensive discussion of methods for exploring the landscape as B is varied of solutions of (22) that are not necessarily local minima of the free energy. In particular, their own uniform sampling method, while necessarily more complex than the iterative method, has the advantages that it searches over all ODFs (not just axisymmetric ODFs), and by design enables construction of the landscape and bifurcation branches of solutions of (22).

The iterative method studied here naturally converges to minima of the free energy as considered over the subset of axisymmetric ODFs. Since, as far as is known at present, minima of the free energy for models having the form of eqn (3) where W is a function of cos(γ₁₂) are axisymmetric,^18–23 the iterative method provides a robust straightforward way of determining minima of the free energy, i.e., the physically observed ODFs, only requiring the user to run the algorithm with a variety of starting values for the iteration to best determine the global minimum. As previously noted below eqn (12), if f is a saddle point but all directions of decrease in the free energy around f are non-axisymmetric ODFs, then f will be a local minimum over the set of axisymmetric ODFs and the iterative method will converge to it when started with a suitable ODF.

3.4.2 Illustrative examples. The figures illustrate the behavior of the iteration and display examples of ODFs that are local minima for a variety of particle interaction kernels W. These include cases where W(γ) is not symmetric about π/2 (Fig. 7c and 8c), and three cases in which there is an applied field (Fig. 8). An example is also given of failure of the iteration when the conditions on the {w_n} being ≤0 for n > 0 are not satisfied (Fig. 9). Although convergence is still possible in the presence of a positive w_n with n > 0, other procedures should be used for finding minima or solutions of the Euler–Lagrange equation in such circumstances.²³

Fig. 2 displays the ODFs that are local minima of the free energy for the hard core reference system (W(γ) = sin(γ) and the applied field V = 0) with B = 14. Graphs of the axial and the planar minima for other values of B are given in ref. 3. The axial ODF is the global minimum of the free energy. The planar ODF is a local minimum with respect to the set of axisymmetric ODFs, but is a saddle point having a two dimensional set of directions of decrease over all ODFs (see Fig. 7 in Yin, Zhang and Zhang²³ where they use the terminology oblate for planar ODFs and prolate for axial ODFs).

Fig. 2 Axial and planar orientation distribution functions that are local minima in the hard core reference system with B = 14. The particle interaction W(γ) = sin(γ) favors particle alignment (γ = 0 or π). The planar local minimum ODF balances unfavorable particle interactions (at θ around π/2) with higher entropy from a more dispersed equatorial orientation distribution. Note that relative particle populations at each value of θ are obtained by multiplying the ODF by sin(θ) to account for the axially symmetric distribution in ϕ.

Fig. 3 illustrates that when all the coefficients w_n for n > 0 in the expansion of the interparticle potential in terms of Legendre polynomials are ≤0, each iterative step decreases the free energy. This has been the case for all calculations that we have carried out (by trapezoidal numerical integration), up to computer finite precision.

Fig. 3 Free energy F vs. iteration number for axial and planar initial ODFs in the hard core reference system with B = 14. Iterations were done until the convergence criterion was satisfied. The free energy strictly decreased at each iteration. The ODF at iteration 0 is the one chosen to start the iteration.

Fig. 4 displays the dependence of particle order on the value of B in the hard core reference system (W(γ) = sin(γ) and V(θ) = 0) as obtained with initial axial (peaks at θ = 0 and θ = π) and planar (peak at θ = π/2) ODFs. The jump at approximately B = 8.883 is indicative of the known first order phase transition at this value of B, e.g., ref. 12 and 33–35.

Fig. 4 Ordering in the hard core reference system. (a) The entropy 〈−ln f(θ)〉 vs. B for converged ODFs resulting from starting the iteration with axial (purple) and planar (green) ODFs. . The jump in this plot for the axial ODFs around B = 8.883 (delineated by dots) corresponds to the location ρ_b of Kayser and Ravaché.¹² (Their ρ* = 4 corresponds to B = 32/π = 10.18592 (marked here with a short vertical black bar) where planar local minima bifurcate from the isotropic ODF.) (b) The order parameter 〈P₂(cos θ)〉 vs. B for converged ODFs, resulting from initial axial (purple) and planar (green) ODFs. . Note the discontinuity around B = 8.883 (delineated by dots) and change in slope around B = 32/π = 10.18592.

Fig. 5 plots the free energy as a function of B for the axial, planar and isotropic solutions of the Euler–Lagrange equation for the hard core reference system. Fig. 5a shows the free energy per particle for all the local minima. Fig. 5b shows the free energy per unit volume for the global minimum, as required to identify the boundaries for coexistence of the isotropic and planar phases by the common tangent construction. The number of θ intervals used here was increased to 1024, as the determination of the coexistence region is quite sensitive to the accuracy of the numerical approximation.³⁴ The coexistence region depicted in Fig. 5b is consistent with results obtained by appropriate thermodynamic calculations in ref. 12, 32, 34 and 36. (Note the scale factor 8/π to go from endpoints of the coexistence region determined in these references to those depicted in Fig. 5b.)

Fig. 5 Free energy vs. B for the hard core reference system. (a) The free energy per particle, F, for the isotropic ODF is compared with F for local minima found using the iterative method, eqn (24). Note that the isotropic ODF is no longer a local minimum of F for B larger than the bifurcation value 32/π. The “unstable” bifurcation branch of ODFs that are not local minima, but lead to the axial ODFs that are local minima,¹² will not be found using the present iteration procedure. (b) The common tangent construction in which the free energy per particle, F, is multiplied by B to give an energy per unit volume, with the zero of energy per particle adjusted to give a horizontal common tangent by choosing C = 6.4387.

Fig. 6 plots results for the Maier–Saupe kernel W(γ) = 1/3 − cos²(γ) = −(2/3)P₂(cos γ) = sin²(γ) − 2/3 with B = 9. The bifurcation point from the isotropic solution is at B = 7.5,^19,20 Appendix 5 in ESI.† Note that some texts omit constants in the kernel because they have no effect on which ODFs are minima of the free energy or solutions of the Euler Lagrange equation. Yin, Zhang and Zhang²³ report (see their Fig. 2) that for the Maier–Saupe kernel, the planar ODF when viewed in the set of all ODFs (not just the axisymmetric ODFs) is a 2-saddle, meaning there is a two dimensional set of descent directions at this ODF. With B = 9, the order parameter r in Fatkullin and Slastikov's eqn (27), where τ = 1/B, is −2.2412853 for the axial ODF and 0.5231352 for the planar ODF. The exact solutions are then given in their eqn (5) with Z the normalizing constant such that eqn (9) is satisfied.

Fig. 6 Planar and axial ODFs from the iterative method for the Maier–Saupe kernel with B = 9, compared with exact values from the analytical solutions²⁰ plotted as black dots. The Maier–Saupe kernel gives results qualitatively similar to the Onsager kernel (W(γ) = sin(γ)) while being more amenable to analysis, and this is reflected in the similarity between Fig. 2 and 6.

The shape of ODFs minimizing the free energy will obviously depend on the form of the particle interaction W(γ). In particular, subject to the effect of entropy, ODFs will tend to have orientations around θ values that give rise to angles between orientations where W has minima. Several examples for V = 0 are shown in Fig. 7. Fig. 7a displays a free energy minimizing ODF obtained when W(γ) = sin(γ) − 8P₆(cos γ). Fig. 7b exhibits a case where W is such that the first B at which there might be a branch of anisotropic solutions of the Euler–Lagrange equation off of the isotropic ODF is attained from 2 of the coefficients in the expansion of W in terms of Legendre polynomials (see Appendix 5 in ESI†). This results in a more complex group of local minima. Fig. 7c shows the pair of minima resulting from W(γ) = −(7π/32)P₃(cos γ) which favors parallel particle alignment and disfavors antiparallel particle alignment.

Fig. 7 Distributions resulting from variations in W. (a) ODF for W(γ) = sin(γ) − 8P₆(cos γ), V = 0, and B = 12, beginning the iteration with a planar ODF. (b) Two of the three ODF's obtained for W(γ) = −(5π/32)P₂(cos γ) − (9π/32)P₄(cos γ), V = 0, and B = 11. For this W, there may be a branch of anisotropic solutions at B = 32/π = 10.186, coming from both its P₂ and its P₄ coefficients. At B = 11, i.e., just beyond this bifurcation point, distinct local minima of the free energy are obtained with an initial planar ODF (green), and an initial ODF having peaks near π/4 and 3π/4 (magenta). Not shown is the strongly axial global minimum (with value 54.8 at θ = 0 and π). (c) The mirror distributions resulting from W(γ) = −(7π/32)P₃(cos γ), V = 0 and B = 12. Recall that if f(θ) is a local minimum of F then g(θ) defined to be f(π − θ) will also be a minimum.

The foregoing behavior can be modulated by an applied field V(θ) superimposed on the particle interaction, as illustrated in Fig. 8a and b where W(γ) = sin(γ). In Fig. 8a, the applied field V(θ) = 2 sin(θ) favors particles having orientations toward θ = 0 or θ = π, relative to other orientations. In Fig. 8b, the applied field V(θ) = −cos(θ) favors particles having orientations around θ = 0, and disfavors those around θ = π. In both cases, the iteration obtains the same converged ODF from initial ODFs that are axial, planar or have one peak at either θ = 0 or θ = π.

Fig. 8 ODF modulation by an external field. (a) Results for W(γ) = sin(γ) with B = 12, and an applied field V(θ) = 2 sin(θ) which favors orientations around both θ = 0 and π. The peak value of this ODF is even larger than that for the hard core reference system with B = 14 (cf.Fig. 2). (b) Results for W(γ) = sin(γ) with B = 12, and an applied field V(θ) = −cos(θ) which favors orientations around θ = 0 and disfavors orientations around θ = π. (c) The pair of minima of the free energy for the dipolar kernel with V = 0 and B = 6 (purple and cyan) were obtained from initial values for the iteration having either a peak at θ = 0 or a peak at θ = π. The dots are exact values.²⁰ These ODFs are mirror images of each other about π/2 and are physically equivalent, since each can be obtained from the other by a rotation of spherical coordinates. In the presence of an applied field V(θ) = −cos(θ) which favors orientation toward θ = 0, for starting values with a peak close to θ = π the iteration converges to an ODF (green) that has attenuated orientation toward θ = π. For other starting values, the iteration converges to an ODF (magenta) that has accentuated orientation toward θ = 0. The ODF with the accentuated peak at θ = 0 has the lower free energy in the presence of the applied field.

As expected, the Maier–Saupe kernel continues to give results qualitatively similar to those of the hard core reference system, under the influence of an applied field. When the fields V(θ) = 2 sin(θ) and V(θ) = −cos(θ) were applied with W(γ) equal the Maier–Saupe kernel and B = 9 (data not shown), the iteration obtained the same converged ODF from initial ODFs that were axial, planar or had one peak at either θ = 0 or θ = π.

Fig. 8c shows results for the dipolar potential with B = 6. The pair of minima when V = 0 determined by the iterative method are displayed along with exact values from the analytical solution.²⁰ (B = 6 corresponds to the order parameter r = ±4.73332.) When the field V(θ) = −cos(θ) is applied, one obtains a minimum with an attenuated peak at θ = π when the iteration is started with an ODF with a peak near θ = π. For other starting values, including an ODF with a peak just to the right of π/2, one obtains the global minimum ODF with an enhanced peak at θ = 0.

When there are w_n > 0 for one or more n > 0, the iteration may not reduce the free energy F at some steps, and may not converge. This is illustrated in Fig. 9 for W(γ) = sin(γ) + P₄(cos γ). When , the iteration does converge (albeit non-monotonically in free energy) from both axial and planar initial ODFs. Nevertheless, whenever there is a positive w_n (n > 0), it is best to use an alternate method for finding local minima.

Fig. 9 The free energy vs. iteration number for W(γ) = sin(γ) + P₄(cos γ), V = 0 and B = 10.8 starting from an axial ODF. This W does not satisfy the conditions in eqn (18) since w₄ = 1 − 9π/256 > 0, and the iteration (24) fails to converge. The oscillation in the free energy that is established by iteration 30 was seen to continue through 3000 iterations.

4 Conclusions

The natural iterative method that the calculus of variations suggests for obtaining orientation distribution functions (ODFs) strictly decreases the free energy at each step (unless already at a fixed point of the iteration) under specific and physically appropriate conditions on the particle interaction term. This has significant implications for convergence. In numerical practice, when these conditions are satisfied, convergence is generally quite rapid. It is important to select a variety of initial ODFs in order for the computations to have a good chance to reach all the physically relevant minima of the free energy. The extent to which the behavior of the iteration considered as a dynamical system will represent all stable branches of bifurcations is an open question.

When there are one or more positive coefficients in the expansion of the interparticle potential in terms of Legendre polynomials, the iterative method can fail to converge, and it is advisable to use alternative computational procedures for obtaining minima of the free energy function and/or solutions of the Euler–Lagrange equation, such as described in, for example.²³

Data availability

There is no experimental data for this paper. Relevant code including examples of running the code and corresponding output is available as indicated in ref. 28 in the paper: A. Berger, R-code-for-calculating-orientation-distribution-functions-for-rodlike-particles, 7 Jan 2024. https://github.com/AlanBerger/R-code-for-calculating-orientation-distribution-functions-for-rodlike-particles.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The author thanks Judith Herzfeld for encouragement to extend the results to a wider class of interactions than the original hard core reference system, and for providing physical insight into the conditions for convergence and the computational results. The author also appreciates the excellent suggestions offered by the referees.

References

1. L. Onsager, The effects of shape on the interaction of colloidal particles, Ann. N. Y. Acad. Sci., 1949, 51, 627–659.
2. M. Franco-Melgar, A. J. Haslam and G. Jackson, Advances in generalized van der Waals approaches for the isotropic-nematic fluid phase equilibria of thermotropic liquid crystals-an algebraic equation of state for attractive anisotropic particles with the Onsager trial function, Mol. Phys., 2009, 107, 2329–2358,  DOI:10.1080/00268970903352335.
3. J. Herzfeld, A. E. Berger and J. W. Wingate, A highly convergent algorithm for computing the orientation distribution functions of rodlike particles, Macromolecules, 1984, 17, 1718–1723,  DOI:10.1021/ma00139a014.
4. R. Aliabadi, S. Nasirimoghadam and H. H. Wensink, Evidence of T-type structures of hard square boards in capillary confinement, Phys. Rev. E, 2023, 107, 054117-1,  DOI:10.1103/PhysRevE.107.054117.
5. R. Aliabadi, S. Nasirimoghadam and H. H. Wensink, Capillary-driven biaxial planar and homeotropic nematization of hard cylinders, Phys. Rev. E, 2022, 105, 064704-1,  DOI:10.1103/PhysRevE.105.064704.
6. V. I. Shmotolokha and M. F. Holovko, Effects of porous media on the phase behaviour of polypeptide solutions, Condens. Matter Phys., 2022, 25, 33602-1,  DOI:10.5488/CMP.25.33602.
7. I. Pihlajamma, R. de Bruijn and P. van der Schoot, Continuum percolation in colloidal dispersions of hard nanorods in external axial and planar fields, Soft Matter, 2021, 17, 10458–10468,  10.1039/D1SM01408K.
8. S. M. Ghazi, F. Behzadi and R. Aliabadi, Second-Virial Onsager Theory and Its Limitations in the Prediction of the Ordering Transitions of Confined Hard Rods between Two Parallel Hard Walls, J. Phys. Soc. Jpn., 2020, 89, 114601-1,  DOI:10.7566/JPSJ.89.114601.
9. M. M. C. Tortora, G. Mishra, D. Prešern and J. P. K. Doye, Chiral shape fluctuations and the origin of chirality in cholesteric phases of DNA origamis, Sci. Adv., 2020, 6, eaaw8331-1,  DOI:10.1126/sciadv.aaw8331.
10. K. May, A. Eremin, R. Stannarius, S. D. Peroukidis, S. H. L. Klapp and S. Klein, Colloidal suspensions of rodlike nanocrystals and magnetic spheres under an external magnetic stimulus: experiment and molecular dynamics simulation, Langmuir, 2016, 32, 5085–5093,  DOI:10.1021/acs.langmuir.6b00739.
11. K. Slyusarenko, V. Reshetnyak and Y. Reznikov, Magnetic field control of the ordering of two-component suspension of hard rods, Philos. Trans. R. Soc., A, 2013, 371, 20120250,  DOI:10.1098/rsta.2012.0250.
12. R. F. Kayser Jr. and H. J. Raveché, Bifurcation in Onsager's model of the isotropic-nematic transition, Phys. Rev. A: At., Mol., Opt. Phys., 1978, 17, 2067–2072,  DOI:10.1103/PhysRevA.17.2067.
13. M. A. Cotter, Hard spherocylinders in an anisotropic mean field: a simple model for a nematic liquid crystal, J. Chem. Phys., 1977, 66, 1098–1106,  DOI:10.1063/1.434044.
14. M. A. Cotter, Generalized van der Waals theory of nematic liquid crystals: An alternative formulation, J. Chem. Phys., 1977, 66, 4710–4711,  DOI:10.1063/1.433686.
15. M. A. Cotter, Generalized van der Waals theory of nematic liquid crystals: Requirements for self-consistency, J. Chem. Phys., 1977, 67, 4268–4270,  DOI:10.1063/1.435364.
16. W. M. Gelbart and B. A. Baron, Generalized van der Waals theory of the isotropic-nematic phase transition, J. Chem. Phys., 1977, 66, 207–213,  DOI:10.1063/1.433665.
17. W. M. Gelbart and A. Gelbart, Mol. Phys., 1977, 33, 1387–1398,  DOI:10.1080/00268977700101151.
18. M. A. C. Vollmer, Critical points and bifurcations of the three-dimensional Onsager model for liquid crystals, Arch. Ration. Mech. Anal., 2017, 226, 851–922,  DOI:10.1007/s00205-017-1146-8.
19. H. Liu, H. Zhang and P. Zhang, Axial symmetry and classification of stationary solutions of Doi-Onsager equation on the sphere with Maier–Saupe potential, Commun. Math. Sci., 2005, 3, 201–218.
20. I. Fatkullin and V. Slastikov, Critical points of the Onsager functional on a sphere, Nonlinearity, 2005, 18, 2565–2580,  DOI:10.1088/0951-7715/18/6/008.
21. H. Zhou, H. Wang, Q. Wang and M. G. Forest, Characterization of stable kinetic equilibria of rigid, dipolar rod ensembles for coupled dipole–dipole and Maier–Saupe potentials, Nonlinearity, 2007, 20, 277–297,  DOI:10.1088/0951-7715/20/2/003.
22. J. M. Ball, Axisymmetry of critical points for the Onsager functional, Philos. Trans. R. Soc., A, 2021, 379, 20200110,  DOI:10.1098/rsta.2020.0110.
23. J. Yin, L. Zhang and P. Zhang, Solution landscape of the Onsager model identifies non-axisymmetric critical points, Phys. D, 2022, 430, 133081,  DOI:10.1016/j.physd.2021.133081.
24. H. J. Weber and G. B. Arfken, Essential Mathematical Methods for Physicists, Elsevier Science, Amsterdam, 2004.
25. P. H. Rabinowitz, Some global results for nonlinear eigenvalue problems, J. Funct. Anal., 1971, 7, 487–513,  DOI:10.1016/0022-1236(71)90030-9.
26. W. Rudin, Principles of Mathematical Analysis, McGraw-Hill, New York, 3rd edn, 1976.
27. J. Scheurle, Selective Iteration and Applications, J. Math. Anal. Appl., 1977, 59, 596–616,  DOI:10.1016/0022-247X(77)90084-1.
28. A. Berger, R-code-for-calculating-orientation-distribution-functions-for-rodlike-particles, 2024. https://github.com/AlanBerger/R-code-for-calculating-orientation-distribution-functions-for-rodlike-particles.
29. G. Lasher, Nematic ordering of hard rods derived from a scaled particle treatment, J. Chem. Phys., 1970, 53, 4141–4146,  DOI:10.1063/1.1673914.
30. B. A. Baron and W. M. Gelbart, Molecular shape and volume effects on the orientational ordering of simple liquid crystals, J. Chem. Phys., 1977, 67, 5795–5801,  DOI:10.1063/1.434786.
31. R. Alben, Pretransition Effects in Nematic Liquid Crystals: Model Calculations, Mol. Cryst. Liq. Cryst., 1971, 13, 193–231,  DOI:10.1080/15421407108083541.
32. H. N. W. Lekkerkerker, P. Coulon, R. van der Haegen and R. Deblieck, On the isotropic-liquid crystal phase separation in a solution of rodlike particles of different lengths, J. Chem. Phys., 1984, 80, 3427–3433,  DOI:10.1063/1.447098.
33. M. J. Freiser, Ordered States of a Nematic Liquid, Phys. Rev. Lett., 1970, 24, 1041–1043,  DOI:10.1103/PhysRevLett.24.1041.
34. G. J. Vroege and H. N. W. Lekkerkerker, Phase transitions in lyotropic colloidal and polymer liquid crystals, Rep. Prog. Phys., 1992, 55, 1241–1309,  DOI:10.1088/0034-4885/55/8/003.
35. C. Erignoux and A. Giuliani, Nematic first order phase transition for liquid crystals in the van der Waals–Kac limit, J. Math. Phys., 2020, 61, 103306-1,  DOI:10.1063/5.0007613.
36. R. van Roij, The isotropic and nematic liquid crystal phase of colloidal rods, Eur. J. Phys., 2005, 26, S57–S67,  DOI:10.1088/0143-0807/26/5/S07.

---

Footnote

† Electronic supplementary information (ESI) available: Appendices for conditions for an iterative method for calculating the orientation distribution of rodlike particles to decrease the free energy. See DOI: https://doi.org/10.1039/d4cp02217c

---