Reduced model of plasma evolution in hydrogen discharge capillary plasmas

A model describing the evolution of the average plasma temperature inside a discharge capillary device including Ohmic heating, heat loss to the capillary wall, and ionization and recombination effects is developed. Key to this approach is an analytic quasistatic description of the radial temperature variation which, under local thermal equilibrium conditions, allows the radial behavior of both the plasma temperature and the electron density to be speciﬁed directly from the average temperature evolution. In this way, the standard set of coupled partial differential equations for magnetohydrodynamic (MHD) simulations is replaced by a single ordinary differential equation, with a corresponding gain in simplicity and computational efﬁciency. The on-axis plasma temperature and electron density calculations are benchmarked against existing one-dimensional MHD simulations for hydrogen plasmas under a range of discharge conditions and initial gas pressures, and good agreement is demonstrated. The success of this simple model indicates that it can serve as a quick and easy tool for evaluating the plasma conditions in discharge capillary devices, particularly for computationally expensive applications such as simulating long-term plasma evolution, performing detailed input parameter scans, or for optimization using machine-learning techniques.

The rapid development of plasma-based accelerator techniques, either laser driven [1,12,[21][22][23][24] or beam driven [3,[25][26][27], is made possible by advances in diagnostics and numerical modeling. Since the laser spot size and/or electron beam radius is small compared to the capillary radius, it is the near-axis plasma properties that are the most important to characterize and are the focus of plasma diagnostic techniques including longitudinal laser interferometry [28][29][30][31] and plasma emission spectroscopy [31][32][33][34]. The purpose of this work is to present a simple numerical model for evaluating the plasma properties on axis in plasma capillary discharges.
Magnetohydrodynamic (MHD) simulations have been successfully used for modeling hydrogen discharge capillary devices [8,11,16,[35][36][37]. MHD models, including the approach developed in this work, are generally applicable to * gregory.boyle@desy.de collisional plasmas with atomic density 10 23 m −3 (i.e., initial gas pressures of greater than approximately a few millibars) and ionization degrees 10 −3 . Simulations usually consist of a system of coupled partial differential equations describing the mass-density, momentum, and energy evolution in one, two, or three dimensions. Reduced geometry and simple equilibrium models have been shown to capture the essential physics for many applications [8,11,16]. These investigations have demonstrated that stable quasistatic conditions are reached during the discharge that can be well described by reduced MHD models.
The creation and subsequent evolution of a capillary plasma due to an electrical discharge is largely dictated by the local plasma temperature. During the discharge, the plasma heats up via Ohmic heating and a radial temperature gradient develops between the on-axis plasma and the cooler wall. In response to the pressure gradient, the plasma density moves away from the axis towards the boundary to reestablish a uniform radial pressure. In quasistatic equilibrium, the balance between Ohmic heating and boundary heat loss results in distinctive temperature and electron density profiles, which can be exploited for guiding high intensity laser pulses [8] and mitigated for active plasma lensing applications [17]. The plasma temperature plays the principle role in specifying the density of ionic states, as well as plasma transport properties, e.g., the thermal and electrical conductivity.
In this work, the plasma dynamics are captured via a model of the average energy evolution, i.e., a single ordinary differential equation. This is achieved through assumptions about the radial variation of the plasma properties based on quasistatic conditions. Section II describes the model of a hydrogen discharge capillary, including the assumptions made about the radial plasma temperature and density profiles during different stages of the discharge. The explicit forms of the various energy input and output mechanisms are also detailed here. In Section III the simulation results are compared against existing one-dimensional (1D) MHD simulations for a variety of discharge and pressure conditions.

II. MODEL DESCRIPTION
The most commonly used gas species in gas-filled capillary discharge devices is hydrogen [1,[9][10][11]24,38]. In this work, the discharge dynamics of a confined axisymmetric cylindrical hydrogen plasma of radius R and length L R are considered. The dynamics of the plasma discharge system are largely dictated by the local plasma temperature, and thus the focus of this work is on the dominant energy exchange processes that can occur. For hydrogen plasmas, these are Ohmic heating, the thermal exchange with the capillary wall, and the reactive energy exchanges due to ionization and recombination. Radiative energy losses are neglected, as the influence of radiation cooling on the plasma dynamics for hydrogen is insignificant for discharge currents I 1.4 MA (the Pease-Braginskii current) [8,39]. Z-pinch effects [40] are also neglected, as the magnetic pressure is small compared to the plasma pressure for the range of discharge parameters considered (see Table I).
The system is treated as a single-fluid plasma that exists in a state of local thermal equilibrium (LTE) between the electrons and ionic species. Since L R, the longitudinal variation of the plasma properties is considered negligible and only the radial variation is included. The radial energy balance equation [41] is where r and t are the radial position and time, respectively, is the total energy, P is the total pressure, v is the radial velocity, and q is the heat flux, all defined for a single-fluid plasma. Q represents the combined remaining sources and sinks of thermal energy, which here is only Ohmic heating. Assumptions underlying Eq. (1), and MHD models more generally, are that (1) the characteristic length scales are much greater than the collisional mean-free-path length, electron and ion gyroradii, and Debye length, and (2) the characteristic time scales are much greater than the collisional mean-free-path time and the inverse of electron and ion gyrofrequencies. A small Debye length implies quasineutrality, and high collisionality implies that the electron and ion velocity distributions are close to a Maxwell-Boltzmann distribution. These conditions are generally satisfied for hydrogen discharge capillary plasmas with atomic density of n a 10 23 m −3 (i.e., initial gas pressures of greater than approximately a few millibars) and ionization degree Z a 10 −3 . The initial breakdown of the plasma, which occurs during the first ≈10 ns, is a complex kinetic phenomenon which cannot be described with MHD models. Instead of modeling the breakdown, an initial temperature (e.g., T 0 = 0.3 eV) is assumed such that the plasma is already slightly ionized. For many applications, the full radial variation is not required, and a single characteristic value representing the plasma conditions, e.g., the average value or on-axis value, is sufficient. Averaging over the radial extent of Eq. (1) yields where it is assumed that there is no net exchange of material with the capillary walls, and where the averaging is defined via φ = 1 πR 2 R 0 2π rφ(r) dr. The particular form of each term in Eq. (2) is detailed in Sec. II B.
A similar expression to Eq. (2), i.e., the average representation of the plasma conditions inside a discharge capillary, was considered in Ref. [42], building upon earlier work in Ref. [43]. The key difference is that, in this work, the radial variation of the plasma properties is considered in evaluating Eq. (2), which will be shown to be critical in accurately describing the average energy evolution. A method for approximating the radial variation of the plasma temperature and electron density is hence required.

A. Radial variation of the temperature and atomic density
This section introduces a method for determining the radial temperature and atomic density, which is the cornerstone of the present work. Specifying the radial behavior directly allows the calculation of the on-axis plasma properties, average plasma properties, and, importantly, the derivative terms at the boundary which control heat flux.
The time evolution is separated into two regimes: (1) the initial uniform regime where the plasma conditions are approximately radially uniform and (2) the final quasistatic regime where the plasma temperature and atomic density vary radially so as to maintain a balance between the energy input and output mechanisms.

Transition from uniform to quasistatic conditions
At early times during the discharge, the weakly-ionizedplasma properties, such as the temperature and atomic density, are essentially uniform radially. As the plasma continues to heat, the axis becomes hotter than the constant-temperature wall, creating a temperature (and hence pressure) gradient. Ionization of the neutral species acts to absorb energy, both slowing the temperature increase and reducing the radial temperature variation. However, once the first level of ionization is near completion the plasma temperature is free to rise rapidly. At this point there is a corresponding rapid rise in pressure gradient causing the plasma density to reorganize towards uniform pressure conditions, i.e., the quasistatic state.
To accurately model the transition from the initial to quasistatic conditions requires the additional calculation of the (radially) spatially resolved density and velocity variables. However, given that the onset of the transition tends to coincide with the rapid rise in temperature nearing full on-axis ionization, the model can be vastly simplified while retaining the important physical phenomena. It is hereafter assumed that the radial pressure is always uniform, and that the plasma temperature and density transition between the uniform and quasistatic regimes occurs instantaneously at time t = t * , which is defined via the on-axis ionization fraction Z a0 (t * ) = 0.9 (see Sec. II B 1). The value of 0.9 has been chosen for its good agreement with previously published 1D simulations [8,11], which are examined in Sec. III. Alternatively, the entire plasma evolution can be simulated assuming either uniform or quasistatic conditions to establish a range of possible values.

Quasistatic conditions
The quasistatic regime is characterized by a uniform radial pressure and a plasma temperature that is described by the steady-state energy balance equation (see Appendix A), where T is the plasma temperature and κ is the plasma thermal conductivity. The precise form of Q depends on the expressions chosen for the Ohmic heating (and radiation losses, when not negligible) as discussed in Sec. II B 3. In principle the exact solution to Eq. (3) could be solved at each time step of the full average energy evolution, consistent with the instantaneous average energy. However, a faster and more efficient method is sought in this section. The thermal conductivity controls the radial redistribution of thermal energy. The total thermal conductivity κ includes contributions from the total charged plasma species (electrons and ions), κ p , and total neutral species (atomic and molecular), κ n , via the simple mixture rule [44], The specific forms of the terms in Eq. (4) are given in Appendix C. The equilibrium method introduced in Ref. [8] employs an approximation to the Spitzer-Härm theory of fully ionized plasmas [45] such that κ ∝ T 5/2 (and Q ∝ T 3/2 ). At low temperature and hence low ionization fractions, collisions with neutral species (as opposed to collisions between charged particles) dominate, resulting in a κ ∝ T 1/2 dependence. For a large proportion of the total discharge time, the plasma temperature near the capillary axis will be multi-eV [8,39], while the (constant) temperature of the capillary wall is sub-eV, indicating the existence of a layer near the boundary dominated by neutral collisions due to the low ionization fraction. The system can then be separated into two distinct regions, i.e., the central plasma-dominated bulk and the neutral-dominated boundary layer near the capillary wall.
To facilitate fast and efficient calculations, an analytic approximation for T (r) is sought. Assuming that (1) Ohmic heating effects ensure that the radial plasma temperature decreases monotonically from a maximum value on axis to the minimum value at the boundary; (2) a "two-region" approach can be employed, differentiating the plasma-dominated bulk from neutral-dominated boundary layer by an internal boundary temperature T b ; and (3) Q is approximately constant with respect to radial position, the exact magnitude of which is chosen such that T (r) in Eq. (3) is consistent with in Eq. (2) at each time step, then an analytic expression for the radial temperature profile in the range [0, R] can be derived (see Appendix A). Treating Q as a uniform energy source under quasistatic conditions in order to analytically define the radial plasma temperature is called the quasistatic uniform-energy-source temperature (QUEST) method. The QUEST method tempera-ture profile is where T 0 , T w , and T b are the temperature on axis, at the wall r = R, and at the internal boundary r = R b , respectively. Clearly when T 0 , that the temperature range spans both the plasma-dominated and neutral-dominated regions, but can be altered easily for other situations.
The plasma-dominated regime is here defined by κ p > κ n , and conversely the neutral-dominated regime by The κ components are weakly dependent on the atomic density, and so in the simulations the T b value corresponding to the initial n a is used. For n a = 10 24 m −3 , T b ≈ 0.9 eV, and this value is used hereafter. An order of magnitude change in n a results in 5% change in the value of T b . The corresponding change in the average plasma temperature calculations is 2%, indicating that the simulation procedure is robust to the choice of T b .
The value of R b can be completely specified by the requirement that the heat fluxes from each region, which obey different temperature power laws, match at the internal bound- The expression for R b in terms of the on-axis temperature T 0 and wall temperature T w is the derivation of which is given in Appendix B. Thus the full radial temperature profile (and R b ) is specified by T 0 , T b , and T w . In the course of a simulation, T b and T w are set as input constants, and only T 0 varies as a function of time.
Example radial temperature profiles, corresponding to select average temperature values, are shown in Fig. 1. Different behavior is demonstrated either side of R b , owing to the different temperature power laws controlling the thermal conductivity in each region. As T 0 increases, the position of the plasma-neutral boundary R b moves towards the capillary wall. It should be noted that R b < R and the heat flux at the capillary boundary is dictated by the neutral-dominated thermal conductivity regardless of how thin the neutral-dominated boundary layer becomes.
The nonuniform plasma temperature described by Eq. (5) implies a nonuniform plasma density under uniform pressure conditions P(r) = P . Assuming uniform total pressure, it follows from the ideal gas law that where Z a is the ionization fraction as defined in Eq. (9), and thus the radial atomic density n a (r) is fully specified by T (r) under LTE conditions. For the trivial case that all properties are radially uniform, i.e., during the initial uniform regime, n a = n a . The radial plasma temperature and electron density profiles resulting from the QUEST method are compared to those from 1D MHD simulations in Sec. III.

B. Energy balance terms
In this section, the particular forms of the energy balance terms necessary to evaluate Eq. (2), i.e., the internal energy , Ohmic heating Q Ohm , and boundary heat flux q(R), are detailed.

Density of ionic states
A single-temperature plasma that exists in a local thermal equilibrium between the electrons and heavy ionic species is assumed. Reference [39] showed that, for a hydrogen discharge capillary, LTE conditions are established in approximately 50 ns, and Refs. [8,11,35] have had success modeling discharge capillaries assuming LTE conditions over the entire discharge lifetime.
By assuming LTE conditions, the density of ionic states is fully specified by the plasma temperature (T ) and atomic density (n a ) via the Saha ionization equation [46]. For a quasineutral hydrogen plasma only single-level ionization is required, and the appropriate Saha ionization equation is where I H is the ionization energy for hydrogen, and Z a = n e /n a is the mean charge per atom which here also represents the ionization fraction. The constants m e , k b , and h are the electron mass, Boltzmann constant, and Planck constant, respectively.
The ion density n i and neutral density n n can be found from the quasineutrality, n i = n e , and particle conservation, n n = n a − n i , conditions, respectively.

Internal energy
The connection between the internal energy and the plasma temperature, accounting for the energy stored in ionization, is given by where C v,a = 3 2 n a k b and C v,e = 3 2 n e k b are the atomic and electronic heat capacities for ideal gases, respectively. The potential energy term U = n e I H represents the amount of ionization energy required to reach the specified density of ionic states from a neutral state.
The time derivative of the internal energy can be rewritten as a function of temperature directly, i.e., where C v then represents an effective heat capacity. The calculation of ∂Z a ∂T is detailed in Appendix D. Note that Z a and ∂Z a ∂T are simply functions of T and n a .

Ohmic heating
The discharge current provides the energy input to the plasma system via Ohmic heating. The Ohmic heating contribution to Q in Eq. (1) is where J is the current density and E is the electric field strength. The connection between the electric field strength and the current density is given by Ohm's law, J = σ E , where σ is the electrical conductivity. The Ohmic heating is driven predominantly by electron interactions, such that the electrical conductivity of a plasma consisting of a mixture of electrons, ions, and neutral species is given by [44] 1 where σ ei and σ en are the electrical conductivity in the fully ionized and weakly ionized limits, respectively. The specific forms of the terms in Eq. (14) are given in Appendix C.
Following the quasistatic approach in Ref. [8], it is assumed that the electric field is homogeneous such that where I = R 0 2π rJ dr is the total current. The current amplitude as a function of time is routinely measured in discharge capillary experiments, and thus I (t ) is treated as an input rather than calculated in an additional coupled-circuit model [42,43].

Boundary heat loss
The dominant energy loss mechanism in (enclosed) hydrogen discharge capillaries is the heat flux through the capillary boundary. The heat flux is given by Fourier's law, q = −κ ∂T ∂r , such that Equation (16) explicitly depends on the radial temperature gradient at the boundary, and thus can be written in terms of T w , T b , and R b via Eq. (5). The main reason to decompose the domain into plasma-dominated and neutral-dominated regions (see Sec. II A) is to capture this term accurately.
At the capillary boundary the thermal conductivity, and hence the energy transfer to the wall, is dominated by neutral collisions due to the low local temperature and ionization fraction. This is in contrast to the plasma bulk where the electron thermal conductivity is generally much larger than the neutral (and ion) species thermal conductivity. The melting point of sapphire capillaries is approximately 2300 K, and in this work T w = 2000 K is used in the simulations. The simulation procedure is very robust to the choice of T w value, with a change of 50% in T w resulting in a 1% change in the calculated average plasma temperature.

C. Numerical solution
Each of the transport properties described in Sec. II B and Appendix C is fully specified by the plasma temperature and atomic density (assuming LTE conditions). Thus the Taylor series approximation of each of the radially varying quantities, f (T (r), n a (r)), in the neighborhood of some reference values, T and n a , is It is assumed that the appropriate reference values, i.e., where the dominant contribution to the average occurs, are the aver-age plasma temperature T = T and average atomic density n a = n a , such that f (T, n a ) ≈ f (T , n a ) and the energy evolution equation, Eq. (2), becomes Note that q(R) only depends on T indirectly through the temperature gradient at the boundary (see Sec. II A). The O(T − T ) and O(n a − n a ) terms, arising from the derivatives in Eq. (18), are more in the nature of correction terms, and have been neglected. Thus it is expected that Eq. (19) works best when the plasma properties are only slowly varying functions of radial position. The validity of this approach is expanded on in Appendix E. An overview of the workflow for the QUEST simulation code is given in Fig. 2. At each time step, the method described in Sec. II A is employed to determine the radial plasma temperature profile consistent with the average temperature. This specifies the remaining transport coefficients and energy balance terms described in Sec. II B. The ordinary differential equation [Eq. (19)] is advanced using a fourth-order explicit Runge-Kutta routine [47]. The radial behavior of the atomic density (and thus the on-axis atomic and electron densities) are determined from Eqs. (7) and (8), in both the uniform and quasistatic regimes.
In comparison to the single ordinary differential equation QUEST algorithm, the 1D MHD simulations of Ref. [8] evolve a system of five coupled partial differential equations. Simulations using the QUEST algorithm typically complete in <1 s on a desktop computer. This indicates that QUEST simulations are particularly useful for computationally expensive problems, such as performing detailed input parameter scans or investigating the long-term (10+ μs) plasma evolution, where full MHD simulations are prohibitively expensive. It also makes optimization of discharge capillary plasma conditions with machine-learning techniques feasible. The simulation results are compared in the following section.

III. SIMULATION BENCHMARKS
The principle goal of the QUEST method is to reproduce the plasma temperature and electron density results of more complex 1D MHD simulations, in a quasianalytic and sig-  [8] and blue crosses represent calculations using the simplified equilibrium model in Ref. [8]. The colored shaded bands in (d) and (e) represent the range of values between assuming the uniform regime and the quasistatic regime. The shaded grey region indicates the times at which the QUEST algorithm assumes uniform conditions, and marks the transition between the uniform and quasistatic regimes corresponding to the discrete jump in the electron temperature and density profiles. nificantly less computationally expensive simulation. In this section, the QUEST simulations are compared to the previous 1D MHD investigations of Refs. [8,11] for a range of discharge current amplitudes and initial gas pressures. After establishing the validity of the QUEST approach, the importance of accurately representing the radial variation of the plasma properties, particularly those close to the capillary boundary, in accurately describing the evolution of the average quantities is demonstrated by comparison with the results of Ref. [42]. The conditions for each simulation are detailed in Table I. A comparison of the plasma temperature and electron density evolution calculated with the QUEST method and with TABLE I. Table of parameters for select plasma discharge capillary simulation literature. R is the capillary radius, P is the initial gas pressure, n a is the initial atomic density, and I p and t p represent the magnitude and time of the discharge current peak. The discharge profiles in simulations B and G1-G6 have the form I (t ) = I p sin (πt/t p ). The discharge profile in simulation C does not have an analytic form, and has been digitized for the comparisons in the present work.  [42] the 1D nonideal MHD simulations of Ref. [8] in a hydrogen discharge waveguide study is shown in Fig. 3. At early times (0-50 ns), represented by the grey shaded region in Figs. 3(c)-3(e), uniform radial temperature and density are assumed. The on-axis plasma temperature shown in Fig. 3(a), and electron density shown in Fig. 3(b), from Ref. [8] are very well reproduced by the QUEST method in the uniform regime. The slow rise in the temperature for the first 50 ns is due to substantial energy being absorbed by the ionization process. The radial profiles in Figs. 3(a) and 3(b) corresponding to 40 ns show good agreement. The 1D MHD profiles exhibit some nonuniformity near the boundary but are predominantly uniform.
At late times (75-150 ns) the results from Ref. [8] are also very well reproduced by the QUEST on-axis temperature and electron density in the quasistatic regime. The radial profiles corresponding to 80 and 100 ns show consistent nonuniform behavior between Ref. [8] and QUEST results. The analytic temperature form in Eq. (5) varies more sharply towards the boundary compared to Ref. [8], resulting in electron temperature profiles that are more sharply peaked. However, the overall agreement is very good. Further radial profile agreement can be expected from using a temperature profile shape that is equivalent to the equilibrium model shape in Ref. [8] but comes at the cost of requiring a numeric rather than analytic solution.
The discrepancy in the intermediate time range of 50-75 ns is due to treating the reorganization of the plasma from uniform to quasistatic regimes as an instantaneous process (see Sec. II A). Although the transition onset time of 50 ns is approximately correct, the transition process takes approximately 20 ns according to MHD simulations, rather than being instantaneous. This is emphasized by the fact that the 1D MHD electron density and temperature results smoothly transition between the QUEST uniform and quasistatic regime bands. The 1D MHD radial profiles corresponding to 60 ns occur during this transition, and hence show behavior that is partway between the uniform and quasistatic regimes, and is thus not well reproduced by the QUEST simulation.
The relationship between the on-axis temperature and the (time-dependent) current discharge amplitude is explicitly shown in Fig. 3(f). There are two distinct temperature "paths" corresponding to heating (lower path) and cooling (upper path) phases, i.e., on which side of the 250 A current peak is being sampled. A simplified equilibrium model from Ref. [8], which is a function of the instantaneous current amplitude, rather than being connected to the average energy evolution, is also included in Fig. 3(f), represented by blue crosses. The equilibrium model provides an identical relationship between T 0 and I during both the heating and cooling phases, and demonstrates good agreement for the cooling phase, particularly near the current peak. However, naturally it does not well represent the heating phase, and cannot describe times after the discharge has turned off (if I = 0, then the equilibrium temperature, etc., are also zero). Although both the equilibrium model of Ref. [8] and the QUEST model are based on a power-law temperature dependence of the transport coefficients, it is clear that the temporal evolution of the average energy must be included to satisfactorily describe the full discharge current lifetime. Figure 4 features on-axis simulation results from Ref. [11], where the authors investigated the effect of significant changes in discharge current magnitude and pressure on the formation of plasma waveguides, and thus represents an ideal range of benchmark conditions for the QUEST method. Many of the comments in the discussion of Fig. 3 apply here too.
In cases G1-G3, G5, and G6, the onset time of the transition is well reproduced by the QUEST method. In the case of G4, the plasma temperature (and ionization fraction) increases very slowly and the transition threshold of Z a0 = 0.9 is not reached until 200 ns. According to MHD simulations, the transition begins approximately 50 ns earlier than predicted using the QUEST method, and it is not clear that quasistatic conditions have been established by the culmination of the discharge. This slow transition between the uniform and quasistatic regimes cannot be accurately modeled by the QUEST approach.
Overall, the QUEST calculations and MHD simulations from Ref. [11] agree very well, particularly in the uniform and quasistatic regimes. The average difference between the QUEST calculations and Ref. [11] over the full discharge profile is 5% for the on-axis plasma temperature, and 10% for the on-axis electron density, for each condition G1-G6. The maximum difference is 40% for both properties, and occurs at the transition between uniform and quasistatic regimes. Better overall agreement is observed for discharge conditions that lead to higher temperatures (i.e., higher currents or lower densities) as these tend to demonstrate sharper transitions.
In Ref. [42] a similar approach to describing the evolution of the average plasma properties was proposed. However, in the formulation of Ref. [42] the treatment of the radial variation of the plasma parameters is substantially different from the present work. A comparison of the average plasma temperature and electron density calculated with the QUEST method and the simulations from Ref. [42] for hydrogen is shown in Fig. 5, and demonstrates considerable disagreement. These differences are significant in both the magnitude and behavior, which indicates an inherent incompatibility between the two approaches.

FIG. 5. Comparison of QUEST simulation results and those from
Ref. [42] for average electron temperature T and electron density n e as a function of time. The simulation conditions are given by "C" in Table I. Solid lines indicate the QUEST method results, while dotted lines indicate the (digitized) simulation results from Ref. [42].
The QUEST radially averaged temperature T is considerably greater over most of the time range. Although Ref. [42] explicitly includes radiative energy losses, the effect is insignificant (less than 0.01% of the dissipated power [39]). The larger peak average temperature indicates a difference in the balance between Ohmic heating and wall heat loss for the two approaches. The energy exchange with the capillary boundary is the dominant heat loss mechanism in hydrogen capillaries, which depends critically on the radial temperature derivative at the boundary, as described in Sec. II B 4. A key component of the QUEST method is the precise representation of this boundary temperature derivative, which differs from the formalism of Ref. [42]. Another important difference is that the effect of ionization and recombination energy exchanges is included in the QUEST model. The energy absorbed during ionization (up to ≈75% of the Ohmic heating power) is responsible for the slow temperature increase at early times, and the release (up to ≈50% of the wall energy loss) during recombination is responsible for the slow temperature decrease at late times.
The peak average electron density from Ref. [42] is the same as the QUEST calculation when assuming uniform regime conditions. However, the transition onset time is predicted to be approximately 35 ns, and the subsequent behavior is calculated in the quasistatic regime. Note that the cooler (and hence less ionized) plasma near the capillary boundary contributes substantially to the averaging due to the high atomic density under quasistatic conditions, reducing the average electron density. The difference in the electron density decrease at late times is due largely to the difference in the plasma temperature evolution predicted by the two methods, as discussed previously.

IV. CONCLUSION
It has been shown that the on-axis plasma temperature and electron density calculated in existing full 1D MHD simulations, which solve a complex system of coupled partial differential equations, can be remarkably well reproduced by the QUEST (quasistatic uniform-energy-source temperature) method, which solves a single, simplified ordinary differential equation for the average plasma temperature evolution. This paves the way for investigations of computationally expensive capillary discharge problems, such as characterizing the longterm plasma evolution, performing detailed input parameter scans, or for employing machine-learning-based optimization techniques, which are infeasible using more complex simulation tools.
The key to the QUEST method is in the assumptions made about the radial temperature behavior, which then specify the remaining plasma properties under local thermal equilibrium conditions. The approach followed here is to split the temporal evolution of the plasma into a "uniform regime," where the plasma temperature is radially uniform, and a "quasistatic regime" where the plasma temperature has a nonuniform but analytic representation under quasistatic conditions. Particular attention is given to the quasistatic radial temperature representation, which is separated into plasma-dominated and neutral-dominated regions, as it determines the heat flux at the capillary boundary-the major energy loss process in these hydrogen discharge capillary systems.
The near-axis plasma properties are the most relevant to many experiments, particularly in beam-driven wakefield acceleration. The on-axis plasma temperature and electron density are compared to the full 1D MHD simulations of Refs. [8,11] for a range of discharge current amplitudes and initial gas pressures. The substantially simpler QUEST method demonstrates good agreement, particularly at early and late times where either uniform or quasistatic conditions have been established. The plasma temperature and electron density are generally within 5% and 10% of Refs. [8,11], respectively. At intermediate times, the 1D MHD results exhibit a mixture of uniform and quasistatic behavior; however, the QUEST method still gives results with differences 40%. When compared to the simplified equilibrium model of Ref. [8], the QUEST method demonstrates that modeling the evolution of the average energy is necessary to adequately describe the plasma conditions over the full discharge current lifetime.
This marks the first time that a model based on the evolution of the average energy in capillary discharge devices has been satisfactorily benchmarked against 1D MHD simulations over the entire discharge profile, and the results here indicate an incompatibility with previous approaches [42,43].
In Ref. [31] it was shown that evaluating the plasma temperature to within a relative error of ≈100% was necessary for agreement between plasma diagnostics based on emission spectroscopy and laser interferometry. The demonstrated success of the QUEST method indicates that it can be used in conjunction with plasma emission spectroscopy techniques to evaluate the electron density from measured emission spectra [48,49]. Future investigations will explore the use of QUEST simulations in plasma cell characterization experiments.
was supported by Helmholtz ARD and the Helmholtz IuVF ZT-0009 program.

APPENDIX A: QUEST TEMPERATURE PROFILE T (r)
Under steady-state conditions, the energy balance equation (1) becomes where Fourier's law for the heat flux q = −κ dT dr has been employed, and κ is the plasma thermal conductivity. If κ = κ 0 T n where κ 0 is a constant, then for radially uniform Q the integration can be performed analytically over the range r ∈ [r L , r R ], yielding where T L = T (r L ) and T R = T (r R ) are the temperatures at each end of the range. Expressions for the heat flux term, and average temperature, follow directly. The average temperature simplifies to T ≈ ( n+1 n+2 )T L when T L T R . In this work, a plasma-dominated region is distinguished from a neutral-dominated region, corresponding to n = 5/2 and n = 1/2, respectively.
where the radial temperature profile from Eq. (5) and heat flux from Eq. (A4) have been employed. Thus the position of the internal boundary R b is specified by T 0 , T b , and T w . These expressions assume that T 0 > T b . When this is not the case, R b /R = 0; i.e., the entire domain is neutral dominated.

APPENDIX C: PLASMA TRANSPORT PROPERTIES
The key transport properties in this work that control the source and redistribution of the thermal energy within the plasma are the electrical conductivity and the thermal conductivity. The electrical conductivity controls the Ohmic heating, which is the main energy input, and the thermal conductivity controls the redistribution of the thermal energy and loss to the capillary wall, which is the main energy output. These both depend on the collision frequency between the electrons, ions, and neutrals.
Here, the neutral species can be atomic or molecular hydrogen. The cross sections for binary interactions involving atomic or molecular hydrogen (including interactions with the charged species) are all of a similar magnitude for the range of plasma temperatures considered in this work [50], and are thus approximated by a single hard-sphere scattering cross section for binary collisions with an "effective" hydrogen neutral. Near the cold capillary wall (T w = 2000 K), under LTE conditions, the neutral species is predominantly molecular [46], and so the hard-sphere scattering cross section for hydrogen gas [51] has been chosen to ensure that the thermal conductivity at the capillary wall is described correctly.
The total thermal conductivity κ includes contributions from the total charged plasma species, κ p , and total neutral species, κ n , via the simple mixture rule [44], where m a is the atomic mass, and ν ab represents the collision rate of species a with b, where e, i, and n represent electrons, ions, and neutrals, respectively, which are detailed below. The coefficients of ν ei and ν ii are taken from Ref. [52] such that Eq. (C2) is consistent in the fully ionized limit. The heavy species-electron collision rates ν ie and ν ne are typically smaller than ν ii and ν nn , respectively, by a factor of √ m a /m e , and are thus not included in this work.
The total electrical conductivity σ of a plasma consisting of a mixture of electrons, ions, and neutral species is given by [44] 1 where σ ei and σ en are the electrical conductivity in the fully ionized and weakly ionized limits, σ ei = 1.96 n e e 2 m e ν ei , σ en = n e e 2 m e ν en , where ν ei and ν en are the electron-ion collision and electronneutral collision rate, respectively, detailed below. Although electron-electron collisions are momentum conserving and do not contribute directly to Eq. (C4), the indirect effect of electron-electron correlations on the electron velocity distribution is included in the numerical coefficient of σ ei , which is taken from Ref. [52]. The electron-ion collision rate ν ei [8] and electron-neutral collision rate ν en [44] are given by 2π m e e 4 n e ln λ ei where a = 145 pm is the kinetic radius for hydrogen gas [51], and ln λ ei is the electron-ion Coulomb logarithm [8] here defined as The Coulomb logarithm is the approximation of a diverging collision integral, and is generally of order 10. In the simulations a floor is applied to the Coulomb logarithm, i.e., max(ln λ ei , 1 2 ln 2), to control the Coulomb collisions at low temperatures [53].
The heavy-species collision rates including ion-ion collisions ν ii [8], ion-neutral collisions ν in , neutral-ion collisions ν ni , and neutral-neutral collisions ν nn are calculated via [44] ν ii = 4 3 π m a e 4 n i ln λ ii where once again the kinetic radius a = 145 pm [51], and where ln λ ii is the ion-ion Coulomb logarithm [8], here defined as Similar to the electron-ion Coulomb logarithm Eq. (C7), a floor is also to the ion-ion Coulomb logarithm, i.e., max(ln λ ii , 1 2 ln 2). Note that, due to the identical masses (ignoring the negligible mass of the electron) of the neutral and ion species, n i ν in = n n ν ni . The use of the hard-sphere scattering model for all neutral collisions, along with a single temperature, results in ν in = ν nn .

For a quasineutral single-level ionization plasma the appropriate Saha ionization equation is
where I H is the ionization energy for hydrogen, and Z a = n e /n a is the mean charge per atom which here also represents the ionization fraction. The constants m e , k b , and h are the electron mass, Boltzmann constant, and Planck constant, respectively. The solution for Z a is then and the derivative with respect to temperature is The ionization state described by Z a and dZ a dT is completely specified by the local plasma temperature T and atomic density n a .

APPENDIX E: VALIDITY OF THE ZERO-ORDER TAYLOR SERIES EXPANSION
The transport properties controlling the plasma dynamics are functions of the local plasma temperature and atomic density. In the zero-order Taylor series approximation [Eq. (18)] it is assumed that the appropriate reference values T and n a are the average plasma temperature T and average atomic density n a , respectively. Thus all radially varying plasma properties are evaluated directly at T and n a to approximate the average value.
In general, the transport properties described in Sec. II B are only weakly dependent on the atomic density, and can be well approximated by plasma temperature power laws. The success of the zero-order Taylor series expansion largely depends on how well the average of these power-law functions can be approximated as a function of the average directly, i.e., how close a parameter ζ (p) = T p T −p is to unity. From Eq. (A3) it follows that, for T L T R , ζ (p, n) ≈ n + 1 n + 1 + p n + 2 n + 1 where now ζ (p, n) is a function of two variables to indicate the dependence on radial temperature power-law index n as well as the power to which the temperature is being raised, p. A plot of ζ (p, n) vs p for plasma-dominated (n = 5/2) and neutral-dominated (n = 1/2) limits is shown in Fig. 6. The  FIG. 6. Variation of ζ (p, n) parameter, defined in Eq. (E1), with power p for two temperature power laws n. The plasma limit corresponds to n = 5/2, while the neutral limit corresponds to n = 1/2. ζ (p, n) is generally close to unity, particularly for plasmadominated conditions, which contributes to the remarkable success of the zero-order Taylor series approximation.
A comparison of the average plasma temperature evolution calculated using (1) the zero-order Taylor series expansion approach (see Sec. II C) and (2) the full radial variation of the plasma temperature and associated plasma parameters when evaluating the quantities in Sec. II B, is shown in Fig. 7 for a select range of discharge conditions. The agreement is remarkably good considering the significant approximation involved in the zero-order Taylor series truncation, with the relative errors being <10%. Simulations using the truncated Taylor series approach are approximately two orders of magnitude faster than with the full radial variation, and thus represent an extremely fast and efficient method of estimating the plasma temperature and density in discharge capillary systems.