Hollow Beam Optical Ponderomotive Trap for Ultracold Neutral Plasma - Dataset

Hollow Beam Optical Ponderomotive Trap for Ultracold Neutral Plasma – Dataset Sergey Saakyan December 15, 2025 Version v2.0 Abstract This dataset contains the numerical data and all information required to reproduce the results reported in the manuscript. 1 Hardware and LAMMPS LAMMPS was compiled and run using the Spack package manager with AMD Optimizing C/C++ and Fortran Compilers (AOCC) following the guidelines on the official AMD web site. The computations were performed on two machines with 16and 64-core AMD CPUs. Both machines run Manjaro Linux. Information about the operating system, CPU model, LAMMPS version, and installed packages is given below. 1.1 First Machine: 64 cores (128 threads) 1.1.1 About OS and CPU $neofetch --stdout OS: Manjaro Linux x86_64 Kernel : 6.5.13 -7 - MANJARO Packages : 1428 ( pacman ) Shell : bash 5.2.32 CPU: AMD EPYC 7742 (128) @ 2.250 GHz GPU: 82:00.0 ASPEED Technology , Inc. ASPEED Graphics Family Memory : 13616 MiB / 64136 MiB 1 1.1.2 About LAMMPS $spack -V 0.22.0. dev0 (2 bcfb72d38daa2d148b5551c21a424633339649c ) $spack load lammps@20231121 $lmp -h Large - scale Atomic / Molecular Massively Parallel Simulator - 21 Nov 2023 <...> Compiler : Clang C++ AMD Clang 16.0.3 ( CLANG : AOCC_4 .1.0 - Build # 270 2023 _07_10 ) with OpenMP 5.0 C++ standard : C ++11 MPI v3 .1: Open MPI v5 .0.1 , package : Open MPI epyc@epyc - server Distribution , ident : 5.0.1 , repo rev : v5 .0.1 , Dec 20, 2023 Accelerator configuration : OPENMP package API : OpenMP OPENMP package precision : double INTEL package API: OpenMP INTEL package precision : single mixed double Active compile time flags : -DLAMMPS_GZIP -DFFT_SINGLE -DLAMMPS_SMALLBIG sizeof ( smallint ): 32 - bit sizeof ( imageint ): 32 - bit sizeof ( tagint ): 32bit sizeof ( bigint ): 64bit <...> Installed packages : FEP INTEL KSPACE LEPTON MANYBODY MOLECULE OPENMP OPT RIGID 2 <...> 1.2 Second Machine: 16 cores (32 threads) 1.2.1 About OS and CPU $neofetch --stdout OS: Manjaro Linux x86_64 Host: MS -7 B86 5.0 Kernel : 6.6.52 -1 - MANJARO Uptime : 1 min Packages : 1520 ( pacman ) Shell : bash 5.2.32 Resolution : 1920 x1080 Terminal : / dev / pts /0 CPU: AMD Ryzen 9 5950 X (32) @ 3.400 GHz GPU: NVIDIA GeForce GT 730 Memory : 516 MiB / 32016 MiB 1.2.2 About LAMMPS $spack -V 0.22.0. dev0 (2 b51980904933f8275e9a60f164a0bc6b2eda070 ) $lmp -h Large - scale Atomic / Molecular Massively Parallel Simulator - 21 Nov 2023 <...> Compiler : Clang C ++ AMD Clang 16.0.3 ( CLANG : AOCC_4 .1.0 - Build # 270 2023 _07_10 ) with OpenMP 5.0 C++ standard : C ++11 MPI v3 .1: Open MPI v5 .0.1 , package : Open MPI Distribution , ident : 5.0.1 , repo rev : v5 .0.1 , Dec 20, 2023 Accelerator configuration : OPENMP package API : OpenMP 3 OPENMP package precision : double INTEL package API: OpenMP INTEL package precision : single mixed double Active compile time flags : -DLAMMPS_GZIP -DFFT_SINGLE -DLAMMPS_SMALLBIG sizeof ( smallint ): 32 - bit sizeof ( imageint ): 32 - bit sizeof ( tagint ): 32bit sizeof ( bigint ): 64bit <...> Installed packages : FEP INTEL KSPACE LEPTON MANYBODY MOLECULE OPENMP OPT RIGID <...> 2 Descriptions of included files toHDF.py – Simple Python script to convert the particle trajectories into MultiIndex DataFrames and save them in HDF format for further processing. This script can be called from the LAMMPS input script immediately after the main simulation. This is a new vectorized version of the script that is substantially faster. This script can be executed in the folder containing the LAMMPS dump file as: 1python toHDF.py 2e-13 50000000 20000 where 2 ·10−13 is the integration time step in seconds, 50000000 the total number of integration steps in simulation, and 20000 is the dump interval (in time steps). 2.1 MD Trajectories The DB.csv file contains the simulation database. The table contents are shown in Fig. 1. The subfolders in the MDtrajectories directory contain independent simulations with the parameters specified in the database file DB.csv. Each folder contains the following four files. 4 Figure 1: DB.csv file content. Here, N=Ni=Neis the number of particles in the simulation, Te is the initial temperature of the electrons, sigma0 is the initial plasma size in µm, w0 is the laser beam waist in µm, lis the Laguerre-Gaussian (LG) beam order, P0 is the LG beam power in kW, duration is the simulation duration in µs, dt is the step size of the simulation in seconds, steps is the number of integration steps, folder name is the name of the subfolder containing the corresponding files. input.lammps – LAMMPS input script. The main idea of this script is to create a frozen dummy particle at the origin of the simulation. The pairwise interaction between this particle and the electrons/ions reproduces the traps ponderomotive potential. In the case of spherical symmetry, this approach, together with the table potential function, substantially improves computational efficiency. log.lammps – The LAMMPS output log. This file is used to extract information about the total energy conservation of the simulation. bound frac.pkl – Precalculated bound-state fraction at each time step. The following code imports the precomputed .pkl file. 1import hickle as hkl 2time, bounded_electrons_id_list,bounded_ions_id_list, bound_el_num_list, bound_ion_num_list = hkl.load('./bound_frac.pkl') ,→ ,→ out.hdf – MD trajectories saved as a MultiIndex DataFrame in a compressed HDF file. The following example shows how to open and read this DataFrame. 1import pandas as pd 5 2df =pd.read_hdf ('./MDtrajectories/sim_num022_l16_N500_P20.0kW/out.hdf') It is also necessary to drop the frozen dummy particle with id 3. 1df.drop(3, level=2, axis=0, inplace=True) Figure 2: MultiIndex DataFrame containing MD trajectories. All quantities are in SI units. A MultiIndex DataFrame enables reasonably fast analysis with a table of over 2 million rows. This dump-step spacing is crucial for bound-state fraction calculations. To reduce the size of the dataset, unused per-particle force columns fx,fy,fz from the initial LAMMPS dump are removed and initial DataFrames are compressed with zlib. After reading the HDF file, the desired data can be selected as follows: 1electrons_in_trap =df.loc[:,:,:,'electron']['r'][df['r']<88e-6].dropna() Here we select particle type ’electron’, then select the position-vector magnitude column rand filter this column for r < 88 µm, then drop all NaN values. Store the unique time values in a separate variable: 1time =electrons_in_trap.index.get_level_values(0).unique() Then we can calculate the number of electrons remaining in the ponderomotive trap volume at each time step: 1number_of_electrons =[] 2for tin time: 3number_of_electrons.append(len(electrons_in_trap.loc[t,:,:])) 4 5#plot time vs number of electrons 6plt.plot(time*1e6, number_of_electrons) 7plt.xlabel(r'Time, $\mu$s') 8plt.ylabel('Number of electrons in PT') 6 and plot the resulting dependence (blue curve in Fig. 3). Next, we subtract the number of bound electrons from the total number of electrons in the trap. 0246810 Time, s 0 100 200 300 400 Number of electrons in PT All electrons = free + bounded Free electrons Bounded electrons Figure 3: The number of electrons in the ponderomotive trap volume. The blue curve is the total number of electrons in the trap (electrons with r < 88 µm) vs. time. The green curve is the number of bounded electrons (Rydberg states) in the trap. The number of trapped electrons equals the number of trapped ions. The orange curve is the number of free electrons in the trap. At low initial electron temperatures, three-body recombination leads to the formation of classical bound states (Rydberg states). It is necessary to subtract the number of bounded electrons from the total number of electrons in the trap to isolate the ultracold-plasma (free-electron) component. The complete code used to plot figure 3 is listed below. 1import pandas as pd 2import matplotlib.pyplot as plt 3import hickle as hkl 4import numpy as np 5 6df =pd.read_hdf ('./MDtrajectories/sim_num022_l16_N500_P20.0kW/out.hdf'),→ 7# drop dummy particle 8df.drop(3, level=2, axis=0, inplace=True) 9# read IDs of bounded electrons from precalculated pkl-file,→ 10 _,bounded_electrons_id_list,_,_,_ =hkl.load ('./MDtrajectories/sim_num022_l16_N500_P20.0kW/bound_frac.pkl'),→ 7 11 # Select electrons inside trap volume r < 88um 12 electrons_in_trap = df.loc[:,:,:,'electron']['r'][df['r']<88e-6].dropna(),→ 13 # Get time column, only unique values 14 time = electrons_in_trap.index.get_level_values(0).unique(),→ 15 # There is no information about bounded fraction in first time interval,→ 16 # get this first time interval 17 n=len(time)-len(bounded_electrons_id_list) 18 # select only time steps from n to the end. 19 electrons_in_trap =electrons_in_trap.loc[time[n:],:,:,:] 20 # Get new time column 21 time = electrons_in_trap.index.get_level_values(0).unique(),→ 22 23 #find number of all electrons and number of bounded electrons in trap volume,→ 24 number_of_electrons, bounded_fraction =[], [] 25 for iin range(len(time)): 26 number_of_electrons.append(len(electrons_in_trap.loc[time[i],:,:])),→ 27 #select only bounded electrons 28 bounded_el_df = df.loc[time[i],:,bounded_electrons_id_list[i],'electron'],→ 29 #select bounded particles in PTrap with r<88um 30 bunded_in_PTrap_df = bounded_el_df[bounded_el_df['r']<88e-6]['r'].dropna(),→ 31 #number of bounded particles in PTrap 32 f_B =len(bunded_in_PTrap_df) 33 bounded_fraction.append(f_B) 34 #get number of free (plasma) electrons 35 free_electrons =np.array(number_of_electrons) - np.array(bounded_fraction),→ 36 # plot the result 37 plt.figure(figsize=(4,3)) 38 plt.plot(time*1e6, number_of_electrons, label ="All electrons = free + bounded"),→ 39 plt.plot(time*1e6, free_electrons, label ="Free electrons"),→ 40 plt.plot(time*1e6, bounded_fraction, label ="Bounded electrons"),→ 41 plt.xlabel(r'Time, $\mu$s') 42 plt.ylabel('Number of electrons in PT') 8 43 44 plt.ylim(0,430) 45 plt.legend() 46 plt.show() 2.2 Bound State Fraction The Recombination folder contains scripts for precalculating bound-state fractions. The output of the pkl recombination id precalc fast.py script is the file bound frac.pkl. The script can be executed as follows: 1python pkl_recombination_id_precalc_fast.py ./sim_num009_l16_N200_P50.0kW --mode memory --nproc 2 --drop-id 3 ,→ ,→ The first argument ../sim num009 l16 N200 P50.0kW is the folder path containing the out.hdf file, two --mode options useful for large files or limited RAM: memory and reopen, number of physical CPU cores --nproc for parallelization, and dummy particle ID --drop-id. In this version of the script, computational performance is optimized. The script recombination fast.py is called from pkl recombination id precalc fast.py and contains the main calculation logic. The simple Bash script recombination pkl calc.sh is used to precalculate bound frac.pkl files for all subfolders of some folder specified within the script. This script also determines the number of available physical CPU cores. The out hdf free calc.sh script iterates over the subfolders in a specified directory and precomputes the out free.hdf file by running out hdf free.py Python script. The out free.hdf contains only the trajectories of free particles (bound particles exluded) which is useful for analysis purposes. It is possible to plot trajectories (e.g. the absolute value of the radius vector) of bounded particles by their IDs. In figure 4 we plot only 10% of the bound particles to get a clearer representation. We select only particles that are bounded at the end of the simulation (at t= 10 µs) and plot their full time evolution. 1import pandas as pd 2import matplotlib.pyplot as plt 3import hickle as hkl 4import numpy as np 5 6df =pd.read_hdf ('./MDtrajectories/sim_num022_l16_N500_P20.0kW/out.hdf'),→ 7# drop dummy particle 8df.drop(3, level=2, axis=0, inplace=True) 9