Complete integrability and equilibrium thermodynamics of biaxial nematic systems with discrete orientational degrees of freedom

We study a discrete version of a biaxial nematic liquid crystal model with external fields via an approach based on the solution of differential identities for the partition function. In the thermodynamic limit, we derive the free energy of the model and the associated closed set of equations of state involving four order parameters, proving the integrability and exact solvability of the model. The equations of state are specified via a suitable representation of the orientational order parameters, which imply two-order parameter reductions in the absence of external fields. A detailed exact analysis of the equations of state reveal a rich phase diagram where isotropic versus uniaxial versus biaxial phase transitions are explicitly described, including the existence of triple and tricritical points. Results on the discrete models are qualitatively consistent with their continuum analog. This observation suggests that, in more general settings, discrete models may be used to capture and describe phenomena that also occur in the continuum for which exact equations of state in closed form are not available.


Introduction
Mean-field models in Statistical Mechanics and Thermodynamics are a powerful tool to explore general qualitative properties of thermodynamic systems that, otherwise, would not be analytically treatable.Both conceptual and historical importance of mean-field models is testified by the celebrated van der Waals and Curie-Weiss models [1], complemented by Maxwell's equal areas rule (see e.g.[2]) which provided the first qualitative description of the mechanisms for the occurrence of phase transitions in fluids and magnetic systems.It is also well established that, in order to obtain accurate quantitative predictions, mean-field models need to be replaced by models with finite range interactions which are generally more challenging, and solvable cases require the use of sophisticated techniques as, for example, the transfer matrix and the renormalisation group, see e.g.[3].Spin models are the archetypal example of models aimed at describing the macroscopic and collective behaviour of systems made up of components with internal degrees of freedom (in the simple case the spin σ = ±1) with pairwise (and also higher order) interactions.Such models, although originally introduced in condensed matter physics to explain magnetic properties of materials are, however, of universal importance, as testified by applications in other disciplines such as Biology, Economics, Social Sciences, see e.g.[4,5,6] and references therein.It is also worth noting that a resurgence of interest, in the last decade, for spin-like mean-field models is due to the studies concerning their deployment for information processing, classification, memory retrieval and, more generally, machine learning purposes [7].These studies, originally inspired by the pioneering work of Hopfield [8], led to the definition of models for neural networks, such as the Boltzmann machines and their variations, based on spin glasses and statistical inference algorithms for training and learning [9].The key idea in this context is that spin particles sit at a node of a graph and possess internal degrees of freedom, i.e. their spin values are interpreted as node states of the neural network associated to the graph.The spin-spin interaction constant corresponds to the weight associated to the links on the network.In this paper, we consider a biaxial version of the discrete Maier-Saupe model for nematic liquid crystals (LCs) as studied in [10], whose structure resembles a multi-partite spin model with spin components subject to suitable constraints.The model consists of a system of particles endowed with an internal assigned geometry and symmetries with only orientational degrees of freedom.Not surprisingly, the exact analytical description of their macroscopic thermodynamic behaviour, phase transitions and emergent properties is, in general, not available and therefore alternative approaches and approximation techniques need to be adopted [11].Numerical simulations [12], Landau's expansion of the free energy [13,14], group representation and bifurcation theory [15] are approaches that allow to explore, at least locally, i.e. in the neighbourhood of specified values for the thermodynamic parameters, the possible occurrence of criticalities and phase transitions, and estimate relevant thermodynamic quantities such as orientational order, specific heat, critical exponents.Mean-field models are effective in providing insights that complement and support the aforementioned methodologies and all together help achieve accurate qualitative description and predictions on key properties of LCs including those that are paramount for technological applications [16].From a physical viewpoint, in the last few decades, the biaxial nematic liquid crystal phase has been the object of much intense study.The story of this phase has its roots back to 1970 [17], when the theoretical physicist Marvin Freiser noted that rather than possessing a rod-like shape (i.e. D ∞h symmetry), as usually assumed, most thermotropic, mesogenic molecules were in fact closer to being board-like, thus intrinsically biaxial (i.e. endowed with D 2h symmetry).Usually, they produce uniaxial nematic phases as a consequence of the rotational disorder around the long molecular axis, which eventually yields the definition of a single macroscopic director.This rotational disorder can be overcome by molecular mutual interactions favoring the molecules to align parallel to one another, thus leading to a thermotropic biaxial nematic phase at sufficiently low temperatures.Accordingly, Freiser understood that mesogens should be expected to exhibit a biaxial nematic phase, in addition to the usual uniaxial one.The prediction of a second nematic phase possessing novel properties and promising potential applications, stimulated considerable interest, as well as not little debate.In fact, on the experimental side, stable biaxial phases have been observed in lyotropic systems since the pioneering work of Yu and Saupe [18].In contrast, the experimental proof in favour of their existence in thermotropic systems has been subject of scrutiny and criticism in [19,20,21].In the period 1986 to 2003, the matter remained controversial with no widely accepted results [19,22,23].However, since 2004, clearer experimental evidence was provided for a few classes of compounds, such as polar bent-core or V-shaped molecules [24,25,26], and organosiloxane tetrapodes or their counterparts with a germanium core [27,28,29,30,31].These compounds have been investigated by several techniques which led to measurements of biaxial order parameters [32].According to these experimental results, an alternative picture of biaxial nematic order has emerged [33,34,35,36,37], based on the idea of biaxial domains reoriented by surface anchoring or external fields.Other researchers [38,39,40] have also pointed out that the biaxial nematic order is related to the onset of smectic fluctuations.Moreover, it has also been remarked that biaxial nematics may be formed from molecules possessing a lower symmetry than the usually assumed D 2h one, as for istance the C 2h symmetry [36,41,42,43,44,45].In addition, quite recently [46,47], low symmetry interaction models have been addressed, involving dipolar contribution, so as to describe polar bent-core molecules.The study of biaxial nematics is not only of theoretical origin, it is also connected with their potential technological applications in displays [33,48,49,50,51,52]: orientation of the secondary director in response to external perturbations is expected to be significantly faster than the primary one [34,41].Biaxial nematic phases have also been produced in colloidal suspensions of inorganic compounds [53,54,55].More recently, in [56], Smalyuhk et al. have considered a hybrid molecular-colloidal soft-matter system with orthorhombic biaxial orientational order and fluidity.This molecular-colloidal complex fluid is made up of only uniaxial rod-like building blocks.In contrast, this complex fluid exhibits a surprising self-assembly into a biaxial nematic liquid crystal with the D 2h point group symmetry.Finally, let us mention that very recently, the emergence of biaxial order upon mechanical strain has been proved experimentally in a nematic liquid crystal elastomer, the first synthetic auxetic material at a molecular level [57].By measuring the order parameters during deformation, the deviation from Maier-Saupe theory was detected for the uniaxial order parameters and the biaxial order parameters were deduced, suggesting the occurrence of biaxiality in the initially uniaxial system.On the theoretical side, after Freiser's first prediction [17], investigations were actively carried on along different approaches such as molecular-field or Landau theories, and later on by computer simulations.By the end of the past century, this collection of theoretical methodologies has shown that single-component models consisting of molecules possessing D 2h symmetry, and interacting via various continuous or hard-core potentials, are capable of producing a biaxial nematic phase under appropriate thermodynamic conditions [33,58,59,60].Theoretical studies usually predict a low-temperature biaxial phase, undergoing a transition to the uniaxial one, which, in turn, finally turns into the isotropic phase.In some cases, the transition takes place directly from the biaxial nematic to the isotropic phase.In the former cases, the ratio between the two transition temperatures (biaxial-to-uniaxial and uniaxial-to-isotropic) often turns out to be rather small in comparison with experimentally known stability ranges of the nematic phase.Both the isotropic-to-biaxial and uniaxial-to-biaxial phase transitions can be either firstor second-order, and, accordingly, the phase diagram exhibits triple and tricritical points.However, in the low temperature range, other phases, such as smectic or solid ones, may become more likely to occur.On the other hand, most theoretical frameworks only allow for isotropic and nematic phases [34], being the positional order not accounted for.Over the years, a rather simple, continuous, biaxial mesogenic pair interaction model has been proposed and investigated by several authors and via several types of techniques.In the literature, this model is known as the generalised Straley interaction [61] and it finds its roots in the celebrated Maier-Saupe model for interacting uniaxial nematic molecules [62,63,64].Actually, over the last two decades, several properties of this model have emerged, such as possible simplifications, additional symmetries and versatility in applications.More precisely, in 2003, new experimental findings on biaxial nematics boosted a renewed theoretical interest by some authors [65,66,67,68,69,70].More precisely, the generalised Straley pair potential model was studied by mean-field, as well as Monte Carlo simulation in the simple-cubic lattice-model version and, correspondingly, the effects produced on the resulting macroscopic behaviour were analysed [58,59,67,68,70,71].Moreover, motivated by the new experimental facts, the single-tensor Landau-de Gennes theory of biaxial nematics has been carefully revisited and a double-tensor Landau theory was put forward and studied [72,73].The hidden link between mean-field and Landau-de Gennes-type treatments has also been studied [74,75,76].The Straley potential model involves three independent parameters, and the aforementioned studies have shown that the model is rather versatile and capable of producing both biaxial and purely uniaxial order.In addition, the effect of strong antinematic terms, i.e. terms promoting misalignment, in the pair potential onto the resulting orientational order has been investigated [59,77].As shown in [77,78], these antinematic terms in the Straley model may destroy biaxiality, producing only uniaxial orientational order, and in some cases show evidence of the existence of a continuous ordering transition, in contrast with the discontinuous phase transition predicted by the simple Maier-Saupe model.Moreover, in [79], the Straley potential only contains antinematic terms, and it is found to produce biaxial order via a mechanism of order by disorder.In [61] the authors investigated the effect of two predominant antinematic couplings of equal strength perturbed by a comparatively weaker calamitic one.The resulting phases are a pure calamitic uniaxial phase, accompanied by an intermediate antinematic uniaxial phase.In this work, we consider a discrete version of the celebrated Maier-Saupe model for nematic LCs as the one considered in [10] and study its biaxial generalisation, i.e.Straley model, further extended to account for the effects of external fields.More specifically, molecules are assumed to be rigid cuboids, with two individual orientational degrees of freedom associated to two of the three principal axes of inertia, as the position of the third axis is automatically determined.It is also assumed that homologous principal axes of inertia interact pairwise for any pairs of molecules in the system.This assumption specifically characterises the mean-field models, where indeed any pair of molecules equally interact independently of their distance, and therefore positional degrees of freedom are not relevant.A further assumption is that orientational degrees of freedom are discrete, namely principal axes can only be parallel to the directions of a pre-defined Cartesian reference frame.The discretisation of orientational degrees of freedom for nematic liquid crystal models was firstly introduced by Zwanzig in [80] and successfully employed in various works, including recent papers [10,81].Although this assumption may seem to be at a glance restrictive, it captures, as observed in [10], with strikingly accuracy, properties of the continuum model.We show, via explicit examples, that the predictions obtained under specific symmetry reductions are consistent with the ones present in the literature for the corresponding continuum models.It is also important to note that, although, on one hand, the above assumptions restrict the model and allow to derive explicit global equations of the thermodynamic order parameters, on the other hand, the model is more general than its continuum analogues and the methodology adopted naturally incorporates external fields interacting with each orientational degree of freedom.Therefore, to the best of our knowledge, we provide the first theoretical study on the equilibrium statistical mechanics of a molecular field theory for biaxial liquid crystals subject to external fields.To solve the model, we show that the partition function Z N of the N −molecules reduced Straley biaxial model with external fields satisfies a remarkable differential identity as a function of the temperature and coupling constants.Using suitably re-scaled independent variables, the differential identity for the partition function of the finite size model is equivalent, up to a linear change of variables, to the heat equation.The required solution is therefore obtained by solving a linear equation with a specific initial condition that is fixed by the value of the partition function for the non interacting model the solution of which is straightforward.The properties of the system in the thermodynamic regime are obtained by studying the behaviour of the free energy in the limit as N → ∞, which corresponds to the semi-classical, or low diffusion, limit, of the heat equation, via a suitable asymptotic expansion of the free energy in powers of N −1 .At the leading order, the problem is solved via a Hamilton-Jacobi equation, which can be explicitly integrated and the solution is given in terms of the orientational order parameters from which the equations of state follow as a stationary point for the free energy functional.We study in detail the solution of the Hamilton-Jacobi equation and, specifically, a related system of quasilinear PDEs for the orientational order parameters.The methodology based on differential identities is applied for the first time in the context of a molecular theory for biaxial nematics described by two tensor fields.We derive explicit expressions for state functions when a finite number of molecules N is considered, as well as a novel system of equations of state, which include the interaction with external fields.We rigorously classify all admissible reductions in the absence of external fields, revealing a rich singularity structure describing transitions between isotropic, uniaxial and biaxial phases.A comparison with results available in the literature shows that our findings are consistent with those obtained with different methods and techniques.
As pointed out in a number of papers [82,83,84,85,86,87,88,10,89,90,91] the nature of the PDEs derived for the orientational order parameters suggests a natural interpretation of the singularities as classical shocks propagating in the space of thermodynamic variables.This allows to explain and, qualitatively, predict some features of the phase diagram based on the general properties of shock waves, as for example the occurrence of tricritical points as a collision and merging mechanism of two shock waves.This example demonstrates how such an interpretation is at the same time intriguing and of practical use.The paper is organised as follows.In Section 2 we introduce the physical model under study, we derive differential identities for the statistical partition function and we provide exact solutions for the model in the finite-size regime.In Section 3, we perform the thermodynamic limit and derive exact equations of state for the full model.Two-parameter reductions are also obtained in the cases of i) zero fields and ii) non-zero fields under special constraints.In Section 4 we present the phase diagram of the model in absence of external fields, and discuss criticality and behavior of the corresponding order parameters.Section 5 is devoted to concluding remarks.

The discrete λ-model for biaxial nematics
Let us consider a system of N interacting Liquid Crystals molecules with D 2h symmetry, whose molecular directors ⃗ m, ⃗ e and ⃗ e ⊥ are mutually orthogonal unit vectors parallel to their principal axes.The orientational state of a given molecule is identified by the directions of its molecular axes.Introducing the tensors (see e.g.[92]) where I is the 3 × 3 identity matrix, we consider the Hamiltonian of the form where q i and b i specify the orientational state of the i−th molecule and the scalar product is a • b := Tr (ab), where Tr is the trace operator.Summation indices i and j run from 1 to N , µ is the non-negative mean-field coupling constant and λ is a parameter weighing the degree of biaxiality.In the present paper, we assume λ ∈ [0, 1].In this range, the ground state for two interacting molecules corresponds to parallel homologous axes, that is e i tend to line up with e j , m i with m j and e ⊥,i with e ⊥,j .When λ = 0, the above Hamiltonian reduces to the classical Maier-Saupe model.The specific choice λ = 1/3 corresponds to the MMM model for liquid crystals with equally nematic interaction among corresponding molecular axes [67], i.e. (2.3) For convenience, we have included self-interaction terms corresponding to i = j.This choice will not affect the result as it corresponds to a shift of the energy reference frame by a constant.
The Hamiltonian (2.2) corresponds to a Straley pair-interaction potential reduced to the case of explicit zero-coupling between q and b tensors [64,65].Such pair-potential, identifying the socalled λ−model, has been studied extensively and its associated phase diagram, in the absence of external fields, has been inferred for specific two-order parameter reductions [65,66,67].It is worth noticing that, although the two tensors q and b are not directly coupled in the Hamiltonian (2.2), they are geometrically related via the constraint q i • b i = 0, thus implying an implicit microscopic coupling.As a result, the macroscopic behaviour will eventually reflect this hidden coupling through an entropic contribution in the free energy, in addition to other possible coupling terms in the order tensors, as we also show in this work.
Assuming that allowed configurations are such that molecular directors are parallel to the axes of a fixed Cartesian reference frame, the Hamiltonian (2.2) can be written as follows where c kl = 1 + δ kl for k, l = 1, 2, and Λ l i , with i = 1, • • • , N , and l = 1, 2, 3, 4, parametrise the components of q i and b i as follows giving six possible orientational states of each molecule.In particular, we have that for the i−th molecule, }, where Upon introducing the quantities M l = i Λ l i /N with l = 1, 2, 3, 4, the Hamiltonian H 0 reads as follows We now proceed with modelling the interaction between the liquid crystal and external fields.Consistently with previous studies on uniaxial [93,94,95,96,97] and biaxial nematics [98,10], we assume that the interaction between an individual biaxial liquid crystal molecule and external fields produces a term that is linear in the molecular tensors.Let ϵ = diag (ϵ 1 , ϵ 2 , ϵ 3 ) and χ = diag (χ 1 , χ 2 , χ 3 ) be two tensors associated with a general external field and let H ex be the Hamiltonian modelling the interaction between the external field and the liquid crystal molecules.Our assumption implies that H ex is of the form (2.9) Hence, the full Hamiltonian for the mean-field model under study in this work is The associated partition function for the Gibbs distribution is given by the expression where the summation refers to all possible configurations of (q i , b i ) and β = 1/T with T denoting the absolute temperature.Upon introducing the rescaled coupling constants t := βµ, x := βϵ 13 , y := βϵ 23 , z := βχ 13 and w := βχ 23 , the partition function reads as In the following, similarly to the case of van der Waals type models [85,87], spin systems [99] and the generalisation of the Maier-Saupe model in [10], we look for a differential identity satisfied by the partition function and calculate the associated initial condition.We observe that the partition function (2.10) satisfies the (4 + 1)-dimensional linear PDE Note that, for λ > 0 , equation (2.11) can be transformed via a linear transformation of the spatial coordinates into the heat equation where x ′ ,y ′ , z ′ , w ′ denote the new coordinates and σ = 1/N is the analogue of the heat conductivity.More precisely, the transformation of coordinates is given by u ′ = P λ u, where The associated initial condition, Z 0,N (x, y, z, w) := Z N (x, y, z, w, t = 0), corresponds to the value of the partition function of the model for non mutual interacting molecules.Given that the exponential is linear in the variables M 1 , M 2 , M 3 and M 4 , the initial condition can be evaluated by recursion and gives the following formula where the index i labels the quadruples The exact solution to the equation (2.11) for a given number of molecules N can be formally obtained by separation of variables using as a basis the set of exponential functions obtained by expanding the N −th power at the r.h.s. of equation (2.12).The solution reads as where Let us define the scalar order parameters m 1 N , m 2 N , m 3 N and m 4 N as the expectation values of, respectively, M 1 , M 2 , M 3 and M 4 , i.e.
Upon introducing the free-energy density as F N := (1/N ) log Z N , the order parameters for an intrinsically biaxial system composed by N molecules can be calculated by direct differentiation as follows where m l N = m l N (x, y, z, w, t; λ), for l = 1, 2, 3, 4. Equation (2.11) implies that the free-energy density satisfies the following differential identity (2.16) In Section 3, we derive the equations of state in the thermodynamic (large N ) regime via a direct asymptotic approximation of the solution to equation (2.16).Before proceeding, it is worth to emphasise that the case λ = 0 implies a reduction of the model (2.16) to the one studied in [10], although the initial condition considered in that work depends on the intrinsic molecular biaxiality parameter ∆, differently from the present case in which the degree of biaxiality of the interaction is entirely contained in the internal energy term.The differences in the two treatments arise as in this paper we are working with two order tensors, while in [10] the socalled geometric approximation on the interaction potential allowed to work with a single order tensor, that is a linear combination of q and b.
3 Thermodynamic limit and equations of state The thermodynamic limit is defined as the regime where the number of particles N is large, i.e.N → ∞.Under the assumption that the free-energy admits the expansion of the form F N = F + O (1/N ) and by using Eq.(2.16) we obtain, at the leading order, the following Hamilton-Jacobi type equation A similar asymptotic expansion for the order parameters m l N = m l +O(1/N ) implies the relations Equation (3.1) is completely integrable and can be solved via the method of characteristics.In particular, the solution can be expressed via the free-energy functional where m 1 , m 2 , m 3 and m 4 are stationary points of the free-energy, i.e.
Equivalently, order parameters are solutions to the following system of equations The term S(m 1 , m 2 , m 3 , m 4 ) represents the entropy of the system and, as discussed below, is uniquely fixed via the initial condition F 0 = F (x, y, z, w, t = 0).
The system (3.3)represents the set of equations of state for the λ-model.Hence, phase transitions can be studied through the analysis of critical points of the equations (3.3).Similarly to the thermodynamic models studied in [85,86,87,90], order parameters m l can be viewed as solutions to a nonlinear integrable system of hydrodynamic type where coupling constants x, y, z, w and t play the role of, respectively, space and time variables.In this framework, state curves within the critical region of a phase transition are the analog of shock waves of the hydrodynamic flow.In order to specify completely equations of state (3.3) we have to determine the function S(m 1 , m 2 , m 3 , m 4 ).We proceed by evaluating Eqs.(3.3) at t = 0, that is where m l 0 = m l (x, y, z, w, t = 0), with l = 1, 2, 3, 4. Equations (3.4) show that the function S(m 1 , m 2 , m 3 , m 4 ) can be obtained, locally, by expressing x, y, z and w as functions of the order parameters m l evaluated at t = 0 and then integrating Eqs.(3.4).Indeed, observing that the initial condition for F is F 0 = F N,0 = (1/N ) log Z 0,N , where Z 0,N is given in (2.12), the required functions can be obtained by inverting the system ∂F 0 ∂w (x, y, z, w) .
(3.5)More explicitly, equations (3.5) read as follows where we have introduced the notation X = exp(x), Y = exp(y), Z = exp(z), W = exp(w).Hence, equations of state (3.3) for the model with external fields are completely determined in terms of the roots of system of equations (3.6).We should also emphasise that system (3.6) is algebraic with respect to the variables X, Y , Z and W .
Remark.The order parameters introduced here are related to the scalar order parameters adopted in [65] by the following linear transformation where S, T, S ′ , T ′ are the scalar order parameters characterising the tensors Q := ⟨q⟩ and B := ⟨b⟩ in their common eigenframe, once the thermodynamic limit is performed.Specifically, by considering the eigenframe (⃗ e x , ⃗ e y , ⃗ e z ), the order tensors can be written as The inverse of the linear transformation (3.7) is In [67] it is claimed that, in the absence of external fields, reductions T = S ′ = 0, or T = ±S and S ′ = ±3T ′ hold, these latter obtained by swapping the axes of the reference frame, ⃗ e x , ⃗ e y , ⃗ e z .These conditions read as m 1 = m 2 = −S/3 and m 3 = −m 4 = T ′ .
In the next section, we will introduce a new parametrisation based on the introduction of the molecular Gibbs weights, which leads to the explicit solutions of the model.

Equations of state
A convenient approach to the evaluation of the entropy of the discrete model and the corresponding equations of state starts from the statistical analysis of the 'initial condition', namely the evaluation of the partition function (2.12) as a function of the external fields at t = 0. Indeed, at t = 0, liquid crystal molecules are mutually independent and expectation values can be evaluated by looking at the one-molecule partition function, e xΛ 1,i +yΛ 2,i +zΛ 3,i +wΛ 4,i . (3.11) The molecular Gibbs weights [2] at t = 0 and as functions of the external fields take the following form p 0,i (x, y, z, w) := e xΛ 1,i +yΛ 2,i +zΛ 3,i +wΛ 4,i Z 0,1 , i = 1, . . ., 6 .

Two-parameter reductions
In this subsection, we will focus on the derivation of two-parameter reductions of the equations of state (3.17a)-(3.17d).Such reductions arise naturally when considering the liquid crystal system constrained to suitable forms of external fields, including the case in which external fields are not present and the phase behaviour is entirely regulated by the mutual interactions among liquid crystal molecules and the temperature.The following holds in the absence of external fields.
Lemma 3.1.In the absence of external fields, that is at x = y = z = w = 0, solutions to the system (3.17a)-(3.17d)are given by one of the following 2-parameter reductions: i) a) p 4 = p 2 and p 3 = p 1 , with Let k > 0 be an arbitrary constant and ϑ be the scale transformation defined by ϑ :  As we prove in Theorem 3.1, Lemma 3.1 has a remarkable implication on the structure of the two order tensors of the theory.In order to proceed, it may be convenient to recall a criterion to characterise the degree of biaxiality of a given order tensor Ω.This will be based on the biaxiality parameter Furthermore, in the extreme case β 2 (Ω) = 1, Ω is said to be maximally biaxial.
The following theorem characterises the allowed forms of the two order tensors of the model.
Theorem 3.1.In the absence of external fields, at all temperatures and values of λ, the order tensors take one of the following two forms i) Q uniaxial and B maximally biaxial; ii) Q and B both uniaxial.
Proof.The result is readily obtained by considering Lemma 3.1 and the transformation φ specified by Eqs.(3.16).The subcases a), b) and c) of Lemma 3.1 in each of the two cases i) and ii), correspond to a particular choice of the principal axis.For instance, the transformation φ evaluated along with case i) a) implies , m 4 = 0 , The statement is proven by evaluating the biaxiality parameter β 2 for Q and B along all cases.Due to the invariance by exchange of principal axes, β 2 for Q and B will take same values for all subcases a), b) and c) of a given case.Without loss of generality, we can consider cases i) a) and ii) a) to get, respectively, i) A direct consequence of the reductions ii) in Theorem 3.1 and the transformation φ is the following corollary.
Corollary 3.1.The equations of state for the model in the case of Q and B both uniaxial, cases ii) in Theorem (3.1), can be written explicitly in terms of the eigenvalues m l .In particular, we have that reductions ii) a), b) and c) can be written as follows ii) a) m 2 = m 1 and m 4 = m 3 with Proof.As shown in the proof of Theorem (3.1), the transformation φ is linear when restricted to the case of Q and B both uniaxial.Hence, the transformation can be easily inverted to get the projection φ −1 : (m 1 , m 2 , m 3 , m 4 ) −→ (p 1 , p 2 , p 3 , p 4 ) along with each particular 2-parameter reduction.The equations in terms of the eigenvalues are then obtained by application of the inverse transformation for the specific reduction to the corresponding set of equations in the p−variables.Taking the case ii) a) as an example, the application of Unlike uniaxial-uniaxial reductions discussed above, uniaxial-maximally biaxial reductions cannot be written in explicit simple form in terms of the m l variables.
In this paper, we will focus our discussion on the phase behaviour in the absence of external fields.We note, however, that the above reductions are also compatible with non-zero fields values subject to suitable constraints.Indeed, the following proposition allows to identify the constraints on the external fields so that the system admits uniaxial-maximally biaxial and uniaxial-uniaxial solutions for Q and B. In such cases we can still consider 2-parameter reductions of the system, with the equations of state also accounting for the action of applied fields.

Order parameters in the 2-parameter reductions
The equations of state (3.17a)-(3.17d)are the critical points of the free-energy functional (3.18).According to the definition of free-energy adopted in this paper, global maxima identify stable states of the system and associated phases.Coexistence curves (hypersurfaces, in general) arise as sets of control parameters for which two or more local maxima are resonant, hence identifying the coexistence of the corresponding phases.In order to proceed, we should therefore first identify all local maxima for each choice of control parameters x, y, z, w, λ and t.In this section, we focus on the complete characterisation of the system in the absence of external fields, i.e.
x = y = z = w = 0, hence relying on Theorem 3.1 and its implications.An immediate consequence is that the onset of phase transitions is determined by analysing the singularities of 2-dimensional maps defined by Eqs.(3.19)-(3.26)(see [10] for an exhaustive treatment).Noticeably, this is a more affordable task compared to the full 4-dimensional problem governing the system when external fields are present, Eqs.(3.17a)- (3.17d).By evaluating the Hessian matrix of the free energy density (3.18), it turns out that none of the critical points in the uniaxial-uniaxial reductions, cases ii)a)-c) in Lemma 3.1) are stable.This result is consistent with what is known from mean-field theories based on a continuum of molecular orientational states [67].Therefore, the uniaxial-maximally biaxial reductions, cases i)a)-c) in Lemma 3.1, are the only ones relevant from the equilibrium thermodynamics viewpoint.
The following subsections focus on the analysis of the resulting phase diagram and the associated order parameters behaviour, with case i)a) being considered for such purpose.Cases i)b) and i)c) can be straightforwardly obtained from case i)a) via suitable linear transformations on m 1 and m 3 , which merely correspond to permutations of axes.The λ-t plane is divided in three regions identifying three distinct macroscopic phases, namely the isotropic (I), the uniaxial nematic (U) and the biaxial nematic (B).The lines separating the different regions are either dotted black lines or solid black lines.The former identify the socalled second-order transition lines, that is the lines associated with phase changes characterised by continuous order parameters but discontinuous derivatives.The latter are instead associated to first-order lines, that is the order parameters and their derivatives experience a discontinuity when the line is crossed.Similarly to the analysis performed in [10], second-order lines are identified by cusp points of two-dimensional maps.The cusp points of the model (red lines in Fig. 1) are given explicitly in terms of the transcendental curve Notice that the cusp set can be seen as the union of two curves intersecting at the point (λ, t) = (1/3, 3).The model admits two tricritical points, (λ = (2/3, 3/2) (blue circle).The three phases coexist at the triple point, (λ tp , t tp ) = (0.234, 2.773) (green circle) identified by the resonance condition for the corresponding maxima.A closer look in the region surrounding the triple point and the uniaxial-biaxial tricritical point is provided in the right column.The cusp points in the λ − T * plane are given by the set The constraint λ = T * identifies the subset of cusp points associated to second-order lines for λ ≥ 2/3.The behaviour for values of λ in the interval λ For values of λ in this range the model predicts two first-order phase transitions, a isotropic-to-uniaxial phase transition at high temperature (low values of t) followed by a uniaxial-to-biaxial phase transition at lower temperatures (higher values of t).While the former is associated to a shock that is static in λ, the latter is originated at the uniaxial-biaxial tricritical point and is identified by a classical shock whose location moves from low temperatures to higher temperatures as the biaxiality parameter λ increases.
As shown in Figs. 1 and 4 (left column), for λ = λ tp = 0.234 the system displays a triple point at which all three phases coexist.This situation is realised as the two shocks associated with the isotropic-to-uniaxial and uniaxial-to-biaxial phase transitions merge at zero external fields, giving rise to a single shock having an amplitude given by the sum of the amplitudes of the two individual shocks.Consistently with the phase diagram in Fig. 1, the uniaxial phase is not energetically accessible for λ > λ tp .For instance, for λ = 0.284 (centre column of Figure 4), the order parameters jump from the isotropic phase to the biaxial phase as the temperature is lowered.This is also the case when λ = 1/3 (right column).The case λ = 1/3, as also discussed in [67], leads to proportionality between the stable branches of the order parameters.Precisely this is 3m 1 ± m 3 = 0, corresponding to T ′ = ±S in the convention adopted by Virga  , λ tp as an example.The right column displays a magnification of the order parameters around the tricritical temperature and triple point temperature.and co-authors in [67].
For values of λ exceeding the value 1/3, the isotropic-to-biaxial phase transition remains first-order until a second tricritical point is disclosed.In Fig. 5, the change in order of the phase isotropic-to-biaxial phase transition is displayed.For λ < λ (IB) tc (left column) both order parameters undergo a discontinuous jump from the isotropic solution to the biaxial one.The shock disappears when λ = λ (IB) tc = 2/3 is considered, and consequently both order parameters experience a gradient catastrophe at t = t (IB) tc = 3/2.Larger values of λ lead to a direct second-order transition from the isotropic to biaxial phase, with the transition value given by t (IB) = 1/λ.
Our results are, both qualitatively and quantitatively, consistent with previous studies [58,67,60].In particular, according to the Monte Carlo simulation results reported in [58,67,60], the values of λ at the first and second tricritical points are ≃ 0.24 and ≃ 2/3, respectively.Moreover, the global Monte Carlo study performed in [60] predicts λ ≃ 0.26 for the triple point.Consistency is also shown with the extended analysis preformed in [73], where a fourth degree Landau potential in the two tensors Q and B is considered.Indeed, the Authors find a zero-fields phase diagram displaying a first-and second order isotropic-to-uniaxial-to-biaxial phase transitions and first-and second-order isotropic-to-biaxial phase transitions, thus disclosing two tricritical points and a single triple point.On the other hand, similar but not totally equivalent phase diagram topologies are described in [72,14], where a sixth degree Landau potential in a single order tensor has been analysed.The phase diagram we obtain (Figure 1) is also qualitatively in agreement with the one obtained in [74], where the Authors derive a Landau expansion of a free energy which is intrinsically linked to a molecular-field theory, and then discuss the Sonnet-Virga-Durand limit.
It goes without saying that a single-tensor theory does not distinguish between intrinsic and phase biaxiality.This separation is made clear and sharp by starting from two molecular tensors q and b in the Hamiltonian (e.g.(2.2)) which, correspondingly, give rise macroscopically to two order tensors Q and B eventually accounting for the phase and intrinsic biaxiality, respectively.In this work, the λ−model in the absence of external fields is shown to admit two-parameter reductions, which account for isotropic, uniaxial and intrinsic biaxial phases.Other models, as for example the Maier-Saupe model (retrieved from the λ−model setting λ = 0), have been shown to produce uniaxiality and phase biaxiality at the macroscopic level and in the presence of external fields [10,95,96,97].

Concluding remarks
In this paper we have analysed in detail a discrete mean-field model for a biaxial nematic liquid crystal subject to external fields, using an approach based on the differential identity (2.11) for the partition function.Upon the introduction of suitable variables, namely the order parameters, the multidimensional linear PDE satisfied by the partition function leads, in the thermodynamic limit, to a set of equations of state involving all four orientational order parameters.The equations are completely solvable by the method of characteristics proving the integrability of the model.
Via the introduction of a novel set of order parameters corresponding to orientational Gibbs weights, we have obtained the equations of state in explicit form.We proved that, in the absence of external fields, the system is fully characterised by two-parameter reductions, and such reductions persist in the case of non-zero external fields subject to suitable constraints.A detailed analysis demonstrates the existence of a rich phase diagram, that is remarkably consistent with the results known in the literature for the standard Maier-Saupe model and its biaxial extensions.Hence, the discrete models of the type studied in this paper capture, at least qualitatively, the most important features of continuum models with external fields for which explicit analytic formulae are not available.These results indeed encourage further studies on integrable biaxial models where the Hamiltonian contains a more general nonlinear dependence on the tensors q and b, as for instance the one implied by the full Straley pair-interaction potential.Moreover, the phenomenology encoded in the equations of state (3.17a)-(3.17d)for general field values, as well as the 2-parameter reductions obtained in Proposition 3.1.1for constrained fields, still need to be further analysed and described.Such cases are currently under investigation and results will be reported in due course.

p 3 p 4 ( 3 . 50 )z + 3 (p 3 − p 4 )
The constraints on the fields follow from looking for either uniaxial-maximally biaxial or uniaxial-uniaxial reductions of the whole set of equations of state, Eqs.(3.17a)-(3.17d).For instance, the uniaxial-maximally biaxial reduction i) a) requires p 3 = p 1 and p 4 = p 2 .Eqs. (3.17a)-(3.17d)restricted to this constraint imply that compatibility of the first two equations restricts fields to y = x, while compatibility of the third and fourth requires w = −z.Elementary algebraic manipulations lead to Eqs. (3.40)-(3.41).The set of field-dependent equations of state in the case i) b) and c), and ii) a), b) and c) are obtained following the same procedure, together with the associated constraints on the fields.

4. 1
Phase behaviour in the absence of external fields The phase diagram of the model in the absence of external fields is shown in Fig 1.The top row shows the phase diagram in the λ − t plane, while the bottom row shows the phase diagram in the λ − T * plane, where T * is the dimensionless temperature defined by T * := 1/t = (k B T )/µ.

Figure 1 :
Figure 1: Zero-fields phase diagram.T op row: phase diagram in the λ − t plane.B ottom row: phase diagram in the λ − T * plane.The figures in the right column represent a magnification of the phase diagram in the region surrounding the triple point (green circle) and the uniaxialbiaxial tricritical point (red circle).The line associated with uniaxial cusp points is indicated in red.

4. 2 ( 3 Figure 2 :
Figure 2: Order parameters in the 2-parameter reduction for small values of λ.Each column shows both order parameters, m 1 (black) and m 3 (red) versus t at a specific value of λ.The solutions corresponding to a global maximum of the free energy are displayed with solid lines, while other solutions are indicated with dotted lines.Left column: λ = 0, that is the uniaxial Maier-Saupe interaction potential.Centre column: λ = 1/6.Right column: λ = λ U B tc = 0.217.

Figure 3 :
Figure 3: Order parameters in the two-parameter reduction for values of λ in the interval λ (U B) tc ≤ λ < λ tp .The two column shows both order parameters versus t for λ = 0.226 ∈