# Phase Separation Driven by Active Flows

Saraswat Bhattacharyya, Julia M. Yeomans  Rudolf Peierls Centre For
Theoretical Physics, University of Oxford, Oxford OX1 3PU, United Kingdom.

###### Abstract

We extend the continuum theories of active nematohydrodynamics to model a two-
fluid mixture with separate velocity fields for each fluid component, coupled
through a viscous drag. The model is used to study an active nematic fluid,
mixed with an isotropic fluid. We find micro-phase separation, and argue that
this results from an interplay between active anchoring and active flows
driven by concentration gradients. The results may be relevant to cell-sorting
and the formation of lipid rafts in cell membranes.

Introduction.—There is increasing evidence that phase ordering plays a
fundamental role in biological processes such as morphogenesis and colony
growth [1, 2, 3, 4]. Phase ordering of cell types (Fig. 1) [5], membrane-less
organelles [6] and RNA-protein mixtures [7] inside cells have been studied
using free-energy prescriptions. However, biological matter is intrinsically
out of thermodynamic equilibrium, and it is unlikely that biological phase
ordering relies entirely on the framework of equilibrium thermodynamics.

Figure 1: (color online.) Cell Sorting: (a), (b), (c) starfish embryo at times
(a) t= 4 hr, (b) t= 8 hr, (c) t = 48 hr showing sorting of embryonic cells
(adapted from Ref [8]). (d) Phase ordering of extensile (magenta) and
contractile (green) epithelial MDCK cells (courtesy of L. Balasubramaniam
[9]).

Self-motile particles form an important class of non-equilibrium systems
called active matter [10]. Active matter exhibits non-equilibrium phase
segregation mechanisms, such as motility-induced phase separation (MIPS) [11,
12, 13] where active particles get trapped in regions of high particle
density. Phase separation has been found in continuum models [14, 15, 16], in
the presence of hydrodynamic interactions [17], or driven by aligning torques
[18, 19]. Active nematics consist of active rod-like particles with
orientational order. The elongated particles generate active dipolar stresses
which destabilise the aligned state and lead to active turbulence
characterised by a chaotic and highly vortical velocity field [20, 21, 22].
Recent work has shown that bacterial biofilms [23, 24], epithelial tissues
[25, 26], and microtubule-motor mixtures [27] can be modelled as active
nematics. There is little work about phase ordering in active nematics [28,
29], but recently microphase separation was predicted to occur in extensile
gels near the nematic-isotropic transition point [30].

It is still not clear whether and how active stresses can lead to phase
ordering in a binary mixture of an active nematic and a passive fluid in the
absence of any thermodynamic driving force. To answer this question, we extend
the continuum theories of active nematohydrodynamics to a two-fluid model with
separate velocity fields, coupled through viscous drag, for each fluid
component [31, 32]. This allows us to model relative motion between the two
fluids which is a requirement for activity-driven phase ordering. Two-fluid
models that have been used in biological contexts include studies of
cytoskeleton dynamics and growing biofilms [33, 34, 35].

Model.—We model a binary mixture of an active nematic fluid (component 1), and
a passive isotropic fluid (component 2), in $d$ dimensions. Each fluid
component, denoted by $i\in{1,2}$, has a density field $\rho^{i}$, velocity
field $\mathbf{u}^{i}$ and a chemical potential $\mathbf{\mu}^{i}$. The active
nematic is associated with a second rank tensor order parameter field [36],
$\mathbf{Q}=\frac{d}{d-1}S\big{(}\mathbf{nn}-\frac{\mathbb{I}}{d}\big{)}$,
where $S$ denotes the degree of nematic ordering, and $\mathbf{n}$ denotes the
orientation of the local nematic director field.

Each component fluid $i$ obeys the mass continuity equation and the momentum
balance equations

$\displaystyle\partial_{t}\rho^{i}+\nabla\cdot(\rho^{i}\mathbf{u}^{i})=$
$\displaystyle 0,$ (1)
$\displaystyle\partial_{t}(\rho^{i}\mathbf{u}^{i})+\nabla\cdot(\rho^{i}\mathbf{u}^{i}\mathbf{u}^{i})=$
$\displaystyle-\rho^{i}\nabla\mu^{i}+\gamma\phi(1-\phi)(\mathbf{u}^{3-i}-\mathbf{u}^{i})$
$\displaystyle+\nabla\cdot\mathbf{\sigma}^{visc,i}+\delta_{i1}\nabla\cdot\mathbf{\sigma}^{nem},$
(2)

where $\phi=\rho^{1}/\rho^{c}$ is the concentration of the active fluid, and
$\rho^{c}=\rho^{1}+\rho^{2}$ is the total fluid density. The forces acting on
the fluids, on the rhs of Eq. (2), are thermodynamic forces, viscous drag
between the component fluids, internal viscous dissipation
$\sigma^{visc,i}=\eta^{i}\big{(}\,(\nabla\mathbf{u}^{i})+(\nabla\mathbf{u}^{i})^{T}-\frac{2}{d}(\nabla\cdot\mathbf{u}^{i})\mathbb{I}\,\big{)}$,
and, for the active nematic fluid, elastic and active forces [20],
respectively. The stresses in the nematic are given by

$\displaystyle\allowdisplaybreaks\mathbf{\sigma}^{el}$
$\displaystyle=2\lambda\big{(}\mathbf{Q}+\mathbb{I}/d\big{)}(\mathbf{Q:H})-\lambda\mathbf{H}\cdot\big{(}\mathbf{Q}+\mathbb{I}/d\big{)}-$
(3)
$\displaystyle\lambda\big{(}\mathbf{Q}+\mathbb{I}/d\big{)}\cdot\mathbf{H}-\big{(}\nabla\mathbf{Q}\big{)}\,.\frac{\partial\mathcal{F}}{\partial\nabla\mathbf{Q}}+\mathbf{Q\cdot
H}-\mathbf{H\cdot Q},$ $\displaystyle\mathbf{\sigma}^{act}$
$\displaystyle=-\zeta\mathbf{Q},$ (4) $\displaystyle\mathbf{\sigma}^{nem}$
$\displaystyle=\phi\,\big{(}\,\mathbf{\sigma}^{el}+\mathbf{\sigma}^{act}\big{)},$
(5)

where $\mathcal{F}$ is the free energy, $\mathbf{H}$ is the molecular field
[20] defined by
$\mathbf{H}=-\frac{\partial\mathcal{F}}{\partial\mathbf{Q}}+\frac{\mathbb{I}}{d}\text{Tr}\bigg{(}\frac{\partial\mathcal{F}}{\partial\mathbf{Q}}\bigg{)}$,
$\lambda$ is the flow aligning parameter, and $\zeta$ is the magnitude of the
activity. Eq. (3) is the (passive) elastic stress [36], while Eq. (4) is the
active stress which produces pusher ($\zeta>0)$ or puller ($\zeta<0$) dipolar
flows.

The nematic tensor field evolves according to [37]

$\displaystyle\partial_{t}\mathbf{Q}+\mathbf{u}^{1}\cdot\nabla\mathbf{Q}-\mathbf{S}=\Gamma\mathbf{H}$
(6)

where
$\mathbf{S}=(\lambda\mathbf{E}^{1}+\mathbf{\Omega}^{1})\cdot\big{(}\mathbf{Q}+\mathbb{I}/d\big{)}+\big{(}\mathbf{Q}+\mathbb{I}/d\big{)}\cdot(\lambda\mathbf{E}^{1}-\mathbf{\Omega}^{1})-2\lambda\big{(}\mathbf{Q}+\mathbb{I}/d\big{)}(\mathbf{Q}:\nabla\mathbf{u}^{1})$
is the co-rotation term describing the response of the orientation field to
strain and vorticity in the flow, while the right hand side describes the
relaxation to a state of minimum free energy. Here, $\mathbf{E}^{1}$ the rate
of strain tensor and $\mathbf{\Omega}^{1}$ the vorticity tensor of the active
fluid.

Equilibrium is described by a free energy [20, 21, 28, 38]

$\displaystyle\mathcal{F}=\int
d^{2}\mathbf{r}\bigg{[}\rho^{c}\Big{(}\frac{1}{3}\ln\rho^{c}+a\,(\phi-1/2)^{2}\Big{)}+\frac{1}{2}\kappa||\mathbf{\nabla}\phi||^{2}$
$\displaystyle+\frac{C}{2}\big{(}S_{nem}^{2}\phi-\frac{1}{2}\text{Tr}(\mathbf{Q}^{2})\big{)}^{2}+\frac{K_{LC}}{2}(\nabla\mathbf{Q})^{2}\bigg{]}$
(7)

where the nematic elastic constant $K_{LC}$, $C$, $a$ and $\kappa$ are
material parameters. The first term represents the pressure contribution to
the free energy from the isothermal equation of state for the fluid [39]. The
rest of the first line describes a Landau-Ginzburg free energy with $a>0$
favouring a homogeneous mixed state ($\phi=1/2$) and $\kappa>0$ penalizing
interfaces. The second line is the Landau-de Gennes free energy describing an
ordered nematic.

We change variables to give one equation for the combined fluid, which is to a
good approximation incompressible, and one for the active fluid. To do this,
we define the centre of mass velocity of the combined fluid
$\mathbf{u}^{c}=\phi\mathbf{u}^{1}+(1-\phi)\mathbf{u}^{2}$, and the relative
slipping velocity between the fluids
$\delta\mathbf{u}=\mathbf{u}^{1}-\mathbf{u}^{2}$. Assuming that the viscous
drag between the components is the fastest relaxation process in the system,
$\delta\mathbf{u}<<\mathbf{u}^{1},\mathbf{u}^{2}$, we neglect all terms of
order $(\delta\mathbf{u})^{2}$, and assume that the drag between the fluids is
linear in $\delta\mathbf{u}$. Moreover, recalling that
$\mu^{i}=\partial\mathcal{F}/\partial\rho^{i}$, the thermodynamic forces can
be written as a stress tensor
$\sum_{i=1}^{2}(-\rho^{i}\mathbf{\nabla}\mu^{i})=-\mathbf{\nabla\cdot\sigma}^{th}$
which, for the free energy in Eq. (7), is
$\mathbf{\sigma}^{th}=p\mathbb{I}+\kappa\big{(}\mathbf{\nabla\nabla}\phi-\frac{1}{2}||\mathbf{\nabla}\phi||^{2}\mathbb{I}\big{)}$
where $p=\rho^{c}/3$ [39, 38, 40]. Assuming that both fluids have the same
viscosity $\eta$, the internal viscous dissipation of the combined fluid is
$\mathbf{F}^{visc,c}=\eta^{c}\nabla^{2}\mathbf{u}^{c}$ where $\eta^{c}=2\eta$.

Adding Eqs. (1), (2) for each component then shows that the combined fluid
satisfies

$\displaystyle\partial_{t}\rho^{c}+\mathbf{\nabla}\cdot(\rho^{c}\mathbf{u}^{c})=0,$
(8)
$\displaystyle\partial_{t}(\rho^{c}\mathbf{u}^{c})+\mathbf{\nabla}\cdot(\rho^{c}\mathbf{u}^{c}\mathbf{u}^{c})=-\mathbf{\nabla}\cdot\sigma^{th}+\mathbf{F}^{visc,c}+\mathbf{\nabla}\cdot\sigma^{nem}.$
(9)

The equations for the first (active) component are

$\displaystyle\partial_{t}\rho^{1}+\mathbf{\nabla}\cdot(\rho^{1}\mathbf{u}^{1})=$
$\displaystyle\,0,$ (10)
$\displaystyle\rho^{1}[\partial_{t}\mathbf{u}^{1}+\mathbf{u}^{1}\cdot\mathbf{\nabla}\mathbf{u}^{1}]=$
$\displaystyle-\phi\mathbf{\nabla}\cdot\mathbf{\sigma}^{th}+\phi\mathbf{G}+\mathbf{F}^{visc,1}$
$\displaystyle+\mathbf{\nabla}\cdot\mathbf{\sigma}^{nem}$ (11)

where
$\mathbf{G}=\gamma(\mathbf{u}^{c}-\mathbf{u}^{1})+(1-\phi)\mathbf{\nabla}\delta\mu$
is the force applied by the passive component, and
$\delta\mu=-\delta\mathcal{F}/\delta\phi$.

Eqs. (8-11) are solved using a lattice Boltzmann (LB) algorithm [38, 40]. For
the equation of the $\mathbf{Q}$ field (6) we use a finite difference (FD)
approach. 5 FD steps are taken for each LB step to optimize simulation speed.
Simulations are run in a periodic box of size 200 x 200 for 20,000 LB time-
steps, with a random initial director configuration. Parameter values are
$\rho^{1}=20,\rho^{2}=20,\gamma=4,\eta_{1}=\eta_{2}=10/3,a=0.0025,\kappa=5,\Gamma=0.1,C=0.1,\zeta=0.1,K_{LC}=0.15,S_{nem}=1,\lambda=0$
unless otherwise indicated.

Phase separation driven by activity.—We solve the equations of motion for a
mixture which comprises equal concentrations of an active nematic fluid and a
passive fluid. The mixture is incompressible so that the total fluid density
remains constant, but the relative concentrations of individual components can
vary in space and time. The two components are initially homogeneously mixed
so that the concentration of the active component $\phi=0.5$ everywhere.
However, above a threshold activity the system quickly self-organises into
dynamically-evolving regions which are, respectively, rich or poor in the
active fluid component. A typical snapshot of the phase-ordered state is shown
in Fig. 2(a), and Movie 1 in the Supplemental Material (SM) [41] shows its
time evolution. The domains continually break up and re-form. They are
elongated by extensile active flows [28] and then pulled apart by the chaotic
active turbulent flows .

We measure the standard deviation of $\phi$, denoted by $\Delta$, to quantify
the degree of phase separation. Fig. 2(b) shows how $\Delta$ evolves with time
towards a dynamical steady state.

Figure 2: (color online.) Active phase ordering: (a) Snapshot of the
concentration profile in the micro-phase-separated state, color bar denotes
concentration of the active phase (b) Standard deviation of the concentration,
$\Delta$, as a function of time.

Mechanism.—We emphasise that the phase separation is entirely due to the
activity, and requires no passive forces. The free energy in Eq. (7) is chosen
to be a single well potential, $a>0$, which favours mixing. To explain the
mechanism causing the phase ordering consider a small fluctuation in
concentration. From Eq. (5) this will lead to an active force tangential to
the interface between the higher and lower activity regions
$\mathbf{F}^{\footnotesize{\mbox{tangential}}}=2\zeta|\nabla(S\phi)|\,(\mathbf{m\cdot
n})(\mathbf{l\cdot n})\,\mathbf{l}$ where $\mathbf{m}$ and $\mathbf{l}$ are
unit vectors normal (pointing away from the more active region) and tangential
to the interface respectively, and $\mathbf{n}$ is a unit vector along the
nematic director (Fig 3(a)) [28], see also SM, Sec. 1 [41]. This force sets up
flows parallel to the interface. The nematic director rotates in the flow
leading to a net alignment parallel (perpendicular) to the interface in
extensile (contractile) systems. This is termed active anchoring [28, 42].

Similarly the active force normal to the interface is
$\mathbf{F}^{\footnotesize{\mbox{normal}}}=-\zeta|\nabla(S\phi)|\,(2(\mathbf{m\cdot
n})^{2}-1)\,\mathbf{m}.$ In the absence of anchoring,
$F^{\footnotesize{\mbox{normal}}}$ averages to zero. However, when there is an
active anchoring of the director, the force drives the active component up the
concentration gradient towards the active nematic region. Hence the relative
flow between the two components enhances concentration fluctuations to form
active regions of concentration $\phi>1/2$. The same mechanism leads to domain
formation in contractile systems, as the change in the direction of the active
anchoring is compensated by the change in sign of the activity. This is
illustrated schematically in Figs. 3(b),(c). A simulation snapshot of an
active droplet is shown in Figs. 3 (d)-(g) .

Figure 3: (color online.) Active anchoring at an interface leads to a net flow
normal to the interface: (a) Vectors introduced in the text. (b) Extensile
nematics with planar active anchoring (inset: extensile flow field); (c)
Contractile nematics with homeotropic active anchoring (inset: contractile
flow field). In both cases the changing concentration of nematic particles
leads to a net flow up the concentration gradient (red arrows). (d)-(g):
Simulation snapshot of an actively phase-separated droplet. (a) Concentration
field. (b) Director field, with planar active anchoring everywhere except at
the tips. (c) Velocity difference between the active and passive fluid
components, which points into the droplet except near $+1/2$ defects. (d) The
combined flow field stretches the droplet, in a way identical to the single
fluid result.

To confirm the relevance of interfacial flows we numerically measured the
dependence of physical quantities on the modulus of the concentration
gradient, averaging over the whole domain. Projecting the nematic director
field onto the concentration gradient yields
$\cos(2\theta_{2})=(\mathbf{\nabla\phi}^{T}\cdot\mathbf{Q}\cdot\mathbf{\nabla\phi})/(S|\nabla\phi|^{2})$,
where $\theta_{2}$ is the angle between the director field and $\nabla\phi$.
This quantity tends to -1 for extensile activity confirming planar active
anchoring (Fig. 4(a)). Similarly, it approaches +1 for contractile activity
corresponding to homeotropic active anchoring (Fig. A1 in SM [41]). The
magnitude of the slip velocity $\mid\\!\\!\delta\mathbf{u}\\!\\!\mid$ is
independent of the concentration gradient $\mid\\!\\!\nabla\phi\\!\\!\mid$
(Fig. 4(b)). However, the cosine of the angle between $\delta\mathbf{u}$ and
$\nabla\phi$,
$\cos(\theta_{1})=(\mathbf{\nabla\phi}\cdot\delta\mathbf{u})/(|\nabla\phi||\delta\mathbf{u}|)$,
increases at large $|\nabla\phi|$, confirming that the direction of the
relative fluid velocity is preferentially up concentration gradients (Fig.
4(c)).

Figure 4: Flow–interface correlations: (a) Active anchoring increases with
concentration gradients. (b) Magnitude of mean slip velocity does not depend
on concentration gradients. (c) Orientation of mean slip velocity points up
concentration gradients. Results are for extensile activity. Shaded regions
show standard deviations of observed values.

Varying the Model Parameters.—We now discuss more detailed numerical results
investigating how the size of the phase separated domains and the difference
in concentration between the two phases, $\Delta$, depend on the model
parameters. Fig. 5(a) shows that $\Delta$ increases, but does not scale
linearly, with the activity coefficient $\zeta$. This is because at higher
activities the flow is more turbulent with many defects which reduce the
strength of the active anchoring.

The formation of concentration gradients is opposed by the bulk free energy
which scales as $a$, the relative drag $\gamma$, and the surface tension
$\kappa$. We consider the balance between the activity and each of these in
turn, choosing parameters where the given relaxation mechanism is dominant.
Results are presented for extensile activity but they also hold for
contractile systems since the same mechanisms are at play (SM, Fig. A2 [41]).

Figure 5: (color online.) Variation of concentration difference $\Delta$ and
average domain size with model parameters: In all the plots above, we vary
activity $\zeta$ $\in[0.03,0.18]$.(a) $\Delta$ scales non-linearly on
increasing activity. Ordering is primarily opposed by (b) bulk free energy
$a\in[0.0010,0.0075]$, (c) drag between the fluids $\gamma\in[4,22]$, (d)
surface tension $\kappa\in[5,25]$. Average domain size as a function of the
(e) active length-scale, $L_{\mbox{act}}$, for $K_{LC}\in[0.05,0.40]$. (f)
Surface tension length-scale, $L_{\mbox{surf}}$, for $\kappa\in[5,25]$. Blue
curves are best-fit power laws for (a)-(d), and best-fit straight lines for
(e), (f). Error bars show standard deviations from 40 simulation values.

In Fig. 5(b) $\Delta$ is plotted as a function of $\zeta/a^{0.5}$ showing
collapse onto a single curve. The collapse is a result of the balance between
the driving force normal to the interface and the thermodynamic force which
restores mixing. The former scales as $\zeta^{2}$ because it depends on the
strength of the anchoring and on the magnitude of the active force, both of
which scale with $\zeta$. The latter scales as $a$.

The data also lies on a single curve when $\Delta$ is plotted against the
ratio $\zeta/\gamma$ (Fig. 5(c)). This shows the expected inverse relation
between the activity, which creates a slipping velocity between the active and
passive fluids, and the viscous drag that dampens the slip.

A similar plot, obtained by changing the surface tension coefficient $\kappa$
at constant $a$ and $\gamma$, shows the best data collapse when $\Delta$ is
plotted against $\zeta/\kappa^{0.8}$ (Fig. 5(d)). We do not have a simple
argument for this dependence as several factors contribute. In addition to the
thermodynamic penalty, an increase in surface tension weakens the activity
gradient by making the interface more diffuse, and hence will alter the
details of the active interfacial forces and active anchoring.

The length-scale of the domains can be calculated numerically from the first
zero of the spatial correlation function of the fluctuations in the
concentration field. For fixed surface tension, Fig. 5e shows that the domain
size is controlled by the usual active nematic length scale
$L_{\footnotesize{\mbox{act}}}=\sqrt{{K_{LC}}/{\zeta}}$. This is, however, not
the only length scale controlling the physics. The domain size can be changed
by varying a second length scale, related to the surface tension,
$L_{\footnotesize{\mbox{surf}}}=\sqrt{{\kappa\eta}/{\zeta}}$, as shown in Fig.
5f. Active droplets break up by forming bend or splay instabilities with
typical lengthscale $L_{act}$, or by self-shearing instabilities [43, 28] with
typical lengthscale $L_{surf}$. Snapshots showing formation and dissolution of
a droplet are shown in SM, Fig. A3 [41].

Threshold for phase ordering.—It is possible to obtain an exact expression for
the activity threshold above which phase ordering occurs, in the limit of
perfect surface anchoring, and ignoring interface curvature. The effects of
curvature are discussed in Sec. 2 of the SM [41].

Consider a concentration gradient around $\phi=1/2$ in the $x$-direction,
given by $\nabla\phi=\mid\\!\\!\nabla\phi\\!\\!\mid\,\\!\\!\hat{x}$. Let the
activity coefficient be $\pm|\zeta|$, where the $+$ $(-)$ sign denotes
extensile (contractile) activity. Assuming strong active anchoring at the
interface, this leads to a director field $Q_{xx}=\mp
S_{\footnotesize{\mbox{nem}}}\sqrt{\phi},Q_{xy}=0$. The normal active force at
the interface is then
$\mathbf{F}_{act}=(3/2)S_{nem}\sqrt{\phi}\mid\\!\\!\zeta\\!\\!\mid\,\mid\\!\\!\nabla\phi\\!\\!\mid\,\,\mathbf{\hat{x}}$,
and the restoring force driven by the free energy is
$\mathbf{F}_{rest}=-2a\rho\phi(1-\phi)|\nabla\phi|\,\,\mathbf{\hat{x}}$. This
leads to an expression for the slip velocity $\delta u_{x}$ in the Stokes
limit:

$\displaystyle\delta
u_{x}=\frac{F_{act}+F_{rest}}{\gamma\phi(1-\phi)}=\frac{(3/2)\mid\\!\\!\zeta\\!\\!\mid
S_{nem}\sqrt{\phi}-2a\rho\phi(1-\phi)}{\gamma\phi(1-\phi)}\mid\\!\\!\nabla\phi\\!\\!\mid\hat{x}.$
(12)

The concentration gradient will increase for activities greater than
$\mid\\!\\!\zeta\\!\\!\mid_{crit}=\frac{4a\rho\phi(1-\phi)}{3\sqrt{\phi}S_{nem}}.$
(see SM, Fig. A4 [41].)

Summary & outlook.— By introducing an active two-fluid model with viscous drag
between the fluids we have found that a homogeneous mixture of an active
nematic fluid and a passive fluid spontaneously phase orders into regions with
different concentrations of each fluid component. This is driven dynamically
by active flows set up at concentration interfaces, and is an example of
liquid-liquid phase separation [44] out of thermodynamic equilibrium.

The two-fluid model reduces to the lyotropic formulation of
nematohydrodynamics that has been successful in defining active and passive
regions within a simulation, in the limit of infinite viscous drag between the
fluids. However, because the lyotropic approach uses a single velocity field
for both phases [28, 45], it cannot account for relative flows between, or
compressibility of, the two components.

Active nematics have been observed to preferentially adhere to a substrate or
container wall, a phenomenon known as active wetting [46, 47, 48, 49]. A
straightforward extension of our argument shows that for an initially
homogeneous active-isotropic mixture, the normal force generated at a boundary
would cause wetting of the boundary by the active nematic for both extensile
and contractile activity, even in the absence of surface thermodynamic forces.
The wetting would then be enhanced by any thermodynamic planar (normal)
anchoring for extensile (contractile) activity.

In future work it will be interesting to study mixtures where both the fluids
are active. In addition to cell sorting this is relevant to the growth and
dynamics of bacterial and algal colonies which contain more than one species
[50]. Moreover, the relaxation of the incompressibility constraint will allow
the study of how compressible flows affect defect behaviour, and whether they
can stabilize structures such as asters or the cellular lumen.

###### Acknowledgements.

Acknowledgements.—We thank L. J. Ruske and I. Hadjifrangiskou for valuable
discussions. SB acknowledges support from the Rhodes Trust and the Crewe
Graduate Scholarship.

## References

  * Halbleib and Nelson [2006] J. M. Halbleib and W. J. Nelson, Cadherins in development: cell adhesion, sorting, and tissue morphogenesis, Genes Dev 20, 3199 (2006).
  * Krens and Heisenberg [2011] S. F. Krens and C.-P. Heisenberg, Chapter six - cell sorting in development,  _Forces and Tension in Development_ , Current Topics in Developmental Biology 95, 189 (2011).
  * Durand [2021] M. Durand, Large-scale simulations of biological cell sorting driven by differential adhesion follow diffusion-limited domain coalescence regime, PLOS Computational Biology 17, 1 (2021).
  * Skamrahl  _et al._ [2022] M. Skamrahl, J. Schünemann, M. Mukenhirn, J. Gottwald, M. Ferle, A. Rübeling, A. Honigmann, and A. Janshoff, Cellular segregation in co-cultures driven by differential adhesion and contractility on distinct time scales, bioRxiv:2022.05.23.492966 (2022).
  * Steinberg [1975] M. S. Steinberg, Adhesion-guided multicellular assembly: a commentary upon the postulates, real and imagined, of the differential adhesion hypothesis, with special attention to computer simulations of cell sorting, Journal of Theoretical Biology 55, 431 (1975).
  * Banani  _et al._ [2017] S. F. Banani, H. O. Lee, A. A. Hyman, and M. K. Rosen, Biomolecular condensates: organizers of cellular biochemistry, Nature Reviews Molecular Cell Biology 18, 285 (2017).
  * Dutagaci  _et al._ [2021] B. Dutagaci, G. Nawrocki, J. Goodluck, A. A. Ashkarran, C. G. Hoogstraten, L. J. Lapidus, and M. Feig, Charge-driven condensation of RNA and proteins suggests broad role of phase separation cytoplasmic environments, eLife 10, e64004 (2021).
  * Suzuki  _et al._ [2021] S. Suzuki, I. Omori, R. Kuraishi, and H. Kaneko, Cell sorting and germ layer formation in reconstructed starfish embryos, Development, Growth & Differentiation 63, 343 (2021).
  * Balasubramaniam  _et al._ [2021] L. Balasubramaniam, A. Doostmohammadi, T. B. Saw, G. H. N. S. Narayana, R. Mueller, T. Dang, M. Thomas, S. Gupta, S. Sonam, A. S. Yap, Y. Toyama, R.-M. Mège, J. M. Yeomans, and B. Ladoux, Investigating the nature of active forces in tissues reveals how contractile cells can form extensile monolayers, Nature Materials 20, 1156 (2021).
  * Ramaswamy [2010] S. Ramaswamy, The mechanics and statistics of active matter, Annual Review of Condensed Matter Physics 1, 323 (2010).
  * Cates and Tailleur [2015] M. E. Cates and J. Tailleur, Motility-induced phase separation, Annual Review of Condensed Matter Physics 6, 219 (2015).
  * Gonnella  _et al._ [2015] G. Gonnella, D. Marenduzzo, A. Suma, and A. Tiribocchi, Motility-induced phase separation and coarsening in active matter, Comptes Rendus Physique 16, 316 (2015).
  * Stenhammar  _et al._ [2013] J. Stenhammar, A. Tiribocchi, R. J. Allen, D. Marenduzzo, and M. E. Cates, Continuum theory of phase separation kinetics for active brownian particles, Phys. Rev. Lett. 111, 145702 (2013).
  * Wittkowski  _et al._ [2014] R. Wittkowski, A. Tiribocchi, J. Stenhammar, R. J. Allen, D. Marenduzzo, and M. E. Cates, Scalar $\phi$4 field theory for active-particle phase separation, Nature Communications 5, 4351 (2014).
  * Redner  _et al._ [2013] G. S. Redner, M. F. Hagan, and A. Baskaran, Structure and dynamics of a phase-separating active colloidal fluid, Phys. Rev. Lett. 110, 055701 (2013).
  * Kuan  _et al._ [2021] H.-S. Kuan, W. Pönisch, F. Jülicher, and V. Zaburdaev, Continuum theory of active phase separation in cellular aggregates, Phys. Rev. Lett. 126, 018102 (2021).
  * Thutupalli  _et al._ [2018] S. Thutupalli, D. Geyer, R. Singh, R. Adhikari, and H. A. Stone, Flow-induced phase separation of active particles is controlled by boundary conditions, Proceedings of the National Academy of Sciences of the United States of America 115, 5403 (2018).
  * Zhang  _et al._ [2021] J. Zhang, R. Alert, J. Yan, N. S. Wingreen, and S. Granick, Active phase separation by turning towards regions of higher density, Nature Physics 17, 961 (2021).
  * Gokhale  _et al._ [2022] S. Gokhale, J. Li, A. Solon, J. Gore, and N. Fakhri, Dynamic clustering of passive colloids in dense suspensions of motile bacteria, Phys. Rev. E 105, 054605 (2022).
  * Doostmohammadi  _et al._ [2018] A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans, and F. Sagués, Active nematics, Nature Communications 9, 3246 (2018).
  * Marchetti  _et al._ [2013] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Hydrodynamics of soft active matter, Rev. Mod. Phys. 85, 1143 (2013).
  * Duclos  _et al._ [2017] G. Duclos, C. Erlenkamper, J.-F. Joanny, and P. Silberzan, Topological defects in confined populations of spindle-shaped cells, Nature Physics 13, 58 (2017).
  * Yaman  _et al._ [2019] Y. I. Yaman, E. Demir, R. Vetter, and A. Kocabas, Emergence of active nematics in chaining bacterial biofilms, Nature Communications 10, 2285 (2019).
  * Dell’Arciprete  _et al._ [2018] D. Dell’Arciprete, M. L. Blow, A. T. Brown, F. D. C. Farrell, J. S. Lintuvuori, A. F. McVey, D. Marenduzzo, and W. C. K. Poon, A growing bacterial colony in two dimensions as an active nematic, Nature Communications 9, 4190 (2018).
  * Saw  _et al._ [2017] T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans, and B. Ladoux, Topological defects in epithelia govern cell death and extrusion, Nature 544, 212 (2017).
  * Balasubramaniam  _et al._ [2022] L. Balasubramaniam, R.-M. Mège, and B. Ladoux, Active nematics across scales from cytoskeleton organization to tissue morphogenesis, Current Opinion in Genetics & Development 73, 101897 (2022).
  * Sanchez  _et al._ [2012] T. Sanchez, D. T. Chen, S. J. DeCamp, M. Heymann, and Z. Dogic, Spontaneous motion in hierarchically assembled active matter, Nature 491, 431 (2012).
  * Blow  _et al._ [2014] M. L. Blow, S. P. Thampi, and J. M. Yeomans, Biphasic, Lyotropic, Active Nematics, Phys. Rev. Lett. 113, 248303 (2014).
  * Caballero and Marchetti [2022] F. Caballero and M. C. Marchetti, Activity-suppressed phase separation, Phys. Rev. Lett. 129, 268002 (2022).
  * Assante  _et al._ [2022] R. Assante, D. Corbett, D. Marenduzzo, and A. Morozov, Active turbulence and spontaneous phase separation in inhomogeneous extensile active gels, arXiv:2209.01869 (2022).
  * Stewart and Wendroff [1984] H. B. Stewart and B. Wendroff, Two-phase flow: models and methods, Journal of Computational Physics 56, 363 (1984).
  * Levy and Lévy [1999] S. Levy and S. Lévy,  _Two-Phase Flow in Complex Systems_, A Wiley-Interscience publication (Wiley, 1999).
  * Mogilner and Manhart [2018] A. Mogilner and A. Manhart, Intracellular fluid mechanics: Coupling cytoplasmic flow with active cytoskeletal gel, Annual Review of Fluid Mechanics 50, 347 (2018).
  * Cogan and Guy [2010] N. G. Cogan and R. D. Guy, Multiphase flow models of biogels from crawling cells to bacterial biofilms, HFSP Journal 4, 11 (2010).
  * Alt and Dembo [1999] W. Alt and M. Dembo, Cytoplasm dynamics and cell motion: two-phase flow models, Mathematical Biosciences 156, 207 (1999).
  * de Gennes and Prost [1993] P. de Gennes and J. Prost,  _The Physics of Liquid Crystals_, International Series of Monogr (Clarendon Press, 1993).
  * Beris and Edwards [1994] A. N. Beris and B. J. Edwards,  _Thermodynamics of Flowing Systems with Internal Microstructure_, Oxford Engineering Science Series (Oxford University Press, New York, 1994).
  * Malevanets and M. Yeomans [1999] A. Malevanets and J. M. Yeomans, Lattice boltzmann simulations of complex fluids: viscoelastic effects under oscillatory shear, Faraday Discuss. 112, 237 (1999).
  * Krueger  _et al._ [2016] T. Krueger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. Viggen,  _The Lattice Boltzmann Method: Principles and Practice_ , Graduate Texts in Physics (Springer, 2016).
  * Malevanets and M. Yeomans [2000] A. Malevanets and J. M. Yeomans, Modeling viscous drag in binary fluid mixtures, Computer Physics Communications 127, 105 (2000).
  * [41] See Supplemental Material at [URL will be inserted by publisher].
  * Ruske and Yeomans [2022] L. J. Ruske and J. M. Yeomans, Activity gradients in two- and three-dimensional active nematics, Soft Matter 18, 5654 (2022).
  * Singh and Cates [2019] R. Singh and M. E. Cates, Hydrodynamically interrupted droplet growth in scalar active matter, Phys. Rev. Lett. 123, 148005 (2019).
  * Hyman  _et al._ [2014] A. A. Hyman, C. A. Weber, and F. Jülicher, Liquid-liquid phase separation in biology, Annual Review of Cell and Developmental Biology 30, 39 (2014).
  * Giomi and DeSimone [2014] L. Giomi and A. DeSimone, Spontaneous division and motility in active nematic droplets, Physical Review Letters 112, 147802 (2014).
  * Adkins  _et al._ [2022] R. Adkins, I. Kolvin, Z. You, S. Witthaus, M. C. Marchetti, and Z. Dogic, Dynamics of active liquid interfaces, Science 377, 768 (2022).
  * Pérez-González  _et al._ [2019] C. Pérez-González, R. Alert, C. Blanch-Mercader, M. Gómez-González, T. Kolodziej, E. Bazellieres, J. Casademunt, and X. Trepat, Active wetting of epithelial tissues, Nature Physics 15, 79 (2019).
  * Alert and Trepat [2020] R. Alert and X. Trepat, Physical models of collective cell migration, Annual Review of Condensed Matter Physics 11, 77 (2020).
  * Banerjee and Marchetti [2019] S. Banerjee and M. C. Marchetti, Continuum models of collective cell migration,  _Cell Migrations: Causes And Functions_ , Advances in Experimental Medicine and Biology 1146, 45 (2019).
  * Ramanan  _et al._ [2016] R. Ramanan, B.-H. Kim, D.-H. Cho, H.-M. Oh, and H.-S. Kim, Algae–bacteria interactions: Evolution, ecology and emerging applications, Biotechnology Advances 34, 14 (2016).


# Disorder-induced local strain distribution in Y-substituted TmVO4

Yuntian Li  Department of Applied Physics, Stanford University, Stanford,
California 94305, USA  Geballe Laboratory for Advanced Materials, Stanford
University, Stanford, California 94305, USA  Mark P. Zic  Geballe Laboratory
for Advanced Materials, Stanford University, Stanford, California 94305, USA
Department of Physics, Stanford University, Stanford, California 94305, USA
Linda Ye  Department of Applied Physics, Stanford University, Stanford,
California 94305, USA  Geballe Laboratory for Advanced Materials, Stanford
University, Stanford, California 94305, USA  W. Joe Meese  School of Physics
and Astronomy, University of Minnesota, Minneapolis, Minnesota 55455, USA
Pierre Massat  Geballe Laboratory for Advanced Materials, Stanford University,
Stanford, California 94305, USA  Department of Physics, Stanford University,
Stanford, California 94305, USA  Yanbing Zhu  Department of Applied Physics,
Stanford University, Stanford, California 94305, USA  Rafael M. Fernandes
School of Physics and Astronomy, University of Minnesota, Minneapolis,
Minnesota 55455, USA  Ian R. Fisher  Department of Applied Physics, Stanford
University, Stanford, California 94305, USA  Geballe Laboratory for Advanced
Materials, Stanford University, Stanford, California 94305, USA

###### Abstract

We report an investigation of the effect of substitution of Y for Tm in
Tm_1-xY_xVO4 via low-temperature heat capacity measurements, with the yttrium
content $x$ varying from $0$ to $0.997$. Because the Tm ions support a local
quadrupolar (nematic) moment, they act as reporters of the local strain state
in the material, with the splitting of the ion’s non-Kramers crystal field
groundstate proportional to the quadrature sum of the in-plane tetragonal
symmetry-breaking transverse and longitudinal strains experienced by each ion
individually. Analysis of the heat capacity therefore provides detailed
insights into the distribution of local strains that arise as a consequence of
the chemical substitution. These local strains suppress long-range quadrupole
order for $x>0.22$, and result in a broad Schottky-like feature for higher
concentrations. Heat capacity data are compared to expectations for a
distribution of uncorrelated (random) strains. For dilute Tm concentrations,
the heat capacity cannot be accounted for by randomly distributed strains,
demonstrating the presence of significant strain correlations between sites.
For intermediate Tm concentrations, these correlations must still exist, but
the data cannot be distinguished from that which would be obtained from a 2D
Gaussian distribution. The cross-over between these limits is discussed in
terms of the interplay of key length scales in the substituted material. The
central result of this work, that local strains arising from chemical
substitution are not uncorrelated, has implications for the range of validity
of theoretical models based on random effective fields that are used to
describe such chemically substituted materials, particularly when electronic
nematic correlations are present.

Heat Capacity, disordered effect

##  I Introduction

Phase diagrams of a wide variety of strongly correlated electron systems can
be traversed via chemical substitution. The tuning effect of different dopants
in different materials can be attributed to a variety of physical effects.
However, a pervasive and inescapable effect of such substitution must always
be the associated inhomogeneous internal strains which are induced by
incorporating inequivalent atoms to a crystal lattice. These local strains can
have a profound effect on various types of emergent electronic order and phase
transitions. Here we focus specifically on the case of electronic nematic
order, whose bilinear coupling with strains of the same symmetry leads to a
particularly strong interplay. The bilinear coupling means that strain plays
the role of an effective conjugate field for the nematic order parameter, and
hence the inhomogeneous strains which arise from chemical substitution have
been conventionally understood as a random field Vojta (2019) both
experimentally Ran  _et al._ (2011); Dioguardi  _et al._ (2016); Forrest  _et
al._ (2016) and theoretically Carlson  _et al._ (2006, 2014); Nie  _et al._
(2014); Ye  _et al._ (2022); Meese  _et al._ (2022); Wiecki  _et al._ (2021).
Despite this appealing perspective, strains in solids relax over long
distances, implying a significant correlation between adjacent sites Landau
and Lifshitz (1970); Mura (1987). Similarly, the elastic response of solids
(described by the elastic stiffness tensor) means that local strains of
different symmetry even on a specific site are not necessarily uncorrelated.
These simple perspectives imply that the degree to which a description of the
effects of inhomogeneous strains arising from chemical substitution as a truly
random effect should be examined in some detail. Indeed, a characterization of
“randomness” associated with the long-range effects of local strains in a real
material system, and the applicability of the random field assumption, remain
respectively an open challenge and an open question.

Figure 1:  Crystal structure of the hypothetical material Y1Tm3V4O16 obtained
via density functional theory (DFT) calculations, with Tm (Y) atoms
represented by purple (yellow) spheres. VO4 units are shown as gray
coordination polyhedra, with O atoms as green spheres. Inset shows a magnified
view of a distorted TmO8 polyhedron overlaid with that of pristine TmVO4 (O
atoms in the pristine structure are shown as grey spheres) to illustrate the
distortions induced by the Y substitution The atomic displacements of the O
sites in the inset are magnified 40 times for clarity. Ionic radii of the
Tm^3+ and Y^3+ ions are indicated.

In this work, we attempt to experimentally address the above issues in a model
system exhibiting ferroquadrupolar order, which may be viewed as a simple
local-moment realization of Ising nematic order Maharaj  _et al._ (2017);
Massat  _et al._ (2022)). The subject of this study, TmVO4, is an insulator
with a tetragonal zircon-type RVO4 crystal structure. The crystal electric
field (CEF) ground state of the Tm^3+ ions is a non-Kramers doublet, with the
first excited CEF state well above the doublet in energy ($77$ K,
$54\mathrm{cm}^{-1}$) Morin and Schmitt (1990), rendering TmVO4 a good
approximation to an ideal pseudospin-$\frac{1}{2}$ system at low-temperature
Wortman  _et al._ (1974).

TmVO4 is known to host considerable quadrupole-lattice interactions, resulting
in complete softening of the $c_{66}$ elastic modulus and a continuous
cooperative Jahn-Teller phase transition at $T_{Q}=2.2$ K Melcher (1976);
Massat  _et al._ (2022). The Tm ions in the low-temperature phase develop a
spontaneous local quadrupole moment (splitting the non-Kramers Jahn-Teller-
active CEF groundstate) while the lattice develops a uniform strain
$\epsilon_{xy}$ of the same symmetry. We introduce disorder to TmVO4 by
substituting Tm with Y. Aside from the overall dilution effect of the magnetic
species and quadrupole moments, the substitution of Y introduces local strains
that can couple to the ferroquadrupolar (Ising-nematic) order parameter. To
illustrate the anticipated effect of substitution, Fig.1 shows a schematic of
the local distorted crystal field environments of the Tm^3+ ion obtained from
Density Functional Theory (DFT) calculations of a hypothetical Y1Tm3V4O16
ordered super-cell (investigated for computational simplicity). Inspection of
the figure reveals how Y sites distort nearby TmO8 clusters, resulting in a
change of the local crystal field environment for the Tm3+ ions relative to
the case of TmVO4. While substitution in Tm_1-xY_x VO4 does not result in an
ordered superlattice, the effect of chemical substitution in terms of inducing
local strains on nearby Tm ions will be comparable. It is the effect of these
substitution-induced strains that we investigate in the present study.

Noting that the non-Kramers doublet can be split by both
$\epsilon_{x^{2}-y^{2}}$ and $\epsilon_{xy}$ strains (which transform as the
$B_{1g}$ and $B_{2g}$ irreducible representations of the tetragonal point
group $D_{4h}$, respectively) Maharaj  _et al._ (2017), these two independent
symmetry channels naturally span a two-dimensional space. According to the
central limit theorem, uncorrelated random strains in the pseudo-spin space
follow a Gaussian distribution, with their means centered at zero (Fig. 2(a)).
The total splitting of the non-Kramers doublet $\Delta$ is determined by the
quadrature sum of the strain field in the two-symmetry channels
$\sqrt{(\eta_{1}\epsilon_{x^{2}-y^{2}})^{2}+(\eta_{2}\epsilon_{xy})^{2}}$,
where $\eta_{1}$ and $\eta_{2}$ are nemato-elastic coupling strengths (Fig.
2(b)). As a result, the probability distribution of $\Delta$ can be obtained
from joint probability density functions (PDF) in both dimensions (Fig. 2(c)).
We will return to the exact functional form of the distribution in the case of
uncorrelated local strains (2D-Gaussian) below.

Experimentally, $\Delta$ is accessible through thermodynamic measurements. In
the following, we map out the entropy landscape of $4f$ quadrupoles in
Tm_1-xY_xVO4 ($x=0\sim 1$) with heat capacity experiments, from which we
extract the evolution of local strain distributions with Tm concentration, to
shed light on the correlation between local strains.

A well-known minimal theoretical model that captures the effect of quenched
random fields on phase transitions is the random field Ising model (RFIM) Imry
and Ma (1975); Aharony  _et al._ (1976); Grinstein (1976), in which even weak
random fields will suppress second-order phase transitions to a lower
temperature for systems with either short-ranged or long-ranged interactions
Vojta (2019); Nattermann (1998); Binder (1983); Gehring  _et al._ (1976).
Indeed, in three-dimensional systems, there exists a threshold random field
strength at which the phase transition will be fully suppressed, whereas in
two dimensions long-range order is completely eliminated even at zero
temperature Zachar and Zaliznyak (2003). Due to its conceptual simplicity, the
RFIM continues to be an important model in understanding the phase behavior of
electronic nematics of real materials in the presence of structural disorder
Carlson  _et al._ (2006, 2014); Nie  _et al._ (2014); Ye  _et al._ (2022);
Meese  _et al._ (2022); Wiecki  _et al._ (2021). .

An important definition one must include when applying the RFIM, however, is
“randomness”. In other words, the characterization of a distribution of random
fields that are coupled to the order parameter. This is somewhat arbitrary
since the only distribution moment that can be clearly seen by experimental
probes is the first, i.e. the mean value of the random field. Higher-order
moments, encoding the spatial correlation in the random fields, can be
postulated ab initio Vojta (2019); Nattermann (1998). A frequent choice,
however, is to ignore spatial correlation entirely, and assume furthermore
that the random fields are all identically drawn from a Gaussian distribution
with zero mean and some standard deviation which characterizes the “disorder
strength”. This choice is common for theoretical treatments of electronic
nematic systems as well, but due to the long-range nature of elasticity, it is
unclear whether one can treat the random fields as uncorrelated without losing
essential physics in actual materials. Our study explores evidence for the
presence of these correlations as a first step towards assessing their
significance.

Figure 2: (a) Schematic of the two-dimensional strain distribution $\rho$ as a
function of strain with $B_{1g}$
($\epsilon_{B_{1g}}=\epsilon_{x^{2}-y^{2}}=\frac{1}{2}(\epsilon_{xx}-\epsilon_{yy})$)
and $B_{2g}$ ($\epsilon_{B_{2g}}=\epsilon_{xy}$) symmetry, expected from the
central limit theorem for uncorrelated strains wherein each of them has a
functional form of a Gaussian distribution centered at zero. (b) The magnitude
of the local strain-induced energy gap, $\Delta$, measured in arbitrary units,
as the quadrature sum of the induced gap in both $\epsilon_{B_{1g}}$ and
$\epsilon_{B_{2g}}$ channels. The color scale follows Eq. 5. (c) The resulting
probability distribution $\rho(\Delta)$ from the double-Gaussian distribution
(see section IV).

##  II Experimental Methods

Single crystals of Tm_1-xY_xVO4 used in the heat capacity measurements were
grown via slow cooling in a flux of Pb2V2O7 using a mixture of Tm2O3 and Y2O3
precursors. More details related to material synthesis can be found in Refs.
Feigelson (1968); Smith and Wanklyn (1974); Oka  _et al._ (2006). The Y
content $x$ of each sample is further determined by two different methods: for
$x\leq 0.315$, the composition is determined by microprobe analysis carried
out at the Stanford Mineral and Microchemical Analysis Facility; for
$x>0.315$, the Tm composition is determined by the magnetic moment extracted
from a.c. magnetic susceptibility measurements in a Quantum Design Magnetic
Property Measurement System (MPMS) in a zero DC magnetic field. More details
of the determination of the doping with the MPMS measurement can be found in
Appendix A.

The high-temperature heat capacity of Tm_1-xY_xVO4 (4 K-35 K) was measured in
a Quantum Design Physical Property Measurement System (PPMS) under high vacuum
conditions with the thermal relaxation technique. The low-temperature heat
capacity measurement was carried out with the same method in a dilution
refrigerator insert from 4 K down to the base temperature. All data were taken
at a zero external magnetic field. The data analysis is performed with the
nonlinear least-squares fitting method.

The super-cell in Fig.1 was created using VESTA from the crystal structure of
TmVO4 obtained from the Crystallographic Information File (CIF) in the
Cambridge Crystallographic Data Center (CCDC) (accession code 2117139)Massat
_et al._ (2022). DFT calculations for the ionic relaxation were performed
using the generalized gradient approximation (GGA) exchange-correlation
functional of Perdew, Burke, and Ernzerhof (PBE) Perdew  _et al._ (1996) as
implemented in the Vienna Ab initio Simulation Package (VASP) Kresse and
Furthmüller (1996). The projected augmented wave (PAW) method Blöchl (1994);
Kresse and Joubert (1999) was used for ion–electron interaction, and the
conjugate-gradient algorithm (IBRION = $2$) was used to relax the ions into
their instantaneous ground state.

##  III Heat Capacity of Tm_1-xY_xVO4

In Fig. 3 we show the heat capacity of Tm_1-xY_xVO4 from $x=0$ to $x=0.997$.
Panel (a) shows the measured total heat capacity, $C_{p}$, normalized by the
molar gas constant $R$, while panels (b) and (c) show the derived $4f$
contribution to the heat capacity, $C_{p}^{4f}$, (the isolation of
$C_{p}^{4f}$ from the phonon contributions to $C_{p}$ is described in detail
in Appendix B). The ferroquadrupolar transition temperature $T_{Q}$, defined
here by taking the average of the maxima and minima of the second derivative
of $C_{p}$, is highlighted with a red triangle for each composition in Fig.
3(a). The phase transition is fully suppressed near $x=0.22$. We note that
$T_{Q}$ is suppressed considerably more rapidly than what would be expected
based solely on the dilution effects of the magnetic ions. This result, even
before an entropy analysis, already suggests that local strains induced by
Y-substitution play a key role in the suppression of $T_{Q}$.

Inspection of Fig. 3(a) reveals the presence of a broad shoulder-like feature
in the heat capacity for temperatures above $T_{Q}$ for compositions above
$x=0.15$ (cyan curve in Fig. 3(a)). This shoulder becomes progressively more
prominent with increasing $x$ as the phase transition is fully suppressed
(Fig. 3(b)). This pervasive feature can be attributed to the Schottky anomaly
associated with the strain-induced splitting of the ground state crystal field
doublet of Tm^3+. However, the anomaly is never described by a single gap –
rather, there is a distribution of gaps. To illustrate this point, in Fig.
3(b), we plot $C_{p}^{4f}$ normalized with its maximum $C_{p}^{4f,max}$
against $T$ normalized by $T_{0}$ ($T_{0}$ refers to the temperature at which
$C_{p}^{4f}$ is maximized). For comparison, we include in this figure a trace
(black dashed line) showing the Schottky anomaly function that would arise
from a single gap $\Delta$. The strong deviation of the observed
$C_{p}^{4f}/C_{p}^{4f,max}$ clearly indicates that there exists a finite
distribution of $\Delta$, which originates from the anticipated inhomogeneous
distribution of substitution-induced local strains. Furthermore, it is clear
from inspection of Fig. 3(b) that the distribution of strain-induced gaps must
change as a function of composition (since the overall shape of the normalized
curves changes as a function of composition). Analysis of this variation
informs the key insights of the present work.

The $T$-evolution of $C_{p}^{4f}/T$ is shown in Fig. 3(c) for all
compositions, providing the raw data from which the integrated entropy can be
extracted. In Fig. 4(a), we extract the total entropy associated with the $4f$
degree of freedom via numerical integration $S^{4f}=\sum(C_{p}/T)\Delta T$.
The $\Delta S^{4f}$ plotted here represents the release of entropy with
respect to that at high temperature: $\Delta
S^{4f}(T)=S^{4f}(T)-S^{4f}\mid_{T=T_{max}}$. We then plot the release of the
entropy with respect to the saturated value at high $T$, from the lowest
accessible $T$. For most compositions (all but the very dilute Tm
concentrations), a large fraction of the total entropy $k_{B}\ln(2)$ per Tm is
recovered upon cooling from 4 K down to below 100 mK, indicating that nearly
all Tm^3+ experience a significant local strain-induced splitting. The entropy
landscape in the $x-T$ plane shown in Fig. 4(b) also highlights that beyond
$x=0.22$ where $T_{Q}$ is fully suppressed, the temperature scale associated
with the entropy release first increases and becomes comparable with $T_{Q}$
of the pristine TmVO4 near $x=0.5$, before progressively decreasing towards 0
in the dilute Tm limit. Such energy scales are consistent with the large
magnetoelastic coupling strength experienced by the non-Kramers doublet ground
state of Tm^3+, as well as the rapid suppression of $T_{Q}$ with $x$.

Figure 3: (a) Temperature and chemical composition dependence of molar heat
capacity $C_{p}$ (normalized by $R$) of Tm_1-xY_xVO4 for $x\leq 0.22$. The
peak in $C_{p}$ that corresponds to the ferroquadrupolar phase transition is
labeled by red arrows. (b) $C_{P}^{4f}$ normalized by the maximum
$C_{p}^{max}$ with respect to $T$ normalized by $T_{0}$ (see text) for
$x>0.22$. (c) $T$-evolution of $C_{p}^{4f}/T$ for $x>0.22$. Figure 4: (a)
Entropy change $\Delta S^{4f}$ normalized by $R$ as a function of $T$ at
selected $x$ (see text). The dashed line shows the value of $-\ln(2)$. (b)
Color map of $\Delta S^{4f}$ in the $T-x$ plane. The solid black line
indicates the extracted ferroquadrupolar transition temperature $T_{Q}$.

##  IV Strain distribution in Tm_1-xY_xVO4

In this section, we focus on a quantitative analysis of the evolution of the
gap distribution $\rho(\Delta)$ with $x$. The distribution of gaps in a
locally disordered two-level system was originally proposed in explaining the
specific heat of spin glasses Phillips (1987). A continuous distribution of
splittings of the ground-state doublet has also been reported in different
pyrochlore compounds van Duijn  _et al._ (2005); Tachibana  _et al._ (2007);
Martin  _et al._ (2017). Our analysis builds on these and similar
perspectives, and aims to uncover significant trends affecting the variation
of $\rho(\Delta)$ with composition in Tm_1-xY_xVO4.

Given a known PDF that represents the distributions of the gaps
$\rho(\Delta)$, the total heat capacity can be obtained from the sum of
individual Schottky anomalies associated with each $\Delta$:

$\bar{C}_{p}^{4f}(T)=N\int_{0}^{\infty}{d\Delta\rho(\Delta)\frac{2\Delta^{2}}{T^{2}\cosh^{2}(\frac{2\Delta}{T})}}$
(1)

The inverse problem of determining $\rho(\Delta)$ from a measured heat
capacity is much harder, as it is an inhomogeneous Fredholm integral equation
of the first kind. The solution $\rho(\Delta)$ is then also subject to the
constraints that it is nonnegative and normalized. Such a problem is ill-posed
generally, with approximate numerical solutions that are sensitive to noise in
the data Markovsky (1988); Yuan and Zhang (2019). In light of these
difficulties, we take an alternative approach by circumventing the inverse
problem altogether. Instead, we present a comparison of the measured heat
capacity data to that which would be obtained from specific distributions,
such as the 2D Gaussian distribution discussed below, which can provide a
meaningful perspective on the actual PDF. Note that $\Delta$ here refers to
half of the energy splitting between the two levels; we use this definition in
the fitting and plotting procedure.

We consider the ground state doublet of $\ce{Tm^{3+}}$ in a Tm_1-xY_xVO4
crystal environment. The disordered local strain environment results in a
finite energy splitting via $B_{1g}$ and $B_{2g}$ symmetry distortions, where
$\epsilon_{B_{1g}}=\epsilon_{x^{2}-y^{2}}=\epsilon_{xx}-\epsilon_{yy}$ and
$\epsilon_{B_{2g}}=2\epsilon_{xy}$, with values that vary throughout the
material. The magnitude of the total splitting $\Delta$ on any given site is
proportional to the quadrature sum of the two strain components:

$\Delta=\sqrt{(\eta_{1}\epsilon_{B_{1g}})^{2}+(\eta_{2}\epsilon_{B_{2g}})^{2}}$
(2)

where $\eta_{1}$ and $\eta_{2}$ are the magneto-elastic coupling strength in
the $B_{1g}$ and $B_{2g}$ channels. For simplicity of annotating the effect of
the symmetry-breaking strain as two mutually transverse field components, we
write:

$\displaystyle x$ $\displaystyle\equiv\eta_{1}\epsilon_{B_{1g}}$ (3)
$\displaystyle y$ $\displaystyle\equiv\eta_{2}\epsilon_{B_{2g}}$
$\displaystyle\Delta$
$\displaystyle\equiv\frac{1}{2}\Delta_{\mathrm{Schottky}}=\sqrt{x^{2}+y^{2}}$

We start with the empirical assumption that $x$ and $y$ are independent random
variables that follow two Gaussian distributions $\rho_{x}(x)$ and
$\rho_{y}(y)$, as one would expect from the central limit theorem. The width
of the $\epsilon_{B_{1g}}$ and $\epsilon_{B_{2g}}$ distributions are
$\sigma_{a}$ and $\sigma_{b}$. The expression for $\rho(\Delta)$ is given by
the integral:

$\rho(\Delta)=\int_{-\infty}^{\infty}{\int_{-\infty}^{\infty}{dxdy\rho_{x}(x)\rho_{y}(y)\delta(\Delta-\sqrt{x^{2}-y^{2}})}}$
(4)

The solution for $\rho(\Delta)$ is then:

$\displaystyle\rho(\Delta)$
$\displaystyle=\int_{-\Delta}^{\Delta}dx\frac{1}{\sigma_{x}\sqrt{2\pi}}e^{-\frac{1}{2}(\frac{x}{\sigma_{x}})^{2}}\frac{\Delta}{\sqrt{\Delta^{2}-x^{2}}}$
(5)
$\displaystyle\left[\frac{1}{\sigma_{y}\sqrt{2\pi}}e^{-\frac{1}{2}(\frac{\sqrt{\Delta^{2}-x^{2}}}{\sigma_{y}})^{2}}+\frac{1}{\sigma_{y}\sqrt{2\pi}}e^{-\frac{1}{2}(\frac{-\sqrt{\Delta^{2}-x^{2}}}{\sigma_{y}})^{2}}\right]$

A detailed derivation of Eq. 5 can be found in the Appendix C. We applied Eq.
5 to Eq. 1 to fit the experimental data. Since the total number of magnetic
species $N$ is solely determined by the size and chemical composition of each
sample, the two remaining free-fitting parameters are just $\sigma_{a}$ and
$\sigma_{b}$.

Figure 5:  (a)-(f) Upper panels plot the low-$T$ $4f$ heat capacity
$C_{p}^{4f}$ of Tm_1-xY_xVO4 with $x=0.25,0.45,0.89,0.97,0.99$ and $0.997$.
The blue dashed line shows a forced single-gap Schottky fitting with the data.
The green dashed line in (a) represents the fitting result with Eq. 1, with
$\rho(\Delta)$ given by a Gaussian PDF. The red dashed lines indicate fit to
the 2D-Gaussian distribution in Eq. 5. Lower panels of (a-f) plot the residual
$\delta^{4f}/R$ (fitting deviations) of each 2D Gaussian fit.

Fig. 5 contrasts the fitting based on Eq. 5 (red dashed lines) and forced
fitting of a single Schottky anomaly (blue dashed lines) with $C^{4f}_{p}/R$
(black empty circles) for a few representative compositions $x$. The residual
of the 2D Gaussian fits and the data $\delta^{4f}$ are also shown for each $x$
underneath each of the main panels. We note that the single Schottky anomaly
scenario poorly fits the data for all $x$. In contrast, the 2D Gaussian
distribution is found to describe the data very well for intermediate
concentrations ($x=0.25-0.97$). However, with further increase in $x$, the
distribution progressively deviates from the 2D Gaussian description, and
becomes worst for the most dilute Tm concentrations, as can be seen in the
large and systematic deviation from zero of $\delta^{4f}$ in the lower panels
of Fig. 5(e) and (f).

Figure 6:  Evolution of the 2D-Gaussian distribution parameters plotted in the
$T-x$ plane. The blue and red symbols represent $\sigma_{a}$ and $\sigma_{b}$
extracted from the fitting to Eq. 5 (against the left axis). $T_{Q}$ is
represented by the green diamond symbols. The gaps of a forced single Schottky
anomaly fit are shown as gray symbols. The bottom axis is shown both as
$(1-x)^{-1/3}$ and as $d_{Tm-Tm}$ (Wigner-Seitz radius of the Tm^3+). The top
axis indicates the $x$ value of each composition. The Root-Mean-Square-Error
(RMSE) of each 2D Gaussian fitting is shown as yellow empty circles against
the right axis.

Fig. 6 summarizes the evolution with composition of the fitting parameters,
and the associated goodness of fitting (expressed as the root-mean-square
error, denoted by RMSE), as a function of the average distance between Tm ions
$d_{Tm-Tm}$. The abscissa is shown using a log scale, and is also expressed
equivalently in terms of the composition as $(1-x)^{-1/3}$. The corresponding
Y concentration $x$ is also labeled on the top axis. The two fitting
parameters in the 2D Gaussian model, $\sigma_{a}$, and $\sigma_{b}$, which
describe the width of the Gaussian distribution along the two-dimensional
space, are shown as blue and red solid circles respectively. To provide a
crude estimate of the characteristic energy scale associated with the strain-
induced splitting, the result of the forced single-gap Schottky fit is also
included (gray squares), though readers are reminded that this single-gap
model provides a very poor fit, strongly suggesting the presence of a
spatially inhomogeneous distribution of gaps. $T_{Q}$ is also shown on this
plot as green diamonds. The yellow circles on the right axis indicate the RMSE
of the respective 2D Gaussian fits.

For compositions $x$ beyond approximately 50%, the values of the gap obtained
from the single Schottky fit and of the widths of the two independent Gaussian
distributions, $\sigma_{a}$ and $\sigma_{b}$, rapidly decline in magnitude.
This is in keeping with the general expectation that as the average distance
between Tm^3+ ions increases (separated by increasingly undistorted YVO4), the
energy scale associated with lattice distortions gets progressively smaller,
and the associated distribution gets narrower. However, $\sigma_{b}$ nearly
vanishes for large enough $x$, resulting in physically unsound gap
distributions $\rho(\Delta)$. Moreover, the systematic increase of RMSE with
$x$ in the dilute limit reveals a significant deviation of the gap
distribution from a 2D Gaussian description. This observation is robust
against experimental uncertainties, including details of the subtraction of
the phonon contribution to the heat capacity. This evidence for a progressive
deviation of the local random strain distribution from a 2D Gaussian is our
primary observation.

##  V Discussion

According to the central limit theorem, for large enough sample sizes, the
normalized sum of independent random variables tends towards a Gaussian PDF.
In the present case, with potentially two independent variables (i.e. the two
relevant components of the strain tensor, $\epsilon_{B_{1g}}$ and
$\epsilon_{B_{2g}}$ experienced at each Tm^3+ site) a truly random
distribution with no correlations would yield a $\rho(\Delta)$ following a 2D
Gaussian distribution (see Fig. 1). Our observations reveal that for a wide
range of intermediate concentrations, the heat capacity in Tm_1-xY_xVO4 cannot
be distinguished from such a double-Gaussian distribution. Consequently, it is
tempting to assume that a random field approach is justified for describing
such a system. However, the observation of substantial deviations from the 2D
Gaussian distribution for dilute concentrations indicates that a more subtle
perspective is appropriate, and indeed that strains arising from substitution
should not be treated as uncorrelated. In the following, we propose a scenario
for the observed crossover of $\rho(\Delta)$ from a double Gaussian
distribution in the concentrated Tm regime to a significant deviation from
that in the dilute Tm regime, which emphasizes the presence of strain
correlations in both limits.

Figure 7: Comparison of randomly scattered strain sources and their impact on
the level-splitting distribution in Tm_1-xY_xVO4. Each row shows a single
realization of random level-splittings, $\Delta$, in a two-dimensional
periodic square lattice of size $256\times 256$ on the left, and the
corresponding measured distribution of $\Delta$, $\rho(\Delta)$, on the right.
These splittings are normalized by each realization’s maximum,
$\Delta_{\max}$. In every case, the $B_{1g}$ and $B_{2g}$ strains are randomly
drawn, and then Eq. 2 is used to compute $\Delta$. The measured histograms are
created with $256^{2}=65\,536$ values of $\Delta$. The first and second rows,
plots (a,b) and (c,d), show the result of randomly placing defects throughout
the simulation and computing the strains from elasticity theory, as explained
in Appendix D. In this case, the defects are edge dislocations. (a,b)
represent a low concentration of defects (high $x$) whereas (c,d) represent a
high concentration of defects (low $x$). The bottom row, (e,f), shows the case
where both strains are Gaussian random and spatially uncorrelated, with a
resulting histogram in the form in Eq. 5. As one introduces more crystalline
defects into the system, one traverses the columns, and the distribution of
level-splittings experienced by the $\ce{Tm^{3+}}$ ions approaches the 2D
Gaussian.

Each Tm ion in Tm_1-xY_xVO4 acts as a local reporter of the strain, with a
splitting of the CEF ground state doublet that depends on the quadrature sum
of the local values of $\epsilon_{B_{1g}}$ and $\epsilon_{B_{2g}}$. These
strains originate from a variety of deviations from perfect crystallinity
throughout the solid, including from extended defects such as dislocations,
from residual strains that the crystal experiences from external stresses, and
of course from the chemical substitution itself. To illustrate an important
point, we briefly consider a thought experiment in which defects are
introduced to TmVO4 without diluting the Tm ions. We also imagine that the
strains are small such that the response is strictly linear, and consequently
that strains from various sources are simply superimposed.

The outcome of this thought experiment is illustrated in Fig. 7. In the
Figure, we modeled the effect of adding edge dislocations randomly throughout
a two-dimensional periodic square lattice. Using well-known expressions for
the strain fields generated by a single disclocation, we compute
$\epsilon_{B_{1g}}$ and $\epsilon_{B_{2g}}$ at each lattice site, from which
we calculate the level-splitting experienced by each $\ce{Tm^{3+}}$ ion. To
avoid the formally diverging strain at the dislocation cores that emerge in
linear elasticity, we place the dislocations at the centers of the square
plaquettes whereas the two-level ions are positioned on each lattice site. For
additional details of the model, please see Appendix D. We choose edge
dislocations as a tractable proof-of-concept to simulate the effects of long-
range strain relaxation arising from a variety of actual defects. The real
situation is of course more complex, with large local strains that presumably
exceed the regime of linear response, but this simplification allows us to
develop an important idea.

The strains associated with each defect relax as a power-law rather than
exponential function of the distance from the defect, and hence there is no
characteristic length scale for the relaxation. Nevertheless, there is an
effective length scale beyond which the strains associated with any given
defect become sufficiently small to be indistinguishable from any slowly
varying background strain that arises not from intrinsic defects but, for
instance, from small external stresses experienced by the material. This
defines an effective radius over which the Tm ions are able to record the
strains arising from defects.

In the dilute limit, most $\ce{Tm^{3+}}$ ions report few strain sources, as
there are only a few defects generating the internal strains[Figs. 7(a,b)].
Thus, the majority of the system is essentially undistorted, as indicated by
the dark region in panel (a). It is clear that the aggregate strains on
adjacent sites are highly correlated, as they are dominated by the strains
generated by the defects. Although the locations of the defects are random,
the number of defects within the “sensing radius” of any given ion is small
enough that the central limit theorem does not apply.

We now imagine increasing the concentration of defects within this
hypothetical system [Figs. 7(c,d)]. The defects are positioned randomly, so
there is no correlation between their locations. Since strains are correlated
over long distances in the dilute limit, there must still be significant
correlations between strains on adjacent sites even in the more concentrated
regime. Nevertheless, as the concentration of defects is progressively
increased, any given Tm site experiences strains arising from an increasingly
large number of defects within its sensing volume, which, in the elastic
limit, are simply added to give the total (local) strain that the site
experiences. Since the dominant effect of the defect locations is random
(uncorrelated), the distribution $\rho(\Delta)$ then starts to crossover into
the 2D Gaussian regime as the central limit theorem locally begins to dominate
the strain statistics in the $B_{1g}$ and $B_{2g}$ strain channels. It is
important to emphasize, meanwhile, that one can obtain a $\rho(\Delta)$
distribution that approximates a 2D Gaussian even if the strains are
correlated between sites. The extreme limit of short-ranged behavior would be
realized only in samples with a high density of defects, as shown in Figs.
7(e,f).

Thus, within this thought experiment, a cross-over is anticipated from a
regime in which clear signatures of strain correlations are apparent (in the
dilute limit) to a regime in which the randomness of the defect locations
dominates the strain distribution, and for which a Gaussian distribution of
gaps is a reasonable approximation. The concentration of defects at which this
crossover occurs is determined by the interplay of several key length scales:
the effective length scale on which strains relax, the distance between
defects, and the distance between Tm “reporters”. Signatures of “randomness”
only become apparent when the number of defects and the number of reporters
within a typical relaxation volume become sufficiently large for the central
limit theorem to apply.

The situation in Tm_1-xY_xVO4 is slightly more complex than the above thought
experiment, but not in any qualitative way. The substitution of Y for Tm
yields a progressive increase in the distance between Tm “reporters”, while
also affecting the strain landscape. Similarly, the length scale over which
strains relax might itself depend on the concentration of impurities. The
interplay between length scales is still apparent, though, and dilute versus
concentrated regimes with different qualitative behaviors can still be
anticipated.

In addition to the above, we note that the constitutive elasticity relations
(i.e. the relations between stress and strain, which at the level of linear
response are governed by the elastic stiffness tensor) necessarily imply that
some arbitrary stress state at a given site is a linear superposition of the
various strain components at that site. Therefore, if the stress state is
fixed by sources of structural disorder in the material, the six strain
components cannot all be taken as independent of one another. Hence,
$\epsilon_{B_{1g}}$ and $\epsilon_{B_{2g}}$ on any given site are correlated
in some specific way determined by the elastic distortion of the crystal. We
emphasize that this is in addition to the correlations that exist for the same
strain component between adjacent sites.

##  VI Conclusions

Substitution of Y for Tm in Tm_1-xY_xVO4 rapidly suppresses the long-range
quadrupole order that is found for the pristine “parent” compound TmVO4.
Splitting of the CEF non-Kramers groundstate doublet of each Tm ion acts as a
local reporter of the quadrature sum of $\epsilon_{B_{1g}}$ and
$\epsilon_{B_{2g}}$ strains on each site. Our heat capacity data reveal that
the splitting is not uniform, but is described by some probability
distribution. While the heat capacity data are not amenable to an inverse
analysis (in which the distribution is determined from the observed heat
capacity), the data can nevertheless be compared to the simple expectations of
a 2D Gaussian distribution that would otherwise be expected if the strains
were truly random (uncorrelated). Significantly, for dilute concentrations the
heat capacity cannot be accounted for by such a distribution, revealing for
the first time clear evidence of strain correlations between sites and between
different strain components at the same site. For larger concentrations,
however, the heat capacity cannot be distinguished from that which would be
obtained from a random distribution. The cross-over between these regimes can
be understood in terms of the progressive change in the number of defects
within the sensing volume of a given Tm ‘reporter’ as the concentration of
impurities $x$ is changed. In the concentrated regime, even though strains are
correlated between sites, the defects that generate the strains themselves are
randomly distributed. Hence, the total strain on each site is reasonably
described by a Gaussian distribution when the number of defects within the
sensing volume of a reporter ion is large enough for the central limit theorem
to apply.

Our results demonstrate that the strains, which arise from chemical
substitution in Tm_1-xY_xVO4, cannot be regarded as uncorrelated random
fields. Consequently, theoretical treatment of this problem (i.e. treatment of
the problem of quadrupole moments that interact in the presence of non-uniform
longitudinal and transverse effective fields that arise from local chemical
substitution) cannot rely solely on the results of the random field Ising
model. Instead, one must consider the subtle effects of strain correlations
that can occur between different sites for the same symmetry channel, as well
as between different symmetry channels, either between different sites or the
same one. This problem is not restricted to just Tm_1-xY_xVO4. Any electronic
nematic system for which chemical substitution is used to traverse the
associated phase diagram falls into the same paradigm. Our results therefore
motivate in-depth theoretical efforts to model the rich problem of correlated
“random” fields on Ising electronic nematic order.

##  VII Acknowledgements

We thank Brad Ramshaw, Patrick Hollister, Jörg Schmalian, Mingde Jiang, and
Evan Reed for fruitful discussions. Low-temperature heat-capacity measurements
performed at Stanford University were supported by the Air Force Office of
Scientific Research Award FA9550-20-1-0252, using a cryostat acquired with
award FA9550-22-1-0084. Crystal-growth experiments were supported by the
Gordon and Betty Moore Foundation Emergent Phenomena in Quantum Systems
Initiative Grant GBMF9068. M.P.Z. was also partially supported by a National
Science Foundation Graduate Research Fellowship under grant number
DGE-1656518. L.Y. acknowledges support from Marvin Chodorow Postdoctoral
Fellowship at the Department of Applied Physics, Stanford University. The
theoretical work at the University of Minnesota (W.J.M. and R.M.F.) was
supported by the U.S. Department of Energy, Office of Science, Basic Energy
Sciences, Materials Science and Engineering Division, under Award No. DE-
SC0020045. DFT calculations by Y.Z. were supported by NSF DMR-1455050.

##  Appendix A Determination of chemical composition via magnetic
susceptibility

The magnetic susceptibility ($\chi$) of TmVO4 along its $c$-axis exhibits a
Curie-Weiss behavior down to $2K$ Cooke  _et al._ (1972), which is not
frequency-dependent in the absence of an external magnetic field Bleaney and
Wells (1980). We observed the same $\chi$ behavior in (Y,Tm)VO4 above $T_{Q}$
(if it exists) owing to the identical crystal field environments throughout
the doping series. Consequently, to accurately identify the chemical
composition of samples that contain the lowest Tm percentages, we compare
their magnetic susceptibility with the calculated specific magnetic
susceptibility. In zero-field conditions, the calculated magnetic
susceptibility can be expressed as:

$\chi=\frac{\frac{1}{4}\mu_{0}Ng_{C}^{2}\mu_{B}^{2}}{k_{B}(T-T_{\theta})}$ (6)

$\mu_{0}$ is the vacuum magnetic permeability. $N$ is number of Tm^3+ ions.
$g_{C}=10.22$ is the g-factor of $\ce{Tm^{3+}}$ along the c-axis. $\mu_{B}$ is
the Bohr magneton. $T_{\theta}$ is the Curie temperature. We fit the equation
above with the measured ac magnetic susceptibility as:

$\chi=\frac{A}{k_{B}(T-T_{\theta})}+\chi_{0}$ (7)

Three fitting parameters are $A=\frac{\mu_{0}g_{c}^{2}\mu_{B}^{2}N}{4k_{B}}$
as the Curie constant, $T_{\theta}$ as the Curie temperature, and $\chi_{0}$
as the measured constant magnetic background from non-magnetic species (mostly
from YVO4 lattices).

The ac magnetic susceptibility measurement was performed with a Magnetic
Properties Measurement System (MPMS), with a drive frequency of 75.7Hz. The
value of $A$ takes the average of all the fitting results from 2K to different
$T_{max}$ values. Finally, we obtain $x$ in the main text calculated from $N$,
molar mass $M$, and the sample mass $m$.

##  Appendix B Specific Heat Model and Determination of $C_{p}^{4f}$

In this section, we discuss the method of finding and subtracting the phonon
background and obtaining the heat capacity per Tm site ($C_{p}^{4f}$) from the
heat capacity measurement. We assume the measured heat capacity is only the
sum of the phonon contribution and the Schottky anomaly of the $4f$ degree of
freedom:

$\displaystyle C_{p}^{m}$ $\displaystyle=C_{p}^{Debye}+C_{p}^{Schottky}$ (8)
$\displaystyle=C_{p}^{Debye}+(1-x)C_{p}^{4f}$

Where $x$ represents the content of $Y^{3+}$. The three-dimensional Debye
model describes the phonon contribution of heat capacity. The fit equation is:

$C_{p}^{Debye}=9Nk_{B}(\frac{T}{T_{eff}})^{3}\int_{0}^{\frac{T}{T_{eff}}}dx\frac{x^{4}e^{x}}{(e^{x}-1)^{2}}$
(9)

The $C_{p}^{Debye}$ refers to the contribution of phonons to the total heat
capacity. $T_{eff}$ is a fitting parameter that represents the effective Debye
temperature. It is important to note that the fitting method used in the Debye
Model with a limited temperature range is empirical and does not take into
account the anisotropic phonon dispersion in the Y_xTm_1-xVO4 series.
Therefore, we use $T_{eff}$ to differentiate it from the experimentally
determined $T_{Debye}$ that was reported elsewhereDenisova  _et al._ (2015).

We applied Eq.3 to heat capacity data from $10$ to $25K$, where the Schottky
contribution is negligible.

To obtain $\Delta C_{p}^{Schottky}$, we calculated the phonon background at
lower temperatures (shown as green and red dashed lines in Fig.8), using
$T_{eff}$ obtained from Debye fitting. The phonon background below $T=0.65K$
was extrapolated using the low-temperature $T^{3}$ limit of the Debye model:

$C_{p}^{Debye}=36Nk_{B}(\frac{T}{T_{eff}})^{3}$ (10)

The fit result based on the specific heat model is plotted in Fig. 8 (j).

Figure 8: (a)-(i) Fitting result of the heat capacity data in the fitting
range from 10K to 25K. The dashed curve plots the fitting equation based on
the Debye Model. The green (red) dashed line region is based on the
calculation of the Debye model (expansion of the low-temperature limit)
extrapolated from the fit result. (j) The effective Debye Temperature
$T_{eff}$ v.s. $x$ obtained from the fitting.

##  Appendix C Derivation of the 2D-Gaussian distribution in (Y,Tm)VO4

The ground state doublet of ${\rm Tm}^{3+}$ affected by local strain is
described by Eq.4 and Eq.5 in the main text. We assume that $\rho_{x}(x)$ and
$\rho_{y}(y)$ are two probability distribution functions corresponding to the
strength of the $\epsilon_{B1g}$ and $\epsilon_{B2g}$ respectively. Qualified
by the uncorrelated random-field assumption, both of the functions are the
Gaussian Distribution with a mean value centered at $0$. The probability
distribution function of the quadrature sum $\Delta$ can be written as:

$\rho(\Delta)=\int_{-\infty}^{\infty}{\int_{-\infty}^{\infty}{dxdy\rho_{x}(x)\rho_{y}(y)\delta(f(\Delta))}}$
(11)

where $f(\Delta)$ is:

$f(\Delta)=\Delta-\sqrt{x^{2}-y^{2}}$ (12)

Therefore, we can now solve for $\rho(\Delta)$ as:

$\displaystyle\rho(\Delta)=$
$\displaystyle\int_{-\infty}^{\infty}{\int_{-\infty}^{\infty}{dxdy\rho_{x}(x)\rho_{y}(y)\delta(\Delta-\sqrt{x^{2}-y^{2}}})}$
(13) $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}dx{\rho_{x}(x)[\int_{-\infty}^{\infty}{dy\rho_{y}(y)\delta(\Delta-\sqrt{x^{2}-y^{2}})}}]$

The identity of a delta function implies:

$\delta[g(y)]=\sum_{i}\frac{\delta(y-y_{i})}{|g^{\prime}(y_{i})|}$ (14)

where $m_{i}$ are roots of the function $g(y)$.

$\displaystyle g(y)$ $\displaystyle=\Delta-\sqrt{x^{2}-y^{2}}$ (15)
$\displaystyle y_{i}$ $\displaystyle=\pm\sqrt{\Delta_{i}^{2}-x^{2}}$
$\displaystyle g^{\prime}(y_{i})$
$\displaystyle=\frac{\partial(\Delta-\sqrt{x^{2}-y^{2}})}{\partial{(\pm\sqrt{\Delta_{i}^{2}-x^{2}})}}=\pm\frac{\sqrt{\Delta_{i}^{2}-x^{2}}}{\Delta}$

Therefore, we now apply Eq.14 to Eq.13 and obtain:

$\displaystyle\rho(\Delta)$
$\displaystyle=\int_{-\Delta}^{\Delta}dx\rho_{x}(x)\frac{\Delta}{\sqrt{\Delta^{2}-x^{2}}}$
(16)
$\displaystyle[\rho_{y}(\sqrt{\Delta^{2}-x^{2}})+\rho_{y}(-\sqrt{\Delta^{2}-x^{2}})]$

Then, we apply two Gaussian distributions with widths $\sigma_{x}$ and
$\sigma_{y}$ respectively, and calculate the analytical solution of the
$\rho(\Delta)$:

$\displaystyle\rho_{x}(x)$
$\displaystyle=\frac{1}{\sigma_{x}\sqrt{2\pi}}e^{-\frac{1}{2}(\frac{x}{\sigma_{x}})^{2}}$
(17) $\displaystyle\rho_{y}(y)$
$\displaystyle=\frac{1}{\sigma_{y}\sqrt{2\pi}}e^{-\frac{1}{2}(\frac{y}{\sigma_{y}})^{2}}$

The analytical solution for $\rho(\Delta)$ is the Eq.5 of the main text.

##  Appendix D Modeling correlated strains from dislocations

In this section, we describe the numerical process used to generate the level-
splitting realizations shown in Fig. 7. We use the same procedure shown in
Meese  _et al._ (2024). To best compare conventional short-ranged quenched
disorder fields with long-ranged elastic strains generated by defects, we use
the longest-ranged fundamental crystalline defect typically observed in
equilibrium, the dislocation Landau and Lifshitz (1970); Mura (1987).
Dislocations are topological defects with a charge known as the Burgers
vector, which we denote as $\boldsymbol{b}$ Chaikin and Lubensky (1995). The
strains are well known to decay as $1/r$ from the core of a straight
dislocation, where $r$ is the radial distance to the dislocation line Landau
and Lifshitz (1970); Mura (1987). In brief, we employ an ensemble of random,
uncorrelated, edge dislocations to generate the $B_{1g}$ and $B_{2g}$
components of the strain tensor. We use linear elasticity theory then to
superimpose the strains from each dislocation within the ensemble at a given
site. From there it follows that the level-splitting, $\Delta$, is the
quadrature sum of the net strains at each site.

In isotropic media, the displacement vector $\boldsymbol{u}$ induced by a
single straight edge dislocation centered at the origin, oriented along the
$z$-axis, and endowed with Burgers vector $\boldsymbol{b}$ is known to be
Landau and Lifshitz (1970); Mura (1987)

$\displaystyle u_{x}(\boldsymbol{r};\boldsymbol{b})$
$\displaystyle=\frac{b_{x}\phi}{2\pi}+\frac{1-2\nu}{4\pi(1-\nu)}b_{y}\ln r$
$\displaystyle\quad+\frac{1}{8\pi(1-\nu)}\left(b_{x}\sin 2\phi-b_{y}\cos
2\phi\right),$ (18) $\displaystyle u_{y}(\boldsymbol{r};\boldsymbol{b})$
$\displaystyle=\frac{b_{y}\phi}{2\pi}-\frac{1-2\nu}{4\pi(1-\nu)}b_{x}\ln r$
$\displaystyle\quad-\frac{1}{8\pi(1-\nu)}\left(b_{x}\cos 2\phi+b_{y}\sin
2\phi\right),$ (19)

where $\boldsymbol{r}=(x,y)=r(\cos\phi,\sin\phi)$ and $\nu$ is the Poisson
ratio. From these expressions, $\epsilon_{B_{1g}}$ and $\epsilon_{B_{2g}}$ are
obtained by differentiating the displacement vector as

$\displaystyle\epsilon_{B_{1g}}(\boldsymbol{r};\boldsymbol{b})$
$\displaystyle=\partial_{x}u_{x}-\partial_{y}u_{y}=-\frac{\left(\boldsymbol{b}\cdot\boldsymbol{\hat{r}}\right)\sin
2\phi}{2\pi(1-\nu)r},$ (20)
$\displaystyle\epsilon_{B_{2g}}(\boldsymbol{r};\boldsymbol{b})$
$\displaystyle=\partial_{x}u_{y}+\partial_{y}u_{x}=\frac{\left(\boldsymbol{b}\cdot\boldsymbol{\hat{r}}\right)\cos
2\phi}{2\pi(1-\nu)r}.$ (21)

These functions are used to compute the strain generated by a dislocation at
every lattice site in a two-dimensional periodic square grid with properly
enforced periodic boundary conditions Meese  _et al._ (2024). To avoid the
singular behavior of the strain at the dislocation cores, we place them at the
centers of the square plaquettes.

The total strain in the $a\in\\{B_{1g},B_{2g}\\}$ irreducible representation
is then computed from a set of $N_{d}$ dislocations as

$\displaystyle\epsilon_{a}(\boldsymbol{r})=\sum_{\gamma=1}^{N_{d}}\epsilon_{a}\left(\boldsymbol{r}-\boldsymbol{r}_{\gamma};\boldsymbol{b}_{\gamma}\right),$
(22)

from which it follows that the level-splittings are

$\displaystyle\Delta(\boldsymbol{r})=\sqrt{\eta_{1}^{2}\epsilon_{B_{1g}}^{2}(\boldsymbol{r})+\eta_{2}^{2}\epsilon_{B_{2g}}^{2}(\boldsymbol{r})}.$
(23)

For simplicity, we have maintained the isotropic approximation, and have set
the magneto-elastic couplings to be equal: $\eta_{1}=\eta_{2}=\lambda$. In
Figure 7, we only consider ratios of $\Delta(\boldsymbol{r})$ to its maximum
value in the lattice, $\Delta_{\max}$, which removes the numerical dependence
on $\lambda$ and $\nu$ without losing the qualitative features of a long-
ranged correlated random strain field.

Each realization of random strain is created via an ensemble of $N_{d}$ random
dislocations. Equilibrium distributions of dislocations must be “dislocation
neutral” Kleinert (1989). Thus, we impose that the dislocations within each
realization appear in pairs of equal and opposite topological charges. This is
the only nonlinear effect we model, however. Thus, dislocations were assumed
to be nearly uncorrelated, up to the condition that the total Burgers vector
vanishes and no two dislocations occupy the same plaquette. Up to these
restrictions, the edge dislocations were created by sampling the Burgers
vectors uniformly from the lattice vectors
$\left\\{\pm\boldsymbol{\hat{x}},\pm\boldsymbol{\hat{y}}\right\\}$ and
uniformly sampling the core locations to be at the center of one of the
$L^{2}$ square plaquettes.

The calculation was performed with the Julia programming language Bezanson
_et al._ (2017) and Fig. 7 was created using the Makie plotting package
Danisch and Krumbiegel (2021).

## References

  * Vojta (2019) T. Vojta, Annual Review of Condensed Matter Physics 10, 233 (2019).
  * Ran  _et al._ (2011) S. Ran, S. L. Bud’ko, D. K. Pratt, A. Kreyssig, M. G. Kim, M. J. Kramer, D. H. Ryan, W. N. Rowan-Weetaluktuk, Y. Furukawa, B. Roy, A. I. Goldman, and P. C. Canfield, Physical Review B 83, 144517 (2011).
  * Dioguardi  _et al._ (2016) A. Dioguardi, T. Kissikov, C. Lin, K. Shirer, M. Lawson, H.-J. Grafe, J.-H. Chu, I. Fisher, R. Fernandes, and N. Curro, Physical Review Letters 116, 107202 (2016).
  * Forrest  _et al._ (2016) T. R. Forrest, P. N. Valdivia, C. R. Rotundu, E. Bourret-Courchesne, and R. J. Birgeneau, Journal of Physics: Condensed Matter 28, 115702 (2016).
  * Carlson  _et al._ (2006) E. W. Carlson, K. A. Dahmen, E. Fradkin, and S. A. Kivelson, Physical Review Letters 96, 097003 (2006).
  * Carlson  _et al._ (2014)
  * Nie  _et al._ (2014) L. Nie, G. Tarjus, and S. A. Kivelson, Proceedings of the National Academy of Sciences 111, 7980 (2014).
  * Ye  _et al._ (2022) X. Ye, R. Narayanan, and T. Vojta, Phys. Rev. B 105, 024201 (2022).
  * Meese  _et al._ (2022) W. J. Meese, T. Vojta, and R. M. Fernandes, Phys. Rev. B 106, 115134 (2022).
  * Wiecki  _et al._ (2021) P. Wiecki, R. Zhou, M.-H. Julien, A. E. Böhmer, and J. Schmalian, Physical Review B 104, 125134 (2021).
  * Landau and Lifshitz (1970) L. D. Landau and E. M. Lifshitz,  _Theory of Elasticity_ , 2nd ed., Course of Theoretical Physics No. 7 (Pergamon Press, 1970).
  * Mura (1987) T. Mura,  _Micromechanics of defects in solids_ , 2nd ed., Mechanics of elastic and inelastic solids; 3 (M. Nijhoff, Dordrecht, The Netherlands, 1987).
  * Maharaj  _et al._ (2017) A. V. Maharaj, E. W. Rosenberg, A. T. Hristov, E. Berg, R. M. Fernandes, I. R. Fisher, and S. A. Kivelson, Proceedings of the National Academy of Sciences 114, 13430 (2017).
  * Massat  _et al._ (2022) P. Massat, J. Wen, J. M. Jiang, A. T. Hristov, Y. Liu, R. W. Smaha, R. S. Feigelson, Y. S. Lee, R. M. Fernandes, and I. R. Fisher, Proceedings of the National Academy of Sciences 119, e2119942119 (2022).
  * Morin and Schmitt (1990) P. Morin and D. Schmitt, “Quadrupolar interactions and magneto-elastic effects in rare earth intermetallic compounds ch 1,” (North-Holland, Netherlands, 1990).
  * Wortman  _et al._ (1974) D. Wortman, R. Leavitt, and C. Morrison, Journal of Physics and Chemistry of Solids 35, 591 (1974).
  * Melcher (1976) R. Melcher (Academic Press, 1976) pp. 1–77.
  * Imry and Ma (1975) Y. Imry and S.-k. Ma, Phys. Rev. Lett. 35, 1399 (1975).
  * Aharony  _et al._ (1976) A. Aharony, Y. Imry, and S.-k. Ma, Phys. Rev. Lett. 37, 1364 (1976).
  * Grinstein (1976) G. Grinstein, Phys. Rev. Lett. 37, 944 (1976).
  * Nattermann (1998) T. Nattermann, in  _Spin glasses and random fields_ (World Scientific, 1998) pp. 277–298.
  * Binder (1983) K. Binder, Zeitschrift für Physik B Condensed Matter 50, 343 (1983).
  * Gehring  _et al._ (1976) G. A. Gehring, S. J. Swithenby, and M. R. Wells, Solid State Communications 18, 31 (1976).
  * Zachar and Zaliznyak (2003) O. Zachar and I. Zaliznyak, Phys. Rev. Lett. 91, 036401 (2003).
  * Feigelson (1968) R. Feigelson, Journal of the American Ceramic Society 51, 538 (1968).
  * Smith and Wanklyn (1974) S. Smith and B. Wanklyn, Journal of Crystal Growth 21, 23 (1974).
  * Oka  _et al._ (2006) K. Oka, H. Unoki, H. Shibata, and H. Eisaki, Journal of crystal growth 286, 288 (2006).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Computational Materials Science 6, 15 (1996).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Phillips (1987) W. A. Phillips, Reports on Progress in Physics 50, 1657 (1987).
  * van Duijn  _et al._ (2005) J. van Duijn, K. H. Kim, N. Hur, D. Adroja, M. A. Adams, Q. Z. Huang, M. Jaime, S.-W. Cheong, C. Broholm, and T. G. Perring, Phys. Rev. Lett. 94, 177201 (2005).
  * Tachibana  _et al._ (2007) M. Tachibana, Y. Kohama, T. Atake, and E. Takayama-Muromachi, Journal of applied physics 101, 09D502 (2007).
  * Martin  _et al._ (2017) N. Martin, P. Bonville, E. Lhotel, S. Guitteny, A. Wildes, C. Decorse, M. C. Hatnean, G. Balakrishnan, I. Mirebeau, and S. Petit, Physical Review X 7, 041028 (2017).
  * Markovsky (1988) A. Markovsky, Applied Mechanics Reviews 41, 247 (1988).
  * Yuan and Zhang (2019) D. Yuan and X. Zhang, SN Applied Sciences 1, 1178 (2019).
  * Cooke  _et al._ (1972) A. Cooke, S. Swithenby, and M. Wells, Solid State Communications 10, 265 (1972).
  * Bleaney and Wells (1980) B. Bleaney and M. R. Wells, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 370, 131 (1980).
  * Denisova  _et al._ (2015) L. T. Denisova, L. G. Chumilina, Y. F. Kargin, A. D. Izotov, and V. M. Denisov, Doklady Physical Chemistry 463, 173 (2015).
  * Meese  _et al._ (2024) W. J. Meese, Y. Li, L. Ye, M. P. Zic, I. R. Fisher, and R. M. Fernandes, In Preparation (2024).
  * Chaikin and Lubensky (1995) P. M. Chaikin and T. C. Lubensky,  _Principles of condensed matter physics_ (Cambridge University Press, Cambridge, 1995).
  * Kleinert (1989) H. Kleinert, in  _Part 3: Gauge Fields in Solids_, Gauge Fields in Condensed Matter, Vol. 2 (WORLD SCIENTIFIC, Germany, 1989).
  * Bezanson  _et al._ (2017) J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, SIAM Review 59, 65 (2017).
  * Danisch and Krumbiegel (2021) S. Danisch and J. Krumbiegel, Journal of Open Source Software 6, 3349 (2021).