Emergence of rapid solidification microstructure in additive manufacturing of a Magnesium alloy
D. Tourret R. Tavakoli A. D. Boccardo A. K. Boukellal M. Li J. Molina-Aldareguia



Rapid solidification microstructure in Mg alloys]
Emergence of rapid solidification microstructure in additive manufacturing of a Magnesium alloy



^1 IMDEA Materials, Getafe, Spain
^2 Sharif University of Technology, Tehran, Iran
^3 Université de Lorraine, IJL/CNRS, Nancy, France
^4 Department of Mechanical Engineering, Universidad Politécnica de Madrid, 28006 Madrid, Spain
damien.tourret@imdea.org



Bioresorbable Mg-based alloys with low density, low elastic modulus, and excellent biocompatibility are outstanding candidates for temporary orthopedic implants.
Coincidentally, metal additive manufacturing (AM) is disrupting the biomedical sector by providing fast access to patient-customized implants.
Due to the high cooling rates associated with fusion-based AM techniques, they are often described as rapid solidification processes.
However, conclusive observations or rapid solidification in metal AM — attested by drastic microstructural changes induced by solute trapping, kinetic undercooling, or morphological transitions of the solid-liquid interface — are scarce.
Here we study the formation of banded microstructures during laser powder-bed fusion (LPBF) of a biomedical-grade Magnesium-rare earth alloy, combining advanced characterization and state-of-the-art thermal and phase-field modeling.
Our experiments unambiguously identify microstructures as the result of an oscillatory banding instability known from other rapid solidification processes.
Our simulations confirm that LPBF-relevant solidification conditions strongly promote the development of banded microstructures in a Mg-Nd alloy.
Simulations also allow us to peer into the sub-micrometer nanosecond-scale details of the solid-liquid interface evolution giving rise to the distinctive banded patterns.
Since rapidly solidified Mg alloys may exhibit significantly different mechanical and corrosion response compared to their cast counterparts, the ability to predict the emergence of rapid solidification microstructures (and to correlate them with local solidification conditions) may open new pathways for the design of bioresorbable orthopedic implants, not only fitted geometrically to each patient, but also optimized with locally-tuned mechanical and corrosion properties.





Keywords: 
Rapid solidification; 
Additive Manufacturing; 
Magnesium alloy;
Microstructure.




§ INTRODUCTION



Metal additive manufacturing (AM) is already making a transformative impact in the biomedical sector, by enabling on-demand cost- and time-effective fabrication of implants customized to the needs of each patient.
In particular, biodegradable and bioabsorbable metals – like magnesium (Mg), iron (Fe), and zinc (Zn) – can degrade gradually in vivo and be either absorbed or eliminated by the body as the surrounding tissue heals, hence circumventing the need for an extra removal surgery <cit.>.
Among these materials, Mg alloys, often combined with rare earth (RE) elements, constitute a prime candidate for orthopedic implants <cit.>, as they offer good biocompatibility, osteopromotive properties <cit.>, and mechanical properties close to that of human bones <cit.>, thus reducing stress shielding during loading of the bone-implant interface <cit.>. 
The cytotoxicity of rare earth (RE) elements has been studied both in vitro and in vivo <cit.>. Some of them, such as lanthanum (La), cerium (Ce) and praseodymium (Pr) exhibit high cytotoxicity — in particular, Ce and Pr chlorides present severe hepatotoxicity <cit.>. In terms of 50% lethal dose (LD_50) <cit.>, Ce shows the highest cytotoxicity (10 mg/kg), followed by La (150 mg/kg), but several other RE elements (> 550mg/kg for Gd, Nd, Dy, and Y) are viable for biomedical applications. For instance, Mg–1.8Zn–0.2Gd alloy was found to exhibit little to no toxicity on several different cell lines <cit.>, and Nd was found to improve mechanical properties and corrosion resistance within acceptable toxicity levels <cit.>.
Moreover, high cooling rates were reported to be efficient in modifying microstructures in a wide range of Mg alloys, resulting in improved mechanical properties <cit.> and corrosion resistance <cit.>.
Therefore, fusion-based AM processes, e.g. laser powder-bed fusion (LPBF), which usually yield high cooling rates, provide a promising route to tune processing parameters and achieve desired microstructures, properties, and performance of 3D-printed Mg-based bioimplants.
Yet, the fundamental understanding of microstructure formation in far-from-equilibrium conditions remains incomplete.

In solidification theory <cit.>, rapid solidification refers to a regime in which the solid-liquid interface exhibits a strong departure from equilibrium when reaching a sufficiently high velocity V. 
One manifestation of rapid solidification is solute trapping, whereby the solid phase grows with a solute concentration in excess of its equilibrium solubility limit.
Thus, the solute partition coefficient k=c_s/c_l, with c_s and c_l the solute concentration on the solid and liquid sides of the interface (respectively), deviates from its equilibrium value and tends toward unity as V approaches and surpasses a characteristic diffusion velocity V_ d, typically of the order of 1 m/s.
Meanwhile, both liquidus and solidus lines of the phase diagram converge toward the so-called T_0-line as V increases.
The velocity-dependence of the interface solute partition coefficient, k, and of the liquidus slope, m, are commonly described using the continuous growth model (CGM) <cit.>

    k(V)     = k_ e+V/V_ d/1+V/V_ d
    m(V)/m_ e   = 1-k(V)+[k(V)+(1-k(V))α]ln[k(V)/k_ e]/1-k_ e  ,

where k_ e and m_ e stand for the respective equilibrium (V≪ V_ d) values of k and m, and α is a solute drag coefficient.
Another rapid solidification effect is the occurrence of kinetic undercooling, associated with a deviation of the interface temperature, T, proportional to the interface velocity. 
The resulting velocity-dependent liquidus and solidus temperatures (for a planar interface, i.e. without curvature undercooling) are thus given by 

    T_ L(V)    = T_ M - m(V)c_∞ - V/μ
    
    T_ S(V)    = T_ M - m(V)/k(V)c_∞ - V/μ

with T_ M the melting temperature of the solvent species, c_∞ the alloy nominal solute concentration, and μ the interface kinetic coefficient.
While solute trapping can be probed by a compositional analysis of rapidly solidified samples, kinetic undercooling is much harder to measure experimentally, since it requires correlating the interface velocity with its temperature during solidification under deep undercooling conditions <cit.>.
Finally, the most evident observable indication of rapid solidification is the transition of the microstructure pattern from dendrites to a planar interface, at velocities beyond the absolute stability threshold V_ a.
An extension of the linear perturbation theory by Mullins and Sekerka <cit.> to high undercooling <cit.> leads to the implicit definition of the threshold velocity above which the planar interface is stable

    V_ a=D_ LΔ T_0(V_ a)/k(V_ a)Γ=D_ Lm(V_ a)c_∞[1-k(V_ a)]/k(V_ a)^2Γ

with D_ L the solute diffusion coefficient in the liquid, Δ T_0 the alloy freezing range between the liquidus (T_ L) and solidus (T_ S) temperatures, and Γ the Gibbs-Thomson coefficient of the solid-liquid interface. 
Equation eq:absstab only provides an order of magnitude of the absolute stability limit, which is better estimated by the maximum of T_ S(V) from eq:TsV <cit.>, while the region with T/ V>0 may lead to an oscillatory instability at the origin of banded microstructures <cit.>.

Indeed, intermediate velocities between the dendritic (V≪ V_ d) and planar partitionless (V≫ V_ d) regimes may lead to an oscillatory instability of the solid-liquid interface. 
Resulting microstructures consist of an alternation of (i) featureless regions of homogeneous supersaturated solute content with (ii) patterned regions exhibiting solute segregation and often decorated with small particles of secondary phases such as intermetallic precipitates.
The resulting bands are nearly parallel to the solidification front (i.e. normal to the main temperature gradient and growth direction).
These banded microstructures have been identified in rapidly solidified samples using a broad range of processing techniques (e.g. splat quenching, melt spinning, laser melting, etc.), notably on several of Al-based alloys (e.g. Al-Cu <cit.>, Al-Pd <cit.>, Al-Zr <cit.>, Al-Fe <cit.>,  Al-Er <cit.>, etc.) but also other alloys, like Ag-Cu <cit.>.
These observations have shown that the width of the patterned portion of the bands decreases with an increase in growth rate <cit.>.
Recently, the formation of these bands in laser-melted thin Al-Cu samples was also observed in situ within a transmission electron microscope (TEM), allowing the measurements and correlation of the morphological transition with an increase of the interface velocity <cit.>.

Analytical and numerical studies have identified the origin of the oscillatory instability as a subtle combination of interface undercooling and solute trapping <cit.>. 
Essentially, within a velocity range just below absolute planar stability, the effect of solute trapping is dominant compared to the effect of atomic attachment kinetics, resulting in an increase in solidus temperature (i.e. the stable temperature for a planar interface) with interface velocity (dT/dV>0). Within a temperature gradient moving at a given velocity V, this range of velocity is understood to be unstable for the growth of a planar interface. Hence, therein, the interface oscillates between two stable branches of solutions with dT/dV<0, namely corresponding to a dendritic pattern at relatively lower V and a planar interface pattern at relatively higher V. This oscillation cycle, first theorized considering instantaneous transition between the two stable branches <cit.> was later confirmed by further analytical and numerical studies <cit.>.
These studies have highlighted the key role of the latent heat diffusion in the selection of the band spacing — rendering it essentially unaffected by the temperature gradient <cit.>, while bands are inversely proportional to the temperature gradient if latent heat is neglected <cit.>.
Recent phase-field simulations have reproduced quantitatively the banding mechanism observed in situ in Al-Cu TEM samples <cit.> and shown that the onset of banding instability was closely approximated by eq:absstab, in spite of the underlying linear stability analysis not incorporating anisotropies of excess free energy and kinetic coefficient essential to the interface pattern selection <cit.>.

In this context, although LPBF is commonly referred to as a “rapid solidification” process, whether or not the solid-liquid interface strongly departs from equilibrium is often ambiguous.
Most additively manufactured microstructures indeed exhibit patterns (dendrites) as well as significant intergranular and interdendritic solute segregation <cit.>, indicative of a relatively weak departure from equilibrium at the solid-liquid interface.
Rather, “banded” microstructures reported in metal AM often stem from a periodic change in local processing conditions, for instance due to gas jets disturbing the melt pool and periodically changing the local cooling rate <cit.> or post-solidification microstructure evolution in the heat affected zone <cit.>.

Interestingly, nearly forty years ago, Mg alloys processed by rapid solidification processes, such as splat quenching and melt spinning, were reported to contain “featureless zones” indicative of partitionless planar solidification beyond absolute stability <cit.>. 
These featureless regions exhibited a greater hardness than both the as-received (primary phase) alloys and the dendritic regions of the rapidly solidified material. 
More recently <cit.>, a bioresorbable Mg alloy (WE43) additively manufactured via LPBF was reported to exhibit intragranular banded microstructures reminiscent of those observed in earlier studies of rapid solidification microstructures in Al alloys.

The underlying objective of the present article is to investigate, using state-of-the-art quantitative modeling tools, whether solidification conditions relevant to LPBF of Mg alloy WE43 are indeed expected to result in the oscillatory banding instability linked to a strong equilibrium departure at the solid-liquid interface. 
To do so, we apply a recently proposed quantitative phase-field model of rapid solidification <cit.>, using local solidification conditions estimated by a thermal model of laser melting and solidification of a powder bed. 
We explore microstructural patterns expected in different regions of the melt pool, as well as their complex formation dynamics.




§ METHODS






 §.§ Experiments


We support our analyses with experimental data already presented in detail elsewhere <cit.>, including electron microscopy, X-ray tomography, mechanical and corrosion testing, as well as in vitro cytocompatibility. 
Hence, here we only summarize the key information relevant to the material and manufacturing conditions leading to the emergence of rapid solidification microstructures and the methods used in their characterization.
The reader is referred to earlier works for further relevant details on additive manufacturing processing <cit.> and sample characterization and testing <cit.>.



  §.§.§ Material and manufacturing.
 

The material feedstock is a powder of Mg alloy WE43 MEO (Meotec GmbH, Aachen, Germany) with a nominal composition of 1.4–4.2 wt% Y, 2.5–3.5 wt% Nd, <1wt% (Al, Fe, Cu, Ni, Mn, Zn, Zr), and balance Mg  with a particle distribution D_10 = 29.1µm, D_50 = 45.8µm and D_90 = 64.4µm <cit.>, naturally covered with a thin Y_2O_3 passivation layer limiting its flammability. 
Lattice structures (open scaffolds) of size (10 mm)^3, with a body-centered cubic unit cell structure and different strut diameters (250 to 750 µm), were additively manufactured by LPBF in an Ar atmosphere with a laser power P≈90W, a beam diameter =90µm, a layer thickness H_ p≈30µm, and a 0^∘/90^∘ scan strategy across successive layers.
In order to prepare TEM samples, an additional bulk cylindrical sample (8 mm in diameter, 30 mm in length) was manufactured with similar processing parameters as the lattices.



  §.§.§ Microstructural characterization.
 

For scanning electron microscope (SEM) observation, samples were cut out of the scaffolds, then manually ground with SiC abrasive paper 4000 grit, followed by mechanical polishing with 3 µm and 0.25 µm diamond paste using an alcohol-based lubricant to avoid oxidation, and finally etching with a mixture of 75 ml ethylene glycol, 24 ml distilled water and 1 ml nitric acid.
For TEM analysis, thin samples of 200 µm thickness were cut from the cylindrical specimen, manually polished down to about 80 µm, and then polished with an etchant of 10.6 g lithium chloride, 22.32 g magnesium perchlorate, 200 ml 2-butoxy-ethanol and 1000 ml methanol at -40^∘C and 95 V. 
TEM observation was carried out in a FEI Talos equipped with a field emission gun operating at 200 kV. 
Chemical composition mapping via energy-dispersive X-ray spectroscopy (EDS) was performed under scanning-transmission electron microscopy (STEM) mode with a SuperX detector, with an acquisition time between 30 and 60 min, and the resulting data was analyzed using Bruker QUANTAX software.
The micrographs appearing in section <ref> (figure <ref>) were further processed (brightness and contrast adjustment) using the software ImageJ <cit.>.




 §.§ Thermal modeling





  §.§.§ Finite element model.
 
We assessed the temperature profile using the finite element (FE) software Abaqus <cit.>.
The main objective of these calculations is to estimate the typical melt pool dimensions, and the temperature gradient G and growth velocity V experienced by the solidification front.
We do not explore the effect of geometry of the printed lattice structures, but rather seek an order of magnitude of the  (G,V) conditions to use in the phase-field simulations of microstructure formation, and hence simply simulate a single linear melting/solidification track.

The heat transfer equation, prescribing the evolution of the temperature field T with time t, follows <cit.>

    ρ c_p∂ T/∂ t= ∇·( κ ∇ T )+Q_ L+Q_ϕ  ,

where ρ is the density, c_p is the specific heat capacity, κ is the (isotropic) conductivity, Q_ L is the absorbed heat due to
the laser, and Q_ϕ is the absorbed/released heat due to phase transformations.

The volumetric heat source from the moving laser is introduced via the  subroutine <cit.> as 

    Q_ L=α_ pG_ LI_z/H_ L

where α_ p is the absorptivity of the powder bed and H_ L is the laser penetration depth. 
The Gaussian distribution along the top (z=z_ L) surface, G_ L, and the parabolic decay along the laser penetration direction z, I_z, are given by  

    G_ L   = 2 P/π^2exp[ -2(x-x_ L)^2+(y-y_ L)^2/^2] 
    
      I_z    = 1/0.75[ -2.25 (z_ L-z/H_ L)^2+1.5  (z_ L-z/H_ L) +0.75 ]   ,

where P is the laser power,  is the laser diameter, and (x_ L, y_ L, z_ L) is the location of the laser heat source center 
within a cartesian coordinate system (x, y, z), with the laser moving at a velocity V_ L in the x+ direction.

The source term related to phase transformation during melting and solidification of the alloy is implemented via the  subroutine and calculated as

    Q_ϕ=ρ L_ f∂ f_ s/∂ t  ,

where L_ f is the latent heat of fusion and f_ s the volume fraction of solid, approximated by a linear interpolation between the alloy liquidus (T_ L) and solidus (T_ S) temperatures.

We apply boundary conditions (BCs) following Newton's law

    q_ c=-h_ c(T-T_∞)   ,

where q_ c is the normal heat flux, h_ c is the interfacial heat transfer coefficient, and T_∞ is the surrounding temperature. 
This condition is applied on all domain boundaries, but the domain is taken large enough laterally to result in a negligible effect of BCs in x and y directions on the melt pool dimensions, while different heat transfer coefficients are applied on the top (z+) alloy/air interface and the bottom (z-) alloy/build plate interface.
Different properties are considered for the powder bed and for the dense (solid or liquid) material. 
The initial powder-bed state transitions to the liquid state when the local temperature exceeds the alloy liquidus temperature, T_ L, and is then considered solid when the temperature drops below the solidus temperature, T_ S.



  §.§.§ Parameters.
 
Material properties and process parameters are listed in Tables <ref> and <ref>, respectively.
Most thermal properties (Table <ref>) are approximated by those of pure Mg using values from classical handbooks <cit.>, including temperature-dependent conductivity κ(T) and specific heat capacity c_p(T) <cit.> — the latter being converted from molar to mass using the molar mass M_ Mg=24.305 g/mol.
For consistency with simulations of microstructure formation, liquidus and solidus temperature are calculated for a Mg-3wt%Nd alloy (see section <ref>).
The powder bed is assumed to have a slightly reduced density (by a factor 0.7, just below that of a close packing of identical spheres, i.e. with a porosity fraction ξ=0.3) but a significantly reduced conductivity.
For the latter, we follow the model proposed by Sih and Barlow <cit.>, derived from the Zehner-Schlünder-Damköhler equation <cit.> assuming spherical particles 


    κ_ p/κ_ g =    (1-√(1-ξ)) (1+ξκ_ r/κ_ g) 
            + √(1-ξ){2/1-κ_ g/κ[1/1-κ_ g/κln(κ/κ_ g) - 1] + κ_ r/κ_ g}

where κ, κ_ p, and κ_ g are the conductivity of the bulk (dense) material, powder bed, and gas, respectively.
The thermal conductivity due to radiation among particles, κ_ r, is


    κ_ r = 4Fdσ T^3

where F is a view factor, d is the average diameter of powder particles, and σ is the Stefan-Boltzmann constant (5.67×10^-8 W/m^2/K^4).
Considering the temperature-dependent thermal conductivity of Mg from <cit.> (Table <ref>) and that of Ar gas measured in <cit.>, a view factor F=1/3 <cit.>, and a particle size d≈ 50µm, we found that the conductivity of the powder-bed was well approximated by a second order polynomial of the temperature and used directly this approximation, listed in Table <ref>, in the simulations.
The resulting conductivity of the powder bed is over two orders of magnitude lower than that of the bulk material, with κ_ p/κ≈0.0058 at 923 K.
Process parameters (Table <ref>) such as laser parameters (P, V_ L, ) and the powder layer thickness (H_ p) are taken according to experimental conditions.
The laser penetration depth H_ L≈1.5H_ p and top/bottom heat transfer coefficients are reasonable orders of magnitude typically used in metal processing simulations <cit.>.











  §.§.§ Simulations.
 

Since we only seek reasonable estimates of the (G,V) conditions in the region surrounding the melt pool during LPBF of our Mg alloy, here we chose to disregard the particular geometry of the lattices and the printing strategy, and rather calculate the temperature distribution along a straight single track at steady state using relevant material properties (Table <ref>) and processing parameters (Table <ref>).
The simulated sample consists of a wide plate of solid Mg alloy (4.5 mm long, 3.75 mm wide, and 470 µm high) with a H_ p=30 µm layer of Mg alloy powder at the top. 
The domain is divided in two regions: (i) the region directly affected by the laser (105 µm deep and 750 µm wide, laterally centered on the laser track path) finely meshed with 1 902 600 linear hexahedral elements (DC3D8) of edge size 2 µm, and (ii) the overall surrounding region meshed with 224 977 linear tetrahedral elements (DC3D4). 
We defined the mesh size by performing a convergence analysis of the temperature field in the melt pool.
The simulated time is 0.9 ms, which is sufficient to reach a stable steady melt pool shape and temperature profile therein.
The simulation completion took approximately 9 hours on 28 CPU cores (Intel Xeon Gold 6130 2.10GHz).



  §.§.§ Post-processing.
 

From the FE-calculated thermal field, we extract local (G,V) solidification conditions along the tail of the melt pool, to be used in the one-dimensional approximation of the temperature field in the PF simulations. 
Since the melt pool is nearly elliptical in shape (see section <ref>, figure <ref>), we approximate the location of T_ L and T_ S isotherms by ellipses <cit.> with  

    r_ L(θ) = √((l_ L d_ L)^2/(d_ Lcos
    (θ))^2 + ( l_ Lsin (θ))^2)

and

    r_ S(θ)  = √((l_ S d_ S)^2/(d_ Scos
    (θ))^2 + ( l_ Ssin (θ))^2)  ,

where r_ L and r_ S are the respective radii of the (T=T_ L) and (T=T_ S) ellipses as a function of the angle measured counterclockwise from the top surface <cit.>

    θ = tan^-1|z-z_0/x-x_0|

with (x_0,z_0) the assumed common center of both ellipses (close but not exactly overlapped with (x_ L,z_ L)).
The dimensions of the melt pool tail appear explicitly in eq:rl and eq:rs as the length (l_ L, l_ S) and depth (d_ L, d_ S) of the solidus (subscript S) and liquidus (L) isotherms.
A linear radial interpolation of the temperature between T_ L and T_ S provides a reasonable approximation of the temperature field in the region where solidification takes place <cit.>.
Hence, since the polar component of the temperature gradient, (1/r)∂ T/∂θ, is usually much lower than its radial component, this approximation also provides a simple and convenient approximation of the amplitude of the temperature gradient 

    G(θ) ≈∂ T/∂ r≈T_ L-T_ S/r_ S(θ)-r_ L(θ)  ,

whereas the local velocity of the isotherms, i.e. the local solidification velocity, is simply 

    V(θ) = V_ Lcos(θ)
      .





 §.§ Microstructure formation modeling





  §.§.§ Phase-field model.
 
To simulate microstructure formation, we use a phase-field (PF) model for rapid solidification recently proposed by Ji et al. <cit.>. 
We adapt the formulation of the anisotropic terms for a hexagonal (hcp) Mg alloy, and apply it in two dimensions (2D), within the basal plane, considering local values of temperature gradient G and growth velocity V estimated via the thermal calculations (Section <ref>).

It is known that using a diffuse interface width W much greater than the actual physical width of the solid-liquid interface W_0 leads to artificial solute trapping effect, even at moderate growth velocities, which is typically corrected by adding an “anti-trapping” term into the solute conservation equation <cit.>.
Meanwhile, when using a physically realistic interface width W=W_0∼ 1nm, PF models may be capable of quantitatively capturing solute trapping <cit.>.
In the model used here <cit.>, the spurious trapping effect due to the use of W/W_0>1 is compensated via an increased solute diffusivity through the interface.
For a given value of S=W/W_0, the interpolation function for the diffusivity across the interface can be chosen to ensure that results for S>1 remain quantitatively close to those at S=1 in terms of partition coefficient k(V) and liquidus slope m(V) across a broad range of interface velocity V, where k varies from its equilibrium value k_ e at moderate V up to k→ 1 at high V.
By enabling upscaled yet quantitative simulations with S>1, the resulting PF simulations allow studying the morphological details of solid-liquid interface pattern evolution at experimentally relevant length and time scales. For instance, the model was recently used to quantitatively reproduce the banding growth instability observed in thin Al-Cu samples melted and observed inside a transmission electron microscope <cit.>.

The evolution of the phase field, ϕ, and solute concentration field, c, follows <cit.>

    τ(𝐧)∂ϕ/∂ t   =    ∇⃗·[ W(𝐧)^2∇⃗ϕ] +ϕ-ϕ^3
    +∑_η=x,y[ ∂_η( |∇⃗ϕ|^2 W(𝐧) ∂ W(𝐧)/∂ (∂_ηϕ)) ] 
    
           -λ g'(ϕ)[ c+T-T_ M/m_ eexp{b(1+g(ϕ))}]  ,
    
        ∂ c/∂ t   =    ∇⃗·{ D_ L q(ϕ)c∇⃗[ln c - bg(ϕ)] } ,

where we use ∂_η to denote ∂/∂η, b=ln(k_ e)/2, g(ϕ)=15(ϕ-2ϕ^3/3+ϕ^5/5)/8, and λ=(2√(2)/3)bm_ eSW_0/[Γ(k_ e-1)], with the Gibbs-Thomson coefficient Γ=γ_0T_ M/(ρ L_ f), where γ_0 is the solid-liquid interface (average) excess free energy, T_ M is the solvent (Mg) melting point, and ρ L_ f is the latent heat of fusion per unit volume. 
The (x,y) cartesian coordinate system used here differs from the (x,y,z) coordinates used in the FE thermal simulations, as the x+ direction of the PF simulation corresponds to the main temperature gradient direction along the fusion line, directed toward the center of the melt pool.
Given a solute diffusivity D=0 in the solid and D=D_ L in the liquid phase, we consider an enhanced diffusivity through the interface, D(ϕ)=D_ Lq(ϕ), using the same function q(ϕ)=A(1-ϕ)/2-(A-1)(1-ϕ)^2/4 as in <cit.>, in order to compensate for the spurious solute trapping effect due to S=W/W_0>1.
For a given value of S, an appropriate value of the coefficient A is estimated similarly as in <cit.>, so as to minimize the departure of k(V) and m(V) from curves calculated at S=1 across the entire velocity range (see section <ref>).
We consider the usual frozen temperature approximation, whereby a one-dimensional temperature field is imposed as T=T_0+G(x-x_0-Vt) with V the velocity of the isotherms (i.e. the growth velocity in the x direction at steady state) with respect to a fixed reference temperature T_0 at a point x_0 at a time t=0.
For the hexagonal close-packed Mg alloy, the anisotropic solid-liquid interface excess free energy, γ( n), and kinetic coefficient, μ( n), here considered within the basal plane, are expressed using spherical harmonics <cit.> as

    a_ s( n)   =γ(n)/γ_0 =  1 - ε^γ_20/4√(5/π) + ε^γ_66/64√(6006/π)cos(6θ)   , 
    
    a_ k( n)   =μ(n)/μ_0 =  1 - ε^μ_20/4√(5/π) + ε^μ_66/64√(6006/π)cos(6θ)   ,
 
where γ_0 and μ_0 are the average values of γ and μ.
In eq:aniso_g-eq:aniso_k, we use the same form for γ(n) and μ(n) and consider that coefficients ε^γ_40=ε^γ_60=ε^μ_40=ε^μ_60=0, which was established by molecular dynamics (MD) calculations for γ( n) <cit.> and is here assumed to be also valid for μ( n).
The resulting anisotropic interface width and relaxation time are W(𝐧) = SW_0a_ s( n) and τ(𝐧) = τ_0 a_ s( n)^2/a_ k( n), where τ_0 = (SW_0)^2/(Γμ_0) and the orientation θ of the interface normal n=(n_x,n_y) is θ=n_y/n_x=tan^-1(∂_y ϕ / ∂_x ϕ ).



  §.§.§ Computational implementation.
 
Equations eq:dpdt-eq:dcdt are solved using an explicit Euler time-stepping scheme with a finite difference spatial discretization, implemented in C-based CUDA (Compute Unified Device Architecture) programming language for massive parallelization on multiple Graphics Processing Units (GPUs).
To minimize artificial grid anisotropy effects, Laplacian and divergence terms in eq:dpdt-eq:dcdt are discretized using operators that are isotropic at order h^2 of the grid spacing h <cit.>.
Anisotropic terms are solved similarly as in <cit.> (see Appendices therein), i.e. by expanding the terms in eq:dpdt containing a_ s( n) or its spatial derivatives, and expressing them as a function of first and second order spatial derivatives of the phase field ϕ.
The expressions used for the sixfold anisotropy considered here are fully developed in the Appendix at the end of this article.
The resulting anisotropy implementation was verified to yield similar results as a different thoroughly validated implementation <cit.> in the low-to-moderate velocity regime (using the model from <cit.>).

The multi-GPU implementation is limited to GPUs hosted on a single node, using a simple layer-wise domain decomposition.
Therein, similarly as in <cit.>, for a parallelization on n GPUs, the domain of total dimension L_x× L_y is divided in n layers of dimension L_x×(L_y/n+2), considering an extra halo layer of points at the top and bottom of each sub-domain for data exchange. 
The time loop has a single global kernel call for the calculation of both ϕ and c at the next time step.
Indeed, compared to a classical antitrapping formulation <cit.>, the absence of a term containing ∂ϕ/∂ t in the equation for the evolution of c (namely the antitrapping term) allows to solve both eq:dpdt and eq:dcdt in the same kernel regardless of parallel thread synchronization, hence slightly improving computational efficiency.
After the execution of the main kernel, the halo grid data is updated via direct GPU-GPU memory copy, thus avoiding costly GPU-to-CPU or CPU-to-GPU data transfer.



  §.§.§ Parameters.
 
Material properties and numerical parameters considered in the phase-field simulations are summarized in Table <ref>.





We approximate the WE43 alloy as a binary Mg-3wt%Nd alloy. 
This is motivated by the fact that only yttrium and neodymium have nominal concentrations above one percent in mass (namely with c_ Y∈[1.4-4.2] wt% and c_ Nd∈[2.5-3.5] wt% <cit.>) and the experimental observation of higher Nd segregation levels <cit.> (see section <ref>, figure <ref>).
In fact, microstructural characterization (see section <ref>) shows that most of the Y is trapped within the protective Y_2O_3 passivation layer of the powder, which is not melted but only broken down into flakes during LPBF, hence not participating directly in solute partitioning during solidification.
Moreover, while both Y and Nd exhibit a significant extension of their solubility range in Mg under rapid solidification conditions <cit.>, the maximum equilibrium solubility (at eutectic temperature) of Nd is substantially lower (4.6 wt%) than than of Y (≈14 wt%) <cit.>, with a lower partition coefficient of the former.
These observations point toward Nd as the most likely species to predominantly participate in the banding instability.
Furthermore, the low nominal alloy concentration also provides confidence in the appropriateness of the PF model, originally intended for relatively dilute binary alloys.

Phase diagram features (T_ M, k_ e, m_ e) are directly taken from the Mg-Nd binary diagram <cit.>.
The Gibbs-Thomson coefficient is calculated as Γ=γ_0T_ M/(ρ L_ f), using an interfacial excess free energy γ_0=0.0899J/m^2 calculated (for pure Mg) via molecular dynamics (MD) calculations <cit.> and ρ and L_ f listed in Table <ref> <cit.>.
The solute (Nd) diffusion coefficient in the liquid state is approximated as D_ L≡ D_ Mg^ Nd≈ (r_ Mg/r_ Nd) D_ Mg^ Mg <cit.>, considering a self-diffusion coefficient D_ Mg^ Mg≈5.95×10^-9m^2/s at T=T_ M (averaged between values calculated in <cit.> and <cit.>) and using Goldschmidt (metallic) radii r_ Mg=160pm and r_ Nd=181pm <cit.>.
Average values and anisotropies of the solid-liquid interface excess free energy, γ(n), and kinetic coefficient, μ(n), were calculated by MD for pure Mg <cit.>, here assuming that coefficients ε^μ_40 and ε^μ_60 in the expression of μ(n) can be neglected, as was shown to be the case for γ(n) <cit.>.

Regarding numerical parameters, we consider a diffuse interface width with S=5 and a corresponding coefficient A=11 assessed by an analysis of k(V) and m(V) curves across the entire velocity range (see section <ref>). 
The grid spacing is set as h=0.8W to allow a reasonable spatial description of the ϕ profile across the interface, and the explicit time step is chosen for stability of the explicit scheme (with a safety factor 0.6).



  §.§.§ Simulations. 
  
With the PF model and parameters described above, we performed rapid solidification simulations considering (G,V) conditions sampled along the tail (i.e. solidifying) part of the melt pool in the calculated thermal field (section <ref>). Before describing the main 2D PF simulations, we comment on the convergence analysis carried out to assess an appropriate value of the interfacial diffusion coefficient A for a given value of S=W/W_0>1, following the same procedure as presented in <cit.> (see section A of the Supplemental Material therein).

In order to identify an appropriate value of A for a given S, we consider a one-dimensional problem, and compute the k(V) and m(V) curves for the reference case at S=A=1 over a broad range of V.
As described in <cit.>, assuming that the phase field ϕ remains close to its theoretical stationary profile ϕ_0(x)=-tanh[x/(√(2)W)], the corresponding concentration profile c(x) can be obtained by numerical integration of


    d c/d x=(c_∞-c)V/D_ Lq(ϕ)+b c d g(ϕ)/d x


with c(±∞)=c_∞.
For a given V, the partition coefficient k(V) can thus be obtained from the ratio of concentrations on the solid (c_s) and liquid (c_l) sides of the interface, respectively approximated as the nominal concentration c_∞ and the maximum value of c(x).
The corresponding slope of the liquidus line is then given by the integration of


    m(V)/m_ e = b/(1-k_ e)c_l∫_-∞^+∞g'(ϕ)cd ϕ/d x  x  .


Once k(V) and m(V) calculated for a given set of (A,S) values, their deviation from the reference case at A=S=1, denoted δ_k and δ_m, are quantified as the (discrete) average of 
√((k_S-k_1)^2)/k_1 
and 
√((m_S-m_1)^2)/m_1, 
respectively, over the calculated velocity range, with k_1 and m_1 the reference values of k and m calculated with A=S=1 and k_S and m_S the values calculated for S≠ 1 <cit.>.
Here, equations eq:dcdt_1D_steady and eq:mV_1D_steady are integrated numerically up to machine precision using Matlab library  <cit.> for 20 different velocities V = 0.024× 1.5^{0, 1, …, 19} m/s. 
An optimal value of A is then chosen such that the averages of δ_k(A) and δ_m(A) are minimized (see <cit.> for further details). 

Then, we perform 2D PF simulations of solidification considering a one-dimensional temperature gradient under various (G,V) conditions listed in Table <ref>.
These conditions are sampled between the solidus and liquidus temperatures along the tail of the melt pool in the FE-predicted thermal field (see section <ref>).
PF simulations start with a planar solid-liquid interface located at the liquidus temperature, with a small random perturbation of amplitude h (i.e. 4 nm) in the x direction.
The phase field is initialized with the 1D stationary solution ϕ(x)=-tanh[(x-x_ i)/(√(2)W)] and the solute concentration field with the corresponding equilibrium c(x)=c_∞exp{ b [g(ϕ(x))+1]}.
The most advanced (highest T) point along the solid-liquid interface is kept at the same location x_ i within the domain during the entire simulation by considering a moving frame in the x direction.
The location x_ i is chosen to leave a sufficient distance to the upper x domain boundary, i.e. a length L_ liquid that allows for a decaying concentration profile toward c_∞ in the occurrence of solute partitioning at the interface.
The liquid length L_ liquid ahead of the interface was set to 1 µm, which represents at least 19 times the characteristic diffusion length D_ L/V for all cases.
The Mg preferred growth direction ⟨112̅0⟩ is aligned with the temperature gradient direction (x+).
The domain size is (10 µm)^2, i.e. a 2500^2 grid, with periodic boundary conditions in the lateral (y) direction and Neumann (no-flux) in the growth (x) direction.
The total simulated times (t_ f), listed in Table <ref>, were verified to be sufficient to yield a steady microstructure behavior.

Each PF simulation was performed on one high-performance computing node.
Each node is equipped with 4 Nvidia RTX 3090 GPUs and an AMD (16-core EPYC 7282) CPU processor.
The completion of each simulation lasted between 100 minutes (for t_ f=0.04 ms at V=1.0 µm/s) and 16.5 hours (t_ f=0.4 ms at V=0.1 µm/s).








§ RESULTS



Figure <ref> shows SEM micrographs in various locations of the samples.
Figure <ref> shows a TEM image, as well as corresponding compositional EDS maps for Mg, Nd, Gd, Y and O elements.









Bright irregular flakes measuring a few micrometers correspond to yttrium oxide particles, with small amount of Gd and Nd, but completely depleted of Mg (Figure <ref>). These are inherited from the passivation layer of the feedstock powder, aimed to mitigate the high flammability of Mg. This high melting temperature oxide layer is fragmented during LPBF. Resulting flakes are transported by Mg fumes and end up mixed with the molten metal. Therefore, they do not directly participate in the solidification process during LPBF, besides acting as solid particles potentially interfering with the advance of the solidification front.

Within the regions solidified during LPBF (i.e. excluding the Y oxide flakes), micrographs show the co-existence of patterned (dendritic or seaweed) regions and banded regions, both of which are decorated with fine Mg-RE intermetallic precipitates (bright fine particles).
Images taken parallel to the build direction, upwards in figure <ref>(a)-(d), mostly exhibit bands normal to the build direction, where one might expect a strong temperature gradient during AM.
Images showing a surface normal to the build direction, i.e. figure <ref>(e)-(f), also exhibit both banded and non-banded regions.
Upon closer look into the banded structures (figure <ref>(b) and (f)), the bright regions of the bands appear composed of fine patterns and intermetallic precipitates. 
The EDS compositional maps in Figure <ref> confirm that bands are composed of chemically homogeneous dark regions alternating with bright regions particularly rich in Nd and also exhibiting significant segregation of Gd, and Y.


Figure <ref> illustrates the results of the FE simulation. 
Panel (b) shows superimposed location of the FE-predicted liquidus (▵ symbols) and solidus (▿ symbols) isotherms and their elliptical approximations eq:rl-eq:rs (solid lines).
The melt pool dimensions are d_ L=55.33 µm, d_ S=60.46 µm, l_ L=119.48 µm, and l_ S=132.06 µm.
The resulting temperature gradient and solidification velocity along the melt pool, following equations eq:Gellip-eq:Vellip appear in figure <ref>(c).
The magnitude of the temperature gradient evolves from 10.3 K/µm at the bottom (θ=90^∘) to 4.2 K/µm at the tail (θ=0^∘) of the melt pool, whereas the solidification velocity jointly goes from 0 to V_ L=1 m/s.
Arrows in figure <ref>(b) and symbols in figure <ref>(c) mark the location and (G,V) conditions considered in the PF simulations.
Black arrows and filled symbols identify configurations leading to banded microstructures (otherwise, arrows are white and symbols are open), illustrating the predominance of the banding phenomenon in most of the melt pool.





Figure <ref> illustrates the results of the PF preliminary 1D convergence analysis. 
Panels (a) and (b) respectively show the velocity dependence of the interface solute partition coefficient, k(V), and liquidus slope, m(V), for different values of (S,A)=(3,6), (5,11), and (8,18.5) (symbols), compared to the reference case with (S,A)=(1,1) (solid gray line) and predictions of the CGM eq:cgm:kv-eq:cgm:mv (dashed black line). 
The latter was calculated considering a diffusion velocity through the interface V_ d = 0.356 V_ d^0ln(1/k_ e)/(1-k_ e) = 4.06 m/s with V_ d^0 = D_ L/W_0 = 5.26 m/s, and a partial solute drag with α=0.645 (see <cit.>).
The selected pairs of (S,A) values were obtained by a minimization of δ_k(A) and δ_m(A) curves, illustrated for S=5 in figure <ref>(c).





Figures <ref> and <ref> show the results of the PF simulations, with a growth direction (x+) represented from left to right. 
Figure <ref> illustrates the final microstructure obtained for (G,V) conditions sampled along the melt pool tail (Table <ref>).
Figure <ref> focuses on the interface region during one banding period for the case at V=1m/s and G=4.2K/µm.
Both figures show the (ϕ=0) location of the interface (black line) superimposed on the solute concentration c(x,y) color map.












§ DISCUSSION



Microstructural characterization (figures <ref> and <ref>) unambiguously identify the origin of banded regions as rapid solidification.
These microstructures are nearly identical to those reported from other rapid solidification processes <cit.>, in particular in Al-Cu alloys <cit.>.
They emerge from an alternation of microsegregation-free (partitionless) planar growth and dendritic (or seaweed) growth with solute partitioning.
Overall, banded microstructures are predominant and ubiquitous throughout our additively manufactured Mg alloy samples.
They do not appear to have been significantly affected by the thermal cycling due to the laser path after solidification. The as-printed microstructures (Figure <ref>) retain features distinctive of rapid solidification banded microstructures, with solute segregation confined to the thin patterned regions  (Figure <ref>). Our interpretation is that the thermal peaks induced by the laser are too fast to trigger any significant diffusion mechanism in the solidified microstructure. However, after a possible post-printing heat treatment, we expect that the fine precipitate, present in both dendritic and banded grain could be fully dissolved into the matrix (but not the Y oxides) <cit.>.
The predominance of banded microstructures is consistent with the results of thermal and microstructure selection simulations, indicating that the majority of the melt pool (black arrows in figure <ref>(b)) solidifies under conditions promoting a banding instability, here found for V≥0.4 m/s (figure <ref>).

The PF convergence study confirms that, for this alloy, PF simulations retain a good accuracy on both k(V) (figure <ref>(a)) and m(V) (figure <ref>(b)) over the entire range of considered interface velocity.
As mentioned in <cit.>, while the asymptotic limits at V≪ V_ d and V≫ V_ d match the CGM predictions, a small deviation appears at intermediate velocities, which is expected when comparing a diffuse interface approach (PF) with a sharp-interface theory (CGM).
Still, the close matching of both k(V) and m(V) curves for S>1 with the reference case at S=1 shows that the interface may be upscaled by up to a factor 8 with negligible deviation, hence ensuring that the calculations performed here with S=5 and A=11 remain very close to those at S=1. 
In 2D with an explicit time stepping scheme, the grid coarsening by a factor 5 leads to a reduction of computational cost by ≈5^4=625.
Results for S=8 even suggest that further upscaling by 8^4=4096 is reasonable, which could possibly allow quantitative simulations at larger scale, e.g. possibly up to the full melt pool scale <cit.>.

The mechanism leading to banded microstructures, illustrated in figure <ref>, is the burgeoning instability described in <cit.>, whereby dendritic tips accelerate, morphologically destabilize, and grow almost (yet not completely) isotropically, until they form a new planar front while decelerating. 
The sudden tip acceleration and morphological destabilization coincides with nearly complete solute trapping, as the solid phase appears to grow at the same solute concentration as the surrounding liquid.
Within its basal plane, the hcp Mg crystal has six ⟨112̅0⟩ preferred growth directions, one of which is aligned with the temperature gradient direction, resulting in secondary sidebranches growing at a ±60^∘ angle from the main (x+) growth direction.
Branches growing at a 60^∘ angle with respect to the x+ direction are almost systematically at the origin of the burgeoning instability. 
This may be attributed to their relatively higher growth velocity in their specific growth direction, effectively faster due to the 60^∘ angle, and hence triggering the instability earlier than dendritic tips oriented along the temperature gradient direction.
PF results show that the dendrite acceleration first leads to significant solute trapping while still dendritic (indicated by the lower shade of blue of dendrite tips in figure <ref>(a)-(b)) before the onset of the burgeoning interface instability.
The instability arises almost simultaneously from several dendritic tips (figure <ref>(b)-(c)).
The locations where different partitionless regions impinge upon one another form slightly segregated channels, and a small deviation from a planar front (figure <ref>(e)-(g)), which often corresponds to the location of the emerging dendrites in the next banding cycle (see attached videos).

The effect of the growth velocity V is apparent in both dendritic and banding regimes (figure <ref>).
In the dendritic regime (V≤0.3 m/s), the interface is mostly composed of dendritic doublon structures, with a notable decrease of primary spacing as V increases.
Doublons consist of two asymmetric fingers, with a mirror symmetry between them, growing cooperatively with a narrow liquid channel running between the two fingers along the symmetry axis.
They constitute a known steady-state growth solution promoted by low anisotropy and high undercooling <cit.>.
Theoretical and computational studies <cit.> of the resulting sharp-interface problem have shown that, when both doublon and dendrite steady-state solutions exists, doublons grow faster than regular symmetric dendrites, thus giving the doublon an advantage in terms of microstructure selection.
For low interface anisotropy, doublons have been described as the “building block” of seaweed structures <cit.>. However, the transition from dendrite to seaweed was shown to occur without the formation of doublon or multiplet when the interface dynamics is strongly dominated by kinetics <cit.>.
Here, the propensity to form doublons may be attributed to a combination of (i) the low interface energy anisotropy, (ii) the 2D simulations (doublons were mostly observed in 2D thin sample experiments <cit.>, while 3D calculations <cit.> and experiments <cit.> demonstrated the possibility of more complex multiplet structures), and (iii) the high temperature gradient G. 
Further PF calculations on a broader range of G (to be reported elsewhere) suggest that a transition toward a seaweed pattern occurs between 10 and 50 K/µm, which seems consistent with the substitution of symmetric dendrites by doublons at G≈10 K/µm.
The transition between dendritic and seaweed growth modes was shown to be discontinuous (i.e. first order), with a range of parameters where both branches of solutions may be (meta)stable <cit.>. It was also shown that seaweed patterns are promoted by a high temperature gradient <cit.>. 
The transition to the seaweed regime and the presence of doublons might also be promoted by the sixfold anisotropy, although this hypothesis remains to be studied in further depth.

Our results suggest a transition velocity between dendrites/doublons and banded patterns between 0.3 and 0.4 m/s.
This is lower than the absolute stability predicted by linear stability analysis eq:absstab, V_ a≈ 2.41 m/s using k(V) and m(V) from figure <ref> (i.e., for simplicity, using the 1D approximation ϕ(x)≈ϕ_0(x) and directly integrating eq:dcdt_1D_steady-eq:mV_1D_steady) or V_ a≈ 1.85 m/s using the CGM eq:cgm:kv-eq:cgm:mv.
However, it is worth noting that the precise maximum stable velocity for a steady dendrite (or doublon) can only be probed by first stabilizing a dendritic pattern, before slowly increasing the V until its morphological destabilization (see <cit.>), and that the resulting threshold V is expected to be highly sensitive to fluctuations/noise. 
Yet, here we start with a noisy quasi-planar interface, which appears more likely to yield an oscillatory instability rather than dendrites or doublons across a velocity range where both steady-state solutions may effectively exist.
As such, the value identified here (0.3∼0.4 m/s) can only be viewed as a lower bound for the maximum dendrite/doublon steady velocity, which may still reasonably approximate non-ideal (noisy) conditions in experiments.
Still, while the order of magnitude is realistic, a deeper study on the effect of noise (and a possible V-ramp approach as in <cit.>) would be required to definitely conclude on the maximum steady dendrite/doublon velocity for this alloy.

At V=0.4 m/s, while most of the domain undergoes the cyclic banding instability, some dendritic/doublon patterns may coexist and grow jointly with the bands during a few (here, 2 or 3) successive banding cycles (see top left and bottom right of figure <ref>(d) and corresponding attached video).
In the banding regime (V>0.4 m/s), an increase of growth velocity leads to a reduction of the width of the patterned/segregated region, consistent with experimental observations in rapidly solidified Al alloys <cit.>.

PF simulations (figure <ref>) predict the emergence of microstructures qualitatively similar to those observed in experimental micrographs (figure <ref>).
Since the selected PF model considers only the formation (solidification) of a primary phase, simulation results show the expected solid-liquid interface pattern, discarding the formation of secondary solid phases that would nucleate in high-Nd segregated liquid regions.
The quantitative comparison between experiments and simulated microstructures is reasonable, with a band spacing of about one to two micrometers.
However, some experimental band spacings (e.g. in figure <ref>(a),(c)) appear lower than predicted by PF simulations.
This discrepancy may be attributed to two main factors. 
First, the fact that simulations are two-dimensional, while the experiments are inherently 3D, with a much more complex temperature history than what we considered here, and the fact that micrographs are only representing a 2D slice of the resulting three-dimensional banded microstructures.
Second, and importantly, the fact current PF simulations neglect the effect of the latent heat release, which was shown to have a central role in the selection of the band spacing <cit.>.
The latent heat was also shown to potentially induce a lateral spreading of the interface instability, resulting in microsegregation-free bands slightly tilted with respect to the temperature gradient <cit.>, which might also affect the near simultaneous burgeoning instability observed throughout the sample in the present simulations.
A deeper investigation on the effect of the latent heat is currently underway and will be reported elsewhere.
Other possible sources of discrepancies could stem from high uncertainties on material parameters, such as the liquid diffusion parameter, here roughly estimated as 5.26×10^-9m^2/s, which seems high but of the typical order of magnitude for liquid metallic alloys (𝒪(10^-9) m^2/s).
The high flow velocities expected in the melt pool may also play a role, in particular resulting in distinct local solidification conditions (e.g. solute concentration, temperature gradient). 
However, given the small volume of the melt pool and the typically low diffusion length for such high solidification rates, we expect that, at the scale of our simulations, the assumption of nearly flat thermal field and a homogeneous (well stirred) liquid remain reasonable.




§ SUMMARY & PERSPECTIVES



We studied the formation of rapid solidification microstructures in 3D-printed (LPBF) samples of biomedical-grade Mg alloy (WE43), combining advanced characterization, finite element thermal simulation at the laser track scale, and phase-field simulations of rapid solidification at the scale of the solid-liquid interface pattern.

Our experiments confirm the prevalence of banded structures composed of pattern-free microsegregation-free regions alternating with patterned (dendrite or seaweed) regions with solute segregation and intermetallic precipitates. 
These microstructures are clearly indicative of a rapid solidification regime at an intermediate velocity V between dendritic and planar growth, previously reported in many rapidly solidified alloys <cit.> — but, to the best of our knowledge, never discussed so far in additively manufactured samples.
Our simulations confirm that the LPBF range of solidification conditions, namely temperature gradient G and growth velocity V, promote the occurrence of banded microstructures throughout the majority of the melt pool (except at its very bottom where the local V≤ 0.3 m/s).

The underlying mechanism of the banding phenomenon was confirmed to be the burgeoning instability of primary (or secondary) dendrite tips recently revealed via PF simulations of rapid solidification in thin Al-Cu TEM samples <cit.>.
For a Mg crystal with sixfold anisotropy in the basal plane, this instability originates principally from sidebranches, growing at a 60^∘ angle from the primary (x+) growth direction, due to their relatively higher velocity respective to their growth direction.
The interface pattern at a V just below the onset of the banding instability exhibit a high density of doublons, attributed to the simulation being two-dimensional, the low capillary anisotropy, and the high temperature gradient G.
As expected, an increase of V leads to a reduction of primary spacing in the dendritic/doublons regime, and to a reduction of the width of the regions exhibiting segregation and microstructural features in the banded growth regime. 
The effect of G was not analyzed in further details, since our PF simulations do not explicitly include the effect of latent heat release and diffusion during solidification, which is known to have a dominant effect over that of the higher-scale temperature gradient <cit.>. 
In spite of this limitation, the present study clearly demonstrates the potential of physics-based simulations in predicting the emergence of rapid solidification microstructures, particularly in the context of metal AM.

Perspectives from this work are multiple. 
The most direct follow-up (and focus of ongoing work), is the incorporation of the latent heat  in the PF model, as well as a deeper study into its effect on the selection of microstructure.
The main associated challenge comes from the different time scales of heat and solute transport, which needs to be addressed computationally, for instance using multiple spatial grids <cit.>.
Another important contribution would be the study of the burgeoning instability and banding mechanism in three dimensions. 
Our results suggest that, with the present PF model, the diffuse interface width W might possibly be coarsened even further than the value S=5 used here and in previous works <cit.> while retaining a reasonable quantitative accuracy on k(V) and m(V).
Therefore, although it would certainly represent a substantial computational effort, three-dimensional simulations might be within reach, considering the ever-increasing size and speed of high-performance computing hardware (in particular GPUs).

Finally, from a longer-term applied viewpoint, as rapidly solidified Mg alloy were reported to exhibit significantly different mechanical <cit.> and corrosion <cit.> behavior than regular casting products, one may expect that locally changing the solidification conditions could enable a local modification of microstructures and properties, and ultimately performance of additively manufactured Mg parts.
Combining microstructural models with track-scale or macroscopic thermal models of additive manufacturing could provide a greater understanding of the effect of various parameters (e.g. laser power/speed/path, part geometry, alloy chemistry) on the local distribution of microstructures within a part, e.g. an implant.
Yet, while rapidly solidified Mg alloys may exhibit enhanced mechanical properties and corrosion resistance, a deeper fundamental understanding at microstructural level remain lacking (e.g., to the best of our knowledge, the effect of banded microstructures on corrosion properties is not clearly established).
Once both process-microstructure and microstructure-properties links are clearly established, one may envision computational strategies based on optimization algorithms, not only for topology, but also for alloy chemistry and local processing parameters in order to achieve localized optimal properties.
In this context, we trust that the quantitative prediction of microstructures and microstructural length scales, and their correlation with local solidification conditions, may play a key role in guiding the design of novel 3D-printed biomedical implants with locally tuned mechanical and corrosion properties.




The authors wish to thank Alain Karma and Kaihua Ji for enlightening discussions and clarifications regarding the phase-field model. 
DT, RT, and AKB gratefully acknowledge support from the Spanish Ministry of Science and the European Union NextGenerationEU, through a Ramón y Cajal Fellowship [RYC2019-028233-I] (DT), the MiMMoSA project [PCI2021-122023-2B] (RT), and the HexaGB project [RTI2018-098245] (AKB).
ADB acknowledges financial support from the European Commission, Horizon Europe, Marie Skłodowska-Curie Actions through the M3TiAM project [Grant agreement 101063099].
ML and JMA acknowledge support by the European Union’s Horizon Europe research and innovation programme [Grant agreement 101047008], the Spanish Government [PID2019-109962RB-100 and PID2022-138417OB-C21], and Madrid Regional Government [Y2020/BIO-6756: i-MPLANTS-CM].




§ SUPPLEMENTARY MATERIAL


The joint multimedia files contains videos corresponding to each panel of figure <ref>.
For cases leading to banding instability (i.e. for V≥0.4 m/s) the last 10% of the simulation time is displayed with a lower time step between frames (1/1000 of the simulated duration, instead of 1/100 for the first 90%) in order to reveal temporal details of the fast growth dynamics during banding.



§ APPENDIX: DEVELOPED EXPRESSIONS OF ANISOTROPIC TERMS



The anisotropic excess free energy of the solid-liquid interface is expressed as

    γ(θ) = γ_0 a_ s(θ) 
      .

We simplify its notations and that of its derivatives with respect to θ, simply writing a_ s instead of a_ s(θ). 
Hence, from eq:aniso_g, we have

    a_ s   = 1+ a_20 + a_66cos(6θ) 
    
    a_ s'    = ∂_θ a_ s = -6 a_66sin(6θ) 
    
    a_ s”   = ∂_θθ a_ s =-36 a_66cos(6θ)

with

    a_20   = -ε^γ_20/4√(5/π)
    
    a_66   = ε^γ_66/64√(6006/π)  .

The function a_ s and its derivatives can be expressed directly as a function first order spatial derivatives of the phase field ϕ using trigonometric identities

    cos(6θ)    = (2cos^2θ -1)[ 4(2cos^2θ -1)^2-3 ] 
    sin(6θ)    = 2sinθcosθ  (3-16sin^2θcos^2θ)

with 

    cosθ   = ∂_xϕ/|∇⃗ϕ|
    sinθ   = ∂_yϕ/|∇⃗ϕ|  .


Then, following the same method as in <cit.> (Appendices therein), the first term on the right-hand side of eq:dpdt may be expanded as 

    ∇⃗·[ W(𝐧)^2∇⃗ϕ] = 
    S^2W_0^2 [ ∇⃗(a_ s^2) ·∇⃗ϕ + a_ s^2∇^2ϕ]

with the first term in brackets expressed as

    ∇⃗(a_ s^2) ·∇⃗ϕ 
    = (∂_xϕ∂_xθ + ∂_yϕ∂_yθ )2 a_ s'a_ s  .

Moreover, the penultimate term on the right-hand side of eq:dpdt may be expressed as 

    ∑_η=x,y[ ∂_η( |∇⃗ϕ|^2 W(𝐧) ∂ W(𝐧)/∂ (∂_ηϕ)) ] 
    = (∂_xϕ∂_yθ - ∂_yϕ∂_xθ )(a_ s”a_ s+a_ s'^2)
      ,

where, since θ=tan^-1(∂_yϕ/∂_xϕ), we have <cit.> 

    ∂_xθ   = ∂_xyϕ∂_xϕ-∂_xxϕ∂_yϕ/|∇⃗ϕ|^2
    ∂_yθ   = ∂_yyϕ∂_xϕ-∂_xyϕ∂_yϕ/|∇⃗ϕ|^2  .

Thus, all terms in eq:a12 and eq:a13, and hence all terms in eq:dpdt may be expressed as local first and second order spatial derivatives of ϕ.




§ REFERENCES
 
iopart-num


