Critical Points and Bifurcations of the Three-Dimensional Onsager Model for Liquid Crystals

We study the bifurcation diagram of the Onsager free-energy functional for liquid crystals with orientation parameter on the sphere. In particular, we concentrate on the bifurcations from the isotropic solution for a general class of two-body interaction potentials including the Onsager kernel. Reformulating the problem as a non-linear eigenvalue problem for the kernel operator, we prove that spherical harmonics are the corresponding eigenfunctions and we present a direct relationship between the coefficients of the Taylor expansion of this class of interaction potentials and their eigenvalues. We find explicit expressions for all bifurcation points corresponding to bifurcations from the isotropic state of the Onsager free-energy functional equipped with the Onsager interaction potential. A substantial amount of our analysis is based on the use of spherical harmonics and a special algorithm for computing expansions of products of spherical harmonics in terms of spherical harmonics is presented. Using a Lyapunov–Schmidt reduction, we derive a bifurcation equation depending on five state variables. The dimension of this state space is further reduced to two dimensions by using the rotational symmetry of the problem and the invariant theory of groups. On the basis of these results, we show that the first bifurcation from the isotropic state of the Onsager interaction potential is a transcritical bifurcation and that the corresponding solution is uniaxial. In addition, we prove some global properties of the bifurcation diagram such as the fact that the trivial solution is the unique local minimiser if the bifurcation parameter is high, that it is not a local minimiser if the bifurcation parameter is small, the boundedness of all equilibria of the functional and that the bifurcation branches are either unbounded or that they meet another bifurcation branch.


Introduction
The isotropic-to-nematic phase transition of rod-like molecules is one of the most studied phenomena in the theory of liquid crystals. It is characterised by an onset of orientational order when either the concentration of molecules is increased or the temperature of the system is decreased. The first model describing such a phase transition was introduced by Onsager in 1949 [25]. Let ρ : S 2 → R be a probability density function characterising the orientation of the molecules, that is ρ( p) 0 for all p ∈ S 2 and S 2 ρ( p) d p = 1. (1) Using a second virial approximation, Onsager derived an expression for the Helmholtz free-energy per molecule which is given as a sum of an entropy term and a second term modelling the two-body interactions The parameter λ which appears in the entropy term incorporates the thermal or athermal effects between the molecules and therefore depends either on the temperature or on their concentration. The parameter ε 0 has been introduced in order to keep the Onsager free-energy functional dimensionless. The operator U : is called interaction operator. It depends on the interaction kernel K (·, ·) : S 2 × S 2 → R which models the two-body interactions between any two hard core molecules. Based on the assumptions we impose on the characteristic of the interactions, different types of interaction kernels arise. One of the first interaction kernels that has been derived is the so called Onsager kernel It describes the excluded volume effect given by the repulsive interactions between any two slender polymer rods in the system. However, due to the non-analytic nature of this particular interaction potential, the equilibria of the energy functional in (2) cannot be computed explicitly. Subsequently, many theories have been formulated. One of the most popular is the so called Maier-Saupe potential [21] K MS ( p, q) which is based on a mean field approach to dispersion forces, attractive interactions between non-polar molecules. A third well-known intermolecular potential is the dipolar potential for polar molecules The equilibrium states of the free-energy functional above equipped with the Maier-Saupe kernel in (5) have been studied extensively in the past. Among the first researchers who became interested in this problem was Freiser [13], who extended the Maier-Saupe interaction potential to the case of asymmetric biaxial molecules. He proved both the existence of a first-order phase transition to a nematic state and a second-order phase transition to a biaxial state. Constantin, Kevrekidis and Titi [10] established the equivalence of the equilibrium states of the Onsager functional with the Maier-Saupe potential in three dimensions to the solutions of a transcendental matrix equation, and thus derived the high concentration asymptotics in terms of the eigenvalues of this particular matrix. A complete classification of the equilibrium states in three dimensions has been provided independently by Fatkullin and Slastikov [12] and by Liu et al. [19]. Using the properties of spherical harmonics, both groups derived explicit formulae for all critical points and proved their axial symmetry. Similarly, Fatkullin and Slastikov also obtained results for the dipolar interaction potential in (6).
One of the first approaches to the original problem involving the Onsager kernel in (4) has been undertaken by Kayser and Raveché [18]. By reformulating the problem as an eigenvalue problem, they derive an iterative scheme that allows them to compute all axially symmetric equilibria of the functional. Using the Fourier coefficients of the Onsager potential in two dimensions, Wang and Zhou [32,33] showed that there exist in fact infinitely many bifurcation branches with different symmetries. Revisiting the same question, Chen et al. [7] also proved that all solutions are axially symmetric in 2D. For a class of interaction potentials involving only a finite number of Fourier modes, Lucia and Vukadinovic [20] verified the existence of continuous branches in two dimensions and they characterised the structure of the bifurcation diagram in terms of the size of the spectral gaps of the interaction operator. A very recent approach to the same problem has also been undertaken by Niksirat and Yu [23] who obtained the local bifurcation structure of all equilibrium states in two dimensions and the uniqueness of the trivial solution if the bifurcation parameter is big.
However, the problem of classifying all critical points of the Onsager freeenergy functional with the Onsager interaction kernel in three dimensions has not yet been addressed. In this paper, we reformulate the problem as an eigenvalue problem and we present a combination of methods that allows us to compute the eigenvalues of a broad class of interaction kernels. On the basis of this result, we obtain a complete set of all bifurcation points corresponding to bifurcations from the isotropic state and we find the local bifurcation structure of the first of these points corresponding to the first phase transition occurring at the biggest value of λ. In particular, we prove the existence of a transcritical bifurcation and we show that all critical points of this bifurcation branch are axially symmetric. Our results about the local bifurcation structure of the Onsager free-energy functional are summarised in the following theorem: and its corresponding solutions are locally uniaxial, that is the solution ρ 2 corresponding to λ 2 is axially symmetric with respect to one axis. Moreover, all other local bifurcations from the trivial solution occur at the bifurcation points λ s = Γ (s/2 + 1 2 )Γ (s/2 − 1/2) 8Γ (s/2 + 1)Γ (s/2 + 2) where s ∈ 2N and the trivial solution ρ 0 = 1 4π is a local minimiser for all λ > π 32 and it is not a local minimiser for λ < π 32 .
The proof of this theorem consists of two crucial steps. The derivation of an explicit expression of all eigenvalues of the interaction operator U corresponding to the Onsager kernel and the exploitation of the symmetry properties inherent to the problem. In particular, the symmetry can be expressed in terms of a group action acting on the state space of solutions. Because we are only interested in solutions up to rotational symmetry, the state space of solutions can be reduced to the orbit space of this group action. Our approach involving invariant theory for groups is generally applicable to bifurcation problems for functions defined on the sphere and is of interest on its own (for details see Sections 4 and 5).
We would like to point out that even though Theorem 1 only refers to the Onsager kernel, our procedure applies to a very general class of intermolecular potentials which is specified by the following assumption: Assumption 2. The interaction potential K (·, ·) : S 2 × S 2 → R satisfies (a) K (·, ·) is continuous in both variables; (b) K (·, ·) is rotationally symmetric, that is K ( p, q) = K (Rq, Rp) for all R ∈ SO (3) and p, q ∈ S 2 .
An important consequence of rotational symmetry is the following proposition which is based on the basic representation theorem of simultaneous invariants of vectors due to Cauchy. Proposition 3. [29, page 29] The interaction potential is rotationally symmetric if and only if it can be written as a function of a scalar K ( p, q) = k( p · q) where k : R → R.

Remark 4.
Within this paper, we use the notations K ( p, q) and k( p · q) interchangeably to denote the interaction kernel of interest.

Remark 5.
Notice that by Proposition 3, rotational symmetry implies the symmetry of the interaction kernel.
In this article we restrict our attention to bifurcations of the Euler-Lagrange equation of the Onsager free-energy functional from the isotropic state ρ 0 = 1 4π equipped with the Onsager kernel. However, our approach is generally applicable to any physical two-body interaction potential that satisfies Assumption 2. The particular combination of methods that we are using allows us a derivation of an expansion of the bifurcation equation which can be repeated equivalently for any other kernel; see Remark 15 for more details. Following the same steps for a dimension reduction, one is faced with a similar recognition problem that needs to be solved on a case to case basis. All in all, a characterisation of the local bifurcation can be achieved using the combination of methods presented in this research article.

Remark 6.
Similar results about the local bifurcation structure of the Onsager freeenergy functional with interaction potential satisfying Assumption 2 can also be found in the unpublished thesis of Wachsmuth [31], which was only drawn to the attention of the author after completion of this work. The major differences between his and our work are as follows: -Instead of working with the set of spherical harmonics, Wachsmuth uses an equivalent description of this set of functions, namely homogeneous polynomials restricted to the sphere. This allows him to perform an elegant dimension reduction to a state space of two dimensions. -Wachsmuth does not derive an explicit expression for the set of eigenvalues of the interaction operator U . He shows that the Laplace-Beltrami operator and U commute which allows him to conclude that the spherical harmonics are also the eigenfunctions of U . In contrast, we derive a procedure that allows us to compute the eigenvalues of U equipped with any kernel admitting an appropriate Taylor expansion, see Theorem 9. -Wachsmuth imposes the assumption that the absolute values of the eigenvalues form a decreasing sequence. By proving that the interaction operator is compact in H 2 (S 2 ) (see Lemma 44), we show that the sequence of eigenvalues corresponding to the interaction operator is always decreasing. -Crucially, we prove local uniaxiality of the solutions occurring at the bifurcation point λ 2 , while Wachsmuth only shows that the solution is infinitesimally uniaxial.
Last but not least, we also prove properties of the global bifurcation diagram which are summarised in the following theorem: Theorem 7. (Global bifurcations of the Onsager free-energy functional) Let K (·, ·) : S 2 × S 2 → R be an interaction kernel satisfying Assumption 2. Then the uniform distribution ρ 0 ( p) = 1 4π is a critical point for all values of λ. If λ is such that The results of Theorems 1 and 7 are summarised in the following picture ( Fig. 1).
This work is structured as follows. In Section 2 we derive a procedure that allows us to compute all eigenvalues and eigenfunctions of interaction operators of the Carrying out a Lyapunov-Schmidt reduction, we reduce the infinite-dimensional problem to a five-dimensional bifurcation equation in Section 3. In Section 4 we reduce the dimension of the problem further by exploiting the rotational symmetry of our setting which results in a simplified two-dimensional bifurcation equation. We solve the corresponding recognition problem in Section 5. Following these results, we restrict our attention to the case of uniaxial solutions in Section 6. In particular, we prove that the transcritical bifurcation also occurs in case of the restricted problem and we therefore deduce that the solutions of the full problem must be uniaxial as well. Section 7 is devoted to a brief analysis of some global properties of the bifurcation diagram, such as the boundedness of all solutions, the uniqueness of the trivial solution as local minimiser if the bifurcation paramter is high and the continuity of all bifurcation branches. We complete this work by summarising our results in Section 8.
An earlier version of this manuscript can also be found on the arXiv [30].

A Complete Description of All Bifurcation Points of the Onsager Free-Energy Functional
The novel contribution of this section is the full characterisation of all bifurcation points of the Euler-Lagrange equation of the three-dimensional Onsager free-energy functional. In addition to writing the interaction potential in terms of associated Legendre polynomials, an idea that goes back to [12], we also use the Taylor expansion of the interaction kernel in three dimensions in order to generalise the applicability of this approach to a general class of two-body interaction potentials. The combination of these steps yields a direct relationship between the coefficients of the Taylor series and the eigenvalues of the interaction operator. This relationship will be made explicit in the case of the Onsager interaction potential.
The Euler-Lagrange equation of the free-energy functional in (2) is given by its derivation is well known and can for example be found in [12] or [19]. The constant c is obtained through rearranging the equation and imposing the constraint that the orientation distribution function ρ integrates to one, see the conditions in (1). In particular, c is given by which denotes the partition function. Introducing the thermodynamic potential φ : it follows that ρ( p) = Z −1 exp(−φ( p)) and we can rewrite the Euler-Lagrange equation as The addition or subtraction of a constant to the free-energy functional in (2) does not change its minimisers. Due to the rotational symmetry of the interaction kernel, S 2 k( p · q) dq = const. Therefore we assume without loss of generality that (See Remark 12 for details in the case of the Onsager interaction potential). Each minimiser of the Onsager free-energy functional must be a solution of the Euler-Lagrange equation. In particular, we prove in Proposition 45 that all critical points are elements of the space C ∞ (S 2 ). The ground state of the Euler-Lagrange equation is given by φ 0 ( p) = 0 which corresponds to the uniform probability distribution ρ 0 ( p) = 1 4π and solves the Euler-Lagrange equation for all λ. It represents the isotropic phase and we are interested in analysing whether other solutions exist. Mathematically speaking, we are interested in those values of λ for which new solutions emerge. Locally, these are the values of λ for which the implicit function theorem is not applicable to the operator E : H 2 S 2 → H 2 S 2 given by More details about the space H 2 S 2 and the differentiability of E and the regularity of critical points can be found in Appendix C.1. In particular, E is infinitely many times Fréchet differentiable as a mapping from H 2 S 2 → H 2 S 2 , see Lemma 46. The implicit function theorem is not applicable if L λ : H 2 S 2 × R → H 2 S 2 , given by is not invertible. Hence we are interested in all non-zero solutions (η, λ) of and therefore in the nullspace of the operator L λ .

Remark 8.
Notice that the constraint that ρ integrates to one translates into the fact that the integration of η as in (9) yields where we used (7).
Rearranging the equation L λ η( p) = 0, it becomes apparent that this problem is in fact equivalent to the eigenvalue problem for the interaction operator in (3) In particular, we observe that for any eigenvalue μ of U , a solution to the Euler-Lagrange equation is given by the corresponding eigenvector of U and Even though this eigenvalue problem involves only a linear integral operator, it is not trivial. In the following theorem, we derive an explicit expression for all eigenvalues of U , and thus we find the set of all possible bifurcation points of the Onsager freeenergy functional locally around φ 0 = 0. In particular we make use of the fact that functions in L 2 (S 2 ) can be expanded in terms of spherical harmonics, see Appendix A. In the case of the function space H 2 S 2 these expansions converge uniformly.
Then the eigenfunctions of the corresponding interaction operator are given by the spherical harmonics Y m l : S 2 → C. The corresponding eigenvalues are given by where s ∈ N.
Proof. By P l (x) we denote the associated Legendre polynomials of order l (for more details see Appendix A). Lemma 34 in Appendix B states that where the sum is taken over all r l such that l − r ≡ 0 mod 2. By assumption, a r ( p · q) r denotes the Taylor expansion of an interaction potential for all ( p, q) ∈ {p, q ∈ S 2 || p · q| < 1}. Hence an application of (13) yields Using the addition theorem for Legendre polynomials [38, page 395] we obtain that with c r l := Integrating the interaction kernel against an arbitrary spherical harmonic Y n s (q), swapping the order of the sum and the integral and using the orthogonality of spherical harmonics, it follows that Notice that the interchange of the infinite sum and the integral needs of course to be verified. Viewing the infinite sum as an integral with respect to the counting measure, Fubini's theorem applies if In particular, using the Cauchy-Schwarz inequality and the fact that all spherical harmonics have unit mass with respect to the L 2 -norm, we obtain Hence, Again we can view the two sums as integrals with respect to the counting measure and we may swap the order of summation according to Fubini's theorem if the absolute value of the underlying function is integrable with respect to any order of integration. Therefore the validity of swapping the order of summation in (18) is simultaneously proved when we establish our original claim in (17). Therefore a basic condition that suffices to be proved in order to guarantee the validity of the interchange of the integrals in (16) is given by ∞ l=0 r =l,l+2,l+4,...
In particular, one needs to show on a case to case basis that the coefficients |c r l | decay faster than l 3/2 . In Lemma 35 we prove this condition in the case of the Onsager kernel. For all other cases this condition is assumed to hold (see the statement of the theorem).
Having established (16), we deduce that From the last equation we may conclude that the eigenfunctions of the interaction operator U η( p) = S 2 k( p · q)η(q) dq associated to the class of two-body interaction potentials k( p · q) described in the statement of this theorem are given by the spherical harmonics Y n k and that their corresponding eigenvalues are given by Let us briefly check this result in the case of the Maier-Saupe interaction operator In terms of Legendre polynomials, the Maier-Saupe kernel can also be written as where P 2 denotes the second Legendre polynomial. Using (14), we have that Following the reasoning as in the proof of Theorem 9, we may conclude that there is only one non-zero eigenvalue in the case of the Maier-Saupe interaction kernel which is given by Each eigenvalue μ O (s) for s > 0 even has multiplicity 2s + 1. The eigenvalue zero has infinite multiplicity.
Proof. Based on the Taylor expansion of the square root function √ 1 − x for all x ∈ (−1, 1), the Taylor expansion of the Onsager kernel is given by for all p, q ∈ S 2 such that | p · q| < 1. We deduce that and thus using (15), Having expressed the Onsager kernel in terms of spherical harmonics, we may apply (12) in Theorem 9 in order to get where we used the orthogonality of the spherical harmonics and Lemma 35 in which we verify the condition in (19) that allows an interchange of the sum and the integral.
In fact, using the theory of generalised hypergeometric functions, it can be shown that this expression for the eigenvalues can be written in a closed form (for details see Lemma 36 in Appendix B). It reads Based on Corollary 10 and (10), we therefore conclude that a set of all possible bifurcation points is given by Remark 11. We will show in Theorem 33 in Section 7 that in fact any point λ O (s) is a bifurcation point.

The Bifurcation Equation of the Onsager Free-Energy Functional
Having found all bifurcation points from the isotropic state of the Onsager freeenergy functional equipped with a general class of interaction operators, we direct our attention to the corresponding bifurcation equations. We would like to mention that this procedure is quite general and may be applied to a very large class of interaction kernels that satisfy Assumption 2. However, in the following analysis we will restrict our attention to the Onsager potential.
The presentation of our results is divided into two parts. We begin by giving a brief introduction to the Lyapunov-Schmidt decomposition, which is our main tool for computing the bifurcation equation, see Section 3.1. In Section 3.2 we focus on the practicalities that are involved when we apply the Lyapunov-Schmidt decomposition to the Euler-Lagrange operator given in (8). In particular, Section 3.2 is divided into four parts: a reformulation of the problem in the language of Section 3.1, see Section 3.2.1, an algorithmic procedure that allows us to fulfil the decomposition, see Sections 3.2.2 and 3.2.4, and the presentation of an algorithm that allows the fast computation of products of spherical harmonics in terms of spherical harmonics, see Section 3.2.3.

The Theory: The Lyapunov-Schmidt Decomposition
This brief introduction to the Lyapunov-Schmidt decomposition is based on [9]. The main idea of the Lyapunov-Schmidt decomposition is the reduction of a possibly infinite-dimensional bifurcation problem to a finite-dimensional one. Let be the equation of interest with w ∈ X , λ ∈ R and F : X × R → R where X denotes a Banach space. Without loss of generality we assume that The linear and non-linear parts of (22) decompose as where R denotes its remainder. We assume that L is a Fredholm operator of index zero, which means that the kernel N (L) has finite dimension d, the range R(L) has finite codimension r and that the dimensions d and r agree. The idea of the Lyapunov-Schmidt decomposition is to reduce the dimension of the equation and its solution by projecting it onto the kernel of L which is, due to the assumptions on L, finite-dimensional. Let P : X → N (L) denote the projection onto the nullspace of L and let Q be the projection onto its complement R(L). Then (22) is equivalent to the system of equations where u := Pw denotes the projection of the solution onto N (L) and v := (1− P)w its projection onto the complement. The second of these equations can uniquely be solved for v using the implicit function theorem and the solution v(u, λ) can then be plugged back into the first equation, which yields the bifurcation equation Note that we dropped the first term appearing in (24a) because (1− Q)L(u+v, λ) = Lv − QLv = 0. The solutions of (25) are equivalent to the solutions of our original problem in (22), see [9], because it only depends on u, which is an element of the finite-dimensional kernel of L, and thus reduces the possibly infinite-dimensional problem to a finite-dimensional one.

Practicalities: An Algorithmic Procedure to Derive the Bifurcation Equation
According to the previous section, the computations for the derivation of the bifurcation equation in (25) are divided into two parts: Step A. An application of the implicit function theorem to (24b) that gives us an explicit expression of v in terms of u and λ.
Step B. Plugging v back into (24a) which results in the bifurcation equation.
By using the implicit function theorem in the first of these two steps, it is often not possible to derive an expression in closed form. Instead, we obtain a Taylor expansion of v(u, λ) (which is justified by the implicit function theorem [3, Theorem 2.3]) and its coefficients are obtained by matching those of the same order that occur in (24a). The resulting expression for v(u, λ) can then be plugged into the expansion of R. An application of the projection P finally yields an expansion of the bifurcation equation in (25).
In order to characterise the bifurcation occurring at λ s , it is often sufficient to obtain an expansion of (25) up to a particular order. We know that a certain order is sufficient by looking at the so-called recognition problem that corresponds to the problem at hand. More details on the particular recognition problem that we need to apply in order to solve the bifurcation equation can be found in Section 5, where we will show that it suffices to compute the bifurcation equation up to third order.
In the following, we will first restate the problem in the language of Section 3.1 before directing our attention to the first of the two main steps mentioned above in Section 3.2.2. In order to make all computations tractable, we write all expressions in terms of spherical harmonics. An essential tool for the comparison of coefficients with matching order is an algorithmic procedure that allows us the fast computation of products of spherical harmonics with a computer. This will be the subject of Section 3.2.3. We conclude this section by performing Step B above which yields the bifurcation equation in five dimensions.

Reformulation of the Problem
In order to find the minimisers of the free-energy functional, we consider its Euler-Lagrange equation and thus the corresponding Euler-Lagrange operator Remark 12. Notice that (26) is an update of the Euler-Lagrange operator in (8). In order to ensure that the assumption E(0, 0) = 0 holds (see also (7)), the interaction operator k(·) had to be translated by the constant (in the case of the Onsager interaction potential, this value is given by π 2 ) and by λ s in λ.
The constant λ s denotes the bifurcation point of interest. According to (10), it is given by λ s = − μ s 4π , where μ s denotes the eigenvalue of order s of the interaction operator associated to the Onsager kernel, see Corollary 10 for details. The variable λ denotes the deviation from this bifurcation point and will act as bifurcation parameter. Having this translation of the problem in mind, we are interested in bifurcations around the point (φ 0 , λ) = (0, 0). In particular, we are looking for the first phase transformation that occurs when the bifurcation parameter is descending. Hence we consider the largest bifurcation point which is, in the case of the Onsager kernel and according to the result of Corollary 10, given by In order to derive the decomposition of the above functional into its linear and nonlinear part, we consider its Taylor expansionÊ up to fourth order in both variables φ and λ around the point (0, 0), namely According to (23), its linear and non-linear parts therefore decompose as respectively. The fact that L is a Fredholm operator of index zero is a consequence of the following corollary in [1].

Corollary 13. [1, Corollary 4.47] Let X and Y be Banach spaces and let S, K : X → Y be linear operators. Assume in addition that S is a Fredholm operator and that K is compact. Then S + K is also a Fredholm operator with
Schematically, L can be written as In Lemma 47 in Appendix C.1 we prove that the operator is compact for all continuous interaction kernels. Hence in order to apply Corollary 13, we only have to show that the identity is a Fredholm operator of index zero. Let us make the following observations Hence the identity is a Fredholm operator of index zero and thus so is L as claimed.
In view of the results of Section 2, it is easy to conclude that the kernel of the operator L, denoted by N (L), is the eigenspace of the spherical harmonics corresponding to the degree of the bifurcation point of interest, which is in our case λ 2 . Thus, which is a five-dimensional space. It follows that the projection P onto N (L) is given by let us further define u m and v l,m such that respectively. We will see that it will be an essential step to write all expressions in terms of spherical harmonics as this simplifies the actions of all operators. In particular, we may deduce from Theorem 9 that Accordingly, we conclude that and similarly, Since the set of all spherical harmonics is an orthonormal basis for the space L 2 (S 2 ) and orthogonal for H 2 S 2 , these operators in fact diagonalises over the set of spherical harmonics. This leads to an enormous simplification of our subsequent computations.

Step A: An Application of the Implicit Function Theorem
It is straightforward to show that the implicit function theorem is applicable to The operator L is invertible in R(L) and hence its inverse is bounded. Since D u R(0, 0) = 0 by definition and the Euler-Lagrange operator E is sufficiently smooth (see Lemmas 46 and 48), the implicit function theorem is applicable. However, writing down the exact solution v(u, λ) is not possible in practice. Instead, an algorithmic procedure is required which has been presented in [9]. It is mainly based on writing v as a Taylor expansion in terms of both variables u and λ which is justified by the differentiability of the implicit function v(u, λ), see the implicit function theorem [3, Theorem 2.3]. In particular, we expand v = (1 − P)φ up to fourth order byv where the termsv i, j denote terms of ith order in u and of jth order in λ. Furthermore, we identify u as before with By definition v ∈ N (L) ⊥ and thus, without loss of generality, we assume that v i, j ∈ (N (L)) ⊥ for all i, j. Plugging this expression into (24b) and matching the terms of the right order on the two sides, we obtain expressions for each of thev i j . We will execute this calculation for the first couple of terms.

Remark 14.
Because the spherical harmonics form an orthonormal basis for L 2 (S 2 ), the v i, j can be written as an infinite series expansion. In fact, when considering the first few steps of the procedure, we will see that this expansion is finite.
Since there does not exist any term that is constant in both variables, we conclude thatv As a second step, we considerv 1,0 which only depends on u. Hence we only consider terms of order one in u on both sides of (24b) taking the expansion ofÊ in (27) and the definition of R in (28) into account. All terms on the right hand side depend on λ or on higher order terms of u, so they can be neglected when we match the coefficients. Thus, we obtain This is equivalent to solving Lv 1,0 = 0. Because we are assuming thatv i, j ∈ N (L) ⊥ , we conclude thatv 1,0 = 0. Similarly, one can deduce that v 0,1 = 0.
As the first non-zero term, we consider 1 2v 2,0 . Matching the order of terms on the two sides of (24b), we obtain This yields an explicit expression forv 2,0 if we bear in mind that the operator L and its inverse L −1 can actually be interpreted as simple multiplications, see (30). However, in order to be able to apply this multiplication to the right hand side, we need to write it in terms of spherical harmonics. Recall that u itself is already given as an expansion in terms of spherical harmonics, see (29), Hence we only have to expand the terms In fact, because the second factor, being the sum of integrals of spherical harmonics of degree l = 2, is zero. In the case of the second term, we observe that we have to compute u 2 which is based on computing products of spherical harmonics. In order to fulfil this task, we have developed an algorithm that is presented in Section 3.2.3. Applying it to the case above yields an expansion of u 2 in terms of spherical harmonics. We are now in the position to apply the interaction operator U , the projection 1 − P and the operator L −1 which, according to (30), reduces to a simple multiplication in this case. We obtain This concludes the computation ofv 2,0 . The computation of the other terms follows the same pattern. However, because they are more complex and tedious, we will omit them. In order to give a good overview of the individual steps for each case, we present a summary of the procedure in form of a general algorithm. We would like to conclude this section by remarking that the expansion ofv resulting from this procedure is in fact finite and that this is always the case when considering a finite expansions of an equation. The reason for this is that u is a linear combination of spherical harmonics that span the kernel of L which is itself finite. Considering products of up to four terms in u andv i, j , we can therefore only reach spherical harmonics up to degree l = 8 4 since the highest degree forv i, j with i + j = 4 can only be l = 8.

Remark 15.
The operator L depends on the bifurcation parameter and the eigenvalues of the interaction operator U . For any kernel satisfying Assumption 2 the expansion of v would only depend on finitely many eigenvalues corresponding to spherical harmonics of low order. That means the resulting bifurcation can then be studied using the combination of methods presented here.

Products of Spherical Harmonics
The crucial step in matching the coefficients in (24b) is to expand products of spherical harmonics again in terms of spherical harmonics. It is known that this can be achieved using Clebsch-Gordan coefficients which arise in angular momentum coupling [17]. However, it is very hard to compute them explicitly because it takes too much computing time. Therefore we devised an algorithm for the computation of products of spherical harmonics with short run time on the computer. The exact algorithm that we implemented on the computer is given in Appendix C.2 while its crucial ideas are the subject of this section.
We first recall the definition of spherical harmonics (see Appendix A) where ϕ and θ denote the polar and the azimuthal angle corresponding to the unit vector p ∈ S 2 . The functions P m l are called associated Legendre polynomials and their exact definition is also given in Appendix A. Thus, a product of two spherical harmonics is given by Notice that we will abuse the notation and drop the (φ, θ ) or p dependence in due course. Due to the factor e i(m+q) , we observe that the order of each spherical harmonic occurring in an expansion of the product in (31) has to be of order m + q while its degree may be arbitrary. On the basis of this observation, we would like to investigate products of associated Legendre polynomials. In fact, we use the following recurrence formula from [34]: Rearranging it gives and thus, Applying the original recurrence rule in (32) again, but this time to the product x P p,q , yields Using the definition of spherical harmonics, we deduce that and thus, we observe that we reduced the degree of the spherical harmonics Y m l by one or two. We repeat this procedure until we reach either a term of the form Y s s+1 or Y s s (which has to happen eventually). Hitting Y s s+1 first, we use another recurrence rule from [34], namely so that we obtain in any case Y s s Y q p . Again we apply a recurrence rule This last product can then be resolved by using repeatedly if necessary. Again this recurrence formula arises from using recurrence formulae for associated Legendre polynomials. The computation of the first couple of steps of the algorithmic procedure above illustrates how an expansion of a product of two spherical harmonics in terms of spherical harmonics can be obtained. The general method involves a couple of other cases to start with, such as spherical harmonics of negative order for example, which also have to be taken into account and for which special rules need to be defined even though the general method works exactly like this. A detailed list of all rules that are necessary to cover all cases that might occur is given in Appendix C.2. In order to implement this algorithm, it is important to use a conditioned replacement rule that makes it faster. This means that the algorithm itself remembers cases that it has computed already so that it will stop once it hits a known target. Remark 17. In fact, we used a fourth order Taylor approximation of the Euler-Lagrange operator and a Taylor approximation of v(u, λ) up to a joint fourth order in both variables. Thus, we obtained an approximation of the bifurcation equation of at least fourth order. More precisely, the bifurcation equationf =R(u+v(u, λ), λ) = 0 is approximated based on the following expansions of v and R:

Step B: Derivation of the Bifurcation Equation
The following calculation shows thatf agrees with f up to fourth order: The bifurcation equation can also be written as an expansion in spherical harmonics of degree l = 2, and therefore admits the form The solutions of this equation are equivariantly equivalent to the solutions of the Euler-Lagrange equation of the Onsager free-energy functional in (2). In order to find the zeroes of (33), each of the coefficientsf i for −2 i 2 needs to be zero. Therefore we seek the solutions to a system of five polynomials, each depending on the six variables u −2 , u −1 , u 0 , u 1 , u 2 and λ which reduces the infinite-dimensional state space H 2 (S 2 ) of our problem to the ten-dimensional one that is C 5 .
However, we are only interested in real solutions to this equation, so we can reduce the dimension further by restricting the equation onto the space of real spherical harmonics with real coefficients. In particular, the real spherical harmonics are defined as Thus, if the function u( p) in (29) is real-valued, its coefficients can be written as where T : R 5 → C 5 for some a := (a −2 , . . . , a 2 ) ∈ R 5 . Using this transformation, we define the following real version of the bifurcation equation which corresponds to the set of real spherical harmonics Equivalently, we will use the expression in order to denote its approximation up to fourth order. Due to its complicated form, we will not state the bifurcation equation in terms of real spherical harmonics explicitly in this section. Instead it can be found in Appendix C.2.

Reducing the Dimension of the State Space of the Bifurcation Equation
The most common approach to investigate the local structure of a system of polynomial equations at a given solution point is to use its Groebner basis. A Groebner basis is the generating set of an ideal in a polynomial ring over a field K [x 1 , x 2 , . . . , x n ]. The advantage of this method is that it can be used in order to find the dimensionality of the solution set at a given point. However, due to the lengthy form of the bifurcation equation in (35), this method is not directly applicable. Instead, we make use of the symmetries of our problem in order to reduce its dimensionality. In particular we consider the Onsager free-energy functional with a rotationally symmetric intermolecular two-body potential. This means that if we can find a rotation that maps one solution of the bifurcation equation (35) into another, we can consider the two solutions to be equivalent.
In order to represent an arbitrary rotation matrix R(ϕ R , ψ R , θ R ) with respect to the Euler angles ϕ R , ψ R and θ R , we use the so-called y-convention. In this convention, a rotation R(ϕ R , ψ R , θ R ) is given by the product where the first rotation is by an angle ϕ R around the z-axis, the second rotation is by an angle θ R around the y-axis and the third rotation is by an angle ψ R around the z-axis. In particular, the three matrices are given by In contrast to the rotation of points in the state space, the rotation of a function u( p) written in terms of complex spherical harmonics of order l = 2 corresponds to a linear transformation of its coefficients u = (u −2 , . . . , u 2 ). This transformation is given by the 5 × 5 Wigner matrix D(ϕ R , ψ R , θ R ). In particular, this relationship is formulated as for all p ∈ S 2 [22]. If we are given two vectors a, b ∈ R 5 each representing the five coefficients of a real-valued solution written in terms of real spherical harmonics, then we say that a ∼ b if and only if there exists a triple of rotation where T denotes the isomorphism given in (34). In other words, M is the real version corresponding to the complex Wigner matrix D(ϕ R , ψ R , θ R ) given in (36).
Having introduced these concepts, we are now in the position to establish symmetry properties of our bifurcation equation.

Proposition 18.
The bifurcation equationf real is equivariant with respect to the Wigner matrices M(ϕ R , ψ R , θ R ), which means that for all M(ϕ R , ψ R , θ R ) and every a ∈ R 5 ,f be the rotation matrix that corresponds to the Wigner matrix M(ϕ R , ψ R , θ R ). For ease of notation, we drop the dependence of R on the Euler angles ϕ R , ψ R and θ R . Defining Rφ( p) := φ(Rp), we see that Hence the operator E is equivariant with respect to the Wigner matrices. By definition where R denotes the non-linear part of E and is therefore also equivariant. The projection P maps all terms onto the spherical harmonics of degree l = 2 and therefore also preserves the equivariance as well as the approximation of the equation up to fourth order. The restriction to the subspace of real solutions is also a linear transformation and therefore we can deduce thatf real is equivariant with respect to the Wigner matrices M(ϕ R , ψ R , θ R ).
Knowing that f is in fact equivariant with respect to rotations, we can restrict our attention to one representative of each class of solutions that can be translated into each other by a rotation. In other words we are interested in the orbit space that corresponds to the action of SO(3) on the space R 5 . To illustrate this idea, consider the example of SO(2) acting on R 2 . The orbit space in this case can then be taken to be the non-negative part of the x-axis and is therefore one-dimensional.
The aim of this section is to show that the orbit space of SO(3) acting on R 5 can in fact be reduced to two dimensions. In particular, using invariant theory for groups, we can prove that the space contains at least one representative of every orbit.
In particular, one can show that the generators of the ring of G-invariant polynomials for any compact Lie-group G separate the orbits; thus for any two distinct orbits Γ and Γ , there exists at least one of the generators taking different values on Γ and Γ [2, Appendix C]. In the case of the group action given in (37) there exist two invariant polynomials I 1 and I 2 generating the ring. The fact that these generators separate the orbits can be reformulated as follows. If then there exist rotation angles ϕ R , ψ R and θ R such that M(ϕ R , ψ R , θ R )a = b and hence a ∼ b. This allows us to check that S contains at least one representative for every orbit by verifying that In other words, we show that the state space can in fact be reduced to two dimensions. In [31] Wachsmuth achieves a somehow similar two-dimensional reduction by the interesting device of considering the set of homogeneous polynomials restricted to the sphere, which is an equivalent description of the set of spherical harmonics.
This section is structured as follows. In Section 4.1 we find an isomorphic representation of the group action, called the Cartan representation, and we derive an explicit way to separate its orbits based on invariant polynomials. Using the isomorphism between the Cartan representation and the group action given in (37) which is based on the Wigner matrix, see (64), we obtain similar expressions for the invariant polynomials with respect to the group action in (37). We conclude this section by proving that the reduced state space S contains a representative of every orbit, see Section 4.2.

The Separation of Orbits of the Group Action SO(3) Acting on R 5
There exists a unique irreducible representation of the group action of SO(3) acting on R 5 [14]. Since the complex and real representations based on the Wigner matrix given in Appendix E.1 and E.2 are irreducible, it is isomorphic to any other irreducible representation of SO(3) acting on the five-dimensional real space.
One of these isomorphic representations is called the Cartan representation, see Appendix E.3. The explicit isomorphism between the two group actions is given in Lemma 49 in Appendix E.4. The Cartan representation is given by (3), (S, X + so(3)) → S X S −1 + so (3) where su(3) and so(3) denote the Lie-algebras associated with the Lie groups SU (3) and SO(3), respectively. The representation above then arises as the adjoint representation of so(3) [6]. In this setting the matrix X ∈ su(3)/so(3) is a traceless skew-Hermitian matrix, that is X = −X . In particular, if we write X = U + i V , with U and V both being real 3 × 3 matrices, then we can make the following observation Thus, in order for X to be skew-Hermitian, U needs to be a traceless skewsymmetric matrix, thus U ∈ so (3), and V needs to be traceless and symmetric.
Thus, any element in su(3)/so (3) can be written as a linear combination of these matrices, that is for any X ∈ su(3)/so(3) there exist (x 1 , . . . , x 5 ) such that x 5 e n .
In particular, we say that any two elements X and Y ∈ su(3)/so(3) are equivalent, and thus are elements of the same orbit, if there exists a matrix S( x n e n S −1 (ϕ R , ψ R , θ R ) = Similarly, this relationship can also be represented by a 5 × 5 matrix mapping x onto y. In particular, x ∼ y if and only if there exist ϕ R , ψ R and θ R such that where the particular form of M C is given in Appendix E.4. Having established an equivalence relation that describes the orbits in the case of the Cartan representation, we are now in a position to describe its orbit space. We recall that for any two real symmetric matrices X and Y , there exists a matrix S ∈ SO(3) such that if and only if the characteristic polynomials of X and Y agree, that is χ(X ) = χ(Y ).
We deduce that the orbits described by the Cartan representation are separated by the coefficients of the characteristic polynomials of the corresponding matrix. In particular, the characteristic polynomial χ(X ) of a matrix X ∈ su(3)/so(3) written in terms of the basis E in (40) is given by Thus, the two invariant polynomials for the group action written in terms of the Cartan representation are In particular, both of these polynomials satisfy the relationship for j = 1, 2, respectively, and for all (ϕ R , ψ R , θ R ) ∈ [0, π) 2 × [0, 2π). Using the isomorphism given in (64), the invariant polynomials that correspond to the group representation based on the Wigner rotation matrix in E.2 are given by In order to follow the programme outlined at the beginning of this section, it would now only be left to show that I 1 and I 2 are generators of the ring of invariant polynomials. However, in this case, we already know that both polynomials separate the orbits by using the explicit isomorphism and the fact stated in (41).

Reduction to a Two-Dimensional Orbit Space
Having found the fundamental invariants corresponding to the group action of SO(3) acting on R 5 , we have an explicit description of all orbits as each orbit corresponds to exactly one constant value in the image of the two invariants. In order to reduce the state space to S, we need to verify that In particular, we can express the left hand side and the right hand side as respectively. Moreover, we observe that {I 1 (a) = r } is a five-dimensional sphere which is connected and compact. Rewriting the problem like this, we observe that we are interested in the extrema of I 2 for points on the sphere with radius r . Since the sphere is compact, we know that these extrema are attained and an application of the intermediate value theorem gives us that I 2 attains in fact all values in between its extrema. Taking these arguments into account, the verification of (39) has thus been reduced to checking that the extrema defining the intervals agree for all r > 0. Because both invariant polynomials I 1 and I 2 are homogeneous, it actually suffices to verify this condition for the case r = 1. Moreover, since −I 2 (a) = I 2 (−a) and a ∈ S 2 if and only if −a ∈ S 2 , it is enough to prove the equality of the maxima. Hence, we only need to show that the following relationship holds: if I 1 (a) = 1, which thus proves the claim, and we conclude that the state space can in fact be reduced to the two-dimensional space S.

Solving the Reduced Bifurcation Equation
Having shown that the equivariant state space of the bifurcation equation can in fact be restricted to the two-dimensional space S in (38), we are now in the position to solve this reduced problem. We denote the reduced form of the real bifurcation equation byf S real which only depends on the two state variables a 0 and a 2 . Its approximation up to fourth order is given bŷ a 2 a 0 − 5 π (56π(17 + 2240π) + 61)a 3 2 a 0 241472(1 + 32π) 2 .
In particular, we observe that this equation can also be rewritten up to third order in a i for i ∈ {0, 2} and λ aŝ where c := √ 5π 448 and d := 9 3136(1+32π) − 5 1792 . However, we need to find conditions ensuring that we consider the equation up to a sufficiently high order to guarantee that the characteristics of the bifurcation do not change. Such a problem is called recognition problem and it will be presented in the context of the two-dimensional bifurcation equation in (43) in the following two sections. First we establish the theory in Section 5.1 before it will be applied to our particular problem in Section 5.2.

The Recognition Problem in Two Dimensions
The aim of this section is to derive non-degeneracy conditions which ensure a one-to-one correspondence between the solutions of the two-dimensional reduced bifurcation equation and the algebraic equation of a simple form of lower order. This problem is called recognition problem. In our presentation we follow [8,Chapter 9].
A one-to-one correspondence between the sets of solutions of two functions f (u, λ) = 0 and g(u, λ) = 0 for u ∈ R n and λ ∈ R exists if f and g are strongly equivalent.

Definition 19.
Two C ∞ -functions f and g are strongly equivalent, denoted by f ∼ g if there exist S : R n × R → L(R n , R n ) and U : R n × R → R n such that det S(u, λ) > 0, However, often it is not straightforward to verify strong equivalence for particular cases. Instead, it is easier if the problem inherits certain symmetry properties that can be taken into account. In the case of the reduced bifurcation equation in (43), the problem is symmetric with respect to the group action of S 3 acting on the space R 2 . In particular, S 3 acts on any element u ∈ R 2 by where g i are defined as Moreover, it is easy to show that the group action is generated by the elements Based on this representation of the group action of S 3 acting on R 2 , we direct our attention to the concept of S 3 -invariance, S 3 -equivariance and S 3 -equivariant equivalence.

Definition 20. A mapping
for all g ∈ S 3 and u ∈ R 2 .
The aim of Section 5.2 will be to show that the reduced bifurcation equation is in fact S 3 -equivariant and equivariantly equivalent to its approximation up to fourth order. In order to verify the equivariant equivalence, we use the following concepts and results.

Theorem 23. [28, Theorem 1]. Any smooth S 3 -invariant mapping h can be represented as
Theorem 24. Any S 3 -equivariant mapping f : R 2 × R → R 2 can be represented by Based on these results, we are now in the position to state the following proposition which yields the existence of a normal form under appropriate conditions.

Existence of a Transcritical Bifurcation of the Onsager Free-Energy Functional up to Equivariance
The aim of this section is to find non-degeneracy conditions that allow us to reduce the problem of solving the bifurcation equation in (35) to an algebraic equation of a simple form and lower order. In particular, we show thatf real is S 3equivariantly equivalent to a mapping G : R 2 × R → R 2 by using the symmetry properties of f real .
In order to prove that f real is S 3 -equivariant, it is sufficient to show that it is equivariant with respect to the generators g 2 and g 4 of the group action, see (44). This in turn is a consequence of Proposition 18. Notice that for two particular choices of Euler angles we have the following two Wigner matrices and Both linear maps leave the subspace S = {(0, 0, x, 0, y)|x, y ∈ R} invariant. Merely viewed on the space S, they correspond to the two elements that generate S 3 . Hence we conclude that f real is S 3 -equivariant and we know from Theorem 24 that it can therefore be written as Moreover, we can further deduce from Proposition 25 that it reduces to the normal form Since it is not trivial to find the explicit form of the coefficients a(·, ·) and b(·, ·), it is not straightforward to verify that these conditions hold in case of the bifurcation equation in (35).
However, this is sufficient in the case of the reduced bifurcation equationf S real in (43) because the first two derivatives of f S real andf S real agree. In this case the coefficients are given bŷ .
Recall that Remark 17 shows thatf real is a fourth order approximation to f real . In f S real we have dropped the fourth order terms so that it agrees with f S up to third order and therefore also verifies the non-degeneracy condition.
Thus, we see that all three non-degeneracy conditions in Proposition 25 hold and we can conclude that there exists a one-to-one correspondence between the set of solutions of f S real in (43) and The set of solutions of G can easily be computed. It is given by the trivial solution a 0 a 2 = 0 0 and the three branches However, looking at the rotational symmetries of these three branches, we see that they in fact coincide after an application of a rotation. Thus, the solutions corresponding to all three branches are equivalent which proves the existence of a simple linear branch and thus the existence of a unique transcritical bifurcation up to rotation.

Uniaxial Solutions
Recall that the Euler-Lagrange operator corresponding to the Onsager freeenergy functional is given by Now, we restrict our attention to the uniaxial solutions of this problem, which means that we assume that any solution ρ s : L 1 (S 2 ) is axially symmetric with respect to the z-axis. Thus, writing ρ in terms of spherical coordinates θ ∈ [0, π) and ϕ ∈ [0, 2π) (see Appendix A for details), we assume that it only depends on θ and is independent of ϕ. In this case, where S 2 dp = 2π 0 dϕ π 0 sin θ dθ, the Onsager free-energy functional in (2) restricted to the set of axially symmetric probability densities can be rewritten as The Euler-Lagrange equation corresponding to this restricted problem does not differ from the one of the full problem. The reason is that the Euler-Lagrange equation is derived by computing In particular, the term within the square brackets does not depend on the variable ϕ p , for details see Lemma 50 in Appendix F. This allows us to apply the fundamental lemma of the calculus of variations with respect to the variable θ p and this gives us the Euler-Lagrange equation sin θ q dθ q dϕ q sin θ p dθ p dϕ p as before. Introducing the thermodynamic potential φ s : it follows that ρ s (θ p ) = Z −1 s exp(−φ s (θ p )) and we can rewrite the Euler-Lagrange equation as Thus, we conclude that the Euler-Lagrange operator for uniaxial probability distributions is equivalently given by A well-known concept that establishes a relationship between the symmetric solutions of a minimisation problem and the symmetrised problem is the principle of symmetric criticality [26]. It states that each critical point of the symmetric problem is a symmetric critical point of the general problem. In the above setting, this principle can be verified explicitly. Since the Euler-Lagrange operators E s in (45) and E in (8) have exactly the same form, we can consider the corresponding Lyapunov-Schmidt decompositions for the general and the symmetric problem simultaneously. In particular, we split φ = u + v and φ s = u s + v s , respectively. The first step consists in solving for v in terms of u and v s in terms of u s , respectively. Taking the derivative with respect to v and v s , respectively, we see that we can use the implicit function theorem in order to solve uniquely for v s and v. Because is axially symmetric. Therefore, solving the bifurcation equation of the full problem with a −2 = a −1 = a 1 = a 2 = 0 yields all uniaxial solutions. In particular, this equation is given by This is now a bifurcation problem with a one-dimensional state variable and we say that it undergoes a transcritical bifurcation at (0, 0) if (0, 0) is a non-hyperbolic fixed point, that is and if the non-degeneracy conditions hold (these conditions can be found in [39, (3.1.65)-(3.1.68)]). It is easy to see that all of these conditions hold for f s (a 0 , λ) and we may therefore deduce that a transcritical bifurcation occurs. Since we do have a unique transcritical solution in the case of the full problem and the symmetric problem, we conclude that these are the same. Therefore we proved the following statement:

Qualitative Behaviour of the Global Bifurcation Diagram
Finally, we complement our bifurcation analysis of the Onsager model by proving some properties of the global bifurcation diagram. We mainly focus on general interaction kernels satisfying Assumption 2. In Section 7.1 we show that any minimiser of the free-energy functional in (2) is bounded and that we can therefore restrict our attention to the set of bounded probability densities. Following this result, we prove in Section 7.2 that the trivial solution is a local minimiser for large λ and that it is not a local minimiser if λ is small. In Section 7.3, we show that the free-energy functional is in fact strictly convex if the bifurcation parameter λ is large and that its trivial solution, ρ 0 = 1 4π is the unique global solution. Finally, we focus again on the Onsager interaction potential and we show in Section 7.4 that any bifurcation branch either meets infinity or that it meets another bifurcation branch, thus proving that all bifurcation branches are continuous and do not end suddenly.

An Upper Bound on all Admissible Probability Density Functions
By constructing a test function ρ that is bounded above by a constant C , we prove that F(ρ, λ) > F(ρ , λ) for all admissible probability densities ρ and thus we show that we can restrict our attention to probability densities that are bounded above. This result holds for all continuous interaction kernels, see Assumption 2 (a).

Lemma 27.
Let ρ ∈ L 1 (S 2 ) be a probability density and suppose that there exists a set A ⊂ S 2 of positive measure such that ρ( p) C := exp(16M/λ) for all p ∈ A where M := max p,q∈S 2 K ( p, q). Then there exists a modification of ρ denoted by ρ , such that ρ is bounded from above by C and F(ρ , λ) < F(ρ, λ).

Proof. Let
. For the ease of presentation, we use the following notation It is easy to see that |C| > 0 since otherwise S 2 ρ( p) d p > 1 contradicting the fact that ρ is a probability density. By the choice of γ , we make sure that ρ still integrates to one. Observe that C = exp(16M/λ) 1 > 1 2π and hence that A ∩ B = ∅. In order to ensure that ρ C on C, we claim that γ < 1 2π . Since ρ is a probability density, it follows that A (ρ( p) − C ) d p 1. Hence and we only have to show that |C| = |{p| 1 2π ρ( p)}| 2π . Suppose the contrary, that is |C| < 2π and thus |A ∪ B| = |{ρ( p) > 1 2π }| 2π since the total measure of the unit sphere is given by 4π and thus |A ∪ B ∪ C| = 4π . Moreover, since ρ is positive, contradicting the fact that ρ integrates to one. Therefore γ < 1 2π . Having defined ρ , we are now in a position to prove that In the following, we deal with the entropic and the interaction term separately. We begin by deriving a lower bound on the entropy term. In particular, where f (x) := x ln(x). Moreover, defining g( p) := ρ( p) − C , we may apply the fundamental theorem of calculus (notice that the following claim holds trivially for g( p) = 0, so we may assume that g( p) = 0 for all p ∈ S 2 ). Hence Since f (x) = ln(x) + 1, which is monotonically increasing, it follows that where the last line follows from the fact that ρ ∈ C and hence ρ( p) < 1 2π for all p ∈ C. Observing that g( p) 0 for all p ∈ A and using the monotonicity of the logarithm as well as the facts that γ < 1 2π , ln 1 π < −1 and λ, γ , |C| > 0, we deduce that This gives us a lower bound on the entropic term and we can direct our attention to a lower bound on the interaction term. Using the definition of ρ , it follows that where we used the symmetry of the kernel and again we defined g( p) := ρ( p)−C and h( p, q) := ρ( p)ρ(q) − C C which are both non-negative for p, q ∈ A. Assuming without loss of generality that the interaction kernel is positive, we observe that all terms with a positive sign are positive while only those with a negative sign are negative. Neglecting the positive terms, we may therefore deduce that where we used the fact that ρ( p) C for all p ∈ A in (47). Since γ < 1 2π we obtain a lower bound on the interaction term.
We are now in the position to derive an overall bound on F. Using (46), (48) and the definition of γ := A g( p)dp |C| , we obtain that Finally, making use of the assumption that C = exp(16M/λ) and using the fact that |C| 2π , the statement follows because

Then any solution of the minimisation problem min A F(ρ, λ) is in fact a solution of the reduced minimisation problem min B F(ρ, λ).
Proof. The statement follows directly from Lemma 27.

Local Properties of the Trivial Solution
We concentrate on proving conditions which ensure when the trivial solution is a local minimum and when it is not.

Proposition 29. Let U be an interaction operator that diagonalises over the set of spherical harmonics, that is
Let μ inf = inf l μ l , then for all λ > − μ inf 4π the trivial solution ρ 0 = 1 4π is a local minimum. Moreover, for all 0 < λ < − μ inf 4π , the trivial solution ρ 0 = 1 4π is not a local minimum.
Proof. The second variation of the Onsager free-energy functional in (2) in H 2 S 2 is given by with z ∈ H 2 S 2 such that S 2 z( p) d p = 0. In particular, we may write and taking into account that z integrates to zero, we may w.l.o.g. assume that a 0 = 0. On this basis, we can rewrite the second variation for ρ 0 = 1 4π as Viewing the spherical harmonics as an orthonormal basis in L 2 (S 2 ), it follows that Hence the second variation is positive, that is I (z, λ) > 0 for all 0 = z ∈ L 2 (S 2 ), if and only if 4πλ > −μ inf . Thus, the trivial solution ρ 0 is a local minimiser of the bifurcation equation for all Moreover, we observe that the following holds: If we choose z close enough to zero, then the higher order terms are negligible and we see that ρ 0 = 1 4π is not a local minimiser for all λ < − μ inf 4π .

Existence of a Unique Solution for High Values of λ
We prove that the Onsager free-energy functional is strictly convex in ρ for all continuous interaction kernels and therefore that the isotropic state, represented by the uniform distribution ρ 0 = 1 4π , is the unique solution if λ is big.
Theorem 31. Let k(·) be a continuous interaction potential and let λ be given such that the following inequality: holds for M := max p,q∈S 2 k( p · q). In particular, the critical value λ can be written as where W denotes the Lambert-W function. Moreover, the Onsager free-energy functional F is strictly convex on the set of all probability densities ρ bounded by C = exp(16M/λ), see Lemma 27.
Proof. Recall that the second variation of F is given by see (49). We know from Lemma 27 that it is enough to consider the minimisation problem among the set of probability densities ρ ∈ L 1 (S 2 ) which are bounded by a constant C = exp(16M/λ) where M denotes an upper bound on the kernel k( p ·q). Therefore an application of Jensen's inequality and imposing the condition that 8π M λ/C yield proving that the Onsager functional is strictly convex. Rearranging this condition and plugging in the constant C = exp(16M/λ) yields the above inequality 8π M exp(16M/λ) λ.
In fact Theorem 31 proves that the Onsager free-energy functional is a strongly convex functional for all probability densities ρ ∈ S as its second derivative is bounded away from zero by a factor which depends on the bifurcation parameter λ. Numerically, this value can be approximated by λ ≈ 38.205 when M = 1 (for example in the case of the Onsager interaction potential).

Corollary 32. Let k(·) be a continuous interaction potential and let λ = 16π
W (2/M) . Then the uniform density ρ 0 = 1 4π is the unique solution to the Onsager free-energy functional for all λ λ .
Proof. The uniform density ρ 0 = 1 4π is indeed a solution to the Euler-Lagrange equation for all λ ∈ R. By Corollary 28 and Theorem 31, the necessary optimality condition is also sufficient [11, Theorem 5.2.4].

Global Behaviour
We conclude this section by showing that every bifurcation branch of the Euler-Lagrange equation associated with the Onsager free-energy functional equipped with the Onsager kernel is either unbounded and meets infinity or that it meets another bifurcation branch.
The Euler-Lagrange equation in terms of the thermodynamic potential is given by which is of the form λId − T where T is compact, see Lemma 44. Therefore we can rewrite the bifurcation problem as and it is thus of the form of bifurcation problems considered by Rabinowitz [27].
where u ∈ E and λ ∈ R. Moreover, let where H (u, λ) is O(||u||) for u near the origin uniformly for bounded λ and let L be a compact and linear map on E. If λ is a real non-zero eigenvalue of L of odd multiplicity, then the solution set {(u, λ) | u = G(u, λ)} possesses a maximal connected and closed subset C u such that (0, λ) ∈ C u and C u either 1. meets infinity in E, or 2. meets (0,λ) where λ =λ is a real non-zero eigenvalue of L.
In particular, Rabinowitz shows that any eigenvalue of odd multiplicity is a bifurcation point and that any branch bifurcating from the trivial solution must either meet infinity or meet another bifurcation branch. In the case of the Onsager free-energy functional with Onsager potential, all eigenvalues have multiplicity s = 2n + 1 for n ∈ N, thus we know that (0, λ s ) is not an isolated solution but that it is instead a member of a non-trivial closed connected set described by Theorem 33.

Conclusion
We presented a procedure to derive the eigenvalues of integral operators of the form depending on kernels that are both continuous and rotationally symmetric. Moreover, we derived an explicit expression for the eigenvalues of the Onsager freeenergy functional involving the Onsager kernel. Based on this result, we know from [27] that all eigenvalues are bifurcation points and we proved that a transcritical bifurcation occurs at λ 2 = π 32 denoting the largest of these eigenvalues. Moreover, we proved that the corresponding solution is uniaxial. Globally, we showed that the trivial solution ρ 0 = 1 4π is the unique solution for high values of λ. Finally, we verified that all bifurcation branches from the isotropic state of the Onsager free-energy functional equipped with the Onsager kernel either tend to infinity or meet another bifurcation branch.
In order to obtain a complete bifurcation diagram, it is left to show that no isolas occur, which means that we do not have any bifurcations that are isolated from all other solutions to our problem at hand. Our approach is not at all limited to the Onsager kernel and not even to quadratic forms given by integral kernels.
Acknowledgements. The author would like to thank John Ball, Yi-Chao Chen, Epifanio Virga and Markus Upmeier for insightful and stimulating discussions relating to the work in this paper. The research leading to these results has received funding from the Eu- where ϕ ∈ [0, 2π) is the polar angle and θ ∈ [0, π) the azimuthal angle. The subscript l denotes its degree while −l m l denotes its order. N lm is a normalisation constant given by and we write P m l for the associated Legendre polynomial given by The Legendre polynomials P n (x) are defined as One of the most important facts about spherical harmonics is that they build an orthonormal basis of L 2 (S 2 ). In general, spherical harmonics are known as eigenfunctions of the Laplacian. However, it seems that they can be understood in a much broader sense as eigenfunctions of rotationally symmetric operators.
where P l denotes the Legendre polynomial [36].
Proof. The claim follows from an induction based on the reccurrence relation for Legendre polynomials with the convention that P −1 = 0 [34]. The base case r = 0 is trivial because P 0 = 1, so let us move on to the actual induction step. Multiplying both sides of (50) by x yields x r +1 = l=r,r −2,...
Changing the index in the two parts of the sum corresponding to P l+1 and P l−1 , we obtain Notice that we excluded the boundary cases from the overall sum. Putting everything together yields which concludes the induction step. (16) is valid in the case of the Onsager kernel. In other words, we can verify condition (19), stating that

Lemma 35. (Supplement to Corollary 10) The change of integral and infinite sum in
Moreover, the eigenvalues of the interaction operator equipped with the Onsager kernel decrease faster than l −3 .
Proof. With regard to (21), the Onsager kernel written in terms of spherical harmonics is given by The coefficients c r l are all negative except for c 0 0 = 4π . Taking this into account and applying Lemma 36, we know that In particular, we know from Section 5.6 in [24] that for b − a ≥ 1, a ≥ 0 and z = x + iy with x > 0, Moreover, for 0 < s < 1 [24]. An application of these two results to the explicit expression of the eigenvalues gives us the following result Claim: Using (52) and (53) with x = l 2 , s = 1 2 , z = l 2 − 1 2 , a = 0 and b = 5 2 yields Inserting this result into (51) yields Comparing this with the harmonic p-series, we observe that for some constants C 1 and C 2 which therefore proves the first part of the claim. The second part follows easily from (54).
In the case of odd integers n = 2k − 1 for k 1, the double factorial n!! can be replaced by the expression n!! = (2k)! 2 k k! [37, page 823]. Hence we can rewrite the terms such that Notice that we also expressed the first factor in the denominator as fraction of two factorials. This sum can now further be manipulated by rewriting it in terms of gamma functions [37, page 1137], Following these results, in C.2 1. we present an algorithm which can be used to express products of spherical harmonics in terms of spherical harmonics (Algorithm 2), 2. and we state the real bifurcation equation up to third order (63).

C.1 Properties of the Euler-Lagrange Operator
First of all, let us clarify the exact definitions and notation that we are using.
Using very similar arguments, one can prove the following more general statement.

Corollary 40. The H k -norm
where a lm denote the coefficients of the expansion of f in spherical harmonics, that is f = ∞ l=0 l m=−l a lm Y m l . Notice that D α acts in tangential direction to S 2 . Lemma 41. Let U be the operator Proof. Suppose that η ∈ H k and let Let U be the integral operator with U Y lm = μ l Y lm where μ l are the corresponding eigenvalues which are given by l even 0 otherwise Rewriting the above yields Due to the fact that η ∈ H 2 (S 2 ), we know that the first part is bounded (see (58)). Hence we are seeking a bound on the term A l which depends on the rate of convergence of the eigenvalues μ l . Using the result of the claim in (54), we know that Thus, ifk k + 3, then implying that Lemma 42. Let T be the operator and assume that k(·) is an interaction kernel satisfying Assumption 2. Then T : Proof. Let φ ∈ H 2 (S 2 ) ⊂ L 2 (S 2 ). Using Lemma 41, we know that Hence where the last step follows by the Sobolev Embedding theorem Lemma 43. The Euler Lagrange operator Proof. Again using the Sobolev Embedding theorem, we know that ||φ|| ∞ C||φ|| H 2 (S 2 ) . Hence This fact together with the reasoning that we used in the proof of Lemma 42, we know that Lemma 44.
Proof. Again, similarly to the proofs of Lemmas 42 and 43, we know that The statement follows because H 3 (S 2 ) is compactly embedded in H 2 (S 2 ).
Proof. We rearrange the Euler-Lagrange equation as Notice that φ is bounded, see Lemma 27. Thus, taking the limit for a sequence p n → p on the right hand side, we see that it is continuous due to the dominated convergence theorem and so is the left hand side. Because k O ( p ·q) is differentiable almost everywhere, except for | p · q| = 1, with an L 1 bounded derivative, we can swap integration and differentiation and conclude that φ is differentiable.
Having proved that φ is continuous and differentiable, we want to show that also the derivative of φ is continuous and differentiable. This will then allow us to proceed by an induction argument.
By [16] one way to define the C 1 -norm is x i cos t − x j sin t, x i+1 , . . . , x j−1 , x i sin t + x j cos t, x j+1 , . . . , x n are the intrinsic derivatives on S n−1 . Here we consider the special case n = 3. An elementary calculation shows that Thus, using the expansion of K O ( p, q) = ∞ l=0 l m=−l c l Y m l ( p)Y m l (q), see (21), we can calculate In particular (λ s + λ)X i j φ( p) is continuous and differentiable. By iteration of this argument, the statement is proved.

Lemma 46. E is infinitely many times Fréchet differentiable as a mapping from
Proof. In order to check Fréchet differentiability of E, we need to verify that We take the first and second directional derivatives of These are given by Thus, by the integral form of the remainder in Taylor's theorem in ε, we have that The · H 2 (S 2 ) -norm of the left hand side can be bounded by dε and our aim is now to establish an overall bound of the form where b : R 2 → R is locally bounded. This would yield that E is Fréchet differentiable. In order to show this, we treat each of the four terms occurring in (61) separately using the triangle inequality and the same reasoning as in the proof of Lemmas 41 and 42. In particular, using the Sobolev embedding theorem yields Bounds for all other terms and higher order derivatives can be obtained similarly and will be omitted. 13: If l = |m|, then define 14: If l > m + 2, then define 15: If l = m or l = m + 1, then define
Since L is a linear operator, the smoothness of E and R is the same. Moreover, the implicit function theorem tells us, that the smoothness of R and v is the same.
Since the smoothness of the bifurcation equation is given by the minimal number of derivatives of R and v (using the chain rule) which are then equal, we conclude that the smoothness of f is the same as the smoothness of E (the Euler-Lagrange equation). The Euler-Lagrange operator is infinitely many times differentiable, see Proposition 46.

E.1 The Wigner Rotation Matrix for Degree l = 2 Spherical Harmonics
The Wigner rotation matrix for complex spherical harmonics of degree l = 2 is given by

E.4 Isomorphism Between the Representation of the Group Action of SO(3) Acting on R 5 in (37) and the Cartan Representation
Lemma 49. (Isomorphism between the group representations). An isomorphism between the two representations of the group action is given by In particular M C = Φ MΦ −1 .
Thus, the expression is invariant with respect to rotations around the z-axis and therefore independent of the variable ϕ p .