# Versatility of type-II van der Waals heterostructures: a case study with
SiH-CdCl2

Achintya Priydarshi  Department of Electrical Engineering, Indian Institute of
Technology Kanpur, Kanpur 208016, India  Abhinav Arora  Department of
Materials Science and Engineering, Indian Institute of Technology Kanpur,
Kanpur 208016, India  Yogesh Singh Chauhan<EMAIL_ADDRESS>Department of
Electrical Engineering, Indian Institute of Technology Kanpur, Kanpur 208016,
India  Amit Agarwal<EMAIL_ADDRESS>Department of Physics, Indian Institute
of Technology Kanpur, Kanpur 208016, India  Somnath Bhowmick
<EMAIL_ADDRESS>Department of Materials Science and Engineering, Indian
Institute of Technology Kanpur, Kanpur 208016, India

###### Abstract

Unlike bilayers or a few layers thick materials, heterostructures are designer
materials formed by assembling different monolayers in any desired sequence.
As a result, heterostructures can be tailor-made for specific functionalities
and applications. Here, we show the potential of heterostructures for several
applications, like piezoelectricity, photocatalytic water splitting, and
tunnel field effect transistor (TFET), by taking SiH-CdCl2 as a representative
system. Our study reveals that the characteristics of the heterostructure
mainly depend on the potential difference between the constituent monolayers.
By leveraging the extensive database of layered materials, many such
combinations with a suitable potential difference can be formed for specific
functionalities. This opens up exciting avenues for designing and engineering
materials with tailored properties, paving the way for next-generation device
applications.

##  I INTRODUCTION

Since the experimental discovery of graphene Novoselov  _et al._ (2004),
numerous two-dimensional (2D) materials have been extensively explored due to
their exceptional electrical, mechanical, optical, and catalytic properties.
These include transition metal dichalcogenides, group-V materials like
phosphorene, and group-IV materials like silicene among others Chhowalla  _et
al._ (2013); Zhang  _et al._ (2017); Liu  _et al._ (2014); Sa  _et al._
(2014); Nahas  _et al._ (2017); Mardanya  _et al._ (2016); Di  _et al._
(2017); Oughaddou  _et al._ (2015). To impart additional functionality to 2D
materials beyond their intrinsic properties, various techniques such as
chemical functionalization, surface adsorption, strain engineering, and
quantum confinement are used. Sofo  _et al._ (2007); Rastogi  _et al._ (2014);
Nahas  _et al._ (2016); Puri and Bhowmick (2018); Priydarshi  _et al._ (2018);
Bhowmick  _et al._ (2013). Owing to their tunable properties, van der Waals
(vdW) heterostructures have emerged as a promising class of 2D materials Geim
and Grigorieva (2013); Novoselov  _et al._ (2016). These heterostructures are
are formed by stacking monolayers in desired sequences, allowing for the
tailoring of material properties for specific applications Deng  _et al._
(2016); Bertolazzi  _et al._ (2013).

Heterostructures can be broadly classified into three types based on their
electronic properties: type I, type II, and type III. Among these, type II
heterostructures exhibit distinct advantages due to the presence of valence
and conduction band edges localized at different monolayers Lee  _et al._
(2015). This band alignment is particularly advantageous for applications such
as solar energy harvesting and photocatalytic water splitting. The latter uses
solar energy to generate oxygen and hydrogen, an excellent way to produce
clean energy in the future. For both these applications, efficient charge
separation is essential to minimize electron-hole recombination. Owing to
this, Type-II heterostructures have been proposed to be potential candidates
for these applications Rawat  _et al._ (2019); Mao  _et al._ (2019); Wang  _et
al._ (2013); He  _et al._ (2019); Li  _et al._ (2021); Mao  _et al._ (2023).
Additionally, type II heterostructures have shown promise in other areas,
including tunneling field effect transistors (TFETs) Zhu  _et al._ (2013) and
piezoelectricity Mohanta  _et al._ (2021).

Although previous studies have highlighted various applications of van der
Waals heterostructures individually, a systematic exploration of their
multifunctional aspects is lacking. In this work, we investigate the
properties of the SiH-CdCl2 heterostructure using ab initio calculations.
While the individual monolayers, namely silicane (SiH) and CdCl2, have been
studied separately, their combination as a heterostructure and the resulting
physical properties are relatively unexplored Zhang  _et al._ (2012); Qiu  _et
al._ (2015); Houssa  _et al._ (2011); Restrepo  _et al._ (2014); Sheng  _et
al._ (2020); Han  _et al._ (2020); Luo  _et al._ (2019). We show that SiH-
CdCl2 heterostructure exhibits a type-II band alignment, making it a potential
candidate for efficient photocatalytic water splitting applications. To
confirm its potential as photocatalyst utilizing solar energy, we calculate
its absorption spectrum and show that its absorption peaks lie in the visible
range. Additionally, we explore the application of the heterostructure in
TFETs, and peizoeletricity. Our findings reveal that the characteristics of
the heterostructure depend on the potential difference.

Our paper is organized in the following manner: Section II lists the details
of the ab initio calculations. In Section III, we present the results,
starting with crystal structure in IIIA, followed by electronic and optical
properties in IIIB and IIIC, respectively. We illustrate the charge
distribution at the heterojunction in IIID, and discuss the resulting
piezoelectricity is discussed in IV. We describe other applications like
photocatalysis and TFET in Sections V and VI, respectively. We conclude the
paper in Section VII.

##  II Methodology

For our ab initio calculations, we used the density functional theory (DFT)
and projector-augmented wave pseudo-potentials, as implemented in the Quantum
Espresso package Giannozzi  _et al._ (2009). Exchange and correlation energies
are calculated using the Perdew-Burke-Ernzerhof (PBE) implementation of the
generalized gradient approximation (GGA). We set the energy cutoff for the
plane wave basis to be 60 Ry. A $k$-mesh of $14\times 14\times 1$ is used for
the 2D Brillouin zone (BZ) integrations. For structural relaxations, we use
force and energy convergence threshold of $10^{-3}$ Ry/a.u and $10^{-4}$ Ry,
respectively. We use Grimme’s DFT-D method to capture the effects of the vdW
interactions Grimme (2006). To eliminate the spurious interactions between
neighboring slabs, a vacuum layer of 20 Å is inserted between adjacent layers
along the $z-$direction (perpendicular to the plane of monolayers). VESTA is
used for visualizing crystal structures Momma and Izumi (2008). We use the
Heyd–Scuseria–Ernzerhof (HSE06) hybrid functional to correct for the known
bandgap underestimation within the standard semilocal DFT approachesHeyd  _et
al._ (2003).

Figure 1:  Top view (top panel) and side view (bottom panel) of isolated monolayers of (a) SiH, (b) CdCl2, and (c) the most stable stacking among the six possible stacking of their heterostructure. We define the interface as located at a distance of $d/2$ from the SiH and CdCl2 monolayers. Table 1: Calculated lattice parameter ($a$), band gap (Eg), and band edge location ($E_{\rm CBM}$ and $E_{\rm VBM}$) for isolated SiH, CdCl2 monolayers and their heterostructure. $E_{\rm CBM}$ and $E_{\rm VBM}$ are determined with reference to the vacuum level. The value of Eg within the brackets and outside corresponds to PBE and HSE-based calculations, respectively. Structure |  $a$ (Å) |  $E_{\rm CBM}$ (eV) |  $E_{\rm VBM}$ (eV) |  $E_{g}$ (eV) |
---|---|---|---|---|---
SiH | 3.86 | -2.59 | -5.54 | 2.95 (2.18) |
CdCl2 | 3.80 | -3.39 | -8.73 | 5.34 (3.98) |
SiH/CdCl2 | 3.84 | -3.34 | -5.77 | 2.43 (1.18) |

##  III Physical properties

###  III.1 Crystal structure

Both SiH and CdCl2 monolayers have a hexagonal unit cell [see Figure 1(a) and
(b)], and the optimized lattice parameters are found to be 3.86 Å and 3.80 Å,
respectively. Other structural parameters for both the monolayers we obtained
[see Table S1 of Supporting Information (SI) for further details] also agree
with the reported values Han  _et al._ (2020, 2021); Yue  _et al._ (2021);
Kistanov  _et al._ (2022). This validates our choice of computational
parameters reported in Section II. Out of several possible stacking sequences
of monolayer SiH and CdCl2 [see Figure S1 in SI for further details], the one
with the lowest energy is shown in Figure 1(c). The most stable structure has
a lattice parameter of 3.83 Å. The parent monolayers and the heterostructure
have very similar lattice parameters, with the difference being less than 1%.
This ensures that none of the layers are under significant strain in the
heterostructure. The interlayer separation in this heterostructure (2.14 Å) is
significantly less than that of bilayer graphene (3.4 Å). This indicates some
electrostatic interaction (in addition to the Van der Waals interaction) among
the parent monolayers of the heterostructure supported by charge
redistribution at the interface.

We calculate the binding energy between the monolayers using Almayyali  _et
al._ (2020)

$E_{b}=(E_{SiH/CdCl_{2}}-E_{SiH}-E_{CdCl_{2}})/S_{0},$ (1)

where E${}_{SiH/CdCl_{2}}$, ESiH, and E${}_{CdCl_{2}}$ are the total energy of
the optimized heterostructure and monolayers, respectively. $S_{0}$ is the
area of the optimized heterostructure. We find the binding energy of the most
stable stacking to be -7.44 meV/Å2, which is comparable to that of other vdW
heterostructures like InSe/SiH (-21.96 meV/Å2) and g-C3N4/MoS2 (-17.8 meV/Å2)
Sheng  _et al._ (2020); Wang  _et al._ (2014). A negative binding energy
suggests sufficient binding between the two constituent monolayers of the
heterostructure.

Figure 2: Band structures for (a) SiH monolayer, (b) CdCl2 monolayer, and (c)
SiH/CdCl2 heterostructure. The black (dashed) lines represent the GGA-PBE
results. Red (solid) lines represent the HSE results. (d) The projected
density of states of the SiH/CdCl2 heterostructure and corresponding highest
occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital
(LUMO). The density of states clearly shows that SiH-CdCl2 heterostructure has
a type-II band alignment with the HOMO and LUMO being confined entirely to SiH
and CdCl2 layer, respectively. Figure 3:  (a) The modification in the
electronic band structure and (b) the band gap variation of SiH/CdCl2
heterostructure as a function of vertical electric field Ez.

###  III.2 Electronic band structure and carrier mobility

The electronic band structure is the most fundamental characteristic for
determining the suitability of a material for electronic and optical
applications. We perform the electronic band structure calculations using the
optimized structures described in the previous section. Among several possible
stackings, the one with the lowest energy is chosen for the heterostructure.
The electronic band structures of isolated SiH monolayer, isolated CdCl2
monolayer, and SiH-CdCl2 heterostructure, calculated by using PBE and HSE06
methods, are shown in Figure 2 (a), (b), and (c), respectively. SiH monolayer
is a semiconductor having an indirect band gap of 2.95 eV. The valence band
maximum (VBM) and conduction band minimum (CBM) is at the $\Gamma$ and
M-point, respectively. CdCl2 monolayer is a direct band gap semiconductor with
a band gap of 5.34 eV. Both CBM and VBM are located at the $\Gamma-$point.
Electronic band structures of SiH and CdCl2 monolayers agree with the reported
theoretical and experimental studies Han  _et al._ (2021); Yue  _et al._
(2021); Kistanov  _et al._ (2022). As shown in Figure 2(c), SiH-CdCl2
heterostructure is a direct band gap semiconductor with a band gap of 2.43 eV.
This is lower than the bandgap of both the parent compounds. Both CBM and VBM
are located at the $\Gamma$ point. Projected density of states (PDOS), as well
as HOMO-LUMO plots [integrated local density of states near the band extremum,
as shown in Figure 2(d)] further confirm that the VBM and CBM states are
confined to the SiH and CdCl2 layer, respectively. Such a configuration is
known as the type-II band alignment, which is very effective for electron-hole
separation. We show below that the type-II band alignment makes SiH-CdCl2
heterostructure a suitable candidate for photocatalysis.

Bandstructure tuning via an out-of-plane electric field is a crucial feature
of layered materials. Our calculations show that an out-of-plane electric
field can also engineer the electronic properties of SiH-CdCl2 heterostructure
[Figure 3]. We find that for an electric field applied in the positive
z-direction, the SiH bands move downwards, and CdCl2 bands move upwards,
leading to an overall bandgap increase. For a very high electric field (say
1.0 V/Å), we see a shift from Type-II to Type-I band alignment in this
heterostructure, with both VBM and CBM belonging to the SiH. Such tunability
in band alignment can be exploited in applications like optoelectronic devices
and thin film transistors (TFTs) Gao  _et al._ (2020); Novoselov  _et al._
(2016). We find that other than the extreme case when the band structure
changes from type-II to type-I, for smaller electric field values, the bandgap
changes linearly [Figure 3 (b)] while maintaining the type-II alignment. For
an electric field applied in the negative (positive) z-direction, we find that
the band edges of SiH and CdCl2 move close to (away from) each other to
decrease (increase) the bandgap. Such a linear change can be attributed to the
formation of the heterostructure, as we do not find any change in bandgap in
isolated monolayers of SiH and CdCl2 by varying electric field from -0.5 V/Å
to 0.5 V/Å[see Figure S2 in the SI]. The externally applied field opposes or
favors the internal field (generated by the intrinsic interlayer potential).
As a result, the bandgap changes depending on the direction of the externally
applied electric field. Later we shall discuss in detail how such an electric
field tunable bandgap is suitable for tunnel field effect transistors.

Table 2: Effective mass, mobility and other related parameters of the heterostructure. Carrier | $m_{x}^{*}\textfractionsolidus m_{0}$ | $m_{y}^{*}\textfractionsolidus m_{0}$ | $C_{x}$ | $C_{y}$ | $E_{x}$ | $E_{y}$ |  $\mu_{x}$ (zigzag) |  $\mu_{y}$ (armchair)
---|---|---|---|---|---|---|---|---
|  |  |  (J/m2) |  (J/m2) | (eV) | (eV) | $cm^{2}/V-s$ | $cm^{2}/V-s$
e | 1.04 | 1.04 | 175.453 | 173.925 | 5.925 | 6.711 | 96.7362 | 75.3983
h | 0.497 | 0.497 | 175.453 | 173.925 | 5.082 | 6.358 | 575.599 | 364.642

Having explored the electronic properties, we now focus on carrier mobility,
an important property for transistors and optoelectronic device applications.
At low energies, the electron-acoustic phonon coupling dominates the
scattering process Kaasbjerg  _et al._ (2013). Acoustic phonon-limited carrier
mobility of the negative and positive charge carriers of SiH/CdCl2
heterostructure is calculated using the deformation potential (DP) theory by
Bardeen and Shockley Kaasbjerg  _et al._ (2013). The phonon-limited mobility
of charge carriers in a 2D crystal along the $~{}\alpha$ ($\alpha=x/y$)
direction is given by,

$\mu_{\alpha}=\dfrac{2e\hbar^{3}C_{\alpha}}{3k_{B}Tm^{*}_{\alpha}\sqrt{m^{*}_{\alpha}m^{*}_{\beta}}(E_{\alpha}^{i})^{2}}~{}.$
(2)

Here, $m^{*}_{\alpha}$ and $m^{*}_{\beta}$ represent the effective mass of the
carriers along the transport and the transverse directions, $T$ is the
temperature, $C_{\alpha}$ is the longitudinal acoustic (LA) phonon’s elastic
modulus in the direction of transport, and $E_{\alpha}$ is the DP constant.

The in-plane elastic modulus ($C_{\alpha}$) along the $\alpha$ direction for
2D systems is given by the equation
$(E-E_{0})/S_{0}=\frac{1}{2}C_{\alpha}(\delta l_{\alpha}/l_{\alpha 0})^{2}$,
where $E-E_{0}$ is the change in total energy owing to applied uniaxial strain
$\delta l_{\alpha}/l_{\alpha 0}$, and $S_{0}$ is the area of the pristine unit
cell. $C_{\alpha}$ is calculated by fitting the strain vs. energy curve
(calculated from first principles) with a parabolic function, with the strain
varied between -0.5% and 0.5% [see Figure S3 (a) in SI]. We determine the DP
constant from the strain-induced shift of the band edges (conduction band
minimum for electrons and valence band maximum for holes). The DP constant is
defined as $E_{\alpha}^{i}=\delta E_{i}/(\delta l_{\alpha}/l_{\alpha 0})$,
where $\delta E_{i}$ is the energy change of the $i^{th}$ band (VBM and CBM
for holes and electrons, respectively) under the applied strain of $\delta
l_{\alpha}/l_{\alpha 0}$. A linear fit of $\delta E_{i}$ vs. strain curve
(obtained from first principle calculations) is used to compute
$E_{\alpha}^{i}$ with the strain varying from $-0.5$ to $0.5$ percent [see
Figure S3 (b) in SI]. The effective mass $m^{*}_{\alpha}\equiv
m^{*}_{\alpha\alpha}=\hbar^{2}/(\partial^{2}E(k)/\partial k_{\alpha}\partial
k_{\alpha})$ of the electron (hole) is calculated using a quadratic fit of the
energy dispersion in the vicinity of the minima (maxima) of the conduction
(valance) band.

We calculate the carrier mobility at room temperature, i.e., $T=300$ K [Table
2]. The hole mobilities are almost six times higher than that of electrons.
This is mainly because of a lower hole effective mass than the electron. Our
calculations reveal significant anisotropy, with both electrons and holes
having 1.3 to 1.5 times higher carrier mobility along the zigzag direction.
This arises from the the lower DP constant of electrons and holes in the
zigzag direction than that of the armchair direction. In the next subsection,
we explore the optical properties of the heterostructure.

###  III.3 Optical absorption

Figure 4:  The variation of the total in-plane absorption coefficient with the
incident photon’s energy.

The optical absorption coefficient is calculated using the following
expression, Gajdoš  _et al._ (2006)

$\alpha(\omega)=\frac{\sqrt{2}\omega}{c}\left[\sqrt{\epsilon_{1}^{2}(\omega)+\epsilon_{2}^{2}(\omega)}-\epsilon_{1}(\omega)\right]^{1/2}.$
(3)

Here $\omega$, is the photon frequency,
$\epsilon_{1}=\mathcal{R}[\epsilon(\omega)]$ and
$\epsilon_{2}=\mathcal{I}[\epsilon(\omega)]$ are the real and imaginary
components of the dielectric function, and c is the speed of light. The real
and imaginary parts of the dielectric function are plotted in Figure S4 in SI,
as a function of incident light’s energy. As shown in Figure 4, for the in-
plane polarized light, the first absorption peak is present in the visible
range at $\sim 2.5$ eV for SiH-CdCl2 heterostructure. Following the first
peak, several peaks are in the range $\sim 3.5-4.5$ eV in the near-ultraviolet
and middle-ultraviolet regions. Calculated absorption coefficients ($\sim
10^{6}$ cm-1) are comparable to perovskite solar cells. A significant
absorption coefficient in the visible and near-ultraviolet regions suggests
the potential applicability of SiH-CdCl2 heterostructures in photocatalytic
and photovoltaic devices.

Figure 5: (a) Plane averaged potential along the height of the SiH/CdCl2
heterostructure, showing a potential drop of 2.36 eV across the interface. (b)
Plane averaged charge density difference ($\Delta\rho$), showing a net
accumulation of negative and positive charge in the CdCl2 and SiH side,
respectively. The Vertical dashed red line marks the SiH/CdCl2 heterostructure
interface as defined in Figure 1 (c).

###  III.4 Heterojunction characteristics

Although the constituent monolayers are held together by weak van der Waals
interaction, inter-layer charge transfer is still possible because of the
electronegativity difference between individual layers. The electrostatic
potential is plotted along the $z-$direction [Figure 5(a)], showing a deeper
potential at SiH monolayer, compared to CdCl2. The difference between the
potential minimum in two constituent monolayers is 2.36 eV, resulting in an
internal electric field along negative $z$ (directed from CdCl2 to SiH). As a
result, a net charge transfer is expected from the SiH to the CdCl2 monolayer.
To confirm this, we plot the planar averaged electron density difference
between the heterostructure and its constituent parent compounds. It is
defined as,

$\displaystyle\Delta\rho(z)=\int
dxdy\left[\rho_{SiH/CdCl_{2}}-\rho_{SiH}-\rho_{CdCl_{2}}\right],$ (4)

where $\rho_{SiH/CdCl_{2}}(x,y,z),~{}\rho_{SiH}(x,y,z)$, and
$\rho_{CdCl_{2}}(x,y,z)$ are the charge densities of the SiH/CdCl2
heterostructure and the isolated SiH and CdCl2 monolayers, respectively. As
shown in Figure 5 (b), we find that $\Delta\rho(z)$ is oscillatory, with
alternating regions of charge accumulation (red) and depletion (blue).
However, a net accumulation of negative (positive) charge can be clearly seen
in the CdCl2 (SiH) side, as expected due to the internal electric field.

Because of its unique type-II electronic band structure, optical absorption
peaks in the visible and near-ultraviolet range, and internal electric field,
the SiH-CdCl2 heterostructure has the potential to be used in several device
applications. We discuss some of these applications in the following sections.

##  IV Application: Piezoelectricity

Figure 6: (a) Plane averaged value of charge transfer ($\Delta Q$), and (b)
dipole moment ($\mu$) in the direction normal to the plane of the SiH-CdCl2
heterostructure. (c)-(d) The change in the out-of-plane polarization as a
function of the vertical and biaxial strain, respectively.

Owing to the built-in electric field, which causes a charge redistribution
across the interface [Figure 5(b)], we expect a net electric dipole moment
along the out-of-the-plane direction. With the heterostructure configuration
shown in Fig. 1(a), the charge transfer from SiH to CdCl2 can be estimated as,

$\Delta Q(z)=\int_{-\infty}^{z}\Delta\rho(z^{\prime})\,dz^{\prime},$ (5)

where $\Delta\rho$ is the planar average charge density difference for the
SiH-CdCl2 heterostructure [Figure 5(b)]. The plot of $\Delta Q$ [Figure 6(a)]
confirms that the CdCl2 monolayer gains electrons (negative charges) from the
SiH monolayer. Using $\Delta Q$ values, interface-induced dipole moment $\mu$
can be calculated as,

$\mu(z)=\int_{-\infty}^{z}z^{\prime}\Delta\rho(z^{\prime})\,dz^{\prime}.$ (6)

We find the value of $\mu(z)$ to be $+0.13\times 10^{-4}$ debye [see Figure 6
(b)], induced due to the formation of vdW heterostructure of SiH and CdCl2
monolayers. The dipole points from the negatively charged CdCl2 layer to the
positively charged SiH layer. The existence of an interfacial dipole moment in
the $z$-direction will undoubtedly result in an out-of-plane polarization and
hence an out-of-plane piezoelectric coefficient, which we will discuss next.

The coefficient of piezoelectricity is described as followsChen  _et al._
(2021):

$e_{ijk}=\dfrac{\partial P_{i}}{\partial\epsilon_{jk}}~{}.$ (7)

Here, $P_{i}$ is the polarization, and $\epsilon_{jk}$ represents the strain
tensor. The in-plane $x$, $y$, and out-of-plane $z$ directions are represented
by the values 1, 2, and 3, respectively. Taking polarization along the out-of-
plane direction ($i=3$) and longitudinal strain along the out-of-plane
direction ($j=k=3$), we calculate the out-of-plane piezoelectric coefficient
$e_{333}=e_{33}$ as follows:

$e_{33}=\dfrac{\partial P_{3}}{\partial\epsilon_{3}}\approx\dfrac{\Delta
P_{3}}{\Delta\epsilon_{3}}.$ (8)

Out-of-plane strain leads to a change in interlayer distance, which can be
achieved in two different ways by the application of vertical strain and
biaxial strain. If we take $\Delta\epsilon_{3}\sim\Delta d/d_{0}$, where
$\Delta d$ is the change of interlayer distance from its equilibrium value
$d_{0}$, the change of polarization is much greater in the case of biaxial
strain. As a result, $e_{33}$ is one order of magnitude higher [Figures 6(c)
and 6(d)] in the case of a biaxial strain (-7.75 pC/m) than that of a vertical
strain (-0.76 pC/m). To understand why biaxial strain has a stronger
influence, we compare the change of the intrinsic interlayer potential in the
heterostructure at finite strain [Figure S5 in SI]. Clearly, the modulation is
much stronger in case of biaxial strain. While applying a perpendicular
tensile strain, we only change the vertical distance between the monolayers,
keeping intralayer bonds unchanged. On the other hand, biaxial strain not only
changes the interlayer distance, but also modifies the intralayer bonds in
individual monolayers, leading to a stronger effect on polarization.
Interestingly, we find that the sign of the piezoelectric coefficients is
negative, which makes this type-II heterostructure suitable for many device
applications, such as flexible transducers and energy harvesters Urbanaviciute
_et al._ (2019). Other van der Waals materials like polyvinylidene fluoride,
CuInP2S6, and BiTeX (X = Cl, Br, and I), are also reported to have a negative
piezoelectric response Katsouras  _et al._ (2016); Kim  _et al._ (2019); Liu
_et al._ (2016). Next we explore the photocatalysis potential of SiH-CdCl2
heterostructure.

##  V Application: photocatalysis

Figure 7: Band diagram of isolated monolayers of SiH, CdCl2, and their
heterostructure. The plot also shows the energy levels for oxidation and
reduction potential for reaction related to water splitting with respect to
vacuum level (0 eV). We also mark the valence band offset (VBO) and the
conduction band offset (CBO) as they are responsible for the separation of
charge carriers produced by photons.

Photocatalytic water splitting requires water to be simultaneously oxidized
(2H2O + 4h+ = 4H+ \+ O2) and reduced (4H+ \+ 4e-=2H2), such that the overall
result is 2H2O = 2H2 \+ O2. The holes (h+) and electrons (e-) required for the
process are generated by the absorption of photons. Because of the type-II
band alignment, the oxidation reaction occurs at the heterostructure’s SiH
layer (VBM), while the reduction reaction occurs at the CdCl2 layer (CBM). A
large absorption peak present at $\sim 2.5$ eV [Figure 4] ensures the
availability of abundant photogenerated charge carriers in the presence of
solar light. The physical separation of the oxidation (induced by holes) and
reduction process (induced by electrons) ensures higher efficiency of the
photocatalytic process.

Additionally, the type-II band alignment of SiH-CdCl2 heterostructure has
other advantages for photocatalytic applications. The band diagram of SiH-
CdCl2 heterostructure is illustrated in Figure 7. The band edges are
calculated with respect to the vacuum level. Band edges (VBM and CBM) of
monolayer SiH are at higher energy levels than those of the monolayer CdCl2.
Even in the heterostructure, respective band edges change marginally. There is
a valence band offset (VBO, between VBM of SiH and CdCl2) of 3.1 eV and a
conduction band offset (CBO, between CBM of SiH and CdCl2) of 0.8 eV. The band
gap between the VBM of SiH and CBM of CdCl2 is 2.43 eV. As shown in Figure 7,
because of the CBO and VBO, electrons and holes are accumulated at the CdCl2
and SiH layer, respectively. With charge separation in different layers,
electrons and holes will have a low recombination rate, thus offering more
efficiency for photocatalysis. A significant difference between electron and
hole mobility [see Table 2] further helps reduce the recombination rate and
increase the lifetime of charge carriers. As discussed earlier, owing to low
effective mass, hole mobility is almost six times higher than that of
electrons. As a result, the photogenerated holes move away faster than the
slower electrons, reducing electron-hole recombination and enhancing the
efficiency of SiH-CdCl2 heterostructure as a photocatalyst.

In addition to the properties mentioned above, energy levels of reduction
(H+/H2) and oxidation (H2O/O2) potential should lie within the band gap of the
photocatalyst for water splitting. As shown in Figure 7, the reduction (-4.44
eV) and oxidation (-5.67 eV) potential lie within the bandgap of CdCl2
monolayer, but the VBM of the SiH monolayer is more positive than the (H2O/O2)
potential, making the latter unsuitable as a photocatalyst for water
splitting. However, the reduction and oxidation potential lies within the
bandgap of SiH-CdCl2 heterostructure, making it suitable for photocatalysis.
As already explained, SiH-CdCl2 heterostructure is a better photocatalyst than
monolayer CdCl2 because the latter has lower efficiency due to electron-hole
recombination. The standard oxidation-reduction potentials are a function of
pH [see Eq. S1 and S2 of the SI]. The values reported in Figure 7 are for a
pH=0. We have verified that the oxidation-reduction potentials lie within the
bandgap of SiH-CdCl2 heterostructure from pH 0 to 14 [see Figure S6 in SI],
ensuring a wide range of applicability for photocatalytic application. Having
demonstrated the photocatalytic abilities of the SiH-CdCl2 heterostructure, we
now explore its potential as a field effect transistor.

##  VI Application: tunnel field effect transistor (TFET)

Figure 8:  (a) Schematic depiction of a 2D crystal p-i-n junction based TFET,
(b) Band diagram of TFET, (c) Interband tunneling current density for
SiH/CdCl2 heterostructure at a bias of V=0.3 V and (d) Current density as a
function of voltage.

The ability of a heterojunction based vertical two-dimensional TFET to
overcome the thermionic limit of subthreshold swing (SS) (60 mV/dec) makes it
one of the most promising options for ultra-thin and low-power integrated
circuits. The analytical model proposed by Jena Ma and Jena (2013), suggests
that the interband tunneling current of a vdW-based TFET can be calculated
using

$J_{T}^{2D}=\frac{qg_{s}g_{v}T_{0}k_{B}T}{\left(2\pi\right)^{2}\hbar}\int\limits_{-K_{max}/\eta}^{+K_{max}/\eta}dk_{y}~{}exp\left[-\frac{E_{y}}{\bar{E}}\right]\\\
\times\ln\left\\{\frac{\left(e^{\beta\left(qV-\eta^{2}E_{y}\right)}+e^{-\beta
E_{y}}\right)\left(1+e^{\beta\left(qV-
E_{y}\right)}\right)}{\left(e^{\beta\left(qV-\eta^{2}E_{y}\right)}+e^{\beta\left(qV-
E_{y}\right)}\right)\left(1+e^{-\beta E_{y}}\right)}\right\\}~{}.$ (9)

Here, WKB tunneling probability
$T_{0}=\exp[-4\sqrt{2m_{R}^{*}}E_{g}^{3/2}/3q\hbar F]$, $\bar{E}=q\hbar
F/2\sqrt{2m^{*}_{R}E_{g}}$, $E_{y}=\hbar^{2}K_{y}^{2}/2m_{v}^{*}$, and
$K_{max}^{2}=2m^{*}_{v}qV/\hbar^{2}$. Other terms needed are as follows: $q$
is the charge, $g_{s}$ and $g_{v}$ are spin and valley degeneracy, reduced
mass $m^{*}_{R}=m_{c}^{*}m^{*}_{v}/(m_{c}^{*}+m^{*}_{v})$, $m_{c}^{*}$ and
$m^{*}_{v}$ are the effective masses of the valence and conduction band
extremum, and $F$ is the vertical field. Further, $\eta$ and $\beta$ are
defined as $\eta^{2}=m^{*}_{v}/m^{*}_{R}$ and $\beta=1/k_{B}T$, where $\hbar$
and $k_{B}$ are the reduced Plank constant and Boltzmann constant,
respectively. For SiH/CdCl2 heterostructure, the electronic band gap $E_{g}$
changes linearly with the applied field [see Figure 3(b)]. We can express it
as $E_{g}=E_{g0}\pm aF$, where $E_{g0}$ is the zero-field band gap of the
SiH/CdCl2 heterostructure. The constant $a=1.56$ is calculated from a linear
fit of the $E_{g}$ vs. field plot in Figure 3 (b). We also verify that
linearity holds good for very small electric field values [see Figure S7 in
SI].

The calculated tunneling current density as a function of vertical field and
voltage (at different temperatures) are shown in Figures 8(c) and 8 (d),
respectively. We find that the interband tunneling happens only when $V$
exceeds zero since the zero drain voltage $V=0$ is defined when the source’s
conduction band edge is aligned with the drain’s valence band edge. Due to the
exponential dependence of the tunneling probability on the field $F$
($T_{0}\propto exp[-4\sqrt{2m_{R}^{*}}E_{g}^{3/2}/3q\hbar F]$), the tunneling
current varies exponentially for low values of F [Figure 8(c)], and begins to
saturate when $F$ exceeds 1.5 V/nm. The tunneling transistors based on 2D
transition metal dichalcogenides (TMDs) show a similar trend Ma and Jena
(2013); Jiang and Li (2014); Balaji  _et al._ (2018); Pal  _et al._ (2020).

The change of bandgap ($E_{g}$) with external field [Figure 3] also affects
the current, as the tunneling probability increases (decreases) due to a
decrease (increase) in $E_{g}$. We carry out a systematic investigation by
calculating the current in three different cases: I) $E_{g}$ decreases with
$F$, II) $E_{g}$ remains unchanged, and III) $E_{g}$ increases with $F$. While
case I and III are observed in SiH/CdCl2 [Figure 3], case II is calculated to
fully understand the extent of the effect of bandgap tuning via an external
electric field. As expected, the lower the bandgap, the higher the tunneling
current, which is the highest in case I, followed by case II, and lowest in
case III. As shown in Figure 8(c), the on-state current = 1300$~{}\mu A/\mu m$
for case I is about 4 and 18 times greater than that in case II and case III,
respectively. Note that case I also has a better SS [Figure 8(c)] than the
rest of the cases. In all cases, SS is less than 60 mV/dec, with the lowest SS
in case I.

The tunneling current as a function of voltage at the room temperature is
presented in Figure 8(d). In contrast to typical metal-oxide-semiconductor
field-effect transistors (MOSFETs), we observe a rectifying curvature. Several
studies conducted on TFETs have demonstrated a similar trend Mallik and
Chattopadhyay (2011); Pal  _et al._ (2011); Lattanzio  _et al._ (2011). Such a
quasi-exponential behavior of current characteristics is typical of band-to-
band tunneling (BTBT) and is due to the drain barrier modulation of the BTBT
window. As depicted in Figure 8(d), at small reverse bias voltages $V\ll
V_{0}=\eta^{2}\bar{E}/q$, the tunneling current changes as $J^{2D}_{T}\sim
V^{1.5}$ Ma and Jena (2013). For larger values of $V\gg V_{0}$, Eq. 9 changes
to $J^{2D}_{T}=G\times V$, where $G\approx[q^{2}/h(g_{s}g_{v}/2\pi)\sqrt{(2\pi
m^{*}_{v}/\hbar^{2})\times(q\hbar F/\sqrt{8m^{*}_{R}E_{g}})}\times T_{0}]$ and
the current thus varies linearly with V Ma and Jena (2013).

##  VII Conclusion

Using ab initio calculations, we demonstrate the potential of a newly
predicted heterostructure formed by SiH and CdCl2 monolayers for different
device applications. The heterostructure has a direct bandgap of 2.43 eV and a
Type-II band alignment in which the VBM and CBM are localized to the SiH and
CdCl2 monolayer, respectively. Such a band alignment is beneficial for
photocatalytic water splitting, as electrons and holes are separated at two
different monolayers, reducing their recombination rate and increasing device
efficiency. We show that the optical absorption spectrum has peaks in the
visible wavelength, confirming photocatalytic device operation using solar
energy. Our calculations further reveal that the heterostructure can be
explored for making TFET devices.

Our study highlights that the characteristics of the heterostructure mainly
depend on the potential difference between the constituent monolayers, leading
to the charge transfer, which is the origin of piezoelectricity. With many
theoretically predicted and experimentally realized layered materials, we can
form many combinations of heterostructures with a suitable potential
difference to engineer specific functionalities. Looking at the potential for
various functionality, we believe that an exciting future lies ahead for
device applications in assembled heterostructures.

## ACKNOWLEDGMENTS

We acknowledge National Super Computing Mission (NSM) for providing computing
resources of “PARAM Sanganak” at IIT Kanpur, which is implemented by C-DAC and
supported by the Ministry of Electronics and Information Technology (MeitY)
and Department of Science and Technology (DST), Government of India. We also
acknowledge the HPC facility provided by CC, IIT Kanpur. We acknowledge the
Science and Engineering Research Board (SERB, file no: EMR/2017/004970) for
financial support.

## References

  * Novoselov  _et al._ (2004) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2004)
  * Chhowalla  _et al._ (2013) M. Chhowalla, H. S. Shin, G. Eda, L.-J. Li, K. P. Loh, and H. Zhang, Nature chemistry 5, 263 (2013)
  * Zhang  _et al._ (2017) K. Zhang, Y. Feng, F. Wang, Z. Yang, and J. Wang, Journal of Materials Chemistry C 5, 11992 (2017)
  * Liu  _et al._ (2014) H. Liu, A. T. Neal, Z. Zhu, Z. Luo, X. Xu, D. Tománek, and P. D. Ye, ACS Nano 8, 4033 (2014), pMID: 24655084
  * Sa  _et al._ (2014) B. Sa, Y.-L. Li, J. Qi, R. Ahuja, and Z. Sun, The Journal of Physical Chemistry C 118, 26560 (2014)
  * Nahas  _et al._ (2017) S. Nahas, A. Bajaj, and S. Bhowmick, Phys. Chem. Chem. Phys. 19, 11282 (2017)
  * Mardanya  _et al._ (2016) S. Mardanya, V. K. Thakur, S. Bhowmick, and A. Agarwal, Phys. Rev. B 94, 035423 (2016)
  * Di  _et al._ (2017) J. Di, J. Xia, H. Li, and Z. Liu, Nano Energy 35, 79 (2017)
  * Oughaddou  _et al._ (2015) H. Oughaddou, H. Enriquez, M. R. Tchalala, H. Yildirim, A. J. Mayne, A. Bendounan, G. Dujardin, M. Ait Ali, and A. Kara, Progress in Surface Science 90, 46 (2015), special Issue on Silicene
  * Sofo  _et al._ (2007) J. O. Sofo, A. S. Chaudhari, and G. D. Barber, Phys. Rev. B 75, 153401 (2007)
  * Rastogi  _et al._ (2014) P. Rastogi, S. Kumar, S. Bhowmick, A. Agarwal, and Y. S. Chauhan, The Journal of Physical Chemistry C 118, 30309 (2014)
  * Nahas  _et al._ (2016) S. Nahas, B. Ghosh, S. Bhowmick, and A. Agarwal, Phys. Rev. B 93, 165413 (2016)
  * Puri and Bhowmick (2018) S. Puri and S. Bhowmick, Phys. Rev. Mater. 2, 044001 (2018)
  * Priydarshi  _et al._ (2018) A. Priydarshi, Y. S. Chauhan, S. Bhowmick, and A. Agarwal, Phys. Rev. B 97, 115434 (2018)
  * Bhowmick  _et al._ (2013) S. Bhowmick, A. Medhi, and V. B. Shenoy, Phys. Rev. B 87, 085412 (2013)
  * Geim and Grigorieva (2013) A. K. Geim and I. V. Grigorieva, Nature 499, 419 (2013)
  * Novoselov  _et al._ (2016) K. S. Novoselov, A. Mishchenko, A. Carvalho, and A. H. Castro Neto, Science 353, 461 (2016)
  * Deng  _et al._ (2016) D. Deng, K. Novoselov, Q. Fu, N. Zheng, Z. Tian, and X. Bao, Nature nanotechnology 11, 218 (2016)
  * Bertolazzi  _et al._ (2013) S. Bertolazzi, D. Krasnozhon, and A. Kis, ACS nano 7, 3246 (2013)
  * Lee  _et al._ (2015) Y. Lee, Y. Hwang, and Y.-C. Chung, ACS Applied Materials & Interfaces 7, 7163 (2015)
  * Rawat  _et al._ (2019) A. Rawat, R. Ahammed, Dimple, N. , M. K. Mohanta, and A. De Sarkar, The Journal of Physical Chemistry C 123, 12666 (2019)
  * Mao  _et al._ (2019) Y. Mao, C. Xu, J. Yuan, and H. Zhao, Journal of Materials Chemistry A 7, 11265 (2019)
  * Wang  _et al._ (2013) Y. Wang, Q. Wang, X. Zhan, F. Wang, M. Safdar, and J. He, Nanoscale 5, 8326 (2013)
  * He  _et al._ (2019) C. He, J. Zhang, W. Zhang, and T. Li, The journal of physical chemistry letters 10, 3122 (2019)
  * Li  _et al._ (2021) X.-H. Li, B.-J. Wang, G.-D. Wang, X.-F. Yang, R.-Q. Zhao, X.-T. Jia, and S.-H. Ke, Sustainable Energy & Fuels 5, 2249 (2021)
  * Mao  _et al._ (2023) Y. Mao, C. Qin, X. Zhou, Z. Zhang, and J. Yuan, Applied Surface Science 617, 156298 (2023)
  * Zhu  _et al._ (2013) Y. Zhu, N. Jain, D. K. Mohata, S. Datta, D. Lubyshev, J. M. Fastenau, A. K. Liu, and M. K. Hudait, Journal of Applied Physics 113, 024319 (2013)
  * Mohanta  _et al._ (2021) M. K. Mohanta, A. Rawat, and A. De Sarkar, Applied Surface Science 568, 150928 (2021)
  * Zhang  _et al._ (2012) P. Zhang, X. Li, C. Hu, S. Wu, and Z. Zhu, Physics Letters A 376, 1230 (2012)
  * Qiu  _et al._ (2015) J. Qiu, H. Fu, Y. Xu, Q. Zhou, S. Meng, H. Li, L. Chen, and K. Wu, ACS nano 9, 11192 (2015)
  * Houssa  _et al._ (2011) M. Houssa, E. Scalise, K. Sankaran, G. Pourtois, V. Afanas’ Ev, and A. Stesmans, Applied Physics Letters 98, 223107 (2011)
  * Restrepo  _et al._ (2014) O. D. Restrepo, R. Mishra, J. E. Goldberger, and W. Windl, Journal of Applied Physics 115, 033711 (2014)
  * Sheng  _et al._ (2020) W. Sheng, Y. Xu, M. Liu, G. Nie, J. Wang, and S. Gong, Physical Chemistry Chemical Physics 22, 21436 (2020)
  * Han  _et al._ (2020) S. Han, Y. Li, and Z. Wang, Physical Chemistry Chemical Physics 22, 17145 (2020)
  * Luo  _et al._ (2019) Y. Luo, Y. Hu, and Y. Xie, Journal of Materials Chemistry A 7, 27503 (2019)
  * Giannozzi  _et al._ (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, Journal of Physics: Condensed Matter 21, 395502 (2009)
  * Grimme (2006) S. Grimme, Journal of computational chemistry 27, 1787 (2006)
  * Momma and Izumi (2008) K. Momma and F. Izumi, Journal of Applied crystallography 41, 653 (2008)
  * Heyd  _et al._ (2003) J. Heyd, G. E. Scuseria, and M. Ernzerhof, The Journal of Chemical Physics 118, 8207 (2003)
  * Han  _et al._ (2021) S. Han, Y. Li, and Z. Wang, Physica E: Low-dimensional Systems and Nanostructures 134, 114869 (2021)
  * Yue  _et al._ (2021) L. Yue, W. Jin, A.-J. Mao, and X. Kuang, The Journal of Physical Chemistry C 125, 12983 (2021)
  * Kistanov  _et al._ (2022) A. A. Kistanov, S. A. Shcherbinin, R. Botella, A. Davletshin, and W. Cao, The journal of physical chemistry letters 13, 2165 (2022)
  * Almayyali  _et al._ (2020) A. O. M. Almayyali, B. B. Kadhim, and H. R. Jappor, Chemical Physics 532, 110679 (2020)
  * Wang  _et al._ (2014) J. Wang, Z. Guan, J. Huang, Q. Li, and J. Yang, Journal of Materials Chemistry A 2, 7960 (2014)
  * Gao  _et al._ (2020) R. Gao, H. Liu, J. Yang, F. Yang, T. Wang, Z. Zhang, X. Liu, H. Jia, B. Xu, and H. Ma, Applied Surface Science 529, 147026 (2020)
  * Kaasbjerg  _et al._ (2013) K. Kaasbjerg, K. S. Thygesen, and A.-P. Jauho, Phys. Rev. B 87, 235312 (2013)
  * Gajdoš  _et al._ (2006) M. Gajdoš, K. Hummer, G. Kresse, J. Furthmüller, and F. Bechstedt, Phys. Rev. B 73, 045112 (2006)
  * Chen  _et al._ (2021) Y. Chen, Z. Tang, H. Shan, B. Jiang, Y. Ding, X. Luo, and Y. Zheng, Physical Review B 104, 075449 (2021)
  * Urbanaviciute  _et al._ (2019) I. Urbanaviciute, X. Meng, M. Biler, Y. Wei, T. D. Cornelissen, S. Bhattacharjee, M. Linares, and M. Kemerink, Materials Horizons 6, 1688 (2019)
  * Katsouras  _et al._ (2016) I. Katsouras, K. Asadi, M. Li, T. B. Van Driel, K. S. Kjaer, D. Zhao, T. Lenz, Y. Gu, P. W. Blom, D. Damjanovic, _et al._ , Nature materials 15, 78 (2016)
  * Kim  _et al._ (2019) J. Kim, K. M. Rabe, and D. Vanderbilt, Phys. Rev. B 100, 104115 (2019)
  * Liu  _et al._ (2016) F. Liu, L. You, K. L. Seyler, X. Li, P. Yu, J. Lin, X. Wang, J. Zhou, H. Wang, H. He, _et al._ , Nature communications 7, 1 (2016)
  * Ma and Jena (2013) N. Ma and D. Jena, Applied Physics Letters 102, 132102 (2013)
  * Jiang and Li (2014) X.-W. Jiang and S.-S. Li, Applied Physics Letters 104, 193510 (2014)
  * Balaji  _et al._ (2018) Y. Balaji, Q. Smets, C. J. L. De La Rosa, A. K. A. Lu, D. Chiappe, T. Agarwal, D. H. Lin, C. Huyghebaert, I. Radu, D. Mocuta, _et al._ , IEEE Journal of the Electron Devices Society 6, 1048 (2018)
  * Pal  _et al._ (2020) A. Pal, W. Cao, and K. Banerjee, IEEE Transactions on Electron Devices 67, 4473 (2020)
  * Mallik and Chattopadhyay (2011) A. Mallik and A. Chattopadhyay, IEEE transactions on electron devices 58, 4250 (2011)
  * Pal  _et al._ (2011) A. Pal, A. B. Sachid, H. Gossner, and V. R. Rao, IEEE Transactions on Electron Devices 58, 1045 (2011)
  * Lattanzio  _et al._ (2011) L. Lattanzio, L. De Michielis, and A. M. Ionescu, IEEE Electron Device Letters 33, 167 (2011)


# Odd Dipole Screening in Radial Inflation

Yang Fu  Beijing National Laboratory for Condensed Matter Physics and
Laboratory of Soft Matter Physics, Institute of Physics, Chinese Academy of
Sciences, Beijing 100190, China  H. George E. Hentschel  Department of
Physics, Emory University, Atlanta, Georgia-30322, USA  Pawandeep Kaur
Department of Chemical & Biological Physics, Weizmann Institute of Science,
Rehovot-7610001, Israel  Avanish Kumar  Department of Chemical & Biological
Physics, Weizmann Institute of Science, Rehovot-7610001, Israel  Itamar
Procaccia  Department of Chemical & Biological Physics, Weizmann Institute of
Science, Rehovot-7610001, Israel  Sino-Europe Complex Science Center, School
of Mathematics, North University of China, Shanxi, Taiyuan 030051, China.

(August 28, 2024)

###### Abstract

The inflation of an inner radial (or spherical) cavity in an amorphous solids
confined in a disk (or a sphere), served as a fruitful case model for studying
the effects of plastic deformations on the mechanical response. It was shown
that when the field associated with Eshelby quadrupolar charges is non-
uniform, the displacement field is riddled with dipole charges that screen
elasticity, reminiscent of Debye monopoles screening in electrostatics. In
this paper we look deeper into the screening phenomenon, taking into account
the consequences of irreversibility that are associated with the breaking of
Chiral symmetry. We consider the equations for the displacement field with the
presence of “Odd Dipole Screening", solve them analytically and compare with
numerical simulations. Suggestions how to test the theory in experiments are
provided.

##  I Introduction

In a series of recent papers the phenomenon of screening in the mechanical
response of amorphous solids to non-uniform stresses or strains was examined
theoretically[1, 2, 3, 4, 5], experimentally [6, 7] and simulationally [8, 9,
10]. The phenomenon at question deals with an amorphous solid that is at
mechanical equilibrium (i.e. with vanishing net force on each constituent
particle), that is subjected to strain, leading to loss of mechanical
equilibrium. Upon relaxation back to mechanical equilibrium, the positions of
constituent particles change, until a new configuration is attained in which
the net force on each particle vanishes. We are interested in the displacement
field ${\bm{d}}({\bm{r}})$ which is the difference between the final and
initial in positions of the centers of mass of each particle.

For a classical elastic solid the equation obeyed by displacement field
${\bm{d}}$ reads [11]:

$\Delta{\bm{d}}+(1+\tilde{\lambda}){\bm{\nabla}}({\bm{\nabla}}\cdot{\bm{d}})=0,\quad\tilde{\lambda}\equiv\lambda/\mu\
.$ (1)

Here $\lambda$ and $\mu$ are the classical Lame’ coefficients. For an
amorphous solids things are different. First, generically any strain on
amorphous solids results in plastic deformation [12, 13]. Typical plastic
events are quadrupolar in nature [14, 15], called sometime “Eshelby
inclusions" due to the resemblance to the Eshelby solution of the displacement
field caused by forcing a circle to change to an ellipse within an elastic
material [16]. When the density of quadrupolar events is finite, one refers to
the resulting quadrupolar field as $Q^{\alpha\beta}({\bm{r}})$. It was shown
the creation of such a field leads to a renormalization of the elastic moduli
[1, 2]. Further, when the quadrupolar field of such events is non-uniform, it
was explained that gradients of the field act as effective dipoles,

${\mathcal{P}}^{\alpha}\equiv\partial_{\beta}Q^{\alpha\beta}\ .$ (2)

When these are present, the equation satisfied by the displacement field
change drastically, reading [1, 2]

$\Delta{\bm{d}}+(1+\tilde{\lambda}){\bm{\nabla}}({\bm{\nabla}}\cdot{\bm{d}})+{\bm{\Gamma}}{\bm{d}}=0.$
(3)

Here ${\bm{\Gamma}}$ is a tensor that needs to be specified, containing
screening parameters. In previous work it was assumed that in isotropic and
homogeneous amorphous solids one can assume that $\Gamma$ is diagonal, leading
to an equation to be solved of the form

$\Delta{\bm{d}}+(1+\tilde{\lambda}){\bm{\nabla}}({\bm{\nabla}}\cdot{\bm{d}})+k^{2}{\bm{d}}=0.$
(4)

The term $k^{2}{\bm{d}}$ is responsible for translational symmetry breaking,
the introduction of a typical length scale $\ell$, $\ell\sim k^{-1}$, and to
screening phenomena that change dramatically the expected displacement field
${\bm{d}}$ from the predictions of Eq. (1). While experiments and simulations
provided support for the predictions of Eq. (4), there was one consequence
that was found to be in only fair agreement. This is a resulting constitutive
equation that relates the dipole and the displacement fields

${\bm{{\mathcal{P}}}}=-\kappa^{2}{\bm{d}}\ .$ (5)

This constitutive relations means that the dipole and the displacement fields
are expected to be co-linear and opposite in direction. This prediction was
not accurately corroborated in simulations, an angle was distinctly existing
between the prediction and the measured direction of the dipole field.

The resolution of this issue was offered in Ref. [7]. It was argued that in
systems that do not conserve energy along a closed loop of straining, the
assumption of diagonal tensor ${\bm{\Gamma}}$ is not tenable. Rather, an odd
tensor is expected, of the form

${\bm{\Gamma}}=\begin{bmatrix}\kappa_{e}^{2}&\mp\kappa_{o}^{2}\\\
\pm\kappa_{o}^{2}&\kappa_{e}^{2}\end{bmatrix}\ .$ (6)

The $\mp$ sign of $\kappa_{0}^{2}$ means that it can be negative or positive,
with an anti-symmetric counterpart $\pm\kappa_{0}^{2}$. The presence of this
odd tensor in the theory has led to the term “odd dipole screening". The
existence of the anti-symmetric tensor ${\bm{\Gamma}}$ leads to breaking
Chiral symmetry. Even for purely radial inflation as discussed below, the
resulting displacement field exhibits non-trivial transverse (angular)
response, as will be shown explicitly in the next Section. Denoting the angle
formed between the (inverse) direction of the displacement field and the
dipole field as $\Phi$, one predicts that

$\tan\Phi=\pm\frac{\kappa_{o}^{2}}{\kappa_{e}^{2}}\ .$ (7)

The angle $\Phi$ can be positive or negative, depending on the sign of
$\kappa_{o}^{2}$, but is expected to be the same for all the angles between
the dipole and displacement fields in a given realization. Tests of this
prediction are offered below.

The structure of this paper is as follows: in Sect. II we discuss the
equations to be solved for two-dimensional radial inflation of an inner disk,
and provide their analytic solutions. We then discuss typical numerical
simulations using classical glass formers, and compare the numerical
measurements to the analytic solutions. A good agreement is reported. In Sect.
III we turn to a deeper analysis of the consequences of Chiral symmetry
breaking. In particular we want to verify the prediction embodied in Eq. (7).
To this aim we will start from the measured displacement field, extract the
quadrupolar field associated with the displacement, and from it the dipole
field. Having done so we can examine the angle between the dipole field and
the inverse direction of the displacement field. The analysis corroborates Eq.
(7) with the actual measured values of the screening parameters $\kappa_{e}$
and $\kappa_{o}$. Section IV offers a summary and discussion of the paper,
including suggestions how to test the theory in experiments.

Figure 1: Schematic diagram of the geometry under discussion. Material is
confined in an annulus between an outer disk of radius $r_{\rm out}$ and an
inner boundary of radius $r_{\rm in}$. The inner radius is then inflated to
$r_{\rm in}+d_{0}$ and then the material is relaxed back to mechanical
equilibrium

##  II Translational and Chiral Symmetry Breaking

As always, the solution of differential equations depend on the boundary
conditions. To make the novel points very clear, we will focus in this paper
on two dimensions, having in mind an amorphous solids confined between an
outer disk of radius $r_{\rm out}$ and an inner boundary of radius $r_{\rm
in}$, cf. Fig. 1.

After equilibration, the inner circle is inflated, $r_{\rm in}\to r_{\rm
in}+d_{0}$, and the system is equilibrated again. As said, the displacement
field studied is the difference between the two configurations after and
before the inflation.

###  II.1 Solving for the displacement field

To solve Eq. (3) with ${\bm{\Gamma}}$ as defined in Eq. (6), the displacement
field ${\bm{d}}$ can be separated into radial and transverse components

${\bm{d}}=d_{r}(r,\theta)\hat{r}+d_{\theta}(r,\theta)\hat{\theta}.$ (8)

For concreteness we will present the solution for the case of $\kappa_{o}^{2}$
positive, and a parallel analysis can be easily done for a negative
$\kappa_{o}^{2}$. We then evaluate the screening term as

${\bm{\Gamma}}{\bm{d}}=\begin{bmatrix}(\kappa_{e}^{2}d_{r}-\kappa_{o}^{2}d_{\theta})\leavevmode\nobreak\
\hat{r}\\\ (\kappa_{o}^{2}d_{r}+\kappa_{e}^{2}d_{\theta})\leavevmode\nobreak\
\hat{\theta}.\end{bmatrix}\ .$ (9)

Now equation (3) can be decomposed into following coupled differential
equations in $r,$ and $\theta$ :

$\frac{(\tilde{\lambda}+2)}{r^{2}}\left[r^{2}d_{r}^{\prime\prime}+rd_{r}^{\prime}-d_{r}\right]-\frac{(\tilde{\lambda}+3)}{r^{2}}\frac{\partial
d_{\theta}}{\partial\theta}+\frac{1}{r^{2}}\frac{\partial^{2}d_{r}}{\partial\theta^{2}}+\frac{(\tilde{\lambda}+1)}{r}\frac{\partial^{2}d_{\theta}}{\partial
r\partial\theta}+\kappa_{e}^{2}d_{r}-\kappa_{o}^{2}d_{\theta}=0,$ (10)
$\frac{1}{r^{2}}\left[r^{2}d_{\theta}^{\prime\prime}+rd_{\theta}^{\prime}+\frac{\partial^{2}d_{\theta}}{\partial\theta^{2}}-d_{\theta}\right]+\frac{(1+\tilde{\lambda})}{r^{2}}\left[\frac{\partial^{2}d_{\theta}}{\partial\theta^{2}}+r\frac{\partial^{2}d_{r}}{\partial
r\partial\theta}\right]+\frac{(\tilde{\lambda}+3)}{r^{2}}\frac{\partial
d_{r}}{\partial\theta}+\kappa_{o}^{2}d_{r}+\kappa_{e}^{2}d_{\theta}=0.$ (11)

Figure 2: Panel (a) and (b): Typical solutions of $a_{0}(r)$ and $b_{0}(r)$
for a fixed $\kappa_{\rm o}$=0.041 and for different values of $\kappa_{\rm
e}=0.0204,0.0612,0.102,0.1429$. Panel (c): Comparison of $a_{0}$ and $b_{0}$
for $\kappa_{\rm o}=0.041$ and $\kappa_{\rm e}=0.204$.

Figure 3: Panel (a) and (b): Typical solutions of $a_{0}(r)$ and $b_{0}(r)$
for a fixed $\kappa_{\rm e}$=0.041 and for different values of $\kappa_{\rm
o}=0.0204,0.0612,0.102,0.1429$. Panel (c): Comparison of $a_{0}$ and $b_{0}$
for $\kappa_{\rm o}=0.204$ and $\kappa_{\rm e}=0.041$.

Taking into account the periodic boundary conditions on the angle $\theta$ we
can now seek a solution of these equations in the form of Fourier series for
$d_{r}(r,\theta)$ and $d_{\theta}(r,\theta)$ respectively:

$\displaystyle
d_{r}(r,\theta)=a_{0}(r)+\sum_{n=1}^{\infty}\left[a_{n}(r)\cos(n\theta)+b_{n}(r)\sin(n\theta)\right],$
$\displaystyle
d_{\theta}(r,\theta)=b_{0}(r)+\sum_{n=1}^{\infty}\left[c_{n}(r)\cos(n\theta)+d_{n}(r)\sin(n\theta)\right].$
(12)

Of course, the more coefficients we keep, the more equations we need to solve.
However, due to the orthogonality of the Fourier coefficients and the
linearity of the equations, different order coefficients do not mix. Thus, for
simplicity, to demonstrate the decoupling, we consider only the first-order
Fourier terms, $n=1$. We have

$\displaystyle
d_{r}(r,\theta)=a_{0}(r)+a_{1}(r)\cos(\theta)+b_{1}(r)\sin(\theta),$
$\displaystyle
d_{\theta}(r,\theta)=b_{0}(r)+c_{1}(r)\cos(\theta)+d_{1}(r)\sin(\theta).$ (13)

After substitution of this ansatz, and after some simplifications, Eqs. (11)
and (II.1) respectively take the following forms:

$\displaystyle\left[r^{2}a_{0}^{\prime\prime}(r)+ra_{0}^{\prime}(r)-a_{0}(r)\right]+\frac{\kappa_{e}^{2}r^{2}}{(\tilde{\lambda}+2)}a_{0}(r)-\frac{\kappa_{0}^{2}r^{2}}{(\tilde{\lambda}+2)}b_{0}(r)$
$\displaystyle+\left[\left[r^{2}a_{1}^{\prime\prime}(r)+ra_{1}^{\prime}(r)\right]+\frac{(\kappa_{e}^{2}r^{2}-\tilde{\lambda}-3)}{(\tilde{\lambda}+2)}a_{1}(r)-\frac{\kappa_{o}^{2}r^{2}}{(\tilde{\lambda}+2)}c_{1}(r)+\frac{(\tilde{\lambda}+1)r}{(\tilde{\lambda}+2)}d_{1}^{\prime}(r)-\frac{(\tilde{\lambda}+3)}{(\tilde{\lambda}+2)}d_{1}(r)\right]\cos(\theta)$
$\displaystyle+\left[\left[r^{2}b_{1}^{\prime\prime}(r)+rb_{1}^{\prime}(r)\right]+\frac{(\kappa_{e}^{2}r^{2}-\tilde{\lambda}-3)}{(\tilde{\lambda}+2)}b_{1}(r)-\frac{\kappa_{o}^{2}r^{2}}{(\tilde{\lambda}+2)}d_{1}(r)-\frac{(\tilde{\lambda}+1)r}{(\tilde{\lambda}+2)}c_{1}^{\prime}(r)+\frac{(\tilde{\lambda}+3)}{(\tilde{\lambda}+2)}c_{1}(r)\right]\sin(\theta)=0,$
(14)
$\displaystyle\left[r^{2}b_{0}^{\prime\prime}(r)+rb_{0}^{\prime}(r)-b_{0}(r)\right]+\kappa_{e}^{2}r^{2}b_{0}(r)+\kappa_{0}^{2}r^{2}a_{0}(r)$
$\displaystyle+\left[\left[r^{2}c_{1}^{\prime\prime}(r)+rc_{1}^{\prime}(r)\right]+(\kappa_{e}^{2}r^{2}-\tilde{\lambda}-3)c_{1}(r)+\kappa_{o}^{2}r^{2}a_{1}(r)+(\tilde{\lambda}+1)rb_{1}^{\prime}(r)+(\tilde{\lambda}+3)b_{1}(r)\right]\cos(\theta)$
$\displaystyle+\left[\left[r^{2}d_{1}^{\prime\prime}(r)+rd_{1}^{\prime}(r)\right]+(\kappa_{e}^{2}r^{2}-\tilde{\lambda}-3)d_{1}(r)+\kappa_{o}^{2}r^{2}b_{1}(r)-(\tilde{\lambda}+1)ra_{1}^{\prime}(r)-(\tilde{\lambda}+3)a_{1}(r)\right]\sin(\theta)=0.$
(15)

Since each line in the above equations has to vanish separately, these two
equations produce a system of six coupled differential equations for the
coefficients $a_{0},b_{0},a_{1},b_{1},c_{1}$ and $d_{1}$. It is important to
realize that keeping higher order terms in Eq. (II.1) would not change these
equations, but will only add more independent equations for higher order
coefficients.

Figure 4: Panel (a): A typical displacement field resulting from an inflation
of the inner boundary. Panel (b): the radial component of the same
displacement field. Panel (c): the transverse component of the same field.

###  II.2 Analytic solutions

For the case of radial inflation these equations simplify further, since the
boundary conditions are $a_{0}(r=r_{\rm in})=d_{0}$ and all the other
coefficients are zero at $r_{\rm in}$. At $r_{\rm out}$ all the coefficients
must vanish. Thus the only remaining nonzero coefficients are $a_{0}(r)$ and
$b_{0}(r)$. All the higher Fourier modes will have zero values unless we
provide different boundary conditions.

The equations for $a_{0}$ and $b_{0}$ can be solved analytically, as is
demonstrated in Appendix A. The final solutions are shown in Eqs. (55) and
(56). Typical solutions are shown in Figs. 2 and 3.

To compare with numerical simulations, we need to extract data for the
coefficients $a_{0}(r)$ and $b_{0}(r)$ from the measured displacement field.
To this aim we compute the angle averages

$\displaystyle a_{0}(r)=\frac{1}{2\pi}\oint_{0}^{2\pi}d_{r}(r,\theta)d\theta,$
(16) $\displaystyle
b_{0}(r)=\frac{1}{2\pi}\oint_{0}^{2\pi}d_{\theta}(r,\theta)d\theta,$

###  II.3 Comparison with simulations

The description of the numerical procedures for the preparation of our
amorphous solids for the inflation of the inner disk is offered in Appendix B.

A typical displacement field and its radial and transverse components are
shown in Fig. 4. Next we compare to the analytic solutions using Eqs. (16).

When we compare the results of the simulations to the exact solutions of the
type displayed in Figs. 2 and 3, we run into the usual difficulties
encountered when one compares continuum theory to granular simulations. While
one expects excellent agreement in the bulk of the system, on the scale of the
grains one can have deviations. In particular the boundary conditions at
$r_{\rm in}$ can be an issue. While we find that the radial boundary condition
$a_{0}(r=r_{\rm in})=d_{0}$ is faithfully obeyed, the other boundary
condition, i.e. $b_{0}(r=r_{\rm in})=0$ is not obeyed precisely. Even though
the inflation is purely radial, due to the disorder and the interaction of the
inner boundary with the bulk, the close particles touching the boundary
receive some transverse displacement. Thus, in order to compare to the theory,
we measured the transverse displacement at the first layer of disks, and used
this as a boundary condition in the analytic solution. We note that having a
transverse displacement will excite also higher order Fourier modes, but due
to the orthogonality of these modes, the solutions for $a_{0}$ and $b_{0}$
will not be affected. Needless to say, the observed displacement field will
look more complex when more Fourier modes are present.

Typical comparisons of numerical and exact solution are shown in Fig. 5.

Figure 5: Comparisons of angle averaged radial ($a_{0}(r)$) and transverse
($b_{0}(r)$) displacement field with that of analytical solutions (solid
line). Upper figure: $r_{in}=1.25$, $r_{out}=80$, $\mu=3.43$,
$\tilde{\lambda}=68.9$, $k_{e}=0.209$ and $k_{o}=0.206$. Here
$a_{0}(r_{in})=0.16$ and $b_{0}(r_{in})=-0.2$. Middle figure: $r_{in}=1.3$,
$r_{out}=80$, $\mu=5.31$, $\tilde{\lambda}=41.3$, $k_{e}=0.099$ and
$k_{o}=0.149$. Here $a_{0}(r_{in})=0.18$ and $b_{0}(r_{in})=-0.02$. Lower
figure: $r_{in}=1.25$, $r_{out}=80$, $\mu=3.73$, $\tilde{\lambda}=39$,
$k_{e}=0.165$ and $k_{o}=0.155$. Here $a_{0}(r_{in})=0.19$ and
$b_{0}(r_{in})=-0.08$.

We note that the analytic solutions depend on the actual elastic moduli $\mu$
and $\lambda$. We have computed these for the actual realizations after
equilibration, from the microscopic theory that is described for example in
Ref. [17]. In contrast, the screening parameters $\kappa_{e}$ and $\kappa_{o}$
were fitted to the data, and were the only free parameters. We find the
agreement between data and theory quite satisfactory, despite some slight
deviation mainly in $b_{0}(r)$ at higher values of $r$. We therefore turn now
to an even more sensitive test for the fit of the screening parameters, i.e.
the angle $\Phi$ defined in Eq. (7).

##  III The Consequences of Chiral Symmetry Breaking

In this section we explore the angle between the (inverse) direction of the
displacement field and the dipole field. This angle is a direct consequence of
the breaking of Chiral symmetry and it must vanish if $\kappa_{o}=0$, i.e. for
a symmetric tensor ${\bm{\Gamma}}$. To this aim we need to extract the dipole
field and find the angle between it and the inverse displacement direction.

###  III.1 Extracting the Quadrupolar and Dipole Fields from displacement data

To accomplish our aim we must extract first the quadrupolar field associated
with the displacement field. Since a quadrupole is represented by the
traceless part of the non-affine strain tensor, so we first compute the total
strain tensor. We follow the procedure advised in Ref. [18], using the
discrete deformation gradient $F$ defined as,

$\displaystyle\varepsilon=\frac{1}{2}\left[I-\left(F^{\text{T}}F\right)^{-1}\right],$
(17)

where $F\equiv\frac{\partial\chi}{\partial\bf{X}}=\frac{\partial
x}{\partial\bf{X}}$ is a second rank tensor. Here $\chi$ maps a particle
position from the initial undeformed configuration $\Omega_{0}$ to final
deformed configuration $\Omega_{D}$ defined by $\bf{x=\chi(X)}$. Here $\bf{X}$
and $\bf{x}$ are the positions of the particles before and after the
deformation in the two configurations.

To calculate $F$ in the atomistic simulation we measure the relative
displacements between two particles. Before we proceed, we would like to
mention that the particle positions have been mapped to grid points.
Therefore, in all our subsequent discussion particle labels and grid points
are synonyms. Let $\Delta X^{ij}=X^{j}-X^{i}$, and $\Delta x^{ij}=x^{j}-x^{i}$
be the relative displacements in $\Omega_{0}$ and $\Omega_{D}$ respectively
between the particles $i$ and $j$. The transformation from configuration
$\Omega_{0}$ to $\Omega_{D}$ is defined by deformation gradient $F$, which
linearly maps the relative position vector $\bf{\Delta X^{ij}}$ to $\bf{\Delta
x^{ij}}$ as follows,

$\displaystyle\Delta x^{ij}=F^{i}\Delta X^{ij},$ (18)

where $F^{i}$ is the deformation gradient measured at the particle $i$. Note
that, instead of $F$, we use a local deformation gradient $\hat{F}$, which is
given by [18],

$\hat{F}D=A\leavevmode\nobreak\ \leavevmode\nobreak\
\text{or}\leavevmode\nobreak\ \leavevmode\nobreak\ \hat{F}=AD^{-1},$ (19)

where $A$ and $D$ are $2\times 2$ matrices, which are obtained from

$\displaystyle{A=\sum_{j=1}^{L}\Delta
x^{ij}\Delta{X^{ij}}^{\text{T}}w_{j}},\leavevmode\nobreak\
{D=\sum_{j=1}^{L}\Delta X^{ij}\Delta{X^{ij}}^{\text{T}}w_{j}}.$ (20)

Here $w_{j}$ is the weight function assigned to particle $j$ depending upon
its relative distance with the central atom $i$, and $L$ denotes the total
neighbors around particle $i$ [18].

The non-affine strain $\varepsilon^{q}$ is obtained after subtracting the
affine strain components $\varepsilon^{a}$ from the total strain $\varepsilon$
measured in Eqn.(17),

$\displaystyle\varepsilon_{\alpha\beta}^{q}=\varepsilon_{\alpha\beta}-\varepsilon_{\alpha\beta}^{a}=\begin{bmatrix}\varepsilon_{11}^{q}&\varepsilon_{12}^{q}\\\
\varepsilon_{21}^{q}&\varepsilon_{22}^{q}\end{bmatrix}.$ (21)

Here
$\varepsilon_{\alpha\beta}^{a}=0.5(\nabla_{\alpha}d_{\beta}+\nabla_{\beta}d_{\alpha})$,
and $\alpha,\beta=1,2$ which are same as $(x,y)$. Note that
$\varepsilon_{\alpha\beta}^{q}$ is calculated at each particle label or each
grid point. Here, $d$ is the quasi-elastic displacement field in purely radial
inflation which is given by [1]:

$d(r)=d_{o}\left[\frac{r^{2}-r_{out}^{2}}{r_{in}^{2}-r_{out}^{2}}\right]\left(\frac{r_{in}}{r}\right).$
(22)

Figure 6: Heat map of the quadrupolar field in typical radial inflation. The
dense regions indicate high values of $Q$, and less dense regions correspond
to the low values of $Q$ as obtained from Eq. (24). The arrows represent the
principal axis or the directions of the $Q$ field, cf. Eq. (26). The
parameters here are $r_{in}=1.25$, $r_{out}=80$, $\mu=3.43$,
$\tilde{\lambda}=68.9$, $k_{e}=0.209$ and $k_{o}=0.206$. Here
$a_{0}(r_{in})=0.16$ and $b_{0}(r_{in})=-0.2$. In the top panel we represent
the full quadrupolar field, whereas in the bottom panel we show the magnified
view of the most dense region.

The non-affine strain $\varepsilon^{q}$ can be decomposed into its trace and
its traceless components [10]:

$\displaystyle\varepsilon^{q}$
$\displaystyle=\frac{1}{2}\begin{bmatrix}\varepsilon_{11}^{q}+\varepsilon_{22}^{q}&0\\\
0&\varepsilon_{11}^{q}+\varepsilon_{22}^{q}\end{bmatrix}+Q\begin{bmatrix}\varepsilon_{11}^{ts}&\varepsilon_{12}^{ts}\\\
\varepsilon_{21}^{ts}&\varepsilon_{22}^{ts}\end{bmatrix},$
$\displaystyle=0.5\times mI+Q\varepsilon^{ts}.$ (23)

Here $I$ is the identity tensor,
$m=(\varepsilon_{11}^{q}+\varepsilon_{22}^{q}$) is the trace of the tensor
$\varepsilon^{q}$, and $\varepsilon^{ts}$ is traceless symmetric tensor
representing a quadrupole, where
$\varepsilon_{22}^{ts}=-\varepsilon_{11}^{ts}$. The quantity $Q$ is the
quadrupolar charge associated with the quadrupolar tensor and is defined as,

$Q^{2}=(\varepsilon_{11}^{ts})^{2}+(\varepsilon_{22}^{ts})^{2}.$ (24)

Eq. (III.1) can also be written in the following form,

$\displaystyle\varepsilon^{q}$
$\displaystyle=\frac{1}{2}\begin{bmatrix}\varepsilon_{11}^{q}+\varepsilon_{22}^{q}&0\\\
0&\varepsilon_{11}^{q}+\varepsilon_{22}^{q}\end{bmatrix}+Q\begin{bmatrix}\cos(2\Theta)&\sin(2\Theta)\\\
\sin(2\Theta)&-\cos(2\Theta)\end{bmatrix}.$ (25)

Therefore the orientation of the quadrupolar field is given by,

$\Theta=\frac{1}{2}\arctan\left(\frac{\varepsilon_{12}^{ts}}{\varepsilon_{11}^{ts}}\right).$
(26)

$\Theta$ defines the principal axis of the quadrupolar field. We show a
typical quadrupolar field in Fig. 6.

Next we calculate the dipole field using the gradient of the quadrupolar
field, ${\mathcal{P}}^{\alpha}\equiv\partial_{\beta}Q^{\alpha\beta}\ $ at each
grid point. The dipole field associated with the quadrupolar field in Fig. 6
is shown in figure 7.

Figure 7: Heat map of the dipole field after a typical radial inflation,
obtained by taking the gradient of the quadrupolar field
${\mathcal{P}}^{\alpha}\equiv\partial_{\beta}Q^{\alpha\beta}\ $. The dense
regions indicate high values of ${\mathcal{P}}\ $, and less dense regions
correspond to the low values of dipole field ${\mathcal{P}}\ $. The arrows
represent the local direction of the dipole field. The parameters here are are
$r_{in}=1.25$, $r_{out}=80$, $\mu=3.43$, $\tilde{\lambda}=68.9$, $k_{e}=0.209$
and $k_{o}=0.206$. Here $a_{0}(r_{in})=0.16$ and $b_{0}(r_{in})=-0.2$. In the
top panel we represent the full dipole field, whereas in the bottom panel we
show the magnified view of the most dense region.

Figure 8: The dipole and displacement vector fields shown together, red arrows
for the dipole field and blue arrows for displacement. The parameters here are
those of Fig. 7. In the top panel we represent the full dipole and
displacement fields, whereas in the bottom panel we show the magnified view of
the most dense region. The averaged angle between the dipole and the inverse
displacement field from the plots is $46.48$ degrees, whereas theoretically we
obtain
$\theta=\arctan\left(\frac{\kappa_{o}^{2}}{\kappa_{e}^{2}}\right)=44.11$
degrees.

Having the dipole field we compute the angle between the inverse displacement
and the dipole field, which is predicted to be the angle $\Phi$ of Eq. (7).
For the present example, we found from simulations an angle of $46.48^{o}\pm
0.55^{o}$ as compared to an angle of $44.11^{o}$ from the direct calculation
of Eq. (7). This close correspondence is typical and similar results were
obtained for other simulations of radial inflation.

##  IV Summary and discussion

The aim of this paper was to explore further the consequences of Chiral
symmetry breaking in the mechanical response of amorphous solids. This
symmetry breaking is additional to the breaking of translational symmetry. The
latter is responsible for the introduction of a typical screening length,
which in the notation of the present paper is $\kappa_{e}^{-1}$. Translational
symmetry breaking due to dipole screening is also the mechanism behind the
hexatic phase transition [19] and the Kosterlitz-Thouless transition [20]. The
appearance of this Chiral symmetry breaking seems novel to the mechanical
response of amorphous solids, and is related to irreversibility and energy
loss (or gain) in a closed loop of straining an amorphous solids. It can be
shown that the presence of the odd tensor ${\bm{\Gamma}}$ in the form of Eq.
(6) is sufficient to render a closed loop of straining dissipative, in the
sense that energy is either lost or gained by the material [7]. It is
worthwhile to remark that a similar notion of “Odd Elasticity" was coined
recently in the context of materials whose mechanics is not derived from a
potential [21]. The present formulation thus may be referred to as “Odd
Plasticity", and it also connects to the lack of energy conservation.

The methods employed above translate naturally to experimental data. All that
one needs is an accurate map of the displacement field, something that was
shown to be very doable in experiments using transparent plastic disks, see
for example Ref. [6]. Once the displacement field is given, the protocol
followed in the previous section is relatively straightforward. It seems quite
obvious that the interest in the presence of Chiral symmetry breaking warrants
a serious experimental attempt to corroborate the simulation results shown
above.

An issue that was not explicitly discussed in the present paper is the
transition from elastic response to anomalous one. It was shown before that
given glass formers or granular system can show both behaviors with an
interesting transition between them, as a function of a control parameter like
pressure, [4]. We now expect that the transition will be accompanied by both
translation and Chiral symmetry breaking, requiring additional scrutiny of the
nature of this transition. This is work in progress in our group that will be
reported in the near future.

Finally, dynamical response to oscillatory or any form of dynamical straining
are known to show dipole screening as well, cf. Ref. [5]. The presence of
Chiral symmetry breaking and it consequences needs to be examined for such
protocols as well.

###### Acknowledgements.

This work had been supported in part by ISF under grant #3492/21
(collaboration with China) and the Minerva Center for “Aging, from physical
materials to human tissues" at the Weizmann Institute.

##  Appendix A Analytic solutions of the equations

In this appendix we solve the Eqs. (10) and (11) for the analytical form of
the radial and transverse displacements functions $d_{r}$ and $d_{\theta}$
respectively. We consider the following forms of the $d_{r}$ and $d_{\theta}$,

$\displaystyle d_{r}(r,\theta)=a_{0}(r)$ $\displaystyle
d_{\theta}(r,\theta)=b_{0}(r).$ (27)

Thus equations (10 and 11) reduce to following coupled equations in $a_{0}(r)$
and $b_{0}(r)$,

$\displaystyle[r^{2}a_{0}^{\prime\prime}+ra_{0}^{\prime}-a_{0}]+\frac{\kappa_{e}^{2}r^{2}}{(\tilde{\lambda}+2)}a_{0}-\frac{\kappa_{o}^{2}r^{2}}{(\tilde{\lambda}+2)}b_{0}$
$\displaystyle=0,$
$\displaystyle[r^{2}b_{0}^{\prime\prime}+rb_{0}^{\prime}-b_{0}]+\kappa_{e}^{2}r^{2}b_{0}+\kappa_{o}^{2}r^{2}a_{0}$
$\displaystyle=0.$ (28)

We can combine the above two coupled equations in $a_{0}(r)$ and $b_{0}(r)$ by
using two Lagrange multipliers $C_{1}$ and $C_{2}$ as follows

$\displaystyle
C_{1}\left\\{[r^{2}a_{0}^{\prime\prime}+ra_{0}^{\prime}-a_{0}]+\frac{\kappa_{e}^{2}r^{2}}{(\tilde{\lambda}+2)}a_{0}-\frac{\kappa_{o}^{2}r^{2}}{(\tilde{\lambda}+2)}b_{0}\right\\}+C_{2}\left\\{[r^{2}b_{0}^{\prime\prime}+rb_{0}^{\prime}-b_{0}]+\kappa_{e}^{2}r^{2}b_{0}+\kappa_{o}^{2}r^{2}a_{0}\right\\}=0,$
(29)

which after some simplification can be written as,

$\displaystyle
r^{2}\left[C_{1}a_{0}^{\prime\prime}+C_{2}b_{0}^{\prime\prime}\right]+r\left[C_{1}a_{0}^{\prime}+C_{2}b_{0}^{\prime}\right]-\left[C_{1}a_{0}+C_{2}b_{0}\right]+\left[C_{1}\frac{\kappa_{e}^{2}r^{2}}{(\tilde{\lambda}+2)}+C_{2}\kappa_{o}^{2}r^{2}\right]a_{0}+\left[C_{2}\kappa_{e}^{2}r^{2}-C_{1}\frac{\kappa_{o}^{2}r^{2}}{(\tilde{\lambda}+2)}\right]b_{0}=0.$
(30)

The above equation can be written in the form of Bessel differential equation.
To show that this is true, we write the above equation in the following form

$r^{2}\left[C_{1}a_{0}^{\prime\prime}+C_{2}b_{0}^{\prime\prime}\right]+r\left[C_{1}a_{0}^{\prime}+C_{2}b_{0}^{\prime}\right]-\left[C_{1}a_{0}+C_{2}b_{0}\right]\\\
+r^{2}[\tilde{C_{1}}a_{0}+\tilde{C_{2}}b_{0}]=0,$ (31)

where,

$\displaystyle\tilde{C_{1}}=\left[C_{1}\frac{\kappa_{e}^{2}}{(\tilde{\lambda}+2)}+C_{2}\kappa_{o}^{2}\right],$
$\displaystyle\tilde{C_{2}}=\left[C_{2}\kappa_{e}^{2}-C_{1}\frac{\kappa_{o}^{2}}{(\tilde{\lambda}+2)}\right].$
(32)

Now let us substitute

$\displaystyle
r^{2}[\tilde{C_{1}}a_{0}+\tilde{C_{2}}b_{0}]=\kappa^{2}r^{2}[C_{1}a_{0}+C_{2}b_{0}],$
(33)

where $\kappa$ is a parameter to be obtained such that it defines the above
transformation. Thus

$\displaystyle[\tilde{C_{1}}a_{0}+\tilde{C_{2}}b_{0}]=[\kappa^{2}C_{1}a_{0}+\kappa^{2}C_{2}b_{0}],$
(34)

or,

$\displaystyle[\tilde{C_{1}}-\kappa^{2}C_{1}]a_{0}+[\tilde{C_{2}}-\kappa^{2}C_{2}]b_{0}=0.$
(35)

Since, $a_{0}(r)$ and $b_{0}(r)$ are arbitrary, therefore to hold the above
equation true, we must have

$\displaystyle\tilde{C_{1}}=\kappa^{2}C_{1},$
$\displaystyle\tilde{C_{2}}=\kappa^{2}C_{2}.$ (36)

Now we can solve these two equations together to determine the value of
$\kappa$. After substituting the values of $\tilde{C_{1}}$ and
$\tilde{C_{2}}$, above equations reduce to

$\displaystyle
C_{1}\frac{\kappa_{e}^{2}}{(\tilde{\lambda}+2)}+C_{2}\kappa_{o}^{2}=\kappa^{2}C_{1},$
$\displaystyle
C_{2}\kappa_{e}^{2}-C_{1}\frac{\kappa_{o}^{2}}{(\tilde{\lambda}+2)}=\kappa^{2}C_{2}.$
(37)

With further simplifications we have

$\displaystyle\frac{\kappa_{e}^{2}}{(\tilde{\lambda}+2)}+\frac{C_{2}}{C_{1}}\kappa_{o}^{2}=\kappa^{2},$
$\displaystyle\kappa_{e}^{2}-\frac{C_{1}}{C_{2}}\frac{\kappa_{o}^{2}}{(\tilde{\lambda}+2)}=\kappa^{2}.$
(38)

Let us take $\frac{C_{1}}{C_{2}}=Z$, then the above equations simplify to

$\displaystyle\frac{\kappa_{e}^{2}}{(\tilde{\lambda}+2)}+\frac{1}{Z}\kappa_{o}^{2}=\kappa^{2},$
$\displaystyle\kappa_{e}^{2}-Z\frac{\kappa_{o}^{2}}{(\tilde{\lambda}+2)}=\kappa^{2}.$
(39)

Now we have two equations giving the same values of $\kappa$. We equate these
two equation to give a quadratic equation in $Z$,

$\displaystyle\frac{\kappa_{e}^{2}}{(\tilde{\lambda}+2)}+\frac{1}{Z}\kappa_{o}^{2}=\kappa_{e}^{2}-Z\frac{\kappa_{o}^{2}}{(\tilde{\lambda}+2)}.$
(40)

Which simplifies to

$\displaystyle
Z^{2}\frac{\kappa_{o}^{2}}{\tilde{\lambda}+2}+\left(\frac{\kappa_{e}^{2}}{\tilde{\lambda}+2}-\kappa_{e}^{2}\right)Z+\kappa_{o}^{2}=0,$
(41)

or

$\displaystyle\kappa_{o}^{2}Z^{2}-\kappa_{e}^{2}(\tilde{\lambda}+1)Z+\kappa_{o}^{2}(\tilde{\lambda}+2)=0.$
(42)

Now if we substitute $Z=\frac{C_{1}}{C_{2}}$, then the above equation can be
written in terms of $C_{1}$ and $C_{2}$ as follows

$\displaystyle\kappa_{o}^{2}C_{1}^{2}-\kappa_{e}^{2}(\tilde{\lambda}+1)C_{1}C_{2}+\kappa_{o}^{2}(\tilde{\lambda}+2)C_{2}^{2}=0.$
(43)

The solutions of $Z$ are

$Z=\frac{\kappa_{e}^{2}(\tilde{\lambda}+1)\pm\sqrt{\kappa_{e}^{4}(\tilde{\lambda}+1)^{2}-4(\tilde{\lambda}+2)\kappa_{o}^{4}}}{2\kappa_{o}^{2}},$
(44)

or

$Z_{1}=\frac{\kappa_{e}^{2}(\tilde{\lambda}+1)+\sqrt{\kappa_{e}^{4}(\tilde{\lambda}+1)^{2}-4(\tilde{\lambda}+2)\kappa_{o}^{4}}}{2\kappa_{o}^{2}},$
(45)

and

$Z_{2}=\frac{\kappa_{e}^{2}(\tilde{\lambda}+1)-\sqrt{\kappa_{e}^{4}(\tilde{\lambda}+1)^{2}-4(\tilde{\lambda}+2)\kappa_{o}^{4}}}{2\kappa_{o}^{2}}.$
(46)

Substituting the values of $Z$ from above equations into Eqn.(A), we obtain
the values of $\kappa$ defining the transformation in Eqn.(33). Therefore,
from equations (31, and 33) we have

$r^{2}\left[C_{1}a_{0}^{\prime\prime}+C_{2}b_{0}^{\prime\prime}\right]+r\left[C_{1}a_{0}^{\prime}+C_{2}b_{0}^{\prime}\right]-\left[C_{1}a_{0}+C_{2}b_{0}\right]\\\
+\kappa^{2}r^{2}[C_{1}a_{0}+C_{2}b_{0}]=0,$ (47)

or

$r^{2}\left[C_{1}a_{0}^{\prime\prime}+C_{2}b_{0}^{\prime\prime}\right]+r\left[C_{1}a_{0}^{\prime}+C_{2}b_{0}^{\prime}\right]\\\
+(\kappa^{2}r^{2}-1)[C_{1}a_{0}+C_{2}b_{0}]=0.$ (48)

We can simplify it little bit more to give

$r^{2}\left[\frac{C_{1}}{C_{2}}a_{0}^{\prime\prime}+b_{0}^{\prime\prime}\right]+r\left[\frac{C_{1}}{C_{2}}a_{0}^{\prime}+b_{0}^{\prime}\right]\\\
+(\kappa^{2}r^{2}-1)\left[\frac{C_{1}}{C_{2}}a_{0}+b_{0}\right]=0.$ (49)

Now, if we substitute

$\displaystyle X(r)$ $\displaystyle=\frac{C_{1}}{C_{2}}a_{0}(r)+b_{0}(r)$
$\displaystyle=\leavevmode\nobreak\ \leavevmode\nobreak\ Za_{0}(r)+b_{0}(r),$
(50)

in the Eqn.(49), we obtain the following differential equation

$\displaystyle r^{2}X^{\prime\prime}+rX^{\prime}+(\kappa^{2}r^{2}-1)X=0.$ (51)

This is a bessel differential equation, where $\kappa$ and $Z$ are already
defined above. A general solution of this equation is

$\displaystyle X(r)=mJ_{1}(\kappa r)+nY_{1}(\kappa r),$ (52)

where $J_{1}$ and $Y_{1}$ are the Bessel functions of first kind, and the
coefficients $m$ and $n$ are the constant parameters to be obtained using the
boundary conditions. Note that, we will have two solutions corresponding to
the two values of $\kappa$ or $Z$. The two values of $Z$ are $Z_{1}$ and
$Z_{2}$, and the corresponding values of $\kappa$ are $\eta$ and $\zeta$
respectively. Then we obtain following two coupled equations in $a_{0}(r)$ and
$b_{0}(r)$ from the equations (A), and (52),

$\displaystyle m_{1}J_{1}(\eta r)+n_{1}Y_{1}(\eta r)=Z_{1}a_{0}(r)+b_{0}(r),$
(53)

and

$\displaystyle m_{2}J_{1}(\zeta r)+n_{2}Y_{1}(\zeta
r)=Z_{2}a_{0}(r)+b_{0}(r).$ (54)

From the above two equations we immediately obtain the analytical forms for
$a_{0}(r)$, and $b_{0}(r)$

$\displaystyle a_{0}(r)=\frac{\left[[m_{1}J_{1}(\eta r)+n_{1}Y_{1}(\eta
r)]-[m_{2}J_{1}(\zeta r)+n_{2}Y_{1}(\zeta r)]\right]}{Z_{1}-Z_{2}},$ (55)
$\displaystyle b_{0}(r)=\frac{Z_{2}[m_{1}J_{1}(\eta r)+n_{1}Y_{1}(\eta
r)]-Z_{1}[m_{2}J_{1}(\zeta r)+n_{2}Y_{1}(\zeta
r)]}{\left(Z_{2}-Z_{1}\right)}.$ (56)

Now we use boundary conditions to determine the coefficients $m_{1},n_{1}$ and
$m_{2},n_{2}$. The boundary conditions on $a_{0}(r)$ and $b_{0}(r)$ are

$\displaystyle a_{0}(r)|_{r=r_{in}}=d_{a},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ a_{0}(r)|_{r=r_{o}}=0,$
$\displaystyle b_{0}(r)|_{r=r_{in}}=d_{b},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ b_{0}(r)|_{r=r_{o}}=0.$ (57)

With the above boundary conditions, we have following four coupled equations
in $m_{1},n_{1}$ and $m_{2},n_{2}$,

$\displaystyle[m_{1}J_{1}(\eta r_{in})+n_{1}Y_{1}(\eta
r_{in})]-[m_{2}J_{1}(\zeta r_{in})+n_{2}Y_{1}(\zeta r_{in})]$
$\displaystyle=(Z_{1}-Z_{2})d_{a},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ (i)$
$\displaystyle[m_{1}J_{1}(\eta r_{o})+n_{1}Y_{1}(\eta
r_{o})]-[m_{2}J_{1}(\zeta r_{o})+n_{2}Y_{1}(\zeta r_{o})]$
$\displaystyle=0,\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ (ii)$ $\displaystyle
Z_{2}[m_{1}J_{1}(\eta r_{in})+n_{1}Y_{1}(\eta r_{in})]-Z_{1}[m_{2}J_{1}(\zeta
r_{in})+n_{2}Y_{1}(\zeta r_{in})]$
$\displaystyle=(Z_{2}-Z_{1})d_{b},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ (iii)$
$\displaystyle Z_{2}[m_{1}J_{1}(\eta r_{o})+n_{1}Y_{1}(\eta
r_{o})]-Z_{1}[m_{2}J_{1}(\zeta r_{o})+n_{2}Y_{1}(\zeta r_{o})]$
$\displaystyle=0.\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ (iv)$

From Eqn. A$(ii)$, we have

$\displaystyle[m_{1}J_{1}(\eta r_{o})+n_{1}Y_{1}(\eta
r_{o})]=[m_{2}J_{1}(\zeta r_{o})+n_{2}Y_{1}(\zeta r_{o})],$ (59)

and using this in equation A$(iv)$, we get

$\displaystyle[m_{2}J_{1}(\zeta r_{o})+n_{2}Y_{1}(\zeta
r_{o})](Z_{2}-Z_{1})=0.$ (60)

Since $Z_{1}\neq Z_{2}$, therefore $(Z_{2}-Z_{1})\neq 0$. Thus we must have

$\displaystyle m_{2}J_{1}(\zeta r_{o})+n_{2}Y_{1}(\zeta r_{o})=0,$ (61)

which gives

$\displaystyle m_{2}=-n_{2}\frac{Y_{1}(\zeta r_{o})}{J_{1}(\zeta r_{o})}.$
(62)

Similarly from Eqns. A$(ii)$, and A$(iv)$, we have

$\displaystyle[m_{1}J_{1}(\eta r_{o})+n_{1}Y_{1}(\eta r_{o})](Z_{2}-Z_{1})=0,$
(63)

which gives

$\displaystyle m_{1}=-n_{1}\frac{Y_{1}(\eta r_{o})}{J_{1}(\eta r_{o})}.$ (64)

Substituting the values of $m_{1}$ and $m_{2}$ into Eqn. A$(i)$, we obtain,

$\displaystyle n_{1}\left[\frac{Y_{1}(\eta r_{in})J_{1}(\eta r_{o})-Y_{1}(\eta
r_{o})J_{1}(\eta r_{in})}{J_{1}(\eta
r_{o})}\right]-n_{2}\left[\frac{Y_{1}(\zeta r_{in})J_{1}(\zeta
r_{o})-Y_{1}(\zeta r_{o})J_{1}(\zeta r_{in})}{J_{1}(\zeta
r_{o})}\right]=(Z_{1}-Z_{2})d_{a}.$ (65)

Similarly substituting the values of $m_{1}$ and $m_{2}$ into Eqn. A$(iii)$,
we obtain,

$\displaystyle n_{1}Z_{2}\left[\frac{Y_{1}(\eta r_{in})J_{1}(\eta
r_{o})-Y_{1}(\eta r_{o})J_{1}(\eta r_{in})}{J_{1}(\eta
r_{o})}\right]-n_{2}Z_{1}\left[\frac{Y_{1}(\zeta r_{in})J_{1}(\zeta
r_{o})-Y_{1}(\zeta r_{o})J_{1}(\zeta r_{in})}{J_{1}(\zeta
r_{o})}\right]=(Z_{2}-Z_{1})d_{b}.$ (66)

For simplicity let us define

$\displaystyle\alpha=\left[\frac{Y_{1}(\eta r_{in})J_{1}(\eta
r_{o})-Y_{1}(\eta r_{o})J_{1}(\eta r_{in})}{J_{1}(\eta r_{o})}\right],$

and

$\displaystyle\beta=\left[\frac{Y_{1}(\zeta r_{in})J_{1}(\zeta
r_{o})-Y_{1}(\zeta r_{o})J_{1}(\zeta r_{in})}{J_{1}(\zeta r_{o})}\right].$

Thus Eqns. (65 and 66) reduce to

$\displaystyle n_{1}\alpha-n_{2}\beta=(Z_{1}-Z_{2})d_{a},$ (69)

and,

$\displaystyle n_{1}\alpha Z_{2}-n_{2}\beta Z_{1}=(Z_{2}-Z_{1})d_{b}.$ (70)

The above two equations can be solved together to give values of $n_{1}$ and
$n_{2}$

$\displaystyle n_{1}$ $\displaystyle=\frac{Z_{1}d_{a}+d_{b}}{\alpha}$
$\displaystyle=\frac{(Z_{1}d_{a}+d_{b})J_{1}(\eta r_{o})}{[Y_{1}(\eta
r_{in})J_{1}(\eta r_{o})-Y_{1}(\eta r_{o})J_{1}(\eta r_{in})]},$ (71)

and

$\displaystyle n_{2}$ $\displaystyle=\frac{Z_{2}d_{a}+d_{b}}{\beta}$
$\displaystyle=\frac{(Z_{2}d_{a}+d_{b})J_{1}(\zeta r_{o})}{[Y_{1}(\zeta
r_{in})J_{1}(\zeta r_{o})-Y_{1}(\zeta r_{o})J_{1}(\zeta r_{in})]}.$ (72)

Now we can write the values of $m_{1}$ and $m_{2}$ from Eqs.(62 and 64)

$\displaystyle m_{1}=-\frac{(Z_{1}d_{a}+d_{b})Y_{1}(\eta r_{o})}{[Y_{1}(\eta
r_{in})J_{1}(\eta r_{o})-Y_{1}(\eta r_{o})J_{1}(\eta r_{in})]},$
$\displaystyle m_{2}=-\frac{(Z_{2}d_{a}+d_{b})Y_{1}(\zeta r_{o})}{[Y_{1}(\zeta
r_{in})J_{1}(\zeta r_{o})-Y_{1}(\zeta r_{o})J_{1}(\zeta r_{in})]}.$ (73)

We list below all the parameters that are needed to determine the functional
forms of $a_{0}(r)$ and $b_{0}(r)$. The values of $Z$ are

$\displaystyle
Z_{1}=\frac{\kappa_{e}^{2}(\tilde{\lambda}+1)+\sqrt{\kappa_{e}^{4}(\tilde{\lambda}+1)^{2}-4(\tilde{\lambda}+2)\kappa_{o}^{4}}}{2\kappa_{o}^{2}},$
$\displaystyle
Z_{2}=\frac{\kappa_{e}^{2}(\tilde{\lambda}+1)-\sqrt{\kappa_{e}^{4}(\tilde{\lambda}+1)^{2}-4(\tilde{\lambda}+2)\kappa_{o}^{4}}}{2\kappa_{o}^{2}}.$
(74)

The two values of $\kappa$, i.e., $\eta$ and $\zeta$ are

$\displaystyle\eta=\sqrt{\kappa_{e}^{2}-Z_{1}\frac{\kappa_{o}^{2}}{(\tilde{\lambda}+2)}},$
$\displaystyle\zeta=\sqrt{\kappa_{e}^{2}-Z_{2}\frac{\kappa_{o}^{2}}{(\tilde{\lambda}+2)}}.$
(75)

and the parameters $m_{1},n_{1},m_{2},n_{2}$ are

$\displaystyle n_{1}=\frac{(Z_{1}d_{a}+d_{b})J_{1}(\eta r_{o})}{[Y_{1}(\eta
r_{in})J_{1}(\eta r_{o})-Y_{1}(\eta r_{o})J_{1}(\eta r_{in})]},$
$\displaystyle m_{1}=-\frac{(Z_{1}d_{a}+d_{b})Y_{1}(\eta r_{o})}{[Y_{1}(\eta
r_{in})J_{1}(\eta r_{o})-Y_{1}(\eta r_{o})J_{1}(\eta r_{in})]},$
$\displaystyle n_{2}=\frac{(Z_{2}d_{a}+d_{b})J_{1}(\zeta r_{o})}{[Y_{1}(\zeta
r_{in})J_{1}(\zeta r_{o})-Y_{1}(\zeta r_{o})J_{1}(\zeta r_{in})]},$
$\displaystyle m_{2}=-\frac{(Z_{2}d_{a}+d_{b})Y_{1}(\zeta r_{o})}{[Y_{1}(\zeta
r_{in})J_{1}(\zeta r_{o})-Y_{1}(\zeta r_{o})J_{1}(\zeta r_{in})]}.$ (76)

##  Appendix B protocols of the preparation of the amorphous solid

Here we outline the details of the numerical simulation to produce the glassy
configurations. We study a two-dimensional poly-dispersed model of $N=20,000$
point particles in an annulus bounded by two rigid walls, with an inner
radius, $r_{in}$ and outer radius, $r_{out}$. The particles are filled in
annular region such that density of the system, $\rho=N/A=1.0$, where
$A=\pi(r_{out}^{2}-r_{in}^{2})$. The binary interactions between point
particles with mass, $m=1$ is given Lennard-Jones (LJ) potential:

$u_{ij}(r)=\begin{cases}u_{ij}^{LJ}+A_{ij}+B_{ij}r+C_{ij}r^{2},\;\;\textbf{if}\;R\leq
R_{ij}^{cut}\\\ 0,\hskip 85.35826pt\text{otherwise}\end{cases}$ (77)

where

$u_{ij}^{LJ}=4\epsilon_{ij}\Big{[}\Big{(}\frac{\sigma_{ij}}{r}\Big{)}^{12}-\Big{(}\frac{\sigma_{ij}}{r}\Big{)}^{6}\Big{]}$
(78)

The parameters $A_{ij}=0.4526\epsilon_{ij}$,
$B=-0.3100\epsilon_{ij}/\sigma_{ij}$, $C=0.0542\epsilon_{ij}/\sigma_{ij}^{2}$
are added to smooth the potential, $u_{ij}(r)$ at cut-off distance,
$R_{ij}^{cut}=2.5\sigma_{ij}$ (upto second derivative). The interaction length
of particles, $\sigma_{i}$ is drawn from a probability distribution,
$P(\sigma)\simeq 1/\sigma^{3}$ in a range between $\sigma_{max}$ and
$\sigma_{min}$ such that the mean $\bar{\sigma}=1$. The mixing rule of
$\sigma_{i}$ for binary interaction are:

$\begin{split}\sigma_{ij}=\frac{\sigma_{i}+\sigma_{j}}{2}\Big{[}1-0.2|\sigma_{i}-\sigma_{j}|\Big{]},\\\
\sigma_{max}=1.61,\sigma_{min}=\sigma_{max}/2.219\end{split}$ (79)

The reduced units for mass, length, energy and time are $m$, $\bar{\sigma}$,
$\epsilon_{ij}=1$ and $\bar{\sigma}\sqrt{m/\epsilon_{ij}}$ respectively. The
interaction between point particles and the walls is of same form as Eq 77,
where $r_{ij}$ is replaced by distance to the wall.

The system is first thermalized at “mother temperature”, $T_{m}=1$ using swap
Monte Carlo and cooled down to $T=0$ using conjugate gradient method. Once the
system is mechanically equilibrated with total force on each particle smaller
than $10^{-8}$, we inflate the inner radius $r_{in}$ such that inner radius
after inflation is $r_{in}+d_{0}$.

After inflation, we equilibrate the system again using conjugate gradient
method and measure the displacement field, $\bm{d}$ by comparing the
equilibrated configurations before and after inflation. To compare with the
theory, we compute the angle-averaged radial component of the radial
displacement field, $d_{r}(r)$.

## References

  * Lemaître  _et al._ [2021] A. Lemaître, C. Mondal, M. Moshe, I. Procaccia, S. Roy, and K. Screiber-Re’em, Anomalous elasticity and plastic screening in amorphous solids, Phys. Rev. E 104, 024904 (2021).
  * Charan  _et al._ [2023] H. Charan, M. Moshe, and I. Procaccia, Anomalous elasticity and emergent dipole screening in three-dimensional amorphous solids, Phys. Rev. E 107, 055005 (2023).
  * Kumar and Procaccia [2024] A. Kumar and I. Procaccia, Elasticity, plasticity and screening in amorphous solids: A short review, Europhysics Letters 145, 26002 (2024).
  * Jin  _et al._ [2024] Y. Jin, I. Procaccia, and T. Samanta, Intermediate phase between jammed and unjammed amorphous solids, Phys. Rev. E 109, 014902 (2024).
  * Hentschel  _et al._ [2024] H. G. E. Hentschel, A. Pomyalov, I. Procaccia, and O. Szachter, Dynamic screening by plasticity in amorphous solids, Phys. Rev. E 109, 044902 (2024).
  * Mondal  _et al._ [2022] C. Mondal, M. Moshe, I. Procaccia, S. Roy, J. Shang, and J. Zhang, Experimental and numerical verification of anomalous screening theory in granular matter, Chaos, Solitons and Fractals 164, 112609 (2022).
  * Cohen  _et al._ [2024] Y. Cohen, A. Schiller, D. Wang, J. Dijksman, and M. Moshe, Odd dipole screening in disordered matter (2024), arXiv:2310.09942 [cond-mat.soft] .
  * Bhowmik  _et al._ [2022] B. P. Bhowmik, M. Moshe, and I. Procaccia, Direct measurement of dipoles in anomalous elasticity of amorphous solids, Phys. Rev. E 105, L043001 (2022).
  * Kumar  _et al._ [2022] A. Kumar, M. Moshe, I. Procaccia, and M. Singh, Anomalous elasticity in classical glass formers, Phys. Rev. E 106, 015001 (2022).
  * Mondal  _et al._ [2023] C. Mondal, M. Moshe, I. Procaccia, and S. Roy, Dipole screening in pure shear strain protocols of amorphous solids (2023), arXiv:2305.11253 [cond-mat.soft Phys.Rev. E, in press] .
  * Landau and Lifshitz [1959] L. D. Landau and E. M. Lifshitz,  _Course of Theoretical Physics Vol 7: Theory of Elasticity_ (Pergamon Press, 1959).
  * Karmakar  _et al._ [2010] S. Karmakar, E. Lerner, and I. Procaccia, Athermal nonlinear elastic constants of amorphous solids, Phys. Rev.E 82, 026105 (2010).
  * Hentschel  _et al._ [2011] H. G. E. Hentschel, S. Karmakar, E. Lerner, and I. Procaccia, Do athermal amorphous solids exist?, Phys. Rev.E 83, 061101 (2011).
  * Malandro and Lacks [1999] D. L. Malandro and D. J. Lacks, Relationships of shear-induced changes in the potential energy landscape to the mechanical properties of ductile glasses, J. Chem. Phys , 4593 (1999).
  * Maloney and Lemaître [2006] C. E. Maloney and A. Lemaître, Amorphous systems in athermal, quasistatic shear, Phys. Rev. E 74, 016118 (2006).
  * Eshelby [1957] J. D. Eshelby, The determination of the elastic field of an ellipsoidal inclusion, and related problems, Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 241, 376 (1957).
  * Lutsko [1989] J. F. Lutsko, Generalized expressions for the calculation of elastic constants by computer simulation, Journal of Applied Physics 65, 2991 (1989).
  * Gullett  _et al._ [2007] P. M. Gullett, M. F. Horstemeyer, M. I. Baskes, and H. Fang, A deformation gradient tensor and strain tensors for atomistic simulations, Modelling and Simulation in Materials Science and Engineering 16, 015001 (2007).
  * Halperin and Nelson [1978] B. Halperin and D. R. Nelson, Theory of two-dimensional melting, Physical Review Letters 41, 121 (1978).
  * Kosterlitz [2016] J. M. Kosterlitz, Kosterlitz–thouless physics: a review of key issues, Reports on Progress in Physics 79, 026001 (2016).
  * Braverman  _et al._ [2021] L. Braverman, C. Scheibner, B. VanSaders, and V. Vitelli, Topological defects in solids with odd elasticity, Phys. Rev. Lett. 127, 268001 (2021).