Collisional ionization problems

[czajah]

Hi,

I'm trying to run 2d simulation of the streamer formation. As first attempt I've created simulation with space filled by neutral species with n = 1e18 particles per cm^3 and at the center I placed "seed" fully ionized plasma with density 1e13 particles per cm^3. What i observed is strange (for me) behaviour of "seed" electrons around newly created ion-electron pairs. It looks like they are blowed out. In the picture below you can see "bubbles" and inside them pairs of newly created ions and electrons pairs marked by red and green pixels respectively (black is for "seed" electrons).

It looks for me like maybe EM fields are calculated respecting newly created electrons but not taking into account changed charges of ions? What do you think?

Dimmensions of the simulation are 1cm x 3cm and the picture show the state after 0.5ns

Also sometimes ionization "explodes". Images below shows newly created ions and electrons (black and red respectively) after 47ps, 79ps and 110ps

47ps

79ps

110ps

What I'm doing wrong?

Simulation parameters:

Main = {
  'EM_boundary_conditions': [['periodic', 'periodic'],
  ['silver-muller', 'silver-muller']],
 'EM_boundary_conditions_k': [],
 'Laser_Envelope_model': False,
 '_fillObjectAndAddToList': <function SmileiComponent._fillObjectAndAddToList(self, cls, obj, **kwargs)>,
 '_list': [<Smilei Main>],
 '_verify': True,
 'apply_rotational_cleaning': False,
 'cell_length': [0.24796961358928027, 0.24796961358928027],
 'cell_sorting': False,
 'clrw': 17,
 'custom_oversize': 2,
 'custom_region_oversize': 2,
 'every_clean_particles_overhead': 100,
 'geometry': '2Dcartesian',
 'global_factor': [],
 'grid_length': [33.723867448142116, 67.44773489628423],
 'interpolation_order': 2,
 'is_pxr': False,
 'is_spectral': False,
 'maxwell_solver': 'Yee',
 'norder': [],
 'number_of_AM': 2,
 'number_of_AM_relativistic_field_initialization': 1,
 'number_of_cells': [136, 272],
 'number_of_damping_cells': [0],
 'number_of_patches': [8, 16],
 'number_of_timesteps': 3169,
 'patch_arrangement': 'hilbertian',
 'poisson_max_error': 1e-14,
 'poisson_max_iteration': 50000,
 'print_every': 100,
 'print_expected_disk_usage': True,
 'pseudo_spectral_guardells': 0,
 'random_seed': 0,
 'reference_angular_frequency_SI': 1000000000000.0,
 'relativistic_poisson_max_error': 1e-22,
 'relativistic_poisson_max_iteration': 50000,
 'save_magnectic_fields_for_SM': True,
 'simulation_time': 500.09005268710973,
 'solve_poisson': True,
 'solve_relativistic_poisson': False,
 'time_fields_frozen': 0.0,
 'timestep': 0.15780689576746915,
 'timestep_over_CFL': None,
 'uncoupled_grids': False
}

Neutral species:

{'_list': [<Smilei Species>,
  <Smilei Species>,
  <Smilei Species>,
  <Smilei Species>],
 '_verify': True,
 'atomic_number': 10,
 'boundary_conditions': [['periodic', 'periodic'], ['remove', 'remove']],
 'c_part_max': 1.0,
 'charge': 0.0,
 'charge_density': None,
 'ionization_electrons': None,
 'ionization_model': 'none',
 'ionization_rate': None,
 'is_test': False,
 'mass': 36720.0,
 'maximum_charge_state': 0,
 'mean_velocity': [0.0, 0.0, 0.0],
 'merge_accumulation_correction': True,
 'merge_discretization_scale': 'linear',
 'merge_every': 0,
 'merge_max_packet_size': 4,
 'merge_min_momentum': 1e-05,
 'merge_min_momentum_cell_length': [1e-10, 1e-10, 1e-10],
 'merge_min_packet_size': 4,
 'merge_min_particles_per_cell': 4,
 'merge_momentum_cell_size': [16, 16, 16],
 'merging_method': 'none',
 'momentum_initialization': 'mj',
 'multiphoton_Breit_Wheeler': [None, None],
 'multiphoton_Breit_Wheeler_sampling': [1, 1],
 'name': 'N_zero',
 'number_density': 3182.6073522665843,
 'particles_per_cell': <function ppc_fun(x, y)>,
 'ponderomotive_dynamics': False,
 'position_initialization': 'random',
 'pusher': 'boris',
 'radiating': False,
 'radiation_model': 'none',
 'radiation_photon_gamma_threshold': 2.0,
 'radiation_photon_sampling': 1,
 'radiation_photon_species': None,
 'regular_number': [],
 'relativistic_field_initialization': False,
 'temperature': [1.9569511835738735e-05],
 'thermal_boundary_temperature': [1.9569511835738735e-05],
 'thermal_boundary_velocity': [0.0, 0.0, 0.0],
 'time_frozen': 0.0,
 'time_relativistic_initialization': 0.0
}

Seed plasma ions:

{'_list': [<Smilei Species>,
  <Smilei Species>,
  <Smilei Species>,
  <Smilei Species>],
 '_verify': True,
 'atomic_number': 10,
 'boundary_conditions': [['periodic', 'periodic'], ['remove', 'remove']],
 'c_part_max': 1.0,
 'charge': 1.0,
 'charge_density': None,
 'ionization_electrons': None,
 'ionization_model': 'none',
 'ionization_rate': None,
 'is_test': False,
 'mass': 36720.0,
 'maximum_charge_state': 0,
 'mean_velocity': [0.0, 0.0, 0.0],
 'merge_accumulation_correction': True,
 'merge_discretization_scale': 'linear',
 'merge_every': 0,
 'merge_max_packet_size': 4,
 'merge_min_momentum': 1e-05,
 'merge_min_momentum_cell_length': [1e-10, 1e-10, 1e-10],
 'merge_min_packet_size': 4,
 'merge_min_particles_per_cell': 4,
 'merge_momentum_cell_size': [16, 16, 16],
 'merging_method': 'none',
 'momentum_initialization': 'mj',
 'multiphoton_Breit_Wheeler': [None, None],
 'multiphoton_Breit_Wheeler_sampling': [1, 1],
 'name': 'N_plus',
 'number_density': 0.031826073522665844,
 'particles_per_cell': <function ion_ppc_fun(x, y)>,
 'ponderomotive_dynamics': False,
 'position_initialization': 'random',
 'pusher': 'boris',
 'radiating': False,
 'radiation_model': 'none',
 'radiation_photon_gamma_threshold': 2.0,
 'radiation_photon_sampling': 1,
 'radiation_photon_species': None,
 'regular_number': [],
 'relativistic_field_initialization': False,
 'temperature': [1.9569511835738735e-05],
 'thermal_boundary_temperature': [1.9569511835738735e-05],
 'thermal_boundary_velocity': [0.0, 0.0, 0.0],
 'time_frozen': 0.0,
 'time_relativistic_initialization': 0.0
}

Seed plasma electrons:

{'_list': [<Smilei Species>,
  <Smilei Species>,
  <Smilei Species>,
  <Smilei Species>],
 '_verify': True,
 'atomic_number': None,
 'boundary_conditions': [['periodic', 'periodic'], ['remove', 'remove']],
 'c_part_max': 1.0,
 'charge': -1.0,
 'charge_density': None,
 'ionization_electrons': None,
 'ionization_model': 'none',
 'ionization_rate': None,
 'is_test': False,
 'mass': 1.0,
 'maximum_charge_state': 0,
 'mean_velocity': [0.0, 0.0, 0.0],
 'merge_accumulation_correction': True,
 'merge_discretization_scale': 'linear',
 'merge_every': 0,
 'merge_max_packet_size': 4,
 'merge_min_momentum': 1e-05,
 'merge_min_momentum_cell_length': [1e-10, 1e-10, 1e-10],
 'merge_min_packet_size': 4,
 'merge_min_particles_per_cell': 4,
 'merge_momentum_cell_size': [16, 16, 16],
 'merging_method': 'none',
 'momentum_initialization': 'mj',
 'multiphoton_Breit_Wheeler': [None, None],
 'multiphoton_Breit_Wheeler_sampling': [1, 1],
 'name': 'eon_prim',
 'number_density': 0.031826073522665844,
 'particles_per_cell': <function ion_ppc_fun(x, y)>,
 'ponderomotive_dynamics': False,
 'position_initialization': 'random',
 'pusher': 'boris',
 'radiating': False,
 'radiation_model': 'none',
 'radiation_photon_gamma_threshold': 2.0,
 'radiation_photon_sampling': 1,
 'radiation_photon_species': None,
 'regular_number': [],
 'relativistic_field_initialization': False,
 'temperature': [0.0019569511835738733],
 'thermal_boundary_temperature': [0.0019569511835738733],
 'thermal_boundary_velocity': [0.0, 0.0, 0.0],
 'time_frozen': 0.0,
 'time_relativistic_initialization': 0.0
}

Species for newly created electrons:

{'_list': [<Smilei Species>,
  <Smilei Species>,
  <Smilei Species>,
  <Smilei Species>],
 '_verify': True,
 'atomic_number': None,
 'boundary_conditions': [['periodic', 'periodic'], ['remove', 'remove']],
 'c_part_max': 1.0,
 'charge': -1.0,
 'charge_density': 0.0,
 'ionization_electrons': None,
 'ionization_model': 'none',
 'ionization_rate': None,
 'is_test': False,
 'mass': 1.0,
 'maximum_charge_state': 0,
 'mean_velocity': [],
 'merge_accumulation_correction': True,
 'merge_discretization_scale': 'linear',
 'merge_every': 0,
 'merge_max_packet_size': 4,
 'merge_min_momentum': 1e-05,
 'merge_min_momentum_cell_length': [1e-10, 1e-10, 1e-10],
 'merge_min_packet_size': 4,
 'merge_min_particles_per_cell': 4,
 'merge_momentum_cell_size': [16, 16, 16],
 'merging_method': 'none',
 'momentum_initialization': 'cold',
 'multiphoton_Breit_Wheeler': [None, None],
 'multiphoton_Breit_Wheeler_sampling': [1, 1],
 'name': 'eon',
 'number_density': None,
 'particles_per_cell': 0,
 'ponderomotive_dynamics': False,
 'position_initialization': 'regular',
 'pusher': 'boris',
 'radiating': False,
 'radiation_model': 'none',
 'radiation_photon_gamma_threshold': 2.0,
 'radiation_photon_sampling': 1,
 'radiation_photon_species': None,
 'regular_number': [],
 'relativistic_field_initialization': False,
 'temperature': [],
 'thermal_boundary_temperature': [],
 'thermal_boundary_velocity': [0.0, 0.0, 0.0],
 'time_frozen': 0.0,
 'time_relativistic_initialization': 0.0}

Collisions setup:

Collisions(
    species1 = ["eon","eon_prim"],
    species2 = ["N_zero"],
    ionizing = "eon",
    debug_every = 1000
)

Smilei version: bv4.4-796-g9cb8dd39-bmaster

Initialization log:

Geometry: 2Dcartesian
 --------------------------------------------------------------------------------
     Interpolation order : 2
     Maxwell solver : Yee
     (Time resolution, Total simulation time) : (6.336859, 500.090053)
     (Total number of iterations,   timestep) : (3169, 0.157807)
                timestep  = 0.900000 * CFL
     dimension 0 - (Spatial resolution, Grid length) : (4.032752, 33.723867)
                 - (Number of cells,    Cell length)  : (136, 0.247970)
                 - Electromagnetic boundary conditions: (periodic, periodic)
                     - Electromagnetic boundary conditions k    : ( [1.00, 0.00] , [-1.00, -0.00] )
     dimension 1 - (Spatial resolution, Grid length) : (4.03, 67.45)
                 - (Number of cells,    Cell length)  : (272, 0.25)
                 - Electromagnetic boundary conditions: (silver-muller, silver-muller)
                     - Electromagnetic boundary conditions k    : ( [0.00, 1.00] , [-0.00, -1.00] )

 Vectorization: 
 --------------------------------------------------------------------------------
     Mode: off

 Patch arrangement : 
 --------------------------------------------------------------------------------

 Initializing MPI
 --------------------------------------------------------------------------------
     applied topology for periodic BCs in x-direction
     MPI_THREAD_MULTIPLE enabled
     Number of MPI process : 1
     Number of patches : 
         dimension 0 - number_of_patches : 8
         dimension 1 - number_of_patches : 16
     Patch size :
         dimension 0 - n_space : 17 cells.
         dimension 1 - n_space : 17 cells.
     Dynamic load balancing: never

 OpenMP
 --------------------------------------------------------------------------------
     Number of thread per MPI process : 16

 Initializing the restart environment
 --------------------------------------------------------------------------------

 Initializing species
 --------------------------------------------------------------------------------

     Creating Species #0: N_zero
         > Pusher: boris
         > Density profile: 2D built-in profile `constant` (value: 3182.607352)

     Creating Species #1: N_plus
         > Pusher: boris
         > Density profile: 2D built-in profile `constant` (value: 0.031826)

     Creating Species #2: eon_prim
         > Pusher: boris
         > Density profile: 2D built-in profile `constant` (value: 0.031826)

     Creating Species #3: eon
         > Pusher: boris
         > Density profile: 2D built-in profile `constant` (value: 0.000000)

 Initializing External fields
 --------------------------------------------------------------------------------
     External field Ey: 2D built-in profile `constant` (value: 0.001760)
     Parameters for collisions #0 :
         Collisions between species (2 3) and (0)
         Coulomb logarithm: 0.00
         Debug every 1000 timesteps
         Collisional ionization with atomic number 10 towards species `eon`

 Initializing Patches
 --------------------------------------------------------------------------------
     First patch created
         Approximately 10% of patches created
         Approximately 20% of patches created
         Approximately 30% of patches created
         Approximately 40% of patches created
         Approximately 50% of patches created
         Approximately 60% of patches created
         Approximately 70% of patches created
         Approximately 80% of patches created
         Approximately 90% of patches created
     All patches created

 Creating Diagnostics, antennas, and external fields
 --------------------------------------------------------------------------------
     Created TrackParticles #0: species eon_prim
         id,x,y,px,py,Ex,Ey,Bx,By,q
     Created TrackParticles #1: species eon
         id,x,y,px,py,Ex,Ey,Bx,By,q
     Created TrackParticles #2: species N_zero
         id,x,y,px,py,Ex,Ey,Bx,By,q
     Created TrackParticles #3: species N_plus
         id,x,y,px,py,Ex,Ey,Bx,By,q

 finalize MPI
 --------------------------------------------------------------------------------
     Done initializing diagnostics, antennas, and external fields

 Solving Poisson at time t = 0
 --------------------------------------------------------------------------------

 Initializing E field through Poisson solver
 --------------------------------------------------------------------------------
     Poisson solver converged at iteration: 118, relative err is ctrl = 0.95 x 1e-14
     Poisson equation solved. Maximum err = 0.00 at i= -1
 Time in Poisson : 0.08

 Applying external fields at time t = 0
 --------------------------------------------------------------------------------

 Applying prescribed fields at time t = 0
 --------------------------------------------------------------------------------

 Initializing diagnostics
 --------------------------------------------------------------------------------

 Running diags at time t = 0
 --------------------------------------------------------------------------------

 Species creation summary
 --------------------------------------------------------------------------------
         Species 0 (N_zero) created with 581040 particles
         Species 1 (N_plus) created with 13936 particles
         Species 2 (eon_prim) created with 13936 particles
         Species 3 (eon) created with 0 particles

[MickaelGrech]

Ciao! The parameters you are using look weird to me. Can you give us what you is your reference angular frequency in SI units, and to what scale it corresponds to in your problem (electron skin depth, wavelength of an electromagnetic wave, or whatever). Also, what would be the characteristic Debye length of the plasma you use as a seed, and/or the plasma you expect to get by the end of your simulation?

[czajah]

@MickaelGrech

Sorry, I gave wrong sim dimmension 1cm x 3cm actually it is 1cm x 2cm.

It's just pilot simulation to check how collisional ionization works. I want to check how the collisional ionization depends on density and temperature of neutrals and maybe external electric field.

reference_angular_frequency_SI: 1000000000000.0 hertz

number_density : 10000000000000.0 / centimeter ** 3 eon temperature : 1 kiloeV ion temperature : 0.01 kiloeV Wp : 178447213442.43182 hertz L_debye : 0.07433941996724053 millimeter cell_length : 0.07433941996724053 millimeter

time_step : 0.15780689576746915 picosecond simulation_time : 5.000000000000001e-07 millisecond simulation_steps : 3169

x_simulation_cells : 136 y_simulation_cells : 272

EM_boundary_conditions : [['periodic', 'periodic'], ['silver-muller', 'silver-muller']] boundary_conditions : [['periodic', 'periodic'], ['remove', 'remove']] macroparticles_per_cell : 16

Ey : 30.0 kilovolt / centimeter

[mccoys]

The charge should be conserved during ionization. Have you verified ?

There can be a number of instabilities if you do not have enough particles per cell for instance.

[czajah]

@mccoys charge is conserved (based on "q" attribute from TrackParticles diagnostics). I'll try with increased number of particles per cell

[mccoys]

When charge is conserved, I mean the total charge. The q of an ion should change.

[czajah]

When charge is conserved, I mean the total charge. The q of an ion should change.

The 'q' of an ion is changing. At the end of simulation the sum of 'q' for N_zero is equal to -1 * sum of 'q' for 'eon'

[MickaelGrech]

I could not go through all the parameters in your namelist, but I've checked a few and they look fine (no CFL violation or obvious mistake in your normalisation). However, your neutral plasma is very dense and if collisional ionization does lead to it being effectively ionized, you will get a very dense quite cold plasma. As a consequence, your resolution will not be sufficient to resolve the Debye length associated to this dense and cold plasma, and numerical heating will kick in. You should check the evolution of the electric field in your simulation, and most important, the total energy and energy balance (from Scalar diagnostics).

[mccoys]

@czajah please let us know if the change of particles per cell or the change of resolution helped.

[czajah]

@mccoys ok, I've been doing something else for the past few days but I'm just sitting down to this problem

[czajah]

@mccoys I've tried 16, 32, 64 particles per cell and cell length 1, 0.5, 0.2 Debye length with no visible change in observed effect (bubbles). Pictures bellow are for 64 particles per cell and cell length 1 Debye length after 0.5ns.

Cells with electrons created by ionnization:

Number of particles per cell:

Maximum electric field at particle position in particular cells in Er units:

Energy is well conserved.

[mccoys]

It looks like new electrons create a space-charge field that removes other electrons. What does the total density look like ?

[czajah]

@mccoys for me it looks like cells in which ionization occured suck in other electrons. Look at the picture with "Number of particles per cell". In the center of bubbles you can see cells with a very high number of particles up to a thousand.

Number density of neutral gas is 1e18 particles/cm3 Number density of seed plasma is 1e13 particles/cm3

[mccoys]

Do you have a plot of total charge density ?

[czajah]

Do you have a plot of total charge density ?

@mccoys not yet. Give me a minute

[czajah]

@mccoys Sum of the charge per cell. Is that what you asked for?

[mccoys]

I asked for the charge density. What did you compute ? Maybe the easier way would be for you to send your input file. Is it possible ?

[czajah]

@mccoys

I asked for the charge density. What did you compute ?

I've summed charges ('q') of all species in particular cells. So it's total charge per cell.

Maybe the easier way would be for you to send your input file. Is it possible ?

I've to clean it a bit first because it is generated from template and uses my local libraries

[mccoys]

I've summed charges ('q') of all species in particular cells. So it's total charge per cell.

It should be the sum of charge*weight. Why don't you directly use the happi tool ? S.Field(0,"Rho").plot. That's it

EDIT: You should actually include either Field or Probe diagnostics in your simulation. It helps a lot.

[czajah]

@mccoys

Input file:

input.py.txt

Another input file for the explosion case from my first post:

explosion.py.txt

You should actually include either Field or Probe diagnostics in your simulation. It helps a lot.

At some point more convienient for me was using TrackParticles diagnostics with chargeless "probe particles" and I used to just checking fields at particle positions.

[czajah]

@mccoys

When I calculated charge*weight charge density start to look normal. I'm confused because I thought that particle weights after initialisation remains constant during whole simulation.

EDIT: ok, now I understand my mistake. I feel so stupid. Anyway I can't wait to understand what I did wrong in case of explosion

[mccoys]

@czajah Thanks to your input file, I have found a small bug in collisional ionization. The new electrons are not set exactly at the ion position so that there is a space charge created. I will push a fix soon.

[mccoys]

@czajah The bug should be fixed in the lastest commit available here on github (master branch). Please report back when you have checked whether it works or not.

[czajah]

@mccoys

I already did something like this fix but with hardcoded idim for my 2d geometry and I found new problem:

#3    Object "./smilei", at 0x558a8b8d255e, in 
#2    Object "./smilei", at 0x558a8b858c19, in VectorPatch::applyCollisions(Params&, int, Timers&)
#1    Object "./smilei", at 0x558a8b67b9a7, in Collisions::collide(Params&, Patch*, int, std::vector<Diagnostic*, std::allocator<Diagnostic*> >&)
#0    Object "./smilei", at 0x558a8b824590, in Patch2D::getPrimalCellVolume(Particles*, unsigned int, Params&)
Segmentation fault (Address not mapped to object [(nil)])

It happens when there is like a dozen new electrons created. When I use number_of_patches=[1,1] it works without crash.

[mccoys]

Well, it looks like you have changed something in the arrays that should not be changed. I recommend you use my patch. Otherwise let me look at your changes.

[czajah]

@mccoys

I'll pull your patch later.

Mine fix:

            new_electrons.Momentum[0].back() *= pr;
            new_electrons.Momentum[1].back() *= pr;
            new_electrons.Momentum[2].back() *= pr;

        cout <<"p"<<endl;
            new_electrons.Position[0].back() = pi->position(0, ii);
        cout<<"p0"<<endl;
            new_electrons.Position[1].back() = pi->position(1, ii);
        cout<<"p1"<<endl;
//            new_electrons.Position[2].back() = pi->position(2, ii);
//      cout<<"p2"<<endl;
        cout<<"popopo"<<endl;

            // Correction for moving back to the lab frame
            pr = w+1. - pr*gamma_s;
            new_electrons.Momentum[0].back() += pr * pi->momentum( 0, ii );

Right now I'm plaing with some other aspects of ionization code.

[czajah]

@mccoys I have your version but unfortunately simulation still blows up :(

[mccoys]

It does not not blow up when I run it. What exactly does it do ?

[czajah]

@mccoys

What exactly does it do ?

Energy rise thousands times.

Electric field at particles positions around newly created ion-electron pairs is still hundreds times bigger than average.

I have questions about code in CollisionalIonization.cpp

1)

F. Pérez et al., Phys. Plasmas 19, 083104 (2012) describes two cases We >= Wi and We < Wi

This (line 280):

        // Lose incident electron energy
        if( U2 < Wi/We ) {

looks like part of We >= Wi case but then where is ionization?

And this (line 255):

        // Ionize the atom and create electron
        if( U2 < We/Wi ) {

looks like part of We < Wi case but then where is slowing down part?

2)

Why Zstar is always incremented (line 299) even if ionization haven't occured?

[mccoys]

I cannot reproduce the blowing up. Are you sure you recompiled and run properly?

1. You have inverted the two cases, but is that the issue you have?

2. Zstar, here, is not the actual ion charge. It is used to test the cumulative ionisation probability for each charge state. The ion is not always ionised (only when U2 < We/Wi) because it is a Monte Carlo approach to particles having different weights. It is true that the order of loops and ifs is not optimal here

[czajah]

@mccoys

I'm sure I recompiled and run properly.

1. when I looked again at this in the morning I realized that that was a stupid question

It looks like I'm doing here something wrong. I think we can close this issue. Thanks for help!