Phase diagram for ortho-para-hydrogen monolayers

The phase diagram for orientational ordering of hydrogen monolayers on graphite and boron nitride is revised in view of current theory and experimental observations from nuclear magnetic resonance (NMR) studies recently reported for ortho-H2 concentrations 0.35 ⩽ c ⩽ 0.92 and temperatures 0.14 ⩽ T ⩽ 1.80 K. The characteristic interaction coupling Γ0 = 0.50 ± 0.03 K and the crystalline field amplitude V0 = 0.70 ± 0.10 K are derived from experimental data, and distinct types of the local orientationally ordered structures are analysed using a proposed model for site-diluted uniaxial quadrupoles on a triangular plane lattice of hexagonal symmetry. The long-range periodic pinwheel structure and the short-range quadrupolar glass (QG) phase are stable above the 2D site-percolation limit, cp = 0.72, and for 0.48 < c < cp, respectively, where quadrupolar-order effects dominate. At very low T, the QG phase shows instability with respect to local dipole-like polarization effects and the ground state changes to a hindered rotor state. Two para-rotational local ordered PRA and PRB structures driven, respectively, by positive and negative crystalline fields are well distinguished, experimentally and theoretically, and are rather different from the unique PR phase suggested earlier by Harris and Berlinsky.


Introduction
The orientational ordering of hydrogen molecules in 3D and 2D is a particularly interesting problem because the ordering is determined by interactions that can be calculated from first principles, and experimentally one can measure the local order parameters directly. The interactions are dominated by short-range, highly anisotropic electrostatic quadrupolequadrupole (EQQ) couplings that are the same as those in the gas state [1]. Because of the incompatibility of the symmetry properties of the EQQ interactions with the geometry of the crystalline lattice, not all nearest-neighbour pairs of molecules can be arranged in the EQQ minimum energy configuration of an isolated pair, which is a 'Tee' structure with mutually perpendicular molecular axes. This frustration is reminiscent of the random-bond exchange of the spin glasses. Studies of the diluted hydrogen systems are particularly valuable compared to other orientational systems and spin glasses in general, because one has clean experimental conditions for which the effects of disorder due to site substitution can be accurately controlled and monitored. Experimental studies therefore allow us to test our understanding of highly frustrated disordered systems for which the relevant timescales are easily accessible by experiment. Furthermore, the relevant degrees of freedom (molecular rotations) are pure quantum states, and are therefore introduced, at least at low temperature and low pressure, as rotor states with a single angular momentum J = 1.
Quadrupolar systems are different from dipolar systems, which exhibit time-reversal symmetry [2]. In the absence of time-reversal symmetry breaking interactions, the local quadrupolar order is described by a second-rank tensor Q αβ = J α J β + J β J α T with α and β referring to Cartesian components in a chosen coordinate system where · · · T refers to a thermodynamic average at temperature T . The tensor Q αβ has five independent components: the three principal local axes, and the alignment and the eccentricity with respect to those axes. Because the order parameters involve second-rank tensors, the free energies will, in the absence of special symmetries, contain third-order invariants and the phase transitions from hightemperature free-rotor states to long-range orientationally ordered phases will be first order 3 . The 35.3 ordered states are periodic for high rotor concentrations, namely a P a3 structure [1] in 3D and a pinwheel (PW) or herring bone (HB) structure [3]- [5] in 2D. The experimental studies for bulk H 2 showed that the long-range order is lost well above the site-percolation threshold for the hcp lattice, c (hcp) p = 0.204. This is similar to the case of the nonmetallic spin glass Eu c Sr 1−c S, which has an fcc lattice with c (f cc) p = 0.208 [6]. This sharp dependence on concentration results from the strong influence of spatial disorder on frustration of interactions with increasing dilution above c (3D) p .
Careful studies at low temperatures [7]- [11] showed that the local ordering at low ortho-H 2 concentrations involves not just a random freezing of local axes, but that the degree of alignment (or quadrupolarization) also varied from site to site, which is the signature of frozen short-range correlations. This new low-concentration, low-temperature state was called a quadrupolar glass (QG) [7]. The QGs are in a way similar to pseudo-spin-1 glasses [12]- [14], Potts glasses [2,14] and structural glasses [15]. One expects that QG states would also occur in 2D films but that the critical concentration for long-range periodic order should be even higher in 2D because of the enhanced effects of disorder for a 2D system with a smaller number of spatial constraints. This surmise was recently verified experimentally in both monolayers (and submonolayers) [16] and in thick films (from two to 12 layers) [16]. The 2D system, however, is appreciably different from the 3D case because of the presence of a crystalline field that is much stronger in 2D than in 3D. The 2D crystalline field results principally from the interaction of each molecule with the substrate. As a consequence, there is a rich variety of different long-range orientational states that were predicted by Harris and Berlinsky [4] for the non-diluted EQQ system. We will see that in the site-diluted case the crystal-field effects will overwhelm those of the short-range quadrupolar interactions at sufficiently low concentrations.
New local-order phases, different to those discussed in [4], have been reported at very low concentrations [16], but up to the present it has been unclear as to whether the pure EQQ coupling or non-direct phonon-mediated effects in combination with crystal-field effects are responsible for the observed NMR line shapes. Furthermore, a classification of the underlying locally ordered structures has not been given. The aim of this study is to provide an in-depth analysis of the observed phase diagram deduced experimentally for solid (o-H 2 ) c (p-H 2 ) 1−c monolayer films on boron nitride (BN) [16]. The analysis is based on the data [16] for the NMR lineshapes and elaborated within the framework of the proposed 2D site-disordered model system (J = 1) c (J = 0) 1−c . A comparative examination of the findings by Harris and Berlinsky for the particular case of the pure J = 1-system [4] is also given.
The paper is organized as follows. After an overview of the fundamental concepts we introduce in section 2 the microscopic Hamiltonian for the 2D random model system and deduce the self-consistent equations for the order parameters and the local conjugate fields. The variation of these order parameters with temperature and concentration for the different types of global and local ordered state are then compared in section 3 with those used to describe the experimental phase diagram. This comparison permits the establishment of a qualitative understanding of the full range of the phase diagram for ortho-hydrogen monolayers that is summarized in section 4.

Model parameters
The phase diagram for the 2D pure J = 1 rotor system was carefully studied by Harris and Berlinsky [4] using two model parameters: (i) the EQQ coupling constant 0 and (ii) the 35.4 crystalline field amplitude V 0 for both signs (see the inset in figure 1). This diagram was proposed [4] for the case of pure o-H 2 , but to the best of the authors' knowledge, no theoretical research on the 2D phase diagram for solid ortho-para-hydrogen mixtures has been given. The order-disorder boundary, separating the para-rotational (PR) and antiferrorotational (AFR) phases in the 3D (J = 1) c (J = 0) 1−c system, and the PR-QG boundary were described using improved mean-field (MF) theories (see, respectively, [10] and [14]). These findings are used as the starting point of our analysis of the 2D PR-long-range and PR-short-range boundaries established experimentally in [3] and [16], respectively. This analysis will be given for the (J = 1) c (J = 0) 1−c 2D model system in terms of the three parameters (T , 0 , and V 0 ) extended to the concentration c (≤1).
A key point to the description of orientational ordering in monolayers of o-p-H 2 is based on the observation that the extrapolation of the order-disorder boundaries on the 2D phase diagram to the c = 1 positions derived for the two distinct experimental systems, i.e., o-p-H 2 on graphite [3] and o-p-H 2 on BN [16], is consistent with the PW-PR boundary found [4] by Harris and Berlinsky for the case of o-H 2 . This implies that both sets of experimental data can be understood through a generalization of a simple model [4]. The values of 0 (=6Q 2 /25R 5 0 , where Q is the molecular electrostatic quadrupole moment) calculated for the nearest-neighbour separations R 0 for the commensurate 2D √ 3 × √ 3 structures on graphite and BN are 0.534 K [4] and 0.470 K [16], respectively (to be compared with the 3D estimate [1,16] 0 = 0.827 K). In general, the magnitude of the crystal-field amplitude is available from NMR studies [16], but additional and independent experiments are needed to establish the sign of the crystalline field amplitude. Thus, one deduces from the observed NMR line shapes |V (exp) 0 | ≈ 0.6 K for H 2 on BN [5], and |V (exp) 0 | = 0.6-0.8 K for H 2 on grafoil [16]. By comparative analysis of the 2D phase diagram predictions [4], it will be seen that the crystal-field amplitude is essentially positive, and we therefore adopt for our estimates V 0 = 0.70 ± 0.10 K and 0 = 0.50 ± 0.03 K.

Order parameter description
The microscopic description of the site-disordered ortho-para-hydrogen system with pure EQQ intermolecular interactions has been explored extensively within the context of the 3D QG problem [8,13,17]. The thermodynamic rotational states of a given ortho-molecule located at a site i are characterized by the aforementioned local tensorial order parameters Q αiβi . This representation of the ordered states can be transformed unitarily to a convenient set of local variables {L i , σ i , η i } related to the principal local order parameters [7,8,12,18]. Here L i is the three-component local quantization axis, σ i = 1 − 3Ĵ 2 zi /2 T is the local alignment and η i = Ĵ 2 xi −Ĵ 2 yi T is the eccentricity, mentioned in the introduction. Following the simplification commonly used in the analyses of the NMR data ('powder approximation' [8]), we employ a reduced two-component order-parameter description, with {L zi , σ i } and η i ≡ 0, based on the truncated 2D Hamiltonian of N quantum rotors with z neighbours, namely, Here c i is a random occupation number whose mean c is given by the configurational average,  [3]; open triangles for commensurate hydrogen monolayers on BN, [16]. The solid symbols refer to NMR studies of [16]: solid circles, transitions to the QG state; diamonds, transitions to the HR state. The inverted triangles refer to the vanishing of the small splitting of the NMR lines in the PR state. Inset: theoretical phase diagram from figure 2 in [4]. a and b are the tricritical points identified by Harris and Berlinsky [4], and c is the minimum in the observed transition temperatures from the PR to the PW state.
where the Legendre function is determined by the angle R 0 of the vector R 0 connecting centres of mass of the nearest rotors. The local-axis tensor ij = 5P 20 (L zi )P 20 (L zj ) (= 00 00 (L i , L j ), see equations (6), (7) in [12] for the general case). All vectors are introduced in a crystal-lattice spherical coordinate system with the z-axis normal to the substrate. In the 2D case, all the intermolecular axes are located in the plane and, thus, ξ 0 = 9/8. The exchange  [3]; open triangles [16].
energy J ij and the crystalline field h i in equation (1) can then be written as 4 with P 20 (L zi ) = (3 cos 2 i − 1)/2 where i is the polar angle of the orientation of the principal molecular axis L i . It is important to note that the signs of the competitive exchange J ij and the field h i are determined by the local-order axes and can lead to either a FR (ferrorotational) or an AFR local structure.
The model free energy F N of the system of N rotors, with both fixed spatial configuration . The thermodynamic average, · · · (MF ) T , is evaluated using the MF Hamiltonian H (MF ) N introduced by the trial local molecular field ε i , namely

35.7
A straightforward calculation for a given configuration {c i } yields where σ i is the local quadrupolarization (or alignment). The local equilibrium conditions (∂F N /∂σ i = 0, ∂F N /∂ε i = 0) provide the local molecular field ε i (T ) = z j =1 J ij σ j c j and the local order parameter σ i (T ) = −∂F (MF ) N /∂ε i . This yields a self-consistent equation for As seen from equation (4), the effective one-site MF description emerges naturally through the macroscopic free energy f (c, T ) obtained in the thermodynamic limit (f (c, T ) = lim N →∞ F N ({c i }, T )/N) that gives rise to the concentration behaviour introduced by a configurational average over all possible configurations Note that our description of the site-disordered (c ≤ 1) system is consistent with the ordered case (c = 1 and z = 6) considered by Harris and Berlinsky. Indeed, one can see that macroscopic equation (17) in [4] follows from the microscopic equation (5) Harris and Berlinsky [4] dealt with the ideal case of a perfect monolayer of o-H 2 (c i = c = 1) adsorbed as a commensurate triangular √ 3 × √ 3 lattice on grafoil. They calculated the free energies using three local order parameters σ i , η i and µ i , and two Euler angles (ϕ i , χ i ) for axially distributed principal axes. In what follows we discuss qualitatively the distinct structures of the 2D phase diagram [4] within a reduced order-parameter set { i , σ i } introduced by equations (1) and (2), which is relevant in description of the NMR data for diluted quadrupolar systems.
Let us review distinct orientational states by Harris and Berlinsky in terms of the set { i , σ i }. In [4] the PR high-T phase was suggested to be a single sublattice formed by rotors with small σ i (T ) (=σ 0 (T )) and oriented perpendicular to the monolayer plane with i = j = 0 in the high-T limit. This phase is a precursor of the PW phase (see the inset in figure 1). As follows from equation (2), both phases are induced by a constant positive crystal field h i = 2V 0 /3 in the presence of a negative exchange interaction with J ij = −3 0 /2. The PR phase is characterized by a positive local orientational order parameter σ 0 (T ) (and η 0 (T ) ≡ 0) with a ground state σ 0 (0) = 1, for sufficiently large crystal fields V 0 > 22.5 0 .

35.8
As shown, the distinct phases of the concentrated 2D system are well distinguished by the sign of the local order parameter, which in turn is determined by the sign of the crystalline field.

Observed phase diagram
Among the structures predicted theoretically (PR, FR, PW, HB-in and HB-out) [4] only a few have been identified experimentally because not all of the conditions required can be realized for a given substrate. As shown in figure 1, the long-range ordered phase inferred from NMR studies at high ortho-H 2 concentrations on graphite [3] is in good accord with similar studies [16] of monolayers on BN for concentrations c > 0.72. Note that the lowest concentration for which the long-range order is deduced (c (exp) p ≈ 0.70) is close to the 2D site-percolation threshold c p = 0.72 established for the honeycomb lattice [6]. Below c p (see figure 1), new short-range-ordered states appear beyond the classification given in [4].
The boundary between the long-range-ordered and the disordered phases, observed [3] for concentrations 0.72 ≤ c ≤ 0.92, presumably the PW-PR boundary, is characterized by experimentally well defined order-disorder transition temperatures T The quantitative analysis given below shows that the difference between the observed PW-PR and the predicted [4] PW-PR boundaries is due to deviations from the short-range FR-like structure in the real PR phase. These deviations are induced by strong local crystalline fields, which eventually lead to either positive or negative macroscopic order parameters σ i (T ) C . The crossover between these two short-range-order states is identified by the characteristic temperature T (1), make a firm basis for numerical analysis of the long-range order observed in the site-disordered system.

Long-range-ordered structures
In general, the MF predictions by Harris and Berlinsky are expected to be valid for diluted systems, at least for small disorder. Thus we expect that the observed long-range boundaries can be described with the help of the reduced coordinates V (c, T )/ (c) and T / (c) that are related to the principal phases predicted [4] by MF for the particular case c = 1. The PW-PR boundary of the 2D J = 1 diagram proposed in [4] is wedged between two multicritical points, (a) (T a / 0 = 3.50, V 0 / 0 = 9.82) and (b) (T b / 0 = 5.61, V 0 / 0 = 2.34), that are shown as points a and b in the inset of figure 1.
First, the question arises as to whether the observed extrapolated critical temperature T (exp) c (1) = 0.84±0.03 K of the order-disorder phase transition is consistent with the EQQ model given in equations (1) and (2) with the parameters V 0 = 0.70 ± 0.10 K and 0 = 0.50 ± 0.03 K.
To answer this question we cast the PW-PR boundary [4], between the points a and b, into an explicit numerical form, namely If one adopts the experimental data on T (exp) c (1) and V 0 , equation (6) results in (exp) 0 = 0.165 ± 0.010 K, i.e., the deduced (exp) 0 = 0 /3. To clarify this discrepancy with the estimate 0 = 0.50 ± 0.03 K obtained for the direct EQQ short-range interaction, we analyse other conceivable EQQ-type interactions that qualitatively would reproduce the PW-PR boundary given in equation (6). We therefore consider the situation where the direct EQQ interaction (that decays with distance as R −α with α = 5) is not dominant in dilute o-p-H 2 systems, but an indirect EQQ interaction, mediated by phonons, such as static (with α = 3) or dynamic (α = 1) intermolecular EQQ-type interactions, is predominant (see discussion of equation (6) in [20]). In this case the numerical fit of equation (6)  c n with n = α/d where d is the space dimension. As seen from the analysis given in figure 2, n (c) deduced from experiments [3,16] is consistent with d = 2 and α = 5. We therefore conclude that the pure direct EQQ interaction governs the site-diluted PW-phase for positive crystalline fields. We infer that the PW-PR boundary predicted in [4] is qualitatively consistent with that observed in the high-concentration case, but the quantitative discrepancy between 0 and where is introduced near point C for the real phase diagram. The proposed form of V (c, T ) accounts for the crystal-field fluctuations and will be justified theoretically below. The solution of equation (7) is analysed in figure 3 for the case of c = 0.70 available from the NMR frequency splitting observed [16] in the PWB phase and is justified by the fitting parameter γ (exp) 1 = −0.8. In figure 3 the ideal case of the PW structure (that assumes [4] i = 0, j = π/2 and therefore γ 1 = −1/2) is also shown. A deviation of the observed nearest-neighbour orientational-order correlation factor γ (exp) 1 from that in the PW phase is associated with the disorder effect of the local crystalline fields of negative sign.

Short-range order effects
For the orientationally disordered, but short-range correlated PRA state by the local order parameter σ A (c, T ), we consider a configurational average of equation (5) using a Gaussian distribution for the local-field fluctuations ε i ≡ ε i −ε 1 (see [14] where the relevant means ε 1A = ε i ≈ σ 2 A . The first moments of the EQQ local field can be introduced on the basis of the 3D results (see equation (4) in [17]), namely (10b) The mean ε 1A of the random local field immediately follows from ε i (T ) given in equation (4). The MF correction δε 1A introduced in equation (10) is due to the EQQ fluctuations and is given in a modified form for the 2D case. The MF variance ε 2 is restricted to the leading Zeeman term, i.e., ε 2 No solutions of equations (9) and (10) with a positive correlation factor γ 1A > 0 are found that would be compatible with the observed quadrupolarization. Again, this indicates that the observed PRA phase is not compatible with the idealized FR structure (γ 1 = 1) suggested [4] for the PR phase. The concentration dependence of the local order parameter σ A (c, T ) is given (see analysis in figure 4 for two different temperatures) by a solution of equations (9) and (10)  ≈ 3/4) is similar to that of the canonical PW phase but with a modification that in-plane rotors show out-of-plane orientations. This picture is consistent with the stability of the HB-out structure, which occurs in positive crystalline fields (see the inset in figure 1). With decrease of temperature, during crossing of the disorder-order PRA-PW boundary of real o-p-H 2 , the short-range disorder of ortho-molecules lying in the plane disappears, that leads to the decrease in the magnitude of γ 2 from 0.7-0.8 to 0.5.
At low temperatures (T < V 0 ) and low concentrations (c < c p ), we expect to find a locally ordered PRB structure that results from short-range EQQ correlations in the presence of large negative random crystalline fields (|V | V 0 ) anticipated in equation (8) below the T (AB) x (c) line. Following the MF prediction [4], this structure is expected to be modulated by competitive local-order structures characteristic of the symmetry of the HB-in (or '2-in' with γ 1 = 1/4) and the HB-out ('2-out', γ 1 ≈ 1/16) phases. As seen from the inset in figure 1, both phases minimize the EQQ interaction and stabilize the long-range order in the concentrated system. In the site-disordered case, the 'mixture' of HB-in and HB-out structures causes high frustration for the EQQ interaction at small c. Indeed, for such a competing 'mixed' HB structure, exemplified  (9) and (10) for the two temperatures. The fitting parameters γ 1 and γ 2 are also shown.
by a rotor 'in plane' ( i = π/2) and its nearest-neighbour 'out of plane' ( j = π/4), the correlation factor is γ 1 = −1/8. The basic condition |γ 1 | γ 2 , required for the formation of the QG phase, seems to be satisfied within this region and consistent with experiment. With further decrease of temperature the QG state becomes unstable and the so-called [16] hindered rotor (HR) phase emerges.
This picture is qualitatively consistent with our analysis of the observed PRB local quadrupolarization σ B calculated from the equation This result follows immediately from the configurational average in equation (5)  ), justifies the analytical form given in equation (8). Analysis of the available experimental data on the basis of equation (11) is shown in figure 5 for the case of c = 0.44. We see that at low temperatures the PRB phase (with γ 1B = 0 and γ 2B = 1) is characterized by a local structure similar to that for the 2D isotropic QG phase. Meanwhile, the first-order fluctuation effects, introduced in equation (11) by the term δε 1B , are insufficient to fit the data for T > 0.6 K and a more complete analysis that takes into account higher-order correlation effects (fluctuation fields) is needed. These effects can be introduced in terms of the fields ε 2η and ε 2σ which are conjugate to the eccentricity and the alignment QG order parameters [14], respectively, η 2 i C = q η and σ 2 i C = q σ . But their description is beyond the scope of the truncated Hamiltonian given in equation (1).  [16]; open triangles, onset of PW state in the anomalous upturn region of the phase transition boundary [16]. The curves designate solutions of equation (13) for an adjustable parameter λ = 2.3. Other parameters shown include the fitting parameter γ 2 . Inset: left, predictions for concentration dependence of the effective crystalline field at distinct temperatures derived from equation (8) using experimental data [16] for T (exp) x ; right, the variance of the crystalline field in bulk ortho-para-H 2 [14]. defined by the observed crossover temperature T (exp) x (c). These fundamental model parameters extend the corresponding characteristics V 0 and 0 of the J = 1 rotor model [4]. The 'bare' interaction parameters 0 = 0.50 ± 0.03 K and V 0 = 0.70 ± 0.10 K are estimated on the basis of the experimental data available for hydrogen monolayers on graphite [3,5] and BN [16]. These scaling relations are justified by analysis of the 2D phase diagram described through the observed temperature-concentration behaviour of the long-range and the short-range order parameters σ (c, T ).

35.13
The orientational monolayer structures are considered using a simplified (truncated) 2D Hamiltonian introduced through the reduced version of the real 3D Hamiltonian elaborated for the model (J = 1) c (J = 0) 1−c 2D system. On one hand, this simplification ignores the eccentricity of the local order (by putting η i ≡ 0), but on the other hand, it permits one to estimate the observed (macroscopic) quadrupolarization σ (c, T ). Analysis of the local structures is developed in terms of a reduced local order parameter set { i , σ i }, where i is the polar angle and σ i is a local alignment with respect to the principal local order axes, and the nearest-neighbour correlation factors γ 1 and γ 2 (defined in equation (10b)).
The observed order-disorder boundary T (exp) c (c) above c p = 0.72 has been compared with the PW-PR instability line described in [4]. We have inferred that instead of the postulated [4] PR phase with all rotors normal to the substrate, the corresponding PRA phase is due to locally 35.16 correlated arrangements formed by nearest neighbours (with γ (exp) 1A = −1/3), which are shortrange precursors of the out-of-plane sublattice (PWA structure with = 0, π/2).
The arrangements of the disordered phases below c p are driven by fluctuations of both the crystalline and the EQQ local fields. In the PRB phase (with σ B < 0), most of the orientationally disordered molecules lie in the plane, but at high temperatures their configuration cannot be established unambiguously within the limited set of order parameters. We have demonstrated that the observed PRA-PRB crossover line is due to the interplay between the EQQ coupling order and the crystalline-field disorder effects, which are both temperature and concentration dependent. It is worth noting that the crystal-field effects, that are much stronger that those in the 3D case, are suppressed by dilution near the site percolation threshold c p . As a result, the short-range QG phase and the long-range PW phases common for, respectively, the disordered and the ordered pure EQQ 3D systems are observed near c p (see figure 1). Note, that according to the predictions given in [4], strong negative crystalline fields tend to align the rotors in the plane but do not prevent establishment of the long-range order in the HB-in phase. This is consistent with the fact that the 2D QG phase appears below the PRB phase. We therefore infer that, unlike the 3D case, the local principal axes for the rotors in the 2D QG state lie mainly in the plane.
The upturn of the phase boundary between the long-range PW structure and the PRA states near the site-percolation threshold is most unusual. This region corresponds to values of concentration where competing positive and negative crystal-field effects are most pronounced. Qualitatively, the MF equations support the stability of the PW phase in the presence of negative fields (see inset in figure 1) and the anomalous upturn is therefore associated with this narrow region of stability of the PWB phase ( figure 3).
Unlike the PRA phase, which is the high-temperature precursor of the PWA phase, the PRB phase on cooling leads to the highly correlated QG state and then to the HR regime, which both have no analogues in the c = 1 diagram [4]. In addition, the stabilization of the PRB phase for high negative crystalline fields (see figure 1 and inset in figure 6) makes it plausible to assume that the local PRB structure is locally similar to the long-range predicted [4] HB-in ('2-in') structure. Furthermore, at moderate negative fields a competition between HB-in and HB-out local structures gives rise to the origin of the 2D QG phase with very small average quadrupolarization |σ G | 1 (see figure 12 in [16]), that also assumes [17] q G σ 2 G . With decrease of temperature, the 2D QG order shows an instability and transforms into a new local order with unspecified σ H R < 0, characteristic of the HR phase driven by strong negative crystalline fields (figure 5). A complete description of the HR phase should include the strong dipole polarization effects [16]. The stabilization of dipolar order was anticipated by the MF prediction of the 'canted ferromagnetic phase' [21].