Grain Boundary Effects in Dealloying Metals: A Multi-Phase Field Study

A multi-phase field model is employed to study the microstructural evolution of an alloy undergoing liquid dealloying. The model proposed extends upon the original approach of Geslin et al. to consider dealloying in the presence of grain boundaries. The model is implemented using a semi-implicit time stepping algorithm using spectral methods, which enables simulating large 2D and 3D domains over long time-scales while still maintaining a realistic interfacial thickness. The model is exercised to demonstrate a mechanism of coupled grain-boundary migration to maintain equilibrium contact angles with this topologically-complex solid-liquid interface during dealloying. This mechanism locally accelerates dealloying by dissolving the less noble alloy metal from (and rejecting the more noble metal into) the migrating grain boundary, thereby enhancing the diffusion-coupled-growth of the liquid channel into the precursor. The deeper corrosion channel at the migrating grain boundary asymmetrically disrupts the ligament connectivity of the final dealloyed structure, in qualitative agreement with published experimental observations. It is shown that these grain boundary migration-assisted corrosion channels form even for precursors with small amounts of the dissolving alloy species, below the so-called \textit{parting limit}


INTRODUCTION
Dealloying of multicomponent metals is known to generate highly topologically complex structures with micron-or nano-scale ligaments and voids [2][3][4] . Interest in this process stems from its use as a method to produce nanoporous metals at large scales. Dealloyed nanoporous metals have been shown to possess remarkable catalytic properties, large capacitance, and high mechanical strength when compared to their bulk counterparts [5][6][7][8] . However, these desirable properties can be deleteriously affected by the presence of grain boundaries in the metal, which are known to affect the ligament size and connectivity of these dealloyed structures 7,[9][10][11][12][13][14] .
In general, systems relevant to dealloying processes can be described as a ternary, wherein a metal alloy precursor made up of species 'A' and 'B' is subjected to corrosion by a fluid dealloying agent 'C.' The thermodynamics of these species' interactions is such that they are all relatively soluble in one another, except species A and C which are relatively insoluble. This results in the selective dissolution of B into C, while A organizes into a topologically complex solid structure. Some examples of porous metals synthesis via dealloying include aqueous solution dealloying (e.g. Ag dealloying from AuAg in acid 15,16 ), liquid metal dealloying (e.g. Ti dealloying from TaTi in liquid Cu 10 ), and, more recently, molten salt dealloying (e.g. Cr dealloying from NiCr in molten salt 17 ).
The key mechanisms by which nanoporous microstructure formation proceeds have been identified by means of numerical and experimental investigations. Early work by Erlebacher et al. 15 illustrated the important role of spinodal driving forces and diffusion processes within the solid-liquid interface in reorganizing the non-dissolving A-species into ligaments. This promotes dissolution of the AB precursor into a bicontinuous structure, made up of two continuous, interpenetrating phases of A-rich solid and C-rich liquid. A phase field model for liquid metal dealloying (LMD) of TaTi precursors in liquid Cu was later developed by Geslin et al. 1 based on a parameterization informed by experimental observations. Their model successfully predicted the early stages of dealloying and A-ligament formation, and elaborated on this mechanism of interface-diffusion-coupled-growth of the A (Ta) rich solid and C (Cu) rich liquid into the AB (TaTi) precursor.
Although the microstructure of the solid alloy precursor is expected to play an important role in this dealloying, the detailed mechanisms underlying how features such as grain boundaries may alter this process are less understood. Experimental studies have shown that grain boundaries locally alter dealloying, resulting in a corroded microstructure that is topologically-and chemically-distinct from the bulk lattice 7,9,[12][13][14] . In one example, liquid metal dealloying (LMD) experiments by McCue et al. 10 on TaTi precursors showed that corrosion at grain boundaries is significantly deeper and generates thick Ta-blocks that asymmetrically disrupt the ligament connectivity of the final porous structure. More recently, experiments by Joo et al. 11 investigated the early stages of this phenomenon in FeNi precursors dealloyed by liquid Mg. They also showed that infiltration by the dealloying agent at the grain boundary generates thick dealloyed plates that are asymmetrically attached to one of the par-ent grains while being separated from the other grain by the liquid melt. Interestingly, at different locations along a corroded grain boundary, the connectivity of these dealloyed plates was seen to switch back and forth between parent grains, resulting in a "wavy" corrosion channel 11 .
Mechanistic understanding of these experimental observations 10,11 remains incomplete. It has been proposed that the excess energy of grain boundaries (GBs) accelerates dealloying through a mechanism of GBwetting by the liquid phase 11 . However, the enhanced kinetics at the GB (both interface mobility and solute diffusivity), should also influence this process by locally altering the diffusion-coupled-growth of the interpenetrating solid and liquid phases. It is critical to understand these mechanisms of grain boundary dealloying corrosion to optimize nanoporous materials synthesis.
The current work aims to understand how the presence of grain boundaries affects local dealloying in a model LMD system. The phase field model by Geslin et al. 1 is extended to model an AB bicrystal precursor corroding in liquid C, using the multi-phase field formalism of Steinbach and Pezzolla 18 . The multi-phase field model uses an advanced semi-implicit time integration procedure by Badalassi et al. 19 , which allows for accessing large lengthand time-scales while still maintaining realistic interfacial thicknesses. 2D and 3D phase field simulations are used to investigate grain boundary effects on dealloying morphology, and numerical predictions are discussed in the context of experimental findings by McCue et al. 10 and Joo et al. 11 . In addition, we discuss how these findings may relate to grain boundary accelerated dealloying seen in metals in molten salt environments 14,20 .

RESULTS
Dealloyed Microstructure at Grain Boundaries Figure 1 presents a 2D multi-phase field simulation of a liquid metal dealloyed bicrystal precursor at various time steps. The simulation uses the thermodynamic and kinetic parameterization of TaTi (hereafter AB) exposed to liquid Cu (hereafter C), and grain boundary energies, mobilities and diffusivities set to realistic values as described in the Methods section. The initial simulation condition and snapshots at 23 and 72 µs are shown in Figs. 1 (a)-(c), respectively. These simulations employ a 2D grid of 2048 × 2048 pixels (409.6 × 409.6 nm 2 ), initially containing a bicrystal precursor (made up of initial A-concentration c A,0 = 40% and B-concentration c B,0 = 1 − c A,0 ) in contact with a liquid phase (made up of pure C). Each image is colored by the local concentration c A , and the gray band tracks the grain boundary. For more details on the simulation cell setup, see the Methods section. Videos of all simulations are provided in the Supplementary Materials.
As shown in Figure 1(b) and 1(c), the liquid dealloying agent penetrates the precursor in a bicontinuous manner leading to the formation of an A-rich topologically-complex solid structure. The A-rich solid forms a highly connected ligament network, except at the grain boundary, where a significantly deeper corrosion channel forms and separates the parent grains. Interestingly, the grain boundary appears to couple with the dealloying process by migrating from its original position and imparting a locally-distinct microstructure. Specifically, the grain boundary begins to migrate towards the left and accelerate local corrosion by dealloying B from deep inside (and rejecting A deep into) the bulk lattice through which it sweeps. An A-rich wake forms from the migrating grain boundary, which locally passivates the growing (right) grain at its solid-liquid interface. We label this flat, Arich structure as the grain boundary migration-assisted dealloying (GBMD) region in Figure 1(c). The resulting morphology is a bicontinuous network of A-rich ligaments that are disrupted by the larger A-rich block formed by the migrating grain boundary.
3D simulations are necessary to explore the dealloying corrosion morphology across the solid-liquid interfacial plane, both near and away from the grain boundaryliquid line. Figure 2 shows successive snapshots from a 3D multi-phase field simulation of a bicrystal precursor with initial A-concentration c A,0 = 20% undergoing liquid metal dealloying. In the top row of Figure 2, 3D views are shown, where only the interfaces (solid-liquid and grain boundary) are filled in and colored by the concentration of species A. In the bottom row of Figure 2, the topological views looking down at the solid-liquid interface for each snapshot are shown, colored by the corrosion depth. In the topological views, the location of the diffuse grain boundary-liquid line is shown by the dark grey band. These simulations employ 256 × 512 × 512 voxels (51.2 × 102.4 × 102.4 nm 3 ) in the X-, Y-, and Z-directions respectively.
In the very early stages of the dealloying in Figure 2(b) and 2(e), A-rich solid nodes begin to set up and the grain boundary begins to migrate, as was seen in early stages of 2D simulations. These 3D simulation snapshots provide additional insight, showing that, away from the grain boundary, these A-rich nodes form as three-dimensional mounds. Meanwhile, along the migrating grain boundary a single two-dimensional ridge forms. Examination of Figure 2(e) reveals that the grain boundary is actually migrating in two opposite directions simultaneously -in some locations it has migrated in the positive Z-direction, while in other locations it has migrated in the negative Z-direction.
As the corrosion continues, the A-rich mounds away from the grain boundary grow and form threedimensional A-rich ligaments. At the grain boundary, the two-dimensional A-rich ridge persists while corrosion channels rapidly penetrate the precursor, separating each grain's ligament network from one another. As was the case in 2D simulations, these GBMD corrosion channels form in the direction the grain boundary is locally migrating.

Grain Boundary Migration-Assisted Dealloying Mechanism
The multi-phase field simulation results herein are used to understand in detail the mechanism underlying the coupling between dealloying and grain boundary migration. This mechanism of grain boundary migrationassisted dealloying, or "GBMD," is driven by the relative energetics of the grain boundary and solid-liquid interface and is reinforced by the solute flux along these interfaces at their junction. Figure 3 illustrates the GBMD mechanism for a bicrystal with c A,0 = 35%. (For a full description of this dealloying process in the absence of grain boundary effects, the reader is referred to Geslin et al. 1 , and only the salient features and grain boundary effects are described herein). Generally throughout the dealloying process, the solid-liquid interface becomes A-rich as B is selectively dissolved from the precursor. Since diffusive transport within the bulk lattice is negligibly small, dealloying is confined to within the solid-liquid interface, except at the grain boundary-liquid junction where diffusion along the grain boundary promotes B-dissolution / A-rejection from / to deeper within the precursor solid alloy. This is shown in Figure 3(a) at an early stage of the dealloying process, where it can be seen that the dissolution is relatively uniform across the solid liquid-interface except at the grain boundary-liquid junction where it is slightly faster. We note here that, in the present simulations, the grain boundary does not simply melt faster than the bulk solid, as one might assume. Rather, since the grain boundary dealloys faster, it builds up in species A, which actually raises the local melting temperature, as it is higher for species A than B.
At this early stage, the dealloying remains symmetric -from left to right -about the vertical grain boundary. Therefore, the dihedral angles in each grain also remain symmetric (θ 1 ≈ θ 2 ), with their magnitude determined by the excess interfacial grand potential of the grain boundary, γ gb , and of the solid-liquid interface, γ sl . This is schematically shown in Figure 3(b). The relative energetics of each interface depend on the prescribed values of the grain-boundary and solid-liquid interfaces for the pure elements (σ gb and σ sl in Eq. 1) used to parameterize the multi-phase field model, as well as the local chemistry that modifies these interfacial free energies relative to their values for the pure elements. With regards to chemistry, the high concentration of solid-A mixing with liquid-C is seen to significantly increase the solidliquid interfacial energy in the model, due to their large miscibility gap. This leads to a large value (near π) for the liquid dihedral angle, θ 3 , as seen in Figure 3(a) and 3(b).
The dealloying front quickly becomes B-depleted / Arich to the point that the solid-liquid interface develops a composition that is within the unstable region of the chemical free energy in Eq. 4. This leads to spinodal instabilities for unmixing within the interface, causing segregation of A-colonies from the penetrating C. Importantly, interfacial solute diffusion processes are ac-tive enough to enable this A-C segregation. As these A-rich structures continuously set-up and coarsen laterally along the solid-liquid interface, the symmetry about the grain boundary eventually breaks. One example of this is shown in Figure 3(c), where the topology of the Arich solid structure has "tilted" slightly so that it slopes down-right at the grain boundary-liquid junction. To reestablish force balance, the grain boundary-liquid junction and grain boundary migrate (towards the right) until θ 1 ≈ θ 2 , as illustrated in Figure 3(d). We note here that the direction of this asymmetry is random, and grain boundary migration is equally likely to initiate towards the right or the left.
Once this migration is initiated, it is generally reenforced through two coupled mechanisms. First, as the grain boundary moves with the triple junction, it migrates through the bulk lattice, continuously providing new, non-dealloyed regions along which to reject A and dealloy B. This results in the A-rich wake following the grain boundary migration path, seen again here in Figure  3(e). A snapshot of this A-rejection into the grain boundary is shown in Figure 3(f), where the local flux of species A is represented by the red streamlines and the gray regions show the interfaces. It can be seen that A -being immiscible with the liquid C and immobile in the bulk lattice -is confined to diffuse along the solid-liquid and grain boundary interfaces. Along the solid-liquid interface, spinodal driving forces drive the "uphill" diffusion of A towards adjacent A-rich structures. As discussed in Geslin et al. 1 , this is essential to the diffusion-coupledgrowth of interpenetrating A-rich solid and C-rich liquid phases into the AB-solid precursor -where the rate of diffusion-coupled-growth is controlled by the rate at which species A can be "evacuated" from the infiltrating C-liquid channel. An essential feature of GBMD is that the grain boundary provides an additional pathway along which to reject species A, thereby locally enhancing this diffusion-coupled-growth.
The second reenforcing mechanism, grain boundary straightening, occurs simultaneously.
As the grain boundary-liquid junction drags the grain boundary, significant curvature is induced at locations deeper inside the solid, as shown in Figure 3(g). If the grain boundary interface mobility is sufficient, capillary forces will move the grain boundary to reduce its curvature. In the simulations, the high interfacial energy of the solidliquid interface that arises when it is enriched by species A enforces the grain boundary-liquid triple point to be constrained to maintain force balance. As a result, the grain boundary can only straighten by migrating into non-dealloyed regions of the metal, allowing it to take on more A from the solid-liquid interface. (It is noted here that the "continuous-gradient" boundary conditions, discussed in Methods, essentially pin the grain boundary to its initial position at the bottom edge, and in a real system the grain boundary would migrate to a greater extent than the simulations show here). shown schematically, (c) grain boundary-liquid junction migration to achieve force balance during in the presence of an asymmetric corrosion front seen in the phase field simulation, and (d) illustrated schematically, (e) later stages of grain boundary migration driven by (f) the flux of species A out of the solid liquid interface down into the grain boundary, and (g) due to curvature forces on the grain boundary.

Grain Boundary Dealloying and the Parting Limit
The degree to which this GBMD mechanism influences the corrosion morphology and depth depends strongly on the initial composition of the precursor. This is illustrated in Figure 4, which shows dealloyed bicrystal precursors with varying A-content, increasing from c A,0 = 20% in Figure 4(a), c A,0 = 30% in Figure 4 Precursors with relatively low A-content, c A,0 = 20% in Figure 4(a) and c A,0 = 30% in Figure 4(b), dealloy to form somewhat lamellar structures with alternating A-rich solid and C-rich liquid phases. The diffusioncoupled-growth of each phase is generally sustained since there is a moderately small amount of A to remove from the infiltrating C-channels. The grain boundary accelerates this effect, via the GBMD mechanism illustrated in As the concentration is increased to c A,0 = 45% in Figure 4(c), a more-connected network of thicker A-rich ligaments forms. Since a higher nominal amount of species A must be removed away from the penetrating C-liquid, diffusion-coupled-growth of these channels is significantly slowed, or even broken-down in some cases. And so, away from the grain boundary, interpenetrating A-rich solid and C-rich liquid phases tend to follow more tortuous pathways into the precursor. However, since the fast-diffusing grain boundary locally assists this diffusioncoupled-growth, a singular C-rich channel rapidly penetrates at the migrating grain boundary, forming a linear pathway into the precursor similar to the lamellar channels forming in precursors with lower A-content.
At the highest A-concentration considered, c A,0 = 50% in Figure 4(d), diffusion-coupled-growth of C-channels into the precursor is not possible away from the grain boundary. Indeed, in the single-crystal simulations in Supplemental Figure SM1(d), the phase field model predicts that, in the absence of the grain boundary, the precursor with c A,0 = 50% will not undergo bicontinuous dealloying. This high nominal A-content saturates the interface after some initial dealloying, thereby forming a passivating A-rich phase that prevents further dissolution of the precursor. These precursors are below the model's predicted parting limit, i.e. the requisite concentration of species B to fully dealloy the precursor.
However, as shown in Figure 4(d), a singular dealloying channel can still invade the precursor by the mechanism of GBMD. The result is a deep dealloying channel at the migrating grain boundary but a relatively passivated corrosion front elsewhere. To clarify, the model does not predict that the grain boundary will lower the parting limit, since this singular infiltration by the dealloying agent does not lead to the formation of a fully-dealloyed A-rich structure. Rather, these results show that for a system below the parting limit, the GBMD mechanism can still lead to the fast infiltration at grain boundaries.

DISCUSSION
The predicted morphologies resulting from GBMD can be compared to previous experimental studies of grain boundary effects in precursors undergoing liquid metal dealloying by McCue et al. 7,10 and Joo et al. 11 (albeit at time-and length-scales much earlier and smaller than the experiments). In these experiments, the dealloying within grains maintains a uniform depth, while at grain boundaries the dealloying is generally much faster/deeper (see Figure 8 of McCue et al. 10 and Figure 1 of Joo et al. 11 ). The thick, plate-like dealloyed structure that results near the migrating grain boundary in Figure 1(c) of the present work qualitatively matches the thick Tarich blocks shown in Figure 8 Figure 1(c), this A-rich plate only connects to the ligament network within one of the parent grains, and is separated from the other grain by the widening liquid melt -which continues to penetrate deeper into the sample along the migrating grain boundary.
The SEM image provided in Figure 1(c) of McCue et al. 7 nicely shows a 3D view of a dealloyed TiTa precursor near a grain boundary, which allows comparison with 3D multi-phase field predictions in Figure 2 of the present work. In both the experimental and the simulation results, dealloying corrosion forms three-dimensional ligament networks within the grains and two-dimensional ridges along grain boundaries. Further, both experiment and simulation show that the transition from threedimensional to two-dimensional corrosion morphology is abrupt.
In Joo et al. 11 , a "waviness" of several of the dealloyed grain boundaries was illustrated, wherein the connectivity of the flat dealloyed plates periodically switches back and forth between parent grains. We also note this waviness is seen along some of the grain boundaries in Figure 8 of McCue et al. 10 . In the present work, Figure  2 shows that the stochastic nature of GBMD can drive grain boundary migration in opposite directions simultaneously. We postulate that this could lead to a wavy dealloyed structure similar to experiments, since in the locations along the Y-axis in Figure 2 where the grain boundary migrates left (-Z direction), GBMD would impart a flat dealloyed plate along the right (+Z) grain, and vice-versa. We further postulate that a migrating grain boundary could periodically switch migration directions along the X-axis in Figure 2, perhaps due to an intermittency in the activation of the GBMD mechanism. Further, potentially larger, simulations would be helpful to explore this phenomenon in more detail.
Although the qualitative agreement between the present simulations and experimental observations described above is noteworthy, detailed quantitative comparisons with experimental results from McCue et al. 7,10 and Joo et al. 11 would require an understanding of how the numerical predictions herein will evolve of time-and length-scales more consistent with experiments. These dealloyed structures coarsen predominantly via surface diffusion processes 1,7,10,17,21 , and, away from the grain boundary, ligaments have been shown to maintain a similar morphology over different dealloying durations, with only the characteristic length scale of ligaments increasing over time 10,22,23 . Regardless of the mechanism, ligaments will coarsen to reduce the overall surface area of the structure via mass transport from surfaces with high curvature to surfaces with lower curvature. The large flat structure formed at the migrating grain boundary should, therefore, grow at the expense of the adjacent highlycurved ligaments. The coarsening rate depends on local curvatures 24 , and so coarsening near the grain boundary should be faster than at locations far away from the grain boundary, given that the features within grains are more similar in size and curvature. Therefore, we expect these numerically predicted corrosion morphologies herein to persist, or become even more exaggerated, with coarsening over longer time-scales.
The GBMD mechanism described in this work is distinct from the diffusion-induced-grain-boundary migration (DIGM) mechanism discussed by Balluffi and Cahn 25 . In the present work, an interdiffusivity is used which assumes fluxes among the solutes sum to zero and there is no explicit accounting for vacancy fluxes. Future work could explicitly include a vacancy species, and model point defect enrichment along a grain boundary that could drive migration through a DIGM process of grain boundary dislocation climb. In any case, the current work illustrates a mechanism of grain boundary migration that occurs absent consideration of vacancy flux.
In the present work, the fast penetration alongside the grain boundary by GBMD occurs in the absence of a grain boundary wetting condition. However, a wetting condition has recently been hypothesized to be the key feature of the heterogeneous dealloying seen experimentally at grain boundaries 11 . Specifically, it has been proposed that the grain boundary is quickly replaced by two solid-liquid interfaces, as would be the case if the appropriate thermodynamic condition is met that the excess energy of the solid-liquid interface is less than half that of the grain boundary. For the present model, the solid-liquid interfacial energy is influenced by the large enthalpy of mixing between species A and C. This mixing enthalpy leads to a spinodal driving force within the interface, which is an essential characteristic of these dealloying systems 1,5 . However, this large mixing enthalpy (Ω AC in Eq. 1) leads to a high excess interfacial grand potential, γ sl , for the dealloyed interface between C-rich liquid and A-rich solid. We note that this excess grand potential (integrated through the solid-liquid interface) is much greater than the solid-liquid interfacial free energy for the pure elements (σ sl ) in Eq. 1. The grain boundary, by comparison, has a much lower excess interfacial grand potential, since it resides in an AB alloy with ideal mixing thermodynamics such that compositional enrichment does not change its interfacial energy. This is evident in the contact angles formed at the grain boundary-liquid junction, where the liquid-phase dihedral angle approaches θ 3 = π.
Quantitative predictions made using the current model stand to be improved by a refinement of thermodynamic and kinetic parameters for the bulk phases and compositions, as well as for interfaces, junctions, and mixtures. For instance, since solute transport along the grain boundary is key to GBMD, predictions by the current work can be extended by atomistically informed models describing how ground boundary diffusion and mobility vary during dealloying. Further, the model highlights the importance of understanding how the mixing enthalpies vary through the interfacial region between the bulk phases, as this would impact the force balance at the triple junction. For example, a higher solubility between species A and C at the solid-liquid interface, which could be captured herein via a ramping-down Ω AC when ϕ 1 ϕ 3 > 0, would lower the predicted excess grand potential of this interface. Similarly, a driving force for solute segregation at the grain boundary would affect the relative interfacial energies and also the migration phenomenon by potentially adding a solute drag effect. The present model could be used in a future study to explore the relative effect of these driving forces through a parameter sensitivity analysis, which could be used to motivate future experiments and lower length scale simulation studies.
We expect that this mechanism of GBMD could also be active in other dealloying systems, e.g. molten salt dealloying (MSD). Indeed, MSD of NiCr precursors has recently been shown to be a viable method of producing three-dimensional bicontinuous metals 17 . Similar to LMD, in MSD the dealloying of the less noble element (Cr) coupled with interfacial diffusion processes leads to the formation of a porous metal made up of the more noble element (Ni). Further, MSD has been shown to be accelerated at, or in some cases confined to, grain boundaries in the NiCr precursor 14,20 . Since MSD is conducted at high homologous temperatures, with respect to the alloy precursor, we expect that grain boundaries are sufficiently mobile to activate this GBMD mechanism. Future work could explore this mechanism in MSD through the incorporation of electrochemical driving forces into the current model.

CONCLUSION
In this work, an advanced multi-phase field model is presented to study metal dealloying in the presence of grain boundaries. This new model expands upon the original work of Geslin et al. 1 to consider grain boundaries by incorporating the multi-phase field formulation by Steinbach and Pezzolla 18 . An advanced semi-implicit time integration algorithm by Badalassi et al. 19 , using spectral methods, is employed, which enables simulating the coupled evolution of an arbitrary number of phases, grains, and conserved species, over large length-and time-scales, while still maintaining a realistic width for interfaces -which is crucial for simulating the fine-scale interfacial curvatures seen in dealloying metals.
Using 2D and 3D simulations, the model is exercised in application to a liquid metal dealloying system (TiTa dealloying in liquid Cu), and we predict phenomena that is qualitatively consistent with the grain boundary effects seen in recent dealloying experiments 7,10,11 . Through a detailed analysis of the simulation results, we identify a mechanism of grain boundary migration-assisted dealloying (GBMD) that controls the morphology of these dealloyed metals. The key feature of this GBMD mechanism is that it assists the diffusion-coupled-growth of the infiltrating liquid dealloying agent, by providing an additional pathway in the metal precursor along which to reject the more noble element (and dealloy the less noble element). For the deliberate synthesis of nanoporous metals by dealloying, this mechanism is predicted to lead to large blocks of the non-dissolving metal species, which asymmetrically disrupt the connectivity of the final bicontinuous structure. For metals designed to withstand dealloying environments, e.g. with high-concentrations of the non-dissolving species, GBMD can still generate a deep dealloying channel that is expected to severely limit the safe operation of these material systems.

Multi-phase field model
The phase field model for LMD from Ref. 1 is extended to model an AB bicrystal precursor corroding in liquid C, using the multi-phase field formalism of Ref. 18 . The bicrystal precursor is made up of two solid phases (one for each grain) with binary composition c A,0 and c B,0 = 1 − c A,0 . The bicrystal is placed in contact with a liquid phase of pure c C,0 = 1. Within one grain, the non-conserved phase field variable ϕ 1 takes on the value ϕ 1 = 1, while in the other grain ϕ 2 = 1. In the liquid phase, the non-conserved phase field variable ϕ 3 = 1 − ϕ 1 − ϕ 2 = 1. The diffuse interface between grains where ϕ 1 ϕ 2 > 0 represents the grain boundary, and interfaces where ϕ 1 ϕ 3 > 0 or ϕ 2 ϕ 3 > 0 represent the solid-liquid interfaces for each grain. The diffuse region where ϕ 1 ϕ 2 ϕ 3 > 0 represents the junction among the two grains the liquid phase.
The phase evolution is given by: for independent phases α = 1, 2. In the above,M αβ and σ αβ are the interfacial mobility and excess free energy for the α-β interface, h ′ = 8 π ϕ α ϕ β is the derivative of an interpolation function between phases, ∆g αβ is the chemical free energy difference between phases, η is the interface thickness, andÑ is the number of phases locally present. The difference in energy between bulk solid and liquid phases is ∆g 13 = ∆g 23 = ∆g sl , and there is no difference in energy between each grain (∆g 12 = 0).
The concentration of species i evolves by the continuity equation: with phase-dependent solute mobility matrix M ij , and chemical potential µ j of species j. In Eq. 2, the jsummation is over independent compositions, c A and c B . The chemical potential is defined by: where f chem is the bulk chemical free energy, and κ is the energy penalty coefficient for composition gradients. Following 1 , the chemical free energy is defined using a regular solution model: The first term is the entropy of mixing, with Boltzmann constant k B , temperature T , and atomic volume V a . The second term is the excess enthalpy associated with mixing elements A and C with factor Ω AC . Following 1 , Ω AC ≫ 0, which leads to a large miscibility gap between A and C, while all other element pairs mix ideally. The third term contains the difference in chemical energy between the liquid and solid phases, ∆g sl ({c i }), weighted by a smooth interpolation function h(ϕ 3 ) where h = 1 in the solid phases (ϕ 1 + ϕ 2 = 1), and h = 0 in the liquid phase (ϕ 3 = 1). The free energy difference between phases is given as: where T i and L i are the melting temperature and latent heat of melting, respectively, for species i. The interpolation function is: The solute mobility is defined as: with a diffusivity that varies among solid, liquid, and grain boundary regions: where D l , and D gb are the diffusivities in the liquid phase and grain boundary region, respectively. Following 1 , diffusivity in the solid phase is assumed to be negligibly small relative to the liquid phase and is set to zero. The form of the solute mobility matrix ensures that, in Eq. 2, Fickian diffusion of Ti is recovered in the liquid Cu. At the solid-liquid interface for each grain, Eqs. 1 and 2 reduce to similar forms as used in the dualphase field model by Geslin et al. 1 . The only major difference lies in the choice of interpolation function, h, between phases. The interpolation function in Eq. 6, following 28 , and the approximation h ′ = 8 π ϕ α ϕ β (where α=1,2 and β=3) in Eq. 1 are selected so that chemical driving forces on ϕ α reach zero when ϕ α = 0. In the present simulations, this is important for preventing the non-physical wetting of one grain along the other grain's liquid interface. This triple-junction delocalization could have also been addressed by increasing the energy of the triple-junction, as discussed by Nestler et al. 29 and Hötzer et al. 30 . During development of the present model, different thermodynamic descriptions for the grain boundary-liquid junction were explored, including using the "anti-symmetric" approximation for the multi-phase field equations 26 and triple-point energy manipulation. In all cases, the grain boundary migrationassisted dealloying mechanism remained active, suggesting that the choice of triple-junction treatment does not qualitatively affect the results presented herein.
The grain boundary mobility isM 12 = 1.2 m/(s·GPa), and grain boundary diffusivity is 7 × 10 −10 m 2 /s. The prescribed interface thickness is η = 4 nm. The grain boundary current (i.e. the 1D integral of the diffusivity across the 1-2 interface) is therefore D gb · 3η/8 = 12.0 m/(s · GPa) D gb 7 × 10 −10 m 2 /s D l 7 × 10 −9 m 2 /s 10 −18 m 3 /s. These kinetic parameters are within realistic ranges for varying grain boundary types [31][32][33][34] . The composition-gradient energy-penalty is κ = 15.0 eV/nm. All the simulations presented in the main text use a grain boundary energy σ 12 = 0.39 J/m 2 , which is chosen to be higher than the solid-liquid interfacial energy (σ 13 = σ 23 = 0.2 J/m 2 ) from 1 , but not so high as to enforce a grain boundary wetting condition. This grain boundary energy is reasonable for the elevated temperature at which liquid metal dealloying is conducted 35,36 . All other material properties are taken from directly from Geslin et al. 1 . All thermodynamic and kinetic parameters for the multi-phase field model are given in Table I. For clarity, in Table I the interfacial parameters for interactions between phases 1 and 2 are denoted "gb" since they correspond to the grain boundary, and interfacial parameters between phases 1 and 3 and between phases 2 and 3 are denoted "sl" for the solid-liquid interface.

Numerical implementation
A representative simulation cell is shown in Figure 5. The simulations are initialized with varying initial alloy compositions c A,0 = 0.2, 0.3, 0.4, 0.45, and 0.5 and, following 1 , each alloy's composition is perturbed around its initial value with a uniform white noise of amplitude ±0.025. Dirichlet boundary conditions are enforced on the top of and bottom edges of each simulation. On the left and right edges, periodic boundary conditions are enforced for the conserved species and the non-conserved liquid phase. Periodic boundary conditions on the right and left edges for the non-conserved grains would result a sharp transition between ϕ 1 and ϕ 2 , which would quickly introduce a second grain boundary. To avoid this, each grain's phase takes on the value of its periodic neighbor at the right and left edges on the simulation cell. That is, ϕ 1 assumes the value of ϕ 2 , and vice-versa, past the left and right edges when solving Eq. 1. These "quasi-periodic" boundary conditions in the horizontal direction are pre- ferred to Neumann (mirror-plane) boundary conditions, so that corrosion channels can freely "wander" across the left and right edges. (In the single crystal simulations, this not an issue, and regular periodic boundary conditions are enforced for all conserved and non-conserved variables). To avoid boundary effects along the bottom edge, "continuous -gradient" boundary conditions are enforced on phases ϕ 1 , ϕ 2 , and ϕ 3 following 27 . This allows the grain boundary to intersect the bottom boundary at non-right angle without penalty.
The phase fields are updated using a forward Euler time integration. However, using forward Euler to update the composition fields would require prohibitively small time-steps to maintain numerical stability. Therefore, the composition fields are updated using the semiimplicit time-stepping procedure following 19 , which uses spectral methods. In the present work, Eq. 2 is computed in real-space, and then c i and ∂c i /∂t are transformed into frequency-space to semi-implicitly update the composition fields. Discrete sine-transforms are used in the vertical directions to maintain Dirichlet boundary conditions, while discrete real-to-complex transforms are used in the horizontal directions for the periodic boundary conditions. The FFTW library is used to compute transforms 37 , and finite difference operators for derivatives in frequency-space are used to avoid spurious Gibbs oscillations in the composition fields 38 . For more details on the numerical implementation in frequency-space, see Supplementary Materials SM1. In the present work, the semi-implicit time-stepping scheme allows taking timesteps of size ∆t = 10 −3 ns, which is approximately 100 times larger when compared to a forward Euler scheme, without loss of stability or accuracy. (Even larger time steps could be used, and still maintain numerical stability, but were found to noticeably change the simulation results). Since numerical integration of Eq. 1 is prone to violating the physical constraint that phases remain on the Gibbs simplex, i.e. 0 ≤ ϕ α ≤ 1 39 , unphysical phases {ϕ α } are projected back onto the Gibbs simplex following 40 . The simulations are discretized with gridpoint spacing ∆x = 0.2 nm, and Eqs. 1 and 2 are computed using the isotropic finite difference scheme by Ji et al. 41 .
Post-processed figures of dealloying precursors were generated using colormaps from Crameri 42 .

DATA AVAILABILITY
The data produced and analysed in this study are available from the corresponding author on reasonable request.

Grain Boundary Effects in Dealloying Metals
where (c i ) n+1 m is the composition of species i at gridpoint m = ⟨m x , m y , m z ⟩, which is being solved for at next time-step, n + 1, based on the composition at the current and previous time-steps, respectively (c i ) n m and (c i ) n−1 m , and their time-evolution, (ċ i ) n m and (ċ i ) n−1 m , which are computed using Eq. 2 of the main text. Along each direction, the gridpoints are indexed m j = 1, 2, ..., N j , where N j is the grid-length in the j-direction. In Eq. SM1, ∆t is the time-step size, and λ and τ are semi-implicit parameters as detailed in [?].
Equation SM1 is solved using spectral methods. We define the following discrete Fourier transform: (ĉ i ) n k = DFT (c i ) n m k (SM2) and associated inverse discrete Fourier transform as: where (ĉ i ) n k is the spectral-space representation of species i at frequency k = ⟨k x , k y , k z ⟩, at a given time-step n. Along each direction, the discrete frequencies are indexed k j = 1, 2, ..., N j . The DFT operator in Eq. SM2 is a multi-dimensional discrete Fourier transform, which is chosen based on the boundary conditions in each direction. In the current work, Dirichlet boundary conditions are enforced along the X-direction, and periodic boundary conditions are enforced for the composition evolution along the Y-and Z-directions. This requires using the sine transform (DST-II) along the X-axis and using real-to-complex transforms along the Y-and Z-axes, and so Eq. SM2 is fully defined as: The associated inverse transform uses the inverse sine transform (DST-III), and complex-to-