Enhanced vibrational stability in glass droplets

Abstract We show through simulations of amorphous solids prepared in open-boundary conditions that they possess significantly fewer low-frequency vibrational modes compared to their periodic boundary counterparts. Specifically, using measurements of the vibrational density of states, we find that the D(ω)∼ω4 law changes to D(ω)∼ωδ with δ≈5 in two dimensions and δ≈4.5 in three dimensions. Crucially, this enhanced stability is achieved when utilizing slow annealing protocols to generate solid configurations. We perform an anharmonic analysis of the minima corresponding to the lowest frequency modes in such open-boundary systems and discuss their correlation with the density of states. A study of various system sizes further reveals that small systems display a higher degree of localization in vibrations. Lastly, we confine open-boundary solids in order to introduce macroscopic stresses in the system, which are absent in the unconfined system and find that the D(ω)∼ω4 behavior is recovered.


Introduction
In contrast to crystalline solids, whose vibrational states are well described by the Debye model, amorphous solids exhibit anomalous mechanical and thermal properties (1)(2)(3)(4)(5)(6).For a system in d dimensions, the Debye model predicts a vibrational density of states (VDoS), D(ω) ∼ ω d−1 arising from phonons (7).On the other hand, disordered and anharmonic systems possess an excess of lowfrequency modes above the Debye prediction, termed the "Boson peak" (8)(9)(10)(11).It has been suggested that the large specific heat of glasses and the plastic failure of amorphous solids are intimately related to these low-frequency nonphononic vibrational modes (12)(13)(14)(15)(16).However, a complete understanding of the structural and statistical properties of these vibrational modes has been elusive and remains an important current topic of interest in the field of disordered solids.
Several theoretical frameworks for the vibrations in disordered systems seek to model amorphous solids as an ensemble of anharmonic oscillators.Such a treatment is motivated by a disordered arrangement of particles containing local "soft spots" where the stiffness associated with a collective vibration is very small.The phenomenological "soft potential model" (SPM) treats these regions as noninteracting oscillators and predicts a VDoS D(ω) ∼ ω 4 for the lowest frequencies of stable inherent structures (17)(18)(19).An extension to the SPM that includes effects of interactions between the anharmonic oscillators continues to retain the ω 4 behavior of the VDoS (20)(21)(22).Other recent theoretical studies also predict a D(ω) ∼ ω 4 at the lowest frequencies (23)(24)(25)."Fluctuating Elasticity Theory" on the other hand predicts a lowfrequency regime composed of extended modes that scale as D(ω) ∼ ω d+1 in d dimensions (26,27).More recent work suggests nonaffine displacements as the source of such non-Debye behavior (28,29).Certain mean field theories such as the "Perceptron Model" (30) and the "Effective Medium Theory" (31) predict nonphononic vibrations with a D(ω) ∼ ω 2 dependence.In this context, numerical investigations of amorphous solids in order to ascertain the nature of the low-lying excitations of such systems are of crucial importance.
Recent numerical studies of amorphous solids have identified a universal D(ω) ∼ ω 4 scaling in the low-frequency regime of the VDoS across a broad class of simulated model systems (32).The universality of this nonphononic power law at low frequencies has been established in two dimensions (2D), three dimensions (3D), and four dimensions (33)(34)(35).In order to extract this behavior, studies have focused on small system sizes, where the lowest frequency quasilocalized vibrations are well separated in energy from the first system-spanning phonon.This suggests that these quasilocalized modes (QLMs) are primarily responsible for the power-law tail in the VDoS (33,36).Further, suppression of system-spanning vibrations using random pinning protocols has been shown to enhance the nonphononic spectra, thereby displaying a pronounced D(ω) ∼ ω 4 behavior (37).Yet other studies have utilized measures such as participation ratios in order to isolate such QLMs that have been shown to contribute to the observed power-law behavior (38).Simulations of various systems including ultrastable glasses (39), silica models (40), long-ranged models (41), finite-temperature systems (42), and random matrix models (43), all provide significant evidence of the ubiquity of the ω 4 regime of the VDoS.However, other studies simulating amorphous solids report deviations from the universal quartic law.High parent temperatures, poor annealing, and small system sizes have each been shown to result in power-law exponents that are less than 4 (44)(45)(46).Recent studies have further identified exponents of 3 and 3.5 in the low-frequency regime (47)(48)(49).Confined three-dimensional thin films have also been shown to possess a low-frequency VDoS of ω 3 (50).
A crucial aspect often overlooked in simulations of amorphous solids is the residual shear stresses arising due to periodic boundary conditions (PBC) (51).This implies that energy-minimized configurations of disordered systems under PBC are unstable to shear deformations.It has recently been shown that the low-frequency regime of the VDoS is modified to D(ω) ∼ ω 5 when considering ensembles stable to simple-shear perturbation (52).Significantly, such an increase in the exponent points to a correlation between the shear stability of the system and a reduction in the propensity for localized vibrations.Given these observations, a natural question that arises pertains to the consequences of stabilization against all possible deformations.Solids formed under open-boundary conditions (OBC) are a suitable candidate since all elements of their pressure tensor are identically zero for each configuration.Such a state is permitted by a lack of confinement at the boundaries.It is important to note that the OBC system is no longer isotropic and introduces possible radial heterogeneity in structure as well as relaxation dynamics.Furthermore, unlike systems under PBC, phonons in OBC are not required to obey the artificial symmetries enforced by the periodicity.Consequently, solids under OBC can accurately capture features of natural solids, including surface as well as system-size effects.
In this paper, we report a characterization of the VDoS of openboundary amorphous solids.The low-frequency vibrational properties of open systems remain relatively unexplored (53), and our study forms the first such computational examination of the localized modes of open-boundary amorphous solids.We observe the low-frequency vibrational spectrum of a simulated model amorphous solid under OBC to be of the form D(ω) ∼ ω δ with δ > 4, both in 2D and 3D.Since an increase in the exponent implies a reduction in the degree of vibrational localization, the corresponding solid ensemble may be said to possess enhanced stability.Notably, the model system displays such a VDoS only when configurations are annealed to their inherent structures.On the other hand, commonly used quenching protocols lead to ensembles with δ ∼ 4. Interestingly, at large enough system sizes, the exponent saturates to a value δ = 5 in 2D and δ = 4.5 in 3D.A detailed investigation of the average stress profile of the solids allows us to identify a surface layer that suggests a source of the system-size effects.Lastly, we also observe that confining open-boundary solids under a harmonic trap recovers the D(ω) ∼ ω 4 behavior of PBC systems, confirming the role of stresses in the stability of solids.
The outline of the paper is as follows.The first section describes the numerical protocol we use to generate stable open-boundary solids.In the second section, we demonstrate the effect of different minimization protocols to generate stable solids.In the third section, we present data to demonstrate the effect of the system size on the vibrational spectrum.In the fourth section, we provide an analysis of the source of the stabilization.In the fifth section, we examine the effect of confinement on the vibrational spectrum of solids.Finally, we conclude and provide directions for future research.

Model and simulation details
In order to simulate the behavior of amorphous solids in OBC, we use models with attractive interactions between the particles.Specifically, we use variants of the canonical Kob-Andersen Lennard-Jones model (54,55) in 2D and 3D.The model consists of binary mixtures of particles in number ratios 65:35 in 2D and 80:20 in 3D, with an interaction potential: where α, β ∈ {A, B} correspond to the two types of particles resulting in three types of interactions.The potential is cut off at a distance r αβ c = 2.5σ αβ .σ AA = 1 is the unit of length, and ϵ AA = 1 is the unit of energy (Boltzmann's constant being unity).The remaining parameters are ϵ AB = 1.5, ϵ BB = 0.5, σ AB = 0.8, and σ AB = 0.88 with all masses set to unity m = 1.0.We equilibrate liquid configurations at high temperatures (T = 0.55 in both 2D and 3D) by performing constant temperature molecular dynamics simulations under PBC.The liquid is then cooled to a low temperature (T = 0.2) at a rate Ṫ = 10 −4 .In order to achieve a density close to that of a system under OBC, we perform zero-pressure molecular dynamics simulations using the isobaric ensemble.
Finally, in order to study solid droplets under OBC, we begin by cutting out a circular (spherical in 3D) region of the zero-pressure liquid.We ensure that all the droplet samples contain a fixed number of particles, say N, by selecting the N closest particles to the center of mass of the corresponding liquid sample referred to henceforth as a "cutout" sample.Solids are configurations that resist deformations, which correspond to the inherent structures of the amorphous configuration of particles.We generate such structures by performing an energy minimization of the particle-position degrees of freedom using three different protocols, namely: (i) Molecular dynamics simulations in the presence of dissipative viscous drag in the athermal limit (damped dynamics, DD) (56), (ii) fast inertial relaxation engine (FIRE) (57), and (iii) the nonlinear conjugate gradient (CG) algorithm (58).Note that we do not equilibrate the cutout sample but instead directly quench or anneal through the energy minimization protocol.
Simulations of the largest system-size three-dimensional zeropressure PBC liquids are performed using LAMMPS (59).The vibrational properties of the solid configurations are probed through the Hessian matrix as defined in Eq. 4. We perform eigenvalue computations using the Intel Math Kernel Library (60) sparse solver routine mkl_sparse_d_ev.Below, we discuss the three protocols we employ in order to find the minimum-energy amorphous configurations.

Damped dynamics
In the DD protocol, a dissipative viscous drag in the absence of a temperature bath or other energy input serves as a method to find the stationary state of the conservative force fields.Such an evolution of the system is described by the equation of motion where  F i is the conservative force on particle i due to interactions with other particles,  v i is the instantaneous velocity of the particle i and γ is the viscous damping parameter.A stable steady-state solution of these equations corresponds to a minimum-energy configuration of the system of particles with  F i = 0 and  v i = 0 for all particles.We employ a velocity Verlet integration scheme modified to incorporate velocity-dependent accelerations.These damped equations of motion allow the system to traverse the basins of multiple inherent states before selecting and settling in a more stable minimum.We employ small damping constants so that the system is able to execute such a relaxation.At large values of γ, the dynamics would be equivalent to a steepest-descent minimization of the energy.Unless explicitly mentioned the value of the damping constant γ is 0.1.

FIRE minimization
The FIRE is an efficient protocol designed to find local minima of multidimensional functions (57).The procedure involves numerically integrating a dynamical equation with variable viscous damping and, additionally, a gradient director.The equation of motion for each particle is where F(t) refers to unit vector along  F(t).These two aspects of the FIRE algorithm that lend it its speed also prevent the system from leaving its original energy basin.Therefore, the inherent structure obtained is primarily a function of the equilibrium sampling and is not dependent on the stability of the minimum of the basin.

CG minimization
We also test the stability of configurations obtained via a CG minimization scheme (58).We implement the nonlinear CG for the open boundary system using the Polak-Ribiere method to obtain the updated direction and the secant minimization method in order to determine the step size of the line search.As in the case of the FIRE algorithm, the CG leads the system to its nearest minimum.

Stable vibrational modes
Harmonic vibrations of a solid can be probed by diagonalizing the "Hessian matrix," defined as where U(r 1 , r 2 , . . .r n ) is the total potential energy of the system and r i denotes the position of particle i.The eigenvectors (ψ k ) of the Hessian matrix correspond to the cooperative displacements of the particles participating in harmonic vibrations (termed normal modes), and the corresponding eigenvalues (λ k ) represent the frequency of that vibration.The system, when perturbed along a normal mode, will perform a pure oscillation with a frequency The Hessian matrix evaluated at an energy minimum of the landscape is positive semi-definite with zero modes corresponding to the global invariants of the Hamiltonian (Goldstone modes).For example, a 3D system under PBC possesses three translational degrees of freedom, and the corresponding Hessian, therefore, has three zero modes.Systems under OBC, on the other hand additionally also possess rotational degrees of freedom.In general, in d dimensions, such a solid has d(d+1) 2 zero modes corresponding to d translational modes and d(d−1) 2 rotational modes.The VDoS of a system of N particles in d dimensions is defined as In Fig. 1a and b, we display a typical open-boundary 2D solid configuration of N = 400 particles and its first nonzero mode corresponding to a frequency ω = 3.24107 × 10 −1 with a participation ratio (defined in Eq. 16) of 1.313714 × 10 −1 .Similarly, in Fig. 1c  and d, we show a typical 3D solid of N = 4, 096 particles and it's first mode corresponding to a frequency ω = 4.039514 × 10 −1 and participation ratio of 6.792794 × 10 −3 .These low-frequency modes are quasilocalized with a small number of particles displaying large participation by forming the core of the vibration.Interestingly, these eigenmodes resemble the quadrupolar modes observed in the low-frequency regime of the vibrational density of states of structural glass formers under PBC (32).

Anharmonic stability analysis
While the Hessian characterizes the curvature of the energy minimum and, thereby, the frequency of vibrations, an anharmonic analysis is necessary in order to understand the vicinity of the energy minimum better.The energy near the minimum may be approximated up to the fourth order as where s is the scalar distance of particles from the minimum along a given direction in the energy landscape.If the particles are displaced along an eigenvector (ψ) of the Hessian (H), is the eigenvalue of the Hessian, and quantify the nonlinearity of the energy landscape along that mode.Given a vibration in the energy landscape as defined by an eigenvector, the corresponding energy minimum is said to be "stable" provided that the minimum is the deepest in its neighborhood.Under the quartic approximation of the energy near the minimum, such stability is achieved when (18)

Effect of minimization protocol
Simulated models studying the properties of amorphous solids have typically been systems under periodic boundaries.It has

Chakraborty et al. | 3
been shown that in such systems, the CG minimization protocol generates configurations with stable minima, provided that configurations are obtained from a sufficiently low parent temperature (32).However, the open-boundary configurations obtained by quenching the cutouts of zero-pressure PBC systems contain, among them, some inherent structures at unstable minima.We, therefore, explore an energy-minimization protocol that anneals the system into an ensemble comprised of primarily stable minima.
We find that different energy-minimization protocols drive the system to different minima.In Fig. 2a, we show the displacements incurred by each particle from a particular liquid configuration when undergoing energy minimization through two different protocols, namely DD and CG.While both procedures yield inherent structure configurations, they do not correspond to the same minimum.Most significantly, DD finds minima in a relatively expanded region in the energy landscape as compared to CG and FIRE.In Fig. 3a, we present corroborating numerical evidence with particles undergoing relatively larger displacements in the case of DD.In order to qualify this difference, we compute the distance between the cutout liquid configuration and the energyminimized solid configuration as where  r L i and  r S i represent the position of particle i in the liquid and solid configurations, respectively.In terms of the radial and tangential displacement incurred by particles at different radial distances, where u r (r) and u θ (r) are the radial and tangential displacement of particles at a radial distance r from the center of mass of the liquid droplet.In Fig. 3, we show the average (b) radial and (c) tangential displacement at different radial distances of the liquid configurations.For the case of DD minimization, particles near the surface undergo comparatively larger displacement in the tangential direction, which results in a stable inherent structure.Since the viscosity of the surrounding medium determines the exact extent of this search space, we use an appropriately chosen damping, optimizing for both improved search area as well as speed.
We now consider the vibrational properties of the minima obtained through the various energy minimization protocols.In Fig. 2, we show the lowest frequency vibrational mode of a solid configuration obtained via (b) CG and (c) DD minimization of the same liquid configuration.Although both modes show signatures of similar quadrupolar vibrations, the inherent structure corresponding to (a) does not satisfy the stability criterion defined in Eq. 10.
An important difference between the ensembles of minima obtained via the various protocols lies in their vibrational stability as defined in Eq. 10.Through Fig. 4 frequency mode in order to extract the stability of the minimum.Interestingly, the CG and FIRE minimization protocols present many unstable inherent structures indicated by the data points above the stability line corresponding to (B 2 3 = 3B 2 B 4 ).At the same time, the DD minimizer produces predominantly stable inherent structures as shown by the data points being well within the regime of stability.This suggests that the DD minimizer allows the system to find a locally deeper minimum.
Since the choice of minimization protocol affects the stability of minima, we next examine its effects on the vibrational spectrum.In Fig. 5, we plot the distribution of the low-lying vibrational frequencies of two-dimensional solid ensembles.While the VDoS of CG and FIRE minimized ensembles display the erstwhile universal lowfrequency power law D(ω) ∼ ω 4 (32,34), ensembles generated via the DD minimizer show a power law of ω 5 in the low-frequency regime, as has been seen previously in shear-stabilized systems (52).Generically the smaller the lowest frequency mode, the closer the system is to instability.An increase in the power of the VDoS from 4 to 5 corresponds to a reduction in the proportion of such smallvalued modes.Therefore, such a change is a signature of a transformation to a more stable solid ensemble enabled by an appropriate choice of an annealing energy minimizer.
In the inset of Fig. 5, we plot the distribution of where states within the dashed lines (|β| < 1) correspond to stable minima.We find that the stability condition is best satisfied by configurations obtained via DD minimization.Thus, we observe that the model system, under OBC, displays a correlation between stable minima and the vibrational stabilization of the corresponding solid ensembles.

System size effects
It is well known that the low-frequency vibrational spectrum of amorphous solids contains system-spanning phonons in addition to localized vibrations (32).The long-wavelength phonons in a solid with linear dimension L, under PBC, possess a frequency that varies with system size as L −1 .This decrease in the frequency of the phonons with increasing system size leads to difficulty in the characterization of quasilocalization at low frequencies.At the same time, the regime of the VDoS that is most amenable to the study of quasilocalized behavior are the frequencies below the first phonon (61).As the phonon frequencies become comparable to the frequencies of QLMs at larger system sizes, it is then necessary to perform a disorder averaging over an increasingly large number of samples in order to distinguish these modes.This becomes computationally infeasible at very large system sizes.Additionally, QLMs also display finite-size effects when the extent of these modes becomes comparable to the size of the system (44).The softness of QLMs, therefore, becomes independent of system size once it surpasses the length scale of the localized vibration.In this context, very small system sizes have been shown to possess an excess of quasilocalized vibrations with D(ω) ∼ ω δ where δ < 4 (44-46), as the finite-size effects dominate in such situations.Solids prepared under OBC can display additional surface effects in their low-frequency vibrational spectrum, as boundary effects not present in PBC systems also play a role.The stiffness of particle motion near the surface of the solid is small compared to the bulk of the material.Therefore, the eigenmodes are softer near the boundary of the solid.It is thus natural to expect that in small open-boundary solids, there will be stronger surface effects than in larger solids.This occurs in addition to the finite-size effects arising from the interplay between the size of the system and the size of the QLMs.It is therefore important to study the variation of the VDoS with system size under OBC in order to identify regimes of validity of the enhanced stability.
A quantitative characterization of the surface effects may be performed through a measurement of the stress distribution within  such solids.Moreover, since the low-lying VDoS is crucially sensitive to the macroscopic stress in the amorphous structure (52), it is reasonable to assume that the finite-size effects displayed by the stress profiles also carry over to the VDoS.In PBC systems, the average stress displays homogeneity over the sample, and therefore the local elastic behavior is also independent of the location in the solid.However, as the macroscopic stress is precisely zero in openboundary solids, any bulk stress is necessarily counterbalanced by the stresses on a "boundary layer."Such a distinct surface region possessing a typical thickness is yet another source of finite-size effects.Additionally, the local elasticity properties of the material are correlated with the local stress distribution, with particles near the boundary displaying lower stiffness in their interaction.In this context, we study the boundary stress layer in solids prepared under OBC, in order to systematically analyze the finite-size effects present in their low-frequency VDoS.

Stress distribution and the boundary layer
In Fig. 6, we describe the stress distributions which highlight the existence of a surface layer in an open-boundary solid.Specifically, in Fig. 6a and b, we plot the spatial distributions of the pressure and shear stress, respectively, as seen in openboundary solids comprised of 900 particles.These stress profiles are constructed by averaging over 10 5 amorphous configurations and using coarse-graining boxes each of dimension 0.28 σ AA × 0.28 σ AA .The typical radius of a configuration is approximately 15 σ AA resulting in about 9, 000 bins. Figure 6c plots the radial distribution of the pressure for amorphous solids of various sizes.The stress distributions display clear indications of a boundary layer, with the effects of the surface spanning over approximately 5 σ AA .We also find that while the pressure in bulk is sensitive to the total size of the droplet, the thickness of the boundary layer is largely independent of the system size.
The microscopic stress field (σ ij (r)) can be expressed as where σ 0 ij (r) represents the stress tensor of the initial liquid configuration and δσ ij (r) is the change in the stress tensor which results from the process of energy minimization.The original distribution of stress σ 0 ij in the liquid is dependent on the preparation protocol.The effect of removing a circular cutout from the zero-pressure liquid configuration gives rise to a nontrivial stress distribution near the boundaries of the cutout.This excess stress causes the particles near the edges to be displaced further inwards (due to the attractive nature of the interaction) near the boundaries.This leads to a new stress-balanced state, with δσ ij displaying larger changes near the boundary.The condition of mechanical equilibrium at the boundary of the solids is satisfied when the normal force acting on the surface due to the internal pressure exactly balances the surface tension force.Our numerical observations suggest that the stress profile within the surface layer is largely independent of the system size, leading to a constant surface tension across different system sizes.This suggests the bulk pressure will decrease as the radius of the solid is increased.In Fig. 6d, we plot the bulk pressure for different system-size solids.We observe a P ∼ 1 R scaling of the bulk pressure with the radius of the solid droplet R, which is consistent with a "Laplace law" for the pressure in a droplet (62).

Vibrational density of states
Having characterized the boundary stress layer, we next turn to the system size dependence of the VDoS.We expect the effects of the boundary to be significant only in system sizes for which the width of the boundary is comparable to the radius of the droplet.In order to discern the lengthscales corresponding to such a crossover, we analyze solids of various system sizes with particle numbers ranging from N = 144 to 1600 in 2D and from N = 216 to 10,000 in 3D.
In Fig. 7, we report the effect of system size on the lowfrequency vibrational spectrum of open-boundary solids in 2D.Boundary layers approximately five particles thick are observed in all systems, independent of size.d) Bulk pressure as a function of system radius.We observe a 1 R scaling of the bulk pressure where R is the radius of the solid droplet.
The small systems (N = 144, 256) display a higher degree of softness in their vibration, as can be seen from their low-frequency behavior of D(ω) ∼ ω 4.5 .The VDoS of systems with particle numbers equal to and larger than 400 appear to exhibit a lowfrequency regime of D(ω) ∼ ω 5 .This shows a crossover as the size of the system is increased beyond a length scale, with the lowfrequency VDoS showing a significant change in the power law, with the exponent varying from 4.5 to 5 in 2D.
The enhanced stability in the VDoS at large system sizes suggests that QLMs in the bulk result in the power law ω 5 , whereas modes localized on the surface are much softer and more unstable.We show further analysis comparing bulk and boundary localization in the next section.Finally, such behavior in the thermodynamic limit of the OBC systems under consideration is consistent with earlier studies in PBC systems (52).
In Fig. 8, we plot the VDoS for different system sizes in 3D.Here we observe a low-frequency behavior of D(ω) ∼ ω 4.5 in contrast to the results in 2D.In the next section, we perform further analyses to examine this behavior.Moreover, as may be seen in the inset of Fig. 8, we further confirm the robustness of this result by examining the VDoS of solids generated at a much higher degree of annealing by utilizing a damping constant that is one order of magnitude smaller (γ = 0.01).

Eigenmode localization-surface vs bulk
In this section, we provide further analysis of the lowest frequency vibrational modes in 2D in order to determine the potential source of the stable vibrations achieved under OBC.
The spatial extent of a mode may be evaluated by measuring the participation ratio (PR), where |ψ k i 〉 denotes the d-dimensional component of the eigenmode ψ k , corresponding to particle i. System-spanning modes possess a PR close to unity, while spatially localized modes display a PR ∼ 1/N.
The contribution of the particles on the surface of the solid to a particular mode may then be estimated through the "Surface Participation" (SP) as defined below: The surface is defined as encompassing all the particles belonging to the boundary layer as determined by the stress analysis described in the previous section.Modes localized on the surface of the solid display an SP close to 1, whereas a vibration that is localized deeper within the bulk of the solid incurs an SP closer to 0. In Fig. 9a, we plot the probability distribution of SP for the lowest frequency modes in two-dimensional solids of system size N = 400.The peak of the distribution corresponds to modes with extended, system-spanning vibrations that all particles participate in.By observation, we consider modes with an SP ≤ 0.7 to be predominantly bulk-localized vibrations and modes with an SP ≥ 0.8 to be predominantly surface-localized.In Fig. 9b and c, we plot histograms of the lowest frequency divided by ω δ min in order to extract the power law.We find that bulk-localized modes (Fig. 9b) display the stabilized δ ≈ 5 while the surface-localized modes (Fig. 9c) present a more unstable δ ≈ 4.5.This points to the possibility that the low-stress environment of the bulk contributes stable vibrations, while the large stresses at the surface result in modes of lower stability.Furthermore, such a separation potentially explains the large degree of softness observed in the vibrational spectrum at small system sizes.
We now provide a characterization of the lowest frequency vibrational modes in three-dimensional open-boundary solids.We study the contribution of surface particles on these modes through their "surface participation" (SP) as defined in Eq. 17.In Fig. 10a, we plot the probability distribution of SP corresponding to the lowest frequency modes in three-dimensional solids of various system sizes.The modes with a value of SP corresponding to Fig. 7. System size dependence of the low-frequency vibrational spectrum of open-boundary solids in 2D.We have sampled the first 100 nonzero eigenvalues of 2 × 10 5 configurations for system sizes N = 144 and 256.The histograms for the larger system sizes are constructed using 5 × 10 5 configurations.At small system sizes, the surface effects are more pronounced leading to more soft modes in their vibrational spectrum.They display a D(ω) ∼ ω δ with δ < 5 in the low-frequency regime.Such effects are reduced upon increasing the size with large systems displaying a D(ω) ∼ ω 5 .Fig. 8. System size dependence of the low-frequency vibrational spectrum of three-dimensional open boundary solids generated by annealing (γ = 0.1).We sample the first 100 nonzero eigenvalues of 2 × 10 5 configurations for system sizes N = 512, 2,197, and 4, 096.For the larger system of size N = 10,000, the histogram is generated by sampling 8 × 10 4 configurations.(Inset) Distribution of low-frequency modes of various system sizes of three-dimensional solids generated via much slower annealing (γ = 0.01).These distributions are drawn by sampling the first 100 low-frequency modes of at least 80,000 configurations for each system size.Solids corresponding to all the system sizes and degrees of annealing studied display a D(ω) ∼ ω 4.5 in the low-frequency regime of their VDoS.

Chakraborty et al. | 7
the peaks of the distributions are extended, system-spanning vibrations.In the case of the distribution corresponding to a system size N = 4, 096, we define the modes with SP ≤ 0.82 as being predominantly bulk-localized vibrations and modes with ≥ 0.9 as predominantly surface-localized.In Fig. 10b and c, we plot histograms of the lowest frequency for bulk and surface-localized modes, respectively.We find that bulk-localized modes display a δ ≈ 4.5 and surface-localized modes show a δ ≈ 4.

Effect of confining stresses
In this section, we show that the vibrational properties of systems under periodic boundaries can by reproduced by imposing macroscopic stresses on open-boundary solids.
The stresses in a system are described through the force moment tensor, defined as where f ij α is the α-component of the force on particle i by particle j, r ij β is the β-component of the vector distance between the particles i and j, σ αβ is the macroscopic stress tensor and A is the area of the two-dimensional system.Solids prepared under OBC have precisely zero macroscopic stresses, i.e. σ xx = σ yy = σ xy = 0.As described in the previous sections, such systems display a D(ω) ∼ ω δ with δ > 4, in contrast to PBC systems where δ = 4.It is, therefore, interesting to analyze the effects of re-introducing stresses on the vibrational spectrum of solids.
In order to accomplish this, we have studied the vibrational modes of solids confined in a radially symmetric harmonic trap centered at the center of mass (CM) of the system, as follows: where κ is the spring constant of the trap, r CM is the location of the CM and V(r 1 , r 2 , . . .r n ) represents the pair-wise particle interaction.We then perform DD to obtain energy-minimized solid configurations.Since the center of the harmonic trap is fixed, it breaks translational symmetry.However, it being radially symmetric, the system retains one zero mode corresponding to global rotations.The translational zero modes are substituted with eigenvalues corresponding to the stiffness of the confining harmonic potential.We, therefore, sample the low-frequency modes while excluding the zero modes as well as this trivial addition to the spectrum.In Fig. 11, we display the effect of confinement on the low-frequency vibrational spectrum at various strengths of the harmonic trap.Remarkably, we find that confinement results in a higher degree of softness in the vibrational spectrum of the solids.(Inset) Radial distribution of pressure for three-dimensional solids of various system sizes.In b) and c), we plot P(ω min ) divided by ω δ min for bulk and surface localized modes, respectively.Modes corresponding to SP ≤ 0.82 are considered bulk-localized vibrations, and modes with SP ≥ 0.9 are considered to be surface-localized.Both histograms are constructed using samples of the lowest frequency modes from 5 × 10 4 configurations.Bulk-localized modes display a δ ≈ 4.5, whereas surface-dominated modes show δ ≈ 4.

Conclusion and discussion
In this paper, we have performed a detailed characterization of the vibrational modes of amorphous solids prepared under OBC in both 2D and 3D.We showed that structures prepared under OBC differ crucially in their stability properties in comparison to solids prepared with PBC.Specifically, we observed that the ∼ ω δ with δ = 4 seen in systems under PBC is modified to a δ ≈ 5 in 2D and δ ≈ 4.5 in 3D, for solids under OBC.This points to the fact that openboundary solids which lack any macroscopic stresses are inherently more stable than their periodic boundary counterparts.These results reinforce the phenomenon observed by Krishnan et al. (52) where shear-stabilized (σ αβ = 0 ∀α ≠ β) configurations under PBC display an increase in the exponent to δ ≈ 5. Further, we probed the nature of energy minima of systems under OBC through an analysis of the anharmonic coefficients associated with their lowest frequency vibrational mode.Surprisingly we have found that the anharmonic stability of the minimum is sensitive to the protocols employed in generating solid configurations.Moreover, the nature of the minima under consideration is also correlated with a change in the exponent of the VDoS.Specifically, the ensembles with unstable minima correspond to a D(ω) ∼ ω 4 , whereas ensembles with stable minima show larger exponents.In particular, we observed clear enhancements in the stability of configurations resulting from dissipative dynamical processes that anneal the system when compared to those obtained via quenching minimization protocols.Next, we also characterized the dependence of the low-frequency modes on the system size.Solids at small systems sizes show a higher degree of softness in their low-frequency vibrations, as evidenced by a decrease in the exponent of the VDoS.We also performed an analysis of eigenmode localizations on the boundary and bulk of the system in order to decipher the source of the enhanced stabilization in systems under OBC.We found strong indications that the bulk of the solid predominantly contains stable vibrations, whereas the surface supports unstable modes.Finally, in order to isolate the cause of instability under PBC, we studied the vibrational spectrum of confined solids using a harmonic trap that introduces finite macroscopic stresses in the solid.We indeed showed that the D(ω) ∼ ω 4 behavior of the VDoS is recovered with an increase in the strength of the confinement.This is consistent with recent work where it has been shown that an increasing strain corresponds to an increase in the propensity for low-frequency modes (63).
Our study highlights aspects of the protocols used to create amorphous structures that have direct consequences on the nature of the sampled energy minima.We showed that various energyminimization procedures populate qualitatively different minima of the landscape, as has also been observed in previous studies (64,65).Our work in the context of open-boundary solids sheds light on the sparse nature of stable minima in such realistic systems.These results further raise important questions about the limits of stability that may be achieved through sufficient annealing.
Through our analysis of the stress profiles of solids under OBC, we found that the crossovers between the stable and unstable behaviors of the VDoS are well described by a bulk-boundary decomposition of the solids.The enhanced stability is derived from the low-stress environment of the bulk of the solid supported by constant boundary stress.These results point to the fact that the frozen-in stresses are crucial in determining the stability of the amorphous solids as well as their VDoS.It would be interesting to use the characterization of the boundary layer and stress distributions that appear in such confined systems to understand the correlations between the macroscopic stress tensor and the distribution of eigenmodes.
Notably, our ensembles exhibit liquid-like surface tensions, as seen by the conformance of bulk pressures to the Laplace law.This observation bares further examination in order to probe any correlation to the enhanced stability of open-boundary solids.Although our results for the cutout protocol reveal striking properties of vibrations in solids under OBC, it would be intriguing to further characterize this behavior with other preparation protocols, such as by evaporation and deposition, which is of direct experimental relevance.Finally, our results in three-dimensional systems regarding their smaller power-law exponent D(ω) ∼ ω 4.5  in comparison to two-dimensional systems is surprising in the context of earlier work displaying a universal ω 5 behavior in shear-stabilized systems under PBC.This points to the possibility that systems in higher dimensions support more localization of vibrations.It would be interesting to explore aspects of this nonuniversality in open-boundary systems.

Fig. 1 .Fig. 2 .
Fig. 1. a) A two-dimensional solid generated from a 400 particle, circular "cutout" of a liquid configuration, via damped dynamics (DD) minimization under OBC.b) The first nonzero vibrational mode of the two-dimensional solid with a participation ratio of approximately 0.26 and a frequency ω ≈ 4.34 × 10 −1 , displaying typical quasilocalized characteristics.c) A three-dimensional solid configuration of 4, 096 particles under OBC generated using DD energy minimization.d) The first nonzero vibrational mode of the three-dimensional solid with a participation ratio of 0.0068 and a frequency ω ≈ 4.04 × 10 −1 .Note that the solid borders of the figures are to aid depth perception, especially in 3D.The systems themselves have no boundaries.

Fig. 3 .
Fig. 3. a) Distribution of the distances between the open-boundary liquid and solid configurations for different minimization protocols.Particles in configurations that underwent DD minimization show relatively larger displacements.(Inset) The distribution plotted in log-scale to highlight the difference at large d c .Average b) radial and c) tangential displacements incurred by particles at different radial distances.DD produces a large displacement of the particles in the tangential direction, further enhanced near the surface of the liquid as compared to other minimization protocols.

Fig. 4 .
Fig. 4. Scatter plot of B 2 3 against B 2 B 4 with a sample size of 2 × 10 5 solid configurations.The dashed line corresponds to B 2 3 = 3B 2 B 4 with the region of stability lying below.CG and FIRE minimization protocols produce some inherent structures at unstable minima, whereas DD primarily generate stable minima.

Fig. 5 .Fig. 6 .
Fig. 5. Distribution of the low-frequency vibrational modes of a 400-particle system in 2D under OBC.The histogram is generated by sampling the first 100 nonzero eigenvalues of 3 × 10 5 configurations.The solid ensembles generated using CG and FIRE minimizers display a low-frequency behavior of ω 4 , whereas using the DD minimizer displays ω 5 .(Inset) Distribution of the stability factor β = B 3 / ������� 3B 2 B 4 √ for solids obtained via the various minimization protocols.The dashed vertical lines correspond to |β| = 1, with the enclosed regions containing stable configurations.

Fig. 9 .Fig. 10 .
Fig.9.a) Distribution of "Surface Participation" (SP) for the lowest frequency modes in two-dimensional solids of system size N = 400.Modes localized on the surface of the solid display an SP of 1, whereas when SP = 0, the vibration is limited to the inner bulk.The peak of the distribution corresponds to extended modes spanning the whole solid.The distribution of the lowest frequency modes with predominantly b) bulk and c) surface localizations.Here, modes corresponding to an SP ≤ 0.7 are considered to be bulk-localized vibrations and modes with an SP ≥ 0.8 are considered surface-localized vibrations.Both histograms are constructed using 10 5 lowest frequency modes.The bulk-localized modes display a δ = 5, whereas surface-localized modes show a δ = 4.5.