Self-excited converging shock structure in a complex plasma medium

The experiments related to this study have been performed in a newly built Dusty Plasma Experimental-II (DPEx-II) device. The extensive study on the experimental set-up and plasma diagnostics can be found in Ref. Arumugam et al. (2021). The device has an asymmetric electrode configuration, which has been used to generate DC plasma and large-sized dusty plasma. The schematic of the setup is shown in Fig.1(a). The uneven electrode configuration with a disc shape serves as an anode, and a wide rectangular tray with upward bents serves as a cathode. The DC glow discharge plasma is formed in an argon gas environment with a pressure of $7.0$ Pa and a DC voltage of 320-360 V. The plasma parameters, such as plasma density and electron temperature, are measured with the help of single and double Langmuir probes and are estimated to be 1-6 $\times~{}10^{15}~{}\rm/m^{3}$ and 2-4 eV, respectively. Mono-dispersive melamine formaldehyde dust particles of diameter 10.66 $\mu$m are introduced into the plasma to create a dusty plasma. The particles levitate above a stainless steel circular ring on the cathode’s right end. This ring also serves as a part of the cathode and creates a radial electrostatic trap due to the radial sheath around it. In the plasma environment, the dust particles get negatively charged as the larger number of electrons, being lighter than ions, attach to their surface. These negatively charged dust particles levitate at a certain height from the cathode by the balance of downward gravitational force and upward electric field of the cathode sheath. The radial expansion of the dust cloud is restricted by the radial sheath imposed by the ring called the confinement ring. Using the collision-enhanced model (CEC) Khrapak and Morfill (2006), the dust charge ($Z_{d}e$) is estimated to be $\sim~{}4.8\times 10^{3}~{}\rm e$. Our estimated charge from CEC theory has also been supported by the experimental verification based on a force balance model in our past study Arumugam et al. (2021) using the same experimental device and under similar discharge conditions. In experiments, the levitation of negatively charged dust particles occurs due to the balance of downward gravitational force ($m_{d}g$) and upward electrostatic force ($QE_{z}$). The mass density of melamine formaldehyde particle (dust particle) is 1.5 g/cm^-3, and the radius ($r_{d}$) is 5.33 $\rm\mu$m which gives the mass of the dust particle to be 9.5$\times~{}10^{-13}$ kg. The electric field strength measured using an emissive probe in the same device without the presence of dust particles is $E_{z}\sim$ 2$\times~{}10^{4}$ V/m Arumugam et al. (2021) for our discharge conditions. The force balance equation $m_{d}g=QE_{z}$ estimates the dust charge to be $Q\sim~{}3\times 10^{3}$e. Although the measurements from the probe are subjected to 20-30 $\%$ error, the estimated charge from the experimentally measured electric field is of the same order as from CEC theory. The 2D dust density is roughly calculated to be $\sim~{}1\times 10^{8}~{}/m^{2}$. The dust particles are illuminated by a horizontally aligned green laser light (XZ plane), and the dynamics are recorded from the top view high-speed camera.

To begin with, an argon pressure of 7 Pa and DC voltage of 320 V is set to create a plasma and dusty plasma. A higher dust density is aimed at introducing more particles in this device than in past experiments. A stable circular dust cloud with high dust density is created at this plasma condition. In this condition, a dust fluid with density $1\times 10^{8}~{}m^{-3}$ is obtained in a strongly coupled liquid state with an estimated coupling parameter $\Gamma\sim 180$. The coupling parameter $\Gamma=(Q^{2}/4\pi\epsilon_{0}ak_{B}T_{d})~{}exp(-a/\lambda_{D})$ is defined as the ratio of inter-particle electrostatic potential energy to the average thermal energy of the particles. Here, $Q,~{}\epsilon_{0},~{}a,~{}k_{B},~{}T_{d},and~{}\lambda_{D}$ are the particle’s charge, the electric permittivity of free space, inter-particle distance, Boltzmann constant, the temperature of the dust particles, and dust Debye length, respectively. The dust Debye length is defined as $\lambda_{D}=(\epsilon_{0}k_{B}T_{i}/n_{i}e^{2})^{1/2}$, where $T_{i},~{}n_{i}$ are the temperature of ions and plasma density, respectively. For our discharge conditions, i.e., pressure, $P=7$ Pa, and discharge voltage $V=320$ V, $T_{i}\sim 0.03~{}eV$, and $n_{i}\sim~{}1\times~{}10^{15}~{}\rm m^{-3}$ Arumugam et al. (2021), we estimated the value of the coupling parameter to be $\Gamma\sim 180$.

In our experiment, the sheath electric field of the horizontal confinement ring has a 3D profile. However, the monodispersive particles mainly form a 2D monolayer levitating at a certain height for our discharge conditions. The vertical component of the electric field, originating from both the cathode and circular ring, is balanced out by the gravitational force acting on the particles. At higher discharge voltages, the particles are compressed in the horizontal plane due to an increase in the horizontal confinement field. As a result, some particles may go down and form a quasi-2D layer. However, a particular monolayer shown in Fig.2 has the highest populated density and a much larger area in the horizontal plane. As a result, particles levitating in this particular layer are subjected to the stronger sheath electric field of the horizontal confinement ring. As the DC voltage is increased from 320 to 340 V, a single-density compressed structure is excited, propagating in the radially inward direction. The snapshots from the fast camera comprising the dust fluid and excited density compressed structure are shown in Fig.2. As time progresses, the density structure propagates radially inward, and the density compression becomes brighter and wider. The last snapshot shows the concentration of energy in the central region where the amplitude and width are higher. It is to be noticed that the luminosity is brighter at the lower edge than the upper edge, as can be seen from Fig.2. We are speculating that this non-uniformity comes from the non-uniform perturbation of the sheath electric field, which arises due to the manufacturing imperfection of the metallic ring responsible for the horizontal confinement or due to the inhomogeneous variation of plasma potential or electric field in the radial direction.

In order to obtain the quantitative parameters of the excited structure, such as amplitude, width, and velocity, we have estimated the radial density distribution of the dust cloud with respect to its center. The spatial profile of dust density is estimated from the scattered light intensity. The dust density is considered to be proportional to the intensity of light scattered from the dust particles, as the used CCD camera has a linear response to the light intensity. Now, to estimate the one-dimensional (1D) density profile as shown in Fig. 3(a), the whole dust cloud is radially binned by concentric circular strips with a width $dr=0.1$ mm. The 1D radial density profile is then obtained from the mean intensity of the scattered light within these circular bins. We define $n_{d}(r)/n_{d0}(r)\sim I(r)/I_{0}(r)$, where $I_{0}(r)$ is the intensity at a particular radial distance $r$ when the cloud is stable with no excited density compression, and $I(r)$ represents the same with the excited density structure, as shown in the images of Fig.2. The amplitude of the crest is defined as the peak of $n_{d}/n_{d0}$, and the width is estimated using the technique mentioned in Ref. Heinrich et al. (2009) by Heinrich et. al..

Fig.3(a) shows the evolution of the radial density profile of the dust cloud at different times. The shift of the peak towards the center ($r=0$) with time indicates the inward propagating nature of the self-excited density crest. Fig.3(b) shows the time evolution of the estimated amplitude and width of the excited density compression. Here, both the amplitude and width of the spontaneously excited density crest increase with time. The obtained amplitude values signify that the structure is nonlinear in nature. The nonlinearity can also be implied from the propagation speed. We have estimated the propagation speed of the density crest by tracking its evolution over time and found it to be $v\sim$ 3 cm/s. We have estimated the dust acoustic speed ($C_{d}$) using the expression, $C_{d}=\omega_{pd}\lambda_{D}$, which has a value $C_{d}\sim$$1.2$ cm/sec for our discharge conditions. Here, $\omega_{pd}=(Q^{2}/2\pi\epsilon_{0}m_{d}a^{3})^{1/2}\approx~{}262$ Hz and $\lambda_{D}\approx 45~{}\rm\mu m$ represent the dust-plasma frequency and dust Debye length, respectively. The estimated Mach number $(M=v/C_{d})$ $\sim$ 2.5 indeed shows that the structure is nonlinear in nature. The overall properties of the excited density pulse, e.g., radial shape, nonlinearity in the amplitude, and $M>1.0$, indicate that the excited density crest is nothing but a shock structure. The increase in width during its propagation is another typical characteristic of a shock structure due to the dissipative losses in the medium. Here, in our medium, dissipation can come from strong coupling effects as well as neutrals streaming. However, the increase in amplitude during its journey is interesting despite its propagation along the increasing dust density towards the center of the cloud. Tadsen et. al. Tadsen et al. (2017) experimentally showed the propagation characteristics of a linear DAW in an inhomogeneous dust cloud and found that the amplitude of linear waves increases with a decrease in dust density. Recently, Arora et. al. Arora et al. (2020) excited the shock waves and studied the propagation characteristics in a linear gradient of dust density. They found that the shock’s amplitude and width increase with a decrease in dust density. In the present study, the increase in the amplitude of the excited shock propagating towards the increasing dust density can be understood as follows. The shock waves have a high energy density distributed in a thin streak, and its converging nature signifies that the area is reducing; hence, energy density increases due to the conservation of energy followed by an acceleration of the shock wave. Thus, in our case, the rise in amplitude of the shock structure is attributed to its converging nature.

As the shock structure in our medium is spontaneously excited by increasing the discharge voltage, we have performed a series of experiments by changing the discharge voltage. Fig.4 shows the snapshots of dust clouds with the excited shock structure for four different discharge voltages at a constant time. Here also, an excitation of shock structure is seen at a discharge voltage of 350 V, similar to 340 V. The only difference is the position of the crest and compression factor. The structure travels more inward over a certain time at a higher voltage. This implies that the propagation speed increases with an increase in discharge voltage. However, if we keep increasing the discharge voltage, the circular symmetry breaks into segments. This is attributed to the non-uniformity in the amplitude of excited density crest at different azimuthal positions, as mentioned in the discussion related to Fig.2. At higher discharge voltages, the sheath electric field and, consequently, the perturbation strength increases. Thus, the scale of the non-uniformity over different azimuthal positions in the amplitude of the excited structure increases at higher voltages. The velocity of the excited structure is proportional to its amplitude. Thus, the non-uniformly excited structure will move with different velocities at different azimuthal locations, resulting in the breakup of the ring shape, forming blob-like structures.

The time evolution of the width and amplitude of the shock for two different discharge voltages are shown in Fig. 5. The subplots show that the amplitude and width of the shock increase as time progresses for each discharge voltage. However, at a particular instant of time, both the amplitude and width of the shock have higher values for higher discharge voltage. Since the exact cause of the excitation of the shock structure observed in our study is still unclear, we hypothesize that the abrupt change in the sheath electric field, triggered by the increasing discharge voltage, could potentially be the reason for the spontaneously excited perturbation. As the outermost particles within the dust cloud predominantly encounter the horizontal sheath electric field, the density perturbation is excited around the periphery of the dust cloud. With the increase in discharge voltage, the electric field strengthens, and the perturbation in the electric field increases. Hence, the density compression of the excited structure increases at higher voltages. We have also estimated the velocity of the inward propagating shock structure for different discharge voltages, as shown in Fig.6. It has been found that the structure excited in higher discharge voltage moves with higher velocity.

We have used an open-source massively parallel classical MD code LAMMPS Plimpton (1995) to perform three-dimensional MD simulations. Initially, two thousand point particles with charge $Q=4.8\times 10^{3}e$ and mass $m_{d}=9.5\times 10^{-13}$ kg are randomly distributed inside a 3D simulation box with length $L=0.05$ cm in all three directions. Here, $e$ represents the charge of an electron. Periodic boundary conditions in all three directions have been considered. Particles are assumed to interact via Yukawa pair potential. The electric field in the form of $\mathbf{E}_{x}(x)=-(m_{d}\omega_{r}^{2}/Q)(x-L/2)\hat{x}$ and $\mathbf{E}_{y}(y)=-(m_{d}\omega_{r}^{2}/Q)(y-L/2)\hat{y}$ has been applied in the x and y-directions, respectively, providing parabolic confinement of particles in the horizontal plane. In the vertical direction, we have considered the force due to gravity $m_{d}g(-\hat{z})$ (i.e., acting vertically downward) and force associated with an externally applied electric field $\mathbf{E}_{z}(z)=QE_{z0}\exp{(-\alpha z})\hat{z}$ (acting vertically upward for the negatively charged particle), mimicking the cathode sheath electric field in experiment Sickafoose et al. (2002); Arumugam et al. (2021). A combined effect of these two forces provides an effective parabolic confinement of particles in the vertical direction Maity and Das (2019). In our simulations, the values of $\omega_{r}$ and $\alpha$ are considered to be $1.0$ Hz and $500$ m^-1, respectively. The equation of motion of any $i$th particle can be expressed as,

where $r_{i}$ and $r_{j}$ are the positions of the $i$th and $j$th particles at a particular time, respectively, and $U(r_{i,j})=(Q/4\pi\epsilon_{0}|\bf{r}_{j}-\mathbf{r}_{i}|)\exp{(-|\mathbf{r}_{j}-\mathbf{r}_{i}|/\lambda_{D})}$ represents Yukawa pair potential between $i$th and $j$th particle. Here, $N$ represents the total number of particles.

Initially, particles were relaxed to a thermal equilibrium state with a desired temperature in the presence of a Nose-Hoover thermostat Nosé (1984); Hoover (1985). Later, the thermostat was disconnected, and the system was allowed to evolve under a micro-canonical ensemble where the system’s total energy remained conserved. In equilibrium, under the influence of all the confining forces, particles levitate at a certain height (i.e., at $z=L/2$ for our chosen values of $E_{z0}=(m_{d}g/Q)\exp{(\alpha L/2)}$), forming a monolayer disc in the horizontal (x-y) plane with radius $R$ with respect to the center of the dust cloud. Now, we introduce a perturbation to the system by applying electric fields in the following forms,

where $t_{0}$ is the time of the perturbation and $\delta$ represents the Dirac’s delta function. The electric perturbation has been applied only to the particles for which the radial locations with respect to the center of the monolayer, $r_{i}>0.9R$, i.e., particles residing at the outer radial locations. This results in an electrostatic force acting on the outer particles of the monolayer along the radially inward direction. As mentioned earlier, MD simulations have been performed only to qualitatively model our experimental observations concerning the shock characteristics, i.e., amplitude, width, velocity, and its converging nature. In simulations, the amplitude of the electric field perturbation is taken to be an arbitrary value so that the initial density perturbation (i.e., amplitude and profile) triggered in the monolayer remains in the same order as observed in the experiment. We have characterized the evolution of dust fluid under the influence of this electrostatic perturbation and presented them in the following subsection.

In equilibrium (without external perturbation), the dust cloud has a radially inhomogeneous density profile with a higher density near its center. The perturbation in the horizontal electric fields in the form of a pulse creates a density hump around the outer boundary of the 2D monolayer. Fig.7 (a)-(d) shows the distribution of particle number density in the x-y plane of the monolayer at times $\omega_{p}dt=8,~{}10,~{}12,~{}14$, respectively, for a particular case when the perturbation in the electric field is $E/E_{0}=0.8$. Here, $E_{0}=Q/4\pi\epsilon_{0}a^{2}\sim 193$ V/m is the average inter-particle electric field where $a=1/(n_{d0}\pi)^{1/2}$ represents the average inter-particle distance and $n_{d0}=N/\pi R^{2}$ is the average $2D$ density of the monolayer. The subplots of Fig.7 show that a single-density crest (red color) is created in the dust fluid in response to the electric perturbation. As time progresses, the density compression moves radially inward and becomes wider, similar to the observation in experiments. We have also estimated the radial density profile of the dust cloud under the influence of the external perturbation at different simulation times and presented them in Fig.8. The figure shows that the peak of density, representing the density crest, moves towards the center, ensuring the radially inward propagation of the crest. We have estimated the velocity of this structure and found that it propagates with a supersonic speed. It is also seen that the density structure has a triangular shape and its peak value increases during its propagation. Thus, the properties of the excited density structure, e.g., its shape, speed, and nonlinearity in the amplitude, ensure that it follows the same characteristics of a converging shock. We have extracted the quantitative parameters, such as amplitude and width, from the radial density profile and shown in Fig.9.

In Fig.9 (a) and (b), we have shown the time evolution of width as well as amplitude for two different electric field perturbations. The amplitude is defined as the maximum of $n_{d}/n_{d0}$, and the width is determined using the technique mentioned in Heinrich et al. (2009). The figure clearly shows that as time evolves, the amplitude of the excited density crest increases, and it becomes wider. Here, the width of the shock structure increases because of dispersion originating from the strong coupling effect of the medium. At the same time, the amplitude increase in its inward propagation is attributed to its converging nature and can be understood from the energy conservation point of view. As mentioned in section II.2, in our experiments, we increased the discharge voltage, which increased the sheath electric field, and found that the density crest is excited with a larger amplitude and width. In the simulation, we have replicated the phenomena by increasing the amplitude of electric field perturbation. Here also, at a given time, the amplitude and width of the excited shock pulse have higher values as compared to lower perturbation, and its amplitude increases as time progresses, as observed in our experiments.

The velocity variation estimated from the simulation runs with the increased electric perturbation has been shown in Fig.10. The figure depicts that velocity increases as we increase the amplitude of perturbation akin to Fig.6 obtained in the experiments. This can be explained as follows. As we increase the amplitude of perturbation, the amplitude of the density hump created at the boundary increases. A typical characteristic of a shock is that its propagation speed increases with an increase in nonlinearity. Since nonlinearity depends upon the amplitude of the excited structure, a density crest with higher amplitude propagates with higher velocity. Henceforth, all the results obtained in the simulation qualitatively explain the experimental observations.