We study the one-dimensional expansion of a thin foil plasma irradiated by a high intensity laser with multi-picosecond (ps) pulse durations by using particle-in-cell simulation. Electrons are found to recirculate around the expanding plasma for many times, which results in stochastic heating leading to increase of the electron temperature in the multi-ps time scale beyond the ponderomotive scaling. The conventional isothermal model cannot describe such an expansion of plasmas in the long time scale. We here developed a non-isothermal plasma expansion theory that takes the time dependence of electron temperature into account for describing the multi-ps interactions 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. The time evolution of ion maximum energy obtained by the non-isothermal theory explains the details of that observed in the simulation.

## I. INTRODUCTION

Laser-matter interactions in the relativistic regime, where the laser intensity exceeds 10^{18} W/cm^{2}, 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,^{1} NIF-ARC,^{2} LMJ-PETAL,^{3} and OMEGA-EP^{4} 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,^{5} the wake field acceleration due to forward Raman scattering and/or self-phase modulation, the betatron acceleration in laser-plasma channels,^{6} 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^{9} and also by experiment^{10} that the electron slope temperature becomes much higher than the conventional ponderomotive scaling^{11} 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

*ω*is the laser angular frequency.

_{L}^{12}The other is that the self-generated magnetic fields on the front surface play an important role in electron acceleration.

^{13}

In interactions of thin solid foils with a relativistic laser in the regime of $I\lambda L2\u22651018\u2009W\u2009\mu 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*is the electron temperature, and

_{e}*M*is the ion mass.

_{i}^{18–22}Hence, the key parameter of the TNSA is the electron temperature

*T*resulting from the laser-plasma interaction. Conventionally in the sub-ps regime, the ponderomotive temperature, $Tp\u2261mec2(1+a02/2)1/2$, has been used to explain the electron temperature observed in experiments.

_{e}^{11}Here,

*m*is the electron rest mass,

_{e}*c*is the speed of light, and $a0=eE0/mec\omega L$ the normalized laser field amplitude, where

*e*is the fundamental charge and

*E*

_{0}is the amplitude of laser electric field. Note that in the conventional self-similar model of TNSA,

*T*is usually assumed to be determined by the peak value of the laser pulse amplitude $a\u03020$, 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.

_{p}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.^{23} It is also known that when the laser pulse duration *t _{p}* and foil width

*L*satisfy the relation 2

*L*<

*ct*, 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.

_{p}^{24}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.

^{25}

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.^{26} As the pulse duration was extended from 1.5 ps to 6 ps while keeping the peak laser intensity constant around $2.5\xd71018\u2009W/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.^{26} 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.

## II. PIC SIMULATION FOR MULTI-PICOSECOND INTERACTIONS OF LASER FIELDS AND SOLID FOILS

### A. Simulation condition

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.^{27} We consider linearly polarized laser fields assuming two pulse functions, i.e., a Gaussian pulse function $f=exp\u2009(\u2212(t\u2212t1)2/tp2)$ with *t _{p}* = 0.9 ps and

*t*

_{1}= 1.2 ps, and a longer pulse that rises and falls with Gaussian functions $f=exp\u2009(\u2212(t\u2212t1)2/tp2)$ and $f=exp\u2009(\u2212(t\u2212t2)2/tp2)$, respectively, while having a flattop region

*f*= 1 for

*t*

_{1}≤

*t*<

*t*

_{2}where

*t*

_{2}= 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

*λ*= 1.05

_{L}*μ*m and the peak normalized amplitude is $a\u03020=1.42$ which corresponds to the intensity of $I=2.5\xd71018\u2009W\u2009cm\u22122$. These parameters correspond to the LFEX experiment.

^{26}

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\xd71022\u2009cm\u22123$ and

*n*=

_{eb}*Zn*, respectively, from the position

_{ib}*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*is

_{ib}*n*= 40

_{eb}*n*where the charge state is

_{c}*Z*= 1 and

*n*= 1.0 × 10

_{c}^{21}cm

^{−3}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*

_{i}_{0}= 0.2 keV and

*T*

_{e}_{0}= 1 keV, respectively. Collisions and ionization processes are not included in the simulations performed in this study.

### B. Time evolution of electron temperature

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\epsilon e/Ne$, where

*ε*is the electron energy. Note that in obtaining the energy distributions, we observed electrons in the rear side region, i.e., from

_{e}*x*= 81

*μ*m to the rear boundary

*x*=

*L*, 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 $\u2212e\varphi f$ at each time is shown by vertical black dashed lines. Here, we set Φ = 0 at the rear surface of the target,

_{x}*x*= 80 μm. The electrons that come out from the target rear with energy below $\u2212e\varphi f$ are reflected by the sheath field of the foil, while energetic electrons whose energy at the target rear exceeds $\u2212e\varphi 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*, by fitting the energy spectrum below $\u2212e\varphi f$. We assume the 1D relativistic Maxwellian as the fitting function and obtain

_{e}*T*= 0.3 MeV, 0.8 MeV, and 1.5 MeV for

_{e}*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+a\u030202/2)1/2=0.2\u2009MeV$. The time evolution of *T _{e}* is fitted by the following exponential function:

Here, $T\u0303e=Te/T0$ is the normalized temperature, *T*_{0} is the initial temperature at time $t=t0,\u2009tmax$ is the time at which *T _{e}* takes the maximum value $Tmax$, and $\alpha \u2261Tmax/T0$. From the initial condition $T\u0303e(t0)=1$, the time scale

*t*is obtained as $tT=(tmax\u2212t0)/(ln\alpha )1/2$. The initial temperature

_{T}*T*

_{0}is assumed to be the ponderomotive temperature

*T*. We fit the parameters

_{p}*α*and $tmax$ in Eq. (1) by the simulation results. The time

*t*

_{0}is set to be

*t*

_{0}= 1 ps, at which

*T*∼

_{e}*T*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

_{p}*α*= 2.3 and 7.2 for 1 and 2 pulse cases, respectively, which represent that $Tmax=0.46\u2009MeV$ in the 1 pulse case and $Tmax=1.43\u2009MeV$ 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*= 0.45 MeV and 1.10 MeV for 1 and 2 pulse cases, respectively.

_{h}### C. Electron heating mechanisms in an expanding plasma

In Figs. 3(a) and 3(b), we present the snapshots of plasma densities, electrostatic potential $\varphi $, 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*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

_{x}*E*shown in (b) represents the spatially averaged value obtained from the neighbouring 30 meshes which corresponds to the length of 0.3

_{x}*μ*m. The value of

*E*is normalized by the pulse peak amplitude of the laser electric field $E\u03020$. The front peak value of $Ex/E\u03020$ at the position of the black arrow is twice larger than the spatial average of $Ex/E\u03020$ in Fig. 3(b) in the rear side where the electrostatic potential $\varphi $ starts to drop (from

_{x}*x*= 100

*μ*m to 130

*μ*m). This relation of

*E*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}*x*∼73

*μ*m, a sharp electron density jump is formed by the laser radiation pressure to create a charge separation and the strong

*E*field. In Fig. 3(a), the electrostatic potential $\varphi $ is shown by the black line, which we obtained by integrating the longitudinal electric field

_{x}*E*along

_{x}*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 $\varphi f\u223c\u22127\u2009MV$. Therefore, the electrons with energy below 7 MeV are confined by the potential $\varphi 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

*ω*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

_{L}*p*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}*x*= 80

*μ*m, i.e., the starting position of plasma expansion,

*n*

_{e}_{0}, is around 3.5

*n*which is higher than the relativistic cutoff density

_{c}*γn*∼ 1.4

_{c}*n*. We also confirmed that

_{c}*n*

_{e}_{0}increases from the subcritical density in time up to

*t*= 2.4 ps and saturates at

*n*

_{e}_{0}∼ 3.5

*n*. 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

_{c}*μ*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*, i.e., $(\gamma x\u22121)mec2$, where $\gamma x\u2261(1+px2/(mec)2)1/2$. We see that the electron recirculates between the rigid front surface

_{x}*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.^{29} 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,^{8} 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 $\varphi $ 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.

## III. PLASMA EXPANSION THEORY WITH TIME DEPENDENT ELECTRON TEMPERATURE

### A. The basic equations

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 _{i}*

_{0}and

*n*

_{e}_{0}=

*Zn*

_{i}_{0}, 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, $\varphi $ is the electrostatic potential, and we neglected the ion pressure term in the equations of motion assuming cold ions

*T*∼ 0. We assume the quasineutral condition

_{i}*n*=

_{e}*Zn*instead of solving the Poisson equation Eq. (5). Note that the quasineutrality breaks up at the expansion front

_{i}*x*=

*x*where the Debye length $\lambda D=(Te/(4\pi nee2))1/2$ exceeds the scale length of plasma expansion, i.e., time integral of the sound velocity

_{f}*C*

_{s.}^{21}In the isothermal case, the condition for breaking of the quasineutrality can also be represented by the ion plasma frequency

*ω*as

_{pi}*ω*< 1, where $\omega pi=(4\pi nie2/Mi)1/2$.

_{pi}tIn 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*is assumed to be constant. Equations (2)–(4) are the basic equations describing the 1D plasma expansion with time-dependent electron temperature.

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

Here, although the coefficient *n _{e}*

_{0}can be a function of time, we assume

*n*

_{e}_{0}= const. for simplicity. As we discussed in Sec. II C,

*n*

_{e}_{0}increases in time due to the recirculation. However, the change of

*n*

_{e}_{0}in time is smaller than that of the temperature

*T*. For the region where the quasineutrality is satisfied, i.e.,

_{e}*x*<

*x*, we obtain the ion density from Eq. (6) as

_{f}### B. Variable transformation

In the case of isothermal expansion, the solution that satisfies Eqs. (2)–(4) and the quasineutral condition *n _{e}* =

*Zn*is known to be described by a single variable, which is defined by

_{i}where the sound velocity *C _{s}* is constant. In this case, the solution for

*n*,

_{i}*v*, and $\varphi $ 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

_{i}*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*.

_{s}tWe here consider a transformation of variables $(x,t)\u2192(\xi ,t)$. Derivative operators are transformed as $\u2202x|t=R\u22121\u2202\xi |t$ and $\u2202t|x=\u2202t|\xi \u2212\xi CsR\u22121\u2202\xi |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*, and $\varphi $ are functions of

_{i}*t*and

*ξ*to be solved from the basic equations [Eqs. (2)–(4)], while

*T*and then

_{e}*C*are given as functions of

_{s}*t*.

### C. The transformed equations

### D. Ordering of each term of equations

We solve Eqs. (11) and (12) on the basis of perturbation expansion assuming that the time scale of temperature evolution $(d(log\u2009Te)/dt)\u22121=(2C\u0307s/Cs)\u22121$ 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\u2202t$, we expand the variables as

where $vi\u2032(1)$ and $\varphi \u2032(1)$ are functions of *t* and *ξ*, and have no explicit dependence on $R=R(t)$. Hereafter, we use the notations $vi(1)\u2261Rvi\u2032(1)$ and $\varphi (1)\u2261R\varphi \u2032(1)$. The ion density is expanded as $ni=ni(0)+\epsilon ni(1)+\cdots $.

and the first order terms are given by

Note that in the zeroth order of *ε*, *T _{e}*, and

*C*are constant as is in the isothermal case, and thus $R=Cst$. Therefore, we assumed $\u2202tvi(0)=0$ and $\u2202t\varphi (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

_{s}*ε*, $vi(0)$ and $\varphi (0)$ are treated as functions of

*t*and

*ξ*with time dependent

*C*and

_{s}*T*.

_{e}### E. The zeroth order self-similar solution

The corresponding electrostatic field *E*^{(0)} is derived as

which is proportional to $R\u22121\u221dt\u22121$ and uniform in *x*. The ion density is then derived from Eq. (7) as

### F. The first order solution

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 $\u2202\xi vi(0)=R\u0307,\u2009\u2202tvi(0)=R\xa8(1+\xi ),\u2009e\u2202\xi \varphi (0)=\u2212Te$, and $e\u2202t\varphi (0)=\u2212T\u0307e(1+\xi )$. Equations (18) and (19) then lead to

where the relation $b=RR\xa8/R\u03072$ 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\u223cO(\epsilon )$. 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., $T\u0307e>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\u0303=R/(Cs0t\u0302)$ where $Cs0=(ZT0/Mi)1/2$ is the sound velocity for the initial temperature *T*_{0}, $\tau =t/t\u0302$ is the time normalized by $t\u0302\u2261(2eN)1/2\omega pi0\u22121,\u2009\omega pi0$ is the ion plasma frequency for the density *n _{i}*

_{0}, i.e., $\omega pi0=(4\pi ni0e2/Mi)1/2$, and

*t*is the time scale of temperature evolution given in Eq. (1). The integrant in Eq. (31), $C\u0303s=Cs/Cs0$, is the normalized sound velocity. In Fig. 5(b), we plot

_{T}*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*

_{0}= 2.3 (red dashed-dotted lines) and

*α*= 7.2 (blue solid lines), where

*T*

_{0}is defined to be the ponderomotive temperature

*T*= 0.2 MeV. Note that the above values of

_{p}*α*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*

_{0}= 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*is merely equivalent to

_{s}*C*

_{s}_{0}

*t*. 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 $C\u0307s$ 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.

### G. Ion maximum energy

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 *ω _{pi}*

_{0}

*t*≫ 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.

^{21}We take the same assumption for the electric field at the ion expansion front

*E*, i.e.,

_{f}*E*is twice larger than Eq. (29), as

_{f}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}* = 2

*E*

^{(0)}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\u2009(1)$, and *λ _{D}*

_{00}denotes the Debye length for density

*n*=

_{e}*n*

_{e}_{0}and initial temperature

*T*=

_{e}*T*

_{0}. 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\u223cE(t=0)$ [Eq. (34)] for *R*/*λ _{D}*

_{0}≪ 1 which corresponds to

*ω*

_{pi}_{0}

*t*≪ 1 and

*E*∼

_{fa}*E*[Eq. (33)] for

_{f}*R*/

*λ*

_{D}_{0}≫ 1 which corresponds to

*ω*

_{pi}_{0}

*t*≫ 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\u0303e=1,\u2009R\u0303=\tau $, *b* = 0, and Θ = 1, the integral in Eq. (36) is solved analytically as $vf(0)=2Cs0ln(\tau +1+\tau 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 $\epsilon 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, $\epsilon imax$ depends on the plasma density at the origin of the expansion, i.e., *n _{e}*

_{0}defined in Eq. (6), through $t\u0302\u221d\omega pi0\u22121\u221d(ne0)\u22121/2$. We show here two results assuming

*n*

_{e}_{0}= 0.2

*n*and 4

_{c}*n*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

_{c}*n*

_{e}_{0}increases in time from subcritical density ∼0.2

*n*to 3.5

_{c}*n*as discussed in Sec. II C. The two lines assuming 0.2

_{c}*n*and 4

_{c}*n*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

_{c}*n*

_{e}_{0}varying from 0.2

*n*to 4

_{c}*n*. 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(\tau +1+\tau 2)$ by black lines. Here, the constant ponderomotive temperature

_{c}*T*= 0.2 MeV is assumed for the black dashed-dotted lines, i.e.,

_{p}*T*=

_{e}*T*= const., and the averaged temperatures ⟨

_{p}*T*⟩ = 0.25 MeV for 1 pulse case and ⟨

_{e}*T*⟩ = 0.68 MeV for 2 pulse case are assumed for the black thin solid lines. Here, ⟨

_{e}*T*⟩ are the averaged values of the temperature for times from

_{e}*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*in the integrant in Eq. (37).

_{e}For both 1 pulse and 2 pulse cases, most of $\epsilon 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 _{e}*

_{0}= 0.2

*n*initially, and it shifts to the first-order theory of

_{c}*n*

_{e}_{0}= 4

*n*later. This reflects the time evolution of

_{c}*n*

_{e}_{0}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 $\epsilon 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*and

_{p}*T*= ⟨

_{e}*T*⟩ 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

_{e}*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, $\epsilon 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 $\epsilon 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).

### H. Discussions and energy conversion to ions

We also investigated the dependence of $Tmax/Tp=\alpha $ on the laser peak amplitude $a\u03020$ 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\xd71018\u2009W/cm2\u2009(a\u03020=1.0)$ to $1.0\xd71019\u2009W/cm2\u2009(a\u03020=2.84)$ for the 2 pulse case. Therefore, $Tmax$ has the same dependence on $a\u03020$ 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 $a\u03020$. Here, under the approximation $(1+R\u03032)1/2\u223cR\u0303$ and Θ ∼1, the integrant on the RHS of Eq. (37) becomes $T\u0303e/R\u0303\u221d\alpha 1/2$, which is independent of $a\u03020$. Therefore, the dependence of the ion maximum energy $\epsilon imax$ on $a\u03020$ is the same as that of the ponderomotive temperature when we assume

*T*

_{0}=

*T*in Eq. (37). In the intensity regime from $a\u03020=1.0$ to 2.84, we obtained a scaling $Tmax\u22436.1\u2009Tp$ for the 2 pulse case. By using this scaling, we can estimate that $\epsilon imax$ reaches to the level of 100 MeV for the laser irradiation with the peak amplitude of $a\u03020\u22653.2$ and the pulse duration of

_{p}*t*≥ 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.

_{p}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^{26} where the energy conversion efficiency to ions increases when the laser pulse duration is extended from 1.5 ps to 6 ps.

## IV. CONCLUSION

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.^{42} 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.

## ACKNOWLEDGMENTS

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.