Inducing mechanical self-healing in polymer glasses

Self-healing refers to the ability of a material to restore functionality after undergoing physical damage. This characteristic is prominent in living materials such as plants¹, mammalian cells², and human skin³. In contrast, synthetic materials lack self-healing mechanisms, and the accumulated damage over time leads to macroscopic fractures and mechanical failure. In recent decades, the search for self-healing mechanisms in synthetic materials has become increasingly important for a sustainable environment^4,5, promising to lower waste production and energy costs by diminishing the need for manufacturing replacements. Chemical modifications introducing dynamic covalent^6,7 and physical bonds^8,9 that break and reform reversibly can add self-healing abilities to synthetic materials. However, these mechanisms strongly rely on a sufficiently high molecular mobility within the material which facilitates mass transport, favoring the reconnection of crack surfaces and subsequent healing^10,11.

Molecular mobility typically shows a dramatic decrease once a material is cooled below its glass-transition temperature T_g. Below T_g, materials exhibit rigidity because molecules are kinetically arrested. Above T_g, materials become flexible and viscous. Although soft materials such as elastomers and hydrogels maintain a solid shape due to covalent bonds, molecular segments have high local mobility at room temperature T_r, i.e. T_g < T_r, which enables the use of intrinsic mechanisms to induce self-healing in soft materials^{12,13,14,15,16}. In contrast, polymer glasses evidence an extremely low local mobility, since T_g >> T_r. The kinetically arrested nature of glassy systems gives rise to mechanical stability and a high elastic modulus, but disables potential intrinsic mechanisms that induce self-healing in soft materials. To overcome this low mobility, alternative approaches rely on the synthesis of polymers with complex architectures or functionalization^8,17,18 that are not suited for large-scale production since the synthesis becomes complex, expensive and may give rise to mechanical properties different from those required^19,20,21,22. Therefore, these strategies apply only to very specific materials and create a critical bottleneck for the development of general self-healing polymer glasses such as everyday plastics.

Here, we demonstrate by numerical simulations a general alternative approach that involves applying controlled oscillatory deformations on damaged polymer glasses to heal micro-cracks. Traditionally, deformations are associated with damage and fracture propagation. Large deformations in glassy materials trigger mechanical rejuvenation, restoring the system to an unaged state, similar to the thermodynamic effect of temperature increase²³. However, experimental studies on glassy polymers undergoing uniaxial deformation have shown enhanced molecular mobility before altering the underlying physical aging^24,25. Recent experiments on metallic glasses have phenomenologically demonstrated cold welding induced by ultrasound vibrations^26,27,28, which was attributed to localized processes enabling the formation of metallic bonds. Additionally, metallic polycrystalline alloys, also characterized by low molecular mobility and amorphous structure at the crystalline ___domain level, can undergo healing under well-controlled tensile high-cycle fatigue²⁹ by merging the crack propagation front at the surface of grain boundaries. Numerical simulations also revealed that stress-driven grain boundary migration occurs as a result of tensile and oscillatory shear deformations in nanocrystalline metals³⁰.

From the theoretical perspective of glasses, recent works by Sastry and coworkers^31,32,33 have shown in the athermal regime that mechanical oscillations can promote glass annealing at the molecular scale, and we recently proved a quantitative equivalence between temperature and shear rate in the molecular relaxation of glasses at T ≳ T_g subjected to steady-state shear flow²³. Building on these results, here we subject glassy systems to volume-preserving oscillatory deformations at varying amplitude and frequency. We identify specific deformation conditions under which molecular mobility can be accelerated without affecting the overall out-of-equilibrium nature of the polymer glasses, and we apply this principle to glasses in which an artificially-generated cylindrical crack can be healed. We find that oscillatory deformations can selectively accelerate the dynamics on the crack surface, promoting crack closure without leading to full melting of the glass in the bulk region. The theoretical understanding of different out-of-equilibrium dynamical regimes is a paradigm shift for the use of mechanical deformations to promote self-healing in polymer glasses.

Our framework involves polymer glasses with short chains to avoid entanglement effects. Beads belonging to the same chain interact via harmonic covalent bonds. Non-consecutive beads and beads from different chains interact attractively, corresponding to non-covalent interactions²³. The bead size is d, setting our length units (see the Methods section for a detailed model description). We explore a range of temperatures T from a fluid \((T > {T}_{g})\) to a glassy \((T < {T}_{g})\) state. We quantify the molecular mobility from the α-relaxation time τ_α obtained from the decay of the intermediate self-scattering function \({F}_{s}\left({q}^{*},t\right)\), with q^* and t corresponding to the length scale of the glassy cage and the simulation time, respectively, see Fig. S1. The structural relaxation time τ_α, in equilibrium and under deformation, is extracted by imposing that \({F}_{s}\left({q}^{*},{\tau }_{\alpha }\right)=1/e\).

Steady shear deformations

Before showing that oscillatory deformations can heal artificially-generated cylindrical cracks in glassy systems, we characterize the bulk α-relaxation dynamics under shear. The oscillatory deformation depends on the strain amplitude γ_A and frequency ω, as \(\gamma \left(t\right)={\gamma }_{A}sin\left(\omega t\right)\), which represents the strain applied to the material as a function of time. The combination of possible values within the set {γ_A, ω} is vast, even when restraining our simulations to a regime where the system temperature remains stable in the canonical NVT ensemble. To restrict and rationalize our design space, we first study the rheological behavior of polymer glasses under steady shear deformations. The stress-strain evolution picture is well established in amorphous solid materials³⁴: there is an initial linear response corresponding to solid elastic behavior at small deformations dominated by sterical cages, while plastic deformations at large strains eventually trigger a viscous-fluid behavior due to cage breaking. This transition is marked by a stress overshoot, placed at the yielding point γ_y,2. This picture is true for particle suspensions interacting sterically. However, in the present system, we consider linear chains interacting through attractive non-covalent bonds. As illustrated in Fig. 1a, where the stress component σ_xy is normalized by k_BT/d³ (with k_B being the Boltzmann constant), a first yielding point γ_y,1 emerges at small strain values as transient bonds break, allowing monomers to escape from interparticle attraction and partially dissipate energy^35,36,37. In addition, Fig. S2 shows that γ_y,1 is shear rate dependent. The existence of two yielding points allows us to identify three different regimes of γ_A, with ω remaining a free parameter. Focusing on systems at T/T_g ≥ 1.0, for which the polymer melt reaches equilibrium, we observe from Fig. 1b that γ_A defined in each regime induces different molecular mobility. Both axes are rescaled by the relaxation time at equilibrium \({\tau }_{\alpha }^{0}\), so that relaxation times are reported in comparison to their equilibrium value, and deformation rates are reported as dimensionless Deborah numbers. This number defines the ratio of the material’s characteristic relaxation time, in our case \({\tau }_{\alpha }^{0}\), to the characteristic flow time, which for oscillatory deformation is defined as γ_Aω. The molecular mobility is extracted from \({F}_{s}\left({q}^{*},t\right)\) under oscillatory deformation, see Fig. S3.

a Shear stress tensor σ_xy as a function of the strain γ for a shear rate \(\dot{\gamma }=0.01\) applied in the xy−plane. The first yielding point corresponding to physical bond breaking is referred to as γ_y,1, whereas γ_y,2 ≃ 2.6 indicates the yielding point corresponding to the transition from elastic solid to viscous-fluid behavior. The existence of these two yielding points delineates three regimes where all beads are frozen (linear regime, blue beads), partially relaxing (intermediate regime), or fully relaxing (nonlinear regime, red beads). b Normalized relaxation time \({\tau }_{\alpha }/{\tau }_{\alpha }^{0}\) as a function of the Deborah number \({\gamma }_{A}\omega {\tau }_{\alpha }^{0}\), for systems at different T subjected to oscillatory deformations. In all plots, the error bars indicate the standard deviation over configurations. c Energy per particle E/N as a function of accumulated strain γ_acc and γ_Aω for a system at T/T_g = 1.00. d General representation of the local mobility under oscillatory shear deformation, \({\tau }_{\alpha }\left({\gamma }_{A},\omega,T\right)\). Black arrows emphasize how the structural relaxation time changes with increasing γ_A and T of the system. e Snapshots at different accumulated strain γ_acc for a polymer melt system at T/T_g  = 1.00, subjected to an oscillatory deformation in the intermediate regime with γ_A = 0.4 and γ_Aω = 0.01. We also represent the displacement of the center of mass for the polymer chains. In this case, chains whose center of mass displaces equal to or more than the radius of gyration R_g are not found.

Full size image

Linear regime (γ _A < γ _y,1)

Oscillatory deformations have no effect on the α-relaxation, as indicated by the flat behavior of \({\tau }_{\alpha }/{\tau }_{\alpha }^{0}\simeq 1\). However, the potential energy per particle E/N as a function of the accumulated strain γ_acc = 4N_cycγ_A, shown in Fig. 1c and where N_cyc represents the deformation cycles, is slightly smaller than the corresponding energy for the system immediately before applying oscillatory deformations. Previous numerical simulations of oscillatory athermal quasistatic shear (AQS) deformations of dense amorphous samples^31,38,39 established the limit at small strain deformations, i.e., γ_A << γ_y,2, as a way to obtain better-annealed glasses. In particular, AQS assumes that the relaxation after deformation occurs on shorter timescales than ω⁻¹. Thus, deformations in the AQS approach depend on the deformation length scale encoded in the amplitude deformation. In our case, the deformation frequency ω introduces a timescale dependence that plays a key role in annealing glasses by accessing deeper energy states at small amplitude deformations^33,40. Although the study of this regime is beyond the scope of this investigation, we could anticipate that in the limit of poorly annealed glasses^33,40, and with decreasing ω (indicating that the deformation timescale is so large that it allows the system to relax the accumulated stress) our findings would align with the observations discussed in AQS simulations⁴¹.

Nonlinear regime (γ _A > γ _y,2)

Oscillatory deformations promote melting similarly to steady shear flow²³. In this regime, each oscillatory cycle induces a strain larger than the characteristic deformation that triggers the viscous-fluid transition. As a result, τ_α decreases below \({\tau }_{\alpha }^{0}\) as γ_Aω increases through ω. Actually, we find that the slope of \({\tau }_{\alpha }/{\tau }_{\alpha }^{0}\) tends to an asymptotic behavior as x^−0.8 by decreasing T and increasing the deformation. The exponent is related to the shear thinning behavior, which is characterized by an apparent viscosity that decreases with increasing shear rate^23,42. Additionally, E/N monotonically increases with γ_Aω, shown in Fig. 1c, whereas the characteristic plateau observed in \({F}_{s}\left({q}^{*},t\right)\) shortened, see Fig. S3.

Intermediate regime (γ _y,1 < γ _A < γ _y,2)

Intermediate strain amplitudes can speed up the molecular mobility without changing the underlying structure of the system. When the characteristic deformation time exceeds the structural relaxation time of the system, identified by Deborah numbers larger than 1, the dynamics accelerates until it reaches a plateau. Maintaining a constant strain amplitude, we detect that the plateau length increases as we approach T_g. On the other hand, τ_α falls to a plateau of lower height as γ_A approaches γ_y,2, see Fig. S4. As γ_Aω increases by increasing ω, the system approaches a critical value after which τ_α resembles that of a polymer fluid. This striking behavior is in agreement with the enhanced molecular mobility observed for polymer glasses under tensile deformations prior to flow^24,25. Similarly, regions of high local mobility have been detected for amorphous materials under oscillatory shear deformations^32,40,41,43. These regions of high local mobility forming shear bands were identified at γ_A > γ_y,2, and we investigate if a similar mechanism is observed in our systems.

One way to detect shear bands is by monitoring the particle energy E/N as a function of γ_acc. As discussed in ref. ⁴¹, E/N would exhibit a sharp upward shift because of the spontaneous onset of shear banding. In the intermediate and nonlinear regime, depicted in Fig. 1c and Fig. S5, we see that E/N evolves smoothly without showing abrupt changes. However, we note that oscillatory deformations with a strain amplitude just above the first yielding point induce strong energy fluctuations, see Fig. S5. Likewise, the evolution of E/N as a function of γ follows nearly identical trajectories in each cycle, reflecting the predominant elastic response of the material (see Fig. S5). Once γ_A > γ_y,1, the elastic behavior gradually weakens due to the breaking of non-covalent bonds, resulting in greater energy dissipation and system reorganization. Consequently, energy fluctuations decrease. Aiming to look for any signs of shear banding, we examine the mean-squared displacement \(\left\langle {r}^{2}\right\rangle\) for single beads and for the center of mass of each chain as a function of γ_acc (see Methods for more details). In Fig. 1e, we identify beads that undergo a displacement \(\left\langle {r}^{2}\right\rangle \ge {d}^{2}\), and chains whose center of mass undergo a displacement equal to or larger than the radius of gyration R_g, for a glassy system subjected to γ_A = 0.4 and γ_Aω = 10⁻¹. Contrary to what has been observed in refs. ^40,41,44, particles with high displacement do not develop a shear band. Furthermore, Fig. S6 shows that higher γ_A promotes the formation of a branched percolating cluster of particles with a high displacement, resembling the scenario observed under steady shear⁴⁵. Therefore, the enhanced molecular mobility we observe does not involve shear banding. This difference can be attributed to the fact that the deformation-induced dynamics are too fast to allow sufficient time for the formation of shear bands. Indeed, in ref. ⁴¹, the maximum shear rate explored is two orders of magnitude lower. Instead, the existence of these plateaus suggests that in the intermediate amplitudes regime, the oscillatory deformation breaks the microscopic cages responsible for the material vitrification. This mechanism, inducing local yielding, allows the system to remain globally arrested but to relax locally. Fig. 1d schematizes molecular dynamics (MDs) behavior under oscillatory deformation, expressed as \({\tau }_{\alpha }\left({\gamma }_{A},\omega,T\right)\).

Promoting self-healing with oscillatory shear deformations

The accelerated dynamics observed in Fig. 1 can be leveraged to heal an artificially-generated cylindrical crack in a glassy system via oscillatory deformations, as shown in Fig 2. As explained in the Methods section, polymer glass systems at T/T_g = 0.5 are prepared with a cylindrical crack of diameter D along the y axis, with walls preventing particles from entering the crack region. Once the glassy state is reached, we remove the cylinder and allow the system to evolve freely. We recognize that these configurations differ in shape and stress distribution from more realistic cracks naturally arising from material aging or damage. Nevertheless, inducing ad-hoc cracks enables precise and rigorous quantification of the mechanisms that trigger healing under mechanical oscillatory deformations, and is a first step to solve the broader problem of fracture healing in glassy materials. The first snapshot in Fig. 2a shows a polymer glass with a cylindrical crack of diameter D/d = 1.0. The evolution of the crack surface area is monitored over time, shown in Fig. S8 and Supplementary Video 1 and 2, respectively. For cylindrical cracks of diameter half the particle size, i.e., D/d = 0.5, the attractive interaction between different chains promotes inward flow. However, for D/d ≥ 1.0, the cracks are stable over time, as shown in Fig. S8a, because of the extremely long relaxation time, which is indicative of the very low molecular mobility that hinders any relaxation time in our observable timescale. Higher temperatures are therefore required, as shown in Figs. 2b and S8a, to close the crack.

a A cylindrical crack of diameter D/d = 1.0 in a polymer glass is repaired after several cycles of oscillatory deformations at γ_A = 1.0 and γ_Aω = 0.1. Beads are colored according to their depth within the material. Dotted circles highlight the initial position of the cylindrical crack. b T−D phase diagram indicating the percentage of closed samples by increasing the overall system temperature. The dashed box highlights the range of cylindrical crack diameters we explore under the action of an oscillatory deformation. c Structural relaxation time τ_α as a function of amplitude γ_A and frequency ω deformation. The error bars indicate the standard deviation over configurations. d γ_A−γ_Aω phase diagram for healed glasses at different D values of cylindrical cracks and varying strain amplitude and shear rate.

Full size image

Keeping T/T_g = 0.5 (glass regime), we focus on the highlighted range shown in Fig. 2b and apply oscillatory deformations. Firstly, we repeat the previous study discussed in Fig. 1b by computing the evolution of τ_α as a function of oscillatory deformations (see Fig. 2c). While we cannot estimate \({\tau }_{\alpha }^{0}\) because the system is out-of-equilibrium (see Fig. S1), we observe qualitatively the same behavior as we approach T_g. Oscillations with γ_A > γ_y,2 lead to the melting of the glass material, whereas for γ_y,1 < γ_A < γ_y,2, the dynamics can be locally accelerated by breaking microscopic cages. Additionally, we also observed that even in the glassy state, no shear band formation occurs (see Fig. S6), which may be attributed to the polymer chain connectivity of our model and the high oscillation frequencies. However, we also note that τ_α develops a stark maximum by increasing ω. This trend is already observed for T ≥ T_g (see intermediate regime in Fig. S4). The presence of this maximum following a sharp decrease in τ_α could indicate non-trivial annealing behavior of the glasses, similar to that observed near the yielding point in Kob-Anderson binary mixtures⁴¹. However, in Fig. S7, we have attempted to establish a connection between τ_α and the degree of molecular mobility by computing the distribution of per-particle Debye-Waller factor, but we do not see signs of any annealing. Instead, we observe that this peak marks the transition from a well-defined unimodal distribution to a distribution with a long tail. Returning to the mechanical self-healing, Fig. 2a displays the closure of a crack of D/d = 1.0 after N_cyc = 72 deformation cycles. Fig. 2d displays a set of γ_A−γ_Aω phase diagrams as a function of the initial D. These diagrams illustrate the percentage of samples in which oscillatory deformations in the intermediate regime successfully closed the crack. While for D/d = 1.0 we consistently achieve crack closure, as D increases, the efficiency of oscillatory deformation decreases. However, we note that this decline begins at high values of ω for D/d = 1.5, and becomes more pronounced for D/d = 2.0. This trend is rationalized by the fact that the oscillation rate of the deformation imposed by ω overcomes the characteristic particle diffusion rate. This, along with the fact that the crack diameter is large enough to suppress particle interactions across the crack, favors particles following the same trajectory in each deformation cycle. Thus, an increase in γ_A is required to bring particles in contact across the crack. Fig. S8b shows the evolution of the crack surface area over time for an initial cylindrical crack of diameter D/d = 1.5 under oscillatory deformations. In addition, we also explore the case in which the system has been previously subjected to an affine oscillatory deformation, taking the system configuration at the maximum deformation and then letting it evolve (see Supplementary information for additional details). Notably, affine deformations and the shape variation of the crack under large deformations are not sufficient to initiate particle flow, and dynamic oscillatory deformations are needed, as we further prove in Fig. S8. At least qualitatively, this healing mechanism seems to also apply to less artificial crack configurations, as we show in Fig. S9 where the glassy material is punctured by indenting it with a sharp object before applying deformations. Still, a systematic extension of our findings to realistic fracture configurations is beyond the scope of this work and should be an important focus of future studies.

One may wonder if we genuinely induce self-healing by accelerating the local dynamics exclusively around the crack. This is a legitimate doubt because deformations are often associated with the breakdown of the material’s structure, inducing mechanical rejuvenation in amorphous materials^23,46,47. In addition to showing that this amplitude regime precedes the material’s yielding point γ_y,2, Fig. 3 shows that indeed cracks can heal purely by accelerated local dynamics, while the bulk remains globally glassy and the final mechanical properties of the healed glass are comparable to those of the pristine equilibrium glass. We identify particles on the surface of the crack and in the bulk, i.e., far away from the crack, and compute their mobility during deformation through the Debye-Waller factor \(\left\langle {u}^{2}\right\rangle\), which corresponds to the value of the mean square displacement at local minima of its logarithmic derivative^48,49, as shown in Fig. S1e. Fig. 3a displays results corresponding to a polymer glass with a crack of diameter D/d = 1.0, subjected to {γ_A, ω} = {1.0, 0.1}, equivalent to γ_Aω = 0.1. The crack surface and particles on the surface are represented in the snapshots in light gray and dark gray colors, respectively, while particles in the bulk are omitted to enhance the clarity of the snapshots. Furthermore, we depict the evolution of \(\left\langle {u}^{2}\right\rangle\) as a function of the time t alongside the corresponding mobility at different T in the absence of deformation, indicated by a horizontal line. While \( < {u}^{2} > \) slightly increases for particles in the bulk, we note that the dynamics for surface particles aligns with that observed for a bulk equilibrium system at T/T_g = 1.0. As the crack closes, the surface dynamics converges to the bulk dynamics, at a value slightly higher than the equilibrium dynamics at T/T_g = 0.50 because of the oscillatory deformations. It is important to note how the particles initially located on the surface remain localized around the region where the crack was. This indicates that local melting induces an inward flow of material, filling the damaged space. In the case where the crack does not close, as depicted in Fig. 3b, the surface and bulk dynamics do not match.

Local mobility \(\left\langle {u}^{2}\right\rangle\) for particles on the crack surface and bulk as a function of time t for polymer glasses subjected to oscillatory deformations of γ_A = 1.0 and γ_Aω = 0.1, with a cylindrical crack of diameter (a) D/d = 1.0 and b D/d = 2.0. Horizontal lines represent the bulk local mobility as a function of T at mechanical rest. Circles highlight the times at which the snapshots were taken, represented on the right. The crack region is represented in light gray, whereas dark gray beads correspond to beads on the crack surface. c The healed system at ω = 0.1 and different γ_A values is subjected to steady shear flow at shear rate \(\dot{\gamma }=0.01\). The resulting stress tensor is compared with the stress curve of the glass at rest, showing that mechanical properties are recovered.

Full size image

Finally, we assess the mechanical properties of the material immediately after closing the crack. In Fig. 3c, we compare the mechanical response under steady shear of the healing glass material with that obtained from a polymer glass in which no crack is present. At low-strain deformations, the first yielding point γ_y,1, associated with physical bond breaking, is less apparent in the stress curves of the healed material. This issue is attributed to the slightly accelerated bulk dynamics, see Figs. 3a and  S10. However, the position of γ_y,2 and the stress overshoot \({\sigma }_{xy}({\gamma }_{y,2})\), indicating the opposing resistance of the material to break and flow, perfectly match. Thus, our findings suggest that the application of small oscillatory deformations has the potential to induce complete self-healing in glassy polymer materials.

In recent decades, oscillatory deformations have emerged as an avenue for exploring new frontiers in the fundamental physics of amorphous materials. This includes enhancing annealing^39,40 and inducing memory^50,51 in amorphous glasses. Here, we used oscillatory deformations to induce self-healing in glassy polymer materials by overcoming their intrinsically low local mobility. The enhanced mobility would enable the use of intrinsic healing mechanisms such as dynamic covalent bonds or supramolecular bonds that strongly depend on the local dynamics⁹.

First, we have studied the local mobility in bulk model glasses without damage. We have identified a range of oscillatory deformations, placed above the strain deformation γ_y,1 corresponding to breaking transient bonds, and below the critical strain γ_y,2 value, which leads to melting, and where the molecular mobility can be meticulously accelerated without globally modifying the underlying structure or mechanical properties of the system. In this range of deformations, which we define as the intermediate regime, particle dynamics are locally accelerated as microscopic cages responsible for the glassy state are broken, while overall dynamics remain arrested. We have shown that the structural relaxation time exhibited a plateau likely due to the acceleration of local dynamics which does not lead to a fully liquid behavior, since the strain amplitude does not exceed γ_y.2. Furthermore, carefully observing the particle displacement, we have found that the presence of a plateau in relaxation time is not accompanied by the formation of shear bands in our simulations, and τ_α dramatically decreases at critical ω values that are strain amplitude-dependent. We have then demonstrated that this range can be exploited to induce mechanical self-healing under highly controlled conditions. In particular, we have studied damaged polymer glasses at T < T_g where the low molecular mobility hinders any relaxation and makes artificially-generated cylindrical cracks stable at equilibrium at any observable timescale. By monitoring the local mobility of the particles, \(\left\langle {u}^{2}\right\rangle\) we were able to observe that the dynamics under oscillatory deformations is more prominently accelerated for the particle located at the crack surface. In oscillatory deformations within the intermediate regime, this leads to a mass flow into the crack within the simulation timescale, while the bulk dynamics is largely unchanged. Furthermore, stress-strain curves for healed samples closely match the pristine material.

Since glassy polymers are out-of-equilibrium materials, their rheological response depends on the aging history and preparation protocol. In the past, it has been studied that more stable glasses are more brittle, making the yielding point a spinodal instability characterized by a sharp discontinuous stress jump⁴⁴. It would be relevant in the future to understand if local mobility controlled by oscillatory deformations depends on the glass stability. Likewise, cracks induced by realistic mechanical deformation are an important next step to investigate since the perfect cylindrical cracks investigated in our study might not represent the structure and distribution of stresses of cracks naturally emerging from material aging or damage. In this scenario, local dynamics and stability would be affected by the highly heterogeneous stress concentrations and irregularity of the surface created by the mechanical deformation^25,52. Furthermore, mechanical deformations could promote the presence of multiple small cracks, giving rise to a yet unexplored scenario that is worth investigating. Nevertheless, our results on polymer glasses, along with previous observations of healing in metal alloys²⁹, nanocrystalline metals³⁰ and metallic glasses^26,28 under mechanical deformations, point to the fact that oscillatory deformations represent a general strategy to induce mechanical self-healing in materials with extremely low local mobility. Previous computational results showing equivalent equilibrium glassy dynamics at the surface of thin films and the curved interface of spherical nanoparticles^49,53,54 suggest that this mechanism might be general regardless of the interface nanostructure, though further studies are needed, given the experimentally known complexity of glassy dynamics at realistic, rough interfaces^55,56. In this respect, ultrasonic vibrations of frequency \(f=\left[0.02-7\right]\times 1{0}^{4}\) Hz and strain amplitude of γ_A = 10⁻²% applied to metallic glasses appear to be particularly promising to study mechanical self-healing²⁶. With the assumption of an experimental structural relaxation time of 100 s at T_g⁵⁷, these deformations correspond to Deborah numbers \(D{e}_{exp}\in \left[1{0}^{2}-1{0}^{4}\right]\), similar to the value De_sim = 10³ of our simulations if we use a fairly common estimate for a computational T_g when τ_α = 10⁴t_LJ. Therefore, we hope that the broader theoretical discussion and the explicit demonstration of molecular self-healing without volume-changing deformations presented here in controlled conditions can stimulate further theoretical studies in more realistic conditions and targeted experimental activities toward the study of self-healing in polymer glassy materials.

Modeling

We perform MDs simulations of fully flexible linear chains of beads linked by harmonic springs. Nonbonded monomers belonging to the same or different chains interact with a truncated LJ potential defined as

where d is the particle diameter, which sets the unit of length, ϵ controls the energy scale, and U_cut ensures that \({U}_{LJ}\left(r=2.5d\right)=0\). Defining m as the mass of the particles, the time units are defined as \(t=\sqrt{m{d}^{2}/\epsilon }\). On the other hand, chemical bonds between connected monomers are modeled by

where k_b = 555.5ϵ/d² is the spring constant and r₀ = 0.97d is the equilibrium bond length. Henceforth, all quantities are expressed in terms of reduced LJ units, i.e., ϵ = 1 and d = 1, with unit monomer mass m and Boltzmann constant. The reduced units can be mapped onto physical units relevant to generic nonequilibrium fluids, by taking MD time, length, and energy units as corresponding roughly to 2 ps, 0.5 nm, and 3.7 kJ/mol, respectively.

We consider systems of N_c = 500 chains of M = 20 monomers, thus avoiding the entanglement, and doing a total of N = 10⁴ monomers. All simulations were performed using LAMMPS simulation package⁵⁸, considering a simulation time step δt = 0.002. Results are averaged over 5 independent simulations.

Sample preparation

Initially, an equilibration process was performed by considering polymer chains enclosed in an orthogonal cubic box of size L with periodic boundary conditions. First, 10⁷ simulation steps were performed in the NPT ensemble by employing a Nosé-Hoover thermostat and barostat with T = 1.0 and \(\left\langle P\right\rangle=0\) allowing full correlation loss of the end-to-end vector of the polymer chains. Next, a quenching was made in the range of temperatures from T/T_g = 0.50 to T/T_g = 1.38, maintaining \(\left\langle P\right\rangle=0\) during 10⁷ simulation steps. In both equilibration procedures, the simulation box was allowed to fluctuate in an isotropic fashion in all three directions of space.

Static and dynamic characterization

Once the system was equilibrated to the required temperature, the volume was fixed by performing simulations in the NVT ensemble by using a Nosé-Hoover thermostat at the required temperature, and static and dynamic properties were computed to characterize the polymer system as a function of T. In particular, the static structure factor was computed as \(S(q)=\langle \frac{1}{N}{\sum }_{ij}{e}^{-i{{{\bf{q}}}}\cdot ({{{{\bf{r}}}}}_{i}-{{{{\bf{r}}}}}_{j})}\rangle \,\), where q is the wave vector, r_i indicates the position of the ith particle, and \(\left\langle \cdots \,\right\rangle\) denotes an ensemble average. Likewise, the dynamic of the system was quantified through the mean-squared displacement \(\langle \Delta {r}^{2}(t)\rangle=\langle \frac{1}{N}\mathop{\sum }_{i=1}^{N}{[{{{{\bf{r}}}}}_{i}(t)-{{{{\bf{r}}}}}_{i}(0)]}^{2}\rangle\) and the self-intermediate scattering function \({F}_{s}({q}^{*},t)=\langle \frac{1}{N}\mathop{\sum }_{i}^{N}{e}^{i{{{{\bf{q}}}}}^{*}\cdot [{{{{\bf{r}}}}}_{i}(t+t{\prime} )-{{{{\bf{r}}}}}_{i}(t{\prime} )]}\rangle\), computed on the wave vector q^* = 2π/l corresponding to the main peak of the static structure factor. The relaxation time τ_α was extracted by imposing \({F}_{s}\left({q}^{*},{\tau }_{\alpha }\right)=1/e\).

Shear deformation

After equilibration, we perform two kinds of simulations: (i) steady shear deformations at a fixed shear rate \(\dot{\gamma }\); (ii) oscillatory shear deformations by applying a sinusoidal deformation in the xy plane defined as \(\gamma \left(t\right)={\gamma }_{A}sin\left(\omega t\right)\), where γ_A is the strain amplitude and ω is the frequency. In both cases, deformations are performed at constant volume and temperature, and hence, the SLLOD equations of motion were used along with the thermostat, in combination with Lees-Edwards boundary conditions. During the steady shear deformation, we monitor the stress tensor component \({\sigma }_{\alpha \beta }=\frac{1}{A}\langle {\sum }_{i\; > \; j}\, {f}_{ij}^{\alpha }{r}_{ij}^{\beta }\rangle\), where \({f}_{ij}^{\alpha }\) is the α-component of the force with respect to interaction defined in eq. (1) and eq. (2), \({r}_{ij}^{\beta }\) corresponds the β-component of the distance vector between particles i and j. In this case, the steady shear deformation was applied in all three directions of space and subsequently averaged. On the other hand, during the oscillatory deformation, we first apply an accumulated strain γ = 4N_cycγ_A = 100%, where N_cyc corresponds to the number of deformation cycles. Then, we compute the self-intermediate scattering function on the wave vector q^* corresponding to the main peak of the static structure factor before applying any deformation. Since oscillatory deformations were applied on the xy plane, \({F}_{s}\left({q}^{*},t\right)\) was computed excluding q_x components, see Fig. S3. Furthermore, we computed the mean-squared displacement per particle, defined as \(\Delta {r}^{2}{\left(t\right)}_{i}={\left[{{{{\bf{r}}}}}_{i}\left(t\right)-{{{{\bf{r}}}}}_{i}\left(0\right)\right]}^{2}\), excluding the x-component.

Inducing self-healing

To study self-healing, we create a crack on polymer glasses with T/T_g = 0.50 and N monomers. The equilibration previously described above is performed in the presence of a cylinder placed at the center of the simulation box of length L_y and diameter D. Monomers interact with the cylinder wall through the Weeks-Chandler-Andersen potential:

Then, the cylinder is removed, and NVT simulations and oscillatory deformation simulations are independently performed during 10⁷ up to 3 × 10⁷ simulation steps. In both cases, the volume of the crack is monitored by employing the surface area algorithm⁵⁹. During both simulations, we compute the mobility \(\left\langle {u}^{2}\right\rangle\) of the monomers as a function of the distance from the center of the simulation box, as well as a function of time. Since the mobility is statistically very noisy, we apply a Savitzky-Golay filter on windows of 21 points and equations of 2 order.

Mechanical properties after the self-healing

Glassy polymer systems, on which the crack was closed, are subjected to steady shear flow at \(\dot{\gamma }=0.01\), and the evolution of σ_xy as a function of \(\gamma=\dot{\gamma }t\) is compared with the corresponding σ_xy for the glass at rest.