We employ classical density-functional theory to investigate the phase diagram of an assembly of mutually penetrable, parallel ellipsoids interacting via the generalized exponential model of index four (GEM-4) pair potential. We show that the crystal phases of the system are obtained from those of the spherically symmetric GEM-4 model by rescaling the lattice vectors. Performing this rescaling in combination with an arbitrary rotation of the lattice leads to infinitely many different structures with the same free energy, thereby implying their infinite degeneracy. These findings generalize to non-zero temperature the results formerly obtained by us [Pini *et al.*, J. Chem. Phys. **153**, 164901 (2020)] for the ground state of a similar system of ellipsoids interacting via a Gaussian potential. According to the mean-field free-energy functional used here, our conclusions apply to soft-core potentials both when they form cluster crystals as the GEM-4 and when they form single-occupancy crystals as the Gaussian itself.

## I. INTRODUCTION

Bounded, repulsive potentials that allow particle overlap are often introduced to model the effective interactions between non-compact macromolecules featuring an open internal structure, such as polymers^{1} and dendrimers.^{2,3} For simplicity sake, those potentials are generally assumed to be spherically symmetric. However, there is strong evidence that the macromolecules which they are meant to represent are often aspherical (see, e.g., Refs. 4–11) and would be then more realistically described by anisotropic interactions.

In an effort to introduce the asphericity of the interaction into the picture without the increased complexity entailed by the need of taking into account the orientational degrees of freedom of each (effective) particle, some time ago, a simplified model potential was proposed,^{12} which in the following will be referred to as soft-core nematic.

According to this model, particles are represented as soft-core ellipsoids of revolution and the pair interaction depends not only on their mutual distance but also on the angle between their orientation and the vector connecting their centers. However, the orientation is assumed from the outset to be the same for all particles. Clearly, this rules out the possibility of employing such a model to study the onset of orientational order. Its purpose is, instead, to describe the effect of asphericity on the behavior of the system, *assuming* that a nematic state has been reached.

In a former publication,^{13} we employed classical density-functional theory to investigate the crystal phases of the soft-core nematic model when the soft-core interaction is given by the generalized exponential model of index four (GEM-4).^{14,15} We showed that the GEM-4 nematic interaction, just like its spherically symmetric counterpart, features crystals with multiple occupancy of the lattice sites and a nearly state-independent lattice constant. The local density displays sharp peaks at the lattice sites, which, not surprisingly, on increasing the density or the degree of asphericity rapidly assume the same ellipsoidal shape and the same orientation as those of the particles. We also found that the system presented a wealth of different crystal structures, which, however, were not extensively discussed.

In a subsequent paper,^{16} we turned to the soft-core nematic model for a pair interaction of Gaussian form, which had already been considered by other authors.^{12,17} In Ref. 16 we were concerned only about the system at zero temperature, when its behavior is determined solely by its energy. In that case, we showed that (i) by rescaling the lattice vectors along the direction of the main axis of the particles, the crystal phases of the Gaussian nematic model are mapped rigorously into those of the spherically symmetric Gaussian potential at a rescaled density and (ii) as a consequence, there exists an infinite number of different lattice structures, which have the same energy. This holds, in particular, for the ground state, which is then infinitely degenerate.

In the light of those findings, in the present paper, we reconsider the GEM-4 nematic model by the same mean-field grand potential functional used in Ref. 13 and show that, according to that functional, at every temperature its crystal phases are related to those of the spherically symmetric GEM-4 by the same mapping described in Ref. 16, thereby implying their infinite degeneracy. On this basis, those phases are classified by the lattice of the spherically symmetric system from which they are obtained, and the phase diagram is determined. These results are fully supported by the direct numerical minimization of the functional.

The generalization at non-zero temperature of the mapping of Ref. 16 can in fact be applied not only to the GEM-4 nematic model but also to the Gaussian nematic model considered there. However, it should be remarked that in both systems, this generalization hinges on the mean-field form of the grand potential functional and therefore might not hold exactly, unlike in the zero-temperature case.

The plan of this paper is as follows: in Sec. II, we briefly describe the soft-core nematic model for GEM-*m* pair interactions. In Sec. III, we introduce the mean-field grand potential functional used in the paper and prove that the ordered configurations of the soft-core nematic model are mapped into those of the corresponding soft-core model with spherically symmetric interactions. In Sec. IV, this result is applied to the GEM-4 nematic model, and its phase diagram is determined. A somewhat technical point is addressed in Sec. V, where the grand potential functional introduced in Sec. III is modified so as to remove the self-interaction term, and it is proved that doing so does not affect the aforementioned mapping. This is relevant for the application to the Gaussian nematic model considered in Sec. VI, where the results of Ref. 16 for the ground-state configurations are extended to non-zero temperature. Our conclusions are drawn in Sec. VII. In Appendix A, evidence of the aforementioned lattice degeneracy is given by reporting the numerical results obtained by the minimization of the grand potential functional for two different thermodynamic states of the GEM-4 nematic model. Appendix B provides a detailed discussion of some technical points related to the Gaussian potential mentioned in Sec. VI, which are not specific to the nematic model.

## II. MODEL

*m*pair potential is given by

*r*is the modulus of the center-to-center vector

**r**of two interacting particles,

*σ*

_{0}sets the unit length,

*ɛ*sets the energy scale,

*m*is a positive integer with

*m*> 2, and

*φ*(

*x*) = exp(−

*x*

^{m}). Since Φ

_{0}(

*r*) assumes a finite value at

*r*= 0, it belongs to the class of soft potentials.

^{14,15}

*m*nematic model considered in this contribution, the GEM-

*m*interaction is generalized to aspherical particles by introducing the pair potential

*λ*is a parameter that characterizes the asphericity of the particles such that

*λ*= 1 corresponds to the spherically symmetric potential (1), and

*ϑ*is the angle enclosed by $r\u0302$ and an arbitrary, but fixed unit vector

**u**(the so-called director), which determines the orientation of the particles, and is assumed to be the same for all particles. Without loss of generality, we identify the

*z*axis with the direction of

**u**; thus, cos

*ϑ*=

*z*/

*r*, and Eq. (2) becomes

*z*axis with aspect ratio

*λ*.

**r**),

*z*axis with aspect ratio 1/

*λ*.

The system is characterized by its number density *ρ* and temperature *T*. As usual, we shall denote by *β* the quantity $(kBT)\u22121$, with *k*_{B} being the Boltzmann constant.

## III. DENSITY-FUNCTIONAL THEORY FOR SOFT, ASPHERICAL PARTICLES

^{18–20}the grand potential of a system can be expressed in the grand-canonical ensemble as a functional Ω[

*ρ*] of the one-particle density (or density profile)

*ρ*(

**r**), whose expression at equilibrium is determined by minimizing Ω[

*ρ*]. Assuming the mean-field approximation, which, according to Refs. 14 and 15, is appropriate for systems with cluster-forming soft interactions, such as the GEM-

*m*with

*m*> 2,

^{15}the grand potential functional of our system is given by

*μ*is the chemical potential and Λ is the thermal wavelength. If

*ρ*(

**r**) is periodic, as we assume throughout this contribution, then one has

**a**

_{i},

*i*= 1, 2, 3, form a set of vectors defining the Bravais lattice.

^{21}For convenience, we collect these vectors and the related vectors of the respective reciprocal lattice,

**b**

_{i},

*i*= 1, 2, 3, in matrices

**A**and

**B**,

**r**

_{m}and

**k**

_{n}, can be expressed as

**r**

_{m}= ∑

_{i}

*m*

_{i}

**a**

_{i}and

**k**

_{n}= ∑

_{i}

*n*

_{i}

**b**

_{i}with the integers

**m**= (

*m*

_{1},

*m*

_{2},

*m*

_{3}) and

**n**= (

*n*

_{1},

*n*

_{2},

*n*

_{3}).

*V*is the volume of the system and

*v*is the volume of the unit cell $C$ of the Bravais lattice. The sum runs over the vectors of the reciprocal lattice

**k**

_{n}, and $\rho \u0302(n)$ is given by

**r**inside the unit cell can be expressed as

**r**=

**A**·

**s**, where

**s**is a vector of components

*s*

_{α}, which vary in a cube $Q$ of unit edge, −1/2 <

*s*

_{α}≤ 1/2, with

*α*= 1, 2, 3. For a given lattice, the density profile is then determined by the function

*ψ*(

**s**) defined in $Q$ such that

*k*

_{λ,n}is given by Eq. (9) with

**k**=

**k**

_{n}, and Eq. (14) becomes

According to Eqs. (16) and (17), both the ideal-gas contribution to the grand potential and the Fourier transform $\psi \u0302(n)$ of the density profile are actually independent of the shape of the lattice: they solely depend on the values taken by the density profile inside the unit cell, whichever the lattice generated by it. The information on the lattice is conveyed *only* via the vectors **k**_{λ,n} at which $\Phi \u03030(k)$ is evaluated.

**r**

_{λ,m}and

**k**

_{λ,n}can also be regarded as vectors of a Bravais lattice with primitive vectors $ai0$ and $bi0$,

*i*= 1, 2, 3,

**A**and

**B**via

**D**

_{λ},

*ψ*

_{0}(

**s**) given by

*λ*. By expressing

*ψ*(

**s**) in terms of

*ψ*

_{0}(

**s**) and taking Eq. (18) into account, we find

*ψ*

_{0}(

**s**) according to Eq. (17).

*β*Ω

_{0}[

*ρ*

_{0}]/

*V*of a system with the spherically symmetric potential (1) and the reciprocal lattice generated by

**B**

_{0}, at the same inverse temperature

*β*as that of the original system and fugacity

*z*

_{0}given by

*z*=

*e*

^{βμ}is the fugacity of the original system. Clearly, Ω is minimized by

*ψ*(

**s**) if and only if Ω

_{0}is minimized by

*ψ*

_{0}(

**s**). Hence, for the grand potentials, density profiles, and average density $\rho \u0304$ at equilibrium, it holds that

*F*, the relationship between the Helmholtz free energy

*F*of the aspherical potential and that of the spherical potential

*F*

_{0}is readily established from Eqs. (23) and (24),

*N*is the number of particles.

*T*, fugacity

*z*, and average density $\rho \u0304$ are mapped

*exactly*into those of the spherical potential at the same temperature

*T*, fugacity

*λz*, and average density $\lambda \rho \u0304$, i.e.,

**A**

_{0}and

**B**

_{0}are known, and the direct and reciprocal lattices of the aspherical case are obtained by Eq. (19).

**A**

_{0}by an arbitrary rotation matrix

**R**. Clearly,

**A**

_{0}and

**RA**

_{0}represent the same kind of Bravais lattice, the rotation amounting just to a change of the reference frame. Hence, the matrices

**A**=

**D**

_{λ}

**A**

_{0}and

**A**′ =

**D**

_{λ}

**RA**

_{0}both describe an equilibrium configuration of the aspherical potential at the same

*T*and

*z*, giving the same grand potential Ω. In both cases, the spherical neighbor shells of the direct and reciprocal lattices of the spherically symmetric system are mapped into ellipsoids,

**A**and

**A**′, in general, do

*not*coincide: their points, even though they lie on the same ellipsoids, have different distances from the centers of the ellipsoids so that their neighbor shells are different. A pictorial representation of this situation was given in Fig. 1 of Ref. 16. Since the rotation

**R**is arbitrary, there are infinitely many different Bravais lattices, which, for a given state, correspond to the same Ω, i.e., the periodic configurations of the system are infinitely degenerate. These degenerate configurations correspond to the same

*ψ*(

**s**), which, however, is attached to different lattices.

**A**. We might then wish to establish if this lattice is actually obtained by applying to the Bravais lattice of the spherical potential corresponding to

**A**

_{0}the rotation-rescaling procedure just described. This is certainly the case if

**A**=

**D**

_{λ}

**RA**

_{0}for some rotation

**R**, but, in fact, this condition should be relaxed because the primitive vectors of a Bravais lattice are not univocally determined. As discussed in Ref. 16, two matrices

**C**and

**C**′ generate the same Bravais lattice if and only if

**C**′ =

**CL**, where

**L**is a matrix with integer elements such that det

**L**= ±1, the identity being of course a special case. Hence,

**A**belongs to the family of degenerate lattices obtained from

**A**

_{0}if and only if there exist a rotation

**R**and a matrix

**L**with the above properties such that

**R**is represented by its Euler angles, the group of rotations can be swept numerically, and it can be easily checked whether Eq. (31) holds for some

**R**and given

**A**and

**A**

_{0}.

## IV. APPLICATION TO THE GEM-4 NEMATIC MODEL

We now specialize the general results discussed in Sec. III to the GEM-4 nematic model considered in Ref. 13. For the isotropic GEM-4 interaction such that $\Phi 0(r)=\epsilon exp[\u2212(r/\sigma 0)4]$, the phase diagram is known^{14,15} and features both bcc and fcc cluster crystals.

**A**

^{bcc}and

**A**

^{fcc}of the primitive vectors of the corresponding lattices of the aspherical system obtained via Eq. (19) are given by

*k** of the reciprocal lattice of the isotropic system. The reason for this choice is that

*k** can be identified, to a very good degree of approximation, with the wave vector

*k*

_{min}of the minimum of $\Phi \u03030(k)$.

^{14,15}Specifically, for the GEM-4 potential, one has

*k** ≃ 5.572/

*σ*

_{0}. However, the above results do

*not*hinge on that approximation.

Equations (32) and (33) represent, in fact, one instance of the infinitely many degenerate lattices obtained from Eq. (30) when **A**_{0} is replaced by either $A0bcc$ or $A0fcc$, respectively. In the following, we shall refer to those classes of lattices as bcc* and fcc*.

The conclusions of the above picture were compared to the numerical results obtained by the direct minimization of the grand potential functional (10). The algorithm is the same one first adopted in Refs. 22 and 23 and subsequently in Ref. 13 for the model studied here. It performs a fully numerical minimization based on the conjugate-gradient method, whereby the optimization is carried out with respect to both the density profile *ρ*(**r**) and the lattice vectors **a**_{i}, with no assumption on the analytical form of *ρ*(**r**), other than it being periodic. To this end, *ρ*(**r**) was discretized by sampling it at 40^{3} = 64 000 grid points. This discretization, while being coarser than those adopted in Refs. 13, 22, and 23, is amply sufficient for the purpose of the present paper. The asphericity parameter *λ* was fixed at *λ* = 1.5, which is the value most often considered in Ref. 13.

As in former studies,^{13,22,23} *μ* was parameterized by introducing the density *ρ*^{⋆} related to *μ* by Eq. (35). Clearly, in the homogeneous phase, *ρ*^{⋆} coincides with the average density $\rho \u0304$, but this is no longer the case if the system is inhomogeneous, when $\rho \u0304$ is obtained *a posteriori*. Even so, the benefit of this parameterization is that the thermal wavelength in Eq. (10) cancels out since it appears both in ln[*ρ*(**r**)Λ^{3}] and in $ln[\rho \u0304\Lambda 3]$. Moreover, it comes handy in the present situation because switching from the aspherical system at fugacity *z* to the spherically symmetric system at fugacity *λz* is equivalent to switching from *ρ*^{⋆} to *λρ*^{⋆}.

The numerical minimization fully confirms the mapping described in Sec. III. In particular, we checked that the numerical results for the grand potential Ω at fugacity *z* are related to the grand potential Ω_{0} at fugacity *λz* by Eq. (24) and that the respective densities $\rho \u0304$ and $\rho \u03040$ are related by Eq. (26). At given *T* and *z*, minimization runs initialized by different trial density profiles gave the same Ω for different Bravais lattices, whose primitive vector matrices **A** rigorously fulfill Eq. (31) for either $A0=A0bcc$ or $A0=A0fcc$. Some of these numerical results are presented in greater detail in Appendix A for two different states.

In the light of the results above, the phase diagram of the GEM-4 nematic model in the temperature–density plane is obtained straightforwardly from that of the spherically symmetric GEM-4 displayed in Fig. 4 of Ref. 22. One needs just to rescale the density axis by the factor 1/*λ* and replace the bcc and fcc domains by the bcc* and fcc* ones.

The phase diagram thus obtained is shown in Fig. 1. As one moves away from the fluid region by increasing the density at constant temperature, one meets a first-order transition from the fluid to (any of the degenerate lattices of) the bcc* phase and subsequently another first-order transition from the bcc* to the fcc* phase. Note that the phase diagram of the spherically symmetric GEM-4 reported in Ref. 14, where the fluid phase was described by the energy route of the random-phase approximation (RPA), presents a fluid–bcc–fcc triple point. Below the triple-point temperature *k*_{B}*T*_{t}/*ɛ* ≃ 0.4, the bcc phase disappears. As stated above, here as well as in Ref. 22, the fluid instead has been described by the (less accurate) van der Waals approximation obtained from Eq. (10) in the homogeneous phase, which is equivalent to the RPA compressibility route. As a consequence, we did not find any triple point, in either the spherical or aspherical case.

Since in the present model the particles have the same orientation by construction, the fluid phase is obviously nematic. At the same time, smectic phases, which would occur in more realistic models of aspherical particles, are not expected here. A smectic phase is periodic along a single direction so that *ρ*(**r**) should represent a one-dimensional, stripe-like modulation of the local density. According to the mapping described in Sec. III, this phase would be necessarily obtained from a similar stripe phase of the spherically symmetric GEM-4. However, at variance with hard-core interactions giving rise to mesophases, such as the short-range attractive, long-range repulsive (SALR) potential,^{23} GEM-4 does not yield a stable stripe phase, nor a stable tubular phase for that matter. As discussed in detail in Ref. 24, although at sufficiently low temperature or high density those phases are favored with respect to the homogeneous fluid, they are always unfavored compared to phases characterized by three-dimensional density modulations, such as the bcc or fcc. This is fully supported by numerical simulation,^{14} which gave no evidence of inhomogeneous phases of lower dimensionality.

*w*

_{0}is the absolute value of $\Phi \u03030(k)$ at its minimum. For the GEM-4 potential, one has $w0\u22430.127\epsilon \sigma 03$. At lower temperature or larger density, the homogeneous phase is no longer a minimum of the grand potential at given

*T*and

*μ*, not even locally. This locus is generally referred to as the

*λ*-line (where the symbol

*λ*holds no relation with the asphericity parameter

*λ*used throughout the paper).

## V. SUBTRACTING OFF THE SELF-INTERACTION TERM

As discussed in Sec. III, the infinite degeneracy of the configurations of the GEM-*m* nematic model is based on the mapping on the configurations of the corresponding spherically symmetric system. Within the mean-field grand potential functional (10), this mapping holds exactly. However, the functional itself is approximate and disregards correlations. Hence, one may wonder whether such a result would still hold, if at least some of the correlation effects neglected in that functional were taken into account.

**r**

_{m}with a function

*f*(

**x**) discussed below. Equation (40) is used to describe the crystal phases of the spherically symmetric GEM-4

^{14,15}and also other so-called

*Q*

^{±}interactions,

^{25}such as the double-Gaussian potential

^{26}introduced in Ref. 3. In those cases, a Gaussian expression for

*f*(

**r**) was assumed, but, in fact, any periodic

*ρ*(

**r**) can be cast in that form by setting $f(r)=\rho (r)\chi C(r)$, where $\chi C(r)$ is the characteristic function of the primitive cell $C$ of the lattice.

**r**

_{m}. Therefore, the number of particles per cell $nc=\rho \u0304v$ is given by

*ρ*(

**r**) by Eq. (40), we can single out the term Ω

_{intra}of the grand potential, which describes the interactions between particles located in the same cells,

*N*/

*n*

_{c}is the number of lattice sites, and the double integral refers to the interaction between particles in a single cell. Because of Eq. (42), the latter scales as $nc2$, whereas in the absence of the self-interaction, it should scale as

*n*

_{c}(

*n*

_{c}− 1). To restore the proper scaling, one should then replace Ω in Eq. (10) with the functional Ω′ obtained by subtracting to Ω the term Ω

_{self}given by

*n*

_{c}. In the above relation, $f\u0303(k)$ is the Fourier transform of

*f*(

**r**) defined as in Eq. (7).

For the GEM-4 interaction, *ρ*(**r**) is strongly localized at the lattice sites, implying that *f*(**r**) has a sharp peak at *r* = 0. As a consequence, Ω_{self} is hardly distinguishable from the state-independent quantity *N*Φ(0)/2 so that subtracting it off amounts to a trivial shift in the chemical potential, which does not affect the phase diagram, provided that the same shift is applied to the fluid phase. However, in fact, the mapping of the nematic model into the spherically symmetric one rigorously holds irrespective of the form of *f*(**r**), as will be shown in the following.

**r**inside the unit cell can be expressed as

**r**=

**A**·

**s**, where

**s**varies in a cube $Q$ of unit edge. Here, we release that constraint on

**s**and adopt the same prescription for any point

**r**in $R3$. Similarly, we define a vector

**u**such that

**k**=

**B**·

**u**with $k,u\u2208R3$. We then introduce the function

*γ*(

**s**) defined by

*γ*(

**s**) given by

*ρ*(

**r**) given by Eq. (40), the function

*ψ*(

**s**) defined in Eq. (15) is expressed as

_{self}/

*V*to the functional Ω/

*V*of Eq. (16) to obtain the functional Ω′,

_{0}(

*r*) is based on the identification of

**k**

_{λ,n}with the Bravais lattice vector $kn0\u2261B0\u22c5n$; see Eqs. (18) and (19). Similarly, in the integral over

**u**of Eq. (51), one has

**k**

_{λ}=

**k**

_{0}≡

**B**

_{0}·

**u**. Following the argument of Sec. III, we introduce the local density

*ψ*

_{0}(

**s**) ≡

*λψ*(

**s**). Within the present representation of

*ρ*(

**r**) via Eq. (40), this amounts to setting

*γ*

_{0}(

**s**) ≡

*λγ*(

**s**). Multiplication of both sides of Eq. (51) by

*λ*then yields

_{0}(

*r*) and the reciprocal lattice generated by

**B**

_{0}at inverse temperature

*β*and fugacity

*z*

_{0}=

*λz*,

*provided*that its number of particles per cell $nc0$ is the same as that of the aspherical system

*n*

_{c}. This is indeed the case since Eqs. (19) and (42) imply

Therefore, the mapping described by Eqs. (24)–(27) relating the grand potential, density profile, and Helmholtz free energy of the nematic model with the corresponding quantities of the spherically symmetric model still holds rigorously even after the self-interaction term is removed. As a consequence, the infinite degeneracy of the ordered configurations of the GEM-4 nematic model stemming from that mapping also holds.

## VI. APPLICATION TO THE GAUSSIAN NEMATIC MODEL

In Ref. 16, the aspherical potential Φ(**r**) of Eqs. (2) and (3) was considered in the special case *m* = 2, i.e., that of the Gaussian nematic potential. The mapping of the periodic phases of the system into those obtained for the spherically symmetric Gaussian potential and their ensuing infinite degeneracy were discussed in detail, casting some light on the findings reported in former investigations of the same model.^{12,17}

*T*= 0, when the Helmholtz free energy

*F*reduces to the internal energy

*E*and Eq. (27) becomes

*T*= 0,

*E*/

*N*is given exactly by the lattice sum

On the basis of the formalism presented in Secs. III and V, we are now in a position to establish that, within the mean-field DFT treatment employed here, degeneracy will persist even at *T* ≠ 0. In order to do so, a preliminary technical point must be addressed, which originates from the fact that the phase diagram of the Gaussian potential^{27–32} is widely different from that of the GEM-*m* potential with *m* > 2.^{14,15}

Both potentials present a soft-core repulsion, but, according to the classification proposed in Ref. 25, the latter belongs to the so-called class of *Q*^{±} interactions, whose Fourier transform has a negative minimum at a non-vanishing wave vector. Systems of particles characterized by such interactions crystallize at every temperature and form cluster crystals featuring multiple occupancy of the lattice sites. By contrast, the Gaussian potential belongs to the class of *Q*^{+} interactions, whose Fourier transform is always positive. Those systems crystallize only at low temperature, form crystals with single occupancy such that *n*_{c} = 1, and display reentrant melting at *T* ≠ 0.

For the Gaussian potential, the mean-field functional (10) *per se* fails to predict the occurrence of crystal phases (see Appendix B). We verified this by checking that the unconstrained numerical minimization of the functional via the same algorithm mentioned in Sec. IV always converges to the uniform phase, irrespective of the temperature and of the trial density profile fed into the algorithm.

As discussed in detail in Appendix B, a closer examination of functional (10) indicates that there are two shortcomings, which prevent it from displaying a minimum at some non-uniform, periodic density profile *ρ*(**r**). The first one is the inclusion of self-interaction in the grand potential, which can be corrected as described in Sec. V. At variance with the case of GEM-*m* interactions with *m* > 2, that correction is important in order to account for freezing. The second shortcoming consists in the fact that the number of particles per cell *n*_{c} is always determined *a posteriori* from the relation $nc=\rho \u0304v$ after the optimization with respect to the lattice constant has been performed. However, for the Gaussian potential, freezing occurs only for *k*_{B}*T*/*ϵ* ≲ 0.01.^{30} At such low temperatures, *n*_{c} is bound to assume integer values because the thermal motion is unable to redistribute the particles among the cells. In this situation, which is, in fact, that commonly found in atomic crystals, at a given density $\rho \u0304$, it is the lattice constant that is determined by *n*_{c}, rather than the other way round. Actually, this effect shows up also in the phase diagram of the GEM-4 potential pertaining to the *Q*^{±}-class at low temperature, where it causes an infinite sequence of isostructural transitions between crystal phases with increasing values of *n*_{c} as the density is increased,^{33,34} which again functional (10) is unable to predict.

The failure to account for the discretization of *n*_{c} at low temperature is in fact the main limit of the present mean-field treatment of soft-core potentials, as evidenced by the comparison with the simulation results of the low-temperature phase diagram of the GEM-4 presented in Ref. 33. For *Q*^{±} interactions, an approach aimed at overcoming such a limit was proposed in Ref. 35, where it was applied to the one-dimensional case. For systems such as the Gaussian model, where *n*_{c} is locked at *n*_{c} = 1, the simplest recipe to correct this deficiency is representing *ρ*(**r**) as in Eq. (40) and treating Eq. (42) with *n*_{c} = 1 as a normalization condition to be enforced from the outset by suitably choosing *f*(**r**). This rather natural prescription is widely adopted, see, e.g., Ref. 30, where *ρ*(**r**) was described as a superposition of normalized Gaussians. In this way, the functional Ω′ = Ω − Ω_{self} with Ω_{self} given by Eq. (44) is made able to deal with the occurrence of crystal phases in the Gaussian model. Since fixing *n*_{c} to unity implies that the lattice constant is determined by $\rho \u0304$, the calculation is most easily performed in the canonical ensemble by minimizing the corresponding free energy *F*′ at fixed $\rho \u0304$.

Turning now to the Gaussian nematic model, the relevant point is that neither of the above modifications affects the possibility to map this system into the spherically symmetric Gaussian model. As for the subtraction of the self-interaction term, this has already been discussed in Sec. V. As for the constraint *n*_{c} = 1 to be enforced on *f*(**r**), we observe that *n*_{c} is conserved by that mapping, as shown by Eq. (53). Therefore, imposing *n*_{c} = 1 in the nematic model entails the same constraint in the spherically symmetric model and vice versa. We are then brought to the conclusion that the infinite degeneracy of the crystal phases of the Gaussian nematic potential found in Ref. 16 at *T* = 0 will not be disrupted by moving to *T* ≠ 0. As the spherically symmetric Gaussian potential features single-occupancy bcc and fcc phases,^{30} single-occupancy bcc* and fcc* phases will be found in the nematic case.

## VII. CONCLUSIONS

The main topic of this paper consists in a study of the ordered configurations of the GEM-4 nematic model. This model is a straightforward generalization of the spherically symmetric GEM-4 potential that has already been extensively investigated^{14,15} such that the pair interaction is assumed to depend not only on the distance between the particles but also on the angle between the vector **r** connecting them and a fixed direction **u**. The iso-surfaces of the potential are no longer spheres, but ellipsoids of revolution around **u** with aspect ratio *λ*.

The present investigation complements a former one of the same model,^{13} which was also focused on the ordered phases, but concerned mostly the effect of the asphericity of the interaction on the shape of the density profile *ρ*(**r**) in the neighborhood of the lattice sites, irrespective of the lattice. Here, instead, we are more interested in the identification of the crystal lattices and in the resulting phase diagram.

To this end, we have resorted to density-functional theory by employing the same mean-field grand potential functional used in Ref. 13 and have shown that, by rescaling the lattice vectors along the direction of **u**, the periodic phases of the GEM-4 nematic model at temperature *T*, fugacity *z*, and average density $\rho \u0304$ are mapped exactly into those of the spherically symmetric GEM-4 at the same temperature *T*, fugacity *λz*, and average density $\lambda \rho \u0304$. Since the phase diagram of the spherically symmetric GEM-4 model is known,^{14,15} the phase diagram of the nematic model is determined straightforwardly on the basis of that mapping.

Obviously, for the spherical potential, an arbitrary rotation of the lattice leads to the same kind of lattice. However, when one turns to the nematic system, the original and the rotated lattices are mapped into intrinsically different structures, which, nevertheless, correspond rigorously to the same grand potential. A rather surprising consequence of this situation is that the periodic phases of the nematic model are infinitely degenerate. It is obvious that the same scenario will apply to any GEM-*m* nematic interaction with *m* > 2 and, more generally, to any soft-core nematic interaction of the *Q*^{±} class.

Both the mapping of the nematic model into the spherically symmetric model and the ensuing degeneracy had, in fact, been discussed in a former study of the zero-temperature configurations of the Gaussian nematic potential,^{16} whereby the asphericity is introduced in the same way as in the GEM-4, but GEM-4 is replaced by a Gaussian. The phase diagram of the Gaussian model is very different from that of the GEM-4 model, and some care must be exerted in order to describe its crystal phases via a mean-field density functional like that employed here:^{30} specifically, the self-interaction term must be subtracted from the functional, and the single-particle occupancy of the lattice sites must be enforced explicitly. Here, we have shown that, after those prescriptions have been taken into account, the mapping still holds. Therefore, within the present mean-field treatment, the results of Ref. 16 concerning the type of crystal phases of the Gaussian nematic model and their infinite degeneracy can be extended to *T* ≠ 0, as for GEM-4. In both cases, the analysis developed here can be easily generalized to an assembly of parallel particles whose shape is described by a generic ellipsoid other than a spheroid and whose particle interaction is soft.

Nevertheless, one should keep in mind that at *T* = 0 the free energy reduces to the energy that is given exactly by Eq. (55) and fulfills Eq. (54). As a consequence, the mapping and the degeneracy of the crystal phases are exact properties of the nematic model, for both the Gaussian and the GEM-4 interactions, whereas at *T* ≠ 0, the functional considered here is just an approximation, and there is no warrant that those properties would still hold if a different, more accurate functional were adopted.

For the Gaussian nematic model, such a possibility is corroborated by the numerical simulation study by Prestipino and Saija,^{12} where they claim that at *T* ≠ 0 the degeneracy disappears. We are not aware of simulation studies of the GEM-4 nematic model, but for the spherically symmetric GEM-4, there is evidence that the mean-field functional (10) considered here is surprisingly accurate,^{14} provided that the temperature is not so low that the effects of the discretization of the site occupancy *n*_{c} considered in Refs. 33 and 34 become important. On the other hand, at *T* = 0, the mapping holds exactly. Hence, it is likely that it will remain valid at all temperatures, either rigorously or nearly so, implying that even a hypothetical exact treatment would still predict infinitely many degenerate or at least almost degenerate structures.

Having said so, it should be acknowledged that the relevance of such a scenario to real systems is highly questionable. The GEM-4 or, more generally, the GEM-*m* nematic model comes about as an attempt to introduce anisotropy so as to represent soft-core particles of ellipsoidal shape. However, its main assumption that all the particles share the same orientation from the outset is clearly highly artificial and has little resemblance with the situation expected for actual aspherical particles, whereby the orientation of each particle is a degree of freedom of the system. First, we notice that, as soon as *T* ≠ 0, the assumption of equal orientation becomes untenable. However, even at *T* = 0, when real aspherical particles may indeed form a perfect nematic crystal, it is by no means granted that they will behave as the present model.

Let us consider, for instance, the two-dimensional system of dumbbells studied in Refs. 36 and 37. Each dumbbell consists of two spherical particles at contact interacting via a GEM-4 potential, to which a small hard core has been added in order to prevent particle superposition. Unlike the GEM-4 nematic model, the orientation of each dumbbell can be chosen at one’s will, and the interaction between two such dumbbells depends on the orientations of both of them with respect to the direction of the line connecting their centers. At *T* = 0, this system forms a nematic crystal, whereby the dumbbells are arranged into a centered rectangular lattice and are aligned along the direction of the longer side of the rectangular unit cell.^{36} At first, the situation may resemble that of the nematic model considered here: in both cases, compared to the lattice expected for the spherically symmetric GEM-4 (the triangular lattice in two dimensions^{38}), the crystal stretches along the direction of the nematic director. However, in the nematic model, the nematic director is *fixed*, and it will remain unchanged upon rotating the crystal, whereas for actual dumbbells, the director is attached to the lattice, and if the crystal is rotated, it will rotate too. Maintaining that it should remain fixed in space upon rotation would lead to sub-optimal arrangements, thereby releasing the degeneracy. Although the study carried out in Refs. 36 and 37 concerns the two-dimensional case, we expect that similar considerations will apply in three dimensions as well.

A different situation in which the infinite degeneracy does instead take place concerns the crystal phases of hard-core spheroids with additional contact interactions at *T* = 0 investigated by Schiller *et al.*,^{39,40} who in fact were the first to point out the reason for the degeneracy. Actually, that system is more complex than the soft-core models studied in this paper because in their case, the occurrence of degeneracy depends on the degree of symmetry of the crystal lattice of the hard-sphere system into which the original lattice is mapped.^{40} However, if the degeneracy is there, the mechanism leading to it is the same as that described here.

## ACKNOWLEDGMENTS

D.P. acknowledges support from Università degli Studi di Milano, Project No. PSR2022_DIP_008-Linea 2. The authors acknowledge TU Wien Bibliothek for financial support through its Open Access Funding Program.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Davide Pini**: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Software (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Markus Weißenhofer**: Data curation (equal); Software (equal). **Gerhard Kahl**: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.

### APPENDIX A: GRAND POTENTIALS AND ROTATION PARAMETERS OF SOME DEGENERATE LATTICES

As an independent test of the degeneracy of the crystal phases of the nematic GEM-4 potential based on the analytical considerations put forth in Sec. III, here, we report the results of the numerical minimization of functional (10) for anisotropy parameter *λ* = 1.5 and two thermodynamic states, corresponding to *k*_{B}*T*/*ɛ* = 0.5, $\rho \u22c6\sigma 03=2$, $\rho \u0304\sigma 03=2.258$ and *k*_{B}*T*/*ɛ* = 0.5, $\rho \u22c6\sigma 03=6$, $\rho \u0304\sigma 03=8.316$, respectively. The minimization runs were initialized by different trial density profiles, resulting in different lattices.

For both states, the values of −Ω/*V*, i.e., of the pressure *P* obtained for three such lattices, are shown in Table I. Those values coincide up to the eleventh decimal digit, which we regard as convincing evidence that they do refer to degenerate structures.

k_{B}T/ɛ
. | $\rho \u0304\sigma 03$ . | ϕ (deg.) . | θ (deg.) . | ψ (deg.) . | $P\sigma 03/\epsilon $ . |
---|---|---|---|---|---|

0.5 | 2.258 | 45.00 | 59.99 | 120.01 | 25.554 893 523 065 |

0.5 | 2.258 | 134.14 | 88.07 | 84.54 | 25.554 893 523 065 |

0.5 | 2.258 | 32.77 | 3.96 | 165.21 | 25.554 893 523 064 |

0.5 | 8.316 | 161.66 | 66.70 | 14.03 | 274.262 904 362 684 |

0.5 | 8.316 | 44.85 | 123.87 | 62.71 | 274.262 904 362 687 |

0.5 | 8.316 | 155.20 | 138.92 | 70.60 | 274.262 904 362 689 |

k_{B}T/ɛ
. | $\rho \u0304\sigma 03$ . | ϕ (deg.) . | θ (deg.) . | ψ (deg.) . | $P\sigma 03/\epsilon $ . |
---|---|---|---|---|---|

0.5 | 2.258 | 45.00 | 59.99 | 120.01 | 25.554 893 523 065 |

0.5 | 2.258 | 134.14 | 88.07 | 84.54 | 25.554 893 523 065 |

0.5 | 2.258 | 32.77 | 3.96 | 165.21 | 25.554 893 523 064 |

0.5 | 8.316 | 161.66 | 66.70 | 14.03 | 274.262 904 362 684 |

0.5 | 8.316 | 44.85 | 123.87 | 62.71 | 274.262 904 362 687 |

0.5 | 8.316 | 155.20 | 138.92 | 70.60 | 274.262 904 362 689 |

Moreover, the matrix **A** of the primitive lattice vectors always fulfills Eq. (31), meaning that these degenerate structures are indeed obtained by applying to the primitive vector matrix **A**_{0} of the spherically symmetric GEM-4 the rotation-rescaling recipe described in Sec. III. For the state at lower density, the lattice generated by **A**_{0} is the bcc, whereas for the state at higher density, it is the fcc. The matrices $A0bcc$ and $A0fcc$ are obtained by setting *λ* = 1 in Eqs. (32) and (33), respectively.

**R**appearing in Eq. (31) is identified by the Euler angles ϕ, θ, ψ via the expression

For each lattice, the Euler angles are also reported in Table I.

### APPENDIX B: MEAN-FIELD DFT FOR GAUSSIAN PARTICLES

In this appendix, we show that, unlike in the GEM-*n* case with *n* > 2, for the Gaussian potential $\Phi 0(r)=\epsilon exp(\u2212r/\sigma 0)2$ the mean-field functional (10) is unable to predict the occurrence of crystal phases and that those phases are retrieved only if (i) the condition *n*_{c} = 1 for the site occupation number is enforced explicitly *and* (ii) the self-interaction term is subtracted off the functional.

*ρ*(

**r**) as in Eq. (40) where, following Refs. 14, 15, and 30, a Gaussian form for

*f*(

**r**) is adopted, giving

*α*is a variational parameter determining the degree of localization of

*ρ*(

**r**), and

*n*

_{c}is related to the nearest-neighbor distance

*d*and the average density $\rho \u0304$ by

*ζ*is a numerical coefficient characteristic of the lattice.

^{14}for the Helmholtz free energy

*F*

*ρ*(

**r**) ln

*ρ*(

**r**) of Eq. (10), the contributions due to the overlap of Gaussians centered at different lattice sites have been disregarded.

_{0}(

*r*) is also a Gaussian, the integrals in Eq. (B3) can be performed analytically. We introduce the dimensionless quantity

^{30}

*d** =

*d*/

*σ*

_{0},

*β** =

*βɛ*, and the summation is performed over the neighbor shells. The numerical constants

*ℓ*

_{i}and

*c*

_{i}are characteristic of the lattice and represent, respectively, the number of lattice points in the

*i*th neighbor shell and the square of the shell radius in units of the square of the nearest-neighbor distance

*d*.

*γ*also

*d*can be regarded as a variational parameter and that

*n*

_{c}is determined

*a posteriori*from Eq. (B2). By substituting Eq. (B2) into Eq. (B5) and differentiating with respect to

*d**, we find

*G*(

*x*) is given by

*x*→ +

*∞*, one has

*G*(

*x*) → 0

^{−}.

*G*(

*x*) at small

*x*, we use the Poisson identity and rewrite

*g*(

*x*) as follows:

*i*th neighbor shell of the reciprocal lattice and the square of the shell radius in units of the square of the reciprocal lattice nearest-neighbor distance.

*η*is another numerical constant such that the nearest-neighbor distance of the reciprocal lattice is expressed as

*η*/

*x*. For

*x*→ 0, Eq. (B9) gives

Since *G*(*x*) is a monotonically increasing function, one has −3 ≤ *G*(*x*) < 0 in the whole region [0, +*∞*). Equation (B6) then shows that the derivative of *F* with respect to *d** is always positive. As a consequence, the optimal value of *F* is assumed for *d* = 0, which, by virtue of Eqs. (B1) and (B2), implies *n*_{c} = 0 and $\rho (r)\u2261\rho \u0304$, i.e., the homogeneous phase. The free energy of the homogeneous phase would become infinitely negative because of ln *n*_{c} in Eq. (B3); however, this is in fact just a shortcoming due to the aforementioned approximation of disregarding the overlap contribution in the term *ρ*(**r**) ln *ρ*(**r**). Clearly, this is not justified for *d* → 0 because in this limit the distance between the Gaussians centered at different lattice sites in Eq. (B1) becomes vanishingly small.

*F*by the free energy

*F*′ given by

Since the above correction does not depend on *d*, the derivative of *F* with respect to *d* is unaffected by this term.

Therefore, in order for the mean-field free energy (B5) or (B12) to predict a crystal phase, it is necessary that the value of the site occupation number *n*_{c} is fixed *a priori*, thereby implying that, at any $\rho \u0304$, the lattice constant *d* is univocally determined by Eq. (B2).

On the basis of the above discussion, both sides of Eq. (B13) are again always positive: the optimal *F* is then assumed for *γ* = 0, i.e., *α* = 0, which, according to Eq. (B1), corresponds once more to the homogeneous phase $\rho (r)\u2261\rho \u0304$. As before, the term $ln\alpha \sigma 02$ in Eq. (B3) would give an infinitely negative value for *F*, but this would be corrected by taking into account the overlap between different Gaussians in the ideal-gas contribution.

*γ*yields

*T** = 1/

*β** =

*k*

_{B}

*T*/

*ɛ*. The above equation always admits solutions, provided that

*T** is low enough.

In the light of the above discussion, we conclude that, in order to describe the crystal phases of the Gaussian potential, both conditions (i) and (ii) stated in the first paragraph of this Appendix must be enforced in the mean-field free energy functional. Having said so, we remark that, since these phases are recovered as minima of *F*′ under the constraint *n*_{c} = 1, their free energy is still higher than that of the homogeneous phase, which represents the *unconstrained* minimum of both *F* and *F*′. As a consequence, the fluid–crystal transition cannot be predicted by using the mean-field *F*′ for both the fluid and crystal phases. If the mean-field *F*′ is used for the crystal, then a different approach should be adopted for the fluid. This is in fact what was done in Ref. 30, where the fluid was described by the virial route of the hypernetted chain theory.

## REFERENCES

*Fundamentals of Inhomogeneous Fluids*

*Solid State Physics*