Non-equilibrium effects in a relativistic plasma sheath model

Plasma sheaths characterized by electrons with relativistic energies and far from thermodynamic equilibrium are governed by a rich and largely unexplored physics. A reliable kinetic description of relativistic non-equilibrium plasma sheaths—besides its interest from a fundamental point of view—is crucial to many application, from controlled nuclear fusion to laser-driven particle acceleration. Sheath models proposed in the literature adopt either relativistic equilibrium distribution functions or non-relativistic non-equilibrium distribution functions, making it impossible to properly capture the physics involved when both relativistic and non-equilibrium effects are important. Here we tackle this issue by solving the electrostatic Vlasov–Poisson equations with a new class of fully-relativistic distribution functions that can describe non-equilibrium features via a real scalar parameter. After having discussed the general properties of the distribution functions and the resulting plasma sheath model, we establish an approach to investigate the effect of non-equilibrium solely. Then, we apply our approach to describe laser–plasma ion acceleration in the target normal sheath acceleration scheme. Results show how different degrees of non-equilibrium lead to the formation of sheaths with significantly different features, thereby having a relevant impact on the ion acceleration process. We believe that this approach can offer a deeper understanding of relativistic plasma sheaths, opening new perspectives in view of their applications.


Introduction
Plasmas out of thermodynamic equilibrium conditions 2 are ubiquitous as much in nature as in the laboratory and have been the focus of countless works since the very beginning of plasma physics [1][2][3][4]. Their great appeal still endures nowadays in numerous research fields, such as astrophysics [5], thermonuclear fusion [6,7], low-temperature plasmas [8], plasma propulsion [9], accelerator physics [10] and laser-plasma interaction [11,12]. An important problem common to all these fields and significant for application purposes is that of plasma sheaths [13,14]: whenever a plasma gets in contact with another physical system (e.g. a solid material, vacuum or another plasma with different macroscopic features), a charged layer-called sheath-sets up across their interface. That is because the plasma electrons usually have a much higher thermal velocity than the plasma ions (the thermal velocity is v t ∼ T 1/2 m −1/2 , where T and m are the species temperature and mass, respectively), so that the electrons will likely charge the adjacent material and leave positively charged the plasma where they come from. In this scenario, whether the plasma electrons are in thermal equilibrium or not may dramatically change the physics near the boundary interface [15], namely the density and fields profiles and their evolution in time.
Non-equilibrium plasma sheaths naturally arise in superintense laser-solid interactions, which are of great interest both for fundamental physics research [16,17] and for the development of promising applications [18,19]. In a common configuration, a high-intensity (I 10 18 W cm −2 ), ultra-short (∼1-10 3 fs) laser pulse is irradiated onto a thin ( 100 μm), solid target. During the interaction, the pulse heats a portion of the target electrons up to relativistic energies through many absorption processes [20]. We will refer to these electrons as hot electrons, as opposed to the cold electrons which do not directly partake in the interaction. In this situation, the hot electrons have a higher mobility than all the other species (i.e. ions and cold electrons) and may be able to travel through and separate from the bulk target, thus generating a sheath at the solid-vacuum interface. In this way, the interaction results in the conversion of the transverse, oscillatory electric field carried by the laser pulse into a very intense longitudinal, quasi-stationary, sheath electric field at the solid-vacuum interface. In this kind of scenario, the sheath field builds up in a relativistic, collisionless and non-equilibrium framework. Certainly, relativistic effects are crucial at such high laser intensities, as it is well-known from both simulations and experimental data [21]. This is by no means surprising since an electron in a laser field of intensity ∼ 10 18 W cm −2 can be brought to relativistic energies within a single optical cycle. Moreover, collisions can be fairly neglected provided that the plasma has a low-enough density n e and a high-enough temperature T e , so that the characteristic time between two electron collisions ν −1 e ∼ n −1 e T 3/2 e is much shorter than the interaction timescale (∼laser duration). This is definitely the case in most experimental situations, even more so as the laser intensity increases so that relativistic effects become more and more important. Besides being relativistic and collisionless, the hot electron generation is a highly nonlinear process due to many overlapping effects driven by the strong laser field. Being the plasma collisionless and the dynamics of hot electron generation ultra-fast (∼ 10 − 100 fs), then the hot electron population does not have enough time to thermalize. Hence, in that case, the assumption of thermal equilibrium is unlikely to be accurate, as is also suggested by other works [22,23].
Within the above-mentioned context, a case in point is that of laser-solid interaction for ion acceleration, especially under the scheme known as Target Normal Sheath Acceleration (TNSA) [24]. Here, the sheath field is exploited to drive the acceleration of the target ions, especially the light contaminants (protons and carbon ions) located on the rear surface as impurities. TNSA has been massively investigated both theoretically and experimentally in the past [25][26][27]. One way to describe the TNSA process relies on a quasi-static sheath model that can take into account kinetic effects in a self-consistent and relativistically-correct way [28]. Nevertheless, in some cases the theory is developed in the non-relativistic [29][30][31] or ultra-relativistic limit [32] and, what is more, in most cases the hot electron population is considered to be in thermal equilibrium [28,30,32,33].
In this paper we address the issue of assessing the role and extent of non-equilibrium effects in a relativistic, self-consistent sheath model and thereafter in the consequential TNSA process. We abandon the hypothesis of thermodynamic equilibrium of the hot electron population by describing it via a non-equilibrium distribution function f that-besides satisfying the relativistic Vlasov equation in presence of a self-consistent electromagnetic potential-depends on an additional parameter α which measures the 'degree of non-equilibrium': the larger α is, the further away from a Jüttner distribution (e.g. in terms of L 2 distance) f is. This modeling approach allows one to directly explore the role of non-equilibrium by letting α vary, something that is not easily done by means of conventional numerical simulations, e.g Particle-In-Cell (PIC) simulations [34]. Among possible reasonable choices for f, we choose a relativistic generalization of the Cairns distribution function [35], also with the goal of extending a previous non-relativistic treatment proposed to model the effects of energetic electrons on TNSA [31]. In section 2 we introduce all the ingredients needed to define our approach and we show that the proposed Cairns-like distribution function fully satisfies the principles of special relativity by writing the main expressions in a manifestly Lorentz-covariant form. Once the distribution function is defined, in section 3 we make use of it in a stationary sheath model and we define a strategy to assess the effects of non-equilibrium. Lastly, in section 4 we show that the Cairns-like function is a reasonable choice to describe the hot electron population generated by super-intense, ultra-short laser-plasma interaction, and we apply these results to TNSA modeling in order to obtain the main quantities of interest, e.g. the accelerating field and the ion maximum energy.

Non-equilibrium, relativistic distribution function
Quantitatively assessing the importance of non-equilibrium of the relativistic hot electron population requires to make a specific choice for its distribution function. We define a scalar distribution function f = f α (x μ , p μ ), which depends upon the time-space coordinates x μ = (ct, x), the energy-momentum four-vector p μ = (p 0 , p) and has an additional dependence on a parameter α 0, such that Here f eq is the equilibrium distribution function in presence of an electromagnetic four-potential A μ = (φ, A) [36][37][38], namely a generalization of the Jüttner distribution: with m e the electron mass, c the speed of light in vacuum, n 0 the unperturbed electron number density (i.e. if A μ = 0), K n the modified Bessel function of the second kind of order n (here n = 2), T the temperature, ζ = m e c 2 /T, e the elementary charge and U μ = γ(u)(c, u) the hydrodynamic four-velocity of the plasma (u is the three-dimensional hydrodynamic velocity and γ(u) its Lorentz factor). Note that, following the argument developed by Hakim [36,38], obtaining (1) as the explicit solution of the collisionless Boltzmann equation in presence of an electromagnetic field that also makes the collisional integral disappear requires two additional assumptions. First, that the (timelike) four-vector U μ does not depend on the space-time variables (and so u is constant in both space and time). Second, that A μ is constant along the direction identified by U μ . The latter hypothesis implies that U μ ∂A ν /∂x μ = 0, which means that the electromagnetic scalar and vector potentials are time-independent in the frame of reference where U μ = (c, 0), known as the Lorentz rest frame. Here, we keep the same hypotheses even if we are looking for a non-equilibrium distribution. Therefore, we will find a non-equilibrium distribution function that is time-independent in the Lorentz rest frame as well. Additionally, f α (x μ , p μ ) together with A μ must solve the Vlasov-Maxwell system, given by (in Gauss units): where we introduced the electromagnetic tensor F μν = ∂A ν ∂x μ − ∂A μ ∂x ν and the current four-vector J μ = (cρ, J). We further assume that any plasma population other than the one described by f (in our case the hot electrons) is part of a cold background, so that the dependence of J μ upon f can be expressed as follows: where Ω ⊆ R 3 is a suitable integration domain and γ(p) = 1 + |p| 2 /(m e c) 2 is the electron Lorentz factor and i is an index that runs through all the cold populations with charges q i and number densities n i . In order to readily satisfy all these requirements, we pick a distribution function that is a function of the variables x μ , p μ only through a constant of motion I, i.e. a quantity such that dI/dτ = 0, where τ is the proper time. In our case, an integral of the motion is U μ p μ − e c A μ [36], that written in the Lorentz rest frame reduces to the three-dimensional single particle Hamiltonian At this point, we need to adopt an explicit form of the distribution function. There are several reasonable options, such as the κ or Vasyliunas distribution function [39][40][41][42], the Tsallis distribution function [43][44][45], a super-Gaussian [46,47], a Maxwellian with a supra-thermal component [48] or a Jüttner distribution with a superimposed shift in the momentum [49]. Here, we choose a Cairns-like distribution function. The Cairns distribution function was first introduced to explain soliton structures observed in cold space plasmas [35] and then employed in a non-relativistic TNSA model [31]. We extend the Cairns distribution function [35] to the relativistic regime, which, in a manifestly covariant form, reads: where N is the normalization factor: The term m e c 2 /T in (6) is added to let f α reduce to the Cairns distribution when taking the non-relativistic limit. The main reason behind the choice of the Cairns distribution lies in its simplicity.. First, f α is built as a sum of two contributions: f eq and a 'perturbation' term proportional to α. This allows us to interpret α as a knob that moves the distribution function from equilibrium (α = 0) to a higher degree of non-equilibrium (α → ∞). Second, when plugging f α into the sheath model, computations in the non-relativistic [31] and ultra-relativistic limits can be partially carried out analytically. In addition, as we will detail in section 4, f α is consistent with 3D PIC simulations of TNSA.
Having written the chosen non-equilibrium distribution function in a manifestly covariant fashion, the hot electrons are automatically described in a relativistically-correct way. A corresponding three-dimensional expression can be retrieved by fixing the frame of reference as the Lorentz rest frame, which leads to: Note that f α (x, p) in (8) only depends on (x, p) through the above-mentioned Hamiltonian H(x, p), thus guaranteeing that it satisfies df/dt = 0. In the Lorentz rest frame, the system of equations (2) and (3) reduces to the familiar electrostatic Vlasov-Poisson system: The expression (8) can be simplified for instance by assuming that all the quantities depend on only one space coordinate, say x, and only one momentum coordinate, say p x . Under these additional hypotheses, the distribution function (8) can be written in the Lorentz rest frame as: where the normalization factor now reads Figure 1 shows the dependence of the distribution defined in equation (8) upon the absolute value of the momentum |p|. On the left, figure 1(a) shows how f changes for varying α, when the electrostatic potential vanishes. As α increases, two symmetric supra-thermal populations arise, while the cold portion of electrons (that around |p| = 0) is reduced. The value α = 1 separates two regimes: one with no maximum points other than zero (α < 1) and one with two additional maximum points (α > 1). At the (pointwise) limit α → ∞, for which f is well-defined, f is that of a two counter-propagating electron bunches with most probable energy given by 2 √ 1 + ζ/ζm e c 2 . On the right, figure 1(b) shows again f defined in equation (8) as a function of |p|, but having fixed α = 1 and for different values of the potential φ. For increasing φ, f consequently increases everywhere. This is because the normalization factor N defined in equation (7) cannot depend on φ, or otherwise f would not satisfy Vlasov equation. This implies that f, as written in equations (8) and (10), is normalized to n 0 only for φ ≡ 0.
Lastly, it is interesting to compute the number density n and kinetic energy density e kin of the hot electron populations described via (10): where ϕ = eφ/m e c 2 . Both n and e kin are functions of the space coordinate x through the potential φ = φ(x) and both depend on the non-equilibrium parameter α. By an additional integration in space, once φ(x) is known, one can compute the total number of electrons N tot and their total energy E tot , which will again depend on α. This implies that two distribution functions of the form (10) that only differ by their α value, will describe two populations with not only different 'non-equilibrium parameter', but also different total electron number and total electron energy. It is clear that a comparison between such populations is not fair if one intends to investigate the importance of non-equilibrium only. Therefore, in order to isolate non-equilibrium effects, one has to find a way to fix N tot and E tot even when varying α. As a consequence, some other quantities (e.g. T, n 0 ) will have to properly vary together with α. We devised a specific procedure to accomplish this goal in the framework of the sheath model for TNSA. We describe the details of our approach in the next section.

Non-equilibrium in a plasma sheath model
We mean to determine the plasma sheath generated by a relativistic hot electron population described by the non-equilibrium distribution function of the form (10). Here, we detail the rationale of our modeling approach, going through the main hypotheses and equations. We consider a one-dimensional plasma consisting of a hot electron population and a cold background. We assume that the hot electron population has already been generated through other previous processes (e.g. laser-plasma interaction) and has reached a stationary non-equilibrium state described by (10). The non-equilibrium feature is attributed to the prior phenomena and to a low collisionality of the relativistic plasma. Under these assumptions, the physics is governed by the one-dimensional version of the Poisson-Vlasov system (9). We describe the build-up of a stationary (in the rest frame), one-dimensional sheath field at the plasma-vacuum interface using an electrostatic and self-consistent approach. Not all the electrons contribute to the sheath formation, rather only those electrons whose kinetic energy does not exceed the potential energy that they experience (we call them trapped electrons). To do so, we adopt the sheath model described in details in references [28,30,32].
From now on, we will use dimensionless quantities in the following units: momenta in m e c, energies and temperatures in m e c 2 , densities in n 0 , lengths in (R e n 0 ) −1/2 (R e = e 2 /m e c 2 is the classical electron radius) and f in n 0 (m e c) −1 . We use the potential energy, rather than the potential itself and define the dimensionless quantity ϕ = eφ/(m e c 2 ). In addition, we use a gauge such that ϕ vanishes at +∞. We remark that these dimensional units, even if seldom adopted, are quite convenient because they do not hide any dependence on the temperature, as it would happen if distances were normalized with respect to the Debye length.
The main equation of the model is the Poisson equation written piecewise inside (x 0) and outside (x 0) the plasma: where n trap (x) is the number density of the trapped electrons and x * < 0 is a point well-inside the plasma where the charge density vanishes by construction, so that local neutrality holds. n trap (x) is given by: where f = f α (x, p) is the Cairns-like distribution (10) that was discussed in section 2. We report it here in dimensionless units for convenience (we now write p instead of p x ): The following 'implicit' boundary conditions are taken. It is assumed that there exist two points x * < 0 and x > 0 such that: where the value ϕ * is a free parameter of the model. Additionally, the continuity of the electric field and of the potential is imposed at the solid-vacuum interface, which implies that there is no charge accumulation on the plasma-vacuum interface. Equation (17) is a relationship of the form ϕ 0 = ϕ 0 (ϕ * ), so that the value of the potential at the interface can be found (as a function of ϕ * ) without having to solve the whole problem. Indeed, this relation represents the actual boundary condition that is used to numerically solve both equations of the system (13) after one spatial integration. To this regard, we point out that the values of the two points x * andx are not determined a priori with the boundary conditions (16). Instead, once the relation ϕ 0 = ϕ 0 (ϕ * ) has been fixed, they can be found by integrating twice equations (13) and expanding up to the first order around ϕ = ϕ * inside the plasma and ϕ = 0 outside the plasma. It results that x * → −∞, so that the plasma region, in this model, is formally semi-infinite, whilex < +∞, so that the sheath field has a finite extension. Equations (13)- (17) mathematically formalize the physical sheath model, which has four degrees of freedom: ζ, n 0 (hidden in the units), ϕ * and α. Figure 2 shows a sketch of the typical solution of the model together with the main quantities used to frame the whole problem. Now, we wish to understand to what extent and how the main quantities are affected by the Cairns-like non-equilibrium. To tackle this issue, the easiest way would be to fix the parameters ζ, n 0 and ϕ * , hence to solve the model for different values of α. However, this procedure does not isolate the effect of non-equilibrium alone because varying α strongly affects the total hot electron number N tot and energy E tot . Indeed, figure 3(a) shows N tot (blue) and E tot (red) as functions of α for fixed ζ, n 0 and ϕ * . Switching the non-equilibrium on (α 1) makes them rapidly increase, while they reach an asymptotic value for α 1. Overall, N tot and E tot can vary by an order of magnitude. In light of this trend, we rather fix the total number N tot and the total energy E tot of the hot electron population, given by: where Σ ⊆ R and A, B, C, D and E are now given by: Because the model is one-dimensional, N tot and E tot are actually the total number and total energy per unit surface area (in units of #e/ R e /n 0 and m e c 2 / R e /n 0 , respectively). Note that we approximated the total electron energy with its kinetic component, assuming that the plasma is ideal. The simplest choice for the integration domain is Σ = [−x d ,x], where x d > 0 is representative of the size of the plasma. Hence, we fix the physical size of the plasma d once and for all, as if we are studying one only plasma configuration, so that x d = d/(R e n 0 ) −1/2 . The only requirement on x d is that it should be large-enough so that the potential approaches its boundary value ϕ * . We choose to keep N tot and E tot constant, equal to their value at equilibrium (α = 0), given by: These expressions of N tot and E tot can be approximated by considering the following observations. First, the typical potential profile, solution of the model, is approximately constant inside the target, where it takes its maximum value, ϕ * . Second, the potential monotonically (≈exponentially) drops to zero outside the target. If it is also possible to assume that ζϕ * 1, then we can approximate the integrals in (19) using the following: which yields These steps let us write N tot and E tot as functions of the free parameters only, at least for the equilibrium case. Note that ζϕ * 1 is a very reasonable hypothesis: indeed ϕ * and ζ −1 represent the maximum and average trapped electron energy, respectively. Figure 3

(b) shows a very good agreement between formulas (19) (solid lines) and their approximations (21) (dashed lines).
In order to keep the values of N tot and E tot constant, equal to their values at equilibrium for every α, we adjust the other parameters (ζ, n 0 and ϕ * ) conveniently for each α. We illustrate the technical details of our approach in appendix A. Now that we have identified a strategy to assess the effects of non-equilibrium, we apply it to the case of laser-ion acceleration driven by the TNSA process.

Non-equilibrium in a laser-driven plasma sheath
We apply the approach discussed so far to describe the sheath field that drives the ion acceleration process in a TNSA scenario with the goal to assess the magnitude of non-equilibrium effects. Now, before going into the details of our results, we discuss the physical meaning of the main hypotheses presented in section 3 within the context of TNSA. First of all, since the very beginning we have been using a static picture, meaning that the electron distribution function in the Lorentz rest frame does not depend on time. This implies that our results will describe the physics on a time frame for which the hot electron population has already been generated, while the ions are still cold and immobile. Second, we assume that f depends on only one space coordinate, x, which is the direction normal to the target. This is reasonable if the laser spot size is large enough ( laser wavelength), so that a one-dimensional (in space) picture can apply. As a consequence, the sheath field always arises along the target-normal direction-even in case of non-normal incidence-as it is well-understood from the features of the accelerated ions [50,51]. The description developed in previous sections is also one-dimensional in the momentum space, meaning that the distribution function depends on only one momentum coordinate. Under the TNSA scenario, it is reasonable to assume such remaining momentum component to be the longitudinal one, p x . This approximation is the more reliable as the electron Lorentz factor γ = 1 + |p| 2 is better approximated with its longitudinal part γ ≈ 1 + p 2 x . In order to check the consistency of these hypotheses in a TNSA scenario, we performed a 3D PIC simulation representative of a standard experimental configuration. This also allows us to highlight those physical features that can be well described by the stationary 1D sheath model, even if the actual physics is time-dependent and 3D.

Comparison with a 3D particle-in-cell simulation
We performed a 3D PIC simulation where a 400 nm-thick solid target is irradiated at normal incidence by a 0.8 μm wavelength, 30 fs, P-polarized laser with a 0 = 4. The relevant details of the simulation are reported in appendix B. In figure 4(a) we show the time evolution of the kinetic energy associated to the electron (green) and ion (purple) populations expressed as a fraction of the initial laser energy (i.e. the total energy of the system). The laser energy is quickly absorbed by the electrons within a few tens of fs (time t = 0 being the beginning of the interaction), until the absorbed energy reaches a peak value ( 10%) while the ions are still at rest. Then, the electrons transfer part of the absorbed energy to both bulk and contaminant ions, of which the latter start to be accelerated under the action of the sheath field. This confirms that the early stages (t < 100 fs) of the acceleration process can be satisfactorily described by means of a static TNSA model. We have also collected the electron distribution function at t = 53 fs (marked by a black vertical line in figure 4(a)) at the plasma-vacuum interface (i.e. at the target rear surface). Figure 4(b) shows the distribution function on the transverse momentum plane (p y , p z ). It can be appreciated that this distribution is symmetric, only depending on the magnitude of the transverse momentum p ⊥ = p 2 y + p 2 z , even in the case of a P-polarized laser. Hence, we compare the distributions along the p x (black) and p ⊥ (blue) axes in figure 4(c). While for relatively cold electrons (p ⊥ , p x < 1m e c) the distributions are similar, once p x , p ⊥ exceed √ 3m e c (to which corresponds a kinetic energy of m e c 2 ) the distribution is dominated by the longitudinal component. This implies that, for hot electrons, the average γ = 1 + |p| 2 ≈ 1 + p 2 x , according to the model hypotheses. Moreover, it can be seen that a clear non-equilibrium feature only appears in the longitudinal momentum distribution. By looking at its shape, three regions can be recognized: a cold part (p x < √ 3m e c) with a steep decreasing slope, a hot part ( √ 3m e c < p x < 3m e c) where the slope becomes milder, and a bump (3m e c < p x < 6m e c).
The distribution in longitudinal momentum is reprinted with black dots in figure 4(d). The cold electron population can be satisfactorily fitted with a Jüttner distribution (blue line, T cold ≈ 40 keV) [29]. On the other hand, the distribution cannot be described in terms of a Jüttner function (which is monotonically decreasing in p x ) at higher energies. This clearly points to some non-equilibrium process. Indeed, we fitted the PIC data for the hot electrons (p √ 3m e c) with either a Jüttner distribution (red line, T hot ≈ 0.2 MeV) and a relativistic Cairns-like distribution (orange line, T Cairns ≈ 0.5 MeV, φ fit ≈ 1.1 MV, α ≈ 4.6). Both functions are in good agreement with the data in the range 1.5m e c < p x < 3m e c. However, only the Cairns-like distribution is able to reproduce the PIC data for p x > 3m e c. A satisfactory fit over the whole spectrum is given by the sum of the cold Jüttner and the Cairns-like distributions (purple dashed line), indicating that this can be a suitable choice to account for non-equilibrium effects in laser-generated relativistic plasma sheaths.
In order to further assess whether the relativistic non-equilibrium sheath model can actually describe a TNSA scenario, we solved it using the Cairns-like distribution function obtained from the fit shown in figure 4(d) (α = 4.6 and T = 0.5 MeV). The fitted distribution function is evaluated at the target-vacuum interface, so that the fitting value φ fit is an estimation of φ(0) = φ 0 . Hence, we fixed the parameter φ * in such a way that φ 0 ≈ φ fit , while n 0 is left as a free parameter. The resulting electric field and potential obtained with n 0 = 5 × 10 20 cm −3 are shown in figures 4(e) and (f), respectively. By comparing the green (PIC) and magenta (model) curves it can be appreciated how the model nicely reproduces the simulation results.

Consequences on target normal sheath acceleration
After having assessed the soundness of the relativistic sheath model based on the Cairns-like distribution to describe TNSA, we investigate the role of non-equilibrium on the TNSA process itself.
Here, we follow the idea of estimating the maximum ion energy through the potential drop at the target-vacuum interface by invoking total energy conservation in an electrostatic picture. Assuming that the ions accelerated by the sheath field are test particles located on the surface x = 0, then the maximum energy that they will acquire is given by Zϕ 0 m e c 2 [32]. Now, ϕ 0 can be written in terms of the parameters ζ, ϕ * and α by expanding equation (17), as was done in reference [28] for the equilibrium case: where β(ϕ) = (1 + ϕ) 2 − 1, In equation (22) the (normalized) maximum ion energy is written as the sum of two contributions: ϕ * − 1/ζ, which does not depend on α, and a combination of other terms that do depend on α. Among the latter, the term proportional to (1 + 2α) would vanish if all hot electrons were trapped (note that, at equilibrium, if all electrons were trapped, i.e. if ζϕ * 1, then ϕ 0 → ϕ * − 1/ζ). On the other hand, the term proportional to α couples the effects of non-equilibrium and trapped electrons, since it would survive even if all electrons were trapped. Now, to numerically solve the model, we need to set the free parameters to values that are relevant in the context of TNSA. To accomplish this goal, it is of particular interest to relate the formal mathematical model to a realistic experimental configuration, suitably described by specific choices of the free parameters. Overall, there are the 4 degrees of freedom (α, n 0 , ζ, ϕ * ) and we wish to study the dependence of the solutions on one of them-α. In order to relate the mathematical model to possible experimental configurations, we express the other free parameters (n 0 , ζ, ϕ * ) in terms of the laser and target properties (e.g. the normalized intensity a 0 , energy E L , spot size σ and the target thickness d), using the following scaling laws: Relation (23) is the well-known ponderomotive scaling for the electron temperature [52], while (24) is an empirical scaling based on experimental results, proposed in reference [32]. Additionally, we introduce equation (25) as a way to relate n 0 to the other two parameters ζ and ϕ * . Here, the laser energy E L is assumed to be converted into kinetic energy of the hot electrons E tot with efficiency η < 1, so that E tot (m e c 2 / R e /n 0 )σ = ηE L . Formally, this leads to (25), where we used the expression of E tot in (21) and the fact that E tot is actually an energy per unit surface (which is why the spot size appears in the equation). We fix the conversion efficiency to ∼ 10%, in agreement with the above-mentioned PIC simulation (see figure 4(a)). Now, once the laser and target configuration is given, we use equations (23), (24) and (25) to find a set of reference parameters. Then we apply the procedure detailed in appendix A in order to find the adjusted values of the parameters n 0 , ζ, φ * that keep N tot and E tot constant for varying α. From now on we fix the plasma thickness to 5 μm and the laser parameters to the following: 0.8 μm wavelength, 30 fs duration and 25 μm 2 spot area. We let the laser energy vary in the range 0.5-10 J. Figure 5(a) shows how the maximum ion energy ϕ 0 varies with α (the parameters of the model are those shown in figure A1). Overall, using this approach, the maximum proton energy can increase by up to ≈40% its value at equilibrium, suggesting that non-equilibrium features of the hot electron population can play a major role in TNSA ion acceleration. Figure 5(b) shows the minimum and maximum predictions for ϕ 0 , obtained with α = 0 (blue) and α = 1 (red) respectively, as functions of the laser energy E L . It can be seen that the scaling with respect to E L is similar for both values of α. The shaded area between the two curves corresponds to all the possible model outcomes in terms of maximum proton energy for every α ∈ [0, ∞).
In figure 6 we show the results of the model obtained with 1 J laser energy. When varying the non-equilibrium parameter α, the spatial profiles of the (a) electrostatic potential, (b) electric field and (c) electron number density change significantly. The further away from equilibrium, i.e. the higher α, the larger the extension of the sheath. This is a consequence of how a higher degree of non-equilibrium rearranges the electrons in the phase space (x, p x ) with respect to the equilibrium case. To this regard, figure 7 shows the phase space (x, p x ) of the hot electrons obtained with 1 J laser energy and three different non-equilibrium parameters: α = 0, 0.1, 1. The superimposed white lines mark the locus of points for which an electron carries as much kinetic energy as the potential energy it experiences: all electrons within this region are trapped. Far away from the target-vacuum interface, namely where the potential is low (ϕ 1), the larger α is, the more electrons are able to escape the potential, so that the potential drop needs  to fall upon a larger spatial region. This leads to a more extended sheath as shown in figure 6. Of course, a larger sheath does not necessarily accelerate the contaminants ions more efficiently. The acceleration of contaminant ions is more efficient as the sheath potential is higher and larger, but not independently. In our scenario, the value α = 1 might be identified as the optimal non-equilibrium parameter. Indeed, lower values of α lead to a lower potential barrier and to a less extended sheath. On the other hand, larger values of α, although generating a longer sheath, lead to lower potential barriers such that the overall result is to accelerate the contaminant ions to lower maximum energies. All in all, α = 1 realizes a convenient trade-off between the height of the potential and its extension. Lastly, note that enabling the non-equilibrium leads to significant differences even for α 1, where small variations of α yield significant variations in the main sheath and acceleration properties, such as the spatial profiles and maximum ion energy. This indicates that even small deviations from equilibrium may be important and, ideally, should be taken into account when modeling the sheath or when interpreting experimental results in light of the considered model. Moreover, even if measuring the electrons distribution function is far beyond experimental capabilities, the electric field and electron density profiles are accessible observables, for example by means of proton radiography [53][54][55]. To this regard, our results offer an interpretative key to relate available experimental observables to a theoretical model, which could help assessing the actual importance of non-equilibrium.

Conclusions
In this paper we have addressed the topic of electron non-equilibrium in the formation of a relativistic sheath. This is of interest in the context of laser-solid interactions, especially when driving the acceleration of ions under the TNSA scheme. We have introduced an explicit form for the distribution function f α of the electron population that generates the sheath, condensing the degree of non-equilibrium in a single parameter α. Moreover, we have written f α in a Lorentz-covariant form in order to manifestly show that it complies with special relativity. The proposed distribution function is also consistent with 3D PIC simulations of laser-solid interaction, which reinforces the soundness of our choice. Electron populations that have the same main macroscopic properties (total number, total energy, hence average energy as well), but different degrees of non-equilibrium, cause the formation of sheaths with significantly different features. These results indicate that the non-equilibrium kinetic properties of the electron population can play an important role and that their macroscopic properties may not be not enough to determine the main sheath quantities and the maximum energy of the ions accelerated via TNSA. Moreover, we find that even a small degree of non-equilibrium can have an important impact on the main observables, i.e. the electric field and the electron density. We stress that, albeit having picked a specific non-equilibrium function, our modeling approach is general. On the one hand the proposed strategy is by no means limited to a specific choice of the distribution function; on the other hand we have determined the model parameters in a closed way, namely without having to introduce any other arbitrary or fitting parameter.
In conclusion, we have identified a new approach to model non-equilibrium in the self-consistent generation of a relativistic sheath. This is a strongly cross-disciplinary topic, being of interest in a variety of fields, such as astrophysics, accelerator physics, nuclear fusion and advanced laser-driven ion acceleration, especially when using non-conventional targets. Our results provides a key for the interpretation of those numerical simulations and experiments that allow to retrieve information on the sheath dynamics.
Here we report the details of our strategy to find the free model parameters n 0 (in the units), ζ and ϕ * that let N tot and E tot constant when varying α, so that we can study the effect of non-equilibrium only.
As a first option, one could obtain an analytical expression for N tot and E tot defined in equations (18) as a function of the free parameters only by adopting the same approximation strategy as done at equilibrium to obtain equations (21) (essentially, evaluate the integrand at ϕ * and multiply by x d ). In this case one has to solve simple algebraic equations to tune the parameters that let N tot and E tot constant, equal to their value at equilibrium. However, using this approach we observed fluctuations in the unapproximated N tot and E tot of the order of ∼ 10%. Hence, we prefer to proceed numerically without any approximation (other than those due to the numerical computations) rather than analytically. We compute N tot and E tot a posteriori, once the model has been solved and the potential profile ϕ(x) is known, via numerical integration of equations (18). Hence, we adopt the following strategy. First, we pick a set of reference parameters n 0 , ζ, ϕ * and fix the product ζϕ * once and for all. We solve the model at equilibrium using these parameters to find N tot (n 0 , ζ, ϕ * , 0) and E tot (n 0 , ζ, ϕ * , 0). Then, for all α = 0, we use the following two-step procedure: • Solve the model with the reference parameters to find N tot (α) and adjust the reference value n 0 to its updated value n 0 (α): n 0 (α) = n 0 N tot (n 0 , ζ, ϕ * , α) N tot (n 0 , ζ, ϕ * , 0) • Solve the model with the adjusted parameter n 0 (α) for many ζ to find the dependence of E tot on ζ; then, in order to find ζ(α) (and so ϕ * (α) according to the fixed product ζϕ * ), solve the implicit equation in the variable ζ: E tot (n 0 (α), ζ, ϕ * , 0) = E tot (n 0 (α), ζ, ϕ * , α).
In this way one obtains the new set of parameter n 0 (α), ζ(α) and ϕ * (α) to be used as inputs in the model when α = 0. Figure A1 shows how the free model parameters n 0 (panel (a)), T (b) and ϕ * (c) need to vary with α, respectively, in order to keep N tot and E tot constant within 1% accuracy (the other parameters are the same as those mentioned in the main text, namely 5 μm thickness, 0.8 μm wavelength, 30 fs duration and 25 μm 2 spot area). According to the chosen distribution function, increasing α implies two main consequences (compare with figure 1): increasing the number of hot electrons at high kinetic energies (p 1) where the potential is low (ϕ ∼ 0); increasing the number of hot electrons at low kinetic energies Figure A1. Main results of our approach to isolate non-equilibrium features in the sheath model. Panels (a), (b) and (c) show how the main free model parameters n 0 , T and eϕ * need to vary with α, respectively, in order to keep the total number N tot and total energy E tot of hot electrons constant within 1% accuracy.
(p ∼ 0) where the potential is high (ϕ 1). Since the potential, solution of the problem considered here, is 1 in a significant portion of its domain and N tot depends exponentially on ϕ(x), then it is easy to see that increasing α would increase N tot . The first step of the procedure compensates this effect by decreasing the parameter n 0 . At this point, simultaneously increasing α and reducing n 0 causes E tot to have a minimum for α ≈ 1. Hence, the second step of the procedure compensates this behavior by increasing ζ −1 up to α = 1 and decreasing it for α > 1. The behavior of the parameter ϕ * is simply a consequence of having fixed the product ζϕ * .

Appendix B. Particle-in-cell simulation
We performed an illustrative 3D particle-in-cell simulation of laser-plasma interaction to support the rationale of our model together with its hypotheses. We used the open source, massively parallel code piccante [56]. In this simulation a λ-wavelength laser pulse hits a thin solid target that is comprised by a bulk layer and a thin contaminant layer on non-illuminated side. The box size was [100λ × 60λ × 60λ] with a resolution of 40 points per wavelength in each direction. The code adopts the Yee algorithm as Maxwell solver with a time resolution given by 0.98 of the Courant condition. The laser is longitudinally cos 2 with duration (full width half maximum of the field) equal to 15λ/c, while it is transversely Gaussian with a waist of 4 μm. The normalized peak vector potential is a 0 = 4. The laser propagates along the x direction and is P-polarized along the y direction. The target, irradiated at normal incidence, is a 0.5λ foil with 40n c electron density and Z/A = 1/2. Its back surface is located at x = 0. The electron population is sampled with 40 macro-electrons and 2 macro-ions per cell. The contaminant layer is 0.05λ-thick with Z = A = 1, 7n c electron density, sampled with 64 macro-electrons and macro-ions per cell. The macro-electrons were initialized with a 5.11 eV temperature, while the ions are initially cold. In our computation we always consider λ = 0.8μm, so that the laser has a 30 fs duration. The interaction with the tail of the pulse begins at time 0λ/c, while the pulse peak hits the front surface at 19λ/c. To obtain the distribution functions shown in figure