1 Introduction
As a quasi-neutral plasma expands into a vacuum, the electrons tend to diffuse much faster the ions, creating a charge separation that accelerates the ions beyond their characteristic sound speed. Dating back to the 1960s, a wealth of literature exists on ambipolar plasma expansion into an infinite vacuum (see early reviews Samir, Wright & Stone Reference Samir, Wright and Stone1983; Sack & Schamel Reference Sack and Schamel1987). The present study addresses a new case where the vacuum is bounded by an electrically insulated surface. The surface could be electrically conducting or insulating, so long as it is an isolated plane. As the surface is charged by the cloud of oncoming electrons, it develops an electric field which repels the electrons and accelerates the ion expansion. The ions, in turn, carry the electron cloud further towards the surface. We derive an analytic theory to explore this feedback process.
We begin with a brief review of the existing literature, emphasizing the prevailing focus on infinite plasma expansion into an infinite vacuum domain. The pioneering work by Gurevich, Pariiskaya & Pitaevskii (Reference Gurevich, Pariiskaya and Pitaevskii1966) introduced the self-similar approach, which assumes quasi-neutrality (neglecting Poisson’s equation) to obtain an analytic model of collisionless isothermal unbounded plasma expansion. Early numerical investigations by Widner, Alexeff & Jones (Reference Widner, Alexeff and Jones1971) and Crow, Auer & Allen (Reference Crow, Auer and Allen1975) established the self-similar model as a reasonable approximation of the bulk plasma behaviour, after a charge-equilibration time and up to a moving boundary where quasi-neutrality breaks down. Beyond the physical boundary, called the ion expansion front, the self-similar solution is not valid. A salient feature of the non-quasi-neutral numerical models is the inclusion of a separate electron cloud region, spanning beyond the expansion front out to many times the Debye length of the quasi-neutral plasma. Semi-analytic models with a more accurate portrayal of the ion front boundary layer (Mora Reference Mora2003; Medvedev Reference Medvedev2011) have shown that the self-similar domain of validity can be extended earlier in time and further in space. Many variations of this problem, beyond the scope of the present discussion, include non-isothermal expansion (Mora & Pellat Reference Mora and Pellat1979; Grismayer et al. Reference Grismayer, Mora, Adam and Héron2008), non-Maxwellian initial distributions (Gurevich & Meshcherkin Reference Gurevich and Meshcherkin1981a; True Reference True1981; Akbari-Moghanjoughi Reference Akbari-Moghanjoughi2015; Bennaceur-Doumaz et al. Reference Bennaceur-Doumaz, Bara, Benkhelifa and Djebli2015), magnetic fields (García-Rubio, Ruocco & Sanz Reference García-Rubio, Ruocco and Sanz2016) and multiple ion species (Elkamash & Kourakis Reference Elkamash and Kourakis2016). Some attention has been given to the case of a finite plasma (Medvedev Reference Medvedev2005; Murakami & Basko Reference Murakami and Basko2006), but still with infinite boundary conditions in the vacuum.
The present study examines the interaction of an expanding plasma with a finite insulated surface. As we demonstrate, our system is quite distinct from the fixed-potential (anode–cathode like) set-up addressed by Gurevich & Pitaevsky (Reference Gurevich and Pitaevsky1975). Our study is motivated by the simulations carried out by Zimmerman et al. (Reference Zimmerman, Farrell, Stubbs, Halekas and Jackson2011), who observed surface feedback on solar-wind plasma expansion into a lunar polar crater. For a steady-state plasma wake structure generated by the uniform horizontal flow of plasma over a large cavity, Zimmerman et al. (Reference Zimmerman, Farrell, Stubbs, Halekas and Jackson2011) find an order of magnitude increase in the steady-state electrostatic potential, compared to the case of free expansion. The identified plasma–surface interaction mechanism – which is the subject of the present study – involves an electron cloud and associated electric field extending from the plasma wake front all the way down to the crater surface. Analytic efforts to describe the solar-wind plasma wake in a lunar crater, by Farrell et al. (Reference Farrell, Stubbs, Halekas, Killen, Delory, Collier and Vondrak2010) and Rhodes & Farrell (Reference Rhodes and Farrell2019), serve as the basis for the present model.
The applicability of the one-dimensional (1-D) approximation depends upon the specific problem of interest. For the local plasma wakes along the lunar surface (Farrell et al. Reference Farrell, Stubbs, Halekas, Killen, Delory, Collier and Vondrak2010; Rhodes & Farrell Reference Rhodes and Farrell2019), neglect of other spatial dimensions is based on two assumptions; (i) a large width to depth ratio leading to a relatively insignificant horizontal field from crater wall charging, and (ii) relatively small surface formations on the crater floor relative to the system size. In essence, we require that the combined ambipolar and surface-generated electric field dominate over transverse field effects. This condition is justified by surface wake simulations (Zimmerman et al. Reference Zimmerman, Farrell, Stubbs, Halekas and Jackson2011) as well as in data describing the analogous macroscopic wake extending about the entire Moon (Farrell et al. Reference Farrell, Kaiser, Steinberg and Bale1998; Xu et al. Reference Xu, Poppe, Halekas, Mitchell, McFadden and Harada2019). Motivated by these lunar phenomena, the present study aims to address the general problem of bounded plasma expansion, as a starting point for similar systems.
The recent analytical model by Rhodes & Farrell (Reference Rhodes and Farrell2019) matches the self-similar solution in the expansion region with a Poisson solution of a bounded electron cloud, to examine charging of an insulated surface in a steady-state plasma wake. We now extend the work of Rhodes & Farrell (Reference Rhodes and Farrell2019) to account for the dynamic feedback effect of surface charge accumulation on the plasma expansion process. The new boundary condition matching is accomplished by an additional solution of Poisson’s equation in a thin layer called the ion front boundary layer, which connects the quasi-neutral layer and the electron cloud. Our four-layer model structure mimics that of Medvedev (Reference Medvedev2011), with the advantage of being analytically derived and self-contained rather than interpolated from simulation data.
In the no-surface limit, we find the same qualitative behaviour as found in Medvedev (Reference Medvedev2011), with the ion front boundary layer expanding linearly in time. Unlike the study by Medvedev (Reference Medvedev2011) and a similar semi-analytic study by Mora (Reference Mora2003), however, the present study focuses on the early expansion phase rather than the long-time asymptotic behaviour. Furthermore, we note that our thin-layer approximation is not intended to improve on the aforementioned ion front models (Mora Reference Mora2003; Medvedev Reference Medvedev2011), but rather to demonstrate the key role of the ion front in plasma–surface coupling.
For the case of bounded expansion, our results show that an insulated surface may be quickly charged by the oncoming electrons. We demonstrate how the presence of the charge-accumulating surface can cause a widening of the ion front boundary layer, acceleration of the ion front, increase of the surface potential magnitude and decrease of the near-surface electron density.
2 Four-layer expansion model
2.1 Model overview
Consider an initially quasi-neutral quiescent hydrogen plasma, with comparable ion and electron temperatures, filling the half-plane defined by 
$-\infty <\unicode[STIX]{x1D701}<0$ at 
$t=0$. The classic 1-D plasma expansion problem supposes an infinite vacuum filling the complementary space 
$0<\unicode[STIX]{x1D701}<\infty$. The entire space is unmagnetized. We now consider the effect of an insulated surface at fixed position 
$\unicode[STIX]{x1D701}_{s}$, assumed much farther than an initial plasma Debye length (
$\unicode[STIX]{x1D701}_{s}\gg \unicode[STIX]{x1D706}_{D}$). We emphasize that our initial condition limits both ions and electrons to the half-plane, in contrast with many existing models which initialize the electrons in an infinite Boltzmann equilibrium. This initial condition is necessary to account for the surface charging and resulting plasma–surface interaction.
Initial plasma expansion process, depicted by the ion density (dashed red) and electron density (solid blue) spatial distributions at two different times. In the initial phase (a) the electrons expand outward while the ions stay relatively stationary. The rightward-moving ion expansion front follows at a later time (b), accompanied by a leftward-moving rarefaction front. The densities are described by (A 4) for the self-similar domain (where 
$n_{e}=n_{i}$), equations (B 1) and (B 8) for the ion front layer, and (C 8) for the electron cloud domain (where 
$n_{i}=0$).
The classic unbounded expansion process is qualitatively described as follows: owing to a relatively high thermal speed, electrons rapidly diffuse into the void while the ions remain approximately fixed. This electron expansion phase is depicted in the first tier of figure 1. Charge separation creates an ambipolar electric field which accelerates the ions into the negative charge region, referred to as the electron cloud (to the right in figure 1). The expansion is accompanied by a rarefaction front, a boundary of plasma depletion, which travels in the opposite direction at the ion acoustic speed (to the left in figure 1). Up to the expansion front, the plasma bulk remains quasi-neutral. At the expansion front there is an ion front boundary layer – a narrow region of positive net charge – beyond which lies the electron cloud. The second tier of figure 1 illustrates the evolution of the expanding plasma density, with the left-moving rarefaction front and right-moving expansion front and ion front boundary layer. Note that in reality the rarefaction front also contains a narrow positive charge layer (Mora Reference Mora2003), not depicted in figure 1, which is neglected in the present model since it has no bearing on the plasma–surface interaction. The spatial density distributions presented in the figure are based on (A 4) for the quasi-neutral plasma domain (where 
$n_{e}=n_{i}$), equations (B 1) and (B 8) for the ion front layer, and (C 8) for the electron cloud domain (where 
$n_{i}=0$).
The plasma expansion process evolves on the time scale of the ion plasma frequency 
$\unicode[STIX]{x1D714}_{i}$. The widely used self-similar expansion model (see appendix A), incorporated as part of the present model, takes 
$\unicode[STIX]{x1D714}_{i}t=1$ as the time step for a quasi-neutral expansion to initialize. For the length scales of interest, 
$\unicode[STIX]{x1D701}_{s}\sim 10{-}100\unicode[STIX]{x1D706}_{D}$, we find that 
$\unicode[STIX]{x1D714}_{i}t\approx 1$ is a reasonable time frame for the initial electron cloud to fill the gap between the plasma and the insulated surface. Thus we take 
$t_{0}=\unicode[STIX]{x1D714}_{i}^{-1}$ to be the initial time for both surface charging and ion front expansion. One oddity of this initial condition is the slightly negative location of the ion front at 
$\unicode[STIX]{x1D714}_{i}t=1$, 
${\sim}-0.6\unicode[STIX]{x1D706}_{D}$, owing to the singular asymptotic behaviour of the self-similar formulation approaching 
$t=0$ (Mora Reference Mora2003). Keeping this caveat in mind, we focus the present work on the qualitative dynamics of plasma–surface interaction, which is largely insensitive to the precise ion front location for a reasonably distant surface.
Both the cold ion plasma and the electron cloud are assumed to expand isothermally, and maintain a Maxwellian velocity distribution via the Boltzmann relation. When the isothermal assumption – for both electrons and ions – is relaxed, numerical computations of unbounded plasma expansion show that the thermal temperature decay (the energy source for ion acceleration) causes a substantial departure from the self-similar solution at the scale where 
$\unicode[STIX]{x1D701}/t\sim (Zm_{i}/m_{e})^{1/2}$ (Denavit Reference Denavit1979; Mora & Pellat Reference Mora and Pellat1979). Here 
$Z$, 
$m_{i}$ and 
$m_{e}$ respectively represent the ion charge number, ion mass and electron mass. Given a hydrogen plasma (
$(Zm_{i}/m_{e})^{1/2}\approx 43$), our range of interest for the initial plasma expansion (
$\unicode[STIX]{x1D701}/t<4$) remains well within the validity of a first-order isothermal approximation. Furthermore, while technically the presence of a charge-accumulating surface causes a truncated velocity distribution that brings into question the application of Boltzmann’s relation, this effect can be shown to be quite small over the length scales of interest. Upon arrival of the electron cloud at the surface, prior to charge accumulation, we find typical potential values of 
$\unicode[STIX]{x1D719}\sim -5kT/e$. Under these conditions only electrons faster the 3 times the thermal speed arrive at the surface, 
$v>3v_{te}$, or approximately 1 % of a Maxwellian velocity distribution. As the surface potential continues to build, the distribution truncation shrinks well below 1 %, and thus a Maxwellian distribution is reasonably justified.
Four solution domains of the plasma expansion process; (i) undisturbed plasma, (ii) quasi-neutral expansion, (iii) ion front boundary layer and (iv) electron cloud, defined by the distribution of the electron density (solid blue) and ion density (dashed red) at a fixed time. The densities are described by (A 4) for the self-similar domain (where 
$n_{e}=n_{i}$), equations (B 1) and (B 8) for the ion front layer and (C 8) for the electron cloud domain (where 
$n_{i}=0$).
The overall plasma expansion structure can be broken into qualitatively distinct spatial regions. We adopt the four-layer framework of Medvedev (Reference Medvedev2011), which explicitly incorporates a positively charged boundary layer at the plasma expansion front. However, unlike the layer model by Medvedev (Reference Medvedev2011), which is interpolated based on simulation data, the present model applies an analytic approximation to describe the initial expansion phase in a self-contained manner. Our analytic formulation facilitates the addition of an insulated surface boundary condition, which is the focus of the present study.
In terms of the spatial variable 
$\unicode[STIX]{x1D701}$, the four layers are the (i) rarefaction front 
$\unicode[STIX]{x1D701}_{r}$, (ii) boundary of quasi-neutrality 
$\unicode[STIX]{x1D701}_{b}$, (iii) ion front 
$\unicode[STIX]{x1D701}_{f}$ and (iv) insulated surface 
$\unicode[STIX]{x1D701}_{s}$. The corresponding four layers, depicted in figure 2 are:
(i) undisturbed plasma:
$-\infty <\unicode[STIX]{x1D701}<\unicode[STIX]{x1D701}_{r}$;(ii) quasi-neutral expanding plasma:
$\unicode[STIX]{x1D701}_{r}<\unicode[STIX]{x1D701}<\unicode[STIX]{x1D701}_{b}$;(iii) ion front boundary layer:
$\unicode[STIX]{x1D701}_{b}<\unicode[STIX]{x1D701}<\unicode[STIX]{x1D701}_{f}$;(iv) electron cloud:
$\unicode[STIX]{x1D701}_{f}<\unicode[STIX]{x1D701}<\unicode[STIX]{x1D701}_{s}$.
We assume continuity of 
$\unicode[STIX]{x1D719}$ and 
$n_{e}$ at all boundaries. The ion front boundary layer is introduced to resolve a continuous electric field at the interface of the quasi-neutral plasma and the electron cloud. The layer width, 
$\unicode[STIX]{x1D6E5}\equiv \unicode[STIX]{x1D701}_{f}-\unicode[STIX]{x1D701}_{b}$, is a free parameter to be uniquely determined in the electric field matching process. This additional boundary condition is necessary to account for the effect of the external field generated by the charge-accumulating surface.
Prior to the application of boundary conditions, independent solutions are obtained in each of the four regions, as described below.
2.2 Independent solutions
We describe the model in general terms and present its independent domain solutions, details of which are found in appendices A–C. Note that domain (i) contains undisturbed plasma with a trivial solution.
A description of domain (ii) – the quasi-neutral expansion region – is described by the well-known self-similar solution (Gurevich et al. Reference Gurevich, Pariiskaya and Pitaevskii1966), detailed in appendix A. The self-similar solution neglects charge separation (dropping Poisson’s equation), assuming a Maxwellian electron fluid that equilibrates with the ion fluid much faster than the ion dynamic time scale 
$\unicode[STIX]{x1D714}_{i}^{-1}$ (defined below). The quasi-neutral plasma layer is considered impervious to the external electric field. This last assumption, we believe appropriate for the relatively short expansion time scales considered in the present study. Field penetration into an expanding plasma is an open problem, beyond the scope of the present paper.
Domain (iii) – the ion front boundary layer – is approximated by a constant ion density, equal to the density at the boundary of quasi-neutrality 
$\unicode[STIX]{x1D701}_{b}$. With the electron density spatially decreasing, this set-up emulates the thin positive charge layer predicted by early numerical computations (Widner et al. Reference Widner, Alexeff and Jones1971; Crow et al. Reference Crow, Auer and Allen1975). Assuming a Maxwellian electron distribution and small potential variation, the electron density and potential are determined by Poisson’s equation (see appendix B). A key feature of this new formulation is that the layer width, rather than being pre-set, is determined by the requirement of electric field continuity. Although the choice of a constant density solution is not unique, it has the benefit of providing a convenient analytic description of the early expansion phase. The resulting thin-layer solution generates a reasonable first-order approximation, providing insight into the physical mechanism by which a charge-accumulating surface can influence plasma expansion.
The solution for domain (iv) – the electron cloud – is determined by solving Poisson’s equation for a Maxwellian electron distribution and zero ion density. The electron density remains continuous at 
$\unicode[STIX]{x1D701}_{f}$, while the ion density jumps to zero. The resulting solution, derived in appendix C, depends on the surface charge density and resulting electric field at the insulated surface 
$\unicode[STIX]{x1D701}_{s}$. In the limit of an infinitely far or uncharged surface, the solution reduces to the known formula for a steady-state Maxwellian electron cloud (Rhodes & Farrell Reference Rhodes and Farrell2019).
The equations below are expressed in terms of the following dimensionless units
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D701}\sim \unicode[STIX]{x1D706}_{D}=\sqrt{{\displaystyle \frac{\unicode[STIX]{x1D716}_{0}kT_{e}}{e^{2}n_{0}}}},\quad t\sim \unicode[STIX]{x1D714}_{i}^{-1}=\sqrt{{\displaystyle \frac{\unicode[STIX]{x1D716}_{0}m_{i}}{e^{2}n_{0}}}},\\ n\sim n_{0},\quad \unicode[STIX]{x1D719}\sim {\displaystyle \frac{kT_{e}}{e}}.\end{array}\right\}\end{eqnarray}$$ Here, 
$n_{0}$ denotes the electron density of the undisturbed plasma. The ion temperature is neglected, and the electron temperature 
$T_{e}$ is assumed constant throughout. In these dimensionless units, the ion acoustic and electron thermal speed are given by 
$c_{i}=1$ and 
$v_{Te}=\sqrt{m_{i}/m_{e}}$. The resulting electric field units are 
$E_{0}=(n_{0}kT_{e}/\unicode[STIX]{x1D716}_{0})^{1/2}$.
Assuming potential and electron density continuity across all boundaries, we find the following solutions for the potential 
$\unicode[STIX]{x1D719}$ in the four domains:
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D719}_{i}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D719}_{ii}=-(\unicode[STIX]{x1D701}/t+1), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D719}_{iii}=\unicode[STIX]{x1D719}_{b}(t)-\sinh [(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D701}_{b}(t))/t], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D719}_{iv}=\unicode[STIX]{x1D719}_{f}(t)-2\log (\unicode[STIX]{x1D6FC}^{-1}\sinh [g_{\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x1D701})]). & \displaystyle\end{eqnarray}$$ Here, the subscript expressions refer to the spatial boundary; 
$\unicode[STIX]{x1D719}_{b}(t)=\unicode[STIX]{x1D719}(t,\unicode[STIX]{x1D701}_{b})$ and 
$\unicode[STIX]{x1D719}_{f}(t)=\unicode[STIX]{x1D719}(t,\unicode[STIX]{x1D701}_{f})$. Functions 
$\unicode[STIX]{x1D701}_{b}(t)$ and 
$g_{\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x1D701})$ are respectively defined in (A 3) and (C 9). Here, the free parameter 
$\unicode[STIX]{x1D6FC}$ characterizes the electron fluid pressure balance
$$\begin{eqnarray}\unicode[STIX]{x0394}P_{e}(\unicode[STIX]{x1D701})={\textstyle \frac{1}{2}}E^{2}(\unicode[STIX]{x1D701})-n_{e}(\unicode[STIX]{x1D701})\equiv n_{f}\unicode[STIX]{x1D6FC}^{2},\end{eqnarray}$$ which is constant throughout the electron cloud region (see appendix C, equation (C 4)). The magnitude of 
$\unicode[STIX]{x1D6FC}$ can be thought of as a measure of the surface electric field contribution. The limit 
$\unicode[STIX]{x1D6FC}\rightarrow 0$ corresponds to an equilibrium pressure balance with an infinitely far surface 
$\unicode[STIX]{x1D701}_{s}\rightarrow \infty$. This special case (
$\unicode[STIX]{x0394}P_{e}=0$) is equivalent to the finite surface floating potential derived by Rhodes & Farrell (Reference Rhodes and Farrell2019), neglecting the kinetic surface charging effect. In the general case, 
$\unicode[STIX]{x1D6FC}\neq 0$, the parameter 
$\unicode[STIX]{x1D6FC}$ defines a family of solutions based on the surface field pressure 
$E_{s}^{2}/2=n_{s}+n_{f}\unicode[STIX]{x1D6FC}^{2}$. The solution in the next section (2.13) relates 
$\unicode[STIX]{x1D6FC}$ to the ion front boundary layer width 
$\unicode[STIX]{x1D6E5}$, through the requirement of electric field continuity.
The surface-generated electric field is assumed to penetrate only into the ion front boundary layer, and not into the quasi-neutral plasma. In the absence of surface charging, the quasi-neutral domain (i) and the electron cloud domain (iii) solutions can be continuously matched at a single expansion boundary (
$\unicode[STIX]{x1D6E5}=0$), permitting a 
$\unicode[STIX]{x1D6FF}$-function ion density layer and corresponding electric field discontinuity (Rhodes & Farrell Reference Rhodes and Farrell2019). With the addition of surface charging, however, a finite-
$\unicode[STIX]{x1D6E5}$ ion front boundary layer (domain (ii)) is needed in order to account for the contribution of the surface-generated electric field.
As the surface field penetrates into the ion front boundary layer, the layer width is shown to expand. An important consequence is that the ions at the front are accelerated faster than in the characteristic ambipolar plasma process. The accelerating ion front then carries more electrons toward the surface, in spite of the surface repulsion of the electron cloud. The present study is the first to establish a theoretical framework for this feedback process between surface charging and plasma expansion, predicted in simulations by Zimmerman et al. (Reference Zimmerman, Farrell, Stubbs, Halekas and Jackson2011).
Next, we describe the evolution of the surface charge and associated field, required for the time-dependent boundary layer matching.
2.3 Surface charge accumulation
The surface electric field, given by the cumulative surface charge, can be computed at each time step based on the velocity distribution of the electron cloud near the surface. In dimensionless units, the electron cloud Maxwellian distribution is written
$$\begin{eqnarray}f(v,x)=\frac{1}{\sqrt{2\unicode[STIX]{x03C0}}v_{Te}}\exp \{\unicode[STIX]{x1D719}(x)-(v/v_{Te})^{2}\}.\end{eqnarray}$$Surface charge accumulation is determined by the electron flux, given by the first velocity moment, evaluated at the surface in terms of the known solution for the electric potential
$$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D70E}_{s}}{\text{d}t}=-\int _{0}^{\infty }vf(v,x=s)\,\text{d}v=-\frac{v_{Te}}{2\sqrt{2\unicode[STIX]{x03C0}}}\text{e}^{\unicode[STIX]{x1D719}_{s}}.\end{eqnarray}$$ Here, 
$\unicode[STIX]{x1D719}_{s}=\unicode[STIX]{x1D719}(t,\unicode[STIX]{x1D701}_{s})$ is the time-dependent surface potential. The corresponding electric field above an insulated surface is given by
$$\begin{eqnarray}E_{s}(t)=-\frac{\unicode[STIX]{x1D70E}_{s}(t)}{2}=E_{s}(t_{0})+\frac{\sqrt{m_{i}/m_{e}}}{4\sqrt{2\unicode[STIX]{x03C0}}}\int _{t_{0}}^{t}\text{e}^{\unicode[STIX]{x1D719}_{s}(t)}\,\text{d}t.\end{eqnarray}$$ Here, the initial charging time is taken to be 
$t_{0}=\unicode[STIX]{x1D714}_{i}^{-1}$, accounting approximately for the initial electron expansion, as discussed earlier in this section. The initial value 
$E_{s}(t_{0})$ is determined by the surface-free electron cloud solution, expressed by (C 11). The time-dependent surface potential 
$\unicode[STIX]{x1D719}_{s}$ is not known a priori, but evaluated consistently with (2.5) and (2.9) along with the boundary layer matching conditions described below.
2.4 Boundary layer matching
The ion front boundary layer width, 
$\unicode[STIX]{x1D6E5}$, is determined by the continuity of both the electron density and the electric field at the ion front, making use of (B 8), equation (B 7) and (C 7). In addition, the electron cloud field at the ion front is characterized by the electric field at the surface, according to (C 3). The combined set of algebraic equations is as follows:
$$\begin{eqnarray}\displaystyle & n_{ef}=t^{-2}\{1-\sinh [\unicode[STIX]{x1D6E5}/t]\}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & E_{f}=t^{-1}\cosh [\unicode[STIX]{x1D6E5}/t]=\sqrt{2n_{ef}(\unicode[STIX]{x1D6E5})}\cosh [\text{arcsinh}\,(\unicode[STIX]{x1D6FC})], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & E_{s}=\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6E5})\sqrt{2n_{ef}(\unicode[STIX]{x1D6E5})}. & \displaystyle\end{eqnarray}$$ Note that (2.12) has been simplified by the assumption that 
$\unicode[STIX]{x1D701}_{s}-\unicode[STIX]{x1D701}_{f}\gg 1$, to avoid removable singularities. The system above is solved for 
$\unicode[STIX]{x1D6E5}$ by isolating the surface parameter 
$\unicode[STIX]{x1D6FC}$ in (2.11)
$$\begin{eqnarray}\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6E5})=\sinh \left[\text{arccosh}\,\left(\frac{\cosh [\unicode[STIX]{x1D6E5}/t]}{\sqrt{2}\{1-\sinh [\unicode[STIX]{x1D6E5}/t]\}^{1/2}}\right)\right].\end{eqnarray}$$ Given a value of 
$E_{s}$, the function 
$\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6E5})$ is inserted into (2.12) to find a unique solution for the boundary layer width 
$\unicode[STIX]{x1D6E5}$.
In the limit of 
$\unicode[STIX]{x1D6FC}\rightarrow 0$, where the surface field is based on the fluid floating potential (as in Rhodes & Farrell (Reference Rhodes and Farrell2019)), equation (2.12) reduces to an analytical expression
$$\begin{eqnarray}\unicode[STIX]{x1D6E5}/t=\text{arcsinh}\,(-1+\sqrt{2})\approx 0.40.\end{eqnarray}$$ This result provides the same qualitative form as in the semi-analytic interpolated model by Medvedev (Reference Medvedev2011), who found that 
$\unicode[STIX]{x1D6E5}/t\approx 1.0$ for 
$\unicode[STIX]{x1D714}_{i}t\gg 1$. Note that Medvedev (Reference Medvedev2011) is concerned with the long-time solution structure, whereas the present model addresses the initial expansion phase 
$\unicode[STIX]{x1D714}_{i}t\sim O(1)$. Owing to the constant ion density approximation which tends to overestimate 
$n_{if}$ at 
$\unicode[STIX]{x1D714}_{i}t\gg 1$, an underestimation of the present boundary layer width is to be expected.
As we show in the following section, the presence of a charge-accumulating surface causes the boundary layer width to expand faster than predicted by (2.14). As the boundary layer width increases, we expect an eventual breakdown of the thin layer approximations, namely that 
$\unicode[STIX]{x1D719}_{f}-\unicode[STIX]{x1D719}_{b}\ll 1$ and 
$n_{i}\approx \text{const.}$ within the layer. These approximations impose the following constraint on our solution validity
$$\begin{eqnarray}\unicode[STIX]{x1D6E5}/t<\text{arcsinh}\,(1)\approx 0.88.\end{eqnarray}$$Beyond the domain of this constraint, the electron cloud density vanishes (in (2.10)) and the solution becomes unphysical. This constraint limits our solution space to a minimum plasma–surface distance and a maximum expansion time.
The presented results below are strictly limited to the constraint in (2.15), upholding the thin-layer approximation 
$\unicode[STIX]{x1D719}_{f}-\unicode[STIX]{x1D719}_{b}\ll 1$. While the model validity time domain depends upon 
$\unicode[STIX]{x1D701}_{s}$, we plot only the domain that is common to all curves. For the smallest presented surface distance of 
$\unicode[STIX]{x1D701}_{s}/\unicode[STIX]{x1D706}_{D}=20$, the validity domain is found to be 
$1<\unicode[STIX]{x1D714}_{i}t<3.6$. Beyond this time the model predicts an electron cloud that becomes entirely depleted (
$n_{e}\rightarrow 0$). Another assumption that may break down at longer times is that the surface-generated field only penetrates the ion front layer, and therefore does not affect the location of the quasi-neutral expansion front. Further work is required to examine the validity of these assumptions for larger times.
3 Results
The results below provide a sample calculation for a fixed surface distance
$\unicode[STIX]{x1D701}_{s}/\unicode[STIX]{x1D706}_{D}=20$, as well as a comparison of the surface effects for different values of 
$\unicode[STIX]{x1D701}_{s}$. Our findings illustrate how the ion front boundary layer grows in time, and how the ion front is accelerated by the presence of the charge-accumulating surface.
3.1 Example with fixed surface distance $\unicode[STIX]{x1D701}_{s}/\unicode[STIX]{x1D706}_{D}=20$
Electrostatic potential spatial distribution (based on (2.2)–(2.5)) for a sequence of three times, with an insulated surface at 
$\unicode[STIX]{x1D701}_{s}/\unicode[STIX]{x1D706}_{D}=20$. The moving domain boundaries are shown for reference: rarefaction front (dashed), quasi-neutral expansion boundary (solid) and ion front (dotted).
We begin with contour plots of the electrostatic potential (figure 3) and associated charge density (figure 4), distributed throughout the region at three time slices. The potential plots are based on (2.2)–(2.5), including a matching parameter determined at each time step. The associated electric field and charge density for each region are detailed in a series of appendices; appendices A–C. Although the results of our 1-D model can be presented by simple curves, the (vertically uniform) contour plots provide a visualization of the plasma expansion and its distinct layers addressed in the following discussion. For example, the density contours shown in figure 4 are equivalent to the simple density curves shown in figures 1 and 2. The contour view places our work in the context of the 2-D plasma wake structures exhibited in the study by Rhodes & Farrell (Reference Rhodes and Farrell2019). Following the contour plots, figure 5 presents the electric field and charge density distributions in the form of curves. Each view highlights different features of the plasma expansion process.
Figure 3 shows the potential distribution, normalized by the constant electron temperature. As discussed previously the initial time is set to 
$\unicode[STIX]{x1D714}_{i}t=1$, when the self-similar solution becomes valid and the surface is set to begin charging. At this initial time the electrons have already escaped to fill the void between the plasma and the surface, leaving behind a positively charged ion front boundary layer (dotted line). At the next shown time steps, 
$\unicode[STIX]{x1D714}_{i}t=2.30$ and 
$\unicode[STIX]{x1D714}_{i}t=3.60$, we observe the ion front moving toward the surface (to the right), followed by the quasi-neutral expansion boundary (solid line). The rarefaction front (dashed line) is observed to move in the opposite direction (to the left). The surface potential, indicated by the colour at the rightmost edge of the domain, becomes more negative in time as the surface is charged by the electron cloud.
Charge density spatial distribution for a sequence of three times, with an insulated surface at 
$\unicode[STIX]{x1D701}_{s}/\unicode[STIX]{x1D706}_{D}=20$. The moving domain boundaries are shown for reference: rarefaction front (dashed), quasi-neutral expansion boundary (solid) and ion front (dotted). The densities are described by (A 4) for the self-similar domain (where 
$n_{e}=n_{i}$), equations (B 1) and (B 8) for the ion front layer and (C 8) for the electron cloud domain (where 
$n_{i}=0$).
Coupled with the potential distribution, figure 4 shows total charge density 
$n=n_{i}-n_{e}$, normalized by the undisturbed electron density 
$n_{e0}$. As in figures 1 and 2, the spatial density distributions are described by (A 4) for the self-similar domain (where 
$n_{e}=n_{i}$), equations (B 1) and (B 8) for the ion front layer and (C 8) for the electron cloud domain (where 
$n_{i}=0$). Beyond the quasi-neutral plasma on the left (white), we observe the positively charged ion front boundary layer (red) followed by the electron cloud (blue). This visualization highlights the double layer which forms at the interface of the ion front boundary layer and the electron cloud, as described in Farrell et al. (Reference Farrell, Stubbs, Halekas, Delory, Collier, Vondrak and Lin2008). The double layer is strongest at 
$\unicode[STIX]{x1D714}_{i}t=1$, when the boundary layer positive charge is most concentrated, and gets weaker over time as the positive charge layer expands. The electron density near the surface is observed to decrease in time, as electrons are absorbed by the surface, which in turn inhibits the arrival of more electrons. As noted in the previous section, at larger times our approximation of constant 
$n_{i}$ within the boundary layer may lead to an overestimation of the ion front ion density in comparison with semi-analytic expansion models (Mora Reference Mora2003; Medvedev Reference Medvedev2011). We also note that in reality the rarefaction front is also associated with a narrow layer of positive charge (Mora Reference Mora2003), which does not interact with the charging surface and is not included in figure 4.
Electric field (solid green) and total charge density (dashed yellow) spatial distribution for three time slices, with 
$\unicode[STIX]{x1D701}_{s}/\unicode[STIX]{x1D706}_{D}=20$. The unique feature is the surface electric field (rightwards of the ion front), which increases in time and depletes the electron density. The field distributions are described by (A 4) for the self-similar domain, equation (B 7) for the ion front layer and (C 7) for the electron cloud domain.
We now re-examine these results, visualized as simple curves. Focusing on a narrow region surrounding the ion front, figure 5 highlights a central element of our model; electric field continuity across the ion front (based on (A 4), (B 7) and (C 7)). The electric field (solid green) and total charge density (dashed yellow) are shown for the three time slices discussed above, 
$\unicode[STIX]{x1D714}_{i}t=1.0$, 
$\unicode[STIX]{x1D714}_{i}t=2.3$ and 
$\unicode[STIX]{x1D714}_{i}t=3.6$. The respective right-moving ion fronts are observed at 
$\unicode[STIX]{x1D701}/\unicode[STIX]{x1D706}_{D}=-0.6$, 
$\unicode[STIX]{x1D701}/\unicode[STIX]{x1D706}_{D}=2.8$, 
$\unicode[STIX]{x1D701}/\unicode[STIX]{x1D706}_{D}=8.7$, and the respective left-moving rarefaction fronts are observed at 
$\unicode[STIX]{x1D701}/\unicode[STIX]{x1D706}_{D}=-1.0$, 
$\unicode[STIX]{x1D701}/\unicode[STIX]{x1D706}_{D}=-2.3$, 
$\unicode[STIX]{x1D701}/\unicode[STIX]{x1D706}_{D}=-3.6$. Each tier of figure 5 emphasizes different features of the expansion process. The values are normalized by the undisturbed plasma electron density 
$n_{0}$ and related electric field units 
$E_{0}=(n_{0}kT_{e}/\unicode[STIX]{x1D716}_{0})^{1/2}$.
In tier (a) of figure 5, when 
$\unicode[STIX]{x1D714}_{i}t=1.0$, the dominant feature is the electron cloud, which fills the space between the plasma and the surface (
$\unicode[STIX]{x1D701}/\unicode[STIX]{x1D706}_{D}=20$). The ion front is marked by a substantial double layer jump in the total charge density, of the order of the undisturbed electron density 
$n_{0}$. Here, by construction, the nearby rarefaction and expansion fronts start out overlapping. Note that the field jump implies a 
$\unicode[STIX]{x1D6FF}$-function layer of positive charge, as discussed in Rhodes & Farrell (Reference Rhodes and Farrell2019). The rarefaction front discontinuity is of little interest to the present study, except in comparison of its magnitude with the peak at the ion front. Mora (Reference Mora2003) further explains how the initial ion front electric field peak starts out relatively small compared to the rarefaction jump – as observed in here – before obtaining the long-time asymptotic relation where they are of comparable magnitude.
In tier (b) of figure 5, 
$\unicode[STIX]{x1D714}_{i}t=2.3$, we clearly observe the characteristic electric field peak which coincides with the ion front, as expected from numerical studies (Crow et al. Reference Crow, Auer and Allen1975; Mora Reference Mora2003; Medvedev Reference Medvedev2011). This continuous electric field peak is notably absent in our previous model (Rhodes & Farrell Reference Rhodes and Farrell2019), which employed a 
$\unicode[STIX]{x1D6FF}$-function approximation for the ion front boundary layer as well as the rarefaction front. Consequentially, this image highlights the finite width of the double layer structure around the ion front, apparent in the shape of the charge density distribution.
Lastly tier (c) of figure 5, 
$\unicode[STIX]{x1D714}_{i}t=3.6$, highlights the build-up of the surface electric field, unique to the present study. The right-hand side of the field distribution, in the electron cloud domain, is observed to approach a large uniform field. This uniform field is approximately given by the vacuum solution, determined by the surface field as the electrons get pushed back into the plasma and the electron cloud density becomes small. By this mechanism the surface field penetrates the ion front and accelerates it beyond the no-surface ambipolar rate. In the context of this feedback mechanism, the ion front boundary layer can be thought of as the mediator between the ambipolar-field-dominated quasi-neutral plasma and the surface-field-dominated electron cloud region.
Additional effects of the surface charge accumulation and associated electric field are demonstrated in the following section by varying the plasma–surface distance.
3.2 Varying surface distance
Having established the general behaviour of the four-layer plasma expansion model, we now explore the surface effects by varying 
$\unicode[STIX]{x1D701}_{s}$, the initial distance between the plasma and insulated surface. We begin with the ion front of the expanding plasma, and then consider the related effects on the insulated surface.
In figure 6, we demonstrate the increased broadening rate of the ion front boundary layer width, 
$\unicode[STIX]{x1D6E5}$, as the surface is set closer to the initial plasma position. The layer width is solved at each time step based upon (2.9), (2.12) and (2.13). The layer broadening effect is caused by penetration of the built up surface field into the boundary layer. Even after a relatively short time, 
$\unicode[STIX]{x1D714}_{i}t=3.5$, the boundary layer width is observed to double in the presence of a floating surface at 
$\unicode[STIX]{x1D701}_{s}/\unicode[STIX]{x1D706}_{D}=20$, relative to the no-surface case. Here 
$\unicode[STIX]{x1D706}_{D}$ is the initial Debye length of the undisturbed plasma.
As a corollary of the accelerated 
$\unicode[STIX]{x1D6E5}$ widening, figure 7 presents the increase in ion front velocity (given by (B 9)) as the surface is brought closer to the plasma. The velocity is calculated by integrating the continuity equation within the ion front boundary layer, resulting in (B 9). At 
$\unicode[STIX]{x1D714}_{i}t=3.5$ and 
$\unicode[STIX]{x1D701}_{s}/\unicode[STIX]{x1D706}_{D}=20$, we observe a velocity boost of approximately 20 % relative to the no-surface case. It is also important to note where the model predicts the surface effects on ion front velocity to become inconsequential. For the time scale shown (
$\unicode[STIX]{x1D714}_{i}t<3.5$), we predict little surface-generated ion front accelerate at distance of 
$\unicode[STIX]{x1D701}_{s}/\unicode[STIX]{x1D706}_{D}>40$. Equivalently, farther surfaces require more time to feel the plasma coupling.
This insulated surface system is shown to be substantially different from that discussed by Gurevich & Pitaevsky (Reference Gurevich and Pitaevsky1975), where a fixed-potential surface causes a deceleration of the ion front. Furthermore, we note that our estimates of the surface effect are conservative, given that the presented model tends underestimate the ion front boundary layer width – and hence the ion front velocity – as time increases. Further study is required for a more precise calculation of the predicted plasma coupling effects of a charge-accumulating surface.
Ion front velocity versus time (given by (B 9)), shown for different distances of the charge-accumulating surface as well as the no-surface limit. Smaller surface distance is shown to increase the ion front velocity.
In addition to the effect of the surface on the expanding plasma, it is interesting to examine how this feedback process changes the conditions at the surface. Figures 8 and 9 show the temporal variations in the surface potential and near-surface electron density. These are respectively calculated based upon (2.5) and (C 8), evaluated at 
$\unicode[STIX]{x1D701}_{s}$, along with the matching (2.9), (2.12) and (2.13), solved at each time step. Each curve is scaled by its respective no-surface value, which depends upon the plasma distance even the absence of a surface. For the shortest surface distance shown, 
$\unicode[STIX]{x1D701}_{s}/\unicode[STIX]{x1D706}_{D}=20$, figure 8 shows a relative surface potential increase of over 50 % (more negative). Figure 9 shows a corresponding order of magnitude reduction in the near-surface electron density.
We note that, even for 
$\unicode[STIX]{x1D701}_{s}/\unicode[STIX]{x1D706}_{D}=40$, where there is little effect on the ion front boundary layer (according to figures 6 and 7), we still observe an order of magnitude reduction in near-surface electron density. Even at much larger distances, e.g. 
$\unicode[STIX]{x1D701}_{s}/\unicode[STIX]{x1D706}_{D}=160$, the surface conditions are found to exhibit non-negligible plasma coupling effects over a relatively short time frame. This localized effect of the surface electric field build-up is relevant to the characterization of ambient plasma conditions for explorers on the surface of lunar craters (Zimmerman et al. Reference Zimmerman, Farrell, Stubbs, Halekas and Jackson2011; Rhodes & Farrell Reference Rhodes and Farrell2019), which motivated the present study. Note that, for application to the lunar wake problem, this dynamic model must be adjusted to the steady-state surface condition. Other applications of plasma expansion theory, such as laboratory inertial fusion (True Reference True1981; Sack & Schamel Reference Sack and Schamel1987) or plasma contactors (Camporeale, Hogan & MacDonald Reference Camporeale, Hogan and MacDonald2015), may also find relevance in this general dimensionless model.
Surface potential versus time (evaluating (2.5) at 
$\unicode[STIX]{x1D701}_{s}$), shown for different distances of the charge-accumulating surface. Each curve is scaled by the no-surface curve at its respective surface–plasma distance.
Near-surface electron density versus time (evaluating (C 8) at 
$\unicode[STIX]{x1D701}_{s}$), shown for different distances of the charge-accumulating surface. Each curve is scaled by the no-surface curve at its respective surface–plasma distance.
4 Summary
We have developed a model of plasma expansion where the vacuum region is bounded by an electrically insulated surface. The model separates the plasma expansion process into four spatial domains; (i) undisturbed plasma, (ii) quasi-neutral expansion, (iii) ion front boundary layer and (iv) electron cloud. Matching the surface-generated electric field with the ambipolar expansion field determines the width of the ion front boundary layer. The model is valid during the early phase of plasma expansion (
$\unicode[STIX]{x1D714}_{i}t\sim 1$), where the surface charging effect is most notable.
As the insulated surface accumulates negative charge, it is shown to repel the electron cloud, expand the ion front boundary layer and accelerate the ion front. Even after a relatively short period 
$3.5\unicode[STIX]{x1D714}_{i}^{-1}$, a surface at 
$20\unicode[STIX]{x1D706}_{D}$ is predicted to accelerate the ion front by 
${\sim}20\,\%$ relative to the no-surface case. No less important, the model is able to predict the parameter space where surface effects do not substantially change the ion front velocity.
Another key result is the increased surface potential magnitude (which becomes more negative) and corresponding reduction in near-surface electron density, in comparison with the no-surface calculation. We find that even in the case where little surface effect is predicted in the ion front (
$t=3.5\unicode[STIX]{x1D714}_{i}^{-1}$ and 
$\unicode[STIX]{x1D701}_{s}=40\unicode[STIX]{x1D706}_{D}$), the surface electric field causes an order of magnitude reduction in the near-surface electron density.
Given the dimensionless scaling in our formulation, these substantial surface effects may be relevant to a variety of plasma systems of different scales, including solar-wind interaction with airless bodies (Rhodes & Farrell Reference Rhodes and Farrell2019), laboratory inertial fusion (True Reference True1981; Sack & Schamel Reference Sack and Schamel1987) or plasma contactors (Camporeale et al. Reference Camporeale, Hogan and MacDonald2015).
Acknowledgements
Research by D.J.R. was supported by an appointment to the NASA Postdoctoral Program at the NASA Goddard Space Flight Center, administered by the Universities Space Research Association under contract with NASA. Research by W.M.F. was supported by NASA SSERVI award DREAM2.
Appendix A. Review of self-similar formulation
This section outlines the solution of domain (ii) from figure 2; the quasi-neutral plasma expansion (bounding the trivial solution in domain (i); the undisturbed plasma). Assuming cold ions and a Maxwellian electron distribution, the combination of the ions’ equation of motion and continuity equation produces a self-similar structure, in terms of a dimensionless parameter 
$\unicode[STIX]{x1D709}\equiv \unicode[STIX]{x1D701}/t$. For a pure hydrogen plasma in the domain where 
$\unicode[STIX]{x1D709}\sim 1$, these assumptions are reasonably accurate (Denavit Reference Denavit1979; Mora & Pellat Reference Mora and Pellat1979).
Under the assumptions above, the self-similar solution in the bulk plasma is stable (Gurevich et al. Reference Gurevich, Pariiskaya and Pitaevskii1966). At the moving boundaries where quasi-neutrality down, inclusion of Poisson’s equation around the discontinuities may lead to localized short wavelength oscillations (Gurevich & Pitaevsky Reference Gurevich and Pitaevsky1975). Additionally, we note that a variety of instabilities may arise with the introduction of multiple ion species (Gurevich, Par & Pitaevski Reference Gurevich, Par and Pitaevski1973; Elkamash & Kourakis Reference Elkamash and Kourakis2016), non-Maxwellian electron distributions (True Reference True1981; Bennaceur-Doumaz et al. Reference Bennaceur-Doumaz, Bara, Benkhelifa and Djebli2015) and magnetic fields (Gurevich & Pitaevsky Reference Gurevich and Pitaevsky1975; García-Rubio et al. Reference García-Rubio, Ruocco and Sanz2016). All of these auxiliary effects are beyond the scope of the present study.
The well-known solution, with the normalization described in (2.1), is
$$\begin{eqnarray}v=\unicode[STIX]{x1D709}+1,\quad n_{e}=\text{e}^{-(\unicode[STIX]{x1D709}+1)},\quad \unicode[STIX]{x1D719}=-(\unicode[STIX]{x1D709}+1),\quad E=\frac{1}{t},\end{eqnarray}$$ where the velocity 
$v$ is defined by the direction of the expansion variable 
$\unicode[STIX]{x1D701}$. Requiring that the expansion length scale exceed the local Debye length (Samir et al. Reference Samir, Wright and Stone1983; Mora Reference Mora2003), this solution is limited to the quasi-neutral wake region where
$$\begin{eqnarray}-1=\unicode[STIX]{x1D709}_{r}<\unicode[STIX]{x1D709}<\unicode[STIX]{x1D709}_{b}=-1+2\ln t,\quad t\geqslant 1.\end{eqnarray}$$ Evaluating the self-similar solution from (A 1) at the quasi-neutrality boundary 
$\unicode[STIX]{x1D709}_{b}$, where
$$\begin{eqnarray}\unicode[STIX]{x1D701}_{b}=(-1+2\ln t)t,\quad t\geqslant 1,\end{eqnarray}$$we obtain the following boundary conditions:
$$\begin{eqnarray}v_{b}=2\ln t,\quad n_{eb}=\frac{1}{t^{2}},\quad \unicode[STIX]{x1D719}_{b}=-2\ln t,\quad E_{b}=\frac{1}{t}.\end{eqnarray}$$ Applying these conditions at the quasi-neutrality boundary, the self-similar solution is analytically continued into the non-neutral ion front layer. Note that the constant electric field in the quasi-neutral expansion zone implies an electric field discontinuity at the rarefaction boundary, since the field vanishes in the undisturbed plasma region (Gurevich & Meshcherkin Reference Gurevich and Meshcherkin1981b). The rarefaction front does not play a role in the surface interaction and is not further addressed in the present paper. The self-similar solution at the quasi-neutrality boundary 
$\unicode[STIX]{x1D709}_{b}$, however, does admit a continuous electric field if properly matched to the non-neutral region as shown in § 2.4.
Appendix B. Ion front boundary layer
This section outlines the solution of domain (iii) from figure 2; the ion front boundary layer. In this domain we approximate the ion density to be constant within a thin layer, given by
$$\begin{eqnarray}n_{i}=n_{eb},\end{eqnarray}$$ the electron density at the quasi-neutral boundary 
$\unicode[STIX]{x1D701}_{b}$. Along with the Maxwellian electrons, the potential is determined by Poisson’s equation
$$\begin{eqnarray}\frac{\text{d}^{2}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D701}^{2}}=n_{b}(\text{e}^{\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{b}}-1)\approx n_{b}(\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{b}).\end{eqnarray}$$ This form uses the normalization convention established in (2.1). The expansion above assumes a reasonably thin layer where 
$\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{b}\ll 1$. The resulting solution, requiring potential and electron density continuity at 
$\unicode[STIX]{x1D701}_{b}$, is given by
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{b}=-\frac{E_{b}}{\sqrt{n_{eb}}}\sinh [\sqrt{n_{eb}}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D701}_{b})], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D701}}=-E_{b}\cosh [\sqrt{n_{eb}}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D701}_{b})], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}^{2}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D701}^{2}}=-\sqrt{n_{eb}}E_{b}\sinh [\sqrt{n_{eb}}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D701}_{b})]. & \displaystyle\end{eqnarray}$$ The particle density, potential and electric field at the quasi-neutrality boundary 
$\unicode[STIX]{x1D709}_{b}$ are assumed to be continuous with the values in the self-similar region given in (A 4). Given 
$\sqrt{n_{eb}}=E_{b}=t^{-1}$, we obtain explicit expressions in terms of 
$\unicode[STIX]{x1D701}$ and 
$t$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{b}=-\!\sinh [(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D701}_{b})/t], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D701}}=-t^{-1}\cosh [(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D701}_{b})/t], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}^{2}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D701}^{2}}=-t^{-2}\sinh [(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D701}_{b})/t], & \displaystyle\end{eqnarray}$$ with the quasi-neutrality boundary located at 
$\unicode[STIX]{x1D701}_{b}=(-1+2\ln t)t$. Owing to the spatially constant ion density approximation (
$n_{i}=n_{eb}$) within the boundary layer, the velocity can be obtained by directly integrating the continuity equation, resulting in
$$\begin{eqnarray}\displaystyle v & = & \displaystyle v_{b}-\left(\frac{\text{d}\ln n_{eb}}{\text{d}t}\right)(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D701}_{b}),\nonumber\\ \displaystyle & = & \displaystyle 2\ln t+\frac{2}{t}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D701}_{b}).\end{eqnarray}$$ The non-neutral ion front 
$\unicode[STIX]{x1D701}_{f}$ and associated boundary layer width, 
$\unicode[STIX]{x1D6E5}=\unicode[STIX]{x1D701}_{f}-\unicode[STIX]{x1D701}_{b}$, is determined by matching the electron density (and corresponding potential) as well as the electric field with the solution in the electron cloud.
Appendix C. Electron cloud
This section outlines the solution of domain (iv) from figure 2; the electron cloud. Here, we calculate the general family of solutions for Poisson’s equation over a Maxwellian electron cloud, allowing for surface potentials greater (more negative) than the fluid floating potential. The family of solutions is parametrized by the surface charge, accumulated over time due to supra-thermal electrons penetrating the electric field.
In dimensionless form, Poisson’s equation for a Maxwellian electron cloud is written as
$$\begin{eqnarray}\frac{\text{d}^{2}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D701}^{2}}=n_{ef}\text{e}^{\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{f}},\end{eqnarray}$$subject to the boundary conditions
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}(\unicode[STIX]{x1D701}=\unicode[STIX]{x1D701}_{f})=\unicode[STIX]{x1D719}_{f}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D701}}(\unicode[STIX]{x1D701}=\unicode[STIX]{x1D701}_{s})=-E_{s}. & \displaystyle\end{eqnarray}$$ Here, 
$\unicode[STIX]{x1D70E}_{s}$ is the surface charge density of the insulated surface. The first integration yields
$$\begin{eqnarray}\frac{1}{2}\left(\frac{\text{d}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D701}}\right)^{2}-n_{ef}\text{e}^{\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{f}}={\displaystyle \frac{1}{2}}E_{s}^{2}-n_{s}\equiv n_{f}\unicode[STIX]{x1D6FC}^{2}.\end{eqnarray}$$The positive sign on the right-hand side of (C 4) above is selected to maintain the boundary conditions, subject to the constraint that the solution is real, single valued and resides in the physical domain where
$$\begin{eqnarray}-\infty <\unicode[STIX]{x1D719}<0,\quad -\infty <\frac{\text{d}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D701}}<0.\end{eqnarray}$$ To obtain the solution for 
$\unicode[STIX]{x1D719}$, we further integrate (C 4), which yields the general solution
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}(\unicode[STIX]{x1D701})-\unicode[STIX]{x1D719}_{f}=-2\log (\unicode[STIX]{x1D6FC}^{-1}\sinh [g_{\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x1D701})]), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D701}}=-\sqrt{2n_{ef}}\unicode[STIX]{x1D6FC}\coth [g_{\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x1D701})], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}^{2}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D701}^{2}}=n_{e}=n_{ef}\unicode[STIX]{x1D6FC}^{2}\text{csch}\,^{2}[g_{\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x1D701})], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle g_{\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x1D701})=\text{arcsinh}\,(\unicode[STIX]{x1D6FC})+\unicode[STIX]{x1D6FC}\sqrt{\frac{n_{ef}}{2}}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D701}_{f}). & \displaystyle\end{eqnarray}$$ Here, the function 
$g_{\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x1D701})$ has been set to satisfy the first boundary condition, equation (C 2), and the positive constant 
$\unicode[STIX]{x1D6FC}$ is to be determined by the second boundary condition, equation (C 3), for a given value of the surface field 
$E_{s}$.
In the limit where 
$\unicode[STIX]{x1D6FC}\rightarrow 0$ we regain the solution of Rhodes & Farrell (Reference Rhodes and Farrell2019), where the surface takes on the Maxwellian fluid floating potential, resulting in the following potential distribution:
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}(\unicode[STIX]{x1D701})-\unicode[STIX]{x1D719}_{f}=-2\log (h(\unicode[STIX]{x1D701})), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D701}}=-\sqrt{2n_{ef}}h(\unicode[STIX]{x1D701})^{-1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}^{2}\unicode[STIX]{x1D719}}{\text{d}\unicode[STIX]{x1D701}^{2}}=n_{e}=n_{ef}h(\unicode[STIX]{x1D701})^{-2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle h(\unicode[STIX]{x1D701})=1+\sqrt{\frac{n_{ef}}{2}}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D701}_{f}). & \displaystyle\end{eqnarray}$$This solution is equivalent to the limit where the plasma–surface distance approaches infinity, i.e. the no-surface case.
Open access
