# Evaluating ion dynamics through Coulomb and Yukawa interaction potentials in one component strongly coupled plasmas

Swati Swagatika Mishra, Sudeep Bhattacharjee, Pascal Brault

► To cite this version:

Swati Swagatika Mishra, Sudeep Bhattacharjee, Pascal Brault. Evaluating ion dynamics through Coulomb and Yukawa interaction potentials in one component strongly coupled plasmas. Physical Review E, 2025, 111 (1), pp.015208. 10.1103/PhysRevE.111.015208. hal-04676799v2

HAL Id: hal-04676799

https://hal.science/hal-04676799v2

Submitted on 11 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Evaluating ion dynamics through Coulomb and Yukawa interaction potentials in one
component strongly coupled plasmas

Swati Swagatika Mishra\* and Sudeep Bhattacharjee\dagger

Department of Physics, Indian Institute of Technology Kanpur, Uttar Pradesh 208016, India

Pascal Brault‡

GREMI UMR 7344 CNRS Université d'Orléans, Cedex 2 France

(Dated: January 9, 2025)

Atmospheric pressure helium plasmas are investigated through molecular dynamics simulations at room temperature (300 K) for various ionization fractions ($\chi_i = 10^{-1} - 10^{-5}$) in the strongly coupled regime (ion coupling parameter, $\Gamma_i \sim 1 - 10$) employing Coulomb and Yukawa interaction potentials. The role of electron screening in ion dynamics and energetics is examined through ion and gas temperatures, mean squared displacement of ions, ion coupling parameter, and radial distribution function of the system. It is found that electron screening in the Yukawa potential significantly limits the disorder-induced heating (DIH) mechanism for strongly ionized plasmas ($\chi_i \ge 10^{-3}$). Whereas, ions show a prominent sub-diffusive behavior associated with the DIH during the non-equilibrium phase for the Coulomb potential. The DIH mechanism is explained using a model based upon the conservation of energy. However, for weakly ionized plasmas ($\chi_i \le 10^{-4}$), the maximum ion temperatures are almost similar for both potentials. Furthermore, electron screening affects the separation distance and arrangement of the ion-neutral pairs for all the values of $\chi_i$. In general, Yukawa potential results in a lower mean potential energy of the interacting particles, which is energetically favourable for the stability of the system.

# I. INTRODUCTION

Plasmas at atmospheric pressure have been studied extensively in a wide range of experiments and found their applications in surface processing and biomedical fields [1–4]. The criticality of their research lies in the fact that it has led to the emergence of a new field known as plasma medicine [5, 6]. These plasmas are often generated in laboratories at room temperatures employing RF sources [7], microwaves [8], DC sources [9], nano-second pulsed discharges [10], and dielectric barrier discharges [11]. Depending on the mode of generation, the value of ionization fraction ($\chi_i = n_i/(n_i + n_n)$, where $n_i$ and $n_n$ are the ion and gas number densities, respectively) in the atmospheric pressure plasmas, ranges from $10^{-1}$ to $10^{-8}$, and divides the plasma into strong, moderate and weakly ionized regimes. Based on their ionization fraction, the plasmas show significant heating and strong coupling among the ions after their generation [9, 12, 13], where the ion coupling parameter ($\Gamma_i$), which is the ratio of mean ion potential energy to the mean ion kinetic energy becomes $\ge 1$. In strongly ionized plasmas, the ion heating is assumed to be caused by the disorder-induced heating (DIH) mechanism [12]. DIH arises as the ions tend to move to a lower potential energy state from a disordered state and release kinetic energy in the process. However, for weakly ionized plasmas, the heating is often attributed to the discharge geometry and type of plasma source [9, 13]. For the moderate and weakly

ionized plasmas, most of the simulation and theoretical efforts are focused on understanding the plasma chemistry, electron temperature ($T_e$), and electron densities ($n_e$) for a given plasma geometry and generation scheme [9, 13, 14]. Furthermore, the contributions of intra and inter-species interactions and internal structures of the bulk plasmas in plasma heating and ion coupling strength are not very well understood, even for the strongly ionized plasmas. In this context, we employ molecular dynamics (MD) simulations to investigate the ion and neutral gas dynamics and their implications in ion energies, ion coupling strength, ion and neutral temperatures, and the internal structures in an atmospheric pressure helium plasma for $\chi_i (= 10^{-1} to 10^{-5})$, at room temperature (300 K), thus covering the strong, moderate and weakly ionized plasma regimes. Here, an instantaneous ionization is considered, unlike the real experiments where the ionization process (generation and growth of the ions) may change depending upon the method of plasma generation, such as DC, RF, microwave, etc. Therefore, we choose the total ionization fraction as the parameter, which may be realized in different plasma systems. MD simulations efficiently track each particle's trajectories and easily convey molecular-level information [15–18]. However, while employing the MD simulations, it is imperative to designate the intra and inter-species interactions with appropriate potentials. Therefore, in this work, we discuss the performance of both Coulomb and Yukawa potentials as the probable ion-ion interaction potential for atmospheric pressure plasmas, which may help us to decide the appropriate potential, depending upon $\chi_i$ and other plasma parameters. Section II of this article discusses the inter-particle potentials in-depth, along with the simulation details, followed by the important

<footnote>\* swatis@iitk.ac.in</footnote>

<footnote>† sudeepb@iitk.ac.in</footnote>

<footnote>‡ pascal.brault@univ-orleans.fr</footnote>
results outlined in section III and a discussion in section IV. The article is concluded in section V.

## II. SIMULATION DETAILS

Typically, the plasmas consist of neutrals, ions, and electrons, and in atmospheric pressure plasmas, the gas temperature ($T_g$) ~ ion temperature ($T_i$), and the electron temperature, $T_e \gg T_i$. Therefore, in most plasma simulations, one component plasma (OCP) model is used [12, 19, 20], and electrons are not explicitly considered in the simulation system. Commonly, for atmospheric pressure plasmas, ionic interactions (which are singly charged He ions in the study) are depicted by the Coulomb potential [12, 21], given by

$$ \phi_C(r) = \frac{1}{4\pi\epsilon_0} \frac{q^2}{r}, \qquad (1) $$

where $\epsilon_0$ is the permittivity of free space, $q$ is the charge of the ions, and $r$ is the separation distance between a pair of ions. This model is seen to work well for quantifying the DIH for strongly ionized argon plasmas with $\chi_i \ge 10^{-2}$ [12]. In the present simulation work, this model will be tested for He plasmas at various ionization fractions. However, in Coulomb potential, the effect of the screening provided by the background electrons is fully neglected, which may have a role in the ion dynamics, as seen in the case of ultracold neutral plasmas [22]. Therefore, the screened Coulomb (Yukawa) potential, given by

$$ \phi_Y(r) = \frac{1}{4\pi\epsilon_0} \frac{q^2 e^{-\kappa r}}{r}, \qquad (2) $$

is also employed to quantify the screening effect at different values of $\chi_i$. Here, $\kappa$ is the screening constant due to the electrons and is calculated assuming the ions and the electrons have the same densities ($n_i = n_e$), using the relation,

$$ \kappa^2 = \frac{n_e q^2}{\epsilon_0 k_B T_e}. \qquad (3) $$

Furthermore, the electron temperature ($T_e$) is chosen to be 1 eV (11600 K) in this work, which is a typical value of $T_e$ in atmospheric pressure plasmas [9]. The values of $\kappa$ are tabulated in Table II. It is observed that for lower values of $\chi_i$, the values of $\kappa \sim 10^{-4}$ Å$^{-1}$. Therefore, it is expected that the effect of electron screening diminishes for weakly ionized plasmas. The inverse of $\kappa$, known as the Debye radius ($r_d$), is $> a_{ii}$ for all the values of $\chi_i$, to ensure effective screening within the Debye sphere [23, 24].

The plasma properties are also expected to be influenced by the neutral gas. Therefore, ion-neutral and neutral-neutral interactions are incorporated in the simulation. The ion-neutral interactions are guided by the

modified charge-induced dipole potential [12], given by

$$ \phi_{ind}(r) = \frac{q^2}{8\pi\epsilon_0} \frac{\alpha_R a_B^3}{r^4} \left( \frac{r_\phi^8}{r^8} - 1 \right), \qquad (4) $$

where, $\alpha_R$ is the relative polarizability for He (= 1.38 [25]), $a_B$ is the Bohr's radius, and $r_\phi$ is the radius at which the repulsive core acts. The parameter $r_\phi$ (= $ca_{in}$, where $a_{in}$ is the interparticle separation distance for an ion and neutral pair, and $c$ is a multiplicative constant) is calibrated using similar techniques followed by Acciarri et al. [12]. For $T_g = 300$ K, $c$ is found to be 0.133.

The neutral-neutral interactions are modeled by Lennard-Jones potential given by,

$$ \phi_{LJ}(r) = 4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right], \qquad (5) $$

where $\epsilon$ is the energy parameter (0.00094 eV for He), and $\sigma$ is the length parameter (2.64 Å for He) [16].

All the three-dimensional MD simulations are carried out using the LAMMPS software package [15]. The maximum number of ions ($N_i$) is 10,000 for $\chi_i = 10^{-1}$, whereas for $\chi_i \le 10^{-4}$, $N_i = 100$ is fixed. The system sizes ($L$) are chosen in such a way that atmospheric gas pressure is maintained and $L \gg r_d$, which ensures an important plasma criterion [23, 24]. The neutral He particles are distributed randomly in the system and assigned with an initial Gaussian velocity distribution corresponding to the mean temperature $T_g = 300$ K. The He atoms are then equilibrated at constant temperature in a canonical (NVT, where particle number ($N$), volume ($V$), and temperature ($T$) are constant) ensemble using a Nose-Hoover thermostat [26]. After the system is thermally equilibrated, a fraction of He particles are ionized, and the ions are introduced randomly in the system depending on the value of $\chi_i$. The initial velocity distribution of the He ions is of the Gaussian type with the mean temperature corresponding to 300 K. The system is then allowed to evolve in a microcanonical (NVE, where $N$, $V$, and energy ($E$) are constant) ensemble, and $T_i$, $T_g$, and ion energies are recorded. As the system evolves, $T_i$ approaches $T_g$, and thermal equilibrium is achieved when $T_i \sim T_g$.

The Coulomb potential calculations in the simulations are carried out using the particle-particle particle-mesh (P3M) solver to incorporate the long-range effects. To employ the P3M method, a cut-off length for the short-range interactions needs to be assigned, for which the force values are calculated in the real space. Furthermore, the ion-neutral potential (cf. eqn. 4) is incorporated in a tabular format in the LAMMPS input script as it is not an inbuilt potential form. To use both the potentials combinedly, their cut-off lengths need to match. Therefore, this common cut-off length is calibrated for all the values of $\chi_i$ so that the maximum possible value can be chosen without letting the simulation system collapse. For $\chi_i = 10^{-1}$, this value is 200 Å, and for $\chi_i \le 10^{-2}$, the value is set at 250 Å($\sim 10 a_{in}$).
<table><caption>TABLE I. Simulation input parameters for different ionization fractions ($\chi_i$). $N_i$ and $N_T$ are the ion numbers and total number of particles, respectively. $L$ is the system size. $n_i$ and $n_n$ are the ion and neutral densities, respectively. $a_{ii}$ and $a_{in}$ are the inter-particle spacings for ion-ion and ion-neutral pairs, respectively. $r_Y$ is the chosen cut-off length for Yukawa potential.</caption><thead><tr><th>χ<sub>i</sub></th><th>N<sub>i</sub></th><th>N<sub>T</sub> × 10<sup>5</sup></th><th>L (Å)</th><th>n<sub>i</sub> (×10<sup>22</sup>m<sup>-3</sup>)</th><th>n<sub>n</sub> (×10<sup>25</sup>m<sup>-3</sup>)</th><th>a<sub>ii</sub> (Å)</th><th>a<sub>in</sub> (Å)</th><th>r<sub>Y</sub> (Å)</th></tr></thead><tbody><tr><td>10<sup>-1</sup></td><td>10000</td><td>1</td><td>1598.7</td><td>244.75</td><td>2.2027</td><td>46.044</td><td>22.136</td><td>500</td></tr><tr><td>10<sup>-2</sup></td><td>1000</td><td>1</td><td>1598.7</td><td>24.475</td><td>2.423</td><td>99.199</td><td>21.444</td><td>500</td></tr><tr><td>10<sup>-3</sup></td><td>1000</td><td>10</td><td>3444.3</td><td>2.4475</td><td>2.445</td><td>213.72</td><td>21.379</td><td>500</td></tr><tr><td>10<sup>-4</sup></td><td>100</td><td>10</td><td>3444.3</td><td>0.24475</td><td>2.4472</td><td>460.43</td><td>21.373</td><td>600</td></tr><tr><td>2 × 10<sup>-5</sup></td><td>100</td><td>50</td><td>5889.7</td><td>0.048949</td><td>2.4474</td><td>787.32</td><td>21.372</td><td>800</td></tr></tbody></table>

<table><caption>TABLE II. Intrinsic quantities of He plasma for different ionization fractions ($\chi_i$). $\kappa$ is the electron screening constant, $\Gamma_i$ is the Coulomb coupling parameter at $t=0$ (i.e., at $T_i=300$ K), $\omega_{pi}$ is the ion-plasma frequency and $dt$ (= 0.001 ps) is the time step chosen for the simulations.</caption><thead><tr><th>χ<sub>i</sub></th><th>κ(×10<sup>-3</sup> Å<sup>-1</sup>)</th><th>Γ<sub>i</sub></th><th>ω<sub>pi</sub> (×10<sup>12</sup>rad s<sup>-1</sup>)</th><th>ω<sub>pi</sub>dt (×10<sup>-3</sup>rad)</th></tr></thead><tbody><tr><td>10<sup>-1</sup></td><td>21.03</td><td>12.1</td><td>1.0336</td><td>1.0336</td></tr><tr><td>10<sup>-2</sup></td><td>6.65</td><td>5.61</td><td>0.3268</td><td>0.3268</td></tr><tr><td>10<sup>-3</sup></td><td>2.1</td><td>2.6</td><td>0.1034</td><td>0.1033</td></tr><tr><td>10<sup>-4</sup></td><td>0.665</td><td>1.21</td><td>0.0326</td><td>0.0326</td></tr><tr><td>2 × 10<sup>-5</sup></td><td>0.297</td><td>0.7</td><td>0.0146</td><td>0.0146</td></tr></tbody></table>

The cut-off length for the Yukawa potential ($r_Y$) is chosen in such a way that for $r > r_Y$, $T_i$ and $T_g$ at equilibrium remain consistent and $r_Y > a_{ii}$. For $\chi_i = 10^{-1}$ and $10^{-2}$, a maximum possible value of $r_Y$ is chosen (= 500 Å) due to smaller values of $L$ (cf. Table I). However, for moderate and weakly ionized systems, ($\chi_i \le 10^{-3}$), the values of $L$ are larger, and therefore, we need to judiciously choose $r_Y$ so that the long-range effects are incorporated without increasing the computational expenses. For $\chi_i \le 10^{-3}$, $T_i$ and $T_g$ at equilibrium remain consistent as $\phi_Y(r)$ reaches $10^{-10}$ times its maximum value. Therefore, $r_Y$ for $\chi_i \le 10^{-3}$ is chosen as the $r$ for which $\phi_Y(r)$ reaches $10^{-10}$ times its maximum value (cf. Table I). The cut-off length for neutral-neutral interactions is chosen to be 50 Å.

$T_i$, $T_g$, and ion potential energy at equilibrium are evaluated for different values of $N_i$ and $L$ for a given value of $\chi_i$ to estimate the system size effect. The quantities are found to be consistent with the reasonable variation of system size. The chosen values of $N_i$, total number of particles ($N_T$), $L$, $n_i$, $n_n$, mean ionic separation length ($a_{ii} = a_{ws} = (3/4\pi n_i)^{1/3}$, where $a_{ws}$ is the Wigner-Seitz radius), $a_{in}$ and $r_Y$ are listed for all the values of $\chi_i$ in Table I. The values of $\kappa$, $\Gamma_i(= q^2/4\pi\epsilon_0 a_{ws} k_B T_i)$, ion plasma frequency ($\omega_{pi} = \sqrt{e^2 n_i / \epsilon_0 m_i}$, where $e$ is the elementary charge and $m_i$ is the mass of the ion), and $\omega_{pi}dt$, where $dt$ is the chosen time step for the simulations are tabulated in Table II for different values of $\chi_i$. $\Gamma_i$ is calculated for Coulomb potential at $t=0$ (at the start of the NVE stage) and $T_i = 300$ K.

The convergence of MD runs is ensured by estimating the potential energy fluctuations in the system. It is observed that the maximum percentage standard deviation

in the ion potential energy after the equilibrium has reached is ~4%. Additionally, all simulation runs are performed for different spatial distributions of neutrals and ions. For all spatial distributions, similar levels of fluctuations are observed. Furthermore, as mentioned in the result section, the maximum error in the mean ion temperature profile is 8.12 %, and the gas temperature profile is 0.2 %, which shows the convergence of the simulations within the reasonable error limits, as can be seen in the gaseous MD simulation works [17].

To validate the simulation protocols, the temperature profiles for $T_i$ and $T_g$, and the ion self-diffusion coefficient in Ar plasma are compared with those obtained from the literature (Acciarri et al. (2022) and (2023)) and are found to be matching well.

Acciarri et al. (2022) [12] focus on the strongly ionized Ar plasma and investigate the temporal evolution of $T_g$, $T_i$, ion kinetic energy, coupling parameter, and radial distribution function employing Coulomb potential for ion-ion interactions, which highlight the disordered induced heating (DIH) as one of the primary heating mechanisms in plasmas such as nanosecond pulsed plasmas. Following the non-equilibrium work, Acciarri et al. (2023) [21] study the ion self-diffusion coefficient for a broad range of $\chi_i$ in the Ar plasmas using MD and MCC simulations, focusing on the diffusion phenomena.

In contrast, our work investigates He plasmas for a broad range of $\chi_i$ using MD simulations. The work attempts to quantify the screening effect due to the electrons in plasma dynamics, energetics, and structure in atmospheric plasmas by employing both Coulomb and Yukawa potentials as the ion-ion interactions. The work investigates the underlying mechanisms of DIH and energy contributions during the heating and cooling phases
<figure><img src="image_2.png" /><figcaption>FIG. 1. Temporal evolution of ion and gas temperatures ($T_i$ and $T_g$) for different ionization fractions. (a) $\chi_i = 10^{-1}$, (b) $\chi_i = 10^{-2}$, (c) $\chi_i = 10^{-3}$ and (d) $\chi_i = 10^{-4}$. The black solid line in Figure 1(a) shows the initial system temperature = 300 K.</figcaption></figure>

of the ion temperature profiles. Furthermore, the consequences of the DIH are also estimated in terms of subdiffusive ion motion and ion coupling parameter. Lastly, the article comments on the structure of the plasma and opens up the scope for future investigations on the one component Yukawa plasmas at atmospheric pressure.

The simulations presented in the current article are run for five different initial spatial distributions of neutrals and ions for each value of $\chi_i$. The values of $T_i$ and $T_g$ are then averaged and presented in the temperature plots (cf. Figure.1).

### III. RESULTS

#### A. Temperature profiles and ion dynamics

Figure 1 shows the spatially averaged ion and neutral temperature profiles during the non-equilibrium condition in the system, i.e., $T_i \neq T_g$, for both potentials (Coulomb and Yukawa). The temperatures are sampled during the NVE run of the simulation. The curves are shown till $T_i$ approaches $T_g$ (~350 ps). The maximum error in the $T_i$ and $T_g$ curves are 8.12% and 0.2%, respectively, which are obtained from the standard deviation errors of the spatially averaged curves. It is observed that the DIH mechanism is prominent for $\chi_i \ge 10^{-3}$, which results in ion heating. $T_g$ shows an increasing trend for $\chi_i \ge 10^{-2}$, which later equilibrates with $T_i$ (cf. Figure 1(a) and (b)). Furthermore, for $\chi_i \ge 10^{-3}$, the $T_i$ curve shows significant oscillations for Coulomb potential. The
<figure><img src="image_3.png" /><figcaption>FIG. 2. Maximum ion temperature ($T_{i_max}$) and equilibrium gas temperature ($T_{g_eq}$) versus the ionization fraction ($\chi_i$).</figcaption></figure>

frequency of the oscillations decreases, whereas the width of each peak increases with a decrease in $\chi_i$. The time period within which the oscillations last is quantified in terms of plasma frequency ($f^{-1} = 2\pi/\omega_{pi}$). For $\chi_i = 10^{-1}$, $10^{-2}$ and $10^{-3}$ the oscillations last for $\sim 13f^{-1}$, $2.6f^{-1}$, and $0.9f^{-1}$, respectively. The oscillations last for a longer period for the higher value of $\chi_i$. It shows that the ions interacting with Coulomb potential suffer rapid displacements from their respective mean positions in the potential well formed by the nearest neighbors [27], resulting in fluctuating repulsive potential energy that gets converted to the oscillating kinetic energy. Furthermore, with a higher number of ions, this phenomenon lasts for longer periods due to higher repulsive energies. In contrast, such oscillations are very much damped for Yukawa potential. For example, for $\chi_i = 10^{-1}$, there is only one secondary peak, which vanishes for lower values of $\chi_i$ ($\le 10^{-3}$). This signifies that the screening due to the electrons weakens the repulsive force between the ions. For weakly ionized plasmas, the peaks gradually vanish, with $\chi_i = 2 \times 10^{-5}$ not showing any peaks for any potential (not shown). Rather, the $T_i$ curve oscillates with significant thermal and statistical fluctuations around a mean temperature. It is also observed that for $\chi_i \ge 10^{-3}$, the ions attain the maximum temperature ($T_{i_max}$) at similar times for both potentials. For example, for $\chi_i = 10^{-1}$, the ions attain $T_{i_max} = 1405.84$ K at $t = 0.57$ ps for Coulomb potential and $T_{i_max} = 792.4$ K at $t = 0.67$ ps for Yukawa potential. However, for $\chi_i \le 10^{-4}$, $T_{i_max}$ is attained faster for Yukawa potential compared to Coulomb potential. Additionally, with decreasing $\chi_i$, the time corresponding to $T_{i_max}$ increases.

Figure 2 shows that the value of $T_{i_max}$ decreases with the decrease in $\chi_i$, demonstrating that the DIH mechanism is prominent for $\chi_i \ge 10^{-3}$. $T_{i_max}$ for a constant $\chi_i$ is higher for Coulomb potential, as seen in Figure 2 for all the values of $\chi_i$, and the values start approaching

for both potentials for $\chi_i \le 10^{-4}$. This signifies that ion heating is directly influenced by electron screening. The equilibrium gas temperature at the start of the equilibrium ($T_{g_eq}$) for both potentials is also plotted in Figure 2. The values are similar for both potentials for $\chi_i < 10^{-2}$.

To study the ion dynamics during the DIH mechanism, the mean squared displacement for the ions (MSD = $\frac{1}{N_i} \sum_{j=1}^{N_i} [r_j(t) - r_j(0)]^2$, where $r_j(0)$ and $r_j(t)$ are the displacements of the $j^{th}$ ion at an initial time and a later time $t$, respectively) is calculated. MSD is calculated independently at different initial times in two scenarios: (1) during the non-equilibrium (NE) phase, as soon as the ions are introduced in the system, and (2) during the equilibrium (EQ) phase ($T_i = T_g$), for both potentials. Figure 3 shows the MSD data for $\chi_i = 10^{-1}$. It is observed that for the NE case, after a short ballistic regime (MSD $\propto t^a$, $a = 2$), the MSD shows a sub-diffusive regime with $a \sim 0.8$, before transitioning to the diffusive regime with $a = 1$ for both potentials. The time corresponding to the end of the ballistic regime also corresponds to $T_{i_max}$. Thereafter, the sub-diffusive regime starts and lasts till $T_i$ approaches $T_g$. Such a sub-diffusive regime is missing for the EQ scenario and for $\chi_i \le 10^{-4}$, even in the NE scenario. This observation confirms the fact that the sub-diffusive behavior is exclusive to the DIH mechanism and it indicates the non-thermal transport in the system. A similar observation has been reported for ultracold neutral plasma (UNP), where the ion velocity distribution during the DIH is non-thermal in nature due to the influence of nearest neighbors [28]. This phenomenon is also observed in biological specimens, where the particle is crowded by its neighbors, which hinders its random motion [29].

## B. Ion coupling parameter

The screening effect due to the electrons also has an effect on the system's energetics. The mean potential energy due to ion-ion pair interaction for Coulomb potential is given by,

$$PE_C = \frac{q^2}{4\pi\epsilon_0 a_{ws}}, \qquad (6)$$

For Yukawa potential, the mean potential energy due to ion-ion pair interaction is given by,

$$PE_Y = \frac{q^2 e^{-\kappa a_{ws}}}{4\pi\epsilon_0 a_{ws}}. \qquad (7)$$

The mean kinetic energy (KE) for both potentials is $k_B T_i$. $T_i$ depends on the chosen potential at the equilibrium condition. Subscripts “C” and “Y” are used to refer to Coulomb and Yukawa potentials, respectively. Figure 4 shows the variation of mean PE and KE with $\chi_i$ for both potentials. It is evident that the mean PE
<figure><img src="image_4.png" /></figure>

FIG. 3. Mean squared displacement for ions (MSD) at $\chi_i = 10^{-1}$ for (a) Coulomb and (b) Yukawa potential. MSD is calculated during the non-equilibrium (NE) and equilibrium (EQ) phases, with two different origin times. The curve for the NE phase shows marking for $T_{i,max}$. The curves extend till the start of the equilibrium ($T_i \sim T_g$). The respective dependence of the MSD curves on $t$ is shown in the graphs. For the sub-diffusive region, MSD $\propto t^{0.8}$, and for the diffusive region, MSD $\propto t$. The solid lines show the fitting lines.

<figure><img src="image_5.png" /></figure>

FIG. 4. Mean potential (PE) and mean kinetic energies (KE) for ions for both potentials.

<figure><img src="image_6.png" /></figure>

FIG. 5. Ion coupling parameter ($\Gamma_i$) versus the ionization fraction ($\chi_i$) for Coulomb and Yukawa potentials. The solid lines are only to guide the eyes.

for Yukawa potential is lower than that for Coulomb potential, which provides an energetically more favorable condition for the plasma system. Meanwhile, the mean KE remains almost constant for both potentials except for the $\chi_i = 10^{-1}$ case.

The screening effect on the system's energetics is quantified through the coupling parameter of the ions ($\Gamma_i$). For Coulomb potential, $\Gamma_{iC} = PE_C/KE_C (=q^2/4\pi\epsilon_0a_{ws}k_B T_i)$ and for Yukawa potential $\Gamma_{iY} = PE_Y/KE_Y (=q^2e^{-\kappa a_{ws}}/4\pi\epsilon_0a_{ws}k_B T_i)$ [30].

It is observed that $\Gamma_i$ decreases with a decrease in $\chi_i$ (cf. Table II and Figure 5) due to the reduced number

of ions for the lower values of $\chi_i$. The values of $\Gamma_i$ for Coulomb potential are higher than those for Yukawa potential due to the screening effect (cf. Figure 5). This finding agrees with the observations for the ultracold neutral plasmas [22, 31].

### C. Internal plasma structure

To probe into the internal structure of the plasma, the radial distribution function ($g(r) = Vn(r)/4N\pi r^2\Delta r$),
<figure><img src="image_7.png" /><figcaption>FIG. 6. Radial distribution function for (a) ion-ion pair and (b) ion-neutral pair, for ionization fraction ($\chi_i$) = $10^{-1}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 7. Change in temperature ($\Delta T = T(t + \Delta t) - T(t)$) versus time for (a) Coulomb potential and (b) Yukawa potential, for ionization fraction ($\chi_i$) = $10^{-1}$. The time step ($\Delta t$) used for the calculations is 0.05 ps.</figcaption></figure>

where $n(r)$ is the number of particles situated from a reference particle at a distance of $r$ and $r + \Delta r$) for various pairs, such as ion-ion and ion-neutral, is computed at equilibrium condition for both potentials. The curves are shown for $\chi_i = 10^{-1}$ in Figure 6. For the ion-ion pair (cf. Figure 6(a)), the $g(r)$ curves for both potentials overlap each other. The values of $g(r)$ are negligible for $r/a_{ii} \le 0.54$. These curves attain a peak ($\sim 1.1$) at $r/a_{ii} = 1.65$ and then remain constant with the value of 1 throughout. The peak height decreases and almost becomes 1 for lower values of $\chi_i$ (not shown). It signifies that the ions do not form any structures in the system [16], and the electron screening does not impact the internal structure of the ions. For the ion-neutral pair (cf. Figure 6(b)), $g(r)$ is negligible for $r < r_\phi$. $r_\phi$ is the radius

at which the repulsive core in the ion-neutral interaction potential acts (cf. eqn. (4)). $g(r)$ for the ion-neutral pair shows multiple peaks for both potentials and attains the primary peak at $r/r_\phi = 1.23$ and $1.10$ for Coulomb and Yukawa potentials, respectively. It is observed that the peak widths are smaller for the Yukawa potential. Furthermore, the multiple secondary peaks indicate the formation of possible shell structures among the ions and neutrals [32]. The position and width of the peaks dictate the ion-neutral interaction energy. The primary peak for the Coulomb potential falls in the region where the potential energy for the ion-neutral interactions is negative with a higher magnitude. A more detailed discussion of the potential energy of an ion due to the ion-neutral interaction can be found in Section IV and Appendix A.
For lower values of $\chi_i$, the primary peak height increases for both potentials, with distinctive secondary peaks (not shown).

## IV. DISCUSSION

To understand the potential energy contribution to the DIH and oscillations in the $T_i$ curve, a model is developed. According to the DIH hypothesis, the increase in ion kinetic energy is attributed to the decrease in their potential energy. Mathematically, for an individual ion, it can be written as

$$ \Delta(k_B T_i) = -(\Delta PE_T), \quad (8) $$

where $PE_T$ is the total pairwise potential energy experienced by an ion due to other ions and neutral particles. Therefore, $PE_T = PE_I + PE_N$, where $PE_I$ is the potential energy due to other ions which can be either due to Coulomb potential ($PE_C$) or Yukawa potential ($PE_Y$) and $PE_N$ is the ion potential energy due to the neutrals.

For two arbitrary initial ($t$) and final ($t+\Delta t = t_f$) time instances, eqn. 8 can be written as,

$$ \begin{aligned} \Delta T &= \underbrace{T_i(t + \Delta t) - T_i(t)}_{\Delta T_i} \\ &= \frac{1}{k_B}((PE_I(t) - PE_I(t + \Delta t)) + (PE_N(t) - PE_N(t + \Delta t))), \quad (9) \end{aligned} $$

where $\Delta T$ is the change in temperature. For Coulomb potential, the expression for $\Delta T$ becomes,

$$ = \underbrace{\frac{q^2}{4\pi\epsilon_0 k_B} \left[ \frac{1}{r(t)} - \frac{1}{r(t + \Delta t)} \right]}_{-\Delta PE_I/k_B} + \underbrace{\frac{q^2 \alpha_R a_B^3}{8\pi\epsilon_0} \left[ \left( \frac{r_\phi^8}{r^{12}(t)} - \frac{1}{r^4(t)} \right) - \left( \frac{r_\phi^8}{r^{12}(t + \Delta t)} - \frac{1}{r^4(t + \Delta t)} \right) \right]}_{-\Delta PE_N/k_B}. \quad (10) $$

Similarly, for Yukawa potential, the expression for $\Delta T$ becomes,

$$ = \underbrace{\frac{q^2}{4\pi\epsilon_0 k_B} \left[ \frac{e^{\kappa r(t)}}{r(t)} - \frac{e^{\kappa r(t+\Delta t)}}{r(t+\Delta t)} \right]}_{-\Delta PE_I/k_B} + \underbrace{\frac{q^2 \alpha_R a_B^3}{8\pi\epsilon_0} \left[ \left( \frac{r_\phi^8}{r^{12}(t)} - \frac{1}{r^4(t)} \right) - \left( \frac{r_\phi^8}{r^{12}(t+\Delta t)} - \frac{1}{r^4(t+\Delta t)} \right) \right]}_{-\Delta PE_N/k_B}. \quad (11) $$

$\Delta T_i$ (cf. eqn. (9)) is the change in ion temperature in the time interval $\Delta t$ and is acquired from the simulations directly, where $\Delta t$ is chosen as 0.05 ps. $\Delta T$ due to the change in potential energies (R.H.S of eqns. (10) and (11)) is acquired from the simulation by taking the collective potential energies of all the ions and then averaging their values by dividing with $N_i$. $\Delta T$ is plotted against $t_f$ in Figure 7 for both potentials for $\chi_i = 10^{-1}$, showing the contribution of different pairwise interactions to $\Delta T$. It is observed that for Coulomb potential, the initial rise of $T_i$ due to DIH is attributed to the $PE_I$ contribution. As the ions start to cool down, the contribution of $PE_I$ diminishes, and $PE_N$ increases. The rapid decrease in $PE_I$ with time is caused by the increased distance between the ions to attain a minimum potential energy state among themselves. For Yukawa potential, $\Delta T_i$ is primarily due to the $PE_I$ contribution and the $PE_N$ contribution is negligible throughout. Furthermore, it is found that, for Coulomb potential, the values of $PE_N$ are negative and large, whereas, for Yukawa potential, these are negative but negligibly small. Upon analyzing the ion-neutral interaction potential, it is understood that the ion-neutral interaction energy is negative and significant only for a small range of separation distance for an ion-neutral pair ($r/r_\phi \sim 1-2.4$) (cf. Figure A.1). For $r/r_\phi > 2.4$, the values of the interaction potential are negative but negligible. This shows that apart from the separation distance of the ion-neutral pair, the number of neutrals and their arrangement around an ion also play a crucial role in determining the values of $PE_N$.

The damped oscillations in the $\Delta T$ curve at long timescales (> 2 ps) are only observed for Coulomb potential. The frequency and the amplitude of these oscillations decrease with a decrease in $\chi_i$ due to lower ion densities (not shown in the graph). Furthermore, for the Yukawa potential, as $\kappa$ decreases and approaches the value 0, the oscillations appear and are comparable to those shown for the Coulomb potential (not shown graphically).

## V. CONCLUSION

Atmospheric pressure He plasma at room temperature is modeled for $\chi_i = 10^{-1} - 10^{-5}$ employing molecular dynamics simulations, where two potentials are used to describe the ion-ion interaction; Coulomb and Yukawa. Yukawa potential incorporates the screening effect caused by the background electrons.

It is found that the DIH, which is prominent for $\chi_i \ge 10^{-3}$, is more significant for the Coulomb potential. The electron screening in the Yukawa potential efficiently limits the ion and gas heating in the system. The maximum ion temperatures for the weakly ionized plasmas ($\chi_i \le 10^{-4}$) are almost similar for both potentials. During the non-equilibrium phase ($T_i \ne T_g$) caused by the DIH, the $T_i$ curve for the Coulomb potential shows an oscillatory behavior for $\chi_i \ge 10^{-3}$. The oscillations are
very much damped for the Yukawa potential. For weakly ionized plasmas ($\chi_i \le 10^{-4}$), such oscillations are absent for both potentials. During this non-equilibrium phase, the ions show a sub-diffusive behavior for $\chi_i \ge 10^{-3}$, which is more prominent for the Coulomb potential. The sub-diffusive behavior is exclusively associated with the DIH mechanism.

The decrease in the potential energy due to the ion-ion interaction contributes to the increase in $T_i$ during the DIH process for both Coulomb and Yukawa potential cases. The contribution to the change in ion temperature by the ion-neutral interactions is dominant for Coulomb potential.

In general, the potential energy due to the ion-neutral interactions is attractive for both potentials and the magnitude is higher for Coulomb potential. Combined with the $g(r)$ data for the ion-neutral pair, it is understood that a higher number of neutrals are arranged near an ion for Coulomb potential, and the ion-neutrals form a shell-like structure for both potentials for all the values of $\chi_i$.

To conclude, it is more scientific to include the screening effect in the simulations for the atmospheric pressure plasma systems, as this effect dictates the plasma dynamics, energetics, and internal structure more realistically. The mean potential energy in a Yukawa system is lower, which is energetically favourable for the stability

of the system. Although the Coulomb potential is easier to implement with long-range effects and is computationally less expensive, the Yukawa potential may be important for ascertaining more physical results for the given experimental conditions. However, the Yukawa potential needs to be investigated further to ascertain how it guides the ion-neutral structural arrangements at room-temperature plasmas and to extend the current finding to high and low gas temperatures at atmospheric pressure, which may be taken up in the future works.

## VI. ACKNOWLEDGEMENT

The authors (S.S. Mishra and S. Bhattacharjee) sincerely thank the HPC facility at IIT Kanpur and the funding agency DST-SERB India for the research grant CRG/2022/000112.

## Appendix A: Ion-neutral interaction potential

Figure A.1 shows the values of the interaction potential ($\phi_{ind}$) and force ($F_{ind}$) for the ion-neutral interactions for $\chi_i = 10^{-1}$. For $r/r_\phi \le 1.14$, the values for $F_{ind}$ are positive (repulsive). $\phi_{ind}$ attains minimum at $r/r_\phi = 1.14$, and the values of $F_{ind}$ becomes negative (attractive).

[1] K. Barman, S. Luhar, R. Rane, D. N. Srivastava, S. K. Nema, and S. Bhattacharjee, Improving electrochemical sensitivity of screen-printed carbon electrodes by atmospheric pressure plasma jet treatment and electrochemical detection of dopamine, Plasma Processes and Polymers **20**, 2200161 (2023).

[2] D. Behmani and S. Bhattacharjee, Plasma parameters and electric field fluctuations in a cold micro-plasma jet interacting with a substrate, Physics of Plasmas **30** (2023).

[3] D. Behmani, M. Keidar, and S. Bhattacharjee, Frequency-locked measurement of floating potential fluctuations in cold micro-plasma jets: Interactive surface feedback effects on reactive species generation, Physics of Plasmas **31** (2024).

[4] E. C. Neyts and P. Brault, Molecular dynamics simulations for plasma-surface interactions, Plasma Processes and Polymers **14**, 1600145 (2017).

[5] T. Bernhardt, M. L. Semmler, M. Schäfer, S. Bekeschus, S. Emmert, and L. Boeckmann, Plasma medicine: Applications of cold atmospheric pressure plasma in dermatology, Oxidative medicine and cellular longevity **2019**, 3873928 (2019).

[6] M. Keidar, D. Yan, I. I. Beilis, B. Trink, and J. H. Sherman, Plasmas for treating cancer: Opportunities for adaptive and self-adaptive approaches, Trends in biotechnology **36**, 586 (2018).

[7] S. Hofmann, A. Van Gessel, T. Verreycken, and P. Bruggeman, Power dissipation, gas temperatures and electron densities of cold atmospheric pressure helium

<figure><img src="image_9.png" /><figcaption>FIG.A. 1. Ion neutral interaction potential and force values for ionization fraction, $\chi_i = 10^{-1}$.</figcaption></figure>

and argon rf plasma jets, Plasma Sources Science and Technology **20**, 065010 (2011).

[8] Z. Ouyang, V. Surla, T. S. Cho, and D. N. Ruzic, Characterization of an atmospheric-pressure helium plasma generated by 2.45-ghz microwave power, IEEE Transactions on Plasma Science **40**, 3476 (2012).
[9] Q. Wang, D. J. Economou, and V. M. Donnelly, Simulation of a direct current microplasma discharge in helium at atmospheric pressure, Journal of Applied Physics **100** (2006).

[10] R. Van der Horst, T. Verreycken, E. Van Veldhuizen, and P. Bruggeman, Time-resolved optical emission spectroscopy of nanosecond pulsed discharges in atmospheric-pressure n2 and n2/h2o mixtures, Journal of Physics D: Applied Physics **45**, 345201 (2012).

[11] K. Barman, D. Behmani, M. Mudgal, S. Bhattacharjee, R. Rane, and S. K. Nema, Characteristics of atmospheric pressure micro-plasma jets in two different modes of excitation depending upon wave amplitude and frequency, Plasma Research Express **2**, 025007 (2020).

[12] M. D. Acciarri, C. Moore, and S. D. Baalrud, Strong coulomb coupling influences ion and neutral temperatures in atmospheric pressure plasmas, Plasma Sources Science and Technology **31**, 125005 (2022).

[13] P. S. Kothnur, X. Yuan, and L. L. Raja, Structure of direct-current microdischarge plasmas in helium, Applied physics letters **82**, 529 (2003).

[14] M. J. Kushner, Modelling of microdischarge devices: plasma and gas dynamics, Journal of Physics D: Applied Physics **38**, 1633 (2005).

[15] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. In't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, et al., Lammmps-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Computer Physics Communications **271**, 108171 (2022).

[16] S. S. Mishra and S. Bhattacharjee, Temperature guided behavioral transitions in confined helium: Gas-wall interaction effects on dynamics and transport in the cryogenic limit, Chemical Thermodynamics and Thermal Analysis **7**, 100073 (2022).

[17] S. S. Mishra and S. Bhattacharjee, Collisional cooling in confined gases at cryogenic temperatures: Time dependent near-wall accumulation and gas dynamics, International Journal of Heat and Mass Transfer **220**, 125008 (2024).

[18] K. Ramkumar, S. S. Mishra, and S. Bhattacharjee, Uncertainties in the transport properties of helium gas at cryogenic temperatures determined using molecular dynamics simulation, Fluid Phase Equilibria **577**, 113979 (2024).

[19] Y. V. Arkhipov, A. Askaruly, A. E. Davletov, D. Y. Dubovtsev, Z. Donkó, P. Hartmann, I. Korolov, L. Conde, and I. M. Tkachenko, Direct determination of dynamic properties of coulomb and yukawa classical one-component plasmas, Physical Review Letters **119**, 045001 (2017).

[20] G. Dharuman, L. G. Stanton, and M. S. Murillo, Controlable non-ideal plasmas from photoionized compressed gases, New Journal of Physics **20**, 103010 (2018).

[21] M. Acciarri, C. Moore, and S. Baalrud, Influence of strong coulomb coupling on diffusion in atmospheric pressure plasmas, Plasma Sources Science and Technology **32**, 115004 (2023).

[22] M. Lyon, S. D. Bergeson, and M. S. Murillo, Limit of strong ion coupling due to electron shielding, Physical Review E **87**, 033101 (2013).

[23] J. A. Bittencourt, *Fundamentals of plasma physics* (Springer Science & Business Media, 2013).

[24] F. F. Chen et al., *Introduction to plasma physics and controlled fusion*, Vol. 1 (Springer, 1984).

[25] M. Puchalski, K. Szalewicz, M. Lesiuk, and B. Jeziorski, Qed calculation of the dipole polarizability of helium atom, Physical Review A **101**, 022505 (2020).

[26] D. J. Evans and B. L. Holian, The nose-hoover thermostat, The Journal of chemical physics **83**, 4069 (1985).

[27] J. Castro, P. McQuillen, H. Gao, and T. Killian, The role of collisions and strong coupling in ultracold plasmas, in *Journal of Physics: Conference Series*, Vol. 194 (IOP Publishing, 2009) p. 012065.

[28] S. D. Bergeson, A. Denning, M. Lyon, and F. Robicheaux, Density and temperature scaling of disorder-induced heating in ultracold plasmas, Physical Review A **83**, 023409 (2011).

[29] M.-B. Luo and D.-Y. Hua, Simulation study on the mechanism of intermediate subdiffusion of diffusive particles in crowded systems, ACS omega **8**, 34188 (2023).

[30] V. Vishnyakov and G. Dragan, Coupling parameter for low-temperature plasma with condensed phase, Condensed Matter Physics **10**, 201 (2007).

[31] T. Ott, M. Bonitz, L. Stanton, and M. S. Murillo, Coupling strength in coulomb and yukawa one-component plasmas, Physics of Plasmas **21** (2014).

[32] M. Yadav, P. Deshwal, S. Maity, and A. Das, Structure formation by electrostatic interactions in strongly coupled medium, Physical Review E **107**, 055214 (2023).