The atmospheric pressure arcs have recently found application in the production of nanoparticles. The distinguishing features of such arcs are small length and hot ablating anode characterized by intensive electron emission and radiation from its surface. We performed a one-dimensional modeling of argon arc, which shows that near-electrode effects of thermal and ionization non-equilibrium play an important role in the operation of a short arc, because the non-equilibrium regions are up to several millimeters long and are comparable to the arc length. The near-anode region is typically longer than the near-cathode region and its length depends more strongly on the current density. The model was extensively verified and validated against previous simulation results and experimental data. The Volt-Ampere characteristic (VAC) of the near-anode region depends on the anode cooling mechanism. The anode voltage is negative. In the case of strong anode cooling (water-cooled anode) when the anode is cold, temperature and plasma density gradients increase with current density, resulting in a decrease of the anode voltage (the absolute value increases). Falling VAC of the near-anode region suggests the arc constriction near the anode. Without anode cooling, the anode temperature increases significantly with the current density, leading to a drastic increase in the thermionic emission current from the anode. Correspondingly, the anode voltage increases to suppress the emission, and the opposite trend in the VAC is observed. The results of simulations were found to be independent of sheath model used: collisional (fluid) or collisionless model gave the same plasma profiles for both near-anode and near-cathode regions.

## I. INTRODUCTION

The atmospheric pressure arcs recently found application in the production of nanoparticles, such as carbon nanotubes^{1–5} and boron-nitride nanotubes.^{6} The distinguishing features of such arcs are typically short length (about several millimeters inter-electrode gap) and hot ablating anode, which provide feedstock material for the growth of nanoparticles. A high anode temperature leads to thermionic electron emission and intensive radiation from its surface. The high electron emission from the hot anode strongly affects the near-anode plasma, especially the space-charge sheath and the heat flux to the anode. Because the arc is short, non-uniform near-electrode regions of thermal and ionization non-equilibrium play an important role in the arc physics.

Electrode ablation significantly increases the complexity of the arc physics, because such arc plasma consists of a mixture of gases with different ionization and transport characteristics. These processes will be considered in the follow-up publications. In this paper, we focus only on the non-equilibrium processes in the arc with non-ablating electrodes, e.g., a short argon arc with tungsten electrodes.

Although arcs have been extensively studied earlier, most of the studies considered only specific aspects of the arc physics: cathodic region,^{7–13} arc column,^{14,15} and anodic region.^{16–24} A detailed review of works on the argon arc modeling can be found in Ref. 25. However, we could not find papers considering arc with a hot emitting anode. Thorough reviews on numerical and experimental studies of the near-anode region of arc discharges can be found in Refs. 26 and 27. Typically, in simulations of the whole arc, thermal or chemical equilibrium plasma model (Saha equation) is used^{28–32} or space-charge sheath effects are not taken into account.^{33–35} The effects of the near-electrode non-equilibrium layers were taken into account in recent arc simulations^{36} making use of specific boundary conditions; bulk plasma was considered equilibrium. The layers were considered infinitely thin, which may be inaccurate for short arcs. Numerical simulation of the whole arc making use of the non-equilibrium plasma model was only recently reported in Ref. 37. The authors of Refs. 36 and 37 focused primarily on plasma–cathode interaction, and for the anode, some simplifications were still used. A thorough non-equilibrium fluid model of plasma of the atmospheric and high-pressure arcs was developed in Ref. 13. The governing equations for species transport and heat transfer, and coefficients therein are derived from kinetic theory.^{38,39} The model is validated by a comparison with the experimental data, in particular, see validation for the cathodic part of argon arc in Ref. 13.

In the present study, this model^{13} was expanded and implemented into a one-dimensional (1D) code for self-consistent simulations of the whole arc, including heat transfer in cylindrical electrodes, and radiation from their surfaces. For short arcs, in which axial gradients are much higher than radial, 1D approximation is rather accurate. The code results were benchmarked against previous simulations^{13} and validated against experimental data.^{40,41} Different models of space-charge sheath regions were implemented into the code: solving the Poisson equation for the collisional sheath-plasma^{13} and the quasineutral plasma with sheath boundary conditions for the collisionless sheath model.^{42,43} The models were applied to both near-anode and near-cathode regions and compared with each other.

Parametric studies of short atmospheric pressure argon arc with tungsten electrodes were performed for various current densities and inter-electrode gap sizes. The non-equilibrium effects in the near-electrode regions were studied. Anodes with and without water cooling were considered. The effect of electron emission on current-voltage characteristic of the near-anode layer was investigated. The analytical formulas for the scaling of non-equilibrium regions' widths and Volt-Ampere's characteristics of these regions and the whole arc are given in the accompanying paper.^{44}

The paper is organized as follows: In Sec. II, the governing equations and boundary conditions for plasma and electrodes are presented. Section III describes the numerical procedure for solving the governing equations. The results of simulations including the validation of the model and parametric studies of the arc are presented and discussed in Sec. IV. Conclusions of this work are summarized in Sec. V.

## II. The 1D ARC MODEL

The model consists of equations for the plasma region, corresponding boundary conditions, and equations for heat transport in the electrodes. The equations for the plasma region are derived from the kinetic theory (see Ref. 13 and references therein for details).

### A. Governing equations for plasma region

#### 1. Momentum balance of species

The momentum balance of individual species in multi-component mixture can be described by Stefan-Maxwell equations,^{39} which accurately account for the species cross-diffusion

where the subscripts *k, j* denote different species: argon atoms $a$, argon ions $i$, and electrons *e*; $nk$, $n$ are the number density of species *k* and plasma and background gas as a whole; $\rho k$, $\rho $ are the mass density of species *k* and plasma and background gas as a whole; $v\u2192$ is the mean (mass-averaged) velocity; $E\u2192$ is the electric field; *k* is the Boltzmann constant; *e* is the elementary charge; $Zk$ is the charge number: −1 for electrons, 1 for ions, 0 for neutrals; $Ckj$ are the coefficients derived from the kinetic theory (typically about unity); $Cia=1$, $Cei=0.506$; definitions of other coefficients are more complicated and can be found in Ref. 13. $R\u2192kth=Ck(e)nkk\u2207Te$ is the thermal-diffusion force, where the term with gradient of temperature of heavy particles is neglected. Kinetic coefficients $Ck(e)$ are defined as follows:

$P$ is the ratio of electron–atom and electron–ion collision frequencies

$Ctdei$, $Ctdea$ are the values of thermal diffusion coefficient for electrons in the limits of strongly and partially ionized plasmas. $Ctdei=0.703$ and $Ctdea$ is defined via integral of electron-atom collision cross-section.^{13,38}$Dkj$, $Tkj$ are the binary diffusion coefficients and temperatures defined as follows:

$mk$, $Tk$ are the mass of particles of sort k, is their temperature. Temperatures and masses of heavy particles are very close and are not distinguished in the model: $Ti=Ta=T$, $mi=ma=mAr$. $\sigma kj$ are the momentum transfer cross-sections. For electron–atom collisions, cross-section data were taken from Ref. 45, and the approximate value of electron-atom momentum transfer cross-section is $3\xd710\u221220\u2009m2$. For ion–atom interactions, the charge-exchange cross-section was used:^{46} $\sigma ia=9.2\xd710\u221219\u2009m2$. Coulomb collisions between electrons and ions can be described by the following cross-section:^{47}

where $\epsilon 0$ is the vacuum permittivity, and $ln\Lambda =ln(8\pi \epsilon 0kTe\epsilon 0kTe/ne/e3)$ is the Coulomb logarithm.

Note that Eq. (1) is written in a general multi-dimensional form. In the code, it was implemented in the 1D form. In 1D approximation for non-ablating electrodes, the mean velocity is zero everywhere and momentum conservation equation for the mixture as a whole takes the form

#### 2. Species conservation equations

For each of the species, the conservation equation can be written in the following form:

where $\Gamma \u2192k=nkv\u2192k$ is the flux of species *k*, and $sk$ is the volumetric source

where $ki$, $kr$ are the reaction rate coefficients. $ki$ is calculated using the approach described in Ref. 48, and $kr$ is calculated making use of the ionization equilibrium condition

where $ne,Saha$ is the equilibrium number density defined by the Saha equation

Here, $gi/ga$ is the ratio of statistical weights of ground state and ionized state [for argon, this ratio is equal to 6 (see Ref. 49)]; *h* is the Planck's constant; and $Eion$ is the ionization energy of argon atoms.

Because the arc is short (shorter than the radii of both electrodes) and the currents are below 100 A (the arc self-magnetic field is small and should not affect plasma profiles), the flow in the inter-electrode region can be considered truly 1D; derivatives are predominantly along the arc axis. Using this assumption, we rewrite equations of conservation of fluxes for heavy and charged particles

Here, $j$ is the current density, which we use as the input parameter of the model. At the boundaries of the plasma domain (walls of the electrodes), the total flux of heavy particles is equal to zero; therefore, it is zero everywhere [Eq. (10)].

#### 3. Energy transport equations

The energy transport of electrons and heavy particles gas is described by equations

Here, $\lambda e=(\lambda e,i\u22121+\lambda e,a\u22121)\u22121$, $\lambda h=\lambda h,i+\lambda h,a$ are the thermal conductivities of electron gas and heavy particles, $\lambda e,i=3.2knenDei/ni$; definitions of $\lambda e,a,$$\lambda h,i$, $\lambda h,a$ can be found in Ref. 13.

$Qe\u2212h=Ae\u2212H(Te\u2212T)$ is the volumetric heat exchange between electrons and heavy particles, $Ae\u2212H=8ne2\sigma ei2kmeTe\pi km$.

$Qion=si\u2009Eion$ is the heat source due to ionization/recombination processes, and

is the volumetric heat loss due to radiation; it is supposed that plasma is optically thin and the radiation from plasma is a net sink of energy^{13}

are kinetic coefficients.^{39}

The left-hand side of Eqs. (12) and (13) represents the convective heat flux carried by electrons or heavy particles, respectively. First terms on the right-hand sides of Eqs. (12) and (13) correspond to the thermal conduction. The second terms on the right-hand sides of these equations represent the effects of the Joule heating.

Note that due to differences in the mobilities of electrons and ions, the ion flux is typically significantly smaller than the electron flux: about ten times smaller in the near-cathode region and about two orders of magnitude smaller elsewhere. Therefore, the second term under “nabla” operator on the left-hand side of Eq. (12) can be neglected.

#### 4. The Poisson equation for the self-consistent electric field

Outside the space-charge sheaths, Gauss's law can be reduced to the quasi-neutrality approximation

#### 5. Equation of state

The plasma is treated as a mixture of ideal gases of electrons and heavy particles

The pressure variation is determined by Eq. (18). Substitution of Eq. (15) into Eq. (18) gives an algebraic relation for pressure

where the second term in the left-hand side part is often referred as the electrostatic pressure. $p0$ is the reference pressure or pressure in a far-away location where there is no electric field. It is an input parameter of the model (1 atm, in the considered case). Note that even in the cathodic space-charge sheath where the electric field is at its highest, it is about 10 V per 1 *μ*m and the electrostatic pressure is about 10^{3} Pa, which is much lower than $p0$. Therefore, the total gas and plasma pressure in Eq. (17) can be approximated by a constant pressure, $p0$ throughout the arc.

### B. Boundary conditions including sheath effects

To define the solution, boundary conditions are required at the plasma-electrode surfaces. The temperature of heavy particles should be equal to the temperature of the electrode's surface

The temperature of the electrode can either be predefined (in the case of water-cooled anode^{50}) or has to be determined from the solution of the heat balance equation as described in Sec. II C.

The boundary conditions for other equations are formulated in terms of fluxes. The net electron flux at the plasma-electrode boundary is a composition of two oppositely directed fluxes: plasma electrons $\Gamma e,bplasma$ and electrons emitted from the electrode surface $\Gamma e,bemiss$

Here and further, subscript *b* relates to the plasma-electrode boundary, $en,b$ denotes the projection to the x-axis of the outward pointing unit normal from the electrode surface (see Fig. 1). The ion flux at the boundary is presented by ions' fall from plasma onto the electrode's surface

Addition of ion current density to the electron current density gives the total current density at the surface

Here, $j$ is the total current density, which is assumed to be constant and known in the 1D set-up. The relations for determining the boundary fluxes $\Gamma e,bplasma$, $\Gamma e,bemiss$, and $\u2009\Gamma i,bplasma$ are given by formulae (24)–(28).

Due to the difference in masses and thermal velocities of electrons and ions, the electron thermal flux is often much higher than the current and the electron-flux-limiting sheaths formed near the electrodes, where the electron and ion densities significantly deviate from the quasi-neutrality assumption. Charge separation results in the voltage drop in the sheath, which can be compared with the total arc voltage. The sheath voltage drop suppresses or accelerates fluxes of charge particles in the sheath, depending on their charge and direction, leading to a positive or negative heat deposition in the sheath. Depending on the direction of the electric current and its density, the wall temperature, and plasma density at the plasma-sheath boundary, the sheath voltage drop can be either positive or negative so that the balance of electric current (22) is satisfied.

We consider two approaches to modeling of the sheaths. In the first approach, we resolve variations in the species densities and the electric field in the sheath. In the second approach, we use the fact that the sheath is typically much narrower than the inter-electrode gap, and can be assumed to be infinitely thin. Then, the effective sheath boundary conditions for fluxes of particles and heat can be applied at the plasma-sheath boundary without resolving the sheath.

#### 1. Effective boundary conditions for collisionless sheath

If the sheath is collisionless (its width is much narrower than the mean free paths of plasma particles), then it is convenient to use effective sheath boundary conditions.^{42,43} The plasma can be modeled using quasi-neutrality assumption (16) instead of using Poisson's equation (15).

Relations given below are in a unified form applicable to both anode and cathode. The sheath voltage drop, $Vsh$, is considered positive if it gives a positive contribution to the total arc voltage. Therefore, the formulas should be used as they are for the cathode, and the sign of the sheath voltage should be reversed for the anode

In the case of positive sheath voltage drop $Vsh$, the ion current is determined by Bohm's condition

For the negative sheath voltage, the ion current is

The electron current from the plasma to the wall is

The electron emission current is

where $jR$ is the current density predicted by Richardson's emission law

$AR$ is the Richardson's constant, $Vw$ is the work function of the electrode material (4.5 V for tungsten), and $ESchott$ is the Schottky correction voltage (about 0.1 V, see Ref. 9, for instance).

Equations (24)–(28) for the fluxes of charged particles in combination with the balance relation for current density (22) allow determining the sheath voltage drop and its sign.

The heat flux to the electron gas at the plasma-electrode boundary can be expressed as

This flux is used as a boundary condition for the electron heat transport equation (12)

The heat flux to the electrode from the gas and the plasma can be expressed as

#### 2. Solving for both plasma and sheath regions

If the sheath is strongly collisional (its width is much higher than the mean free paths of plasma particles), then it can be described by the fluid governing equations (1)–(17). The boundary conditions in this case are as follows. The heat flux to the electron gas is defined by

A zero number density condition can be used for the ion density (see Ref. 13)

The equation of current conservation (22) can, in this case, be used to determine the electron number density at the boundary

### C. 1D Model of heat transfer in the electrodes

To account for the cooling of electrodes by radiation from their sides, in a quasi one-dimensional approximation, we can neglect the temperature variation in the radial direction. The heat flux through the electrode is decreasing due to losses at its side surfaces and increasing due to Joule heating

Here, $\lambda el$ is the thermal conductivity of the electrode material (assumed to be constant, 170 W/m/K for tungsten), $rel$ is the electrode radius, $Tamb=300\u2009K$ is the ambient gas temperature, and $\lambda gas$ is the thermal conductivity of gas surrounding the electrode; the constant value of 0.1 W/m/K for $\lambda gas$ was used; it corresponds to the gas temperature of about 3000 K (close to the temperature of the cathode tip where major heat losses take place). There is no need for a more accurate resolution of the gas thermal conduction term, because it appeared to be of minor importance compared with the radiation term. $Nu$ is the Nusselt number taken to be equal to 1.1 (see Ref. 51), $\sigma $ is the Stefan-Boltzmann constant, $\epsilon $ is the emissivity of the electrode surface, taken to be equal to 1, $\rho el$ is the electrical resistivity of the electrode material (assumed to be constant, small for metallic electrodes).

In the 1D model, Eq. (35) is solved numerically with the heat flux coupled to heat flux from plasma at the arc side minus the surface radiation flux

Here, $qto\u2009electrode$ term is defined in Eq. (31), $qrad,tip$ term takes into account the radiation from the front surface of the electrodes including mutual radiation.^{50} The solution of Eq. (35) is used to define the temperature of the electrode's surface, which is used as a boundary condition in Eq. (19).

## III. NUMERICAL PROCEDURE

As mentioned earlier, a set of governing equations is different depending on a sheath modeling approach utilized. In the case of full sheath and plasma modeling approach, when the sheath and plasma regions are resolved with the same collisional model, one needs to take into account the deviation between the densities of electrons and ions, and treat them as two different variables and solve Poisson's equation (15) throughout the whole computational domain comprising of plasma and sheaths. In the case of effective sheath boundary conditions, the sheath regions are not resolved by the plasma model and quasi-neutrality condition (16) can be utilized. Hence, one can reduce the number of independent variables by exclusion of electron density from the governing equations. Two slightly different solution procedures were utilized depending on the sheath modeling approach.

In general case, Eqs. (1) and (6) for ions can be transformed to one second order differential equation describing the ion transport. From Eq. (1) for ions, using relations between the fluxes of particles (10) and (11)

where $\nu k,j=438kTkj\pi mkjCkj\sigma kjnj\u2009\u2009$ is the effective collision frequency of species *k* with species *j*.

where $Di=k(T+Te)/(\nu i,amAr)$, $V\u2192i=(k\u2207T+kCi(e)\u2207Te+eE\u2192)/(\nu i,amAr)$.

The electric field can be determined from Eq. (1) for the electrons

In the case of full plasma-sheath modeling, differential equations (38), (12), and (13) are solved with the boundary conditions defined in Subsection II B 2 together with closures given by Eqs. (10), (11), (15), (17), and (39) to form a fully closed system of governing equations. Due to non-linearity, differential equations (38), (12), and (13) were solved iteratively in the 1D code; and an implicit scheme was used. At each iteration, the coefficients of these equations were computed using independent variables from the previous iteration, and the source terms were linearized around values of independent variables from the previous iteration. The electric field can be obtained from Eq. (39), the electron density—from Eq. (15). Linearized equations are discretized using second order schemes resulting in tridiagonal matrices solved using tridiagonal matrix algorithm. For numerical stability, the implicit coupling between the electric field and electron number density was used. The electron density in a numerator of the second term on the right-hand side of (39) was expressed via the ion number density and electric field using Poisson's law (15). Therefore, the algebraic equation (39) for the electric field was transformed into the second order differential equation

The Neumann boundary conditions were used for this equation: the gradient of electric field at the boundary is determined from a difference between the ion number density (33) and electron number density (34) using Poisson's law (16).

When using the effective sheath boundary conditions, quasi-neutrality approximation for plasma is used ($ne=ni$) and the electric field can be excluded from the relation for the ion flux resulting in implicit and more numerically stable coupling of the variables. Substitution of Eq. (39) into Eq. (37) gives a relation for the ion flux (a similar approach was used in Ref. 52)

where $D=\mu k(T+Te)$ is the ambipolar diffusion coefficient.

$DT=\mu kT$, $DTe=\mu k(1+Ce(e)+Ci(e))Te$ are the thermal diffusion coefficients

Substitution of the ion flux given by Eq. (41) into Eq. (6) and taking into account that $Ae\u226a1$ yields transport equation for ions

where $V\u2192=\u2212DT\u2207lnT\u2212DTe\u2207lnTe$.

A set of governing Eqs. (42), (12), and (13) with closures (10), (11), and (39) for the electric field is solved iteratively in an implicit manner.

Note that in 1D formulation, all derivatives are taken in primary coordinate along the axis of symmetry of the arc. Diffusion and convection in radial direction are assumed to be negligible. This assumption should be valid because the arc is short (the distance between the electrodes is smaller than their diameters). The radial energy loss from the plasma is taken into account by the radiation term in Eq. (12). The radial energy losses from the side surfaces of electrodes are also taken into account and results show that for long and thin electrodes considered here, these losses are of major importance.

## IV. RESULTS AND DISCUSSION

### A. Benchmarking of 1D code and sheath models

For benchmarking of the code, 1D simulations of the near-cathode region of atmospheric pressure argon arc were performed for the same conditions as in Ref. 13. The current density was 5 × 10^{6} A/m^{2}, and the fixed temperature 3500 K of the cathode surface was assumed. The full plasma-sheath modeling approach was employed; the computational grid was strongly refined near the cathode to resolve variation in the plasma parameters in the sheath. A comparison of our results with the simulations of Ref. 13 is given in Fig. 2; notice the logarithmic scale of x-axis. A very nice agreement on the profiles of various plasma parameters was achieved—the profiles obtained in two simulations almost coincide. The regions of deviation from the ionization equilibrium, thermal equilibrium, and quasi-neutrality are clearly visible and have different thicknesses. The sheath width appeared to be about one micron. This is less or comparable to the mean free paths of all plasma species; hence, the sheath is weakly collisional and the applicability of full plasma-sheath modeling approach is questionable in this case.

Benchmarking of two different sheath modeling approaches: (i) full plasma-sheath fluid modeling; and (ii) instead of sheath modeling, using the effective boundary conditions for collisionless sheath (23)–(31) was also performed. The results of simulations of the near-cathode region with the effective boundary conditions are plotted in Fig. 2 with blue lines.

The results obtained with two sheath modeling approaches are in very good agreement. This agreement is surprising at the first sight because the sheath models compared were supposed to be applicable to the opposite conditions: applicability of the full sheath-plasma model is justified in case of strictly collisional sheath (which is typical for pressures of tens of atmospheres), and the effective boundary conditions were designed for the collisionless sheath. The similarity of the heat fluxes to the cathode surface using these two approaches was already noted in Refs. 13 and 42. The new results show that not only the heat fluxes at the surface, but actually the whole profiles of plasma parameters appeared to be almost the same.

A comparison of performance of the sheath modeling approaches was also performed for the near-anode region. Pressure and current density were the same as for the near-cathode region. The anode temperature was 3900 K. The results are plotted in Fig. 3; notice that the anode is on the right side of the plots, and logarithmic scale is used for the x-axis. Note that the electric potential profile appeared to be non-monotonic due to intensive electron emission from the anode surface: a negative electric field in the quasi-neutrality region and a positive electric field in the sheath region. One can define such a near-electrode region as a double layer. However, despite the complexity of the layer structure, the results obtained with the two different approaches were in very good agreement with each other.

The independence of the results on the sheath modeling approach can be explained by the fact that both approaches capture the major effect determining voltage and current composition in the sheath and correspondingly plasma density at the sheath edge. This effect is the reduction of fluxes of charged species by the Boltzmann exponent due to the sheath voltage. For the boundary conditions (23)–(31), this Boltzmann factor is implemented explicitly. In a full modeling approach, the Boltzmann factor is manifested in an exponential decrease in densities of charged particles with the electric potential profile according to Eq. (1) because the density gradients are the major effect in the sheath regions.

Note that the comparison of sheath modeling approaches was also performed in a wide range of current densities, 2 × 10^{6} A/m^{2}–2 × 10^{7} A/m^{2} (typical for atmospheric pressure argon arc, see, for example, Refs. 37 and 53). In all the considered cases, a good agreement was observed.

According to the results discussed earlier, both approaches were cross-verified against each other and can be claimed applicable for modeling of atmospheric pressure arcs.

### B. Parametric studies of argon arc

#### 1. Discussion of the arc structure: Arc column and near-electrode regions

The plasma density and temperature profiles for arcs of 2.5 mm and 5 mm lengths obtained in the simulations for various current densities are displayed in Figs. 4 and 5. The self-consistent heat transfer between the plasma and the electrodes, 6 mm in diameter and 10 cm in length, was solved. The anodes of the arcs of different lengths are aligned with each other in the figures.

In the middle part of the long arc, ionization and thermal equilibria take place: the plasma density profile is determined by the Saha equation (9), and the temperatures of electrons and heavy particles are equal. Two sub-regions can be distinguished inside of the arc column: the complete equilibrium region where all plasma parameters are uniform (marked with green background in Figs. 4 and 5) and the local equilibrium region where plasma is non-uniform (white background in the Figs. 4 and 5). Near the electrodes, the ionization and thermal equilibria are not maintained; the non-equilibrium regions are marked with blue background in the figures.

As can be seen from Figs. 4 and 5, the near-electrode non-equilibrium regions are rather autonomic: profiles of plasma parameters are almost the same in the near-anode region of short and long arcs at the same current density. This is valid for all current densities considered, except for the lowest one (2.5 × 10^{6} A/m^{2}) at which the near-anode region tends to be longer than the short arc itself. But even in this case, the profiles in the near-electrode non-equilibrium regions for short and long arcs are rather close. In the near-cathode region, profiles of plasma parameters are close in cases of long and short arc for all current densities considered.

Note that the autonomic behavior of the near-electrode regions would be expected in case of very long arc with an area of well-established plasma equilibrium between the near-electrode regions. However, in the presented results, an absolute equilibrium with no variation in the plasma parameters is established at most current densities. Therefore, the plasma parameters in the arc column have rather small influence on the near-electrode regions.

At higher current densities, the equilibrium region occupies a significant part of the long arc, but is not observed in the short arc at any of the current densities considered. With a decrease in the current density, the non-equilibrium and local equilibrium regions significantly elongate, especially near the anode, reaching several millimeters in width and resulting in the vanishing of the equilibrium region in the 5 mm-long arc.

This transformation results in a qualitative change in the electric potential profile in the long arc (see Fig. 6). The plasma equilibrium region in the central part of the arc is characterized by a constant electric field and linear growth of the electric potential (see Fig. 6). This region is commonly referred as a positive arc column. With an increase in the current density, the arc column becomes more pronounced: it elongates due to the narrowing of the near-electrode regions, and the electric field in the column becomes stronger. The positive column cannot be distinguished at the lowest current density considered (2.5 × 10^{6} A/m^{2}).

In the near-cathode region, the electric potential significantly jumps from a reference zero value at the wall. A major part of the arc voltage is typically gained in this region; however, in the case of higher current densities, the role of the arc column increases and can become dominating at a larger gap size. Near the anode, the temperature and plasma concentration decrease, causing electron diffusion toward the anode. It results in a decrease in the electric field in order to preserve a constant electron current density. In order to describe this variation in the electric field, it is convenient to express the electron current density as a composition of electron mobility in the electric field, electron diffusion, and thermal diffusion. This relation can be obtained from Eq. (39) neglecting the last term

Here, $\sigma $ is the electrical conductivity (not to confuse with cross-sections).

The contributions of different electron current density components are displayed in Fig. 7 for the arc 2.5 mm length and current density 5 × 10^{6} A/m^{2}. As mentioned earlier, no equilibrium region is present in the short arc. As seen from Fig. 7, the diffusion component is the largest throughout the arc, and the electric field does not play a dominant role anywhere.

Near the electrodes, at distances about 0.5 mm, due to significant gradients of plasma density, the magnitude of the diffusion component becomes several times higher than the total current density. It results in a significant increase in the absolute value of the electric field in these areas in order to keep the total current constant. Near the cathode, the diffusion component is negative, the electric field is positive, and vice versa near the anode. Therefore, significant heat source and heat sink are located in the near-cathode and the near-anode plasma, respectively.

The thermal diffusion plays a significant role throughout the arc. It is comparable to other electric field components and should be taken into account in the accurate modeling; however, it is typically several times smaller than the diffusion component and does not lead to any qualitative effects. Note that the same conclusions are valid for longer arcs and other current densities not displayed in Fig. 7.

Note that the electron–ion collision frequency is much higher than the electron–atom one in most of the arc [$\nu e,i\u226b\nu e,a$ in the denominator of parameter $\sigma $, Eq. (43)]. The ratio of the collision frequencies, or parameter $P$ defined in Eq. (2), is about $(\sigma eana)/(\sigma eine)$ or $(\sigma ea/\sigma ei)(1\u2212\alpha i)/\alpha i$ where $\alpha i=ne/(ne+na)$ is the ionization degree. The ratio of cross-sections $\sigma ea/\sigma ei$ is about $0.01$ at temperatures below 18 000 K. According to Fig. 8, the ionization degree is much higher than $0.01$ almost everywhere in the arc for all current densities considered. Hence, the parameter $P$ is small compared with unity throughout the arc, including major parts of the near-electrode regions. This statement will be used in the analytical arc model presented in the accompanying paper.^{44}

#### 2. Discussion of the near-cathode region

A significant electric potential jump near the cathode (see Fig. 6) results in heating the electrons, as can be seen in Fig. 4. On the contrary, the temperature of heavy particles decreases toward the cathode to reach the temperature of the electrode at its surface. Provided the electron density is in ionization equilibrium (n_{e, Saha}, dotted lines in Fig. 5), an increase in the electron temperature near the cathode would lead to an increase in the plasma density, whereas the actual plasma density decreases due to non-equilibrium fast acceleration of ions onto the cathode surface as described by Bohm's condition (24) for positive sheath. The deviation between the actual and equilibrium number densities exceeds an order of magnitude and results in high net production of ions that move toward the cathode due to electric field and density gradient. The ion current in the near-cathode region grows from zero value outside the region to about 15% to 25% of the arc current at the cathode surface (see Fig. 9). Therefore, according to Eq. (22), the electric current at the cathode surface is mostly carried by emitted electrons [75%–85%; the current of electrons from the plasma is negligible because it is suppressed by the voltage drop in the sheath; see Eq. (26)]. This qualitative picture is in accordance with predictions of the previous theoretical models of the cathodic region;^{7,8} however, a quantitative disagreement in the values of ion current ratio is observed with these studies. An opposite trend for ion current fraction was observed in simulations^{13} where a fixed cathode temperature was used. Therefore, it is important to take the heat transfer in the electrode into account in order to accurately describe the near-cathode plasma.

Electron emission requires certain temperature of the electrode surface, as specified by Eq. (28). The cathode temperature slightly goes up with an increase in the total current density (see Fig. 10), but in general, it remains in the vicinity of 3500 K. The temperature of the electrode is determined by the energy balance between the electrode and the arc plasma; therefore, heat transfer through the electrode should be taken into account, and interfacial conditions should be used.

Heat fluxes in the plasma transferred by electrons and heavy particles are plotted in Figs. 11(a) and 11(b) correspondingly; positive values correspond to the heat fluxes directed away from the cathode (having positive projection to the x-axis). The conductive heat fluxes of the electrons and heavy particles correspond to the first terms on the right-hand sides of Eqs. (12) and (13) correspondingly. The electron convective heat flux corresponds to the left-hand side of Eq. (12). Ions carry their ionization energy; the corresponding heat flux, $\Gamma iEi$, represents the convective heat flux transferred by heavy particles. As mentioned earlier, a heavy component as a whole does not move [see Eq. (10)]; therefore, its thermal energy does not contribute to the convective heat transfer.

As seen from Fig. 11(a), the electron convective heat flux is significant throughout the arc. Electrons gain significant thermal energy in the vicinity of the cathode and transfer it as they drift toward the anode. The conductive heat transferred by the electrons is significant only at lower current densities and in the near-cathode region, where electron temperature gradients are high.

It should be noted that electrons gain a significant amount of energy in the cathodic space-charge sheath: as can be seen from Fig. 12, about 2/3 of the near-cathode voltage drop is attributed to the sheath. Electrons emitted from the cathode's surface accelerate significantly in the weakly collisional sheath and then exchange the energy in collisions with plasma species, mostly with plasma electrons rather than with heavy particles, due to difference in masses. High heat deposition in a very thin near-electrode region manifests in a sharp increase in the electron temperature at the cathode's surface (see Fig. 4).

Generally speaking, there is a region in the vicinity of the cathode where the emitted electrons are not thermalized with plasma electrons. However, the thickness of this region is about several electron mean free paths, i.e., several microns, much smaller than the length scale at which variation in the electron temperature occurs. Therefore, the description of the near-cathode plasma with fluid model should be valid.

Hot plasma electrons in the near-cathode region transfer a significant portion of their thermal energy to heavy particles in elastic and inelastic collisions corresponding to ionization and heating; see the last two terms on the right-hand side of Eq. (12). Heavy particles bring this energy back to the cathode [see Fig. 11(b)] via thermal conduction and convection. These two constitutive parts of the heat flux are comparable in the near-cathode region and are both small outside the region; for illustrative purposes, a logarithmic scale was used for the x-axis of Fig. 11(b). As mentioned earlier, convective heat flux to the cathode is associated with the ionization energy ions bring to the cathode during recombination at the surface.

In the previous theoretical studies,^{7–9} conductive heat flux to the cathode was neglected, leading to overestimation of the ion current fraction. This question is addressed in more details in the accompanying paper.^{44}

The statement of low conductive heat transfer at the plasma boundary of the near-cathode region was previously used in a theoretical study,^{9} but in that paper, it was used as an assumption, without a proof. The validity of this assumption was especially questionable for short arcs. In the present study, we provide a justification based on the results of simulations of the whole arc. This simplification will be used in the arc model presented in the accompanying paper.^{44} Asymptotic relations for the near-cathode region voltage drop, electron temperature, ion current, and size of the region are presented in the accompanying paper.^{44}

#### 3. Discussion of the near-anode region

The temperature of heavy particles decreases and reaches the temperature of the anode at the electrode surface. The electron temperature is close to the temperature of heavy particles due to collisional heat exchange [term $Qe\u2212h$ in Eqs. (12) and (13)]. Because the net electron flux is directed from plasma toward the anode, electron gas is losing energy at the boundary. It is easy to explain by rewriting Eq. (29) for the electron heat flux at the anode boundary in the following form:

where $\Gamma e,b=\Gamma e,bplasma\u2212\Gamma e,bemiss$ is net electron flux to the electrode. Apparently, the first term on the right-hand side of Eq. (44) is negative. The second term is also negative if the electron temperature at the boundary is higher than the temperature of the electrode.

Due to negative heat flux of the electron gas, the deviation between the temperatures of electrons and heavy particles is not strong, unlike in the near-cathode region. The electron temperature generally follows the downtrend of heavy particles temperature in the near-anode region (see Fig. 4), though it does not reach the temperature of the electrode; the simulations predicted the electron temperature at the anode surface of about 5500 K, weakly depending on the current density.

Electrons bring energy to the anode and heat it up. The energy flux brought to the anode by the electrons can be expressed from Eq. (31)

Here, $je,a=e\Gamma e,a$ is the net electron current density to the anode. It is very close to the total current density, because the ion current at the anode is small. The largest term in Eq. (45) is the first term on the right-hand side due to a large value of work function compared with the temperature.^{27} The convective heat flux is not high due to relatively low temperature; the sheath voltage drop typically does not exceed 1 V, as will be shown below. Therefore, the heat flux to the anode is roughly proportional to current density, and as a result, the temperature of the anode tip significantly increases with the current density (see Fig. 10).

The low electron temperature near the anode results in a significant decrease in the equilibrium plasma density (n_{e, Saha}; see Fig. 5). The actual plasma density also decreases, though it is higher than the equilibrium density (a net recombination of ions takes place in the near-anode region). Due to a significant decrease in the plasma density, electron diffusion plays an important role in the near-anode region, driving electrons toward the anode. It results in a strongly negative electric field (see Figs. 6 and 7) and in negative voltage in the near-anode layer.

Electric potential profiles in the vicinity of the anode are plotted in Fig. 13 for two different current densities and for anodes with and without water cooling. In case of cooled anodes, two temperatures of the anode surface were considered: 1 000 K and 3 000 K. For both temperatures, the electron emission from the anode surface [Eq. (28)] is negligible and, as can be seen from the figure, the potential decreases fast in the vicinity of the anode, and the sheath voltage drop is negative to suppress the excess of electrons from plasma moving to the anode. Anode without additional cooling reaches temperatures above 4000 K at the current densities considered. The electron emission is very high at such conditions, and the sheath voltage drop is positive in order to suppress the emission.

In Fig. 14, the voltage drop in the anode sheath and whole near-anode non-equilibrium layer are plotted as functions of current density. Anodes with and without water cooling are considered. Because the voltages do not depend significantly on a particular value of the anode temperature in the case of cooled anode, the results are plotted for only one temperature (1000 K).

As can be seen from Fig. 14, in the case of cold anode, the voltage drop in the sheath is negative for all current densities considered. Voltages in the space-charge sheath and in the non-equilibrium layer as a whole decrease with current density (the absolute values of the voltages increase), whereas in the case of the hot anode, the sheath voltage drop is positive and increases up to 0.7 V; the voltage in the non-equilibrium layer increases only at low current densities, but generally a positive trend is observed, which is in qualitative agreement with anode voltage measurements^{54} in a carbon arc with the hot anode. It leads to a conclusion that in the case of the cold anode, arc constriction in the near-anode layer is energetically advantageous contrary to the hot anode case.

Accurate quantitative asymptotic relations for voltage and temperature in the non-equilibrium region and its size can be found in the second paper of the series.^{44}

#### 4. Significance of thermal and ionization non-equilibrium effects

Additional computational runs were performed to study the effects of thermal and ionization non-equilibrium on integral arc characteristics and profiles of plasma parameters. Local thermal equilibrium (LTE) was assumed in these computational runs: the electron temperature was taken equal to the temperature of heavy particles. One energy balance equation was solved instead of Eqs. (12) and (13)

Boundary condition (19) was used with a simple relation for heat flux to the electrode

In addition to LTE assumption, the local ionization equilibrium (also often referred as chemical equilibrium, LCE) was assumed in one of the computational runs: Saha equation (9) was used to determine plasma density instead of the ion transport equation (42); ion production rate $si$ in the ionization heat term $Qion$ was calculated using Eq. (6). Effective collision-less sheath boundary conditions were used for quasi-neutral plasma equations.

In Fig. 15, results of non-equilibrium simulations for 2.5 mm arc at 5 × 10^{6} A/m^{2} current density (dark lines) are compared with the results obtained with an LTE assumption (red lines) and both LTE + LCE assumptions (green lines).

As can be seen from the figure, in the non-equilibrium case and the LTE case, the voltages in the near-cathode region and in the whole arc as well as temperatures of the electrodes are almost the same, but the plasma profiles are significantly different, especially in the near-cathode region where significant deviation between temperatures of electrons and heavy particles should take place. The LTE simulations give lower electron temperature in the near-cathode region resulting in significantly less steep plasma density and electric potential profiles and prediction of drastically low value for the ion current.

In the “LTE + LCE” case, the plasma density near the electrodes is strongly underestimated due to applying the Saha equation at low temperatures. As a result, higher voltages in the near-electrode regions are required to heat up and maintain the plasma: about 5 V extra in each near-electrode region, i.e., 10 V in the arc total. Correspondingly, significantly higher temperatures of the electrodes are observed.

### C. Validation against experimental data

The arc model was validated against experimental data of Refs. 40 and 41 in which the atmospheric pressure argon arc with cylindrical tungsten electrodes 3 mm in diameter was run at arc currents of 30 A, 50 A, and 100 A. The inter-electrode gap width was varied from 0.3 mm to 3.5 mm.

In Fig. 16, the arc voltage is plotted as a function of current density. Results of non-equilibrium 1D simulations with self-consistent heat transfer in the plasma and the electrodes are compared with the experimental data. Reasonable qualitative and quantitative agreement is observed. At larger inter-electrode gaps, the arc voltage linearly increases with the gap size. This behavior can be explained by the elongation of equilibrium region of the arc column (see Figs. 4–6). At smaller gaps (below 0.5–2 mm, depending on the arc current), the near-electrode non-equilibrium regions overlap and the dependence of the arc voltage on current is different. Simulations show that the arc voltage decreases with current if the arc width is comparable to the non-uniform near-electrode layers. The same non-monotonic behavior is observed in the experimental data as well.

## V. CONCLUSIONS

The 1D model of argon arc with cylindrical tungsten electrodes was implemented in a numerical code. The model features coupled heat transfer in the plasma and electrodes and an accurate account of non-equilibrium effects in plasma including thermal, ionization non-equilibrium, electron diffusion, thermal diffusion, and effects of space-charge sheaths. The model was benchmarked against previous simulations^{13} and validated against experimental data.^{40} Good agreement was obtained.

The role of non-equilibrium plasma effects in arc self-organization was studied. It was shown that the electron diffusion significantly affects the electric field and leads to its reversal near the anode. The thermal diffusion of electrons also plays a significant role. However, its effect is smaller than that of classical diffusion.

The thermal non-equilibrium effect ($Te\u2260Ta$) is important, especially in the near-cathode region where a significant deviation between the temperatures of electrons and heavy particles takes place (see Fig. 4). Disregarding the thermal non-equilibrium yields a low electron temperature in the near-cathode region affecting plasma density and electric potential profiles and resulting in the prediction of drastically low value for the ion current (see Fig. 15). Similarly, assuming thermal and ionization equilibria results in the prediction of drastically low value for the plasma density near the electrodes. Correspondingly, a higher voltage is required to maintain the arc plasma under the assumption of full thermal and ionization equilibriums: additional 5 V in each near-electrode region, i.e., 10 V in the total arc voltage. Moreover, the temperatures of the electrodes are overestimated as well.

Two space-charge sheath modeling approaches were implemented in the code: (i) in the first approach, the sheath is assumed collisional and variations in the species densities in the sheath are resolved directly using the fluid model,^{13} and the Poisson equation is solved for the electric field; (ii) in the second approach, the sheath is assumed collisionless and effective sheath boundary conditions for fluxes of heat and charged particles^{42,43} are applied at the plasma-sheath boundary without resolving the sheath, and quasi-neutrality assumption is utilized for the plasma.

Both collisional and collisionless space-charge sheath modeling approaches yielded the same plasma profiles for the near-electrode regions and entire arc (see Figs. 2 and 3). Independence of the results on the sheath modeling approach can be explained by the fact that both approaches capture the major effect determining voltage and current composition in the sheath and correspondingly plasma density at the sheath edge. This effect can be explained by the fact that the electron density profile obeys the Boltzmann law due to high gradients of the electron density in the anode and cathode regions. For the collisionless boundary conditions (23)–(31), this Boltzmann factor is implemented explicitly. In addition, ionization in the thin sheath regions does not contribute to the ion flux toward electrodes.

Parametric studies of short atmospheric pressure argon arc with tungsten electrodes were performed for various current densities and inter-electrode gap sizes (see Figs. 4 and 5). It was shown that the near-electrode regions can become up to several millimeters long, significantly affecting the arc characteristics. The near-anode region is typically wider and its size more strongly depends on the current density.

A major part of the arc voltage typically falls in the near-cathode region; however, the voltage drop in the arc column (where the Joule heating is locally balanced by radiation losses) increases with arc current and may become a considerable part of the arc voltage (see Fig. 6). The voltage drop in the cathode sheath constitutes a major fraction of the voltage in the region in the near-cathode region. The cathode temperature is about 3500 K and does not change significantly with current. A significant part of the heat flux to the cathode is due to the thermal conduction of heavy particles.

The ion current fraction at the cathode is about 20% of the total current and decreases with an increase in the total current density. This is in qualitative agreement with predictions of the previous models of the cathode region.^{7,8}

The Volt-Ampere characteristic (VAC) of the near-anode region depends on the anode cooling mechanism (see Fig. 14). The anode voltage is negative. In the case of strong anode cooling (water-cooled anode) when the anode temperature is sufficiently low so that the thermionic emission from the anode is negligible, the anode voltage decreases with the current density (the absolute value of the negative anode voltage drop increases) due to stronger gradients of plasma density and temperature (see Figs. 4 and 5) at higher current densities. Falling VAC of the near-anode region suggests the arc constriction near the anode. Without anode cooling, the anode temperature increases significantly with the current density, leading to a drastic increase in the thermionic emission current from the anode. Correspondingly, the potential barrier near the anode forms to limit the thermionic emission, and the total anode voltage increases with the current density—the opposite trend in the VAC is observed as compared with the case of strong cooled anode.

## ACKNOWLEDGMENTS

The authors are grateful to Vlad Vekselman (PPPL, NJ), Yevgeny Raitses (PPPL, NJ), Mikhail Shneider (Princeton University, NJ), Nelson Almeida (Universidade da Madeira, Portugal), Mikhail Benilov (Universidade da Madeira), Ken Hara (Texas A&M University, TX), and Marina Lisnyak (Université d'Orléans, France) for fruitful discussions and valuable input.

The research was funded by the U.S. Department of Energy (DOE), Office of Fusion Energy Sciences.