pyplasma 0.4.3

Python module for plasma formation modelling.

pyplasma

Yet another 3D FDTD solver in Python. The focus of this package is the modelling of plasma formation in intense laser irradiated dielectric medium. Field ionization (Keldysh model) is implemented, along with three models for avalanche ionization trigger by impact ionization. Single rate equation (SRE), multiple rate equations (MRE) and the novel delayed rate equations (DRE) model which paper is available here: https://arxiv.org/abs/1906.08338. All figures of the DRE paper can be reproduced using the scripts under examples/paper_figures/

CPU calculations optimized with the numpy library as well as GPU calculations with Pytorch are implemented, based upon the work of http://github.com/flaport/fdtd.

Notable features are:

• 0D, 1D, 2D and 3D solvers
• CPU or GPU calculations
• Easily deployable real time data visualization while running simulation
• Multiple modes for fast calculation of Keldysh ionization rates (Pre-calculated interpolation tables, polynomial fits), alog with brute force
• Non-linear optics (2nd and 3rd order) based upon https://arxiv.org/abs/1603.09410
• Surface roughness (white noise, perlin noise or fractal noise)

Installation

Pyplasma is available on pip

pip install pyplasma

For manual installation of the developpement version,

git clone https://github.com/jldez/pyplasma.git 
cd pyplasma
python setup.py develop --user

Get started

A few examples are available in the examples/ directory.

cd pyplasma/examples
python one_dimension.py

If Pytorch is installed, you can try the three_dimensions.py script, or the self_reconfiguration.py script that reproduces the results of https://arxiv.org/abs/1702.02480

API description

Import

In a python 3.6+ environnement, start by importing pyplasma

from pyplasma import *

Backend

Then, set the backend. For 0D, 1D and small 2D simulations, CPU calculations are prefered:

set_backend('numpy')

For 3D calculations, if a GPU with CUDA is available and Pytorch is installed, use instead:

set_backend('torch.cuda')

Simulation domain

The simulation domain is created using the Domain class. For 0D, no arguments are needed. A N dimensions domain is created with N elements lists in the arguments grid and size. For example, a 3D domain can be created with

dom = Domain(grid=[50,128,128], size=[1*um,3*um,3*um], pml_width=200*nm)

For now, the boundaries can only be PMLs at both the x-axis ends and periodic along the y and z axes. The posibility to customize this is planned for a future update.

Laser source

The laser source is created with the Laser class

laser = Laser(wavelength=800*nm, pulse_duration=10*fs, fluence=7e4, t0=10*fs)

Its amplitude can be set directly with the E0 argument or indirectly set by a pulse fluence if the pulse duration (FWHM) is set and finite. The t0 argument is the time correspoding to the peak of the gaussian enveloppe. The laser source is attached to the simulation domain with

dom.add_laser(laser, position='default', source_mode='tfsf', ramp=True)

The default position is right next to the PML at the low x-axis boundary. Propagation is directed towards +x and the polarization is along the z axis. This will be made customizable in a future update. A ramp can be added to smoothly turn on the laser source and prevent numerical artefacts.

Material

Optics and non-linear optics

The dielectric medium parameters can be set using the Material class, for example

material = Material(index=1.45, resonance=120*nm)

The refractive index and the resonance are the only two mandatory parameters. If resonance is not a wanted element of the simulation, simply choose a value far from the laser wavelength. Non-linear optics is included with the 2nd and 3rd order susceptibilities chi2 and chi3 arguments.

Drude model

The free currents from charge carriers are calculated based upon the Drude model. Add this to the simulation with

material = Material(index=1.45, resonance=120*nm, drude_params={'damping':1e15, 'rho':3e26, 'm_CB':1, 'm_VB':2})

The drude_params arguments takes a dictionary. The plasma damping is set with the damping key, the plasma density (electron-hole pairs density) with the rho key, the effective mass of the electrons in the conduction band with the m_CB key (relative to the electron's free mass) and the effective mass of the holes in the valence band with the m_VB key (relative to the electron's free mass).

Ionization

Plasma formation is activated by adding the ionization_params argument when creating the Material instance:

material = Material(index=1.45, resonance=120*nm, drude_params={'damping':1e15},
                    ionization_params={'rate_equation':'dre','bandgap':9*eV,'density':2e28,'cross_section':1e-19})

The rate_equation key is used to set the ionization model. It can be set to

• fi for field ionization (Keldysh) only
• sre for the single rate equation model (the key alpha_sre is then necessary to set the impact ionization rate in m^2/J)
• mre for the multiple rate equations model (the key cross_section is then necessary to set the carrier-neutral collisional cross-section in 1/m^2)
• dre for the delayed rate equations model (the key cross_section is then necessary to set the carrier-neutral collisional cross-section in 1/m^2)

For all cases, the bandgap (in J) and the molecular density (in 1/m^3) of the material has to be indicated. The plasma density will saturate at density, as only single ionization per molecule is accounted for (for now).

A electron-hole pair recombination rate can be set with the recombination_rate key.

An optional correction to the Keldysh model to account for plasma damping can be toggle on with the key fi_damping_correction set to True. This is unpublished work and should be kept off while it is still under investigation. Or be stolen and published by anyone, be my guest!

The Keldysh field ionization model is extremely computationaly expansive. While it can be calculated by brute force at every FDTD grid cells and time steps with fi_mode key set to brute, it is not reasonable. By default, the fi_mode is set to linear, which will pre-calculate an interpolation table and perform linear interpolation instead. The size of the table is set with the fi_table_size key (default is 1000) and will work for electric field amplitudes between 10^3 V/m and 3*E0 (E0 is the peak amplitude of the laser). However, the linear interpolation is performed using the CPU and may cause a speed bottleneck for GPU calculations. It becomes less of a problem for large 3D grids, as the rest of the calculations catch up in computationnal cost. For small grids, CPU is prefered anyway, but for medium grid sizes, one may want to try to set fi_mode to nearest. This is less accurate, but fully works on GPU. However, it is extremely memory expansive and the fi_table_size will probably have to be reduced considerably (which hurts accuracy even more). Finaly, to fix all these issues at the cost of approximating the Keldysh model to a polynomial fit, the fi_mode can be set to fit, which is quite fast, memory efficient and also works on GPU.

Add material to the simulation

The material instance just created is added to the simulation (after adding the laser source to the domain) with

dom.add_material(material, boundaries={'xmin':300*nm})

If the domain is 1+ dimensions, the material is added everywhere, unless boundaries are specified with the keys xmin, xmax, ymin, ymax, zmin and zmax.

Surface roughness

To add surface roughness to the xmin boundary (other boundaries not possible yet), use (after adding material to the domain)

surface_roughness(material, boundary='xmin', amplitude=20*nm, noise='fractal', feature_size=100*nm, show=True)

The amplitude is the maximum thickness added to the surface. The roughness is generated from a random 2D map (that can be visualized at the start of the simulation with show set to True). The randomness of that map can be tuned with the noise argument set to

• white for a random height between 0 and amplitude at each grid cell
• perlin for Perlin noise, for which the characteristic bump sizes is set with feature_size
• fractal for layered Perlin noise maps, for which the characteristic bump sizes is set with feature_size

Data extraction from simulations

Under contruction. The API is in its early phase and will ultimately change. However, one can play around with what is available based upon the examples.

Run simulation

The simulation is ran using

dom.run(15*fs)

The total time of the simulation is the only required argument. The time steps are automatically set, but stability can be compromised by various elements added to the simulation. If stability issues occur, the time steps can be reduced with the argument stability_factor that multiplies the time steps (default is 0.95).