Unconventional Continuous Structural Disorder at the Order-Disorder Phase Transition in the Hexagonal Manganites

Sandra H. Skjærvø, Quintin N. Meier, Mikhail Feygenson, Nicola A. Spaldin, Simon J. L. Billinge, Emil S. Bozin, and Sverre M. Selbach NTNU Norwegian University of Science and Technology, Department of Materials Science and Engineering, NO-7491 Trondheim, Norway ETH Zürich, Materials Theory, Wolfgang Pauli Strasse 27, CH-8093 Zürich, Switzerland Forschungszentrum Jülich, JCNS, D-52425 Jülich, Germany Chemical and Engineering Materials Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA Brookhaven National Laboratory, Condensed Matter Physics and Materials Science Department, Upton, New York 11973, USA Columbia University, Department of Applied Physics and Applied Mathematics, New York, New York 10027, USA


I. INTRODUCTION
The multiferroic hexagonal manganites, h-RMnO 3 (R ¼ Dy-Lu, In, Y, or Sc), are improper ferroelectrics where the polarization emerges as a secondary effect due to an improper coupling to the primary distortion mode. This results in unusual ferroelectric domain structures in which topological protection of the domainwall intersections causes fundamentally and technologically interesting physical properties ranging from early universe analogues [1][2][3] to nanoscale conducting channels [4][5][6][7][8][9][10][11]. In spite of multiple studies, the evolution of the polarization with temperature has not been explained on a microscopic level. In particular, studies based on powder neutron [12,13] and x-ray [13][14][15][16][17][18] diffraction show good agreement with a polar model describing the average structure at low and intermediate temperatures, while structural anomalies have been reported between 800 K and the Curie temperature T C . These findings led to a range of reported values for T C and proposals of two distinct structural phase transitions [12][13][14][15][16]19], although it is now understood that there is in fact only one phase transition at T C , with the polarization slowly emerging as a secondary effect [20][21][22][23][24].
The high-symmetry nonpolar average structure of the prototypical hexagonal manganite h-YMnO 3 above T C displays P6 3 =mmc symmetry. This structure features corner-sharing MnO 5 trigonal bipyramids separated by layers of Y 3þ ions, each coordinated by eight oxygens [ Fig. 1(a)]. The unit cell triples across T C at around 1250 K, as a noncentrosymmetric, but zero-polarization, zone-boundary K 3 mode condenses [25]. This K 3 distortion consists of z-axis displacements of the Y atoms and a corresponding tilt of the trigonal bipyramids. It is described by a two-component order parameter (Q, Φ) with amplitude Q and angle Φ [21]. Here, Q is related to the magnitude of bipyramidal tilting and Y displacements, and Φ is related to the direction of the bipyramidal tilting and the displacement pattern of the Y, e.g., updown-down.
The resulting energy landscape resembles a Mexican hat [ Fig. 1(b)] with three lower-symmetry space groups in the brim. The group-subgroup relationships are illustrated in Fig. 1(c). The subgroups are characterized by the angle of the order parameter Φ. A general angle leads to P3c1 symmetry [IM, Fig. 1(e)], but special values of Φ lead to higher symmetries. For Φ ¼ nðπ=3Þ, with n ¼ 0; 1; ...; 5, the system is polar with space group symmetry P6 3 cm The order-parameter observables for the polar subgroup are the bipyramidal tilt amplitude (angles α A and α P ) and corrugation of Y cations (ΔY). Note that ΔY denotes the distance in the c direction between Y1 and Y2. Green arrows indicate the directions of the bipyramidal tilts. The angle between the O1-O2 line and the c axis defines the apical tilt α A and is a robust measure of the order-parameter amplitude Q irrespective of the value of Φ. The plane through the three in-plane oxygens (one O3 and two O4) relative to the ab plane defines the planar tilt angle α P and is related to both the order-parameter amplitude Q and angle Φ. (e)-(g), Structures of the subgroups at different order-parameter angles Φ: Y off-centering pattern (top) and bipyramidal tilting directions, indicated by green arrows (bottom). Atomic sites for the three subgroups are labeled and coded with colored circles to emphasize which positions are symmetry related in each phase. Atom positions with the same color markings are aligned along the c axis and have the same multiplicity. A detailed overview of the atomic positions for the space groups is given in Fig. S1 of the Supplemental Material [26]. system becomes antipolar with space group symmetry P3c1 [AP, Fig. 1(g)].
Landau theory analysis of the K 3 mode following Artyukhin et al. [21], in combination with first-principles calculations of hexagonal YMnO 3 , suggests an insignificantly small energy difference between the polar and antipolar symmetry when only the K 3 mode is included. Stabilization of the polar state occurs only when the K 3 mode couples to a polar Γ − 2 mode that causes a shift of the Y atoms towards the Mn-O layer. Because of this improper coupling, nearly all Φ angles result in a net polarization that reaches its maximum in the P structure (around 6 μC cm −2 at room temperature) but vanishes for the AP structure. The P structure becomes favored by about 100 meV per unit cell in comparison with the AP structure and by more than 600 meV per unit cell in comparison with the NP structure at 0 K [21]. The exact values of the relative energies are highly temperature dependent. The resulting Z 6 symmetry of the configuration space causes the unusual sixfold ferroelectric domain patterns characteristic of h-YMnO 3 [27].
This established model [21] of the Mexican-hat Landau free energy [ Fig. 1(b)] describes the average symmetry evolution of the system reduced to the degrees of freedom given by the order parameter, but it does not address the underlying microscopics. In particular, whether the transition mechanism is closer to the displacive limit [the order parameter (Q, Φ) goes to zero both locally and on average at T C ] or the order-disorder limit (local order parameter is conserved) is not known.
Here, we provide such a local-structure description of the atomic structure of h-YMnO 3 from ambient temperature to 1273 K, across the ferroelectric transition T C . We combine pair distribution function (PDF) analysis of neutron total scattering data with conventional Rietveld refinement to probe structural coherence and to distinguish short-range from average long-range order. This analysis reveals a surprising and unconventional behavior of the local structure as a function of temperature that cannot be explained either by a conventional order-disorder or by a displacive transition picture. In a conventional order-disorder transition, the low-temperature distortions persist above the transition in the local structure but are not evident in the long-range ordered structure due to averaging over variants of the distorted structure [28][29][30]. Our PDF results show that below about 800 K, the average and local structures evolve consistently, with smoothly decreasing distortions related to a decreased order-parameter amplitude. However, between about 800 K and T C , the average and local structures progressively diverge from each other. Upon crossing T C , the local structure below 16-20 Å changes significantly less than the average structure, suggesting an order-disorder transition. The fit of the polar model to the local structure gradually deteriorates, while fits with the antipolar structure gradually improve, consistent with increasing fluctuations of the order-parameter angle Φ on heating. At T C , where long-range cell tripling disappears in the average structure, the fits of the polar and antipolar models become equivalent, though neither model fits ideally. The transition at T C is therefore an order-disorder transition in the sense that the average structure transitions to the undistorted high-symmetry nonpolar structure, while the structural distortions persist in the local structure. In other words, the amplitude Q of the order parameter never goes to zero in the local structure. However, it is unconventional in the sense that the disordering does not occur only between the local polar variants but between all possible angles of the order parameter. These total scattering results are corroborated by a first-principles-based effective Hamiltonian predicting the same behavior.
These two key discoveries-the unconventional nature of the order-disorder transition and the extensive temperature region of symmetry-lowering fluctuations below T C where the material accesses a wide range of structures, intermediate between the polar and antipolar subgroups, at the local scale-reconcile previous inconsistencies in the literature. Such hidden disorder, revealed by total scattering, is also anticipated in other materials where competing interactions give rise to structural frustration, or where energy barriers of different magnitudes give rise to different thermal evolution of order-parameter components.

II. STRUCTURAL DESCRIPTION OF THE ORDER PARAMETER
The order parameter (Q, Φ) is related to three observable atomic displacements [ Fig. 1(d)]: The first, α A , is the bipyramidal tilt angle calculated from the line connecting the apical oxygens (O1, O2) relative to the c axis; the second, α P , is the bipyramidal tilt of the plane through the three planar oxygens (O3, O4) relative to the ab plane; and the third is the out-of-plane off-centering (corrugation) of the yttrium ions, ΔY ¼ cðz Y1 − z Y2 Þ. The latter two are strictly defined only for the polar ground-state symmetry [ Fig. 1(d)]. Shifting the order-parameter angle Φ away from the polar symmetry and towards the antipolar or intermediate space groups breaks the symmetry of the Y2 and O4 sites such that the corresponding order-parameter observables α P and ΔY have to be calculated by including the additional Wyckoff sites (see Fig. S1 in the Supplemental Material [26]).

A. Average structure
We first consider the average structure behavior based on qualitative and quantitative assessment of our neutron powder diffraction data. On the qualitative side, our reciprocal space data show that the (102) and (202) super-reflections of the polar structure disappear above 1223 K, corresponding to the ferroelectric Curie temperature T C (see Fig. S2 in the Supplemental Material [26]), in line with previous reports [12,15,18,24].
Quantitative analysis of the average structure of h-YMnO 3 was carried out by Rietveld refinements of the measured diffraction patterns with the three special space group models NP, P, and AP, for which representative fits are given in Figs. S3-S5 of the Supplemental Material [26]. In Fig. 2(a), we have presented the fit residuals r w for the P and AP models, found in the brim of the Mexican-hat energy landscape, below T C , and for the NP model, found at the top of the Mexican hat, above T C . Refinement parameters from fitting with the P model are presented in Figs. 2(b)-2(d). The models are described and visualized in Figs. 1(e) and 1(g).
The models give adequate fits to the data with weighted profile agreement factors (r w ) below 10%. The polar model gives better agreement than the antipolar model for all temperatures below T C [ Fig. 2(a)]. Despite having fewer variables, the NP model gives a comparable r w to that of the P model above T C . The most surprising observation is a discontinuous jump in r w extracted from the P and AP models at 800 K. The lattice parameters [ Fig. 2(b)] vary smoothly with temperature, also through T C . Above T C , the lattice parameter a of the nonpolar model is multiplied by ffiffi ffi 3 p for direct comparison with the tripled unit cell of the polar model.
There is no discontinuity of the lattice parameters at 800 K that might explain the change in r w there. However, there is a smooth drop-off from linear temperature dependence of the c-axis parameter with an onset at around 800 K. This result is highlighted in the figure by the gray dashed line extrapolated from the low-temperature linear behavior. Such a drop-off has been noted previously [12], as has the extended region of zero thermal expansion of the c axis above T C .
We now turn to the three order-parameter observables α A , α P , and ΔY, whose temperature dependencies are shown in Fig. 2(c). All of these parameters have, by definition, a value of 0 in the high-symmetry nonpolar structure and are thus expected to decrease smoothly to zero on heating to T C , as they are all direct observables of the order-parameter amplitude Q (Fig. 1). The refinements of the α parameters to the polar model smoothly decrease with increasing temperature, becoming close to zero at T C , consistent with a continuous transformation. On the other hand, the corrugation parameter ΔY keeps a large value right up to T C , which would be expected for a discontinuous transition. From about 800-880 K, the fits become unstable and the α parameters start to decrease faster than linearly, in agreement with previous studies [12], possibly explaining the jump in r w in panel (a).
The vanishing polyhedral tilting and persistent ΔY could indicate a scenario in which the structure has off-centered Y ions combined with untilted Mn-O 5 polyhedra, which would give very long out-of-plane Y-O distances. We thus performed density functional calculations to check if this scenario is plausible and found that it is highly unfavorable (see Fig. S6 in the Supplemental Material [26]). We note that the Y and O isotropic atomic displacement parameters (ADPs) U shown in Fig. 2(d) are anomalously large compared to the Mn U, which is nonlinear upon approaching T C . There is also a discontinuous jump in the ADPs on moving into the high-symmetry nonpolar structure, suggesting that broken local symmetry persists above T C [31][32][33][34]. Hence, we conclude that the structural behavior obtained from the reciprocal space refinements is physically unfeasible at the local scale, motivating the use of a local structure-sensitive method.

B. Local structure
In order to investigate the local structure, we perform a PDF analysis on the same neutron scattering data as analyzed in the previous section. Representative PDFs at different temperatures are shown in Fig. 3(a), plotted over a wide range of r (see detailed plots in Figs. S7 and S8 in the Supplemental Material [26].
PDF analysis utilizes both Bragg and diffuse scattering information, revealing time-averaged snapshots of the atomic structure on multiple length scales. In contrast to a conventional crystallographic approach that seeks the highest possible symmetry model consistent with the Bragg data component only, the PDF analysis explores whether such symmetry is broken on a nanometer length scale and, if so, how. The peak positions correspond to interatomic distances, while the widths of the peaks indicate distributions of interatomic distances due to thermal motion and also lower symmetry.
A characteristic of neutron PDFs of materials containing manganese (and other negative neutron scattering length materials) is that peaks corresponding to Mn-non-Mn pairs are negative. This characteristic is most clearly seen for the Mn-O nearest-neighbor bond length at around 2 Å. This peak is a doublet, as is clearly seen in the 298 K data Apart from thermal broadening effects, several other qualitative changes can be observed in the experimental PDFs with increasing temperature, as shown in Fig. 3(a). While we observe qualitative changes in PDF features at all length scales, these seem to become much more pronounced in the high-r region (above 20 Å), suggesting length-scale-dependent structural changes on heating. The apparent order-disorder behavior of the system across T C , inferred from unphysical behavior of structural parameters derived from the Rietveld analysis, can also be further explored by comparing measured and simulated PDFs corresponding to temperatures straddling T C . If the transition were displacive, the PDFs would show clear changes of similar magnitude over the entire r range. We illustrate this in Fig. 3(c) by comparing simulated PDFs using the  [26]. The data temperature is color coded, as indicated by the horizontal color bar. Comparison of PDFs for temperatures bracing the structural transition temperature T C are shown in panels (b) and (c) as follows. Experimental PDFs from our measurements are displayed in panel (b). Simulated PDFs, calculated using the average structure information from Rietveld refined structures reported by Gibbs et al. [12] across the transition, are shown in panel (c). Corresponding difference curves are shown underneath. The shaded green and red areas behind the difference curve indicate, respectively, regions of low and high difference. The presence of length scales exhibiting different behavior can be readily evidenced by comparing the difference curves in panels (b) and (c). At low r, the measured PDFs change appreciably less than what would be expected if the average crystallographic picture were applicable to the local structure, as portrayed by the simulated PDFs. In contrast, the changes across the transition at high r are of similar amplitude in both panels. Note that the structural phase transition temperature in Ref. [12] is 1258 K, 35 K higher than in our measurements. Comparison of simulated PDFs based on structural parameters from the Rietveld refinements of Gibbs et al. to that of this work is provided in Fig. S9 of the Supplemental Material [26]. Rietveld refined structures reported by Gibbs et al. [12] that portray such a displacive picture. In this case, substantial structural changes are evident in the difference curve over the entire r range considered [see Fig. 3

(c)]. However, this
is not what is observed when experimental PDFs measured at temperatures just below and just above T C are compared, as can be seen in Fig. 3(b). It is evident that the observed structural changes are much smaller in the low-r region up to about 15 Å, as compared to those in the r range between 20 and 40 Å. This behavior is characteristic for an orderdisorder mechanism, where the local structure changes insignificantly across the transition, while the global symmetry change is evident in the measures of the average structure. The r values of 16-20 Å separating the local regime from the average regime correspond to approximately three unit cells in plane and 1.5-2 out of plane.
Quantitative information about the local structure can be further extracted by fitting the structural models to the PDF data. When symmetry breaking occurs only in the local structure but not in the average structure, it is customary to fit the PDF data with lower-than-crystallographic symmetry models over the relevant r range [28]. We fitted the PDF data with models corresponding to the four space groups found in the Mexican-hat energy landscape: the P, AP, IM, and high-symmetry NP models over an r range of 12 Å. This length scale corresponds to 1-2 unit cells and is well within the r range implicated by the direct data comparisons, as discussed in relation to Fig. 3(b). Note that the IM model does not represent a single configuration of atoms, but it can describe any angle Φ of the order parameter in the brim of the Mexican hat, including the P and AP models.
In Fig. 4, we show representative PDFs at 298 K, fitted with the polar and antipolar space groups, and at 1273 K, fitted with the polar, antipolar, and nonpolar space groups. Importantly, and in agreement with other evidence of an order-disorder transition presented so far, the lowsymmetry polar model gives much better agreement with the 1273 K data set than the high-symmetry nonpolar model, as is apparent from the appreciably lower χ 2 and clear appearance of misfit regions in the difference curve plotted below. While the antipolar structure fits the measured PDF significantly worse at room temperature [ Fig. 4(b)], it provides a comparable fit to that of the polar model at 1273 K. We will return to this observation later.
Next, we consider the PDF-derived lattice parameters and order parameter observables α A , α P , and ΔY, shown in Fig. 5. The refined lattice parameters are in very good agreement with those obtained from Rietveld refinement, and they show the same behavior. The order-parameter observables are all expected to retain finite values above T C for the now-established order-disorder transition. Indeed, we observe finite values of all three order-parameter observables. Surprisingly, however, they do not decrease by the same amount: α A decreases to about half its low-T value, while α P retains about 80% of its low-T value. Note that ΔY decreases to about 60% of its low-T value. While these observations support the suggested order-disorder behavior, they also indicate that the P model used to fit the PDF data is not able to fully represent the structural changes upon heating.
This case is further substantiated by the uncertainties of ΔY, which become significantly larger than those of α A and α P above 1000 K, suggesting that Y displacements become progressively more difficult to determine within the polar model with increasing temperature. We note that, as with the average structure refinements, the atomic displacement parameters and/or their uncertainties experience an anomalous increase above about 800 K, indicative of the increasing inadequacy of this model for describing the local structure as temperature increases.
Finally, we compare the quality of the fits with the P model to fits with the AP, IM, and high-symmetry NP models. The comparison of the fit residuals r w for the four . 4. Fits of the crystallographic models to the PDF data. The GðrÞ at 298 K fitted over the range r ¼ 1.6-12 Å by the (a) lowsymmetry P and (b) AP models. Corresponding difference curves are shown beneath the PDFs, and the overall fit residuals χ 2 are as indicated. Fits to the 1273 K GðrÞ data by the low-symmetry P and AP models over the same r range are shown in panels (c) and (d), respectively. (e) Fit of the high-symmetry NP model to the 1273 K PDF data: Arrows mark regions of particularly bad fits. High fit residual and eventful difference curves demonstrate that the nonpolar model does not provide an adequate description of the local structure of the high-temperature phase. Fits of the PDF data for additional temperatures are shown in Figs. S10 and S11 of the Supplemental Material [26]. models [ Fig. 5(d)] clearly shows that the AP model gives distinctly worse agreement than the P model at room temperature. As temperature increases, the AP model fits progressively improve, while the polar model fits deteriorate. This observation of stagnating or even increasing fitting error is anomalous, as for systems that do not exhibit structural changes the fit residuals are expected to decrease with heating as thermal fluctuations broaden the data (see, for example, the behavior of bulk FCC Ni in Fig. S14 of the Supplemental Material [26]). At around 1100 K, the two models provide a similar, but not ideal, fit quality, suggesting that neither one of the two provides a complete description of the underlying local structure. Notably, anomalous jumps in U are not fully reconciled by either the P, AP, or IM model of the local structure despite reasonably low fit residuals, indicating that the underlying local structure is of greater complexity than that described by these low-symmetry models (see Supplemental Material, Fig. S13). We observe that the IM model is comparable to the P model over the entire temperature range, with a slightly lower r w , which we attribute to more degrees of freedom (16 for IM, 11 for P). This observation suggests that there is no other minimum at a different order-parameter angle Φ since the degrees of freedom in the intermediate model cover all angles in the brim of the Mexican hat.
We further examine the observation that the antipolar and polar models seem to fit equally well at higher temperatures by fitting the data to constructed models in which we manually rotate the order-parameter angle Φ. The models are constructed starting from the refined structure for the P model and rotating Φ while keeping the order-parameter amplitude Q constant. This method is justified as the Q for refinements with all three subgroup models obtains comparable values (see Fig. S11 of the Supplemental Material [26]). We then fit the models at different temperatures, allowing the polar degrees of freedom as well as the ADPs to refined as presented in Fig. 5(e). As is evident from the figure, the roomtemperature data are much better fitted by the polar than the antipolar model, while as the temperature rises, the fit quality becomes independent of the order-parameter angle and the absolute fitting error increases overall. We note here that the invariance of the fitting error with respect to the order-parameter angle cannot be explained by the order-parameter amplitude Q going to zero, as the constructed models used here have the same amplitude of structural distortions as shown in panel (b), with the amplitude never falling below 50%.
Thus, from our inspection of the fitting errors, we make four key observations on heating: (1) for the P model, r w increases; (2) for the P and AP models, r w becomes similar; (3) for the NP model, r w , while decreasing, remains much larger; and (4) the fitting error becomes independent of the order-parameter angle when approaching T C .
These observations can be explained by a scenario in which the local structure in the high-symmetry phase is neither described by the NP model nor does it take any well-defined average low-symmetry structure. Instead,  [26] for details on the calculation of error bars). (d) Fit residual r w for fitting GðrÞ between 1.6 and 12 Å for the P model found at the minima and the AP model found at the local maxima in the brim of the Mexicanhat potential using fixed lattice parameters found from refinements between 1.6 and 22 Å. The grey vertical line at 1223 K shows the Curie temperature T C , and the grey shaded area between 800 K and 880 K shows the temperature interval where refinements become incoherent. (e) Fitting error r w as a function of the order-parameter angle for a rotation of the angle between the P and the AP structures at constant amplitude. The amplitude for each temperature is shown in panel (b). Comparisons of the refined order-parameter observables and ADPs of the polar model with those of the antipolar and intermediate space group models are shown in Figs. S12 and S13, respectively, of the Supplemental Material [26]. multiple different local configurations in the brim of the hat become represented simultaneously upon heating. Such a scenario is achieved by increasingly large fluctuations of Φ, which increase the sampled phase space and thereby reduce the local symmetry. This kind of continuous disorder is the most likely explanation for the observed anomalies in the fits.

C. Effective Hamiltonian analysis
To better understand the unusual behavior indicated by our local structural analysis, we next analyze an effective Hamiltonian that describes the local structure and energetics of YMnO 3 . We focus specifically on whether our PDF analysis can be rationalized in the context of a competition between displacive and order-disorder behaviors. In the case of a conventional single-component order-parameter system described by a double-well potential, the two energy scales that determine whether displacive or order-disorder behavior dominates are (i) the energy cost per site, V, of removing the distortion that deforms the structure from the maximum to the minimum of the double well (this is the energy difference per unit cell between the high-and low-symmetry structures) and (ii) the intersite energy cost W associated with changes in the distortion from one unit cell to the next (this is the additional energy of forming one locally undistorted unit cell in a background that is homogeneously distorted to the double-well minimum [35]). The displacive limit is reached when V ≪ W since the energy cost for removing the distortion of the units is less than the cost of introducing disorder between locally distorted units. In contrast, the order-disorder limit occurs when W ≪ V, where locally ordered units are strongly favored and the cost of orienting them in different directions relative to each other is small [35][36][37][38]. The situation is more complex for YMnO 3 with its two-component order parameter, and we consider the following effective Hamiltonian, which is a lattice version of the well-established Landau free energy of Refs. [20,21] (details in Appendix C). The resulting effective Hamiltonian is analogous to an XY model with variable spin sizes. We separate the Hamiltonian into the on-site Mexican-hat-potential contribution H V and the intersite interaction energy contribution H W , which describes the energy cost of distorting the tilt pattern: where Here, x i is the coordinate of the ith lattice site, where each lattice site corresponds to a 30-atom unit cell, which is the smallest unit necessary to contain the trimerization distortions. The two-component order parameter ðQ; ΦÞ is a generalized coordinate representing the local trimerization amplitude Qðx i Þ and angle Φðx i Þ. The polarization Pðx i Þ is improperly coupled to Qðx i Þ and is the main cause of the energy barriers in the brim of the Mexican-hat potential. The parameters in H V are taken directly from Ref. [21], and J is extracted from the parametrized gradient terms in the same work (see Appendix C for details and a table of values). Our goal is to determine whether this Hamiltonian, with the parameters appropriate for YMnO 3 , describes a system in an order-disorder or displacive limit. We begin by analyzing whether the local distortion amplitudes Qðx i Þ remain finite or tend to zero at T C . We derive the following expressions for the cost of undistorting a site, V Q , as well as the interaction energy associated with undistorting a site within a homogeneously trimerized background, W Q , using the parameters from Appendix C. Here, n is the number of nearest neighbors, N is the number of lattice sites, and Q 0 ¼ 0.96 Å is the Q value that minimizes the local energy (it finds the bottom of the Mexican hat). We find that V Q =W Q ≈ 1.7. As explained above, the limit of a fully displacive system is given by V Q =W Q ≪ 1, with the small V Q energy cost favoring the local units undistorting when approaching T C , whereas in the opposite orderdisorder limit, V Q =W Q ≫ 1, the local amplitude of the order parameter is preserved on approaching T C . Our value of V Q =W Q ≈ 1.7 therefore points towards an orderdisorder-like behavior of the order-parameter amplitude.
Given this finding that the local order-parameter amplitude remains finite at T C , we next explore the displacive versus order-disorder behavior of the order-parameter angle, corresponding to changes of the angle around the brim of the Mexican hat. In this case, we compare the on-site energy barrier between two local minima within the brim of the hat, V Φ , and the intersite energy cost W Φ arising from a change of the trimerization angle by 2π=6 (between two local minima). We obtain the values giving V Φ =W Φ ≈ 0.4. The small energy barriers between neighboring states in the brim of the hat allow the system to readily access a continuum of angles, leading to a predominantly displacive-like behavior of the orderparameter angle.
In summary, our analysis of the effective lattice Hamiltonian suggests a mixed scenario in which the order-parameter amplitude shows order-disorder behavior, and the order-parameter angle shows displacive behavior at T C . This case is consistent with the model that we proposed above to explain our PDF analysis, in that it predicts that the local order-parameter amplitude should not vanish [Figs. 5(b) and 5(d)] and that the P and AP models, which differ in their local order-parameter angle, should fit equally well at T C [ Fig. 5(d)] but still not fully represent the mix of IM phases present.
Finally, we performed a series of Monte Carlo (MC) simulations on the effective lattice Hamiltonian including only nearest-neighbor interactions. These simulations were performed using the Metropolis algorithm [39] on a 10 × 10 × 10 unit cell with periodic boundary conditions, averaging six different runs over Oð10 6 Þ sweeps. To improve the convergence, we discretized the phase space in 36 equally spaced states in ½0; 2π for Φ and 11 different states in the range ½0; 1.25Q 0 for the amplitude Qðx i Þ. We then analyzed the behavior of the expectation value of the order parameter hQi and the local amplitude of the order parameter ffiffiffiffiffiffiffiffiffiffi hQ 2 i p (see Appendix C for details). We found a T C of about 2200 K, which overestimates the experimental T C by about 60%. This finding is not unreasonable considering our use of the long-wavelength limit for calculating J and our neglect of beyond-nearestneighbor and nonlocal phononic effects. The latter, in particular, will cause the depth of the local Mexican-hat potential to decrease with increasing temperature and will therefore decrease the transition temperature.
The observed behavior, presented in Fig. 6, is consistent with our analysis of the effective Hamiltonian above, as well as with our interpretation of the PDF measurements. When the system crosses T C , the local order-parameter amplitude ffiffiffiffiffiffiffiffiffiffi hQ 2 i p is reduced by only around 20%, while its macroscopic expectation value hQi goes to zero [ Fig. 6(a)]. Considering the distribution of angles, at low temperatures the system is localized at one angle, while the distribution becomes broader around the expectation value upon heating, and disorders between all angles at T > T C [ Fig. 6(b)] [40]. We present a statistical analysis in Fig. 6(c), which shows that far away from T C only the minimum energy states are occupied. The number of occupations of sites with intermediate order-parameter angles increases progressively with increasing temperature. The system ends up having an almost continuous ffiffiffiffiffiffiffiffiffiffi hQ 2 i p and its macroscopic expectation value hQi, as a function of reduced temperature. Note that the local order parameter does not go to zero above T C , whereas its macroscopic expectation value vanishes above T C due to the disordering of the order-parameter angle. (b) Histograms of the order parameter angle at a temperature well below T C (left), just below T C (center), and above T C (right). We see that the angle is well defined at low temperatures, broadens upon heating, and can have any value above T C . (c) Probability P that the states in the minima (Φ ¼ ½ð2πnÞ=6) (brown circles) and the maxima (Φ ¼ f½ð2n þ 1Þπ=6g) (green circles) within the brim of the hat are occupied. At low T, only the minima are occupied [PfΦ ¼ ½ð2πnÞ=6g ¼ 100%], but at higher temperatures, the fractions converge towards the continuous limit marked by the line denoted "C," which corresponds to an equal distribution between all angles. The probabilities were extracted by creating 36 bins in ½0; 2π, and as a result, C has a 16.7% fraction. distribution of order parameter angles over a substantial region below T C . Again, this result is consistent with the fits of the PDF data, with the AP and P phases having similar errors in a large region below T C [Fig. 5(c)].

IV. DISCUSSION
We have studied the structural evolution with temperature of hexagonal YMnO 3 at both local and average length scales. The average structure from reciprocal-space Rietveld refinements is captured by the conventional P model below T C and changes to the NP structure across the phase transition, in agreement with previous studies. However, on the local scale, both the symmetry and the structural evolution are very different from their average behavior, with the two components of the order parameter (Q, Φ) even showing distinct behavior. Combining our experimental (total scattering þ PDF) and theoretical (DFT þ MC) results detailed above, we present the following conceptual order-disorder model for the structural phase transition (Fig. 7).
At low temperatures, the local structure corresponds to the P model of one of the six minima, consistent with the average structure [ Fig. 7(a)] [41]. Upon heating above 800 K, fluctuations of the angle Φ allow the local structure to access configurations of lower symmetry corresponding to the IM model, as illustrated Fig. 7(b). This region of the crossover corresponds to the reduced temperature starting at around ϵ ¼ 0.5 in the MC simulations (Fig. 6), where hQ 2 i and hQi 2 start to deviate from each other and where PfΦ ¼ ð2n þ 1Þπ=6g approaches PfΦ ¼ ð2nπÞ=6g. This case explains the anomalous increase in fit residuals r w for the P model above 800 K and the increasingly similar r w for the different order-parameter angles [ Fig. 5(e)]. This inadequacy of a single structural model is the reason for the artifacts observed in the Rietveld fits that we described in Sec. III A.
On approaching T C , these local fluctuations of Φ from its mean value become more pronounced. Close to, but below, T C , the measured local structure is a dynamic superposition of structural configurations that fluctuate within a continuum of values of Φ [ Fig. 7(c)] while the average order parameter retains nonzero value. This scenario explains our r w fits (Fig. 5), where we saw a similar error of P and AP already in a region rather far below T C . It is consistent with our results from the effective Hamiltonian, where the distribution of order-parameter angles reveals an effectively flat free-energy landscape at high temperatures [ Fig. 6(c)].
Above T C [ Fig. 7(d)], long-range order is lost and the system disorders between a continuum of angles Φ, corresponding to all possible local IM structures with P3c1 symmetry. Since the IM structures are polar, the structure remains polar above T C on the local scale. This special order-disorder transition is unconventional in the sense that it does not disorder between the six degenerate ground states of the low-temperature ferroelectric P model but between a continuum of structures with all possible angles Φ of the order parameter. This continuous disorder above T C distinguishes the ferroelectric transition in YMnO 3 from conventional orderdisorder transitions, e.g., BaTiO 3 exhibiting disorder between different, but discrete, [111]-oriented ground states [28,30]. The continuous positional degree of freedom in YMnO 3 is, in some ways, reminiscent of α-AlF 3 [33] and cristobalite [31], in which the structures are disordered with arbitrarily bent bonds between rigid polyhedral units. Here, however, because of the complexity of the structure, the disorder involves many atoms and bonds on a much larger length scale.
The high T C of YMnO 3 could also be important for the observed disorder. Notably, the corresponding P6 3 =mmc-P6 3 cm improper ferroelectric transition in 2H-BaMnO 3 [42] at 130 K is reported to be predominantly displacive [43]. However, the crystal structure of 2H-BaMnO 3 is fundamentally different in terms of polyhedral units and connectivity from that of YMnO 3 , making a direct comparison difficult. We hope that our work will motivate a detailed PDF study of the phase transition in 2H-BaMnO 3 .
It is well established that there is only one phase transition in the hexagonal manganites at T C with the polarization slowly emerging as a secondary effect [20][21][22][23][24], despite the many reported anomalies observed between 800 K and T C . Until now, these anomalies have not been fully explained. The model we have presented here, describing how the structural order parameter behaves in the local structure, provides an explanation for these anomalies. Entering the regime where fluctuations of the order parameter become increasingly important around 800 K, the temperature dependence of the macroscopic order parameter is renormalized when approaching T C [35,44]. This result will affect related macroscopic properties such as thermal expansion and specific heat [15], and can explain the thermal evolution of the lattice parameters: The reported change of slope 300-400 K below T C and the flattening of the c parameter just above T C [12,[15][16][17]45] can be expected from the coupling of the order parameter to strain. The temperature range of these anomalies also coincides with the reported Ginzburg temperature for YMnO 3 [13], which is a measure of how far below T C one must be to describe the system by mean-field theory. The accessible phase space for local order-parameter fluctuations far below T C becomes limited to a small region around the local minimum. This freeze-in of fluctuations is in compliance with the reported glassy behavior of the system [46].
The gradual loss of long-range structural coherence upon approaching T C , as illustrated in Fig. 7, is not easily detectable by conventional diffraction probing the average structure. Observations of α P and α A decreasing continuously to zero at T C while ΔY goes abruptly to zero are demonstrated both by our own Rietveld refinements and others [12,47]. While discontinuity in one orderparameter observable across T C has been suggested as a signature of a weakly first-order transition [47], the continuous molar volume (lattice parameters) and very subtle changes in the PDFs across T C do not support this interpretation.
Finally, we note that order-disorder and displacive transitions show different features in their dynamical properties. For a purely displacive transition, the frequency of the soft phonon mode goes to zero at T C . On the other hand, the soft-mode frequency always remains nonzero for a purely order-disorder transition. Instead, one observes the emergence of an additional anharmonic relaxational mode at ω ¼ 0, called the central mode [48][49][50][51]. While Gupta et al. [52] did not observe softening of any phonon branches below T C , Bouyanfif et al. [19] and Bansal et al. [53] observed a partial softening upon heating towards T C . However, it is not possible to identify a central mode from their data. These results suggest a mixed displacive and order-disorder character of the phase transition, and our model presented in Fig. 7 reconciles and explains these unusual spectroscopic observations.
We anticipate that the Mexican energy landscape will lead to unusual dynamics for the central mode. The influence of the central mode on properties like dielectric susceptibility have been studied for proper order-disorder ferroelectrics [54,55] but may be very different in the case of this improper ferroelectric with unconventional disorder. We hope our work motivates future studies in this direction.

V. CONCLUSIONS
We discovered unconventional behavior across the improper ferroelectric transition in the hexagonal manganites where the amplitude of the two-dimensional order parameter is order-disorder-like, but the order-parameter angle behaves displacively. This unusual behavior is indicated by our finding that the local structure across T C differs strongly from the established average structure symmetry and that the local low-symmetry distortions are conserved in the high-symmetry phase, both evidence of an order-disorder mechanism. Second, we observe that the low-symmetry distortions at high temperature are not restricted to specific minima but span a continuum of structures found in the brim of the Mexican-hat-shaped energy landscape, effectively lowering the local symmetry from P6 3 cm to P3c1. Our findings reconcile decades of ambiguous reports on hexagonal manganites by showing that the reported anomalies in macroscopic data are caused by the onset of a continuum of local structures.
The only requirement for the presented behavior is the condensation of a disordered multicomponent order parameter with a small barrier between the different orderparameter directions. While such a continuous nature of the order parameter on the local scale is commonly found in magnetic structures, it is highly unusual in crystalline solids where the directionality of bonds and coupling to lattice strain normally favor discrete minima in the potential energy landscape. To the best of our knowledge, this work is the first report of such fundamentally different behavior of two components of a multicomponent order parameterone component showing order-disorder behavior and the other displacive.
Although reports of continuous structural disorder in crystalline solids are rare, the underlying structural and energetic requirements are not. Competing interactions and two characteristic energy barriers of different magnitude could lead to a frustrated structure where all bonds cannot be optimized simultaneously and where fluctuations of one order-parameter component can easily be masked by thermal vibrations. In addition to the improper ferroelectric hexagonal manganite YMnO 3 studied here, similar structural and energetic features can also be found in magnetoelectric multiferroics [56], piezoelectrics at morphotropic phase boundaries [57], and relaxor ferroelectrics [58].
Hidden disorder is thus anticipated in a wide range of material systems with multicomponent order parameters and characteristic energy barriers of different magnitudes. Specific examples of systems where the continuous structural disorder introduced here could be relevant are the lowtemperature tetragonal phase of superconducting cuprates [59], which have been shown to exhibit order-disorder characteristics [60], the pyrochlore Cd 2 Re 2 O 7 where Goldstone-like modes and order-disorder features have been reported [61][62][63], the shallow energy landscape of the layered ferroelectric perovskite PbSr 2 Ti 2 O 7 [64], and (111)-strained perovskite oxides like LaAlO 3 and SrMnO 3 , where the curvature of the Mexican-hat-shaped potential energy surface can be tuned by strain through the choice of substrate [65,66].
This work illustrates that the atomic-scale mechanisms and underlying physics of phase transitions in complex crystalline solids with multicomponent order parameters and competing interactions can easily be overlooked if the local structure is not adequately probed, e.g., by total scattering methods. Neutron total scattering was performed at the Nanoscale-Ordered Materials Diffractometer (NOMAD) [67] at the Spallation Neutron Source at Oak Ridge National Laboratory. A powder sample with mass of about 1.5 g was sealed in a 6-mm-diameter vanadium container and measured in an ILL-type vacuum furnace. The NOMAD detectors were calibrated using scattering from diamond standard powder, and Si standard powder was used to obtain the instrument parameter file for Rietveld refinements. The data were collected at room temperature and between 473 and 1273 K in steps of 20 K for 60 min at each temperature. Measurements of two subsequent temperature cycles between 1113 K and 1373 K on the same sample were performed to investigate chemical expansion as a function of oxygen loss and were found to be insignificant (see Fig. S15 of the Supplemental Material [26]). The structure factor SðQÞ was obtained by normalizing the scattering intensity to the scattering from a solid vanadium rod, and the background was subtracted using an identical, empty vanadium can. Pair distribution functions (PDF) were obtained by Fourier transform of SðQÞ with Q min ¼ 0

ACKNOWLEDGMENTS
The average crystal structure over the whole temperature range was determined by Rietveld refinements with the space groups P6 3 =mmc and P6 3 cm [68], using TOPAS Academic v.5 [69]. PDFs were fitted to the same space groups at different ranges of r using PDFgui [70]. Isotropic atomic displacement factors (U) were used [26], and more details on calculation of error bars for ΔY, α P , and α A can be found in the Supplemental Methods [26]. We note that error bars for the real space and reciprocal space refinements cannot be quantitatively compared. Details of the PDF method can be found elsewhere [32]. Simulated PDFs were calculated using the PDFgui software with lattice parameters, atomic coordinates, and atomic displacement parameters from the Rietveld refinements, both with our data and with the data from Gibbs et al. [12]; they were scaled such that direct visual comparison could be made. The parameter δ 1 was set to the value obtained for the relevant space group in the PDF refinements. For the assessment of fitting errors, we used the following quantities as implemented in PDFgui [70]: where N is the number of data points, p is the degrees of freedom in the fit, G calc and G obs are the calculated and measured PDFs, and

APPENDIX B: COMPUTATIONAL DETAILS
For our density functional calculations, we used the LDA þ U approximation [71,72] as implemented in the abinit PAW plane-wave code [73][74][75][76], with a cutoff energy of 30 Hartree and a k-point mesh of 6 × 6 × 2. A collinear A-type magnetic ordering with a U of 8 eV on the Mn d orbitals was applied. While this magnetic ordering underestimates the band gap, it does not break the space group symmetry, which is important for this study. Here, 30-atom unit cells with lattice parameters from our PDF refinements were relaxed with respect to atomic positions such that the forces were converged within 1E-6 Ha/Bohr (see Supplemental Methods [26]).

APPENDIX C: EFFECTIVE HAMILTONIAN
The conventional expression for the Landau free energy is based on a continuous field model and gives the evolution of the average of the order parameter anywhere in the system. From first principles, we are only able to access properties at 0 K.
In order to access information about the finitetemperature properties and order-parameter fluctuations, we map the Landau free energy into a discretized form.
We take the 30-atom unit cell of YMnO 3 as the discretized unit since a local order parameter cannot be defined for a smaller unit. The terms in the Mexican-hat-potential contribution H V are taken directly from the parameters calculated by Artyukhin et al. in Ref. [21] and reproduced in Table I. To extract the interaction energy terms H W , we mapped the gradient terms onto the XY model Hamiltonian, in order to enforce the necessary 2π periodicity. Note that this expression reproduces the Landau free energy in the continuum limit. We extract the coupling parameter J by dividing each gradient term (s Q in the notation of Ref. [21]) by the square of the relevant lattice constant J i ¼ s i =a 2 i (to approximate the energy cost of changing the trimerization between two unit cells) and then averaging the in-plane and out-of-plane values using J ¼ ðJ xy J z Þ 1=3 [3].
We neglect any gradient contribution from the polar mode.
In our Monte Carlo simulations, we assume that the polarization is minimized on each lattice site, giving ∂H ∂Pðx i Þ ¼ 0; leading to which we insert in Eq. (1). To calculate the expectation values of the order parameter, we need to project it onto an arbitrary axis, for example, From this calculation, the macroscopic expectation value for the order parameter can be calculated using and the expectation value of the local order-parameter amplitude is given by