Thermal and field-induced transitions in ferroquadrupolar Kondo systems

Recent experiments have examined the impact of a magnetic field on ferroquadrupolar orders in the intermetallic Kondo material PrTi$_2$Al$_{20}$. Motivated by this, we use extensive Monte Carlo simulations to study a diamond lattice XY model of non-Kramers pseudospin-$1/2$ Pr$^{3+}$ moments which crucially incorporates three-spin interactions. This model supports a thermal $Z_3$ Potts ordering transition upon cooling from the paramagnetic phase into the ferroquadrupolar phase. An applied magnetic field along the [110] direction leads to a thermal Ising transition out of the quadrupolar ordered phase. A magnetic field along the [001] direction leads to only thermal crossovers, but supports a spinodal transition out of metastable domains which could be strongly pinned by coupling to elastic lattice deformations. We propose noise measurements as a potential probe to"hear"the spinodal transition. Our work highlights the importance of multispin interactions in Kondo materials near the small-to-large Fermi surface transition.

One basic question which arises when one confronts the plethora of broken symmetry states in Kondo materials is whether one needs to go beyond the simple two-spin RKKY model in modelling the effective interaction between local moments. Indeed, as the Kondo coupling in heavy fermion materials increases, we expect multispin interactions can arise from higher-order perturbation theory beyond the simple RKKY limit. One setting in which such multispin interactions have been investigated extensively is in the vicinity of Mott transitions in quasi-two-dimensional organic materials [37][38][39][40][41] where it has been shown to potentially stabilize exotic quantum * arunp@physics.utoronto.ca spin liquids. From this viewpoint, we expect multispin interactions to also emerge naturally in Kondo materials if we view the the large-to-small FS transition as an "orbital selective Mott transition" of the f -electrons [42,43]. The impact of such couplings has only recently been investigated in multipolar Kondo systems [44][45][46][47], although there has been some suggestive previous work in dipolar Kondo materials [48,49]. Given this, we ask the following questions. Are there any heavy fermion multipolar systems where multispin couplings play a role? Can such interactions lead to observable signatures?
We address these questions in the context of recent experiments on the Pr(TM) 2 X 20 family of materials, where the Pr 3+ ions feature a non-Kramers ground state doublet, which acts as a pseudospin-1/2 degree of freedom on the diamond lattice [22,23]. As discussed in the literature, two components of this pseudospin carry a quadrupolar moment while the third component describes an octupolar moment [23,50]. In this paper, we focus on PrTi 2 Al 20 , which has been proposed to host a ferroquadrupolar (FQ) ordered ground state [19,20,22,28,33] below T Q ≈ 2K, well before the system enters a low temperature superconducting state with transition temperature T SC ≈ 0.2K.
Recent experiments have studied the non-trivial impact of a magnetic field on this ferroquadrupolar ordered state, discovering a strong dependence of the response on the field direction [33,51,52]. For a magnetic field along the [111] direction, the quadrupolar transition appears nearly unaffected, while there appear to be distinct fieldinduced transitions for fields along [001] and [110] directions. We argue here that an appropriate low-energy microscopic model for this material must necessarily include three-spin interactions, and that it reveals itself via the impact of a magnetic field.
Our key results are the following. We show that the model pseudospin-1/2 Hamiltonian for local Pr 3+ moments must include crucial symmetry-allowed three-spin couplings. We use classical Monte Carlo (MC) simulations to study the ordered states, thermal fluctuations, and the impact of a magnetic field in this model. We uncover thermal and field-induced phase transitions and crossovers which are qualitatively consistent with experimental observations. However, our model does not display a subset of field-induced transitions which have been inferred from certain experiments [51]. We argue here that such experiments could potentially probe spinodal transitions out of metastable ground states; such metastable states do exist in the model and may be rendered visible by the strong coupling between the quadrupolar order and elastic lattice deformations. We propose that noise measurements could be used to "hear" such spinodal transitions. While our work here focusses on PrTi 2 Al 20 , our main results are broadly applicable to ferro-quadrupolar orders in diverse materials.

II. MODEL
We consider a simple low energy diamond lattice model for the pseudospin-1/2 non-Kramers doublets, ignoring higher crystal field levels which are split off by a relatively big energy scale ∼ 50 K, which is much larger than the observed ferroquadrupolar transition temperature. The pseudospin Hamiltonian we propose takes the form where τ ⊥ ≡ (τ x , τ y ) denotes pseudospin-1/2 Pauli matrices and H.c. refers to the Hermitian conjugate. We denote nearest neighbor pairs by ij , while the notation ijk refers to shortest site-triplets on the diamond lattice as illustrated in Fig. 1 An easy-plane interaction with J>0 is appropriate to describe ferroquadrupolar XY order in PrTi 2 Al 20 . The magnetic field B = (B x , B y , B z ) couples to the pseudospin at O(B 2 ). This arises within second-order perturbation theory [53,54] via intermediate states involving higher crystal field multiplets, with α > 0. Most importantly, the term Γ > 0 is the simplest symmetry allowed three-spin interaction which breaks the XY symmetry and leads to a Z 3 clock anisotropy. While such clock terms have been previously discussed within Landau theory [45,53,54], there can be no such single-site clock anisotropy term for pseudospin-1/2 models. Hence, the clock-like anisotropy for pseudospin-1/2 cases must necessarily arise from multi-site couplings at the lattice scale. We note that this multispin interaction allows for the τ z eigenvalue to change in steps of ±3, which cannot arise in any RKKY-type two-spin exchange model.
Our motivation here is to understand the ordered phases and thermal transitions of such quadrupolar spin models. We will thus focus on a mean-field theory and large scale classical MC simulations of this model replacing τ ⊥ by a classical XY vector spin. It would be interesting in the future to examine the impact of quantum spin fluctuations in this model.

III. MEAN-FIELD THEORY
At zero temperature and at mean-field level, we replace uniformly τ + i = e iθ which leads to an energy per spin The magnetic field thus competes with the Γ-term, allowing us to probe the impact of the reduction of symmetry from U (1) to Z 3 . We see that applying a field in the [111] direction will not couple at all to the quadrupolar field. A magnetic field along [110] direction gives αB 2 sin θ, while a field along the [001] direction gives −2αB 2 sin θ. Fig. 2 plots the energy landscape as a function of θ and B for these latter two field directions. As seen from Fig. 2(a), a magnetic field along [001] favors θ = π/2, while the two other zero-field minima become metastable minima at θ = (7π/6 − δ, 11π/6 + δ), where δ ≈ √ 3αB 2 /(54Γ − αB 2 ) for weak fields. These metastable minima lie at an energy ≈ 3αB 2 above the ground state, vanishing at a field B sp = 9Γ/α which marks a field-induced spinodal transition.
For a [110] field, as seen from Fig. 2(b), there are two energy minima which lie at θ = (7π/6 + δ , 11π/6 − δ ), where we find that δ ≈ √ 3αB 2 /(108Γ + αB 2 ) for small B. The third zero-field minimum becomes a metastable minimum at θ = π/2, which lies at an energy ≈ 3αB 2 /2 above the global minima. This will convert the thermal Z 3 clock transition into an Ising transition since the three ground states of the Z 3 clock model have been reduced to just two degenerate ground states. Eventually, the two minima merge at B which marks the end point of the Ising transition, where B = 54Γ/α. The metastable minimum at θ = π/2 persists until B .
We thus expect that for the [001] field direction, the field should immediately round off the Z 3 thermal transition into a crossover by selecting one of the three ground states. For the [110] direction on the other hand, we expect the Z 3 thermal transition to convert into an Ising transition for arbitrarily weak fields, with the Ising transition vanishing above a critical field B .

IV. METASTABLE MINIMA
The reason why the metastable minima in Fig. 2 may be important to explore in this system is the following. Imagine we consider starting from a paramagnetic state at high temperature. When we cool below the ferroquadrupolar transition T Q at zero field, we would end up having different Z 3 domains of a typical size L(T ) at some temperature T < T c . Ferroquadrupolar order in this system induces a lattice distortion, which arises from coupling an elastic distortion field ϕ to the quadrupolar degree of freedom, which can be understood via a coarsegrained Hamiltonian where τ ⊥ (r) is the coarse grained quadrupolar order parameter, and λ denotes the magnetoelastic coupling. The two-component elastic distortion field ϕ may be written in terms of the elastic strain tensor ε, as ϕ x = ε xx − ε yy and ϕ y = (2ε zz −ε xx −ε yy )/ √ 3. The impact of quadrupolar order will thus be to produce a small nonzero lattice distortion ϕ. This elastic deformation along different directions in the different domains will tend to collectively pin the local order. Thus, we see that while an applied [001] field will favor a single domain, we have to thermally excite the system out of the metastable domains in order to get to the true equilibrium state. If thermal fluctuations are not significant at low temperature, then such domains might get stuck until we reach a threshold field corresponding to a mean-field spinodal transition; this effect may reveal itself in certain experiments.

V. MONTE CARLO SIMULATIONS
We have carried out extensive classical MC simulations of the Hamiltonian H 0 from Eq. (1). While standard MC updates sufficed to explore the equilibrium phase diagram via measurements of the specific heat and ferroquadrupolar XY order parameter M FQ = i τ ⊥ i , exploring the metastable transitions required us to choose a special update engineered to probe the free energy as function of the angle θ of M FQ . The update involves a local update conserving the direction of M FQ , and a global update jumping between two angles. Combining multiple such simulations at slowly varying angles (typically ∆θ = 2π/1080) we recover the relative weights between them, and ultimately estimate the free energy landscape. Further details on this procedure are provided in the Appendix. Simulations were typically done with two million thermalization and eight million measurement sweeps for a linear system size of L = 9 (corresponding to 2L 3 = 1458 spins) in mapping out the phase diagram and for L = 6 when studying the metastable regions.

A. Zero field phase diagram
The phase diagram of the model with J 1 = 1 in the absence of any magnetic field is shown in Fig. 3(a). Based on a finite-size scaling analysis of specific heat data (for For Γ = 0, we expect this transition to be in the universality class of the 3D XY model, and the corresponding ferroquadrupolar order parameter m FQ = | M FQ | /N indeed continuously vanishes when we heat above the transition temperature T Q ≈ 1.3J 1 as seen from Fig. 3(b). When we turn on Γ = 0, the clock anisotropy suppresses fluctuations and enhances T Q ; furthermore, the transition becomes first order, as is confirmed by the increasingly sharp and discontinuous drop of m FQ across T Q .
B. Impact of nonzero magnetic field Fig. 4 shows the impact of a magnetic field on the phase diagram for fixed J 1 = 1 and Γ = 0.01J 1 . The upper and lower halves in this diagram correspond to fields along the [001] and [110] directions, respectively. We will discuss in the following section that Γ/J 1 ∼ 10 −3 for experiments on PrTi 2 Al 20 ; however, the numerical simulations are more challenging for such small Γ. We thus choose to work with a larger Γ in the MC simulations. The magnetic field required to induce the relevant transitions or crossovers scales as ∝ √ Γ as indicated by mean field theory. We can thus use our MC results, with suitable scaling, to make useful comparisons with experiment.
As expected, a sufficiently large magnetic field leads to a crossover temperature scale since it favors a single free energy minimum as seen from the free energy plots for (I) and (IV) in the left panel, where the color at the bottom depicts the favored angle θ. This crossover temperature T * , indicated by the dotted line, is detected in our MC simulations as a broad hump in the specific heat which does not scale with system size (based on simulations done for linear system sizes L = 6,9,12).
At low field, the [110] direction leads to an Ising transition, denoted by the solid black line, into a phase where there are two degenerate minima as seen from the free energy plot (III) in the left panel. Different MC runs (initialized with a random state) in this regime lead to the system ending up in one or the other minimum, which is depicted by the colors in phase (III) with corresponding θ values shown in the left panel.
For the [001] field direction, even at low field, a single free energy minimum is selected as seen from the left panel (II). The dashed line indicates the crossover field beyond which the metastable free energy minima in the left panel (II) disappear; this corresponds to the spinodal transition discussed from the perspective of mean-field theory above.
The results from our extensive MC simulations are thus broadly consistent with expectations based on mean field theory, but with thermal fluctuations suppressing the magnetic field scale required to induce the observed phase transitions and crossovers. We next turn to the experimental implications of this phase diagram.

VI. EXPERIMENTAL IMPLICATIONS
Our classical MC simulations show that the zero field ferroquadrupolar transition for Γ = 0 occurs at T Q ≈ 1.3J 1 . As Γ increases, T Q increases and the transition becomes more visibly first order, consistent with the behavior of the 3D Z 3 clock (or equivalently 3-state Potts) model. Since the experiments [19,22] see what appears to be a nearly continuous thermal transition at T Q ≈ 2.2K, we assume Γ J 1 . We thus use the value of T Q at Γ = 0, to roughly estimate J 1 ∼ 1.7 K. Microscopic calcu-lations [53] using the measured crystal field levels [19,22] yield α = (gµ B ) 2 (7/3E 4 − 1/E 5 ), where g = 4/5. For PrTi 2 Al 20 , the relevant excited crystal field levels [22] lie at E 4 ≈ 65 K and E 5 ≈ 107 K. This yields α ≈ 0.008 K/T 2 . Assuming the Ising transition for the [110] field direction [51,52] disappears at B ∼ 3 T, we are led to estimate Γ = αB 2 /54 ≈ 10 −3 K, so that indeed Γ J 1 . The spinodal transition for the [001] field direction is then expected to occur around B sp ∼ 1.2 T.
Recently, transport, magnetization, and 27 Al nuclear magnetic resonance (NMR) experiments [51,52] have been used to further explore the phase diagram of PrTi 2 Al 20 . At high fields, B 4 T, for both [001] and [110] directions, there is a significant enhancement of the magnetization [52] upon cooling below T 3 K. However, for low fields, B 2 T, this strong enhancement is absent. Within our theoretical framework, the bulk magnetization is given by the field derivative of the free energy, m = −∂F/∂ B (note that this is not m FQ ). For both field directions, this is given (up to a sign) by | m| ∝ B τ y . For sufficiently high fields, we expect the system to evolve from τ y ∼ 0 for high temperature, to a nearly polarized value | τ y | ∼ 1 at low temperature, so that there would be a significant increase in | m| below a crossover temperature. By contrast, in the presence of metastable domains which we expect at low fields, τ y would be greatly reduced via averaging over the domains, since τ x will also be nonzero in some domains. This leads to the suppression of the bulk magnetization in low fields, so that the sharp increase upon cooling seen at higher fields will now be absent, in qualitative agreement with the data. Furthermore, NMR measurements of the Knight shift [52] are consistent with the bulk susceptibility from the magnetization measurements at high field, but in disagreement at low fields; this disagreement might also indirectly signal the presence of an inhomogeneous domain structure at low fields. The presence of domains is predicted to lead to NMR line splittings, or to inhomogeneous line broadening if the splitting is weak. This expectation is qualitatively borne out from the experimental data [52], but a detailed theoretical understanding needs further analysis using the microscopic hyperfine couplings. Finally, scattering from such an inhomogeneous domain structure could partially contribute to the experimentally observed resistivity anomalies [51].
In order to estimate the typical linear dimension L D of Z 3 domains, we ask when the system with an average order parameter pointing along an XY angle θ would rather break up into domains of the discrete Z 3 order to save bulk anisotropy energy, governed by Γ, at the expense of a domain wall cost arising from J 1 . Assuming a lat-tice constant a, we thus equate 6Γ(L D /a) 3 ∼ J 1 (L D /a) 2 which, for Γ/J 1 ∼ 10 −3 , leads to L D ∼ 160a. This might be the size of typical domains we expect to get pinned by elastic lattice deformations.
One possible experimental route to further exploring such a spinodal origin of the magnetization and transport anomalies could be noise spectroscopy. For instance, resistivity measurements in nanowires of high temperature cuprate superconductors exhibit a telegraph noise, which has been attributed to fluctuating nematic domains or charge stripe domains [55,56]. Similar field and temperature dependent resistivity noise measurements might be valuable in PrTi 2 Al 20 . Another possible experiment might be to detect the actual sound associated with the avalanche of domain rotations one expects near these metamagnetic transitions.
Our proposal of strong spin-lattice coupling leading to field-induced anomalies is distinct from, but not entirely at odds with, a previously proposed explanation [52], which has considered the impact of additional fielddependent quadrupolar exchange couplings within an effective Landau theory. While the microscopic origin of this effect has been attributed to field-induced changes in the Fermi surface [52], and thereby the RKKY Kondo couplings, such terms may also occur if we incorporate field dependent spin-phonon coupling and integrate out the phonons. The microscopic details of such a mechanism, and its connection with the metastable domain picture discussed here, remains a topic for future study.
In summary, understanding the nature of the fielddependent phase transitions and anomalies in PrTi 2 Al 20 may help deepen our understanding of multipolar orders in heavy fermion materials. Finally, our work suggests that multispin interactions must play a broadly important role in Kondo materials.

Computation of Free Energy
The pair-sampling method allows us to fix two angles of the XY order parameter {θ, θ + ∆θ}. By counting the number of sweeps spent at each angle we can determine the relative weight between them Z(θ + ∆θ) Z(θ) = N (θ + ∆θ) N (θ) .
By setting an initial value for Z(0) we can compute successive Z(θ > 0). From this we can derive the angle resolved free energy F (θ) = − log(Z(θ))/β. Note that this process becomes increasing expensive at low temperatures, requiring small ∆θ and a large number of sweeps to get finite counts N (θ) > 0. The free energy per site F (θ)/N site from the two algorithms is compared in Fig. 5(c) and 5(d), and they show very similar angle dependence, although there is some difference is in the amplitude of the free energy variation for a field applied along the [001] direction. The second algorithm with a smoother distribution of proposed updates is likely to be a better representation of the true free energy curve.
Comparing the computed result to the mean-field free energy, we find that the angle dependence is nearly identical; the Monte Carlo and mean-field curves match closely up to an overall ∼ 2.4 scale factor, which reflects a renormalization of Γ due to thermal fluctuation effects beyond mean field theory.