**2. Materials and Methods**

As a model system, we consider void nucleation at *θ*-particles in an FCC Al matrix. *θ* is the thermodynamically stable intermetallic phase of the Al-Cu system and is commonly observed in Al-Cu-copper alloys (e.g., 2xxx series) in the overaged condition [34]. *θ*-particles have a composition of Al2Cu and a body-centered tetragonal C16 crystal structure. They are incoherent with the matrix and typically adopt plate-like geometries [35]. However, for

simplicity in this work we will use a spherical precipitate geometry. While there is evidence that voids may nucleate at *θ*-particles [36], we emphasize that we are using the *θ*-Al system as a model incoherent precipitate system with the goal of gaining general insight into the micromechanics and kinetics of void nucleation.

MD simulations were performed using LAMMPS [37] with the Al-Cu angulardependent interatomic potential of Apostol and Mishin [38]. The angular-dependent potential framework is an enriched version of the embedded atom method that enables incorporation of angular-dependent interactions. These interactions enable the potential to capture the lattice constants, anisotropic elastic moduli, and formation energy of the *θ*-phase with reasonable accuracy. Our simulation cell geometry is shown in Figure 1a; we initially inserted an incoherent spherical particle with a radius of *R* = 50 Å into a pure Al lattice using the zero-temperature lattice parameters predicted by the potentials. The Al lattice and *θ*-particle are oriented so that their unit cell axes are aligned with the simulation box. The *c*-axis of the *θ*-lattice is oriented in the *z*-direction of the simulation cell. Periodic boundary conditions were used in all directions with a 200 Å cubic simulation cell. The sequence of each simulation is shown in Figure 1b. During the relaxation stage, we used an *NPT* (constant number of atoms *N*, pressure *P*, and temperature *T*) ensemble and simulated 2 ps at the chosen temperature and zero hydrostatic stress. Subsequently, during the ramping stage we ramped the hydrostatic stress up to the target value *σ<sup>H</sup>* over a duration of 23 ps. This duration was chosen empirically; if the stress was ramped too quickly we observed "premature" void nucleation, likely because of stress spikes resulting from imperfect performance of the barostat. See Appendix A for additional information. Finally, during the holding stage, the hydrostatic stress was held constant for the duration of the simulation until the particle completely debonded from the matrix or the simulation terminated after 1-week of wall time. All simulations used a thermostat damping parameter value of 0.01 ps, barostat damping parameter value of 1 ps, and a time step size of 0.001 ps.

**Figure 1.** Simulation details. (**a**) Snapshot showing periodic simulation cell with a spherical *θ*-particle loaded hydrostatically; (**b**) time history of hydrostatic stress from a sample simulation at *T* = 400 K and *σ<sup>H</sup>* = 6.0 GPa with simulation stages marked. When the particle completely delaminates the applied stress can no longer be sustained, causing the precipitous drop at the end.

We note that while the MD barostat controls the average (virial) stress state in the simulation cell, the local stress state may vary. In fact, we expect there to be variation because the elastic constants between the matrix and particle differ, i.e., this is an Eshelby inhomogeneity problem [39]. Furthermore, the particle images resulting from periodic boundary conditions will interact with each other, further complicating the stress field. These effects are quite complex, especially given the anisotropic nature of the C16 *θ*-phase. For simplicity, in our analysis of the data we assume that the applied hydrostatic stress *σ<sup>H</sup>*

dominates the delamination behavior; our successful thermal activation analysis below justifies this assumption. We note that Pogorelko and Mayer have analyzed the spatially varying stress field near a second-phase particle under uniaxial loading and its influence on the delamination process [26,27,29].

Results from a total of 290 MD simulations are reported in this study at temperatures ranging from 200 to 400 K and stresses in the range of 5.7 to 7.2 GPa (the precise stress range differed for each temperature). In most cases, 10 simulations with different thermalization histories (initial atomic velocities) were performed at each stress-temperature condition. In many simulations, we observed nucleation of dislocations at the particle interface. To enable efficient detection of the appearance of Shockley partial dislocations in our simulation cell, we exploited the fact that atoms situated in stacking faults (e.g., produced by a Shockley partial dislocation) appear as HCP atoms when analyzed via common neighbor analysis (CNA) [40]. Hence, by simply monitoring the number of "HCP" atoms in the simulation cell *NHCP*, we could identify when a dislocation appeared. We note that there was always a small, non-zero number of HCP atoms detected, due to thermal noise in the lattice and the imperfect detection capacity of CNA. To prevent false detection of a dislocation, we established a threshold value for the appearance of a dislocation, *N<sup>d</sup> HCP*, based on empirical analysis of our data. We set this threshold at *N<sup>d</sup> HCP* = 40 HCP atoms, so if *NHCP* > 40 we "detected" appearance of a dislocation. Furthermore, we applied a moving average to the raw *NHCP* vs. time data using a window of width 0.005 ps. This served to smooth out the data a bit and remove spurious spikes in *NHCP* which did not lead to a sustained increase in *NHCP* over time (as was expected if a dislocation had nucleated and remained in the system). This approach to dislocation detection was validated by manually analyzing several datasets in OVITO Pro [41]. All simulation snapshots were produced using OVITO Pro [41].
