Transverse fields to tune an Ising-nematic quantum critical transition

TL;DR

This paper demonstrates how transverse fields, such as strain or magnetic fields, can tune orbital-nematic order to a quantum critical point, revealing intertwined order parameters through non-trivial quantum relations.

Contribution

It introduces a practical method to tune orbital-nematic order to quantum criticality using transverse fields, linking different order parameters via quantum commutation relations.

Findings

Transverse strain or magnetic field can tune nematic order to a quantum critical point.

Non-trivial commutation relations intertwine different order parameters.

Berry-phase terms describe the effective field theory of the transition.

Abstract

The paradigmatic example of a continuous quantum phase transition is the transverse field Ising ferromagnet. In contrast to classical critical systems, whose properties depend only on symmetry and the dimension of space, the nature of a quantum phase transition also depends on the dynamics. In the transverse field Ising model, the order parameter is not conserved and increasing the transverse field enhances quantum fluctuations until they become strong enough to restore the symmetry of the ground state. Ising pseudo-spins can represent the order parameter of any system with a two-fold degenerate broken-symmetry phase, including electronic nematic order associated with spontaneous point-group symmetry breaking. Here, we show for the representative example of orbital-nematic ordering of a non-Kramers doublet that an orthogonal strain or a perpendicular magnetic field plays the role of the transverse field, thereby providing a practical route for tuning appropriate materials to a quantum critical point. While the transverse fields are conjugate to seemingly unrelated order parameters, their non-trivial commutation relations with the nematic order parameter, which can be represented by a Berry-phase term in an effective field theory, intrinsically intertwines the different order parameters.

Figures

Tables

Table 1. Table S1: Transformation laws of ( Φ 1 , Φ 2 , Φ 3 ) subscript Φ 1 subscript Φ 2 subscript Φ 3 (\Phi_{1},\Phi_{2},\Phi_{3}) under point group transformations and time reversal.

Symmetry	${\Phi }_1$	${\Phi }_2$	${\Phi }_3$
$M_h$	$-$	$-$	$+$
$M_d$	$+$	$-$	$-$
$R_{\pi /2}=M_h{M}_d$	$-$	$+$	$-$
$𝒯$	$+$	$-$	$+$

Table 2. Table S2: Character Table for D 4 ​ h subscript 𝐷 4 ℎ D_{4h} , with the reducible representations Γ 6 superscript Γ 6 \Gamma^{6} and for general J 𝐽 J , Γ J superscript Γ 𝐽 \Gamma^{J} shown at the bottom. Determining the crystal field splitting of the J = 6 𝐽 6 J=6 Hund’s rule ground state amounts to determining how the Γ 6 superscript Γ 6 \Gamma^{6} representation (which is reducible in D 4 ​ h subscript 𝐷 4 ℎ D_{4h} ) is decomposed into irreducible representations of D 4 ​ h subscript 𝐷 4 ℎ D_{4h} .

Irrep	BSW notation111due to Bouckaert, Smoluchowski & Wigner (1936)	$E$	$2C_4$	$C_2$	$2C_2^{\prime }$	$2C_2^{\prime \prime }$	$i$	$2S_4$	${\sigma }_h$	$2{\sigma }_v$	$2{\sigma }_d$
$A_{1g}$	${\Gamma }_1^{+}$	1	1	1	1	1	1	1	1	1	1
$A_{2g}$	${\Gamma }_2^{+}$	1	1	1	-1	-1	1	1	1	-1	-1
$B_{1g}$	${\Gamma }_3^{+}$	1	-1	1	1	-1	1	-1	1	1	-1
$B_{2g}$	${\Gamma }_4^{+}$	1	-1	1	-1	1	1	-1	1	-1	1
$E_g$	${\Gamma }_5^{+}$	2	0	-2	0	0	2	0	-2	0	0
$A_{1u}$	${\Gamma }_1^{-}$	1	1	1	1	1	-1	-1	-1	-1	-1
$A_{2u}$	${\Gamma }_2^{-}$	1	1	1	-1	-1	-1	-1	-1	1	1
$B_{1u}$	${\Gamma }_2^{-}$	1	-1	1	1	-1	-1	1	-1	-1	1
$B_{2u}$	${\Gamma }_4^{-}$	1	-1	1	-1	1	-1	1	-1	1	-1
$E_u$	${\Gamma }_5^{-}$	2	0	-2	0	0	-2	0	2	0	0
${\Gamma }^{J=6}$		13	-1	1	1	1	13	-1	1	1	1
${\Gamma }^J$		$\left(2J+1\right)$	${\left(-1\right)}^{\text{floor}\left(J/2\right)}$	${\left(-1\right)}^J$	${\left(-1\right)}^J$	${\left(-1\right)}^J$	${\left(-1\right)}^J\left(2J+1\right)$	${\left(-1\right)}^{J+\text{floor}\left(J/2\right)}$	1	1	1

Equations

$${\Phi }_{\alpha }(\vec{R})\equiv \sum _{a,a^{\prime }}{c}_{a,\vec{R}}^{†}{\tau }_{a,a^{\prime }}^{(\alpha )}{c}_{a^{\prime },\vec{R}}^{},$$

$$\left[{\Phi }_{\alpha }(\vec{R}),{\Phi }_{\beta }({\vec{R}}^{\prime })\right]=i{\delta }_{\vec{R},{\vec{R}}^{\prime }}{ϵ}_{\alpha \beta \gamma }{\Phi }_{\gamma }(\vec{R})$$

$${\Phi }_3(\vec{R})\equiv {\hat{n}}_x(\vec{R})-{\hat{n}}_y(\vec{R}).$$

$$\mathcal{H}\to \mathcal{H}-\sum _{\vec{R}}\mathbf{h}\cdot \mathbf{\Phi }(\vec{R})$$

$$h_3$$

$$h_1$$

$$h_2$$

$${\Phi }_3(\vec{R})\equiv \sum _{{\vec{R}}^{\prime }}f(\vec{R}-{\vec{R}}^{\prime })\left[c_{\sigma ,\vec{R}}^{†}{c}_{\sigma ,{\vec{R}}^{\prime }}+H.C.\right]$$

$$\left[{\Phi }_3(\vec{R}),\sum _{{\vec{R}}^{\prime }}{\Phi }_1({\vec{R}}^{\prime })\right]=0.$$

$$S_{\mathrm{top}}[\mathbf{\Phi }]={\int }_0^{1/T}d\tau \mathbf{A}(\mathbf{\Phi })\cdot \dot{\mathbf{\Phi }}.$$

$$S[\mathbf{\Phi }]=i\sum _{\vec{R}}\Theta [{\mathbf{\Phi }}_{\vec{R}}]+S_{LGW}[\mathbf{\Phi }]$$

$$\Theta [\mathbf{\Phi }]={\int }_0^{\beta }d\tau \mathbf{A}(\mathbf{\Phi })\cdot \dot{\mathbf{\Phi }}.$$

$${\nabla }_{\mathbf{\Phi }}\cdot \mathcal{B}=\sum _i\frac{q_i}{2}\delta (\mathbf{\Phi }-{\mathbf{\Phi }}_i),$$

$$\mathbf{A}(\mathbf{\Phi })=q_0\left[\frac{\Phi -{\Phi }_1}{4\pi \Phi }\right]\left[\frac{{\Phi }_2{\mathbf{e}}_3-{\Phi }_3{\mathbf{e}}_2}{{\Phi }_2^2+{\Phi }_3^2}\right],$$

$$\mathbf{A}(\mathbf{\Phi })=q_0\left[\frac{1}{8\pi {\Phi }_1^2}\right]\left[{\Phi }_2{\mathbf{e}}_3-{\Phi }_3{\mathbf{e}}_2\right]\left[1+\mathcal{O}\left(\frac{{\Phi }_2^2+{\Phi }_3^2}{{\Phi }^2}\right)\right].$$

$$S^{\mathrm{eff}}[{\Phi }_3]=S_0[{\Phi }_3]+\frac{q_0^2{\chi }_2}{128|{\bar{\Phi }}_1{|}^4}\sum _{\vec{R}}{\int }_0^{\beta }d\tau |{\dot{\Phi }}_3{|}^2$$

$$H_1$$

$$H_2$$

$${ϵ}_{1\mathbf{k}}$$

$${ϵ}_{2x\mathbf{k}}$$

$${ϵ}_{2y\mathbf{k}}$$

$$s_{\mathbf{k}}$$

$$d_{\mathbf{k}}$$

$$H_{1\text{int}}=\frac{-V}{N}\sum _{\mathbf{k},{\mathbf{k}}^{\prime }}{d}_{\mathbf{k}}{d}_{{\mathbf{k}}^{\prime }}{c}_{\mathbf{k}}^{†}{c}_{\mathbf{k}}{c}_{{\mathbf{k}}^{\prime }}^{†}{c}_{{\mathbf{k}}^{\prime }}$$

$$H_{2\text{int}}=-\frac{V}{N}\sum _{\mathbf{k},{\mathbf{k}}^{\prime }}\left({\Psi }_{\mathbf{k}}^{†}{\sigma }_z{\Psi }_{\mathbf{k}}\right)\left({\Psi }_{{\mathbf{k}}^{\prime }}^{†}{\sigma }_z{\Psi }_{{\mathbf{k}}^{\prime }}\right)$$

$${\varphi }_1$$

$${\varphi }_2$$

$$H_{1,\text{MF}}$$

$$H_{2,\text{MF}}$$

$${\varphi }_1=\frac{V}{N}\sum _{\mathbf{k}}\left(\cos k_x-\cos k_y\right)f\left({\epsilon }_{1\mathbf{k},\text{MF}}\right)$$

$${\varphi }_2=\frac{V}{N}\sum _{\mathbf{k}}\frac{2\left({\varphi }_2+t_2{d}_{\mathbf{k}}\right)}{\sqrt{4(t_2{d}_{\mathbf{k}}+{\varphi }_2{)}^2+{\epsilon }_{xy}^2}}\left[f({\epsilon }_{-\mathbf{k}})-f({\epsilon }_{+\mathbf{k}})\right],$$

$${\epsilon }_{\pm \mathbf{k}}=-2t_1{s}_{\mathbf{k}}\pm \sqrt{4(t_2{d}_{\mathbf{k}}+{\varphi }_2{)}^2+{\epsilon }_{xy}^2}$$

$${\Gamma }^{J=6}=\sum _i{a}_i{\Gamma }_i$$

Full text

Transverse fields to tune an Ising-nematic quantum

critical transition

Akash V. Maharaj1,2, Elliott W. Rosenberg3, Alexander T. Hristov1,2, Erez Berg4,5, Rafael M. Fernandes6, Ian R. Fisher2,3, and Steven A. Kivelson1,2

1Department of Physics, Stanford University, Stanford, CA 94305, USA

2Stanford Institute for Materials and Energy Sciences, SLAC National Accelerator Laboratory, Menlo Park, California 94025, USA

3Department of Applied Physics, Stanford University, Stanford, CA 94305, USA

4Department of Condensed Matter Physics, Weizmann Institute of Science, Rehovot, Israel 761001

5Department of Physics, University of Chicago, Chicago, IL 60637, USA

6School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA

The paradigmatic example of a continuous quantum phase transition is the transverse field Ising ferromagnet. In contrast to classical critical systems, whose properties depend only on symmetry and the dimension of space, the nature of a quantum phase transition also depends on the dynamicsSachdev (2007); Sondhi et al. (1997). In the transverse field Ising model, the order parameter is not conserved and increasing the transverse field enhances quantum fluctuations until they become strong enough to restore the symmetry of the ground state. Ising pseudo-spins can represent the order parameter of any system with a two-fold degenerate broken-symmetry phase, including electronic nematic order associated with spontaneous point-group symmetry breakingKivelson et al. (1998). Here, we show for the representative example of orbital-nematic ordering of a non-Kramers doublet that an orthogonal strain or a perpendicular magnetic field plays the role of the transverse field, thereby providing a practical route for tuning appropriate materials to a quantum critical point. While the transverse fields are conjugate to seemingly unrelated order parameters, their non-trivial commutation relations with the nematic order parameter, which can be represented by a Berry-phase term in an effective field theorySenthil and Fisher (2006), intrinsically intertwines the different order parametersFradkin et al. (2015).

Ising nematic states are increasingly recognized as important members of the phase diagrams of several families of strongly correlated electronic materials. In the iron-based superconductors, electronic nematic phases are unambiguously presentChuang et al. (2010); Chu et al. (2010); Dusza et al. (2011); Yi et al. (2011); Lu et al. (2014); Fernandes et al. (2014), and the maximum superconducting transition temperature often correlates with the location of a putative quantum critical point with a nematic character. In the cuprate superconductors, there is a wealth of indirect evidenceAndo et al. (2002); Howald et al. (2003); Hinkov et al. (2008); Lawler et al. (2010); Cyr-Choinière et al. (2015) that electronic anisotropy is enhanced in the underdoped regions of the phase diagram, with a possible quantum critical point near optimal doping. As it is breaks a discrete rotational symmetry, Ising-nematic order is relatively robust to disorder and low dimensionality. Moreover, nematic critical fluctuations have zero momentum, and so couple to (nearly) all low energy fermionic quasiparticles. Consequently, there has been considerable theoretical interest in the role of nematic quantum critical fluctuations as a route to non-Fermi liquid metallic behavior and even superconductivityKim and Kee (2004); Maier and Scalapino (2014); Lederer et al. (2015); Metlitski et al. (2015); Kang and Fernandes (2016).

Identifying appropriate means to tune systems through a continuous nematic quantum phase transition is therefore of considerable importance. Previous studies have tuned nematicity through magnetic fields, hydrostatic pressure, or chemical composition. The latter two parameters, however, cannot be varied continuously during an experiment. Moreover, chemical substitution frequently changes the band filling and inevitably affects the degree of disorder, potentially confounding attempts to delineate the roles played by individual variables. Here, we demonstrate that orthogonal strain or perpendicular magnetic fields promote quantum fluctuations that suppress a nematic phase transition, opening a promising avenue to tune to a nematic quantum critical point.

This effect bears a close resemblance to the physics of the transverse-field Ising model. To construct the nematic pseudo-spins, consider a tetragonal material with electronic nematic order associated with the splitting of a non-Kramers doublet. The corresponding Wannier functions of the orbital doublet at each site transform according to the $E_{g}$ irreducible representation, and therefore have either $xz$ symmetry (represented by the orbital index $a=x$) or $yz$ symmetry (with orbital index $a=y$). To simplify the discussion, and without loss of generality, we consider spinless electrons. In the Hilbert space corresponding to a single site $\vec{R}$ with an orbital doublet, any operator can be expressed as a linear combination of the total number operator and the vector pseudo-spin operators

[TABLE]

where $c_{a,\vec{R}}^{\dagger}$ creates an electron in the Wannier orbital of symmetry $a$, $\tau_{a,a^{\prime}}^{(\alpha)}$ are the Pauli matrices, and $\alpha=1,$ 2, 3. The components $\Phi_{\alpha}(\vec{R})$ satisfy the same commutation relations as a spin operator,

[TABLE]

where $\epsilon_{\alpha\beta\gamma}$ is the Levi-Civita symbol. In terms of the orbital density $\hat{n}_{a}\equiv c_{a,\vec{R}}^{\dagger}c_{a,\vec{R}}$, $\Phi_{3}$ can be identified as the difference in occupancy of $x$ and $y$ orbitals (Figure 1)

[TABLE]

We identify $\Phi_{3}$ as the nematic order parameter that breaks the equivalence between the $x$ and $y$ axis; in group-theory language, it has $B_{1g}$ ($x^{2}-y^{2}$) symmetry. The phase in which $\left\langle\Phi_{3}\right\rangle\neq 0$ breaks tetragonal symmetry but preserves horizontal and vertical mirror symmetries. By changing the orbital basis, it is easy to see that $\Phi_{1}$ is associated with breaking the equivalence between the two diagonals $x\pm y$, and corresponds to a nematic order parameter with $B_{2g}$ ($xy$) symmetry, which breaks tetragonal symmetry but preserves diagonal mirror symmetries. Finally, $\Phi_{2}$ is identified with the difference in occupancy of $x\pm iy$ orbitals, and corresponds to an orbital magnetic moment that preserves tetragonal symmetry but breaks time-reversal symmetry (see Figure 1).

Expressed in terms of local multipole moments, $\Phi_{3}$ and $\Phi_{1}$ correspond to electric quadrupole moments, while $\Phi_{2}$ corresponds to a magnetic dipole moment, immediately delineating the appropriate conjugate fields for each of the three components of the vector operator, as shown in Figure 1. The combination of time-reversal and point group symmetries forbid terms in the Hamiltonian, $\mathcal{H}$, that are odd functions of $\Phi_{\alpha}$. However, application of symmetry breaking strains $\varepsilon_{xx}-\varepsilon_{yy}$ (corresponding to unequal lattice distortions along $x$ and $y$) and $\varepsilon_{xy}$ (corresponding to a shear distortion of the lattice) or a uniform magnetic field perpendicular to the plane $H_{z}$ induce linear terms of the form

[TABLE]

where, to leading order:

[TABLE]

The analogy with the transverse field Ising model is now apparent. The “longitudinal” field, $h_{3}$, is a symmetry breaking field: in a phase with $x^{2}-y^{2}$ symmetry, even an infinitesimal $h_{3}$ lifts the two-fold degeneracy of the state, selecting the phase in which $h_{3}\langle\Phi_{3}\rangle>0$, so a finite $h_{3}$ necessarily smears the nematic transition. On the other hand, the other two components of ${\bm{h}}$ behave as “transverse” fields. They reduce the symmetry of the Hamiltonian while preserving the symmetry $\Phi_{3}\to-\Phi_{3}$, which permits a well-defined nematic transition in the presence of non-zero $h_{1}$ and/or $h_{2}$. However, large enough values of the transverse fields will preclude nematic order, since the commutation relations in Eq. 2 imply that a pseudo-spin with a well-defined value of $\Phi_{1}$ or $\Phi_{2}$ must be highly uncertain in $\Phi_{3}$. This makes it possible to induce a $B_{1g}$ nematic quantum phase transition by applying shear strain or magnetic field (see Fig. 2).

An informative realization of Ising nematic order, which clearly illustrates the effect of transverse fields, is ferroquadrupolar order in $4f$ intermetallic compoundsGehring and Gehring (1975); Morin and Schmitt (1990); Santini et al. (2009). In their trivalent state, the rare earth elements Ce - Yb have a partially filled $4f$ orbital whose small extent implies that its electronic states are effectively localized. A hierarchy of energy scales then determines the character of the ground state. Focusing first on a single site, Hund’s rules in the limit of strong spin-orbit coupling dictate that the total angular momentum $J$ is the only good quantum number. The local crystal electric field (CEF) then acts as a perturbation, splitting the $(2J+1)$-fold degenerate Hund’s rule ground state. The CEF eigenstates are linear combinations of these states with their character determined by the point group symmetryGehring and Gehring (1975). For the specific case of intermetallic systems the additional conduction electrons mediate not only an effective interaction between the local magnetic moments but also between the local charge distributions via a generalization of the Ruderman-Kittel-Kasuya-Yosida (RKKY) exchange mechanismSchmitt and Levy (1985); Morin, P. and Schmitt, D. (1988); Aléonard and Morin (1990); Morin and Schmitt (1990). At the two “ends” of the $4f$ series (corresponding to Ce and Pr on the left, or Tm and Yb on the right), Hund’s rules necessarily imply that these elements have relatively large orbital angular momenta (and hence large electric quadrupole moments), but relatively small total spin (and hence small magnetic exchange energies). For these elements, quadrupole-quadrupole interactions can exceed spin-spin interactions. Materials with such constituents can exhibit ordered phases in which the local $4f$ orbitals develop a spontaneous quadrupole moment at a higher temperature than any long range magnetic order.

This is precisely the case for the tetragonal intermetallic compounds TmAg2 and TmAu2, which display ferroquadrupolar order with $B_{1g}$ symmetry, below a critical temperature of approximately 5K and 7K respectivelyMorin and Rouchy (1993); Kosaka et al. (1998); Morin et al. (1999). Near the ferroquadrupolar transition, the only states which are both thermally populated and significant to the transition are the two CEF eigenstates belonging to the non-Kramers doublet, labeled “$E_{g}$ doublet” in Fig. 3. The states of this doublet can be treated as a pseudo-spin $1/2$; projecting the Stevens operators $O_{2}^{2}=J_{x}^{2}-J_{y}^{2}$, $P_{xy}=(J_{x}J_{y}+J_{y}J_{x})$, and $J_{z}$ to this doublet yields operators with the same commutation relations defined in Eq. 2 (see supplementary material for more details). Consequently, $\varepsilon_{xy}$ and $H_{z}$ are appropriate transverse fields to tune the ferroquadrupolar transition to a quantum critical point. Indeed, a magnetic field oriented along the crystalline $c$-axis has been shown to suppress quadrupole orderMorin and Rouchy (1993); Morin et al. (1999), though the effect of shear strain predicted here has yet to be demonstrated.

Our preceding discussion focused on nematic systems whose relevant low-energy states are a non-Kramers doublet. In this case, we emphasized the commutation relations in Eq. 2. These imply the transverse field enhances quantum fluctuations of the order parameter and, when sufficiently large, drives the nematic transition temperature $T_{c}\to 0$. However, this result does not generalize to all possible types of Ising nematic order, as symmetry does not guarantee the availability of an external, experimentally applicable transverse field.

Consider, for example, the case of a tetragonal, single band metal undergoing a Pomeranchuk transition at which the Fermi surface undergoes a symmetry-breaking nematic distortion. In this case, an $x^{2}-y^{2}$ nematic order parameter must involve at least nearest-neighbor sites,

[TABLE]

where $f$ is the “$d$-wave” form factor: $f(\pm\hat{x})=-f(\pm\hat{y})=1$ and $f(\vec{R})=0$ otherwise. From a symmetry perspective, this is no different than the situation analyzed above. Indeed, a field with the symmetry of $\Phi_{1}$ can be constructed by choosing $f(\vec{R}-\vec{R}^{\prime})$ to have $xy$ character in Eq. 6.

However, because the system has only a single $s$-orbital, it is not possible to write down an orbital magnetic moment, i.e. there is no $\Phi_{2}(\vec{R})$. More importantly, while the commutation relations of the nematic operators are moderately complicated , they commute when averaged over sites.

[TABLE]

Consequently, although external strain still couples to $\vec{\Phi}$ as in Eqs. 4 and Transverse fields to tune an Ising-nematic quantum critical transition (with $h_{2}=0$), Eq. 7 implies that $h_{1}$ does not enhance quantum fluctuations of $\Phi_{3}$. Of course, changing the value of any term in the microscopic Hamiltonian which does not explicitly break a relevant symmetry will generally result in a shift in $T_{c}$. However, while $T_{c}$ will be a parametric function of $h_{1}$, it does not necessarily decrease for increasing strain.

To illustrate this crucial point, we calculated the Hartree-Fock phase diagram of such a one-band model undergoing a Pomeranchuk transition characterized by a non-zero $\Phi_{3}$ (see Supplementary Material). Application of shear strain (i.e. a finite $h_{1}$) decreases the second-order transition temperature. However, as it approaches $T=0$, the nematic transition becomes generically first-order, preempting a quantum critical point. By contrast, when the system undergoing the Pomeranchuk transition has two bands arising from $x$ and $y$ orbitals, the Hartree-Fock phase diagram gives a nematic transition that remains second-order all the way to $T=0$.

As a first step toward a more general formulation of the problem, which highlights the role of purely quantum effects, it is instructive to consider the effective field theory of the pseudo-spin $\vec{\Phi}$. We focus on an insulating tetragonal system with an odd number of sites. If there is a (non-Kramers) orbital doublet that transforms according to a two-dimensional representation ($E_{g}$ or $E_{u}$) of the tetragonal point group, then there is a local pseudo-spin-1/2 associated with each unit cell of the crystal. In this case, all the states – including the ground-state – are at least doubly degenerate. It is illuminating to draw an analogy with the insulating spin-1/2 antiferromagnet: in that case, the two-fold degeneracy of every state is reflected in a topological term in the effective field theory. Correspondingly, the effective field theory for $\vec{\Phi}$ also contains a quantized Berry-phase term:

[TABLE]

where $\tau$ is the imaginary time and ${\bm{A}}({\bm{\Phi}})$ is the Berry connection associated with a monopole in “nematic space” of topological charge $q$. The two-fold degeneracy, absent for the case of singlet orbitals, is encoded in the parity of $q$, which has important consequences for the system. While in the insulating antiferromagnet this term distinguishes integer from half integer spins, in our ferroquadrupolar (nematic) system it encodes the non-trivial commutation relations between $\Phi_{1}$, $\Phi_{2}$, and $\Phi_{3}$. Given that $\Phi_{1}$ and $\Phi_{2}$ break entirely different symmetries than $\Phi_{3}$, in a usual Landau-Ginzburg-Wilson treatment of the problem, explicit reference to these other forms of order would be supressed except near a fine-tuned multi-critical point. In the present problem, these different orders are implicitly intertwined by the commutation relations in Eq. 2, such that anything that increases one component of $\Phi$ increases the amplitude of the quantum fluctuations of the others.

In summary, our work introduces a powerful new tuning parameter, namely shear strain, that enhances quantum nematic fluctuations and which can be varied continuously, without introducing disorder, and without breaking time reversal symmetry. Longitudinal strain has already been extensively used to probe the nematic susceptibility in various systemsChu et al. (2010). Transverse strain provides a new way to tune through the nematic phase diagram, and especially to access the quantum critical regime.

For many materials of current interest, the magnitude of the coupling between nematic and elastic degrees of freedom is not known quantitatively. However, for the specific examples of TmAu2 and TmAg2 discussed here, estimates of the magnetoelastic coupling coefficient combined with the measured elastic modulus $c_{66}$ indicate that shear strains of order 1% would be sufficient to access the quantum phase transition (see Sec. IV D of Supplementary Material). Such strains are certainly not inconceivable, and they might be smaller in other materials depending on microscopic details.

A metal near a time-reversal symmetric Ising-nematic quantum critical point is expected to have an enhanced instability towards superconductivityKim and Kee (2004); Maier and Scalapino (2014); Lederer et al. (2015); Metlitski et al. (2015); Kang and Fernandes (2016). Thus, tuning the ferroquadrupolar transition of TmAu2 and TmAg2 to a quantum critical point using shear strain may drive these compounds into a superconducting phase. In contrast, a $c$-axis magnetic field naturally suppresses superconductivity. Finally, shear strain does not only couple to $\mathbf{Q}=0$ ferroquadrupolar order; $\varepsilon_{xy}$ may also act as a transverse field for antiferroquadrupolar order by promoting quantum fluctuations between components of each local pseudo-spin.

We acknowledge useful conversations with Daniel Agterberg and Andrey Chubukov. AVM, ATH, IRF, and SAK were supported by the Department of Energy, Office of Basic Energy Sciences under contract DE-AC02-76SF00515. EWR was supported by the Gordon and Betty Moore Foundation EPiQS Initiative through grant GBMF4414. RMF was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award number DE-SC0012336. EB was supported by the Israel Science Foundation under Grant No. 1291/12 and by the US-Israel BSF under Grant No. 2014209.

I Role of Berry phase terms in field theoretic formulation

To achieve a more general understanding of the quantum aspect of the problem, it is useful to consider the problem from the more abstract – less microscopic – perspective of an effective field theory. As a first step in that direction, we can integrate out the microscopic degrees of freedom, leaving us with an effective action, $S[{\bm{\Phi}}]$, which determines the quantum statistical mechanics of a three-component collective field, ${\bm{\Phi}}$, such that $\Phi_{3}$ is the nematic order parameter of interest. For simplicity, we will assume we are dealing with in insulator, such that $S[{\bm{\Phi}}]$ can be sufficiently well approximated by a local expression that a conventional Landau-Ginzburg-Wilson (LGW) approach to the problem is reasonable. In contrast, in a metallic system, $S$ might be a complicated non-local functional; a full analysis of this case is beyond the scope of the present work. We will comment on the implications for the metallic case at the end.

For a classical problem, we would simply take $S=S_{LGW}$ where the usual symmetry analysis would be applied to determine the allowed terms in powers of the fields and their derivatives. Indeed, if we are not near a fine-tuned multicritical point, we would simply drop (or more formally, integrate out) the remaining components of ${\bm{\Phi}}$; certainly, if there were no sign of orbital ferromagnetism, we would drop $\Phi_{2}$.

But for the quantum problem, there is another possible Berry-phase term that is allowed:

[TABLE]

where

[TABLE]

Here, ${\bm{A}}({\bm{\Phi}})$ is the Berry connection.

The form of ${\bm{A}}({\bm{\Phi}})$ is constrained by symmetry: $A_{1,3}$ must be odd under time reversal, while $A_{2}$ is even. Under the point group symmetry transformation, the components of $A_{1,2,3}$ transform in the same way as $\Phi_{1,2,3}$, respectively. In addition, ${\bm{A}}({\bm{\Phi}})$ can have singularities at a discrete set of points. Namely, the Berry curvature $\mathcal{B}\equiv{\bm{\nabla}}_{\mathbf{\Phi}}\times{\bm{A}}$ satisfies

[TABLE]

where ${\bm{\Phi}_{i}}$ is a discrete set of points in ${\bm{\Phi}}$ space and $q_{i}$ are a set of integers. (The quantization of $q_{i}$ comes from the requirement that the amplitude of the path integral is single valued.)

A “monopole” of ${\mathcal{B}}$ at ${\bm{\Phi}}_{i}$ must be accompanied by symmetry-related partners under mirror reflections relative to the horizontal and diagonal directions of the square lattice, and by time reversal. Hence, if there is a monopole at $\mathbf{\Phi}_{i}=(\Phi_{1,i},\Phi_{2,i},\Phi_{3,i})$, it must be accompanied by monopoles at $(\pm\Phi_{1,i},\pm\Phi_{2,i},\pm\Phi_{3,i})$. A monopole at the origin, $\mathbf{\Phi}_{i}=0$, does not have symmetry relatives. Note that the components of $\mathcal{B}$ transform in the same way as those of $\mathbf{\Phi}$ under the point group operations and time reversal (see Section II of this Supplement).

Imagine deforming our theory, e.g., by changing the microscopic Hamiltonian, while maintaining all the symmetries of the problem. An isolated monopole of $\mathcal{B}$ at $\mathbf{\Phi}_{0}=0$ with strength $q_{0}$ can change into a monopole of strength $q_{0}-2$, by splitting into three monopoles along any of the three $\Phi$ axes, of strengths $1$, $q_{0}-2$, and $1$, respectively. Importantly, the parity of $q_{0}$ is invariant under such deformations. We can thus classify the allowed Berry phase terms according to the parity of $q_{0}$ of the monopole at the origin.

To understand the two classes of Berry phase terms, we examine their microscopic origin. Consider the low-lying states in the Hilbert space of a single unit cell. These states transform under the point group symmetry transformations as one of the irreducible representations of D4h. If the low-energy states form a two-dimensional representation (Eg or Eu), then these states are doubly degenerate. In this situation, all the low-energy states of an $L\times L$ unit cell system (with odd $L$) form doubly degenerate multiplets. This can be seen by the fact that, acting on a single site, the mirror transformations along the diagonal and the horizontal directions anti-commute with each other: $M_{d}M_{h}=-M_{h}M_{d}$. This property is preserved for any finite-size lattice with an odd number of sites that is symmetric under the point group.

If a symmetry breaking field is applied (either of the $\Phi_{1,2,3}$ character), this degeneracy is split. In the path integral formulation, a monopole at $\mathbf{\Phi}=0$ of strength $q_{0}=1$ encodes the double degeneracy at the origin of $\mathbf{\Phi}$ space. This is similar to the Berry phase term in a coherent-state path integral of a single spin-$\frac{1}{2}$ particle. For example, in the model consisting of a pseudospin arising from an $E_{g}$ doublet, $M_{h}M_{d}=-M_{d}M_{h}$ when acting on a single unit cell (or any finite lattice of size $L\times L$ unit cells with odd $L$), and thus $q_{0}=1$.

Note that this argument relies on having a well-defined transformation law of the low-lying states of each unit cell. It does not apply if the local Hilbert space contains multiplets transforming as different irreducible representations of the point group. In particular, this complicates the analysis if the system is metallic, since adding a single electron to a unit cell can change the transformation properties of the low-lying states.

To understand the effect of the Berry-phase term, let us consider the case in which there is a field conjugate to one component of the three, so we can count on this component always having a substantial magnitude. For convenience, let us consider the problem in the presence of a transverse field $h_{1}>0$, so that in doing the path integral we can always assume that $\Phi_{1}$ is substantial. For an appropriate gauge choice and rescaling of the components of $\mathbf{\Phi}$ such that

[TABLE]

where $\Phi=\sqrt{\Phi_{1}^{2}+\Phi_{2}^{2}+\Phi_{3}^{2}}$. This gauge corresponds to having a quantized flux coming in through an infinitely narrow solenoid along the $-{\bm{e}}_{1}$ axis and then having the magnetic flux spread out from the origin where it terminates. This corresponds to the $\mathcal{B}$ field of a nematic monopole, as described above. Notice that for $\Phi_{2,3}\ll\Phi$

[TABLE]

We are interested in the behavior of $\Phi_{3}$ in the disordered state or in the ordered state not too far from the QCP, so that the condition $|\Phi_{3}|\ll|\Phi_{1}|$ remains valid. We can safely approximate $\Phi_{1}$ by a non-zero constant value, $\bar{\Phi}_{1}$, which is an increasing function of the transverse field, $h_{1}$. Assuming that we are not near a time reversal broken phase, we can integrate out $\Phi_{2}$, which, since its fluctuations are small, can be treated in a Gaussian approximation. The resulting effective action for $\Phi_{3}$ is

[TABLE]

where $\chi_{2}$ is the susceptibility with respect to the other transverse field, $h_{2}$, and $S_{0}$ has all the $\Phi_{3}$ dependent terms from $S_{LGW}$.

In short, there is an additive contribution to the effective mass for $\Phi_{3}$ that is proportional to $q_{0}^{2}$, i.e. it comes from the quantum character of the fluctuations. This term is, in turn, a strongly decreasing function of $|\bar{\Phi}_{1}|$; an increasing transverse field thus causes a decrease in the effective mass, thereby increasing the quantum fluctuations of $\Phi_{3}$.

Before concluding this section, let us comment on the case of a metallic system undergoing a ferroquadrupolar transition. As mentioned above, in this case, the effective action for the order parameter field $\mathbf{\Phi}$ is complicated by the presence of non-local terms. Therefore, strictly speaking, the analysis presented in this section does not apply. However, if our system is composed of localized quadrupoles coupled to itinerant electrons, as in the $4f$ systems TmAg2 and TmAu2 discussed in the main text, we may imagine artificially setting the coupling between the two sets of degrees of freedom to zero, while keeping the exchange coupling between the local quadrupoles finite. Then, applying a diagonal strain tunes through a ferroquadrupolar QCP. Turning the coupling between the itinerant electrons and the local quadrupoles back on, we expect the QCP to survive over a finite range of coupling. Therefore, one may expect that a diagonal strain can tune through a QCP in the metallic case as well, although the properties of the QCP are very different from the case of an insulator.

II Transformation of $\mathbf{\Phi}$ under symmetries

The symmetries of the components of $\mathbf{\Phi}$ under the horizontal and diagonal mirror symmetries, $M_{h}$ and $M_{d}$, and under time reversal, $\mathcal{T}$ are summarized in Table S1. We can construct the combination of the symmetry operations, $\tilde{M}_{h}=M_{h}\mathcal{T}$, $\tilde{M}_{d}=M_{d}\mathcal{T}$, and $\mathcal{T}$, such that under each of these symmetries, one component of $\mathbf{\Phi}$ is odd, while the other two are even. Thus, these symmetries act as reflections in $\mathbf{\Phi}$ space about the $\Phi_{2}-\Phi_{3}$, $\Phi_{3}-\Phi_{1}$, and $\Phi_{1}-\Phi_{2}$ planes, respectively.

The components of $\mathbf{A}(\mathbf{\Phi})$ transform under $M_{d,h}$ in the same way as those of $\mathbf{\Phi}$. Under $\mathcal{T}$, they transform oppositely to the components of $\mathbf{\Phi}$. From this, we can derive the transformation law of the Berry curvature. $\mathcal{B}_{a}=\varepsilon_{abc}\partial A_{b}/\partial\Phi_{c}$ transforms as $A_{a}$ (and $\Phi_{a}$) under the point group operations, and in the same way as $\Phi_{a}$ (i.e., oppositely to $A_{a}$) under time reversal. For example, consider $\mathcal{B}_{3}=\partial A_{1}/\partial\Phi_{2}-\partial A_{2}/\partial\Phi_{1}$;$A_{1}$ and $\Phi_{2}$ are both $\mathcal{T}$ odd, while $A_{2}$ and $\Phi_{1}$ are both even, so $\mathcal{B}_{3}$ is $\mathcal{T}$ even, as is $\Phi_{3}$.

III Mean field theory in a metallic system - dichotomy between 1 and

2 band models

III.1 The models

Let us try to examine the dichotomy between single band and multi-band versions of the transverse field Ising model, at the level of a mean-field field theory. We consider starting from the 1- and 2- band Hamiltonians:

[TABLE]

In the second Hamiltonian, the orbital degrees of freedom have $x$ and $y$ symmetries (e.g. $d_{xz}$ and $d_{yz}$ orbitals). We imagine that the interaction terms $H_{1\text{int}}$ and $H_{2\text{int}}$ serve to drive both systems into a nematic phase, with $x^{2}-y^{2}$ symmetry.

For the single band model, we choose the dispersion to be that of a square lattice, nearest neighbor tight binding model, which in the presence of the $xy$ strain has the form :

[TABLE]

(we will set $t=1$), while for the two band model, symmetry under $x\rightarrow y$ means we can choose:

[TABLE]

where I have defined the $s$-wave and $d$-wave form factors

[TABLE]

Here, the $xy$ strain term has the form $\varepsilon_{xy}\Psi_{\bm{k}}^{{\dagger}}\sigma_{x}\Psi_{\bm{k}}$ where $\Psi_{\bm{k}}^{{\dagger}}=(c_{x\bm{k}}^{{\dagger}}\,\,c_{y\bm{k}}^{{\dagger}})$ is a two component spinor as usual.

III.2 Self consistency equations

To derive the appropriate mean-field self consistency equations, we assume that the dominant interaction in the case of the 1 band model is a forward scattering type interaction:

[TABLE]

where $N$ is the system size. Meanwhile in the two band model, the dominant interaction is

[TABLE]

We will introduce mean field order parameters in both models, of the form

[TABLE]

where expectation values are taken in the quadratic Hamiltonians:

[TABLE]

and

[TABLE]

The self consistency equation for the 1 band model takes the form:

[TABLE]

where $f$ is the Fermi distribution function, while for the two band model it is:

[TABLE]

where

[TABLE]

III.3 Results

We now numerically solve these self consistency equations for a model with a single orbital per site (Fig  S1), and then for a 2 orbital model (Fig.  S2). We have solved the mean field equations on a 2 dimensional lattice with $256\times 256$ points by iterating the self consistency equations until convergence is achieved.

Generically we find that for the single band metal, the nematic phase terminates at a first order quantum phase transition. This is in fact consistent with previous studies Khavkine et al. (2004); Yamase et al. (2005)

On the other hand, continuous zero temperature transitions (i.e. quantum critical points) are possible for the two band model (Fig.  S2), provided the hopping anisotropy parameter $t_{2}$ is not too large.

IV Quadrupolar order in $4f$ intermetallics

As we discussed in the main text, the $4f$ intermetallic compounds TmAg2 and TmAu2 provide a specific realization of Ising nematic order arising from the splitting of a non-Kramers doublet. Here, we show how crystal field effects give rise to an $E_{g}$ doublet as the lowest energy eigenstates, and then discuss how Ising nematic (ferroquadrupolar) order arises using a mean field treatment of the quadrupole-quadrupole interactions.

IV.1 Crystal field effects in a $4f$ system: how a doublet arises from a $J=6$ state.

In the rare earth intermetallic compounds TmAg2 and TmAu2, the Tm ion takes the Tm3+ ($4f^{12}$) state. From Hund’s rules the $4f$ electronic orbitals are filled such that $L=5$ and $S=1$, giving the electronic multiplet a total angular momentum state of $J=6$. Thus for a local Tm site the spherical harmonics $Y^{m}_{6}$ form the natural basis to construct wavefunctions. The surrounding Au ions create a crystalline electric potential which obeys the tetragonal point group symmetry $D_{4h}$, and so the degeneracy of the 13 $Y^{m}_{6}$ states is removed . The nature of the resulting eigenstates (but not their relative energies) can be inferred purely from symmetry arguments. Here we demonstrate this analysis explicitly, following closely the manipulations described in Dresselhaus et al. (2007).

The general goal is to examine how the $\Gamma^{J=6}$ representation of $SO(3)$ is reduced in the tetragonal environment (with point group $D_{4h}$), i.e. we would like to find the coefficients $a_{i}$ in

[TABLE]

where $\Gamma_{i}$ are the irreducible representations of $D_{4h}$. To do this, we must follow these steps:

• •

First, calculate the character $\chi^{(\text{reducible})}(\mathcal{C}_{k})$ of the (reducible) $\Gamma^{J=6}$ representation under the group operations $\mathcal{C}_{k}$ of $D_{4h}$. Using the 13 $Y^{m}_{6}$ spherical harmonics as a basis, the $\chi^{(\text{reducible})}(\mathcal{C}_{k})$ is the trace of the matrix that corresponds to the group element $\mathcal{C}_{k}$.

• •

We then use the orthogonality of characters to obtain the coefficients $a_{j}$:

[TABLE]

where $h=16$ is the order of the group $D_{4h}$, $N_{k}$ is the number of elements in class $k$, and $\chi^{(\Gamma_{n})}(\mathcal{C}_{k})$ is the character of the irreducible representation $k$ in $D_{4h}$.

For example, it can be shown that under a rotation about the $z$ axis by $\alpha$, the spherical harmonic $Y^{m}_{J}$ is transformed as

[TABLE]

so that, under the group element corresponding to four fold rotations ($C_{4}$), the character is

[TABLE]

The characters of the remaining symmetry operations in $D_{4h}$ have been appended to the character table (Table S2), both for arbitrary $J$, and then specifically for $J=6$. Having written the reducible representation $\Gamma^{J=6}$, we can now reduce it on the irreducible representations of $D_{4h}$, using the orthogonality of characters. This gives

[TABLE]

This implies that the CEF Hamiltonian breaks the degeneracy of the 13 $J=6$ states in such a way that there are 3 $E_{g}$ doublets, and 7 singlets. Two of the singlets have wavefunctions which have $A_{1g}$ symmetry, one that has $A_{2g}$ symmetry, two have $B_{1g}$ symmetry, and two have $B_{2g}$ symmetry. We stress that group theory does not determine the relative energies of these representations. For TmAg2, the ground state is determined by the eigenvalues of the CEF Hamiltonian exhibiting the $D_{4h}$ point group symmetry surrounding the Tm sites. This Hamiltonian can be written as a sum of the Stevens operators $\hat{O}^{m}_{l}$ :

[TABLE]

The coefficients $B^{m}_{l}$ have been experimentally inferred through inelastic neutron scattering and magnetic susceptibility to be: $B^{0}_{2}=0.0909$K, $B^{0}_{4}=0.3265$mK, $B^{4}_{4}=-64.5$mK, $B^{0}_{6}=0.0897$mK, $B^{4}_{6}=-1.30$mK Morin and Rouchy (1993). Diagonalizing this CEF Hamiltonian gives the energy spectrum displayed below.

The ground state of TmAg2 is an $E_{g}$ doublet with the closest excited state, an $A_{1g}$ singlet, at 14K. These states determine much of its low temperature behavior.

IV.2 Projecting to the non Kramer’s $E_{g}$ doublet

The basis vectors corresponding to different irreducible representations can be found by use of the projection operator:

[TABLE]

where $\Gamma_{n}$ is the representation of interest, $l_{n}$ is the dimension of the representation, $h$ is the order of the group, $\chi^{(\Gamma_{n})}\left(\mathcal{C}_{k}\right)$ is the character of $\Gamma_{n}$ for group element $\mathcal{C}_{k}$, and $\hat{P}_{\mathcal{C}_{k}}$ is corresponding operator. When the projection operator is applied to an arbitrary function, it projects a linear combination of the basis vectors $\ket{\Gamma_{n}k}$ for that particular irreducible representation.

[TABLE]

For the case of the $E_{g}$ doublet, $l_{n}=2$, $h=16$ for $D_{4h}$, and we can look at the character table above to find the traces of the matrices. Assuming an arbitrary function consisting of all 13 degenerate states:

[TABLE]

Thus the only possible states that can compose the $E_{g}$ doublet states are the odd states. Furthermore, given the behavior of spherical harmonics under the rotation $C_{4}$ one can show that only states that combine angular momenta that differ by $\Delta m=4$ are allowed. To see this, consider the action of $\hat{C}_{4}$ on a state $\ket{\psi}=A\ket{m_{1}}+B\ket{m_{2}}$. We have $\hat{C}_{4}\psi=e^{im_{1}\pi/2}\left(A\ket{m_{1}}+e^{i(m_{2}-m_{1})\pi/2}\ket{m_{2}}\right)$. For $\psi$ to be an eigenstate, we must have $(m_{2}-m_{1})\pi/2=2n\pi$, i.e. $m_{2}-m_{1}=4n$. This, in addition to the fact that the doublet must preserve time reversal symmetry forces us to conclude that the $E_{g}$ representation can be expressed in the spherical harmonic basis as

[TABLE]

Now that we know the form of the $E_{g}$ doublet, let us examine the effects of applying a magnetic field as well as applying strains. In a vector space consisting of basis vectors $\ket{\psi^{\Gamma_{5}}_{1}}$ and $\ket{\psi^{\Gamma_{5}}_{2}}$, a magnetic field $H_{z}$ applied in the $z$-direction measures the $z$ component of angular momentum for each state, and thus is proportional to the operator $J_{z}$:

[TABLE]

where $C=5|e|^{2}+|f|^{2}-3|g|^{2}$. If we were to apply the $B_{1g}$ strain $\varepsilon_{xx}-\varepsilon_{yy}$ the corresponding quadrupole operator $\hat{O}^{2}_{2}=J_{x}^{2}-J_{y}^{2}$ takes the form:

[TABLE]

where $D=21f^{2}+12\sqrt{10}fg+2\sqrt{165}eg$ assuming real coefficients. Finally, the operator for $B_{2g}$ nematic strain $\varepsilon_{xy}$ is the quadrupole operator $\hat{P}_{xy}=\frac{1}{2}(J_{x}J_{y}+J_{y}J_{x})$ which takes the form:

[TABLE]

where $E=\frac{21}{2}f^{2}-6\sqrt{10}fg+\sqrt{165}eg$.

These operators can clearly be identified with the Pauli matrices (Pseudospin operators). Note that a change of basis will rotate these operators into the canonical orbital basis that was used in the main text. We have therefore shown how projecting onto the ground state $E_{g}$ doublet of this system gives rise to the pseudo-spin 1/2 operators.

IV.3 Quadrupole-Quadrupole interactions and the spin 1/2 transverse field Ising model

As discussed in the main text, a generalized version of RKKY interactions Schmitt and Levy (1985); Morin, P. and Schmitt, D. (1988); Aléonard and Morin (1990); Morin and Schmitt (1990) gives rise to Quadrupole-Quadrupole (QQ) interactions. In the presence of a $B_{2g}$ strain, the full Hamiltonian of the $J=6$ system is then:

[TABLE]

Here, $J_{ij}$ is the strength of (attractive) QQ interactions (including strain renormalizations) between sites $i$ and $j$, $\hat{H}_{CEF}$ is the crystal field Hamiltonian outlined in the previous sections, and $B_{B2g}$ is the strength of the coupling to external strains $\varepsilon_{xy}$. Here, $\hat{O}^{2}_{2}$ and $P_{xy}$ are the full $13\times 13$ matrices representing Stevens operators. In Section IV.4 we describe our mean field treatment of this full Hamiltonian.

Before doing so, let us note that the spin 1/2 version of this Hamiltonian may be obtained by projecting these operators onto the ground state $E_{g}$ doublet of $\hat{H}_{CEF}$. Applying this projecting operator, we find the projected version of this Hamiltonian to be

[TABLE]

where $\tilde{J}$ and $\tilde{B}$ are proportional to the original constants $J_{ij}$ and $B_{B2g}$ of Eq. S36. Thus, our mapping to the spin $1/2$ transverse field Ising model is complete.

IV.4 Mean field equations

In the previous section we described how group theoretic arguments guarantee the existence of three $E_{g}$ doublets when a $J=6$ system is embedded in a tetrgonal environment. The specific crystal field Hamiltonian took the form

[TABLE]

where the $\hat{O}^{m}_{l}$ are the Stevens operators and $B^{m}_{l}$ are the corresponding coefficients, which in the case of TmAu2 and TmAg2 result in a ground state $E_{g}$ doublet. This crystal field Hamiltonian ignores quadrupole-Quadrupole interactions and quadrupole-lattice (magneto-elastic) interactions. In this section, we treat these terms within the framework of mean field theory. The corresponding mean field Hamiltonian takes the form

[TABLE]

Here, we have adopted the conventional notation $K_{B1g}=zJ_{ij}$ where $z$ is the number of neighbors for which $J_{ij}$ in Eq. S36 is non-zero. The terms containing strains $\varepsilon_{x^{2}-y^{2}}$ and $\varepsilon_{xy}$ are the bilinear magneto-elastic terms, and the terms containing the quadrupole expectation values are the mean-field version of the quadrupole-quadrupole interaction terms, which arise due to a generalized version of the RKKY interactionSchmitt and Levy (1985); Morin, P. and Schmitt, D. (1988); Aléonard and Morin (1990); Morin and Schmitt (1990). $\hat{O}^{2}_{2}$ is the quadrupole operator, whose expectation value is the $B_{1g}$ order parameter, and is written in terms of angular momentum operators as $J_{x}^{2}-J_{y}^{2}$. Similarly, $P_{xy}$ is the quadrupole operator whose expectation value is the $B_{2g}$ order parameter, and is written as $\frac{1}{2}(J_{x}J_{y}+J_{y}J_{x})$. Assuming the material is free to relax, i.e. the strain reaches the value that minimizes the free energy, the total Hamiltonian can be written as :

[TABLE]

The corresponding mean field equation for the $B_{1g}$ order parameter has the form

[TABLE]

which we iterate numerically until self consistency is achieved. The results of this process for temperatures from 10K to 1K and transverse strains from 0 to 1.5% is shown in the contour plot of the order parameter $\braket{\hat{O}^{2}_{2}}$ vs. temperature and strain (Fig. S4). Figure 3 of the main text shows the lowest 5 eigenenergies of the mean field Hamiltonian as a function of transverse strain $\varepsilon_{xy}$, for a fixed temperature of $T=2$K.

IV.5 Accuracy of transverse field Ising approximation for a $J=6$ system

In Section IV.3 we described how the low energy physics of the full $J=6$ system is that of a (pseudo) spin 1/2 transverse field Ising model. However, the presence of other CEF states in TmAg2 and TmAu2 means that this statement is approximate. This is because the ground state $E_{g}$ doublet will mix with the excited states in both the ordered phase and under large applied strains. An important question which we must therefore address is to what extent can we treat the low energy physics full $J=6$ system as that of the transverse field Ising model? In this section, we show:

• •

The nature of the symmetry breaking (and hence the Ising character of the phase transition) is unaffected by higher energy states that admix with the ground state $E_{g}$ doublet.

• •

The transverse field nature of $B_{2g}$ strain is still appropriate when we consider the full $J=6$ system.

• •

The low energy physics is effectively that of a two level system, (i.e. an $S=1/2$ pseudo-spin), as the character of the ground state wave function is predominantly that of the lowest energy CEF doublet ($\Gamma^{(1)}_{5}$).

Each of these points is discussed below.

IV.5.1 Ising universality in the $J=6$ system

It is important to stress that regardless of whether we treat the essential degrees of freedom as the full $J=6$ system (which is already an approximation, since this describes just the lowest energy state of the spin orbit interaction), or as the projected pseudo-spin 1/2, the phase transition from tetragonal to quadrupolar order remains an Ising phase transition, and so all the critical phenomena are in the Ising universality class. The essential point is that the transition from a system with point group $D_{4h}$, to one with point group $D_{2h}$ is necessarily Ising in character. That is, a finite expectation value for the operator $\langle\hat{O}^{2}_{2}\rangle$ implies that single $Z_{2}$ symmetry is being broken ($C_{4}$ rotations down to $C_{2}$ rotations).

IV.5.2 Transverse field nature of $B_{2g}$ strains

While the projection to the ground state $E_{g}$ doublet makes manifest the relationship between operators $\hat{P}_{xy}$, $\hat{O}^{2}_{2}$, $\hat{H}_{z}$, and the Pauli matrices, we note that the essential physics of the transverse field Ising model is apparent even in the full $J=6$ description. In particular, the quadrupole operator $\hat{O}^{2}_{2}$ and the $B_{1g}$ operator $\hat{P}_{xy}$ do not commute. Including the operator magnetic field operator ($\hat{J}_{z}$) and using the definitions of Stevens operators in terms of angular momentum operators, we have the commutation relations:

[TABLE]

The fact that $\hat{O}^{2}_{2}$ and $\hat{P}_{xy}$ do not commute means that these two operators do not share the same ground state. Increasing the $B_{2g}$ strain enhances quantum fluctuations between these differing ground states, and must eventually lead to a quantum phase transition, i.e. $B_{2g}$ strain acts in the same way as a transverse field in an Ising magnet. This is exactly the same argument that was presented in the main text for an orbital doublet (pseudo-spin 1/2); here we have shown that it still applies even when the full $J=6$ manifold is considered. For completeness, we note that upon projecting these operator equations into the low energy doublet, we recover the typical spin $1/2$ commutation relations ($\mathfrak{su}(2)$ algebra).

IV.5.3 Validity of low energy limit

Our approximation of treating the system as a spin $1/2$ Ising model refers to projecting all operators and eigenstates into the basis of the lowest energy $E_{g}$ doublet (i.e. the $\Gamma^{(1)}_{5}$ doublet) of the CEF Hamiltonian. This is an approximation due to the mixing of higher energy states with this lowest energy $\Gamma^{(1)}_{5}$ doublet under the perturbations of the quadrupole and transverse field operators $\hat{O}^{2}_{2}$ and $\hat{P}_{xy}$ respectively. However, one can show that $\hat{O}^{2}_{2}$ and $\hat{P}_{xy}$, only mix states of $E_{g}$ character, which are well separated in energy from the ground state (see Fig. S5). Therefore, any deviation from spin $S=1/2$ physics is suppressed by this energy gap.

The lower panel of Fig. S5 illustrates this point in a quantitative manner. We plot the overlap of the ground state with its projection in to the lowest energy $\Gamma^{(1)}_{5}$ doublet of the CEF Hamiltonian. In particular, denoting the eigenstates of low energy doublet of the CEF Hamiltonian as $|\psi^{\Gamma_{5}}_{1}\rangle$ and $|\psi^{\Gamma_{5}}_{2}\rangle$, the projection operator into these states is

[TABLE]

We then consider the evolution of the ground state $|\psi_{1}\rangle$ of a mean field Hamiltonian in the full $J=6$ basis,

[TABLE]

Figure S5 (middle panel) shows the evolution of the spectrum of this Hamiltonian as various terms are turned on sequentially (upper panel). The maximum strength of each term is chosen to correspond to the maximum magnitudes of each operator in our self consistent mean field calculation (see Fig. S4). The lower panel shows the magnitude squared of the overlap

[TABLE]

which is essentially the $\Gamma^{(1)}_{5}$ content of the ground state. As is clear from the figure, the “ spin $1/2$ ” character of the ground state is very close to 1 throughout the evolution of the quadrupole and transverse field operators, indicating the validity of treating this system as a pseudo-spin $1/2$.