Flash sintering incubation kinetics

The microstructural mechanisms leading to onset of the flash sintering are demonstrated experimentally and theoretically for Yttria Stabilized Zirconia, YSZ. Three regimes leading to flash event are identified: (1) Radiation-dominated regime, where the oven controls the heating of the sintered sample, and a small subset of particle-particle contacts and surfaces of the green body define percolative paths for the charge to flow along and across the interfaces; (2) Transition regime, where charge transport is suppressed across particle contact misorientations and deflects to surficial and small angle particle contact misorientations. As a result, internal Joule heating takes over externally-driven radiation heating. Finally, (3) Percolative regime, where the concentration of oxygen vacancies drastically increases at particle contacts, surfaces, and triple junctions, and enables charge to flow through multiple paths, generating large amounts of Joule heating, resulting in the onset of a flash event. The validated theory sets the stage to rationalize the microstructural evolution and charge transport on a ceramic green body during flash sintering.


INTRODUCTION
The state-of-the-art of energy storage, biomedical, electronic, and high-temperature materials and devices demand a broad palette of ceramic materials with high-performance properties. Typically, these ceramic materials are manufactured using powder compacts fired at high temperatures to enhance mass transport and densification. This classic technique called sintering 1 , consumes large amounts of thermal energy, requires from hours to days of processing, and contributes a great deal to CO 2 emissions, increasing its manufacturing cost and making it detrimental to the environment 2,3 .
Modifications to this basic approach include the application of external mechanical pressure 4 , electric 5 , and magnetic fields 6 to decrease manufacturing cost, applied energy, and processing time. These new techniques enabled an unprecedented degree of control of grain size, porosity, and crystallographic texture. Specifically, flash sintering has emerged as a very promising processing technique 5,[7][8][9][10][11] , for the fabrication of electronic [12][13][14] , semiconducting [15][16][17] , and composite solids 18,19 . Flash sintering was pioneered by Cologna, Rashkova, and Raj, and provides remarkably rapid densification, with processing times in the range of minutes to hours, at furnace temperatures significantly lower than those used by traditional sintering techniques 5 . 3YSZ was the first ceramic material processed using flash sintering 5 , see experimental results shown in Fig. 1. At a heating rate of 10 K/min, for a constant applied electric field of 60 V/cm, the onset furnace temperature for flash is measured as 1310 K. The increase in applied electric field from 60 to 120 V/cm reduces the onset temperature for flash to 1123 K 5 .
Flash sintered specimens are rapidly densified with a sudden increase in electrical conductivity, and is followed by a bright light emission event which is referred to as the flash event 20 . Also, a nonlinear increase in power density is experimentally reported in the 10 to 50 mW/mm 3 range for a wide variety of ceramics, for different temperatures and electric fields 21 . The range of observed power density is attributed to Joule heating, and correlates to the onset of flash 21,22 . As the sample builds up to the flash event, the total conductivity follows an Arrhenius-type behavior as a function of the sample temperature. The onset condition for flash was first predicted using a coarse-grain thermal runaway model developed by Zhang, Jung, and Luo for ZnO 23 , Todd et al. for YSZ 24 , and also Dong and Chen for a portfolio materials 25 . Recent rapid thermal annealing experiments on ZnO 26 , and 3YSZ 27 , without an electric field achieved rapid densification rates analogous to flash sintered samples indicating that the heating rate of the sample is an important controlling factor for fast densification and grain growth rates.
To explain the underlying mechanisms that lead to the flash event, including thermal annealing approaches 28 , a wide range of mechanisms have been brought forward. Raj proposed the possibility of field-induced Frenkel pairs formation at particle contacts to facilitate the mass transport for flash sintering 29 . Naik et al. proposed the nucleation of Frenkel pairs at grain boundaries, but required an unphysically large dielectric constant contrast between grains and grain boundaries in order to contribute an inhomogeneous volumetric dipolar free energy that would drive the sintering process, while the grain boundaries experience large Joule heating 30 . Most recently, Jongmanns et al. suggested through molecular dynamics simulations on aluminum the possibility of large concentrations of Frenkel pairs for sample temperatures above the Debye temperature 31 . However, Schie et al. 32 reported that a substantially large field of 10 GV/m is required to generate Frenkel pairs in ionic ceramics, using HfO 2 as a model material.
In contrast, Narayan proposed that electric fields alter the atomic mobility and enthalpy of migration of the diffusion process. Also, he proposed that the electric field would generate anionic vacancies as a result of the partial reduction of the material, which are preferentially segregated at grain boundaries and dislocation loops 33 . The excess vacancies are proposed to lead to localized Joule heating at grain boundaries, including the possibility of premelting, potentially leading to a flash event 34 .
Another mechanism that has been proposed to rationalize the onset of the flash event is based on the hypothetical development of significant temperature gradients where overheating/premelting of particle contacts define the abrupt increase in the electrical conductivity. Here, particles and particle contacts soften and transform into liquid, potentially leading to a flash event 35 . However, because of the high thermal diffusivity of 3YSZ, heat transfers rapidly from local Joule heating hot spots, large thermal gradients for time-scales greater than~1 μs cannot be physically justified 24,36 .
Experimentally, the detailed microstructural analysis of flash sintered samples of 3YSZ by M'Peko et al. indicates that there is a sudden increase of the concentration of oxygen vacancies near interfaces 37 . Additional Transmission Electron Microscopy studies on YSZ bicrystals by Yoshida et al. suggested that the increase in crystallographic misorientation facilitates yttrium segregation as a result of local disorder imposed by the interfacial energy 38,39 . Further, atomistic studies by Lee et al. highlight the segregation of oxygen vacancies and yttrium defects at interfaces followed by the development of depletion regions of defects in the immediate neighborhood of grain boundaries 40,41 .
Chen and Khachaturyan pioneered diffuse interface models to study the electrostatic and polarization effects in ionic solids 42 . Bishop and coworkers developed a thermodynamically consistent phase-field model for spinodal decomposition of charged domains in ionic ceramics 43 , and García et al. developed a generalized variational formulation which naturally includes Maxwell's Equations to describe the thermodynamics and kinetics of electrically active systems 44 . Separately, Guyer and coworkers proposed a phase-field theory to describe the local equilibrium of a double layer associated to an electrochemical interface 45 . Recently, Vikrant et al. introduced a variational formulation, which includes electrochemical and chemomechanical effects to describe the equilibrium space charge layers at a grain boundary and their associated transport properties 46 . Most recently, Vikrant and García applied this thermodynamic framework to understand the effect of misorientation and chemical dopant on the interfacial structural and electrochemical phase transitions and their overall impact on ionic transport properties 47 .
In spite of the rapid progress to harness the potential of flash sintering, a fundamental understanding of the mechanisms leading to the flash event remain unavailable. Further, the microstructural mechanisms controlling the incubation time and the critical power density necessary to flash are not known. In this paper, the very first thermodynamically consistent theory to explain the mechanisms that lead to the flash event is presented by considering N chemical species, f½V Zj j g ¼ f½V Z1 1 ; ; ½V ZN N g, of an ionic ceramic green body for three distinct microstructural regions, Here, each single crystal particle in a two-dimensional polycrystalline ionic ceramic green body is described using a crystallographic orientation order parameter, θ. This description sets the stage to rationalize the concurrent microstructural evolution and charge transport on a ceramic green body during electric-field assisted sintering. For Yttria Stabilized Zirconia and following recent work on charged interfaces 46,47 , the sum of all the contributions from chemical, electrical, and interfacial energy to the total free energy functional is given by where f is the Helmholtz free energy per unit volume, α × is the gradient energy coefficient of order-disorder interface, α d is the gradient energy coefficient of disorder-porous interface, α v is the gradient energy coefficient of porous-order interface, ρ is the electrostatic charge per unit volume, ϕ is local electrostatic potential, gðη Þ ¼ η 2 is a coupling function, and s 1 and s 2 are structural coupling parameters, as presented by Kobayashi et al. 48 .
Locally, for small deviations away from equilibrium, the variational derivatives associated to Eq. (1) constitute the structural and electrochemical driving forces for microstructural evolution 44,49 . The resultant kinetic equations are Here, the sintering phase is able to undergo an η × → η d transition, and defines a non-conserved order-disorder phase transformation 48,50 . Similarly, the crystallographic orientation of every granular particle is a non-conserved quantity 48,50 . In contrast, the porous phase, η v , and the concentration of solute and point defects in the condensed phase are conserved quantities 43,44,49 .
The structural discontinuity at the pore-condensed phase interface favors the accumulation of charged defects and enables the possibility of a charged surface. The concentration of solute and point defects in the condensed phase are conserved, in agreement with several authors 43,44,49 . In the absence of free surfaces or porosity, and in the limit of a single component system, Eq. (2) reduces to the work by Kobayashi et al. 48 .
For an externally applied electric field, the spatial inhomogeneity in conductivity of the green body due to particle-particle contacts and pore surfaces 47 induces local voltage gradients, which constitute electrostatic deviations, δϕ, away from the equilibrium potential. The total local electrical potential, ϕ T , is defined as ϕ T = ϕ + δϕ, is a solution of Coulomb's Equation (Eq. (7) in the "Methods" section). Here, the homogeneous solution, δϕ is a deviation away from the equilibrium electrostatic potential, ϕ. The solution to the homogeneous part of Coulomb's Equation is calculated using the charge continuity equation: , is the total electrical conductivity, and The temperature evolution of a green body due to ovenimposed heating is described by the homogenized heat transport equation 51 : where the right-hand side of Eq. (4) embodies the contributions from Joule heating, as the statistical ensemble of microstructural features generated by the ceramic green body during the period leading to the flash event, defined herein as the flash incubation period. Here, 〈κ〉 is average electrical conductivity, ∇ 〈δϕ〉 is externally applied macroscopic voltage drop across a width of the sample, σ e eAðT 4 f À hTi 4 Þ for radiation losses, and ∂hη i ∂t accounts for the coarse grained rate of change of order-disorder phase transitions. In the absence of bulk defects and order-disorder phase transitions, Eq. (4) reduces to existing thermal runaway models for flash sintering, e.g., refs [23][24][25] . See Supplementary Numerical Implementation section for numerical details and Supplementary Table 1 for material properties. Figure 2 summarizes the effect of Eqs. (1) through (4) on particle contacts, free surfaces, and triple junctions of a statistically representative green body on the charge transport and Joule heating subjected to a constant heating rate of 10 K/min and an applied electric field of 60 V/cm. Here, when the furnace reaches a temperature of T f = 1000K (see Fig. 2a), the local particle wetting is accompanied by interfacial structural disorder, which in turn induces the formation of a spatial distribution of contact areas that enable charge to percolate across a random distribution of particles. In agreement with previous work 46 , the increase in temperature at the particle-particle contacts provides vibrational energy to the oxygen vacancies, kinetically enabling these defects to segregate to surfaces and interfaces, as a result of the driving force imposed by the high interfacial energy value. The resultant electrochemical state induces an interfacial attraction of ½Y 0 Zr due to its opposite charge polarity. This results in the development of a positive electrostatic interfacial potential as the sample is heated up towards the flash event. Calculations demonstrate that only those particle contacts and surfaces that spatially percolate with locally higher ionic conductivity will control the transport of ions. Figure 2b shows that for a set of stream lines seeded at representative starting points on the far left edge of the sample, the resultant charge transport paths favor contact areas with locally higher ionic conductivity. Further, note that in very tortuously connected regions the most efficient path for charge flux corresponds to surfaces and triple junctions.

RESULTS
As the furnace temperature increases, e.g., to T f = 1200K, the concentration of oxygen vacancies at the large-angle particle contact interfaces increase resulting in the spontaneous appearance of fast ionic conducting paths, which compete for efficient charge flow with the already established particle contacts. As a result, when a local ion transport path becomes shorter than an already established route, charge flux is locally rerouted, microstructurally improving the Joule heating process (see Fig. 2c).
When the furnace increases its temperature to T f = 1300K, the number of ionic conducting paths that favorably contribute to the Joule heating generation increases. Large sections of favorable transport path are dominated by charge flow along particle surfaces. The increase in the number of particle network paths as the temperature increases enables the granular microstructure to transition from being externally heated by the furnace to increasing its internal temperature through self-generated Joule heating. Here, particle-particle contacts and triple junctions generate a lot of Joule heating as a result of the high selfinduced electrical conductivity in the interfaces and surfaces, in 3YSZ, defining a random distribution of hot spots. Further, the localized Joule heating favors particle-particle wetting for low angle misorientation locations, while decreasing it for those particle contacts with high angle misorientations. Thus, heating of the sample induced by the imposed electrical current favors necking formation and particle contacts, much before the flash event takes place, in agreement with experiments 18 . Finally, the rapid heating of a green body induced by the current density can initiate a faster densification process, as experimentally observed 26,27 . Overall, the performed analysis demonstrates that the incubation of the flash event is a result of the accumulation of point defects at particle contacts and surfaces in the 3YSZ system, effectively controlling the flash event by heating up the sample through the internal percolating charge conducting paths, which in turn increases the density of efficient Joule heating generating centers.
The compounding effects from the different surfaces and particle-particle contacts, including particle size polydispersity, long-range particle interactions, curvature effects, and the effects of the applied macroscopic voltage, define microstructural deviations from ideality and add inhomogeneities to the spatial distribution of ½Y 0 Zr , and ½V ÁÁ O , and the corresponding charge flux, during the flash event. Figure 3 shows the spatial distribution of yttrium and oxygen vacancies for a region in a randomly distributed particle population possessing a lognormal size distribution, as one would expect in a commercial powder 52 , for a furnace temperature T f = 1200K. The internal particle contacts display charge accumulation and depletion regions, as described herein; however, contact areas and interfacial curvature develop as a result of the misorientation-controlled wetting of the particles. These developing interfaces induce local microstructural ½Y 0 Zr and ½V ÁÁ O inhomogeneities that directly impact the positiondependent ionic conductivity. Simulations demonstrate that small contact areas develop a higher concentration of yttrium and oxygen vacancies; however, corner effects (junctions) dominate the segregation of the interface and shift the ½Y 0 Zr to a preferential corner. In addition, the vicinity of contact areas with negative radius of curvature become more depleted. Finally, particles whose diameter is comparable to the width of the electrochemically active area become fully depleted of yttrium and oxygen vacancies, making their core ionically very insulating. Figure 3d directly compares the predicted and experimental ½Y 0 Zr spatial distribution in the vicinity of a particle-particle contact, for a tilt misorientation, Δθ = 30°, demonstrating a very good agreement. The predicted extent of depletion zone is~3.3 nm, in good agreement with experimental value,~3.1 nm. Observed differences are a result of the limitations to measure grain boundary core thickness, as spatially resolved by EDS. Further, free surfaces also display a nonuniform distribution of charge. The vicinity of surfaces close to triple junctions display a larger amount of ½Y 0 Zr and ½V ÁÁ O , as point charges tend to segregate to edges and corners. These surficial heterogeneities are also a function of curvature differences between the particles in contact. Thus, surfaces with smaller radius of curvature will prefer to segregate less charge. The triple junctions are structurally more disordered than surfaces and particle contacts and accumulate larger amounts of ½V ÁÁ O and ½Y 0 Zr in agreement with EDS results in Fig.  1e, f, and g. A direct comparison against the existing isolated internal pores demonstrates that ½Y 0 Zr segregates at surfaces and interfaces, and that ½V ÁÁ O correlates directly with its accumulation and depletion. Figure 4a shows the predicted macroscopic sample temperature, T s , and sintering time as a function of furnace temperature, T f , for an externally applied electric field of 60 V/cm. The sample temperature linearly increases during the first 3000 s as a result of the oven-induced radiation heating. Here, the initial particle contacts and surfaces of percolated particles are formed, as demonstrated herein. The charge flow that is established through the percolating paths contributes to increase the Joule heating and electric power density of the sample. At T f = T s = 1120K, a kinetic microstructural transition is induced and enables the internal heating of the sample to be self sustained (see Fig. 4b), radiating heat back to the furnace. In this regime, the sample temperature becomes increasingly higher than the furnace temperature. The sample undergoes a thermal and electrical runaway at T on = 1312K, which is defined herein as the onset of flash sintering. While the sample internally heats up at temperatures that are higher than the oven temperature, surfaces and interfaces become increasingly disordered, but interfacial melting is never favored energetically.
The corresponding average predicted current density, 〈I〉, as a function of furnace temperature is shown in Fig. 4c, and highlights that the total conductivity of the heated sample follows an Arrhenius-type response as a function of sample temperature, with an activation energy corresponding to the majority charge  Supplementary Fig. 1). Inset (b) shows ½V ÁÁ O distribution, (c) shows ½Y 0 Zr distribution. Inset (d) directly compares the EDS measurement (red) and phase field prediction (blue) of the ½Y 0 Zr spatial distribution in the vicinity of a particle-particle contact, for a tilt misorientation, Δθ = 30°. The segregation of ½Y 0 Zr at the triple junctions, particle contacts, and surfaces are in good agreement with EDS results in Fig. 1e-g. The particle-particle contacts and free surfaces accumulate ½V ÁÁ O due to thermally enabled chemical driving forces, which naturally attracts the ½Y 0 Zr due to its opposite charge polarity.
carrier, the oxygen vacancies energy of migration, E m = 1.35 eV. The predicted critical current density is~1.0 A/cm 2 in excellent agreement with experimental results 18 . The compounding effects demonstrated in Figs. 2-4 result in the predicted macroscopic power density response, as shown in Fig. 1a as a function of furnace temperature and time for 3YSZ for a constant heating rate of 10 K/min and an applied electric field of 60 V/cm. A direct comparison to experimental results as reported by Cologna et al. 5 shows an excellent agreement. For these conditions, and as a result of the underlying microstructural fields, the flash temperature is identified herein as T s~1 312K, a few degrees away,~1300K, as compared to existing coarse-grain thermal runaway models [23][24][25] . The predicted incubation time for the flash event is~4230s. The increase in applied electric field to 120 V/cm reduces the onset flash temperature to~1125K and incubation time to~3040s. The average current density shows that the total conductivity follows an Arrhenius-type behavior as a function of the sample temperature, as the sample builds up to the flash event. At the microstructural level, HAADF-STEM (see Fig.  1b-d) and EDS mapping (see Fig. 1e-g) of flash sintered 3YSZ samples subjected to an applied electric field of 60 V/cm, demonstrates the ½Y 0 Zr segregation at pore surfaces, interfaces, and triple junctions (see Supplementary Experimental Setup for more details). The ½Y 0 Zr segregation in polycrystalline 3YSZ is a result of positively charged interfaces due to excess ½V ÁÁ O . Also, ½Y 0 Zr segregation alters the ½V ÁÁ O distribution at interfaces and enhances the ionic conductivity along the interfaces could potentially lead to a flash event 47 . Therefore, Fig. 1a highlights that the critical power density necessary to reach the flash event corresponds to a microstructural self-induced state favored by the Joule heating advent of percolative ionically conductive particle networks. Thus, as the electric field increases, the power density necessary to incubate this state decreases hyperbolically.

DISCUSSION
A thermodynamically consistent phase-field theory was demonstrated to rationalize the underlying microstructural mechanisms leading to the flash event in great agreement with experimental results. The structural, chemical, and electrostatic interactions of point defects are directly correlated with the intrinsic surfaces and interfaces which control the incubation kinetics of electric fieldassisted sintering. Three regimes leading to the onset of flash were identified: (1) Radiation dominated regime, where the initial contacts and surfaces of particles are established and enable charge to start to flow as a result of the heat gained from the oven thermal radiation; (2) Transition regime, where wide depletion zones of ½V ÁÁ O and ½Y 0 Zr form in the neighborhood of large-angle contacts between particles, impeding the charge flow across particles where charge was previously favored, deflecting it to new surficial and small-angle interfaces. Here, the contribution to the Joule heating of the sample continuously increases and dominates to increase the temperature of the sample over external radiation heating. Finally, (3) Percolating regime, where ½V ÁÁ O increases drastically at the core of the particle contacts, surfaces, and triple junctions allows charge along flow in multiple percolative paths, transitioning from interfacially-driven Joule heating to volumetric-driven Joule heating, resulting in a flash event.
The developed theory mainly focuses on explaining the incubation kinetics alone, i.e., only the build-up process to the flash event. Important effects occurring after the incubation event has resulted in a flash event, such as time-dependent elastic and plastic flow kinetics, densification, and grain growth during and after the flash event are beyond the scope of this work. The current variational formulation sets the stage to understand the densification and grain growth kinetics during the flash event by introducing the rates of generation, reaction, and sink of oxygen vacancies and cation vacancies, which can be readily incorporated in densification models, as reported by several authors, e.g., refs [53][54][55] . Also, rigid body motion of powder particles are not considered but can be introduced in the existing model to describe the mass transport during densification and grain coarsening kinetics, as reported by several authors, e.g., refs 55,56 . Fig. 4 Macroscopic Joule heating response of a ceramic green body during flash sintering. Macroscopic thermal and electrical behavior of 3YSZ sample as described by microstructurally averaging the Joule heating power density response of the simulated ceramic green body (see Fig. 2). Inset (a) shows the predicted sample temperature, T s , as a function of furnace temperature, T f , and sintering time. An external electric field of 60 V/cm and a constant heating rate of furnace, 10 K/min is applied. Inset (b) shows the different heating rate contributions to heating the sample: radiation heating (red), Joule heating (blue), and enthalpy of melting (green). Inset (c) shows predicted current density, 〈I〉, as a function of furnace temperature. The current density follows an Arrhenius relationship with the sample temperature till incubation time. The abrupt increase in conductivity after the incubation time defines a thermal and electrical runaway, which in turn defines the onset of the flash event.
The current formulation is demonstrated on 3YSZ and can be readily extended to other chemistries, such as 8YSZ and ZnO by incorporating the correct thermodynamic and kinetic data. Recent work, e.g., ref. 47 , and well established first principles studies 57 , demonstrate that the grain boundary ionic transport properties for 3YSZ can be enhanced, while the attractive interactions between ½V ÁÁ O and ½Y 0 Zr drop the ionic mobility of oxygen vacancies in 8YSZ, and lead to a significant decrease in the local particle-particle interfacial ionic conductivity. Finally, the application of the theoretical approach proposed herein to other chemistries will enable a detailed understanding on the generality, limits, and applications of flash sintering of single crystals and powder specimens 28 , including different particle sizes, morphologies and textures of powder specimens.
Overall, the developed theory provides a fundamental basis to engineer the microstructural evolution and charge transport of ionic ceramic green bodies before and during electric-field assisted sintering for a wide variety of technologically advanced ceramics, such as solid oxide fuel cells, solid state-based lithiumion batteries, and high-temperature gas sensors.

Theoretical framework
Define the Helmholtz free energy density, f, of an ionic ceramic green body for three distinct microstructural regions: the solid crystalline lattice, η × , the structurally disordered, interfacial disordered regions, η d , and those regions describing the open and closed porosity, η v . Site conservation is imposed on the system through the expression, η × + η d + η v = 1. f is a function of N chemical species, f½V Z1 1 ; ; ½V ZN N g 50 . The thermochemical volumetric free energy based on equilibrium defect principles 58,59 is defined herein as The free energy of formation of the ith species is where f i corresponds to the formation energy of ith species in the crystalline lattice, f d i to the disorder phase, and f v i to the pore phase. pðη i Þ ¼ 1 2 η 2 i ð15 À 25η i þ 15η 2 i À 3η 3 i À 15ð1 À η i Þ P j≠i η 2 j Þ is an interpolation function 60 , hðη i ; η j Þ ¼ η 2 i η 2 j is a pair wise function to prevent spontaneous phase transformations, and w ij is the energy barrier height. In the present study, two major species are considered, i.e., oxygen vacancies, ½V ÁÁ O , and yttrium defects, ½Y 0 Zr . The concentration of cation vacancies and cation interstitials are between two to nine orders of magnitude smaller than that of ½V ÁÁ O and ½Y 0 Zr , based on the existing experimental and DFT studies [61][62][63] . Other impurities, precipitates, and defects can be incorporated into this model through well established thermodynamic principles, i.e., by editing Eq. (5), and through the introduction of phase-field descriptions that capture the formation of new phases 57,[64][65][66] . In the absence of electrical effects for a binary system, Eq. (5) reduces to neutral grain boundary complexions formulation by Tang et al. 58 . In the absence of electrostatic interactions in a single component system, Eq. (5) reduces to the multiphase field model, as pioneered by Steinbach et al. 67 .
In agreement with recent work 44,46,47 , the extended free energy density, f ∘ , is a function of the electrostatic potential, ϕ, the electrostatic charge per unit volume, ρ, the electric displacement field vector, D ! , and the electric field vector, E ! , which allows to formulate contributions from the electrostatic charge energy density, ρϕ, and the dipolar energy density, By performing the Legendre transform, on f ∘ , the electrochemical free energy density, f ec , is minimized 59,68 .
Physical constrains, such as E ! ¼ À∇ϕ (a solution to Faraday's Law under constant magnetic field), and D ! ¼ ϵ E ! ¼ Àϵ∇ϕ (a constitutive law for electrically active systems in the absence of pyroelectricity and ferroelectricity), are directly substituted in Eq. (6). The crystallographic orientation in each single-crystal particle in the ionic ceramic green body is described using the order parameter, θ. The interfacial energy from particle contacts due to local misorientation, | ∇ θ|, is incorporated by extending the KWC model 48 . Two terms contributing to particle contacts interfacial energy are: s 1 g(η × )| ∇ θ|, and s2 2 gðη Þj∇θj 2 . Here, gðη Þ ¼ η 2 is a coupling function, and s 1 and s 2 are structural parameters, as proposed by Warren et al 69 . The sum of the contributions from thermochemical, electrical, and interfacial energy to the total free energy functional is given by Eq. (1).
In agreement with recent work 46,47 , Eq. (1) defines the equilibrium of a polycrystalline ceramic with electrostatically active grain boundaries and free surfaces by minimizing the free energy functional with respect to the controlling variables, f½V Zj j g ¼ f½V Z1 1 ; ; ½V ZN N g; fη i g ¼ fη ; η d ; η v g; θ; ρ; ϕ. For ease in the description, except for solid ↔ liquid phase transitions, and crystallographic order ↔ disorder volumetric, conserved or non-conserved, solid-state phase transformations such as chemical phase separation are not included; however, they can be easily incorporated, e.g., see refs 49,58,70,71 . Additional physical constraint on the volumetric charge density to the local composition of the ith species is ρ ¼ P N i¼1 eZ i ½V Zi i 43-47 The variational derivatives of Eq. (1) after substituting The first three rows of Eq. (7) define the structural state of every granular volume element of material, including crystalline grains, interfacial particleparticle contact, and internal and external surfaces. ξ i in the fourth row of Eq. (7) corresponds to the structural and electrochemical potential 46,47 , and is a local thermodynamic driving force for mass and charge accumulation of the ith chemical species at particle contacts and surfaces due to local structural, chemical, and electrical inhomogeneities. ξ i reduces to electrochemical potential in the absence of local structural inhomogeneities and chemical potential gradients in the absence of electrical effects, in agreement with multiple authors [43][44][45][46][47] . Lastly, the fifth row of Eq. (7) represents Coulomb's equation in its differential form, in agreement with previous work [43][44][45][46][47] .
Locally, for small deviations away from equilibrium, the structural and electrochemical variational derivatives associated to Eq. (1) and Coulomb's Equation constitute the driving forces for microstructural evolution 44,49 . The resultant kinetic equations are summarized in Eq. (2).
The effect of heating imposed by the oven is described by considering the contributions to the enthalpy of the system, h ¼ u À σ $ Á ε $ , for a volume element of material, where σ $ is the mechanical stress tensor and ε $ is the elastic strain tensor. Following Boettinger et al. 49 , and in agreement with Eq. (1), the change in enthalpy density is expressed as For an isobaric, or constant stress system, Eq. (8) reduces to Here, dS $ c σ dT T , where c σ is the heat capacity per unit volume at constant stress of 1 atm, and ΔH ×→d is the latent heat of melting per unit volume. The generalized change in enthalpy for M order-disorder phases, P M i¼1 ∂H ∂η i dη i , reduce to −ΔH ×→d dη × , for the case of only one orderdisorder phase transition. By applying the enthalpy method 72 , heat balance is readily accounted for as Here, Stefan-Boltzmann radiation exchange between the sintered solid and the surrounding oven walls is described by, _ Q rad ¼ σ e eAðT 4 f À T 4 Þ, as a boundary condition at the solid-pore interface. A is the effective surface K.S.N. Vikrant et al. area between the sample and the oven walls, and e is emissivity of the sintered gray body 73 .
In the limit of a sample whose characteristic physical dimensions are much smaller than the representative thermal diffusion distance, The radiation heat transfer of the internal surfaces of the green body is set to zero since the local thermal gradients across the porous microstructure are effectively zero. Equation (10) results in the homogenized heat transport equation 51 , Eq. (4).

DATA AVAILABILITY
The simulation data from this study are available upon request.

CODE AVAILABILITY
The partial differential equations were implemented in FiPy 3.1. The theoretical data from this study is available upon request.