A semi-analytic model of magnetized liner inertial fusion

Presented is a semi-analytic model of magnetized liner inertial fusion (MagLIF). This model accounts for several key aspects of MagLIF, including: (1) preheat of the fuel (optionally via laser absorption); (2) pulsed-power-driven liner implosion; (3) liner compressibility with an analytic equation of state, artificial viscosity, internal magnetic pressure, and ohmic heating; (4) adiabatic compression and heating of the fuel; (5) radiative losses and fuel opacity; (6) magnetic flux compression with Nernst thermoelectric losses; (7) magnetized electron and ion thermal conduction losses; (8) end losses; (9) enhanced losses due to prescribed dopant concentrations and contaminant mix; (10) deuterium-deuterium and deuterium-tritium primary fusion reactions for arbitrary deuterium to tritium fuel ratios; and (11) magnetized alpha-particle fuel heating. We show that this simplified model, with its transparent and accessible physics, can be used to reproduce the general 1D behavior presented throughout the original MagLIF paper [S. A. Slutz et al., Phys. Plasmas 17, 056303 (2010)]. We also discuss some important physics insights gained as a result of developing this model, such as the dependence of radiative loss rates on the radial fraction of the fuel that is preheated.

The MagLIF concept at Sandia uses the electromagnetic pulse supplied by the Z accelerator to radially implode an initially solid cylindrical metal tube (liner) filled with preheated and premagnetized fusion fuel (deuterium or deuterium-tritium). The implosion is a result of the fast z-pinch process, where a large gradient in the applied magnetic field pressure operates near the liner's outer surface. 5,45 One-and two-dimensional simulations of MagLIF using the LASNEX radiation magnetohydrodynamics code 46 predict that if sufficient liner integrity can be maintained throughout the implosion, then significant fusion yield (>100 kJ) can be attained on the Z accelerator when deuterium-tritium fuel is used and the accelerator's Marx generators are charged to 95 kV to obtain a peak drive current of about 27 MA. 1,5 To maintain liner integrity throughout the implosion, liners with thick walls have been proposed 1,2 and used experimentally. 3,4,7,8,13,14,16 Thick walls mitigate the deleterious effects of the Magneto-Rayleigh-Taylor (MRT) instability, [3][4][5]7,8,13,14,45,[47][48][49][50][51][52] which first develops on the liner's outer surface, and then works its way inward, toward the liner's inner/fuel-confining surface, throughout the implosion. The parameter that is typically used to describe the robustness of a liner to the MRT instability is the liner's initial aspect ratio, where r l0 is the liner's initial outer radius and δ l0 is the liner's initial wall thickness. Lower A r0 liners are more robust to the MRT instability, but their use results in slower implosion velocities; thus, there is a tradeoff between liner robustness and implosion efficiency. Two-dimensional LASNEX simulations predict that liners with A r0 < 10 should be robust enough to keep the MRT instability from overly disrupting the fusion burn at stagnation. 1 These preliminary LASNEX simulations predict a broad optimum in the fusion yield surrounding A r0 ≈ 6 (see Fig. 10 in Ref. 1). Experiments 7,8 on the Z accelerator with no initial axial magnetic field have found the implosion dynamics of A r0 = 6 Be liners to be in good agreement with the 2D LASNEX predictions found in the original MagLIF paper 1 (e.g., compare Fig. 2 in Ref. 8 with Fig. 9 in Ref. 1). By contrast, A r0 = 6 Be liner experiments on the Z accelerator that included an initial axial magnetic field of about 10 T have revealed the presence of a helical instability structure that is fundamentally 3D in nature, 13,14 and thus could not be captured by the 2D LASNEX simulations; regardless, the inclusion of the initial axial field seems only to have improved overall implosion stability. 14 Moreover, the first fully-integrated experimental tests of MagLIF have shown that a fuel implosion driven by an A r0 = 6 Be liner, combined with preheat and premagnetization, does indeed result in fusion relevant conditions on the Z accelerator. [15][16][17] Therefore, we assume that the qualitative 1D versus 2D behavior shown in Fig. 10 of Ref. 1 holds fairly true for A r0 6, and thus that MagLIF can be modeled reasonably well in 1D for A r0 6.
The leading hypothesis for what seeds the MRT development is the electrothermal instability (ETI). [10][11][12] Recent experiments have shown that ETI growth can be tamped by applying a thick dielectric coating to the outer surface of the metal liner, thereby delaying the onset of nonlinear MRT growth. 12 Thus, it is hopeful that the application of dielectric coatings may enable MagLIF's 2D and/or 3D ≈ 1D performance to be extended to liners with aspect ratios significantly greater than 6. Note that these assumptions do not account for possible azimuthal variations; the effects of azimuthal variations on fuel confinement are beyond the scope of this paper, but they are presently being investigated in other studies. 53 Presumably, the deleterious effects of azimuthal variations are also mitigated by the use of low A r0 liners, and perhaps by dielectric coatings as well.
Since MagLIF calls for low A r0 liners to mitigate the MRT instability, the resulting liner implosions are too slow to shock heat the fuel (i.e., the fuel heating due to compression is approximately adiabatic). Thus, to heat the fuel to fusion-relevant temperatures (>1 keV) while simultaneously lowering fuel convergence requirements, the MagLIF concept calls for preheating the fuel just prior to the implosion. The parameter typically used to describe fuel convergence is the convergence ratio, where r g0 and r g (t) are the initial and time-dependent radii of the fuel-liner interface, respectively. The convergence ratio is often evaluated at particular times of interest, such as the time of peak fusion power, C rp , and the time of minimum r g , i.e., the time of "bounce", C rb . In many cases, C rp and C rb have fairly similar values, but differences between these two parameters can indicate important phenomena such as fuel collapse or fuel ignition. To reduce C rp and C rb requirements to 30, the MagLIF concept calls for preheat temperatures of 100 eV. At the Z facility, fuel preheating has been accomplished using the Z beamlet laser (ZBL), 16,54 which was originally coupled to the Z pulsed-power accelerator for diagnostic purposes (e.g., radiography 55,56 ). In the first fullyintegrated MagLIF experiments, 16 ZBL provided about 2 kJ of 532-nm light in a 2-ns pulse. Recently, ZBL has been upgraded to deliver about 4 kJ in a 4-ns pulse, and plans are in place to further increase ZBL's delivered pulse energy to 6 kJ and beyond over the next several years.
To keep the fuel hot during the relatively slow implosion of a MagLIF liner, the concept requires premagnetization of the fuel using an axially-aligned B z field. About 10-50 T are to be supplied by external B z coils. This initial seed field is to be amplified by a factor of 10 2 -10 3 within the fuel-filled volume of the imploding liner via magnetic flux compression. This large axial field is required to mitigate energy loss from the fuel due to elec-tron and ion thermal conduction. Additionally, the axial field should enhance α-particle confinement and heating of the fuel, and thus increase the overall fusion yield.
In the first fully-integrated MagLIF experiments, 16 as well as in liner implosion dynamics experiments, 13,14 initial seed fields of up to 10 T have been supplied by a newly developed applied B on Z (ABZ) subsystem at the Z facility. 57 Analysis of the neutron diagnostics data collected during the first fully-integrated MagLIF experiments indicates that significant flux compression did indeed occur. 17,18 Also, the ABZ capabilities at the Z facility have recently been upgraded to 15 T, and plans are in place to further increase these fields to 30 T in the next couple of years. 57 To elucidate some of the key physics issues relevant to MagLIF, we have developed a new semi-analytic model of the concept. This model is formulated as a system of ordinary differential equations (ODEs) that are straightforward to solve, particularly with standard software tools such as MATLAB ® , IDL ® , Mathematica ® , etc. This model accounts for: (1) preheat of the fuel (optionally via laser absorption); (2) pulsed-power-driven liner implosion; (3) liner compressibility with an analytic equation of state, artificial viscosity, internal magnetic pressure, and ohmic heating; (4) adiabatic compression and heating of the fuel; (5) radiative losses and fuel opacity; (6) magnetic flux compression with Nernst thermoelectric losses; (7) magnetized electron and ion thermal conduction losses; (8) end losses; (9) enhanced losses due to prescribed dopant concentrations and contaminant mix; (10) deuterium-deuterium and deuterium-tritium primary fusion reactions for arbitrary deuterium to tritium fuel ratios; and (11) magnetized α-particle heating. This model has been implemented in a code called SAMM (Semi-Analytic MagLIF Model). Simulations using SAMM typically take about 30 seconds to run on a laptop using the ode23 solver in MATLAB ® . Using the parallel computing cluster at Sandia, parameter scans of about 2000 simulations can be completed in as little as 10 minutes.
Semi-analytic models have proven themselves useful in the past. 26,58 Their fast run times allow the system parameter space to be explored rapidly. Also, the physics included in semi-analytic models are often transparent and thus accessible by anyone, which is not always the case for more sophisticated codes (e.g., the LASNEX code 46 that the original MagLIF paper 1 was based on). Finally, one of the most important aspects of developing simplified models such as SAMM, is that it forces those involved to verify and understand the results generated by more sophisticated codes. For example, when differences are observed between a simple model and a more advanced code, can those differences be explained in terms of the simplifying assumptions made in the reduced model? This type of questioning often leads to new physics insights that would otherwise be overlooked. For example, during the development of SAMM, we found an important relationship between radiative loss rates and the initial radial fraction of the fuel that is preheated; this finding and others are discussed further in this paper.
The remainder of this paper is organized as follows. In Sec. II, we present our semi-analytic model of MagLIF. In Sec. III, we verify this model relative to the 1D results presented in the original MagLIF paper (i.e., Ref. 1). In Sec. IV, we summarize this work. Unless otherwise specified, all units are SI.

II. A SEMI-ANALYTIC MODEL OF MAGLIF
A. Overview of the model and the radial distribution of the driving azimuthal magnetic field An overview of the semi-analytic MagLIF model is provided in Fig. 1. The liner implosion, and thus the fuel implosion, is driven by the pressure associated with the azimuthal magnetic field, which is supplied by the pulsedpower driver (e.g., the Z accelerator). Because of the cylindrical symmetry, the azimuthal field in the vacuum region is where µ 0 = 4π × 10 −7 H/m is the permeability of free space and I l is the liner current. We assume that B θ is partially diffused into the liner wall, and that its distribution in the liner region can be described by where the constant power β is found by forcing the amplitude of B θl (r) to drop by one e-folding within one skin depth, δ skin , of the liner's outer surface; that is, we set where e = 2.71828 is the base of the natural logarithm, and solve for β, which gives The skin depth is given by where ρ e is the initial electrical resistivity of the liner material, and τ r is the rise time of the driving B θ pulse (e.g., ∼ 130 ns for the Z accelerator). Some parameters of interest are provided in Table I. The purpose of assuming a simple power-law dependence for B θl (r), rather than an exponential dependence, is that it allows us to integrate both the total magnetic flux and the total magnetic field energy analytically, which enables us to obtain simple analytic expressions for the inductance, circuit response, and ohmic dissipation due to B θl .  1. Schematic overview of the semi-analytic MagLIF model. There are three primary regions: a fuel region, a liner region, and a vacuum region. The system height is h; the thermally-insulating axial magnetic field, which is initially distributed uniformly over all regions, is Bz; the radius of the fuel-liner interface is rg; the liner's outer radius is r l ; the return current radius is rrc; the azimuthal magnetic field, which drives the cylindrical implosion, is B θ . Normalized profiles are shown for B θ in the vacuum region, B θv (magenta), and in the liner region, B θl (orange); their analytic expressions are given by Eqs. 3 and 4 in the text (B θ is assumed to be zero in the fuel region). The liner region is further divided into multiple concentric liner shells; this discretization is necessary to avoid overdriving the fuel. Within the fuel region, normalized profiles are shown for the gas pressure, pg (blue), the gas temperature, Tg (red), the gas density, ρg (green), and the radiation temperature, Tr (cyan). The pressure profile is flat throughout the fuel (i.e., we have made an isobaric assumption due to the subsonic nature of MagLIF implosions). The gas temperature and density profiles thus have an inverse dependence to one another; their analytic expressions are given by Eqs. 105-109 in the text. The radiation temperature is nearly constant across the fuel region. The fuel region is further divided into a hot spot region from r = 0 to r h and a cold dense shelf region from r h to rg. The gas temperature in the shelf region is equal to the radiation temperature (i.e., the fuel material and the radiation field are in thermodynamic equilibrium in the shelf region). The shelf region erodes away throughout the implosion, until r h = rg, due to thermal transport from the hot spot to the shelf. The shelf region is only present if the fuel is preheated from r = 0 to r < rg.  Fig. 3), Z0=0.18 Ω, L=8.34 nH, C=8.41 nH, and L0≈5 nH. The elements Lv(t) and L lc (t) are given by Eqs. 16 and 20 in the text and they represent the coupling of the circuit to the dynamic volume of the model illustrated in Fig. 1. to adequately represent the Z accelerator coupled to zpinch loads. 59 With this model, simulations are driven by an open-circuit voltage, ϕ oc (t). This voltage is twice the forward-going voltage at the vacuum-insulator stack on Z, and it can be constructed from the data of previous Z experiments using Eq. 4 in Ref. 59. An example ϕ oc (t) waveform for Z is shown in Fig. 3. This circuit model also requires R loss (t) to be specified. This is a resistance used to model shunt current losses on Z, and it can be derived from previous Z experiment data using the methods discussed in Refs. 59-61. For simplicity, we will take R loss (t) to be large and constant (>10 Ω) throughout this paper, which essentially eliminates current loss from the results presented herein.
To solve the circuit shown in Fig. 2, we begin by writing the voltage across the capacitor and thus,İ The liner current is and thusφ The voltage required to supply B θ flux to the model's volume illustrated in Fig. 1 is Since the axial current density in the liner wall goes to zero at r = r g , which is implied by our assumption for B θl (r), then by the integral form of Faraday's law, where the path integral curve c is around the model's volume illustrated in Fig. 1, E is the electric field vector, dl is an infinitesimal path element vector along curve c, and Φ θ is the total B θ flux in the model's volume illustrated in Fig. 1, we know that we can write ϕ rc also as where Φ θv and Φ θl are the total B θ fluxes in the vacuum and liner regions, respectively. For the vacuum region, we haveΦ where L v is the standard coaxial vacuum inductance 62 and thus,L For the liner region, we have and thus, It is convenient to define an effective inductance that can be used to describe the circuit response due to the B θ flux in the liner region. From Eq. 19, we define this effective inductance as and thus,L With these expressions, Eq. 19 can be rewritten aṡ Finally, combining Eqs. 12, 14, 15, and 22 and solving forİ l givesİ

C. Drive energetics and ohmic liner heating
It is important to realize that the definition for L lc above is not the same as the definition for standard inductance; the standard inductance is found from the total magnetic field energy in the liner region 62 which gives L l (t) = µ 0 h(2βr l + r l + r g )(r l − r g ) 4πr 2 l (β + 1)(2β + 1) , and thusL l (t) = − µ 0 h(r g + βr l )(ṙ g r l − r gṙl ) 2πr 3 l (β + 1)(2β + 1) .
Note that with this standard inductance,Φ θl = L lİl + L l I l . Differences exist between L lc and L l because the liner current is distributed radially across the liner region. The distributed current is a result of our assumed B θl (r), which we chose in an attempt to model the magnetic diffusion and ohmic dissipation known to occur in thickwalled liner implosions. 63 These two definitions for inductance provide us with a clean and simple way to express the following powers being delivered to the volume of Fig. 1: (a) the total electromagnetic power, P EM ; (b) the power going into just the azimuthal magnetic field, P B θ ; (c) the power going into the kinetic motion and compression of the fuel, liner, and axial magnetic field, P kin ; and (d) the ohmic heating power, P Ω .
The total electromagnetic power is the sum P EM = P B θ + P kin + P Ω ; it can also be expressed in terms of the drive circuit as From Eq. 24, the power going into just the B θ field is To obtain an expression for P kin , we consider the case of a perfectly conducting liner, where P Ω = 0, and L lc = L l = 0 because the liner current is a surface current at r = r l (i.e., β → ∞). 64 This gives which is an upper bound for P kin . From numerical tests with finite conductivity, we find that to within a few percent over most of the implosion Finally, the power associated with ohmic dissipation is found from Substituting Eqs. 27, 28, and 30 into Eq. 31 gives the approximate (and lower bound) solution Note that Eqs. 27, 31, and 32 are valid only if the pulsedpower generator can supply azimuthal flux to the liner wall faster than the flux can be dissipated ohmically, and thus δ skin < r l − r g . In typical MagLIF experiments, this condition is essentially always met due to the use of electrically-conductive, thick-walled liners, and because the pulsed liner current changes rapidly with time throughout almost the entire experiment. 65 In numerical tests, we have found that Eqs. 31 and 32 generate ohmic dissipation rates that agree well with those calculated by full radiation magnetohydrodynamics simulations. Accounting for P Ω is important because without it, the liner region remains cooler and more compressible than expected throughout the implosion (i.e., on a "lower adiabat"). Upon stagnation, an overly compressed liner region will drive the fuel to a pressure that is too high, and thus result in an overly optimistic fusion yield.

D. Liner dynamics and compression
To account for liner compressibility, we first divide the liner region into N ls 20 concentric thin liner shells and write an equation of motion for the interface between each shell (internal interfaces), as well as for the fuel-liner interface and for the liner-vacuum interface (external interfaces). To each internal interface, we assign a mass of m ls = m l /N ls , where m l is the total liner mass; to each external interface, we assign a mass of m ls /2. There are N li = N ls + 1 liner interfaces in total. The radial positions of the liner interfaces, r l,i , are distributed from r g ≡ r l,i=1 to r l ≡ r l,i=N li . The initial radial position of interface i + 1 is given by where ρ l0 is the initial mass density of the liner material. The center of mass radial position for liner shell s is given by Since there are N li interfaces, there are N li equations of motion. The equations of motion for the internal liner interfaces are given bÿ where p l,s is the effective pressure in liner shell s. Each p l,s is comprised of material pressure, p ml,s , azimuthal magnetic field pressure, p B θl ,s , axial magnetic field pressure, p B zl ,s , and pseudo-pressure due to artificial viscosity, q l,s . The equation of motion for the fuel-liner interface isr where p g is the gas pressure in the fuel region and pB zg is the pressure due to the average axial magnetic field in the fuel region. The equation of motion for the liner-vacuum interface is where p B θlv is the pressure of the azimuthal magnetic field evaluated at the liner-vacuum interface (r = r l ) and pB zv is the pressure due to the average axial magnetic field in the vacuum region. The magnetic field pressures are simply where C is the component (either θ or z) and R is the region (either g for the fuel/gas region, l for the liner region, or v for the vacuum region). We also note the following. Since we have made an isobaric assumption in the fuel, and since, for MagLIF, the gas pressure typically dominates the magnetic field pressure in the fuel, we simply use the average axial magnetic field in the fuel for pB zg , where Φ zg is the total axial magnetic flux in the fuel. By contrast, in the liner region, we assume B zl (r)/ρ l (r) = const., and thus for p B zl ,s , we use is the average axial magnetic field in the liner region, Φ zl is the total axial magnetic flux in the liner region, is the average mass density in the liner region, is the mass density of shell s, and is the volume of shell s. For p B θl , we evaluate our analytic expression for B θl (r) (see Eq. 4) at the center of mass of each liner shell, r ls (t) (see Eq. 34). For p Bzv , we use the average axial magnetic field in the vacuum, where Φ zv is the total axial magnetic flux in the vacuum region of Fig. 1. Finally, for p B θlv , we use the azimuthal magnetic field evaluated at r = r l , and thus the resulting azimuthal field pressure at r=r l is With the magnetic field pressures and liner interface equations of motion thus defined, we can revisit Eqs. 27-32. Specifically, we can now obtain a more exact expression for the kinetic power, where δ(p B θ ) i is defined as the difference in the azimuthal magnetic field pressures on either side of interface i. This expression for P kin can be used to replace the approximate expression given in Eq. 30, and then from Eq. 31, we can obtain a more exact expression for P Ω , thus replacing Eq. 32. However, in practice, both Eqs. 30 and 48 work fine, since typically, they are within a few percent of each other over most of the implosion. To model the internal material pressure of the liner shells, we use the expression The first term, p 0 (ρ l,s ), is the zero temperature "cold curve" for the liner material, which is a function of only the mass density, ρ l,s . To model the cold curves of various liner materials, we use the functional form of Birch-Murnaghan, 66-68 where ρ l0 is the liner material's mass density at zero temperature and zero pressure, and where A 1 , A 2 , γ 1 , and γ 2 are fitting parameters. For lithium (Li), beryllium (Be), and aluminum (Al), we find reasonably good fits to SESAME equation of state data 69 using the values provided in Table II. The second term on the right hand side of Eq. 49 is the ideal gas thermal contribution to the equation of state, where E l,s is the total thermal energy of the liner material in shell s. With finite liner temperature, E l,s can change due to adiabatic compression and/or expansion. Additionally, E l,s can increase because of artificial viscosity and ohmic heating. We have found that including artificial viscosity is useful because it dampens the spring-like reverberations that occur in the liner's wall thickness throughout the implosion. The artificial viscosity model that we use is similar to that discussed in Ref. 70. Our pseudo-viscous pressure for shell s is given by where a q is an arbitrary coefficient of O(1) that specifies the length scale of the artificial viscosity in multiples of the radial thickness of shell s, The overall fusion calculations are not very sensitive to the exact value chosen for a q . We find that a q in the range of 1-3 works well for our purposes. Note Eq. 51 states that the artificial viscosity contributes to the shell pressure only during times of compression. The rate of change of the internal thermal energy of liner shell s is thus given bẏ The first term on the right hand side of Eq. 53 is due to adiabatic compression/expansion, while the second term is due to artificial viscous heating during compression only. The quantities P r and P c represent the radiative and thermal conduction loss powers from the fuel region, respectively; they will be described in more detail below. For simplicity, the powers P r , P c , and P Ω are all assumed to be fully absorbed by the liner and evenly distributed throughout the liner, hence they are each divided by N ls in Eq. 53. The liner is permitted to cool through adiabatic expansion as well as blackbody radiation from the liner's outer surface. The blackbody power is given by where σ = 5.67 × 10 −8 W/(m 2 ·K 4 ) is the Stefan-Boltzmann constant. For simplicity, energy losses from the liner due to P BB are evenly distributed throughout the liner, hence P BB is also divided by N ls in Eq. 53. Note that in Eq. 55, we use the average liner temperature, which is given bȳ where k = 1.38 × 10 −23 J/K is the Boltzmann constant, E l is the total thermal energy of the liner, N l = N at (1 +Z l ) is the total number of particles in the liner (ions and electrons), N at = m l /(uA) is the total number of atomic nuclei in the liner, u = 1.66 × 10 −27 kg is the unified atomic mass unit, A is the atomic mass number of the liner material (see Table II), andZ l is the average ionization state of the liner. We approximate the average ionization state using 71 where Z nuc is the atomic number (nuclear charge) of the atoms comprising the liner material (see Table II). The average liner temperature in keV is simplȳ where q e = 1.6 × 10 −19 C (or J/eV) is the charge of an electron. For reasonable fusion calculations, we have found it necessary to discretize the liner region into concentric thin shells, as described above, because even a compressible cylindrical slab-like model for the liner region (i.e., no internal interfaces, just the fuel-liner and liner-vacuum interfaces connected by the equation of state described above) results in very optimistic fusion yields due to very optimistic fuel compression. This occurs because in the cylindrical slab-like model, too much mass is assigned to the fuel-liner interface. The only way to reduce the mass at the fuel-liner interface is to use a finer discretization. We find that solutions tend to converge when approximately 20 liner shells are used.

F. Fuel ionization, energetics, and adiabatic heating
Within the fuel region, we assume adiabatic compression and a radially constant pressure profile (an isobaric assumption). In reality, there will be small radially dependent pressure waves that reverberate within the fuel; however, the slow liner implosion, which on average is subsonic relative to the fuel's sound speed, 1 implies that our adiabatic approximation is well justified. This also allows us to assume equal ion and electron temperatures, which we do throughout this semi-analytic model.
To model the energetics in the fuel region, we assume that the fuel is an ideal gas, and thus the isobaric particle pressure is related to the total thermal energy of the fuel, E thg , by where is the volume of the gas. Furthermore, the sum of the thermal energy and the ionization energy gives the total energy of the fuel We have found that accounting for the ionization energy is important for MagLIF in cases where the available preheat energy is low (∼ 100 J). For example, about 50 J would be required to fully ionize the deuterium fuel used in recent experiments at the Z facility. 16 This is a substantial fraction of the 100-300 J of preheat energy that is thought to have coupled to the fuel in these experiments. 15 For the ionization energy of particle species s in the fuel, we have found good fits to tabulated data 72 using the analytic expression where Z nuc,s is the atomic number (nuclear charge) of atomic species s, andZ s is the average ionization state of species s. Here we again make use of the approximate formula 71Z s = min 20 T g,keV , Z nuc,s , whereT and whereT g,keV andT g are the mean temperature of the entire fuel region in units of keV and K, respectively; N i and N e are the total number of ions and electrons in the fuel, respectively; N s is the total number of ions of atomic species s in the fuel; andZ g is the average ionization state of the fuel. Finally, summing over all species gives The purpose of accounting for various particle species, s, is so that we can study the effects of prescribed levels of dopants and/or contaminants ("mix") in the fuel region. Also note that since ionization depends on temperature and temperature depends on ionization, this would require an iterative solution within our overall ODE system solve. However, since we are only roughly estimating the ionization, and since it changes relatively slowly, we avoid having to implement an iterative solution by holding the ionization fixed during a system solve and updating it only between successive system solves (which is typically every ∼ 100 ps).
With ionization established, we note the following simple relationships, which will be used throughout the remainder of this manuscript: where m g is the total mass of the fuel; A s is the atomic mass number of species s; m d = 3.34 × 10 −27 kg is the mass of a deuteron; m t = 5.01 × 10 −27 kg is the mass of a triton; N d and N t are the number of deuterons and tritons in the fuel, respectively;m i is the average mass of an ion in the fuel;ρ g is the average mass density in the fuel;n i andn e are the average ion and electron densities in the fuel, respectively; andn s is the average ion density in the fuel for species s. The dynamics of the total fuel energy are described bẏ where P pdV is the adiabatic heating rate, P ph is the fuel preheating rate, P α is the heating rate due to α-particle energy deposition, P r is the radiative cooling rate, P c is the cooling rate due to electron and ion thermal conduction, andĖ ends is the energy loss rate due to fuel escaping out of the top and/or bottom ends of the imploding cylinder. The adiabatic heating rate is given by G. Fuel preheating (optionally via laser absorption) The fuel preheating rate can be described simply using a square pulse where E ph , τ ph , and t ph are the preheating energy, pulse length, and turn-on time, respectively. The experimental MagLIF program at Sandia National Laboratories has been investigating the approach of preheating the fuel using the Z beamlet laser. 16 In general, MagLIF does not require that the preheating be done with a laser. One could imagine finding a way to preheat the fuel using some of the energy supplied by the pulsed-power driver, which may enable more total preheat energy, as well as more efficient preheating of the fuel. Nevertheless, should the preheating be done with a laser, we can account for the beam propagation and deposition analytically.
An analytic description of laser propagation and deposition is particularly useful for quickly evaluating the various MagLIF target designs presently being considered for use at the Z facility. For example, MagLIF targets (see Fig. 2 in Ref. 16) have consisted of a cylindrical laser entrance channel (LEC), filled will fuel, that resides axially between the bottom of the laser entrance window (LEW) and the top of the imploding region (z = h in Fig. 1). We define the axial length of this channel as ∆z LEC . Furthermore, a beam dump has been used, residing below the imploding region (below z = 0 in Fig. 1). This beam dump is simply a reservoir of fuel that is present in case the laser propagates completely through the imploding region. If the fuel reservoir were not present, the laser could hit and ablate high-Z material from various target and/or electrode surfaces, which could spray back into the fuel of the imploding region, providing a source of contaminant mix.
To account for the energy absorbed in each region (i.e., the LEC region, the imploding region, and the beam dump region), we use the inverse bremsstrahlung model outlined in Ref. 1. By assuming an ionization state for the fuel while it is interacting with the beam; for example, we use the fully-ionized condition given bȳ and by ignoring thermal conduction and hydrodynamic motion (i.e., we assume fast isochoric laser heating), the heating of the fuel by laser light can be described by, 73 where ε b (z, t) is the energy density in the fuel that is interacting with the beam, I b (z, t) is the beam intensity in the fuel, κ(z, t) is the absorption coefficient, ν ei is the electron-ion collision frequency, c = 3 × 10 8 m/s is the speed of light in vacuum, ω p is the plasma frequency in the fuel, and ω b is the laser frequency. The plasma frequency is given by 74 where m e = 9.1 × 10 −31 kg is the mass of an electron, and ǫ 0 = 1/(µ 0 c 2 ) is the permittivity of free space. The laser frequency is given by where λ b is the laser wavelength. For preheating the fuel with the Z beamlet laser 54 at the Z facility, λ b = 532 nm. The electron-ion collision frequency is given by 75 where we have accounted for multiple ions species, s; T b is the temperature of the fuel that is interacting with the beam; and ln Λ is the Coulomb logarithm (see Appendix A). Next, using we factor out the ε b (z, t) dependence from Eq. 82 and evaluate the remaining factors at t = t ph ; we define the result asκ This implies that the only ε(z, t) dependence in κ(z, t) is that shown explicitly by substituting Eq. 86 into Eq. 85 and then Eq. 85 into Eq. 82. However, there is also an implicit ε(z, t) dependence in the calculation of ln Λ in Eq. 85, which we neglect for simplicity. For the kT dependence in ln Λ, we simply use the average energy per particle required to meet our full ionization assumption. That is, we use Eq. 62 withZ s = Z nuc,s , and thus we use kT = 13.6 · qe · s Z ζ nuc,s in the ln Λ calculation of Eq. 85.
The purpose of factoring out ε b (z, t) in Eqs. 81-87 is to make it clear that, for this subroutine of calculating laser absorption and propagation, we are simply freezing all dynamic variables [other than ε b (z, t) and I b (z, t)] to their values at the time when the laser is first applied (i.e., our assumption of fast isochoric laser heating). With this done, Eq. 81 becomes which has an exact solution given by 1 where and where z LEW = h+∆z LEC is the axial location of the LEW (where the beam first interacts with the fuel), and z f (t) is the axial location of the laser heating front as the beam "bleaches" its way through the fuel. Note that we have written this solution for a downward propagating beam, i.e., in the −ẑ direction, and that this solution is applicable only for t ph ≤ t ≤ (t ph + τ ph ) and z f (t) ≤ z ≤ z LEW . The energy deposition rates in the LEC region, in the imploding/fusion region, and in the beam dump region, are given by respectively, where r b is the radius of the beam and I b (z LEW , t) = P in /(πr 2 b ). Here, P in is the laser power that makes it through the LEW and enters the fuel.
For fusion calculations, we only need the energy deposited in the imploding region from z = 0 to h, hence we only need P ph as given by Eq. 94 (or by Eq. 79 if ignoring laser absorption). We assume that the energy deposited in the imploding region is instantly distributed uniformly in both the axial and radial directions. In the axial direction, this is a reasonable assumption since thermal conduction is uninhibited in this direction (i.e., the applied magnetic field that undergoes flux compression, B zg , is axially aligned, and thus inhibits thermal transport only in the radial direction). Furthermore, fullyintegrated 2D radiation magnetohydrodynamics simulations have shown that the preheat energy redistributes axially on a timescale that is fast compared to the implosion time. 15 In the radial direction, full radiation magnetohydrodynamics simulations show that a pressure wave (a "blast wave") is generated by the rapid preheating of the fuel, which expands radially outward from approximately the beam radius until it hits and reflects off of the liner's inner surface. 15 Subsequent pressure waves then reverberate within the fuel, but settle to a reasonably isobaric state on a timescale that is short relative to the overall implosion time.
Note that even though the pressure (overall energy density) in the fuel is assumed to be distributed uniformly in both the axial and radial directions, the fuel temperature and density profiles are assumed to be uniform only in the axial direction. Moreover, the radial temperature and density profiles are significantly affected by the radial extent of the fuel that is initially preheated relative to the radial extent of the overall fuel region. This is discussed in more detail below.
For the radius of the beam, r b , we ignore laser-plasma interactions and use the beam radius at the axial midpoint of the imploding region. Furthermore, we assume that the beam follows the focusing and defocusing cones described by 76 where ∆z b is the distance from the beam's focal plane to the axial position of interest, f /# is the beam's "fnumber", bf is the diameter of the beam in the beam's plane of best focus, and c b = 428.6 m −1 is a fitting parameter. Note that the second term in Eq. 96 is significant only in regions close to the focal plane; it is a first-order corrective term for describing the "waist" of the beam. As an example of using Eq. 96, we consider the experiments of Ref. 16, where the f /10 Z beamlet laser was focused to a 250-µm spot size, roughly 3.5 mm above the LEW; thus, with ∆z b = 3.5 mm, Eq. 96 gives a spot size on the LEW of 450 µm, as quoted in Ref. 16. Additionally, the length of the LEC was about 2.1 mm, and the length of the imploding region was about 7.5 mm; thus, with ∆z b = (3.5 + 2.1 + 7.5/2) = 9.35 mm, Eq. 96 gives r b = 492 µm. For reference, r g0 was 2.325 mm. Note that we have not accounted for absorption or scattering due to the laser entrance window. For this, we rely on separate laser-only experiments that measured transmission through foils similar to those used for LEWs on MagLIF experiments; this was done for various spot sizes on the foil and will be presented in a future publication. 77,78 For now, we simply note that where T LEW is the transmission through the LEW and P laser is the laser power on the vacuum side of the LEW. We comment that modeling, simulating, and experimentally diagnosing laser transmission, propagation, and absorption in such a system is a very challenging problem, even for much more sophisticated simulation codes; these are presently very active areas of research within the MagLIF program and elsewhere.

H. Fuel heating via magnetized α-particle energy deposition
The next term to calculate in Eq. 77 is the heating rate due to α-particle energy deposition, P α . To do so, we note that the α particles are produced only by DT reactions and that the energy carried by each α particle is Following Ref. 33, the fraction of Q α that is deposited in our magnetized cylindrical fuel is calculated using where and where l α is the mean free path of an α particle, m α = 6.64 × 10 −27 is the mass of an α particle, v α0 = 2Q α /m α is the birth velocity of an α particle, Z α = 2 is the charge of an α particle, r αL is the Larmor radius of an α particle, and the ln Λ calculation is described in Appendix A. The α particle heating rate is then given by whereṄ dt is the primary DT reaction rate, which is given below (see Eq. 158).

I. Model of fuel hot spot and dense outer shelf regions
For the remaining terms in Eq. 77, which all describe energy losses from the fuel, and for describing magnetic flux loss due to the Nernst thermoelectric effect, we need to assume something about the temperature and density gradients in the fuel. For this semi-analytic model of MagLIF, we assume temperature and density profiles that are comprised of two parts: a low-density hot spot region and a cold dense shelf region (see Fig. 1). This separation into two regions occurs when the preheating is done rapidly across a region from r = 0 to r ph0 , where r ph0 < r g (t ph ). 79 This separation occurs because the rapid preheating causes a blast wave to propagate radially outward from r ph0 toward the outer radius of the overall fuel region, r g (t). This blast wave redistributes mass by boring out a low-density hot spot and pushing colder non-preheated fuel up against the liner's inner surface.
For the low density hot spot region, we assume a fuel (gas) temperature profile of the form and a radiation temperature given by where T c is the central (peak) temperature at r = 0, r h is the outer radius of the hot spot region (see Fig. 1), ξ is the power-law parameter that specifies the curvature of the profile, and T B is the brightness temperature, which is found iteratively while calculating the radiative cooling rate, P r (see discussion about P r below, in Sec. II J). By comparing to detailed calculations (see for example Fig. 5 in Ref. 1), we find the hot spot profile to be reasonably well modeled using ξ ≈ 2 with Nernst effects, and ξ ≈ 3.5 without Nernst effects. The results are not overly sensitive to this choice. We have found good agreement using anything from ξ = 2 to ξ = 8. Because of our isobaric assumption, specifying the temperature or the pressure at any point in the fuel automatically determines the other, i.e., ρ g (r) · T g (r) = const. Thus, Eq. 105 means that the density profile in the hot spot region is given by where ρ c =ρ gTg /T c is the density at r = 0. In the cold dense shelf region, the fuel temperature and the radiation temperature are given by and thus by our isobaric assumption, the fuel density in the shelf region is given by The temperature gradients in each region are thus given by for r h < r ≤ r g .
(110) Note that the temperature and density profiles are nearly flat in the shelf region and that T g (r), T r (r), ∇T g , ρ g (r), and ∇ρ g are all finite throughout the entire fuel region. The reasons for using these functional forms, as well as how to find T c , ρ c , r h , and T B , are deferred to the discussion below on radiative losses from the fuel (Sec. II J).
With these profiles, we can now define the following radially dependent number densities in the fuel: n e ≡ n e (r) =Z g n i (112) n s ≡ n s (r) = N s · ρ g (r)/m g , which are for the ions, electrons, and the ions of species s, respectively. Also, by assuming a "frozen-in" condition for the axial magnetic field in the fuel plasma, we have B zg (r)/ρ g (r) = const. =B zg /ρ g , and thus Note that the Nernst effect (discussed below) actually breaks the frozen-in condition (cf. Fig. 5 in Ref. 1). Nevertheless, this expression for B zg (r) provides us with a simple first-order approximation of the radial dependence of the axial magnetic field that is accurate enough for the purposes of this model.

J. Radiative losses
The next term to calculate in Eq. 77 is the radiative cooling rate, P r . This calculation determines the fuel temperature and density profiles described in Eqs. 105-109 by finding T B iteratively while applying isobaric and conservation of mass arguments. The formulation is derived from a two-temperature gray model of nonequilibrium radiative diffusion. For more details on this type of model, see the discussion on pages 261-262 in Ref. 71.
Using a two-temperature model was motivated after studying the results of full radiation magnetohydrodynamics simulations of MagLIF and observing the following: (1) the radiation temperature is nearly constant throughout the fuel, dropping only slightly in the shelf region as r approaches r g ; (2) the radiation temperature is significantly lower than the fuel temperature in the hot spot region; (3) the radiation temperature is roughly equal to the fuel temperature in the shelf region; (4) the radiative flux across the fuel-liner interface is well approximated by (1 − α s )σT 4 s , where α s 0.9 is the albedo of the liner's inner surface and T s (r) is the temperature in the shelf region (both the fuel and radiation temperatures since they are equal in the shelf region).
These phenomena occur because the liner's inner surface heats and compresses rapidly to a point where it can no longer absorb or transmit radiation efficiently. Thus, this surface begins to reflect/reemit radiation back into the fuel region, effectively trapping a significant portion of the radiation in the fuel region. Moreover, thermal diffusion from the liner's hot inner surface to the colder material deeper within the liner region is slow and on a timescale that is long compared to the implosion time. This slow diffusion is due in part to the high heat capacity of the liner material, and it helps to enable the radiation trapping.
These observations inspired the following twotemperature gray model for radiative losses over the volume of the fuel, where opacity effects become stronger in and near the cold dense shelf region. From Eq. 6.62 in Ref. 71, we haveε where ε g is the energy density of the fuel, κ r is an averaged opacity described by Eq. 6.60 in Ref. 71, and κ p is the Planck mean opacity. The first term on the right hand side of Eq. 115 describes the rate at which the fuel absorbs energy from the radiation field, while the second term describes the rate at which the fuel loses energy to the radiation field. Since the radiation is trapped in the fuel region by the liner's hot inner surface, and since, in the shelf region, the fuel mass and radiation fields are in thermal equilibrium, we assume that κ r ≈ κ p , and thuṡ Note that a consequence of Eq. 116 is that there is no net energy exchange between the fuel and the radiation field when the two are in thermodynamic equilibrium (T r = T g ). Next, from Eq. 5.22 in Ref. 73, we have where J is the frequency-integrated emission coefficient for the bremsstrahlung mechanism and A br = 1.57 × 10 −40 m 3 · K − 1 2 · J/s. 80 Substituting Eq. 117 into Eq. 116 and integrating over the volume of the fuel from r ′ = 0 to r ′ = r gives the cumulative radiative cooling power from the fuel as a function of r P rv (r) = A br · 2πh ·Z 2 g r 0 n i n e T g 1 − Note that since we have set T r = T g in the shelf region, the contribution of the shelf region to the P rv integral is zero; this is consistent with the fuel's mass and radiation field being in thermodynamic equilibrium in the shelf region. The first term in Eq. 118 represents a radiation source and the second term represents a radiation sink due to opacity and absorption. Note that for the case of T r = 0 and constant density and temperature profiles, integrating Eq. 118 over the entire fuel region gives the familiar bremsstrahlung radiation power from an optically-thin volume emitter: Next, to find T r , we invoke the additional constraint that, in the shelf region (r h ≤ r ≤ r g ), the volume radiation described by Eq. 118 must be consistent with the gray body surface radiation described by To do this, we first note that as long as P rv (r h ) = P rs (r h ), then P rv (r) = P rs (r) throughout the shelf region because of the functional form assumed for the shelf temperature in Eq. 108, i.e., T g (r) = T r (r) ∝ r −1/4 . This functional form was chosen to ensure that Eq. 120 would return a constant value, and thus enforce our assumption of no radiative loss or gain within the shelf region. We also note that, in the hot spot region (0 < r ≤ r h ), we have assumed T r = const. = T B , and thus T r (r h ) = T B . Therefore, to ensure consistency in the shelf region, we must find the values of r h and T B that satisfy where To find the values of r h and T B that satisfy Eqs. 121-122, we use a simple bisection method and evaluate the expressions for P rv (r h ) and P rs (r h ) iteratively until P rv (r h ) = P rs (r h ) to within 1%. 81 Also, before we can test our iterative guesses for r h and T B in Eqs. 121-122, we need to find ρ c and T c to complete the profile definitions given in Eqs. 105-109. To do this, we use our isobaric assumption and conservation of mass arguments. Additional computational details for finding r h , T B , ρ c , and T c are provided in Appendix B.
Throughout a given simulation, α s = 0.9 is held fixed. However, the results are not very sensitive to this choice. In practice, we find that any reasonable value, say in the range of 0.5-0.95, works fine.
Since, in the shelf region, P rv (r) = P rs (r) = const. = P rv (r h ) = P rv (r g ), the total radiative cooling rate for Eq. 77 can be taken from anywhere in the shelf region; we take P r = P rv (r g ). (123) Note that in this model, if the entire fuel region is preheated uniformly, then the cold dense shelf region will never exist (only a hot spot region will exist), and thus all of the fuel mass will contribute to the radiation losses. Also, due to the isobaric assumption and the conservation of fuel mass, the absence of a shelf region means that the central peak temperature, T c , will be lower for the same amount of energy deposited in the fuel. For these reasons, uniformly preheating the entire fuel is not optimal. Moreover, preheating only a central portion of the fuel, say from r = 0 to r ph0 ≈ 0.5 · r g (t ph ), broadens the optimal MagLIF operating space by enabling the use of higher initial fuel densities (i.e., if these higher initial fuel densities were used in cases where all of the fuel was preheated, then the radiation losses would cool the preheated fuel too rapidly to obtain good fusion yield). This is fortunate since the experimental MagLIF program at Sandia is presently investigating the approach of preheating the fuel via a laser with a small beam radius relative to r g (t ph ).
Finally, we comment that assuming this shelf region is present in experiments may be idealistic. Azimuthal or other 3D asymmetries (particularly during laser preheating) might make it impossible to establish or maintain this shelf region. If a shelf region cannot be maintained (which would be difficult to diagnose experimentally), then the optimal operating space for MagLIF might not be as broad as that suggested by simulations where the initial fuel density is scanned over and r ph0 < r g (t).

K. Magnetized electron and ion thermal conduction losses
The next term to calculate in Eq. 77 is the energy loss rate due to electron and ion thermal conduction. To do this, we use the Epperlein-Haines transport equations for the electrons 82 and the Braginskii transport equations for the ions. 83 The energy loss rate due to electron thermal conduction is given as a function of r bỹ where κ e = n e kT g τ ei m e · 6.18 + 4.66x e 1.93 + 2.31x e + 5.
is the coefficient for electron thermal conduction perpendicular to B zg , x e ≡ ω ce τ ei is the electron Hall parameter, ω ce = q e B zg /m e is the electron cyclotron frequency, and τ ei = 1/ν ei is the average time between electron-ion collisions. Here, ν ei is the electron-ion collision frequency given by 75 where we have accounted for multiple ion species, s, and where the ln Λ calculations are described in Appendix A. Similar to electron thermal conduction, the energy loss rate due to ion thermal conduction is given as a function of r byP where is the coefficient for ion thermal conduction perpendicular to B zg , x i ≡ ω ci τ ii is the ion Hall parameter, ω ci = q eZg B zg /m i is the ion cyclotron frequency, and τ ii = 1/ν ii is the average time between ion collisions.
Here, ν ii is the ion-ion collision frequency given by 75 where n h ≡ n d + n t ,Z h ≡Z d =Z t , and where we have assumed that the plasma is comprised mainly of fuel deuterons and tritons (subscript h for hydrogen isotopes), with an average dominant ion mass of that is much less than the mass m s of any of the dopant/contaminant particles. The ln Λ calculations are described in Appendix A. Finally, the total radially dependent energy loss rate due to thermal conduction is given bỹ P c (r) =P ce (r) +P ci (r).
For this calculation, we use radially dependent parameters n i (r), n e (r), x e (r), and x i (r) [and thus B zg (r), ν ei (r), ν ii (r), n h (r), n s (r), and ln Λ(r)], since they are all readily calculable from our given expressions for T g (r) and ρ g (r). While studying thermal conduction using full radiation magnetohydrodynamics simulations (as well as our hot spot model), we observed that ion thermal conduction dominates over electron thermal conduction in the inner regions of the fuel hot spot. By contrast, near the edge of the hot spot (i.e., r ≈ r h ), the thermal transport is dominated by electron conduction. Thus there is a handoff from ions to electrons at some radius in the fuel hot spot region. Because of this, and because our hot spot profile is prescribed and determined by the radiation model, we wish to be conservative with regards to the expected losses due to thermal conduction, and therefore we look for a maximum inP c (r). That is, for the thermal conduction loss power from the hot spot region to the shelf or liner region, we use For the thermal conduction loss power from the shelf region to the liner region, we simply use The overall thermal conduction loss power from the fuel to the liner depends on whether or not the shelf region exists, and thus, for the thermal conduction power in Eq. 77, we use P c = P cs for r h < r g P ch for r h = r g .
Note that even if a shelf region exists, we still need to calculate P ch in order to estimate the shelf erosion rate (see Eq. 150 below).

L. End losses
The final term to calculate in Eq. 77 is the energy loss rate due to fuel losses out of the ends of the imploding cylinder,Ė ends . We estimate this loss by assuming that the fuel (with its associated mass and energy densities) flows across the top (z = h) and bottom (z = 0) planes of the imploding region at the hydrodynamic sound speed, which is determined by the pressure and mass density of the fuel. Therefore, the energy losses out of the top and bottom planes, respectively, are given bẏ where r top (t) and r bot (t) are the radii for the top and bottom apertures that the fuel can escape through, and c g (r) is the hydrodynamic sound speed, given by where γ g = 5/3 is the ratio of specific heats for an ideal gas. The factor of (3/4) 4 is a result of the derivation presented in Ref.
1. Finally, the total energy loss due to end losses is 84Ė Similarly, for mass losses from the fuel, we assume ion density flows across the apertures, and thus, for ion species s, we havė The number of deuterons and tritons in the fuel change because of fusion reactions as well as end losses (see Eqs. 162 and 163 below). By contrast, end losses are the only way for prescribed dopant and/or contaminant particle numbers to change. Thus, we havė where N mix (t) = s =d,t N s (t) is the total number of dopant and contaminant particles. Note that we only need one differential equation to describe the evolution of all of the "mix" particles because their number distribution amongst themselves remains constant in time, i.e., and therefore, Finally, should fuel preheating be done with a laser entering from above the imploding region, then these end loss calculations should be made for at least the top plane, where the radius of the top aperture is given by and where r LEC ≥ r b is the radius of the laser entrance channel. The calculation for the bottom plane, however, depends on whether the bottom plane is completely closed off with electrode material (i.e., an implosion "glide plane") or if it has an opening for a beam dump. Should a beam dump be used, then the radius of the bottom aperture is given by where r dump is the radius of the beam dump.

M. Magnetic flux loss due to the Nernst thermoelectric effect
With all of the terms in Eq. 77 thus defined, the next important quantity to calculate is the flux loss of the axial magnetic field from the fuel region to the liner region due to the Nernst thermoelectric effect. For MagLIF with a hot fuel region, flux losses due to the Nernst effect dominate over flux losses due to resistive diffusion alone. 1 To model flux losses due to the Nernst effect, we follow Ref. 83 to writė F (x e ) = 1.5x 3 e + 3.053x e x 4 e + 14.79x 2 e + 3.7703 .
Note that, like thermal conduction, Nernst losses also depend on a function of the Hall parameter, F (x e ≡ ω ce τ ei ), as well as gradients in the fuel temperature. Also, we assume that the magnetic flux lost by the fuel is transported to the liner, and thus As the axial flux is transported into the liner, it is likely that it is dissipated rapidly by resistive diffusion, which would heat the liner material somewhat. 85 However, the energy density in the liner for either case (either pure magnetic energy or ohmic dissipation leading to thermal energy) is roughly the same, and thus will affect the liner dynamics similarly given the other simplifying assumptions made in our model. Thus, for simplicity, we ignore ohmic dissipation of the axial magnetic field. Finally, we assume that the axial magnetic flux in the vacuum region is conserved, thus Φ zv (t) = Φ zv0 .

N. Erosion of the fuel's dense outer shelf region
The cold dense shelf region, generated by preheating only a central portion of the fuel, provides a buffer region in the fuel between the hot spot and the cold liner wall. This buffer region significantly reduces radiation losses by effectively reducing the amount of fuel mass that contributes to the radiation losses. It also reduces thermal conduction losses and axial magnetic flux losses from the fuel to the liner because the temperature gradient at the edge of the nearly flat shelf region is much less than that at the edge of the hot spot region. However, this buffer region is not static. It begins to erode away immediately after its formation due to thermal conduction from the hot spot region to the shelf region. To describe this erosion, we differentiate E = 3 2 N kT with respect to time, rearrange the terms, and make simple substitutions to approximate the mass transfer rate from the shelf region to the hot spot region aṡ whereT h is the average temperature in the hot spot region andT s is the average temperature in the shelf region. This states that thermal energy must be deposited in the shelf for a particle to make a "quantum" jump from having an average thermal energy of 3 2 kT s in the shelf region to an average thermal energy of 3 2 kT h in the hot spot region; the larger the temperature jump and/or mass transfer rate, the larger P ch − P cs needs to be. The average temperatures in each region are found by invoking our isobaric assumption, i.e., are the mean fuel densities in the hot spot and shelf regions, respectively, and f mh is the fraction of the total fuel mass that is part of the hot spot. We keep track of f mh rather than the absolute fuel mass in the hot spot since the absolute mass changes (due to fusion reactions and end losses). Thus we havė Note that initially (at the time of preheat), we have where if laser absorption is being calculated, r ph0 = r b . From our example above using the laser preheating parameters of recent MagLIF experiments, 16 we get Thus, initially, only about 4% of the fuel mass is part of the hot spot region, though, radially, the hot spot region occupies 21% of the overall fuel radius r g , and quickly expands to occupy ∼ 90% of r g before the shelf begins eroding away with any significance.

O. Primary fusion reaction rates
For an arbitrary deuterium to tritium fuel ratio, the DT reaction rate is given by 86 and the two dominant DD reaction rates are given by 86 where σv dt , σv dd, 3 He , and σv dd,t are the reactivity parameters for the DT, DD-3 He, and DD-T reactions, respectively. The reactivity parameters are calculated using 87,88 where the coefficients C 0 -C 7 are the fitting parameters provided in Table III.
The reactivity rates, along with end losses, provide the rate of change of the total number of deuterons and tritons in the fuel, respectively, aṡ The reactivity rates also provide the total fusion power, as given by where the energy yields per DT, DD-3 He, and DD-T reaction are, respectively, Q dt = 17.6 × 10 6 · q e (165a) Q dd, 3 He = 3.27 × 10 6 · q e (165b) Since one neutron is released per DT reaction and one neutron is released per DD-3 He reaction, the primary DT and primary DD neutron yields are given by and thus the total primary neutron yield is Finally, the total fusion energy yield is P. Summary of the semi-analytic MagLIF model The dynamics of our model are described by Eqs. 9, 11, 23, 35, 36, 37, 53, 77, 142, 147, 155, 158, 159, 160, 162, and 163. These 2N ls + 13 ordinary differential equations are repeated here for emphasis and clarity: (i = 2, 3, . . . , N li − 1) + (P r + P c + P Ω − P BB ) /N ls (s = 1, 2, . . . , N ls ) Note that if the liner current is prescribed, rather than derived from the voltage-driven circuit model, then Eqs. 170a-170c are discarded, and only 2N ls + 10 ordinary differential equations are required to describe the dynamics of our model. Furthermore, if we ignore the Nernst effect, preheat all of the fuel uniformly, use either pure D 2 fuel or a 50/50% DT mix, and we ignore dopants and contaminants, then the number of required differential equations drops to only 2N ls + 5.

III. MODEL VERIFICATION
In Fig. 4, we present SAMM simulation results where the intent was to reproduce the 1D results presented in the original MagLIF paper. 1 These SAMM results show that our simplified model does indeed capture the 1D behavior presented throughout Ref. 1. All of the simulations presented in Fig. 4 used DT fuel.
In Fig. 4(a), we have plotted the fuel radius, liner radius, drive current, and average fuel temperature as a In Fig. 4(b), we have plotted the normalized fuel temperature, magnetic field strength, and fuel density at stagnation (at peak burn rate) as a function of the normalized fuel radius for two SAMM simulations, one with the Nernst effect included (solid lines), and one without the Nernst effect included (dashed lines). All quantities are normalized to the simulation without the Nernst effect. 90  solutions were found as a function of initial fuel temperature. Using these initial conditions directly with SAMM, shown here in Fig. 4(d), we generated the results shown in Fig. 4(c). The resulting convergence ratios were 16.3±1.4, 25.7±3.4, and 32±5.5, and thus are only approximately equal to 10, 20, and 30.
In Fig. 4(e), we have plotted the optimum fusion energy yields per unit liner length as a function of the initial axial magnetic field strength. Here we have held the maximum convergence ratio fixed at either 10, 20, or 30. To find these optimum constant convergence ratio solutions, we had to scan over the input parameter space that consists of the initial preheat energy and the initial fuel density. These results should be compared with Fig. 7 Fig. 4(f), we have plotted the optimum fusion energy yields per unit liner length as a function of the peak drive current for two cases, one with α-deposition turned on, and one with α-deposition turned off. Here, as in Ref. 1, for the case with α-heating, we scanned over the input parameter space that consists of the initial preheat energy and the initial fuel density to find the optimum solutions that had a maximum convergence ratio of 20 (for the case with no α-heating, the same input parameters were used, but α-deposition was simply  turned off, thus the convergence ratios varied somewhat). In Fig. 4(g), we have plotted the ratio of the peak fuel temperature with α-heating turned on to the peak fuel temperature with α-heating turned off. We have also plotted the ratio of the total fusion energy yield to the total energy absorbed by the target (liner and fuel), i.e., the "target gain". In Fig. 4(h), we have plotted the optimum fusion energy yields per unit liner length as a function of the initial liner aspect ratio (AR). This figure should be compared with Fig. 10 in Ref. 1. In the SAMM results, the transition to higher yields with higher aspect ratios is a bit more abrupt, and then rolls off sooner; however, the general trend is the same as that shown in Ref. 1. Also, although not explicitly shown in Ref. 1, the eventual rolloff in yield with higher aspect ratio is expected. That is, to keep the implosion time constant with higher aspect ratio, the overall liner mass must be decreased while the liner's initial outer radius is increased (and while the liner's initial wall thickness is decreased). Eventually, the liner's mass decreases to a point where it cannot provide sufficient tamping and inertial confinement at stagnation.
In Fig. 4(i), we have plotted the normalized yield as a function of the premixed fraction of beryllium in the fuel. The yields are normalized to the case with no mix/dopants. This figure should be compared with Fig. 11 in Ref. 1. In Fig. 4(j), we consider laser propagation through the fuel using the analytic treatment discussed in Sec. II G.
Here we plot the laser heating front propagation distance as a function of time for the preliminary point design of Ref. 1. The laser pulse duration is 10 ns. At 5 ns (red), the front of the "bleaching wave" has reached 0.56 cm into the fuel, which is consistent with Fig. 13(b) of Ref. 1. By 8 ns (blue) and 11 ns (green), the bleaching wave has propagated completely through the 0.65-cm-long imploding fuel region, which is consistent with Figs. 13(c) and 13(d) of Ref. 1, respectively.
In Fig. 4(k), we present a color-filled contour plot of total fusion energy yield as a function of both the fuel preheat temperature and the laser beam spot size (radius). Overlaid on top of the color-filled contour plot are the contours of constant preheat energy required to obtain the preheat temperature at the given spot sizes. This figure should be compared with Fig. 15 in Ref. 1. In Ref. 1, this figure was generated by additionally scanning over a range of initial fuel densities to obtain only the solutions where the convergence ratios were 25. For simplicity, we used the input densities and preheat energies found by this previous effort to generate the SAMM results presented here in Fig. 4(k), and thus the convergence ratios were allowed to vary somewhat.
Finally, in Fig. 4(l), we have plotted results for mass end losses out of the fuel region through the laser entrance hole using the analytic treatment discussed in Sec. II L. Here we are plotting the ratio of mass remaining after fusion burn to the mass of the fuel at the start of the simulation as a function of the axial length of the liner. We have done this for three cases, where the ratio of the radius of the laser entrance hole to the initial radius of the fuel is either 1.0, 0.5, or 0.25. This figure should be compared with Fig. 16 in Ref. 1.

IV. SUMMARY
In this article, we have presented a new semi-analytic model of MagLIF. This model, implemented in a code called SAMM, accounts for: (1) preheat of the fuel (optionally via laser absorption); (2) pulsed-power-driven liner implosion; (3) liner compressibility with an analytic equation of state, artificial viscosity, internal magnetic pressure, and ohmic heating; (4) adiabatic compression and heating of the fuel; (5) radiative losses and fuel opacity; (6) magnetic flux compression with Nernst thermoelectric losses; (7) magnetized electron and ion thermal conduction losses; (8) end losses; (9) enhanced losses due to prescribed dopant concentrations and contaminant mix; (10) deuterium-deuterium and deuteriumtritium primary fusion reactions for arbitrary deuterium to tritium fuel ratios; and (11) magnetized α-particle heating. We have also shown that SAMM is capable of reproducing the general 1D behavior presented throughout the original MagLIF paper (i.e., Ref. 1).
In closing, we note that the development of this model has led to several physical insights that were not fully appreciated previously (e.g., the dependence of radiative loss rates on the radial fraction of the fuel that is preheated). Additionally, this model's accessible physics and fast run times (∼ 30 seconds/simulation) make it a useful pedagogical tool, especially for students, experimentalists, or any researcher interested in MagLIF.
search groups in general for many useful technical discussions regarding various aspects of MagLIF; in particular, we thank A. B. Sefkow, M. R. Martin, C. A. Jennings, R. W. Lemke, T. J. Awe, D. C. Rovang, D. C. Lamppa, and W. A. Stygar. Finally, we thank the personnel of the Pulsed-Power Sciences Center (including the Z and ZBL facilities) at Sandia National Laboratories; without their hard work and dedication to excellence, this work would not have been possible.
Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed-Martin company, for the United States Department of Energy's National Nuclear Security Administration, under Contract No. DE-AC04-94AL85000. is the electron thermal velocity. For the ln Λ calculation in Eq. 85, we use the following assignments:Z s → Z nuc,s , kT → 13.6 · qe · s Z ζ nuc,s , n e →n e , and n s →n s .
For the ln Λ calculation in Eq. 101, we use the following assignments: kT → kT g , n e →n e , and n s →n s .
For the ln Λ calculations in Eqs. 126 and 129, we use the following assignments: kT → kT g (r), n e → n e (r), and n s → n s (r). and zero. Conversely, if P rs (r h ) > P rv (r h ), then T B is too high and r h is too small; the correct value for r h will be found between this value that is too small and r g .

Updating ρc and Tc
Before we can test our iterative guesses for r h and T B in Eqs. 121-122, we first need to find ρ c and T c to complete the profile definitions given in Eqs. 105-109. To do this, we use our isobaric assumption to first update ρ g (r h ) =ρ gTg /T B ≡ ρ B . Next, we find the value for ρ c that makes the total mass in the hot spot, as defined by Eq. 107, consistent with m gh = f mh m g (see Appendix B 4 for additional details). Then, we again apply our isobaric assumption to get T c =ρ gTg /ρ c . Finally, with all profile parameters defined and updated, we evaluate P rv (r h ) and P rs (r h ) and test whether the results meet our desired solution accuracy.