Dynamics of baryon ejection in magnetar giant flares: implications for radio afterglows, r-process nucleosynthesis, and fast radio bursts

We explore the impact of a magnetar giant flare (GF) on the neutron star (NS) crust, and the associated baryon mass ejection. We consider that sudden magnetic energy dissipation creates a thin high-pressure shell above a portion of the NS surface, which drives a relativistic shockwave into the crust, heating a fraction of these layers sufficiently to become unbound along directions unconfined by the magnetic field. We explore this process using spherically-symmetric relativistic hydrodynamical simulations. For an initial shell pressure $P_{\rm GF}$ we find the total unbound ejecta mass roughly obeys the relation $M_{\rm{ej}}\sim4-9\times10^{24}\:\rm{g}\:(P_{\rm GF}/10^{30}\:\rm{ergs}\:\rm{cm}^{-3})^{1.43}$. For $P_{\rm{GF}}\sim10^{30}-10^{31}\:\rm{ergs}\:\rm{cm}^{-3}$ corresponding to the dissipation of a magnetic field of strength $\sim10^{15.5}-10^{16}\:\rm{G}$, we find $M_{\rm{ej}}\sim10^{25}-10^{26}\:\rm{g}$ with asymptotic velocities $v_{\rm{ej}}/c\sim0.3-0.6$ compatible with the ejecta properties inferred from the afterglow of the December 2004 GF from SGR 1806-20. Because the flare excavates crustal material to a depth characterized by an electron fraction $Y_e\approx0.40-0.46$, and is ejected with high entropy and rapid expansion timescale, the conditions are met for heavy element $r$-process nucleosynthesis via the alpha-rich freeze-out mechanism. Given an energetic GF rate of roughly once per century in the Milky Way, we find that magnetar GFs could be an appreciable heavy $r$-process source that tracks star formation. We predict that GFs are accompanied by short $\sim$minutes long, luminous $\sim10^{39}\:\rm{ergs}\:\rm{s}^{-1}$ optical transients powered by $r$-process decay ("nova brevis"), akin to scaled-down kilonovae. Our findings also have implications for the synchrotron nebulae surrounding some repeating fast radio burst sources.


INTRODUCTION
Magnetars are neutron stars (NSs) with surface magnetic field strengths ≳ 10 14 G (e.g., Duncan & Thompson 1992;Usov 1992;Kouveliotou et al. 1998;Kaspi & Beloborodov 2017) and relatively slow rotation periods of several seconds or longer.Although the physical origin of their strong magnetic fields remains uncertain, Galactic magnetars are likely produced in at least tens of percent of core-collapse supernovae (e.g., Beniamini et al. 2019).
One of the hallmarks of magnetars is their transient outbursts, as manifested through a spectrum of hard X-ray and soft gamma-ray bursts (e.g., Göǧüş et al. 1999Göǧüş et al. , 2000)).Based on the total energy released, the bursts can be divided into three categories: short bursts (isotropic energies  ≲ 10 41 ergs), intermediates flares ( ∼ 10 41 − 10 43 ergs), and giant flares (GFs) ( ∼ 10 44 − 10 47 ergs).See, e.g., the reviews by Turolla et al. (2015); Kaspi & Beloborodov (2017).Thus far, only a handful of GFs have been observed in the Local Group (e.g., Mazets et al. 1979;Hurley et al. 1999Hurley et al. , 2005)), with some ★ E-mail: jakub.cehula@mff.cuni.czevidence for GFs outside of the Local Group (e.g., Tanvir et al. 2005;Ofek et al. 2006;Svinkin et al. 2021;Minaev & Pozanenko 2020;Trigg et al. 2023).Perhaps the most notable example is the GF that occurred in December 2004 from SGR 1806-20 (e.g., Palmer et al. 2005;Hurley et al. 2005), an event which released ∼ 2 − 4 × 10 46 ergs in hard X-rays and soft-gamma rays during the initial bright spike that lasted for ≈ 0.2 − 0.5 s, was followed by a decaying minutes-long tail modulated at the NS rotation period, and a ∼ hour-long hard X-ray afterglow that peaked ∼ 600 − 800 s after the initial spike (e.g., Mereghetti et al. 2005;Frederiks et al. 2007).The GF produced a radio afterglow observed with multiple facilities in the following weeks and months (e.g., Cameron et al. 2005).The light curve shows a steep decay (∝  −2.7 ) around 9 days after the initial spike (Gaensler et al. 2005), followed by a rebrightening around day 25 lasting for about a week (Gelfand et al. 2005).Motion of the radio centroid and its polarization properties point to an asymmetric predominantly one-sided outflow (Taylor et al. 2005).A dynamical model capable of explaining the radio afterglow was proposed by Gelfand et al. (2005); Granot et al. (2006), in which the early steep decline phase arises from baryon-rich ejecta colliding with a shell surrounding a pre-existing cavity, while the rebrightening occurs as the comoving ejecta and shell decelerate upon sweeping up more material in the ambient medium.Most of the ejecta energy resides in material expanding at mildly relativistic speeds, with an initial velocity ∼ 0.7, kinetic energy ∼ 10 44.5 − 10 46 ergs, and total mass ∼ 10 24.5 − 10 26 g, for an assumed distance of 15 kpc (Granot et al. 2006).This scenario is also broadly consistent with the observed hard X-ray afterglow (Mereghetti et al. 2005).Additional modeling of the optical and radio afterglows generated as magnetar flares interact with their environments was presented in (Margalit et al. 2020;Wei et al. 2023), with Wei et al. (2023) emphasizing the potential to detect extragalactic flares with current and future radio telescopes.
An alternative interpretation of the observations is that GFs are relativistic coronal mass ejections, in rough analogy with solar flares (e.g., Lyutikov 2006Lyutikov , 2015;;Mehta et al. 2021).In this picture, GFs are purely magnetospheric events that produce strongly relativistic, strongly magnetized, and baryon-poor ejecta.Lyutikov (2006) constrains the mass ejected during the December 2004 GF to be ≲ 10 22 g.This model can reproduce the late-time behaviour of the radio light curve, but (at least thus far) provides no explanation for the radio rebrightening (Mehta et al. 2021).
In the baryon-ejection scenario of Gelfand et al. (2005); Granot et al. (2006), the sudden energy release responsible for the GF results from the gradual build-up of magnetic stresses in the NS crust, released in a sudden starquake/crustquake (e.g., Thompson & Duncan 1995, 2001;Perna & Pons 2011;Lander et al. 2015).Note, however, that the dichotomy between the solar flare and the starquake/crustquake paradigm is not clearly defined (e.g., Sharma et al. 2023).Recently, Demidov & Lyubarsky (2023), building on earlier work by Thompson & Duncan (1995), considered that a portion of the energy released during the initial bright spike, is trapped in a radiatively cooling electron-positron pair fireball in the NS magnetosphere.In such a strong magnetic field, the vacuum becomes birefringent and the photons are split into ordinary and extraordinary modes.The ordinary mode interacts strongly with matter while the extraordinary mode interacts only weakly.This allows for material to be ablated from the NS surface by irradiation, with Demidov & Lyubarsky (2023) finding that ∼ 10 18 g can be ejected this way over ∼ 100 s.However, this mass-loss phase, which occurs during the minutes-long tail/fireball stage of the GF, is not sufficient to explain the mildly relativistic ejecta of mass ≳ 10 24.5 g inferred by Granot et al. (2006), which instead is more likely to have been ejected during the initial stages of the flare when the magnetosphere is still open.We propose a model for such a prompt ejecta phase in the present paper.
The ejection of baryonic material during magnetar GFs would have consequences for a number of topics beyond just the phenomenology of Galactic magnetars.In particular, NS crust material is neutronrich, such that the decompression of GF ejecta into space may give rise to the conditions necessary for the creation of heavy elements via the rapid neutron-capture process (-process).The origin of the process is a long-standing mystery in nuclear astrophysics (Burbidge et al. 1957;Cameron 1957;Lattimer & Schramm 1974).Additionally, magnetars are considered the most likely central engines for fast radio bursts (FRBs;Lorimer et al. 2007;Lyubarsky 2014;Metzger et al. 2017;Kumar et al. 2017;Beloborodov 2017;Bochenek et al. 2020;CHIME/FRB Collaboration et al. 2020).These models require large-quantities of baryon-rich transrelativistic ejecta to explain the high rotation measures inferred from the synchrotron nebulae that surround particularly active FRB sources (Margalit & Metzger 2018;Zhao & Wang 2021;Sridhar & Metzger 2022), such as FRB 121102 (e.g., Michilli et al. 2018) and FRB 190520B (e.g., Niu et al. 2022).

Physical Picture for Baryon Ejection in Magnetar Giant Flares
All of the above motivates us to consider the hypothesis that the sudden energy release above the NS surface in a GF, drives a strong relativistic shock into the outer layers of the star that unbinds baryonic material from the NS crust.
We illustrate this process schematically in Fig. 1.We envision that one side of the magnetar magnetosphere undergoes reconnection or strong magnetic field dissipation (panel 1).The sudden release of energy that accompanies the beginning of the GF produces a high-pressure region above the NS surface comprised of photons and electron/positron pairs (red region in panel 2).We denote this pressure as  GF .The geometric thickness of this region is of order the initial magnetospheric energy density scaleheight.We denote the radial thickness of the higher-pressure region as Δ GF throughout this paper and we envision that it is of order or less than the NS radius (Δ GF ≲  NS ).
The high pressure of the GF shell  GF drives two shockwaves.The first is outwards into the low-density region outside the high-pressure shell.The second shockwave is into the NS atmosphere (denoted by the teal region in panel 2).The outwards-propagating shock meets no resistance and expands at ultra-relativistic speeds (outward arrows with  sh ≃  in panel 3).Such an ultra-relativistic shock, as it collides and shocks the gaseous environment surrounding the magnetar on much larger scales ≳ 10 12 cm, may produce a coherent radio burst akin to an FRB through the synchrotron maser mechanism (e.g., Lyubarsky 2014;Beloborodov 2017;Metzger et al. 2019).The inwards-propagating shockwave heats the NS crust while simultaneously encountering the very steep density profile of the NS.The sudden release of pressure from the outer edge of the high-pressure shell is communicated to smaller radii by a rarefaction wave, which moves inwards and eventually reduces the pressure being applied to the star on the characteristic timescale ∼ Δ GF / s ≲ 10 over which the shell radius appreciably expands, where  s ≃ / √ 3 is the sound speed of the relativistically-hot plasma.
Meanwhile, the inwards-propagating shock moves deeper into the NS crust.The region behind it is separated from the high-pressure shell by a contact discontinuity (teal-red transition in panel 3).As the shock moves inwards, it encounters regions of the crust with higher and higher pressure, slowing the shock down, and weakening it until its speed eventually drops below that of the rarefaction   = / √ 3. The shockwave effectively terminates where the initial pressure of the high-pressure region ( GF in panel 2) becomes approximately equal to the pressure at some depth within the NS crust  cr ≃  GF .The shocked NS crust, which is originally pushed inwards by the inwardpropagating shock, has a much higher specific energy than before the GF eruption, and now starts to expand outwards, transforming its thermal energy back into kinetic energy (akin to a spring uncoiling or a thermal blastwave).
During these phases, at least a portion of the magnetosphere becomes open because of the magnetosphere dissipation that initiated the GF.The open magnetosphere allows this super-heated material to escape (panel 4).As the flare subsides and the pressure is relieved, decompression of the shock-heated surface layers reconverts a fraction of the deposited thermal energy into bulk kinetic energy, allowing matter to escape (see Sec. 4.1).
The magnitude of the pressure in the GF determines the depth to which the NS is shocked, as follows from the approximate equality  GF ≃  cr .While the outer crust may be composed of lighter nuclei like iron with electron fraction 0.464,   decreases with depth in the NS crust (see Sec. 2).As we show below, and as can be confirmed (1) A sudden release of magnetic energy in the magnetosphere results in (2) the creation of a thin high-pressure layer ( GF ; red shell) above a portion of the NS surface (right hand side in the diagram); (3) The high pressure drives two strong shocks, one into the surrounding low-density magnetosphere propagating at  sh ≃  and a second into the crust of the NS (dark blue).Eventually, the inward-propagating shock stalls (typically meters to hundreds of meters below the original NS surface), and the pressure drop associated with the outward-expanding shock wave is communicated at the sound speed  s to the inner edge of the high-pressure shell; (4) The shocked neutron-rich crustal layers (teal shell), photodissociated into free nucleons by the shock, re-expands into space.A fraction of this material now possesses enough energy to reach the escape speed and become unbound.These layers possess sufficiently high entropy for an -rich freeze-out during seed nucleus formation, enabling a large enough neutron-to-seed ratio for the synthesis of heavy -process elements in the escaping ejecta (Sec.5.2) and powering a short-lived kilonova-like transient (Sec.5.3).
analytically, for large enough  GF , the shock-heated material should dissociate into free nucleons, leading to the ejection of high-entropy and relatively low-  material, again depending on the initial value of  GF .As it expands, the free nucleons will recombine, first to  particles and, via subsequent  and neutron captures to heavier elements via the -process (Sec.5.2; as in other contexts; e.g., Qian & Woosley 1996).An optical transient is expected from the expanding ejecta, as in models of kilonovae (Sec.5.3).
The amount of material unbound  ej from the NS crust is easily estimated at order-of-magnitude using several different analytic criteria based on either the total energy of the GF or the total energy density of the shocked crustal material.For example, an absolute upper limit on  ej arises from the global criterion that the total energy of the GF exceeds the minimal kinetic energy required to unbind the ejecta from the gravitational well of the NS.A more stringent, but more realistic local upper limit on the ejecta mass arises by consid-ering only those layers of the shocked crust which separately achieve positive energy.Considering that the NS crust starts from rest, this criterion is effectively that the internal energy of the shocked material exceeds the gravitational binding energy.We expand on these estimates of the ejected mass in Sec.4.3 and compare with our numerical calculations in Sec.4.4.
The material will only escape if it is not radiative.For most of the parameter space we explore, by post-processing our numerical results, we find that neutrino losses are of negligible importance to the dynamics of the outflowing ejecta.For sufficiently large  GF , corresponding to extremely powerful GF not yet observed, neutrino cooling may become dynamically important (Sec.4).

This Paper
Although the physical processes that give rise to a GF and its aftermath are undoubtedly a complex multi-dimensional magnetohydrodynamic event, here we present a simplified model as a preliminary exploration.We perform a suite of 1D spherically-symmetric special-relativistic hydrodynamic simulations, with a simplified Γlaw equation of state and ignoring radiative losses.As described above and sketched in Fig. 1, we initialize our simulations by setting up a high-pressure shell of defined width above the NS surface, and investigate the outcome for different values of the pressure  GF and the high-pressure shell width Δ GF , under the admittedly strong assumption that the effects of magnetar-strength magnetic fields can be neglected after a significant fraction of the magnetic field has been dissipated.Let us stress again here that the unbound ejecta from the Dec. 2004 GF could not have been spherically symmetric because the observed X-ray emission would have been obscured during the first minutes following the flare due to large scattering optical depth (Granot et al. 2006).Thus, our line of sight had to be along a relatively baryon-poor and radiation-rich portion of the outflow.Multi-dimensional models to understand the ejecta asymmetry will be the subject of future work.
This paper is organised as follows.We start by building a simple analytical NS crust model in Sec. 2. We describe the setup of our simulations and diagnostics of numerical results in Sec. 3. We present our simulation results in Sec. 4, starting with an in-depth analysis of the Fiducial model, and then moving on to explore a wider suite of simulations covering a wide range of initial shell pressures and widths (∼ GF energies), comparing our results for the ejecta properties to the analytic estimates.In Sec. 5 we summarize our results and discuss their consequences for a variety of topics, including -process nucleosynthesis, kilonova-like transients resulting from heavy element ejection, and the environments surrounding repeating FRB sources.Finally, we show our analytic estimates for mass-weighted ejecta distributions in Appendix A and the comparison of ejecta distributions for different simulation setups in Appendix B.

CRUST MODEL
We begin by describing the initial pre-GF radial structure of the NS crust, which we assume to be cold (i.e., non-accreting; though see Chamel & Haensel 2008 for a more general discussion).The crustal material resides in its absolute ground state in nuclear equilibrium, i.e. "cold catalysed material" (though we note this may not always be a good approximation if multiple GFs occur in rapid succession, ejecting matter faster than the crust can re-establish −equilibrium; see Sec. 5.4 for further discussion).Below the lowdensity NS atmosphere, resides the "outer crust", which starts at densities  ≈ 10 4 g cm −3 and extends to a depth corresponding to the neutron drip line at  ND ≈ 4.3 × 10 11 g cm −3 (Rüster et al. 2006).Below the outer crust, the "inner crust" extends to nuclear density at  nuc = 2.8 × 10 14 g cm −3 (Chamel & Haensel 2008).The total baryonic mass of the outer crust for a typical 1.4 M ⊙ NS is of order 10 28 − 10 29 g (10 −5 -10 −4 M ⊙ ) (Pearson et al. 2011).Considering that the inferred baryonic ejecta from the December 2004 GF was ≲ 10 26 g (Granot et al. 2006), we conclude that GFs of even significantly greater energy than those observed thus far, will mainly affect the outer crust.The equation of state (EOS) of cold catalysed matter for  <  ND should therefore be sufficient to describe those layers of the star affected by the GF-driven shockwave dynamics.

Structure of the Pre-Flare Crust
The EOS of the outer crust of nonaccreting cold NSs was intensively studied by Rüster et al. (2006), who utilized the theory of Baym et al. (1971), experimentally measured atomic masses from Audi et al. (2003), and made different assumptions about the nuclear mass model at densities higher than those accessible to laboratory experiment.Rüster et al. (2006) show that the EOS for  <  ND is well-established, with only small variations in the pressure of a few percent of pressure at a given  for different theoretical nuclear mass models.Hence, we only show their results here for two different theoretical nuclear mass models, namely, the nonrelativistic Skyrme force model BSk8 of Samyn et al. (2004) and the relativistic mean field TMA model of Geng et al. (2005), see the top and the middle panel in Fig. 2. For our purposes we only need the crustal pressure as a function of the crustal density  cr ( cr ) (top panel) and the crustal electron fraction  e,cr ( cr ) (middle panel).
We approximate the pressure of the crust as a polytrope, where log 10  * / ergs cm −3 = 19, log 10  * / g cm −3 = 4,  * is the polytropic constant and Γ * = 1.43 is the approximate effective polytropic exponent (fit to the crust pressure-density relationship).
To approximate the electron fraction  e,cr ( cr ) we use the constant  e,Fe = 26/(26 + 30) ≈ 0.464 for log 10  cr / g cm −3 < 6.90 and the quadratic fit of the BSk8 model with the least square method for log 10  cr / g cm −3 > 6.90.The top panel of Fig. 2 shows that the polytrope (1) approximates the true pressure  cr ( cr ) reasonably well, while the middle panel shows the constant  e,Fe + the quadratic fit is a good approximation of  e,cr ( cr ) for log 10  cr / g cm −3 ∈ ⟨4, 12⟩.

Density and Mass Coordinates
Neglecting general-relativistic corrections, the equation of hydrostatic equilibrium reads where  is the radius, () the gravitational acceleration,  the gravitational constant and  NS is the total NS mass (we neglect the self-gravity of the crust).Equations ( 1), (2) then define the crust density profile where  surf =  *  Γ * surf ,  NS is the NS radius, and  surf = ( NS ),  surf = ( NS ) are the surface density and pressure, respectively.The density scale-height is defined by   ≡ − cr /(d cr /d), as shown by the solid line in the bottom panel of Fig. 2.
The solid lines indicate our approximations, the polytrope (green; equation 1), the constant  e,Fe + the quadratic fit of  e,cr ( cr ) (red), and the density scale-height   (purple).The purple dashed line shows the resolution of our fiducial model (see Table 1).The black dotted vertical line indicates the neutron drip line at  ND .
12 km as a second (top) horizontal axis in Fig. 2. The total mass of the outer crust is ∼ 10 29 g when integrated up to the neutron drip density  ND (vertical dotted line), consistent with the generalrelativistic calculations of Pearson et al. (2011).The bottom panel of Fig. 2 shows that the fiducial resolution of our simulations (Sec.3) is sufficient to resolve   down to mass-depth  cr ∼ 10 22 g but not the surface layers further out.We use these pressure, density, and mass profiles (equations 1,3,4) as initial conditions for our numerical experiments (see Sec. 3.1) and in our analytic estimates of the ejecta mass (Sec.4.3).

NUMERICAL APPROACH AND DIAGNOSTICS
We perform 1D spherically symmetric-relativistic hydrodynamics (RHD) simulations of GF shockwaves with the Pluto code (Mignone et al. 2007).We describe the Pluto code and the simulation setup in Sec.3.1 and introduce the quantities needed for diagnostics of the simulation data in Sec.3.2.

Simulation Setup
Pluto is a multiphysics code designed for the treatment of discontinuous astrophysical flows (Mignone et al. 2007).The discontinuities are dealt with utilizing different Godunov-type high-resolution shockcapturing schemes for integrating a system of conservation laws.A general system of conservation laws can be written in the following form (Mignone et al. 2007) where  is a vector of conservative quantities, T denotes a tensor of fluxes of the components of , and  denotes the source term.Pluto includes different physics modules.We use the RHD module.
We use an ideal gas constant-Γ EOS (Mignone et al. 2021), such that the specific enthalpy obeys, where  is the thermal pressure,  is the rest-mass density, and Γ (corresponding to the hot shocked gas) is different from the index Γ * which defined the initial cold crust (see below).We initiate the simulation by setting up a shell of uniform highpressure  GF and width Δ GF ≲  NS above the NS surface, with the inner boundary at  =  NS and the outer boundary at  =  NS +Δ GF .Insofar that GFs are powered by the large-scale reconfiguration and dissipation of the magnetic field, we motivate the chosen values of the pressure assuming it should roughly equal the pressure of the magnetar's pre-GF magnetic field, i.e.
where  is the surface magnetic field strength and need not be identical to the polar magnetic field strength inferred from vacuum dipole spindown.We do not self-consistently include the effects of the magnetic field itself in our 1D simulations; rather, we implicitly assume that following a strong magnetic dissipation event the dynamical effects of the magnetic field can be neglected to first order across at least a portion of the NS surface, with material thus able to freely escape along the torn open magnetic field lines (Thompson & Duncan 1995) after the magnetic field reconfiguration (Fig. 1; we discuss the limitations of this assumption in Sec.5.5).
The energy in the high-pressure GF shell is given by, where  GF,30 =  GF /(10 30 ergs cm −3 ), Δ GF,1 = Δ GF /(1 km), ΔΩ is the solid angle subtended by the shell (Fig. 1), and the final line assumes ΔΩ = 4 and  NS = 12 km.For these parameters, the energy scales  ∼ 10 44 − 10 47 ergs of Galactic magnetar GF correspond to  GF ∼ 10 25 − 10 28 ergs cm −3 for Δ GF = 1 km.We assume a NS of mass  NS = 1.4 M ⊙ .As our default assumption, we place the inner boundary of the simulation region at  in = 10 km, the NS surface at  NS = 12 km, and the outer boundary at  out = 20000 km.The inner boundary condition is reflective, while the outer boundary employs an outflow condition.The initial density profile inside the star ( <  NS ) is given by our assumed crustal model (equation 3), which for numerical purposes we attach onto a power-law profile () =  surf (/ NS ) −3 above the surface  >  NS (containing negligible mass relative to the crust).The initial pressure profile follows equations (1,3), for  <  NS , we assume roughly constant pressure in the high-pressure shell  =  GF , for  ∈ ⟨ NS ,  NS + Δ GF ⟩, and the same power-law as for the density profile in the outer region () =  surf [()/ surf ], for  >  NS +Δ GF .We set the minimum density and pressure in the initial profiles to be  floor ≃ 2×10 −24  ND ,  floor ≃ 1×10 −24  cr ( ND ) to avoid numerical instabilities.
We take Γ = 4/3 in our simulations (cf.Γ * = 1.43, equation 1), motivated by the fact that the pressure and energy density of the high-pressure shell and shocked crustal material is dominated by radiation and relativistic non-degenerate electron/positron pairs.By adopting a Γ−law EOS, we also neglect the energy released by the shock-dissociation of alpha particles into free nucleons and their eventual recombination in the outflow (Sec.3.2.2);this is a reasonable approximation because the energy released ≃ 7 MeV per nucleon is modest compared to the asymptotic specific kinetic energy of the unbound ejecta, which is typically comparable to the NS gravitational binding energy ≈ 150 MeV per nucleon.For our default simulation B15.0_1km we set  surf = 1 g cm −3 ,  GF = 3.98 × 10 28 ergs cm −3 (corresponding to  = 10 15 G), and Δ GF = 1 km.However, we also run a suite of simulations spanning a wide range of  GF ∼ 10 25 − 10 32 ergs cm −3 and Δ GF = 0.1 − 10 km that corresponds to  GF ∼ 10 43 − 10 51 ergs.
As our default assumption, we employ a uniform radial grid of  u = 1500 points between  in and  u−s = 13 km and stretched grid of  s = 10000 points up to the outer boundary at  out .Our default uniform resolution of 2 m near the NS surface allows us to resolve the density scale height interior to a mass coordinate of ∼ 10 22 g, as shown in the bottom panel of Fig. 2. The length of the -th stretched grid element is   s Δ u , where  s is the stretching ratio and Δ u is the length of a uniform grid element.By default, we use a piecewise total variation diminishing (TVD) linear reconstruction, 2nd-order TVD Runge-Kutta time stepping, CFL condition of 0.5, and a simple TVD Lax-Friedrichs Riemann solver (Mignone et al. 2021).However, as we shall discuss in Sec.4.2, we also explore the sensitivity of our results to the resolution, reconstruction technique (from linear to parabolic; Mignone 2014) and Riemann solver (TVD Lax-Friedrichs to two-shock HLLC; Mignone & Bodo 2006), finding moderate quantitative differences in key quantities such as ejecta mass and entropy distribution relative to those obtained for the default setup.Let us also mention that a test simulation performed without a high-pressure shell ( GF = 0) gives rise to no ejecta.Each simulation is run for approximately one light-crossing time of the domain, i.e. 67 ms, which we find is sufficient for convergence of the unbound ejecta properties.

Diagnostics and Analysis
To analyze the effects of neutrino cooling and nucleosynthesis, we must extract the fluid temperature  from our simulation data.Under the conditions of high temperatures  ≫    2 ≈ 0.5 MeV and high entropies which characterize the shock-heated NS crust material, electrons/positrons are relativistic and form a non-degenerate ideal gas.The total pressure is therefore comprised of baryons and radiation (photons, electron, positrons) in thermal equilibrium, for which: where  u is the atomic mass unit,  is the Boltzmann constant,  is the radiation constant, and  = 1 is the mean molecular weight of the free nucleons.The entropy per baryon in relativistic particles in units of  can be approximated as (Qian & Woosley 1996) where  MeV ≡ /MeV and  8 = /(10 8 g cm −3 ).Radiation pressure dominates over ion gas pressure for  b ≫ 1.
The importance of degeneracy effects on the populations of electrons/positrons can be assessed by the ratio  e ≡  e /(), where  e is the electron chemical potential parameter.Using the approximate expression (Qian & Woosley 1996;their Eq. 6), we shall find  e ≪ 1 (equivalently,  b ≫ 2.4) for the shocked layers ultimately unbound from the star.This illustrates that degeneracy effects on the EOS and pair-capture weak interaction rates can be neglected to first order.We now introduce several of the key diagnostic quantities associated with neutrino cooling (Sec.3.2.1),nucleosynthesis (Sec.3.2.2),and for estimating the unbound ejecta layers (Sec.3.2.3).

Neutrino Cooling and Leptonization
Although the photon optical depths are sufficiently high surrounding the star during the post-flare evolution that photons are essentially trapped in the fluid, the shocked layers are transparent to neutrinos.As we do not include radiative cooling in our simulations, we must therefore check that neutrino cooling has a negligible impact on the dynamics of the system.An approximate condition for this to be justified is that the expansion timescale of the shocked gas,  exp , be shorter than the neutrino cooling time-scale,   .We estimate the local expansion time-scale as where we use the absolute value || because the velocity can become negative close to the NS surface in our simulations.An alternative but closely related global dynamical time-scale,  dyn , related to the unbound ejecta layers, will be defined in Sec.3.2.3.For the neutrino cooling timescale, we estimate where   =  pp +  cc is the specific neutrino cooling rate.For the latter, we consider the two most relevant processes: pair-production,   ν ↔ e + e − , and charged-current electron/positron capture onto nucleons,  e n ↔ e − p, νe p ↔ e + n.The combined neutrino cooling rate per unit volume for pairproduction by  e νe ,   ν , and   ν is approximately given by (e.g., Thompson 2002) Likewise, the cooling rate associated with the charged-current processes is (e.g., Thompson et al. 2001) We also define the total energy radiated away in neutrinos: Except for the most powerful magnetar GF, we shall find   ≫  dyn is satisfied in the layers ultimately unbound from the star; this indicates that (1) neutrino losses can be neglected on the ejecta dynamics; (2) after being shock-heated, the entropy  b in the unbound layers remains largely conserved as they decompress and undergo nucleosynthesis.
Although the material ultimately ejected from the star originates from the neutron-rich crust (  < 0.5), weak interactions in the hot shocked ejecta can in principle also change   from these initial values.Neglecting neutrino absorption reactions, the electron fraction evolves as (e.g., Qian & Woosley 1996) where  e−p is the rate of e − p → n e and  e+n is the rate of e + n → p νe .
The electron fraction therefore evolves on a timescale (e.g., Beloborodov 2003), where the final line assumes the high-temperature non-degenerate limit,  e+n ≈  e−p ≈  ∝  5 .We shall find that    ≫  exp is always satisfied in the unbound ejecta layers; this indicates that weak interactions are slow and hence   of the ejecta layers will retain their original values from the NS crust (Fig. 2).

Nucleosynthesis
We first must address whether the shocked crustal material will be dissociated into free nucleons (protons and neutrons).The free nucleon mass fraction can be approximated by the following expression: (Woosley & Baron 1992) As we shall show, most of the crustal layers ultimately unbound from the star are shock heated to sufficiently high temperatures to be initially dissociated ( N ≃ 1).As these shocked layers cool and decompress adiabatically away from the star, heavy element nucleosynthesis can in principle occur.In particular, once the temperature drops below ∼ 0.5 − 1 MeV, free nucleons will recombine into alpha particles, which themselves can undergo a neutron-aided version of the triple alpha process to create 12 C and ultimately (after additional  captures) heavier "seed nuclei" (Woosley & Baron 1992).When the ejecta is neutron-rich (  < 0.5), any remaining free neutrons can then capture onto these seed nuclei, creating -process elements (e.g., Hoffman et al. 1997).Neutron captures proceed to create nuclei up to a maximum mass number determined by the ratio of neutrons to seed nuclei.
For the moderately neutron-rich composition   ≈ 0.4 − 0.5 of the unbound ejecta in our simulations, heavy -process production beyond the second or third peak is only possible by suppressing the formation of seed nuclei (so-called "alpha-rich freeze-out"), as occurs for low densities (high entropy) or rapid expansion timescale.In particular, whether the -process proceeds up to the third peak (nuclear mass number  ∼ 195) can be characterized by a single quantity (Hoffman et al. 1997;Thompson & ud-Doula 2018) where  ,0 =  exp (,   )/(1 s) is the expansion time (equation 12) measured at the radius   where seed nuclei begin to form (typically,  MeV ≈  ,MeV = 0.5).We define the -particle formation radius according to where we take the minimum because the temperature profile is not monotonic in radius.To approximate the expansion time  ,0 of the ejecta at a given radius , we assume homologous expansion, i.e.  ∝  −3 .Further, assuming the entropy per baryon  b (equation 10) in the ejecta remains constant during its expansion, we have  ∝  −1 .This allows us to estimate  ,0 as follows where  exp,0 =  exp /(1 s) and we have assumed   ≫  NS ,  ≫  NS .If  >  crit , then the third -process peak is achieved, while if instead  <  crit the -process terminates at lower  ≪ 195.
The general trend of the analytical expression agrees well with the numerical simulations and thus serves as a useful guideline for the third -process peak detection (Thompson & ud-Doula 2018).

Unbound Ejecta
In analogy to equation (4) we can define a time-dependent mass coordinate  (, > ) as the mass above a given radius , this time calculated directly from the simulation data, where again  (, >  out ) = 0.The time-dependent ejecta mass  ej, (, > ) above a given radius  can thus be defined by where is the total energy density, comprised of kinetic, internal, and gravitational potential components, respectively.Mass shells with  > 0 are considered to form part of the unbound ejecta.Insofar that our simulations are only run slightly longer than the light-crossing time over the domain, we are justified to neglect any mass leaving the outer boundary of the simulation domain.We define the minimum radius  ej,min () of the unbound ejecta according to Thus, the total ejecta mass at a given time  becomes We define a global dynamical timescale of the ejecta as where  0.5 () = (,  ej,0.5 ()) and  ej,0.5 () is the half-mass radius defined by This estimate is meant to approximate, in a single quantity, the characteristic expansion time experienced by the bulk of the unbound ejecta (as distinct from the local expansion timescale defined in equation 12).Finally, to estimate the electron fraction of the ejected material,  e (e.g., as enters equation 20) we use the mapping  e (, ) =  e,cr (( (, > ))) (see Fig. 2), where ( (, > )) is obtained by inverting equations ( 1) and ( 4), respectively, i.e.  =  cr ( cr ( cr )).

RESULTS
Table 1 summarizes our suite of simulations and several key properties of the ejecta, while Table 2 shows convergence tests performed for a single model.We begin in Sec.4.1 by describing in detail the Fiducial model with  = 10 16 G ( GF = 10 30.6 ergs cm −3 ) and Δ GF = 1km.In Sec.4.2 we perform convergence tests which assess the sensitivity of our results to things like the the numerical scheme and adopted resolution.Then in Sec.4.3 we derive analytic constraints on the unbound ejecta mass, following the physical picture outlined in Sec.1.1.In Sec.4.4 we describe the broader simulation results covering a wide range of shock strengths (different  GF ) and shell-thicknesses (Δ GF ) and describe the summary trends and compare them to the analytic estimates.

Fiducial Model
Figures 3 and 5 show early-time snapshots zoomed-in on the NS surface, of the radial profiles of several quantities relevant to the ejecta dynamics, nucleosynthesis, and neutrino emission, while Fig. 7 shows the time evolution of several ejecta properties during this phase.By "early-times" we mean shortly after the high-pressure shell has been initialized, when the NS material is being pushed inwards and shocked, then released due to decrease of pressure communicated from the outside by the rarefaction wave, and, subsequently, undergoes outwards re-expansion as the high thermal energy is transformed back into kinetic energy (analogy to an uncoiling spring).This process takes few × 10 s in our Fiducial simulation.Later-time snapshots of the same quantities on a large radial grid, are shown in Figs. 4 and 6, illustrating the expansion of the unbound ejecta away from the NS.Mass-weighted distributions of several key properties of the unbound ejecta are shown in Fig. 8.
The shell-pressure of the Fiducial model  GF = 4×10 30 ergscm −3 corresponds to a uniform temperature  GF ≈ 2.0 × 10 11 K in the radiation-dominated medium.The second snapshot in Fig. 3 ( = 10.0 s; orange lines) corresponds to the point of the greatest compression of the NS crust.The density profile at this point is compressed inwards to around 0.4 km below the initial NS radius, with  s = / √ 3 ≈ 0.58 in the shocked region.Considering that 10 s • 0.58c ≈ 1.7 km, the rarefaction wave launched inwards from the outer edge of the high-pressure region has by this time caught up with the shock being driven into the NS crust, reducing the pressure sufficiently to allow material to reverse its radial direction and to start rebounding outwards.
The third snapshot in Fig. 3 at  = 33.4s shows the expansion of the previously shocked material, as its thermal energy is converted back into the kinetic energy of outwards motion (the uncoiling of the spring; thin green lines).This behavior is perhaps most easily visible in the density profile, which now extends hundreds of meters beyond the initial NS surface, as does the half-ejecta mass radius  ej,0.5 (equation 29).The positive outwards velocity  > 0 of the shocked material is also apparent.The fact that  < 0 at small  <  NS by this point shows that no more unbound ejecta originating from the NS is expected.Indeed, the ejecta shell is on the verge of being sonically disconnected from the NS surface (see also Fig. 4).
The second snapshot in the late-time profiles shown in Fig. 4 corresponds to  = 0.67ms (orange lines).As can be seen most clearly in the density profile, the ultra-relativistic shock driven into the region outside the high-pressure shell is followed by a low-density cavity region and, subsequently, the shocked NS crust ejecta.The latter has become sonically disconnected from the NS star already at this point, with the total ejecta mass  ej nearly saturated at its final value ≈ 10 25.5 g (see also Table 1).By the third snapshot at  = 50 ms, the velocity profile decreases moving to smaller radii in the ejecta, before passing through zero becoming negative at  < few × 100 km.This position marks the location of the now-weakened shock, separating the subsonic shocked NS crust still tightly bound to the NS from the supersonically-expanding largely unbound ejecta layers.The shock is also clearly visible in the density, pressure, and temperature profiles.
Figs. 5 and 6 show three radial profile snapshots, this time of quantities relevant to the ejecta thermodynamics, nucleosynthesis and neutrino emission.The early-time snapshots in Fig. 5 follow those in Fig. 3.We first note that all of the shock-heated material (the two later snapshots; orange and green lines) that becomes unbound is dissociated into free nucleons ( N = 1).These layers are also non-degenerate ( e ≪ 1 at  ∼  ej,0.5 ), justifying our neglect of electron/positron degeneracy effects on the EOS and weak interaction rates.We further see that while the shock reaches down to layers in the NS crust where  e ≈ 0.30, the unbound ejecta has a higher electron fraction  e ≈ 0.43 because it originates further out.The hierarchy of timescales   / dyn ∼   e / dyn ∼ 10 2 at around  ej,0.5 for the second snapshot (orange lines) and   / dyn ∼   e / dyn ∼ 10 3 near  ej,0.5 for the third snapshot (green lines), have the implications that: (1) neutrino losses on the ejecta dynamics can be safely neglected (see also Fig. 7); (2)  e in the ejecta changes sufficiently slowly due to weak interactions that our assumption of using the original crustal  e is justified.Finally, the fact that  exp / dyn ∼ 1 around  ej,0.5 implies that our conclusions are insensitive to which definition of the expansion/dynamical timescale we use.
The late-time snapshots shown in Fig. 6 follow those in Fig. 4. Again, the unbound ejecta is seen to remain non-degenerate ( e ≪ 1) throughout the entire simulation, and that its electron fraction saturates at  e ≈ 0.44 (see also Table 1).By the third snapshot (green lines) we see that   / exp and   e / exp become quite small ∼ 10 −2 close to the NS surface, indicating that our neglect of neutrino losses and weak interactions are not self-consistent at small radii and late times.However, because these inner bound layers are by now casually disconnected from the unbound ejecta shell, our predicted ejecta properties are not sensitive to changes in the late-time behavior that would result if neutrino cooling or weak interactions had been self-consistently included.Finally, the entropy of the unbound ejecta is observed to achieve values  b ∼ few × 100, corresponding to an -process figure of merit parameter  ∼ 10 11 (see also Table 1); we return to the implications of this finding in Sec.5.2.Fig. 7 shows the time evolution of several key quantities over the first ∼0.7 ms of evolution.The unbound ejecta mass  ej increases rapidly initially, before changing course to decrease after ∼ 10 s, and then increases again at ∼ 400 s, ultimately saturating at a value Table 1.Simulation results.The superscript * denotes simulations artificially terminated before 67 ms.The subscript _two-sh indicates simulations in which the Riemann solver is changed from TVD Lax-Friedrichs to two-shock HLLC (Mignone et al. 2021).Note the short-hand notation log 10 ≡ lg and that quantity ej denotes a mass-weighted average of the unbound ejecta as defined by equation ( 27) as measured at the end of the simulation.

Model name
Model ≃ 3 − 4 × 10 25 g by  ∼ 600 s.The initial "overshooting" of  ej occurs due to energy exchange between mass shells (i.e., layers which temporarily acquire  > 0 from the shock but later again become bound).The later-time increase at 400 s is a result of a supersonic shock originating in the NS crust due to complex interplay between the original shock, the NS surface, and the inner reflecting boundary condition.Subsequent weaker shocks originating in the NS crust manifest in Fig. 4 as discrete "steps" in various profiles at few×10 km as shown with the orange line; however, these additional shocks do not significantly increase the unbound ejecta mass.The behaviour of  e,ej,min follows that of  ej as expected from the  e,cr ( cr ) profile of the NS crust.While the total radiated neutrino energy   increases sharply during the early shock phase (≲ 10s), its final value reaches only a few percent of  GF (achieved when  ej saturates).This again illustrates that neutrino losses are not important for the ejecta dynamics, at least for the Fiducial simulation (see Sec. 4.4).After ∼ 200 s, the alpha-formation radius   in the ejecta saturates around a characteristic value ≈ 30 km a couple NS radii above the surface.The radii of the unbound ejecta  ej,min ,  ej,0.5 , and  ej,0.1 decrease slightly during the initial NS compression phase at  ≲ 10 s, (see Fig. 3), before increasing monotonically during the subsequent ejecta expansion and escape.The radii  ej,min ,  ej,0.5 ,  ej,0.1 intersect   at Table 2. Convergence tests of the model B15.0_1km with  = 10 15 G ( GF = 10 28.6 ergs cm −3 ) and Δ GF = 1 km.The default simulation setup is described in Sec.3.1.The notation of the ejecta properties is the same as in Table 1.Explanations of the changes to the default setup: 1 we change the NS surface density  surf from 1 g to 10 2 g, 2 we change the scaling of the initial density (pressure) profile outside of the NS (high-pressure shell) from ∝  −3 to ∝  −2 , 3 we change the inner boundary  in from 10 km to 11 km, 4 we change the outer boundary  out from 20000 km to 10000 km and the simulation duration from 66.7 ms to 33.3 ms, 5 we change the values of minimum density  floor / ND (pressure  floor / cr ( ND )) in the initial profile from 2 × 10 −24 (1 × 10 −24 ) to 2 × 10 −22 (1 × 10 −22 ), 6 we change the Γ-constant from 4/3 to that of the pre-flare NS crust Γ = Γ * = 1.43 (Sec.2.1), 7 we switch from 2nd-order to 3rd-order TVD Runge-Kutta time stepping (Mignone et al. 2021), 8 we switch from TVD Lax-Friedrichs to two-shock HLLC Riemann solver (Mignone et al. 2021), 9 we switch from piecewise TVD linear to parabolic reconstruction scheme (Mignone et al. 2021). ≈ 210 s, 227 s, 243 s, respectively, causing the mass-weighted distribution of log 10  ,0 to peak sharply at −3.64.Fig. 8 shows mass-weighted distributions of quantities relevant to the ejecta dynamics and nucleosynthesis, at three late-time snapshots.A comparison of the final two snapshots (orange and green lines, respectively) shows the mass-weighted distributions have largely saturated, i.e. the unbound ejecta is freely-expanding material with fixed properties.This terminal electron fraction distribution samples values from the NS surface down to  e ≈ 0.435.The expansion time distribution, log 10  exp,0 , is narrowly concentrated at the simulation duration , i.e. the time since the onset of the GF.This is consistent with our assumption of free expansion when calculating  ,0 , because for  = const.we have  exp = const.(equation 12; see also Fig. 6).The unbound ejecta layers obey  >  crit , thereby satisfying the conditions for third-peak -process nucleosynthesis (Sec.5.2).The similarity between the distributions of  and Be implies that a very similar result for  ej is obtained whether one adopts an energy criterion ( > 0) or an alternative Bernoulli criterion (Be > 0).The qualitative shapes of the mass-weighted distributions are derived analytically in Appendix A.
Model B15.0_1km_0.5x_1.43shows that increasing Γ from 4/3 to Γ * = 1.43 (corresponding to the polytropic index of the cold crust) results in a slightly lower ejecta mass  ej , leaves  ej practically untouched, lowers  b,ej , and thus  ej , significantly by a factor of about 4 and 50, respectively.However, physically Γ ≃ 4/3 is satisfied to high accuracy in the radiation-dominated shocked ejecta anywhere that  b ≫ 1.Indeed, we have checked that Γ = 4/3 is satisfied to better than 0.05 per cent by post-processing model B15.0_1km.
Switching the Riemann solver from TVD Lax-Friedrichs to twoshock HLLC (B15.0_1km_0.5x_two-sh)increases  ej by a factor of 3, increases  ej by around 0.06, but leaves  b,ej and  ej almost unchanged.On the other hand, models run using the two-shock HLLC solver are not more sensitive to resolution changes (B15.0_1km_twosh)than for those using TVD Lax-Friedrichs.Switching from the linear to parabolic reconstruction scheme (B15.0_1km_0.5x_par)increases  ej by a factor of around 5 − 6, increases  ej by around 0.22, and increases  ej by around one order of magnitude.Also, the change to parabolic reconstruction makes the model more sensitive to resolution changes (B15.0_1km_2x_par).The massweighted distributions of the key ejecta properties for the four different models (B15.0_1km,B15.0_1km_two-sh, B15.0_1km_2x, B15.0_1km_2x_par) in Fig. B1 (Appendix B) show that changing the resolution (B15.0_1km_2x) or solver (B15.0_1km_two-sh)gives behavior qualitatively similar to the default setup (B15.0_1km),but changing the reconstruction (B15.0_1km_2x_par)results in some apparently aphysical results, such as a double-peaked velocity distribution, leading us to disfavor this scheme.In summary, the degree to which our numerical results for the ejecta properties are conserved can be conservatively estimated based on differences that result from different assumed solvers.As such, we conclude that  ej is reliable to a factor ≲ 3,  ej to ≲ 0.05, log 10  b,ej to few hundredths, and log 10  ej to few tenths in the case of the B15.0_1km model.As a final sanity check, we note that a test simulation performed without a high-pressure shell ( GF = 0) gives rise to zero ejecta, i.e.  ej = 0. Since we have not set up the atmosphere in strict hydrostatic equilibrium with the adopted Γ * , we do see a time-dependent hydrodynamical re-arrangement of the mass.The simulation runs until the NS atmosphere develops a region of very low density where the outer power-law profile is attached to the steep density gradient of the NS atmosphere.This happens after about ∼ 30 ms of evolution and is to be expected, because we use a value of Γ = 4/3 that is not consistent in the deep NS atmosphere.

Analytic Constraints
We now provide analytic estimates of the unbound ejecta mass,  ej , following the physical picture outlined in Sec.1.1.Consistent with the 1D nature of our simulations, we assume spherically symmetric outflows occuring from all sides of the neutron star surface (solid  angle ΔΩ = 4); however, in the more realistic case that the outflows are confined to only a fraction of the total surface (Fig. 1), these may be considered as estimates of the isotropic mass-loss rate and rescaled accordingly (see Sec. 5.5 for further discussion).
An absolute upper limit on  ej arises from the global criterion that the total energy of the GF exceed the minimal kinetic energy required to unbind the ejecta from the gravitational well of the NS.Equating  GF (equation 8) to  ej  2 esc /2, where  esc = √︁ 2  NS / NS is the escape speed, we obtain: where  NS,12 =  NS /(12 km) and  NS, This global criterion corresponds to the assumption of a homogeneous energy distribution, such that each mass shell possesses just enough energy to escape.A more stringent, but more realistic local upper limit on the ejecta mass arises by considering only those layers of the shocked crust which separately achieve positive energy.Considering that the NS crust starts effectively from rest, this criterion reads  +  ≳ 0 (equation 25).Relating this condition on the post-shock density  sh of the marginally-ejected layer,  sh ≈ 7 ej ≲  GF /(Γ − 1)/(  NS / NS ), to its original pressure (equation 1) and associated mass-coordinate (equation 4) in the crust, we obtain the local upper limit on the ejecta mass: where we have used  sh ≈ 7 ej , corresponding to the jump conditions for a strong shock for a radiation-dominated Γ = 4/3 plasma in the Newtonian limit (Thompson 1986;their Eq. 24).
Even  ej,max,loc represents only an upper limit on the ejecta mass, because the shock does not have infinite time to drive into the crust before the driving external pressure is relieved by the outwards expansion of the high-pressure shell, which thus limits the amount of material shock-heated to the positive energy.As the shock propagates to higher densities deeper inside the crust, it weakens and propagates slower, enabling the rarefaction wave launched by the outwards expansion of the shell to catch up to the shock, relieving the pressure and allowing material to escape to infinity.For a thin initial high-pressure region (Δ GF ≪  NS ) this occurs roughly when  sh ≲ / where  sh is the shock velocity in the rest frame of the unshocked NS crust, see also Fig. 1.Substituting this condition into the relativistic shock jump condition (Blandford & McKee 1976;their Eq. 8), where  sh is the post-shock pressure and  s the shock's Lorentz factor in the NS frame, yields  ej ≳  sh / 2 as the criterion for a layer to be ejected.Finally, assuming an approximate equality of pressure across the contact discontinuity which separates the shocked crust from the high-pressure shell, where  ∈ ⟨0.1, 1⟩ is some constant of order unity, we obtain a local  4) on a zoomed-out radial grid.The meaning of all quantities is the same as in Fig. 5.
lower limit on the ejecta mass: Taking the ratio of equations ( 31) and (34) yields for  NS,12 = 1,  NS,1.4 = 1, and  = 0.5.As we shall show in Sec.4.4,  ej,min,loc assuming  = 0.5 provides a good fit to the normalization and functional form of the ejecta masses  ej ( GF ) measured from our suite of numerical simulations (equation 27) for Δ GF = 2 km (see Fig. 10).
One of the key quantities relevant to nucleosynthesis in the ejecta is the entropy it acquires from the shock (Sec.5.2).The entropy of the shocked layer can be estimated as (equation 10)  b,sh ≃ 5.21 where  sh,8 = 7 cr,8 and  sh ≃ (12 GF /11) 1/4 are the density and temperature of the post-shock gas, respectively.In Fig. 9 we show  b,sh as a function of  cr obtained from  cr ( cr ) (Fig. 2) for different values of  GF .We see that for sufficiently powerful flares ( GF > 10 30 ergs cm −3 ) to unbind  cr ≳ 10 24.5 g according to even our conservative local criteria, entropies of  b,sh ≳ 100 can be achieved in these ejecta layers.Such entropy values can be sufficient for third-peak -process nucleosynthesis (Sec.5.2).
. Time evolution of several quantities for the Fiducial model, including: unbound ejecta mass  ej (defined as all layers with  > 0, equation 27); minimum electron fraction of the unbound ejecta  e,ej,min ; ratio of the total energy radiated in neutrinos   (based on post-processing, equation 16) to the initial energy of the high-pressure sure  GF (equation 8); and various critical radii defining the unbound ejecta layers,  ej,min (equation 26),  ej,0.5 (equation 29),  ej,0.1 , and an estimate of where -particle formation occurs   (equation 21).Here,  ej,0.1 is defined in analogy to  ej,0.5 , as the radius exterior to which 10% of the unbound ejecta mass resides.

Parameter Study
In this section we compare the analytic constraints derived in Sec.4.3 with simulation results shown in Table 1.
In Fig. 10 we show our results for the total ejecta mass  ej from four sets of simulations with Δ GF,1 ∈ {0.1, 1, 2, 10} and other model parameters according to Table 1.For comparison we show analytic estimates of the maximal mass based the global energy criterion  ej,max,gl (equation 30) and the local energy criterion  ej,max,loc (equation 31) together with the minimal mass  ej,min,loc (equation 34).Red and blue shaded regions correspond to regions of parameter space where the shock energy is too weak to dissociate the ejecta into free nucleons (i.e. N ≪ 1; equation 19), or too strong for neutrino cooling to negligibly influence the ejecta dynamics (i.e.  / dyn < 1; equations 13,28), respectively.The quantities  N and   / dyn are measured for the models with Δ Gf = 1 km at times close to the point of the greatest compression of the NS surface ( = 10.0 s) and at the half-mass radius  =  ej,0.5 .The measured values  N (B13.5_1km)∼ 4 × 10 −2 ,  N (B14.0_1km)≈ 1,   / dyn (B16.0_1km)∼ 50, and   / dyn (B16.5_1km)∼ 2 × 10 −1 inform the painted regions in Fig. 10.Thus, for  ≲ 10 13.5 G ( GF ≲ 10 25.5 ergs cm −3 ) and  ≳ 10 16.5 ( GF ≳ 10 31.5 ergs cm −3 ), our physical assumptions break down and the accuracy of the simulation results are questionable.
As discussed in Sec.4.3, the ejecta mass  ej,max,gl ∝  GF based on the global energy criterion (equation 30) for Δ GF = 1km is much less restrictive than the local energy constraint  ej,max,loc ∝  Γ * GF (equation 31).The dependence of our numerically calculated  ej (equation 27) on  GF more or less follows the trend predicted by the minimum ejecta mass  ej,min,loc ∝  Γ * GF estimate derived from the local energy criterion (equation 34).We return to the implications of our results for the December 2004 GF of SGR 1806-20 in Sec.5.1.

SUMMARY AND DISCUSSION
Using 1D RHD simulations and analytic estimates, we have investigated the hypothesis that the sudden dissipation of magnetic energy above the NS surface in a magnetar GF drives a relativistic shockwave into the NS crust, ultimately giving rise to the ejection of neutron-rich material (Fig. 1).Our results for the system evolution and ejecta properties are summarized in Figs.7-10.The videos produced from the Fiducial model are publicly available in Zenodo, on YouTube, and alongside the paper as supplementary material on the journal website.We now describe the implications of our results for GF radio afterglows, -process nucleosynthesis, and FRBs.We conclude by highlighting some caveats and limitations of our numerical calculations and considering the prospects for future directions.  .Shown for comparison with symbols are analytic estimates for the total unbound ejecta mass (equations 31, 34) as well as the ejecta masses  ej measured from our simulations for a range of initial shell-thicknesses, Δ GF .For models capable of reproducing the baryon ejecta mass inferred from the Dec. 2004 GF of SGR1806-20 (vertical gray-shaded region), we should expect ejecta entropies  b,sh ≳ 100 in the range necessary for -process nucleosynthesis (Sec.5.2).

Radio Afterglows of Magnetar GFs
Assuming a spherical or nearly spherical outflow covering a large fraction of the NS surface, the ejecta masses  ej ≳ 10 24.5 g and kinetic energies we find are consistent with those inferred from the radio afterglow of the Dec. 2004SGR 1806-20 GF (Granot et al. 2006;Gelfand et al. 2005) for an assumed initial shell pressure  GF ≳ 10 30 ergscm −3 (Fig. 10), corresponding to a dissipated magnetic field of strength  ≳ 10 15.5 G (equation 7).The latter is indeed similar to the dipole field strength  dip ≈ 2 × 10 15 G measured from the spin-down rate of SGR 1806-20 (Olausen & Kaspi 2014).However, matching the thermal energy contained in the highpressure shell (equation 8) to the ≳ 2 − 4 × 10 46 ergs of gamma-ray emission observed from the 2004 GF (Palmer et al. 2005;Hurley et al. 2005) would appear to require a very thin shell ≲ 0.01 km, possibly pointing to some tension with the model or the observations.For example, the true energy budget of the 2004 GF might be considerably higher than implied by the observed (isotropic) gammaray luminosity.We further discuss this tension and another possible resolution involving an aspherical outflow in Sec.5.5.

r-Process Nucleosynthesis
Our finding that magnetar GFs can give rise to the ejection of significant quantities of NS crust material has implications for the synthesis of -process elements.Expanding initially cold NS matter was suggested as an -process site by Lattimer & Schramm (1974); Lattimer et al. (1977).Here the ejecta is far from cold: the GFdriven shock heats up the NS crustal material sufficiently to dissociate the heavy nuclei into free nucleons ( N ≃ 1; e.g., Fig. 5) for  GF ≳ 10 25.5 ergs cm −3 ( ≳ 10 13.5 G; see Fig. 10).Nevertheless, the temperatures achieved are not sufficiently high for weak interactions, particularly positron captures on neutrons, to appreciably increase the electron fraction from the neutron-rich composition   ∼ 0.4 − 0.5 of the pre-flare crust.We furthermore find that the entropy attained by the ejecta is sufficiently high (Fig. 9), and its expansion rate through the seed-nuclei formation region sufficiently fast ( ≳  crit ; equation 20) for an alpha-rich freeze-out (Hoffman et al. 1997), allowing neutron captures to proceed up to and even beyond the third -process peak (e.g., Fig. 6).Our results therefore support magnetar GFs as a new astrophysical site for the -process.
How important are magnetar GF to Galactic chemical evolution as a whole?The average production rates in the Milky Way for all -process elements, as well as just those above the second ( > 130) and third ( > 195) -process peaks, are roughly given by  r,all ≈ 1.9 × 10 −6 M ⊙ yr −1 ,  r, >130 ≈ 2.5 × 10 −7 M ⊙ yr −1 , and  r, >195 ≈ 5 × 10 −8 M ⊙ yr −1 , respectively (Rosswog et al. 2017).The rate of Galactic magnetar GFs is poorly constrained observationally, but if we assume that SGR 1806-20-like flares of energy ≳ 10 46 ergs occur at a rate of Γ GF ∼ 10 −2 yr −1 (Beniamini et al. 2019), then the corresponding GF -process production rate is  r,GF ∼ 10 −9 M ⊙ yr −1 [Γ GF /(10 −2 yr −1 )], assuming  ej ∼ 10 26 g of -process ejecta per flare (Granot et al. 2006).While this corresponds to a negligible fraction of  r,all , it represents a few percent of  r,A>195 , i.e. of those elements particularly difficult to synthesize in neutrino-driven core-collapse supernovae (e.g., Thompson et al. 2001).The yield estimated based on Galactic magnetars could furthermore represent a lower limit if a rare sub-population of magnetars (e.g., particularly young or highly-magnetized sources) are even more active than Galactic magnetars, as hinted by the existence of very active repeating FRB sources (see Sec. 5.4 for further discussion).
Even a subdominant heavy -process production which occurs promptly after star formation (as satisfied by GF given the young ages of Galactic magnetars; ≲ 10 4 yr) could play an important role in the Galactic chemical evolution through enrichment of the most metal-poor stars in the Milky Way or its satellite dwarf galaxies (e.g., Ji et al. 2016).Several studies have suggested that another channel (in addition to the established site of NS binary mergers) with short delay time relative to star-formation is needed to explain the [Eu/Fe]-[Fe/H] observations (e.g., Hotokezaka et al. 2018;Côté et al. 2019;Siegel et al. 2019;Van der Swaelmen et al. 2023;Lian et al. 2023).Channels that might satisfy this requirement include "magnetorotational" supernovae (e.g., Winteler et al. 2012;Nishimura et al. 2015; see however Mösta et al. 2018), disk outflows from hyper-accreting black holes (e.g., MacFadyen & Woosley 1999;Siegel et al. 2019), and proto-magnetar neutrino-driven winds (Thompson & ud-Doula 2018).We suggest magnetar GFs be added to the list of potential prompt enrichment sites.

Kilonova-like optical transients from magnetar GFs
As in the case of NS binary mergers, the ejection of -process nuclei from magnetar GF flares could power a short-lived kilonova-like transient powered by radioactive decay (Li & Paczyński 1998;Metzger et al. 2010).The transient peaks on the diffusion timescale, where  is the opacity of -process nuclei (e.g., Tanaka et al. 2020).The peak luminosity can be roughly estimated from the radioactive The effective temperature of the emission,  eff ≃ ( pk /4 SB ( ej  pk ) 2 ) 1/4 ≈ 2 × 10 4 K, where  SB is the Stefan-Boltzmann constant, corresponds to optical/UV wavelength emission.At an 15 kpc assumed distance of SGR 1806-20,  pk corresponds to a peak apparent g-band AB magnitude  AB ≈ 7. Thus, we predict that substantial baryon ejection in Galactic GF would be accompanied by short ≲ minutes-long bright optical flare, a transient event we dub "nova brevis" (a brief/short nova).Robotic telescopes with very large instantaneous fields of view covering much of the entire night sky (e.g., Evryscope; Law et al. 2014) appear best equipped to detect such bright, but very rare flares.Alternatively, for extragalactic magnetar GFs, a rapidly slewing space-based optical/UV telescope, which triggers on the short gamma-ray spike of the GF, could potentially detect the nova brevis signal for events in the relatively nearby universe ≲ 10 Mpc.Several UV satellite missions of this type, which can slew on a timescale of minutes with a sufficiently large field of view to encompass the gamma-ray error region, are currently planned, including ULTRASAT (Sagiv et al. 2014), QUVIK (Werner et al. 2023), or UVEX (Kulkarni et al. 2021).

Fast Radio Burst Environments
Even prior to the discovery of an FRB from a flaring Galactic magnetar (Bochenek et al. 2020;CHIME/FRB Collaboration et al. 2020), extragalactic magnetars were considered the most promising models for the central engines of FRBs (e.g., Lyubarsky 2014;Metzger et al. 2017;Beloborodov 2017; see Lyubarsky 2021 for a review).Two of the best studied and most active repeating sources, FRB 121102 (Chatterjee et al. 2017) andFRB 190520B (Niu et al. 2022), are spatially coincident with luminous synchrotron point sources.The sub-parsec sizes and high luminosities of these sources are consistent with them being young (∼ decades-centuries old) magnetized nebulae filled with relativistic electrons, the latter accelerated at the termination shock of a quasi-steady trans-relativistic "wind" created by the accumulation of ejecta from magnetar flares (Metzger et al. 2017;Beloborodov 2017;Margalit & Metzger 2018).The extremely high and time-variable rotation measure of the bursts from FRB 121102 (Michilli et al. 2018) can then be explained by the bursts passing through the same magnetized turbulent nebula responsible for the persistent synchrotron emission, but only provided that the flare ejecta has an electron-ion composition (a nebula of electron/positron pairs would impart no net rotation measure).Margalit & Metzger (2018) showed that all the basic properties of the FRB 121102 persistent source (size, flux, self-absorption constraints) and the large but decreasing rotation measure (RM) of the bursts, can be explained by a magnetar injecting a transrelativistic outflow with a velocity  ∼ 0.5 and time-averaged mass flux  ∼ 10 19 − 10 21 g s −1 .Given that the rarest, most powerful FRBs from FRB 121102 (perhaps analogs to magnetar GF) occur roughly once per day on average (e.g., Nicholl et al. 2017), this  corresponds to a per-flare baryon ejection of ∼ 10 24 − 10 26 g, similar to that inferred for for SGR 1806-20 GF (Granot et al. 2006).This scenario predicts that repeating FRB sources may go "dark" for a period of hours to days after powerful bursts (those generated by GF-energy events) as a result of any subsequent FRBs undergoing absorption or induced Compton scattering by the baryon shell (Metzger et al. 2019).
Taken together, our results suggest that young active FRB sources like FRB 121102 may be generating a high mass-flux  of heavy process nuclei.Over the Δ ≳ 10 yr active lifetime of FRB 121102, this would correspond to an -process yield of  r ∼ Δ ∼ 10 −6 − 10 −4  ⊙ , coincidentally similar to the potential -process yield of the proto-neutron star wind phase (e.g., Thompson et al. 2001).However, the estimated birthrate of repeating FRB sources suggest that only a small fraction ≲ 1% of all magnetars are as active as FRB 121102 (e.g., Margalit et al. 2020).Within a scenario for -process production from magnetar GF, a key open question is the relative contributions from the hyper-active magnetars responsible for repeating FRB sources and comparatively inactive magnetars like those in our Galaxy.The possibility of a high rate of magnetar crust removal due to repeated GF also raises the question of whether the crustal composition will have time to maintain −equilibrium (as assumed in our crust profile; Fig. 2), or whether the layers being excavated are even more neutron-rich than we have assumed.Future work should explore the long-term evolution of the NS crust subject to rapid mass-loss, using nuclear reaction network calculations similar to those used to study the opposite process of mass accretion (e.g., Schatz et al. 1999).

Simplifications, Caveats, and Future Work
Let us now address some of the simplifications and idealizations employed in this work.Firstly, we have neglected radiation losses; while a good approximation for photons in the highly opticallythick shock and outflow, neutrino losses can become relevant at high temperatures and densities such as those achieved in the shocked crust.We have explored neutrino losses in post-processing, finding that they are only important for the flare dynamics for extremely high-pressures  GF ≳ 10 31.5 ergs cm −3 ( ≳ 10 16.5 G; blue shaded region in Fig. 10).
A more significant caveat is that we do not include magnetic fields in our simulations.Strong magnetic fields (magnetization  m > 1) qualitatively change the dynamics of the plasma as a result of flux freezing and the strong Lorentz force (Lyutikov 2022;Barkov et al. 2022).Although the energy source of magnetar flares is clearly the magnetic field, we are implicitly assuming that immediately following the flare, the efficiency of magnetic dissipation is high enough that the plasma covering at least a portion of the NS surface is comparatively weakly magnetized.If the magnetic field in the crust remains sufficiently high that  m ≳ 1, then the transverse component of the field may suppress the formation of strong shocks, resulting in weaker heating of the plasma than we have assumed.Future RMHD simulations are necessary to explore the effects of strong magnetic fields on our results.
Our simulations are only 1D and hence implicitly assume spherical symmetry, despite the fact that the radio observations of the December 2004 GF indicate one-sided outflow on large scales with a 2:1 axis ratio (Taylor et al. 2005).In the physical situation, the high-pressure region, unconfined by the magnetic field, may cover only a small fraction of the NS surface  open = ΔΩ/4 ≪ 1.While for a fixed shell-pressure  GF , our 1D analytic estimates predict the total ejecta mass to scale as  ej ∝  open  1.43 GF (equations 31,34; Fig. 10), if one instead fixes the flare energy  GF ∝  open  GF (equation 8), then  ej ∝  1.43  GF  −0.43 open .Thus, confining a flare of the same total energy to a smaller and smaller covering fraction of the NS surface (smaller  open ≪ 1), increases the total ejecta mass, potentially bringing the ejecta mass implied by the radio afterglow of the Dec. 2004 GF from SGR 1806-20 into better accord with the energy budget implied by the radiated gamma-rays (Sec.5.1).How realistic is the sudden appearance of a high-pressure region above a localized region of the NS surface?At least the "sudden" aspect of such a picture is supported by the rapid rise-time of the gamma-ray light curves ≲ 0.25 ms (e.g., Palmer et al. 2005), comparable to the light crossing time of the magnetosphere (e.g., Lyutikov 2006).The source of dissipated energy may be Alfvén waves launched by a substantial crustal failure (e.g., Perna & Pons 2011;Lander et al. 2015) or an otherwise source of large-scale reconnection of magnetic field (e.g., Parfrey et al. 2013;Lander 2016).
While a 1D calculation may be a reasonable approximation to the physical picture if  open is not too small, confining the dissipated energy source to a narrower and narrower patch of the NS surface (in essence, envisioning a highly-localized "magnetic bomb") must eventually violate a 1D picture.Indeed, a natural expansion of our approach would be to increase the number of spatial dimensions of our numerical simulations to two or three.The inclusion of additional physical effects−magnetic fields, neutrino cooling, radiation transport, and general relativity−would also add more realism to the set-up.Predictions for the ejecta nucleosynthesis could also be improved using nuclear reaction network calculations.We leave these important extensions to future work.

Figure 1 .
Figure 1.Schematic illustration of the envisioned process of baryonic mass ejection following a magnetar GF (see alsoSec.4.3).(1) A sudden release of magnetic energy in the magnetosphere results in (2) the creation of a thin high-pressure layer ( GF ; red shell) above a portion of the NS surface (right hand side in the diagram); (3) The high pressure drives two strong shocks, one into the surrounding low-density magnetosphere propagating at  sh ≃  and a second into the crust of the NS (dark blue).Eventually, the inward-propagating shock stalls (typically meters to hundreds of meters below the original NS surface), and the pressure drop associated with the outward-expanding shock wave is communicated at the sound speed  s to the inner edge of the high-pressure shell; (4) The shocked neutron-rich crustal layers (teal shell), photodissociated into free nucleons by the shock, re-expands into space.A fraction of this material now possesses enough energy to reach the escape speed and become unbound.These layers possess sufficiently high entropy for an -rich freeze-out during seed nucleus formation, enabling a large enough neutron-to-seed ratio for the synthesis of heavy -process elements in the escaping ejecta (Sec.5.2) and powering a short-lived kilonova-like transient (Sec.5.3).

Figure 3 .
Figure 3. Radial profiles for the Fiducial model of quantities relevant to the ejecta dynamics, shown at three different early-time snapshots in time relative to the initialization of the high-pressure shell as marked.These quantities include the density , pressure , temperature , velocity  (in comparison to the local escape velocity  esc and sound speed  s ), total energy density  (equation 25), and the mass coordinate  (>  ) (equation 23); the latter matches the analytic coordinate  cr (equation 4) almost perfectly.The radius exterior to which half the unbound ejecta mass (equation 29) is located are shown on each profile by solid dots.On top of the cumulative mass profiles  (>  ) we show the original NS crust profile and two estimates of the total unbound ejecta mass  ej (equation 27).

Figure 4 .
Figure 4. Radial profiles for the Fiducial model of quantities relevant to the ejecta dynamics, but now shown at three later snapshots on a zoomed-out radial grid.The meaning of all quantities is the same as in Fig. 3.

↓Figure 5 .
Figure 5. Radial profiles for the Fiducial model of quantities relevant to nucleosynthesis and neutrino emission, shown at the same three early-time snapshots shown in Fig. 3.These quantities include: the free nucleon mass fraction  N (equation 19); degeneracy parameter  e =   /( ) (equation 11); electron fraction  e ; the ratios of different key timescales   / dyn , where   =   ,  e , or  exp (equations 28,13,18,12, respectively); entropy per baryon  b (equation 10); and the -process figure-of-merit parameter  (equation 20).

↓Figure 6 .
Figure 6.Radial profiles for the Fiducial model of quantities relevant to nucleosynthesis and neutrino emission, shown at three different late-time snapshots (same as in Fig.4) on a zoomed-out radial grid.The meaning of all quantities is the same as in Fig.5.

MFigure 9 .
Figure9.Analytic estimate of entropy obtained behind the shock driven into the NS crust,  b,sh (equation 36), as a function of the density, or, equivalently, mass-depth  cr , of the crust, for different values of the externally applied shell-pressure  GF (different colored lines, as marked).Shown for comparison with symbols are analytic estimates for the total unbound ejecta mass (equations 31, 34) as well as the ejecta masses  ej measured from our simulations for a range of initial shell-thicknesses, Δ GF .For models capable of reproducing the baryon ejecta mass inferred from the Dec. 2004 GF of SGR1806-20 (vertical gray-shaded region), we should expect ejecta entropies  b,sh ≳ 100 in the range necessary for -process nucleosynthesis (Sec.5.2).

Figure B1 .
Figure A1.Analytically derived mass distributions for different unbound ejecta properties for three different values of the pressure of the shock-heated material  sh as labeled.We show the classical  , (equation A5) and the relativistic  ,,rel (equation A9) velocity distributions, the distribution of the expansion timescale   at the -formation radius  ,exp,  (equation A14), and the entropy distribution  , (equation A17).We also show the numerical  -distribution  , together with an approximate analytic solution  , , (equation A22).The distributions are normalized so that ∫ ∞ −∞  , (  )d = 1.
2NS in the crust and  cr () ≫  surf .We show  cr for  NS = 1.4 M ⊙ and  NS = Unbound ejecta mass  ej from our suite of simulations as a function of the initial shell-pressure  GF (equivalently, strength of the dissipated magnetic field ; top horizontal axis) for different initial shell thicknesses Δ GF as indicated.Shown for comparison are analytic estimates bounding the ejecta mass (see Sec. 4.3).The range of ejecta masses inferred from the radio afterglow of theDec.2004GF of SGR 1806-20 (Granot et al. 2006) is shown with gray shading (Sec.5.1).