Structure formation by electrostatic interactions in strongly coupled medium

Electrostatic interactions play an overwhelmingly important role and form the basis of many processes, involving plasmas Tsytovich (2007); Filippov et al. (2015); Ratynskaia et al. (2006), soft matter Naji et al. (2005), biological medium Holm et al. (2001); Levin (2002), etc. In these cases, the system is often modelled by variants of Langevin and Fokker Planck equations involving electrostatic interactions (e.g. the Poisson Nernst Planck equations Bazant et al. (2004); Kilic et al. (2007)). Formations of colloidal suspensions, polymers, and clusters in these contexts have always attracted attention and often have significant technological and scientific relevance. Our attempt here is to understand the process of structure formation with the help of Molecular Dynamics simulations in the simplest scenario of a collection of uniformly distributed charged particles of both signs, interacting via long-range Coulomb force along with Lennard Jones (LJ) interaction. The LJ potential helps in switching off the attractive Coulomb interaction amidst oppositely charged particles at short distances preventing their collapse. One is thus seeking the formation of clusters in this medium, wherein two or more particles join together to form a stable bound configuration. It is clear that the thermal effects will let the particles fly away inhibiting the survival of any such binding configuration. However, if the medium is in a strongly coupled state for which the mutual interaction energy exceeds the thermal energy, the possibility of the formation of such a bound state exists. We, therefore, explore the dynamics in the strongly coupled regime for this medium.

We have chosen a collection of equal and oppositely charged particles interacting with electrostatic Coulomb and LJ potential in the strongly coupled regime for our model. This is essentially like a quasi-neutral strongly coupled plasma state of matter Fortov et al. (2006); Ichimaru (1982). Preparation of a strongly coupled plasma medium is usually a challenge as it is typically created by doing violence to matter and hence is seething with thermal energy. The strongly coupled state requires that the parameter $\Gamma=Q^{2}/aK_{b}T>1$, where $Q$ is the charge state of the particle, $a$ is the inter-particle separation and $T$ is the temperature and $K_{b}$ is the Boltzmann constant. Thus the strong coupling condition can be achieved by (i) enhancing the charge $Q$, (ii) reducing the inter-particle distance $a$, and (iii) decreasing the temperature $T$.

Working with high $Q$ particles has led to the observation of strong coupling effects even at room temperatures in dusty plasmas Chu and Lin (1994); Thomas et al. (1994) in the form of crystallization of highly charged micron-sized dust particles against a background of normal electron-ion plasma. The strongly coupled dusty plasma medium has garnered a lot of attention in the plasma community for a couple of decades now. The strong coupling effects playing a distinct role in the collective modes have been studied extensively both theoretically Kaw and Sen (1998) and experimentally Murillo (2004). For the strongly coupled dusty plasma medium properties like collective structures Kumar et al. (2017); Das et al. (2014), linear and nonlinear waves Kumar et al. (2018), instabilities Tiwari et al. (2012); Das and Kaw (2014), phase transitions Maity and Das (2019); Maity and Arora (2022) have been observed. The interplay of single particle and collective dusty plasma dynamics has also been studied Maity et al. (2018) The dynamics of small dust clusters have illustrated dynamical equilibrium states with chaotic motion Deshwal et al. (2022); Maity et al. (2020).

In super-dense stars Van Horn (1991), the interior of planets, inertial confinement plasmas, etc., the reduction of inter-particle separation $a$ can lead to strong coupling behavior. In some of these contexts, however, the dynamics might involve electromagnetic effects. There has also been tremendous progress in the creation of Ultracold plasma experimentally using laser-cooling techniques. For instance, xenon atoms via a photo-ionization process are confined by a magneto-optical trap Killian et al. (1999); Kulin et al. (2000); Killian et al. (2001) with typical density in the range of $10^{5}cm^{-3}$ to $10^{10}cm^{-3}$ with either both or one of the species in a strongly coupled state. Many properties of ultracold plasma have been studied like the formation of Rydberg atoms Killian et al. (2001). Any heating in such a system needs to be avoided to keep the system in a strongly coupled state. For this, studies have been conducted to identify and restrict the inherent heating mechanisms that are prevalent in these systems. Mechanisms such as disorder-induced heating Cummings et al. (2005); Guo et al. (2010) and three-body recombination(3BR) Kuzmin and O’Neil (2002) have been identified as some prominent factors causing the heating of the medium. A study of thermodynamics and transport properties such as pressure and internal energy Tiwari et al. (2017) etc., has also attracted attention. It is also interesting to note that studies in magnetically confined plasma Gilbert et al. (1988) have shown ions in the strong coupling regime arranged in concentric shells and demonstrate solid and liquid-like behavior. In another work by T. Pohl et. al., Pohl et al. (2004) the crystallized structure was observed at the center of expanding laser-cooled ultracold plasma in which short-range concentric ion shells are formed with appropriate initial conditions.

Our objective here is to explore the possibility and nature of structure formation in a neutral plasma. The study also has relevance in other contexts of soft matter, and biological systems in which structure formation is an important issue, and electrostatic interactions are believed to be responsible for the same. Here, one is seeking the formation of structures with the most elementary coulomb interaction with Lennard Jones potential at short distances to avoid the collapse of unlike charged particles. We make use of the open-source molecular dynamics code of LAMMPS for carrying out these studies numerically.

This article has been organized as follows. In Section II, we discuss the simulation setup and the choice of parameters. In section III, we investigate the formation of classical bound states. In sections IV and V, we study the effect of an external perturbation introduced in the 2-dimensional and 3-dimensional simulation systems. The perturbation is in the form of inserting a massive and heavily charged particle into the medium. This could be equivalent to applying a biased probe in experiments. This inserted particle remains fixed in space while the particles in the medium respond to its potential. The charge in it is chosen to be quite high, about two orders higher than the charge of other particles in the system. The evolution of the system is studied in both two and three dimensions. Oppositely charged particles assemble around the externally placed charge particle to screen its potential. However, when the underlying plasma is strongly coupled, the oppositely charged particles assemble in a crystalline pattern around it. The behavior of the crystalline pattern is understood with the help of the Radial distribution function and Voronoi diagrams. As expected the potential due to this crystalline charge cloud typically resembles the Debye shielding profile. However, there are some deviations, arising due to the crystalline and discrete nature of the shielding cloud. The value of $\sigma$, which typically defines the distance where the LJ repulsive potential goes to zero, is also found to play a role. The role of external perturbation on the formation and existence of the seemingly fragile classical bound state is studied. It is observed that the disturbance induced by the insertion of the external highly charged particle leads to the formation of higher numbers of bigger-sized bound state clusters. This happens as the particles placed far apart now get an opportunity for close encounters. However, in all these simulations no three-dimensional bound state structure formation was observed. In section VI it is shown that when some two-dimensional bound states are evolved in isolation (absence of other plasma particles) they reorganize as three-dimensional bound structures. In fact, it is shown that the 3-d structures that form are energetically favorable. Yet such structures do not form in the backdrop of plasma medium. It clearly shows that the plasma medium plays an important role in the formation of the clusters. In Section VII we summarize and conclude our studies.

We have carried out 2D and 3D Molecular Dynamics simulations using an open-source classical MD simulation software LAMMPS Plimpton (1995). We have typically used the parameters associated with electron-ion plasma in an ultra-cold regime. The particles in the simulation box is interacting with each other through the long-range pair Coulomb potential ($V_{pC}$) and an additional short-range repulsive LJ potential ($V_{LJ}$) as given by

$$V(r_{kl})=V_{pC}+V_{LJ}$$

(1)

where,

$${V}_{pC}=\frac{Q_{k}Q_{l}}{4\pi\epsilon_{0}r_{kl}}$$

(2)

$${V}_{LJ}=4\epsilon\left[\left(\frac{\sigma}{r_{kl}}\right)^{12}-\left(\frac{\sigma}{r_{kl}}\right)^{6}\right]$$

(3)

This short range LJ potential avoids the recombination of electrons and ions. Here $Q_{k}$ and $Q_{l}$ define the $k^{th}$ and $l^{th}$ particles respectively. Also, $r_{kl}$ represents the distance between these two particles. Here, $\epsilon$ and $\sigma$ are the usual Lennard-Jones potential parameters where $\epsilon$ defines the strength of the LJ potential well and $\sigma$ describes the distance at which inter-particle LJ potential becomes zero. The PPPM (particle-particle particle-mesh) method Heath (1997) is used for the calculation of long-range Coulomb interactions. In our case, we have assumed both electrons and ions to be in a strongly coupled regime which means the coulomb coupling parameter ($\Gamma$) for both species is greater than unity.

In the two-dimensional (2D) study, we have chosen a rectangular simulation box of length $L_{x}=L_{y}=1.4\times 10^{-3}$m in the $\hat{x}$ and $\hat{y}$ directions, respectively. Periodic boundary conditions are considered in all directions. Initially, $10^{4}$ electrons and $10^{4}$ ions with densities $1.0\times 10^{10}$ m^-2 are distributed randomly inside the simulation box. This ensures that there exist no initial spatial pre-correlations between the charged particles. The charge $(Q_{i})$ and mass $(m_{i})$ of the ion are taken to be $1.6\times 10^{-19}$ C and $100m_{e}$, respectively. Here, $m_{e}$ is the mass of an electron, i.e., $m_{e}=9.109\times 10^{-31}$Kg. The mass of the ion in simulation is chosen smaller compared to even the proton mass, merely to reduce the simulation time. In our simulations, the temperature of electrons and ions are both considered to be equal and at $T_{e}=T_{i}=0.1$ K. All the length scales are normalized by the average inter-particle distance $a$. In 2D, the average inter-particle distance $a$ and the areal density $n$ are related as $a=(n\pi)^{-1/2}$. Thus for our choice of $n$, the value of average inter-particle distance turns out to be $a=5.64\times 10^{-6}$m. In our simulations, the time scales are normalized by $\omega_{pe}^{-1}$, where $\omega_{pe}=(ne^{2}/m_{e}\epsilon_{o})^{1/2}$ is electron-plasma frequency.

In our 3D simulations, we have taken the length of the simulation box to be $L_{x}=L_{y}=L_{z}=3.7\times 10^{-5}$ $m$, in $\hat{x}$, $\hat{y}$, and $\hat{z}$ directions, respectively. Again, $10^{4}$ electrons and $10^{4}$ ions with densities $3.7\times 10^{17}$ m^-3 are distributed randomly inside the simulation box. The relation between density $n$ and average inter-particle distance is $a=(3/(4n\pi))^{1/3}$. The temperature of electrons and ions is chosen to be $T_{e}=T_{i}=1$ K, which corresponds to the Coulomb coupling parameter $\Gamma=20$ i.e. in the strongly coupled regime.

As mentioned earlier, particles (both, electrons and ions) are interacting with long-range pair Coulomb potential. Additionally, a short-range LJ pair potential is present in the pair interactions between particles. The parameters associated with these interaction potentials are as follows. The cutoff distance for Coulomb pair potential is chosen to be $r_{c}=20a$. The parameter $\sigma$ associated with the LJ potential is chosen to be $10^{-2}a$ for most of the simulation studies. We have, however, also varied the value of $\sigma$. The normalized parameter $\epsilon$ is chosen to be $558K_{B}T$ and $75K_{B}T$ for 2D and 3D, respectively.

In both 2D and 3D simulations, Nose-Hoover thermostats Nosé (1984); Hoover (1985) have been used to achieve the thermal equilibrium associated with the desired value of Coulomb coupling parameter $\Gamma$.

Once the system achieves thermodynamic equilibrium with the attached thermostat, its evolution shows the formation of several distinct varieties of bound clusters. In these structures, a few electrons and ions combine with each other and form a stable configuration. Some such structures have been displayed in Fig.1. There are linear structures as shown in the subplot (a) and (b) of this figure. In fact in subplot (b) about 9 particles (4 electrons and 5 ions) have joined together to form a linear chain. Such linear chains have alternating particles with opposite signs attached together. Subplot(c) shows a 2-D structure in which two linear chains join adjacently to form a 2-d structure. In subplot(d) we have shown the formation of a square structure in which the diagonally opposite particles are of the same sign. Two such square structures can join together as shown in subplot(e) or a square can have a linear chain attached to itself as shown in subplot(f) resembling a kite. Subplot (g) is again a square structure with $9$ particles, each side of the square having $3$ particles. The circular ringed structure in subplot (h) forms when the particles at the end of the linear chain join together. Subplot (i) is again a kite-like structure like that shown in (f). It should be noted that the bound structures that form do not necessarily have an equal number of ions and electrons. A mismatch of one unit charge (but not higher than this) has been observed.

We illustrate the dynamical process of formation of some of these structures in Fig.(2) by showing the snapshot of a particular zoomed region from the entire simulation region which has been shown in the subplot(a) of the figure at $t=0$, where the randomly distributed particles have been shown. From subplot (b) to (t) we show the snapshots at various times for a zoomed region from the simulation box containing $6$ particles. The entire dynamics is covered within a fraction of the plasma period. During this period no other particle entered this region hence the particle number remained at $6$. From subplot (b) to (f) one can observe the attachment of an electron to an ion-electron pair placed at the top of the box. A $3$ particle linear chain can be observed to be rotating at the bottom right corner of the subplot. From subplot(g) to (k) one can observe how the two $3$ particle linear chains approach each other and join together. After combining, they wiggle a lot going through a variety of phases shown in subplots (l) to (p). Finally, as they relax they form a ringed structure which is found to be stable.

We would like to emphasize here that, while the simulations have been carried out both in 2-D and 3-D, we observe no formation of three-dimensional clusters. In the next section, we observe what happens to these structures when the medium is perturbed by inserting a highly charged heavy particle having no dynamics of its own. We also study how the plasma particles accumulated around it to shield the potential of this externally placed particle in the medium.

When a weakly coupled plasma is disturbed by an external potential, the plasma tries to shield this potential within a distance of Debye length. It is, however, not clear how a strongly coupled plasma would react to an externally applied field. In this section, we try to explore this particular issue in a two-dimensional simulation set-up. Once the system gets settled in thermal equilibrium an external point charge is inserted at the center of the simulation box. The mass ($M_{P}$) and charge ($Q_{P}$) of external perturbation are chosen to be large as $100m_{i}$ and $100Q_{i}$, respectively. Since the externally inserted charge has a very high mass it remains stationary and is not allowed to move. The electrons get attracted towards this externally applied positive charge whilst the ions are repelled from it. The insertion thus initiates a dynamical response from the plasma. We investigate the behavior of the shielding cloud that accumulates around this charge and also the influence of such a dynamic response from the plasma on the classical bound pairs. It is observed that the shielding occurs in the form of an electron cloud arranged in a crystalline pattern due to the strongly coupled nature. It is interesting to note that despite the plasma getting perturbed, even bigger-sized bound state clusters form in the bulk region. We discuss these two effects in the subsequent subsections.

We have also carried out our simulation studies in three dimensions. We distribute $10^{4}$ electrons and $10^{4}$ ions randomly in the simulation box and seek the attainment of thermodynamic equilibrium at the desired temperature. Thereafter, we perturb the system by inserting a massive particle having a charge of $100Q_{i}$ at the center of the simulation box. Since the inserted particle has a positive charge again the electrons in the system rush to shield its charge. This time, however, the system is three-dimensional so the electrons form a three-dimensional shielding cloud. We now investigate whether the shielding cloud has a crystalline form and if so what kind of the arrangement takes place in the cloud. As the structure cannot be visualised in 3-D we have employed certain diagnostics to perceive it in the best possible manner.

In the simulations that we have carried out so far, we have observed only 2-D bound state formation. The bound state shows considerable stability and survives for a long time in the plasma where it faces constant bombardments from other particles. In this section, we explore the stability of a couple of chosen bound structures in isolation by tracking their evolution.

We first study the stability of the ring structure shown in subplot (a) of Fig(16) as an initial configuration and track its evolution. It quickly (within a fraction of the plasma period) attains the form shown in subplot(b). We also placed a pair of rings on top of each other (subplot(c) ) and a linked pair of rings (subplot(e). Both these configurations relax immediately to a cuboid form as shown in subplots (d) and (f) of Fig(16). It appears as a puzzle that the ring structure which was stable in the presence of background plasma and survive for several plasma periods, is unstable and quickly relax to the cuboid structure in the absence of background plasma. It should also be noted that the final configuration is a 3-D bound state which was never observed in the plasma.

We now look at the energetics of the ringed and the cuboid structure. The interparticle spacings of the two-ringed and the cuboid structure are shown in Fig(17). The total potential energy of the ringed structure and the cuboid structure in terms of the interparticle spacings shown in Fig(17) has been evaluated.

$$PE_{pC}^{ring}=6\frac{Q_{i}Q_{j}}{4\pi\epsilon_{0}}\left(\frac{-3}{a}+\frac{2}{b}-\frac{1}{c}+\frac{2}{d}-\frac{2}{e}+\frac{1}{f}\right)$$

(4)

$$PE_{lj}^{ring}=6(2V_{lj}(a)+2V_{lj}(b)+V_{lj}(c)+2V_{lj}(d)+V_{lj}(e)+V_{lj}(f))$$

(5)

whereas, all the parameters in terms of $a$ are

$b=2asin(60^{\circ})$, $c=2a$

$d=\sqrt{2}a$, $e=2a$, $f=\sqrt{5}a$

The total potential energy of ringed structure is

$PE^{ring}(a)=PE_{pC}^{ring}+PE_{lj}^{ring}$

whereas,

$$PE_{pC}^{cuboid}=\frac{Q_{i}Q_{j}}{4\pi\epsilon_{0}}\left(\frac{-18}{a}+\frac{22}{b}-\frac{8}{c}-\frac{8}{d}+\frac{4}{e}\right)$$

(6)

$$\begin{split}PE_{lj}^{cuboid}=20V_{lj}(a)+4V_{lj}(2a)+22V_{lj}(b)+8V_{lj}(c)+\\ 8V_{lj}(d)+4V_{lj}(e)\end{split}$$

(7)

whereas, all the parameters in terms of $a$ are

The total potential energy for cuboid structure is

$PE^{cuboid}(a)=PE_{pC}^{cuboid}+PE_{lj}^{cuboid}$

We have then plotted the expression for the potential energy for the two cases as a function of $a/\sigma$ as shown in Fig.18. The green and magenta curve corresponds to overlapping ring and cuboid structure, respectively. The potential energy of the cuboid structure is found to be less than the overlapped ring structure. We believe this is the reason for favoring their formation in isolation. The presence of plasma clearly seems to modify this.

We have simulated an ultracold collection of an equal number of negative and positively charged particles (e.g. like electrons and ions, albeit their mass ratio is chosen to be different from the realistic case for ease of computation) in 2D and 3D by MD simulation in which both species are strongly coupled. These electrons and ions-like particles are interacting with long-range Coulomb potential, and a short-range LJ pair potential is present in the pair interaction between particles. At equilibrium, we have observed many types of bound structure formation. With the introduction of external perturbation (highly charged and massive point particle) we observe the expected phenomena of shielding. The particles are in a strongly coupled regime, therefore, the shielding cloud arranges in specific patterns. The shielding potential profile unlike the weakly coupled cases is determined by the choice of the LJ parameter of $\sigma$ defining the distance at which the interparticle potential is minimum. In the crystalline pattern, the inter-particle distance gets determined by this parameter. Thus when $\sigma$ is larger than the Debye length, the shielding is broader than expected from the Debye screening process. On the other hand, when it is smaller than the Debye length the screening is sharper. Thus one needs to be cautious while drawing inferences of certain plasma properties while simulating the Ultra Cold strongly coupled plasma using the short-range repulsive form to overcome the blowing up of the attractive potential amidst oppositely charged particles.

This research work was supported by a J.C. Bose fellowship grant (Grant No. JCB/2017/000055/SSC) from the Department of Science and Technology (DST), Government of India and the Core Research Grant (Grant No. CRG/2018/000624) of the science and Engineering Research Board (SERB), Government of India. The authors thank IIT Delhi HPC facility for computational resources. M.Y. is thankful to the University Grants Commission [Grant No. 1316/CSIR-UGC NET DEC.2017] for funding the research.