#  Divergent Thermal Expansion and Grüneisen Ratio in a Quadrupolar Kondo
Metal

A. Wörl<EMAIL_ADDRESS>Experimental Physics VI, Center
for Electronic Correlations and Magnetism, University of Augsburg, 86159
Augsburg, Germany  M. Garst  Institute for Theoretical Solid State Physics,
Karlsruhe Institute of Technology, 76131 Karlsruhe, Germany  Institute for
Quantum Materials and Technology, Karlsruhe Institute of Technology, 76131
Karlsruhe, Germany  Y. Yamane  Current affiliation: Department of Material
Science, Graduate School of Science, University of Hyogo, Kamigori, Hyogo
678-1297, Japan  Department of Quantum Matter, Graduate School of Advanced
Science and Engineering, Hiroshima University, Higashi-Hiroshima 739-8530,
Japan  S. Bachus  Experimental Physics VI, Center for Electronic Correlations
and Magnetism, University of Augsburg, 86159 Augsburg, Germany  T. Onimaru
Department of Quantum Matter, Graduate School of Advanced Science and
Engineering, Hiroshima University, Higashi-Hiroshima 739-8530, Japan  P.
Gegenwart<EMAIL_ADDRESS>Experimental Physics VI,
Center for Electronic Correlations and Magnetism, University of Augsburg,
86159 Augsburg, Germany

###### Abstract

We report on the low-temperature thermal expansion and magnetostriction of the
single-impurity quadrupolar Kondo candidate Y1-xPrxIr2Zn20. In the dilute
limit, we find a quadrupolar strain that possesses a singular dependence on
temperature $T$, $\varepsilon_{\mathrm{u}}\sim H^{2}\log 1/T$, for a small but
finite magnetic field $H$. Together with the previously reported anomalous
specific heat $C$, this implies a quadrupolar Grüneisen ratio
$\Gamma_{u}=\partial_{T}\varepsilon_{\mathrm{u}}/C\sim H^{2}/(T^{2}\log 1/T)$
whose divergence for finite $H$ is consistent with the scenario of a
quadrupolar Kondo effect. In addition, we find a singular behavior of the
isotropic strain $\varepsilon_{\mathrm{B}}$ in zero magnetic field resulting
in a divergence of both the volume thermal expansion and the volume Grüneisen
parameter. We speculate that this behavior might be also induced by putative
Kondo correlations via elastic anharmonicities or static strain disorder.

In solid state physics, non-Fermi liquid phases describe unconventional
metallic states of matter. Such exotic phases have been extensively studied in
the framework of heavy fermion (HF) quantum criticality [1]. A magnetic
quantum critical point typically forms in Ce- and Yb-based intermetallic
systems that are located in between a Kondo screened Fermi-liquid and a long
range antiferromagnetically ordered state. The divergence of the Grüneisen
parameter, defined as the ratio of volume thermal expansion to specific heat,
is a universal signature of pressure sensitive quantum critical points [2, 3].

In addition, there has been a particular focus on exploring non-Fermi liquid
states that are related to electric quadrupole moments and on identifying
potential links between these states and magnetic HF quantum criticality [4,
5, 6, 7]. Quadrupolar ground states form, for instance, in non-Kramers Pr- and
U-based intermetallics, given that certain symmetry constraints are fulfilled.
Possible cause for unconventional metallic behavior in these materials is the
quadrupolar Kondo effect, which was postulated by Cox in 1987 [8]. Here, the
simultaneous overscreening of a localized quadrupole moment by two-channels of
conduction electrons, which are related to their spin degree of freedom, leads
to non-Fermi liquid behavior in specific heat ($C/T\propto\log 1/T$),
electrical resisitivity ($\rho/\rho_{0}\propto 1+A\sqrt{T}$) and quadrupole
susceptibility ($\chi_{\mathrm{Q}}\propto\log{1/T}$) as well as an
unconventional residual entropy of $S=R\log\sqrt{2}$ [9].

Cubic Pr-based 1-2-20 systems are prototypical to study such novel quadrupolar
[10, 11] as well as higher multipolar related correlation effects [12, 13,
14]. PrIr2Zn20, for instance, has a well defined quadrupolar non-Kramers
$\Gamma_{3}$ ground state doublet (point group $T_{d}$) [15] and displays
clear signatures of the quadrupolar Kondo lattice effect [4], which are cut
off by antiferroquadrupolar order at 0.11 K [16]. Recent studies on highly
diluted Y1-xPrxIr2Zn20 provided direct evidence of the single-impurity
quadrupolar Kondo effect, based on characteristic non-Fermi liquid behaviors
found in the specific heat, electrical resistivity and elastic constant [5,
17]. Nevertheless, evidence for the residual entropy $S=R\log\sqrt{2}$ remains
elusive so far [5].

The concept of the Grüneisen parameter can be generalized to all irreducible
representations of elastic strains. Such a generalized Grüneisen parameter is
expected to diverge close to a quantum critical point provided that the
associated stress couples to a relevant operator of the critical fixed-point
[2, 18]. For the single-impurity quadrupolar Kondo fixed-point, this is the
case for the quadrupolar stress $\sigma_{\mathrm{u}}$, that breaks the cubic
symmetry, and, as a consequence, destabilizes the non-Fermi liquid physics and
quenches the residual entropy. This implies that the quadrupolar Grüneisen
ratio $\Gamma_{\mathrm{u}}=(\partial_{T}\varepsilon_{\mathrm{u}})/C$ diverges
in a characteristic manner. Here,
$\varepsilon_{\mathrm{u}}=(2\varepsilon_{zz}-\varepsilon_{xx}-\varepsilon_{yy})/\sqrt{3}$
is the strain component of the $\Gamma_{3}$ doublet that is conjugate to
$\sigma_{\mathrm{u}}$ with the corresponding thermal expansion
$\alpha_{\mathrm{u}}=\partial_{T}\varepsilon_{\mathrm{u}}$ and
$C=C_{\mathrm{m}}/V_{\mathrm{m}}$ the specific heat with the molar volume
$V_{\mathrm{m}}$ and the molar specific heat $C_{\mathrm{m}}$. Using a
thermodynamic identity,
$\Gamma_{\mathrm{u}}=1/T(\mathrm{d}T/\mathrm{d}\sigma_{\mathrm{u}})_{S}$ also
quantifies the adiabatic change of temperature upon the variation of
$\sigma_{\mathrm{u}}$ and thus describes an elastocaloric effect [19].

The application of a magnetic field $\boldsymbol{H}$ will split the non-
Kramers $\Gamma_{3}$ doublet resulting in a finite quadrupole moment $\langle
Q\rangle\sim\chi_{\mathrm{Q}}H^{2}$ of each Pr3+ ion, that is proportional to
$H^{2}$ for small fields due to time-reversal symmetry. This induces a local
quadrupolar strain via the elastic coupling $g_{\Gamma_{3}}$, and after
averaging over Pr disorder results in a homogeneous strain
$\varepsilon_{\mathrm{u}}\sim
g_{\Gamma_{3}}n_{\mathrm{Pr}}\chi_{\mathrm{Q}}H^{2}/c_{\mathrm{u}}$ where
$n_{\mathrm{Pr}}$ is the Pr3+ density and $c_{\mathrm{u}}=(c_{11}-c_{12})/2$
the corresponding elastic constant. The quadrupolar Kondo effect predicts
$\chi_{\mathrm{Q}}\sim\log 1/T$ resulting in a singular temperature dependence
of the quadrupolar strain $\varepsilon_{\mathrm{u}}\sim
n_{\mathrm{Pr}}H^{2}\log 1/T$. Finally, using that the specific heat at low
temperatures is dominated by the contribution of the Pr3+ ions,
$C_{\mathrm{m}}\sim n_{\mathrm{Pr}}T\log 1/T$ for small magnetic fields [5],
one expects for the quadrupolar Grüneisen parameter $\Gamma_{\mathrm{u}}\sim
H^{2}/(T^{2}\log{1/T})$. Consequently, thermal expansion and magnetostriction
experiments are ideally suited for the investigation of the non-Fermi liquid
behavior associated with the quadrupolar Kondo effect.

In this work, we perform such measurements on highly diluted single
crystalline Y1-xPrxIr2Zn20. As our key finding, we confirm experimentally that
the temperature dependence for the quadrupolar strain
$\varepsilon_{\mathrm{u}}$ and the associated quadrupolar Grüneisen parameter
$\Gamma_{\mathrm{u}}$ is indeed consistent with the suggested quadrupolar
Kondo scenario [5, 17]. Remarkably, we also find that the volume strain
$\varepsilon_{\mathrm{B}}$ shows a similar singular behavior in zero magnetic
field although isotropic stress does not directly couple to a relevant
operator of the quadrupolar Kondo fixed-point. We speculate that an indirect
coupling might be generated either by elastic anharmonicities or by static
strain disorder accounting for the experimentally observed divergence in the
volume thermal expansion.

Ultrahigh-resolution thermal expansion and magnetostriction measurements were
carried out in a dilution refrigerator using a miniaturized capacitive
dilatometer [20]. Central to this study is a highly diluted Y1-xPrxIr2Zn20
single crystal with $x=0.036$, for which relative length changes were measured
along a cubic $\langle 100\rangle$ direction, with magnetic fields applied
either parallel or perpendicular to the measurement direction. The magnetic
field direction is defined as [001] in the following, so that
$\varepsilon_{zz}$ and $\varepsilon_{xx}$ correspond to longitudinal
$\varepsilon_{\parallel}$ and transverse strain $\varepsilon_{\perp}$,
respectively. As the applied field does not break the symmetry within the
$(x,y)$ plane, we can assume that $\varepsilon_{xx}=\varepsilon_{yy}$. The
measurement of longitudinal and transverse strain in magnetic field then
allows to infer the values of the isotropic strain
$\varepsilon_{\mathrm{B}}=\varepsilon_{xx}+\varepsilon_{yy}+\varepsilon_{zz}$
and the quadrupolar strain $\varepsilon_{\mathrm{u}}$, as detailed in the
Supplemental Material (SM) [21]. An unavoidable side effect of the
experimental technique is a small force of approximately $4$ N [20] acting on
the sample along the measurement direction, which corresponds to a tiny
unixial stress of a few MPa. To clarify whether this effect has an impact on
the deduced relative length changes, we performed complementary measurements
on a $[111]$ oriented single crystal with a comparable Pr-concentration of
$x=0.033$. For further characterization, a special uniaxial stress capacitive
dilatometer [27], that exerts a roughly fifteen times larger force on the
sample than the miniaturized dilatometer, was employed. For details on the
single crystalline samples examined in this study, see SM [21].

First, we discuss the quadrupolar thermal expansion coefficient
$\alpha_{\mathrm{u}}=\partial_{T}\varepsilon_{\mathrm{u}}$ of Y1-xPrxIr2Zn20
with $x=0.036$. The respective data is derived from the measurement of the
longitudinal $\alpha_{\parallel}$ and the volume thermal expansion $\beta$ for
$\boldsymbol{H}\parallel[001]$ via the the relation
$\alpha_{\mathrm{u}}=\sqrt{3}(\alpha_{\parallel}-\beta/3)$. A detailed
derivation of this relation and an overview of the data of
$\alpha_{\parallel}$ and $\beta$ is provided in the SM [21].

Figure 1: (a) Temperature dependence of the quadrupolar thermal expansion
coefficient $\alpha_{\mathrm{u}}$ at various magnetic fields
$\boldsymbol{H}\parallel[001]$. The dashed lines are crystal electric field
(CEF) calculations. (b) Quadrupolar Grüneisen parameter normalized to magnetic
field $\Gamma_{\mathrm{u}}/H^{2}$ as a function of temperature for
$\boldsymbol{H}\parallel[001]$. The black solid line denotes the theoretically
expected temperature dependence of $\Gamma_{\mathrm{u}}$ for the quadrupolar
Kondo effect. The inset shows $\alpha_{\mathrm{u}}/H^{2}$ vs. $T$ on a log-log
scale, together with power law divergences as dashed and solid lines.

The quadrupolar thermal expansion coefficient $\alpha_{\mathrm{u}}$ is shown
in Fig. 1(a) on a logarithmic temperature scale ranging from 0.05 K to 4 K for
various magnetic fields up to 10 T. For small fields, $\alpha_{\mathrm{u}}$
increases down to lowest temperatures. At around 1.5 T a maximum develops
whose position shifts to higher temperatures as a function of increasing
field. At the same time, the height of the maximum continuously decreases with
$H$. For $\mu_{0}H\gtrsim 8$ T, this peak is quantitatively captured by
crystal electric field (CEF) calculations [21] as indicated by the dashed
lines.

The ratio $\alpha_{\mathrm{u}}/(\mu_{0}H)^{2}$ is shown in the inset of Fig.
1(b). For temperatures $k_{\mathrm{B}}T\gg\mu_{\mathrm{B}}\mu_{0}H$, the data
collapses onto a single curve that exhibits a crossover from a
$\alpha_{\mathrm{u}}/H^{2}\sim\partial_{T}\chi_{\mathrm{Q}}\sim 1/T$ behavior
at low temperature, consistent with the quadrupolar Kondo effect, to a
$1/T^{2}$ dependence at high temperature as expected for a Curie
susceptibility $\chi_{\mathrm{Q}}\sim 1/T$ of a fully localized $\Gamma_{3}$
doublet[21]. The crossover temperature of $\sim 0.6$ K is consistent with
previous studies [5, 17]. By using the molar 4f specific heat $C_{\mathrm{m}}$
[21], measured on a single crystal from the same batch with a comparable Pr-
concentration of $x=0.044$, we evaluate the quadrupolar Grüneisen parameter
that follows the expected singular behavior $\Gamma_{\mathrm{u}}\sim
H^{2}/(T^{2}\log 1/T)$ for low magnetic fields as shown by the black solid
line in Fig. 1(b). The crossover field of approximately 2 T is comparable to
the crossover field as suggested by the measurement of $\alpha_{\mathrm{u}}$.

Figure 2: Quadrupolar thermal expansion coefficient normalized to magnetic
field $\alpha_{\mathrm{u}}/(\mu_{0}H)^{2}$ as a function of temperature for
magnetic field $\boldsymbol{H}\parallel[001]$. The dashed solid line indicates
the $T^{-1}$ single-impurity quadrupole Kondo dependence. Arrows denote a
characteristic temperature $T^{*}$, below which deviations from the universal
single-impurity quadrupole Kondo behavior arise. The inset displays $T^{*}$ as
a function of magnetic field, with the red solid line indicating a $H^{2}$
field dependence.

Within the framework of the quadrupolar Kondo model, the deviation from the
characteristic quadrupole Kondo behavior for larger magnetic fields can be
attributed to both a channel- and a quadrupolar-asymmetry, induced in linear
and quadratic order in $H$, respectively. At small fields, the first effect is
expected to dominate and results in a crossover temperature $T^{*}\sim H^{2}$,
see Ref. [9], separating non-Fermi and Fermi-liquid behavior. In order to
specify the origin of the magnetic field induced crossover in highly diluted
Y1-xPrxIr2Zn20, Fig. 2 displays the data of
$\alpha_{\mathrm{u}}/(\mu_{0}H)^{2}$ on a logarithmic temperature scale at
small magnetic fields up to 2.5 T. The crossover temperature $T^{*}$, which is
determined as the temperature at which deviations from the universal single-
impurity quadrupole Kondo behavior arise, is indicated by an arrow for each
magnetic field. The inset displays the characteristic temperature $T^{*}$ as a
function of magnetic field, whereby the red solid line indicates a quadratic
magnetic field dependence. Indeed, the characteristic temperature $T^{*}$
estimated from the quadrupolar thermal expansion data can be well scaled with
$T^{*}\sim H^{2}$. This indicates that the field induced channel-asymmetry is
the dominating perturbation at low magnetic field leading to the
experimentally found deviations from the quadrupole Kondo behavior at low
temperature, which is in very good agreement with the theoretical expectation.

Figure 3: Magnetic field variation of the quadrupolar magnetostriction
$\varepsilon_{\mathrm{u}}$ for $\boldsymbol{H}\parallel[001]$ at different
temperatures. The dashed lines are CEF calculations. The inset shows the
quadrupole susceptibility $\chi_{\mathrm{Q}}$ extracted from the quadratic
dependence $\varepsilon_{\mathrm{u}}\sim
n_{\mathrm{Pr}}\chi_{\mathrm{Q}}H^{2}$ at small $H$ for three different Pr
concentrations $x$. The CEF prediction is shown as a solid blue line, whose
temperature dependence is characterized by the superposition of a constant Van
Vleck, $\chi_{\mathrm{VV}}$, and a $C/T$ Curie contribution which cannot
describe the experimental data at low $T$.

The results on the quadrupolar thermal expansion coefficient
$\alpha_{\mathrm{u}}$ and the quadrupolar Grüneisen parameter
$\Gamma_{\mathrm{u}}$, which are indicative of the single-impurity quadrupole
Kondo effect, are further corroborated by the quadrupolar magnetostriction
$\varepsilon_{\mathrm{u}}$ displayed in Fig. 3. The quadrupolar
magnetostriction coefficient $\varepsilon_{\mathrm{u}}$ is derived from the
measurement of the longitudinal $\varepsilon_{\parallel}$ and the volume
magnetostriction coefficient $\varepsilon_{\mathrm{B}}$ by using the relation
$\varepsilon_{\mathrm{u}}=\sqrt{3}(\varepsilon_{\parallel}-\varepsilon_{\mathrm{B}}/3)$
for $\boldsymbol{H}\parallel{[001]}$. The data of $\varepsilon_{\parallel}$
and $\varepsilon_{\mathrm{B}}$ used for the calculation of
$\varepsilon_{\mathrm{u}}$ is provided in the SM [21]. By analyzing the
initial quadratic field dependence of $\varepsilon_{\mathrm{u}}\sim n_{\rm
Pr}\chi_{\mathrm{Q}}H^{2}$ [21], we extract the quadrupolar susceptibility
$\chi_{\mathrm{Q}}$ that is shown for various Pr doping $x$ in the inset of
Fig. 3. The CEF calculation, that predicts a Curie-like $1/T$ temperature
dependence on top of a constant Van Vleck contribution, can only capture the
behavior of $\chi_{\mathrm{Q}}$ at elevated temperature. At low $T$ and for
the lowest doping concentration $x=0.036$ our results are consistent with a
logarithmic temperature dependence of $\chi_{\mathrm{Q}}$. For higher doping,
the susceptibility is suppressed at low temperatures indicating a quenching of
the quadrupolar Kondo effect, likely induced by interactions between Pr3+
ions.

We now turn to the discussion of the volume thermal expansion coefficient
$\beta=\partial_{T}\varepsilon_{\mathrm{B}}$. In Fig. 4(a), $\beta$ is shown
for the highly diluted single crystal with $x=0.036$. Remarkably, the thermal
expansion coefficient at zero field increases down to lowest temperatures.
This continuous increase is quenched by the application of a field and, for
$\mu_{0}H\gtrsim 6$ T, the volume thermal expansion practically vanishes.
Similarly, $\beta$ is also suppressed by larger doping as shown in the inset
for zero field.

Figure 4:  (a) Volume thermal expansion coefficient $\beta$ as a function of
temperature at various magnetic fields $\boldsymbol{H}\parallel[001]$. The
inset shows $\beta$ at zero magnetic field for three different Pr-
concentrations $x$. (b) Temperature dependencies of the uniaxial thermal
expansion coefficients $\alpha_{[001]}$ and $\alpha_{[111]}$ measured with the
miniaturized dilatometer at zero magnetic field. (c) Temperature dependencies
of $\alpha_{[001]}$ measured with the miniaturized (closed circles) and the
uniaxial stress dilatometer (open circles). (d) Temperature dependencies of
$\alpha_{[111]}$ obtained by the miniaturized (filled circles) and the
uniaxial stress dilatometer (open circles) together with specific heat data,
partially taken from Ref. [5].The blue solid line is a fit to $\alpha_{[111]}$
using $\alpha_{\mathrm{fit}}=\gamma C_{\mathrm{\mathrm{m}}}+bT^{-1}$. The
inset displays the temperature dependence of the bulk Grüneisen parameter
$\Gamma_{\mathrm{B}}$.

The continuous increase of $\beta$ on cooling at zero field and low doping is
unexpected at first sight as the influence of the quadrupolar moments should
in lowest order average out in the bulk thermal expansion. A parasitic signal
could potentially arise from the presence of a small uniaxial stress from the
sample holder that is on the order of a few MPa in case of the miniaturized
dilatometer and breaks the isotropy of the environment. In order to clarify
this influence, we compare in Fig. 4(b) the uniaxial thermal expansion along
crystallographic $[001]$ and $[111]$ directions measured on two single
crystals with almost the same doping level. While a uniaxial stress along
$[111]$ only activates the bulk strain $\varepsilon_{\mathrm{B}}$, a uniaxial
stress along $[001]$ induces both bulk and quadrupolar strains where the
latter couples linearly to the $\Gamma_{3}$ ground state doublet of the Pr3+
ions. At very low temperature, we find indeed a difference between
$\alpha_{[001]}$ and $\alpha_{[111]}$ that we attribute to the presence of a
small uniaxial stress exerted by the dilatometer.

To demonstrate the tremendous effect of uniaxial stress along [001], we
performed a complementary measurement using a uniaxial stress dilatometer with
an applied uniaxial stress approximately fifteen times larger than in the
miniaturized dilatometer, see Fig. 4(c). Under higher uniaxial stress, the
divergence of $\alpha_{[001]}$ is suppressed and a maximum forms at around
0.11 K, suggesting the breakdown of the two-channel Kondo effect due to the
linear in strain splitting of the $\Gamma_{3}$ ground state doublet.

Figure 4(d) displays the temperature dependence of $\alpha_{[111]}$ measured
with both the miniaturized and the uniaxial stress dilatometer at zero
magnetic field. As there is no notable difference between the two data sets,
we conclude that the divergent temperature dependence of $\alpha_{[111]}$ is
in fact an intrinsic volume effect and not artificially induced by the
uniaxial stress applied externally along the measurement direction. This
result is striking as it is markedly different from the molar specific heat of
the 4f electrons $C_{\mathrm{m}}$ as shown by the red symbols. The thermal
expansion is well described by a fit $\alpha_{\mathrm{fit}}=\gamma
C_{\mathrm{m}}+b/T$ (blue line) that assumes two contributions. The first
obeys Grüneisen scaling with a constant $\gamma$ and the second contribution
diverges like $1/T$, which is reminiscent of the behavior of
$\alpha_{\mathrm{u}}/H^{2}$ shown in the inset of Fig. 1(b). Consequently, the
bulk Grüneisen ratio
$\Gamma_{\mathrm{B}}=V_{\mathrm{m}}3\alpha_{[111]}/C_{\mathrm{m}}\approx
3V_{\mathrm{m}}\gamma+3V_{\mathrm{m}}b/(TC_{\mathrm{m}})$ exhibits a
divergence as a function of temperature similar to $\Gamma_{\mathrm{u}}/H^{2}$
shown in Fig. 1(b).

In the following, we speculate on the origin of this surprising finding of an
intrinsic bulk thermal expansion that increases down to lowest temperatures.
The elastic coupling $g_{\Gamma_{3}}$ to the non-Kramers $\Gamma_{3}$ doublet
leads in second order perturbation theory to the following correction to the
elastic Hamiltonian

$\displaystyle\delta\mathcal{H}=-\frac{g_{\Gamma_{3}}^{2}\chi_{\mathrm{Q}}}{2}\sum_{\alpha}(\varepsilon^{2}_{\mathrm{u}}(\vec{r}_{\alpha})+\varepsilon^{2}_{\mathrm{v}}(\vec{r}_{\alpha})).$
(1)

Here, $\varepsilon_{\mathrm{u/v}}(\vec{r}_{\alpha})$ are the two components of
the local quadrupolar strain doublet at the position $\vec{r}_{\alpha}$ of the
Pr3+ ion with index $\alpha$. After averaging over Pr disorder this
Hamiltonian recovers the renormalization of the quadrupolar elastic constant
$\delta c_{\mathrm{u}}=-n_{\mathrm{Pr}}g_{\Gamma_{3}}^{2}\chi_{\mathrm{Q}}$
whose temperature dependence was confirmed by elastic constant measurements
[17]. Treating the Hamiltonian of Eq. (1) in perturbation theory, one obtains
a correction to the free energy density of the form $\delta f\sim\frac{\delta
c_{u}}{2}\langle\varepsilon_{\mathrm{u/v}}^{2}\rangle$. Here, a finite
expectation value $\langle\varepsilon_{\mathrm{u/v}}^{2}\rangle$ might arise
either from dynamic zero-point fluctuations of the quadrupolar strain or from
static strain disorder [21]. This results in a contribution to the bulk
thermal expansion $\delta\beta=\partial_{p}\partial_{T}\delta
f=-\frac{1}{2}(\partial_{p}n_{\mathrm{Pr}}g_{\Gamma_{3}}^{2}\langle\varepsilon_{\mathrm{u/v}}^{2}\rangle)\partial_{T}\chi_{\mathrm{Q}}$.
The divergence $\delta\beta\sim 1/T$ originating from $\chi_{Q}\sim\log 1/T$
is consistent with the observations in Fig. 4(d).

Hence, dynamic as well as static strain fluctuations can in principle explain
the observed singular behavior in $\beta$. However, both scenarios also imply
a contribution to the specific heat $\delta C_{\mathrm{m}}\sim 1/T$, which so
far has not been successfully identified. A possible reason is a small
prefactor of $\delta C_{\mathrm{m}}$, or, equivalently, a large effective
Grüneisen parameter $\delta\Gamma_{\mathrm{B}}=\delta\beta/\delta
C_{\mathrm{m}}$ characterizing the hydrostatic pressure dependence of
$n_{\mathrm{Pr}}g_{\Gamma_{3}}^{2}\langle\varepsilon_{\mathrm{u/v}}^{2}\rangle$.
This dependence can arise from elastic anharmonicities, e.g., via the elastic
coupling $g_{\Gamma_{3}}(p)$ or, in case of strain disorder, from a response
of $n_{\mathrm{Pr}}\langle\varepsilon_{\mathrm{u/v}}^{2}\rangle$ to pressure
changes. From Fig. 4(d) we can estimate a lower bound
$\Gamma_{\mathrm{B}}>2.15$ GPa-1, and with the isothermal compressibility
$\kappa_{T}\approx 0.01$ GPa-1, estimated from elastic constant measurements
on PrIr2Zn20 [28], this implies a dimensionless Grüneisen parameter
$\delta\Gamma_{\mathrm{B}}\kappa_{T}$ at least of the order of 200 to be
compared with a typical value of 1 for metals [29]. As elastic anharmonicities
are in general unlikely to account for such large values, this suggests static
strain disorder to be at its origin. Hydrostatic pressure will influence the
local strain distribution and might thus affect sensitively the quadrupolar
environment of each Pr ions resulting in a large response. Local strain fields
were also invoked to explain the temperature and field dependence of the
specific heat [5] and the elastic constant [17], respectively. Future nuclear
quadrupolar resonance (NQR) measurements might help clarifying this issue.

In this work, we established the symmetrized quadrupolar expansivity, i.e.,
thermal expansion and magnetostriction, as well as the quadrupolar Grüneisen
parameter as a highly sensitive probe to uncover quadrupolar fluctuations. For
the example of highly diluted Y1-xPrxIr2Zn20 we demonstrated that it gives
access to the quadrupolar susceptibility. Our experimental results are in
excellent agreement with expectations for the single-impurity quadrupolar
Kondo effect and corroborate its emergence in this material. On a more general
note, we anticipate the quadrupolar expansivity and quadrupolar Grüneisen
parameter to become an important tool also for the investigation of other
materials with pronounced quadrupolar correlations, in particular, for
materials close to a nematic quantum critical point as discussed for certain
cuprate and iron pnictide superconductors [30, 31].

###### Acknowledgements.

Long term collaboration with R. Küchler as well as helpful discussions with Y.
Tokiwa, R. Yamamoto, T. Yanagaisawa and H. Kusunose are gratefully
acknowledged. The work in Augsburg was supported by the German Research
Foundation (DFG) via the Project No. 107745057 (TRR80). The work in Hiroshima
was supported by the Center for Emergent Condensed-Matter Physics (ECMP),
Hiroshima University and financially supported by Grants-in-Aid from MEXT/JSPS
of Japan [Grant Nos. JP26707017, JP15H05886 (J-Physics), JP15KK0169,
JP18H01182, and JP18KK0078]. M.G. is supported by the DFG via the Project No.
422213477 (TRR 288, project A11).

## References

  * Gegenwart  _et al._ [2008] P. Gegenwart, Q. Si, and F. Steglich, Nat. Phys. 4, 186 (2008).
  * Zhu  _et al._ [2003] L. Zhu, M. Garst, A. Rosch, and Q. Si, Phys. Rev. Lett. 91, 066404 (2003).
  * Gegenwart [2016] P. Gegenwart, Rep. Prog. Phys. 79, 114502 (2016).
  * Onimaru  _et al._ [2016] T. Onimaru, K. Izawa, K. T. Matsumoto, T. Yoshida, Y. Machida, T. Ikeura, K. Wakiya, K. Umeo, S. Kittaka, K. Araki, T. Sakakibara, and T. Takabatake, Phys. Rev. B 94, 075134 (2016).
  * Yamane  _et al._ [2018a] Y. Yamane, T. Onimaru, K. Wakiya, K. T. Matsumoto, K. Umeo, and T. Takabatake, Phys. Rev. Lett. 121, 077206 (2018a).
  * Matsubayashi  _et al._ [2012] K. Matsubayashi, T. Tanaka, A. Sakai, S. Nakatsuji, Y. Kubo, and Y. Uwatoko, Phys. Rev. Lett. 109, 187004 (2012).
  * Shimura  _et al._ [2015] Y. Shimura, M. Tsujimoto, B. Zeng, L. Balicas, A. Sakai, and S. Nakatsuji, Phys. Rev. B 91, 241102(R) (2015).
  * Cox [1987] D. L. Cox, Phys. Rev. Lett. 59, 1240 (1987).
  * Cox and Zawadowski [1998] D. L. Cox and A. Zawadowski, Adv. Phys. 47, 599 (1998).
  * Onimaru and Kusunose [2016] T. Onimaru and H. Kusunose, J. Phys. Soc. Jpn. 85, 082002 (2016).
  * Sakai and Nakatsuji [2011] A. Sakai and S. Nakatsuji, J. Phys. Soc. Jpn. 80, 063701 (2011).
  * Tsujimoto  _et al._ [2014] M. Tsujimoto, Y. Matsumoto, T. Tomita, A. Sakai, and S. Nakatsuji, Phys. Rev. Lett. 113, 267001 (2014).
  * Patri  _et al._ [2019] A. S. Patri, A. Sakai, S. B. Lee, A. Paramekanti, S. Nakatsuji, and Y. B. Kim, Nat. Commun. 10, 4092 (2019).
  * Patri and Kim [2020] A. S. Patri and Y. B. Kim, Phys. Rev. X 10, 041021 (2020).
  * Iwasa  _et al._ [2013] K. Iwasa, H. Kobayashi, T. Onimaru, K. T. Matsumoto, N. Nagasawa, T. Takabatake, S. Ohira-Kawamura, T. Kikuchi, Y. Inamura, and K. Nakajima, J. Phys. Soc. Jpn. 82, 043707 (2013).
  * Onimaru  _et al._ [2011] T. Onimaru, K. T. Matsumoto, Y. F. Inoue, K. Umeo, T. Sakakibara, Y. Karaki, M. Kubota, and T. Takabatake, Phys. Rev. Lett. 106, 177001 (2011).
  * Yanagisawa  _et al._ [2019] T. Yanagisawa, H. Hidaka, H. Amitsuka, S. Zherlitsyn, J. Wosnitza, Y. Yamane, and T. Onimaru, Phys. Rev. Lett. 123, 067201 (2019).
  * Zacharias  _et al._ [2015] M. Zacharias, I. Paul, and M. Garst, Physical Review Letters 115, 025703 (2015).
  * Ikeda  _et al._ [2019] M. S. Ikeda, J. A. W. Straquadine, A. T. Hristov, T. Worasaran, J. C. Palmstrom, M. Sorensen, P. Walmsley, and I. R. Fisher, Review of Scientific Instruments 90, 083902 (2019).
  * Küchler  _et al._ [2017] R. Küchler, A. Wörl, P. Gegenwart, M. Berben, B. Bryant, and S. Wiedmann, Rev. Sci. Instr. 88, 083903 (2017).
  * [21] See Supplemental Material for details on the experimental methods, crystal electric field calculations, theory on the volume thermal expanison, data of the longitudinal, transverse and volume thermal expansion and magnetostriction and the quadrupole field susceptibility analysis, which includes Ref. [22, 23, 24, 25, 26].
  * Yamane  _et al._ [2018b] Y. Yamane, T. Onimaru, K. Uenishi, K. Wakiya, K. T. Matsumoto, K. Umeo, and T. Takabatake, Physica B: Condensed Matter 536, 40 (2018b).
  * Yamane  _et al._ [2018c] Y. Yamane, T. Onimaru, K. Wakiya, K. T. Matsumoto, K. Umeo, and T. Takabatake, AIP Adv. 8, 101338 (2018c).
  * Pott and Schefzyk [1983] R. Pott and R. Schefzyk, Journal of Physics E: Scientific Instruments 16, 444 (1983).
  * Wörl  _et al._ [2019] A. Wörl, T. Onimaru, Y. Tokiwa, Y. Yamane, K. T. Matsumoto, T. Takabatake, and P. Gegenwart, Phys. Rev. B 99, 081117(R) (2019).
  * Lea  _et al._ [1962] K. Lea, M. Leask, and W. Wolf, J. Phys. Chem. Solids 23, 1381 (1962).
  * Küchler  _et al._ [2016] R. Küchler, C. Stingl, and P. Gegenwart, Rev. Sci. Instr. 87, 073903 (2016).
  * Ishii  _et al._ [2011] I. Ishii, H. Muneshige, Y. Suetomi, T. K. Fujita, T. Onimaru, K. T. Matsumoto, T. Takabatake, K. Araki, M. Akatsu, Y. Nemoto, T. Goto, and T. Suzuki, J. Phys. Soc. Jpn. 80, 093601 (2011).
  * Ashcroft and Mermin [1976] N. Ashcroft and N. Mermin,  _Solid State Physics_ (Saunders College Publishing, 1976) Chapt. 25.
  * Auvray  _et al._ [2019] N. Auvray, B. Loret, S. Benhabib, M. Cazayous, R. Zhong, J. Schneeloch, G. Gu, A. Forget, D. Colson, I. Paul, A. Sacuto, and Y. Gallais, Nat. Commun. 10, 5209 (2019).
  * Worasaran  _et al._ [2021] T. Worasaran, M. S. Ikeda, J. C. Palmstrom, J. A. W. Straquadine, S. A. Kivelson, and I. R. Fisher, Science 372, 973 (2021).

Supplemental Material for
“Divergent Thermal Expansion and Grüneisen Ratio in a Quadrupolar Kondo Metal”
A. Wörl, M. Garst, Y. Yamane, S. Bachus, T. Onimaru, and P. Gegenwart

##  I Experimental methods

This section provides detailed information on the investigated single
crystalline samples of Y1-xPrxIr2Zn20. Furthermore, the two different
capacitive dilatometer employed in this study are briefly introduced.

###  I.1 Y1-xPrxIr2Zn20 Single Crystals

The Y1-xPrxIr2Zn20 single crystals studied in this work were grown by use of a
Zn self-flux technique. For more details on the single crystal growth process,
we refer to Ref. [22]. In this study, we focus on four differently doped
Y1-xPrxIr2Zn20 single crystals with Pr-concentrations of $x=0.033,0.036,0.09$
and $0.49$. The shape, crystallographic directions and dimensions of each
single crystal are illustrated in Fig. S1, whereby the given length values
were determined with an accuracy of $\pm 0.001$ mm.

Figure S1: Sketches of the studied Y1-xPrxIr2Zn20 single crystals which
illustrate the shape, crystallographic directions and dimensions for each
single crystal.

The Pr concentrations $x$ of the single crystalline sample with $x=0.49$ was
determined by use of electron-probe microanalysis (EPMA), whereby the
uncertainty in the concentration is $\Delta x\pm 0.06$. As EPMA turned out as
not sensitive enough for characterizing the higher diluted samples with
$x=0.033$, $x=0.036$ and $x=0.09$, their magnetization was measured at 1.8 K
and 1 T and the respective Pr-concentration estimated by comparing the
experimental value with the theoretical one suggested by the crystal electric
field (CEF) calculation. These are the same characterization procedures as
already applied in previous studies [5, 22, 17]. The single crystal with
$x=0.044$ on which the specific heat presented in Fig. 4(d) of the main paper
was measured, is from the identical sample batch as the two highly diluted
single crystals with $x=0.033$ and $x=0.036$, on which the thermal expansion
was obtained. The uncertainty in the Pr concentration of the three highly
diluted single crystals with $x=0.033$, $x=0.036$ and $x=0.044$ is estimated
at $\Delta x=\pm 0.01$. The Pr content of these three single crystalline
samples can therefore be considered as the same within the estimated error of
the experimental technique used for the determination of the Pr concentration.
The uncertainty in the Pr concentration of the single crystalline sample with
$x=0.09$ is estimated at $\Delta x=\pm 0.02$.

The previously published values for the residual resistance ratio ($RRR$) of
Y1-xPrxIr2Zn20 imply that the highly diluted single crystals are of very high
quality. The $RRR$ of the highly diluted sample with $x=0.036$ is estimated at
160, which is the value determined previously for a sample ($x=0.044$) from
the same batch [22]. By increasing the Pr-concentration the $RRR$ decreases
continuously, which can be explained by increasing disorder. The $RRR$ of the
$x=0.49$ sample is estimated at only 8, which is the previously reported value
of a single crystal $(x=0.44)$ from the same batch [22].

###  I.2 Dilatometry

Ultrahigh resolution relative length change measurements were carried out by
employing a miniaturized capacitive dilatometer [20]. The dilatometer is made
of a copper beryllium alloy and has super compact dimensions of
$15\,\mathrm{mm}\times 14\,\mathrm{mm}\times 14.7\,\mathrm{mm}$ [20]. This
allows for its rotation inside the small sample space of a dilution
refrigerator and enables relative length change measurements both parallel and
perpendicular to magnetic field.

The relative length change $\Delta L$ of the sample was determined by
measuring the capacitance of the dilatometer as a function of temperature or
magnetic field. To convert the respective capacitance values into relative
length changes, we used the formula by Pott and Schefzyk [24], which takes
account of the maximal adjustable capacitance value $C_{\mathrm{max}}$ of the
dilatometer. This corrects a small error which originates from non perfectly
parallel capacitor plates. For the miniaturized capacitive dilatometer, the
maximal adjustable capacitance value is $C_{\mathrm{max}}=50$ pF. The
capacitor plate radius of the miniaturized dilatometer is 5 mm [20].

As the sample is clamped by four leaf springs into the dilatometer, a small
inevitable force acts on the sample in measurement direction. This force
varies with the adjusted capacitance value and equals $4$ N when the
capacitive dilatometer is operated at a capacitance value of $20\,$pF [20]. As
measurements on the different samples were performed at capacitance values
ranging between $19.8$ pF and $20.7$ pF, we simply assume that the same force
of $4$ N acts on all investigated samples. The uniaxial stress on the sample
is estimated as $\sigma=F/A$, where $F=4$ N is the force exerted on the sample
by the miniaturized capacitive dilatometer and $A$ and the surface area of the
sample. As the sample shapes are close to rectangular (cf. Fig. S1), we
approximately calculate the surface area by the lengths of the two sides
perpendicular to the measurement direction.

To quantify the impact of uniaxial stress on the measurement results, we
carried out additional thermal expansion measurements on the $x=0.036$ (along
[001]) and on the $x=0.033$ (along [111]) single crystals by use of a uniaxial
stress dilatometer [27]. This dilatometer is also manufactured from a copper
beryllium alloy and its working principle is identical to the previously
described miniaturized capacitive dilatometer, with the important difference
that its leaf springs are much stronger, therefore exerting a higher force on
the sample. Furthermore, it has a larger capacitor plate radius of 7 mm [27].
Relative length changes were calculated by the same procedure as discussed
above, whereby the maximal adjustable capacitance value of the uniaxial stress
dilatometer is $C_{\mathrm{max}}=90$ pF. The force acting on the sample in
measurement was estimated by means of the previously published relation
between capacitance and force [27], which yields approximately 59.1 N and 58.7
N for the measurements on the $x=0.033$ and $x=0.036$ samples, respectively.
After the measurement under high uniaxial stress the appearance of the two
investigated samples was checked. The $x=0.036$ sample did not show any
indication of damage. The $x=0.033$ sample exhibited small cracks.

The thermal expansion and magnetostriction curves presented in the main text
are background corrected, whereby the cell effect of the dilatometer was
determined by measuring relative length changes of the capacitive dilatometer
with a Cu specimen placed inside.

##  II Theoretical considerations

###  II.1 Crystal Electric Field Hamiltonian

The fully localized quadrupole moments of each Pr3+ ion ($J=4$) is described
by the crystal electric field (CEF) Hamiltonian [26] for the $T_{d}$ point
group of the cubic crystal symmetry

$\mathcal{H}_{\mathrm{CEF}}=W\left[\frac{x}{60}(O_{4}^{0}+5O_{4}^{4})+\frac{1-|x|}{1260}(O_{6}^{0}-21O_{6}^{4})\right]-g_{J}\mu_{\mathrm{B}}\boldsymbol{J}\boldsymbol{H}$
(S1)

where $O_{k}^{q}$ are the Stevens operators given in terms of the angular
momentum operator $\boldsymbol{J}$. The CEF parameters used to describe the
thermal expansion and magnetostriction of Y1-xPrxIr2Zn20 $(x=0.036)$ are
$W=-1.50$ K and $x=0.537$. The same parameters were used previously in order
to account for the temperature and field dependence of entropy for a sample
with a comparable Pr concentration of $x=0.044$ [23]. We note that the value
$W=-1.50$ K differs slightly from the value $W=-1.22$ K [15], which was
suggested by inelastic neutron scattering on pure PrIr2Zn20. Possible
explanation for the difference between the Y diluted system and pure system is
a strengthening of the CEF effect due the smaller ionic radius of Y as
compared to Pr. The second term of the Hamiltonian is the Zeeman term, where
$g_{J}$ is the Landé factor and $\mu_{\mathrm{B}}$ the Bohr magneton.

In the absence of a magnetic field $\boldsymbol{H}$, the groundstate of Eq.
(S1) is given by a degenerate non-Kramers doublet with quadrupolar
$\Gamma_{3}$ symmetry,

$\displaystyle|\Gamma_{3},+\rangle$
$\displaystyle=\sqrt{\frac{7}{24}}\Big{(}|4\rangle-\sqrt{\frac{10}{7}}|0\rangle+|-4\rangle\Big{)},$
(S2) $\displaystyle|\Gamma_{3},-\rangle$
$\displaystyle=\sqrt{\frac{1}{2}}\Big{(}|2\rangle+|-2\rangle\Big{)},$ (S3)

with the spin eigenstates $|m\rangle$ and $m\in\\{-4,-3,...,3,4\\}$. These
states can be associated with the Hilbert-space of a pseudospin operator
$\boldsymbol{\sigma}_{\Gamma_{3}}$. An important role is played by the Stevens
operators $O^{0}_{2}\propto 2J_{z}^{2}-J_{x}^{2}-J_{y}^{2}$ and
$O^{2}_{2}\propto J_{x}^{2}-J_{y}^{2}$. After projection onto the degenerate
groundstate doublet, these operators can be identified with components of the
pseudospin operator $\boldsymbol{\sigma}_{\Gamma_{3}}$. Consequently, their
thermal expectation value with respect to the non-Kramers doublet vanishes at
zero magnetic field, $\langle O^{q}_{2}\rangle=0$ for $q=0,2$. We will be
particularly interested in a finite magnetic field applied along the cubic
crystal axis $[001]$ that splits the doublet and leads to a finite expectation
value of $\langle O^{0}_{2}\rangle\sim\mathcal{O}(H^{2})$.

###  II.2 Quadrupolar Kondo coupling to the conduction electrons

The local pseudospin $\boldsymbol{\sigma}_{\Gamma_{3}}$ will couple with a
Kondo exchange coupling $J_{K}$ to the conduction electrons via two channels
corresponding to the two electronic spin-configurations potentially realizing
the two-channel Kondo (2CK) effect as suggested by Cox [8]. In the dilute
limit of a few Pr3+ ions the interaction of each pseudospin with the itinerant
electrons might then be governed at low temperature by the non-Fermi liquid
fixed point of the two-channel Kondo model characterized by a residual entropy
of $\frac{k_{B}}{2}\log 2$ per ion. Close to this fixed-point, the quadrupolar
pseudospin susceptibility $\chi_{Q}$, that characterizes the response of the
pseudospin to an infinitesimally small pseudomagnetic field
$\boldsymbol{H}_{\Gamma_{3}}$, is logarithmically divergent with temperature,
$\chi_{Q}\sim\log\frac{1}{T}$. A pseudomagnetic field is generated when the
quadrupolar symmetry of the Pr3+ environment gets broken. This is, for
example, achieved by either a suitably applied elastic strain or by the
application of a magnetic field, where $H_{\Gamma_{3}}\sim H^{2}$ for small
fields $H$.

###  II.3 Elastic coupling to the local quadrupole moment

The Stevens operators can locally couple to the strain tensor of the crystal
lattice $\varepsilon_{ij}(\vec{r})$,

$\displaystyle\mathcal{H}_{\rm strain}=\int
d\vec{r}\frac{1}{2}\varepsilon_{ij}(\vec{r})C_{ijkl}\varepsilon_{kl}(\vec{r})-g_{\Gamma_{3}}\sum_{\alpha}\left[O^{0}_{2,\alpha}\varepsilon_{u}(\vec{r}_{\alpha})+O_{2,\alpha}^{2}\varepsilon_{\mathrm{v}}(\vec{r}_{\alpha})\right]-g_{\Gamma_{1}}\sum_{\alpha}O^{0}_{0,\alpha}\varepsilon_{B}(\vec{r}_{\alpha})$
(S4)

where $C_{ijkl}$ is the elastic constant matrix and $O^{q}_{k,\alpha}$ are the
Stevens operator associated with the Pr3+ ion enumerated by the index $\alpha$
and located at position $\vec{r}_{\alpha}$. We have assumed for simplicity
that the elastic coupling constants, $g_{\Gamma_{3}}$ and $g_{\Gamma_{1}}$,
are the same for all $\alpha$. The Stevens operators $O^{0}_{2,\alpha}$ and
$O^{2}_{2,\alpha}$ generate a local quadrupolar stress, whereas
$O^{0}_{0,\alpha}=1$ generates a local isotropic stress. The crystal responds
to the local stresses leading to a homogeneous and inhomogeneous strain
distribution.

####  II.3.1 Macroscopic strains induced by the Pr ions

In the following, we focus on the homogeneous macroscopic strain denoted by
$\varepsilon_{ij}=\frac{1}{V}\int d\vec{r}\varepsilon_{ij}(\vec{r})$. The
macroscopic eigenstrains and the corresponding elastic constants in a cubic
system are listed in Table S1 where we adopt the convention of the previous
literature on Y1-xPrxIr2Zn20.

Symmetry | Symmetry Strains | Symmetry Elastic Constants
---|---|---
$\Gamma_{1}$ | $\varepsilon_{\mathrm{B}}=\varepsilon_{xx}+\varepsilon_{yy}+\varepsilon_{zz}$ | $(c_{11}+2c_{12})/3$
$\Gamma_{3}$ | $\varepsilon_{\mathrm{u}}=(2\varepsilon_{zz}-\varepsilon_{xx}-\varepsilon_{yy})/\sqrt{3}$ | $(c_{11}-c_{12})/2$
| $\varepsilon_{\mathrm{v}}=\varepsilon_{xx}-\varepsilon_{yy}$ | $(c_{11}-c_{12})/2$
$\Gamma_{5}$ | $\varepsilon_{xy},\varepsilon_{xz},\varepsilon_{yz}$ | $4c_{44}$
Table S1: Symmetries, respective symmetry strains and symmetry elastic
constants for cubic symmetry.

Minimizing Eq. (S4) with respect to the macroscopic strains we obtain

$\displaystyle\varepsilon_{\mathrm{B}}$
$\displaystyle=\frac{g_{\Gamma_{1}}n_{\rm
Pr}}{(c^{0}_{11}+2c^{0}_{12})/3},\qquad\varepsilon_{\mathrm{u}}=\frac{g_{\Gamma_{3}}}{(c^{0}_{11}-c^{0}_{12})/2}\frac{1}{V}\sum_{\alpha}\left\langle
O^{0}_{2,\alpha}\right\rangle,\qquad\varepsilon_{\mathrm{v}}=\frac{g_{\Gamma_{3}}}{(c^{0}_{11}-c^{0}_{12})/2}\frac{1}{V}\sum_{\alpha}\left\langle
O^{2}_{2,\alpha}\right\rangle,$ (S5)

where $n_{\mathrm{Pr}}$ is the Pr density. The macroscopic strain induced by
the Pr3+ ions is related to the thermal expectation values $\langle
O^{q}_{k,\alpha}\rangle$, and $(c^{0}_{11}+2c^{0}_{12})/3$ and
$(c^{0}_{11}-c^{0}_{12})/2$ are the background elastic constants with
$\Gamma_{1}$ and $\Gamma_{3}$ symmetry, respectively. Whereas the expectation
value $\left\langle O^{0}_{0,\alpha}\right\rangle=1$ is trivial, the ones in
the $\Gamma_{3}$ channel are governed by the thermal occupation of the non-
Kramers doublet. Importantly, these expectation values $\langle
O^{q}_{2,\alpha}\rangle$ will in general depend themselves on both elastic
coupling $g_{\Gamma_{3}}$ and the Kondo coupling $J_{K}$. Moreover, they will
also depend on the local strain distribution at the site $\vec{r}_{\alpha}$.
Only after averaging over this strain disorder one obtains

$\displaystyle\varepsilon_{\mathrm{B}}$
$\displaystyle\approx\frac{n_{\mathrm{Pr}}g_{\Gamma_{1}}}{(c^{0}_{11}+2c^{0}_{12})/3}\langle
O^{0}_{0}\rangle,\qquad\varepsilon_{\mathrm{u}}\approx\frac{n_{\mathrm{Pr}}{g_{\Gamma_{3}}}}{(c^{0}_{11}-c^{0}_{12})/2}\langle
O^{0}_{2}\rangle,\qquad\varepsilon_{\mathrm{v}}\approx\frac{n_{\mathrm{Pr}}{g_{\Gamma_{3}}}}{(c^{0}_{11}-c^{0}_{12})/2}\langle
O^{2}_{2}\rangle,$ (S6)

with the Pr density $n_{\mathrm{Pr}}$.

In general, a macroscopic strain along a certain direction specified by the
unit vector $\hat{n}$ is given by $\hat{n}_{i}\varepsilon_{ij}\hat{n}_{j}$. In
particular, the strains along the cubic directions [100], [010] and [001] as
well as along the diagonal [111] are related to the eigenstrains via

$\displaystyle\frac{\Delta L}{L}\bigg{\rvert}_{[001]}$
$\displaystyle=\varepsilon_{zz}=\frac{1}{3}\varepsilon_{\mathrm{B}}+\frac{1}{\sqrt{3}}\varepsilon_{\mathrm{u}},$
(S7) $\displaystyle\frac{\Delta L}{L}\bigg{\rvert}_{[100]}$
$\displaystyle=\varepsilon_{xx}=\frac{1}{3}\varepsilon_{\mathrm{B}}-\frac{1}{2\sqrt{3}}\varepsilon_{\mathrm{u}}+\frac{1}{2}\varepsilon_{\mathrm{v}},$
(S8) $\displaystyle\frac{\Delta L}{L}\bigg{\rvert}_{[010]}$
$\displaystyle=\varepsilon_{yy}=\frac{1}{3}\varepsilon_{\mathrm{B}}-\frac{1}{2\sqrt{3}}\varepsilon_{\mathrm{u}}-\frac{1}{2}\varepsilon_{\mathrm{v}},$
(S9) $\displaystyle\frac{\Delta L}{L}\bigg{\rvert}_{[111]}$
$\displaystyle=\frac{1}{3}\varepsilon_{\mathrm{B}}+\frac{2}{3}\Big{(}\varepsilon_{xy}+\varepsilon_{zx}+\varepsilon_{yz}\Big{)}.$
(S10)

We singled out the [001] direction along which the magnetic field will be
applied. In the absence of any internal or external strains that break
additional symmetries, we expect that the shear strains vanish
$\varepsilon_{xy}=\varepsilon_{zx}=\varepsilon_{yz}=0$. In this case the
strain along [111] is only sensitive to the bulk strain $\varepsilon_{B}$. In
case the cubic symmetry is maintained for zero magnetic field, the
$\Gamma_{3}$ doublet remains degenerate and $\langle O^{q}_{2}\rangle=0$. In
turn, this implies vanishing quadrupolar strains
$\varepsilon_{\mathrm{u}}=\varepsilon_{\mathrm{v}}=0$, so that only the
isotropic bulk strain $\varepsilon_{\mathrm{B}}$ will contribute to a finite
expansion. A magnetic field along $\boldsymbol{H}\parallel[001]$ leads to a
tetragonal distortion via the elastic coupling resulting in a finite value for
$\varepsilon_{\mathrm{u}}$. The strain $\varepsilon_{\mathrm{v}}$ still
remains zero implying equivalent strains $\varepsilon_{xx}$ and
$\varepsilon_{yy}$ within the plane perpendicular to the applied magnetic
field.

Evaluating the thermal average in Eqs. (S6) with the CEF Hamiltonian Eq. (S1)
in zeroth order in the Kondo coupling $J_{K}$ and treating the elastic
coupling $g_{\Gamma_{3}}$ in Eq. (S4) in a mean-field approximation, one
obtains a characteristic dependence of the uniaxial strains on temperature and
magnetic field that are shown as dashed lines in Fig. 1 of the main text.
Here, we used $(c_{11}^{0}-c_{12}^{0})/2=52.771$ GPa [17] and
$N_{\mathrm{Pr}}=0.101\times 10^{27}$ m-3.

In Fig. S2 the CEF prediction for the quadrupolar thermal expansion,
$\alpha_{u}=\partial_{T}\varepsilon_{u}$, is shown indicating that at small
fields its temperature dependence is governed by the Curie-Weiss behavior,
$\alpha_{u}/H^{2}\sim\partial_{T}\chi_{Q}\sim 1/T^{2}$.

Figure S2: Quadrupolar thermal expansion
$\alpha_{u}=\partial_{T}\varepsilon_{u}$ divided by $(\mu_{0}H)^{2}$ versus
temperature as predicted by the crystal field Hamiltonian of Eq. (S1). The
simulation parameters are the ones of single crystal with Y1-xPrxIr2Zn20
$(x=0.036$).

####  II.3.2 Influence of local strains

In the presence of finite local strains induced by crystal defects or
inhomogeneities, the symmetry of the local Pr environment will be broken which
splits the non-Kramers doublet. This results in a finite expectation value of
the quadrupole moments, that for small strains is given by $\langle
O^{0}_{2,\alpha}\rangle\approx\chi_{Q}g_{\Gamma_{3}}Q^{u}_{ij}\varepsilon^{\rm
loc}_{ij}(\vec{r}_{\alpha})$ and $\langle
O^{2}_{2,\alpha}\rangle\approx\chi_{Q}g_{\Gamma_{3}}Q^{v}_{ij}\varepsilon^{\rm
loc}_{ij}(\vec{r}_{\alpha})$. Using Eq. (S5) we then obtain a contribution to
the macroscopic quadrupolar strains

$\displaystyle\varepsilon_{u/v}=\frac{g_{\Gamma_{3}}\chi_{Q}g_{\Gamma_{3}}Q^{u/v}_{ij}}{(c^{0}_{11}-c^{0}_{12})/2}\frac{1}{V}\sum_{\alpha}\varepsilon^{\rm
loc}_{ij}(\vec{r}_{\alpha})$ (S11)

where we abbreviated
$Q^{v}_{ij}=\delta_{ix}\delta_{jx}-\delta_{iy}\delta_{jy}$, and
$Q^{u}_{ij}=\frac{1}{\sqrt{3}}(2\delta_{iz}\delta_{jz}-\delta_{ix}\delta_{jx}-\delta_{iy}\delta_{jy})$.
After disorder averaging, this contribution vanishes by definition because
$\int d\vec{r}\varepsilon^{\rm loc}_{ij}(\vec{r})=0$. However, for a given
strain distribution this might yield a finite contribution that scales
linearly with the Pr density $\varepsilon_{u/v}\sim n_{\rm Pr}\chi_{Q}$ and
breaks the cubic crystal symmetry even at zero magnetic field.

The Pr ions themselves will also give rise to local strain fields. In general,
the local strains, $\varepsilon^{\rm
loc}_{ij}=\frac{1}{2}(\partial_{i}u_{j}+\partial_{j}u_{i})$, are represented
in terms of the displacement vector $u_{i}$. Minimizing Eq. (S4) with respect
to $u_{i}$ we obtain for the displacement vector in Fourier space
$u_{i}(\vec{k})$

$\displaystyle-
C_{ijkl}k_{j}k_{k}u_{l}(\vec{k})=\sum_{\alpha}e^{-i\vec{k}\vec{r}_{\alpha}}ik_{j}\Big{[}g_{\Gamma_{3}}Q^{v}_{ij}\langle
O^{0}_{2,\alpha}\rangle+g_{\Gamma_{3}}Q^{u}_{ij}\langle
O^{2}_{2,\alpha}\rangle+g_{\Gamma_{1}}\delta_{ij}\Big{]}$ (S12)

where we used that $\langle O^{0}_{0,\alpha}\rangle=1$. In case the last term
dominates and the quadrupolar expectation values $\langle
O^{q}_{2,\alpha}\rangle$ can be neglected, the Pr3+ ions induce a local strain
field given by

$\displaystyle
u_{l}(\vec{k})=-ig_{\Gamma_{1}}D^{-1}_{ln}(\vec{k})k_{n}\sum_{\alpha}e^{-i\vec{k}\vec{r}_{\alpha}}$
(S13)

where we introduced the dynamical matrix $D_{il}(\vec{k})=C_{ijkl}k_{j}k_{k}$.
The resulting local strain field is given by

$\displaystyle\varepsilon^{\rm
loc}_{ij}(\vec{r})=g_{\Gamma_{1}}\int\frac{d\vec{k}}{(2\pi)^{3}}\mathcal{K}_{ij}(\vec{k})\sum_{\alpha}e^{i\vec{k}(\vec{r}-\vec{r}_{\alpha})}=g_{\Gamma_{1}}\sum_{\alpha}\tilde{\mathcal{K}}_{ij}(\vec{r}-\vec{r}_{\alpha})$
(S14)

with the kernel
$\mathcal{K}_{ij}(\vec{k})=\frac{1}{2}(k_{i}D^{-1}_{jn}(\vec{k})k_{n}+k_{j}D^{-1}_{in}(\vec{k})k_{n})$
and its Fourier transform $\tilde{\mathcal{K}}_{ij}$. It captures the long-
range influence of the Pr impurities that is mediated by static strain fields.

Plugging this into Eq. (S11) we then obtain for the macroscopic quadrupolar
strain

$\displaystyle\varepsilon_{u/v}=\frac{g_{\Gamma_{1}}g_{\Gamma_{3}}^{2}\chi_{Q}}{(c^{0}_{11}-c^{0}_{12})/2}\frac{1}{V}\sum_{\alpha,\beta}Q^{u/v}_{ij}\tilde{\mathcal{K}}_{ij}(\vec{r}_{\alpha}-\vec{r}_{\beta})=n_{\rm
Pr}\frac{g_{\Gamma_{1}}g^{2}_{\Gamma_{3}}\chi_{Q}}{(c^{0}_{11}-c^{0}_{12})/2}\int\frac{d\vec{k}}{(2\pi)^{3}}Q^{u/v}_{ij}\mathcal{K}_{ij}(\vec{k})|z(\vec{k})|^{2}.$
(S15)

where we abbreviated
$z(\vec{k})=\frac{1}{\sqrt{N_{\mathrm{Pr}}}}\sum_{\alpha}e^{i\vec{k}\vec{r}_{\alpha}}$
with $N_{\mathrm{Pr}}$ being the number of Pr ions.

For wavevectors $|\vec{k}|L\gg 1$ much larger than the inverse linear system
size and uniformly distributed positions of Pr ions, $z(\vec{k})$ is a random
complex number with a phase uniformly distributed on $[0,2\pi)$. After
disorder averaging $|z|^{2}\to 1$, the above expression vanishes because the
integral over wavevectors contains the quadrupolar matrix $Q^{u/v}_{ij}$. For
a specific disorder configuration, however, the integral might be finite and
for dimensional reasons it will scale as $1/a_{\mathrm{Pr}}^{3}$ where
$a_{\mathrm{Pr}}$ is a characteristic length scale, that is expected to obey
$a_{\mathrm{Pr}}^{3}\sim 1/n_{\mathrm{Pr}}$. This suggests that the static
local strains induced by Pr ions give rise to a macroscopic quadrupolar strain
$\varepsilon_{u/v}\sim n_{\mathrm{Pr}}^{2}\chi_{Q}$ but with a prefactor
suppressed by two powers of the Pr density. For this reasons, the local
strains induced by the Pr ions themselves represents a sub-leading correction
in the dilute limit.

These disorder induced corrections, that break the cubic symmetry even in zero
field, could in principle lead to systematic errors in our data analysis. In
particular, care must be taken when identifying the bulk thermal expansion
from uniaxial thermal expansion measurements along a single cubic direction
$\langle 100\rangle$ at zero field, see the discussion of Fig. 4 in the main
text. Not only a small external stress from the sample holder but also the
disorder might break the cubic symmetry inducing a finite quadrupolar
macroscopic strain with an anomalous $T$-dependence induced by the response of
the non-Kramers doublet. However, thermal expansion measurements along
$\langle 111\rangle$ will not be affected by $\varepsilon_{u/v}$, see Eq.
(S10). Disorder might induce a finite shear, e.g. $\varepsilon_{xy}$, but
without coupling to the non-Kramers doublet and, therefore, without the
anomalous $T$-dependence due to the susceptibility $\chi_{Q}$.

###  II.4 Anomalous bulk thermal expansion induced by quadrupolar Kondo
correlations

In the following subsections, we focus on the limit of both vanishing applied
stress and magnetic fields. The linear coupling of the strain to the
quadrupolar moment realizes an effective spin-boson model, and it is expected
that this interaction eventually quenches the residual entropy of the 2CK
fixed-point. There are two candidates for this quenching process: the
degeneracy of the non-Kramers doublet is either lifted by dynamical strain
fields, i.e., by phonons, or by static local strains. Both lead to a
perturbative correction to the free energy that resembles a high-temperature
tail of a Schottky anomaly, however, here for the quench of the anomalous
$\frac{k_{B}}{2}\log 2$ entropy.

This process might be influenced by the application of a hydrostatic pressure.
We will discuss possible thermodynamic signatures in the bulk thermal
expansion and the specific heat. Whereas both quantities are expected to show
an anomalous $T$-dependence, Grüneisen scaling is actually expected. Possible
reasons for a large Grüneisen parameter is discussed.

####  II.4.1 Renormalization of the elastic constant

Integrating out the electronic quadrupolar degrees of freedom, one obtains a
correction to the Hamiltonian in second order in the elastic coupling that
reads for small strains

$\displaystyle\delta\mathcal{H}_{\rm
strain}=-\frac{g^{2}_{\Gamma_{3}}\chi_{Q}}{2}\sum_{\alpha}\Big{(}\varepsilon_{u}(\vec{r}_{\alpha})^{2}+\varepsilon_{v}(\vec{r}_{\alpha})^{2}\Big{)}$
(S16)

where we assumed that the quadrupolar moments of distinct Pr3+ ions are
uncorrelated. After disorder averaging,

$\displaystyle\overline{\delta\mathcal{H}_{\rm strain}}=\int
d\vec{r}\Big{(}-\frac{g^{2}_{\Gamma_{3}}n_{\rm
Pr}\chi_{Q}}{2}\Big{(}\varepsilon_{u}(\vec{r})^{2}+\varepsilon_{v}(\vec{r})^{2}\Big{)}\Big{)}$
(S17)

this correction reduces to a renormalization of the elastic constant

$\displaystyle\frac{c_{11}-c_{12}}{2}=\frac{c^{0}_{11}-c^{0}_{12}}{2}-g^{2}_{\Gamma_{3}}n_{\rm
Pr}\chi_{Q}.$ (S18)

The disorder averaging should be justified for strains that vary on a length
scale that is much larger than the distance between Pr3+ ions. As the
susceptibility $\chi_{Q}\sim\log\frac{1}{T}$ diverges for low temperatures at
the 2CK fixed-point, the renormalized elastic constant is reduced and might
eventually vanishes at ultralow temperatures indicating the tendency towards a
structural instability. The logarithmic $T$-dependence of this elastic
constant has been confirmed by ultrasound experiments [17].

####  II.4.2 Correction to the thermodynamics

The effective Hamiltonian (S16) will lead to a correction to the free energy
density

$\displaystyle\delta f=-\frac{1}{2}g^{2}_{\Gamma_{3}}n_{\rm
Pr}\chi_{Q}\langle\varepsilon^{2}_{E}\rangle_{0}$ (S19)

The finite expectation value $\langle\varepsilon^{2}_{E}\rangle_{0}$ might
possess two different origins. First, it might arise from dynamic strain
fields, i.e., phonon excitations. For long-wavelength phonons we can use the
disorder averaged Hamiltonian (S17) and
$\langle\varepsilon^{2}_{E}\rangle_{0}\equiv\frac{1}{V}\int
d\vec{r}\langle\varepsilon_{u}(\vec{r})^{2}+\varepsilon_{v}(\vec{r})^{2}\rangle_{0}$.
At lowest temperatures, this expectation value then corresponds to the zero-
point fluctuations of phonons evaluated with the bare elastic constants.
Second, it might arise from static local strain fields,
$\langle\varepsilon^{2}_{E}\rangle_{0}=\frac{1}{N_{\rm
Pr}}\sum_{\alpha}\Big{(}\varepsilon_{u}(\vec{r}_{\alpha})^{2}+\varepsilon_{v}(\vec{r}_{\alpha})^{2}\Big{)}$
induced by crystal defects and inhomogeneities, where $N_{\rm Pr}$ is the
number of Pr ions. In both cases, this yields a negative correction to the
entropy

$\displaystyle\delta S=-\partial_{T}\delta
f=\frac{1}{2}g^{2}_{\Gamma_{3}}n_{\rm
Pr}\langle\varepsilon^{2}_{E}\rangle_{0}\partial_{T}\chi_{Q}=-\frac{1}{2}g^{2}_{\Gamma_{3}}n_{\rm
Pr}\langle\varepsilon^{2}_{E}\rangle_{0}A_{Q}\frac{1}{T}.$ (S20)

where we used $\chi_{Q}\approx A_{Q}\log\frac{1}{T}$ with a constant $A_{Q}$.
This should be compared to the $-1/T^{2}$ behavior of the high-temperature
tail of a standard Schottky anomaly. This entropy entails a singular
correction to both the specific heat $\delta C=T\partial_{T}\delta S$ and the
volume thermal expansion $\delta\beta=-\partial_{p}\delta S$ as a function of
$T$,

$\displaystyle\delta C$ $\displaystyle=\frac{1}{2}g^{2}_{\Gamma_{3}}n_{\rm
Pr}\langle\varepsilon^{2}_{E}\rangle_{0}A_{Q}\frac{1}{T},\qquad\delta\beta=\frac{1}{2}\Big{[}\partial_{p}\Big{(}g^{2}_{\Gamma_{3}}n_{\rm
Pr}\langle\varepsilon^{2}_{E}\rangle_{0}A_{Q}\Big{)}\Big{]}\frac{1}{T}$ (S21)

where $\partial_{p}$ is the derivative with respect to hydrostatic pressure.
The ratio of both yield a reduced Grüneisen parameter

$\displaystyle\delta\Gamma=\frac{\delta\beta}{\delta
C}=\partial_{p}\log(g^{2}_{\Gamma_{3}}n_{\rm
Pr}\langle\varepsilon^{2}_{E}\rangle_{0}A_{Q})$ (S22)

that is determined by the pressure dependence of the product of characteristic
parameters $g^{2}_{\Gamma_{3}}n_{\rm
Pr}\langle\varepsilon^{2}_{E}\rangle_{0}A_{Q}$. In case that, e.g., the main
$p$-dependence arises from the elastic coupling $g_{\Gamma_{3}}$, the reduced
Grüneisen parameter reduces to $\delta\Gamma=\partial_{p}\log
g^{2}_{\Gamma_{3}}$.

##  III Analysis of the experimental data

###  III.1 Thermal expansion and Grüneisen parameter for
$\boldsymbol{H\parallel[001]}$

This section discusses the longitudinal and transverse thermal expansion
coefficients that were measured for $\boldsymbol{H}\parallel[001]$. Based on
these two measurements, the volume and quadrupolar thermal expansion
coefficients can be deduced. In addition, the specific heat data, used for the
calculation of the quadrupolar Grüneisen ratio, is presented.

Figure S3: Temperature dependence of the (a) longitudinal thermal expansion
coefficient $\alpha_{\parallel}$ and (b) transverse thermal expansion
coefficient $\alpha_{\perp}$ together with CEF calculations for
$\boldsymbol{H}\parallel[001]$. (c) Temperature dependence of the volume
thermal expansion coefficient $\beta=\alpha_{\parallel}+2\alpha_{\perp}$ at
different magnetic fields. The inset shows the temperature dependence of
$\beta=3\alpha_{[001]}$ normalized to the Pr-concentration $x$ for three
differently doped samples at zero magnetic field.

####  III.1.1 Longitudinal, transverse, volume and quadrupolar thermal
expansion

In order to deduce the quadrupolar and volume thermal expansion coefficients
$\alpha_{\mathrm{u}}$ and $\beta$, which are shown in Fig. 1(a) and Fig. 4(a)
of the main text, longitudinal $(\alpha_{\parallel}=\alpha_{[001]}\parallel
H)$ and transverse $(\alpha_{\perp}=\alpha_{[100],[010]}\parallel H)$ thermal
expansion measurements were carried out for $\boldsymbol{H}\parallel[001]$.

Figure S3(a) and (b) shows the longitudinal ($\alpha_{\parallel}$) and
transverse thermal expansion coefficients ($\alpha_{\perp}$) at various
magnetic fields together with CEF simulations, which were calculated by
following the approach explicitly discussed in section II. The quadrupole-
strain coupling constant $g_{\Gamma_{3}}$ was determined at -28.9 K by fitting
the simulated data to the high field data at 10 T where the system is expected
to be in the fully localized state. This value is comparable to
$g_{\mathrm{\Gamma_{3}}}=-38.0$ K [25] of pure PrIr2Zn20.

By using the longitudinal and transverse thermal expansion coefficients, the
volume thermal expansion coefficient, as shown in Fig. S3(c), is obtained

$\beta=\alpha_{\parallel}+2\alpha_{\perp}.$ (S23)

By using $\alpha_{\parallel}$ and $\beta$, the quadrupolar thermal expansion
coefficient $\alpha_{\mathrm{u}}$ then calculates with the help of Eq. (S7) as

$\alpha_{\mathrm{u}}=\sqrt{3}\left(\alpha_{\parallel}-\frac{1}{3}\beta\right),$
(S24)

which is presented in Fig. 1(a) of the main text.

####  III.1.2 Grüneisen Parameter

In order to calculated the quadrupolar Grüneisen parameter plotted in Fig.
1(b) of the main text, the specific heat data measured on a Y1-xPrxIr2Zn20
single crystal from the same batch as the one used for the thermal expansion
measurement, having a comparable Pr-concentration of $x=0.044$, was used. The
specific heat of this particular sample has already been published in Refs.
[5, 23]. Subsequently shown data was partially taken from Refs. [5, 23]
($T\gtrsim 0.8)$ and partially remeasured with a different setup to obtain
higher resolved data down to lower temperatures ($T<0.8$ K). Figure S4 shows
the 4f contribution to the molar specific heat as a function of temperature of
Y1-xPrxIr2Zn20 ($x=0.044$) at various magnetic fields
$\boldsymbol{H}\parallel[001]$.

Figure S4: Temperature dependence of the 4f contribution to the molar specific
heat $C_{\mathrm{m}}$ at magnetic fields $H\parallel[001]$ up to 4 T. The
specific heat of this particular single crystal has already been shown in Ref.
[5, 23] (except data at 1 T). The here presented data was partially taken from
Refs. [5, 23] ($T\geq 0.8$ K) and partially remeasured ($T<0.8$ K).

In zero magnetic field a quadrupolar nuclear contribution and a lattice
contribution, determined from a reference measurement on YIr2Zn20, were
subtracted. For the data obtained in magnetic field, a magnetic nuclear
Schottky contribution was subtracted additionally. As outlined previously by
Yamane et al. [5], this particular sample shows a behavior $C/T\sim\log 1/T$
in zero magnetic field consistent with the single-impurity quadrupolar Kondo
scenario.

The molar volume, used for the calculation of the Grüneisen parameter, was
determined by using the specific heat sample’s lattice constant
$a=14.197\times 10^{-10}$ m and the number of formula units per unit volume
$Z=8$, as

$V_{\mathrm{m}}(x=0.044)=\frac{V_{\mathrm{m}}\times a^{3}}{Z}=2.154024\times
10^{-4}\,\frac{\mathrm{m}^{3}}{\mathrm{mol}}$ (S25)

As no specific heat measurement was performed for $\mu_{0}H=0.5$ T, the
specific heat data at zero magnetic field was used to calculate the
quadrupolar Grüneisen parameter at 0.5 T. Fig. S4 shows that the specific heat
at 0 T and 1 T are nearly identical, which suggests that $C_{\mathrm{m}}$ can
be considered as nearly field independent for $\mu_{0}H\leq 1$ T. This
justifies the use of the zero field specific heat for the calculation of
$\Gamma_{\mathrm{u}}$ at 0.5 T. For all other $\Gamma_{\mathrm{u}}$ curves
shown in Fig. 1(b) in the main text, thermal expansion and specific heat were
measured at the same magnetic field. As pointed out in Section I.A., by taking
the error in the Pr content of the highly diluted samples into account, the Pr
content of the two single crystals employed for the dilatometry ($x=0.033$ and
$x=0.036$) can be considered as the same as the Pr content of the single
crystal employed for the specific heat measurement ($x=0.044$). Consequently,
a normalization of the thermal expansion and the specific heat to the Pr
content was not necessary in order to calculate the Grüneisen parameter.

###  III.2 Magnetostriction for $\boldsymbol{H\parallel[001]}$

In this section we present the longitudinal and transverse magnetostriction
coefficients measured for $\boldsymbol{H}\parallel[001]$ on highly diluted
Y1-xPrxIr2Zn20 $(x=0.036)$. From these measurements, the volume
magnetostriction and the quadrupolar magnetostriction can be determined.
Finally, we discuss how the low temperature magnetostriction can be employed
to estimate the quadrupole susceptibility.

Figure S5: (a) Longitudinal and transverse magnetostriction coefficients for
$\boldsymbol{H}\parallel[001]$, measured on the highly diluted Y1-xPrxIr2Zn20
single crystal with $x=0.036$. (b) Volume magnetostriction determined from the
longitudinal and transverse magnetostriction shown in (a) via the relation
$\Delta V/V=\Delta L/L_{\parallel}+2\Delta L/L_{\perp}$.

####  III.2.1 Longitudinal, transverse, volume and quadrupolar
magnetostriction

To determine the volume and quadrupolar magnetostriction coefficients
$\varepsilon_{\mathrm{B}}$ and $\varepsilon_{\mathrm{u}}$ of the the highly
diluted Y1-xPrxIr2Zn20, longitudinal ($(\Delta
L/L)_{\parallel}=\varepsilon_{\parallel}=\varepsilon_{[001]}\parallel H$) and
transverse ($(\Delta
L/L)_{\perp}=\varepsilon_{\perp}=\varepsilon_{[100],[010]}\perp H$)
magnetostriction measurements were performed which are shown in Fig. S5. CEF
calculation are shown as dashed lines. The curves were calculated by use of
the same CEF parameter $W=-1.50$ K and $x=0.537$ by which also the thermal
expansion was simulated. For the background elastic modulus the value reported
for highly diluted Y1-xPrxIr2Zn20 in the Supplemental material of Ref. [17],
namely $(C_{11}-C_{12})/2=52.771$ GPa. The magnetostriction analysis suggests
a slightly larger quadrupole strain coupling constant of
$g_{\Gamma_{3}}=-35.0$ K than the thermal expansion ($g_{\Gamma_{3}}=-28.9$
K). The Pr density per unit volume is $N_{\mathrm{Pr}}(x=0.036)=0.101\times
10^{27}$ m-3. $\varepsilon_{\parallel}$ and $\varepsilon_{\perp}$ allow for
the calculation of the bulk strain $\varepsilon_{\mathrm{B}}=\Delta V/V$ via
the relation

$\varepsilon_{\mathrm{B}}=\varepsilon_{\parallel}+2\varepsilon_{\perp}.$ (S26)

The bulk strain determined in this manner is shown in Fig. S5(b) as a function
of magnetic field $\boldsymbol{H}\parallel[001]$ for various temperatures. As
compared to the longitudinal and transverse magnetostriction, the volume
change is markedly small, which is in line with the previous results on pure
PrIr2Zn20 [25]. At the lowest measured temperatures, $(\Delta V/V)$ exhibits a
small magnetic field dependence which is gradually suppressed as temperature
increases. Thereby, $(\Delta V/V)$ changes the most at intermediate field
ranging from 1 T to 3 T, which can be likely connected to the suppression of
the divergent volume thermal expansion in magnetic field. The field variation
of $(\Delta V/V)$ at different temperatures is similar to the one of the
previously reported ($c_{11}-c_{12}$)/2 elastic constant [17], providing
further indication for a possible coupling between volume strain and the
quadrupole susceptibility via static or dynamic strain fields.

The quadrupolar magnetostriction $\varepsilon_{\mathrm{u}}$, which is shown in
Fig. 3 of the main text, can simply be calculated by using Eq. (S7) as

$\varepsilon_{\mathrm{u}}=\sqrt{3}\left(\varepsilon_{\parallel}-\frac{\varepsilon_{\mathrm{B}}}{3}\right).$
(S27)

####  III.2.2 Quadrupole susceptibility

Next, we briefly discuss the estimation of the quadrupole susceptibility from
the magnetostriction data. As outlined in the previous section, the
expectation value of the quadrupole operator $\left\langle
O^{0}_{2}\right\rangle$ depends for small $H$ quadratically on magnetic field

$\left\langle O^{0}_{2}\right\rangle=\chi^{F}_{\mathrm{Q}}(\mu_{0}H)^{2},$
(S28)

where $\chi^{F}_{Q}\sim\chi_{Q}$. This initial quadratic field dependence is
thereby mainly generated by the mixing of $\Gamma_{3}$ ground state and the
first excited $\Gamma_{4}$ state in magnetic field. Plugging this into Eq.
(S6) one obtains

$\varepsilon_{\mathrm{u}}=\frac{N_{\mathrm{Pr}}{g_{\Gamma_{3}}}}{(c^{0}_{11}-c^{0}_{12})/2}\chi_{\mathrm{Q}}(\mu_{0}H)^{2}.$
(S29)

This suggests that the initial quadratic field dependence of
$\varepsilon_{\mathrm{u}}$ measured at various temperatures, provides direct
access to the quadrupole susceptibility $\chi_{\mathrm{Q}}$. Fig. S6 shows
$\varepsilon_{\mathrm{u}}$ versus $(\mu_{0}H)^{2}$.

Figure S6: $\varepsilon_{\mathrm{u}}$ as a function of $(\mu_{0}H)^{2}$ at
different temperatures to illustrate the quadratic field dependence of
$\epsilon_{\mathrm{u}}$ for small $\boldsymbol{H}\parallel[001]$.

By fitting the low field magnetostriction with a linear regression that yields
the initial slope $m$, we can then determine $\chi_{\mathrm{Q}}$ at selected
temperatures. For $\chi_{\mathrm{Q}}$ it follows

$\chi_{\mathrm{Q}}=\frac{(c_{11}^{0}-c_{12}^{0})/2}{N_{\mathrm{Pr}}g_{\Gamma_{3}}}m$
(S30)

In the inset of Fig. 3 of the main text, we also present $\chi_{\mathrm{Q}}$
of two higher doped Y1-xPrxIr2Zn20 single crystals with $x=0.09$ and $x=0.49$.
As for these two samples the magnetostriction was only measured parallel to
magnetic field, $\chi_{\mathrm{Q}}$ is determined by evaluating the strain
parallel to magnetic field

$\varepsilon_{\parallel}=\frac{1}{\sqrt{3}}\frac{N_{\mathrm{Pr}}{g_{\Gamma_{3}}}}{(C^{0}_{11}-C^{0}_{12})/2}\chi_{\mathrm{Q}}^{\mathrm{F}}H^{2},$
(S31)

under the assumption that $\varepsilon_{\mathrm{v}}=0$ and that the field
dependent contribution of the isotropic bulk strain $\varepsilon_{\mathrm{B}}$
is small as compared to the uniaxial quadrupolar contribution
$\varepsilon_{\mathrm{u}}$. The former assumption is justified, as
$\epsilon_{\mathrm{u}}\gg\epsilon_{\mathrm{B}}$ for both the highly diluted
system with $x=0.036$ and the pure system PrIr2Zn20 [25].