Tailoring magnetic hysteresis of Fe-Ni permalloy by additive manufacturing: Multiphysics-multiscale simulations of process-property relationships

Designing the microstructure of Fe-Ni permalloy by additive manufacturing (AM) opens new avenues to tailor the materials' magnetic properties. Yet, AM-produced parts suffer from spatially inhomogeneous thermal-mechanical and magnetic responses, which are less investigated in terms of process simulation and modeling schemes. Here we present a powder-resolved multiphysics-multiscale simulation scheme for describing magnetic hysteresis in materials produced via AM. The underlying physical processes are explicitly considered, including the coupled thermal-structural evolution, chemical order-disorder transitions, and associated thermo-elasto-plastic behaviors. The residual stress is identified as the key thread in connecting the physical processes and in-process phenomena across scales. By employing this scheme, we investigate the dependence of the fusion zone size, the residual stress and plastic strain, and the magnetic hysteresis of AM-produced Fe21.5Ni78.5 permalloy on beam power and scan speed. Simulation results also suggest a phenomenological relation between magnetic coercivity and average residual stress, which can guide the magnetic hysteresis design of soft magnetic materials by choosing appropriate AM-process parameters.


Introduction
The Fe-Ni permalloy has been widely studied in recent decades owing to its extraordinary magnetic permeability, low coercivity, high saturation magnetization, mechanical strength, and magneto-electric characteristics.The material has been widely employed in conventional electromagnetic devices, such as sensors and actuators, transformers, electrical motors, and magnetoelectric inductive elements.Fe-Ni-based permalloys modified with additives are also promising candidate materials for multiple novel applications, such as wind turbines, all-electric vehicles, rapid powder-conversion electronics, electrocatalysts, and magnetic refrigeration [1,2,3,4].
Due to the increasing importance of additive manufacturing (AM) technologies, the possibilities of designing soft magnetic materials by AM have been explored in a number of studies [5,6,7,8,9,10,11,12].However, due to the delicate interplay of process conditions and resulting properties, there are several open questions that need to be answered in order to obtain AM-produced Fe-Ni permalloy parts with the desired property profile.Magnetic properties of the Fe-Ni system depend on the chemical composition, as depicted in Fig. 1a.Fe 21.5 Ni 78.5 is typically selected targeting the chemical-ordered low-temperature FCC phase (also known as awaruite, L1 2 , or γ phase, as the phase-diagram shown in Fig. S1a), which possesses a minimized coercivity and peaked magnetic permeability [13,8].The main problem is to increase the generation of γ phase using AM, as the phase transition kinetics from the chemical-disordered high-temperature FCC phase (also known as austenite, A1, or γ phase) to the chemical-ordered γ phases is extremely restricted [14].On the laboratory timescale, growing the γ phase into considerable size requires annealing times on the order of days [15,16,17,18,19].Due to the rapid heating and cooling periods, only the γ phase exists in AM-processed parts [9].Combining in-situ alloying with AM methods allows to stabilize the γ phase in printed parts [8], but the restricted kinetics still limits the growth of a long-range ordered phase [14].Another question concerns the influences of associated phenomena on magnetic hysteresis.Although there are studies on the effects of crystallographic texture and orientation [11,20].Unfortunately, a pivotal discussion regarding interactions among residual stress, microstructures, and physical processes leading to the magnetic hysteresis behavior is missing.
Recently, it has been shown that soft magnetic properties can be designed by controlling magneto-elastic coupling [21,22,23].Along these lines, the residual stress caused by AM and the underlying phase transitions could be key for tuning the coercivity in an AM-processed Fe-Ni permalloy.Micromagnetic simulations by Balakrishna et al. [22] presented that magneto-elastic coupling plays an important role in governing magnetic hysteresis, as the existence of the pre-stress shifts the minimized coercivity from the composition Fe 25 Ni 75 to Fe 21.5 Ni 78.5 , while the magneto-crystalline anisotropy is zero at Fe 25 Ni 75 but non-zero at Fe 21.5 Ni 78. 5 .Based on a vast number of calculations, the dimensionless constant (C 11 − C 12 )λ 2 100 /2K u ≈ 81 was proposed as a condition for low coercivity along the 100 crystalline direction for cubic materials (incl.Fe-Ni permalloy) [23].Here, C 11 and C 12 are the components of the stiffness tensor, λ 100 is the magnetostrictive constant, and K u is the magneto-crystalline constant.Nevertheless, the influence of the magnitude and the states of the residual stress were not comprehensively analyzed.Yi et al. discussed magneto-elastic coupling in the context of AM-processed Fe-Ni permalloy for the first time [21].The simulations were performed with positive, zero, and negative magnetostrictive constants under varying beam power.The results showed that the coercivity of the Fe-Ni permalloy with both positive and negative magnetostrictive constants rises with increasing beam power.In contrast, the permalloy with zero magnetostrictive constant showed no dependence on the beam power.This demonstrates the necessity of magneto-elastic coupling in tuning coercivity in AM-produced permalloy and illustrates the potential of tailoring properties of permalloys via controlling the residual stress during AM.
Understanding the residual stress in AM and its interactions with other physical processes, such as thermal and mass transfer, grain coarsening, and phase transition, is never the oak that felled at one stroke.Taking the popular selective laser melting/sintering (SLM/SLS) method as an instance, the temperature gradient mechanism (TGM) explains the generation of residual stress by considering the heating mode and the cooling mode [24,25,26].The heating mode presents a counterbending with respect to the building direction (BD) of newly fused layers (Fig. 1b).This is because the thermal expansion in an upper overheated region gets restricted by the lower old layer/substrate.Plastic strain can also be generated due to the activated plasticity of the material and compensate the local stress around the heat-affected zone.The cooling mode, in contrast, presents a bending towards the BD of the newly fused layer due to the thermal contraction between the fusion zone and the old layer/substrate.
It should be noted that TGM only provides a phenomenological aspect by employing the idealized homogeneous layers.In the practical SLM/SLS, varying morphology and porosity of the powder bed create inhomogeneity in not only the temperature field but also the on-site thermal history on the mesoscale (10-100 µm), inciting varying degrees of thermal expansion, and eventually leading to the development of the thermal stress in various degree.Stochastic inter-particle voids and lack-offusion pores also create evolving inhomogeneity in material properties on the mesoscale [27,28], leading to the shifted local conditions for developing the residual stress as interpreted by TGM.In other words, mesoscopic inhomogeneity and coupled thermo-structural evolution should act as a long-range factor in residual stress development.On the other hand, due to the relatively smaller lattice parameter, the continuous growth of the γ phase would also result in the increasing misfit stress between itself and the chemical-disordered γ matrix [19,18] (as also presented in Fig. 1c).Therefore, the nanoscopic solid-state phase transition also contributes to the development of residual stress as a short-range fluctuation, which is almost effectless to the mesoscopic phenomena yet still influences the local magnetic behavior [22].To sum up, the residual stress from AM processes that eventually affects the magnetization reversal via magneto-elastic coupling should already reflect such long-range (morphology and morphology-induced chronological-spatial thermal inhomogeneity) and short-range factors (misfit-induced fluctuation).This is the central challenge that this work addresses.
In this work, we developed a powder-resolved multiphysics-multiscale simulation scheme to investigate the hysteresis tailoring of Fe-Ni permalloy by AM under a scenario close to practical experiments.This means that the underlying physical processes, including the coupled thermal-structural evolution, chemical order-disorder transitions, and associated thermo-elastoplastic behaviors, are explicitly considered and bridged by accounting for their chronological-spatial differences.The influences of processing parameters (notably the beam power and scan speed) are analyzed and discussed on distinctive aspects, including the size of the fusion zone, the development of residual stress and accumulated plastic strain, the γ/γ transition under the residual stress, and the resulting magnetic coercivity of manufactured parts.It is anticipated that the presented work could provide transferable insights in selecting processing parameters and optimizing routine for producing permalloy using AM, and deliver a comprehensive understanding of tailoring the hysteresis of soft magnetic materials in unconventional processing.

Multiphysics-multiscale simulation scheme
In this work, we consider SLS as the AM approach due to its relatively low energy input as compared to other methods, like SLM.SLS allows us to obtain a stable fusion zone and thus to gain better control of the residual stress development since we don't need to consider melting and evaporating processes and the associated effects, such as the keyholing and Marangoni convection.The microstructure of SLS-processed parts is porous, which allows us to explore the effects of lack-of-fusion pores on the development of residual stress and plastic strain.Based on an overall consideration of all possible phenomena involved in the SLS of the Fe-Ni permalloy, two chronological-spatial scales are integrated in this work: (i) On the mesoscale, with the characteristic length of several 100 µm, powders are fused/sintered around the laser spot, creating the fusion zone.Featured phenomena such as partial/full melting, necking, and shrinkage among powders can be observed.High gradients in the temperature field are also expected due to laser scanning and rapid cooling of the post-fusion region.By choosing a typical scan speed of 100 mm s −1 , this stage lasts only 10 ms.
(ii) On the nanoscale with a characteristic length well below 1 µm, the chemical order-disorder (γ/γ ) transition can be observed once the on-site temperature is below the transition temperature, as presented in Fig. 1c.Owing to the difference in thermodynamic stability, redistribution of the chemical constituents by inter-diffusion between γ and γ phases is coupled to the phase transition.Due to the extremely restricted kinetics, it would cost several hundred hours of annealing to have the γ phase formation in a mesoscopic size [15,29,18].
We explicitly consider a three-stage processing route consisting of an SLS, a cooling, and an annealing stage.As shown in the inset of Fig. 2, the domain temperature will rise from a pre-heating temperature (T 0 ) during the SLS stage.After that, the whole powder bed would gradually cool down to T 0 .Finally, the processed powder bed enters the annealing stage at T 0 where the γ/γ transition continues.The cooling stage lasts three times longer than the SLS stage, and the sequential annealing stage takes far longer than the two other stages.Taking a 500 µm scan section with a typical scan speed of 100 mm s −1 as an example, the SLS and cooling stages would last 5 and 20 ms, respectively, and the annealing time is on the order of 100 hours.Remarkably, the inhomogeneous and time-varying temperature field in the mesoscopic powder bed can be treated as uniform and nearly constant for the nanoscopic γ/γ transition, as the local heating and cooling stages induced by laser scan are negligible compared to the time required for the γ/γ transition.Nonetheless, the long-term mechanical response (notably the residual stress) remains after the first two stages.It would further influence the γ/γ transition during the annealing stage and the resultant magnetic hysteresis behavior by electro-magnetic coupling [30,21,22].
The simulations are arranged in a subsequent scheme to recapitulate the aforementioned characteristics on different scales while balancing the computational cost-efficiency, as shown in Fig. 2. Accepting that heat transfer is only strongly coupled with microstructure evolution (driven by diffusion and underlying grain growth) but weakly coupled with mechanical response during the SLS-process stage, we employ the non-isothermal phase-field model proposed in our former work [27] to simulate the coupled thermo-structural evolution, and perform the subsequent thermo-elasto-plastic calculations based on the resulting transient mesoscopic structure (hereinafter called mesostructure) and temperature field from the SLS simulations.In other words, mechanical stress and strain are developed under the quasi-static microstructure and temperature field.This is based on the fact that the thermo-mechanical coupling strength is negligible for most metals [31], unlike the strong inter-coupling among mass as well as heat transfer and grain growth [32,33].From a kinetic point of view, the propagation of elastic waves is generally faster than thermal conduction and diffusion-based mechanisms, like grain coarsening and solid-state phase transition.Next, taking the nanoscopic domains that are sufficiently small and can be regarded as "homogenized points" on the mesostructure, we transfer the historical quantities on the sampled coordinates, notably the temperature and stress histories, to the subdomains as the transient uniform fields and perform the non-isothermal phase-field simulations of the γ/γ transitions.This also means that mesoscopic temperature gradients are disregarded in the nanoscopic simulations in this work.Finally, we connect the nanoscopic Ni concentration and stress field to the magnetic properties, incl.the saturation magnetization M s (X Ni ), magnetocrystalline anisotropic strength K u (X Ni ), and magnetostriction constants λ 100 (X Ni ) and λ 111 (X Ni ), and perform the magnetoelastic coupled micromagnetic simulations for the local hysteresis.Both non-isothermal phase-field and thermo-elasto-plastic models are numerically implemented by the finite element method (FEM), which allows handling the geometric complexity and adaptive meshing at considerable numerical accuracy.The micromagnetic models, on the other hand, are implemented by the finite difference method (FDM) to allow for GPU-accelerated high-throughput calculations [34].Simulation domains are collectively illustrated in Fig. S2.Apart from the main workflow, the proposed scheme also involves other methods, such as the discrete element method (DEM) and CALculation of PHAse Diagrams (CALPHAD) approach, to deliver information such as the powder size (R i ) and center (O i ) distributions and the thermodynamic/kinetic parameters that required in the simulations.Details regarding the modeling and the simulation setup are explicitly given in the Method section.

SLS single scan simulations and coupled thermal-microstructural evolution
Here we present the results of SLS single-scan simulations of a Fe 21.5 Ni 78.5 powder bed in Argon atmosphere.A powder bed with an average thickness of h = 25 µm is placed on a substrate with the same composition and thickness of 240 µm.The powder size distribution is presented in Fig. S1b.The simulation domain has the geometry of 250 × 500 × 300 µm.The melting point of Fe 21.5 Ni 78.5 is T M = 1709 K, and the initial temperature of the powder bed is set as the pre-heating/annealing temperature T 0 = 0.351T M = 600 K.The temperature at the substrate bottom is set as T 0 throughout the simulations.D FWE2 = 200 µm (i.e., full-width at 1/e 2 ) is adopted as the nominal diameter of the laser spot, within which around 86.5% of the power is concentrated.The full width at half maximum intensity (FWHM) is then calculated as D FWHM = 0.588D FWE2 = 117.6 µm, characterizing 50% power concentration within the spot.
Fig. 3a shows the evolution of simulated microstructure for a single scan of P = 30 W and v = 100 mm s −1 .In the overheated region, particles may be fully/partially melted.The tendency to reduce the total surface energy leads to the motion of the localized melt flowing from convex to concave points, which contributes to the fusion of the powders.In regions with T ≤ T M , no melting occurs.However, the temperature of the particles is sufficiently high to induce diffusion, evidenced by the formation of necking between adjacent particles.Since the local temperature is well above the γ/γ transition temperature during the SLS processes, there is no ordering at this stage.
The temperature profiles of the powder bed for different beam powers and scan speeds are presented in Fig. 3b, where the temperature field strongly depends on particle morphology, as the isotherms are concentrated around the surface concaves and sintering necks among particles.Apart from this, one can observe relatively dense isotherms at the front and the bottom of the overheated region, indicating a large temperature gradient.While the laser spot is moving, this temperature gradient becomes smaller, as the isotherms tend to be sparser.This indicates a fast heating process followed by slow cooling.Comparing the beam spot for different processing parameters presented in Fig. 3b shows that increasing the beam power and/or decreasing the scan speed enhance the heat accumulation at the beam spot, resulting in the overheated region with increasing size.For the same scan speed v = 100 mm s −1 , increasing the beam power from P = 30 W (Fig. 3a) to 35 W (Fig. 3b) leads to more significant overheated region.On the other hand, for the same power of P = 30 W, increasing the scan speed from v = 100 mm s −1 (Fig. 3a) to 150 mm s −1 (Fig. 3b) leads to reduced overheated region.

Development of stress and plastic strain during SLS single scan
In order to analyze the stress evolution during SLS, we use isotropic hardening plasticity to describe Fe 21.5 Ni 78.5 with temperaturedependent mechanical properties, including thermal expansion coefficient α, Young's modulus E, yield stress σ y , and hardening tangent modulus E t .Spatial interpolation of the mechanical properties according to the order parameter ρ, which is ρ = 1 in the materials and ρ = 0 in the atmosphere/pores, is also performed to consider structural inhomogeneities due to pore formation.Details are described in the section on Methods.The domain-average quantities are defined as σD e = Ω ρσ e dΩ Ω ρΩ , pD e = Ω ρp e dΩ Ω ρΩ where Ω is the simulation domain volume, and ρ is the order parameter indicating the substance.σ e is the von Mises stress and p e is the accumulated (effective) plastic strain.To eliminate the boundary effects, which lead to heat accumulation on the boundary that intersects with the scan direction, the simulation domain with a geometry 250 × 250 × 250 µm is selected from the center of the domain for processing with the transient T and ρ fields mapped, as shown in Fig. 2. Fig. 4a presents the evolution of the domain-average von Mises stress during the SLS and cooling stages.σ e develops with the temperature rise due to the generation of laser-induced heat in the powder bed.However, once the laser spot moves into the domain, followed by the overheated region, the σD e drops along with TD continuing rising to the maximum.This is because the points inside the overheated region lose their stiffness, as the material is fully/partially melted, thereby presenting zero stress (Fig. 4b 1 ).Surroundings of the overheated region also present relatively low stress due to the sufficient reduction of stiffness at high temperatures.When the laser spot moves out of the domain, the thermal stress starts to develop along with the cooling of the domain (Fig. 4b 2 -b 4 ).The stress around the concave morphologies on the powder bed, incl.concaves on the surface and sintering necks among particles, rises faster than the one around the traction-free convex morphologies on the surface and unfused powders away from the fusion zone.This may attribute to the locally high temperature gradient around the concave morphologies during the SLS process and, thereby, the strong thermal traction.At the end of the cooling, the convex morphologies and unfused powders have relatively lower stress developed, while the locally concentrated stress can be observed around the concave morphologies on the powder bed, like the surface concaves and sintering necks among particles (Fig. 4b 4 ).Due to the convergence of the σ e vs. time in the cooling stage, the stress at the end of the cooling stage (t = 20 ms) will be regarded as the residual stress in the following discussions.
The development of the accumulated plastic strain p e presents an overall increasing tendency vs. time during the SLS and cooling stages, as shown in Fig. 4c.To emulate the effects of full/partial melting, we reset the p e in the overheated region (T > T M ).This reset does not, however, influence the accumulation of p e outside of the overheated region (Fig. 4d 1 ), distinguished from σ e that suffers reduction not only inside but also outside of the overheated region due to loss of stiffness at high temperature.As a result, the growth of the pD e slows down only when the laser spot moves into the domain rather than tends to reduce as σD e (t).The continuous accumulation of p e also results in the distinctively concentrated p e at the fusion zone's outer boundary (Fig. 4d 2 -d 4 ), attributing to the high temperature gradient at the front and bottom of the overheated region during the SLS where the existing high thermal stress locally activates the plastic deformation and then contribute to the rise of the p e .The same reason can be used to explain the concentrated p e around the pores and concave morphologies like sintering neckings near the fusion zone.In contrast, unfused powders and substrate away from the fusion zone present nearly no accumulation of p e , as the onsite thermal stress due to relatively low local temperature is not high enough to initiate the plastification of the material.

Nanoscopic γ/γ transition with residual stress
Before sampling many points inside the fusion zone for studying the subsequent γ/γ transition and carrying out micromagnetic simulations, we first examined five selected points from the middle section of the SLS simulation domain, counting its repeatability along the x-direction (SD).Profiles of σ e and p e are presented on this selected middle section in Fig. 5a.These five points are taken from the profiling path along z-direction (BD) and y-direction, as σ e and p e along these two paths are explicitly presented in Fig. 5b 1 -b 2 .σ e gradually rises along z-profiling path, as the gap between the normal stresses (specifically between σ xx and σ zz , and between σ yy and σ zz since there are almost no differences between σ xx and σ yy .σ e reaches a peak around the fusion zone boundary (FZB) and then decreases.Notably, there is a reverse of σ zz from positive (tension) to negative (compression) across the FZB.Similarly, p e increases to a peak around the FZB and then decreases along z-profiling path, and the reversed components of ε pl are observed across the FZB.This implies the bending of the SLS-processed mesostructure caused by the thermal contraction between the fusion zone and substrate, where the high temperature gradient is depicted.Along the y-profiling path, both σ e and p e present no monotonic tendencies, receiving the influences from the morphologies, yet still reach peaks around the FZB, respectively.The Lamé's stress ellipsoids are illustrated in Fig. 5c for visualizing the stress state of selected points with the directions of the principal stress denoted.Points near the surface and in the fusion zone (P 1 , P 2 , P 4 , and P 5 ) are under tensile stress states, while the point P 3 near the FBZ has one negative principle stress (compressive) along BD due to the bending caused by the thermal contraction in the fusion zone.The Ni concentration and nanoscopic stress redistribution due to the γ/γ transition are simulated in the middle section perpendicular to SD as well (Fig. S2a 3 ).It is also coherent with the diffuse-controlled 2D growth of γ phase experimentally examined by the Johnson-Mehl-Avrami-Kolmogorov (JMAK) theory [29].We initiate the γ nuclei randomly but with a minimum spacing of 100 nm according to the experimental observation in Fig. 1c using Poisson disk sampling [35].A relatively longer annealing time of 1200 h to obtain sufficient γ phase formation, as presented in Fig. 5d.The Ni concentration (X Ni ) at the centers of the grown γ phase is relatively low, close to the equilibrium value 0.764 at T γ/γ = 766 K.With the growth of the γ phase, Ni is accumulated at the interface due to the relatively large γ/γ interface mobility compared with the inter-diffusive mobility of Ni species, agreeing with its diffuse-controlled growth examined experimentally [29].Among the points, P 3 with one principal compressive stress has relatively larger γ grown after annealing.This is due to the growing γ phase with relatively smaller lattice parameters (in other words, inciting shrinkage eigenstrain inside γ phase) being mechanically preferred with compressive stress.As P 2 and P 5 with similar stress states, similar γ phase formations are shown.P 1 has relatively less γ phase grown, as it has the highest normal stresses as tensile (σ xx , σ yy , and σ zz ) compared to other points, as shown in Fig. 5b 1 .

Magnetic hysteresis behavior of stressed nanostructures
The nanoscopic distributions of stress and Ni concentration are imported from the results of the γ/γ transition, in which both the long-range and short-range factors are embodied.The magneto-elastic coupling is implemented by considering an extra term f em in the magnetic free energy density along with the exchanging, magneto-crystalline anisotropy, magnetostatic, and Zeeman contributions [30,36,37], as described in the section on Methods.To consider contributions from normal and shearing stresses, a homogeneous in-plane configuration for the normalized magnetization m is oriented in an angle of ϑ = 45 • relative to the crystalline orientation u, which is assumed to be the z-direction.Meanwhile, an effective coupling field B em = −∂ f em /∂(M s m) is calculated for analyzing the magneto-elastic coupling effects in vectorial aspect (Fig. 5e).It should be noticed that within the range of segregated X Ni from 0.781 to 0.810 as shown in Fig. 5d, λ 100 ranges from 1.274 ×10 −5 to 5.249 ×10 −6 , presenting a positive shearing magnetostriction.On the other hand, λ 111 inside of the γ phase (X Ni = 0.781 ∼ 0.785) and in the γ matrix (X Ni = 0.785) has positive values, ranging from 2.41 ×10 −6 to 1.96 ×10 −6 , while X Ni shifts from 0.781 to 0.785.On the γ/γ interface with X Ni = 0.810 , however, a negative value of λ 111 = −8.429×10 −7 is obtained.This implies that there is a negative normal magnetostriction on the γ/γ interfaces and a positive normal magnetostriction in the bulk of γ and γ phases.Due to the existing shrinkage eigenstrain, there are locally high f em contributions inside the γ phase compared to the γ matrix, which causes strong magneto-elastic coupling effects.Among all points P 1 -P 5 , P 1 has comparably lower magneto-elastic coupling in both γ and γ phases.Owing to similar stress states, P 2 , P 4 , and P 5 have similar profiles of f em and B em , where P 5 has slightly stronger coupling effects.Remarkably, P 2 , P 4 , and P 5 all have the local B em lying at an angle about 135 • , as emphasized by dashed-dotted circle Fig. 5e.As for P 3 , however, the angle between m and B em further increases to 142 to 165 • together with larger magnitude (|B em |) comparing to other points, demonstrating enhanced reversal effects to the m.
The hysteresis curves of the nanostructures also reflect the X Ni -dependence of K u , λ 100 and λ 111 .For comparison, the hysteresis curves simulated on a stress-free reference with a homogeneous X Ni = 0.785 are also plotted.In order to take numerical fluctuations into account, ten cycles of the hysteresis were examined for each nanostructure/reference with the averaged one presented in Fig. 5f.Notice that K u ranges from −0.137 to −0.364 kJ m −3 according to Fig. 1a, which can be regarded as the easy-plane anisotropy as the magnetization prefer to orient in the plane perpendicular to BD. Comparing selected points, P 2 and P 5 have similar coercivity, which is 0.55 mT for P 2 and 0.58 mT for P 5 , as these two points have similar stress states.P 4 and P 1 have the coercivity of 0.26 mT and 0.09 mT, respectively.However, the coercivity at P 3 is infinitesimal.This is phenomenologically due to the strong coupling field B em that reverses m at the relatively low external field, as shown in Fig. 5e.As the coupling field B em is directly related to the stress state, the relatively strong shear stress implied by the Lamé's stress ellipsoids at P 3 (Fig. 5c) may be one of the reasons for the infinitesimal coercivity.This should be further examined in future studies.

Discussion
In the following, we discuss the relation between fusion zone geometries and the processing parameters, notably the beam power P and scan speed v.It should notice that the size of the fusion is directly linked to the overheated region, which varies with the different combinations of P and v, as shown in Fig. 3b.The usage of the conserved OP ρ in representing the powder bed morphologies and the coupled kinetics between ρ and local temperature allows us to simulate the formation of the fusion zone in a way close to the realistic setup of SLS.Here the indicator ξ is utilized to mark the fusion zone, which is initialized as zero and would irreversibly turn to one once the temperature is above T M [38].Two characteristic sizes, namely the fusion zone width b and the fusion zone depth d, are defined by the maximum width and depth of the fusion zone, as shown in the inset of Fig. 6.Fig. 6a and Fig. 6b present the maps of b and d vs. P and v, with the isolines indicating the identical volumetric specific energy input U V that is calculated as where the width of the laser scan track w takes D FWHM = 117.6 µm and the average thickness of the powder bed h = 25 µm.Since 50% total power is concentrated in the spot with D FWHM , the efficiency takes η = 0.5.According to the observed geometry of the fusion zone, the maps can be divided into three regions: P and v located in the region (R1) would result in a continuous fusion zone, as shown in Fig. 6c 1 , c 3 -c 6 .The depth of the fusion zone in (R1) normally penetrates the substrate, implying the formation of a considerable size of the overheated region, within which the melting-resolidification would take the dominant role.P and v located in the region (R2) generate small and discontinuous fusion zones, as two classical geometries as shown in Fig. 6c 2 and c 7 .Its limited size implies the typical partial melting and liquid-state sintering mechanism, as the melt flow is highly localized and can only help the bonding among a subset of powder particles.It is worth noting that these discontinuous fusion zones should attribute to the thermal inhomogeneity induced by local stochastic morphology rather than the mechanisms like the Plateau-Rayleigh instability and balling, where a significant melting phenomenon is required [39,40,41].In the region labeled as (R3), no fusion zone is generated under the chosen P and v, and the solid-state sintering process remains dominant in the powder bed.Fig. 7a and b present the maps of average residual stress σe and plastic strain pe inside the fusion zone vs. P and v, with the specific energy input U V .The average residual stress and plastic strain inside the fusion zone are calculated with the fusion zone indicator ξ from the relations σe = Ω ξσ e dΩ Ω ξ dΩ , pe = Ω ξ p e dΩ Ω ξ dΩ . ( Since the properties inside the fusion zone are focused, only the results with P and v located in the continuous fusion zones (R1) and discontinuous ones (R2) are selected and discussed.Generally, the increase of the σe follows the direction of increasing specific energy input U V , leading to the enlargement of the fusion zone.In other words, increasing the size of the fusion zone receives larger contractions between itself and the substrate, leading to the rise of the residual stress inside the fusion zone, comparing Fig. 7c 1 -c 3 and c 1 , c 4 -c 5 .It worth noting that σe presents a rapid increase on U V below 120 J mm −3 , as shown in Fig. S3.When U V > 60 J mm −3 , there is almost no increasing of σe , which is reflected as the sparse contours beyound the isoline U V = 60 J mm −3 .This implies the saturation of residual stress in the fusion zone when U V is sufficiently large, despite the fusion zone would continue enlarging along with the further increase of U V .The map of pe presents a ridge at P = 30 W, which is different from the monotonic dependence of σe on U V .This means the pe reaches a minimum at P = 30 W for every selected v.This may reflect the competition between enlarging fusion zone and the accumulation of plastic strain in the fusion zone.In the low-power range, the accumulation of plastic strain is slower than the growth of the fusion zone, as the pe reduces along with increasing P. When P > 30 W, the accumulation of plastic strain is faster than the growth of the fusion zone, as the pe increases along with increasing P. Interestingly, such a tendency is not observed in reducing v with fixing P, as the pe always grows monotonically with decreasing v at every selected P.
Moreover, σ e drastically rises across the former boundary between powder bed and substrate for all selected combinations of P and v, as shown in Fig. 7c 1 -c 5 .This is due to the different mechanical responses between the porous powder bed and homogeneous substrate to the thermal stress formed together with the overheated zone.As the size of the fusion zone enlarged, the high concentration of σ e extends further into the interior of the substrate, as the thermal contraction becomes enhanced between the fusion zone and the substrate.Meanwhile, a relatively low concentration of p e is observed inside the fusion zone.
And a relatively higher concentration is located at the fusion zone's outer boundary, reflecting the continuing accumulation of certain regions, as shown in Fig. 7d 1 -d 5 .
Many points located on the mid-section of the fusion zones were sampled, considering their repeatability along the xdirection (SD), as shown in Fig. S4.It can also eliminate the boundary effects and distractions from the quantities outside the fusion zone.Nanoscopic γ/γ transition and hysteresis simulations were performed on each point subsequently.The average coercivity of the fusion zone Hc is then calculated directly as the point-wise average of the resulting local hysteresis.At least three hysteresis cycles were performed on each point to reduce the fluctuations due to the numerical scheme.In Fig. 8a, we present the average coercivity Hc map with respect to the P and v.The resulting coercivities under all examined P and v locate in the experimentally measured range (from 0.06 to 4 mT).In the map, the lower-right region shows a smaller coercivity, and the upper-left region shows a higher coercivity.Increasing P from 27.5 to 35 W results in the rise of Hc from 0.35 to 0.41 mT (by 17%) when v = 100 mm s −1 , and decreasing v from 125 to 50 mm s −1 results in the rise of Hc from 0.34 to 0.44 mT (by 29%) when P = 30 W. This is mainly due to the increasing fraction of region with high local coercivity in the fusion zone for increasing P and decreasing v (comparing Fig. 8c 1 , c 2 , c 3 , and c 1 , c 4 , c 5 , respectively).It is also evident that the high local H c points also possess a high local volume fraction of γ phase Ψ γ (comparing Fig. 8c 1 -c 2 with d 1 -d 5 ).This implies an enhanced magneto-elastic coupling effect at high Ψ γ , as discussed in Fig. 5e.However, the effect from the local residual stress on the pattern of high local H c region should also be stressed, as the points with low local H c around where the local σ e has a drastic change, e.g., the former boundary between powder bed and substrate, comparing (Fig. 8c 1 -c 2 with e 1 -e 5 , especially c 3 with e 3 and c 4 and e 4 ).Moreover, it remains unclear for the existing "islands" of low local H c located in the high local σ e and Ψ γ region.They may be subject to a stress state similar to the P 3 in Fig. 5 where nearly zero coercivity is obtained.Nonetheless, a data-driven investigation should be conducted as one upcoming work to connect the local stress state to the coercivity.
To examine the dependence of Hc from the phenomenological aspect, we firstly performed the nonlinear regression analysis of Hc on the specific energy input U V .The results are presented in Fig. 8b.Notably, Hc relates to U V by allometric scaling rule, i.e., Hc = C H (U V ) I H with C H and I H the parameters.The analysis gives the correlation coefficient R 2 = 86.54%with relatively large uncertainty located on low and high U V regions.Regressed I H = 0.31 ± 0.03 also implies a diminishing scaling of Hc by U V , as the increment of Hc decreases with the increase of U V .However, similar to the one for σe , such a simple scaling rule between the specific energy input and coercivity might be challenged since U V may not be able to identify the uniquely H c for the SLS-processed part, as the isoline of U V in Fig. 8 an evidently intersects with the contour of Hc .
Regression analysis of Hc on σ e was also performed.The exponential growth rule was chosen based on the tendency of Hc ( σe ), i.e., Hc = A σ exp( σe S σ ) + H σ with A σ , S σ and H 0 as the parameters.Here we adopt the A σ as the growth pre-factor and S σ as the stress scale.When σe = 0, we have H c | σe =0 = A σ + H 0 ≈ H 0 , meaning H 0 can be regarded as the stress-free coercivity.The analysis gives a relatively higher correlation coefficient R 2 = 91.31%cf. the one on U V , with relatively large uncertainty located on low U V region.Compared to the homogeneous stress-free reference in Fig. 5f that is 0.45 mT, this H 0 is around 24% smaller, owing to the contributions from the infinitesimal-coercivity points.It also presents the rapid growth of Hc after σe = 206 MPa with around 40% increment, demonstrating evident effects on Hc with increasing residual stress.

Conclusion
In summary, processing-property relationship in tailoring magnetic hysteresis of Fe 21.5 Ni 78.5 has been demonstrated in this work by conducting multiphysics-multiscale simulation.The residual stress is unveiled to be the key thread since it readily carries both long-range (morphology and morphology-induced chronological-spatial thermal inhomogeneity) and short-range information (misfit-induced fluctuations) after the processes.Influences of beam power and scan speed have been investigated and presented on distinctive phenomena, including the geometry of fusion zones, the residual stress and accumulated strain, and the resultant coercivity of the manufactured parts.The following conclusions can be drawn from the present work: (i) The simulated mesoscopic residual stress states are coherent with TGM interpretation.Further details like the concentrated stress around the concave morphologies (surface concave, sintering necks, etc) beyond the TGM interpretation are also delivered.The accumulated plastic strain is evidently observed at the fusion zone's outer boundary.
(iii) Large magneto-elastic coupling energy is observed inside γ phase with the corresponding effective field imposing rotating effects on the magnetization.These effects vary point-wisely according to residual stress states and γ phase formation, and eventually lead to different local coercivity.Remarkably, the point around the bottom of the fusion zone is examined to have nearly zero coercivity, which may attribute to the on-site stress state with a principal compressive stress along the building direction with another two tensile ones, i.e., implying a relatively significant shear stress.
(iv) The relation between the average residual stress and coercivity of the fusion zone is examed to follow the exponential growth rule with a correlation coefficient of 91.31%.The stress-free coercivity derived from the exponential law is 0.34 mT, and the rapid growth of average coercivity is observed when average residual stress exceeds 206 MPa, implying a potential threshold of average residual stress for restricting coercivity of the manufactured parts around the stress-free value.On the other hand, average residual stress beyond the threshold can be considered in the context of effectively increasing the resultant coercivity of the manufactured permalloy parts.
Despite the present findings, several points should be further examined and discussed in future works: (i) The conditions deriving the near-zero coercivity should be extensively investigated in the sense of residual stress states rather than its effective value, noting that residual stress is a 2 nd -order tensor.Magneto-elastic coupled micromagnetic simulations should be performed under the classified residual stress states to rationalize the factors leading to the vanishing of coercivity as in the present findings.
(ii) The present findings are only examined at relatively low specific energy input and, correspondingly, low generated residual stress, as the SLS is chosen in this work.It is anticipated to conduct the simulations with relatively high energy input, like SLM, and examine the influences of magneto-elastic coupling on the magnetic hysteresis with comparably higher residual stress cases.Influences on residual stress development and, eventually, the coercivity of manufactured permalloy from multilayer and multitrack AM strategies should also be examined.

Thermodynamic framework
In order to describe the microstructure of an SLS-manufactured Fe-Ni alloy, a conserved order parameter (OP) ρ is employed to represent the substance and atmosphere/pores, and a set of non-conserved OPs {φ ϕ χ } are employed to represent the grains with the superscript ϕ = γ, γ representing the phases and the subscript χ = 1, 2, . . ., N representing the orientations, extended from our former works [27,32].Counting the thermal, chemical, and mechanical contributions, the temperature field T, the strain field ε, and sets of local chemical molar fraction {X ϕ A } of the chemical constituents A = Fe, Ni are considered.On the other hand, as a ferromagnetic material under the Curie temperature T C = 880 K for the composition Fe 21.5 Ni 78.5 , the thermodynamic contribution due to the existing spontaneous magnetization m is also counted.The framework of the free energy density functional of the system is then formulated as follows where f ch represents the contributions from the chemical constituents.f loc and f grad , together as the f intf , presents the contributions from the surface and interfaces (incl.grain boundaries and phase boundaries) [42].f el is the contribution from the elastic deformation, and f mag is the contribution from spontaneous magnetization and magnetic-coupled effects.
It is worth noting that this uniform thermodynamic framework does not imply that a single vast inter-coupled problem with all underlining physics should be solved interactively and simultaneously crossing all involved scales.As sufficiently elaborated in the subsection Multiphysics-multiscale simulation scheme, it is more practical and effective to conduct the multiphysics-multiscale simulations in a subsequent scheme and concentrate on rationalizing and bridging the physical quantities and processes among problems and scales.In that sense, the free energy density functional, originated from Eq. (4), should be sufficiently simplified regarding the distinctiveness of each problem at the corresponding scale.This will be explicitly introduced in the following sections.

Mesoscopic processingssing simulations
Here we consider the SLS processing on a mesoscopic powder bed by using ρ to differentiate pore-substance and {φ ϕ χ } to differentiate polycrystalline orientations.According to the high-temperature phase diagram of the Fe-Ni system [43], the γ phase exists within a relatively large temperature range (from T M = 1709 K to the transition starting temperature T γ/γ = 766 K) for the composition Fe 21.5 Ni 78.5 .On the other hand, since the SLS together with cooling stages would only last a relatively infinitesimal time (on the order of 10 ms) to the following annealing stage (more than 10 h), there is almost no change for γ phase to grow into mesoscopic size.In this regard, we treat all existing polycrystals during the SLS stage as the γ phase.Considering Fe-Ni as a binary system where the constraints X Ni + X Fe = 1 and φ γ χ + φ γ χ = ρ always holds, we then only take the OP set {φ γ χ } (χ = 1, 2, ..., N) as well as ρ for the mesoscopic simulations due to the absence of the γ phase on the mesostructures.The profiles of ρ and {φ γ χ } across the surface and grain boundary between two adjacent γ-grains are illustrated in Fig. 9b.We also take simplified notations X = X Ni in this subsection as the independent concentration indicators, while X Fe = 1 − X.
Due to the co-existing of substance (γ-grains with Ni composition X 0 = 0.785) and pores/atmosphere, the chemical free energy density should be formulated as where h ss and h at are monotonic interpolation functions with subscripts "ss" and "at" representing the substance and pore/atmosphere and are assumed to have the polynomial forms as The temperature-dependent chemical free energy f γ ch is modeled by the CALPHAD approach with , where f γ ref is the term corresponding to the mechanical mixture of the chemical constituents (in this case, Fe and Ni), f γ id is the contribution from the configurational entropy for an ideal mixture, f γ mix is the excess contribution due to mixing, and f γ mag is the contribution due to the magnetic moment.The parameters fed in Eq. (6), including the atom magnetic moment β γ , the Curie temperature T γ C , and the interaction coefficient L γ , are described in the way of Redlish-Kister polynomials [44] which is generally formulated for a binary system as ) represents the Inden polynomial, obtained by expanding the magnetic specific heat onto a power series of the normalized temperature T/T γ C [45,46].R is the ideal gas constant.V sys m is the molar volume of the system.All the thermodynamic parameters for the CALPHAD are obtained from Ref. [43] while the molar volume of the system is obtained from the database TCFE8 from the commercial software Thermo-Calc [47].
Since the variation of Ni composition is negligible in between T M and T γ/γ , we pursue a simple but robust way of implementing f ch under a drastically varying T during the SLS stage.Taking T M as referencing temperature, Eq. ( 5) is then re-written as where f T M ref is a referencing chemical free energy density at T M , which can be omitted in the following calculations.c r is a relative specific heat landscape, i.e., c r (ρ, T) = h ss (ρ)c γ v (T) + h at (ρ)c at v (T) with the volumetric specific heat for γ grains and pores/atmosphere.Notably, c γ v can be thermodynamically calculated as follows at a fixing pressure p 0 and composition X 0 .
It should be noticed that the c γ v obtained by Eq. (8) has a discontinuous point at T C , which is due to the 2 nd -order Curie transition, as shown in Fig. S5a.L is the latent heat due to the partial/full melting, which is mapped by the interpolation function h ml .Here h ml adopts a sigmoid form with a finite temperature band ∆ T which reaches unity once T → T M and is smooth enough to ease the drastic change in f ch .
On the other hand, to explain the free energy landscape across the surface and γ/γ interface (or γ grain boundary) under varying temperatures, we adopt the non-isothermal multi-well Landau polynomial and gradient terms from our former works [27,32], i.e., with W γ/γ and κ γ/γ are temperature-independent parameters obtained from the surface and γ/γ interface energy Γ sf , Γ γ/γ and diffuse-interface width γ/γ , and τ sf (T) and τ γ/γ (T) are the dimensionless tendencies inherited from the temperature dependency of Γ sf and Γ γ/γ , i.e., along with the constraint among parameters for having the sample profile of ρ and φ γ χ across the surface [27], i.e., In this work, we give γ/γ = 2 µm, the temperature-dependent Γ sf and Γ γ/γ are presented in Fig. S5c.The total free energy density landscape at stress-free condition ( f tot = f ch + f intf ) is illustrated in Fig. 9c.We can tell that the term f ch modifies the relative thermodynamic stability of the substance by shifting the free energy minima via temperature changes.In contrast, γ grains at the same temperature do not show a difference in stability until the on-site temperature of one is changed.The governing equations for the coupled thermo-structural evolution are formulated as follows [27,28] where Eq. (11a) is the Cahn-Hilliard equation with the mobility tensor M specifically considering various mass transfer paths, incl.the mobilities for the mass transfer through the substance (ss), atmosphere (at), surface (sf) and grain boundary (gb).As elaborated in our former work [27], the localized melt flow driven by the local curvature is also modeled by one effective surface mobility M eff ml .M is then formulated as [32] with the 2 nd -order identity tensor I and projection tensors T sf and T gb for surface and grain boundary, respectively [32,48].The T-dependent values for M sf , M gb , M ss , and M eff ml are presented in Fig. S5d.The interpolation functions on the surface and γ grain boundaries are defined as Eq. (11b) is the Allen-Cahn equation with the scalar mobility L, which is derived from the γ grain boundary mobility G γ/γ as [49,50] which is also presented in Fig. S5d.Eq. (11c) is the heat transfer equation that considers the laser-induced thermal effect as a volumetric heat source q v .
in which p xy indicates the in-plane Gaussian distribution with a moving center r O (v, t).P is the beam power and v is the scan velocity with its magnitude v = |v| as the scan speed.The absorptivity profile function along depth da/dz is calculated based on Refs.[27,51].The phase-dependent thermal conductivity tensor is formulated in a form considering the continuity of the thermal flux along the normal/tangential direction of the surface [48,52], i.e., with where K ss and K at are the thermal conductivity of the substance and pore/atmosphere.N is the 2 nd -order normal tensor of the surface [48,33].Thermal resistance on the surface and γ/γ interface are disregarded, and will be presented in the upcoming works.While temperature-dependent K ss takes the linear form in this work (Fig. S5b), K at specifically considers the radiation contribution via pore/atmosphere and is formulated as where K 0 ≈ 0.06 W m −1 K −1 is the thermal conductivity of the Argon gas, F = 1/3 is the Damk öhler view factor [53], and is the Stefan-Boltzmann constant.rad is the effective radiation path between particles, which usually takes the average diameter of the powders [54].
As the boundary conditions (BC), no mass transfer is allowed on all the boundaries of the mesoscopic domain, which is achieved by setting Neumann BC on ρ as zero.The temperature at the bottom of the substrate mesh (z min ) is fixed at T 0 via Dirichlet BC on T. Heat transfer is allowed only via the pore/atmosphere, achieved by the combined BC of convection and radiation and masked by h at ), as illustrated in Fig. S2b 1 .

Mesoscopic thermo-elasto-plastic simulations
As elaborated in the subsection Multiphysics-multiscale simulation scheme and our former work [28], the subsequent thermoelasto-plastic simulation was carried out for the calculation of the thermal stress and deformation of the mesostructures from the non-isothermal phase-field simulations of SLS processing.The transient temperature field and substance field ρ are imported into the quasi-static elasto-plastic model as the thermal load and the phase indicator for interpolating mechanical properties.Adopting small deformation and quasi-static assumptions, the mechanical equilibrium reads where σ is the 2 nd -order stress tensor.The top boundary is set to be traction free, and the other boundaries adopt rigid support BCs, which only restrict the displacement component in the normal direction of the boundary Fig. S2b 2 .
Taking the Voigt-Taylor interpolation scheme (VTS), where the total stress is interpolated according to the amount of the substance and pore/substance across the interface, i.e., σ = h γ σ γ + h γ σ γ , while assuming identical strain among phases [55,56,57].In this regard, the stress can be eventually formulated by the linear constitutive equation where the 4 th -order elastic tensor is interpolated from the substance one C ss and the pores/atmosphere one C at , i.e., In this work, isotropic mechanical properties are considered.C ss is calculated from the Youngs' modulus E(T) and Poisson's ratio ν.In contrast, C at is assigned with a sufficiently small value to guarantee the numerical convergence.The thermal eigenstrain ε th is calculated using the interpolated coefficient of thermal expansion, i.e., where α ss is obtained from temperature-dependent molar volume of the binary system V sys m , i.e., Meanwhile, for plastic strain ε pl , the isotropic hardening model with the von Mises yield criterion is employed.The yield condition is determined as with where σ e is the von Mises stress.s is the deviatoric stress, σ y is the yield stress when no plastic strain is present.The isotropic plastic modulus H can be calculated from the isotropic hardening tangent modulus E t and Young's modulus E as H = EE t /(E − E t ).These temperature-dependent mechanical properties are collectively presented in Fig. S6.p e is the accumulated plastic strain, which is integrated implicitly from the plastic strain increment ε pl obtained from the radial return method [58,59].It is worth noting that the plastic strain is reset as zero once T > T M in emulating the effect of partial/full melting.

Nanoscopic chemical order-disorder (γ/γ ) transition
Once the temperature drops below T γ/γ , the historical quantities (notably T(t) and σ(t)) will be sampled and imported to the nanoscopic domain for the subsequential γ/γ simulations.Here we consider the nanoscopic (γ/γ ) transition that occurs localized inside one γ grain, where ρ and only one φ γ χ are unity while other OPs are zero.This means the orientation information of the γ grain and the influences from the surface and γ/γ interface have been omitted on this scale.The growing γ is also known to be orientation-coherent based on the experimental observations [19,18].In this regard, the subscript χ, indicating the different grain orientations, is dropped in the following discussions.Magnetic contribution f mag is also dropped since the magnetic-field-free γ/γ transition is scoped in this work.The profiles of φ γ and φ γ across the γ/γ interface are illustrated in Fig. 9b.Similarly, we take simplified notations in this subsection, such as φ = φ γ and X = X Ni as the independent phase and concentration indicators, while X Fe = 1 − X and φ γ = 1 − φ.
Due to the co-existing of both γ and γ once the temperature is below T γ/γ = 766 K, the chemical free energy density should be formulated as where h γ and h γ are monotonic interpolation functions and can adopt the polynomial formulation as Similarly the elastic contribution is formulated as [60,61] f el (T, {ε where On top of the CALPHAD modeling of the free energy density of the chemical-disordered γ phase, the four-sublattice model is employed to describe the phase with chemical ordering.The model takes the element fractions Y (s) i (i = Fe or Ni indicating the chemical constituents, s = 1, 2, 3, 4 indicating the sublattice site, see inset of Fig. 9c) in each sublattice as the inner degree-offreedoms, representing with where f FCC ijkl are the free energies of the stoichiometric compounds with only one constituent (i, j, k, l =Ni or Fe) occupied on each site [43].L FCC (s)ijkl is the interaction parameter, corresponding to the mixing of constituents on the s-th site while others (u, v, w = s) are with the fractions Fe = 1 should be applied to guarantee the conservation of atom.It is also worth noting that due to the thermodynamic equivalence of four sublattice sites, the equivalence of f FCC ijkl and L FCC (s)ijkl regarding the combination of sublattice constituents must be considered, as explicitly explained in Ref. [43].All the thermodynamic parameters for the CALPHAD are obtained from Ref. [43].In Fig. S8.we present the calculated f γ ch and f γ ch from the T γ/γ to the pre-heating temperature T 0 = 600 K with the varying equilibrium concentration of X γ e and X γ e , the site element fraction Y (s) i , and the calculated phase fraction ψ γ and ψ γ under the equilibrium.On the other hand, since there is only the γ/γ coherent interface, the non-isothermal local and gradient free energy density are then formulated in the typical double-well fashion, i.e., with adapting the same non-isothermal form as the one used in the SLS simulations with the dimensionless temperature tendency τ γ/γ (T).The temperature-independent parameters W γ/γ and κ γ/γ are obtained from the interface energy Γ γ/γ and diffuse- interface width γ/γ , i.e., noting that the relation of γ/γ here corresponds to the case when adjusting parameter as two in Eq. ( 53) of the Ref. [62].In this work, we tentatively take τ γ/γ as one and estimate Γ γ/γ = 0.025 J m −2 , which is a commonly estimated value for the coherent interfaces and lies between the experimental range from 0.008 to 0.080 J m −2 for the Ni-base alloys [63].The diffuse-interface width γ/γ is given as 5 nm.The total free energy at stress-free condition ( f tot = f ch + f intf ) is illustrated in Fig. 9d, where the free energy density path obeying the mixing rule is also illustrated across the γ/γ interface between two equilibrium phases (i.e., with X γ e and X γ e ).The governing equations of the nanoscopic γ/γ transition is formulated as follows [62,60] Notably, Eq. (28a) embodies the mixing rule of the local Ni concentration X from the phase ones X γ and X γ , considering the γ/γ interface as a two-phase mixture with φ as the local phase fraction.This detaches the chemical and local contributions to the interface energy to allow rescalability of the diffuse-interface width.Eq. (28b) is the constraint to the phase concentration X γ and X γ to obtain the maximum driving force for the interface migration, as briefly elaborated in Fig. S9.In return, the drag effect might be eliminated along with the vanishing of the driving force for trans-interface diffusion [64,65,66,67], which should be specifically evaluated and discussed for the γ/γ transition in the Fe-Ni system.The diffusive mobility M X here is directly formulated by the atom mobilities M Fe and M Ni in the FCC lattice considering the inter-diffusion phenomena [68], i.e., and the interface migration mobility M φ is derived by considering the thin-interface limit of the model and the interface migration rate that was originally derived by Turnbull [65,50], i.e., where the dimensionless Cahn number Ch = γ/γ / ˆ γ/γ , characterizing the degree of the rescaling of the diffuse-interface width γ/γ from the realistic interface width ˆ γ/γ , which is estimated as ˆ γ/γ = 5 ā with ā the average lattice parameter of the γ and γ phases based on the experimental observation [69].Length of the burgers vector is also calculated from ā by |b| = ā/ √ 2. Φ γ/γ is a newly defined thermodynamic factor with the estimated interface concentration X γ/γ as The detailed derivations are shown in the Supplementary Note 2. The temperature-dependent atom mobility M Ni (T) and M Fe (T) are obtained from the mobility database MOBFE3 from the commercial software Thermo-Calc [47].
In this work, it should be highlighted that a temperature-dependent dimensionless calibration factor ω(T) is additionally associated with the atom mobilities, which is utilized to be calibrated from the experimentally measured γ/γ transition with respect to time at various temperatures.The calibrated atom mobility is then shown as (noting A = Fe, Ni) Based on the Arrhenius relation on temperature for M A (T) and M * A (T) [68,70], this ω(T) is postulated to follow the Arrhenius relation as well, i.e., ω(T) = ω 0 exp(−Q ω /RT) with the pre-factor ω 0 and the activation energy Q ω .We implemented a simple calibration algorithm by iteratively performing the regression on ω as the time scaling factor to the simulated transient volume fraction of γ phase, i.e., Ψ γ (ωt) with respect to the experimental measurements obtained from [29], as shown in Fig. 10a.The IC of the γ nuclei was generated using Poisson disk sampling [35] with the prescribed minimum nuclei distance according to the observation shown in Fig. 1b.The calibrated ω(T) indeed shows consistency to the Arrhenius relation, confirming our postulate.
As for momentum balance in Eq. (28e), we have to explicitly consider both long-range (morphology and morphologyinduced chronological-spatial thermal inhomogeneity) and short-range factors (misfit-induced fluctuation) factors of the mechanical response on the current scale.In that sense, the stress should be considered in the following form where σ ms comes from the mesoscale and σ is incited due to the misfit of growing γ phase.Assuming the stiffness tensor C ss has no differences between the two phases, we then take a uniform elastic strain that attributes to the mesoscopic stress, i.e., σ ms = C ss : ε ms el .The constitutive relation can then be represented as where ε is the total strain calculated on the nanoscopic domain, and ε γ mis is the misfit strain induced by growing γ phase.ε γ mis is the relative difference between lattice parameters of the γ and γ phases, i.e., ε γ mis = (a γ − a γ )/a γ with a γ and a γ obtained from the temperature-dependent molar volume V γ m and V γ m , respectively.This is presented in Fig. S7b.At T 0 = 600 K, this ε γ mis = −1.32×10 −3 .Alongside with ε ms el as an eigenstrain, the periodic displacement BC are applied to the nanoscopic domain, as shown in Fig. S2b

Micromagnetic hysteresis simulations
Below the Curie temperature, the magnetization of most ferromagnetic materials saturates with constant magnitude (M s ).Therefore in micromagnetics, it is important to have a normalized magnetization vector that is position-dependent, i.e., m.This vector field can be physically interpreted as the mean field of the local atom magnetic moments, but yet sufficiently small in scale to resolve the magnetization transition across the domain wall.However, variation of m(r) across the γ/γ interface is tentatively disregarded as an ideal exchange coupling between two phases.Magnetic properties in the ferromagnetic γ phase are also tentatively assumed to be identical to the ferromagnetic γ at the same Ni-concentration due to the lack of experimental/theoretical investigations on the magnetic properties of individual phases.In other words, only the Ni-concentration dependency of magnetic parameters is explicitly considered in this work, while the exchange constant A ex takes constant as 13 pJ/m [71].In that sense, superscript ϕ, indicating the phase differences, is dropped in the following explanation.We let the orientation u of the nanoscopic subdomain align on the z-direction (BD), and the magnetic free energy density is eventually formulated as with and the magnetostrictive strain ε em on the cubic basis as follows [36,37] Here, f ex is the exchange contribution, recapitulating the parallel-aligning tendency among neighboring magnetic moments due to the Heisenberg exchange interaction.The norm ∇m here represents ∑ j |∇m j | 2 with j = x, y, z and m = [m x , m y , m z ]. f ani represents the contribution due to the magneto-crystalline anisotropy.It provides the energetically preferred orientation to local magnetizations with respect to the crystalline orientation u according to the sign of the K u .f ani represents the contribution due to the magneto-crystalline anisotropy.It provides the energetically preferred orientation to local magnetizations with respect to the crystalline orientation u, concerning the sign of the K u .Defining an orientation angle by ϑ = arccos u • m, the case when K u > 0 leads two energetic minima at ϑ = 0 and π, that is when the magnetization lies along the positive or negative u direction with no preferential orientation, i.e., the easy-axis anisotropy.When K u < 0, the energy is minimized for ϑ = π/2, meaning that any direction in the plane perpendicular to u is thermodynamically preferred, i.e., the easy-plane anisotropy [30], as shown in Fig. 9e.As the resulting X Ni varies from 0.781 to 0.810 as presented in Fig. 5, local K u always takes the negative value in this work.The magnetostatic term f ms counts the energy of each local magnetization under the demagnetizing field created by the surrounding magnetization.The Zeeman term f zm counts the energy of each local magnetization under an extrinsic magnetic field H ext .f em is the contribution due to the magneto-elastic coupling effects.
To simulate the hysteresis behavior of the structure during a cycling H ext , we calculate the magnetization configuration m(r) under every incremental H ext change by conducting the constrained optimization of a stationary Landau-Lifshitz-Gilbert equation, which is mathematically formulated as where α d is the damping coefficient, taking α d = 0.02 [72].This also means that the magnetic hysteresis is evaluated under the quasi-static condition.The simulation domains with the FD grids have the same construction as the FE meshes used in the γ/γ transition simulations to ease the quantity mapping in-between.Periodic BC was applied on the boundaries perpendicular to z-direction by macro geometry approach [73], while Neumann BC was applied on the other boundaries [34].

Implementations and parallel computations
Both non-isothermal phase-field and thermo-elasto-plastic models are numerically implemented by the finite element method within the program NIsoS [27,32], developed by the authors based on the MOOSE framework (Idaho National Laboratory, ID, USA) [75,76].The 8-node hexahedron Lagrangian elements were chosen to mesh the geometry.A transient solver with preconditioned Jacobian-Free Newton-Krylov method (PJFNK) was employed in both models.Each simulation was executed with 96 AVX512 processors and 3.6 GByte RAM per processor based on MPI parallelization.The associated CALPHAD calculations were conducted by open-sourced package PyCALPHAD [77], and the thermodynamic data intercommunication was carried out by customized Python and C++ codes.The DEM-based powder bed generation is conducted by the open-sourced package YADE [27,78].
For SLS simulations, the Cahn-Hilliard equation in Eq. (11a) was solved in a split way.The constraint of the order parameters was enforced by the penalty method.To reduce computation costs, h-adaptive meshing and time-stepping schemes are used.The initial structured mesh is presented in Fig. S2a 1 .The additive Schwarz method (ASM) preconditioner with the incomplete LU-decomposition sub-preconditioner was also employed for parallel computation of the vast linear system, seeking the balance between memory consumption per core and computation speed [79].The backward Euler method was employed for the time differentials, and the constraint of the order parameters was fulfilled using the penalty method.Due to the usage of h-adaptive meshes, the computational costs vary from case to case.The peak DOF number is on order 10,000,000 for both the nonlinear system and the auxiliary system.The peak computational consumption is on the order of 10,000 core-hour.More details about the FEM implementation are shown in the supplementary information of Ref. [27].
For thermo-elasto-plastic simulations, a static structured mesh was utilized Fig. S2a 2 to avoid the hanging nodes generated from the h-adaptive meshing scheme.In that sense, the transient fields T and ρ of each calculation step were uni-directionally mapped from the non-isothermal phase-field results (with h-adaptive meshes) into the static meshes.This is achieved by the MOOSE-embedded SolutionUserObject class and associated functions.The parallel algebraic multigrid preconditioner BoomerAMG was utilized with the Eisenstat-Walker (EW) method to determine linear system convergence.It is worth noting that a vibrating residual of non-linear iterations would show without employing the EW method for this work.The DOF number of each simulation is on the order of 1,000,000 for the nonlinear system and 10,000,000 for the auxiliary system.The computational consumption is on the order of 1,000 CPU core-hour.
For γ/γ transition simulations, a static uniform mesh was utilized Fig. S2a 3 .2 nd backward Euler method was employed.The additive Schwarz method (ASM) preconditioner with the complete LU-decomposition sub-preconditioner was also employed for parallel computation.The simulations were performed in a high-throughput fashion with 100∼1,000 transition simulations as a batch for one set of processing parameters.The DOF number of each simulation is on the order of 1,000,000 for the nonlinear system and 10,000,000 for the auxiliary system.The computational consumption of each simulation is 500 CPU core-hour by average.

Fig. 1 .
Fig. 1.(a) Ni concentration dependent magnetic properties, incl.magneto-crystalline anisotropic constant K u , magnetostriction constants λ 100 and λ 111 , saturation magnetization M s and initial relative permeability µ r , modified with permission from Balakrishna et al. [22].(b) Schematics of temperature gradient mechanism (TGM) in explaining the generation of residual stress in heating and cooling modes, where the tensile σ (t) and compressive σ (c) stress states are denoted.Inset: A bent AM-processed part due to residual stress.Image is reprinted with permission from Takezawa et al. [26] under the terms of the Creative Commons CC-BY 4.0 license.(c) Bright-field image of a Fe-Ni permalloy microstructure after annealing at 723 K for 50 h.Inset: electron diffraction pattern of the microstructure.SEM and electron diffraction images are reprinted with permission from Ustinovshikov et al. [19].
K u (X Ni ) M s (X Ni ) l 1nn (X Ni )s+s Fig.

Fig. 3 .Fig. 4 .
Fig. 3. (a 1 )-(a 4 ) Simulation results of SLS processing of a Fe 21.5 Ni 78.5 powder bed and substrate with beam power 30 W and scan speed 100 mm s −1 at different time points; Figs.(b 1 )-(b 4 ) show results with varying beam power and scan speed at the time point when the laser center locates at x = 310 µm.Overheated regions, where T > T M , are drawn with a continuous color map, while areas with T ≤ T M are shown as isotherms.The laser spot haracterized by D L and D FWHM is also indicated.

Fig. 5 .
Fig. 5. Profiles of (a) the von Mises residual stress σ e and the accumulated plastic strain p e on the middle section perpendicular to the x-direction (SD).The selected points P 1 -P 5 are also indicated.The components and effective values of (b 1 ) the residual stress and (b 2 ) the plastic strain along the z-and y-profiling paths as indicated in (a).(c) The Lamé's stress ellipsoids, representing the stress state at P 1 -P 5 , with the principle stresses denoted and colored (marine: tension; red: compression).(d) Ni concentration profile after annealing for 1200 h at T 0 = 600 K. (e) the magneto-elastic coupling energy f em calculated under a homogeneous in-plane m configuration with ϑ = 45 • w.r.t. the easy axis u.The corresponding effective coupling field B em is also illustrated.(f) The average hysteresis loop of 10 cycles each performed on the nanostructures (NS) at P 1 -P 5 .The reference curve (Ref.) is performed on the stress-free homogeneous nanostructure with only γ phase and the composition Fe 21.5 Ni 78.5 .

1 (Fig. 7 .
Fig. 7. Contour maps of (a) average residual stress σe and (b) average plastic strain pe in the fusion zone w.r.t. the beam power and the scan speed.The dotted lines represent isolines of different specific energy inputs U V .Profiles of (c 1 )-(c 5 ) residual stress and (d 1 )-(d 5 ) plastic strain are plotted for different processing parameters.The boundary of the fusion zone is also indicated.

Fig. 8 .
Fig. 8. (a) Contour map of the average coercivity Hc in the fusion zone w.r.t. the beam power and the scan speed.The dotted lines represent different specific energy input isolines U V .(b) Nonlinear regression of Hc on specific energy input U V and average residual stress σe , with the regression parameters indicated correspondingly.It presents that Hc subjects to scaling rule on U V and to exponential growth law on σe with the correlation coefficient R 2 = 86.54%and 91.31%, respectively.Profiles of (c 1 )-(c 5 ) the local coercivity; (d 1 )-(d 5 ) the local volume fraction of γ phase; and (e 1 )-(e 5 ) the local residual stress on the mid-section of the fusion zone for different processing parameters.

FeFig. 10 .
Fig. 10.(a) Workflow for the calibration of the atom mobility by iteratively performing regression on the factor ω from Ψ γ (ωt) w.r.t. the experimental measured Ψ exp γ (t).The ω that has a difference of less than 0.5% to the last interaction is identified as the converged value for the current temperature.(b) Simulated Ψ γ (t) at 784 K with the mobilities before (ω = 1) and after (ω = 3) calibration cf.Ψ exp γ (t) from [29], with the nanostructures denoted at corresponding time point.(c) Simulated Ψ γ (t) vs Ψ exp γ (t) at different temperatures.Inset: Regression of calibrated ω(T) to the Arrhenius equation, which presents consistency.Notice that both experiments and simulations are performed with composition Fe 25 Ni 75 .