Stress correlations in near-crystalline packings

Roshan Maharana$\star$ and Kabir Ramola$\dagger$

Tata Institute of Fundamental Research, Hyderabad 500107, India

⋆<EMAIL_ADDRESS>†<EMAIL_ADDRESS>
## Abstract

We derive exact results for stress correlations in near-crystalline systems in
two and three dimensions. We study energy minimized configurations of
particles interacting through Harmonic as well as Lennard-Jones potentials,
for varying degrees of microscopic disorder and quenched forces on grains. Our
findings demonstrate that the macroscopic elastic properties of such near-
crystalline packings remain unchanged within a certain disorder threshold, yet
they can be influenced by various factors, including packing density,
pressure, and the strength of inter-particle interactions. We show that the
stress correlations in such systems display anisotropic behavior at large
lengthscales and are significantly influenced by the pre-stress of the system.
The anisotropic nature of these correlations remains unaffected as we increase
the strength of the disorder. Additionally, we derive the large lengthscale
behavior for the change in the local stress components that shows a $1/r^{d}$
radial decay for the case of particle size disorder and a $1/r^{d-1}$ behavior
for quenched forces introduced into a crystalline network. Finally, we verify
our theoretical results numerically using energy-minimised static particle
configurations.

###### Contents

  1. 1 Introduction
  2. 2 Models
    1. 2.1 Short-ranged repulsive Harmonic interaction
    2. 2.2 Attractive Lennard-Jones interaction with cut-off
  3. 3 Numerical Simulations
    1. 3.1 Harmonic Model
    2. 3.2 LJ Model
  4. 4 Elastic properties of near-crystalline packings
  5. 5 Displacement fields induced by microscopic disorder
    1. 5.1 Short-ranged repulsive harmonic interaction
    2. 5.2 Attractive Lennard-Jones interaction with cut-off
  6. 6 Stress correlations induced by microscopic disorder
    1. 6.1 Local stress fluctuations
    2. 6.2 Stress correlations in two dimensional systems
    3. 6.3 Comparision with amorphous systems at large lengthscales
    4. 6.4 Continuum limit
    5. 6.5 Stress correlations in three dimensions
  7. 7 Response to a point force
  8. 8 Distribution of Stresses
  9. 9 Discussion and Conclusion
  10. A Green’s functions in the Harmonic model
    1. A.1 Green’s functions for displacement fields
    2. A.2 Green’s functions for change in local stresses in Fourier space
  11. B Angular variation of stress correlations at large lengthscales
  12. C Hammer projection

##  1 Introduction

Jammed athermal materials find relevance in various fields such as soft
condensed matter physics, material science, civil engineering and metallurgy
[1, 2]. Additionally, jammed packings also arise in fields such as biophysics,
where cellular tissues are well described by soft potential models [3, 4].
They arise when it is not feasible to achieve true thermodynamic equilibrium.
The stability of athermal solids against mechanical disturbances can be
attributed to the macroscopic rigidity arising from the network of constituent
particles [5, 6, 7, 8, 9]. This collective elasticity arises in any system of
interacting particles at low temperatures and is observed universally in both
crystalline and amorphous structures of athermal solids [10, 11, 12, 13, 14,
15, 16, 17, 18, 19, 20, 21, 22, 23, 24]. The reference states that make up the
collection of amorphous solids are highly dependent on the preparation method,
and each configuration satisfies the conditions of local equilibrium i.e., the
force and torque balance of each constituent. Amorphous structures, while
stable in a local sense, are typically not the lowest energy states of their
individual components [25, 26, 27]. Consequently, jammed packings of soft
particles can exhibit both amorphous as well as crystalline structures.

Near-crystalline materials demonstrate a range of unique properties and serve
as a bridge between the physics of crystals and amorphous materials, providing
valuable insights into the behavior of athermal ensembles [28, 29, 30, 31, 32,
33, 34]. The large-scale elasticity properties exhibited in both amorphous and
crystalline solids link these two typically distinct branches of condensed
matter physics [35, 30, 36]. Recent studies on near-crystalline materials have
revealed various characteristics similar to those found in fully amorphous
materials, including the presence of quasi-localized modes [37, 38, 39]. Such
near-crystalline materials therefore help in establishing a connection between
the well studied physics of crystals and that of amorphous solids by
introducing disorder gradually into athermal crystalline packings. On the one
hand, while crystals have ordered structures, amorphous solids are
characterized by random and inflexible structures that arise from the
competing interactions between constituent particles. Despite their distinct
local structures [40], crystalline and amorphous packings exhibit many similar
elastic properties [41, 18]. Given the long-range displacement correlations in
these systems [42, 43], one may reasonably question whether the microscopic
structure affects the large-scale elasticity properties [44]. Therefore a
crucial question to address is how global rigidity is manifested in distinct
networks and whether this can be detected in the local stress tensor
fluctuations [45, 46]. Although some properties of athermal ensembles can be
described by temperature-like variables [47, 48], and despite several attempts
at a unifying framework, there is still a lack of understanding of these
properties in non-isotropic materials and near-crystalline systems. It is
therefore of interest to generate ensembles that can be precisely
characterized theoretically. In this paper, we develop a microscopic theory
for stress correlations in near-crystalline systems arising from various types
of disorders such as particle size disorder or due to quenched forces.

Stress correlations are a key ingredient in understanding the physics of
disordered systems, and there have been several recent studies that establish
their importance in amorphous materials [49, 50]. Such correlations provide
valuable insights into the collective behavior of interacting particles and
are widely used in fields such as material science, fluid dynamics, and
geophysics [51, 19, 52, 53]. Understanding stress correlations can help us
predict the strength and stability of materials under various external
conditions [54, 55] as well as give insights into the rheology of particulate
packings, such as their ability to flow or resist deformation. Stress
correlations provide a deeper insight into the degree of rigidity or
floppiness within particle packings and how they react to external influences
like shear or compression [17, 41]. Recently, stress correlations have also
been studied in various types of systems such as glasses, granular packings,
and gels amongst others, and has become a question that has attracted
considerable interest [17, 18, 19, 20, 21, 22, 23]. There have been several
theoretical studies that use material isotropy and homogeneity in amorphous
materials to derive the large length scale anisotropic behavior of the stress
correlations [56, 57, 58]. Several numerical studies have also explored stress
correlations in computer-simulated disordered packings [22, 23].

The main results of this paper can be summarized as follows. We derive the
displacement fields due to the introduction of particle size disorder or due
to external quenched force in a crystalline system through a microscopic
disorder perturbation expansion. Using the linear order displacement and force
fields, we derive the components of the change in the local stress tensor on
each grain. At large lengthscales, the local stresses show anisotropic
$1/r^{d}$ radial behavior for particle size disorder and $1/r^{d-1}$ behavior
for external force quenching. We analyze the local pressure fluctuation which
shows similar radial behavior yet isotropic at large lengthscales. We also
measure the global bulk and shear modulus for such a near-crystalline system
which show excellent match with simulations for finite small disorder. We then
derive the configurational averaged correlations of the local stress
fluctuations which are verified through numerical simulations in two different
models in both two and three dimensions. We show that the stress correlations
in disordered crystals show different behavior to that of an isotropic
amorphous material at a high packing fraction or high pressure limit.

The outline of the paper is as follows. In Section 2 we introduce the
microscopic models, while in Section 3, we present the corresponding
preparation protocols. In Section 4 we employ the microscopic approach to
derive macroscopic properties, specifically the bulk and shear moduli of near-
crystalline systems. In Section 5 we present a detailed derivation of the
displacement fields due to the introduction of microscopic disorder in the
crystalline packing. In section 6 we derive the change in the local stress
tensor components and their correlations through this microscopic approach and
compare these predictions against direct numerical simulations. Additionally,
we draw parallels between these results and those obtained in the recently
developed VCTG framework for amorphous systems. In Section 6.5, we extend our
theory to three dimensional fcc arrangement of particles. In Section 7, we
derive the displacement and force fields due to point forces. In Section 8, we
derive the distributions of change in local stresses from the linear relations
between the stress and microscopic disorder in a near-crystalline system.
Finally, we conclude and provide directions for future investigations in
section 9.

Figure 1: Schematic of particle arrangement in a hexagonal close packing for
two different models considered in this paper. (a) Soft repulsive interactions
with a Harmonic potential, where the neighboring particles overlap with the
central particle. (b) Lennard-Jones (LJ) interaction with a cutoff. Here, the
outermost circle corresponds to the interaction range with respect to the
central particle (red).

##  2 Models

We study two well-known canonical glass-forming model potentials that can be
used to create amorphous, as well as near-crystalline structures: short-ranged
Harmonic interactions, and an attractive Lennard Jones interaction with a
cutoff.

###  2.1 Short-ranged repulsive Harmonic interaction

For the case of the short-ranged Harmonic model, we examine systems consisting
of frictionless soft disks in two dimensions and spheres in three dimensions
with different levels of overcompression. These particles interact with one
another through a one-sided pairwise potential [59, 14, 31] which takes the
following form:

$\displaystyle V_{a_{ij}}\left(\vec{r}_{ij}\right)$
$\displaystyle=\frac{K}{\alpha}\left(1-\frac{\left|\vec{r}_{ij}\right|}{a_{ij}}\right)^{\alpha}\Theta\left(1-\frac{\left|\vec{r}_{ij}\right|}{a_{ij}}\right),$
(1)

where, $\vec{r}_{ij}$ represents the displacement vector between particles $i$
and $j$, situated at positions $\vec{r}_{i}$ and $\vec{r}_{j}$, respectively.
Here $a_{ij}$ are called the quenched interaction lengths that are defined as
the sum of individual radii, denoted as $a_{ij}=a_{i}+a_{j}$. In this study,
we select $\alpha=2$ to establish a harmonic pairwise potential between the
particles. The length parameters $a_{ij}$ are then set as follows [60, 42, 43,
32, 61, 62]

$a_{ij}=2a_{0}+\eta a_{0}(\zeta_{i}+\zeta_{j}),$ (2)

where $a_{0}$ is the radius of each particle in the crystalline state. The
variables $\zeta_{i}$ are independent and identically distributed random
numbers drawn from a uniform distribution ranging between $-1/2$ and $1/2$,
that are individually assigned to each particle within the system. The
parameter $\eta$ (polydispersity) controls the magnitude of the disorder.

###  2.2 Attractive Lennard-Jones interaction with cut-off

We investigate particles interacting via long-ranged power-law potentials,
which are smoothened up to the second order at a specified cutoff interaction
length ($r_{ij}^{c}$), set at $2.5a_{ij}$ [29, 39]. This cut-off is set to
speed up the numerical simulations and for most purposes, a cut-off greater
than $1.5$ yields similar mechanical properties [63]. The smoothened LJ
potential for a cut-off interaction length $2.5a_{ij}$ can be represented as

$V_{a_{ij}}\left(\vec{r}_{ij}\right)=4K\left[\left(\frac{a_{ij}}{\left|\vec{r}_{ij}\right|}\right)^{12}-\left(\frac{a_{ij}}{\left|\vec{r}_{ij}\right|}\right)^{6}+\sum_{l=0}^{2}c_{2l}\left(\frac{\left|\vec{r}_{ij}\right|}{a_{ij}}\right)^{2l}\right]\Theta\left(2.5-\frac{\left|\vec{r}_{ij}\right|}{a_{ij}}\right),$
(3)

The disorder is introduced through the length parameters $a_{ij}$ and can be
represented as [29, 39],

$a_{ij}=\left\\{\begin{array}[]{cc}\lambda_{\mathrm{SS}}&\text{ both
}i,\text{and }j\text{ are unlabeled },\\\
\eta\left(\lambda_{\mathrm{SL}}-\lambda_{\mathrm{SS}}\right)+\lambda_{\mathrm{SS}}&\text{
either }i\text{ or }j\text{ is labeled },\\\
\eta\left(\lambda_{\mathrm{LL}}-\lambda_{\mathrm{SS}}\right)+\lambda_{\mathrm{SS}}&\text{
both }i\text{ and }j\text{ are labeled }.\end{array}\right.$ (4)

This model corresponds to a bidisperse system where $\eta$ controls the
strength of the disorder. For theoretical simplicity, instead of using length
parameter $a_{ij}$, we define an onsite parameter $t_{i}$. The variable
$t_{i}$ takes a value of either $1$ or $0$, depending on whether the particle
at position $\vec{r}_{i}$ is labeled or not. Using $t_{i}$ we can express
$a_{ij}$ as follows

$\displaystyle a_{ij}=\lambda_{SS}+$
$\displaystyle\eta\left[(t_{i}+t_{j})\left(\lambda_{SL}-\lambda_{SS}\right)+t_{i}t_{j}\left(\lambda_{LL}+\lambda_{SS}-2\lambda_{SL}\right)\right].$
(5)

##  3 Numerical Simulations

To test our theoretical predictions, we conduct simulations of an athermal,
over-compressed triangular lattice (hcp) in two dimensions (2D) and a face-
centered cubic (fcc) lattice in three dimensions (3D) with soft, frictionless
particles with varying levels of particle size disorder. We employ periodic
boundary conditions to account for boundary effects. Our focus is on states in
which every particle achieves force balance, i.e., configurations
corresponding to energy minima. To achieve this, we utilize the Fast Inertial
Relaxation Engine (FIRE) algorithm as described in Ref. [64] to minimize the
energy of the system.

Here the inter-particle separation is kept fixed at $R_{0}$ in the initial
crystalline state. We initially consider a rectangular (2D) grid spacings of
$R_{0}/2$ and $\sqrt{3}R_{0}/2$ along $x$ and $y$ directions respectively and
cubic (3D) lattice with a grid spacing of $R_{0}$. To create a triangular/fcc
arrangement, particles are placed on alternate grid points that satisfy the
respective conditions, i.e., $n_{x}+n_{y}=2n$ for a triangular lattice and
$n_{x}+n_{y}+n_{z}=2n$ for an fcc lattice, where $n$ is an integer. This
technique is an extension of the one used in reference [65] for generating a
hexagonal close packing in two dimensions. A square/cubic lattice has
$4L^{2}/8L^{3}$ grid points in total (since there are $2L$ grid points along
each coordinate axis), of which only half are occupied by particles, yielding
the total number of particles, $N=2L^{2}$ in 2D and $4L^{3}$ in 3D.

###  3.1 Harmonic Model

We first consider the Harmonic model, where particles experience exclusively
repulsive interactions and each grain interacts solely with its closest
neighbors. In a two dimensional triangular lattice, this corresponds to six
neighboring grains, while in a three dimensional fcc lattice, it corresponds
to twelve neighboring grains. The degree of compression in the lattice is
indicated by the packing fraction, which is set to $\phi=0.92/0.96/0.98$ in 2D
and $\phi=0.80$ in 3D. This is in comparison to the marginally jammed
triangular/fcc lattice with a packing fraction of $\phi_{c}\approx 0.9069$ in
2D and $\phi_{c}\approx 0.74$ in 3D. The interparticle spacing ($R_{0}$) is
defined by the initial particle radius, which is set to $a_{0}=0.5$, and the
packing fraction, calculated using the formula
$R_{0}=2a_{0}(\phi_{c}/\phi)^{1/d}$, in the absence of any disorder. Here we
have chosen bond stiffness/interaction strength as $K=0.5$. The numerical
results presented in this study are averaged over $200$ different realizations
of disordered states. These simulations were performed for system sizes of
$N=6400$ and $10000$ particles in two dimensions, and $N=32000$ and $250000$
in three dimensions, with different strengths of particle size disorder
($\eta$).

###  3.2 LJ Model

For the case of Lennard-Jones (LJ) interactions, every particle interacts with
its neighboring particles located within a cutoff radius defined as
$\left|\vec{r}_{ij}\right|/a_{ij}\leq 2.5$. In the absence of perturbation
($\eta=0$), this condition implies that there are a total of 18 neighboring
particles within the interaction range corresponding to each grain in the
triangular lattice. Here we have chosen $K=0.5$, with
$\lambda_{SS}=1,\lambda_{LL}=1.4\lambda_{SS}$ and
$\lambda_{SL}=1.2\lambda_{SS}$ for the interaction potential. Similar to the
Harmonic model, we have performed simulations for systems of size $N=6400$ in
$2D$. Since the volume associated with each particle is not defined, we define
a number density ($\rho_{N}$) as our initial parameter instead of a packing
fraction. Each number density, $\rho_{N}$ corresponds to different values of
initial pressure ($P$) in the system. Here the results are presented for
$P=0/0.27/4.24$. To achieve the initial pressure we use a Berendsen barostat
[66] which is implemented in the FIRE algorithm during the energy
minimisation.

##  4 Elastic properties of near-crystalline packings

In this section, we derive the results that can be used to compute the
relevant macroscopic elastic properties in a near-crystalline granular packing
composed of soft particles. The fundamental property we are interested in is
the global pressure, denoted by the symbol $P$ which can be written as,

$P=d^{-1}\sum_{\mu}\Sigma_{\mu\mu}=(dV)^{-1}\sum_{\mu,\langle
ij\rangle}r^{\mu}_{ij}f_{ij}^{\mu},$ (6)

where $\Sigma_{\mu\mu}$ are the diagonal components of the global stress
tensor where $d$ and $V$ represent the spatial dimension and total volume of
the system respectively. Here $r^{\mu}_{ij}$ and $f_{ij}^{\mu}$ are the $\mu-
th$ component of the relative displacement and force between the particles $i$
and $j$. The configurational averaged total pressure can be defined as

$\displaystyle\left\langle
P\right\rangle=d^{-1}\sum_{\mu}\langle\Sigma_{\mu\mu}\rangle=$ $\displaystyle
d^{-1}\sum_{\mu}\langle\Sigma_{\mu\mu}^{(0)}\rangle+\langle\delta\Sigma_{\mu\mu}^{(0)}\rangle,$
(7)

where, $\Sigma_{\alpha\beta}^{0}=V^{-1}\sum_{\langle
ij\rangle}r^{\alpha(0)}_{ij}f_{ij}^{\beta(0)}$ are the components of the
global stress tensor of the crystalline system without the disorder. For a
small magnitude of the disorder strength $\eta$, the average change in the
global pressure is zero (i.e., $\langle\delta\Sigma_{\mu\mu}^{(0)}\rangle=0$),
which we will show in the following section. For any regular arrangement of
particles in $d$-dimension, we can write the relative distance between the
neighboring particles and their corresponding forces as

$\vec{r}_{ij}^{(0)}=R_{0}\hat{r}_{ij}^{0},\hskip
14.22636pt\vec{f}_{ij}^{(0)}=\frac{K}{a_{0}}\left(1-\frac{R_{0}}{2a_{0}}\right)\hat{r}_{ij}^{0},$
(8)

where $R_{0}$ is the relative distance between two neighboring particles in
the crystalline arrangement. Inserting these values, one can obtain the final
form of the averaged net pressure as,

$\displaystyle\left\langle P\right\rangle=$
$\displaystyle\frac{z_{0}NKR_{0}}{2dVa_{0}}\left(1-\frac{R_{0}}{2a_{0}}\right)=\frac{z_{0}K\rho_{N}^{0}}{d}\left(\frac{\phi_{0}}{\phi}\right)^{-1+1/d}\left(1-\left(\frac{\phi_{0}}{\phi}\right)^{1/d}\right).$
(9)

Here $z_{0}$ is the coordination number of each grain and $\rho_{N}^{0}$ is
the number density in the marginally jammed crystal. We can also represent the
number density of an over-compressed crystal as,
$\rho_{N}=N/V=\rho_{N}^{0}(\phi/\phi_{0})^{1/d}$. In the two dimensional
triangular lattice, $z_{0}=6$ and $\rho_{N}^{0}=1/2\sqrt{3}a_{0}^{2}$. The
system is subjected to an isotropic strain, where the box dimensions along all
the Cartesian directions increase by a factor of $(1+\epsilon)$. Consequently,
the packing fraction $\phi$ changes to $\phi/(1+\epsilon)^{d}$. The pressure
$P$ of this isotropically strained system can be expressed as

$\displaystyle\left\langle P^{\prime}\right\rangle=$
$\displaystyle\frac{z_{0}NkR_{0}}{2dVa_{0}}\left(1-\frac{R_{0}}{2a_{0}}\right)=\frac{z_{0}k\rho_{N}^{0}}{d(1+\epsilon)^{d-1}}\left(\frac{\phi_{0}}{\phi}\right)^{-1+1/d}\left(1-(1+\epsilon)\left(\frac{\phi_{0}}{\phi}\right)^{1/d}\right).$
(10)

Figure 2:  $(a)$ Variation of configurational averaged bulk ($\langle
B\rangle$) and shear modulus ($\langle G\rangle$) with polydispersity for a
near-crystalline packing of soft particles in a system of size $N=256$ with
packing fraction, $\phi=\phi_{0}+0.02$ in two dimensions. From the above plot,
we can notice that the elastic properties of these disordered crystals are
similar to that of a perfect crystal for a disorder strength, $\eta\leq 0.03$
which also changes with the initial packing fraction. $(b)$ Variation of
$\langle B\rangle$ and $\langle G\rangle$ with initial packing fraction for a
fixed polydispersity, $\eta=0.005$.

Here, $\epsilon\sim\delta V/2V$ is proportional to the volumetric strain
applied to the system. The average bulk modulus ($\langle B\rangle$) for these
near-crystalline systems can be obtained by finding the ratio of change in
bulk pressure to the volumetric strain,

$\langle B\rangle=\left|\frac{\delta\left\langle P\right\rangle}{\delta
V/V}\right|=\lim_{\epsilon\to 0}\left|\frac{\left\langle
P^{\prime}\right\rangle-\left\langle
P\right\rangle}{2\epsilon}\right|=\frac{z_{0}k\rho_{N}^{0}}{2d}\left(\frac{\phi_{0}}{\phi}\right)^{-1+1/d}\left((d-2)\left(\frac{\phi_{0}}{\phi}\right)^{1/d}-(d-1)\right).$
(11)

In two dimensions, the average bulk modulus for a small magnitude of disorder
in a crystalline packing can be written as $\langle
B\rangle=\frac{\sqrt{3}K}{4a_{0}^{2}}(\phi/\phi_{0})^{1/2}$, which we have
also verified numerically for $\phi=\phi_{0}+0.02$ as shown in Fig. 2. Next,
we consider a near-crystalline system that is sheared along
$\alpha\beta$-plane with a shear amplitude of $\gamma$. Only the affine part
of the displacements is considered for small shear amplitude, leading to the
following expressions for the relative displacement and force components

$\displaystyle r_{ij}^{\alpha}$
$\displaystyle=R_{0}(\cos{\theta_{ij}^{0}}+\gamma\sin{\theta_{ij}^{0}}),\hskip
14.22636ptr_{ij}^{\beta}=R_{0}\sin{\theta_{ij}^{0}},\hskip
14.22636ptR_{ij}=\left(\sum_{\alpha=1}^{d}(r_{ij}^{\alpha})^{2}\right)^{1/2}.$
(12)

where $\theta_{ij}$ is the angle made by the projection of $r_{ij}$ on the
$\alpha\beta$-plane to the $\alpha$-axis. Using the above relations, we can
write the $\alpha\beta$-component of the stress tensor as

$\displaystyle\langle\Sigma_{\alpha\beta}\rangle$
$\displaystyle=\frac{1}{2V}\sum_{i}\left(\sum_{j}r^{\alpha}_{ij}f^{\beta}_{ij}\right)=k\rho_{N}^{0}\left(\frac{\phi_{0}}{\phi}\right)^{-1+1/d}\left(\sum_{j=1}^{z}\left(1-\frac{R_{ij}}{2a_{0}}\right)(\cos{\theta_{ij}^{0}}\sin{\theta_{ij}^{0}}+\gamma\sin^{2}{\theta_{ij}^{0}})\right).$
(13)

For a two dimensional triangular lattice, the above equation simplifies to

$\displaystyle\langle\Sigma_{xy}\rangle$
$\displaystyle=\frac{3NkR_{0}}{2Va_{0}}\left[-\frac{R_{0}}{2a_{0}}\gamma+\frac{1}{\sqrt{3}}\left(\frac{\gamma-1/\sqrt{3}}{\sqrt{1+(\gamma-1/\sqrt{3})^{2}}}+\frac{\gamma+1/\sqrt{3}}{\sqrt{1+(\gamma+1/\sqrt{3})^{2}}}\right)\right].$
(14)

Considering a small shear amplitude, $\gamma$, we can do a linear
approximation of the shear stress and take the ratio of change in the shear-
stress to shear-strain to obtain the shear-modulus

$G=\lim_{\gamma\to
0}\frac{\langle\Sigma_{xy}\rangle}{\gamma}\sim\frac{3NkR_{0}}{2Va_{0}}\left(\frac{R_{0}}{2a_{0}}-\frac{3}{4}\right)=\frac{\sqrt{3}k}{2a_{0}^{2}}\left(1-\frac{3}{4}\left(\frac{\phi}{\phi_{0}}\right)^{1/2}\right).$
(15)

Both the expressions for bulk and shear modulus for various packing fraction
and particle size disorder are validated through numerical simulations in
near-crystalline packings of soft particles as shown in Fig. 2. Any local
fluctuations of pressure and shear stresses giving rise to local fluctuations
in bulk and shear modulus are discussed in the later section. Given the planar
bulk and shear moduli, one can obtain the expressions for planar Young’s
modulus and Poisson’s ratio as

$\displaystyle E$
$\displaystyle=\frac{4BG}{B+G}=\frac{\sqrt{3}k}{a_{0}^{2}}\left(\frac{\phi}{\phi_{0}}\right)^{1/2}\left(\frac{4-3\left(\phi/\phi_{0}\right)^{1/2}}{4-\left(\phi/\phi_{0}\right)^{1/2}}\right),$
(16) $\displaystyle\nu$
$\displaystyle=\frac{B-G}{B+G}=\frac{5\left(\phi/\phi_{0}\right)^{1/2}-4}{4-\left(\phi/\phi_{0}\right)^{1/2}}.$

Similar techniques can be used in a system of particles interacting via long-
ranged Lennard-Jones interactions in a near-crystalline packing with average
particle separation between the nearest neighbor being $R_{0}$. The bulk and
shear modulus for such a system can be computed as

$\displaystyle B$
$\displaystyle=\frac{24\sqrt{3}k}{R_{0}^{8}}\left(a-\frac{4b}{R_{0}^{6}}\right),$
(17) $\displaystyle G$
$\displaystyle=\frac{4\sqrt{3}k}{R_{0}^{8}}\left(a-\frac{5b}{R_{0}^{6}}\right),$

where,

$a=\sum_{i=1}^{N}\frac{n_{i}}{m_{i}^{6}}\text{, and
}b=\sum_{i=1}^{N}\frac{n_{i}}{m_{i}^{12}}.$ (18)

Here $n_{i}$ is the number of particles on the $i^{th}$ spherical cell and
$m_{i}$ is the ratio of the distance of the $i^{th}$ cell from the central
particle to $R_{0}$, i.e., $m_{i}=R_{i}/R_{0}$. In the $N\to\infty$ limit, we
arrive at $a=6.37588$ and $b=6.00981$.

##  5 Displacement fields induced by microscopic disorder

In the previous section, we presented theoretical results related to the
average stress tensor components in a nearly crystalline arrangement of soft
particles, where we made the assumption that the average changes in local
stress are negligible. The local stress components are proportional to the
square of the interparticle distances of all neighboring particles within the
cut-off distance. Consequently, in order to formulate the expressions for
local stress, it is necessary to derive the displacements of individual
particles in a disordered configuration. Below, we derive the displacement and
force fields resulting from the introduction of disorder into a crystalline
network.

In both Harmonic and LJ model, we begin with a crystalline packing of
monodisperse particles in a fixed volume. In the short-ranged repulsive model
we start with a finite overcompression whereas in the attractive LJ
interaction model, the initial volume is fixed such that the initial pressure
is set to $P=0$. We then introduce disorder in the effective particle sizes
i.e., $a_{ij}$ as given in Eqs. (2), (5). As a response to this disorder, the
particles are displaced from their crystalline positions to maintain force
balance, as

$r^{\mu}_{i}=r^{\mu(0)}_{i}+\delta r^{\mu}_{i}.$ (19)

Given that particle $j$ is one of the neighboring particles of particle $i$ in
the initial crystalline lattice, the relative displacement between their
positions can be expressed using the basis lattice vectors of the crystalline
lattice as
$\vec{r}_{ij}^{(0)}=\vec{r}_{j}^{(0)}-\vec{r}_{i}^{(0)}=\vec{\Delta}_{j}$. The
discrete Fourier transform of the change in the relative displacement $\delta
r_{ij}^{\mu}$ can be expressed as:

$\displaystyle\mathcal{F}\left[\delta
r^{\mu}_{ij}\right]=\sum_{i}e^{i\vec{r}_{i}^{(0)}.\vec{k}}\delta
r^{\mu}_{ij}=$
$\displaystyle\sum_{i}e^{i\vec{r}_{i}^{(0)}.\vec{k}}\left(\delta
r_{j}^{\mu}-\delta
r_{i}^{\mu}\right)=\left[e^{-i\vec{\Delta}_{j}.\vec{k}}-1\right]\delta\tilde{r}^{\mu}(\vec{k}),$
(20)

where, $\delta\tilde{r}^{\mu}(\vec{k})=\mathcal{F}\left[\delta
r^{\mu}(\vec{r})\right]$, corresponds to the discrete Fourier transform of
particle displacements from their crystalline positions. As a response to the
disorder as well as the displacements in the particle positions, the forces
$\vec{f}_{ij}$ acting between adjacent particles $i$ and $j$ also change. This
variation can be expressed as a perturbation relative to the forces between
these particles in the initial crystalline state, represented as

$f_{ij}^{\mu}=f_{ij}^{\mu(0)}+\delta f_{ij}^{\mu}.$ (21)

Every individual component of the excess force $\delta f_{ij}^{\mu}$ acting
between particles $i$ and $j$ can be Taylor expanded up to first order about
its value in the crystalline ground state, in terms of $\delta r^{\mu}_{ij}$
($=\delta r^{\mu}_{j}-\delta r^{\mu}_{i}$) and $\delta a_{ij}$ ($=\delta
a_{i}+\delta a_{j}$) as

$\delta f_{ij}^{\mu}=\sum_{\nu}C_{ij}^{\mu\nu}\delta r^{\nu}_{ij}+C_{ij}^{\mu
a}\delta a_{ij},$ (22)

where the first-order Taylor coefficients $C_{ij}^{\mu\nu},C_{ij}^{\mu a}$
depend only on the form of the potential between the interacting particles and
the initial crystalline arrangement. These coefficients can be represented as,
$\left.C_{ij}^{\mu\nu}=\left(\partial f_{ij}^{\mu}/\partial
r_{ij}^{\nu}\right)\right|_{\\{\vec{r}_{ij}^{(0)},a_{ij}^{(0)}\\}}$ and
$\left.C_{ij}^{\mu a}=\left(\partial f_{ij}^{\mu}/\partial
a_{ij}\right)\right|_{\\{\vec{r}_{ij}^{(0)},a_{ij}^{(0)}\\}}$. For energy-
minimized configurations, the force balance condition dictates that the net
force acting on each particle $i$ is zero. This means that for all interacting
neighbors $j$, the sum of the force deviations, denoted as
$\delta\vec{f}_{ij}$, is equal to zero, expressed as
$\sum_{j}\delta\vec{f}_{ij}=0$. By applying this condition in the linear order
expression for forces as given in Eq. (22), we obtain $d$-equations, where $d$
represents the dimension of the system, for each particle $i$. These equations
are given as follows:

$\displaystyle\sum_{j}\sum_{\nu}C_{ij}^{\mu\nu}\delta r^{\nu}_{ij}=$
$\displaystyle-\sum_{j}C_{ij}^{\mu a}\delta a_{ij},$ (23)
$\displaystyle\mathcal{P}_{1}\ket{\delta r}=$
$\displaystyle\mathcal{P}_{2}\ket{\delta a}.$

Here, we have $Nd$ such equations corresponding to $Nd$-variables
(displacement components). Since the system is translationally invariant,
performing a discrete Fourier transform on the equation can convert the
$Nd$-equations of $Nd$-variables into $d$-equations of $d$ variables. This
simplification reduces the complexity of the problem and diagonalizes the
large matrices $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$. So the Fourier
transform of the Eq. (23) leads to

$\displaystyle\sum_{\nu}A^{\mu\nu}(\vec{k})\delta\tilde{r}^{\nu}(\vec{k})=B^{\mu}(\vec{k}),$
(24)

where,

$\displaystyle A^{\mu\nu}(\vec{k})=$
$\displaystyle\sum_{j}\left(1-e^{-i\vec{k}.\vec{\Delta}_{j}}\right)C_{ij}^{\mu\nu},\hskip
14.22636ptB^{\mu}(\vec{k})=-\mathcal{F}\left[\sum_{j}C_{ij}^{\mu a}\delta
a_{ij}\right].$ (25)

In $d$-dimensions, $A^{\mu\nu}$ would be a $d\times d$ symmetric matrix. Here
the expression for $A^{\mu\nu}$ and $B^{\mu}$ have the same form for both
short and long-range models where the only difference lies in the number of
interacting neighbor particles. Now we can obtain the displacement fields in
Fourier space by inverting Eq. (24),

$\delta\tilde{r}^{\mu}(\vec{k})=\sum_{\nu}(A^{-1})^{\mu\nu}(\vec{k})B^{\nu}(\vec{k}).$
(26)

Since $\delta\tilde{r}$ is expressed as the product of $A^{-1}$ and $B$, its
inverse Fourier transform can be written as a convolution resulting in the
displacement fields in real space, as shown below:

$\delta
r^{\mu}(\vec{r})=\mathcal{F}^{-1}\left[\delta\tilde{r}^{\mu}(\vec{k})\right]=\frac{1}{N}\sum_{\vec{k}}\exp(-i\vec{k}.\vec{r})\delta\tilde{r}^{\mu}(\vec{k}).$
(27)

Here Eqs. (27) and (26) correspond to the displacement fields and their
Fourier transform in the presence of particle size disorder. The exact
expressions of these displacements are model dependent, and we discuss the two
different scenarios in detail below.

###  5.1 Short-ranged repulsive harmonic interaction

The displacement fields for the Harmonic model has been studied extensively
[60, 32, 43, 42] where the disorder is introduced in the particle radius,
$\delta a_{ij}=\eta a_{0}(\zeta_{i}+\zeta_{j})$. Next, putting this in Eq.
(25) we get the expression for $B^{\mu}$,

$\displaystyle B^{\mu}(\vec{k})=$
$\displaystyle-D^{\mu}(\vec{k})\delta\tilde{a}(\vec{k})=-\eta
a_{0}D^{\mu}(\vec{k})\tilde{\zeta}(\vec{k}),$ (28)
$\displaystyle\text{where,}\hskip 14.22636ptD^{\mu}(\vec{k})=$
$\displaystyle\sum_{j}\left(1+e^{-i\vec{k}.\vec{\Delta}_{j}}\right)C^{\mu
a}_{ij}.$

Here $\delta\tilde{a}(\vec{k})=\mathcal{F}\left[\delta a(\vec{r})\right]=\eta
a_{0}\tilde{\zeta}(\vec{k})$, is the Fourier transform of $\delta a_{i}$. We
have defined $|\vec{\Delta}_{j}|=R_{0}$ as the magnitude of the relative
distance between grains in an over-compressed crystalline system. Now putting
the expression of $B^{\nu}(\vec{k})$ in the expression of
$\delta\tilde{r}^{\mu}(\vec{k})$ in Eq. (26), we obtain

$\delta\tilde{r}^{\mu}(\vec{k})=\underbrace{\left[-\sum_{\nu}(A^{-1})^{\mu\nu}(\vec{k})D^{\nu}(\vec{k})\right]}_{\tilde{G}^{\mu}(\vec{k})}\delta\tilde{a}(\vec{k})=\tilde{G}^{\mu}(\vec{k})\delta\tilde{a}(\vec{k}).$
(29)

We can get the displacement fields by taking an inverse discrete Fourier
transform as given in Eq. (27) as,

$\displaystyle\delta r^{\mu}(\vec{r})=$
$\displaystyle\sum_{\vec{r}^{\prime}}G^{\mu}(\vec{r}-\vec{r}^{\prime})\delta
a(\vec{r}^{\prime}),$ (30) $\displaystyle\text{where,\hskip
14.22636pt}G^{\mu}(\vec{r})=$
$\displaystyle\mathcal{F}^{-1}\left[\tilde{G}^{\mu}(\vec{k})\right].$

###  5.2 Attractive Lennard-Jones interaction with cut-off

The above-mentioned formulation for displacement fields can also be extended
to any sort of interaction where every grain can interact with 6 or more
neighbors depending on the interaction cut-off. For example, this cutoff is
$|\vec{r}_{ij}|/a_{ij}\leq 2.5$ for the LJ model, where the microscopic
disorders are incorporated into the bond distances as

$\displaystyle\delta a_{ij}=$
$\displaystyle\eta\left[(t_{i}+t_{j})\left(\lambda_{SL}-\lambda_{SS}\right)+t_{i}t_{j}\left(\lambda_{LL}+\lambda_{SS}-2\lambda_{SL}\right)\right].$
(31)

Therefore, $B^{\mu}(\vec{k})$ has the form

$\displaystyle B^{\mu}(\vec{k})=$
$\displaystyle\eta\left[\left(\lambda_{SS}-\lambda_{SL}\right)\tilde{D}^{\mu}(\vec{k})\tilde{t}(\vec{k})-\frac{\left(\lambda_{LL}+\lambda_{SS}-2\lambda_{SL}\right)}{N}\sum_{\vec{k}^{\prime}}\left[\tilde{t}(\vec{k}^{\prime})\tilde{t}(\vec{k}-\vec{k}^{\prime})\tilde{D}^{\mu}(\vec{k}-\vec{k}^{\prime})\right]\right].$
(32)

In the above expression, $\tilde{D}^{\mu}(\vec{k})$ has the same expression as
given in Eq. (28) with $j$ going from $1-18$ for all the interacting neighbors
within the range $|\vec{r}_{ij}|/a_{ij}\leq 2.5$. The only difference in the
expression of $B^{\mu}$ in LJ model to that of the Harmonic model is that the
magnitude of $\vec{\Delta}_{j}$s are not constant for all the interacting
neighbors. In this study, we have chosen,
$\lambda_{SL}=(\lambda_{SS}+\lambda_{LL})/2$ for numerical simulations, which
simplifies the problem by removing the nonlinear term in the above expression
for $B^{\mu}(\vec{k})$. In this approximation we can write, $\delta
a_{i}=\eta(\lambda_{SL}-\lambda_{SS})t_{i}$. Now the expressions for
$B^{\mu}(\vec{k})$ and $\delta\tilde{r}^{\mu}(\vec{k})$ in Eq. (26) can be
written as,

$\displaystyle B^{\mu}(\vec{k})=$
$\displaystyle-D^{\mu}(\vec{k})\delta\tilde{a}(\vec{k})=-\eta(\lambda_{SL}-\lambda_{SS})D^{\mu}(\vec{k})\tilde{t}(\vec{k}),$
(33) $\displaystyle\delta\tilde{r}^{\mu}(\vec{k})=$
$\displaystyle\underbrace{\left(-\sum_{\nu}(A^{-1})^{\mu\nu}(\vec{k})D^{\nu}(\vec{k})\right)}_{\tilde{G}^{\mu}(\vec{k})}\delta\tilde{a}(\vec{k}).$

where,
$\delta\tilde{a}(\vec{k})=\eta(\lambda_{SL}-\lambda_{SS})\tilde{t}(\vec{k})$.

##  6 Stress correlations induced by microscopic disorder

In this section, we focus on the fluctuations and correlations of the local
stress tensor components. Using the linear order displacement fields derived
in the earlier section, we can compute the stress correlations using a similar
perturbation expansion for the minimally polydisperse system. For any athermal
jammed system, the components of the global stress tensor can be represented
as

$\displaystyle\Sigma_{\alpha\beta}$ $\displaystyle=V^{-1}\sum_{\langle
ij\rangle}r_{ij}^{\alpha}f_{ij}^{\beta},$ (34)

where, $r_{ij}^{\alpha}=r_{ij}^{\alpha(0)}+\delta r_{ij}^{\alpha}$ and
$f_{ij}^{\beta}=f_{ij}^{\beta(0)}+\delta f_{ij}^{\beta}$. Here
$r_{ij}^{\mu(0)}$ refer to the $\mu$-th component of the interparticle
distance whereas $f_{ij}^{\mu(0)}$ denote the $\mu$-th component of the force
acting between particle $i$ and $j$ in the crystalline lattice without any
microscopic disorder. In a system with small particle size polydispersity
($\eta$), both the deviation of the particle positions ($|\delta\vec{r}|$)
from the crystalline positions and the change in inter-particle forces
($|\delta\vec{f}|$) are in the order of $\delta a\sim\eta$. For small values
of $\eta$, we can neglect higher-order terms in the expression for global
stress which leads to

$\displaystyle\Sigma_{\alpha\beta}\sim V^{-1}\sum_{\langle
ij\rangle}\left(r^{\alpha(0)}_{ij}f_{ij}^{\beta(0)}+r^{\alpha(0)}_{ij}\delta
f_{ij}^{\beta}+\delta r^{\alpha}_{ij}f_{ij}^{\beta(0)}\right).$ (35)

Therefore the incremental change in the global stress (i.e
$\delta\Sigma=\Sigma-\Sigma^{(0)}$) due to the introduction of disorder in the
crystalline system can be written as

$\displaystyle\delta\Sigma_{\alpha\beta}=V^{-1}\sum_{i}\underbrace{\left[\sum_{j}(r^{\alpha(0)}_{ij}\delta
f_{ij}^{\beta}+\delta
r^{\alpha}_{ij}f_{ij}^{\beta(0)})\right]}_{\delta\sigma_{\alpha\beta}(\vec{r}_{i})}.$
(36)

The above expression represents the net change in global stress as a linear
combination of $\delta\sigma(\vec{r}_{i}^{0})$, which we define as the change
in local stress at the lattice position $\vec{r}_{i}^{(0)}$, and can be
expressed as follows:

$\displaystyle\delta\sigma_{\alpha\beta}(\vec{r}_{i}^{0})=$
$\displaystyle\sum_{j}(r^{\alpha(0)}_{ij}\delta f_{ij}^{\beta}+\delta
r^{\alpha}_{ij}f_{ij}^{\beta(0)})=\sum_{j}\left[\Delta_{j}^{\alpha}\sum_{\nu}C_{ij}^{\beta\nu}\delta
r_{ij}^{\nu}+\Delta_{j}^{\alpha}C_{ij}^{\beta a}\delta a_{ij}+\delta
r^{\alpha}_{ij}f_{ij}^{\beta(0)}\right].$ (37)

In the above expression, $\delta\sigma_{\alpha\beta}(\vec{r}_{i})$ is a linear
summation of particle displacements and change in radii with coefficients. As
we have demonstrated earlier, the Fourier transform of $\delta
r^{\alpha}_{ij}$ and $\delta a_{ij}$ have simple relationships due to the
translational invariance of the system, as shown in Eqs. (29) and (33).
Therefore, we can simplify the problem by performing a discrete Fourier
transform of Eq. (37) which leads to

$\displaystyle\delta\tilde{\sigma}_{\alpha\beta}(\vec{k})=\mathcal{F}\left[\delta\sigma_{\alpha\beta}(\vec{r}_{i}^{0})\right]$
$\displaystyle=\sum_{i}e^{i\vec{r}_{i}^{0}.\vec{k}}\delta\sigma_{\alpha\beta}(\vec{r}_{i}^{0})=\sum_{j}\left[\Delta_{j}^{\alpha}\sum_{\nu}C_{j}^{\beta\nu}\left[-1+F_{j}(\vec{k})\right]\delta\tilde{r}^{\nu}(\vec{k})\right.$
(38) $\displaystyle+\left.\Delta_{j}^{\alpha}C_{j}^{\beta
a}\left[1+F_{j}(\vec{k})\right]\delta\tilde{a}(\vec{k})+\left[-1+F_{j}(\vec{k})\right]\delta\tilde{r}^{\alpha}(\vec{k})f_{j}^{\beta(0)}\right].$

where, $F_{j}(\vec{k})=\exp\\{-i\vec{\Delta}_{j}.\vec{k}\\}$. We can further
simplify the above expression by replacing $\delta\tilde{r}^{\nu}(\vec{k})$ by
$\tilde{G}^{\nu}(\vec{k})\delta\tilde{a}(\vec{k})$ as given in Eqs. (29) and
(33), to arrive at

$\displaystyle\delta\tilde{\sigma}_{\alpha\beta}(\vec{k})=$ $\displaystyle
S_{\alpha\beta}(\vec{k})\delta\tilde{a}(\vec{k})\hskip
11.38092pt\text{where,}$ (39) $\displaystyle S_{\alpha\beta}(\vec{k})=$
$\displaystyle\sum_{j}\left[\left[1+F_{j}(\vec{k})\right]C_{j}^{\beta
a}\Delta^{\alpha(0)}_{j}+\left[-1+F_{j}(\vec{k})\right]\left(\sum_{\nu}\Delta^{\alpha(0)}_{j}C_{j}^{\beta\nu}\tilde{G}^{\nu}(\vec{k})+f_{j}^{\beta(0)}\tilde{G}^{\alpha}(\vec{k})\right)\right].$

The sum over $j$ pertains to all neighboring particles, i.e., all the
particles that interact with the central particle in the crystalline state
without the disorder. Here, $S_{\alpha\beta}$ represents the Fourier transform
of the Green’s function for the change in local stress components. Next, we
can obtain the change in local stresses in real space as a convolution by
performing an inverse Fourier transform of Eq. (74). This yields

$\delta\sigma_{\alpha\beta}(\vec{r})=\mathcal{F}^{-1}\left[\delta\tilde{\sigma}_{\alpha\beta}(\vec{k})\right]=\sum_{\vec{r}^{\prime}}S_{\alpha\beta}(\vec{r}-\vec{r}^{\prime})\delta
a(\vec{r}^{\prime}).$ (40)

where $S_{\alpha\beta}$ is Green’s function for the change in local stresses.
Next, we can write the form for the Fourier transform of the change in the
local pressure as

$\delta\tilde{P}(\vec{k})=d^{-1}\sum_{\alpha=1}^{d}\delta\tilde{\sigma}_{\alpha\alpha}(\vec{k})=d^{-1}\delta\tilde{a}(\vec{k})\sum_{\alpha}S_{\alpha\alpha}(\vec{k}).$
(41)

Since the Fourier transform of the change in the local stresses is linearly
proportional to $\delta\tilde{a}(\vec{k})$, we can also derive the
configurational average of the local stress and pressure correlations as

$\displaystyle\langle\delta\tilde{\sigma}_{\alpha\beta}(\vec{k})\delta\tilde{\sigma}_{\mu\nu}$
$\displaystyle(\vec{k}^{\prime})\rangle=\left<\delta\tilde{a}(\vec{k}).\delta\tilde{a}(\vec{k}^{\prime})\right>S_{\alpha\beta}(\vec{k})S_{\mu\nu}(\vec{k}^{\prime}),$
(42) $\displaystyle\langle\delta\tilde{P}(\vec{k})\delta\tilde{P}$
$\displaystyle(\vec{k}^{\prime})\rangle=\frac{\left<\delta\tilde{a}(\vec{k}).\delta\tilde{a}(\vec{k}^{\prime})\right>}{d^{2}}\sum_{\alpha,\beta}S_{\alpha\alpha}(\vec{k})S_{\beta\beta}(\vec{k}^{\prime}).$

Using the translational invariance of the system, the configurational average
of the microscopic correlations of $\delta\tilde{a}(\vec{k})$ between two
points $\vec{k}$ and $\vec{k}^{\prime}$ in Fourier space can be written as,

$\displaystyle\left<\delta\tilde{a}(\vec{k}).\delta\tilde{a}(\vec{k}^{\prime})\right>=$
$\displaystyle\frac{N\eta^{2}}{48}\delta_{\vec{k},-\vec{k}^{\prime}},\hskip
42.67912pt\text{(Harmonic)},$ (43) $\displaystyle=$
$\displaystyle\frac{N\eta^{2}(\lambda_{SL}-\lambda_{SS})^{2}}{4}\delta_{\vec{k},-\vec{k}^{\prime}}\hskip
14.22636pt\text{(LJ)}.$

In the $|\vec{k}|\to 0$ limit,
$\langle\delta\tilde{\sigma}_{\alpha\beta}(\vec{k})\delta\tilde{\sigma}_{\mu\nu}(-\vec{k})\rangle$
becomes independent of the magnitude of $|\vec{k}|$ and only has an angular
dependence which we represent as

$C_{\alpha\beta\mu\nu}(\theta)=\lim_{|\vec{k}|\to
0}\left<\delta\tilde{\sigma}_{\alpha\beta}(\vec{k}).\delta\tilde{\sigma}_{\mu\nu}(-\vec{k})\right>.$
(44)

The observed stress correlations are therefore anisotropic in the $k\to 0$
limit, corresponding to a pinch-point singularity at $k=0$ [41]. Due to the
finite system size, and to avoid effects introduced by the periodic
boundaries, we have integrated the stress correlations in Fourier space in the
narrow window of $0.5\leq|\vec{k}|\leq 1.5$. The integrated stress
correlations in a small window of $k\in[k_{min},k_{max}]$ near $k\to 0$ can be
expressed as

$\bar{C}_{\alpha\beta\mu\nu}(\theta)=\int_{k_{min}}^{k_{max}}dk\langle\delta\tilde{\sigma}_{\alpha\beta}(k,\theta)\delta\tilde{\sigma}_{\mu\nu}(k,\pi+\theta)\rangle.$
(45)

In real space, this translates to integrating the stress correlations at
intermediate to large lengthscales. The angular dependence of the integrated
correlations are plotted in Fig. 5 for both Harmonic and LJ model.

###  6.1 Local stress fluctuations

We can express the correlation of the excess local stress between two points
$\vec{r}$ and $\vec{r}^{\prime}$ in real space using the expression for stress
correlation in $k$-space given in Eq. (42) as,

$\displaystyle\langle\delta\sigma_{\mu\nu}(\vec{r})$
$\displaystyle\delta\sigma_{\mu\nu}(\vec{r}^{\prime})\rangle=\sum_{\vec{k},\vec{k}^{\prime}}\frac{\langle\delta\tilde{\sigma}_{\mu\nu}(\vec{k})\delta\tilde{\sigma}_{\mu\nu}(\vec{k}^{\prime})\rangle}{N^{2}}e^{-i(\vec{k}.\vec{r}+\vec{k}^{\prime}.\vec{r}^{\prime})}=\frac{\eta^{2}}{48N}\sum_{\vec{k}}S_{\mu\nu}(\vec{k})S_{\mu\nu}(-\vec{k})e^{-i(\vec{r}-\vec{r}^{\prime}).\vec{k}}.$
(46)

Figure 3:  Local stress fluctuation with increasing polydispersity for $(a)$
disordered triangular lattice (2d) with repulsive harmonic particles of system
size, $N=256$ and initial packing fraction, $\phi=0.92$ and $(b)$ disordered
fcc lattice (3d) with $N=4000$ and $\phi=0.80$. Here all the numerical stress
fluctuation varies as $c_{\alpha\beta}\eta^{2}$, where,
$c_{\alpha\beta}=\frac{1}{48N}\sum_{\vec{k}}S_{\alpha\beta}(\vec{k})S_{\alpha\beta}(-\vec{k})$.

Therefore the local stress fluctuations at the same site can be written as

$\langle\delta\sigma^{2}_{\mu\nu}(\vec{r})\rangle=\frac{\eta^{2}}{48N}\sum_{\vec{k}}S_{\alpha\beta}(\vec{k})S_{\mu\nu}(-\vec{k}).$
(47)

The theoretical prediction for local stress fluctuations as a function of
increasing polydispersity is presented in Figure 3$(a)$ and $(b)$ for a two
dimensional triangular lattice and a three dimensional fcc lattice,
respectively. These theoretical stress fluctuations align perfectly with the
numerical results.

###  6.2 Stress correlations in two dimensional systems

In the case of the two dimensional Harmonic model, every grain has six
neighbors in the near-crystalline system. The expressions for the Green’s
functions $S_{\alpha\beta}$, as defined in Eq. (74), depend solely on the
nearest neighbor arrangement. Our analytic results demonstrate that these
Green’s functions in Fourier space have no radial dependence at small values
of $|\vec{k}|$, corresponding to larger lengthscales in real space. Keeping
only the first term in the Taylor expansion we can write these Green’s
functions as

$\displaystyle S_{xx}=\frac{\mathcal{C}(R_{0},K)}{|\vec{k}|^{2}}$
$\displaystyle\left[\left(\frac{R_{0}}{a_{0}}-1\right)k_{y}^{2}+\left(2-\frac{R_{0}}{a_{0}}\right)k_{x}^{2}\right]+\mathcal{O}(k^{2}),$
(48) $\displaystyle S_{yy}=\frac{\mathcal{C}(R_{0},K)}{|\vec{k}|^{2}}$
$\displaystyle\left[\left(\frac{R_{0}}{a_{0}}-1\right)k_{x}^{2}+\left(2-\frac{R_{0}}{a_{0}}\right)k_{y}^{2}\right]+\mathcal{O}(k^{2}),$
$\displaystyle S_{xy}=\frac{\mathcal{C}(R_{0},K)}{|\vec{k}|^{2}}$
$\displaystyle\left[\left(3-\frac{2R_{0}}{a_{0}}\right)k_{x}k_{y}\right]+\mathcal{O}(k^{2}),$
$\displaystyle\delta\tilde{P}=\frac{\mathcal{C}(R_{0},K)}{2}$
$\displaystyle\delta\tilde{a}(\vec{k})\left[1+f_{1}|\vec{k}|^{2}+f_{2}|\vec{k}|^{2}+\ldots\right],$

where,

$\mathcal{C}(R_{0},K)=\frac{6KR_{0}(R_{0}-a_{0})}{(2R_{0}-a_{0})a_{0}^{2}},\hskip
8.5359ptf_{1}=-\left(a_{0}/3\right)^{2}\cos^{2}{3\theta}\text{ and
}f_{2}=f_{1}^{2}.$ (49)

Here the functions $f_{1}$ and $f_{2}$ in the expression for $\delta\tilde{P}$
are influenced by the initial packing fraction and the orientation of the
reciprocal lattice vector $\vec{k}$. The above expressions for Green’s
functions can be reformulated in a polar coordinate system to describe
behavior at large lengthscales as follows.

$\displaystyle\lim_{|\vec{k}|\to
0}S_{xx}(|\vec{k}|,\theta)\sim\frac{\mathcal{C}(R_{0},K)}{2}\left(1-\left(\frac{2R_{0}}{a_{0}}-3\right)\cos(2\theta)\right),$
(50) $\displaystyle\lim_{|\vec{k}|\to
0}S_{yy}(|\vec{k}|,\theta)\sim\frac{\mathcal{C}(R_{0},K)}{2}\left(1+\left(\frac{2R_{0}}{a_{0}}-3\right)\cos(2\theta)\right),$
$\displaystyle\lim_{|\vec{k}|\to
0}S_{xy}(|\vec{k}|,\theta)\sim-\frac{\mathcal{C}(R_{0},K)}{2}\left(\frac{2R_{0}}{a_{0}}-3\right)\sin(2\theta).$

Figure 4: Change in local stress due to the introduction of disorder in a
single particle (i.e., the Green’s function for change in the local stress) in
a 2D near-crystalline packing (HCP) of soft particles. Here figure (a), (b),
and (c) correspond to $\delta\sigma_{xx}(\vec{r})$,
$\delta\sigma_{xy}(\vec{r})$ and $\delta\sigma_{yy}(\vec{r})$ for $\delta
a(\vec{r}=0)=1$ and $\delta a(\vec{r}\neq 0)=0$ for every other grain.

All the Green’s function above and their correlations show anisotropic
behavior with an angular periodicity of $\pi$. Performing an inverse Fourier
transform on the above expression yields the behavior of these Green’s
functions at large lengthscales in real space, which is presented below as

$\displaystyle S_{xx}(\vec{r})=-S_{yy}(\vec{r})$
$\displaystyle\sim-\frac{\mathcal{C}(R_{0},K)}{2}\left(\frac{2R_{0}}{a_{0}}-3\right)\frac{\cos(2\theta)}{|\vec{r}|^{2}},$
(51) $\displaystyle S_{xy}(\vec{r})=S_{yx}(\vec{r})$
$\displaystyle\sim-\frac{\mathcal{C}(R_{0},K)}{2}\left(\frac{2R_{0}}{a_{0}}-3\right)\frac{\sin(2\theta)}{|\vec{r}|^{2}}.$

Fig. 4 displays all the components of the Green’s functions at large
lengthscales whose functional forms are given in Eq. (51). The long-ranged two
dimensional LJ model also exhibits similar behavior at large lengthscales. The
Green’s functions as well as stress correlations in LJ model have equal
angular behavior as that of the Harmonic model as plotted in Figs. 4 and 5,
only difference lying in the magnitude of the stress fluctuations which
depends on the initial macroscopic properties of these systems like global
pressure, box size. In the two dimensional Harmonic model, all six unique
stress correlations as well as the pressure correlation in Fourier space at
small magnitudes of $\vec{k}$ can be represented using the expressions given
in Eqs. (42) and (48). For example, the correlation of change in local
pressure at large lengthscales ($k\to 0$) can be expressed as,

Figure 5: Angular dependence of stress correlations in Fourier space
integrated in a small window of $|\vec{k}|$ i.e
$\bar{C}_{\alpha\beta\mu\nu}(\theta)=\int_{k_{min}}^{k_{max}}dk\langle\delta\tilde{\sigma}_{\alpha\beta}(k,\theta)\delta\tilde{\sigma}_{\mu\nu}(k,\pi+\theta)\rangle$
for different initial over-compression (Harmonic)/pressure (LJ) with disorder
strength of $\eta=0.005$. The first row represents all the six distinct
correlations for the Harmonic model and the second row corresponds to the LJ
model. Here the solid and dashed lines correspond to the theoretical results
where points correspond to the numerical data for the two different models.
Here we have chosen $k_{min}=0.1$ and $k_{max}=1.0$ and system size $N=6400$
for both the models.

$\displaystyle\mathcal{P}(R_{0},\eta)$ $\displaystyle=\lim_{k\to
0}\left<\delta\tilde{P}(\vec{k}).\delta\tilde{P}(-\vec{k})\right>\sim
N\langle\delta a^{2}\rangle\frac{\mathcal{C}^{2}(R_{0},K)}{4},$ (52)

which is a constant for a given packing of particles. The higher-order terms
in the Taylor expansion of pressure correlation in Fourier space reveal the
anisotropic crystalline nature of the inherent lattice. However, the
significance of higher-order terms becomes apparent at smaller length scales,
where lattice symmetry becomes a dominant factor, as illustrated in Fig. 6. In
Fig. 5 we show an exact match between the stress correlations obtained from
numerical simulations and the predictions from the microscopic theory for both
short-ranged Harmonic and long-ranged LJ model in two dimensions. The
aforementioned correlations are applicable to systems that possess a finite
average pressure ($R_{0}<2a_{0}$). However, in the limit where the average
pressure is zero ($\langle P\rangle\to 0$), i.e., as $R_{0}$ approaches
$2a_{0}$, the results obtained from the VCTG framework [41, 18] for amorphous
materials lacking any crystalline symmetries are reproduced i.e.

$\displaystyle\lim_{R_{0}\to 2a_{0}}\left(\lim_{k\to 0}S_{\alpha\beta}\right)$
$\displaystyle=\frac{\mathcal{C}(2a_{0},K)}{(-1)^{1-\delta_{\alpha\beta}}k^{2}}\prod_{i=\alpha,\beta}\left(\sqrt{k^{2}-k^{2}_{i}}\right),$
(53)

which can be used to write the correlations between the different components
of the stresses as,

$\displaystyle C_{\alpha\beta\mu\nu}=$
$\displaystyle\frac{N\eta^{2}}{48}\left[\lim_{R_{0}\to 2a_{0}}\left(\lim_{k\to
0}S_{\alpha\beta}(\vec{k})S_{\mu\nu}(-\vec{k})\right)\right]=\frac{4\mathcal{P}(2a_{0},\eta)}{(-1)^{2-\delta_{\alpha\beta}-\delta_{\mu\nu}}k^{4}}\prod_{i=\alpha,\beta,\mu,\nu}\left(\sqrt{k^{2}-k^{2}_{i}}\right).$
(54)

where, $\mathcal{C}(2a_{0},K)=4K/a_{0}$ and
$\mathcal{P}(2a_{0},\eta)=\eta^{2}K^{2}/12a_{0}^{2}$. In the same $R_{0}\to
2a_{0}$ limit, the functions in Eq. (48) for change in local pressure have the
following form i.e., $f_{1}=-\left(a_{0}/3\right)^{2}\cos^{2}{3\theta}$,
$f_{2}=f_{1}^{2}$ and so on. For a finite system with particle radius $a_{0}$,
the maximum magnitude of $k$ is $\pi/a_{0}$. So $\left|f_{1}k^{2}\right|\leq
1$ for all values of $k$, except near the boundary of the $1^{st}$ Brillouin
zone. So the change in the local pressure due to particle size defect for
$k\ll\pi/a_{0}$ can be written as,

$\displaystyle\lim_{R_{0}\to 2a_{0}}$
$\displaystyle\delta\tilde{P}(\vec{k})=\frac{\mathcal{C}(2a_{0},K)\delta\tilde{a}(\vec{k})}{2}(1-f_{1}k^{2}+f_{1}^{2}k^{4}+...)$
(55)
$\displaystyle\sim\frac{\mathcal{C}(2a_{0},K)\delta\tilde{a}(\vec{k})}{2(1+f_{1}k^{2})}=\frac{\mathcal{C}(2a_{0},K)\delta\tilde{a}(\vec{k})}{2\left(1+\left(a_{0}/3\right)^{2}k^{2}\cos^{2}{3\theta}\right)}.$

Using the above approximation, we can calculate the correlation of local
pressure which have the sixfold symmetry (i.e $\cos^{2}{3\theta}$) for
$|\vec{k}|>0$, which we have shown in Fig. 6 for a single particle disorder.
In the zero average pressure limit, the findings for the stress correlations
exhibit unexpected universal behavior. This universal behavior is
characterized by the observation of similar anisotropic stress correlations at
large length scales across various jammed athermal packings [41, 18, 17, 56,
57].

Figure 6: Correlation in the change in local pressure in Fourier space
($\langle\delta\tilde{P}(\vec{k})\delta\tilde{P}(-\vec{k})\rangle/\langle\delta
a^{2}\rangle$) due to introduction of a single particle disorder in a
crystalline packing of soft particles with $\phi=0.92$.

###  6.3 Comparision with amorphous systems at large lengthscales

Recently the stress correlations in fully amorphous packings have been
successfully predicted within a field-theoretic framework [41, 18]. This
Vector Charge Theory of “emergent” elasticity is defined by the following
equations:

$\displaystyle\partial_{i}\Sigma_{ij}$ $\displaystyle=$ $\displaystyle f_{j},$
(56) $\displaystyle E_{ij}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\delta_{i}\psi_{j}+\delta_{j}\psi_{i}),$ (57)
$\displaystyle\sigma_{ij}$ $\displaystyle=$
$\displaystyle(\delta_{ijkl}+\chi_{ijkl})E_{kl}=\Lambda^{-1}_{ijkl}E_{kl}.$
(58)

The stress tensor field is represented by $\sigma$, and in this context, $E$
serves a role that is similar to that of the strain field in canonical
elasticity theory. The emergent elasticity modulus tensor is defined as
$\Lambda$, and the equations bear a notable resemblance to those of canonical
linear elasticity. The components of the $\Lambda$ tensor can be interpreted
as "emergent" elastic moduli. By utilizing the tensor gauge theory for
polarizable isotropic media, all 6 distinct stress correlations at small
magnitudes of $\vec{k}$ can be obtained which has the same form as the ones in
Eq. (54) for near-crystalline systems in the $P\to 0$ limit with
$\mathcal{P}(2a_{0},\eta)$ replaced by a constant $K_{2D}$ that depends on the
elastic properties of the system [41, 18]. So for the finite pressure, the
stress correlations for the near-crystalline systems can be represented as a
summation of the stress correlation for an isotropic amorphous system (from
the VCTG framework) and the non-isotropic part of the stress correlation which
contains the information about the crystalline symmetry i.e.

$C_{\alpha\beta\mu\nu}=C_{\alpha\beta\mu\nu}^{t}+d_{\alpha\beta\mu\nu},$ (59)

where $C_{\alpha\beta\mu\nu}^{t}$, represents the stress correlations obtained
from the VCTG framework. This finite shift $d_{\alpha\beta\mu\nu}$, can also
be seen from the numerically obtained correlations as plotted in Fig. 5 which
can not be explained in the VCTG framework. Therefore for an overcompressed
disordered crystal with finite average pressure, the angular behavior of
stress correlations show an additional anisotropic term which depend on the
system pre-stress. For small pre-stress i.e., $R_{0}\to 2a_{0}$, we have
$|d_{\alpha\beta\mu\nu}/C_{\alpha\beta\mu\nu}|\to 0$, and consequently
numerically the stress correlations are indistinguishable to that of the
amorphous systems. The exact expressions for the angular dependence of the
above correlation functions are detailed in Appendix B.

###  6.4 Continuum limit

In the small $|\vec{k}|$ limit, using the expression for
$S_{\alpha\beta}(\vec{k})$ as given in Eq. (48) we can rewrite the change in
the local stresses in the Fourier space as,

$\displaystyle\delta\tilde{\sigma}_{\alpha\beta}=\left[k_{y}^{2}\phi_{\alpha\beta
1}+k_{x}^{2}\phi_{\alpha\beta 2}-k_{x}k_{y}\phi_{\alpha\beta
3}\right]|\vec{k}|^{-2}\delta\tilde{a}(\vec{k}).$ (60)

Using Voigt notation [67], we can replace $xx=1,yy=2,xy=3$. Since the stress
tensor has three independent components in two dimensions, we can represent
the Fourier transform of their change in the following matrix form,

$\displaystyle\lim_{k\to 0}\begin{bmatrix}\delta\tilde{\sigma}_{1}(\vec{k})\\\
\delta\tilde{\sigma}_{2}(\vec{k})\\\
\delta\tilde{\sigma}_{3}(\vec{k})\end{bmatrix}=$
$\displaystyle\underbrace{\left(\frac{\delta\tilde{a}(\vec{k})}{|\vec{k}|^{2}}\right)}_{\tilde{\Psi}(\vec{k})}\underbrace{\begin{bmatrix}\phi_{11}&\phi_{12}&\phi_{13}\\\
\phi_{21}&\phi_{22}&\phi_{23}\\\
\phi_{31}&\phi_{32}&\phi_{33}\end{bmatrix}}_{\hat{\Phi}}.\underbrace{\begin{bmatrix}k_{y}^{2}\\\
k_{x}^{2}\\\ -k_{x}k_{y}\end{bmatrix}}_{\ket{A(\vec{k})}},$ (61)
$\displaystyle\ket{\delta\tilde{\sigma}(\vec{k})}=$
$\displaystyle\hat{\Phi}\left(\tilde{\Psi}(\vec{k})\ket{A(\vec{k})}\right),$

where,

$\displaystyle\phi_{11}$
$\displaystyle=\phi_{22}=\mathcal{C}(R_{0},K)\left(R_{0}/a_{0}-1\right),$ (62)
$\displaystyle\phi_{33}$
$\displaystyle=\mathcal{C}(R_{0},K)\left(2R_{0}/a_{0}-3\right),$
$\displaystyle\phi_{12}$
$\displaystyle=\phi_{21}=\mathcal{C}(R_{0},K)\left(2-R_{0}/a_{0}\right),$
$\displaystyle\phi_{13}$ $\displaystyle=\phi_{23}=\phi_{31}=\phi_{32}=0.$

Given that $\delta\tilde{\sigma}_{i}(\vec{k})$ characterizes the Fourier
transform of the change in local stresses as $|\vec{k}|\to 0$, its inverse-
Fourier transform reveals the changes in local stresses at large lengthscales
due to the defect at the origin. Specifically, this pertains to the coarse-
grained local stress fluctuations in real space, which can be expressed as:

$\delta\sigma_{i}(\vec{r})=\left[\phi_{i1}\partial_{y}^{2}+\phi_{i2}\partial_{x}^{2}-\phi_{i3}\partial_{x}\partial_{y}\right]\Psi(\vec{r}),$
(63)

where,

$\Psi(\vec{r})=\frac{1}{(2\pi)^{2}}\int
d^{2}ke^{i\vec{r}.\vec{k}}\left(\frac{\delta\tilde{a}(\vec{k})}{|\vec{k}|^{2}}\right).$
(64)

The summation in the above equation is performed over reciprocal lattice
points in a triangular lattice configuration. The function $\Psi(\vec{r})$ is
a non-isotropic field contingent upon the symmetry of the lattice. In an
isotropic system, the force balance criterion is:

$\partial_{i}\sigma_{ij}(\vec{r})=f^{j}(\vec{r})=0.$ (65)

In the case of an amorphous system, we can establish an isotropic field
$\Psi^{{}^{\prime}}(\vec{r})$, which lacks any lattice symmetry, yet meets the
force balance condition as expressed in Eq. (65), in a manner such that

$\displaystyle\delta\sigma_{ij}^{{}^{\prime}}(\vec{r})=$
$\displaystyle\epsilon_{ia}\epsilon_{jb}\partial_{a}\partial_{b}\Psi^{{}^{\prime}}(\vec{r}).$
(66)

In the Voigt notation [67],

$\displaystyle\delta\sigma_{i}^{{}^{\prime}}(\vec{r})=$
$\displaystyle\left[\phi_{i1}^{{}^{\prime}}\partial_{y}^{2}+\phi_{i2}^{{}^{\prime}}\partial_{x}^{2}-\phi_{i3}^{{}^{\prime}}\partial_{x}\partial_{y}\right]\Psi^{{}^{\prime}}(\vec{r}),$
(67) $\displaystyle\text{where, }\phi_{ij}^{{}^{\prime}}=$
$\displaystyle\lim_{R_{0}\to
2a_{0}}\phi_{ij}=\mathcal{C}(2a_{0},K)\delta_{ij}.$

So in large length scales the $\hat{\Phi}$ tensor gives a difference between
the coarse-grained local stresses for an amorphous system and a disordered-
crystal. For an amorphous system $\hat{\Phi}$ is an identity matrix whereas in
near-crystalline systems $\hat{\Phi}$ also contains off-diagonal elements.
This difference vanishes for a marginally jammed disordered crystal i.e.,
$R_{0}\to 2a_{0}$.

Figure 7: Angular variation of the change in local stress correlation in
Fourier space due to particle size disorder in a 3D near-crystalline (fcc)
system as $|\vec{k}|\to 0$. Here
$C_{\alpha\beta\mu\nu}(\theta,\phi)=\lim_{|\vec{k}|\to
0}\left<\delta\tilde{\sigma}_{\alpha\beta}(\vec{k}).\delta\tilde{\sigma}_{\mu\nu}(-\vec{k})\right>=\frac{N\eta^{2}}{48}S_{\alpha\beta}S_{\mu\nu}$,
where exact expressions for $S_{\alpha\beta}$ and its angular variations are
given in Eqs. (86), (87) and (88). Here we have plotted correlations the in
Hammer projection as given in eq. (92) for
$(H_{x}/2\sqrt{2})^{2}+(H_{y}/\sqrt{2})^{2}\leq 1$.

###  6.5 Stress correlations in three dimensions

In this section, we derive the stress correlations in a three dimensional
system induced by microscopic disorder. The method developed earlier for the
displacement fields and the change in local stresses is still valid for the 3D
fcc lattice with the only difference arising from the number of nearest
neighbors and their arrangement in space. We can also write the displacement
fields of the particles as,
$\delta\tilde{r}^{\alpha}(\vec{k})=G^{\alpha}(\vec{k})\delta\tilde{a}(\vec{k})$
[61]. Using this we can write down the expression for change in the local
stress components in an fcc lattice in Fourier space as,
$\delta\tilde{\sigma}_{\alpha\beta}(\vec{k})=S_{\alpha\beta}(\vec{k})\delta\tilde{a}(\vec{k})$.
The exact expressions for Green’s functions for displacement fields and the
change in local stress components due to both particle size disorder and force
pinning in fcc lattice are detailed in the appendix A.1 and appendix A.2
respectively. At large lengthscales, the Green’s function in real space for
change in local stress has the following radial behavior,

$S_{\alpha\beta}(\vec{r})=\mathcal{F}\left[S_{\alpha\beta}(\vec{k})\right]\sim\frac{B_{\alpha\beta}(\theta,\phi)}{r^{3}},\hskip
14.22636pt\forall r\gg R_{0}.$ (68)

Similar to two dimensional systems, we can write the correlation between
different components of stress using the Green’s functions defined above in
$|\vec{k}|\to 0$ limit as

$\displaystyle C_{\alpha\beta\mu\nu}(\theta,\phi)=$
$\displaystyle\langle\delta a^{2}\rangle\lim_{|\vec{k}|\to
0}S_{\alpha\beta}(\vec{k})S_{\mu\nu}(-\vec{k}).$ (69)

Here the preliminary theoretical results for the Green’s functions for local
stress and the stress correlations in Fourier space are given in the Appendix
A.2 and all the distinct components of these correlations are plotted in
Figure. 7.

##  7 Response to a point force

Subsequently, we examine the change in the local stress components within a
crystalline system caused by finite quenched forces. We start by introducing
finite external quenched forces, represented as $\vec{f}_{a}(\vec{r}_{i})$, to
each grain $i$ in the crystalline system. The sum of these forces satisfies
the condition $\sum_{i}\vec{f}_{a}(\vec{r}_{i})=0$. To balance the forces on
each grain in the system, particles shift from their original lattice
positions. In the case of an ideal crystalline system with forces attached to
particles, the force balance requirement for each grain $i$ can be expressed
as follows:

$\displaystyle\sum_{j}\left(f^{(0)\mu}_{ij}+\delta f^{\mu}_{ij}\right)=$
$\displaystyle-(f_{a})^{\mu}_{i},$ (70)

where $f_{ij}^{(0)}$ are the forces along the bond between particle $i$ and
$j$ in the initial crystalline system whereas $\delta f_{ij}$ correspond to
the change in the bond forces due to the external pinning. Here $\delta
f_{ij}$ can be approximated using their first-order Taylor series expansion,
resulting in the following linear expression:

$\displaystyle\sum_{j}\sum_{\nu}C_{ij}^{\mu\nu}\delta r^{\nu}_{ij}=$
$\displaystyle-(f_{a})^{\mu}_{i}.$ (71)

The above expression is similar to the linear equation for the displacement
fields due to particle size disorder. By employing a similar approach as
explained earlier for particle size disorder, we can obtain the expression for
displacement fields resulting from pinned forces, which is presented as
follows:

$\delta\tilde{r}^{\mu}(\vec{k})=-\sum_{\nu}\underbrace{\left(A^{-1}\right)^{\mu\nu}(\vec{k})}_{\tilde{\mathcal{G}}^{\mu\nu}(\vec{k})}\tilde{f}^{\nu}_{a}(\vec{k}),$
(72)

whose inverse Fourier transform gives the displacement fields in real space
due to force quench. The large lengthscale behavior of these displacements has
been studied thoroughly [32]. We have given a brief description of these above
expressions in the Appendix A.1. Using these displacement fields and the
expression for the change in local stresses, one can write the components of
change in local stresses in Fourier space as,

$\delta\tilde{\sigma}_{\alpha\beta}(\vec{k})=\sum_{\mu}\mathcal{S}^{\mu}_{\alpha\beta}(\vec{k})\tilde{f}_{a}^{\mu}(\vec{k}),\text{
where},$ (73)

$\displaystyle\mathcal{S}_{\alpha\beta}^{\mu}(\vec{k})=$
$\displaystyle\sum_{j}\left[\left(-1+F_{j}(\vec{k})\right)\left(\sum_{\nu}\Delta^{\alpha(0)}_{j}C_{j}^{\beta\nu}\tilde{\mathcal{G}}^{\mu\nu}(\vec{k})+f_{j}^{\beta(0)}\tilde{\mathcal{G}}^{\mu\alpha}(\vec{k})\right)\right].$
(74)

Figure 8:  Change in local stress components per unit force
($\delta\sigma_{\alpha\beta}(\vec{r})/f_{a}$) applied at the origin in 2D
near-crystalline systems of $N=10000$ particles with periodic boundary
conditions. Here the top panel $(a)$ corresponds to direct numerical
simulations and $(b)$ represents the theoretical results as given in Eq. (78).

We can connect the Green’s functions for the displacement and stress fields
produced due to force quench to that of the microscopic particle size disorder
which can be represented as follows,

$\displaystyle\tilde{G}^{\mu}(\vec{k})=$
$\displaystyle\sum_{\nu}\tilde{\mathcal{G}}^{\mu\nu}(\vec{k})\left(\sum_{j}C_{j}^{\nu
a}\left[1+F_{j}(\vec{k})\right]\right),$ (75) $\displaystyle
S_{\alpha\beta}(\vec{k})=$
$\displaystyle\sum_{\nu}\mathcal{S}^{\nu}_{\alpha\beta}(\vec{k})\left(\sum_{j}C_{j}^{\nu
a}\left[1+F_{j}(\vec{k})\right]\right).$

In real space, the change in the local stress components can be represented as

$\delta\sigma_{\alpha\beta}(\vec{r})=\sum_{\vec{r}^{\prime}}\sum_{\mu}\mathcal{S}^{\mu}_{\alpha\beta}(\vec{r}-\vec{r}^{\prime})f^{\mu}_{a}(\vec{r}^{\prime}),$
(76)

where $\mathcal{S}^{\mu}_{\alpha\beta}(\vec{r})$ corresponds to the Green’s
function for the change in local stress components due to pinned force in real
space. These Green’s functions in the above expression can be understood as
the change in local stresses due to a unit force applied to a single particle
at the origin. The exact expressions for these Green’s functions in Fourier
space (in $k\to 0$ limit) are detailed in Appendix. A.2. We can perform
Fourier transform to obtain the large lengthscale behavior of these Green’s
functions which shows $1/r$ radial behavior in two dimensions and
$\mathcal{S}_{\alpha\beta}^{\gamma}\sim 1/r^{2}$ and
$\mathcal{S}_{\alpha\alpha}^{\alpha}\sim 1/r$ in three dimensions. In two
dimensions we can write these Green’s functions explicitly as

$\displaystyle\mathcal{S}^{x}_{xx}(\vec{r})=$
$\displaystyle\left(4\left(R_{0}-a_{0}\right)\cos{\theta}+a_{0}\cos{3\theta}\right)/(2R_{0}-a_{0})r,$
(77) $\displaystyle\mathcal{S}^{x}_{xy}(\vec{r})=$
$\displaystyle\left(2R_{0}-a_{0}+2a_{0}\cos{2\theta}\right)\sin{\theta}/(2R_{0}-a_{0})r,$
$\displaystyle\mathcal{S}^{x}_{yy}(\vec{r})=$ $\displaystyle
a_{0}\cos{3\theta}/(2R_{0}-a_{0})r,$
$\displaystyle\mathcal{S}^{y}_{xx}(\vec{r})=$ $\displaystyle
a_{0}\sin{3\theta}/(2R_{0}-a_{0})r,$
$\displaystyle\mathcal{S}^{y}_{xy}(\vec{r})=$
$\displaystyle\left(2R_{0}-a_{0}-2a_{0}\cos{2\theta}\right)\cos{\theta}/(2R_{0}-a_{0})r,$
$\displaystyle\mathcal{S}^{y}_{yy}(\vec{r})=$
$\displaystyle\left(4\left(R_{0}-a_{0}\right)\sin{\theta}-a_{0}\sin{3\theta}\right)/(2R_{0}-a_{0})r.$

To verify our results we have taken a system of $N=10000$ particles in a
triangular lattice arrangement and assigned a point force ($f_{a}\hat{y}$)
directed along the $y-$axis to a particle located at the origin ($0,0$). To
make the net force in the system zero we apply an additional $-f_{a}\hat{y}/N$
to all the particles in the system. All the particles will rearrange to
balance the point force, which results in a change in the local stress
profile. So the displacements and the change in local stresses can be written
as,

$\displaystyle\delta x(\vec{r})=-f_{a}\mathcal{G}^{xy}(\vec{r}),\hskip
14.22636pt\delta y(\vec{r})=-f_{a}\mathcal{G}^{yy}(\vec{r}),$ (78)
$\displaystyle\delta\sigma_{\alpha\beta}(\vec{r})=f_{a}\mathcal{S}^{y}_{\alpha\beta}(\vec{r}).$

The results above can be verified in Fig. 8 where we have shown the match
between the analytically and numerically obtained results for change in local
stress profile due to a quenched force along $y$-axis. We have also done a
preliminary study on the effect of pinning in 3d systems which we have
detailed in the appendix A.2.

##  8 Distribution of Stresses

Computations of stress correlations in athermal amorphous materials implicitly
assume a partition function description [41, 18, 57, 58]. We show below that
such an ensemble indeed emerges from the fluctuations of the individual radii.
We obtain the joint probability distribution of changes in local stresses in
Fourier space by utilizing the linear relationships between the change in
local stress and disorder, as described in equations (40) and (76).

$\displaystyle P($
$\displaystyle\delta\tilde{\sigma}_{1}(\vec{k}),\delta\tilde{\sigma}_{2}(\vec{k}),\delta\tilde{\sigma}_{3}(\vec{k}))=\int_{-\infty}^{\infty}d\left(\delta\tilde{a}_{k}\right)p(\delta\tilde{a}_{k})\prod_{i=1}^{3}\delta\left(\delta\tilde{\sigma}_{i}^{R}-S_{i}\delta\tilde{a}(\vec{k})^{R}\right)\delta\left(\delta\tilde{\sigma}_{i}^{I}-S_{i}\delta\tilde{a}(\vec{k})^{I}\right)$
(79)
$\displaystyle=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\prod_{j=1}^{3}df_{j}dg_{j}e^{i\left(f_{j}\delta\tilde{\sigma}_{j}^{R}+g_{j}\delta\tilde{\sigma}_{j}^{I}\right)}\underbrace{\int_{-\infty}^{\infty}d\left(\delta\tilde{a}_{k}\right)p(\delta\tilde{a}_{k})e^{-i\left(\sum_{m=1}^{3}\left(\delta\tilde{a}(\vec{k})^{R}f_{m}+\delta\tilde{a}(\vec{k})^{I}g_{m}\right)S_{m}\right)}}_{h\left(\\{f,g,S\\}\right)}.$

In the above equation, the superscripts $R$ and $I$ correspond to the real and
imaginary parts respectively. Since the $\delta a$s are drawn from a uniform
distribution we can obtain the distribution of its Fourier transformed
variable $\delta\tilde{a}(\vec{k})$ as $p(\delta\tilde{a}_{k})=(48/\pi
N\eta^{2})^{1/2}\exp{-48|\delta\tilde{a}(\vec{k})|^{2}/N\eta^{2}}.$ So using
this distribution we can rewrite the above equation as,

$\displaystyle h$
$\displaystyle\left(\\{f,g,S\\}\right)=\exp{-\frac{N\eta^{2}}{192}\sum_{m,n=1}^{3}(f_{m}f_{n}+g_{m}g_{n})\left(S_{m}S_{n}\right)}=\exp{-\frac{N\eta^{2}}{192}(\bra{f}\hat{S}\ket{f}+\bra{g}\hat{S}\ket{g})},$
(80)

where, $\bra{f}=(f_{1}\text{\hskip 2.84544pt}f_{2}\text{\hskip
2.84544pt}f_{3})$, $\bra{g}=(g_{1}\text{\hskip 2.84544pt}g_{2}\text{\hskip
2.84544pt}g_{3})$ and $\hat{\mathcal{S}}_{mn}=S_{m}S_{n}$. Therefore, using
the above expression we can rewrite the joint probability distribution of the
stress components in Fourier space as

$\displaystyle
P(\delta\tilde{\sigma}_{xx}(\vec{k}),\delta\tilde{\sigma}_{yy}(\vec{k}),\delta\tilde{\sigma}_{xy}(\vec{k}))$
$\displaystyle=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\prod_{j=1}^{3}df_{j}dg_{j}e^{i\left(\bra{f}\ket{\delta\tilde{\sigma}^{R}}+\bra{g}\ket{\delta\tilde{\sigma}^{I}}\right)}e^{-\frac{N\eta^{2}}{192}\left(\bra{f}\hat{S}\ket{f}+\bra{g}\hat{S}\ket{g}\right)}$
(81) $\displaystyle=\left(\frac{48}{\pi
N\eta^{2}}\right)^{1/2}\exp{-\frac{48}{N\eta^{2}}\bra{\delta\tilde{\sigma}^{R}}\mathcal{\hat{S}}^{-1}\ket{\delta\tilde{\sigma}^{R}}}.$

where,

$\bra{\delta\tilde{\sigma}^{R}}=\begin{bmatrix}\delta\tilde{\sigma}_{xx}(\vec{k})&\delta\tilde{\sigma}_{yy}(\vec{k})&\delta\tilde{\sigma}_{xy}(\vec{k})\end{bmatrix},\text{
and
}\mathcal{\hat{S}}=\begin{bmatrix}S_{xx}S_{xx}&S_{xx}S_{yy}&S_{yy}S_{xy}\\\
S_{yy}S_{xx}&S_{yy}S_{yy}&S_{xy}S_{xy}\\\
S_{xy}S_{xx}&S_{xy}S_{yy}&S_{xx}S_{xy}\end{bmatrix},$ (82)

where the form of the $S_{\alpha\beta}$ are model dependent and their exact
forms are provided in the appendix A.2.

Earlier studies [41, 18] have shown that the generalized elastic constants in
amorphous packings are directly related to the correlations in components of
the stress tensor. This explicitly does not depend on the strength of the
disorder in the system. However, as we have shown, in near-crystalline
packings the distributions depend on the strength of the disorder but the
elastic constants are independent of them. So it is not evident that the
formulations developed earlier [41, 18] for generalised elastic constants in
amorphous systems can be directly implemented in the context of near-
crystalline packings.

##  9 Discussion and Conclusion

In this study, we have examined the elastic properties, stress fluctuations as
well as spatial stress correlations in near-crystalline athermal systems. We
have obtained exact theoretical results for the macroscopic elastic properties
by utilizing the fact that the average change in local stresses due to
microscopic disorder in crystalline athermal solids is negligible. Our
findings reveal that these elastic properties remain unaffected by the degree
of disorder within a crystalline packing but are influenced by various initial
conditions, such as packing fraction, pressure, and the strength of particle
interactions. Furthermore, we have presented both numerical and theoretical
results for local stress fluctuations and their spatial correlations within
energy-minimized configurations of soft particles in both two and three
dimensions. Notably, all these fluctuations and correlations exhibit a
quadratic variation with the strength of the disorder. We have shown that the
components of the stress tensor display anisotropic long range decay in both
two and three dimensional near-crystalline packings irrespective of the
interaction potential. For particle size disorder we observe a $1/r^{d}$
radial decay of the change in the stress-tensor components whereas for
external quenched forces, we have established a slower $1/r^{d-1}$ radial
decay. The stress correlations in disordered crystals differ significantly
from those observed in isotropic amorphous materials at high packing fractions
or under high-pressure conditions [57, 58]. Notably, we have observed
additional non-isotropic angular behavior in the stress correlations, which
becomes prominent at higher packing fractions.

Crucially, for the case of near-crystalline packings, we have found that the
correlations depend on the strength of the disorder (proportional to
$\eta^{2}$) introduced into the system, whereas the macroscopic elastic
coefficients are largely independent of disorder. This is in contrast to
stress correlations in amorphous materials, where the magnitude of the
correlations have been related to elastic constants which are independent of
the degree of disorder [41, 18]. It would therefore be very interesting to
examine the crossover between this near-crystalline and amorphous behaviour as
the degree of disorder increased beyond a critical threshold.

Several interesting questions still remain for further research. For example,
the behavior of these stress correlations as we increase the disorder in the
system may not vary quadratically to the strength of the disorder as linear
perturbation expansion is not valid at a high enough disorder. It would
therefore be interesting to study these properties across the crystalline to
amorphous transition. It would also be intriguing to extend our analysis to
dynamical systems where particles obey local force balance constraints only on
average.

Recent studies [68] have connected the fluctuating elastic constants to the
quasilocalised vibrational modes in amorphous materials. This hypothesis would
be interesting to examine microscopically within the context of near-
crystalline materials. Our study also highlights the importance of prestress
in the stress correlations of athermal materials. Increasingly, studies have
shown that accounting for the impact of stresses or prestress is crucial in
understanding mechanical properties of amorphous materials [69, 70]. Our
techniques could therefore be used to test the impact of frozen-in stresses on
the elasticity characteristics of both crystalline and amorphous materials.

## Acknowledgements

We thank Surajit Chakraborty, Pinaki Chaudhuri, Debankur Das, Bulbul
Chakraborty, Subhro Bhattacharjee, Jishnu Nampoothiri and Palash Bera for
useful discussions. The work of K. R. was partially supported by the SERB-
MATRICS grant MTR/2022/000966. This project was funded by intramural funds at
TIFR Hyderabad from the Department of Atomic Energy (DAE), Government of
India.

##  Appendix A Green’s functions in the Harmonic model

###  A.1 Green’s functions for displacement fields

The displacements of every grain (in Fourier space) due to pinned forces can
be written using the translational invariance of the system as detailed in the
recent publications [60, 32] which have the form

$\delta\tilde{r}^{\mu}(\vec{k})=-\sum_{\nu}\underbrace{\left(A^{-1}\right)^{\mu\nu}(\vec{k})}_{\tilde{\mathcal{G}}^{\mu\nu}(\vec{k})}\tilde{f}^{\nu}_{a}(\vec{k}).$
(83)

Similarly, the linear order displacement fields due to disorder in the
particle sizes can be expressed as,

$\delta\tilde{r}^{\alpha}(\vec{k})=\sum_{\nu}\tilde{G}^{\alpha}(\vec{k})\delta\tilde{a}^{\alpha}(\vec{k}).$
(84)

Below we have given the general form of these Green’s functions with various
initial packing fractions, particle size disorder and random quenched forces
for both two and three dimensional harmonic soft particle systems.

Force pinning (2d)
---
$\tilde{\mathcal{G}}^{xx}(\vec{k})$ |  $\frac{R_{0}a_{0}^{2}}{K}\frac{m_{1}}{\left(m_{1}m_{2}-m_{3}^{2}\right)}$ |  where, $m_{1}=-3(R_{0}-a_{0})+(R_{0}-2a_{0})\cos{2k_{x}}+(2R_{0}-a_{0})\cos{k_{x}}\cos{k_{y}},$
$\tilde{\mathcal{G}}^{yy}(\vec{k})$ |  $\frac{R_{0}a_{0}^{2}}{K}\frac{m_{2}}{\left(m_{1}m_{2}-m_{3}^{2}\right)}$ |  $m_{2}=-3(R_{0}-a_{0})+R_{0}\cos{2k_{x}}+(2R_{0}-3a_{0})\cos{k_{x}}\cos{k_{y}},$
$\tilde{\mathcal{G}}^{xy}(\vec{k})$ |  $\frac{R_{0}a_{0}^{2}}{K}\frac{m_{3}}{\left(m_{1}m_{2}-m_{3}^{2}\right)}$ |  $m_{3}=\sqrt{3}a_{0}\sin{k_{x}}\sin{k_{y}}.$
Force pinning (3d)
$\tilde{\mathcal{G}}^{\alpha\alpha}(\vec{k})$ |  $\frac{R_{0}a_{0}^{2}}{K}\frac{p_{\beta}p_{\gamma}-q_{\alpha}^{2}}{h}$ |  where, $h=p_{x}p_{y}p_{z}+2q_{x}q_{y}q_{z}-p_{x}q_{x}^{2}-p_{y}q_{y}^{2}-p_{z}q_{z}^{2},$
$\tilde{\mathcal{G}}^{\alpha\beta}(\vec{k})$ |  $\frac{R_{0}a_{0}^{2}}{K}\frac{p_{\gamma}q_{\gamma}-q_{\alpha}q_{\beta}}{h}$ |  $p_{\alpha}=(-3R_{0}+4a_{0}+(R_{0}-2a_{0})\cos{k_{\beta}}\cos{k_{\gamma}}$ $+(R_{0}-a_{0})\cos{k_{\alpha}}(\cos{k_{\beta}}+\cos{k_{\gamma}})),$ $q_{\alpha}=a_{0}\sin{k_{\beta}}\sin{k_{\gamma}},\hskip 28.45274pt\forall\alpha\neq\beta\neq\gamma\in\\{x,y,z\\}.$
Particle size disorder (2d)
$\tilde{G}^{\alpha}(\vec{k})$ |  $\sum_{\beta}\tilde{\mathcal{G}}^{\alpha\beta}(\vec{k})D^{\beta}(\vec{k})$ |  where, $D^{x}(\vec{k})=i\frac{K(R_{0}-a_{0})}{a_{0}^{3}}\left(2\cos{k_{x}}+\cos{k_{y}}\right)\sin{k_{x}},$ $D^{y}(\vec{k})=i\frac{\sqrt{3}K(R_{0}-a_{0})}{a_{0}^{3}}\cos{k_{x}}\sin{k_{y}}.$
Particle size disorder (3d)
$\tilde{G}^{\alpha}(\vec{k})$ |  $\sum_{\beta}\tilde{\mathcal{G}}^{\alpha\beta}(\vec{k})D^{\beta}(\vec{k})$ |  where, $D^{\alpha}(\vec{k})=i\frac{K(R_{0}-a_{0})}{\sqrt{2}a_{0}^{3}}\left(\cos{k_{\beta}}+\cos{k_{\gamma}}\right)\sin{k_{\alpha}},\hskip 2.84544pt\forall\alpha\neq\beta\neq\gamma\in\\{x,y,z\\}.$

###  A.2 Green’s functions for change in local stresses in Fourier space

In the linear approximation, the change in the local stress due to external
force pinning in both two and three dimensional systems can be expressed in
Fourier space as,

$\delta\tilde{\sigma}_{\mu\nu}(\vec{k})=\sum_{\nu}\mathcal{S}^{\alpha}_{\mu\nu}(\vec{k})\tilde{f}^{\alpha}_{a}(\vec{k}).$
(85)

For particle size disorder the expressions for Green’s functions
$S_{\alpha\beta}$ as given in Eq. (74) which only depend on the nearest
neighbor arrangement exhibit simple relationships at small magnitudes of
$|\vec{k}|$ (corresponding to large lengthscales in real space), which are
expressed below as

Force Pinning (2d)
---
$\mathcal{S}^{\alpha}_{\alpha\alpha}(\vec{k})$ |  $-2ik_{\alpha}\left(k_{\beta}^{2}(4R_{0}-a_{0})+k_{\alpha}^{2}(4R_{0}-5a_{0})\right)/\left(k^{4}(2R_{0}-a_{0})\right),$
$\mathcal{S}^{\alpha}_{\beta\beta}(\vec{k})$ |  $-2ik_{\alpha}(k_{\alpha}^{2}-3k_{\beta}^{2})a_{0}/\left(k^{4}(2R_{0}-a_{0})\right),$
$\mathcal{S}^{\alpha}_{\alpha\beta}(\vec{k})=\mathcal{S}^{\alpha}_{\beta\alpha}(\vec{k})$ |  $-2ik_{\beta}\left(k_{\alpha}^{2}(2R_{0}-5a_{0})+k_{\beta}^{2}(2R_{0}-a_{0})\right)/\left(k^{4}(2R_{0}-a_{0})\right).$
Force Pinning (3d)
$\mathcal{S}^{\alpha}_{\mu\nu}(\vec{k})$ |  $i\sqrt{2}R_{0}n^{\alpha}_{\mu\nu}(\vec{k})/d(\vec{k})$ where, $d(\vec{k})=c_{0}^{2}(c_{0}-a_{0})k^{6}-\frac{3a_{0}^{2}}{2}c_{0}k^{2}\sum_{\alpha\beta\in\\{x,y,z\\}}k_{\alpha}^{2}k_{\beta}^{2}+5a_{0}^{3}(k_{x}k_{y}k_{z})^{2},$ $n^{\alpha}_{\alpha\alpha}(\vec{k})=c_{0}k^{2}\left[c_{0}^{2}k^{2}+2a_{0}c_{1}(k^{2}-k_{\alpha}^{2})\right]-3a_{0}^{2}c_{2}k_{\beta}^{2}k_{\gamma}^{2},$ $n^{\alpha}_{\beta\beta}(\vec{k})=a_{0}k_{\alpha}\left[2c_{0}c_{1}k^{4}+3a_{0}(8c_{3}-a_{0})k_{\beta}^{2}(k_{\alpha}^{2}+k_{\beta}^{2})+a_{0}c_{0}k^{2}k_{\alpha}^{2}\right.$ $\left.\hskip 28.45274pt-(16c_{3}^{2}+a_{0}(2c_{2}+a_{0}))k^{2}k_{\beta}^{2}\right],$ $n^{\alpha}_{\alpha\beta}(\vec{k})=c_{0}k_{\beta}\left[2c_{0}c_{3}k^{4}-12R_{0}a_{0}k^{2}k_{\alpha}^{2}\right.$ $\hskip 28.45274pt\left.-3a_{0}^{2}\left((4k_{\alpha}^{2}-k_{\beta}^{2})(k_{\alpha}^{2}+k_{\beta}^{2})-k^{2}(k_{\alpha}^{2}-k_{\beta}^{2})\right)\right],$ $n^{\alpha}_{\beta\gamma}(\vec{k})=2a_{0}c_{0}k_{\alpha}k_{\beta}k_{\gamma}\left[2c_{0}k_{\alpha}^{2}+(2R_{0}+c_{0})(k_{\beta}^{2}+k_{\gamma}^{2})\right],\hskip 14.22636pt$ $\forall\alpha\neq\beta\neq\gamma\in\\{x,y,z\\}.$ with $c_{0}=(2R_{0}-3a_{0}),c_{1}=(R_{0}-2a_{0})$, $c_{2}=(2R_{0}-a_{0})$, $c_{3}=(R_{0}-a_{0}).$
Particle size disorder (2d)
$S_{\alpha\alpha}(\vec{k})$ |  $\frac{\mathcal{C}(R_{0},K)}{|\vec{k}|^{2}}\left[\left(\frac{R_{0}}{a_{0}}-1\right)k_{\beta}^{2}+\left(2-\frac{R_{0}}{a_{0}}\right)k_{\alpha}^{2}\right],$
$S_{\alpha\beta}(\vec{k})$ |  $\frac{\mathcal{C}(R_{0},K)}{|\vec{k}|^{2}}\left[\left(3-\frac{2R_{0}}{a_{0}}\right)k_{\alpha}k_{\beta}\right],\hskip 14.22636pt\forall\alpha\neq\beta\in\\{x,y\\},$ with $\mathcal{C}(R_{0},K)=\frac{6KR_{0}(R_{0}-a_{0})}{(2R_{0}-a_{0})a_{0}^{2}}.$

Particle size disorder (3d)
---
$S_{\mu\nu}(\vec{k})$ |  $\left(2c_{3}R_{0}/a_{0}^{3}\right)\left(m_{\mu\nu}(\vec{k})/d(\vec{k})\right)$ where, $d(\vec{k})=c_{0}^{2}(c_{0}-a_{0})k^{6}-\frac{3a_{0}^{2}}{2}c_{0}k^{2}\sum_{\alpha\beta\in\\{x,y,z\\}}k_{\alpha}^{2}k_{\beta}^{2}+5a_{0}^{3}(k_{x}k_{y}k_{z})^{2},$ $m_{\alpha\alpha}(\vec{k})=c_{0}^{3}k^{6}-a_{0}^{2}(4c_{1}-a_{0})k_{\alpha}^{2}k_{\beta}^{2}k_{\gamma}^{2}-c_{0}a_{0}k^{2}\left(4c_{1}k_{\alpha}^{4}+a_{0}^{2}k_{\beta}^{2}k_{\gamma}^{2}\right)$ $\hskip 42.67912pt+c_{0}(8c_{1}^{2}+a_{0}c_{0})k^{4}k_{\alpha}^{2},$ $m_{\alpha\beta}(\vec{k})=c_{0}k_{\alpha}k_{\beta}\left(2c_{0}k^{2}-a_{0}(k_{\alpha}^{2}+k_{\beta}^{2})\right)\left(2c_{1}k^{2}+a_{0}(k_{\alpha}^{2}+k_{\beta}^{2})\right),$ $\hskip 42.67912pt\forall\alpha\neq\beta\neq\gamma\in\\{x,y,z\\}.$ with $c_{0}=(2R_{0}-3a_{0}),c_{1}=(R_{0}-2a_{0})$, $c_{2}=(2R_{0}-a_{0})$, $c_{3}=(R_{0}-a_{0}).$

Figure 9: Green’s function for change in the local stress components in $k\to
0$ limit which only has anisotropic angular behavior. Here the Green’s
functions are represented in Hammer projection i.e., $\\{\theta,\phi\to
H_{x}(\theta,\phi),H_{y}(\theta,\phi)\\}$ as given in Eq. (92) for
$(H_{x}/2\sqrt{2})^{2}+(H_{y}/\sqrt{2})^{2}\leq 1$.

In the 3D Harmonic model, when we approach the limit of $|\vec{k}|\to 0$ and
$R_{0}\to 2a_{0}$ with $a_{0}=1/2$, the components of $S_{\alpha\beta}$ take
on specific forms:

$\displaystyle S_{\alpha\alpha}=$
$\displaystyle\frac{8(k_{\beta}^{2}+k_{\gamma}^{2})\left(|k|^{4}-k_{\beta}^{2}k_{\gamma}^{2}\right)}{|k|^{6}+(k_{x}^{6}+k_{y}^{6}+k_{z}^{6}+2k_{x}^{2}k_{y}^{2}k_{z}^{2})},\hskip
5.69046pt\text{and}\hskip
5.69046ptS_{\alpha\beta}=\frac{-8(k_{\alpha}k_{\beta})\left(|k|^{4}-k_{\gamma}^{4}\right)}{|k|^{6}+(k_{x}^{6}+k_{y}^{6}+k_{z}^{6}+2k_{x}^{2}k_{y}^{2}k_{z}^{2})},\hskip
14.22636pt$ (86)
$\displaystyle\forall\alpha\neq\beta\neq\gamma\in\\{x,y,z\\}.$

Now, in a spherical polar coordinate system, where
$k_{x}=k\sin{\theta}\cos{\phi}$, $k_{y}=k\sin{\theta}\sin{\phi}$, and
$k_{z}=k\cos{\theta}$, substituting these values into the equations above
yields the angular behavior of $S_{\alpha\beta}(\vec{k})$ as $|\vec{k}|\to 0$,

$\lim_{|\vec{k}|\to
0}S_{\alpha\beta}(|\vec{k}|,\theta,\phi)=\frac{g_{\alpha\beta}(\theta,\phi)}{h(\theta,\phi)}.$
(87)

In the specific case of $R_{0}\to 2a_{0}$ with $a_{0}=1/2$, the components of
$g_{\alpha\beta}(\theta,\phi)$ and $h(\theta,\phi)$ are given by,

$\displaystyle g_{xx}=$ $\displaystyle
8(\sin^{2}{\theta}\sin^{2}{\phi}+\cos^{2}{\theta})\left(1-\sin^{2}{\theta}\sin^{2}{\phi}\cos^{2}{\theta}\right),$
(88) $\displaystyle g_{yy}=$ $\displaystyle
8(\sin^{2}{\theta}\cos^{2}{\phi}+\cos^{2}{\theta})\left(1-\sin^{2}{\theta}\cos^{2}{\phi}\cos^{2}{\theta}\right),$
$\displaystyle g_{zz}=$ $\displaystyle
8\sin^{2}{\theta}\left(1-\sin^{4}{\theta}\sin^{2}{\phi}\cos^{2}{\phi}\right),$
$\displaystyle g_{xy}=$ $\displaystyle
2\sin^{4}{\theta}\sin^{2}{2\phi}\left(1-\cos^{2}{\theta}\right),$
$\displaystyle g_{yz}=$ $\displaystyle
2\sin^{2}{2\theta}\sin^{2}{\phi}\left(1-\sin^{2}{\theta}\cos^{2}{\phi}\right),$
$\displaystyle g_{zx}=$ $\displaystyle
2\sin^{2}{2\theta}\cos^{2}{\phi}\left(1-\sin^{2}{\theta}\sin^{2}{\phi}\right),$
$\displaystyle h(\theta,\phi$
$\displaystyle)=1+\sin^{6}{\theta}(\sin^{6}{\phi}+\cos^{6}{\phi})+\cos^{6}{\theta}+2\sin^{4}{\theta}\sin^{2}{\phi}\cos^{2}{\phi}.$

##  Appendix B Angular variation of stress correlations at large lengthscales

All the six distinct stress correlations in a two dimensional Harmonic soft
particle system where the underlying arrangement is a triangular lattice have
the following forms,

$\displaystyle C_{xxxx}=4\mathcal{P}(R_{0},\eta)$
$\displaystyle\left[\left(\frac{R_{0}}{a_{0}}-1\right)\sin^{2}{\theta}+\left(2-\frac{R_{0}}{a_{0}}\right)\cos^{2}{\theta}\right]^{2}\xrightarrow{R_{0}\to
2a_{0}}4\mathcal{P}(2a_{0},\eta)\sin^{4}{\theta},$ (89) $\displaystyle
C_{yyyy}=4\mathcal{P}(R_{0},\eta)$
$\displaystyle\left[\left(\frac{R_{0}}{a_{0}}-1\right)\cos^{2}{\theta}+\left(2-\frac{R_{0}}{a_{0}}\right)\sin^{2}{\theta}\right]^{2}\xrightarrow{R_{0}\to
2a_{0}}4\mathcal{P}(2a_{0},\eta)\cos^{4}{\theta},$ $\displaystyle
C_{xyxy}=4\mathcal{P}(R_{0},\eta)$
$\displaystyle\left[\left(3-\frac{2R_{0}}{a_{0}}\right)\sin{\theta}\cos{\theta}\right]^{2}\xrightarrow{R_{0}\to
2a_{0}}4\mathcal{P}(2a_{0},\eta)\sin^{2}{\theta}\cos^{2}{\theta},$
$\displaystyle C_{yyxx}=4\mathcal{P}(R_{0},\eta)$
$\displaystyle\left[\sin^{2}{\theta}\cos^{2}{\theta}+\left(\frac{R_{0}}{a_{0}}-1\right)\left(2-\frac{R_{0}}{a_{0}}\right)(\sin^{4}{\theta}+\cos^{4}{\theta})\right]$
$\displaystyle\xrightarrow{R_{0}\to
2a_{0}}4\mathcal{P}(2a_{0},\eta)\sin^{2}{\theta}\cos^{2}{\theta},$
$\displaystyle C_{xxxy}=4\mathcal{P}(R_{0},\eta)$
$\displaystyle\left(3-\frac{2R_{0}}{a_{0}}\right)\left[\left(\frac{R_{0}}{a_{0}}-1\right)\sin^{3}{\theta}\cos{\theta}+\left(2-\frac{R_{0}}{a_{0}}\right)\sin{\theta}\cos^{3}{\theta}\right]$
$\displaystyle\xrightarrow{R_{0}\to
2a_{0}}-4\mathcal{P}(2a_{0},\eta)\sin^{3}{\theta}\cos{\theta},$ $\displaystyle
C_{yyyx}=4\mathcal{P}(R_{0},\eta)$
$\displaystyle\left(3-\frac{2R_{0}}{a_{0}}\right)\left[\left(\frac{R_{0}}{a_{0}}-1\right)\sin{\theta}\cos^{3}{\theta}+\left(2-\frac{R_{0}}{a_{0}}\right)\sin^{3}{\theta}\cos{\theta}\right]$
$\displaystyle\xrightarrow{R_{0}\to
2a_{0}}-4\mathcal{P}(2a_{0},\eta)\sin{\theta}\cos^{3}{\theta}.$

The stress correlations shown above for marginally jammed crystals i.e., for
$R_{0}\to 2a_{0}$ have similar angular behavior as that of an isotropic
amorphous material which have been shown earlier in several studies [56, 57,
58, 41, 18]. But for an overcompressed disordered crystal with finite average
pressure, the angular behavior of stress correlations show an additional
anisotropic term which depend on the system overcompression. Therefore for
finite pressure, we can rewrite the above stress correlations of a disordered
crystal as a summation of stress correlations of an isotropic amorphous
material with the addition of a finite shift i.e.,
$C_{\alpha\beta\mu\nu}=C^{t}_{\alpha\beta\mu\nu}+d_{\alpha\beta\mu\nu}$, with

$\displaystyle C_{xxxx}^{t}=$ $\displaystyle 4K_{2D}\sin^{4}{\theta},\hskip
56.9055ptd_{xxxx}(\theta)=4K_{2D}t_{0}\left(t_{0}+2\sin^{2}{\theta}\right),$
(90) $\displaystyle C_{yyyy}^{t}=$ $\displaystyle
4K_{2D}\cos^{4}{\theta},\hskip
56.9055ptd_{yyyy}(\theta)=4K_{2D}t_{0}\left(t_{0}+2\cos^{2}{\theta}\right),$
$\displaystyle C_{xxxy}^{t}=$
$\displaystyle-4K_{2D}\sin^{3}{\theta}\cos{\theta},\hskip
22.76228ptd_{xxxy}(\theta)=-4K_{2D}t_{0}\sin(\theta)\cos(\theta),$
$\displaystyle C_{yyxy}^{t}=$
$\displaystyle-4K_{2D}\sin{\theta}\cos^{3}{\theta},\hskip
22.76228ptd_{yyxy}(\theta)=4K_{2D}t_{0}\sin(\theta)\cos(\theta),$
$\displaystyle C_{xyxy}^{t}=$ $\displaystyle
4K_{2D}\sin^{2}{\theta}\cos^{2}{\theta},\hskip 31.2982ptd_{xyxy}(\theta)=0,$
$\displaystyle C_{yyxx}^{t}=$ $\displaystyle
4K_{2D}\sin^{2}{\theta}\cos^{2}{\theta},\hskip
31.2982ptd_{yyxx}(\theta)=4K_{2D}t_{0}(1+t_{0}),$

where,

$K_{2D}=\mathcal{P}(R_{0},\eta)\left(\frac{2R_{0}}{a_{0}}-3\right)^{2},\hskip
5.69046pt\text{and}\hskip
5.69046ptt_{0}=\frac{2a_{0}-R_{0}}{2R_{0}-3a_{0}}\sim
2\left(1-\left(\frac{\phi_{0}}{\phi}\right)^{1/2}\right).$ (91)

For an isotropic media with average pressure close to zero i.e $R_{0}\to
2a_{0}$ both the theories for stress correlations dealing with two completely
different scenarios give the same results as $d_{\alpha\beta\mu\nu}\to 0$.

##  Appendix C Hammer projection

The Hammer projection is a useful technique for visualizing functions that
have a fixed radial coordinate, meaning they depend only on angular variables.
This projection is particularly valuable in two dimensional visualizations.
Any function with this characteristic can be effectively represented using
Hammer projection. In Hammer projection, we transform from spherical
coordinates $(\theta,\phi)$ to Hammer coordinates $(H_{x},H_{y})$ using the
following equations:

$\displaystyle H_{x}=$
$\displaystyle\frac{2\sqrt{2}\cos{\left(\theta-\pi/2\right)}\sin{\left(\phi/2\right)}}{\sqrt{1+\cos{\left(\theta-\pi/2\right)}\cos{\left(\phi/2\right)}}},\hskip
28.45274ptH_{y}=\frac{\sqrt{2}\sin{\left(\theta-\pi/2\right)}}{\sqrt{1+\cos{\left(\theta-\pi/2\right)}\cos{\left(\phi/2\right)}}}.$
(92)

These equations allow us to map spherical coordinates to Hammer coordinates,
facilitating the visualization of angular-dependent functions in a two
dimensional space.

## References

  * [1] A. Boromand, A. Signoriello, F. Ye, C. S. O’Hern and M. D. Shattuck,  _Jamming of deformable polygons_ ,  Phys. Rev. Lett. 121, 248003 (2018),  10.1103/PhysRevLett.121.248003.
  * [2] C. P. Broedersz, X. Mao, T. C. Lubensky and F. C. MacKintosh,  _Criticality and isostaticity in fibre networks_ ,  Nature Physics 7(12), 983 (2011),  https://doi.org/10.1038/nphys2127.
  * [3] D. Bi, X. Yang, M. C. Marchetti and M. L. Manning,  _Motility-driven glass and jamming transitions in biological tissues_ ,  Phys. Rev. X 6, 021011 (2016),  10.1103/PhysRevX.6.021011.
  * [4] A. J. Licup, S. Münster, A. Sharma, M. Sheinman, L. M. Jawerth, B. Fabry, D. A. Weitz and F. C. MacKintosh,  _Stress controls the mechanics of collagen networks_ ,  Proceedings of the National Academy of Sciences 112(31), 9573 (2015),  10.1073/pnas.1504258112.
  * [5] M. Wyart,  _On the rigidity of amorphous solids_ ,  In _Annales de Physique_ , vol. 30, pp. 1–96. EDP Sciences (2005).
  * [6] P. R. Onck, T. Koeman, T. van Dillen and E. van der Giessen,  _Alternative explanation of stiffening in cross-linked semiflexible networks_ ,  Phys. Rev. Lett. 95, 178102 (2005),  10.1103/PhysRevLett.95.178102.
  * [7] M. Wyart, H. Liang, A. Kabla and L. Mahadevan,  _Elasticity of floppy and stiff random networks_ ,  Phys. Rev. Lett. 101, 215501 (2008),  10.1103/PhysRevLett.101.215501.
  * [8] A. Sharma, A. Licup, K. Jansen, R. Rens, M. Sheinman, G. Koenderink and F. MacKintosh,  _Strain-controlled criticality governs the nonlinear mechanics of fibre networks_ ,  Nature Physics 12(6), 584 (2016),  10.1038/nphys3628.
  * [9] M. F. J. Vermeulen, A. Bose, C. Storm and W. G. Ellenbroek,  _Geometry and the onset of rigidity in a disordered network_ ,  Phys. Rev. E 96, 053003 (2017),  10.1103/PhysRevE.96.053003.
  * [10] A. Baule, F. Morone, H. J. Herrmann and H. A. Makse,  _Edwards statistical mechanics for jammed granular matter_ ,  Rev. Mod. Phys. 90, 015006 (2018),  10.1103/RevModPhys.90.015006.
  * [11] S. Henkes, C. S. O’Hern and B. Chakraborty,  _Entropy and temperature of a static granular assembly: An ab initio approach_ ,  Phys. Rev. Lett. 99, 038002 (2007),  10.1103/PhysRevLett.99.038002.
  * [12] D. Bi, S. Henkes, K. E. Daniels and B. Chakraborty,  _The statistical physics of athermal materials_ ,  Annu. Rev. Condens. Matter Phys. 6(1), 63 (2015),  10.1146/annurev-conmatphys-031214-014336.
  * [13] M. E. Cates, J. P. Wittmer, J.-P. Bouchaud and P. Claudin,  _Jamming, force chains, and fragile matter_ ,  Phys. Rev. Lett. 81, 1841 (1998),  10.1103/PhysRevLett.81.1841.
  * [14] C. S. O’Hern, S. A. Langer, A. J. Liu and S. R. Nagel,  _Random packings of frictionless particles_ ,  Phys. Rev. Lett. 88, 075507 (2002),  10.1103/PhysRevLett.88.075507.
  * [15] H. Yoshino,  _Replica theory of the rigidity of structural glasses_ ,  The Journal of Chemical Physics 136(21), 214108 (2012),  https://doi.org/10.1063/1.4722343.
  * [16] J. Geng, D. Howell, E. Longhi, R. P. Behringer, G. Reydellet, L. Vanel, E. Clément and S. Luding,  _Footprints in sand: The response of a granular material to local perturbations_ ,  Phys. Rev. Lett. 87, 035506 (2001),  10.1103/PhysRevLett.87.035506.
  * [17] H. Vinutha, F. D. Diaz Ruiz, X. Mao, B. Chakraborty and E. Del Gado,  _Stress–stress correlations reveal force chains in gels_ ,  The Journal of chemical physics 158(11) (2023),  https://doi.org/10.1063/5.0131473.
  * [18] J. N. Nampoothiri, M. D’Eon, K. Ramola, B. Chakraborty and S. Bhattacharjee,  _Tensor electromagnetism and emergent elasticity in jammed solids_ ,  Phys. Rev. E 106, 065004 (2022),  10.1103/PhysRevE.106.065004.
  * [19] A. Rida, E. Martinez, D. Rodney and P.-A. Geslin,  _Influence of stress correlations on dislocation glide in random alloys_ ,  Phys. Rev. Mater. 6, 033605 (2022),  10.1103/PhysRevMaterials.6.033605.
  * [20] F. Vogel, A. Zippelius and M. Fuchs,  _Emergence of goldstone excitations in stress correlations of glass-forming colloidal dispersions_ ,  Europhysics Letters 125(6), 68003 (2019),  10.1209/0295-5075/125/68003.
  * [21] D. S. Dagur, C. Mondal and S. Roy,  _Spatial stress correlations in strong colloidal gel systems_ ,  Phys. Rev. B 108, 024106 (2023),  10.1103/PhysRevB.108.024106.
  * [22] J. P. Wittmer, A. N. Semenov and J. Baschnagel,  _Correlations of tensor field components in isotropic systems with an application to stress correlations in elastic bodies_ ,  Phys. Rev. E 108, 015002 (2023),  10.1103/PhysRevE.108.015002.
  * [23] S. Mahajan, J. Chattoraj and M. P. Ciamarra,  _Emergence of linear isotropic elasticity in amorphous and polycrystalline materials_ ,  Phys. Rev. E 103, 052606 (2021),  10.1103/PhysRevE.103.052606.
  * [24] T. Yanagishima, J. Russo and H. Tanaka,  _Common mechanism of thermodynamic and mechanical origin for ageing and crystallization of glasses_ ,  Nature communications 8(1), 15954 (2017),  https://doi.org/10.1038/ncomms15954.
  * [25] L. D. Landau and E. Lifshitz,  _Theoretical physics, vol. 7, theory of elasticity_ ,  Science, Moscow, Main editorial board for physical and mathematical literature (1987),  https://doi.org/10.1063/1.3057037.
  * [26] B. Cui, G. Ruocco and A. Zaccone,  _Theory of elastic constants of athermal amorphous solids with internal stresses_ ,  Granular Matter 21(3), 1 (2019),  https://doi.org/10.1007/s10035-019-0916-4.
  * [27] G. Biroli and P. Urbani,  _Breakdown of elasticity in amorphous solids_ ,  Nature physics 12(12), 1130 (2016),  https://doi.org/10.1038/nphys3845.
  * [28] P. Chaudhuri, S. Karmakar, C. Dasgupta, H. R. Krishnamurthy and A. K. Sood,  _Equilibrium glassy phase in a polydisperse hard-sphere system_ ,  Phys. Rev. Lett. 95, 248301 (2005),  10.1103/PhysRevLett.95.248301.
  * [29] H. Mizuno, S. Mossa and J.-L. Barrat,  _Elastic heterogeneity, vibrational states, and thermal conductivity across an amorphisation transition_ ,  EPL (Europhysics Letters) 104(5), 56001 (2013),  10.1209/0295-5075/104/56001.
  * [30] C. P. Goodrich, A. J. Liu and S. R. Nagel,  _Solids between the mechanical extremes of order and disorder_ ,  Nature Physics 10(8), 578 (2014),  https://doi.org/10.1038/nphys3006.
  * [31] H. Tong, P. Tan and N. Xu,  _From crystals to disordered crystals: A hidden order-disorder transition_ ,  Scientific reports 5, 15378 (2015),  https://doi.org/10.1038/srep15378.
  * [32] P. Acharya, D. Das and K. Ramola,  _Disorder perturbation expansion for athermal crystals_ ,  Phys. Rev. E 104, 034608 (2021),  10.1103/PhysRevE.104.034608.
  * [33] G. Tsekenis,  _Jamming criticality of near-crystals_ ,  EPL (Europhysics Letters) 135(3), 36001 (2021),  10.1209/0295-5075/ac0ffc.
  * [34] P. Charbonneau, E. I. Corwin, L. Fu, G. Tsekenis and M. van der Naald,  _Glassy, gardner-like phenomenology in minimally polydisperse crystalline systems_ ,  Phys. Rev. E 99, 020901 (2019),  10.1103/PhysRevE.99.020901.
  * [35] M. Otto, J.-P. Bouchaud, P. Claudin and J. E. S. Socolar,  _Anisotropy in granular media: Classical elasticity and directed-force chain network_ ,  Phys. Rev. E 67, 031302 (2003),  10.1103/PhysRevE.67.031302.
  * [36] S. Gelin, H. Tanaka and A. Lemaître,  _Anomalous phonon scattering and elastic correlations in amorphous solids_ ,  Nature materials 15(11), 1177 (2016),  https://doi.org/10.1038/nmat4736.
  * [37] N. Xu, V. Vitelli, A. J. Liu and S. R. Nagel,  _Anharmonic and quasi-localized vibrations in jammed solids—modes for mechanical failure_ ,  EPL (Europhysics Letters) 90(5), 56001 (2010),  10.1209/0295-5075/90/56001.
  * [38] M. Shimada, H. Mizuno, M. Wyart and A. Ikeda,  _Spatial structure of quasilocalized vibrations in nearly jammed amorphous solids_ ,  Phys. Rev. E 98, 060901 (2018),  10.1103/PhysRevE.98.060901.
  * [39] E. Lerner and E. Bouchbinder,  _Disordered crystals reveal soft quasilocalized glassy excitations_ ,  Phys. Rev. Lett. 129, 095501 (2022),  10.1103/PhysRevLett.129.095501.
  * [40] W. A. Phillips and A. Anderson,  _Amorphous solids: low-temperature properties_ , vol. 24,  Springer,  https://doi.org/10.1007/978-3-642-81534-8 (1981).
  * [41] J. N. Nampoothiri, Y. Wang, K. Ramola, J. Zhang, S. Bhattacharjee and B. Chakraborty,  _Emergent elasticity in amorphous solids_ ,  Physical review letters 125(11), 118002 (2020),  10.1103/PhysRevLett.125.118002.
  * [42] D. Das, P. Acharya and K. Ramola,  _Long-range correlations in pinned athermal networks_ ,  Phys. Rev. E 104, 014503 (2021),  10.1103/PhysRevE.104.014503.
  * [43] D. Das, P. Acharya and K. Ramola,  _Displacement correlations in disordered athermal networks_ ,  Journal of Statistical Physics 189(2), 1 (2022),  https://doi.org/10.1007/s10955-022-02981-9.
  * [44] I. Goldhirsch and C. Goldenberg,  _On the microscopic foundations of elasticity_ ,  The European Physical Journal E 9(3), 245 (2002),  https://doi.org/10.1140/epje/i2002-10073-5.
  * [45] V. V. Vasisht, P. Chaudhuri and K. Martens,  _Residual stress in athermal soft disordered solids: insights from microscopic and mesoscale models_ ,  Soft Matter 18(34), 6426 (2022),  10.1039/D2SM00615D.
  * [46] A. Barbot, M. Lerbinger, A. Hernandez-Garcia, R. García-García, M. L. Falk, D. Vandembroucq and S. Patinet,  _Local yield stress statistics in model amorphous solids_ ,  Phys. Rev. E 97, 033001 (2018),  10.1103/PhysRevE.97.033001.
  * [47] S. F. Edwards and R. Oakeshott,  _Theory of powders_ ,  Physica A: Statistical Mechanics and its Applications 157(3), 1080 (1989),  https://doi.org/10.1016/0378-4371(89)90034-4.
  * [48] R. Blumenfeld and S. F. Edwards,  _On granular stress statistics: Compactivity, angoricity, and some open issues_ ,  The Journal of Physical Chemistry B 113(12), 3981 (2009),  https://doi.org/10.1021/jp809768y.
  * [49] M. Maier, A. Zippelius and M. Fuchs,  _Emergence of long-ranged stress correlations at the liquid to glass transition_ ,  Phys. Rev. Lett. 119, 265701 (2017),  10.1103/PhysRevLett.119.265701.
  * [50] L. Klochko, J. Baschnagel, J. Wittmer and A. Semenov,  _Long-range stress correlations in viscoelastic and glass-forming fluids_ ,  Soft Matter 14(33), 6835 (2018),  https://doi.org/10.1039/C8SM01055B.
  * [51] H. Tong, S. Sengupta and H. Tanaka,  _Emergent solidity of amorphous materials as a consequence of mechanical self-organisation_ ,  Nature communications 11(1), 1 (2020),  https://doi.org/10.1038/s41467-020-18663-7.
  * [52] M. Schindler,  _A numerical test of stress correlations in fluctuating hydrodynamics_ ,  Chemical Physics 375(2-3), 327 (2010),  https://doi.org/10.1016/j.chemphys.2010.05.008.
  * [53] P. Mora and D. Place,  _Stress correlation function evolution in lattice solid elasto-dynamic models of shear and fracture zones and earthquake prediction_ ,  Earthquake Processes: Physical Modelling, Numerical Simulation and Data Analysis Part II pp. 2413–2427 (2002),  https://doi.org/10.1007/978-3-0348-8197-5_13.
  * [54] T. S. Majmudar and R. P. Behringer,  _Contact force measurements and stress-induced anisotropy in granular materials_ ,  nature 435(7045), 1079 (2005),  https://doi.org/10.1038/nature03805.
  * [55] G. Lois, J. Zhang, T. S. Majmudar, S. Henkes, B. Chakraborty, C. S. O’Hern and R. P. Behringer,  _Stress correlations in granular materials: An entropic formulation_ ,  Phys. Rev. E 80, 060303 (2009),  10.1103/PhysRevE.80.060303.
  * [56] A. Lemaître,  _Structural relaxation is a scale-free process_ ,  Physical review letters 113(24), 245702 (2014),  10.1103/PhysRevLett.113.245702.
  * [57] A. Lemaître,  _Stress correlations in glasses_ ,  The Journal of Chemical Physics 149(10), 104107 (2018),  https://doi.org/10.1063/1.5041461.
  * [58] A. Lemaître, C. Mondal, I. Procaccia and S. Roy,  _Stress correlations in frictional granular media_ ,  Phys. Rev. B 103, 054110 (2021),  10.1103/PhysRevB.103.054110.
  * [59] D. J. Durian,  _Foam mechanics at the bubble scale_ ,  Phys. Rev. Lett. 75, 4780 (1995),  10.1103/PhysRevLett.75.4780.
  * [60] P. Acharya, S. Sengupta, B. Chakraborty and K. Ramola,  _Athermal fluctuations in disordered crystals_ ,  Phys. Rev. Lett. 124, 168004 (2020),  10.1103/PhysRevLett.124.168004.
  * [61] R. Maharana,  _Athermal fluctuations in three dimensional disordered crystals_ ,  Journal of Statistical Mechanics: Theory and Experiment 2022(10), 103201 (2022),  10.1088/1742-5468/ac9466.
  * [62] R. Maharana, J. N. Nampoothiri and K. Ramola,  _First-contact-breaking distributions in strained disordered crystals_ ,  Physical Review E 106(6), 064901 (2022),  10.1103/PhysRevE.106.064901.
  * [63] S. Toxvaerd and J. C. Dyre,  _Communication: Shifted forces in molecular dynamics_ ,  The Journal of Chemical Physics 134(8), 081102 (2011),  10.1063/1.3558787.
  * [64] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler and P. Gumbsch,  _Structural relaxation made simple_ ,  Phys. Rev. Lett. 97, 170201 (2006),  10.1103/PhysRevLett.97.170201.
  * [65] T. Horiguchi,  _Lattice green’s functions for the triangular and honeycomb lattices_ ,  Journal of Mathematical Physics 13(9), 1411 (1972),  https://doi.org/10.1063/1.1666155.
  * [66] H. J. Berendsen, J. v. Postma, W. F. Van Gunsteren, A. DiNola and J. R. Haak,  _Molecular dynamics with coupling to an external bath_ ,  The Journal of chemical physics 81(8), 3684 (1984),  https://doi.org/10.1063/1.448118.
  * [67] T. Belytschko, W. Liu and B. Moran,  _Nonlinear Finite Elements for Continua and Structures_ ,  John Wiley and Sons, Ltd (2000).
  * [68] W. Schirmacher,  _Some comments on fluctuating-elasticity and local oscillator models for anomalous vibrational excitations in glasses_ ,  Journal of Non-Crystalline Solids 357(2), 518 (2011),  https://doi.org/10.1016/j.jnoncrysol.2010.07.052,  6th International Discussion Meeting on Relaxation in Complex Systems.
  * [69] S. Zhang, E. Stanifer, V. V. Vasisht, L. Zhang, E. Del Gado and X. Mao,  _Prestressed elasticity of amorphous solids_ ,  Physical Review Research 4(4), 043181 (2022),  10.1103/PhysRevResearch.4.043181.
  * [70] J. Liu, G. Nian, Q. Cao, S. Qu, X. Wang, D. Zhang and J. Jiang,  _Effect of pre-stress on the onset of yielding in bulk metallic glasses_ ,  Journal of Non-Crystalline Solids 503, 44 (2019),  https://doi.org/10.1016/j.jnoncrysol.2018.09.020.