{'fraction_non_alphanumeric': 0.04408911629890203, 'fraction_numerical': 0.015819209039548022, 'mean_word_length': 4.295326258749153, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Onsager's missing steps retraced

Onsager's paper on phase transition and phase coexistence in anisotropic colloidal systems is a landmark in the theory of lyotropic liquid crystals. However, an uncompromising scrutiny of Onsager's original derivation reveals that it would be rigorously valid only for ludicrous values of the system's number density (of the order of the reciprocal of the number of particles) Based on Penrose's tree identity and an appropriate variant of the mean-field approach for purely repulsive, hard-core interactions, our theory shows that Onsager's theory is indeed valid for a reasonable range of densities.


Introduction
The most successful theory of the formation of a nematic liquid crystal phase due to steric interactions is Onsager's celebrated work on the effects of shape on the interaction of colloidal particles [1]. The model successfully describes not only the nematic-isotropic transition [2], but also the phase separation into an isotropic and nematic phase [3].
The motivation of this paper is the desire to understand the success of Onsager's classical theory for liquid crystals based on purely repulsive, hard-core interactions. This theory, appropriately credited with being the first incarnation of modern density functional theory, was highly influenced by Mayer's theory for real gases and by the cluster expansion formalism [4]. However, when this putative parental connection is scrutinized, as done in some detail in Sect. 2, one readily realizes that the system's number densities that would rigorously grant compatibility for the two theories are ludicrously small, of the order of 1/N, where N is the total number of particles. This at least is the disheartening conclusion one arrives at by following closely and uncompromisingly the steps of Onsager's original derivation.
Then, where does the success of Onsager's theory come from? Here we attempt to answer this question by retracing the steps that were missed in Onsager's derivation. Our answer will not be trivial, but we hope that it will be reassuring: Onsager's theory is well rooted, its domain of validity is limited, but meaningful.
The major mathematical tool of our inquiry will be Penrose's tree identity [5], combined with a mean-field approach to purely steric interactions within the canonical ensemble. In this perspective, Onsager's theory, though unchanged in content, will emerge in a different light. This should not be a surprise: of Onsager we could say what Brydges [6] said of Mayer: "He was a pioneer and his derivations twist and turn". This is the structure of the paper. In Sect. 2, we reproduce and criticize Onsager's original derivation. In Sect. 3, we set the scene where our developments are placed by introducing the Helmholtz free energy in a canonical ensemble of complex particles with multiple degrees of freedom. Sect. 4 is devoted to Penrose's tree identity and its relevance to our estimate of the free energy. In Sect. 5, we combine Penrose's identity with a variant of the mean-field approximation that we find suitable for repulsive, hardcore interactions, which are at the heart of our theory. We shall also arrive there at an explicit and reasonable bound on the number density that marks the range of validity of Onsager's theory. In Sect. 6, we connect our approach with the density functional theory for liquid crystals, of which Onsager's model is perhaps the most rightful ancestor. Finally, in Sect. 7, we collect the conclusions of our paper and comment on the possible avenues that it opens for the study of dense nematic liquid crystals, well above the range of validity of Onsager's theory.
The paper is closed by an appendix, which is rather technical in nature: it presents the major mathematical tool employed in this paper, namely, the graph-combinatoric, asymptotic estimates of Britikov [7] that afford corroborating the validity of Onsager's conclusion, if not his derivation.

Onsager's missing steps
In this section, we condense our critique of Onsager's seminal paper [1]. We adopt Onsager's own point of view, even his language and notation whenever possible. [We refer to his equation (xx) in [1] as (Oxx)].
Onsager's model expresses the irreducible Mayer cluster integrals [4] β 1 and β 2 in terms of particle shape parameters, and estimates the partition function B(T ) using these. The free energy is then readily obtained from the partition function; the equilibrium state of the system is that which minimizes the free energy. Onsager's model has been widely used, with volume fractions up to and in excess of 0.70 (see, for example, [8]). Onsager himself expressed the hope that the results would describe "fairly well" isotropic solutions with volume fractions up to 0.25. We shall thus focus on the range of densities for which his free energy expression can be expected to hold.
Onsager's derivation proceeds as follows. The number of states (measured in units of generalized volume) are given by (O1) where u stands for the total interaction potential and dτ is the element of volume in configuration space. The resulting Helmholtz free energy is given by For repulsive forces, the hard-particle pair-interaction energy is (O17) where and for simplicity, we ignore orientation and regard the q i s' as only coordinates of position. The Mayer function is defined as in (O18) by We note here that Φ ij ≦ 0. Expanding the integrand of (1) gives (O19), Substitution into the integrand of (1) and rearranging terms produces where we have kept, in accordance with [1], the lowest order terms (according to the number of variables τ i which appear in the integrand). We note that, by its very definition, I ≧ 0. Integration gives Now, omitting the higher order terms (not shown) and noting that w ij only depends on the difference of particle coordinates, since there are N(N − 1)/2 ≃ N 2 /2 pairs and N(N − 1)(N − 2)/3! ≃ N 3 /6 triplets, we have, Using the cluster integral notation (O20), in (9) we can set with the units of volume, and with the units of volume squared. Since Φ ij ≦ 0, it follows that β 1 ≦ 0 and β 2 ≦ 0. In terms of these, we have We now consider β 1 and β 2 in more detail. Since w 1j depends only on the difference of coordinates, we can integrate over τ 1 , and obtain where τ 12 indicates the difference of the coordinates of particles 1 and 2. It is clear that β 1 is just the excluded volume of particles 1 and 2; we therefore write where v 0 is the volume of one particle, and note that, for spherical particles in 3D, a = 8. Similarly, for β 2 , we can integrate out τ 1 and obtain which we write as For spherical particles b = 15. We then write equation (12) in terms of a and b, and obtain Defining the volume fraction φ as we rewrite (17) as where I must be non-negative. Assuming that φ is small, we neglect the last term, and note that non-negativity requires that (If the last term is not neglected, the constraint on φ is even stronger.) This is our key criticism. Since the number of states must be non-negative, Onsager's approach requires that the volume fraction occupied by particles be of the order of 1/N. Given that even modest systems have N ≈ 10 20 particles, clearly the region of validity of the model is so small as to be of little practical value. We shall return to equation (19) and (20) subsequently.
Ignoring the constraint on the volume fraction, substitution of (17) into (1) and (2), gives where also Stirling's approximation has been used. Assuming that 1 2 aNφ + 1 3 bNφ 2 ≪ 1 and expanding the logarithm in Taylor's series, we obtain after having substituted for φ in (21) and made use of both (14) and (16). Equation (22) is precisely (O21). It seems therefore that, although not at all obvious from (22) and not pointed out in his paper, Onsager's expression for the free energy is valid only in the regime where the volume fraction of particles is no more than of order 1/N, as prescribed by the hidden, but stringent requirement of positivity for I at the chosen level of approximation. If this region of validity is observed, the theory is essentially useless for any practical system.
One might be temped to argue that the positivity of I in (19) may be maintained by the terms which have been neglected; indeed, this would have to be the case. However, these terms cannot then be neglected as they have been in (21) and (22).
Alternately, one might be tempted to argue, as in [9], that the theory is valid in the limit as the aspect ratio of particles goes to infinity, in which case terms of order greater than 2 may indeed be neglected. However, even in that case, the constraint in (20), which is on density rather than shape, must still hold.
Historically, Onsager's starting point (22) has naively been regarded as a truncation of the density expansion (13.43) in Mayer and Mayer's book [4]. That too, however, had not been fully justified before Onsager's times. * * The proof given in [4] is based on the estimate of the logarithm of a sum with a large number of terms by the logarithm of the maximum term. The Mayers duly warn the reader that this is acceptable only when all terms in the sum are positive, which unfortunately is not the case here. They refer to the work of Born and Fuchs [10] for a rigorous proof of the very same statement, which does not call on that restrictive (and unverified) assumption. However, a critical reading of Born and Fuchs' paper reveals a number of technical assumptions that must be verified for the result to be proven rigorously: Onsager's beautiful and seminal theory of the isotropic-nematic transition is in impressive agreement with experiment and simulation; it must therefore be either correct for realistic densities, or at least nearly correct. However, justification for the theory in this regime is absent in [1]. We hope to provide one here, independent of either Mayers' expansion and its convergence.

Helmholtz free energy
We consider a classical system of N particles, for which x i , i = 1, . . . , N, denote the positions in space of the individual centres of mass and ω i ∈ Ω represents the ith particle's generalized orientation. In general, we think of Ω as an m-dimensional space, on which the action of the rotation group SO(3) in ordinary three-dimensional space is specified by a group G of mappings with G ∋ g R : Ω → Ω corresponding to the rotation R ∈ SO(3), and enjoying the property g R 1 R 2 = g R 1 g R 2 , for all rotations R 1 and R 2 . We shall denote by dx the volume measure in three-dimensional space, and by dω a G-invariant measure over Ω, for which |Ω| is the total measure. Ω is the orientation space. For rigid particles, which represent the most common form of colloids, we can identify Ω with SO(3); the appropriate measure on Ω is then the Haar measure on SO(3), when this latter is viewed as a topological group, and if we identify ω with the triple (ϑ, ϕ, ψ) of Euler's angles then dω is most conveniently expressed as dω = sin ϑdϑdϕdψ (23) and |Ω| = 8π 2 . In the following, whenever needed, we shall denote by q i = (x i , ω i ) the configuration of the ith particle and by q = (q 1 , . . . , q N ) is the collective configuration of all particles. For a given rotation R ∈ SO(3), we shall denote Rq := (Rq 1 , . . . , Rq N ), where Rq i := (Rx 1 , g R (ω i )). Likewise, for any translation a, we shall write aq := (aq 1 , . . . , aq N ), where aq i := (x i + a, ω i ). To shorten our formulae, we shall also write dq i = dx i dω i and, occasionally, dq = dq 1 . . . dq N . We shall further assume that all particles are confined within a region B of volume V in three-dimensional space, so that Q := B × Ω will be their configuration space. Similar general treatments of systems of complex particles can also be found in [12] and [13]. * In a given collective configuration q ∈ Q N , particles interact through a pair-wise potential U = U(q i , q j ), which has the following properties: † one is the hypothesis about the existence of a single maximum point in their application of the method of steepest descent to estimate the partition function; another assumption is stated in their inequality (3.12), which involves a whole range of functions defined as sums of series (which one fails to see how can be justified in practical terms). The first rigorous proof of the convergence of the Mayer expansion in the canonical ensemble came only a few years ago [11].
* For liquid crystal molecular theories, this formalism was also adopted in [14] and Chap. 1 of [15]. †As we learn from standard textbooks, such as [16] and [17], these in general guarantee a good statistical mechanical behaviour of the system. We therefore assume them, though they would not be strictly needed for the applications to hard-core colloidal particles that we envision here.
Hard-core potentials between extended, deformable particles are encompassed in this class: for them, U(q i , q j ) = +∞ whenever the pairs (x i , ω i ) and (x j , ω j ) are such that particle i and particle j would partially overlap in space. We are especially interested in steric interactions, for which U(q i , q j ) is either +∞ or 0, depending on whether particles in configurations q i and q j do overlap or not, respectively. Clearly, for these interactions, b = 0 in (24).
The purely configurational Helmholtz free energy is defined as where k is Boltzmann's constant, T is the absolute temperature, and Z N is the canonical partition function, Here and U is the total potential defined in (24). In (26), the factor 1/N! accounts for the indistinguishability of particles. As is customary, kinetic terms have been omitted in (27), under the assumption that they contribute quadratically to the Hamiltonian and so can be integrated exactly and then disposed of as coefficients of the integral of states in (27); these coefficients depend only on the temperature and shift the free energy in (25) by an inessential additive constant. The case of steric interactions between rigid particles falls in the category for which this assumption is valid. Following Mayer [18], we introduce the definition and rewrite (27) in the following equivalent form The Mayer function Φ ij is particularly meaningful in the case of steric interactions, since then its values are either 0, when particles i and j do not overlap, or −1, when they do; for these interactions, Φ ij is the anti-characteristic function of the overlapping set (in configuration space) for two particles. Clearly, if any two rigid particles overlap, it follows from (29) that that configuration does not contribute to the integral of states G N . The introduction of Mayer's functions is an ingenious device that allows the extension of the integral in (29) to the whole available space without overcounting the number of admissible states. As already done by Mayer and his collaborators in a well-known series of papers [18,19,20,21,22,23], following an earlier suggestion of Ursell [24], (see also the textbook [4], for a systematic exposition of this method * ) we give (29) yet another equivalent form, phrased in the language of graphs: where G N is the collection of all (undirected) graphs on N labelled vertices and the pair (i, j) denotes the edge connecting vertices i and j in a graph G. Here vertices represent particles and edges their interactions. Again, in the paradigmatic case of steric interactions, an edge links any two overlapping particles, whose shared space is inaccessible in reality and must be subtracted from the total count of configurational states. In this view, the Mayer function Φ ij becomes an edge function defined on all possible graphs of G N . It should be noticed that G N also contains the empty graph, which has no edges, but only vertices: when in (30) G is the empty graph, the product over the non-existing edges is conventionally set equal to 1 to warrant accord between (29) and (30).
We shall see in the following section how an identity valid for general edge functions reduces (30) to a form that readily suggests an approximate evaluation of G N , appropriate if only graphs with a limited connectivity could be selected from G N .

Penrose's tree identity
A general, quite unexpected identity proved by Penrose [5] can be used to write (30) in an equivalent and rather suggestive form. In recent years, Penrose's identity has been variously re-interpreted and extended [33,34]. Here, for completeness, we shall present its original formulation, after recalling some elementary graph terminology.
A graph G is said to be connected if for any two vertices in G there is at least one path of edges linking them. The collection of all connected graphs on ν vertices will be denoted by G is called a tree if for any two vertices in T there is precisely one path of edges linking them. The collection of all trees on ν vertices will be denoted by T ν . A forest is a disjoint union of trees. The collection of all forests on ν vertices will be denotes by F ν . * According to [25], Yvon [26] was another precursor of Mayer in writing cluster expansions. He is also credited with the stipulation that cluster integrals in the more general case of a non-uniform system should be identical in their topological form with the cluster expansions valid in the uniform case. This conclusion was implicit in the iterative procedure on which his method is based. For this extension of the cluster expansion formalism, which will also be recalled later in Sect. 6 of this paper, see [27,28,29], whose general ideas can be retraced in Chapter 5 of the book [30]. A later extension, but along different lines, was proposed in [31] (see also [32] for a broad review). All these formulations are in the grand-canonical ensemble, whereas here we are only interested in the canonical ensemble counterpart, which is far less common. Penrose's identity is based on a reduction algorithm, often also called a partition scheme, which associates a tree T ∈ T ν with any connected graph G ∈ G (1) ν and, when reversed, starting from any tree T ∈ T ν , it identifies the maximal graph T * ∈ G (1) ν that reduces to T under the direct partition scheme.
We shall not reproduce Penrose's proof here (see [35], for more pedagogical details), but for completeness we shall illustrate Penrose's partition scheme in a simple case. This scheme consists in the following steps: (a) Choose one vertex in a connected graph G, say vertex 1; (b) Assign a weight to every vertex i = 1, defined as the number w i of edges in the shortest path linking i to 1 in G (by definition, w 1 = 0); (c) Delete all edges between vertices of equal weight; (d) For every vertex i = 1, delete all edges joining vertex i to a vertex with weight w i − 1, apart from the one that among these has the least index.
Each of these steps leaves the weight w i unchanged. The result of the construction is independent of the choice of vertex 1. By repeatedly applying this scheme, any graph G ∈ G (1) ν is reduced to a tree T ∈ T ν . Figure 1 shows a simple example of this reduction.
To construct the maximal graph T * originated by a tree T , Penrose's reverse scheme prescribes the following steps: (a ′ ) Choose one vertex in a tree T , say vertex 1; (b ′ ) Assign a weight w i to every vertex i, as in the direct scheme above; (c ′ ) Join all pairs of vertices with the same weight; (d ′ ) Join every vertex i = 1 of weight w i to all vertices of weight w i − 1 with index greater than the largest index of the vertices of weight w i − 1 to which is already joined in T .  These steps, precisely as the direct steps (a)-(d) illustrated above, leave the weights w i unchanged and produce a graph independent of the choice of vertex 1. By repeated application of these steps to a tree T , we obtain the maximal graph T * that Penrose's direct partition scheme reduces to T . Figure 2 shows the application of the reverse construction to the tree T in Fig. 1(c); as expected in general, the maximal graph T * in Fig. 2(c) is different from the parent graph of T in Fig. 1(a).
We are now in a position to state Penrose's tree identity. Given an edge function Φ defined on all connected graphs G on ν vertices, In a less formal, but perhaps more effective language, identity (31) says that "the sum over all connected graphs on ν vertices is equal to the sum over just tree graphs on ν vertices with more complicated weights" [6]. In (31), such correcting weights are expressed as products over the edges of the graph difference T * \ T , which for the simple examples in Fig. 2 is shown in Fig. 3; simpler expressions can be given by use of a clever hard-sphere construction in two higher spatial dimensions, as shown in [36]. It is worth remarking that for purely steric interactions and so all correcting weights in (31) are bounded in the same interval. Equation (31) is valid for connected graphs, but it is easily extended to all graphs on N vertices in G N , which are of interest here. Let F N denote the collection of all forests on N vertices. For any forest F ∈ F N , we construct the maximal graph F * ∈ G N reducible to F by applying Penrose's reverse scheme to all connected components of F , that is, to all trees of F , and taking the union of the maximal graphs T * l obtained from each individual tree T l ⊂ F . Thus the integral of states in (30) can be rewritten as where the edge function Φ is now meant to be the Mayer function (28). In the following section, making use of (33), we shall introduce a mean-field approximation to G N , which will be proved to be equivalent to Onsager's approximation to the free energy of colloidal orientable particles, within a range of validity which can be quantitatively estimated in a general fashion.

Mean-field approximation
Though exact, formula (33) still offers formidable difficulties to the explicit computation of the total number of states G N , even in the simplest case, that of a gas of hard spheres. For this reason, we resort to a mean-field approximation, which at this level of detail is embodied by replacing Φ in (33) with its average Φ over Q 2 . The reader should heed that with this approximation we lose all information concerning particle orientations and that the resulting free energy will thus be that of an isotropic phase. To recover this lost information, in Sect. 6, we shall revisit Onsager's original "mixture argument", which we still find compelling.
We find it more convenient to introduce the function f := −Φ, which for steric interactions represents the characteristic function of the overlapping set in Q of two particles. Thus, setting f = − Φ , we replace G N with its mean-field approximation G N , for which all integrals in (33) are immediately evaluated: where e(G) denotes the number of edges in a graph G. Though immensely simpler than (33), even (34) is not explicitly computable in terms of f and N. Before introducing a further and possibly drastic approximation in the mean-field formula for G N , which will be shown to be equivalent to Onsager's unjustified approximation, we comment briefly on the physical meaning of f .
For a pair of particles at positions x and x ′ with configurations ω and ω ′ interacting through a potential U(x, ω; x ′ , ω ′ ) enjoying the properties described in Sect. 3, f is formally defined as where Since for all translations a, assuming that B is a macroscopic domain of volume V and that the range of f (where it does not vanish) is much smaller than the linear size of B, as is certainly the case for purely steric interactions, we may write f in (36) as where ∆B := B − B is the Minkowski difference * between B and itself, and the function M is defined as having assumed that the origin 0 belongs to B. Said differently, (38) follows from assuming that the fluid is uniform, so that any particle can be taken as the origin of coordinates.
For purely steric interactions M takes only values 1 or 0; we say that it is the characteristic function of the Mayer region in ∆B × Ω 2 . Equation (38) clearly reveals f as the excluded volume fraction. For rigid particles, Ω can be identified with the rotation group SO(3) and the invariance under rotations of U entails that where the reduced Mayer function M R depends only on the relative displacement r and the relative rotation R of two particles. Formally, where I is the identity in SO(3). Our further approximation, admittedly severe, but meaningful-as we shall see, neglects the excluded volume fraction f in the non-tree component of (34), where we simply replace 1 − f with 1: (42) * The Minkowski difference A − B between sets of points A and B in a Euclidean space represented as vectors issued from a given origin is the set of all possible differences a − b of vectors a ∈ A and b ∈ B. Thus, for example, for B a ball of radius R, B − B is a ball of radius 2R.
It is clear what we renounce with this approximation for G N and where this can go badly wrong: (34) shows that as f grows the admissible states become scant, but this is not captured by (42). Qualitatively, this latter should become inadequate as the number density increases; our approximation should possess the character of a mean-field for dilute systems. We shall see below that our method is capable of determining quantitatively the bound for ρ 0 within which the approximation in (42) is expected to be valid. One may also wonder whether (42) can be interpreted in terms of neglected correlations. Properly speaking, no bound is imposed on the order of correlations being considered: it is more the topology of correlations than their order to be limited. Said differently, the permitted level of correlation among interacting particles does not allow for closed loops in their connectivity network, but it is otherwise unlimited. Properly speaking, G N in (42) should not be called an approximation for the real number of states G N because we can offer no estimate of the error |G N − G N |; rather, at best G N mimics G N instead of approximating it.
No matter how simple (42) may look, computing explicitly the sum over all forests F on N vertices is not a trivial matter. We shall accomplish this task by rearranging the terms in (42) and estimating this sum in the limit of large N. Such an estimate will be performed with a view on obtaining the mean-field free energy per particle, which according to (25) and (26) is defined as This will be especially instrumental to the purpose of computing the thermodynamic limit, where both N and V diverge, while keeping ρ 0 fixed. We first rewrite (42) by grouping all forests F on N vertices by the number n of edges they feature (alternatively, they could be grouped by the number K = N − n of trees they possess): where C(N, n) is the cardinality of the collection F (n) N ⊂ F N of all forests with n edges on N vertices. In the limit of large N, C(N, n) was variously estimated by Britikov [7]. * Full details on Britikov's estimates are given in Appendix A. It is proved there that where x, which is defined as is a quantity of order O(1) in the thermodynamic limit. However, equation (46), which will soon lend itself to the retracing of Onsager's missing steps, is valid only for 0 < x < x 0 , where (see Appendix A) and W is the Lambert function ( [43], § 4.13.10). For rigid particles with average pair-excluded volume v e , and the requirement for the validity of (46) reads as This upper bound for the density is to be compared with the least upper bound found by Lebowitz and Penrose [44] to ensure convergence of the virial series (for the pressure in terms of density) within Mayer's theory in the grand-canonical ensemble. In the case of steric interactions, Lebowitz and Penrose's sufficient condition for the convergence of the virial series reads as * The bound in (50) was derived to ensure validity of our approximate evaluation of G N in the canonical ensemble, and so it has no connection with the bound in (51) which was derived to ensure the analyticity of the equation of state in the grand-canonical ensemble. Though dissimilar in origin, (50) and (51) are both similar requirements of diluteness.
Combining (46) and (44), by use of (47) and the classical Stirling's approximation, we arrive at the following formula for the free energy per particle in the mean-field approximation, which readily reduces to to express the free energy per unit volume for a uniform system. Computing the Legendre transform of F N /V as a function of ρ 0 delivers the pressure P as (54) * Lebowitz and Penrose's condition has recently been extended by Procacci and Yuhjtman [34], but for steric interactions their extended condition reproduces that of Lebowitz and Penrose. The latter was also shown to be sufficient for the convergence of the virial series in the canonical ensemble by Pulvirenti and Tsagkarogiannis [11] and for the analyticity of the Helmholtz free energy as a function of density in the canonical ensemble by Morais and Procacci [45]. Expressing the upper bound in (51) in terms of the Lambert function W was suggested in [46].
For rigid particles such as hard spheres (54) takes the more familiar form P kT ρ 0 which is the van der Waals equation for a purely repulsive interaction potential. In our theory, all these formulae are subject to the same bound on density (50) that makes (46) valid.
Though (54) and (55) have the appearance of Mayer's virial series truncated at the first correction to the ideal gas pressure, here they have a completely different justification: they are originated from neglecting all non-tree graphs in (34) and summing over all possible forests of trees in (42). No hypothesis on the existence of a virial series was ever made, and so no truncation has occurred in (54) or (55). Instead, the term linear in ρ 0 featuring there results from a laborious re-summation on all trees in the average total number of states. Had we been able to include graphs containing cycles in (34) (and to count properly them; combinatorics is the main obstruction here), we would have further corrected the ideal gas equation of state by adding a function of ρ 0 , which need not be a power of the density.
There is one special instance where (54) and (55) are exact: as shown in [47,48] for a fluid of hard-spheres, this is the limiting case of infinite space dimensionality (see also [49] and [50] for the debate that this hypothetical limit has sparked). The proof of (55) in [47] is phrased in the traditional framework of Mayer's cluster expansion: it was shown that all higher order correction terms in the virial expansions vanish in the high dimensionality limit. In the language of our paper, this result is not surprising, since in infinite spatial dimensions all graphs entering (34) are bound to be forests: loosely stated, in infinite dimensions tentative cycles miss closure most of the times.
In our approach, firmly rooted in three-dimensional space, (46) fails to hold for x > x 0 because Britikov's estimates show that the average total number of states G N , as given by (45), then alternates in sign as N increases, in contrast with its very meaning. This is a failure that we ascribe to the minimal connectivity hypothesis, as we call the approximation that led us to (42), since it neglects all the cycles present in the maximal graphs F * generated from forests by Penrose's reverse algorithm. Such a minimal connectivity approximation is bound to break down as the density increases; the merit of Britikov's estimates is to have made explicit the upper bound x 0 of x above which this approximation cannot be trusted. We shall indicate in the following section how this approximation is related to Onsager's free-energy functional for liquid crystals.

Onsager's approximation in density-functional form
Here, building upon (52), which Onsager [1] had taken with little, or no comment from Mayer and Mayer's book [4] (see Sect. 2 for a discussion of such missing steps as they emerge from a critical reading of Onsager's paper), * we shall arrive at Onsager's * Mayer and Mayer's book [4] had been published just a couple of years before Onsager's first announcement of his theory [51]. celebrated functional for systems of rigid, anisomeric particles non-uniformly distributed in orientation space Ω. This is perhaps the first manifestation of the general approach that is now known as the density functional theory.
Onsager's original derivation is ingenious; rigid particles are regarded to be distributed in their orientation as if they were members of different coexisting species. This is the essence of the famous multi-species argument that can be retraced in textbooks (see also [35] for an exposition that parallels closely Onsager's original paper). Rather than reproducing here the classical derivation, we shall present a synthetic, formal argument, which could be extended to more general settings.
We start from assuming that particles are distributed in configuration space Q with a density ρ, which obeys the normalization constraint With this in mind, we regard (53) as providing the local free energy per unit volume in Q. Thus the total free energy functional F [ρ] associated with the distribution density ρ will follow by integrating the local free-energy density over Q. While the integration of the first two terms in the right side of (53) is rather unambiguous, the integration of the last term requires discernment. This is a term quadratic in the density where the extra volume |Ω|V entails a further integration. Combining these formal prescriptions with (35), we justify the following definition for F , which is the standard form of density functional theory [52]. An alternative and possibly more convincing derivation is given in [53], where (57) is obtained through a partition of Q in interacting cells individually described by (53). In both ways, we replace the unweighted average of (minus) the Mayer function in equation (53) with its ρ(q)-weighted average, thus reintroducing the information on particle orientations that had been lost in the mean-field approximation. The functional in (57) is more general than that employed by Onsager, which was intended for homogeneous systems of rigid particles interacting only via hard-core repulsion. We imagine such particles to be distributed in orientation space Ω, which here we identify with SO(3), through a probability density ̺, which is normalized to unity, Therefore we represent the density ρ as where ρ 0 is the number density in (43). By letting ρ in (57) be expressed as in (59) and making use of (56) and (58), we readily arrive at a functional of ̺ expressing the dimensionless free energy per particle as 1 NkT where V e (ω, ω ′ ) is the pair-excluded volume defined as and M is the characteristic function of the Mayer region defined by (39). The average pair-excluded volume v e introduced in (49) is thus related to V e through v e = 1 |Ω| 2 F O is Onsager's free-energy functional for colloidal systems: it is a functional of the distribution density ̺ whose minimizers, representing equilibria, depend on the number density ρ 0 , which features as a bifurcation parameter. The more general functional F in (57) has so far lived, as it were, two independent lives. One life springs from the extension of Mayer's cluster expansions to nonhomogeneous systems first obtained by Stillinger and Buff [28], building on earlier work of Yvon [26] and Morita and Hiroike [27]. In this approach, (57) is simply viewed as a truncated functional in the same spirit as (54) is seen (but not by us) as a truncated virial expansion.
The other life springs from the density functional theory originally put forward by Hohenberg and Kohn [54] to describe the ground state of an electron gas and further extended to encompass atomic fluids by Evans [55], the fluid-solid transition and interfaces by Ramakrishnan and Yussouff [56] followed by Haymet and Oxtoby [57], and liquid crystals by Stecki and Kloczkowski [58] followed by Sluckin and Shukla [59]. * As lucidly phrased in [59], "[w]e may regard the theory expounded in this paper as a mesoscopic theory of the isotropic-nematic transition. In spirit it lies between a fully microscopic theory-in which the liquid partition function is evaluated, and a macroscopic theory-which depends totally on a number of phenomenological parameters." A central role is played in this theory by a functional of the form (57), where f (x, ω; x ′ , ω ′ ) is replaced by the direct correlation function (first introduced by Orstein and Zernike [62]), which in Onsager's case simply reduces to the pair-excluded volume V e (ω, ω ′ ). In this approach, the Mayer function is actually the lowest-order approximation to the direct correlation function in a density-functional expansion of the free energy for a non-uniform fluid (as shown, for example in Sects. 3.8 and 3.9 of [63]).
In our approach, we aspire to give (57) a third life: we see it as mean-field mimic of the free energy generated by a minimal connectivity microscopic network. * It is still a mesoscopic theory, but with a clear microscopic basis and a well-defined range of validity.

Conclusion
Although we were sceptical of the rigour of Onsager's original derivation of his celebrated model for lyotropic liquid crystals, we were confident in its validity, which had been confirmed (at least qualitatively) by both experiments and simulations. Yet the fundamental questions remained as to why and in what range Onsager's model is valid. This paper offers an answer to both issues.
The validity of Onsager's model does not rest on the neglect of all but the first irreducible cluster integral of Mayer's theory. Rather, it rests on the (admittedly complicated) resummation of all possible forests of tree graphs that contribute to the total number of states. Surprisingly, the net result is the same, but the theoretical significance is quite different. It is still an approximation, but its range of validity, given explicitly by the upper bound in (48), is small but reasonable. In this approximation, correlations are not quantitatively limited, only their topology is.
Onsager's bifurcation analysis of the minimizers of the free-energy functional F O in (60) for rigid rods was performed within a one-parameter family of distribution functions ̺. In his analysis, he introduced a dimensionless concentration c, which turns out to be precisely our x in (49). The critical value x c for the first-order isotropicnematic transition is x c . = 3.681 and the isotropic phase encounters a transcritical bifurcation at x = x t . = 4 (see also [35]). These values were confirmed by the combined analytic-numerical strategy adopted in [66] and by a deeper bifurcation analysis recently performed in [67]. Both x c and x t fall considerably above the upper bound x 0 . = 0.139 in (48). So a leap of faith is still required to trust quantitatively the results of Onsager's analysis; we have shown that it is at least qualitatively justifiable.
In graph theoretic language, cycles are to be added in the minimal connectivity strategy, which here only includes trees, to improve the approximate evaluation of the total number of states. They might be classified, for example, as are Husimi trees [68], of which the (ordinary) trees considered in this paper are special cases. Once higher order Husimi trees are considered, the functional F O is expected to change, and not necessarily by supplementing the quadratic form in ̺ now in (60) with a cubic one, but we have no clue yet as to what form this functional may have.
Pursuing our interest in the behaviour of dense nematic liquid crystals, we have recently taken a completely different avenue [69]: instead of building the average total number of states G N through higher levels of connectivity, we have estimated it by recursion on the number of particles N. That approach and the one presented in this paper are completely different, we hope they will be shown to be complementary.

Acknowledgments
We are grateful to P. Hohenberg and R. Petschek for inspiring discussions on the subject of this paper. Moreover, E.G.V. would like to express his gratitude to D. Tsagkarogiannis for having introduced him to the charming intricacies of cluster expansions, to D. Brydges for his learned and suggestive encouragements to pursue the theoretical approach presented here, to A. Procacci for explaining the subtleties of Penrose's tree identity, and to O. Penrose himself for his patience and thoughtfulness in listening to a preliminary presentation of the contents of this paper while he was visiting Oxford. Finally, E.G.V. acknowledges the kind hospitality of the Oxford Centre for Nonlinear PDE, where this work was completed while he was visiting the Mathematical Institute at the University of Oxford. This manuscript has benefited from the comments and suggestions of two Reviewers who have contributed (especially one of them) to improve considerably the presentation of our work. We are grateful to both.

Appendix A. Britikov's estimates
Our estimate for G N in Sect. 5 relies heavily on a paper by Britikov [7] that proves various asymptotic estimates for the number C(N, n) of forests with n edges on N vertices as N → ∞. Here we transliterate in full Britikov's results, * even though we shall essentially use only two.
Let ν(N, n) be the function defined by It is not difficult to show that ν obeys the bound Britikov [7] proved that for N → ∞ dt if ν is bounded, (A.5) * Britikov's original estimates were phrased in terms of the number K of trees in a forest with n edges on N vertices, which is given by K = N − n.
It should be noted that by (A.2) the estimate (A.6) is only valid for n → ∞. Moreover, by setting n = N − 1, which corresponds to the case of a forest with a single tree on N vertices, in the limit as N → ∞, (A.6) delivers Cayley's exact classical result [70], To approximate the sum in (45) for large values of N, we set where ⌊·⌋ denotes the floor of the enclosed numerical quantity and γ is a parameter in the interval (0, 1]. * That is, we cover with straight lines the admissible domain n ≦ N − 1, where C(N, n) is defined. Therefore, for given γ, ν becomes a function of N: and so only fomulas (A.4), (A.5), and (A.6) apply to our case, depending on whether 0 < γ < 1 2 , γ = 1 2 , or 1 2 < γ ≦ 1, respectively. For γ = 1 2 , (A.5), which is to be computed for ν = 0 and applies only if N 2 = N 2 , gives which was obtained by use of Stirling's approximation in the form see § 5.11.1 of [43]. In (A.11), c 0 is a constant that can be evaluated explicitly, Though we have carefully computed c 0 , as we shall do below for the similar constants c 1 and c 2 , we are aware that any multiplicative constant scaling G N does not affect the free energy per particle F N /N in the thermodynamic limit. We start by considering 0 < γ < 1 2 . Making use of (A.4) in (A.8), we compute (−x) n , (A.14) * Similarly, we denote by ⌈·⌉ the ceiling of a numerical quantity. We recall that, for a real number a, ⌊a⌋ is the greatest integer that is smaller than or equal to a, whereas ⌈a⌉ is the smallest integer that is greater than or equal to a.
where we have used (47). If we note that for large N the factor 1 − 2n N 1 2 in (A.14) varies much more slowly than N n n! , as n increases, we can replace it by its average in the range of variability of n. A simple computation shows that and with the aid of maple (A.14) can be rewritten as where Q is the incomplete Gamma function ratio defined as ( [43], § 8. We are interested in the behaviour of the sum in (A. 16) for large N and fixed x. As shown in § 8.11 of [43], a number of classic asymptotic formulae are know for Q(a, z), but in them, invariably, only one argument diverges while the other is kept fixed. Temme [71] proved instead a theorem that describes the asymptotic behaviour of Q(a, z) when both arguments diverge, while their ratio is kept fixed, which is precisely what is needed in (A. 16 is kept fixed as a → ∞, and η = 2(λ − 1 − ln λ).
(A. 20) In the case that interests us λ = −2x, and letting η = α + iβ we convert (A.20) into where µ(x) = 1 + 2x + ln 2x. (A.22) * Temme, who had obtained the leading term of (A.18) in [72], proved in [71] a recursive formula for the remainder, of which (A.18) is the first occurrence. Here, we follow, in particular, Sect. 5 of [71], where a convention about the choice of the root of (A.20) is given for λ ranging in the complex plane.
There are two distinct representations of the sum in (A.16) via (A.18), depending on the sign of µ. There is a single zero of µ(x) for x ≧ 0; this is precisely x 0 , as defined in (48). Also, by use of the asymptotic expression erfc z ≈ 1 √ π 1 z e −z 2 , z → ∞, |arg z| < 3π A key result we obtain from (A.26) is that for x > x 0 the partial sum in (A.25) alternates in sign as N grows. Since the average number of states G N is a positive quantity, this fact is to be interpreted as a signal that the approximation of setting equal to its maximum the contribution to G N from the non-tree graphs in (34) breaks down for x > x 0 . This alone should suffice to confine our attention to the interval 0 < x < x 0 . Nevertheless, we continue our analysis for completeness.
Completing the sum in (A.8) is more complicated, but fortunately it can be estimated via Temme's formula (A.18). Here we sketch the essential steps of our reduction strategy. First, we obtain from (A.6) that C(N, n)(−1) n f n ≈ * Again, as for c 0 and c 1 , the constant c 2 is estimated here for completeness; its value will be inessential to the following development.
where U is Kummer's confluent hypergeometric function. * Exploiting relationships between U and Whittaker's confluent hypergeometric function ( [43], § 13.14.5) and between this latter and the incomplete Gamma function ( [43], § 13.18.5), one can show that which combined with (A.29) finally yields Reduced to this form, the above sum can be evaluated asymptotically for large N by use of (A.18), which here leads us to where x 0 is the same as in (48). Both forms of (A.32) are oscillating functions. To estimate G N , since x 0 < 1 4x 0 , we need to put together (A.11), (A.26), and (A.32) in the interval 0 < x < x 0 , where they may constitute a positive function. Since by (47) and (A.11) the contribution given to the total number of states in (A.8) by the term with n = N 2 (for N even) can be written, to leading order, as (−2ex) N 2 , we arrive at the following estimate where an inessential multiplicative constant has been dropped, A(x, N) is a bounded function which can easily be read off from the first line of (A.32), and µ is defined as in (A.22). Since µ(x) < 0 for 0 < x < x 0 , it follows from (A.33) that for large N G N ≈ V N |Ω| N e −N x , (A. 34) provided that 0 < x < x 0 . This is precisely formula (46) of Sect. 5. Our attempts to simplify this lengthy proof of (A.34) have so far failed.