Reorganising of Mironov's for merge pull request into templates

doc/Sphinx/Understand/ionization.rst

@@ -72,6 +72,66 @@ over a period :math:`2\pi/\omega` leads to the well-known cycle-averaged ionizat
   \,I_p\,\left( \frac{2 (2 I_p)^{3/2}}{\vert E\vert} \right)^{2n^\star-\vert m \vert -3/2}\,
   \exp\!\left( -\frac{2 (2 I_p)^{3/2}}{3 \vert E\vert}  \right)\,.
 
+
+In :program:`Smilei`, four models are available to compute the ionization rate of :eq:`ionizationRate1`.
+
+**Tunnelling Ionization (PPT-ADK) for :math:`m=0`**
+
+In the classical model, the ionization rate of :eq:`ionizationRate1`
+is computed for :math:`\vert m \vert=0` only.
+Indeed, as shown in [Ammosov1986]_, the ratio :math:`R` of the ionization rate
+computed for :math:`\vert m\vert=0` by the rate computed for :math:`\vert m\vert=1` is:
+
+.. math::
+
+  R = \frac{\Gamma_{{\rm qs},\vert m \vert = 0}}{\Gamma_{{\rm qs},\vert m \vert = 1}}
+  =  2\frac{(2\,I_p)^{3/2}}{\vert E\vert}
+  \simeq 7.91\,10^{-3} \,\,\frac{(I_p[\rm eV])^{3/2}}{a_0\,\hbar\omega_0[\rm eV]}\,,
+
+where, in the practical units formulation, we have considered ionization
+by a laser with normalized vector potential :math:`a_0=e\vert E\vert /(m_e c \omega_0)`,
+and photon energy :math:`\hbar\omega_0` in eV.
+
+
+
+**Tunnelling Ionization (PPT-ADK) for :math:`|m|>0`**
+In this model,d ependence on the magnetic quantum number :math:`m` is added. 
+
+:math:`m` is attributed to eac electron in accordance with the following rulse:
+1. Since :math:`\Gamma_z(m=0)>\Gamma_z(m=1)>\Gamma_z(m=2)>...` we assume that for electrons
+with the same azimuthal quantum number :math:`l`, the states with the lowest value of
+:math:`|m|` are ionized first.
+2. Electrons with the same azimuthal qunatum number :math:`l` occupy the sub-shells in the
+order of increasing :math:`|m|` and for the same :math:`|m|` in the order of increasing :math:`m`. 
+
+With this algorithm, by knowing the atomic number A, we can assign a unique set of
+quantum numbers :math:`nlm` to each electron on the atomic sub-shells and identify their extraction
+order during successive ionization. 
+
+
+**Barrier Suppression Ionization (Tong&Ling)**
+The formula proposed by Tong and Lin [Tong2005]_ extends the tunnelling ionization rate to the barrier-suppression
+regime. This is achieved by introducing the empirical factor in :eq:`ionizationRate1`:
+
+.. math::
+
+  \Gamma_{Z^\star}^{TL} = \Gamma_{Z^\star} \times \exp (-2\alpha_{TL}n^{\star2}\frac{E}{(2I_p)^{3/2}}),
+
+where :math:`\alpha_TL` is an emprirical constant with value typically from 6 to 9. The actual value
+should be guessed from empirical data. When such data is
+not available, the formula can be used for qualitative analysis of the barrier-suppression
+ionization (BSI), e.g. see [Ciappina2020]_. The module was tested to reproduce the results from this paper.
+
+
+**Barrier Suppression IOnization (Ouatu)**
+Ouatu implemented a different model of BSI initially proposed in [Kostyukov2018]_ which was used in the
+study [Ouatu2022]_.
+
+
+
+To Be removed
+""""""""""""""""""""
+
 In :program:`Smilei`, following [Nuter2011]_, the ionization rate of :eq:`ionizationRate1`
 is computed for :math:`\vert m \vert=0` only.
 Indeed, as shown in [Ammosov1986]_, the ratio :math:`R` of the ionization rate
@@ -340,3 +400,11 @@ References
 .. [Schroeder2014] `C. B. Schroeder, J.-L. Vay, E. Esarey, S. S. Bulanov, C. Benedetti, L.-L. Yu, M. Chen, C. G. R. Geddes, and W. P. Leemans, Phys. Rev. ST Accel. Beams 17, 101301 <https://journals.aps.org/prab/abstract/10.1103/PhysRevSTAB.17.101301>`_
 
 .. [Gibbon] P. Gibbon, Short Pulse Laser Interactions with Matter - An Introduction, Imperial College Press (2005)
+
+.. [Tong2005] `Tong X. M., Lin C. D., J. Phys. B: At. Mol. Opt. Phys. 38 2593 (2005) <https://iopscience.iop.org/article/10.1088/0953-4075/38/15/001>`
+
+.. [Ciappina2020] `M. F. Ciappina, S. V. Popruzhenko., Laser Phys. Lett. 17 025301 (2020) <https://iopscience.iop.org/article/10.1088/1612-202X/ab6559>`
+
+.. [Kostyukov2018] `I. Yu. Kostyukov, A. A. Golovanov, Phys. Rev. A 98, 043407 (2018) <https://journals.aps.org/pra/abstract/10.1103/PhysRevA.98.043407>`
+
+.. [Ouatu2022] `I. Ouatu et al, Phys. Rev. E 106, 015205 (2022) <https://journals.aps.org/pre/abstract/10.1103/PhysRevE.106.015205>`

ionization_module_tests_and_man/cmap_map.py

@@ -0,0 +1,33 @@
+import matplotlib
+import numpy as np
+import matplotlib.pyplot as plt
+
+def cmap_map(function, cmap):
+    """ Applies function (which should operate on vectors of shape 3: [r, g, b]), on colormap cmap.
+    This routine will break any discontinuous points in a colormap.
+    """
+    cdict = cmap._segmentdata
+    step_dict = {}
+    # Firt get the list of points where the segments start or end
+    for key in ('red', 'green', 'blue'):
+        step_dict[key] = list(map(lambda x: x[0], cdict[key]))
+    step_list = sum(step_dict.values(), [])
+    step_list = np.array(list(set(step_list)))
+    # Then compute the LUT, and apply the function to the LUT
+    reduced_cmap = lambda step : np.array(cmap(step)[0:3])
+    old_LUT = np.array(list(map(reduced_cmap, step_list)))
+    new_LUT = np.array(list(map(function, old_LUT)))
+    # Now try to make a minimal segment definition of the new LUT
+    cdict = {}
+    for i, key in enumerate(['red','green','blue']):
+        this_cdict = {}
+        for j, step in enumerate(step_list):
+            if step in step_dict[key]:
+                this_cdict[step] = new_LUT[j, i]
+            elif new_LUT[j,i] != old_LUT[j, i]:
+                this_cdict[step] = new_LUT[j, i]
+        colorvector = list(map(lambda x: x + (x[1], ), this_cdict.items()))
+        colorvector.sort()
+        cdict[key] = colorvector
+
+    return matplotlib.colors.LinearSegmentedColormap('colormap',cdict,1024)

ionization_module_tests_and_man/input.py

@@ -0,0 +1,210 @@
+import numpy as np
+# import sys
+# sys.path.append('/home/ent/Smilei_ionization/simulations/')
+
+
+l0 = 2.0*np.pi                # laser wavelength
+t0 = l0                       # optical cicle
+Lsim = [12.*l0]               # length of the simulation
+Tsim = 24.*t0                 # duration of the simulation               
+resx = 128.                    # nb of cells in one laser wavelength
+dx = dy = dz = l0/resx
+cell_length  = [dx]
+dt  = 0.95*dx/np.sqrt(1.)
+rest = t0/dt
+timesteps = int(Tsim/dt)
+
+Z = 18
+mass = 39.948
+
+
+# laser pulse input parameters
+a0 = 100.
+omega = 1.
+N_cycles = 10.
+xf = Lsim[0]/2.
+tcenter = N_cycles*l0
+
+Lmu  = 0.8                              # mu-m
+
+c = 299792458.                          # m/sec
+omega_ref = 2*np.pi*c/(Lmu*1.e-6)       # 1/sec
+eps0 = 8.8541878128e-12                 # F⋅m^−1
+charge_SI = 1.60217663e-19              # C
+m_e_SI = 9.1093837e-31                  # kg
+N_r = eps0*m_e_SI*omega_ref**2/charge_SI**2
+print('Reference density = ',N_r*10**(-6), ' cm-3')
+
+I_S = 4.65e29                #W/cm^2
+l_C =  3.8615901e-13         #m
+I_ref = 0.5*(a0*omega_ref*l_C/c)**2.*4.65e29
+print('Reference intensity = ',I_ref, ' W/cm^2')
+
+
+# n0 = 5.e13/(N_r*1.e-6)  # initial density 5.e13 cm^(-3)
+# nppc = 16                               # nb of particle per cells
+
+
+n0 = 5.e13/(N_r*1.e-6)
+thickness = 0.25*l0
+position = Lsim[0]/2.-thickness/2.
+def nf(x):
+    return n0*np.heaviside(x-position, 1.)*np.heaviside(position+thickness-x, 1.)
+nppc = 128
+
+Main(
+    geometry = "1Dcartesian",
+
+    interpolation_order = 2 ,
+
+    cell_length = cell_length,
+    grid_length  = Lsim,
+
+    number_of_patches = [1],
+
+    timestep = dt,
+    simulation_time = Tsim,
+
+    reference_angular_frequency_SI = omega_ref,
+
+    EM_boundary_conditions = [ ["silver-muller", "silver-muller"] ],
+
+)
+
+def sin2_envelope(t, tcenter, N):
+    if (-np.pi*N <= t-tcenter) and (t-tcenter <= np.pi*N):
+        return np.sin((t-tcenter+np.pi*N)/2./N)**2
+    else:
+        return 0.
+
+
+
+Laser(
+    box_side = "xmin", 
+    omega = omega,
+    space_time_profile = [ lambda t: 0., 
+                           lambda t: a0*sin2_envelope(t, tcenter, N_cycles)\
+                                       *np.cos(omega*((t-tcenter)-xf)) ],  
+)
+
+
+Species(
+    name = 'atom_tunnel',
+    ionization_model = 'tunnel',
+    ionization_electrons = 'electron',
+    atomic_number = Z,
+    position_initialization = 'regular',
+    momentum_initialization = 'cold',
+    particles_per_cell = nppc,
+    mass = mass*1836.0,
+    charge = 0.,
+    # number_density = n0,
+    number_density = nf,
+    boundary_conditions = [['remove','remove']]
+)
+
+
+Species(
+    name = 'atom_TL',
+    ionization_model = 'tunnel_TL',
+    ionization_tl_parameter = 6,
+    ionization_electrons = 'electron',
+    atomic_number = Z,
+    position_initialization = 'regular',
+    momentum_initialization = 'cold',
+    particles_per_cell = nppc,
+    mass = mass*1836.0,
+    charge = 0.,
+    # number_density = n0,
+    number_density = nf,
+    boundary_conditions = [['remove','remove']]
+)
+
+Species(
+    name = 'atom_BSI',
+    ionization_model = 'tunnel_BSI',
+    ionization_electrons = 'electron',
+    atomic_number = Z,
+    position_initialization = 'regular',
+    momentum_initialization = 'cold',
+    particles_per_cell = nppc,
+    mass = mass*1836.0,
+    charge = 0.,
+    # number_density = n0,
+    number_density = nf,
+    boundary_conditions = [['remove','remove']]
+)
+
+Species(
+    name = 'atom_full_PPT',
+    ionization_model = 'tunnel_full_PPT',
+    ionization_electrons = 'electron',
+    atomic_number = Z,
+    position_initialization = 'regular',
+    momentum_initialization = 'cold',
+    particles_per_cell = nppc,
+    mass = mass*1836.0,
+    charge = 0.,
+    # number_density = n0,
+    number_density = nf,
+    boundary_conditions = [['remove','remove']]
+)
+
+Species(
+    name = 'electron',
+    position_initialization = 'regular',
+    momentum_initialization = 'cold',
+    particles_per_cell = 0,
+    mass = 1.0,
+    charge = -1.0,
+    charge_density = 0.0,
+    boundary_conditions = [['remove','remove']],
+#   time_frozen = 2.*Tsim
+)
+
+
+# DiagScalar(
+#     every = rest
+# )
+
+
+DiagParticleBinning(
+    deposited_quantity = "weight",
+    every = 1,
+    species = ["atom_tunnel"],
+    axes = [
+        ["charge", -0.5, Z+0.5, Z+1]
+    ]
+)
+
+DiagParticleBinning(
+    deposited_quantity = "weight",
+    every = 1,
+    species = ["atom_TL"],
+    axes = [
+        ["charge", -0.5, Z+0.5, Z+1]
+    ]
+)
+
+DiagParticleBinning(
+    deposited_quantity = "weight",
+    every = 1,
+    species = ["atom_BSI"],
+    axes = [
+        ["charge", -0.5, Z+0.5, Z+1]
+    ]
+)
+
+DiagParticleBinning(
+    deposited_quantity = "weight",
+    every = 1,
+    species = ["atom_full_PPT"],
+    axes = [
+        ["charge", -0.5, Z+0.5, Z+1]
+    ]
+)
+
+DiagFields(
+    every = 100,
+    fields = ['Ex','Ey','Ez','Bx','By','Bz']
+)