Further author information: (Send correspondence to J.R.M.)
J.R.M.: E-mail<EMAIL_ADDRESS>
# The conceptual design of GMagAO-X: visible wavelength high contrast imaging
with GMT

Jared R. Males  Steward Observatory, University of Arizona  Laird M. Close
Steward Observatory, University of Arizona  Sebastiaan Y. Haffert  Steward
Observatory, University of Arizona  Olivier Guyon  Steward Observatory,
University of Arizona  Victor Gasho  Steward Observatory, University of
Arizona  Fernando Coronado  Steward Observatory, University of Arizona
Olivier Durney  Steward Observatory, University of Arizona  Alexander Hedglen
James C. Wyant College of Optical Sciences, University of Arizona  Maggie
Kautz  James C. Wyant College of Optical Sciences, University of Arizona
Jamison Noenickx  Steward Observatory, University of Arizona  John Ford
Steward Observatory, University of Arizona  Tom Connors  Steward Observatory,
University of Arizona  Doug Kelly  Steward Observatory, University of Arizona

###### Abstract

We present the conceptual design of GMagAO-X, an extreme adaptive optics
system for the 25 m Giant Magellan Telescope (GMT). We are developing GMagAO-X
to be available at or shortly after first-light of the GMT, to enable early
high contrast exoplanet science in response to the Astro2020 recommendations.
A key science goal is the characterization of nearby potentially habitable
terrestrial worlds. GMagAO-Xis a woofer-tweeter system, with integrated
segment phasing control. The tweeter is a 21,000 actuator segmented deformable
mirror, composed of seven 3000 actuator segments. A multi-stage wavefront
sensing system provides for bootstrapping, phasing, and high order sensing.
The entire instrument is mounted in a rotator to provide gravity invariance.
After the main AO system, visible (g to y) and near-IR (Y to H) science
channels contain integrated coronagraphic wavefront control systems. The fully
corrected and, optionally, coronagraphically filtered beams will then be fed
to a suite of focal plane instrumentation including imagers and spectrographs.
This will include existing facility instruments at GMT via fiber feeds. To
assess the design we have developed an end-to-end frequency-domain modeling
framework for assessing the performance of GMagAO-X. The dynamics of the many
closed-loop feedback control systems are then modeled. Finally, we employ a
frequency-domain model of post-processing algorithms to analyze the final
post-processed sensitivity. The CoDR for GMagAO-X was held in September, 2021.
Here we present an overview of the science cases, instrument design, expected
performance, and concept of operations for GMagAO-X.

###### keywords:

adaptive optics

##  1 INTRODUCTION

The coming generation of 25-40 m Extremely Large Telescopes (ELTs) promises to
transform the study of extrasolar planets. The improved sensitivity due to
larger collection areas compared to today’s 5-10 m class telescopes will
enable highly sensitive radial velocity searches, as well as transit
spectroscopy on a range of nearby planets. Perhaps the area of exoplanet
science which stands to benefit the most is direct imaging. With adaptive
optics (AO) delivering diffraction limited resolution, the improvement in
sensitivity scales as diameter $D^{4}$. Just as importantly, the improvement
in spatial resolution with $\propto 1/D$ will allow direct imaging at smaller
separations than currently possible with today’s telescopes. If we equip them
with “extreme” AO (ExAO[1]) systems and coronagraphs, the ELTs will be able to
characterize large numbers of exoplanets. Most excitingly, such systems will
enable us to characterize large numbers of temperate, mature exoplanets for
the first time[2, 3]. Basic characterization (e.g. broadband albedo
measurement, molecule detection) of a few nearby exoplanets with such
characteristics may be possible with current telescopes in reflected visible
and near-IR light [4] and the thermal infrared[5]. But it is the ELTs which
will finally enable detailed studies of solar system like planets.

The most profound result of this transformative leap in capabilities will be
the ability to search for the signatures of life on planets orbiting other
stars. The Astro2020 Decadal Survey[6] established “Pathways to Habitable
Worlds” placing the search for life on planets other than Earth as a Priority
Area for astrophysics over the next decade. One of the four key capabilities
needed to achieve this is: Ground-based extremely large telescopes equipped
with high-resolution spectroscopy, high-performance adaptive optics, and high-
contrast imaging.[6]

The Giant Magellan Telescope (GMT) presents a unique opportunity to realize
the vision of Astro2020. With the smallest diameter of the 3 planned ELTs, it
can achieve higher Strehl with fewer deformable mirror (DM) actuators – making
it an ideal facility for short wavelength ExAO [7]. Furthermore, the large-
segment design of the GMT provides a straightforward way to harness existing
technologies to assemble the required DM. Finally, the current instrument
suite of the GMT[8] has space for an ExAO instrument, which could then be
coupled to the already in-progress high resolution spectrographs at GMT. Here
we present the results of the conceptual design study for GMagAO-X, an
instrument in development to fill this niche and be ready at or shortly after
first light of the GMT. At the conclusion of the conceptual design study,
GMagAO-X underwent a conceptual design review (CoDR) in September, 2021, with
an external panel, and is now in the preliminary design phase.

##  2 SCIENCE OBJECTIVES

The diffraction limited resolution resolution of the GMT at 500 nm (V band)
will be 4.1 mas, and at 1600 nm (H band) it will be 13.2 mas. These
unprecedented spatial resolutions, when coupled with $D^{4}$ exposure time
scaling, will make the GMT a truly transformative scientific instrument.
Taking advantage of this, GMagAO-X will enable high angular resolution
investigations of many compelling science cases. Here we list the top priority
science objectives identified for CoDR, all of which address the Astro2020
theme of “Worlds and Suns in Context” which “…encompasses the interlinked
studies of stars, planetary systems, and the solar system.”[6].

  1. 1.

Search For Life [new science]

     * •

Search for life on terrestrial exoplanets orbiting nearby stars

  2. 2.

Characterize Older Temperate Exoplanets [new science]

     * •

Reflected light characterization of wide range of radii and mass

  3. 3.

Measure Orbit and Mass of Exoplanets [today’s science, better]

     * •

Feed to G-CLEF for Precision-RV

     * •

Identify new candidates

     * •

Obtain precise ephemerides for imaging targets

  4. 4.

Study Planet Formation at Low Mass and Small Separation [today’s science,
better]

     * •

Search for and characterize thermally self-luminous young planets (YJH)

     * •

Characterize forming planets through the H-alpha accretion signature

  5. 5.

Circumstellar Disk Structure and Disk-Planet Interactions [today’s science,
better]

     * •

Structure in scattered light at high spatial resolution

     * •

Study disk-planet interactions

     * •

Jets and outflows

  6. 6.

Refine Stellar Evolution Models [today’s science, better]

     * •

High spatial and spectral resolution characterization of benchmark binaries

  7. 7.

High Spatial Resolution Mapping [today’s science, better]

     * •

Moons, asteroids, stellar surfaces

The key rationale for the conceptual design of GMagAO-X is that the first two
science objectives, which require reflected light characterization of
temperate exoplanets, are the most demanding. An ExAO instrument capable of
meeting these top science objectives will therefore be capable of meeting the
others. Placing our emphasis on these objectives is validated by Astro2020’s
“Pathways to Habitable Worlds”.

##  3 INSTRUMENT OVERVIEW

###  3.1 Optical and Mechanical Design

The optical and mechanical design of GMagAO-X is presented in detail in these
proceedings[9], as are the details of our wavefront sensing strategy[10], and
the plans for segment phasing control[11] and the 21,000 actuator “parallel
DM”[12]. In brief, the main components of GMagAO-X consist of

  * •

Woofer DM: ALPAO DM 3228 (64 actuator diameter) [9]

  * •

1st Stage WFS: 64x64 infrared Pyramid WFS [10]

  * •

Phase Control: 7x piezo T/T and Piston stages [9, 11]

  * •

Tweeter DM: seven 3000 actuator MEMS DMs[9, 12]

  * •

Coarse Phase Sensing: Holographic Dispersed Fringe Sensor[13, 14] (HDFS,
probably IR using 10% of J-H)

  * •

2nd Stage WFS & Fine Phase Sensing: 212x212 PWFS (I-J bands). Also considering
Zernike WFS [10]

  * •

Coronagraph NCP DM: 3000 actuator MEMS DM

  * •

Coronagraph: Lyot-architecture, supporting PIAACMC

  * •

Science: SDI with dual cameras, fiber feeds to other GMT instruments (namely
G-CLEF)

GMagAO-X is designed to be build-able today, with existing technology. All of
the above components can be purchased or fabricated today, making this
instrument well suited to the early-science phase at GMT.

Figure 1 shows an iso view of GMagAO-X, which will be an air-isolated two-
level floating optical table mounted on the Gregorian instrument rotator (GIR)
in a folded port (FP). The optical table system will be housed inside a barrel
which rotates to maintain gravity down as the GMT elevation changes during an
observation. We plan to use a deployable tertiary mirror, rather than the
facility M3, to allow much higher optical quality on a smaller optic. See
Close et al in these proceedings for a more complete description.[9]

Figure 1:  Iso view of GMagAO-X. See text for details.

###  3.2 Electronics

Shown in Figure 1, non-rotating electronics racks house control equipment for
the rotator, as well as the AO and coronagraph components. Deformable mirror
controllers will consume one of the racks shown. Computing equipment will be
housed off the telescope mount with low-latency fiber communications in
between.

The on-instrument electronic devices, such as filter wheels and focus stages,
are based on similar devices in MagAO-X[15, 16, 17].

###  3.3 Software

The instrument rotation control software, as well as safety and interlock
systems, will be based on standard GMT components. This will ensure that
GMagAO-X can be kept safe during any telescope operation when not observing.

Instrument control software will be based on the now well-tested MagAO-X
instrument control software. This includes the CACAO real-time software
system.[18]

##  4 AO ERROR BUDGET

Based on the chosen AO components and parameters, we first calculated a simple
error budget. The results are shown in Figure 2 for the WFS operating at 800
nm and science observations at 850 nm. For each guide star magnitude, the
d_opt column shows the optimum actuator spacing, which is an integer multiple
of 14 cm. This is to account for WFS detector binning. WFS exposure time is
optimized at each magnitude. The resulting wavefront error terms are shown as
a function of guide star magnitude in units of radians at 800 nm. Included are
measurement (WFS noise), time-delay, fitting (a function of d_opt),
chromaticity of scintillation affecting the OPD, chromaticity of the index of
refraction, dispersive anisoplanatism, and static non-common path errors. The
resulting Strehl ratio is shown in the last column.

Figure 3 is the same except for a science observing wavelength of 1600 nm.

Figure 2: Simple error budget for GMagAO-X for WFS operating at 800 nm and
observations at 850 nm.  Figure 3: Simple error budget for GMagAO-X for WFS
operating at 800 nm and observations at 1600 nm.

##  5 E2E PERFORMANCE MODELING

We next developed an end-to-end (E2E) model of GMagAO-X. We have been
developing a semi-analytic frequency domain framework for such analyses which
models coronagraph performance[19], and we have extended it to include focal
plane intensity dynamics and post-processing [20]. The model explicitly
includes modal gain optimization, and for this study we treated only frozen-
flow turbulence and the high-order AO loop. Figure 4 illustrates the E2E
model.

Figure 4: Illustration of the end-to-end semi-analytic model under development
for GMagAO-X. The model treats control system dynamics and produces the post-
coronagraph contrast[19] and intensity power spectrum[20] at any location in
the focal plane.

We used the E2E model to perform a trade study to compare our chosen 7x3k
“parallel DM”[9, 12] to a plausible monolithic DM with 128x128 actuators. The
7x3k has 190 effective actuators across. As shown in Figure 5 the fitting
error at 800 nm makes the 128x128 monolithic DM a secondary choice for
GMagAO-X.

We also compared using an EMCCD to a CMOS detector. The EMCCD was assigned
readout noise of 0.5 electrons, and the CMOS 2 electrons. The EMCCD excess
noise factor was included. The CMOS WFS performs better on $\sim 4$ to $\sim
10$ magnitude stars where the excess noise is important. The lower readnoise
of the EMCCD provides better performance on fainter stars.

Figure 5: Comparison between the 7x3k “Parallel DM” and a monolithic 128x128
actuator DM in terms of Strehl ratio at 800 nm. The fitting error term is
significant for the largest currently planned monolithic DM. Such a DM with
roughly 190 actuators across will be required to match the Parallel DM, and
such devices may be possible in the future but we do not expect them to be
available on the timescales of GMagAO-X.

We next computed the optimum Strehl (Figure 6) and contrast at 3 $\lambda/D$
(Figure 7). Here we compare the simple integrator (SI) and the Linear
Predictor (LP)[19]. The results depend modestly on whether Strehl or contrast
is used as the criterion for optimization. The predictive controller shows a
significant gain in post-coronagraph contrast as expected.

Figure 6: Predicted Strehl ratio when the end-to-end model of GMagAO-X is
optimized for best Strehl ratio and for best contrast at 3 $\lambda/D$. The
linear predictor (LP) controller has a significant impact over the simple
integrator (SI).  Figure 7: Predicted raw PSF contrast ratio when the end-to-
end model of GMagAO-X is optimized for Strehl ratio and for best contrast at 3
$\lambda/D$. The linear predictor (LP) controller has a significant impact
over the simple integrator (SI).

##  6 EXOPLANET YIELD

To assess the suitability of the conceptual design for meeting the science
objectives, we next analyzed the temperate exoplanet yield of GMagAO-X. This
was based on the catalog of known exoplanets maintained by NExScI. The catalog
was queried to identify exoplanet host stars within the declination range of
GMT. The parameters of these stars were first normalized using standard main
sequence relationships, filling in any missing values such as luminosity or
colors. Next, the exoplanets were analyzed for detectability. Minimum mass was
converted to mass assuming $\sin(i)=60^{o}$, and this was then converted to
radius using a mass-to-radius relationship we developed for use on these
lightly-irradiated planets. The Lambertian phase function was assumed. We
consider a range of models for geometric albedo $A_{g}$. For the results
presented here, we assume a constant $A_{g}=0.4$, which is more conservative
than the mean value of models such as the Earth, Venus, and gas-giant models
such as those of[21]. The planet:star flux ratio of each planet as a function
of its orbit was then computed. The output of the E2E model described above
was then converted to a photon noise and a speckle noise contribution, and the
signal-to-noise ratio as a function of time was computed. Here we assume that
quasi-static speckles are controlled by the focal-plane and low-order WFS
strategies designed into GMagAO-X sufficiently that the main contributor to
the photon and speckle noise terms is residual atmospheric turbulence. We
initially assume a coronagraph with an inner working angle (IWA) of
$1\lambda/D$, an assumption we will later relax.

Figure 8 displays the planet:star flux ratio vs separation for all of the
exoplanets we assess to be detectable by GMagAO-X. In Figure 9 we show the key
physical characteristics of these planets, radius vs. luminosity normalized
separation (a proxy for temperature). Most of the planets have planet:star
flux ratios $>10^{-8}$ (meaning they are brighter than this). The majority are
temperate gas-giants, and a range of host spectral types are represented.
There are $\sim$10 currently known terrestrial planets in the liquid water
habitable zones of their stars. It is these planets which we will search for
the signatures of life.

Figure 8: Planet:Star flux ratio vs. separation in $\lambda/D$ at 800 nm for
the GMagAO-X RV-detected planet sample.  Figure 9: Characteristics of the
GMagAO-X RV-detected planet sample. Against the y-axis we plot planet radius,
inferred by assuming $\sin(i)=60^{o}$ and using a custom mass-to-radius
conversion. The x-axis shows semi-major axis normalized by the square root of
the stellar luminosity. The gray region denotes the liquid water habitable
zone.

In Figure 10 we show the spatial distribution, in right ascension and
declination, of these targets. This is a key input to the potential concept of
operations for an instrument like GMagAO-X, since it motivates a queue mode
campaign where anytime conditions are favorable and a known planet hosting
star is near transit it can be observed. This will allow efficiently building
exposure time for detailed characterization of these exoplanets.

Figure 10: Spatial distribution of the GMagAO-X RV-detected planet sample.
There will be several planets observable on any given night.

Finally, we analyze the impact of performance worse than the optimistic
assumptions made. We compared coronagraph IWA of 1, 2, 3, 4, and 5
$\lambda/D$, along with post-coronagraph contrast worse by multiplicative
factors of 1, 10, 100, and 1000, respectively. Figure 11 shows the impact of
these degradations assuming a 20 msec atmospheric speckle lifetime[20].

Figure 11: Number of planets characterized with an albedo measurement for
various performance degradations assuming a 20 msec residual speckle lifetime.
Left shows the total number at each combination of degradations, and right
shown the number of planets vs. total exposure time for the survey.

We also consider post-processing, specifically using WFS telemetry. This will
lower the speckle lifetime sufficiently to make such observations photon-noise
dominated, making the effective residual speckle lifetime 0[20]. See also
Guyon et. al. in these proceedings[22]. The dramatic improvement in yield is
illustrated by comparing Figure 12 to Figure 11.

Figure 12: Number of planets characterized with an albedo measurement for
various performance degradations assuming a 0 msec effective residual speckle
lifetime. Left shows the total number at each combination of degradations, and
right shown the number of planets vs. total exposure time for the survey. We
expect to achieve this effective residual speckle lifetime through telemetry-
based post-processing[20]. Compare the dramatic improvement from Figure 11.

##  7 CONCLUSION

We have completed the conceptual design of GMagAO-X, a 21,000 actuator ExAO
coronagraph for the GMT, and successfully completed CoDR in September, 2021.
GMagAO-X takes advantage of GMT’s unique design, enabling rapid development
and deployment using currently available technologies. When on-sky at, or
shortly after, first-light of the GMT, GMagAO-X will enable the
characterization of hundreds of nearby, known-from-RV, temperate, mature
extrasolar planets. Roughly 10 of these are terrestrial planets in the liquid
water habitable zone of their host stars. Such planets are prime targets to
search for life to meet the challenge of Astro2020’s “Pathways to Habitable
Worlds” priority area.

###### Acknowledgements.

We are very grateful for support from the University of Arizona Space
Institute (UASI). This research has made use of the NASA Exoplanet Archive,
which is operated by the California Institute of Technology, under contract
with the National Aeronautics and Space Administration under the Exoplanet
Exploration Program.

## References

  * [1] Guyon, O., “Extreme Adaptive Optics,” ARA&A 56, 315–355 (Sept. 2018).
  * [2] Guyon, O., Martinache, F., Cady, E. J., Belikov, R., Balasubramanian, K., Wilson, D., Clergeon, C. S., and Mateen, M., “How ELTs will acquire the first spectra of rocky habitable planets,” in [Adaptive Optics Systems III ], Ellerbroek, B. L., Marchetti, E., and Véran, J.-P., eds., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 8447, 84471X (July 2012).
  * [3] Males, J. R., Close, L. M., Guyon, O., Morzinski, K., Puglisi, A., Hinz, P., Follette, K. B., Monnier, J. D., Tolls, V., Rodigas, T. J., Weinberger, A., Boss, A., Kopon, D., Wu, Y.-l., Esposito, S., Riccardi, A., Xompero, M., Briguglio, R., and Pinna, E., “Direct imaging of exoplanets in the habitable zone with adaptive optics,” in [Adaptive Optics Systems IV ], Marchetti, E., Close, L. M., and Vran, J.-P., eds., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 9148, 914820 (July 2014).
  * [4] Lovis, C., Snellen, I., Mouillet, D., Pepe, F., Wildi, F., Astudillo-Defru, N., Beuzit, J. L., Bonfils, X., Cheetham, A., Conod, U., Delfosse, X., Ehrenreich, D., Figueira, P., Forveille, T., Martins, J. H. C., Quanz, S. P., Santos, N. C., Schmid, H. M., Ségransan, D., and Udry, S., “Atmospheric characterization of Proxima b by coupling the SPHERE high-contrast imager to the ESPRESSO spectrograph,” A&A 599, A16 (Mar. 2017).
  * [5] Wagner, K., Boehle, A., Pathak, P., Kasper, M., Arsenault, R., Jakob, G., Käufl, U., Leveratto, S., Maire, A. L., Pantin, E., Siebenmorgen, R., Zins, G., Absil, O., Ageorges, N., Apai, D., Carlotti, A., Choquet, É., Delacroix, C., Dohlen, K., Duhoux, P., Forsberg, P., Fuenteseca, E., Gutruf, S., Guyon, O., Huby, E., Kampf, D., Karlsson, M., Kervella, P., Kirchbauer, J. P., Klupar, P., Kolb, J., Mawet, D., N’Diaye, M., Orban de Xivry, G., Quanz, S. P., Reutlinger, A., Ruane, G., Riquelme, M., Soenke, C., Sterzik, M., Vigan, A., and de Zeeuw, T., “Imaging low-mass planets within the habitable zone of $\alpha$ Centauri,” Nature Communications 12, 922 (Jan. 2021).
  * [6] National Academies of Sciences, E. and Medicine, [Pathways to Discovery in Astronomy and Astrophysics for the 2020s ] (2021).
  * [7] Close, L. M., Males, J. R., Follette, K. B., Hinz, P., Morzinski, K., Wu, Y.-L., Kopon, D., Riccardi, A., Esposito, S., Puglisi, A., Pinna, E., Xompero, M., Briguglio, R., and Quiros-Pacheco, F., “Into the blue: AO science with MagAO in the visible,” in [Adaptive Optics Systems IV ], Marchetti, E., Close, L. M., and Vran, J.-P., eds., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 9148, 91481M (Aug. 2014).
  * [8] Sitarski, B. N. et al., “Science instrumentation progress at the Giant Magellan Telescope,” Proc. SPIE 12184, #68 (2022).
  * [9] Close, L. M. et al., “The Optical and Mechanical Design for the 21,000 Actuator ExAO System for the Giant Magellan Telescope: GMagAO-X,” Proc. SPIE 12185, #76 (2022).
  * [10] Haffert, S. Y. et al., “Visible extreme adaptive optics for GMagAO-X with the triple-stage AO architecture (TSAO),” Proc. SPIE 12185, #516 (2022).
  * [11] Hedglen, A. et al., “First lab results of segment/petal phasing with a pyramid wavefront sensor and a novel holographic dispersed fringe sensor (HDFS) from the Giant Magellan Telescope high contrast adaptive optics phasing testbed,” Proc. SPIE 12185, #42 (2022).
  * [12] Kautz, A. et al., “A Novel Hexpyramid Pupil Slicer for an ExAO Parallel DM for the Giant Magellan Telescope,” Proc. SPIE 12185, #167 (2022).
  * [13] Haffert, S. Y., Close, L. M., Hedglen, A. D., Males, J. R., Kautz, M., Bouchez, A. H., Demers, R., Quirós-Pacheco, F., Sitarski, B. N., Gorkom, K. V., Long, J. D., Guyon, O., Schatz, L., Miller, K., Lumbres, J., Rodack, A., and Knight, J. M., “Phasing the Giant Magellan Telescope with the holographic dispersed fringe sensor,” Journal of Astronomical Telescopes, Instruments, and Systems 8(2), 1 – 21 (2022).
  * [14] Hedglen, A. D., Close, L. M., Haffert, S. Y., Males, J. R., Kautz, M. Y., Bouchez, A. H., Demers, R. T., Quirós-Pacheco, F., Sitarski, B. N., Guyon, O., Gorkom, K. V., Long, J. D., Lumbres, J. R., Schatz, L. H., Miller, K. L., Rodack, A. T., and Knight, J. M., “Lab tests of segment/petal phasing with a pyramid wavefront sensor and a holographic dispersed fringe sensor in turbulence with the Giant Magellan Telescope high contrast adaptive optics phasing testbed,” Journal of Astronomical Telescopes, Instruments, and Systems 8(2), 1 – 26 (2022).
  * [15] Males, J. R., Close, L. M., Miller, K., Schatz, L., Doelman, D., Lumbres, J., Snik, F., Rodack, A., Knight, J., Van Gorkom, K., Long, J. D., Hedglen, A., Kautz, M., Jovanovic, N., Morzinski, K., Guyon, O., Douglas, E., Follette, K. B., Lozi, J., Bohlman, C., Durney, O., Gasho, V., Hinz, P., Ireland, M., Jean, M., Keller, C., Kenworthy, M., Mazin, B., Noenickx, J., Alfred, D., Perez, K., Sanchez, A., Sauve, C., Weinberger, A., and Conrad, A., “MagAO-X: project status and first laboratory results,” in [Adaptive Optics Systems VI ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 10703, 1070309 (July 2018).
  * [16] Males, J. R., Close, L. M., Guyon, O., Hedglen, A. D., Van Gorkom, K., Long, J. D., Kautz, M., Lumbres, J., Schatz, L., Rodack, A., Miller, K., Doelman, D., Snik, F., Bos, S., Knight, J. M., Morzinski, K., Gasho, V., Keller, C., Haffert, S., and Pearce, L., “MagAO-X first light,” in [Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 11448, 114484L (Dec. 2020).
  * [17] Males, J. R. et al., “MagAO-X: current status and plans for Phase II,” Proc. SPIE 12185, #8 (2022).
  * [18] Guyon, O., Sevin, A., Gratadour, D., Bernard, J., Ltaief, H., Sukkari, D., Cetre, S., Skaf, N., Lozi, J., Martinache, F., Clergeon, C., Norris, B., Wong, A., and Males, J., “The compute and control for adaptive optics (CACAO) real-time control software package,” in [Adaptive Optics Systems VI ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 10703, 107031E (July 2018).
  * [19] Males, J. R. and Guyon, O., “Ground-based adaptive optics coronagraphic performance under closed-loop predictive control,” Journal of Astronomical Telescopes, Instruments, and Systems 4, 019001 (Jan. 2018).
  * [20] Males, J. R., Fitzgerald, M. P., Belikov, R., and Guyon, O., “The Mysterious Lives of Speckles. I. Residual Atmospheric Speckle Lifetimes in Ground-based Coronagraphs,” PASP 133, 104504 (Oct. 2021).
  * [21] Cahoy, K. L., Marley, M. S., and Fortney, J. J., “Exoplanet Albedo Spectra and Colors as a Function of Planet Phase, Separation, and Metallicity,” ApJ 724, 189–214 (Nov. 2010).
  * [22] Guyon, O. et al., “High Contrast Imaging at the Photon Noise Limit with WFS-based PSF Calibration,” Proc. SPIE 12185, #13 (2022).


# Enabling Resource-Aware Compilation of Spiking Neural Networks to
Neuromorphic Hardware

Adarsha Balaji, Shihao Song, Anup Das, Jeffrey Krichmar, Nikil Dutt,
James Shackleford, Nagarajan Kandasamy, and Francky Catthoor A. Balaji, S.
Song, A. Das, J. Shackleford, and N. Kandasamy are with the Department of
Electrical and Computer Engineering, Drexel University, Philadelphia, PA,
19104 USA e-mail<EMAIL_ADDRESS>Krichmar and
N. Dutt are with the Department of Computer Science, University of California,
Irvine, CA, USA.F. Catthoor is with Imec, Belgium and KU Leuven,
Belgium.Manuscript received Month Day, Year; revised Month Day, Year.

# Enabling Resource-Aware Mapping of Spiking Neural Networks via Spatial
Decomposition

Adarsha Balaji, Shihao Song, Anup Das, Jeffrey Krichmar, Nikil Dutt,
James Shackleford, Nagarajan Kandasamy, and Francky Catthoor A. Balaji, S.
Song, A. Das, J. Shackleford, and N. Kandasamy are with the Department of
Electrical and Computer Engineering, Drexel University, Philadelphia, PA,
19104 USA e-mail<EMAIL_ADDRESS>Krichmar and
N. Dutt are with the Department of Computer Science, University of California,
Irvine, CA, USA.F. Catthoor is with Imec, Belgium and KU Leuven,
Belgium.Manuscript received Month Day, Year; revised Month Day, Year.

###### Abstract

With growing model complexity, mapping Spiking Neural Network (SNN)-based
applications to tile-based neuromorphic hardware is becoming increasingly
challenging. This is because the synaptic storage resources on a tile, viz. a
crossbar, can accommodate only a fixed number of pre-synaptic connections per
post-synaptic neuron. For complex SNN models that have many pre-synaptic
connections per neuron, some connections may need to be pruned after training
to fit onto the tile resources, leading to a loss in model quality, e.g.,
accuracy. In this work, we propose a novel unrolling technique that decomposes
a neuron function with many pre-synaptic connections into a sequence of
homogeneous neural units, where each neural unit is a function computation
node, with two pre-synaptic connections.  This spatial decomposition technique
significantly improves crossbar utilization and retains all pre-synaptic
connections, resulting in no loss of the model quality derived from connection
pruning.  We integrate the proposed technique within an existing SNN mapping
framework and evaluate it using machine learning applications on the DYNAP-SE
state-of-the-art neuromorphic hardware. Our results demonstrate an average 60%
lower crossbar requirement, 9x higher synapse utilization, 62% lower wasted
energy on the hardware, and between 0.8% and 4.6% increase in model quality.

###### Index Terms:

Neuromorphic Computing, Spiking Neural Networks (SNNs), Machine Learning,
Computation Graph.

##  I Introduction

Spiking Neural Networks (SNNs) are machine learning approaches designed using
spike-based computation and bio-inspired learning algorithms [1]. In an SNN,
pre-synaptic neurons communicate information encoded in spike train to post-
synaptic neurons, via the synapses (see Fig. 1).  One of the key quality
metrics in SNNs is the inter-spike interval (ISI) defined as the inverse of
the mean firing rate of a neuron, which directly impacts model quality, e.g.,
accuracy.

Figure 1: An example of spiking neural network.

SNNs are executed on tile-based neuromorphic hardware such as DYNAP-SE [2]. A
tile consists of a crossbar [3], which is illustrated in Fig. 2. At each
crosspoint of a crossbar there is a synaptic element, which can be implemented
using Non-Volatile Memory (NVM) [4]. Beside neuromorphic computing, NVMs such
as Phase-Change Memory (PCM) are also used as main memory in conventional
computers to improve performance and energy efficiency [5, 6, 7, 8].

Within each crossbar, synaptic weights between pre- and post-synaptic neurons
are programmed as conductance of NVM cells. A pre-synaptic neuron’s voltage
(applied on a row) is multiplied by the conductance to generate current
(according to Ohm’s Law). Current summation occurs on each column according to
Kirchoff’s Current Law, when integrating excitation from its pre-synaptic
neurons. A neuron circuit takes as its input this current and generates at its
output a train of spikes. A spike is generated only if the current is higher
than a threshold. The spike firing frequency increases with input current,
saturating at a frequency determined by the refractory period of the neuron.

Figure 2: A crossbar in a neuromorphic hardware.

A $n\times n$ crossbar in a tile can accommodate only $n$ pre-synaptic
connections per post-synaptic neuron ($n=128$ for the crossbars in DYNAP-SE).
Fig. 3 reports the number of neurons with more than 128 pre-synaptic
connections, i.e., fanins as a fraction of the total number of neurons in a
few standard machine learning models, which are pruned iteratively to
eliminate the near-zero weights without loss in accuracy [9] (see our design
framework in Fig. 4).

Figure 3: Fraction of neurons with 128 or more pre-synaptic connections
(fanin).

We observe that, even after model pruning, on average 2.6% of neurons in these
models cannot be mapped to the crossbars in DYNAP-SE. There are two currently-
used solutions to this problem – 1) implement larger crossbars, which
increases the power consumption exponentially, and 2) remove synaptic
connections, which reduces the model quality.

Instead, we propose an unrolling technique which decomposes each neuron having
many pre-synaptic inputs or fanin into a sequence of homogeneous neural units,
where each neural unit is a function computation node having a maximum fanin-
of-two (or FIT). Our technique ensures information integrity as well as the
quality of these unrolled machine learning models. Furthermore, the unrolling
technique allows denser packing of the homogeneous neural units per crossbar,
significantly improving the crossbar utilization. We evaluate our technique
using standard machine learning applications and demonstrate an average 60%
lower resource usage and 62% lower wasted energy. The proposed spatial
decomposition technique retains all pre-synaptic connections for each neuron
in an SNN, improving model quality between 0.8% and 4.6% compared to an
existing SNN mapping approach.

##  II Proposed Technique

Fig. 4 illustrates the integration of the proposed technique (shown in the two
colored boxes) inside an existing design flow SpiNeMap [10] for mapping SNNs
to the hardware.111The proposed approach can be integrated with many other
mapping approaches such as the data flow-based performance-oriented SNN
compilation technique of [11, 12, 13, 14], the circuit aging oriented SNN
mapping technique of [15, 16, 17], endurance-aware mapping technique of [18,
19], and the SNN compiler of [20, 21]. This design flow incorporates 1)
Artificial Neural Network (ANN) models written in PyTorch [22] and Tensorflow
[23], and 2) SNN models written in PyCARL [24] and Brian [25]. For ANN models,
analog operations are first converted to spike-based operations using the
ANN2SNN converter [26].  These models are first pruned to eliminate near-zero
weights [9]. This is to keep the model size small for embedded platforms. The
pruned models are then presented to the proposed unrolling technique, which
compiles these models to crossbar-based hardware.

Figure 4: Integration of the proposed approach in an existing design flow.

Our proposed technique works in two steps – 1) the unrolling step, which
decomposes the SNN to limit the number of pre-synaptic connections for each
neuron, and 2) the weight normalization step, which ensures the quality of the
decomposed model. The latter is then clustered and mapped to the neuromorphic
hardware using SpiNeMap [10].

###  II-A Spatial Decomposition using Model Unrolling

We propose an unrolling approach, which decomposes a neuron function
computation with many fanins into a sequence of homogeneous neural units,
where each neural unit is a computation node with a maximum fanin-of-two
(FIT). Fig. 5 illustrates the decomposition (  ❷ ) of the neuron function
shown in  ❶ . Here, one $m$-input neuron function is decomposed into $(m-1)$
two-input neural units connected in sequence.

Figure 5: Unrolling a neuron functionality.

The neuron function $y_{o}=f\left(\sum_{i=1}^{m}n_{i}\cdot w_{i}\right)$ is
represented as

$\footnotesize y_{o}=f(u_{m-1}),\text{ where }u_{i}=\begin{cases}f(n_{1}\cdot
w_{1}+n_{2}\cdot w_{2}),\text{ for }i=1\\\ f(u_{i-1}+n_{i+1}\cdot
w_{i+1}),\text{ otherwise}.\end{cases}$ (1)

where $f$ represents the neuron functionality of generating spike trains with
a mean firing rate proportional to its input excitation,
$n_{1},n_{2},\cdots,n_{m}$ are the $m$ pre-synaptic neurons of the post-
synaptic neuron $n_{o}$, and $w_{1},w_{2},\cdots,w_{m}$ are the corresponding
synaptic weights.  The total number of FIT neural units generated from a
neural network with $N$ neurons is $\mathcal{N}=\sum_{i=1}^{N}(m_{i}-1)$,
where $m_{i}$ is the fanin of neuron $n_{i}$.

Figure 6: Demonstration of SNN clustering using SpiNeMap [10] (

❷

) and the proposed approach (

❹

).

###  II-B Weight Normalization

We apply weight normalization to optimize the synaptic weights $w_{i}$ of the
unrolled model (Fig. 5 ❶ ) such that the firing rate of a neuron in this model
is proportional to its input activation in the original model (Fig. 5 ❷ ). The
weight normalization is performed for each decomposed neural unit and the
weight normalization factor is applied to all its pre-synaptic weights. Using
Fig 5, the weight updates are

$\footnotesize w_{i}=\begin{cases}w_{i}/S_{\text{norm}}^{1},\text{ for
}i=1,2\\\ w_{i}/S_{\text{norm}}^{i-1},\text{
otherwise}.\end{cases}\text{~{}and~{}}w_{u_{i},u_{j}}=1/S_{\text{norm}}^{j}$
(2)

The normalization factor is computed as the maximum activation on the
corresponding synaptic weight in the original model using a batch from the
training set, i.e.,

$\footnotesize
S_{\text{norm}}^{i}=\begin{cases}k\cdot\max\\{a_{1}+a_{2}\\}\text{ for }i=1\\\
k\cdot\max\\{a_{i+1}\\}\text{ otherwise}.\end{cases},$ (3)

where $a_{i}$ is the activation on the synaptic weight $w_{i}$ in the original
model and the scaling factor $k$ is used to limit the mean firing rate of a
neuron to lower energy consumption on the hardware [26].  The weight
normalization overhead can be reduced by allowing non-uniform decomposition of
neuron functions. This is part of our future exploration.

###  II-C Motivating Example

Fig. 6 provides a motivating example demonstrating 1) how the proposed
technique results in a difference in the SNN mapping to hardware when compared
to SpiNeMap and 2) the improvements in model quality.

This example illustrates the mapping of three neuron functions $y_{1},y_{2}$
and $y_{3}$, shown in  ❶ to a neuromorphic hardware consisting of $4\times 4$
crossbars. SpiNeMap and other similar techniques will use three crossbars to
implement the three functions, as shown in  ❷ . Since a crossbar can
accommodate only a limited number of pre-synaptic connections per output
neuron (4 in this example), the component $x_{6}\cdot w_{g}$ of neuron
function $y_{2}$ cannot be mapped to the crossbar. This may result in a
degradation of the model quality.

In the proposed technique, neuron functions $y_{2}$, with fanin of 5 and
$y_{3}$, with fanin of 4 are decomposed to generate homogeneous neural
computation units $u_{1},u_{2},\cdots,u_{7}$ as shown in  ❸ . The neuron
function $y_{1}$ is not decomposed because the proposed technique unrolls only
those neuron functions that have fanin greater than 2.  During packing of
these units to a crossbar with a limited number of input ports, some of the
decomposed units may need to be combined to generate larger units. In this
example, the new computations are represented using functions
$f_{a},f_{b},f_{c}$ & $f_{d}$, where $f_{a}$ and $f_{c}$ can be implemented on
the first crossbar alongside $y_{1}$, while $f_{b}(y_{2})$ and $f_{d}(y_{3})$
are implemented on a second cluster (  ❹ ).  Finally, the two generated
clusters are shown in  ❺ .

We make the following three key observations. First, the proposed technique
only uses two crossbars to implement the three neuron functionalities
$y_{1},y_{2}$ and $y_{3}$, one less than that used by SpiNeMap. This reduces
hardware requirement and improves energy consumption. Second, it does not
eliminate any component of the neuron functionality, which improves the model
quality. Third, it increases the crossbar utilization.

##  III Results and Discussion

###  III-A Evaluation Methodology

We evaluated 10 machine learning applications that are representative of the
three most commonly used neural network classes — convolutional neural network
(CNN), multi-layer perceptron (MLP), and recurrent neural network (RNN).
Table I summarizes the original models and their baseline accuracy, which is
assessed from the inter-spike interval [10].

TABLE I: Applications used to evaluate our approach. Class | Applications | Synapses | Neurons | Topology | Accuracy
---|---|---|---|---|---
CNN | LeNet | 282,936 | 20,602 | CNN | 85.1%
AlexNet | 38,730,222 | 230,443 | CNN | 90.7%
VGG16 | 99,080,704 | 554,059 | CNN | 69.8 %
HeartClass [26, 27] | 1,049,249 | 153,730 | CNN | 63.7%
MLP | DigitRecogMLP | 79,400 | 884 | FeedForward (784, 100, 10) | 91.6%
EdgeDet [28] | 114,057 | 6,120 | FeedForward (4096, 1024, 1024, 1024) | 100%
ImgSmooth [28] | 9,025 | 4,096 | FeedForward (4096, 1024) | 100%
RNN | HeartEstm [29] | 66,406 | 166 | Recurrent Reservoir | 100%
VisualPursuit [30] | 163,880 | 205 | Recurrent Reservoir | 47.3%
DigitRecogSTDP [31] | 11,442 | 567 | Recurrent Reservoir | 83.6%

We evaluated these applications on the DYNAP-SE neuromorphic hardware [2],
where each tile has one $128\times 128$ crossbar. Table II reports the
hardware parameters of DYNAP-SE.

TABLE II: Major simulation parameters extracted from [2]. Neuron technology |  28nm FD-SOI
---|---
Synapse technology |  Hfo2-based OxRRAM
Supply voltage |  1.0V
Energy per spike |  50pJ at 30Hz spike frequency
Energy per routing |  147pJ
Switch bandwidth |  1.8G. Events/s

###  III-B Hardware Requirement

Figure 7 compares the number of crossbars required by SpiNeMap and the
proposed technique for each of the evaluated applications. Crossbar numbers
are normalized to SpiNeMap. We make the following four observations.

Figure 7: Crossbars needed for the evaluated applications.

First, on average, the proposed technique requires 60% fewer crossbars than
SpiNeMap to execute these applications. This reduction is because the
unrolling technique allows to densely pack the fanin synapses of different
neurons into the same crossbar, reducing the overall crossbar count. This
reduces the hardware energy consumption. Second, the number of crossbars
required in the proposed technique is 10x lower than SpiNeMap for EdgeDet,
even though EdgeDet has no neurons with more than 128 fanin synapses. This is
because many neurons in EdgeDet have fanin close to the maximum limit that a
crossbar in DYNAP-SE can accommodate. Therefore, these high-fanin neurons
cannot be packed in the same crossbar by SpiNeMap, needing a separate crossbar
for each of these neurons. The proposed technique, on the other hand, unrolls
a high-fanin neuron to create a structure with a maximum-fanin-of-two. This
results in a denser packing of crossbars and a correspondingly lower crossbar
requirement. Third, for CNN-based applications (LeNet, AlexNet, VGG16, and
HeartClass), the proposed technique requires 37% fewer crossbars than
SpiNeMap. We observe that the reduction over SpiNeMap is higher for networks
with more layers; the reduction is 16% for LeNet, compared to 48% for VGG16.
This is because, with more layers, the proposed technique has more freedom to
improve crossbar utilization (see Sec. III-C). Finally, for RNN-based
applications (HeartEstm, VisualPursuit, and DigitRecogSTDP), the proposed
approach requires 80% fewer crossbars than SpiNeMap. This is due to the cyclic
nature of connections between neurons in these applications. SpiNeMap cannot
optimally pack neurons in cyclic connections to the crossbars. The proposed
technique decomposes a cyclic connection to pack its neurons densely into
crossbars, reducing the crossbar requirement.

Figure 8: Neuron utilization.

Figure 9: Synapse utilization.

###  III-C Crossbar Utilization

Figures 8 and 9 compare respectively, the neuron and synapse utilization of
SpiNeMap and the proposed technique on DYNAP-SE for each evaluated
application. We observe that the average neuron utilization of SpiNeMap is 53%
and of the proposed technique is 62% for these applications. The average
synapse utilization of SpiNeMap is only 9% and of the proposed technique is
25% for these applications. The reason for the high crossbar utilization is
because of the proposed spatial decomposition technique, which allows to
densely pack fanin and fanout synapses from different neurons in the same
crossbar, improving utilization. Higher utilization leads to lower wasted
energy as reported next.

###  III-D Wasted Energy

Figure 10 reports the energy wasted for each evaluated application on DYNAP-SE
using the proposed technique, normalized to SpiNeMap. The wasted energy
incorporates the neurons and synapses in each crossbar that are not utilized
during the execution of these applications. We observe that the energy wasted
using the proposed technique is on average 62% lower than SpiNeMap. This
significant reduction is due to 1) the reduction in the use of crossbars and
2) the increase in the utilization of neurons and synapses in each crossbar.

Figure 10: Wasted energy on DYNAP-SE.

###  III-E Model Quality

Figure 11 plots the model quality i.e., the accuracy of each evaluated
application on DYNAP-SE using SpiNeMap and the proposed technique. For
comparison, the quality using software simulation is also reported in the
figure as Baseline. We make the following two key observations.

Figure 11: Model quality for the evaluated applications.

First, the accuracy loss for some applications such as ImgSmooth, EdgeDet,
HeartEstm, and VisualPursuit are the same for SpiNeMap and the proposed
technique. This is because no neuron in these applications has more than 128
fanin synapses (see Fig. 3) and therefore, no pre-synaptic connections are
eliminated by SpiNeMap. So, the quality of the two techniques are the same.
However, for these applications, the proposed technique is significantly
better than SpiNeMap in terms of crossbar usage, their utilization, and the
wasted energy. Second, for all other applications that do not initially fit on
the available crossbar space, the quality of the proposed technique is better
than SpiNeMap by average 3% (between 0.8% and 4.6%). For these applications,
the proposed technique is better in terms of all the evaluated metrics.

##  IV Discussion

We note that, the decomposition technique proposed here negatively impacts
accuracy and energy in following two aspects. First, by decomposing a large
neuron function, the proposed technique maps some of the decomposed synaptic
connections on the shared interconnect of a neuromorphic hardware, negatively
impacting spike latency. This can lower the model accuracy. However, the
SpiNeMap framework, which integrates the decomposition technique minimizes
such impact by intelligent cluster mapping and placement on the hardware [10].
Furthermore, by not eliminating any synaptic connections, the proposed
technique, in fact, improves accuracy compared to SpiNeMap (see Sec. III-E).
Second, with additional synapses mapped to the shared interconnect, the energy
consumption on the interconnect increases. However, we find this increase in
energy is much lower than the energy savings obtained by reducing the crossbar
usage (see Sec. III-D).

##  V Conclusion

We present a technique to map SNN-based applications to crossbar-based
neuromorphic hardware. The proposed technique involves unrolling of neurons,
which decomposes a complex neuron functionality into a sequence of homogeneous
neural units, where each neural unit is a fanin-of-two (FIT) neuron.  The
unrolling technique significantly improves crossbar utilization and ensures
information integrity, resulting in no loss of model quality derived from
connection pruning.  We integrate this unrolling technique inside an existing
SNN mapping framework and evaluate it using machine learning applications for
a state-of-the-art neuromorphic hardware. Our results demonstrate an average
60% lower crossbar requirement, 9x higher synapse utilization, 62% lower
wasted energy, and 3% increase in model quality compared to an existing SNN
mapping approach. In the future, we will explore the trade-offs involved in
non-uniform tree decomposition of neural function.

## Acknowledgment

This work is supported by the National Science Foundation Award CCF-1937419
(RTML: Small: Design of System Software to Facilitate Real-Time Neuromorphic
Computing).

## References

  * [1] W. Maass, “Networks of spiking neurons: The third generation of neural network models,” _Neural Networks_ , 1997.
  * [2] S. Moradi _et al._ , “A scalable multicore architecture with heterogeneous memory structures for dynamic neuromorphic asynchronous processors (DYNAPs),” _TBCAS_ , 2017.
  * [3] F. Catthoor _et al._ , “Very large-scale neuromorphic systems for biological signal processing,” in _CMOS Circuits for Biological Sensing and Processing_ , 2018.
  * [4] A. Mallik _et al._ , “Design-technology co-optimization for OxRRAM-based synaptic processing unit,” in _VLSIT_ , 2017.
  * [5] S. Song _et al._ , “Enabling and exploiting partition-level parallelism (PALP) in phase change memories,” _TECS_ , 2019.
  * [6] S. Song _et al._ , “Exploiting inter-and intra-memory asymmetries for data mapping in hybrid tiered-memories,” in _ISMM_ , 2020.
  * [7] S. Song _et al._ , “Improving phase change memory performance with data content aware access,” in _ISMM_ , 2020.
  * [8] S. Song _et al._ , “Aging aware request scheduling for non-volatile main memory,” in _ASP-DAC_ , 2021.
  * [9] S. Han _et al._ , “Deep compression: Compressing deep neural networks with pruning, trained quantization and huffman coding,” _arXiv_ , 2015.
  * [10] A. Balaji _et al._ , “Mapping spiking neural networks to neuromorphic hardware,” _TVLSI_ , 2020.
  * [11] S. Song _et al._ , “Compiling spiking neural networks to neuromorphic hardware,” in _LCTES_ , 2020.
  * [12] A. Das _et al._ , “Mapping of local and global synapses on spiking neuromorphic hardware,” in _DATE_ , 2018.
  * [13] A. Balaji _et al._ , “A framework for the analysis of throughput-constraints of snns on neuromorphic hardware,” in _ISVLSI_ , 2019\.
  * [14] A. Das _et al._ , “Dataflow-based mapping of spiking neural networks on neuromorphic hardware,” in _GLSVLSI_ , 2018.
  * [15] S. Song _et al._ , “Improving dependability of neuromorphic computing with non-volatile memory,” in _EDCC_ , 2020.
  * [16] A. Balaji _et al._ , “A framework to explore workload-specific performance and lifetime trade-offs in neuromorphic computing,” _CAL_ , 2019\.
  * [17] S. Song _et al._ , “A case for lifetime reliability-aware neuromorphic computing,” in _MWSCAS_ , 2020.
  * [18] T. Titirsha _et al._ , “Thermal-aware compilation of spiking neural networks to neuromorphic hardware,” in _LCPC_ , 2020.
  * [19] T. Titirsha _et al._ , “Reliability-performance trade-offs in neuromorphic computing,” in _CUT_ , 2020.
  * [20] Y. Ji _et al._ , “Bridge the gap between neural networks and neuromorphic hardware with a neural network compiler,” in _ASPLOS_ , 2018.
  * [21] A. Balaji _et al._ , “Run-time mapping of spiking neural networks to neuromorphic hardware,” _JSPS_ , 2020.
  * [22] A. Paszke _et al._ , “Pytorch: An imperative style, high-performance deep learning library,” in _NeurIPS_ , 2019.
  * [23] M. Abadi _et al._ , “Tensorflow: A system for large-scale machine learning,” in _OSDI_ , 2016.
  * [24] A. Balaji _et al._ , “PyCARL: A PyNN interface for hardware-software co-simulation of spiking neural network,” in _IJCNN_ , 2020.
  * [25] D. F. Goodman _et al._ , “The brian simulator,” _Front. in Neurosc._ , 2009.
  * [26] A. Balaji _et al._ , “Power-accuracy trade-offs for heartbeat classification on neural networks hardware,” _JOLPE_ , 2018.
  * [27] A. Das _et al._ , “Heartbeat classification in wearables using multi-layer perceptron and time-frequency joint distribution of ECG,” in _CHASE_ , 2018\.
  * [28] T. S. Chou _et al._ , “CARLsim 4: An open source library for large scale, biologically detailed spiking neural network simulation using heterogeneous clusters,” in _IJCNN_ , 2018.
  * [29] A. Das _et al._ , “Unsupervised heart-rate estimation in wearables with Liquid states and a probabilistic readout,” _Neural Networks_ , 2018.
  * [30] H. J. Kashyap _et al._ , “A recurrent neural network based model of predictive smooth pursuit eye movement in primates,” in _IJCNN_ , 2018.
  * [31] P. U. Diehl _et al._ , “Unsupervised learning of digit recognition using spike-timing-dependent plasticity,” _Front. in Comp. Neuroscience_ , 2015\.


# Athermal quasistatic cavitation in amorphous solids: effect of random
pinning

Umang A. Dattani<EMAIL_ADDRESS>The Institute of Mathematical Sciences,
C.I.T. Campus, Taramani, Chennai 600113, India  Homi Bhabha National
Institute, Anushakti Nagar, Mumbai, 400094, India  Smarajit Karmakar
<EMAIL_ADDRESS>Tata Institute of Fundamental Research, 36/P, Gopanpally
Village, Serilingampally Mandal,Ranga Reddy District, Hyderabad, 500046,
Telangana, India  Pinaki Chaudhuri<EMAIL_ADDRESS>The Institute of
Mathematical Sciences, C.I.T. Campus, Taramani, Chennai 600113, India  Homi
Bhabha National Institute, Anushakti Nagar, Mumbai, 400094, India

(August 28, 2024)

###### Abstract

Amorphous solids are known to fail catastrophically via fracture, wherein
cavitation at nano-metric scales is known to play a significant role. Micro-
alloying via inclusions is often used as a means to increase the fracture
toughness of amorphous solids. Modeling such inclusions as randomly pinned
particles that move only affinely and do not participate in plastic
relaxation, we study how the pinning influences the process of cavitation-
driven fracture in an amorphous solid. Using extensive numerical simulations
and probing in the athermal quasistatic limit, we show that just by pinning a
very small fraction of particles, the tensile strength is increased and also
the cavitation is delayed. Further, the cavitation that is expected to be
spatially heterogeneous becomes spatially homogeneous by forming a large
number of small cavities instead of a dominant cavity.

##  I Introduction

The mechanical properties of amorphous solids are utilized in diverse
applications in industries and our daily lives Schuh  _et al._ (2007); Bonn
_et al._ (2017). Therefore, mechanical failure of these materials is an area
of concern. Hence, understanding the physical processes that lead to the
failure of these structurally disordered solids is a domain of current active
research, with the primary goal being to figure out design pathways that can
sustain against such failures Rodney  _et al._ (2011); Liu  _et al._ (2021).
Cavitation, i.e., the formation of nano-cavities within the solid, has been
identified as a precursor to eventual failure via fracture Bouchaud  _et al._
(2008); Murali  _et al._ (2011). In experiments, the fracture in amorphous
solids has been shown to propagate via the coalescence of cavities in the
solid along the direction of the crack Shen  _et al._ (2021), which has
motivated numerical investigations to analyze the mechanisms underlying the
cavitation process Guan  _et al._ (2013); Chaudhuri and Horbach (2016);
Altabet  _et al._ (2018a, b); Paul  _et al._ (2020). Recently Dattani  _et
al._ (2022), it has been demonstrated that the plasticity associated with
cavitation has the same universal characteristic elastoplastic response of
amorphous solids undergoing failure via cavitation under uniform
expansionDattani  _et al._ (2022) and on exploring a combination of loading
scenarios where cavitation can occurDattani  _et al._ (2023). In Ref.Dattani
_et al._ (2023), we demonstrated that a combination of deformation processes,
for example, uniform expansion followed by oscillatory shear of a certain
amplitude, can enhance the formation of cavities at densities that are much
higher than the cavitation density observed in uniform expansion only. These
results suggest that in a natural deformation process in which various forms
of deformations will be coupled together, the material can show unpredictable
failure behavior, which will be difficult to control. Thus a systematic study
of how such a cavity-dominated failure process in an amorphous solid can be
mitigated effectively will be of significant interest for practical
applications.

In the quest of making glasses with higher fracture-toughness, the seeding of
the amorphous solids with micro-alloyed inclusions has gained a lot of
popularity in the last few years Hofmann  _et al._ (2008); Hardin and Homer
(2015); Torquato and Haslach Jr (2002); Tyukodi  _et al._ (2016); Garrett  _et
al._ (2012); González (2016). In the numerical modeling, a minimal model of
these micro-alloyed inclusions describes these inclusions as pinned/frozen
particles Gendelman  _et al._ (2014); Bhowmik  _et al._ (2019a, b); Liu  _et
al._ (2023). In the context of probing their mechanical behavior, it is
assumed that they only move affinely during the deformation and do not
actually undergo non-affine motion. Using such a model system, systematic
studies investigating the response of a pinned amorphous solid to a shear-
deformation have explored the microscopic theoriesGendelman  _et al._ (2014),
yielding mechanismsBhowmik  _et al._ (2019a), suppression of shear-banding
Bhowmik  _et al._ (2019b), development of intrinsic length-scales Bhowmik  _et
al._ (2019a); Liu  _et al._ (2023) etc. As these studies focus on the shear-
deformation of high-density amorphous solids, they do not access the region
where cavitation instabilities occur, i.e., under axial tension Altabet  _et
al._ (2018a); Dattani  _et al._ (2022, 2023). Due to the importance of
cavitation instabilities in fracture, it becomes important to study the
response of an amorphous solid under deformation modes where cavitation can
occur.

In this work, we, therefore, study the response of amorphous solids with
micro-alloyed inclusions modeled as pinned particles under uniform expansive
deformation under athermal quasistaic conditions. We find that, even in
presence of a very small fraction of pinned particles, cavitation becomes more
spatially homogeneous and the cavitation point shift to lower densities and
lower pressures, implying a higher load-bearing capacity of the pinned solid.
The sharp brittle-yielding-like transition seen in unpinned solids becomes
more gradual with significantly smaller sizes of plastic events due to which
system size dependence becomes very weak. On tracking the eigenvalues of the
Hessian near cavitation instabilities, we find that, on the potential energy
landscape, cavitation occurs via a saddle-node bifurcation and the average
spatial decay of displacements in the plastic eigenmodes from the plastic
center reveal a length-scale in pinned solids that parallels the length-scale
set by average distance between two pinned sites. The presence of a length-
scale explains the absence of system-size effects and the drastic decrease in
the mean sizes of plastic events. Our findings thus reveal how micro-alloyed
inclusions can suppress the cavitation and how the presence of a length-scale
of plasticity controls the deformation response of such a pinned solid.

The manuscript is organised as follows. After initial introductory discussion
in Section I, we provide in Section II a brief overview of the model amorphous
solid that we consider for our study and the methodology of our simulations as
well as analysis. In Section III, we discuss the detailed findings of our
investigations regarding the presence of random pinning and its influence on
the cavitation process. Finally, we provide a concluding discussion in Section
IV.

##  II Model and Methods

###  II.1 Model details and initial states

We use the well-characterized two-dimensional model consisting of two species
labelled $A$ and $B$ at $65:35$ concentration ratio, interacting via pairwise
Lennard-Jones potential. The interaction parameters are – $\sigma_{AA}=1.0$,
$\sigma_{BB}=0.88$, $\sigma_{AB}=0.8$, $\epsilon_{AA}=1.0$,
$\epsilon_{BB}=0.5$, $\epsilon_{AB}=1.5$ Brüning  _et al._ (2008). With this
model, we smoothen the interaction potential up to first two derivatives. The
form of the interactions between $i^{th}$ and $j^{th}$ particle becomes:

$\phi(r_{ij})=4\epsilon_{\alpha\beta}\left[\left(\frac{\sigma_{\alpha\beta}}{r_{ij}}\right)^{12}-\left(\frac{\sigma_{\alpha\beta}}{r_{ij}}\right)^{6}\right]+u(r_{ij})$
(1)

where,

$u(r_{ij})=C_{0}+C_{2}\left(\frac{r_{ij}}{\sigma_{\alpha\beta}}\right)^{2}+C_{4}\left(\frac{r_{ij}}{\sigma_{\alpha\beta}}\right)^{4}$
(2)

Here, $\alpha$ and $\beta$ correspond to either of the labels $A$ or $B$. The
constants $C_{0}$, $C_{2}$ and $C_{4}$ are determined by requiring the
potential and its first two derivatives to be zero at the cutoff
$r=2.5\sigma_{ij}$. The simulations have been performed for a variety of
system sizes ranging from $N=10^{3}$ to $N=10^{5}$.

###  II.2 Initial states

To prepare initial states for our study, we first equilibriate the system at
$T=1.0$ (in LJ units), which is in the liquid regime, followed by cooling at a
constant rate of $10^{-4}$ per MD timestep to a final temperature of $T=0.01$
Bhowmik  _et al._ (2019a), which is in the glassy regime. The corresponding
glass transition temperature of the model system is at $T=0.44$ Brüning  _et
al._ (2008). The athermal states used in our study are generated by obtaining
inherent structure states corresponding to the glassy configurations at
$T=0.01$, via conjugate gradient (CG) minimizationPolak and Ribiere (1969).

###  II.3 Athermal Quasistatic Expansion

Starting from a spatially homogeneous high density state ($\rho=1.2$ for
KABLJ) having positive barostatic pressure, we study the athermal quasi-static
response (i.e. in the absence of any thermal effects and in the limit of
vanishing driving rates) of this system to isotropic expansion Altabet  _et
al._ (2018a). In each expansion step, a constant volume strain is applied on
the system by rescaling the length of the box by a factor $(1+\epsilon)$ along
with affine transformation of particle coordinates, followed by minimization
of the energy of this strained configuration using the conjugate gradient
algorithm Polak and Ribiere (1969). The values of $\epsilon$ are varied from
$\epsilon=10^{-4}$ to $\epsilon=10^{-9}$. The AQE simulations are done using
LAMMPSPlimpton (1995).

###  II.4 Pinning

We choose a small fraction of particles $c=0.01$ to $c=0.05$ in the generated
solid and freeze their motion. The particles are chosen randomly as long as no
two pinned particles lie within the cut off of the interaction between
themBhowmik  _et al._ (2019a). This helps avoiding the scenario where two
close-by pinned sites increase the energy of the system. The pinned particles
only move affinely when the strain is applied. During the energy minimization,
these pinned particles are not allowed to move.

###  II.5 Hessian of potential energy

LAPACKEAnderson  _et al._ (1999) is used for doing the stability analysis of
the local minima states, by computing eigenvalues and eigenvectors of the
Hessian matrix $\mathcal{H}^{\alpha\beta}_{ij}$, which is defined as

$\mathcal{H}_{ij}^{\alpha\beta}=\frac{\partial^{2}U\left(\\{\mathbf{r}_{i}\\}\right)}{\partial
r_{i}^{\alpha}\partial r_{j}^{\beta}},$ (3)

where $U\left(\\{\mathbf{r}_{i}\\}\right)$ is the potential energy of the
system and $\mathbf{r}_{i}$ is the position vector of particle $i$. The
indices $\alpha,\beta\in\\{x,y\\}$ whereas $i,j\in\\{1,\ldots,N\\}$.

If we now consider a system of $N$ particles, where particle numbers,
$i=1,\cdots,m$ are free and particle numbers $i=m+1,\cdots,N$ are pinned, then
the potential energy of such a system can be expressed as

$U(r)=\frac{1}{2}\left[\sum_{i,j=1;i\neq
j}^{m}\phi_{ij}+2\cdot\sum_{i=1}^{m}\sum_{j=m+1}^{N}\phi_{ij}\right],$ (4)

where the first term comes from the interactions between unpinned particles,
the second term comes due to the interactions between the pinned and the
unpinned particles. Note that the term due to interactions between pinned
sites is set to zero because of our pinning protocol

By substituting Eq.(4) in Eq.(3) , it can be shown that the first term in the
sum of Eq.(4) gives a contribution,

$\small H^{ij}_{\alpha\beta}=-\sum_{k,i;k\neq
i}\left[\left(\frac{\phi_{r}}{(r^{ik})^{3}}-\frac{\phi_{rr}}{(r^{ki})^{2}}\right)r^{ki}_{\alpha}r^{ki}_{\beta}-\delta_{\alpha\beta}\frac{\phi_{r}}{r^{ki}}\right](\delta^{ji}-\delta^{jk}).$
(5)

The second term of the sum in Eq.(4) gives,

$H^{ij}_{\alpha\beta}=-2\,\sum_{k=0}^{m}\,\sum_{l=m+1}^{N}\left[\frac{\phi_{r}\,r^{kl}_{\alpha}\,r^{kl}_{\beta}}{(r^{kl})^{3}}-\frac{\phi_{rr}\,r^{kl}_{\alpha}\,r^{kl}_{\beta}}{(r^{kl})^{2}}-\delta_{\alpha\beta}\frac{\phi_{r}}{r^{kl}}\right]\delta^{ij},$
(6)

where $\phi_{r}$ and $\phi_{rr}$ are first and second derivatives of the pair-
potential with respect to variable $r$ respectively.

Figure 1: For $N={10}^{5}$ and different pinning concentrations, $c$, as
marked: Variation with density $\rho$ of (a) Pressure $P$ (b) Energy per
particle $U/N$ (c) Pressure susceptibility $\chi_{P}$ (d) Energy
susceptibility $\chi_{U}$. Change in (e) average size of pressure jump
$\Delta{P}$, and (f) average size of energy drops $\Delta{U}$, with pinning
concentration.  Figure 2: Density field for different pinning concentrations,
$c=$ 0.00 (a), 0.01 (b), 0.02 (c), 0.05 (d), measured at $\rho=0.982$. Figure
3: (Top) Variation of ensemble averaged pressure ($P$) with density ($\rho$)
for different system sizes, as labelled, at $c=$ 0.00 (a), 0.01 (b), 0.02 (c),
0.05. (d) (Bottom) Corresponding variation of pressure susceptibility
$\chi_{p}$ with density $\rho$, in each case. Figure 4: (a) Pressure $P$ vs
density $\rho$ for an expansion trajectory, using $N=2500$, corresponding to
$c=0.05$. (b) Demonstration of square-root singularity for lowest eignevalue
of the Hessian, $\lambda_{\rm min}$ (appropriately scaled by fit parameter
$A$), computed at the density points in (a) corresponding to occurrence of
plastic instabilities; $\rho_{c}$ is the estimate density location of the
event in each case. (c) Eigenmode, (d) displacement field just before the drop
and (e) displacement fields across the plastic drop, occurring at one such
plastic event, near $\rho\approx 1.03277$. (f)-(j) Evolution of the density
field across the expansion trajectory shown in (a).

##  III Results

###  III.1 Yielding and spatial ramifications

Upon expanding the amorphous solid isotropically under quasistatic loading, as
discussed in previous worksAltabet  _et al._ (2018a); Dattani  _et al._
(2022), the pressure of the solid decreases monotonically, eventually reaching
negative values. After a certain threshold, a sharp jump in the pressure
accompanies the cavitation of the solid. Upon further expansion, the cavities
grow and merge, leading to system-spanning fracture of the solid Dattani  _et
al._ (2022). Here, we expand the same amorphous solid isotropically but with a
small fraction of particles pinned(frozen) and compare it with the case
without pinningAltabet  _et al._ (2018a); Dattani  _et al._ (2022). As shown
in Fig. 1(a), compared to the unpinned solid, the pinned solid, on average,
does not show a large pressure jump that is seen around the first cavitation
event. The location of the yield point, which is usually marked by a turn in
the pressure-density curves, shifts to lower and lower densities with
increasing concentration of pinned particles, $c$, demonstrating an increase
in the load-bearing capacity of pinned solids. In the pinned solids ($c\neq
0$), the turning of the curve occurs due to small pressure jumps that occur
gradually, leading to a smooth average pressure vs density ($P-\rho$) curve as
opposed to the ensemble-averaged trajectories of an unpinned solid ($c=0$),
which show an abrupt jump. A similar trend reflects in the per particle energy
vs density ($U/N-\rho$) plots in Fig. 1(b) i.e. no large energy drop shows up
in the pinned solid unlike the case of unpinned solid and in fact, the energy
per particle of the pinned solid keeps increasing with increasing pinning
concentration. The nature of critical-like behaviour in the ($P-\rho$) or
($U/N-\rho$) plots can be studied by measuring the fluctuations of pressure
and energy at a given value of density across ensembles as characterized by
the following susceptibilitiesDattani  _et al._ (2022),

$\displaystyle\chi_{p}(\rho)$ $\displaystyle=$ $\displaystyle N\,\left(\langle
P^{2}(\rho)\rangle-\langle P(\rho)\rangle^{2}\right),\quad\mbox{and}$ (7)
$\displaystyle\chi_{u}(\rho)$ $\displaystyle=$
$\displaystyle(1/N)\,\left(\langle U^{2}(\rho)\rangle-\langle
U(\rho)\rangle^{2}\right)$ (8)

and shown in Fig. 1(c) & (d) respectively. A sharp susceptibility peak around
the big pressure jump / energy drop seen in the unpinned solid Dattani  _et
al._ (2022) is tamed down due to pinning. For higher pinning concentrations, a
clear peak around the yielding is not seen implying that the yielding of the
solid is more localised and gradual. The gradual nature of cavitation in
pinned solids is also echoed in the average size of pressure jumps
$\langle\Delta P\rangle$ and energy drops $\langle\Delta U\rangle$ encountered
during expansion (see 1(e) & (f)) for different pinning concentrations $c$.
The average size of avalanches, $\langle\Delta P\rangle$ & $\langle\Delta
U\rangle$ decreases drastically with increasing values of $c$ implying
suppression of large pressure jumps and energy drops.

To probe the spatial-nature of cavitation in pinned solids, we look at the
coarse-grained spatial density maps for different values of $c$ at a same
value of density ( shown in Fig. 2). They suggest that the cavitation occurs
at multiple sites in the solid gradually over the course of expansion for
$c\neq 0$ instead of a heterogeneous cavitation starting with a big cavity for
$c=0$. The creation of a large number of cavities in the system also explains
the increase in energy with expansion seen in Fig. 1(b) as presence of a large
number of particles on the surfaces is expected to increase the energy of
pinned solids. So, to summarise, pinning smooths out the sharp brittle
yielding-like cavitation transition seen in amorphous solids without any
micro-alloying and causes less heterogeneous cavitation. This scenario is
consistent with a previous study on effect of pinning on yielding transition
under simple shearBhowmik  _et al._ (2019a) where pinning made the yielding
more spatially homogeneous.

###  III.2 System-size effects

In Ref. Dattani  _et al._ (2022), strong system size effects were observed in
the $P-\rho$ the $\chi_{p}$ curves for athermal quasistatic expansion of
amorphous solid. Hence it is important to probe the dependence on system sizes
for pinned solids as well. The $P-\rho$ plots and $\chi_{p}-\rho$ plots for
different system sizes across the pinning concentrations are shown in Fig.3
top panel and bottom panels, respectively. Unlike the case of unpinned solids,
the $P-\rho$ and $\chi_{p}-\rho$ curves for pinned solids show little to no
dependence on the system size. This occurs due to the emergence of an
intrinsic length scale of plasticity $\xi<<L$ in these systems due to the
imposed random pinning constraints, which are discussed in the subsequent
paragraphs.

###  III.3 Irreversible plastic events on the potential energy landscape

Under athermal quasistatic shear, near a plastic instability, the lowest non-
zero eigenvalue (apart from the two zero-modes in 2D) of the hessian matrix of
the potential energy Eq. 3 is known to vanish as a square root of strain
difference from the point of instability Malandro and Lacks (1999); Maloney
and Lemaître (2006). i.e. $\lambda_{min}\sim\sqrt{\gamma_{c}-\gamma}$, where
$\lambda_{min}$ refers to the minimum eigenvalue. This occurs due to a saddle-
node bifurcation on a potential energy landscape (PEL) where the local minima
in which the system resides becomes unstable in one direction. One of the
possible ways to arrive at such a square-root power law has been discussed in
Ref. Karmakar  _et al._ (2010a), where it has been shown how the nature of
non-affine displacements in amorphous solids near a plastic instability with
only $\lambda_{min}\to 0$ gives rise to the square-root singularity. Vanishing
of only the lowest non-zero eigenvalue near the plastic instability ensures
that the spatial map of eigenmode corresponding to the vanishing eigenvalue
dominates the displacement field on the approach to such a plastic instability
Maloney and Lemaître (2006). The same scenario of a square-root singularity
with only one vanishing eigenvalue was shown to hold under athermal
quasistatic expansion on approach to a cavitation instability as well Dattani
_et al._ (2022).

In the current context of the athermal quasistatic expansion of pinned solids,
we find that the same square-root singularity scenario holds. In Fig. 4(a), we
show one such trajectory corresponding to $c=0.05$ & $N=2500$. Fig.4(b) shows
the square root singularity at which the lowest-eigenvalue of the Hessian
vanishes as $\lambda_{min}\sim\sqrt{(\rho-\rho_{c})/\rho_{c}}$ on approach to
the plastic instability at the points marked in Fig.4(a), where $\rho_{c}$ is
the point at which the plastic instability occurs. The Fig. 4(d)-(e) show an
eigenmode on the approach to the plastic instability, displacement field on
the approach to the plastic instability, and the displacement field across the
pressure drop, respectively, for one of the pressure-jumps in Fig.4(a)& (b).
As evident from the vector-field maps, the displacement fields on the approach
to instability are predicted by the eigenvector of $\lambda_{min}$, but the
displacement fields across the pressure jump do not have a high overlap with
the eigenmode on the approach to the instability. This occurs because of the
cascade/avalanche nature of the plastic jump. These avalanches are also more
localized spatially due to pinning, unlike those seen in unpinned solids which
can be system spanningLemaître and Caroli (2009); Karmakar  _et al._ (2010b);
Dattani  _et al._ (2022). All in all, the mechanism of plasticity on PEL being
analogous to unpinned solids does not come as a surprise because pinning only
blocks certain pathways of relaxation on the PEL of unpinned solids Niblett
_et al._ (2018).

Figure 5: (a) Spatial decay profile of magnitude of displacements in the
eigenmode before the yielding, $e(r)$, leading to the plastic events, shown
for $c=0.00,0.01,0.02,0.05,0.08$, with averaging done over $10$ plastic events
for each case. $r$ is the distance from the plastic center defined by a
particle undergoing maximum displacement in a mode. (b)Dependence of length
scale ($\xi_{\rm fit}$) extracted from fits to $e(r)$ (see text for
discussion) with pinning concentration; also shown is the behaviour of
$\xi_{\rm pin}$, extracted from pinning density, with $c$.

###  III.4 Spatial decay of plastic modes and a length scale of plasticity

The plasticity under shear deformation is known to occur via localised
rearrangements of particles in a shear transformation zone Argon and Kuo
(1979); Falk and Langer (1998); Richard  _et al._ (2021a). These localised
rearrangements(and their eigenmode) have a quadrupolar shape Eshelby (1957);
Gartner and Lerner (2016) and the radial part of the displacements of the
medium decay as $r^{-(d-1)}$; where $r$ is the distance from the center of
rearrangement, and $d$ is the spatial dimensions. Even though the
rearrangements are local, the displacement fields have a long-range character,
i.e., in 2d, they decay as $1/r$ from the center. The eigenmodes, $\vec{e}$,
corresponding to the vanishing eigenvalue $\lambda_{min}$ are known to decide
the direction of failure on the PEL for both shear and expansion Maloney and
Lemaître (2006); Tanguy  _et al._ (2010); Manning and Liu (2011); Dattani  _et
al._ (2022).

In the context of pinning, we, too, look at the radial decay profile of
displacements from the plastic center 111defined as particle undergoing
highest displacement in an eigenmode of the eigenmode $\vec{e}$, just before
the plastic event, i.e., at a distance of $\delta\rho\approx 10^{-7}$, in the
pre-yield regime. We choose to study the spatial profiles in the pre-yield
regime, where density inhomogeneities (which can interfere with the profile
shapes) are largely absent. The decay profiles averaged over $10$ plastic
events for each pinning concentration are shown in Fig.5(a). For the unpinned
case, as one would expect, the decay profile $e(r)\sim 1/r$. For the pinned
solids, interestingly, we find two regimes. For the smaller pinning
concentrations, viz. $c=0.01$, $c=0.02$ and $c=0.05$, we see that the decay
profile fits a screened power-law function $e(r)\sim\exp(-r/\xi)/r^{\eta}$
with $\eta\approx 0.42$ for $c=0.01$ and $c=0.02$ whereas $\eta=0.13$ for
$c=0.05$. For larger pinning concentration, viz. $c=0.08$, we find that the
decay profile is well-fitted by an exponential form $e(r)\sim\exp(-r/\xi)$
suggesting a crossover from a screening-like decay to a purely exponential
decay. The scale-dependent exponential part of the decay profile allows us to
extract a length scale $\xi$ from the fit parameters. The values of
$\xi_{fit}$ for different pinning concentrations are shown using red color in
Fig.5(b).

For a fixed pinning concentration $c$, one naturally obtains a length scale
$\xi_{pin}=\sqrt{1/(c\cdot\rho)}$, which denotes an average distance between
any two pinned sites. But, since the eigenmode decay profiles shown in
Fig.5(a) are sampled at different densities from the expansion trajectory, the
blue curve showing $\xi_{pin}=\sqrt{1/(c\cdot\rho)}$ for the sampled densities
in Fig.5(b). The length scale extracted from the fits in Fig. 5(a) thus shows
reasonable parallels with the length scale introduced by pinning as seen in
Fig.5(b). The decay exponents of length scales with pinning concentrations for
$\xi_{pin}$ is slightly different from exponent $0.5$ primarily due to the
changes in density where these decay profiles are obtained for each pinning
concentration.

The exponential nature of the decay profiles and the extracted length scales
thus suggest that, with increasing pinning concentrations, the displacements
in the eigenmode become more and more localized. This implies that the
consequence of plastic events in a pinned solid has reduced spatial effects,
unlike the unpinned solid, which has a non-local power-law character. Hence,
this will restrict the size of cascades/avalanches, which are caused due to
triggering of multiple such plastic modes at different spatial
locationsStanifer and Manning (2022). This is also consistent with the data in
Fig.1(e) & (f), where the magnitude of the drops in pressure/energy decrease
with an increase in pinning concentration. The presence of a plasticity length
scale in pinned solids also explains why there are negligible system size
effects in the data shown in Fig. 3 as long as $\xi<<L$, which is the case for
the system sizes that we report. As the length-scale $\xi$ prevents the system
from acting as a whole beyond the lengths $l>\xi$.

##  IV Summary & Discussion

To summarize, using extensive numerical simulations, we have investigated how
random inclusions, at the level of atoms/particles, influence the cavitation
failure in amorphous solids. Our studies are done in the athermal quasistatic
limit, i.e., where fluctuations originating from thermal fluctuations or
finite deformation rates are absent. Like in previous studies, we model the
inclusions as randomly pinned particles that only undergo affine motion during
the mechanical deformation process but do not undergo any non-affine motion
during the plastic relaxation.

We demonstrate how an amorphous solid with a small concentration of micro-
inclusions, intended to make the solid stronger, can delay the cavitation and
increase the tensile strength. The resulting glass obtained from the
inclusions is thus not only stronger glass but also less susceptible to
catastrophic fracture failure via cavitation. From a thermodynamic point of
view, the delay in cavitation and, thus, the increase in tensile strength is
likely a result of a change in the densities of coexisting solid and gas
phases in the temperature-density Sastry (2000); Testard  _et al._ (2011)
phase diagram of the amorphous solid. A sharp yielding-like cavitation
transition observed in amorphous solids Altabet  _et al._ (2018a); Dattani
_et al._ (2022) accompanied by an extensive peak in the pressure fluctuations
characterized via a susceptibility $\chi_{p}$ is suppressed by particle-
pinning. Particle pinning not only increases the tensile strength but also
decreases the size of avalanches that lead to cavitation in amorphous solids
by restricting them spatially; as a result, the cavitation occurs more
homogeneously in this pinned solid. The expansion response of a pinned solid
is somewhat analogous to the shear response of a ductile solid, as the sharp
macroscopic drop leading to yielding is absent. But the analogy with a ductile
solid doesn’t hold completely as the yielding via cavitation is delayed in
pinned solids. In a ductile solid, the yielding under shear would occur at
lower strains and bulk stress Yeh  _et al._ (2020). Also, with increasing
system sizes, the ductile solid shows brittle behaviour under shear via shear-
banding-like events, as claimed in a recent work Richard  _et al._ (2021b). In
our range of probing, the finite size effects in the pressure $P$ vs $\rho$
curves and the pressure fluctuation $\chi_{p}$ vs $\rho$ curves do not show
such a brittle-like character for the large system sizes. The absence of
finite-size effects in the $P-\rho$ and $\chi_{p}-\rho$ curves, along with
decreasing sizes of avalanches, is explained using a length-scale extracted
from the spatial decay of displacements in plastic-eigenmodes.

We have modeled these micro-alloyed inclusions via random pinning, which is a
drastic simplifying approximation. Future studies can perhaps create more
realistic models of inclusions in amorphous solids and study their deformation
response. In Ref. Das  _et al._ (2017), it was found that a particle that has
a larger size compared to the typical size of the particles of the host medium
can act as temporary pinning sites (termed as “soft pinning”) over a timescale
relevant to the medium relaxation time. Similarly, particles having different
geometric shapes, like rods, can have significantly lower mobility due to
their enhanced caging effect in a crowded environment and can also act as
possible micro-alloying agents. Recent experimental realization of this “soft
pinning” concept in colloidal glasses Ganapathi  _et al._ (2018), molecular
glasses Das  _et al._ (2021), and glassy polymer mixture Kikumoto  _et al._
(2020) is indeed very encouraging, indicating a clear possibility to test some
of our results in experiments.

##  V Acknowledgements

We thank the HPC facility at IMSc-Chennai for providing computational
resources. S.K. would like to acknowledge the support from Swarna Jayanti
Fellowship Grants No. DST/SJF/PSA-01/2018-19 and No. SB/SFJ/2019-20/05 and
Core Research Grant from SERB via grant CRG/2019/005373. PC acknowledges
funding from SERB via grant MTR/2022/001034.

## References

  * Schuh  _et al._ (2007) C. A. Schuh, T. C. Hufnagel, and U. Ramamurty, Acta Materialia 55, 4067 (2007).
  * Bonn  _et al._ (2017) D. Bonn, M. M. Denn, L. Berthier, T. Divoux, and S. Manneville, Reviews of Modern Physics 89, 035005 (2017).
  * Rodney  _et al._ (2011) D. Rodney, A. Tanguy, and D. Vandembroucq, Modelling and Simulation in Materials Science and Engineering 19, 083001 (2011).
  * Liu  _et al._ (2021) C. Liu, S. Dutta, P. Chaudhuri, and K. Martens, Physical Review Letters 126, 138005 (2021).
  * Bouchaud  _et al._ (2008) E. Bouchaud, D. Boivin, J.-L. Pouchou, D. Bonamy, B. Poon, and G. Ravichandran, Europhysics Letters 83, 66006 (2008).
  * Murali  _et al._ (2011) P. Murali, T. Guo, Y. Zhang, R. Narasimhan, Y. Li, and H. Gao, Physical Review Letters 107, 215501 (2011).
  * Shen  _et al._ (2021) L.-Q. Shen, J.-H. Yu, X.-C. Tang, B.-A. Sun, Y.-H. Liu, H.-Y. Bai, and W.-H. Wang, Science Advances 7, eabf7293 (2021).
  * Guan  _et al._ (2013) P. Guan, S. Lu, M. J. B. Spector, P. K. Valavala, and M. L. Falk, Phys. Rev. Lett. 110, 185502 (2013).
  * Chaudhuri and Horbach (2016) P. Chaudhuri and J. Horbach, Physical Review B 94, 094203 (2016).
  * Altabet  _et al._ (2018a) Y. E. Altabet, A. L. Fenley, F. H. Stillinger, and P. G. Debenedetti, The Journal of Chemical Physics 148, 114501 (2018a), https://doi.org/10.1063/1.5019274 .
  * Altabet  _et al._ (2018b) Y. E. Altabet, A. L. Fenley, F. H. Stillinger, and P. G. Debenedetti, The Journal of Chemical Physics 148, 114501 (2018b).
  * Paul  _et al._ (2020) K. Paul, R. Dasgupta, J. Horbach, and S. Karmakar, Phys. Rev. Research 2, 042012 (2020).
  * Dattani  _et al._ (2022) U. A. Dattani, S. Karmakar, and P. Chaudhuri, Phys. Rev. E 106, 055004 (2022).
  * Dattani  _et al._ (2023) U. A. Dattani, R. Sharma, S. Karmakar, and P. Chaudhuri, arXiv preprint arXiv:2303.04529 (2023).
  * Hofmann  _et al._ (2008) D. C. Hofmann, J.-Y. Suh, A. Wiest, G. Duan, M.-L. Lind, M. D. Demetriou, and W. L. Johnson, Nature 451, 1085 (2008).
  * Hardin and Homer (2015) T. J. Hardin and E. R. Homer, Acta Materialia 83, 203 (2015).
  * Torquato and Haslach Jr (2002) S. Torquato and H. Haslach Jr, Appl. Mech. Rev. 55, B62 (2002).
  * Tyukodi  _et al._ (2016) B. Tyukodi, C. A. Lemarchand, J. S. Hansen, and D. Vandembroucq, Physical Review E 93, 023004 (2016).
  * Garrett  _et al._ (2012) G. R. Garrett, M. D. Demetriou, J. Chen, and W. L. Johnson, Applied Physics Letters 101, 241913 (2012).
  * González (2016) S. González, Journal of Materials Research 31, 76 (2016).
  * Gendelman  _et al._ (2014) O. Gendelman, A. Joy, P. Mishra, I. Procaccia, and K. Samwer, Acta materialia 63, 209 (2014).
  * Bhowmik  _et al._ (2019a) B. P. Bhowmik, P. Chaudhuri, and S. Karmakar, Phys. Rev. Lett. 123, 185501 (2019a).
  * Bhowmik  _et al._ (2019b) B. P. Bhowmik, S. Karmakar, I. Procaccia, and C. Rainone, Phys. Rev. E 100, 052110 (2019b).
  * Liu  _et al._ (2023) Y. Liu, H. Liu, and H. Peng, Journal of Non-Crystalline Solids 601, 122052 (2023).
  * Brüning  _et al._ (2008) R. Brüning, D. A. St-Onge, S. Patterson, and W. Kob, Journal of Physics: Condensed Matter 21, 035117 (2008).
  * Polak and Ribiere (1969) E. Polak and G. Ribiere, Revue française d’informatique et de recherche opérationnelle. Série rouge 3, 35 (1969).
  * Plimpton (1995) S. Plimpton, Journal of Computational Physics 117, 1 (1995).
  * Anderson  _et al._ (1999) E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen,  _LAPACK Users’ Guide_ , 3rd ed. (Society for Industrial and Applied Mathematics, Philadelphia, PA, 1999).
  * Malandro and Lacks (1999) D. L. Malandro and D. J. Lacks, The Journal of Chemical Physics 110, 4593 (1999), https://doi.org/10.1063/1.478340 .
  * Maloney and Lemaître (2006) C. E. Maloney and A. Lemaître, Phys. Rev. E 74, 016118 (2006).
  * Karmakar  _et al._ (2010a) S. Karmakar, E. Lerner, and I. Procaccia, Phys. Rev. E 82, 026105 (2010a).
  * Lemaître and Caroli (2009) A. Lemaître and C. Caroli, Phys. Rev. Lett. 103, 065501 (2009).
  * Karmakar  _et al._ (2010b) S. Karmakar, E. Lerner, I. Procaccia, and J. Zylberg, Phys. Rev. E 82, 031301 (2010b).
  * Niblett  _et al._ (2018) S. P. Niblett, V. K. de Souza, R. L. Jack, and D. J. Wales, The Journal of Chemical Physics 149, 114503 (2018), https://doi.org/10.1063/1.5042140 .
  * Argon and Kuo (1979) A. Argon and H. Kuo, Materials Science and Engineering 39, 101 (1979).
  * Falk and Langer (1998) M. L. Falk and J. S. Langer, Phys. Rev. E 57, 7192 (1998).
  * Richard  _et al._ (2021a) D. Richard, G. Kapteijns, J. A. Giannini, M. L. Manning, and E. Lerner, Phys. Rev. Lett. 126, 015501 (2021a).
  * Eshelby (1957) J. D. Eshelby, Proceedings of the royal society of London. Series A. Mathematical and physical sciences 241, 376 (1957).
  * Gartner and Lerner (2016) L. Gartner and E. Lerner, Phys. Rev. E 93, 011001 (2016).
  * Tanguy  _et al._ (2010) A. Tanguy, B. Mantisi, and M. Tsamados, Europhysics Letters 90, 16004 (2010).
  * Manning and Liu (2011) M. L. Manning and A. J. Liu, Phys. Rev. Lett. 107, 108302 (2011).
  * Note (1) Defined as particle undergoing highest displacement in an eigenmode.
  * Stanifer and Manning (2022) E. Stanifer and M. L. Manning, Soft Matter 18, 2394 (2022).
  * Sastry (2000) S. Sastry, Physical Review Letters 85, 590 (2000).
  * Testard  _et al._ (2011) V. Testard, L. Berthier, and W. Kob, Physical review letters 106, 125702 (2011).
  * Yeh  _et al._ (2020) W.-T. Yeh, M. Ozawa, K. Miyazaki, T. Kawasaki, and L. Berthier, Phys. Rev. Lett. 124, 225502 (2020).
  * Richard  _et al._ (2021b) D. Richard, C. Rainone, and E. Lerner, The Journal of Chemical Physics 155, 056101 (2021b).
  * Das  _et al._ (2017) R. Das, S. Chakrabarty, and S. Karmakar, Soft matter 13, 6929 (2017).
  * Ganapathi  _et al._ (2018) D. Ganapathi, K. H. Nagamanasa, A. Sood, and R. Ganapathy, Nature communications 9, 397 (2018).
  * Das  _et al._ (2021) R. Das, B. P. Bhowmik, A. B. Puthirath, T. N. Narayanan, and S. Karmakar, arXiv preprint arXiv:2106.06325 (2021).
  * Kikumoto  _et al._ (2020) G. Kikumoto, N. Torii, K. Fukao, C. P. Royall, H. Yao, Y. Saruyama, and S. Tatsumi, arXiv preprint arXiv:2003.06089 (2020).