Fast ion acceleration in a foil plasma heated by a multi-picosecond high intensity laser

Laser-matter interactions in the relativistic regime, where the laser intensity exceeds 10¹⁸ W/cm², have opened up various innovative applications such as intense x-ray and neutron sources, compact particle accelerators, and fast ignition-based laser fusion. Recently, high power lasers with kilo-Joule (kJ) energy such as LFEX,¹ NIF-ARC,² LMJ-PETAL,³ and OMEGA-EP⁴ have been operated, which enable to investigate the relativistic plasma interactions with multi-picosecond (ps) pulse durations.

In the ps relativistic laser-matter interactions, the plasma on the front surface, i.e., the laser-irradiated surface, expands significantly, which changes the laser absorption and hot electron dynamics. The plasma expands with the density ranging from a few percent to near-critical density in a spatial scale of the order of 10 μm at the front side of the target, and forms a large scale corona plasma. In this corona plasma, the intense laser accelerates/heats electrons through various mechanisms such as stimulated Raman back scattering,⁵ the wake field acceleration due to forward Raman scattering and/or self-phase modulation, the betatron acceleration in laser-plasma channels,⁶ stochastic acceleration due to the overlapping of incident and reflected laser pulses,^7,8 and so on. Recently, long time laser interactions toward 10 ps with thick overdense plasmas have been studied. It is demonstrated by particle-in-cell (PIC) simulation⁹ and also by experiment¹⁰ that the electron slope temperature becomes much higher than the conventional ponderomotive scaling¹¹ during the multi-ps interactions. Several new electron acceleration mechanisms in a multi-ps regime have been proposed. One is the electron acceleration by the interplay between the 2ω_L ponderomotive force and sheath electrostatic potential, where ω_L is the laser angular frequency.¹² The other is that the self-generated magnetic fields on the front surface play an important role in electron acceleration.¹³ 

In interactions of thin solid foils with a relativistic laser in the regime of $IλL2≥1018 W μm2/cm2$⁠, where I is the laser intensity and λ_L is the laser wavelength, the target-normal sheath acceleration (TNSA) of ions can take place at the rear side of the foil.^14–17 In the TNSA, ions are accelerated by the sheath electric field created at the expansion front of the foil plasma, and the characteristic expansion speed of the ions is represented by the sound velocity $Cs=(ZTe/Mi)1/2$⁠, where Z is the ion charge state, T_e is the electron temperature, and M_i is the ion mass.^18–22 Hence, the key parameter of the TNSA is the electron temperature T_e resulting from the laser-plasma interaction. Conventionally in the sub-ps regime, the ponderomotive temperature, $Tp≡mec2(1+a02/2)1/2$⁠, has been used to explain the electron temperature observed in experiments.¹¹ Here, m_e is the electron rest mass, c is the speed of light, and $a0=eE0/mecωL$ the normalized laser field amplitude, where e is the fundamental charge and E₀ is the amplitude of laser electric field. Note that in the conventional self-similar model of TNSA, T_p is usually assumed to be determined by the peak value of the laser pulse amplitude $â0$⁠, i.e., the electron temperature is constant in time. However, in the multi-ps regime, such an isothermal TNSA model will not be applicable since the change of the electron temperature from the ponderomotive scaling is unignorable.

Electron heating beyond the ponderomotive temperature takes place also in interactions of thin solid foils with high-contrast femtosecond (fs) short-pulse lasers. In this case, the interaction of the strong laser fields at a rigid foil surface with little plasma expansion can be realized, which yields the stochastic vacuum heating.²³ It is also known that when the laser pulse duration t_p and foil width L satisfy the relation 2L < ct_p, relativistic electrons are accumulated by the recirculation between the front and rear surfaces of the foil, since electrons are trapped by the sheath electric fields.²⁴ As a result, the hot electron density in the foil plasma increases. Such a high energy absorption from laser to hot electrons in thin foils leads to an efficient TNSA for ions.²⁵ 

Recently, an ion acceleration experiment using kJ, ps high intensity laser LFEX and thin aluminum (Al) foil target of thickness L = 5 μm was conducted.²⁶ As the pulse duration was extended from 1.5 ps to 6 ps while keeping the peak laser intensity constant around $2.5×1018 W/cm2$ and laser spot size with a full width at half maximum (FWHM) of 70 μm, the proton maximum energy was increased beyond the one-dimensional (1D) isothermal TNSA model. It is observed that the slope temperature of hot electrons increased with pulse duration up to 3 ps and saturated for a longer pulse. In this system, plasma expansion during the multi-ps pulse duration and electron recirculation around the foil of the width L ≪ ct_p/2 are expected to incorporate. Note that the above recirculation is effective when the laser spot size is large enough in comparison with the plasma expansion distance.

In this paper, we investigate the multi-ps laser interaction with a solid thin foil plasma in 1D geometry based on the PIC simulation, and revisit the plasma expansion theory to study the ion acceleration including the electron temperature evolution during the multi-ps laser irradiation. In Sec. II, we perform PIC simulations for interactions between a solid foil plasma and multi-ps laser pulses which are similar to the LFEX experiment conditions.²⁶ We here obtain the time evolution of the electron temperature. Based on the simulation results, in Sec. III, we develop a quasineutral plasma expansion theory that takes the time dependence of electron temperature during the laser irradiation into account. By using the developed theoretical model, we derive the ion maximum energy achieved in the interaction and compare it with the simulation results. Conclusions are given in Sec. IV.

We investigate interactions between a thin foil plasma and intense laser fields with the pulse durations of ps and multi-ps based on the relativistic one-dimensional PIC simulation using EPIC code.²⁷ We consider linearly polarized laser fields assuming two pulse functions, i.e., a Gaussian pulse function $f=exp (−(t−t1)2/tp2)$ with t_p = 0.9 ps and t₁ = 1.2 ps, and a longer pulse that rises and falls with Gaussian functions $f=exp (−(t−t1)2/tp2)$ and $f=exp (−(t−t2)2/tp2)$⁠, respectively, while having a flattop region f = 1 for t₁ ≤ t < t₂ where t₂ = 2.7 ps, as shown in Fig. 1(a). The FWHM of the single Gaussian pulse is 1.5 ps. The flattop pulse is assumed to represent a pulse train consisting of two Gaussian pulses. Hereafter, we refer the above Gaussian and flattop pulses to as 1 pulse and 2 pulse cases, respectively. The laser wave length is λ_L = 1.05 μm and the peak normalized amplitude is $â0=1.42$ which corresponds to the intensity of $I=2.5×1018 W cm−2$⁠. These parameters correspond to the LFEX experiment.²⁶ 

The length of the simulation box is L_x = 154 μm and 256 μm in the x direction for the 1 and 2 pulse cases, respectively, with the mesh size of 10 nm. The laser field is excited by an antenna located at the left boundary. We set a fully-ionized uniform plasma of L = 5 μm thickness with initial ion and electron densities of $nib=4.0×1022 cm−3$ and n_eb = Zn_ib, respectively, from the position x = 75 μm, as shown in Fig. 1(b). On the front side of the target, we put a 1 μm thick preformed plasma, the ion density of which increases linearly from 0 to $nib$⁠. The corresponding electron density n_eb = Zn_ib is n_eb = 40 n_c where the charge state is Z = 1 and n_c = 1.0 × 10²¹ cm⁻³ is the cutoff density for the incident laser light. The charge-to-mass ratio of the ions is Z/A = 1/2. The initial temperatures of ions and electrons are T_i0 = 0.2 keV and T_e0 = 1 keV, respectively. Collisions and ionization processes are not included in the simulations performed in this study.

Figures 2(a)–2(c) show the temporal evolution of the electron energy distribution obtained in the simulation for the 2 pulse case at (a) t = 1.2 ps, (b) 2.4 ps, and (c) 3.6 ps. The vertical axis denotes the distribution function normalized by the total counted number of electrons N_e, i.e., $dNe/dεe/Ne$⁠, where ε_e is the electron energy. Note that in obtaining the energy distributions, we observed electrons in the rear side region, i.e., from x = 81 μm to the rear boundary x = L_x, so that bulk cold electrons in the solid foil are not included in the spectrum. The absolute value of the electrostatic potential energy at the ion expansion front in the rear side $−eϕf$ at each time is shown by vertical black dashed lines. Here, we set Φ = 0 at the rear surface of the target, x = 80 μm. The electrons that come out from the target rear with energy below $−eϕf$ are reflected by the sheath field of the foil, while energetic electrons whose energy at the target rear exceeds $−eϕf$ can go beyond the ion front. Since the ion acoustic speed related to the plasma expansion depends upon the bulk hot electron temperature, we determine the electron temperature that is effective for the plasma expansion, T_e, by fitting the energy spectrum below $−eϕf$⁠. We assume the 1D relativistic Maxwellian as the fitting function and obtain T_e = 0.3 MeV, 0.8 MeV, and 1.5 MeV for t = 1.2 ps, 2.4 ps, and 3.6 ps, respectively, as shown by the blue dotted lines.

In Fig. 2(d), the time evolution of electron temperature for 1 and 2 pulse cases are shown in red and blue points, respectively. Pulse shape functions for 1 and 2 pulse cases are shown on the top of the figure for reference. We see that for both cases, the temperature increases in time exceeding the ponderomotive temperature for the corresponding peak intensity, i.e., $Tp=mec2(1+â02/2)1/2=0.2 MeV$⁠. The time evolution of T_e is fitted by the following exponential function:

Here, $T̃e=Te/T0$ is the normalized temperature, T₀ is the initial temperature at time $t=t0, tmax$ is the time at which T_e takes the maximum value $Tmax$⁠, and $α≡Tmax/T0$⁠. From the initial condition $T̃e(t0)=1$⁠, the time scale t_T is obtained as $tT=(tmax−t0)/(lnα)1/2$⁠. The initial temperature T₀ is assumed to be the ponderomotive temperature T_p. We fit the parameters α and $tmax$ in Eq. (1) by the simulation results. The time t₀ is set to be t₀ = 1 ps, at which T_e ∼ T_p is satisfied in both 1 and 2 pulse cases in the simulation. The fitted functions are shown by solid lines in Fig. 2(d). We find α  =  2.3 and 7.2 for 1 and 2 pulse cases, respectively, which represent that $Tmax=0.46 MeV$ in the 1 pulse case and $Tmax=1.43 MeV$ in the 2 pulse case. This tendency of increase of the temperature is similar to that obtained in the LFEX experiment given in Ref. 26 where the slope temperatures of electron are measured as T_h = 0.45 MeV and 1.10 MeV for 1 and 2 pulse cases, respectively.

In Figs. 3(a) and 3(b), we present the snapshots of plasma densities, electrostatic potential $ϕ$⁠, and longitudinal electric field in the 2 pulse case at t = 3.46 ps, around which the electron temperature exhibits the maximum value as seen in Fig. 2(d). The laser pulse still continues at this time with the amplitude at around half maximum. Ion and electron charge densities, normalized by the cutoff density n_c, are shown in Fig. 3(a) by the red line and in (b) by the blue line, respectively. At this time, the ion front is positioned at around x = 143 μm as pointed by the red arrow in Fig. 3(a), while electrons distribute exceeding the position of the ion expansion front. The charge separation generates the sheath electric field E_x as shown by the black line in Fig. 3(b), whose peak moves with the ion expansion front as pointed by the black arrow. Note that the field E_x shown in (b) represents the spatially averaged value obtained from the neighbouring 30 meshes which corresponds to the length of 0.3 μm. The value of E_x is normalized by the pulse peak amplitude of the laser electric field $Ê0$⁠. The front peak value of $Ex/Ê0$ at the position of the black arrow is twice larger than the spatial average of $Ex/Ê0$ in Fig. 3(b) in the rear side where the electrostatic potential $ϕ$ starts to drop (from x = 100 μm to 130 μm). This relation of E_x fields between inside of and at the expansion front of the plasma is consistent with the result in Ref. 21. On the other hand, at the front side x ∼73 μm, a sharp electron density jump is formed by the laser radiation pressure to create a charge separation and the strong E_x field. In Fig. 3(a), the electrostatic potential $ϕ$ is shown by the black line, which we obtained by integrating the longitudinal electric field E_x along x from x = 75 μm to the –x direction in the front side and from x = 80 μm to the +x direction in the rear side. The potential is assumed to be zero at x = 75 μm and 80 μm. We see that the potential at the ion expansion front reaches $ϕf∼−7 MV$⁠. Therefore, the electrons with energy below 7 MeV are confined by the potential $ϕf$⁠.

Figures 3(c)–3(e) represent electron distributions in phase space x-p_x for times t = 0.98 ps, 2.52 ps, and 3.46 ps, respectively. The red dotted line corresponds to the position of the ion expansion front at each time. As seen in Fig. 3(c), electron bunches are launched from the front side in the 2ω_L frequency. In Figs. 3(d) and 3(e), it is seen that during the multi-ps interaction, most electrons are trapped by the plasma sheath and recirculate around the target, resulting in the distribution which is almost symmetric with respect to the p_x axis. The electrons directing to the rear side of the plasma lose their momentum during passing through the rear sheath potential. A small number of electrons go beyond the ion expansion front, i.e., beyond the red dotted line, which causes the breaking of charge neutrality of the bulk plasma. Notice that density of the hot electrons at x = 80 μm, i.e., the starting position of plasma expansion, n_e0, is around 3.5 n_c which is higher than the relativistic cutoff density γn_c ∼ 1.4 n_c. We also confirmed that n_e0 increases from the subcritical density in time up to t = 2.4 ps and saturates at n_e0 ∼ 3.5 n_c. This represents that the hot electrons accumulate around the foil due to the recirculation as is also mentioned in Ref. 24. The expansion front moves about 30 μm in the 1 ps from (d) t = 2.52 ps to (e) t = 3.46 ps, so that the expansion speed is about 10% of the speed of light c. In Fig. 3(f), the trajectory of a selected electron in phase space is shown for the period t = 0 ps–2 ps (blue line), for t = 2 ps–3 ps (light blue line), and for t = 3 ps–4 ps (black line). The inversed triangle and a circle indicate the position at times t = 0 ps and 3.46 ps, respectively. The right side axis shows the kinetic energy associated with p_x, i.e., $(γx−1)mec2$⁠, where $γx≡(1+px2/(mec)2)1/2$⁠. We see that the electron recirculates between the rigid front surface x = 73 μm and the expanding rear plasma region during t = 0 ps–3 ps with increasing its maximum energy by the multiple kick from the laser field at the front side. In t = 3 ps–4 ps, the electron gains a large momentum after excursion at the peripheral of the front surface, i.e., y = 67 μm–73 μm, and the recirculation distance increases drastically. In this phase, the electron is accelerated to the energy above 1.5 MeV. Such recirculating electrons contribute dominantly in creating the expanding plasma tail.

To see the electron behavior in detail, the trajectories in t–x space and time-energy space of the same electron in Fig. 3(f) are shown in Figs. 4(a) and 4(b), respectively. The black thin dashed-dotted lines at x = 74 μm and 80 μm represent the initial surface positions of the plasma. The recirculation length in the x direction increases gradually in time due to the expansion of the plasma and increase of the electron energy. During the laser pulse time, the electron recirculates more than 20 times. Kinetic energy of the same electron shown in Fig. 4(b) reaches its maxima when the electron is kicked by the fluctuating electric potential and the ponderomotive force of the laser on the front surface. The trajectory in (a) and energy in (b) for t = 2.2 ps–2.6 ps are closed up in (c) and (d), respectively. We see in Fig. 4(d) that the electron energy fluctuates when the electron is positioned on the front surface, e.g., during t = 2.54 ps–2.55 ps. During this time, the electron is trapped by a potential well near the rigid cutoff surface as indicated by Ref. 28. In t = 1.2 ps–3.3 ps in Fig. 4(b), the maximum energy of the recirculating electron varies randomly. This indicates that the phase of the acceleration field is stochastic when the electron is injected into the laser-plasma interaction region.²⁹ Therefore, as seen in Fig. 4(d), the electron energy after the kick at the front side decreases (t = 2.27 ps), increases (t = 2.39 ps), or keeps the same value before the kick (t = 2.55 ps). This stochasticity originates from the non-harmonic nature of the electron recirculation. At t = 3.1 ps–3.3 ps [shaded area in (a) and (b)], the electron energy oscillates rapidly between 0 to 1 MeV around the mean position of x = 70 μm. This is due to the electron interaction with laser and fluctuating electrostatic fields in the underdense plasma. In such an underdense regime, electron interaction with incident and reflected laser waves leads to stochastic heating,⁸ which contributes in generating the hot electron components beyond the ponderomotive temperature. In Fig. 4(e), the phase plot of the same electron during t = 3.1 ps–3.3 ps is shown. Note that this figure corresponds to a close-up picture of Fig. 3(f). As indicated by the left arrow in Fig. 4(e), the electron enters into the underdense region and suffers from a strong acceleration at x = 73.5 μm to p_x ∼0. However, the electron continues to move toward the –x direction with a fluctuating momentum. After a complex excursion, the electron is finally injected back to the foil plasma. After gaining an energy of about 1.6 MeV at t = 3.3 ps, the electron reaches to the position of x = 109 μm at t = 3.46 ps where the electron energy becomes almost zero. Here, it is worth noting that the absolute value of the potential $ϕ$ at x = 109 μm at time t = 3.46 ps is 1.6 MV, as shown in Fig. 3(a). After t = 3.5 ps, i.e., after the laser peak, in Fig. 4(a), the turning point of the electron in the front side moves gradually outward, i.e., to the negative x-direction, due to the decrease of the strong laser radiation pressure. The electron continues recirculation with a long time period in the adiabatically-expanding plasma after the laser pulse.

One can also see in Fig. 4(a) that the frequency of the recirculation increases after t = 1.2 ps and gradually decreases from around t = 2.2 ps. The frequency increases when the velocity of the recirculating electron increases, and decreases when the recirculation length increases as the plasma expands. In the end, when the recirculation period becomes comparable with the laser pulse time, the effect of recirculation on the electron heating will be inefficient. Note that Fig. 4(b) indicates the energy gain and the loss of the recirculating electron. Before t = 3.4 ps, the electron energy increase due to the multiple laser heating with reflection by the rear-side sheath similar to the “Fermi acceleration,” while after t = 3.4 ps, the electron loses the energy due to the adiabatic plasma expansion as the “Fermi deceleration,” i.e., adiabatic cooling during multiple reflection by the sheath.

We here consider a theoretical model of 1D plasma expansion that can include the time evolution of electron temperature which we obtained in the PIC simulation in Sec. II. We assume that a plasma occupies the area of x ≤ 0 in 1D geometry at initial time t = 0 with ion and electron densities n_i0 and n_e0 = Zn_i0, respectively. The expansion of the plasma into the vacuum region x > 0 can be described by the equation of motion and the equation of continuity for ions,

together with the force balance equation for electrons

and also with the Poisson equation

Here, v_i is the ion velocity, $ϕ$ is the electrostatic potential, and we neglected the ion pressure term in the equations of motion assuming cold ions T_i ∼ 0. We assume the quasineutral condition n_e = Zn_i instead of solving the Poisson equation Eq. (5). Note that the quasineutrality breaks up at the expansion front x = x_f where the Debye length $λD=(Te/(4πnee2))1/2$ exceeds the scale length of plasma expansion, i.e., time integral of the sound velocity C_s.²¹ In the isothermal case, the condition for breaking of the quasineutrality can also be represented by the ion plasma frequency ω_pi as ω_pit < 1, where $ωpi=(4πnie2/Mi)1/2$⁠.

In this study, we assume that the electron temperature appearing in Eq. (4) is time dependent, i.e., T_e = T_e(t), which is different from the isothermal condition where T_e is assumed to be constant. Equations (2)–(4) are the basic equations describing the 1D plasma expansion with time-dependent electron temperature.

From Eq. (4), the Boltzmann distribution is obtained as usual but with the time-dependent temperature as

Here, although the coefficient n_e0 can be a function of time, we assume n_e0 = const. for simplicity. As we discussed in Sec. II C, n_e0 increases in time due to the recirculation. However, the change of n_e0 in time is smaller than that of the temperature T_e. For the region where the quasineutrality is satisfied, i.e., x < x_f, we obtain the ion density from Eq. (6) as

In the case of isothermal expansion, the solution that satisfies Eqs. (2)–(4) and the quasineutral condition n_e = Zn_i is known to be described by a single variable, which is defined by

where the sound velocity C_s is constant. In this case, the solution for n_i, v_i, and $ϕ$ is self-similar where no characteristic spatial and time scales appear. On the other hand, in this study, in order to consider the time-dependent electron temperature and the corresponding sound velocity C_s(t), we introduce a new variable on the analogy of Eq. (8) as

where R is defined by

which corresponds to the distance of plasma expansion during time t with a time dependent sound velocity. Note that for the isothermal case C_s = const., Eq. (9) reduces to Eq. (8) with R = C_st.

We here consider a transformation of variables $(x,t)→(ξ,t)$⁠. Derivative operators are transformed as $∂x|t=R−1∂ξ|t$ and $∂t|x=∂t|ξ−ξCsR−1∂ξ|t$ where each derivative is executed keeping the subscript variable written on the right of the vertical line unchanged. Hereafter, we do not write the subscripts explicitly. In the new coordinate (ξ, t), we consider that n_i, v_i, and $ϕ$ are functions of t and ξ to be solved from the basic equations [Eqs. (2)–(4)], while T_e and then C_s are given as functions of t.

Under the transformation of variables $(x,t)→(ξ,t)$⁠, Eqs. (2) and (3) are transformed to

where the dot represents time derivative d/dt. Note that $Ṙ=Cs$⁠. Equations (11) and (12) are the simultaneous equations to solve v_i and $ϕ$⁠. The ion density n_i is then determined according to Eq. (7).

We solve Eqs. (11) and (12) on the basis of perturbation expansion assuming that the time scale of temperature evolution $(d(log Te)/dt)−1=(2Ċs/Cs)−1$ is slow compared to that of plasma expansion R/C_s. Namely, we introduce a dimensionless parameter b as

where ε is a smallness parameter, and expand Eqs. (11) and (12) up to the first order of ε. Considering that each Eqs. (11) and (12) has a term that contains $R∂t$⁠, we expand the variables as

where $vi′(1)$ and $ϕ′(1)$ are functions of t and ξ, and have no explicit dependence on $R=R(t)$⁠. Hereafter, we use the notations $vi(1)≡Rvi′(1)$ and $ϕ(1)≡Rϕ′(1)$⁠. The ion density is expanded as $ni=ni(0)+εni(1)+⋯$⁠.

On the basis of the above ordering in ε, the zeroth order terms in Eqs. (11) and (12) are given by

and the first order terms are given by

Note that in the zeroth order of ε, T_e, and C_s are constant as is in the isothermal case, and thus $R=Cst$⁠. Therefore, we assumed $∂tvi(0)=0$ and $∂tϕ(0)=0$ in Eqs. (16) and (17), considering that there is no explicit time dependence in both Eqs. (11) and (12). In the first order of ε, $vi(0)$ and $ϕ(0)$ are treated as functions of t and ξ with time dependent C_s and T_e.

The solution of the zeroth order equations Eqs. (16) and (17) is obtained as

The corresponding electrostatic field E⁽⁰⁾ is derived as

which is proportional to $R−1∝t−1$ and uniform in x. The ion density is then derived from Eq. (7) as

Equations (20), (21), and (23) are the self-similar solutions for the conventional isothermal plasma expansion.¹⁸ 

In deriving the first order solution, we substitute the zeroth order solutions Eqs. (20) and (21) to the first order Eqs. (18) and (19). Derivatives of the zeroth order variables in Eqs. (18) and (19) are calculated as $∂ξvi(0)=Ṙ, ∂tvi(0)=R¨(1+ξ), e∂ξϕ(0)=−Te$⁠, and $e∂tϕ(0)=−Ṫe(1+ξ)$⁠. Equations (18) and (19) then lead to

where the relation $b=RR¨/Ṙ2$ is used. Finally, the solution up to the first order of ε is obtained as

Equations (26)–(28) represent the ion velocity, static potential, and ion density, respectively, for quasi-neutral 1D plasma expansion with time-dependent electron temperature under the first order approximation with respect to $b∼O(ε)$⁠. The corresponding electrostatic field is given as

where the coefficient Θ, which includes the first-order correction, is given by

Note that Θ  =  1 for the zeroth order of ε, i.e., for b = 0. The first-order correction means that the plasma scale length is shorter by a factor of 1 + b/2 than the time integration of the ion acoustic speed, R, when one assumes an increasing temperature, i.e., $Ṫe>0$ and thus b > 0. Hence, during the pulse time where the electron temperature increases, the electrostatic field is larger by 1 + b/2 than that obtained in the zeroth-order approximation, i.e., Eq. (22).

We here calculate the expansion length R and the expansion parameter b for the temperature evolution obtained in the PIC simulation in Sec. II B. We obtain the expansion length R by integrating the sound velocity C_s, which is proportional to T_e^1/2, in time according to Eq. (10) as

Here, we normalized the expansion length as $R̃=R/(Cs0t̂)$ where $Cs0=(ZT0/Mi)1/2$ is the sound velocity for the initial temperature T₀, $τ=t/t̂$ is the time normalized by $t̂≡(2eN)1/2ωpi0−1, ωpi0$ is the ion plasma frequency for the density n_i0, i.e., $ωpi0=(4πni0e2/Mi)1/2$⁠, and t_T is the time scale of temperature evolution given in Eq. (1). The integrant in Eq. (31), $C̃s=Cs/Cs0$⁠, is the normalized sound velocity. In Fig. 5(b), we plot R given by Eq. (31) for 1 pulse and 2 pulse cases by red dashed-dotted and blue solid lines, respectively. Here, we assumed the temperature function given in Eq. (1) assuming α = T_max/T₀ = 2.3 (red dashed-dotted lines) and α  =  7.2 (blue solid lines), where T₀ is defined to be the ponderomotive temperature T_p = 0.2 MeV. Note that the above values of α correspond to the 1 pulse and 2 pulse cases of PIC simulation, respectively. The corresponding temperature functions are shown in Fig. 5(a). We also show the isothermal case where T_e = T₀ = const. by black dashed lines. Pulse shape functions for 1 and 2 pulse cases are shown on the top of the figures for reference. We see for the isothermal case (black dashed line) that R increases linearly with time, since the time integration of C_s is merely equivalent to C_s0t. On the other hand, when the temperature is time-dependent, R has a finite curvature and deviates from the linear dependence on time. For the temperature function Eq. (1), the ordering parameter b is obtained as

Notice that b becomes zero at t = 0 and $t=tmax$⁠. This is because the expansion length R is zero at the initial time t = 0 and the time-derivative of the sound velocity $Ċs$ is zero at the temperature peak $t=tmax$⁠. In Fig. 5(c), we plot b given by Eq. (32) for the temperature function in Fig. 5(a). At first, the value of b increases as the temperature grows, and begins to decrease from t ∼ 0.8 ps and 1.6 ps for 1 pulse and 2 pulse cases, respectively. The parameter b becomes zero at the time $t=tmax$⁠, i.e., t = 2.1 ps and 3.6 ps for 1 and 2 pulse cases, respectively, since the temperature function is flat at $t=tmax$ as seen in Fig. 5(a). We see that the absolute value of b is around 0.5 at maximum up to the time t = 2.5 ps and 4.2 ps in 1 pulse and 2 pulse cases, respectively, which correspond to the time after the pulse peak. Therefore, the present perturbation analysis is considered to be applicable during the laser irradiation. Note that for the isothermal case T_e = const., the parameter b is zero as shown by the black dashed line.

In Fig. 5(d), we plot Θ given in Eq. (30). For the isothermal case, Θ  =  1 for all the time as shown by the black dashed line. In non-isothermal cases (red and blue), Θ becomes larger than unity during the laser pulse time where the electron temperature increases. Θ  >  1 indicates that the electric field Eq. (29), is larger than that for the zeroth order approximation Θ  =  1. After the pulse time, Θ decreases toward –∞. This is due to the fact that b on the RHS of Eq. (30) decreases to –∞ at the corresponding time, which denotes that the present perturbation analysis becomes not applicable after that time.

We here calculate the maximum energy of ions that are accelerated by the sheath field generated by the multi-ps scale plasma expansion. The first order solution for the plasma expansion with time-dependent electron temperature [Eqs. (26)–(29)] is obtained for the region ω_pi0t ≫ 1 where the quasi-neutral assumption is applicable. At the ion expansion front where the charge neutrality breaks, it is known that the electric field becomes twice larger than that in the quasi-neutral region.²¹ We take the same assumption for the electric field at the ion expansion front E_f, i.e., E_f is twice larger than Eq. (29), as

As is discussed for Fig. 3(b) in Sec. II C, the front electric field observed in the PIC simulation is about two times larger than the electric field inside the expanding plasma. Note that in the zeroth order of ε, i.e., for b = 0, Θ is unity, and thus we obtain E_f = 2E⁽⁰⁾ as in Ref. 21.

On the other hand, the electric field at time t = 0 is obtained by integrating the Poisson equation, Eq. (5), from x = 0 to ∞ assuming the Boltzmann distribution for electrons, Eq. (6), and n_i = 0 for ions as

where $eN=exp (1)$⁠, and λ_D00 denotes the Debye length for density n_e = n_e0 and initial temperature T_e = T₀. In order to see the acceleration of the expansion front in whole time of the laser interaction, one can consider continuing the electric field Eq. (33) to that at the time t = 0, Eq. (34), as discussed in Ref. 21. We here define the electric field

This satisfies $Efa∼E(t=0)$ [Eq. (34)] for R/λ_D0 ≪ 1 which corresponds to ω_pi0t ≪ 1 and E_fa ∼ E_f [Eq. (33)] for R/λ_D0 ≫ 1 which corresponds to ω_pi0t ≫ 1. Equation (35) thus represents the electrostatic field which is valid for all t.

Now, by using the global expression Eq. (35), we can obtain the expression for the front velocity that is valid for any time t. We integrate the equation of motion at the ion front, i.e., $Midvf(t)/dt=ZeEfa(t)$⁠, where E_fa on the RHS is given by Eq. (35), as

For the isothermal case where $T̃e=1, R̃=τ$⁠, b = 0, and Θ = 1, the integral in Eq. (36) is solved analytically as $vf(0)=2Cs0ln(τ+1+τ2)$⁠, which is consistent with Eq. (10) in Ref. 21. Since the ion velocity is larger for large ξ as indicated by Eq. (26), the maximum energy of ions can be obtained from the front velocity as $εimax=Mivf2/2$⁠, namely,

Equation (37) gives the time evolution of ion maximum energy in non-isothermal conditions in the first order approximation with respect to ε.

We integrate Eq. (37) numerically assuming the time dependent temperature given by Eq. (1), where the parameters α and $tmax$ are obtained from the fitting of the PIC results as in Sec. II B. The ion species is assumed as Z = 1 and Z/A = 1/2. The resultant ion maximum energy for 1 pulse and 2 pulse cases is shown in Figs. 6(a) and 6(b), respectively, by red bold solid lines. Here, $εimax$ depends on the plasma density at the origin of the expansion, i.e., n_e0 defined in Eq. (6), through $t̂∝ωpi0−1∝(ne0)−1/2$⁠. We show here two results assuming n_e0 = 0.2 n_c and 4 n_c in the integration by the lower and upper red bold solid lines, respectively. This is according to the PIC result that the hot electron density n_e0 increases in time from subcritical density ∼0.2 n_c to 3.5 n_c as discussed in Sec. II C. The two lines assuming 0.2 n_c and 4 n_c thus correspond to the lower and higher limits of the theoretical prediction for the ion maximum energy, respectively. The shaded areas in Fig. 6 represent the allowable areas for the ion maximum energy for the density n_e0 varying from 0.2 n_c to 4 n_c. The black circles show time evolution of the ion maximum energy observed in the PIC simulations. We also plot Eq. (37) with Θ  =  1 by blue dashed lines, and the isothermal theory where the front velocity is given by $vf(0)=2Cs0ln(τ+1+τ2)$ by black lines. Here, the constant ponderomotive temperature T_p = 0.2 MeV is assumed for the black dashed-dotted lines, i.e., T_e = T_p = const., and the averaged temperatures ⟨T_e⟩ = 0.25 MeV for 1 pulse case and ⟨T_e⟩ = 0.68 MeV for 2 pulse case are assumed for the black thin solid lines. Here, ⟨T_e⟩ are the averaged values of the temperature for times from t = 0 ps to 4 ps in the 1 pulse case and from t = 0 ps to 7 ps in the 2 pulse case. Note that the integration of Eq. (37) assuming Θ  =  1 (blue dashed lines), which corresponds to the zeroth order theory with respect to ε, is different from the isothermal model since the time dependent temperature is used for T_e in the integrant in Eq. (37).

For both 1 pulse and 2 pulse cases, most of $εimax$ of the simulation is seen inside the allowable area of the first order non-isothermal theory until the laser pulse is almost finished. We see that the PIC result follows the first-order theory of n_e0 = 0.2 n_c initially, and it shifts to the first-order theory of n_e0 = 4 n_c later. This reflects the time evolution of n_e0 due to the accumulation of hot electrons by recirculation. The difference of zeroth-order and first-order of the non-isothermal theory becomes obvious in the 2 pulse case where the zeroth order theory predicts smaller values of $εimax$ in the latter half of the pulse time. The difference between the zeroth-order and the first-order originates from the change of the value of Θ from unity in the first order as we have seen in Fig. 5(d). On the other hand, the isothermal theory exhibits a different time evolution with no inflection point. The PIC results are seen outside the allowable area of the isothermal models for both T_e = T_p and T_e = ⟨T_e⟩ during most of the time. Note that after the time when the laser incidence is almost finished, the first order theoretical lines are changed to dotted lines. This is because the absolute value of expansion parameter b increases rapidly toward unity, which corresponds to the limit of applicability of the present perturbation analysis, after t = 2.5 ps and 4.2 ps in 1 pulse and 2 pulse cases, respectively. As indicated by the vertical dashed-dotted-dotted lines in Fig. 6, the adiabatic theory is to be applied after the time where the energy input by laser is almost finished.^30–33 Therefore, the non-isothermal theory obtained here can be regarded as a theory that mediates between isothermal and adiabatic regimes. Note that in the 1 pulse case, $εimax$ of the PIC simulation in the adiabatic regime is seen inside the allowable area of the isothermal models. This denotes that for the TNSA accelerations with picosecond laser pulses, the final value of $εimax$ may be predicted by the conventional isothermal model using a time-averaged electron temperature instead of the ponderomotive temperature, if one defines a final time of acceleration as the empirical model in Ref. 34. However, for multi-ps laser interactions, the isothermal model is not applicable even for the prediction of the final maximum energy of ions as indicated by Fig. 6(b).

We also investigated the dependence of $Tmax/Tp=α$ on the laser peak amplitude $â0$ for a given pulse duration. We confirmed that the ratio of the maximum temperature T_max and the ponderomotive temperature T_p is almost constant, α  = 6.4 ± 1.0, in the intensity regime from $I=1.2×1018 W/cm2 (â0=1.0)$ to $1.0×1019 W/cm2 (â0=2.84)$ for the 2 pulse case. Therefore, $Tmax$ has the same dependence on $â0$ as the ponderomotive temperature but with a different factor α ∼6.4 in the 2 pulse case. Note that α is likely to increase with respect to the pulse duration, since the ratio depends on the recirculation, namely, the number of interactions with the laser. The function shape of T_e(t) and the peak time $tmax$ are similar among the different $â0$⁠. Here, under the approximation $(1+R̃2)1/2∼R̃$ and Θ ∼1, the integrant on the RHS of Eq. (37) becomes $T̃e/R̃∝α1/2$⁠, which is independent of $â0$⁠. Therefore, the dependence of the ion maximum energy $εimax$ on $â0$ is the same as that of the ponderomotive temperature when we assume T₀ = T_p in Eq. (37). In the intensity regime from $â0=1.0$ to 2.84, we obtained a scaling $Tmax≃6.1 Tp$ for the 2 pulse case. By using this scaling, we can estimate that $εimax$ reaches to the level of 100 MeV for the laser irradiation with the peak amplitude of $â0≥3.2$ and the pulse duration of t_p ≥ 3 ps. Note that to achieve such an energetic ion acceleration by the long pulse TNSA, not only the laser intensity (or amplitude) but also the spot size have to be increased in order to sustain the quasi 1D condition during the faster plasma expansion by the acoustic speed with higher electron temperature.

Finally, we show the time evolutions of energy partition in the system in Fig. 7 for electrons (blue solid), ions (red dash-dot), and the sum of electrons and ions (black bold). The plotted values are the ratio of each energy to the total energy input by the laser pulse. Figures 7(a) and 7(b) correspond to the 1 and 2 pulse cases, respectively. The pulse shape function for each case is shown on the top of the figure for reference. In both 1 and 2 pulse cases, the electron energy increases initially, and then it is transferred to the ions during the plasma expansion phase. After the pulse peak, the sum of the particle energies shown by the black solid line saturates which indicates that the expansion is adiabatic after the laser incidence. Namely, the electron thermal energy is converted to the ion expansion energy adiabatically. The saturation times are t ∼ 2.2 ps and 3.8 ps for 1 pulse and 2 pulse cases, respectively, which are consistent with the adiabatic regime we estimated in Fig. 6. The energy partition rate of ions increases from 1 pulse case, i.e., 6%, to 2 pulse case, i.e., 9%. This tendency is consistent with the recent experiment by the LFEX laser²⁶ where the energy conversion efficiency to ions increases when the laser pulse duration is extended from 1.5 ps to 6 ps.

We studied the one-dimensional expansion of a thin foil plasma irradiated by a high intensity laser with multi-picosecond pulse durations by using PIC simulation. During the interaction, the hot electron temperature increases beyond the ponderomotive scaling. The key mechanism leading to such an electron heating is found to be the recirculation of electrons around the expanding foil plasma. An electron trajectory in the plasma shows that electrons are accelerated stochastically in the front side region by fluctuating electrostatic and laser ponderomotive forces many times during the recirculation. Such electron heating associated with the recirculation process becomes important in multi-ps time scale laser-foil interactions.

Based on the numerical results, we developed a plasma expansion theory that takes the time dependence of electron temperature into account in one-dimensional geometry. By assuming that the time scale of electron temperature evolution is slow compared with the plasma expansion time scale, we derived a non-self-similar solution in the first order approximation. We further obtained a formula for the time evolution of ion maximum energy applicable in the non-isothermal situations. By substituting the time evolution of electron temperature obtained in the PIC simulation, the time evolution of ion maximum energy obtained by the theory agrees well with that observed in the PIC simulation until the laser pulse incidence is almost finished. Therefore, the non-isothermal theory obtained here can be regarded as a theory that mediates between isothermal and adiabatic regimes.

So far, various advanced models, based on the isothermal plasma expansion theories, have been studied for describing complex laser-plasma interactions, e.g., the two temperature model assuming hot and cold electron components,^35,36 the model for plasmas consisting of multi-ion species with different charge-to-mass ratios,^37–39 and the model assuming a truncated Maxwellian velocity distribution for electrons.^40,41 Our present study provides a new basis for extending such theories to the regime of multi-ps interactions of thin foils with high intensity lasers. Multi-dimensional effects in the multi-ps regime, e.g., effects of self-generated magnetic fields on the heating and confinement of electrons, are also important issues, which will be studied separately. Recently, a 100 Tesla level magnetic field can be generated in laboratory.⁴² By applying such a strong magnetic field along the laser axis, electron lateral diffusion can be suppressed, so that the electron recirculation will become more 1D like. Hence, our 1D model could be applicable even in tightly focused laser experiments. In the current study, the time evolution of electron temperature is obtained on the basis of numerical results. It is important to develop a theoretical scaling law for the electron heating dominated by the electron recirculation in the multi-ps laser-plasma interactions, which we will study in a future paper.

This study was supported by a Grant-in-Aid for the Japan Society for the Promotion of Science (JSPS) Fellows (Grant No. JP26-6405) and the JSPS KAKENHI Grant No. JP15K21767. This work was also supported by the JST A-STEP Grant No. AS27721002c. We appreciate Professor T. Johzaki and Professor Y. Kishimoto for fruitful discussions.