Quasi-static equations for the plasma in the presence of an intense laser pulse

Assuming the plasma response and the laser pulse evolution to be quasi-static i.e. changes in plasma/ laser occur on a time scale much longer than the typical time scale for electrons to move across the plasma, a set of quasi-static equation can be derived, valid for arbitrary laser intensities. For details of the quasi-static equations check P. Sprangle, E. Esarey, A. Ting , Phys. Rev. Lett 64, 002011 (1990).

Here we examine only the plasma response i.e. we solve numerically the equation for the plasma potential ϕ, the the rpesence of the laser vector potential a:
(((∂^2)ψ)/(∂(ζ^2)))=-(1/2)(1-((1+(a^2))/(((1+ϕ))^2)))

There are only (implicit) analytical solutions for a flat top profile: profile very far from what is used in simulations


Normalized units assumed throughout the notebook: 
space normalized to (c/(Subscript[ω, pe])), time normalized to (1/(Subscript[ω, pe])), charge to e, amd mass to (Subscript[m, e]). (same as numerical simulation)

Speed of light coordinates [check P. Sprangle, E. Esarey, A. Ting , Phys. Rev. Lett 64, 002011 (1990)]: 
ζ=z - ct and τ=t (attention: ζ has meaning of distance to the front of the laser pulse, thus behind the laser pulse leading edge ζ <0)

In this notebook, x plays the role of  ζ

Clear[L,a0,u,a2]; (* reset variables *)

(*------ Laser features -----*)

a0 = 4; (* normalized laser vector potential *)
xrise=  2 Pi /2;   (* laser rise time *) 
xfall=xrise ; (* laser fall time *)
k0 = 10; (* ratio of teh alser wavenumber to the plasma wavenumber *)

(* aux variables *)
x0=-(xrise+xfall);  (* position of trailing edge of laser pulse *)
L = (xrise+xfall); (* pulse length *)


(*------ Laser pulse longitudinal profile -------*)

(* wave laser profile - same as used in osiris *)
f[x_]=a0 If[x<xrise,If[x>0,10 (x/xrise)^3-15  (x/xrise)^4+6  (x/xrise)^5,0],If[x<xrise+xfall,10 ((xrise+xfall-x)/xfall)^3-15 ((xrise+xfall-x)/xfall)^4+6 ((xrise+xfall-x)/xfall)^5,0]]  Sin[k0 x]; 

 (* Gaussian Laser profile *)
(* f[x_]=a0 Exp[-2 (x-x0/4)^2/(xrise^2)] Sin[k0 x]; *) 

 (* Half sine Laser profile *)
(* f[x_]=a0 If[x<xrise+xfall,If[x>0,Sin[Pi x/(xrise+xfall)],0],0] Sin[k0 x];
*)

(* Flat top Laser profile *)
(* f[x_]=a0 If[x<xrise+xfall,If[x>0,1,0],0] Sin[k0 x]; *)

(* Plot laser pulse profile *)

xmax=-2 L;
Plot[{f[x-x0]},{x,xmax,0},PlotLabel->"Laser Pulse Profile",PlotRange->All,PlotPoints->100,Frame->True,PlotStyle->{{Hue[0.4],Thickness[0.015]}}, Background->GrayLevel[1],TextStyle->{FontFamily->"Helvetica",FontSize->14}, FrameTicks->{{{0,0},{-Pi,-Pi},{-2Pi ,-2Pi},{-3Pi ,-3 Pi},{-4Pi,-4Pi}},Automatic,Automatic,Automatic},AspectRatio->1];

(* Solution of the quasi-static equation for the potencial u[x] *)

(* size of the box where we are solving the equation, in units of L - laser pulse length *)
nbox =4; 
xmax=- nbox L;

solution=NDSolve[{(u')'[x]==-0.5(1-(1+f[x-x0]^2)/(1+u[x])^2),u'[0]==0,u[0]==0},u,{x,xmax,0},MaxSteps->10000];

Plot[{Evaluate[u[x]/.solution]},{x,xmax,0},PlotRange->All,PlotPoints->100,Frame->True,PlotLabel->"Plasma Potential",Background->GrayLevel[1],PlotStyle->{{Hue[0.8],Thickness[0.010]}},TextStyle->{FontFamily->"Helvetica",FontSize->14},FrameTicks->{{{0,0},{-Pi,-Pi},{-2Pi ,-2Pi},{-3Pi ,-3 Pi},{-4Pi,-4Pi}},Automatic,Automatic,Automatic},AspectRatio->1];


Plot[{Evaluate[-u'[x]/.solution]},{x,xmax,0},PlotRange->All,PlotPoints->100,Frame->True,PlotLabel->"Longitudinal Electric Field",Background->GrayLevel[1],PlotStyle->{{Hue[0.1],Thickness[0.010]}},TextStyle->{FontFamily->"Helvetica",FontSize->14},FrameTicks->{{{0,0},{-Pi,-Pi},{-2Pi ,-2Pi},{-3Pi ,-3 Pi},{-4Pi,-4Pi}},Automatic,Automatic,Automatic},AspectRatio->1];

Plot[{Evaluate[-(1-(1+f[x-x0]^2)/(1+u[x])^2)/.solution]},{x,xmax,0},PlotRange->All,PlotPoints->100,PlotLabel->"Electron Density",Background->GrayLevel[1],PlotStyle->{{Hue[0.95],Thickness[0.010]}},TextStyle->{FontFamily->"Helvetica",FontSize->14},FrameTicks->{{{0,0},{-Pi,-Pi},{-2Pi ,-2Pi},{-3Pi ,-3 Pi},{-4Pi,-4Pi}},Automatic,Automatic,Automatic},AspectRatio->1,Frame->True];