Unveiling hidden multipolar orders with magnetostriction

Broken symmetries in solids involving higher order multipolar degrees of freedom are historically referred to as “hidden orders” due to the formidable task of detecting them with conventional probes. In this work, we theoretically propose that magnetostriction provides a powerful and novel tool to directly detect higher-order multipolar symmetry breaking—such as the elusive octupolar order—by examining scaling behaviour of length change with respect to an applied magnetic field h. Employing a symmetry-based Landau theory, we focus on the family of Pr-based cage compounds with strongly correlated f-electrons, Pr(Ti,V,Ir)2(Al,Zn)20, whose low energy degrees of freedom are purely higher-order multipoles: quadrupoles \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\cal{O}}_{20,22}$$\end{document}O20,22 and octupole \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\cal{T}}_{xyz}$$\end{document}Txyz. We demonstrate that a magnetic field along the [111] direction induces a distinct linear-in-h length change below the octupolar ordering temperature. The resulting “magnetostriction coefficient” is directly proportional to the octupolar order parameter, thus providing clear access to such subtle order parameters.

I n crystalline solids, the combination of spin-orbit coupling and crystal electric fields places strong constraints on the shape of localized electronic wavefunctions 1 . The quantum mechanically defined multipole moments provide a useful measure of the resulting complex angular distribution of the magnetization and charge densities 2,3 . Most conventional broken symmetry phases in solids involve the magnetic dipole moment of the electron. Remarkably, a large class of strongly correlated electron materials display non-trivial higher order multipolar moments, e.g., quadrupolar or octupolar moments, whose fluctuations and ordering leads to a rich variety of phases, such as quadrupolar heavy Fermi liquids 4-6 , superconductivity [7][8][9] , and unusual multipolar symmetry-broken phases 2,3,10-13 . While multipolar ordered phases fall under the purview of the celebrated Landau paradigm of symmetry-broken phases, they have been termed as so-called 'hidden orders': mysterious phases of matter whose orderings are invisible to conventional local probes (such as neutron scattering or magnetic resonance), but are remarkably still known to exist as their onset triggers non-analytic signatures in thermodynamic measurements 4,[14][15][16][17] . Studying the mysterious ordering patterns of higher order multipoles is also often rendered challenging since they typically coexist with conventional dipolar moments. Examples of such symmetry breaking which are of great interest include spin-nematic order 18 in spin S ≥ 1 quantum magnets, quadrupolar charge order in transition metal oxides, and higher multipolar order in actinide dioxides, such as NpO 2 19 , and f-electron heavy fermion materials 20 , such as URu 2 Si 2 21-29 and UBe 13 [30][31][32] . The quest to probe such orders has led to novel experimental techniques, e.g., elasto-resistivity [33][34][35] to elucidate the quadrupolar order associated with orbital nematicity in the iron pnictides. A broad understanding of the nature of these symmetry-broken phases, and means to definitively demonstrate their existence, has proven to be a challenging, yet stimulating, endeavour for both theory and experiment.
Our work is motivated by a recent series of experiments on the Pr-based cage compounds Pr(Ti,V,Ir) 2 (Al,Zn) 20 which form an ideal setting to study multipolar moments and associated hidden orders 8,17,[36][37][38] . In these systems, the 4f 2 electrons of Pr 3+ ions subject to CEFs host a ground non-Kramers doublet with solely higher-order moments: quadrupoles (O 20 and O 22 ) and octupole (T xyz ) 14,39 . Uncovering and understanding the pattern of multipolar ordering across this family of materials remains an important open problem.
The nature of the quadrupolar ordering in these cage compounds has been indirectly examined with a few techniques 40,41 such as ultrasound experiments [42][43][44][45] (indicating softening of elastic modulus at quadrupolar ordering temperature, T Q ), as well as NMR measurements (where the magnetic field-induced dipole moment is strongly dependent on the underlying quadrupolar phase 46 ). More recently, magnetostriction and thermal expansion strain experiments 47 have also lent themselves as possible probes to study the transitions and the underlying quadrupolar phase. By contrast, the octupolar ordered state has continued to remain an elusive phase of matter, with only indirect hints of its existence from NMR 48 and μSR 49 measurements, but as yet no direct probe to reveal its existence 50 . More recently, some of us (A.S. and S.N., unpublished) have begun experiments to study angle-dependent magnetostriction, the change in sample length induced by a magnetic field which can point along various crystalline directions, in a wide class of materials with multipolar degrees of freedom.
In this work, motivated by these experiments, we theoretically discuss how magnetostriction provides a novel means to directly probe multipolar order parameters. The central observation of this paper is that an applied magnetic field allows for a linear coupling between lattice strain fields and a uniform octupole moment which depends strongly on the applied field direction. In the absence of a dipolar moment, this enables measurements of the magnetostriction to directly reveal the hidden octupolar order parameter. We investigate such field-scaling behaviour of the magnetostriction for various magnetic field directions by employing a symmetry-based Landau theory, which allows us to highlight the universal aspects of the physics and to reveal its applicability to a broader class of materials. Specifically, our Landau theory permits both antiferro-quadrupolar ordering (AFQ) and ferro-octupolar ordering (FO), and we examine our theory along different field directions in three temperature regimes. Denoting the quadrupolar and octupolar transition temperatures as T Q and T O , respectively, we consider the regimes (i) the paramagnetic phase above both transition temperatures where the system exhibits pure quadrupolar order, and (iii) below both ordering temperatures (T < T Q ; T O ) where the system features coexisting quadrupolar and octupolar orders. We make definite predictions for all possible combinations of length change and magnetic field directions, which can be tested in future experiments.

Results
Magnetostriction as a probe of multipolar ordering. Our studies predict a linear-in-h scaling behaviour for particular length changes, for T < T O . The coefficient of the linear-in-h term, i.e. the "magnetostriction coefficient", is directly proportional to the ordered ferro-octupolar moment, thus providing a clear and distinct means to directly probe this order parameter. This linearin-h behaviour appears for magnetic fields applied along the [111] and [100] directions. For other magnetic field (and length change) directions, we predict the signature of quadrupolar ordering as a constant plus quadratic-in-h scaling behaviour in the length change; although the scaling behaviour explicitly involves the FQ order, the AFQ order parameter can be inferred from the FQ, as it scales (to leading order) as the square root of the FQ order parameter. We present our theoretical predictions for the scaling behaviours in Table 1 for a variety of magnetic field and length change directions. A quick way to see this linearin-h result is to note that the elastic energy of a strained cubic crystal is given by 51,52 where the crystal's deformation is described by the components of the strain tensor ϵ ik , and c ij is the elastic modulus tensor describing the stiffness of the crystal. Note that we use the normal modes of the cubic lattice to write the elastic energy in this elegant form. Here c B is the bulk modulus, ϵ B ϵ xx þ ϵ yy þ ϵ zz is the volume expansion of the crystal, ϵ ν ð2ϵ zz À ϵ xx À ϵ yy Þ= ffiffi ffi 3 p and ϵ μ ðϵ xx À ϵ yy Þ are lattice strains that transform as the Γ 3g irreducible representation (irrep.) of the O h group, and the offdiagonal strain components transform as the Γ 5g irrep. of O h group. The subscript g indicates even under time-reversal and spatial inversion (parity). Knowing ϵ ij determines the fractional length change along thel-axis via ðΔL=LÞ ℓ ¼ P ij ϵ ij'i'j , where' i is the i th component of unit vectorl; Supplementary Note 1 provides a more detailed discussion on this expression. As discussed below, an applied magnetic field enables a linear coupling between the strain field and the time-reversal-odd ferro-octupolar moment, m, via a term in the free energy , this leads to ðΔL=LÞ ð1;1;1Þ ¼ ðϵ xy þ ϵ yz þ ϵ xz Þ=3 and so (ΔL/L) (1,1,1) ∝ (g O /c 44 ) mh. This direct relation between the linear-in-h magnetostriction coefficient and the ferro-octupolar order parameter for a magnetic field along the [111] direction is one of the central results of our paper. Furthermore, we predict a characteristic hysteresis in the octupolar moment and the associated parallel length change as a function of magnetic field, arising from the symmetryallowed cubic-in-h coupling of the magnetic field to the octupolar moment. Very recent (unpublished) experiments on PrV 2 Al 20 indeed appear to find a hysteretic linear-in-field magnetostriction, for a [111] magnetic field, below a transition at T* ≈ 0.65 K. Our theoretical results for magnetostriction in the presence of octupolar order thus lend strong support to the idea that this approach, pursued in recent experiments performed by some of us (A.S. and S.N., unpublished), herald the unambiguous discovery of octupolar order.
The theoretical roadmap which gives rise to this striking result requires (i) the Landau free energy of the multipolar moments, and (ii) coupling between the multipolar moments and lattice strain. We present these key ingredients below.
Landau theory of multipolar order. We present in this section, for the sake of self-containedness and to specify our notation, the Landau theory of multipolar order first introduced in ref. 51 . We focus on key aspects of the model here, and relegate the symmetry-based derivation as well as the complete form of the free energy to Methods; the symmetry transformations of the multipolar moments are given in Supplementary Note 2.
The 4f 2 electrons of Pr 3+ ions in the family of rare-earth metallic compounds Pr(Ti,V,Ir) 2 (Al,Zn) 20 reside on a diamond lattice of cubic space group Fd 3m. Surrounding each Pr 3+ ion is a Frank-Kasper (FK) cage (16 Al atom polyhedra). The crystalline electric field (CEF) of this FK cage, with T d point group symmetry, splits the J = 4 multiplet of the 4f 2 electrons. The ground states are experimentally found to be a non-Kramers doublet, and they transform as the basis states of the Γ 3g irrep. of T d ; here the subscript g(erade) and u(ngerade) denote even and odd under time-reversal, respectively. Moreover, this doublet is energetically well separated from the excited states, and so for energies much lower than this gap (≳50 K 4 ), the Γ 3g doublets form an ideal basis to describe the low energy degrees of freedom. The Γ 3g doublets can give rise to time-reversal even quadrupolar , as well as a time-reversal odd octupolar moment 6 J x J y J z which transforms as Γ 2u (where the overline represents the fully symmetrized product). Using the constructed pseudospin basis ({|↑〉, |↓〉}) from the Γ 3g doublets, allows the multipolar moments to be neatly denoted by an effective pseudospin-1/2 operator τ = (τ x , τ y , τ z ) The perpendicular component of the pseudospin vector τ ⊥ ≡ (τ x , τ y ) denotes the quadrupole moments, while τ z denotes the octupolar moment. We also define the raising/lowering pseudospin operators τ ± = τ x ± iτ y . The ordering of these multipolar degrees of freedom acts as a mean field on the pseudospins, and breaks the degeneracy of the non-Kramers doublet. In order to describe these pseudospinsymmetry-broken phases, we resort to a Landau theory approach, focussing on the following order parameters, Here, angular brackets 〈...〉 denote thermal averages, while the A, B subscripts denote the two sublattices of the diamond lattice. The complex scalars ϕ andφ describe ferroquadrupolar (FQ) and anti-ferroquadrupolar (AFQ) orders, respectively, while the real scalars m andm denote the ferro-octupolar (FO) and antiferrooctupolar (AFO) order parameters. We henceforth use the convention ofφ ¼ jφje iα and ϕ = |ϕ|e iα for the complex order parameters.
In this work, we focus on a system where the primary order parameters are AFQ and FO. As discussed in previous works 53, 54 , the Landau theory of a system with AFQ order necessarily admits a 'parasitic' secondary order parameter FQ. Such mixing does not occur for the octupolar order parameter; motivated by explaining the experiments performed by some of us (A.S. and S.N., unpublished) on PrV 2 Al 20 , we choose to work with only FO order and ignore the AFO order parameter. We can thus construct our Landau theory using the order parameters ϕ,φ, and m, based on the local T d symmetry instilled by the FK cage, Here, the free energies Fφ, F m , and F ϕ denote the independent free energies of the AFQ, FO, and FQ orders, and Fφ ;ϕ;m describes the interactions between the different multipolar order parameters. Figure 1 shows the zero magnetic field phase diagram, depicting both quadrupolar and octupolar transitions; with two primary order parameters AFQ (and its accompanying parasitic FQ moment) and FO ordering at Table 1 Scaling relation for relative length change of system ΔL=L ' along direction ' for magnetic field applied along b n direction For each b n, we present the length change parallel and (the two) perpendicular directions with respect to b n. The FQ moment term is expressed as g Q jϕj Φ 1;2 þ κ 1;2 h 2 due to the even-in-h behaviour of the quadrupolar moment, where Φ1,2 is the zero magnetic field quadrupolar moment which arises due to AFQ spontaneously ordering. Here, the two types of g Q (and κ1,2, Φ1,2) include the complex angledependent parts (α) and the quadrupolar-lattice strain coupling; as described in Supplementary Note 6, there are two possible combinations of the complex angle dependency in g Q , which we denote by the subscripts 1, 2. Since Φ1,2,κ1,2 arise from the parasitic FQ moment, they are diminutive, as compared with the conduction electrons' m is a re-defined octupolar moment, including the octupolar-lattice strain coupling.
critical temperatures of T Q and T O , respectively. The 'kink' in the AFQ (as well as FQ) at the octupolar ordering temperature reflects the non-analytic behaviour of the octupolar moment at its critical temperature. We present in Supplementary Note 3 the values of the Landau parameters used for Fig. 1.
In order to study magnetostriction, it is important to understand how the magnetic field couples to the multipole moments. Due to the lack of magnetic dipole moment supported by the Γ 3g doublet, the magnetic field does not couple linearly to the states. One can derive the low energy magnetic field Hamiltonian by performing second-order perturbation theory in h · J, where the low energy subspace is spanned by the Γ 3g doublet, and the high energy subspace is spanned by the excited triplets Γ 4,5 . This leads to F mag ½ϕ;φ, which involves the quadrupolar moments coupling quadratically to the magnetic field,~h 2 τ x,y . The coupling of the magnetic field to the octupole moment (after performing third-order perturbation theory) is of the form~h x h y h z τ z . The Oðh 3 Þ term is neglected at this stage, and its role is revived in the discussion of hysteresis.
Symmetry allowed coupling of multipoles to lattice modes. We now turn our attention to the problem of coupling the lattice normal modes of the cubic crystal to the multipolar moments. We recall that the cubic crystal structure supports macroscopic normal modes that transform as irreps. of O h , while the Landau free energy of the multipolar moments (F) is constructed subject to symmetries of the local T d environment. The symmetry constraints on F ensure that in principle only select normal modes of the crystal that transform as the irreps. of T d are permitted to couple to the multipolar moments. In the present case, all the cubic normal modes presented in Eq. (1) also transform as irreps. under T d (as can be explicitly verified), and so all of the aforementioned strain modes can participate in the coupling.
Coupling of quadrupolar moment to lattice strain. Coupling between the quadrupolar moments and the lattice normal modes appears as a natural choice, as the quadrupolar moments and the lattice strains are both even under time-reversal. Moreover, both the normal modes fϵ μ ; ϵ ν g and the quadrupolar moments fO 22 ; O 20 g transform as Γ 3g irreps. of T d (the aforementioned lattice normal modes also transform as Γ 3g in O h , as T d is a subgroup of O h ). This similarity in how they transform under T d allows a linear coupling between the aforesaid lattice normal modes and quadrupolar moments. Thus, the Landau free energy of the multipolar moments gets augmented by the following lattice elastic energy and coupling terms to quadrupolar moments, where g Q is the coefficient of coupling between the quadrupolar moments and lattice strain tensors. Note that we include the coupling of the lattice strain to the quadrupole moment on each sublattice. Using the definition of the order parameter ϕ from Eq.
(3), and minimizing F strain;Q with respect to ϵ μ ; ϵ ν yields the total strain for each normal mode Substituting these expressions back into Eq. (4), we find that the strain-optimized F strain;Q ½ϕ ¼ À g 2 Q 2ðc 11 Àc 12 Þ jϕj 2 renormalizes the mass term of the FQ order.
Coupling of octupolar moment to lattice strain. A direct linear coupling between the octupolar moment T xyz and the lattice normal modes is not permitted, as the octupolar moment is odd under time-reversal. However, this potential difficulty can be alleviated by the introduction of the time-reversal odd magnetic field h which assists in the coupling between the lattice degrees of freedom and octupolar moment. Thus, the Landau free energy of the multipolar moments gets augmented by the following lattice elastic energy and the coupling terms to the octupolar moments, where we use the definition of m from Eq. (3), and g O is the coefficient of coupling between the octupolar moment and lattice strain. We also include another symmetry-allowed direct coupling between the magnetic field and the same lattice normal modes (with proportionality constant γ c , equivalent on both sublattices). Physically, this kind of term could arise from the independent coupling of the magnetic field and lattice strain to the conduction electrons (and after integrating out the conduction electrons). We discuss in Supplementary Note 5 how the numerical values of these coupling constants can be obtained from experimental observations in conjunction with our theoretical predictions.
Minimizing with respect to the lattice degrees of freedom yields the following expressions for the (total) lattice strains Substituting the expression for the minimized lattice strains from Eq. (7) into Eq. (6) yields F strain;O ½m, where the mass term of the octupolar moment is renormalized by a term quadratic in h; it also introduces an Oðh 3 Þ coupling term between the octupolar moment and the magnetic field, which renormalizes the coefficient of the already present h x h y h z m from third-order in perturbation theory in h ⋅ J. Length change under magnetic field along various directions. Equipped with the necessary ingredients in the previous subsections, we can now examine the relative length change, ΔL/L, for magnetic fields applied along [100], [110], [111] directions and examine the scaling in magnetic field strength, h. For the sake of clarity, we stress that we consider here the complete Landau theory of multipolar moments coupled to lattice strain fields (after having integrated out the lattice degrees of freedom): F½ϕ;φ; m ¼ F Q;O ½ϕ;φ; m þ F mag ½ϕ;φ þ F strain;Q ½ϕþ F strain;O ½m. The scaling relations can be inferred by substituting the expressions for the (extremized) strain in Eqs. (5) and (7) This equation has a striking conclusion as it pertains to observing hidden order. The mysterious octupolar moment can now be determined (up to a proportionality constant) by measuring the slope of the linear-in-h behaviour of the length change both parallel and perpendicular to magnetic fields applied along the [111] direction. This provides a clear signature for the onset of the octupolar ordering as well as a means to study the general behaviour of the octupolar moment (up to a proportionality constant) with respect to other external variables such as temperature, T. Moreover, we discover that length change parallel to the magnetic field along [111] has (negative) twice the slope of the linear-in-h term and (negative) twice the quadratic background as the length changes perpendicular ℓ ¼ ð1; À1; 0Þ; ð1; 1; À2Þ to the field [111]. Furthermore, from Table 1, the octupolar moment analogously appears in the length change perpendicular to the magnetic field along the [100] direction. Indeed, the sign of linear-in-h coefficient flips for the two presented orthogonal directions. All of these provide distinct means to validate the theory. Next, for magnetic fields along the [110] direction, we find that the length changes parallel, ℓ ¼ ð1; 1; 0Þ, and perpendicular, ℓ ¼ ð1; À1; 1Þ; ðÀ1; 1; 2Þ, to the field are purely quadratic-in-h and do not possess a linear-in-h scaling behaviour. Thus, these length changes (for this choice of magnetic field) do not provide information about the octupolar moment; the quadratic in h behaviour arises from the conduction electrons and/or the quadrupolar moment. We provide in Supplementary Note 4 a justification of the scaling behaviours of the multipolar moments, and in Supplementary Note 6 the corresponding general length change expressions. We note that the scaling behaviours presented here and in Supplementary Note 6 neglect the cubicin-h coupling, which breaks the ℤ 2 symmetry (m ↔ −m) of the octupolar moment. This introduces a 'flip' in the octupolar moment at h = 0 (and at T<T O where m has spontaneously ordered, i.e. m ≠ 0): for h = 0 + , the +|m| solution is 'chosen', and as we crossover to h = 0 − , the now physically distinct −|m| solution is 'chosen' (this is seen in Fig. 2). A similar phenomena is observed in usual ferromagnetism, below the ordering temperature. The neglect of this term is due to the consideration of weak, perturbative magnetic fields in this study. It is likely that this term could become more important (with regard to the scaling behaviour) for larger magnetic fields, but this is beyond the field ranges considered in this work.
Hysteretic behaviour of octupolar ordering. We are motivated in this section by preliminary experimental observations found by some of us (A.S. and S.N., unpublished) of hysteretic behaviour in the length change along the [111] direction below the supposedoctupolar temperature. Hysteresis arises from the existence of domains and the motion of domain walls in the presence of obstructing 'pinning sites', which have not been taken into account in the Landau theory we have studied. In order to incorporate such effects, we adapt the phenomenological approach due to Jiles and Atherton 55,56 which has been used to study hysteresis loops in ferromagnetic and ferroelastic materials. This approach identifies the order parameter (obtained by minimizing the Landau free energy) as its ideal bulk value, where the Landau theory includes a direct coupling u f mh 3 of the ferrooctupolar moment m and the external [111] magnetic field. The total macroscopic octupolar moment (m exp ) is obtained by solving the constructed Jiles and Atherton model, which is heuristically derived in Supplementary Note 7. The key point to note is that the hysteresis loop arises from having a time-reversal odd moment (and domains) coupling to the magnetic field.
We present the calculated hysteresis for T < T O in Fig. 2a. The initial condition chosen to obtain the hysteresis loop is such that at h = 0, the ideal configuration is not being met (i.e. m exp ≠ m); this depicts the realistic scenario of having not all domains aligned in the same direction at h = 0. The depicted hysteresis is reminiscent of the hysteresis in ferromagnets. We obtain the corresponding length  1, 1, 1) direction as shown in Fig. 2b, which for small magnetic fields displays the linear-in-h scaling behaviour. To better observe this linear-in-h scaling in the length change, we present the derivative of the length change with respect to the magnetic field in the inset of Fig. 2b. The linear-in-h scaling behaviour of the length change is more clearly apparent as a constant y-intercept in the inset; the further linear scaling in the inset is due to the background quadratic-in-h scaling behaviour of ΔL/L from the conduction electrons (~γ term).
Furthermore, we note that the field strength h* corresponding to the minimum of the length change [i.e. d(ΔL/L)/dh = 0] provides a threshold above which the conduction electron background dominates over the linear-in-h scaling behaviour. For this particular length change direction, h Ã ¼ . We note that dimensionally h has units of energy (as we have set the Bohr magneton, μ B = 1 here) and the strain tensor ϵ is dimensionless; this implies that the composite quantity g O jmj is dimensionless, while the conduction electron term γ c (which scales like an off-diagonal magnetic susceptibility, from Eq. (6)) has units of (energy) −1 . Dimensional analysis thus suggests γ c~D OS, where (DOS) is the conduction electron density of states at the Fermi level. To proceed further with the other quantities, we note that m itself is a dimensionless quantity; subsequently g O is also dimensionless. This follows from the convention used in Landau theory where m $ T O ÀT T O β , and as such for low temperatures (T → 0) we can take m as an O(1) number. If we also take g O $ Oð1Þ, then from above h*~DOS −1 . Thus, the location of the minimum field is inversely dependent on the conduction electron DOS at the Fermi level: when the conduction electron DOS is small, the minimum field h* is correspondingly large.

Discussion
In this work, motivated by recent and ongoing experiments on Pr (Ti,V,Ir) 2 (Al,Zn) 20 , we have used Landau theory of multipolar orders coupled to lattice strain fields to study magnetostriction in systems with quadrupolar and octupolar orders. Our theoretical results for magnetostriction in the presence of octupolar order appear consistent with recent magnetostriction experiments performed by some of us (A.S. and S.N., unpublished) on PrV 2 Al 20 where the onset of unusual linear-in-field and hysteretic magnetostriction is observed for fields along the [111] direction for T < 0.65 K; in particular, we predict linear-in-h scaling of the length change for length changes (both parallel and perpendicular) to magnetic fields applied along the [111] direction, and also for length changes perpendicular to [100], below T O . Moreover, the coefficient of the linear-in-h term is directly proportional to the octupolar moment, thus giving a distinct signature for the onset of octupolar ordering as well as a means to detect/measure the octupolar moment. In addition, we can qualitatively understand the quadratic-in-field background magnetostriction observed in these experiments; we predict that this scaling arises from the quadrupolar moments and/or direct coupling of the conduction electrons to the external magnetic field and the appropriate lattice normal modes. The summary of all the scaling behaviours is presented succinctly in Table 1.
Our results are broadly applicable to a variety of multipolar orders in cubic systems. For instance, the conclusions here are extendable to the cluster octupolar moments suggested in antiferro-magnetically interacting magnetic moments in pyrochlore iridates 57 . Furthermore, the results presented here can be extended to other more 'typical' probes of multipolar ordering/ fluctuations. For instance, due to the permitted octupolar-strain coupling, we expect to observe elastic constant softening in the elastic constant c 44 at the ordering temperature, T O , under the application of a magnetic field. We note that it is the c 44 elastic constant that softens, as it is the associated elastic constant with the off-diagonal strain normal modes. Similarly, we expect elastoresistivity experiments 58 to be a probe for octupolar susceptibility. The application of an elastic strain with T 2g symmetry (such as xy, xz, or yz) in the presence of a magnetic field would result in an associated anisotropic resistivity (ρ xy , ρ xz , ρ yz ), which will be proportional to the octupolar susceptibility. Finally, we expect that Pr(Ti,V,Ir) 2 (Al,Zn) 20 compounds will possess the characteristic low frequency Raman quasielastic peak 59 , associated with octupolar fluctuations; specifically, under the application of a magnetic field along the [001] direction, we expect the quasielastic peak to appear in the xy symmetry Raman spectra.
In terms of future work, an interesting avenue to explore is that of the coupling of the conduction electrons to the multipolar moments, as well as to the lattice strain and magnetic field. In particular, the origin of the conduction electron term in Eq. (6), introduced in our phenomenological model from symmetry arguments, is a fascinating direction to explore (as well as potential other terms arising from conduction electrons). We discuss in Supplementary Note 8 a possible origin of the conduction electron term of Eq. (6). Understanding the nature and role of the conduction electrons will also help shed light on the quantum critical behaviour and superconductivity in such multipolar Kondo lattice systems 4,7,8,32,[60][61][62][63][64] .

Methods
Non-Kramers ground states of Pr 3+ ions. The ground states are experimentally found to form a non-Kramers doublet written in |J z 〉 basis as Constructing a pseudospin basis ({|↑〉, |↓〉}) from the Γ 3g doublets as allows the multipolar moments to be succinctly written as the effective pseudospin-1/2 operator τ = (τ x , τ y , τ z ) in Eq. (2) in the main text. The local T d symmetry instilled by the FK cage provides a constraint on the possible terms permitted in the Landau theory. The generating elements of T d are S 4z (improper rotation of π/2 about the b z-axis) and C 31 (rotation of 2π/3 about the body diagonal [111] axis). In addition to these point group symmetries, we also require that the terms in the Landau theory be invariant under spatial inversion about the diamond bond centre I (which swaps the A and B sublattices), as well as time-reversal Θ. The behaviour of the multipolar moments under these symmetry constraints is detailed in Supplementary Table 1. As described in the main text, we construct our Landau theory using the order parameters ϕ,φ, and m.
Interacting multipolar orders. Equipped with the symmetry knowledge from Supplementary Table 1 we can now write down the Landau free energy for this particular multipolar ordered system as Here, the free energies Fφ, F m , and F ϕ denote the independent free energies of the AFQ, FO, and FQ orders. Settingφ ¼ jφje iα and ϕ = |ϕ|e iα , we get The first two terms in Eqs. (12)- (14), in square brackets, are the usual mass and quartic interaction terms for AFQ, FO and FQ order parameters. We choose tφ ¼ ðT À T Q Þ=T Q , and t m ¼ ðT À T O < T Q , where T denotes the temperature. Focussing on the mass term alone, decreasing T will thus lead to an anti-ferroquadrupolar order for T < T Q , and a lower temperature transition into a state with coexisting ferro-octupolar order when T < T temperatures will be affected by the interplay of the two order parameters; in particular, the true octupolar transition T O will be renormalized from its bare value T ð0Þ O due to the onset of quadrupolar order (besides fluctuation effects which we do not consider here). A measure of how close the two transition temperatures are to each other is provided by the ratio (T Q − T O )/(T Q + T O ). Finally, since FQ is not considered to be a primary order parameter, we choose a large positive mass term, t ϕ . The remaining non-trivial terms in Eqs. (12) and (14) are the unusual sixth order and cubic "clock" terms, with respective coefficients wφ and v ϕ , which fix the phases of the AFQ and FQ order parameters. We set lφ > jwφj to ensure that the free energy is bounded from below.
The couplings between the different multipolar order parameters are encapsulated in Fφ ;ϕ;m , namely between AFQ and FQ moments (g 1 , g 2 ), and between the quadrupolar and the octupolar moments ðu ϕm ; uφ ;m Þ where the term g 1 is a symmetry-allowed cubic term.
Coupling of magnetic field to multipolar moments. Due to the lack of magnetic dipole moment supported by the Γ 3g doublet, the magnetic field does not couple linearly to the states. One can derive the low energy magnetic field Hamiltonian by performing second-order perturbation theory in h ⋅ J, where the low energy subspace is spanned by the Γ 3g doublet, and the high energy subspace is spanned by the excited triplets Γ 4,5 ; here h has units of energy as we have set the Bohr magneton, μ B = 1. This leads to In the above Eq. (16), h = (h x , h y , h z ) with |h| = h, and γ 0 À14 3ΔðΓ 4 Þ þ 2 ΔðΓ 5 Þ , where Δ (Γ 4 ), Δ(Γ 5 ) are the gaps between the low energy doublets and the corresponding triplet states at zero magnetic field. The effective coupling to the ferroquadrupolar Based on the form of the coupling in Eq. (16), we infer that ψ H transforms identically to ϕ under the relevant symmetries. Going to third-order in perturbation theory leads to a further O(h 3 ) coupling of the magnetic field to octupole moment of the form~h x h y h z τ z .
Thus, the symmetry-allowed effective magnetic field coupling to the quadrupolar moments is where jψ H j ¼ γ 0 4 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3ðh 2 x À h 2 y Þ 2 þ ð3h 2 z À h 2 Þ 2 q , and tanðθ H Þ ¼ 1 ffiffi 3 p 3h 2 z Àh 2 ðh 2 x Àh 2 y Þ . The first (second) line in Eq. (17) is the symmetry-allowed coupling to the AFQ (FQ). The third line involves couplings permitted due to pure symmetry reasons that renormalize the mass terms of the AFQ and FQ. Physically they arise from conduction electron mediated magnetic couplings (having integrated out the conduction electrons); similar coupling to the octupolar moment is also permitted [~h 2 m 2 ], which is formally introduced via the magnetic field assisted coupling of the octupolar moment to the lattice strain. In the main text, we discuss magnetic fields applied along the [100], [110] and [111] directions. For clarity, we present the value for |ψ H | and θ H for the magnetic field directions discussed in subsequent sections in Table 2. In the presence of the magnetic field, it is possible for additional couplings between the quadrupolar and octupolar moments to be induced, such as $ m 2 cosðα À θ H Þjφjjψ H j ; ð18Þ $ m 2 sinðθ H þ 2αÞjφj 2 jψ H j ; ð19Þ These terms are merely the usual quadratic-in-field coupling to the quadrupolar moment (Eqs. (16) and (17)) with m 2 multiplied into it. Due to symmetry constraints, we cannot have terms which are linear in the octupolar, quadrupolar and magnetic field (breaks C 31 symmetry). These above terms do not affect the leading scaling behaviour of the magnetostriction, as they have the same order of h as previous terms in the free energy. Specifically, the terms are quadratic-in-h and can be thought of as renormalizing the mass term of the octupolar moment. We recall that the octupolar mass term already contains a quadratic-in-h term, which arose from integrating out the elastic strain in Eq. (6), and so these new terms merely modify the coefficient of the previous quadratic-in-h expressions/terms. The Landau theory is numerically minimized using standard minimization/ optimization schemes. The hysteresis differential equation is numerically solved using Runge-Kutta 4th order methods. We use the initial condition of m ir = 0 for h = 0 to obtain the depicted solution, with k = 100, α = 10 −3 , c = 0.01.

Data availability
All relevant data are available upon reasonable request to the corresponding author.

Code availability
All relevant codes are available upon reasonable request to the corresponding author.
, the magnetic field does not directly couple to the quadrupolar moments, but can do so vias H and s H .