Global formation of topological defects in the multiferroic hexagonal manganites

The spontaneous transformations associated with symmetry-breaking phase transitions generate domain structures and defects that may be topological in nature. The formation of these defects can be described according to the Kibble-Zurek mechanism, which provides a generic relation that applies from cosmological to interatomic lengthscales. Its verification is challenging, however, in particular at the cosmological scale where experiments are impractical. While it has been demonstrated for selected condensed-matter systems, major questions remain regarding e.g. its degree of universality. Here we develop a global Kibble-Zurek picture from the condensed-matter level. We show theoretically that a transition between two fluctuation regimes (Ginzburg and mean-field) can lead to an intermediate region with reversed scaling, and we verify experimentally this behavior for the structural transition in the series of multiferroic hexagonal manganites. Trends across the series allow us to identify additional intrinsic features of the defect formation beyond the original Kibble-Zurek paradigm.

Topological defects are ubiquitous in nature, emerging in various forms in a large variety of physical systems from atomic to cosmic length scales. In the context of cosmology, Kibble first inspected the link between the possible topology of the corresponding defects and gauge symmetry breaking [1]. Subsequently, Zurek derived a scaling law relating the density of defects and the speed at which the transition point is crossed [2]. Their combined theory is known as the Kibble-Zurek mechanism. Under the appropriate conditions, this mechanism is expected to describe the formation of topological defects in a system that is driven through a continuous phase transition at a finite cooling rate. Since Kibble-Zurek scaling is determined by the critical behavior and should be the same for all systems in the same universality class, Zurek proposed the study of condensed-matter analogues to cosmic systems for its verification.
A variety of condensed matter systems have been investigated to date in an effort to verify the Kibble-Zurek mechanism. Early attempts were carried out on liquid crystals [3,4], superfluid 4 He and 3 He [5][6][7][8] and superconducting rings [9,10]. More recent studies have been conducted on multiferroics [11][12][13], Bose-Einstein condensates [14,15], ionic crystals [16,17], Landau-Zener setups [18], and colloidal monolayers [19]; for a review see Ref. [20]. The case of multiferroics is particularly interesting, as they have provided the first experimental setting clearly compatible with a Kibble-Zurek scal-ing beyond mean-field [11]. On the other hand, for the same system a drastic reversal of this scaling (termed "anti-Kibble-Zurek scaling" [11]) has been reported for fast quenches, although its origin is not understood and its existence has been questioned [13].
In this work we combine first-principles calculations and the theory of critical phenomena to provide a global picture of the Kibble-Zurek mechanism in which, by increasing the cooling rate, defect formation evolves from the fluctuation-dominated Ginzburg region to the meanfield regime. This picture naturally encompasses features of anti-Kibble-Zurek behavior, which can emerge from the crossover between these two distinct regimes. Our model system for this investigation is the series of hexagonal multiferroic manganites, RMnO 3 , here with R = Y, Dy, Er and Tm.
In our scanning probe measurements, both the Kibble-Zurek scaling and the anti-Kibble-Zurek behavior are demonstrated unequivocally as a general feature in hexagonal manganites. In addition, trends that we uncover by studying the RMnO 3 series as a whole reveal additional quantitative features suggesting that the topological defect formation is affected by supplementary ingredients beyond the original Kibble-Zurek theory.
We discuss emergence of additional time-and lengthscales as likely candidates for these extra features, which can appear naturally from the propagation of the phasetransition front, the vortex-growth process, or directly from the eventual discrete nature of the corresponding symmetry breaking.
We prepared single crystals of YMnO 3 , DyMnO 3 , and ErMnO 3 using the floating-zone (FZ) technique as described in Methods (see also Supplementary Information). To complete our analysis, we also consider data for TmMnO 3 reported in [13]. These RMnO 3 compounds undergo a high-temperature lattice-distortive unit-celltrimerizing transition at T Y c 1259 K , T Dy c 1223 K, T Er c 1429 K [12], and T Tm c 1514 K [12,21]. By driving the systems through this structural transition, topological defects are created in the peculiar form of discrete vortices as sketched in Fig. 1 [11][12][13][21][22][23][24]. These vortices correspond to particular solutions of the Landau free energy [24] (1) where Q = (Q cos Φ, Q sin Φ) is the primary order parameter associated with the condensation of a zoneboundary K 3 phonon. This condensation induces the spontaneous polarization P ∼ Q 3 cos 3Φ (Φ = nπ/3, with n = 0, 1, . . . , 5). The polarization alternation of the resulting six trimerization-polarization domain states around the vortices (see Fig. 1c) enables their real space imaging by piezoresponse force microscopy (PFM). The parameters of Eq. 1 for the series of hexagonal manganites calculated in this work using density functional theory (see Methods) are given in Table I. Our samples were annealed above the transition temperature and then quenched at cooling rates ranging from 10 −2 K/min to ≈ 10 5 K/min in a temperature interval around T c of at least ±100 K. To exclude surface-related effects, samples were then thinned by about 100 µm before polishing them and resolving the ferroelectric domain structure by PFM at room temperature (see Methods). The resulting images are shown for DyMnO 3 and Here a0 denotes the zero-temperature value of the parameter a (in the simplest case a = −a0Tcε, where ε = (T − Tc)/Tc is the reduced temperature). The parameter s in (1) is the mean value of the anisotropic stiffness s = (s 2 xy sz) 1/3 . The zero-temperature correlation length ξ0, renormalized zero-temperature correlation lengthξ0, and Ginzburg number Gi are derived from these values. At the ordering temperature Tc, sets of three bipyramids tilt towards or (in this case) away from a common center. This trimerizes the unit cell and induces a spontaneous polarization ±P . The order parameter of the trimerization-polarization is (Q cos Φ, Q sin Φ) with Q and Φ as sketched. b, Topological defects are lines (points on the sample surface) around which Φ changes monotonically in a clockwise or counterclockwise fashion. Close to Tc this change is gradual due to the "dangerously irrelevant" character of the Z6 anisotropy. At lower temperature the Z6 anisotropy becomes fully relevant and six domain states with discrete values Φ = n · 60 • , n = 0, 1, . . . , 5, emerge. c, Possible arrangement of the six domain states around the vortex-like topological defect. For each domain state the bipyramidal tilt pattern is indicated with arrows representing Q and Φ. Fig. 2; the behavior of YMnO 3 is similar. Both compounds exhibit the characteristic domain pattern in which meeting points of the six trimerizationpolarization domains identify the location of the topological vortex defects. DyMnO 3 shows an increase of vortex density n up to a cooling rate of 1 K/min, followed by a striking decrease of n by two orders of magnitude upon further increase of the cooling rate. ErMnO 3 displays qualitatively the same behaviour, but the decrease of n sets in at higher cooling rates than in DyMnO 3 . Note that previous experiments [12,13]  at slower cooling rates, and so did not discover the turn around in the slope of n.
The vortex density as a function of cooling rate for all four RMnO 3 compounds is displayed in Fig. 3, clearly revealing that both the Kibble-Zurek and the anti-Kibble- Data for TmMnO3 were taken from Ref. [13]. Starlike red symbols indicate data points taken on flux-grown samples in Ref. [12]. Insets show PFM images of the domain structure after quenching at 1 K/min for 5 × 5 µm 2 sections. The TmMnO3 image is sketched as its domain size would exceed the shown section.
Zurek behaviour are generic features of the vortex formation in the hexagonal magnanites. We also find an intriguing dependence on the crystal chemistry, with larger R 3+ radius correlating with higher n at a given cooling rate within the Kibble-Zurek region, as well as with a decrease of the cooling rate at which the turn around occurs. In addition, we observe deviations from Kibble-Zurek scaling in the ultraslow-cooling regime. We return to these deviations later, after discussing first the origin of the turn around between Kibble-Zurek and anti-Kibble-Zurek behavior.
According to the Kibble-Zurek mechanism, the vortices are expected to emerge from critical fluctuations with a density that is essentially determined by the rate of cooling through the phase transition and the critical slowing down of the system [25]. When the relaxation of the order parameter becomes slower than the changes introduced by the quenching the thermal fluctuations freeze out and give rise to a non-equilibrated order-parameter distribution.
We first revise this picture by making a distinction between two fluctuation regimes, namely the Gaussian approximation about mean-field and the Ginzburg (or scaling) regimes [26][27][28]. Roughly speaking, fluctuations are assumed to be non-interacting fluctuations in the Gaus-sian approximation while their interaction becomes crucial and controls the critical properties in the Ginzburg regime. These fluctuations determine the vortex density n expected according to the Kibble-Zurek hypothesis as Here ξ(t * ) is the coherence length at the freeze-out time t * , which in turn is the time at which the quenching process becomes faster than the relaxation of the order parameter, τ . For a linear quench with T (t) = (1 + t/τ q )T c , where τ q is the characteristic time set by the cooling rate r = T c τ −1 q , the freeze-out time is given by t * ∼ (τ 0 τ zν q ) 1/(1+zν) [2,20]. Here z and ν are the dynamical critical exponent and the critical exponent for the correlation length according to τ (ε) = τ 0 /|ε| zν and ξ(ε) = ξ 0 /|ε| ν , where ε = (T − T c )/T c , is the (timedependent) reduced temperature [29][30][31].
We see in Eq. (2) that the precise form of the Kibble-Zurek scaling is fundamentally related to the nature of the critical fluctuations at the freeze-out. In the meanfield regime the critical exponents are ν = 1/2 and z = 2, which leads to a mean-field Kibble-Zurek exponent 2ν/(1 + zν) = 1/2. In addition, the microscopic correlation length is ξ 0 = s/a 0 in terms of the Landau free-energy parameters in Eq. (1) with the expansion a = −a 0 T c ε. On the other hand, the critical exponents for Eq. (1) in the Ginzburg regime take the values ν = 0.672 and z 2 of the 3D XY model [13,21,32,33]. As a result, the Kibble-Zurek exponent becomes 2ν/(1 + zν) = 0.58. The onset of interaction between fluctuations characterizing the transition to the Ginzburg regime modifies not only the critical exponents, but also "microscopic" parameters such as ξ 0 , which becomesξ wheres = (s 2 xy s z ) 1/3 represents an averaged gradient coefficient (see Supplemental Material).
The crossover between the mean-field and Ginzburg regimes is defined by the Ginzburg-Levanyuk criterion as the point at which the order-parameter fluctuations in a correlation volume of size ξ 3 reach the magnitude of the order parameter itself [26][27][28]. This crossover can be estimated from observables such as the specific heat. Thus, with Eq. (1) one obtains the so-called the Ginzburg- This has to be compared with the reduced freeze-out temperature ε(t * ) = (τ 0 /τ q ) 1/(1+zν) to determine the formation of vortices in the Kibble-Zurek picture. At slow cooling rates, where ε(t * ) < Gi, the Kibble-Zurek mechanism probes the fluctuation-dominated Ginzburg region. However, as the cooling rate increases, the freeze-out eventually occurs so far from T c that ε(t * ) > Gi and hence the Kibble-Zurek physics emerges from Gaussian fluctuations in the mean-field regime. The critical exponents and the microscopic parameters then change accordingly.
The expected crossover for a system described by Eq.
(1) is illustrated in Fig. 4. Since the Kibble-Zurek exponent 2ν/(1 + νz) is larger in the Ginzburg regime than in the mean-field regime an offset between the related lines in the double-logarithmic plot is expected. The density of defects is then expected to show a dropdown (red line in Fig. 4) when the transition from the Ginzburg to the mean-field behaviour occurs. The Gi values that we obtain from our density functional theory (DFT) calculations, given in Table I, are compatible with our observed crossovers, further justifying our use of the Ginzburg regime value 0.58 for the Kibble-Zurek exponent in the scaling regime at slow cooling. We also see in Table I, however, that the hexagonal manganites have a consistently larger correlation length ξ 0 in the Ginzburg than in the mean-field regime, ξ 0 . This should tend to reduce the offset shown in Fig. 4 and could even lead to a crossing of the Ginzburg and mean-field graphs (Fig. 4 inset), along with a smooth transfer between the regimes as indicated by the orange line. A central factor discriminating between the two scenarios, however, is the renormalization of the relaxation time τ 0 . Equation (3) shows that the correlation length increases in the Ginzburg regime, and with this we expect the relaxation time τ 0 to go up. Hence, the dropdown scenario is the more likely one, providing an explanation for our experimental results in Figs. 2 and 3. We expect that experiments at higher quench rates, which we have not yet been able to access, should show a second turn around, as vortex densities begin once again to increase with quench rate following a mean-field scaling.
It is worth noting that, in the seemingly unrelated case of a quantum phase transition driven by a noisy control parameter, the possibility of a dropdown has also been pointed out [34]. This possibility is in fact analogous in the sense that the departure from the Kibble-Zurek scaling is also related to the running of the critical exponents -from their nominal value to the noise-fluctuation limit in that case.
Next we address the deviations observed in our ErMnO 3 samples in the ultraslow-cooling regime. As shown in Fig. 3, here the vortex density increases but without reaching expected value. Such behaviour was already discussed for other systems [4,14,35] and attributed to vortex-antivortex annihilation. This explanation is consistent with our finding that the fall-off is largest at very slow cooling, when the time spent close to T c is large. It is also consistent with its absence in the flux-grown samples of Ref. [12]; our samples are grown using the floating zone method which we expect to have distinctly different stoichiometric modifications. We indeed observed that the mobility of vortices below T c can vary substantially among samples depending on the annealing gas atmosphere.
We now analyze the pronounced trends with the chemistry across the RMnO 3 series that are clear in Figure  3. In particular, the vortex density for a given cooling rate in the slow-cooling regime increases by three orders of magnitude from Tm to Dy, i.e. with decreasing R 3+ radius. In the Kibble-Zurek picture this is directly related to the microscopic parameters ξ 0 and τ 0 . Our first-principles calculations (see Table I) indicate that the mean-field value for ξ 0 is essentially the same in the four systems, as is the renormalized valueξ 0 ( 1.3ξ 0 ) in the fluctuation regime. Therefore, to obtain a factor ∼ 10 3 in the vortex density in either the conventional mean-field or fluctuation-dominated Kibble-Zurek picture, τ 0 would need to change by an unphysical factor of ∼ 10 6 . Strictly speaking, the renormalization of ξ 0 should be computed at the freeze-out temperature T * = [1 + (τ 0 /τ q ) 1/(1+zν) ]T c . However, this does not help us to reconcile the differences, since to increase the vortex density by the observed factor 10 3 , the freeze-out temperature would have to reach T * ∼ 100 T c > ∼ 10 5 K. This is again unphysical. We thus see that the vortex formation in the RMnO 3 systems displays quantitative features challenging the interpretation in terms of the original Kibble-Zurek mechanism. Specifically, the trend with chemistry does not fit with the standard scenario, even with the role of critical fluctuations beyond the meanfield description fully taken into account.
Another trend in Fig. 3 is the several-orders-of-magnitude shift of the cooling rate at which the departure from the Kibble-Zurek scaling occurs and which is not predicted by the calculated values for Gi in Table I. It is striking, however, that the huge spread of the vortex density in the Kibble-Zurek range is substantially reduced in the anti-Kibble-Zurek regime (from 1000 to 4), and in fact almost restores the approximate R-independence expected from Table I. At least in part, extrinsic factors may be responsible for these unexpected trends with chemistry. Temperaturedependent processes related to chemical or mechanical impurities could affect the vortex formation differently across the RMnO 3 series, simply because of the substantial increase of T c from Dy to Tm. In addition, the chemical and mechanical quality of the samples may change with R. DyMnO 3 , in particular, is much harder to grow hexagonally than e.g. TmMnO 3 because of the closer proximity of the competing perovskite phase. In fact, we find that batches of the same material grown under slightly different conditions can show, at a specific cooling rate, a vortex-density spread up to factor four with overall trends with cooling rate like in Fig. 4 (see Supplementary Information).
These factors can modify, in particular, the thermal conductivity of the samples. This conductivity is another ingredient limiting the overall thermalization of the system, and is such that the instantaneous temperature can additionally be position dependent across the sample. If this happens, the local transitions will not be simultaneous during the experiments. Instead, there will be a phase-transition front propagating at the speed v f = dx dt | T (x,t)=Tc , which needs to be compared with the characteristic speed v c = ξ(ε)/τ (ε) = (ξ 0 /τ 0 )|ε| z of the order-parameter relaxation [36][37][38][39]. Thus, if the front propagates much slower than the order-parameter perturbations, the formation defects can be heavily suppressed in a composition-dependent fashion.
Finally, we briefly discuss other possible origins of the anti-Kibble-Zurek behavior. The thermal conductivity could be one of them since, according to the above, a maximum in the vortex density can be expected at v f = v c , followed by quick decrease as the cooling rate (and hence v c ) increases.
On the other hand, the reversal from Kibble-Zurek to anti-Kibble-Zurek scaling has also been demonstrated in a recent Bose-Einstein-condensate transition experiment [15]. In that case, it has been ascribed to vortex-vortex repulsive interactions. Such interactions could also play a role in our RMnO 3 systems, which certainly display a rather dense vortex pattern in the crossover region (see Fig. 2). In addition, there are additional time scales that can become relevant in the problem. The most obvious one is the time-scale needed to obtain a fully developed vortex from the initial seed. Indeed, by increasing the cooling rate, this process can be expected to be gradually suppressed -simply because the phonons will not have time to propagate the order-parameter perturbations [40]-which will eventually lead to a decrease in the vortex density.
Thus the overall trend in Fig. 3 may still be largely determined by intrinsic properties, and hence directly related to the corresponding universality class. In this regard, the discrete nature of the symmetry breaking can play a crucial role, especially if, as has been repeatedly underlined (see e.g. [41][42][43]), the defect formation is decided after crossing the transition temperature. In fact, a generic feature of Z 6 models like Eq. (1) is the emergence of a second correlation length below T c associated with the phase of the order parameter supplementing that of its amplitude [21]. This, however, is missed in the usual Kibble-Zurek mechanism formulated for U (1) models. This second correlation length diverges faster than the standard correlation length [21]. Thus, if the critical slowing down of the phase itself becomes a dominant process, the vortex density might become affected in a substantial way.
In summary, we have shown theoretically that the Kibble-Zurek picture inherently includes a transfer from the asymptotic scaling, where the system is well inside the Ginzburg region, to mean-field scaling, where fluctuations represent small perturbations. This transfer is general since all phase transitions, from the inter-atomic to the cosmological scale, are formally expected to undergo such a crossover. We matched our predictions with the experimental behavior in the series of hexagonal manganites, and showed that this crossover can be behind the striking reversal of the topological defect density as a function of the cooling rate (anti-Kibble-Zurek scaling). In addition, using density-functional-theory, we quantified the expected Kibble-Zurek behaviour across the RMnO 3 series. This quantification revealed the presence of sizeable chemical trends that call for a vital upgrade of the original Kibble-Zurek considerations for the systems in this class.

Methods
Sample preparation Samples were grown at ETH (DyMnO3), PSI (YMnO3) and Lawrence Berkeley National Laboratory (ErMnO3) by the floating-zone technique as described in the supplementary information and elsewhere [44][45][46]. Rods were oriented by Laue diffraction and cut into zoriented platelets of a few mm lateral size and a thickness of about 0.7 mm.
For minimizing sample-dependent drifts in measurements of the cooling-rate-dependent ferroelectric domain vortex density, all samples of a batch were pre-annealed under identical conditions by heating above Tc for several hours (YMnO3: 1400 K, DyMnO3: 1270 K, ErMnO3: 1470 K) and cooling through Tc at a rate of 5 K/min. To reveal the intrinsic 3D bulk domain structure and suppress surface-dependent effects, samples were then thinned by at least 100 µm by lapping with Al2O3 powder. This was followed by etch-polishing with a silica slurry, revealing shiny surfaces with a rms-roughness below 1 nm before determining the vortex density by PFM.
We found that this density showed systematic variations by a factor of about two, depending on the location of a sample in the original floating-zone rod. For our experiment we selected a sub-set of samples displaying the same vortex density after the pre-annealing within the statistical error of the vortex count.
The quench experiments were performed like the preannealing experiments but with a dwell time above Tc of few tens of minutes and varying cooling rate through Tc. In order to avoid accumulation of chemical drift occurring during the quench cycles, each data point was gained from a different specimen of the preselected set.
First-principles calculations For our density functional calculations we use the projector-augmented wave method as implemented in the abinit code [47][48][49][50]. We use a plane-wave cutoff of 30 Ry and a 6 × 6 × 2 k-point grid. To take into account correlation effects on the Mn atoms we use the LDA+U method within the fully localized limit method as introduced by Lichtenstein et al. [51,52]. We choose a value of U of 8 eV and J = 0.88 eV. For all our calculations we adopt an A-type magnetic ordering of the Mn ions and freeze the rare-earth f electrons in the pseudopotential cores. To extract the parameters in the Landau free energy, Eqn. 1, we first fully relax the P 63/mmc structure to an accuracy of 10 −6 Ry/Bohr. We find unit cell volumes of 364.73 Å 3 , 364.42 Å 3 , 367.50 Å 3 and 353.00 Å 3 for YMnO3, DyMnO3, ErMnO3 and TmMnO3 respectively. We then calculate the force constants using the finite displacement method and extract the eigenvectors of the force constant matrix. We then gradually freeze in the eigenvector of the unstable K3 mode and fit a sixth order polynomial to extract the a0 and b terms. To calculate the gradient term s we exploit the fact that in q-space s(∇Q) 2 reduces to sq 2 |Qq| 2 and s can then be obtained by fitting a parabola to the corresponding branch of the force constant dispersion as shown in [24]. We note that the values of the parameters are dependent on the choice of exchange-correlation functional because a scales quadratically and b to the fourth power with the lattice constant. The trends across the series, however, are robust to the computational details.