Critical points and bifurcations of the three-dimensional Onsager model for liquid crystals

We study the bifurcation diagram of the Onsager free-energy func- tional for liquid crystals with orientation parameter on the sphere in three dimensions. In particular, we concentrate on a very general class of two- body interaction potentials including the Onsager kernel. The problem is reformulated as a non-linear eigenvalue problem for the kernel operator, and a general method to find the corresponding eigenvalues and eigenfunctions is presented. Our main tools for this analysis are spherical harmonics and a special algorithm for computing expansions of products of spherical harmon- ics in terms of spherical harmonics. We find an explicit expression for the set of all bifurcation points. 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 theory of invariant polynomials. On the basis of these result, we show that the first bifurcation occurring in case of the Onsager interaction potential is a transcritical bifurcation and that the corresponding solution is uniaxial. In addition, we prove some global proper- ties of the bifurcation diagram such as the fact that the trivial solution is the unique local minimiser for high temperatures, that it is not a local minimiser if the temperature is low, 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 due to an excluded volume effect 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) dp = 1. (1) Using a second virial approximation, Onsager derived the free-energy functional F (ρ) := S 2 k B τ ρ(p) ln(ρ(p)) + 1 2 U (ρ)(p)ρ(p) dp (2) where U : L 1 (S 2 ) → L ∞ (S 2 ) denotes the interaction operator given by The interaction kernel K(·, ·) : S 2 × S 2 → R describes the excluded volume effect given by the interaction of one test rod with all other polymer rods in the system. In particular Onsager was interested in the kernel which has been named after him. 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 and variations of Onsager's model have been formulated. One of the most popular is the so called Maier-Saupe potential [22] K MS (p, q) = 1 3 − (p · q) 2 .
It is based on an approximation of the mean field approach used by Onsager and assumes that the interactions between the rods are represented by the second moment of the probability density function for the orientation of the rods. A third well-known intermolecular potential is the asymmetric interaction potential which is also called the dipolar potential 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 [12], who extended the Maier-Saupe interaction potential to the case of asymmetric 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 [9] 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 by Fatkullin and Slastikov [11] and independently by Liu, Zhang and Zhang [20]. 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 [31,32] showed that there exist in fact infinitely many bifurcation branches with different symmetries. Revisiting the same question, Chen, Li and Wang [6] 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 [21] 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 [24] who obtained the local bifurcation structure of all equilibrium states in two dimensions and the uniqueness of the trivial solution for high temperatures.
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 derive a method that allows us to compute the eigenvalues of general interaction kernels. On the basis of this result, we obtain a complete set of all bifurcation points and we find the local bifurcation structure of the first of these points corresponding to the first phase transition occurring at the highest temperature. In particular, we prove the existence of a transcritical bifurcation and we show that all critical points of this bifurcation branch are axially symmetric. These results about the local bifurcation structure of the Onsager free-energy functional are summarised in the following theorem.
Theorem 1 (Local bifurcations of the Onsager free-energy functional) Let the Onsager free-energy functional F in (2) be equipped with the Onsager kernel K O (p, q) = 1 − (p · q) 2 . Then a transcritical bifurcation from the trivial solution ρ = 1 4π occurs locally at the temperature and it is uniaxial. Moreover, all other local bifurcations from the trivial solution occur at the temperatures τ s = Γ (s/2 + 1 2 )Γ (s/2 − 1/2) 8k B Γ (s/2 + 1)Γ (s/2 + 2) where s ∈ 2N and the trivial solution ρ = 1 4π is a local minimiser for all τ > π 32kB and it is not a local minimiser for τ < π 32kB . 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 methods apply to a very general class of intermolecular potentials. One critical property of most physical intermolecular potentials is rotational symmetry.
For a brief introduction to spherical harmonics see Appendix A.
In this article we restrict our attention to bifurcations of the Euler-Lagrange equation of the Onsager free-energy functional at ρ = 1 4π equipped with the Onsager kernel. However, our approach is generally applicable as long as sufficient information of the eigenvalues λ l is available. On this basis our methods allow a derivation of an expansion of the bifurcation equation, see Remark 15 for more details. Applying a similar 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 methods presented in this research article. Remark 6. Similar results about the local bifurcation structure of the Onsager free-energy functional with interaction potential satisfying Assumptions 5 (a)-(c) can also be found in the unpublished thesis of Jakob Wachsmuth [30], 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 indeed the eigenvectors of U as well. In contrast, we derive a method that allows us to compute the eigenvalues of U equipped with any kernel admitting an appropriate Taylor expansion, see Theorem 8. -Wachsmuth imposes the assumption that the absolute values of the eigenvalues form a decreasing sequence. Assumption 5 (d) is more general and allows the sequence of eigenvalues to fluctuate if it converges to zero. -Crucially, we prove local uniaxiality of the solutions occurring at the bifurcation point τ ⋆ , 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 5 (a)-(c). Then the uniform distribution ρ 0 (p) = 1 4π is a critical point for all temperature values τ . If τ is such that where M = max p,q∈S 2 K(p, q), then the uniform distribution is in fact the unique solution of the Onsager free-energy functional in three dimensions. Moreover, any non-trivial solution of the Onsager free-energy functional is bounded and all bifurcation branches either meet infinity or they meet another bifurcation branch.
This work is structured as follows. In Section 2 we derive a method that allows us to compute all eigenvalues and eigenfunctions of interaction operators of the form given in (3) which correspond to interaction kernels satisfying Assumption 5. 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 a 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 for high temperatures and the continuity of all bifurcation branches. We complete this work by summarising our results in Section 8.

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 Onsager free-energy functional. In particular, we use the Taylor expansion of the interaction kernel in three dimensions.
The Euler-Lagrange equation of the free-energy functional in (2) is given by where λ := k B τ . Its derivation is well known and can for example be found in [11] or [20]. 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 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. Therefore we assume without loss of generality that Each minimiser of the Onsager free-energy functional must be a solution of the Euler-Lagrange equation. In particular, we prove in Proposition 41 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 temperatures λ = k B τ . 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 is not invertible. Hence we are interested in all non-zero solutions (φ, λ) of and therefore in the nullspace of the operator L λ . Rearranging the equation L λ η(p) = 0, it becomes apparent that this 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 free-energy functional locally around φ 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 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.
Remark 9. Observe that we have not used Assumption 5 (d) in the statement of Theorem 8. Instead, we provide a sufficient condition implying Assumption 5 (d) by assuming that K(p, q) can be written as a convergent Taylor series and that its coefficients a r satisfy condition (9) (if the sum of the absolute values of all eigenvalues in (9) is finite, the sequence of eigenvalues converges to zero).
Proof. By P l (x) we denote the associated Legendre polynomials of order l (for more details see Appendix A). Lemma 35 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 (11) yields Using the addition theorem for Legendre polynomials [36, page 395] we obtain that 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 ∞ r=0 S 2 l=r,r−2,...
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 ∞ r=0 S 2 l=r,r−2,...
Moreover, it follows from [14, Proposition 7.0.1] that any spherical harmonic Y m l is bounded by .
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 (15) is simultaneously proved when we establish our original claim in (14). Therefore a basic condition that suffices to be proved in order to guarantee the validity of the interchange of the integrals in (13) 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 36 we prove this condition in case of the Onsager kernel. For all other cases this condition is assumed to hold (see the statement of the theorem).
Having established (13), 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 4πa s+2r (s + 2r)! 2 r r!(2s + 2r + 1)!! with s ∈ N.
Corollary 10. The eigenfunctions and eigenvalues of the interaction operator U associated to the Onsager kernel k O (p · q) = 1 − (p · q) 2 are given by the spherical harmonics {Y n s ∈ L 2 (S 2 ) : s ∈ N even , −s ≤ n ≤ s} and Proof. Based on the Taylor expansion of the square root function √ 1 − x for all x ∈ (−1, 1), it is easy to see that the Taylor expansion of the Onsager kernel is given by for all p, q ∈ S 2 such that |p · q| < 1. We deduce that if r is even 0 if r is odd and thus using (12) Having expressed the Onsager kernel in terms of spherical harmonics, we may apply (10) in Theorem 8 in order to get where we used the orthogonality of the spherical harmonics and Lemma 36 in which we verify the condition in (16) 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 37 in Appendix B). It reads Based on Corollary 10 and (8), we therefore conclude that a set of all possible bifurcation points is given by Remark 11. We will show in Theorem 34 in Section 7 that in fact any point λ O (s) is a bifurcation point.

The bifurcation equation of the Onsager free-energy functional
Having established a general method for deriving the bifurcation points of the interaction operator, 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 Assumptions 5 (a)-(c). However, in the following analysis we will restrict our attention again 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 (7). 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 [8]. The main idea of the Lyapunov-Schmidt decomposition is the reduction of a possibly infinite-dimensional bifurcation problem to a finite-dimensional one. Let F (w, λ) = 0 (19) 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 (19) decompose as where L(η) := ∂F (0 + ǫη, 0) ∂ǫ ǫ=0 (20) and where R denotes its remainder. We assume that L is a Fredholm operator of index zero, that means that the kernel N (L) has finite dimension d, the range R(L) has finite codimension r and that both 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 (1 − P ) be the projection onto its complement. Then (19) is equivalent to the system of equations where u := P w 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 (21a) because L is linear, and thus it follows that P L(u + v) = 0. Similarly, we can drop the linear part of (21b). The solutions of (22) are equivalent to the solutions of our original problem in (19), see [8], 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.
3.2 Practicalities: An algorithmic procedure to derive the bifurcation equation According to the previous section, the computations for the derivation of the bifurcation equation in (22) are divided into two parts: Step A. An application of the implicit function theorem to (21b) that gives us an explicit expression of v in terms of u and λ.
Step B. Plugging v back into (21a) 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 (21a). 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 (22).
In order to characterise the bifurcation occurring at λ 2 , it is often sufficient to obtain an expansion of (22) 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 (23) is an update of the Euler-Lagrange operator in (7) which has been translated by a factor λ s in order to ensure that the assumption E(0, 0) = 0 holds.
The constant λ s denotes the bifurcation point of interest. According to (8), it is given by 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). In particular, we are looking for the first phase transformation that occurs when the temperature is descending. Hence we consider the largest bifurcation point which is, in 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 non-linear part, we consider its Taylor expansionÊ up to fourth order in both variables φ and λ around the point (0, 0), namelŷ According to (20), 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 T, K : X → Y be linear operators and assume in addition that T is a Fredholm operator and that K is compact. Then T + K is also a Fredholm operator with index(T + K) = index(T ). In Lemma 40 in Appendix C.1 we prove that the interaction operator is a compact operator for all interaction kernels satisfying Assumption 5. In order to prove compactness, it is crucial that the operator satisfies Assumption 5 (d). Schematically, L can be written as 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 and 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 8 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 can in fact be viewed as multiplication operators on the space of functions L 2 (S 2 ), such that Lw(p) = λ 2 + µ l 4π w(p) for all w ∈ L 2 (S 2 ). 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, 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 [8]. 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, see implicit function theorem [3, Theorem 2.3]. In particular, we where the termsv i,j denote terms of i th order in u and of j th 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 (21b) and matching the terms of the right order on both sides, we obtain expressions for each of thev ij . 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 0,0 = 0. 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 (21b) taking the expansion ofÊ in (24) and the definition of R in (25) 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 both sides of (21b), 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 (27). 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 (26), 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 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 (27), 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 5 the expansion of v would only depend on finitely many eigenvalues corresponding to spherical harmonics with low order. That means the resulting bifurcation can then be studied with the methods presented here.

Products of spherical harmonics
The crucial step in matching the coefficients in (21b) 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 subject of this section.
We first recall the definition of spherical harmonics (see Appendix A) where ϕ and θ denote the polar and the azimuth 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 (28) 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 [33] Rearranging it gives and thus, Applying the original recurrence rule in (29) again, but this time to the product xP 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 [33], 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. That means that the algorithm itself remembers cases that it has computed already so that it will stop once it hits a known target.

Step B: Derivation of the bifurcation equation
Having found an explicit expression forv(u, λ) in Section 3.2.2, we are now in the position to plug this expression into (21a) which finally gives us an approximation of the bifurcation equation up to fourth order f (u, λ) := PR(u +v(u, λ), λ).
Remark 16. One can show that the bifurcation equation corresponding to the Onsager free-energy functional is smooth, see Lemma 43.
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 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 equivalent to the solutions of the Euler-Lagrange equation of the Onsager free-energy functional in (2). In order to find the zeroes of (30), 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 i.e. 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 (26) 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.3.
Remark 18. From now on, we will refrain from referring to (32) as bifurcation equation in terms of real spherical harmonics but we will also simply refer to it as bifurcation equation.
4 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 (32), 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. That means, that if we can find a rotation that translates one solution of the bifurcation equation (32) into another, we can consider both 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 [23]. 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 tuple of rotation angles (ϕ R , ψ R , θ R ) such that where T denotes the isomorphism given in (31). In other words, M is the real version corresponding to the complex Wigner matrix D(ϕ R , ψ R , θ R ) given in (33).
Having introduced these concepts, we are now in the position to establish symmetry properties of our bifurcation equation.
Proposition 19. The bifurcation equationf real is equivariant with respect to the Wigner matrices M (ϕ R , ψ R , θ R ), that means that for all M (ϕ R , ψ R , θ R ) and every a ∈ R 5 , Proof. Let R(ϕ R , ψ R , θ R ) 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 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 that f 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 (34) 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 (I 1 (a), I 2 (a)) = (I 1 (b), I 2 (b)), 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 [30] 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 (34) which is based on the Wigner matrix, see (51), we obtain similar expressions for the invariant polynomials with respect to the group action in (34). We conclude this section by proving that the reduced state space S contains a representative of every orbit, see Section 4.2.  (3), respectively. The representation above then arises as the adjoint representation of so(3) [5]. 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 + iV , 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. It follows that 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 x n e n S −1 (ϕ R , ψ R , θ R ) = Proof. Omitted.
Based on Lemma 20, 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 (37) 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 (51), the invariant polynomials that correspond to the group representation based on the Wigner rotation matrix in D.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 isomorphism and Lemma 20.

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 I 1 (a), I 2 (a) | a ∈ R 5 = {I 1 (a), I 2 (a) | a ∈ S} . (36) 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 in fact 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 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 (36) 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 (35), 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ŷ 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, in 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 (39) in the following two sections. First we establish some theory in Section 5.1 before this theory 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 an algebraic equation of a simple form of lower order. This problem is called recognition problem. In our presentation we follow [7, 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 21. Two C ∞ -functions f and g are strongly equivalent, denoted by f ∼ g if there exist S : R n × R m → L(R n , R n ) and U : R n × R m → 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 case of the reduced bifurcation equation in (38), 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 g i u = u ′ for all i ∈ {1, . . . , 6}, 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 3equivariant equivalence.
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 26. Any S 3 -equivariant mapping f : R 2 × R → R 2 can be represented by with h 1 (u) = u 2 1 + u 2 2 and h 2 (u) = u 3 1 − 3u 1 u 2 2 . 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 (32) to an algebraic equation of a simple form and lower order. In particular, we show thatf real is S 3 -equivariantly 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 (40). This in turn is a consequence of Proposition 19. 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 26 that it can therefore be written as Moreover, we can further deduce 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 (32). However, it is sufficient to show that they hold in case of the reduced bifurcation equationf S real in (39) because the first two derivatives agree. In this case the coefficients are given bŷ .
Recall that Remark 17 shows thatf real is a fourth order approximation to f real . Inf 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 27 hold and we can conclude that there exists a one-to-one correspondence between the set of solutions of f S real in (39) 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, that 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 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 where z s : [0, π) → R such that π 0 z s (θ p ) sin θ p dθ p = 0. In particular, the term within the square brackets does not depend on the variable ϕ p , for details see Lemma 45 in Appendix E. This allows us to apply the fundamental lemma of the calculus of variations with respect to the variable θ p and gives us the Euler-Lagrange equation as before. Introducing the thermodynamic potential φ s : S 2 → R as 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 (41) and E in (7) 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 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 case of the full problem and the symmetric problem, we conclude that these are the same. Therefore we proved the following statement.
Theorem 28. All solutions of the Onsager free-energy functional are uniaxial in a neighbourhood of the trivial solution φ = 0 and locally around λ = λ 2 .

Qualitative behaviour of the global bifurcation diagram
Finally, we complete our bifurcation analysis of the Onsager model by proving some properties of the global bifurcation diagram. 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.
Using this result, we show that the free-energy functional is in fact strictly convex for high temperatures and that its trivial solution, ρ 0 = 1 4π is the unique global solution. In Section 7.2 we prove that the trivial solution is a local minimiser for high temperatures and that it is not a local minimiser for low temperatures. Finally, 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 interaction kernels that are continuous and symmetric, see Assumption 5 (a)-(b).
Lemma 29. 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) dp > 1. By the choice of γ, we make sure that ρ ⋆ is a probability density that still integrates to one. Observe that C ⋆ = exp(16M/λ) ≥ 1 > 1 2π and hence A ∩ B = ∅. In order to prove that 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 where ξ 1 (p) ∈ [C ⋆ , g(p) + C ⋆ ] and ξ 2 (p) ∈ [ρ(p), ρ(p) + γ] with ρ ∈ C. Since f ′ (x) = ln(x) + 1, which is monotonically increasing, it follows that ≥λ(ln C ⋆ + 1) A g(p) dp − γλ C (ln(ρ(p) + γ) + 1) dp ≥λ(ln C ⋆ + 1) A g(p) dp − γλ C ln 1 2π + γ + 1 dp 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 K(p, q)ρ(p) dq dp + A×B K(p, q)ρ(q)g(p) dp dq + A×C K(p, q)ρ(q)g(p) dq dp − γC ⋆ A×C K(p, q)ρ(q) dq dp − γ B×C K(p, q)ρ(q) dq dp 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 (42). Since γ < 1 2π we obtain a lower bound on the interaction term.
Thus, we are now in the position to derive an overall bound on F . Putting both lower bounds together and using 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 Corollary 30. Let two sets A and B be given by 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 29.

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.
Proof. The second variation of the Onsager free-energy functional in (2) in L ∞ is given by with z ∈ L ∞ (S 2 ) such that 2 S z(p) dp = 0. The Stone-Weierstraß theorem guarantees that we can find a m,l such that for some ǫ > 0. Using the fact that the interaction operator U turns into a multiplication operator when it is applied to a spherical harmonic, we can rewrite the second variation as where we skip the term l = 0 so that ρ + z still integrates to 1. 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πλ > −µ 2 . 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 ρ = 1 4π is not a local minimiser for all λ < λ 2 .

Existence of a unique solution for high temperatures
We prove that the Onsager free-energy functional is strictly convex in ρ and therefore that the isotropic state, represented by the uniform distribution ρ 0 = 1 4π , is the unique solution for high temperatures. Theorem 32. Let k(·) be a continuous interaction potential and let λ be given such that the following inequality holds 8πM exp(16M/λ) ≤ λ.
Then the Onsager free-energy functional F is strictly convex on the set of all probability densities ρ bounded by C ⋆ , see Lemma 29.
Proof. Recall that the second variation of F is given by k(p · q)z(p)z(q) dp dq ≥ 0, see (43). We know from Lemma 29 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 λ S 2 z 2 (p) ρ(p) dp + S 2 ×S 2 k(p · q)z(p)z(q) dp dq 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/λ) ≤ λ.
Theorem 32 proves that the Onsager free-energy functional is a strongly convex functional for all probability densities ρ ∈ S whenever the temperature parameter satisfies the inequality stated in Theorem 32. Numerically, this value can be approximated by λ ⋆ ≈ 38.205 when M = 1 (for example in case of the Onsager interaction potential).
Corollary 33. Let k(·) be a continuous interaction potential and let λ ⋆ be chosen such that 8πM exp(16M/λ) ≤ λ. 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 30 and Theorem 32, the necessary optimality condition is also sufficient [10, Theorem 5.2.4].

Continuity of all bifurcation branches
We conclude this section by showing that every bifurcation branch of the Onsager functional 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 − U where U is compact, see Lemma 40. Therefore we can rewrite the bifurcation problem as and it is thus of the form of bifurcation problems considered by Rabinowitz [27].
Theorem 34 [27,Theorem 1.3]. Let G : E × R → E be a compact and continuous non-linear operator and consider the problem 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 case of the Onsager functional, all eigenvalues have multiplicity 2n + 1, 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 34.

Conclusion
We found a general method to derive the eigenvalues of integral operators of the form U (ρ) = depending on kernels satisfying Assumption 5 and we derived an explicit expression for the eigenvalues of the Onsager free-energy 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 the temperature τ ⋆ = λ 2 /k B where λ 2 = π 32 denotes the largest of these eigenvalues. Moreover, we showed that the corresponding solution is uniaxial. Globally, we showed that the trivial solution ρ = 1 4π is the unique solution for high temperatures. Finally, we verified that all bifurcation branches either meet infinity or they meet another bifurcation branch.
In order to obtain a complete bifurcation diagram, it is left to show that no isolas occur, that 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. This will be demonstrated in forthcoming work. A A brief introduction to spherical harmonics and associated Legendre polynomials The spherical harmonics are defined as where ϕ ∈ [0, 2π) denotes 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 P m l is the set of associated Legendre polynomials given by One of the most important facts of this set of 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. x r = l=r,r−2,...
Notice that we excluded the boundary cases from the overall sum. Putting everything together yields which concludes the induction step.
Lemma 36 (Supplement to Corollary 10). The change of integral and infinite sum in (13) is valid in case of the Onsager kernel. In other words, we can verify condition (16), stating that Moreover, the eigenvalues of the interaction operator equipped with the Onsager kernel decrease faster than l −3 .
Comparing this with the harmonic p-series, we observe that for some constants C1 and C2 which therefore proves the first part of the claim. The second part follows easily from (46).
C Technical steps in the derivation of the bifurcation equation

C.1 Regularity of critical points and differentiability of E
First we show that the interaction operator is a compact operator on the space H s (S 2 ) for all s ∈ R.
Lemma 40. The interaction operator U is a compact operator H s (S 2 ) → H s (S 2 ) for all s ∈ R.
Proof. By (17) we see that U acts as a multiplication operator applied to the orthogonal basis of spherical harmonics. Due to Assumption 5 (d), we know that its eigenvalues converge to zero. Thus, U is a compact operator.
The following Proposition shows that solutions to the Euler-Lagrange equation are smooth.
Proof. We rearrange the Euler-Lagrange equation as Notice that φ is bounded, see Lemma 29. Thus, taking the limit for a sequence pn → 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 kO(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 [15] one way to define the C 1 -norm is xi cos t − xj sin t, xi+1, . . . , xj−1, xi sin t + xj cos t, xj+1, . . . , xn 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 KO(p, q) = ∞ l=0 l m=−l c l Y m l (p)Y ⋆m l (q), see (18), we can calculate In particular (λs + λ)Xijφ(p) is continuous and differentiable. By iteration of this argument, the statement is proved.
The Sobolev spaces on the sphere can be defined in terms of several equivalent ways through the coefficient in the expansion in spherical harmonics or in local coordinates ϕ and θ. In terms of f = ∞ In order to write down the equivalent norm in local coordinates, we follow [16] and use spherical coordinates ϕ and θ on the sphere, see Appendix A. Hence the Riemannian metric in terms of local coordinates is given by g(ϕ, θ) = sin 2 θ 0 0 1 and the Christoffel symbols are given by The H s (S 2 ) norm is then given by [16] f H s (S 2 ) = s i=0 |∇ s f | 2 dp(g) 1 2 where ∇ s denotes the covariant derivative. The corresponding norm is given by Subsequently, we are only interested in H 2 S 2 and thus we only need ∇f and ∇ 2 f . These are given by Hence explicitly, the norm can be expressed as sin 2θ∂ θ f sin θdϕdθ 1 2 By the Sobolev embedding theorem [15], Lemma 42. 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 First we compute the first and second directional derivatives of E 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 (50) separately using the triangle inequality. In the following we sketch the calculation for the first two terms. In case of the first term, we have .
We now use the definition of the · H 2 (S 2 ) -norm given in (48) on the last of these factors η 2 e −φ−ǫη in order to show that this term can be bounded by an expression of the form b( φ H 2 (S 2 ) , η H 2 (S 2 ) ) η 2 H 2 (S 2 ) . Because the expression in (48) is quite long, we will also break it up in its individual summands and call them T1,i for i = 1, . . . , 3. We start with the first one A f 2 sin θdϕdθ for f = η 2 (φ, θ) exp (−φ − ǫη) (notice that we will drop the (ϕ, θ)-dependence from now on). In particular, This gives us a bound on the term T1,1. In case of the second term, we have Using (a + b) 2 ≤ 2 a 2 + b 2 , we obtain that 14: If l > m + 2, then define 15: If l = m or l = m + 1, then define

C.4 Properties of the Bifurcation Equation
Lemma 43. The bifurcation equation f (u, λ) associated to the Euler-Lagrange equation for the Onsager free-energy functional is smooth.
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.

D.4 Isomorphism between the representation of the group action of SO(3)
acting on R 5 in (34) and the Cartan representation