Evidence for a nematic component to the Hidden Order parameter in URu2Si2 from differential elastoresistance measurements

Measurements of the differential elastoresistance of URu$_2$Si$_2$ reveal that the fluctuations associated with the 17 K Hidden Order phase transition have a nematic component. Approaching the"Hidden Order"phase transition from above, the nematic susceptibility abruptly changes sign, indicating that while the Hidden Order phase has a nematic component, it breaks additional symmetries.

For materials that harbor a continuous phase transition, the susceptibility of the material to various "fields" can be used to understand the nature of the fluctuating order and hence the nature of the ordered state. In the present work we use anisotropic bi-axial strain to probe the nematic susceptibility of URu 2 Si 2 , a heavy fermion material for which the nature of the low temperature "Hidden Order" state has defied comprehensive understanding for over 30 years [1][2][3][4] despite considerable theoretical attention [5][6][7][8][9][10][11][12]. Our measurements reveal that the Hidden Order phase has a nematic component, confirming earlier torque measurements that reveal a broken four-fold symmetry [13] and strongly constraining theoretical models describing the Hidden Order state [14,15].
For a tetragonal material, the conjugate field to electronic nematic order is anisotropic biaxial in-plane strain (either aniso = ( xx − yy ) or aniso = γ xy for spontaneous order in the [100] or [110] directions respectively) [16]. The nematic susceptibility for each of these orientations is then χ N ∝ dψ d aniso , where ψ represents any thermodynamic quantity measuring the induced anisotropy in mutually orthogonal directions ( [17,18]. For a continuous electronic nematic phase transition, χ N diverges towards the phase transition, signaling an instability towards nematic order. Such behavior was recently observed for the representative underdoped Fe-pnictide Ba(Fe 1−x Co x ) 2 As 2 , demonstrating that the structural phase transition that precedes the onset of colinear an-tiferromagnetic order in that material is driven by electronic nematic order [18][19][20].
There are also classes of order for which there is a single transition to an ordered state that has a nematic component, while simultaneously breaking additional symmetries. A simple example would be a unidirectional incommensurate density wave with order parameter ∆ i = ∆(Q i ), for which both translational and tetragonal symmetries are simultaneously broken at the same transition (Q i corresponds to the ordering wavevector of the density wave, and i = x or y). In such a case, anisotropic strain is not conjugate to the order parameter (i.e., the coupling term in the free energy expansion is not bilinear). Nevertheless, the inherent nematicity associated with the unidirectional density wave motivates introduction of a nematic order parameter N which couples linearly to (|∆ x | 2 − |∆ y | 2 ) and to which anisotropic strain is linearly coupled. For such a situation, an analysis based on standard Ginzburg-Landau theory for coupled order parameters (see Supplemental Material [21]) reveals that the nematic susceptibility follows a Curie-Weiss temperature dependence at high temperatures (as the density wave fluctuations orient in the anisotropic strain field) but then develops an abrupt additional divergence at the transition temperature. The additional feature arises from a symmetry-allowed contribution to the free energy of the form where λ is a coupling constant [21,22]. Within the assumption of homogeneity, such a coupling can only occur if the density wave order is described by a multicomponent order parameter [23]. As we describe in detail below, we find that URu 2 Si 2 appears to fall into the latter class of material, having a nematic component but also breaking additional symmetries at the phase transition. Our measurements are based on a novel technique that probes the differential response in the electrical resistivity to anisotropic biaxial strain in mutually perpendicular directions (Figure 1). We measure the induced resistivity anisotropy, N = ρ yy − ρ xx 1 2 (ρ yy + ρ xx ) ∼ ((∆R/R) yy − (∆R/R) xx ), which, by symmetry, is proportional to ψ in the limit of asymptotically small values. The induced anisotropy is characterized by specific terms in the associated elastoresistivity tensor m ij that relate changes in the resistivity to strains experienced by the material [20,21,24]. In the regime of infinitesimal strains (linear response), the elastoresistivity coefficients (m 11 − m 12 ) and 2m 66 are linearly proportional to the bare (unrenormalized) nematic susceptibility, χ N [100] and χ N [110] , for strains in the [100] and [110] directions respectively [20]. The proportionality constant relating the resistivity anisotropy to the nematic order parameter is governed by microscopic physics, but does not contain any singular contributions. Consequently, any divergence of the induced resistivity anisotropy is directly related to divergence of the nematic susceptibility. Since this is ultimately a derivative of a transport measurement, the technique is especially sensitive to static or fluctuating order which affects the Fermi surface.
Of particular relevance to the current experiments, recent magnetic torque measurements (measured with a small applied magnetic field in the ab-plane) indicated the onset of a two-fold anisotropy in the Hidden Order phase of URu 2 Si 2 [13]. This result is apparently supported by the observation of a weak orthorhombicity that onsets at the same temperature in high quality crystals [15]. These subtle effects indicate that the Hidden Order phase has a nematic component, but the degree of anisotropy that is observed depends heavily on the crystal size [13] (for torque) and possibly also on the crystal quality [15] (for X-ray diffraction), leading to some contention as to whether or not these are intrinsic effects. By probing the temperature dependence of the nematic susceptibility in the tetragonal state, our current experiments confirm that the Hidden Order phase does indeed have a nematic component. In addition, our observation of a sharp downward deviation from mean field Curie-Weiss behavior close to T HO implies that the Hidden Order phase is not a pure nematic and must break additional symmetries.
Anisotropic strain is achieved by gluing thin crystals of URu 2 Si 2 to the surface of a PZT piezoelectric stack. The  elastoresistivity coefficients m 11 − m 12 and 2m 66 are determined from the difference of (∆R/R) yy and (∆R/R) xx for a given set of strains for [100] and [110] oriented crystals respectively [20,21]. To first order this is also equal to the induced resistivity anisotropy N , expressed below in terms of the anisotropic biaxial strain ( yy − xx ): where x and y refer to the axes along which the piezo stack induces biaxial strain, indicated in Figure 1(a) and following the description given in previous work [19]. For each temperature, five voltage cycles were performed. The strain was measured using mutually perpendicular strain gauges glued to the back surface of the piezo stack. Representative measurements of ( ∆R R ) yy for five different temperatures are shown in Figure 1 (b) and (c) for [110] and [100] oriented crystals respectively, revealing a linear response. This procedure was performed for both the ( ∆R R ) xx and ( ∆R R ) yy orientations by removing the crystal from the piezo with acetone between each measurement. The slopes ∆R/R versus strain for each orientation were found using a linear fit and the difference was taken. The resulting elastoresistivity coefficients 2m 66 and m 11 −m 12 are plotted in Figure 2 as a function of temperature.
The data shown in Figures 1 and 2 reveal a striking anisotropy in the elastoresistivity coefficients, with 2m 66 values considerably larger, and also more strongly temperature-dependent, than m 11 − m 12 . For comparison, the elastoresistivity coefficients of simple metals are small (of order one) and essentially isotropic [20]. A sharp downward anomaly at T HO = 17.15 K marks the phase transition to the Hidden Order state. We return shortly to the detailed behavior in the vicinity of T HO , first focusing on the temperature dependence of the differential elastoresistivity further above the transition. Before doing so, we note that for temperatures below approximately 15 K, the resistivity drops rapidly with decreasing temperature due to the high quality of the samples used [21], leading to a decrease in the signalto-noise ratio of the measured elastoresistance (evident in the increased scatter in the 2m 66 and m 11 − m 12 coefficients shown in Figure 2) and obscuring detailed analysis of the elastoresistance in the Hidden Order phase.
Recalling that 2m 66 is proportional to the nematic susceptibility in the [110] direction (χ N [110] ), the large values of 2m 66 (red data points in Figure 2) and the strong temperature dependence both indicate a diverging nematic susceptibility, and hence motivate fitting the data to the Curie-Weiss form (2m 66 = C/[(T − θ)] + 2m 0 66 ). Since the elastoresistance only develops significant values below approximately 60 K, the range of temperatures over which this fit can be applied is necessarily small. Nevertheless, if we attempt fits over different temperature ranges, the best fit, defined by a reduced chi-squared statistic closest to 1, occurs for a range from approximately 30 K up to the maximum of 60 K, with a Weiss temperature θ = 15.2 ± 0.3 K (see supporting material for a discussion of the fitting procedures) [21]. Extrapolation of the best fit function to lower temperatures (shown as a green line in Figure 2) reveals significant deviation from Curie-Weiss-like behavior for temperatures below approximately 30 K. The onset of this deviation at 30 K is coincident with a suppression in (T 1 T ) −1 seen in recent NMR experiments [25], suggestive of a common origin, possibly associated with the progressive development of strongly fluctuating order. Perhaps coincidentally, a finite polar Kerr effect signal also develops around this temperature [26].
In contrast, m 11 −m 12 , which is proportional to χ N [100] , is small and exhibits only a weak temperature dependence (cyan data points in Figure 2). Given the large values of 2m 66 , it is likely that this weak temperature dependence derives from slight sample misalignment (i.e., a small amount of contamination from 2m 66 in the [100] The green curve is a fit of 2m66 to a Curie-Weiss temperature dependence 2m66 = C T −θ + 2m 0 66 from 60 K to 30 K, as described in the main text. The fit is extrapolated below 30 K to emphasize deviations from Curie-Weiss behavior below this temperature. The values for the fitting parameters are; C = 375 K, θ = 15.2 K, and 2m 0 66 = −7.9. The inset shows (2m66 − 2m 0 66 ) −1 versus temperature, with 2m 0 66 determined from the green fit in the main figure. orientation). For the remaining discussion we focus solely on the 2m 66 data.
The most striking aspect of the temperature dependence of the elastoresistivity coefficients is the sharp downward anomaly observed at T HO . The 2m 66 data are replotted in Figure 3(a) on an expanded scale for the temperature window from 15 K to 19 K, spanning T HO . The anomaly exactly aligns with the sharp peak in the temperature derivative of the resistivity dρ/dT ( Figure  3(b,c)), indicating that the effect is associated with the thermodynamic phase transition at T HO . For comparison, heat capacity measurements were also performed for one of the crystals cleaved from the same rod that was used for the elastoresistance measurements. As has been previously established for URu 2 Si 2 and as anticipated for a metallic system [27], dρ dT follows the heat capacity close to T HO , with slight differences in the peak position attributed to small differences in thermometry and residual resistivity values between the two samples.
Elastoresistance measurements for temperatures greater than ∼150 mK away from T HO reveal a linear response to strain, and the eleastoresistivity coefficients are accurately obtained from a linear fit over the FIG. 3: (Color online) Temperature dependence of (a) the 2m66 elastoresistivity coefficient, (b) the resistivity, ρ, (c) the temperature derivative of the resistivity, dρ dT , and (d) the electronic specific heat [28] in the immediate vicinity of THO. THO obtained from the transport measurements is shown by a black vertical dashed line. Differences between THO obtained from transport and electronic specific heat data (gray dashed lines) are attributed to differences in thermometry and/or RRR of the samples used for the two measurements. entire range of strain. However, for data points within ∼150 mK of T HO , the elastoresistance becomes highly nonlinear [21]. The origin of this effect is at present unclear, but there are two natural candidates. On the one hand, the effect could arise from strain-induced changes to T HO [4]. Within a simple mean field model of a two-component order parameter, such as would give rise to a Curie-Weiss-like nematic susceptibility, anisotropic strain results in a linear increase/decrease in T HO for the two distinct components of the order parameter [21]. Estimates of ρ(T ) for fixed strain obtained from the elastoresistance measurements indicate a shift in T HO of order ± ∼100 mK for strains of order ±10 −3 , consistent with such an interpretation. Equally, as the phase transition is approached and the free energy landscape becomes progressively flatter and more strongly determined by the fluctuating order, the range of strains applied will at some point exceed the regime of linear response. In both cases, the susceptibility can still be measured by the slope of the elastoresistance at zero strain, though in practice signal-to-noise limits the range over which fits can be made. Linear fits to the elastoresistance over the full range of strain in this regime underestimate the slope at zero strain, rounding the sharp anomaly seen in Figure 3 (a) and prohibiting a critical scaling analysis close to T HO [21].
The temperature dependence of the 2m 66 elastoresistivity coefficient of URu 2 Si 2 , which is proportional to χ N [110] , establishes that the fluctuating order associated with the Hidden Order phase has a nematic component. The progressive alignment of the fluctuating order in the anisotropic strain field yields the observed Curie-Weiss-like temperature dependence. Significantly, since 2m 66 does not continue to diverge monotonically towards the eventual phase transition, but rather diverges in the opposite direction as the phase transition is finally reached, it is clear that the Hidden Order phase is not a pure nematic and must break additional symmetries, as The same coupling to the crystal lattice that yields the resistive response to strain must also cause an orthorhombic distortion below T HO , consistent with reports of high resolution X-ray diffraction measurements [15] and also with the observation of two-fold anisotropy in the magnetic torque [13]. Whether or not this anisotropy is observed in macroscopic crystals is then a matter of the resolution of the specific measurement and the relative population of the two orthogonal orthorhombic domain orientations (which will be influenced by both crystal size and quality, via effects associated with pinning of domain walls). Both aspects illustrate the distinct advantage of probing the nematic susceptibility in the tetragonal phase (i.e., for T > T HO ), for which there are no domains (one is simply probing the susceptibility of the symmetric phase) and for which the resistivity anisotropy is extremely sensitive to subtle perturbations of the Fermi surface.
Despite these advances, important questions remain. In particular, since nematic fluctuations in the B 2g channel couple linearly to the shear modulus, one anticipates a softening of the C 66 elastic constant approaching the Hidden Order phase transition. However, earlier measurements of the elastic moduli [29,30] do not show such a softening. At least in principle, weak coupling between the nematic order parameter and the crystal lattice (consistent with the small magnitude of the orthorhombic distortion observed in X-ray diffraction experiments [15]) would limit any appreciable renormalization of C 66 to small values of the reduced temperature close to the phase transition [31]. Nevertheless, the present results, in conjunction with other measurements that reveal a twofold anisotropy in the Hidden Order phase [13], clearly motivate a careful reinvestigation of the elastic properties of URu 2 Si 2 .
In conclusion, we have shown that the 2m 66 elastoresistivity coefficient of URu 2 Si 2 follows a Curie-Weiss-like temperature dependence from 60 K (the highest temperature at which a finite elastoresistance can be detected by the current experimental method) down to 30 K, indicating that the flucuating order has a nematic character. As the Hidden Order phase transition is approached, the differential elastoresistance progressively deviates from the high temperature mean field behavior, and then, very close to the phase transition, develops a sharp downward divergence ( Figure 2). Since the differential elastoresistance for this configuration is proportional to the nematic susceptibility χ N [110] , this behavior directly reveals that anisotropic strain is not a conjugate field for the order parameter, and hence the phase transition must break additional symmetries beyond four-fold rotational symmetry. Taken together these data demonstrate that the Hidden Order phase has a multicomponent order parameter [23]. These results put tight theoretical constraints on any proposed theory of the Hidden Order phase transition; moreover, our measurements reveal the utility of differential elastoresistance measurements in determining the nature of subtle electronic phase transitions. The

Electronic Nematic Transition
To model a second order electronic nematic phase transition with coupling to the tetragonal lattice, one can Taylor expand the free energy to quartic order in terms of the nematic order parameter N and the structural order parameter ε: In this expansion, the bilinear term λN ε is the lowest order coupling allowed by symmetry, T N is the bare nematic temperature which would characterize this transition in the absence of coupling to the lattice, and a 0 , b, c, d, λ are positive coefficients. Due to the bilinear coupling, a nonzero N induces a finite ε, which compels the C 4 → C 2 point group symmetry breaking of the lattice. In order to compute the nematic susceptibility from this free energy, we find the values of N which minimize f by differentiating with respect to N : This defines a cubic equation for N which can be solved in terms of ε; however, we are interested in the linear response ∂N ∂ε , which we can find by implicit differentiation: The nematic susceptbility is given as χ N ≡ lim ε→0 ∂N ∂ε , and since N = 0 for T > T N , above the transition we have the Curie-Weiss form:

Nematic × (?) Symmetry
A nematic response is also possible if the Hidden Order state breaks four-fold spatial rotation symmetry but does not couple to strain as a bilinear term in the free energy (unlike the purely electronic nematic transition in Section A 1). In this scenario, the Hidden Order state can be described by a two-component vector order parameter ∆, the nematic response can be modeled in mean field theory by introducing a separate nematic order parameter N , and the two are coupled in the form of Equation (1) of the main text. In this section, we show that within such a theory, a sharp jump in the nematic susceptibility occurs at the Hidden Order transition and is proportional to the singular heat capacity shift of Hidden Order. This correction to the nematic susceptibility can be seen from both mean field considerations-where the jump occurs infinitesimally below the transition-and upon including Gaussian fluctuations (Sections A 2 a and A 2 b, respectively).
A natural question here concerns the origin of N . While mean field theory simply treats ∆ and N as coupled order parameters, and so is agnostic about their origin, it is possible that N is an avatar of the difference in fluctuations of the two components of ∆. This has been argued to be the case in the iron pnictide superconductors [32], and in Section A 2 b we outline how the nematic order parameter arises in such a scenario.
Another possibility is that there may be some intrinsic nematic ordering tendency which is preceded by the Hidden Order transition, in which case it must be included in the free energy from the start. As we show below, the free energies for these two possibilities are identical, and so the sharp discontinuity in the susceptibility occurs in either case, provided the Hidden Order transition occurs first.

a. Mean Field Intuition
We consider here the mean field theory of a twocomponent vector order parameter ∆ = (∆ x , ∆ y ) coupled to a nematic order parameter N . The minimal symmetry-allowed mean field free energy expansion takes the form In order that the free energy be bounded from below, we take the coefficients a 0 , b, r 0 , and u to be positive, but the parameters α, γ, and λ are of indeterminate sign. T HO and T N are the bare Hidden Order and nematic transition temperatures (respectively) which would characterize those phase transitions in the absence of both the external strain field h, and the coupling term −λN |∆ x | 2 − |∆ y | 2 , and we assume that T HO > T N so that the nematic order onsets parasitically at the Hidden Order transition. The coefficient w determines whether the Z 2 symmetry is broken at the transition-for w > 0 only one component orders below T HO , while for w < 0 both components become non-zero and there is no Z 2 symmetry breaking. Finally, note that we will formally take the external strain field h to zero in computing the nematic susceptibility. We will henceforth assume that w > 0 so that the system spontaneously breaks Z 2 symmetry at the transition, and we also assume that ∆ x is the component which orders (i.e. ∆ y = 0 throughout). The free energy in (A5) then reduces to DefiningT HO (h, N ) ≡ T HO + 2γh r0 + 2λN r0 , we can recast (A6) in the more suggestive form This form emphasizes that the effect of the γh|∆ x | 2 and λN |∆ x | 2 terms is essentially to shift the bare Hidden Order temperature T HO . The phase transition occurs when ∆ x acquires a nonzero expectation value, which is found in mean field theory by computing the stationary points of f : Plugging these back into (A7), we find that αhN is the free energy for a second order phase transition characterized by a single nematic order parameter. This would be the appropriate free energy to model a nematic transition in a liquid crystal, for instance, as it does not include coupling to a crystal lattice.
N will assume values that minimize the free energy, which we find by differentiation: r0 . (A10) defines a cubic equation in N that we could solve in terms of h; treating N as an implicit function of h, we can take the derivative of (A10) with respect to h and solve for the linear response ∂N ∂h : where we have used that ∂T HO ∂h = 2γ r0 . Above and at T = T HO , the system has not transitioned to the Hidden Order state, and hence ∆ x = N = 0; furthermore, the nematic susceptibility is defined in the limit of vanishing h. In these limits,T HO (h, N ) = T HO and the nematic susceptibility is given by (A12) Note that T − HO indicates a temperature infinitesimally below T HO , and so there is an abrupt change in the nematic susceptibility across the transition temperature. Above T HO , the nematic susceptbility is exactly as for the case of an electronic nematic transition, but there are additional contributions below the ordering temperature. This abrupt change occurs because of the change in the free energy at T HO , which in mean field manifests as a discontinuous second derivative of the energy (or, equivalently, a discontinuity in the heat capacity).
We can relate the parameter u to the magnitude of the heat capacity singularity at T HO by focusing on the Hidden Order contribution to the free energy on cooling through the transition and invoking the thermodynamic (A13) Using this result, we find that the nematic susceptibility in terms of the heat capacity anomaly at T HO is Thus, across the Hidden Order transition (where ∆ x orders), the nematic susceptibility receives a correction that is proportional to the specific heat capacity jump of the Hidden Order. In mean field theory, this only arises just below the Hidden Order transition temperature, but fluctuations above the transition temperature have the same effect.

b. Functional Form
In this section, we illustrate how a nematic order parameter can arise from fluctuations of an underlying order with Z 2 symmetry. Here we follow closely the procedure outlined in [32], where we begin with a twocomponent vector order parameter with the free energy density with r ∝ (T − T c )/T c , v > g, and ∂ = ∇ the spatial gradient operator. This is a simple re-writing of (A5) with v = (2w + u)/4 and g = (2w − u)/4. The nematic order parameter in this situation can in fact correspond to the difference in fluctuations of the two components of ∆. More formally, we can study such a nematic state by employing a Hubbard-Stratanovich transformation in which we decouple both fourth order terms above. We begin with the original partition function, which is a functional integral over the order parameter fields ∆ x , ∆ y : In writing S eff in this way, we have assumed that ∆ x and ∆ y are frequency and momentum dependent.
To study the nematic phase, we introduce auxiliary scalar fields N for |∆ x | 2 − |∆ y | 2 and ζ for |∆ x | 2 + |∆ y | 2 and perform a Hubbard-Stratonovich transformation to decouple the quartic order (in ∆ i ) terms in the free energy. The Hubbard-Stratonovich transformation is a generalization of the identity This represents a major simplification because we have reduced the theory to one which includes only Gaussian integrals, but this comes at the cost of two additional integrals over the extra auxiliary order parameter fields. We would like to compute the nematic susceptibility above the temperature at which both ∆ and N order (i.e., ∆ x = ∆ y = N = 0), and so, in addition to transforming the free energy in (A5), we also add in the effect of a symmetry breaking strain field h. Proceeding to Gaussian level in the order parameter fields and the fluctuations, this corresponds to the theory where r ∝ (T − T c )/T c and a = 1 g ∝ (T − T N )/T N (Note: we assume that T c > T N so that nematic order develops parasitically when ∆ orders). The coefficient α reflects the strength and sign of the coupling of N to the external field, relative to that of ∆.
The field ζ (with quadratic coefficient b = −1/v) essentially renormalizes the coefficient r in (A19), and so we neglect its fluctuations and analyze the theory by integrating out ∆ x and ∆ y : with |∆ x | 2 given by and |∆ y | 2 given by Thus, taking the saddle point of the partition function in (A20), we have for small values of the applied field This in turn implies that the nematic susceptibility above T c is given by What (A24) shows is that in the disordered phase of ∆, the nematic susceptibility above the phase transition is the sum of a Curie-Weiss term (the same as for the electronic nematic case) and an additional contribution which, as will be shown below, has a specific heat-like dependence on temperature.

c. Specific Heat-like Singularity
To see why ∂G(r) ∂r behaves like the specific heat singularity due to the phase transition, we will express the partition function in terms of the mean field order parameters (including spatial fluctuations), derive the free energy f , and then compute the heat capacity singularity due to the phase transition as c sing Although we focus only on the heat capacity contribution due to the phase transition c sing , there are many other degrees of freedom that contribute to the total heat capacity; however, these other degrees of freedom do not change abruptly at T c .
Following this procedure, we consider the full partition function (with fluctuations) up to quartic order in the ∆ fields: We begin by expanding ∆ about its mean field value, i.e., by letting where the mean field values of the components of the order parameter are given by and |∆ y | = 0 for all temperatures (since the state is unidirectional). Expanding the terms in (A25) to quadratic order in the fluctuations ψ x and ψ y , Observing that the linear terms in ψ x cancel, the new theory including fluctuations is given by where we have introduced the length scales Using the formula ln det G −1 = Tr ln G −1 , we can integrate out the fluctuation fields ψ x , ψ y in (A29) to obtain the free energy density: For T > T c , the two terms in the integral are the same. We can now compute the heat capacity singularity c sing = ∂ Applying this result to (A31) and substituting the r dependence of |∆ x | , we find that for T > T c c sing = ∂ ∂r where G(r) is the correlation function that we wrote down in the previous section. Hence, ∂G/∂r is equal to the fluctuational contribution to the heat capacity and proportional to the nematic susceptibility correction in the disordered phase.

Appendix B: Elastoresistivity Measurements
For a tetragonal material, the elastoresistivity tensor (which relates the normalized resistivity change to the strain in different directions) is defined by In writing this equation, we have employed the Voigt notation, which allows us to express the rank 4 elastoresistivity tensor as a matrix and the rank 2 strain and resistivity change tensors as 6-component vectors. Crucially for our experiment, we define our Cartesian axes relative to the piezo stack itself, with the x-(y-)direction parallel to the short (long) axis of the piezo. The elastoresistivity coefficients are determined by measuring the resistance change in response to an applied strain. A sample is mounted to a piezo stack with a thin layer of Devcon 5-min epoxy. Applying a positive voltage to the piezo stack causes it to expand in the longitudinal direction (long axis of the piezo) and contract in the transverse direction (short axis of the piezo). The ratio of the longitudinal expansion, which applies positive yy strain on the sample, to the transverse contraction which applies − xx strain to the sample, is the Poisson ratio of the piezo stack, ν p = − xx yy . Note that this convention for the piezo stack Poisson ratio is the inverse of that defined in our previous work [20]. For thin samples mounted on the piezo stack the in-plane sample strain is equal to the piezo strain. The out-ofplane sample strain is free to expand or contract and is controlled by the Poisson ratio of the material ν s . The measured change in resistance to applied strain from the piezo stack is therefore an admixture of several elastoresistivity coefficients. To disentangle particular elastoresistivity coefficients, measurements are made in the longitudinal ( ∆R R ) yy and transverse ( ∆R R ) xx configurations along high symmetry directions. Specifically, measurements are made for crystals oriented such that the current flows along the [100] direction and the [110] directions [20]. Differential elastoresistance measurements for these two configurations yields: = yy (1 + ν p )2m 66 (B2) The above two equations are used to extract m 11 −m 12 and 2m 66 , which correspond to the B 1g and B 2g symmetry channels respectively. As explained in the main text, these coefficients are proportional to the components of the nematic susceptibility tensor χ N [100] and χ N [110] , relevant to nematic order in the [100] and [110] directions, respectively.

Crystal Growth
Single crystals of URu 2 Si 2 were grown by the Czochralski technique and then electro-refined to improve the purity. The crystals were oriented by Laue diffraction in back-scattering geometry and cut with a wire saw. In some cases, the crystals were cleaved with a razor blade and bar shaped samples were extracted. After elastoresitivity measurements the in-plane orientations of the URu 2 Si 2 samples were confirmed using a commercial four-circle single crystal diffractometer in an Eulerian cradle geometry, with CuK α1 radiation, of wavelength λ = 1.54Å. Fine alignment of the diffraction geometry was done using the (004) out-of-plane orientation peak, which is the one with the strongest intensity for this family of planes. Subsequently, the [112] peaks were identified by moving to the corresponding ω and ξ (azimuthal) positions, and performing φ scans (rotation about the out-of-plane axis).  Table I provides details for the three samples measured, where RRR = R(300K) R(4.2K) . The resistivity as a function of temperature for the three samples is plotted in Figure 4.  [110] were used for the elastoresistivity measurements in the main text. Sample 2 [100] was also measured and showed the same temperature dependence in m11 − m12 as Sample 1.

Proximity to Strain-induced Transition
It has previously been established that hydrostatic pressure in excess of approximately 5 kBar stabilizes the competing local moment antiferromagnetic (LMA) groundstate of URu 2 Si 2 (for a relatively recent review of key experimental results see [4] and references therein). In addition, recent theoretical treatment employing a Ginzburg-Landau model based on the lowest-lying crystal field states predicts a similar stabilization of the LMA state due to uniaxial stress [10]. In this theory, the HOto-LMA transition is predicted to occur for uniaxial stress σ ∼ 0.6 GPa directed along either the [100] or [110] directions. Although this effect has yet to be established experimentally, it is important to compare the strains experienced by the crystals in our experiment, with those that would be anticipated for such large stresses. Since the elastic moduli of URu 2 Si 2 are known [29], we can readily calculate the strain that would result from uniaxial stress of this magnitude. In the following analysis we show that our experiment employs strains that are considerably smaller than those corresponding to the critical stress boundaries.
For tetragonal symmetry, strain and stress are related by the tensor equation In writing this tensor equation we have employed the Voigt notation convention, whereby xx → 1, yy → 2, zz → 3, yz → 4, zx → 5, and xy → 6. Answering our main question requires inverting the elastic moduli tensor, C ij , which for tetragonal symmetry gives where α ≡ C 33 (C 11 + C 12 ) − 2C 2 13 . Explicit inversion of this matrix was done in Mathematica, but the exact form is also given in [33].
In this form, we can now consider the cases of a uniaxial stress of magnitude σ applied in the [100] and [110] directions. Starting with the simpler case, a uniaxial stress of magnitude σ applied in the [100] direction amounts to the case where σ xx = σ and all other σ ij = 0, or equivalently that For the more complicated case of a uniaxial stress applied in the [110] direction, σ xx = σ yy = σ xy = σ 2 and all other σ ij = 0. This result is taken from [34], but it can be seen from two different vantage points: 1) by resolving the applied stress into its x-and y-components while accounting for the enlarged cross sectional area of the plane normal to the [110] direction, or 2) by applying a 45 • orthogonal coordinate transformation to the stress tensor. Using this result, we obtain that At this stage, we can now plug in numbers to estimate the strain required to move across the HO-LMA phase boundary. For (C 11 , C 12 , C 13 , C 33 , C 44 , C 66 ) = (255, 48, 86, 313, 133, 188) GPa (taken from [29], which are extrapolated to T = 0 K), a stress of magnitude σ = 0.6 GPa applied in the [100] direction induces a strain of while a stress of magnitude σ = 0.6 GPa applied in the [110] direction induces a strain of This analysis reveals that the strain in the [100] ([110]) direction associated with a uniaxial stress of 0.6 GPa in the [100] ([110]) direction is of order 10 −3 . In contrast, the largest strain experienced by the samples in our experiments corresponding to a full voltage sweep of the piezo stack, is of order ∆l l ∼ 10 −4 . Some of the transverse strains associated with the uniaxial stress are nevertheless smaller than 10 −3 , so the analysis described in the main text to extract the elastoresistivity coefficients was repeated for a smaller range of strains. The temperature dependence of the elastoresistivity coefficients was the same for both small (∼ 10 −5 ) and large (∼ 10 −4 ) strain regimes, consistent with the linear response observed for the resistance as a function of strain (see for example Figure 1 of the main text). The only difference in the temperature dependence of the elastoresistivity coefficients for the two strain ranges was that the signal-tonoise became worse for the smaller strain range. These effects are discussed in more detail in Section C 4.
Another source of strain comes from mounting the sample to the piezo stack, which creates a "built-in" strain from the epoxy [35]. To test for the effects of built-in strain, the resistivity as a function of temperature was measured with the sample free-standing and when glued to the piezo stack ( Figure 5). The two transport curves are essentially identical, showing no discernible shift in the Hidden Order transition temperature. Therefore, we conclude that any built-in strain from the epoxy does not contribute to the resistive response, i.e., R(V )−R0 , where R 0 is the free standing resistance and R(V=0) is the resistance of the sample mounted on the piezo at zero applied voltage.

Curve Fitting
Neglecting the sharp anomaly at T HO , the elastoresistivity coefficient 2m 66 exhibits a monotonic increase with decreasing temperature from 60 K (the highest temperature at which there is a measureable elastoresistance) down to just above T HO . In the following two sections, we briefly describe the procedures used to fit the data first to Curie-Weiss temperature-dependence (motivated by the Ginzburg-Landau models described in Section I), and then to a more general power law.

a. Curve Fitting -Curie-Weiss
Motivated by the Ginzburg-Landau theory described in Section I and in the main text, the 2m 66 elastoresistance data were fit to a Curie-Weiss temperature dependence 2m 66 = C T −θ +2m 0 66 . When the data are fit over the full temperature range from 60 K down to 18 K (i.e. to just above the anomaly centered at 17.2 K), a clear systematic deviation can be seen by inspection (red curve of Figure 6), which is confirmed by the reduced χ 2 value (inset Figure 6). Similar fits were performed between 60 K and different lower-bound cut-off temperatures to find a reduced χ 2 which was closest to 1, which is how we define the best fit (inset to Figure 6). The best fit determined by this method is found to be for the range between 60 K and 30 K (green curve, Figure 6), with fit parameters given in the main text.

b. Curve Fitting -Power Law
In addition to fitting the data to a Curie-Weiss temperature dependence, as described above, a more general fit to an arbitrary power law was also attempted. However, with less than a decade in reduced temperature, the data are not well-suited to such a scaling analysis, and the results are not conclusive. We include this analysis here for completeness. 18  We start by assuming that the elastoresistance is given by a single power law: where 2m 0 66 , C, θ and α are treated as fit parameters. Importantly, α is a nonlinear fit parameter, which greatly complicates the problem. One can linearize α (the parameter we are most interested in) by subtracting off 2m 0 66 , taking the natural log of both sides, and then taking a temperature derivative. If we could precisely estimate 2m 0 66 and sensibly smooth the data in taking the temperature derivative, then we would be left with the simple task of fitting The problem with the procedure leading to (C2) is that if we misestimate 2m 0 66 by some amount δ, then we no longer have a linearized fit form for α: Furthermore, there is no sensible way to linearize (C3) via a Taylor expansion since we expect the temperature dependent term to dominate near criticality but to progressively approach zero at higher temperatures. These issues and others have been addressed in the literature [36]. Instead, we chose a graphical method. Assuming the functional form (C1) and fixing a particular α, one can unbiasedly estimate the other three linear parameters (and in particular the temperature independent term 2m 0 66 ) through a normal fitting routine. Since (C1) is equivalent to it is straightforward to compute 1 (2m66−2m 0 66 ) 1/α and to graph the result versus temperature. Doing this for many values of α yields the critical exponent, as the right choice will result in a quantity linearly proportional to temperature.
Importantly, this method avoids the pitfalls of nonlinear parameter estimation and numerical differentiation; however, it is still quite susceptible to errors in 2m 0 66 , as any error δ results in the highly nonlinear expression The nonlinearity from misestimating 2m 0 66 expressed in (C5) is worth keeping in mind since 1 (2m66−2m 0 66 ) 1/α is not linear in temperature for any α spanning the range 0.5 ≤ α ≤ 3 (Figure 7). The fits are done over the temperature range 20 K ≤ T ≤ 60 K, where we are limited on the low end by fluctuations above T HO and on the high end by our ability to resolve the elastoresistive response. Note that the estimate of 2m 0 66 is itself dependent on the particular α chosen, and so varies between curves. Higher values of α yield highly unphysical estimates of θ and are hence neglected.
From the temperature dependence of the 2m 66 data (Figure 7), we do not unequivocally observe a power law divergence with an extractable critical exponent. The task of fitting critical exponents is potentially complicated by our limited range in temperature (less than a decade), and this might be why it is difficult to distinguish between different α's. The results are also potentially skewed by incorrect estimations in the temperature independent parameter 2m 0 66 . Consequently, the most physically meaningful fit to the data is provided by the Curie-Weiss model described in the main text (i.e. fixing the exponent α = 1).

Nonlinear Elastoresistance near THO
For temperatures more than ∼150 mK away from T HO , the elastoresistance is linear with strain over the entire range of strain experienced by the material in our experiments (which, in the vicinity of T HO is approximately −0.5 × 10 −4 < yy < +1 × 10 −4 ); however, very close to T HO = 17.15 K, the elastoresistance develops a nonlinear response to strain. This effect can be seen in Figure  8(a), which shows successive measurements of ( ∆R R ) yy for a [110] oriented crystal as a function of yy for temperatures spanning T HO . In the nonlinear regime, the magnitude of the elastoresistivity coefficents are slightly underestimated if the data are fit to a linear function over the entire data range. This effect can be seen in Figure 8(b), which shows ( ∆R R ) yy / yy evaluated for different strain ranges (see Figure 8 caption for further detail). It is important to note, however, that the sign change of the slope of the elastoresistance is a robust feature, observed for large and small strains and for temperatures that extend beyond the nonlinear regime.
The physical origin of the nonlinear response very close to T HO is unclear. The effect could be associated with critical behavior, or possibly with strain-induced changes in T HO . To investigate the latter idea, we estimate the rate at which T HO is affected by anisotropic strain. We do this by extracting the temperature dependence of the resistivity for fixed values of the strain. This is obtained by performing a linear fit to the elastoresistance over different strain ranges (see Figure 8 caption for further detail) and using the fit parameter to calculate the resistance for a specific strain (we choose ±5 × 10 −4 , but the exact value is immaterial). The resulting curves are shown in Figure 8(c), together with the temperature de-pendence of the resistivity for the same unstrained crystal (mounted on the piezo, with 0 V applied to the piezoelectric stack). For all temperatures more than ∼150 mK from T HO , the resistance of the strained sample is independent of the fit range used to calculate the elastoresistivity coefficients, consistent with the linear elastoresistive response. Inspection of the curves in Figure  8(c) indicates that T HO increases (decreases) for negative (positive) strains yy − xx . This effect is consistent with expectations for a phase transition involving a multicomponent order parameter and with the model developed in Section I, as can be appreciated by inspection of Equation (A7). The analysis indicates a change in T HO of approximately 100 mK per 10 −3 strain. Close to T HO , such a strain-induced change in T HO would yield a nonlinear response, possibly accounting for some of the nonlinearity at T HO . Significantly, for smaller strains, the change in T HO would be correspondingly smaller. Since the anomaly in 2m 66 is observed for large and small strain ranges alike (Figure 8(b)), and since the response is linear for all temperatures more than ∼150 mK from T HO (in comparison to the width of the anomaly, which extends for over 1000 mK), we conclude that the anomaly itself is not primarily associated with strain-induced changes in T HO . Nevertheless, the nonlinear effects described above preclude a critical scaling analysis, at least for the strains employed in the current study. ). There is a clear sign change in the slope of ( ∆R R )yy at THO = 17.15 K. In a small temperature range (±0.15 K) about THO, the resistive response to applied strain is nonlinear. The slopes in (a) are fit and plotted in (b) for different strain regions. Black points refer to the full strain region of −50 V to +150 V (which correspond to strains ∼ −0.25×10 −4 to 1 × 10 −4 ). Orange points in (b) are from slopes fit over a 0 V to 150 V strain region in (a) (strains of ∼ 0 to 1 × 10 −4 ). Blue points in (b) are from slopes fit over a −50 V to 0 V strain region in (a) (strains of ∼ −0.25 × 10 −4 to 0). Red points in (b) are from slopes fit over a −50 V to 50 V region in (a) (strains of ∼ −0.25 × 10 −4 to 0.25 × 10 −4 ). Comparing the black points and the red points in (b), it can be seen that the magnitude of the elastoresistivity coefficients are slightly underestimated if the data are fit to a linear function over the entire strain range and overestimated if fitting only the negative strain region. In (c) we plot resistance as a function of temperature for a fixed strain in the [110] sample (see main text for further details). From the shift in the point of inflection in the resistance, application of +5 × 10 −4 strain (red curve) and of −5 × 10 −4 strain (green curve) shifts the transition temperature THO by approximately 100 mK per 10 −3 of strain.