Anomalous collisional absorption of laser light in plasma using particle-in-cell simulations

TL;DR

This study uses particle-in-cell simulations with Monte Carlo collisions to investigate anomalous collisional laser absorption in under-dense plasma, revealing a non-monotonic dependence on laser intensity at low electron temperatures.

Contribution

First to demonstrate anomalous collisional laser absorption using PIC simulations with Monte Carlo modeling, connecting experimental observations with simulation results.

Findings

Absorption increases with laser intensity up to a critical point at low temperatures.

Absorption then decreases beyond the critical intensity, showing non-monotonic behavior.

The results bridge the gap between experimental findings and theoretical models.

Abstract

Collisional absorption of laser light in a homogeneous, under-dense plasma is studied by a new particle-in-cell (PIC) simulation code considering one-dimensional slab-plasma geometry. Coulomb collisions between charge particles in plasma are modeled by a Monte Carlo scheme. %[J. Comput. Phys. {\bf 25}, 205 (1977)]. %Both PIC and MC parts are individually benchmarked. For a given target thickness of a few times the wavelength of 800~nm laser of intensity $\I0$, fractional absorption ($\alpha$) of light due to Coulomb collisions (mainly between electrons and ions) is calculated at different electron temperature $\Te$ by introducing a total velocity $v = \sqrt{\vth^2 + \v0^2}$ dependent Coulomb logarithm $\ln\Lambda(v)$, where $\vth$, and $\v0$ are thermal and ponderomotive velocity of an electron. It is found that, in the low temperature regime ($\Te\lesssim15$~eV), fractional absorption of light anomalously increases with increasing $I_0$ up to a maximum corresponding to an intensity $I_c$, and then it drops when $I_0>I_c$. %(approximately) obeying the conventional scaling, i.e., %$\alpha \propto I_0^{-3/2}$ when $I_0>I_c$. Such an anomalous variation of $\alpha$ with $I_0$ in the low intensity regime was demonstrated earlier in experiments, and recently explained by classical and quantum models [Phys. Plasmas {\bf 21}, 013302 (2014); Phys. Rev. E {\bf 91}, 043102 (2015)]. % using the total velocity dependent cut-offs. %Here, we report anomalous nature of laser absorption by Here, for the first time, we report anomalous collisional laser absorption by %PIC simulations assisted by Monte Carlo collisions, PIC simulations, thus bridging the gap between models, simulations, and experimental findings.

Figures

Equations

$$\frac{\partial \bar{\mathbf{B}}}{\partial t}$$

$$\frac{\partial \mathbf{E}}{\partial t}$$

$$\dot{\mathbf{p}}$$

$$\frac{\partial B_x}{\partial t}$$

$$\frac{\partial E_z}{\partial t}$$

$$\frac{\partial E_y}{\partial t}$$

$$\frac{\text{mbox}\mathrm{d}{\psi }_{\pm }}{\text{mbox}\mathrm{d}t}=-4\pi J_z(t,y).$$

$$\frac{B_x^{n+1}(k+1/2)-B_x^n(k+1/2)}{\Delta t}$$

$$\frac{E_z^{n+1/2}(k+1)-E_z^{n-1/2}(k+1)}{\Delta t}$$

$$\frac{E_y^{n+1/2}(k+1)-E_y^{n-1/2}(k+1)}{\Delta t}$$

$$⟨{\delta }^2⟩=\frac{e_{\alpha }^2{e}_{\beta }^2{n}_L\ln \Lambda }{8\pi {ϵ}_0^2{m}_{\alpha \beta }^2{u}^3}\Delta t_c.$$

$$\ln \Lambda =\ln [{\lambda }_D(T_{\alpha }+T_{\beta })/2e_{\alpha }{e}_{\beta }]$$

$$\Theta =2\arctan \sqrt{-2⟨{\delta }^2⟩\ln (1-R_1)},$$

$$\Delta u_x$$

$$\Delta u_y$$

$$\Delta u_z$$

$$\Delta u_x$$

$$\Delta u_y$$

$$\Delta u_z$$

$${\mathbf{v}}_{\alpha }(t+\Delta t_c)$$

$${\mathbf{v}}_{\beta }(t+\Delta t_c)$$

$$⟨\frac{\text{mbox}\mathrm{d}{v}_{\alpha }}{\text{mbox}\mathrm{d}t}⟩$$

$$⟨\frac{\text{mbox}\mathrm{d}{\mathcal{E}}_{\alpha }}{\text{mbox}\mathrm{d}t}⟩$$

$$⟨\frac{\text{mbox}\mathrm{d}{v}_{\perp }^2}{\text{mbox}\mathrm{d}t}⟩$$

$${\nu }_{\mathrm{s}}$$

$${\nu }_{\mathrm{t}}$$

$${\nu }_{\mathrm{d}}$$

$$E_{z}(t,y_{l})=E_{0}\left\{\begin{array}[]{ll}\sin^{2}(\omega t/2n_{c})\cos(\omega t);&0<t<n_{c}T\\ 0&t>n_{c}T,\\ \end{array}\right.$$

$$\ln \Lambda =0.5\ln \left[1+{\left(b_{\max }/b_{\perp }\right)}^2\right],$$

$${\nu }_{\mathrm{ei}}(t)=\frac{{\omega }_p^2\ln \Lambda }{v_{\mathrm{os}}^3(t)}\left[\mathrm{erf}(u(t))-\frac{2}{\sqrt{\pi }}u(t)e^{-u(t{)}^2}\right].$$

$$\alpha =1-\exp (-2{\kappa }_i{L}_p)$$

Full text

I introduction

One of the main objectives of the researchers working in the field of laser-plasma interaction (LPI) is to couple more laser energy with the plasma (or matter) so as to obtain more energetic end-products, e.g., energetic charge particles or intense radiations. Therefore, it is of prime importance to know the underlying physical process (collisional and collisionless) by which laser energy is coupled to the plasma during the interaction. Earlier experiments bach ; teubner ; price ; cerchez and theoretical studies rozmus ; dong ; kruerBook ; mulserBook ; shalomBook ; gibbonBook have already reported various absorption processes, e.g., linear resonance manes77 , anharmonic resonance mulserx0 ; mulserx1 ; mulserx2 ; kundu2 ; kundu3 ; kundu4 ; kot03 ; cerchez , Brunel heating brunel ; mulserx3 , skin layer absorption bauer07 , $\bm{J}\times\bm{B}$ heating hong06 etc., which often depend on parameters of the laser, and the plasma gibbonBook ; kruerBook ; shalomBook ; mulserBook . For example, while passing through under-dense plasma (where plasma frequency $\omega_{\mathrm{p}}$ is less than the laser frequency $\omega$) an intense $p$-polarized short laser pulse can be absorbed by exciting wake-fields and instabilities gibbonBook ; kruerBook ; shalomBook ; mulserBook . On the other hand, in an overdense plasma with an under-dense pedestal, linear resonance absorption (LR) of $p$-polarized light may occur by meeting the resonance condition $\omega_{\mathrm{p}}=\omega$ in a specific location of the density gradient. Most often $p$-polarized light is used by experimentalists because of its ability to drive the plasma particles more efficiently, and relatively less attention is paid in LPI using $s$-polarized light. However, absorption of both $s$- and $p$- polarized light in plasma may happen through the electron-ion collision bornath ; bor01 ; kruerBook ; shalomBook known as inverse bremsstrahlung (IB) if laser intensity is below 10^{17}\,\mbox{\rm Wcm{}^{-2}}.

In this work, we concentrate on the absorption of a $s$-polarized laser light in a homogeneous, under-dense plasma-slab due to IB, since collisional and collisionless absorption processes are coupled together for a $p$-polarized light where it is difficult to know what fraction of the collisional absorption contributes to the total absorption. We are also motivated by some earlier experimental results riley93 ; shalomBook of collisional absorption with $s$-polarized light which shows that fractional absorption $\alpha$ of light (i.e., ratio of absorbed energy to the incident laser energy) anomalously increases initially with increasing laser intensity $I_{0}$ up to a maximum value about an intensity $I_{c}$, and then it drops nearly obeying the conventional scaling kruerBook ; shalomBook , i.e., $\alpha\propto I_{0}^{-3/2}$. Although there are numerous analytical models bor01 ; bornath ; sch97 ; hil05 ; men13 ; mol12 ; pert72a ; pert75 ; rand64 ; weng95 ; silin ; mulser1 ; mulser2 ; mulser3 ; kull01 ; wesson ; pert95 ; rae92 ; catto77 ; kremp01 ; brantov03 ; skupsky ; riley93 which directly or indirectly describe conventional (standard) collisional absorption (CA) without the effect of background plasma, less attempts were made to examine above mentioned anomalous collisional absorption (ACA). Recently, ACA process has been explained in the low temperature regime ($T_{\mathrm{e}}\lesssim 15$ eV) by postulating a total velocity dependent Coulomb logarithm $\ln\Lambda(v)$ (where $v=\sqrt{v_{\mathrm{th}}^{2}+v_{0}^{2}}$, $v_{\mathrm{th}}$ and $v_{0}$ are the thermal and ponderomotive velocity) in an analytical model of electron-ion collision frequency $\nu_{\mathrm{ei}}$ mulser1 ; kundu1 and more rigorous kinetic treatment kunduPRE . However, to what extent this analytical approach is valid can only be answered numerically with self-consistent dynamics of plasma background under the laser irradiation. To this end, we have developed an one-dimensional electromagnetic particle-in-cell code (henceforth we call it EMPIC1D) where variation of physical quantities (charge density, current density, electro-magnetic fields) depend only on the one spatial coordinate along the laser propagation direction while considering all three velocity components of charge particles. In a particle-in-cell (PIC) simulation a reduced number of computational particles is used to represent plasma, instead of a large number of actual physical particles birdsall ; hockney . This technique reduces the computational load, and enables to study the dynamics of an actual physical system of large number of charge particles. Sizes of these PIC particles (computational particles) are typically on the order of a numerical grid, they can pass through each other during the interaction, and Coulomb collisions do not naturally happen sentoku ; cadjan ; takizuka ; takizuka2 ; sma ; manheimer ; oliphant . For this reason, Coulomb collisions are explicitly added in all PIC codes. To include Coulomb collision in our EMPIC1D code a Monte Carlo (MC) technique proposed by Takizuka and Abe takizuka is adopted. Recently, this scheme is used in the PARASOL electrostatic PIC code takizuka2 to study kinetic effects in tokamak plasmas. It conserves total energy and total linear momentum before and after a collision event in the velocity space.

In this work, with the Monte Carlo collision assisted EMPIC1D code, for the first time we show the aforementioned anomalous absorption in an under-dense plasma in the low temperature regime similar to our earlier analytical works kundu1 ; kunduPRE . Thus we bridge the gap between the experimental findings, analytical models, and PIC simulations.

Plasma is assumed to be pre-ionized. Laser intensity is kept below 10^{18}\,\mbox{\rm Wcm{}^{-2}} so that relativistic effects are less important. For convenience, atomic units (a.u.) are used unless mentioned explicitly, i.e., $-e=m=4\pi\epsilon_{0}=\hbar=1$, where $|e|$ is the electronic charge and mass, $\epsilon_{0}$ is the permittivity of free space, and $\hbar$ is the reduced Planck constant.

This article is organized in the following manner. Details of the EMPIC1D code is given in Sec.II, with appropriate benchmarking in the absence of collision. In Sec.III we independently benchmark the Monte Carlo (MC) binary collision scheme. Then we combine the collision module with the EMPIC1D code, and study collisional absorption of a $s$-polarized light in an under-dense plasma slab in Sec.IV where a comparison is also made between simulation and theoretical results. A summary is given in Sec.V.

II Details of the PIC code

Here we give only necessary details of the PIC code. In PIC a collection of physical particles is represented by a computational particle so that the charge to mass ratio $q/m$ of the computational particle remains same as that of a physical particle. The following Maxwell-Lorentz system of equations (in the normalized form) is solved numerically after the discretization in space and time :

[TABLE]

Here, ${\bm{E}},{\bar{\bm{B}}}$ are the electric and magnetic part of the electromagnetic field, $\bm{p}$ is the particle momentum corresponding to its velocity $\bm{v}$ and position $\bm{r}$ at a time $t$. ${\bm{J}}$ is the current density vector, $c$ is the speed of light in the free space. The scaling ${\bar{\bm{B}}}=c{\bm{B}}$ connects actual magnetic field ${\bm{B}}$ with the scaled magnetic field ${\bar{\bm{B}}}$. The other equations, namely, ${\bm{\nabla}}\cdot{\bar{\bm{B}}}=0$ and the Gauss’s law ${\bm{\nabla}}\cdot{{\bm{E}}}=4\pi\rho$ are not explicitly solved in a standard multi-dimensional PIC scheme, thus saving a substantial amount of computer time. However, ${\bm{\nabla}}\cdot{\bar{\bm{B}}}=0$ is ensured by choosing a staggered grid, called Yee mesh. The charge and current conservation follows from $\partial{\rho}/\partial{t}+{\bm{\nabla}}\cdot{\bm{J}}=0$, thus ensuring ${\bm{\nabla}}\cdot{{\bm{E}}}=4\pi\rho$. Note that ${\bm{E}},{\bar{\bm{B}}},{\bm{J}},\rho$ are calculated on the grid points. Therefore ${\bm{E}},{\bar{\bm{B}}}$ field are interpolated to obtain corresponding fields ${\bm{E}_{p}},{\bar{\bm{B}}}_{p}$ at the particle (particle charge $q$ and mass $m$) position ${\bm{r}}$ using linear weighting scheme, and the Lorentz equation (3) is solved using the standard leap-frog method. The advantage of the scaling ${\bar{\bm{B}}}=c{\bm{B}}$ is that, it reduces equations (1) and (2) identical in form in the free space (i.e., when $\bm{J}=0$) and the amplitudes of $\bm{E},\bm{\bar{B}}$ becomes comparable. From now onward, for convenience, we shall write $\bm{B}$ in stead of $\bm{\bar{B}}$ unless mentioned explicitly. A typical cycle of the PIC simulation is shown in Fig.1 which clearly depicts where the binary collision module (the MC part) should be incorporated.

II.1 Simplification in one-dimension

Let us consider a s-polarized light (propagating in $y$-direction) with transverse field components $E_{z},B_{x}$. The physical quantities (e.g., charge density, current density, electro-magnetic fields) are assumed to depend only on the space co-ordinate $y$, while retaining all three velocity components ($v_{x},v_{y},v_{z}$) of particles. Components of Eqn.(1)-(2) reads

[TABLE]

Equation (6) gives longitudinal component of the electric field $E_{y}$ in our case. It is important to mention that, the numerical implementation of our PIC code is little different from some of the traditional 1D-PIC codes, namely, EM1BND birdsall , LPIC++ gibbonBook ; lichters , but closely follow the implementation in PIC codes PSC harmut , VLPL pukhov , VPIC bowersVPIC , OSIRIS foncecaOSIRIS , and VORPAL nieterVORPAL . In EM1BND birdsall and LPIC++ gibbonBook ; lichters , by performing addition and subtraction of Eqn.(4) and (5), and writing $\psi_{\pm}=E_{z}\pm B_{x}$ one finds $(\partial_{t}\pm c\partial_{y})\psi_{\pm}=-4\pi J_{z}(t,y)$; where $\psi_{\pm}$ can be recognized as the two propagating solutions of the wave equation. The advantage in this traditional procedure birdsall ; gibbonBook ; lichters is that the partial derivative $(\partial_{t}\pm c\partial_{y})$ can be written in terms of the total derivative in time w.r.t. an observer moving at a speed $\pm c$, leading to

[TABLE]

For a given $J_{z}$, Eq.(7) is solved as an ordinary differential equation (ODE) to obtain transverse fields $E_{z}=(\psi_{+}+\psi_{-})/2$, $B_{x}=(\psi_{+}-\psi_{-})/2$ on the grid. It also allows larger time step $\Delta t=\Delta y/c$. The disadvantage is that, the longitudinal component $E_{y}$ is obtained by solving the Poisson’s equation ${\partial E_{y}}/{\partial y}=4\pi\rho$ explicitly (not from Eq.(6)), and it needs a separate algorithm than solving the transverse components. Moreover, this traditional scheme is hard to extend in multi-dimensional case. In our EMPIC1D code we use FDTD (finite difference in time domain) method for the solution of all field components (both traverse and longitudinal), in stead of the traditional addition-subtraction method mentioned above. Thus we use only one kind of algorithm for $E_{z},B_{x},E_{y}$ which is extendable to PIC simulations in higher dimensions as in Refs.harmut ; bowersVPIC ; pukhov ; foncecaOSIRIS ; nieterVORPAL . Using FDTD procedure on the Yee mesh, and assuming $t=n\Delta t$, $y=k\Delta y$, Eqn.(4),(5),(6) can be written as

[TABLE]

To ensure numerical stability we take $c\Delta t/\Delta y=1/2$, which decides the time step $\Delta t$ for a chosen grid size $\Delta y$. The dispersion due to the FDTD discretization is minimized by choosing sufficient number of spatial grids (minimum 40 is taken) per wavelength of light. VLPL like pukhov improvement, namely, dispersion free scheme we plan to include in future. The current density $\bm{J}$ due to the motion of charge particles is computed using the “explicit current conserving scheme” by Umeda et al. umeda which satisfies $\partial{\rho}/\partial{t}+{\bm{\nabla}}\cdot{\bm{J}}=0$. Thus we avoid explicit solution of the Poisson’s equation to obtain $E_{y}$. We use “perfectly matched layer” (PML) absorbing boundary condition sullivan for the electromagnetic fields. For charge particles, however, depending upon physical situations, absorbing, periodic, and reflecting boundary conditions are used.

II.2 Benchmarking of the PIC code

The EMPIC1D code is verified for different cases (where analytical solutions exist) : (i) plasma oscillation and energy conservation (without external light source), (ii) reflection of the laser light interacting with an over-dense plasma and corresponding field in the skin layer, (iii) transmission and reflection coefficient of laser light in an under-dense plasma, and (iv) interaction of light with inhomogeneous plasma and formation of standing waves. As a representative case, for the purpose of benchmarking, here we consider only the last case.

II.2.1 Interaction of s-polarized light with inhomogeneous plasma

Consider the interaction of a laser light with an inhomogeneous plasma which consists of an under-dense region and an over-dense region separated by a critical density ($\rho_{c}=\omega^{2}/4\pi$) surface, illustrated in Fig.2. This problem was investigated analytically in Refs.kruerBook ; shalomBook ; Ginzburg for time independent case. However its numerical verification with PIC method is still rare in the literature. Benchmarking of this problem alone can justify the correctness of interaction of laser field with under-dense and over-dense plasma in a single simulation run using EMPIC1D. We assume that density variation of plasma is linear $\rho_{0}(y)=\rho_{\mathrm{c}}y/L$ along the propagation direction $y$ of a $s$-polarized light described by $E_{z},B_{x}$. $L$ is the distance of the critical density surface (vertical dashed line) from the vacuum plasma interface on the left side. A continuous laser light $E_{0}\sin(\omega t)$ of amplitude $E_{0}$ is sourced through the vacuum plasma interface. After propagating through the under-dense region $\rho_{0}(y)/\rho_{c}<1$, light is reflected from the critical surface. The reflected light and the incident light are superimposed to form a standing wave after a long time when steady state is reached, and the resultant field $E_{z}$ satisfies the wave equation kruerBook ; shalomBook $\frac{\mbox{\rm d}^{2}E_{z}}{\mbox{\rm d}\eta^{2}}-\eta E_{z}=0$ with $\eta(y)=(\omega^{2}/c^{2}L)^{1/3}(y-L)$. Analytical solution $E_{z}(\eta)=c_{1}A\mathrm{i}(\eta)+c_{2}B\mathrm{i}(\eta)$ is found in terms of Airy functions $A\mathrm{i},B\mathrm{i}$ with $c_{1},c_{2}$ as constants (independent of $\eta$). For the standing wave solution $c_{2}=0$ is chosen, since $B\mathrm{i}\rightarrow\infty$ as $\eta\rightarrow\infty$, leading to $E_{z}(\eta)=c_{1}A\mathrm{i}(\eta)$. In reality, the problem is time dependent. Desired steady state analytical solution may not be reached during the early time of interaction. The constant $c_{1}$ depends upon the amplitude and phase of the resultant wave at the vacuum plasma interface at a given instant of time. We adjust $c_{1}$ with the effective magnitude of the field $E_{s}^{t_{n}}$ at a time $t=t_{n}$ retrieved from the PIC simulation at the vacuum plasma interface (at $y_{l}$) to match the analytical solution with the PIC simulation such that $c_{1}(t_{n})=E_{s}^{t_{n}}(y_{l})/A\mathrm{i}(\eta(y_{l}))$. It leads to the analytical solution as $E_{z}^{t_{n}}(y)=\left[E_{s}^{t_{n}}(y_{l})/A\mathrm{i}(\eta(y_{l}))\right]A\mathrm{i}(\eta(y))$.

We take laser field parameters as $E_{0}=0.2$ a.u. , and $\omega=0.057$ a.u. corresponding to I_{0}\approx 7\times 10^{15}\,\mbox{\rm Wcm{}^{-2}} and $\lambda=800$ nm respectively. The length of the simulation box is $L_{b}=1000\Delta$ ($\Delta=\Delta y=\lambda/40$ is the grid size) with total number of grids $N_{g}=1000$. Out of these 1000 grids, plasma particles centrally occupy 990 grids while 5 grids are left on each boundaries to separate the plasma from the vacuum. The source is located at the left plasma vacuum interface at $y=y_{l}=5\Delta$. The right boundary of the plasma is at $y=y_{r}=995\Delta$. There are $N=7920$ PIC ions (assumed to be stationary) each of charge $q_{s}=0.00325$ a.u. and an equal number of overlapping PIC electrons with mass $m_{s}=q_{s}$, and charge $-q_{s}$. Particles are loaded according to the scheme given in Ref.birdsall to obtain a linear density profile as shown in Figs.2(a)-(d) with $\rho_{0}/\rho_{\mathrm{c}}$ rising linearly from 0 to 2 (dark solid line, blue). The field profiles $E_{z}/E_{0}$ from PIC simulation (dashed oscillatory, red) and the analytical solution (solid oscillatory, green) are plotted against $y/\Delta$ in Figs.2(a-d) at different times $t_{n}=5,10,15,20$ respectively. At an early time $t_{n}=5$, PIC solution for $E_{z}$ does not match with the analytical solution (in Fig.2a) since at this early time superposition of the reflected wave with the incident wave is still incomplete in the PIC simulation to form the desired standing wave. As time increases, difference between the PIC profile and the analytical profile gradually decreases which is evident from the comparison between Fig.2a and Fig.2b. After a sufficiently longer time $t_{n}=15$ and beyond, the reflected wave meets the incident wave with required amplitude and phase so that the PIC solution matches (in Figs.2c,d) the analytical solution (see oscillatory dashed line coincides the solid line). The match between analytical and the numerical solution clearly ensures the correctness of our PIC simulation without collision.

III Binary collision: model and simulation

In the PIC simulation, as discussed above, particles may pass through each other during their close encounter and collision effects are omitted. To implement binary collision in the EMPIC1D code we have followed the Monte Carlo scheme given by Takizuka and Abe takizuka and also the work by Ma et al. sma , and Sentoku et al. sentoku . For conciseness we only show the validation of our implementation with a minimum detail. The main approximation of binary collision is that at a given instant only two particles will collide, and the effect of collision arises due to the cumulative effect of many small angle binary collisions. Within a computational cell, particles are paired randomly (ion-ion, ion-electron, electron-electron) and then collision is performed between every pair. The maximum impact parameter in a fully ionized quasi-neutral plasma being of the order of the Debye length truvnikov , the maximum size of the collision grid is also restricted to the Debye length. Collision event takes place in the velocity space, meaning that the velocity components of the particles changes but the co-ordinates are not influenced at that time. The post collision velocities are obtained by going to the center of mass (COM) frame of the respective collision pairs and then back to the laboratory frame. Due to collision, during a small time interval $\Delta t_{c}$ (which is sufficiently small compared to the mean relaxation time), the direction of velocities of the colliding particles changes but not their magnitudes. For instance, we consider a system of two particles from two species $\alpha$ and $\beta$ having velocities $v_{\alpha}$ and $v_{\beta}$, masses $m_{\alpha}$ and $m_{\beta}$, densities $n_{\alpha}$ and $n_{\beta}$, and charges $e_{\alpha}$ and $e_{\beta}$. At a given instant $t$, the effect of collision leads to the rotation of the relative velocity ${\bm{u}}={\bm{v}}_{\alpha}-{\bm{v}}_{\beta}$ in the COM frame of the two particles. The relative velocity ${\bm{u}^{\prime}}={\bm{v}^{\prime}}_{\alpha}-{\bm{v}^{\prime}}_{\beta}$ after the collision (primes represents quantities after collisions) in the COM frame, and the rotation of the velocity vector ${\bm{u}}\rightarrow{\bm{u}^{\prime}}$ can be described by the scattering angle $\Theta$ and the azimuth angle $\Phi$ (see Fig.3) which are chosen randomly for a given pair ($\alpha,\beta$). In order to find $\Theta$, a parameter $\delta$ is introduced such that $\delta=\tan(\Theta/2)$.

The variable $\delta$ is chosen randomly from a Gaussian distribution such that its mean is zero, and corresponding variance $\left<\delta^{2}\right>$ is takizuka ; takizuka2 ; sma

[TABLE]

Here $n_{L}$ is the minimum density between $n_{\alpha}$ and $n_{\beta}$, $m_{\alpha\beta}=m_{\alpha}m_{\beta}/(m_{\alpha}+m_{\beta})$ is the reduced mass, and $\ln\Lambda$ is the Coulomb logarithm. $\ln\Lambda$ can be calculated as truvnikov

[TABLE]

where $\lambda_{D}$ is the Debye length, and $T_{\alpha},T_{\beta}$ are the temperatures of the respective species. Above expression of $\ln{\Lambda}$ is valid in the absence of external laser field. Otherwise, $\ln\Lambda$ should include the response of electrons to the laser field strength $E_{0}$ and the frequency $\omega$. We shall use Eq.(12) only for the validation of binary collision event when there is no external force. The necessary modification of $\ln{\Lambda}$ with laser field will be discussed later. Deflection angle $\Theta$ is calculated by using Box-Muller method with distribution $p(\delta)d\delta=(1/\langle\delta^{2}\rangle)\exp(-\delta^{2}/2\langle\delta^{2}\rangle)\delta d\delta$ as given in Refs.pert99 ; cohen2013 ,

[TABLE]

where $R_{1}$ is an uniform random number between 0 and 1. The azimuth angle $\Phi$ is chosen as $\Phi=2\pi R_{2}$, with $R_{2}$ as an uniform random number between 0 and 1. The change in velocity components in the laboratory frame can be calculated as takizuka

[TABLE]

where $u_{\perp}=\sqrt{u_{x}^{2}+u_{y}^{2}}$. When $u_{\perp}=0$, we take

[TABLE]

Final post collision velocities in the laboratory frame reads

[TABLE]

III.1 Validation of the binary collision model

In order to validate our implementation of binary collision, we consider particle collisions in a two component plasma. The conservation of energy, and momentum for each colliding pair of particles are rigorously checked. We benchmark different relaxation rates (collision frequencies): (i) mean rate of change of velocity of the electrons (called slowing down frequency $\nu_{\mathrm{s}}$), (ii) mean rate of change of energy of the test electrons (energy transfer frequency, $\nu_{\mathrm{t}}$), and (iii) the rate of spread of the test particle velocity transverse to its original direction (deflection frequency, $\nu_{\mathrm{d}}$). Analytically, these frequencies are obtained by the test particle theory truvnikov ; spitzer in which a test particle (electron, designated by $\alpha$) is assumed to move through a medium of field particles (ions, designated by $\beta$) having Maxwellian velocity distribution. These frequencies are related to the transport equations truvnikov ; spitzer

[TABLE]

and can be expressed in terms of the integral $\mu(x)=\frac{2}{\pi}\int_{0}^{x}\sqrt{t}\exp(-t)dt$ and its derivative $\mu^{\prime}(x)=\frac{2}{\pi}\sqrt{x}\exp(-x)$ as

[TABLE]

Here $\nu_{0}={e_{\alpha}^{2}e_{\beta}^{2}n_{\beta}\ln{\Lambda}}/\left({8\pi\epsilon_{0}^{2}\sqrt{2m_{\alpha}\mathcal{E}_{\alpha}^{3}}}\right)$ is the basic collision frequency, $\mathcal{E}_{\alpha}=m_{\alpha}v_{\alpha}^{2}/2$ is the kinetic energy of the projectile, and $x=m_{\beta}\mathcal{E}_{\alpha}/{m_{\alpha}T_{\beta}}$. To obtain $\nu_{\mathrm{s}},\nu_{\mathrm{t}},\nu_{\mathrm{d}}$ an average over a group of test particles (electrons) having the same velocity are taken for the better statistics.

We consider 10000 test particles (electrons) and an equal number of field particles (ions). Test particles are assumed to be composed of 40 groups, and each of 250 particles in a group has the same initial velocity. Their velocities are normalized by the thermal velocity, $v_{\mathrm{th}}=\sqrt{T_{\alpha}/m_{\alpha}}$. Here we assume $m_{\beta}/m_{\alpha}=10000$ (i.e. the field particles are highly massive), $T_{\alpha}=T_{\beta}=T=5$ eV, $e_{\alpha}=e_{\beta}=e$ (charge of the electron), $n_{\alpha}=n_{\beta}=n_{d}=10^{15}\,\mathrm{cm}^{-3}$, and $\nu_{0}\Delta t_{c}=10^{-3}$. The value of $\nu_{0}\Delta t_{c}$ is chosen small to ensure small angle collisions sma . In the Monte Carlo simulation, $\nu_{\mathrm{s}},\nu_{\mathrm{t}},\nu_{\mathrm{d}}$ are calculated using the Eqn.(22),(23), and (24) with the velocities from Eqn.(20), and (21) before and after collisions.

Figure 4 shows the normalized frequencies ($\nu_{\mathrm{s}}/\nu_{0},\nu_{\mathrm{t}}/\nu_{0}$, and $\nu_{\mathrm{d}}/\nu_{0}$) versus the energy of the test particles $\mathcal{E}_{\alpha}/T_{\alpha}$ after $n_{t}=20$ time steps. Each solid circle (numerical) represents the average over a group of test particles. An excellent agreement between the numerical results and the analytical values (solid lines) from Eqn.(25),(26), and (27) ensures the correctness of our implementation.

IV Absorption in an under-dense plasma

We have integrated the above collision module with the EMPIC1D code described in Sec.II to study collisional absorption of light incident normally on an under-dense plasma slab of uniform density. The simulation domain consists of $N_{g}=500$ computational cells with the plasma slab at the center. Initially each computational cell contains equal number of electrons and ions so that plasma is charge neutral. The temporal profile of the laser pulse (at the left boundary, $y_{l}$) is chosen as

[TABLE]

with $n_{c}$ as the number of cycles, and $T=2\pi/\omega$ as the laser period. The pulse is numerically excited at $y=y_{l}$, propagates in free space, then strikes the plasma slab. The intensity, wavelength, number of cycles, the duration of pulse, the width $L_{p}$ of the plasma slab can be varied as desired. Accordingly the length of the computational domain, and the number of computational cells $N_{g}$ are also adjusted. We choose the laser wavelength $\lambda=800$ nm with $n_{c}=4$-cycles, and the total pulse duration $\approx 30$ fs. The size of a computational cell is chosen as $\Delta=200$ a.u. which yields the PIC time step ${\Delta t}_{\mathrm{PIC}}=0.729248$ a.u. Length of the plasma is chosen as $L_{p}\approx 1.32\lambda$ with a plasma density $\rho/\rho_{c}\approx 0.136$. Temperature of ions are kept fixed in all simulations at $T_{\mathrm{i}}=5$ eV while temperature of electrons are kept fixed at different values, e.g., $T_{\mathrm{e}}=5,10,...,50,100$ eV for a given laser intensity. The chosen value of $T_{\mathrm{i}}$, however, is found to have negligible effect on the overall results of collisional absorption. Above parameters are kept fixed during a simulation run unless mentioned explicitly. To simulate inverse bremsstrahlung absorption in presence of a laser field the Coulomb logarithm should not be same as Eq.(12) for ordinary collisions, since it does not include the laser field parameters. Because there is no unique model of $\ln\Lambda$ in presence of laser field, we use a modified Coulomb logarithm kundu1 ; mulser1 ; pert95

[TABLE]

where $b_{max}=V_{t}/\max(\omega,\omega_{p})$, $b_{\perp}=e^{2}/(4\pi\epsilon_{0}m_{e}V_{t}^{2})$, with $V_{t}^{2}=v_{th,e}^{2}+v_{0}^{2}$ as the total velocity. The effect of laser field is incorporated through the ponderomotive velocity $v_{0}=eE_{0}/m_{e}\omega$. In the absence of laser field, $V_{t}=v_{th,e}$, $b_{max}=V_{t}/\omega_{p}=\lambda_{D}$ and $\ln{\Lambda}$ in Eq.(29) becomes nearly equal to that given in Eq.(12).

From the simulation, we record total kinetic energy $ke$ gained by the particles, the electric part $ee=\sum_{1}^{N_{g}}E_{j}^{2}\Delta/8\pi$ and the magnetic part $me=\sum_{1}^{N_{g}}B_{j}^{2}\Delta/8\pi$ of the electromagnetic energy at every time step, giving the total energy $te=ke+ee+me$. Figure 5 shows temporal variation of various energies at a given laser intensity I_{0}=5\times 10^{14}\,\mbox{\rm Wcm{}^{-2}} for the two cases: (a) without collision, and (b) with collision between electrons and ions. In Figs.5 (a)-(b) for the initial time upto $t/T\approx 3$, all energies $ee,me$, and $te$ increase sharply since the laser pulse is entering the simulation domain. For $t/T\approx 3-4$ the value of $ee,me$ remains almost constant and $te$ reaches a maximum because entire pulse has appeared in the simulation domain. The pulse strikes the plasma slab about $t/T\approx 4$, and only after this time, for $t/T=4-7.5$, $ke$ first increases and then drops with the corresponding drop and increase in $ee$, and $me$ while $te$ remain conserved at the highest value. After $t/T\approx 7.5$, values of $ee$, $me$ and $te$ sharply drop since the laser pulse is leaving the finite simulation domain, and gets absorbed (artificially) in the right boundary. The constant value of total energy, when the entire pulse is inside the computational box (for $t/T=3-7.5$), indicates conservation of energy in the simulation. In Fig.5(a), without collision, $ke$ reaches a maximum value, and finally drops to zero before $t/T\approx 7.5$. This is expected, because particles can not retain this energy, and finally give back to the electromagnetic fields (which is also evident from the corresponding drop and rise of $ee$ and $me$ between $t/T=4-7.5$), resulting no net absorption. However, when collision is taken into account [in Fig.5(b)] $ke$ increases monotonically in time starting at $t/T\approx 4$ (with corresponding drop in $ee$ and $me$, meaning absorption of the pulse), and reaches a non-zero saturation value around $t/T=6.8$ much before the pulse has left the simulation box. $ke$ does not drop to zero even after the pulse is over which clearly shows that s-polarized light can be absorbed due to collisions and the laser energy can be transferred to the charge particles.

We now find nature of collisional absorption by varying the intensity of the laser pulse for a given initial temperature $T_{\mathrm{e}}$. In reality, $T_{\mathrm{e}}$ should also vary during the interaction. But our pulse being very short we assume it to be unchanged. The other parameters, such as plasma thickness $L_{p}$, plasma density $\rho$, ion temperature $T_{\mathrm{i}}$ are kept constant as above.

Figure 6 shows fractional absorption $\alpha$, defined as the ratio of the final kinetic energy retained in the particles to the maximum of $te$ (which is actually the total energy in the laser pulse), versus the peak intensity for $T_{\mathrm{e}}=5,10,50,100$ eV (solid lines). It is seen that, for higher temperatures $T_{\mathrm{e}}>20$ eV, $\alpha$ initially remains almost constant (or vary slowly) upto a certain value I_{c}\approx 6\times 10^{13}\,\mbox{\rm Wcm{}^{-2}} of the peak intensity, and then decreases gradually for intensities $I_{0}>I_{c}$. This is the conventional result of collisional absorption reported in earlier works bornath ; bor01 ; decker with a $\ln\Lambda$ independent of the ponderomotive velocity $v_{0}=E_{0}/\omega$. However, at a lower $T_{\mathrm{e}}<20$ eV, it is found that $\alpha$ initially increases with the intensity, reaches a maximum value about an intensity $I_{c}$, then drops similar to the high temperature case. Such an anomalous behavior (initial increase followed by a drop) of fractional absorption versus the laser intensity was reported experimentally with normally incident $s$-polarized light (of wavelengths 800 nm shalomBook and 268 nm riley93 on an under-dense plasma with the peak absorption more than 30%. Incorporating a total velocity dependent $\ln\Lambda$, in our EMPIC1D code assisted by Monte Carlo binary collision we reproduce similar anomalous nature of collisional absorption in the low temperature regime. Our results indicate that fractional absorption due to collisional processes can be as high as 40% or even more for different plasma and laser parameters. For the shake of completeness, numerical results are compared with analytical estimates (dashed line in Fig.6) using a modified $\ln\Lambda$ as in Eq.(29) in the ballistic model mulser1 ; kundu1 of time-dependent electron-ion collision frequency

[TABLE]

Where $u(t)=v_{\mathrm{os}}(t)/\sqrt{2}v_{\mathrm{th}}$, and $v_{\mathrm{os}}(t)$ is oscillation velocity of the electron in the laser field. Averaging $\nu_{\mathrm{ei}}(t)$ over a laser period leads to average $\overline{\nu}_{\mathrm{ei}}$ and fractional absorption

[TABLE]

of a continuous light of frequency $\omega$ in an under-dense plasma slab kruerBook ; shalomBook at normal incidence. Here $\kappa_{i}=(\rho/\rho_{\mathrm{c}}){\overline{\nu}_{\mathrm{ei}}}/{v_{\mathrm{g}}}$, and $v_{\mathrm{g}}=c\sqrt{1-\rho/\rho_{\mathrm{c}}}$ is the group velocity of light. Analytical result using Eq.(31) (dashed line) at a lower temperature $T_{\mathrm{e}}=5$ eV shows qualitative agreement with the EMPIC1D result, and confirms the anomalous nature of collisional absorption which was also reported by quantum and classical kinetic models kundu1 ; kunduPRE . However, there are discrepancies between numerical and analytical results at higher temperatures, which may be due to (i) time varying field experienced by particles, (ii) movement of ion back-ground to conserve momenta and energy during binary collisions in the numerical simulations as opposed to the analytical model where all particles experience same peak laser field $E_{0}$, and ions are considered stationary.

V Summary

Collisional absorption of s-polarized laser light in a homogeneous, under-dense plasma is studied by a new particle-in-cell (PIC) simulation code considering one-dimensional slab-plasma geometry. To account for Coulomb collisions between charge particles a Monte Carlo (MC) binary collision scheme is used in the PIC code. For a given target thickness of a few times the wavelength of 800 nm laser fractional absorption of light due to Coulomb collisions is calculated at different electron temperature $T_{\mathrm{e}}$ using a total velocity $v=\sqrt{v_{\mathrm{th}}^{2}+v_{0}^{2}}$ dependent Coulomb logarithm $\ln\Lambda(v)$. In the low temperature regime ($T_{\mathrm{e}}\lesssim 15$ eV) it is found that fractional absorption ($\alpha$) of light anomalously increases initially with increasing intensity $I_{0}$ up to a maximum value corresponding to an intensity $I_{c}$, and then it drops approximately obeying the conventional scenario, i.e., $\alpha\propto I_{0}^{-3/2}$ when $I_{0}>I_{c}$. Anomalous increase of $\alpha$ with $I_{0}$ was demonstrated in some earlier experiments shalomBook ; riley93 , and recently explained by various models kundu1 ; kunduPRE using total velocity dependent cut-offs. Here, we report anomalous nature of laser absorption by self-consistent PIC simulations assisted by Monte-Carlo collisions, thus bridging the gap between the models, simulations, and experiments.

Acknowledgements.

The author would like to thank Anshuman Borthakur for the initial help in the Monte-Carlo simulations, the Plasma Science Society of India (PSSI) for providing partial financial support as a PSSI fellowship to carry out this work and Sudip Sengupta for valuable suggestions.