Electron transpiration cooling for the leading edges (LE) of hypersonic aircraft utilizes thermionic emission; however, space-charge effects limit the electron emission rate, potentially diminishing the efficiency of this cooling mechanism. We develop a variational weak form of the Poisson equation that describes the sheath potential and then numerically solve it using the finite element method. This formulation has two main benefits: (1) the space-charge limit condition can be incorporated as a constraint and (2) it allows for the analysis of three-dimensional geometries with complex boundary conditions. We demonstrate that the current emitted from the surface of an LE is generally a small fraction of the Child–Langmuir limit due to space charge. We then propose several methods to enhance the emitted current from the surface and to boost the cooling effect of thermionic emission. These include increasing the plasma density, applying a negative surface potential, and using fringe fields under suitable geometric conditions. For a LaB_{6} emitting LE, the total emitted current is shown to be minimal and independent of the temperature of a surface with floating potential. However, when a negative potential is applied and the surface is heated, the emitted current follows the Richardson–Dushman relationship up to a critical temperature, beyond which it remains constant. At an applied surface potential of −5 V, the critical temperature is around 1700 K.

## I. INTRODUCTION

^{1}and Air Force Research Laboratory X-51A.

^{2}However, the widespread adoption of these vehicles has been constrained by conflicting design criteria. One such criterion is the preference for a small radius of the leading edge (LE), represented as $RLE$, to achieve an improved lift-to-drag ratio. However, this design decision raises concerns regarding thermal management. The surface heat flux $q\u2033$ at the leading edge can pose challenges, as described by the following relationship:

^{3}

Hypersonic vehicles have used various thermal protection systems, which can be categorized as passive, semi-passive, and active operation.^{4} In initial designs, ultrahigh temperature ceramics have been used extensively.^{5–7} Passive systems employ insulating materials to restrict temperature rise, heat sinks to conduct heat away from high heat flux regions, or hot structures made of high melting point materials.^{8–11} Semi-passive cooling techniques include embedded heat pipes that improve effective conductance^{12–15} and ablative materials that remove heat through thermochemical decomposition.^{16–18} On the other hand, active cooling systems use pumped coolants that are injected directly into the flow field or in a closed flow loop to cool the internal structure.^{19–23} However, semi-passive and active cooling methods are generally more expensive and complex compared to passive systems, and they also increase the size and weight of the vehicle.

Electron transpiration cooling (ETC) is an innovative thermal management technique that uses thermionic emission to cool extremely hot surfaces.^{25–27} When a suitable material is subjected to high temperatures, some of its electrons gain enough energy to overcome the material's work function^{27} and escape. This electron emission carries away significant amounts of heat energy, actively cooling the surface. ETC can potentially be a passive and highly effective cooling system for the leading edges of hypersonic vehicles, although challenges remain to finding optimal materials and designing practical and robust systems. In some advanced concepts, emitted electrons can be collected downstream on a cooler surface. Recombination at this surface releases heat but in a more manageable location. ETC has attracted considerable interest, leading to a multitude of recent computational and analytical investigations.^{27–31} However, experimental research in this domain is still in its infancy, with several recent publications that address this topic.^{32–37} Initial computational fluid dynamics studies emphasized the importance of space charges resulting from emitted electrons in predicting the cooling efficacy of ETC.^{27} More recent modeling efforts have integrated cesium seeding to alleviate space charge buildup and enhance electron emission, thereby refining ETC performance.^{31} Easily ionized cesium vapor has been shown to decrease the space charge in thermionic energy converters.^{38–40}

In the more common method of cooling via ablative structures, composites, resins, and ceramics are used in reentry vehicles and missiles. The main advantage is the ability to handle extreme heat fluxes effectively. However, the main disadvantage is that the ablated material is consumed, and thus, the mechanism works for a limited duration. In addition, the aerodynamic drag and in-flight stability can be greatly altered by changes in the geometry due to ablation. On the other hand, ETC relies on electron emission from heated surfaces by using metals with low work functions (e.g., tungsten). The main advantage of this mechanism is its ability to provide continuous passive cooling without material loss over longer lifespans. It is suitable for applications requiring precise geometry of the vehicle during flight; however, it involves operation at very high temperatures and, as we will see later, is greatly limited by space-charge accumulation close to the emitting surface.

The concept of space charge and its effect on electron emission is a subject of great interest. The classic result of Child^{41} and Langmuir^{42} gives the maximum current density of an emitter in a vacuum and has been understood by plasma physicists for more than a century. Extensions accounting for a finite velocity of emitted electrons,^{43} a Maxwellian distribution of emitted electrons,^{42,44} relativity,^{45} and quantum effects^{46} would follow. The emissive character of a surface is altered in the presence of a quasineutral plasma. Naturally, the potential very close to the surface is modified when plasma ions and electrons are present, forming a Langmuir sheath.^{47} Electrons emitted from the surface result in a modified space charge solution with the flow of ions and electrons emitted analogously to a double sheath.^{47,48}

Although there are several models to describe the electrostatic potential in this situation,^{49,50} analytical solutions are difficult to obtain without assumptions and are limited to simple geometries. More recently, an analysis of the transition between temperature-limited and space-charge-limited electron flow for a planar diode with a thermionic cathode that has spatially nonuniform emission properties was presented by Chernin *et al.*^{51} and Jassem *et al.*^{52} Their key finding was that the shape of the Miram or “rollover” curve is significantly affected by non-uniformities in the space-charge density. Chernin *et al.* analyzed the steady-state limiting current in an infinite periodic array of thin electron sheets spaced by period, *p*, in a planar diode of gap voltage, *V*, and gap separation, *d.*^{53} Using an iterative solution of a non-linear integral equation,^{53,54} they determined that the average anode current density is bounded by the classical 1D Child–Langmuir limit.

Here, we present a method to determine the potential distribution directly by solving the weak form of Poisson's equation using the finite element method (FEM). Space-charge constraints are enforced at the boundary using a constrained variational technique. First, we validate the method by applying this technique to the classic Child–Langmuir (CL) solution^{41,42} in 1D and then extending it to a numerical 2D model of an electron beam that emits electrons in vacuum.^{55–57} We will then apply the technique to various plasma conditions relevant to hypersonic flight and determine plasma parameters and the electrostatic potential. The goal is to understand the influence of space charge on the limiting of the electron current at the LE surface.

The constrained variational formulation of the Poisson equation is derived in Sec. II. The FEM is then used to solve the problem numerically, and its precision is confirmed by comparing it with the analytical solution of the 1D CL equation and the existing 2D numerical results in Sec. III. The results of the model under different conditions of the LE, with and without surface electron emission, are detailed in Sec. IV. Finally, a summary of the study and the conclusions of key findings are presented in Sec. V.

## II. VARIATIONAL FORM OF POISSON'S EQUATION

*ρ*represents the charge density,

*Z*is the charge per ion,

*n*is the ion density,

_{i}*n*is the electron density,

_{e}*ε*

_{0}is the permittivity of free space (8.854 × 10

^{−12}F m

^{−1}), $e\u22481.602\xd710\u221219\u2009C$, the elementary charge, and Δ is the Laplacian operator.

^{58}Poisson's equation can be reformulated as

*trial*function, $\varphi \u0303$ (sometimes called “test”), and integrating over the volume, Ω, gives the following result:

*dA*and the unit normal

**. The test or trial function is selected from a linear function space**

*n**V*

_{0}that contains functions of the type $\varphi \u0303(r)$ and $\lambda \u0303(r)$, which vanish on the domain boundary,

*unconstrained weak form*of Poisson's equation,

*shape functions*within the elements that make up the domain functions $\varphi $ and $\varphi \u0303$.

*action*integral, $S$, for fields (and not single particles), is obtained by integration over the four-vector spacetime of the element $dx4=dtdx1dx2dx3$. Since the problem is independent of time, the action integral is then given by integration over three-vector space only $d\Omega =dx1dx2dx3$,

*λ*can still be a full field variable. The entire domain form of the augmented action integral is then

*λ*can now be obtained if we minimize the function $\Psi (\epsilon 1,\epsilon 2)$ w.r.t the two parameters $\epsilon 1,\epsilon 2$. More precisely, the Gateaux derivatives of the functional $F$ vanish when $F$ is minimum,

The weak form of second order (PDE) offers several significant advantages for numerical solutions, mainly due to how it handles gradients, boundary conditions, and function spaces. The order of derivatives applied to the solution function is reduced compared to the strong form. This allows for the use of lower-order basis functions in numerical methods such as the FEM. It naturally incorporates certain types of boundary conditions (especially Neumann) into the formulation through integration by parts, making it easier to enforce these conditions within numerical methods. The weak form is particularly advantageous for problems defined in complex or irregular domains. For example, it allows the use of unstructured meshes, which can more accurately represent complex geometries. The weak form allows for the discretization of the domain into elements and the approximation of the solution with basis functions in a FEM framework. It can also offer better stability and convergence properties for numerical solutions because the choice of function spaces and basis functions can be tailored to ensure these properties.

## III. MODEL VALIDATION

In this section, we assess the accuracy of the variational solution for the sheath potential by comparing it with the existing analytical and numerical solutions for standard problems.

### A. Child–Langmuir

^{41,42}is first considered for two infinite parallel plates separated by a distance

*d*, as shown in Fig. 1(a). An electron current is emitted from the cathode at

*x*= 0 where the potential is $\varphi =0$. A constant potential,

*V*, is applied to the anode at

*x*=

*d*. The maximum current emitted from the cathode is achieved when the electric field at

*x*= 0 is equal to zero (that is, $\u2207\varphi =0$). Poisson's equation, Eq. (2), describes the potential distribution within the volume subject to the boundary conditions: $\varphi (0)=\u2207\varphi (0)=0\u2009\u2009and\u2009\u2009\varphi (d)=V$. An electron in position

*x*has potential $\varphi $ and its kinetic energy is $12mev2=e\varphi $, where

*m*is the electron mass and

_{e}*v*is the velocity. The emitted current density is $J=enev$, where the magnitude of the elementary charge is used, and

*n*= 0 because there are only electrons in this case. Solving for

_{i}*n*and then plugging into Eq. (2) gives

_{e}*J*, is analytically available and is known as the CL current,

*J*,

_{CL}^{41,42}

*x*= 0 as a space charge constraint at the boundary $\u2202\Omega $,

*E*is a constant surface electric field (

_{s}*E*= 0 under space-charge-limited conditions) and

_{s}**is the unit surface normal. The associated null function is $g(x,\varphi ,\u2207\varphi )=n\xb7\u2207\varphi +Es=0$. Examining the constraint**

*n**g*, the second term in Eq. (16) disappears since $\u2202g\u2202\varphi =0$. Differentiating the functional $S$ gives $\u2202S\u2202\varphi =\u2212\rho $ and $\u2202S\u2202\u2207\varphi =\epsilon \u2207\varphi $. Inserting in Eq. (16), along with the expression for

*Q*, and where $\kappa =1\epsilon 0me2e$ leads to

Figure 1(a) shows a comparison between the analytical and constrained weak-form solution. The normalized potential and electric field agree well with the root mean square errors (RMS) of 0.9% and 3%, respectively. The higher RMS error for the electric field is caused by the large gradient in the electric field near the cathode, *x* = 0, where the boundary constraint is applied. Figure 1(b) shows details of the solution near the cathode region for various numbers of elements. It is clear that the solution for the potential field is more accurate than the electric field, which is expected because of the high electric field gradients near the cathode. What is interesting here is the rapid convergence of the electric field solution for a grid size less than $10\u22122$, which corresponds to only 100 elements. More complex geometry requires a much smaller grid size, as we will see later.

### B. Two-dimensional Child–Langmuir

^{55–57}featuring two parallel plates of length

*L*, separated by a distance

*d*, as shown in Fig. 2. The cathode at

*y*= 0 has a potential $\varphi =0$, and a constant potential $\varphi =V$ is applied to the anode at

*y*=

*d*. Electrons are emitted from a strip of width

*w*on the cathode. Outside of the emitting surface at the edge of the parallel plates ( $|x|=L/2$), the electric field corresponds to a vacuum field. Again, the maximum emitted current density is achieved when the electric field normal to the surface throughout the width strip

*w*equals zero. The problem is described by the following equations:

^{55,56,59}Let us rewrite the unknown current density across the width strip

*w*as $J(x)=\lambda (x)JCL$, where the Lagrange multiplier, $\lambda (x)$, now varies with

*x*. The right-hand side of Eq. (22) is then rewritten as

The results of the variational model are almost identical to the solutions in Ref. 57. The Lagrange multiplier, *λ*, is plotted in Fig. 3(a) for different widths of the emitting strip, *w*, and two different plate separation distances, *d*. The width of the emitted strip ranges from *w* = 0.5 to 6 cm. For the smallest width, *w* = 0.5 cm, *λ* > 1 in the center of the beam, representing an increase in the current emitted relative to the CL limit. An enhancement is observed for a larger plate separation distance. As the width of the emitting strip is increased, this enhancement is less pronounced at the beam center. However, electrons are able to escape into the vacuum region at the edge of the emitting strip, resulting in an increase in the emitted current. Figure 3(b) shows the normalized electric field on the cathode (*y* = 0). Across the emitting strip, this value is constrained because of the space-charge condition, but it increases rapidly to the vacuum value near the edge. Figure 4 illustrates the normalized potential distribution for different widths of the emitting strip and two separation distances between the plates. At the edges of the emitting strips, the potential distribution forms steep gradients, leading to the significant increase and direction change in the electric field depicted in Fig. 3(b).

The spatial distribution of the normalized potential between the cathode and the anode for a plate length, *L* = 8 cm, separated by distances of *d* = 1 and 2 cm, and beam widths of *w* = 0.5, 2, and 6 cm, are shown in Fig. 4. It is noted that the spatial distortion of the potential field in the space-charge region is close to the cathode. The normalized potential (solid) and the electric field (dashed) along the centerline (*x* = 0) of the 2D geometry are plotted in Fig. 5(a) for various emitting strip widths and a plate separation distance, *d* = 1 cm. The result resembles the CL solution in Fig. 1(b). The peak potential and electric field occur for the smallest width of the emitting strip (*w* = 0.5 cm), although the difference in potential is very small. The difference in the normalized electric field was greater when the separation distance was increased to 2 cm. The Lagrange multiplier, *λ*, along the centerline (*x* = 0) changes when the width of the emitting strip is kept constant and only the plate separation distance is changed. As noted in Ref. 57, *λ* decreases as the ratio between the width of the emitting strip and the separation distance, *w*/*d*, increases. Figure 5(b) shows the normalized total current, $I/ICL$, on the emitting strip as a function of this ratio for different widths of the emitting strip. The solid lines are the results of our numerical simulations, and the markers are data from Watrous *et al.*^{56}

The enhancement of the emitted current over the CL value, corresponding to $I/ICL$ > 1, is more pronounced when both the widths of the emitting strips and the ratio *w*/*d* are decreased. As *w* increases, this enhancement becomes less pronounced, and $I/ICL$ approaches a constant value of 1. Although these results are in agreement with those of Ref. 56 they reveal the explicit dependence of the current on the width of the strip itself, *w*, as can be seen in Fig. 5(b). This particular aspect was not examined in the study conducted by Watrous *et al.*^{56} Their investigation concluded that the current is solely a function of a single parameter (*w*/*d*), a conclusion that our findings do not corroborate. The agreement with the results shown in Fig. 13 of Ref. 56 is notable within the specific conditions they analyzed; however, their analysis was restricted to particular ranges of *w*/*d* for each respective value of *w*.

## IV. HYPERSONIC LEADING EDGE

The density of electrons and ions in a quasi-neutral plasma is such that negligible electric fields are present, thus maintaining the potential distribution in the bulk of the plasma. When an electrode (e.g., the surface of a hypersonic vehicle) is inserted, the potential distribution near the surface is disturbed by the random motion of electrons and ions that impinge on the surface. The flux of electrons to the surface exceeds the flux of ions because of their higher energy and greater mobility. A negative charge develops on the surface, while corrective fields allow the potential distribution in the bulk plasma to remain undisturbed.^{60} Plasma electrons in thermal equilibrium have a Boltzmann distribution with local potential: $ne=n0\u2009exp(e\varphi Te)$, where *n*_{0} is the plasma density far from the surface and *T _{e}* is the plasma electron temperature (eV).

*floating*potential,

^{61}A similar relationship was derived by Hobbs and Wesson

^{49}when electrons are emitted from a floating electrode,

The plasma conditions used in our analysis are based on the experimental data from Kang *et al.*^{62} and are listed in Table I. Although not explicitly mentioned, the altitude and flight conditions in their research correspond to Mach 26–27 flight conditions. The speed of sound, *a*, was calculated using the expression $a=\gamma RT\u221e$, where *γ* and *R* are the specific heat ratio (*γ* = 1.4) and the gas constant for air (*R* = 287 J kg^{−1} K^{−1}). The free stream temperature, $T\u221e$, at different altitudes is found in the US Standard Atmosphere.^{63} The electron temperature, *T _{e}*, and the electron density (i.e., plasma), which we use to set the value of

*n*

_{0}, were based on measurements of the electrostatic probe taken at the junction between a hemispherical nose of a radius of 6 in. (152.4 mm) and the 9° half-angle cone that followed.

^{62}The ion mass corresponds to an oxygen plasma. The electron temperature and plasma density are based on the measurements of Kang

*et al.*

^{62}

Altitude . | Speed of sound, . | Mach number, . | Electron temp. . | Plasma density, . | Debye length, . |
---|---|---|---|---|---|

(km) . | a (m s^{−1})
. | $M\u221e$ . | T (eV)
. _{e} | n_{0} (m^{−3})
. | λ (_{D}μm)
. |

84 | 278 | 27.6 | 0.43 | 1 × 10^{15} | 154 |

81 | 281 | 27.2 | 0.43 | 5 × 10^{15} | 69 |

76 | 288 | 26.6 | 0.52 | 1 × 10^{16} | 54 |

71 | 295 | 25.9 | 0.60 | 9 × 10^{16} | 19 |

Altitude . | Speed of sound, . | Mach number, . | Electron temp. . | Plasma density, . | Debye length, . |
---|---|---|---|---|---|

(km) . | a (m s^{−1})
. | $M\u221e$ . | T (eV)
. _{e} | n_{0} (m^{−3})
. | λ (_{D}μm)
. |

84 | 278 | 27.6 | 0.43 | 1 × 10^{15} | 154 |

81 | 281 | 27.2 | 0.43 | 5 × 10^{15} | 69 |

76 | 288 | 26.6 | 0.52 | 1 × 10^{16} | 54 |

71 | 295 | 25.9 | 0.60 | 9 × 10^{16} | 19 |

### A. Plasma sheath without electron emission

^{58}

^{,}

*v*and using ion flux conservation $n0v0=ni(r)vi(r)$, leads to $ni(r)=n0\varphi 0\varphi 0\u2212\varphi (r)$. Assuming $\varphi (r)\u226a\varphi 0$ and expanding in a Taylor series, $ni\u2248n0(1\u2212\varphi 2\varphi 0)$ and $ne\u2248n0(1\u2212e\varphi Te)$. Poisson's equation then becomes

_{i}*n*and

_{i}*n*are the plasma ion and electron densities, respectively. The positive curvature of the sheath potential guarantees that the ions are not reflected, which gives the Bohm acoustic velocity criterion.

_{e}^{60}The corresponding

*potential drop*in the presheath = - $\varphi 0=\u2212Te2e$ relative to the plasma.

^{58}Since the plasma core is quasineutral, this will correspond to the value of $\varphi 0$, which we will use as a boundary condition. The Bohm criterion is that the speed of the ions as they enter the sheath must exceed the sound speed of the ions $cs=Temi$. The electron density at the edge of the sheath is $ne=n0\u2009exp(\u2212e\varphi 0Te)=n0\u2009exp[(eTe)(\u2212Te2e)]\u22480.6n0$.

*n*

_{0}= 9 × 10

^{16}m

^{−3}to ensure the potential increased monotonically across the sheath. The FEM method is now used to solve the corresponding weak form of (29) using Eqs. (7) and (8), with the source term $Q(\varphi )$ defined in Eq. (30). The assumed shape, dimensions, and mesh details are shown in Fig. 6. The geometry of the leading edge in Fig. 6(a) has a 1 mm radius of curvature and a 5° wedge angle. We modeled the sheath region ahead of the leading edge as the 200

*μ*m annular region shown in Fig. 6(b), where the floating potential, Eq. (25), is applied at the inner boundary [i.e., $\varphi (RLE)=\varphi s$] representing the surface of the vehicle. The total mesh had 1.01 × 10

^{6}elements.

The normalized number density, $n/n0$, normal to the stagnation point is plotted in Fig. 7(a). As expected, the normalized electron density is $n/n0$ = 0.6 when $\varphi (Rpre)=\varphi 0$ is the boundary condition on the presheath surface. When $\varphi (Rpre)$ = 0 for *n*_{0} = 9 × 10^{16} m^{−3}, the edge of the sheath is closer to the surface and $n/n0$ = 0.6 at $\u2248$ 0.09 mm from the stagnation point. Moving further away from the surface, the ion and electron densities reach equality, satisfying the quasineutral plasma condition. The potential (solid lines) and the electric field (dashed lines) across the sheath are shown in Fig. 7(b). The potential increases monotonically from the surface potential [Eq. (25)] to the pre-sheath potential, $\varphi (Rpre)$. The magnitude of the surface potential increases as the plasma density increases. When *n*_{0} = 9 × 10^{16} m^{−3}, $\varphi (Rpre)$ = 0 and $\varphi $ reaches $\varphi 0$ much closer to the surface. The maximum value of the electric field occurs at the surface, and the magnitude increases with plasma density. The magnitude of the electric field at the presheath boundary decreases with higher plasma densities. The normalized potential, $\varphi /\varphi s$, across the sheath for the plasma conditions in Table I is shown in Fig. 8. In comparison of the images, the potential gradient across the sheath increases as the plasma density increases, indicating that the sheath thickness is decreasing.

### B. Plasma sheath with electron emission

A space potential well is formed near the LE surface when electrons are emitted and the surface at which $E=0$ moves into the plasma region from the LE surface and is called the “virtual cathode” (VC). We will extend the 1D model proposed by Takamura *et al.*^{50} and summarize the key equations that lead to a weak variational solution, while detailed derivations are given in Ref. 50. The schematic shown in Fig. 9 divides the region between the wall and the plasma into two regimes, *α* and *β*, where the region *α* extends from the plasma to the virtual cathode, while the thin layer *β* describes the region between the VC and the solid surface.^{50} In this work, we model the region *α* using the FEM numerical method, while the region *β* is modeled using Eqs. (41) and (42) of Ref. 50.

*A*= 120 A cm

_{RD}^{−2}K

^{−2}is the R-D constant, which is the theoretical limit. Here, we use a material-specific value for LaB

^{6}equal to 30 A cm

^{−2}K

^{−2}.

^{64}The work function of the LE surface is $\varphi work$ (V), and its temperature is

*T*(K). Once space charge accumulates ahead of the LE surface, the VC potential develops and electrons must overcome an energy barrier of ( $\varphi VC\u2212\varphi s$), and the current that goes through the VC surface into the plasma region is now

_{w}*λ*, field will be determined by the numerical FEM solution. Physically, it represents a local “enhancement” (“damping”) factor that describes the increase (decrease) in electron emission from the LE surface above

*J*. We can find an expression for the unknown virtual potential in terms of the Lagrange multiplier field,

_{CL}*d*,

*et al.*

^{65}used the normalized sheath thickness, $k=d/\lambda D$ = 2.2, for a floating electrode. Hagino

*et al.*

^{66}derived a modified sheath thickness that was a function of the sheath voltage (i.e., $\varphi s$). In practice, the sheath thickness may be several multiples of the Debye length and vary with the plasma conditions. The Debye length is used here for consistency as it is not expected to greatly impact the trends in our results nor detract from the primary thrust of the paper, which is to demonstrate the efficacy of the variational scheme.

^{50}in terms of the ion Mach number, $M\u2261v0/Cs$, and the plasma ion density,

*n*

_{0}. The ion Mach number,

*M*, differs from the freestream Mach number, $M\u221e$, and satisfies the Bohm criterion in Sec. IV A,

^{50}we finally determine the split of the electron density between the plasma, $nep$, and the surface emission, $nes$, and give the following expressions:

*λ*as

*n*

_{0}= 9 × 10

^{16}m

^{−3}. In Eq. (21), the electric field at $RVC$ is

*E*= 0, which is the required condition for space-charge limitation. Note that, according to Fig. 9, the VC is located just ahead of the LE so that $\varphi VC\u2260\varphi s$. The separation distance in region

_{s}*β*is assumed to be much less than the thickness of the sheath (i.e., the annular region of 200

*μ*m of the computational domain).

^{50}Figure 10 shows the problem setup for the FEM variational solution with electron emission. The assumed leading-edge geometry is identical to Fig. 6(a). The meshed model geometry for the case with electron emission had 0.96 × 10

^{6}total elements. Electrons are emitted from the curved section of the leading edge, highlighted in Fig. 10(a), and the electric field is constrained due to the space-charge condition. Figure 10(b) shows the increase in the mesh density near the emitting surface. This is necessary to capture the large electric field gradient due to the space-charge constraint and the transition from a curved to a flat surface.

The normalized densities of the ions (solid lines), plasma electrons (dashed lines), and emitted electrons (dashed-dotted lines) are shown in Fig. 11(a), where the abscissa is the distance relative to the VC. The plasma electron density is depleted at the VC and increases to the plasma density, *n*_{0}, near the sheath edge. The emitted electron density is maximized at the VC but is a small fraction of the plasma density. In particular, the emitted electron density increases with the plasma density, indicating that a more dense plasma can alleviate the space-charge limitation, resulting in a larger emitted electron density at the VC. The potential distribution and electric field across the sheath are shown in Fig. 11(b). The electric field at the virtual cathode is constrained because of the space-charge condition. We note here that the electric field decreases sharply very near the VC surface and vanishes exactly at the VC surface. This behavior is clearly seen in the figure inset, where details of the electric field distribution near the VC are shown with a resolution of 0.5 *μ*m. The magnitude of the potential in the VC and the potential gradient (i.e., the electric field) across the sheath increases with the plasma density with an exception because when *n*_{0} = 5 × 10^{15} m^{−3}, the plasma temperature, and thus the assumed surface temperature are identical to the *n*_{0} = 1 × 10^{15} m^{−3} case. Examining Eqs. (34) and (35), the product of $\lambda JCL$ will determine the value of $Jes$ and $\varphi VC$. The space charge should be more pronounced for a lower plasma density [Eq. (36)], and the emitted current (that is, $\lambda JCL$) should therefore be smaller, resulting in a larger magnitude of $\varphi VC$.

The normalized virtual cathode potential is shown in Fig. 12(a). Electrons are only emitted from the curved section of the leading edge, which corresponds to the central portion of the graph. The normalized virtual cathode potential increases with plasma density, with the exception of the case *n*_{0} = 5 × 10^{15} m^{−3}, as discussed in the previous paragraph. The electric field across the curved section of the leading edge in Fig. 12(b) is constrained to zero, satisfying the space-charge condition. Figure 12(c) shows the emitted current, $Jes$, across the leading edge. Large gradients at the transition from the curved section to the flat surface are caused by the electron's ability to emit toward a more favorable electric field in this region, similar to the 2D validation case [see Fig. 3(a) and discussion in Sec. III B]. The emitted current is zero on the flat surface (i.e., no emission). The Debye length, *λ _{D}*, decreases as the plasma density increases, and therefore, the modified CL current increases ( $\u221d\varphi s3/2/\lambda D2$). In the limit of very small plasma density (i.e., vacuum), a very large value of the Debye length occurs, resulting in negligible emitted current. Thus, the thermionic emission becomes stronger at lower altitudes and increases the plasma density.

The normalized potential distribution, $\varphi /\varphi s$, in Figs. 13(a) and 13(b) shows the spatial structure of the sheath around the curved and flat sections of the LE for *n*_{0} = 1 × 10^{1} and 9 × 10^{16} m^{−3}, respectively. As expected, the sheath thickness decreases as the plasma density increases. The transition from the curved region to the flat surface is shown in detail in Figs. 13(c) and 13(d). The magnitude and direction of the electrostatic force on the emitted electrons ( $Fe=\u2212eE$) in this region highlight the effects of the strong fringe field on the acceleration of electrons away from the region of the curved surface of the surface. Such strong fringe field effects are a consequence of the imposed boundary conditions on the surface potential, forcing high gradients close to the transition regions from an emitting to a nonemitting surface. This observation is consistent with the earlier work of Watrous and coworkers,^{56} as can be seen in Fig. 3(a), and may be utilized as a method to enhance the *average* value of the emitted current.

By introducing extra electrical bias, we can regulate the surface potential, thereby altering the sheath and emitted electron current. The potential distribution across the sheath is shown in Fig. 14(a) for $\varphi app$ = 0 (floating surface), −1, −2, −3, and −10 V when *n*_{0} = 1 × 10^{15} (solid lines) and 9 × 10^{16} m^{−3} (dashed lines), the two extreme plasma densities analyzed so far. The distance is measured from the virtual cathode surface in the normal direction to the tip of the LE. Again, the potential gradient increases with the plasma density. The electric field in Fig. 14(b) is constrained at the virtual cathode due to the space-charge condition, and the magnitude adjacent to the virtual cathode increases with the plasma density and the applied potential. The normalized virtual cathode potential, $\varphi VC/\varphi s$, is shown as a function of the arc length distance in Fig. 15(a). The normalized virtual cathode potential decreases as $\varphi app$ increases for both plasma densities. In Fig. 15(b), $Jes$ also increases, showing that the application of additional electrical bias can be used to modify the virtual cathode potential and enhance electron emission, an important observation for the design of thermionic systems for thermal management and power generation.

The emitted current, $Jes$, is also a function of the wall temperature through the dependence on the Richardson–Dushman current in Eq. (34). Up to this point, we have considered the surface temperature to be a constant value of $Te/3$. Various heat fluxes will be important in the temperature ranges of interest, most notably radiative cooling due to its dependence on $Tw4$, the combination of which will determine *T _{w}*. Next, we vary

*T*at a constant plasma density of

_{w}*n*

_{0}= 9 × 10

^{16}m

^{−3}and show the variation in the emitted current, $Jes$. Figure 16(a) shows $Jes$ vs

*T*for $\varphi app$ = 0 (floating), −3, and −5 V, along with Eq. (33). The numerical values of $Jes$ are less than the Richardson–Dushman current for the temperature range of 1700–2500 K. It is also a weak function of

_{w}*T*, but a much larger dependence on $\varphi app$ is still observed. Although our results are not for vacuum conditions, the dependence of the emitted current on surface temperature is consistent with recent literature on the Miram curves, where the transition from thermionic emission to space-charge limit in the vacuum case has been significantly advanced.

_{w}^{51,53,67}This builds on the Longo scaling for the shape of the transition as a function of current.

^{67,68}

We showed previously that we can modify $Jes$ by applying an additional potential, $\varphi app$, to the surface above the floating potential [Eq. (26)]. We extend this part of the study in Fig. 16(b), which shows the relationship between $Jes$ and the total surface potential, $\varphi s=\varphi f+\varphi app$. Initially, $Jes$ increases quickly with $\varphi s$ but levels off, suggesting some limit to the enhancement in current achieved by applying an electric potential to the surface.

## V. SUMMARY AND CONCLUSIONS

We have introduced an innovative variational formulation for the constrained Poisson equation, optimally tailored for numerical FEM solutions of the potential distribution within a plasma sheath under conditions of surface electron emission. The efficacy of this approach stems from a fundamental aspect of weak formulations of partial differential equations: the ability to augment the functional description with additional constraints through the use of Lagrange multipliers. This methodology is advantageous for various applications where surface conditions are imposed, such as surface charging, biasing, secondary electron emission, and surface ion injection. The computational method supports a wide range of applications involving highly complex three-dimensional shapes, using a segmented application of Dirichlet and Neumann boundary conditions. In our particular case, which deals with enforcing space-charge limitations on a portion of the leading edge of a hypersonic aircraft, the Lagrange multiplier is defined as the ratio of emitted current to the space-charge limited current. This formulation guarantees the enforcement of the space-charge constraint at the electron-emitting section (curved area) of the surface, thus allowing the use of both Dirichlet and Neumann boundary conditions.

At first, we matched our numerical solutions with the analytical solution of the classic 1-D CL problem, attaining an accuracy of 0.9% for the potential and 3% for the electric field. Next, we extended the validation to a 2D CL-type scenario, where electrons are emitted from a narrow strip on the cathode in vacuum. We observed an increase in the density of emitted current beyond the CL limit and demonstrated that the total emitted current decreases with the ratio *w*/*d*, which is the width of the strip to the distance from the cathode to the anode. The numerical results of our model were within a few percent of those obtained by Watrous and co-workers,^{56,57} who used a finite-volume solution combined with iterative methods to satisfy the space-charge condition. Moreover, we revealed an additional dependence of the ratio of the total emitted current to the C-L value on the width of the emitting strip.

The variational approach to Poisson's equation was utilized to address a practical issue, specifically to ascertain the spatial configuration of the sheath potential around the critical leading edge of a hypersonic aircraft. The example considered here involves a leading edge with a 1 mm radius of curvature and a 5° wedge angle. It was assumed that the curved section is composed of a material capable of emitting a thermionic electron current to assist in cooling, with the emitted current being absorbed into the plasma region. The primary challenge is that such emission is constrained by the buildup of space charges near the emitting surface, which can hinder this cooling mechanism.

We demonstrated that the emitted surface current is typically a small fraction of the CL current because of the space charge, and hence, surface cooling via thermionic emission would have a small effect. However, there are several possible ways to enhance the emitted current from the surface and increase the cooling contribution of thermionic emission. First, the emitted current was found to increase as the plasma density increases ahead of the LE, which is the case during lower-altitude flights. Second, a larger negative surface potential was found to boost the thermionic emission. Third, the total emitted current can be further enhanced through geometry, as emitted electrons are likely to escape the space-charge repulsion near the termination regions of the emitting surface, similar to observations under vacuum conditions.^{69} Fringe electric fields occur because the electric field lines do not terminate abruptly at the edges of conductive surfaces. Instead, they spread out or “fringe” around the edges, leading to a less uniform field in these regions, and hence allowing more electrons to escape the space charge ahead of the emitting strip.

We have identified a relationship between surface temperature and the magnitude of the emitted current specifically for LaB_{6}. In the absence of an applied surface potential, we showed that the total emitted current remains minimal and exhibits no dependence on the surface temperature. Consequently, elevating the surface temperature does not facilitate an increase in the thermionically emitted current. Conversely, when a negative potential is applied to the emitting surface, the emitted current adheres to the Richardson–Dushman equation up to a critical temperature, beyond which it becomes invariant with respect to the surface temperature. For example, at an applied surface potential of −5 V, the critical temperature is approximately 1700 K.

## ACKNOWLEDGMENTS

The support of the U.S. Department of Energy under Project No. DE-SC0018410:0008 and the Defense Advanced Research Projects Agency, DARPA, MACH program Project No. DOD-ARO W911NF2010016 at UCLA is acknowledged.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Nasr M. Ghoniem:** Conceptualization (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). **David B. Brown:** Investigation (equal); Software (equal); Writing – original draft (equal). **Yue Huang:** Software (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.

### APPENDIX: NUMERICAL METHOD FOR CONSTRAINED WEAK FORMS OF PDEs

^{70,71}Consider a general PDE of the form: $L\varphi =f\u2002in\u2002\Omega $ with boundary conditions: $\varphi =g\u2002on\u2002\Gamma D$ and $\u2202\varphi \u2202n=h on \Gamma N$ and an additional constraint over the domain: $\u222b\Omega \varphi \u2009d\Omega =C$, where $L$ is a differential operator, $\Omega $ is the domain, $\Gamma D$ is the Dirichlet boundary, $\Gamma N$ is the Neumann boundary, and $C$ is a given constant. Multiply the PDE by a test function $\varphi \u0303$ and integrate over the domain $\Omega $: $\u222b\Omega \varphi \u0303L\varphi \u2009d\Omega =\u222b\Omega \varphi \u0303f\u2009d\Omega $. Now, applying integration by parts to transfer derivatives from $\varphi $ to $\varphi \u0303$ (for second-order PDEs),

*C*is a constant. The augmented weak form becomes

## REFERENCES

*Entry Vehicle Heating and Thermal Protection Systems: Space Shuttle, Solar Starprobe, Jupiter Galileo Probe*

*Structures Technology for Future Aerospace Systems*

*The Emission of Electricity from Hot Bodies*

*Fundamentals of Electric Propulsion: Ion and Hall Thrusters*

*Minimum Ionic Kinetic Energy for a Stable Sheath*

*Handbook of Thermionic Properties: Electronic Work Functions and Richardson Constants of Elements and Compounds*

*The Finite Element Method: Linear Static and Dynamic Finite Element Analysis*