Ionization waves in the PK-4 direct current neon discharge

The PK-4 system is a micro-gravity dusty plasma experiment currently in operation on-board the International Space Station. The experiment utilizes a long DC discharge in neon or argon gases. We apply our 2D particle-in-cell with Monte Carlo collisions discharge simulation to compute local plasma parameters that serve as input data for future dust dynamics models. The simulation includes electrons, Ne+ ions, and Ne m metastable atoms in neon gas and their collisions at solid surfaces including secondary electron emission and glass wall charging. On the time scale of the on-board optical imaging, the positive column appears stable and homogeneous. On the other hand, our simulations show that on microsecond time scales the positive column is highly inhomogeneous: ionization waves with phase velocities in the range between 500 m s−1 and 1200 m s−1 dominate the structure. In these waves, the electric field and charged particle densities can reach amplitudes up to 10 times of their average value. Our experiments on ground-based PK-4 replica systems fully support the numerical findings. In the experiment, the direction of the DC current can be alternated, which has been found to favor dust particle chain formation. We discuss possible mechanisms for how the highly oscillatory plasma environment contributes to the dust particle chain formation.


Introduction
The aim of this investigation is to contribute to the understanding of the collective dust particle chain formation observed in the recent PK-4 dusty plasma experiment. While approaching this goal we have developed a new, GPU accelerated cylindrical particle-in-cell simulation, built two ground-based replicas of the original micro-gravity experiment with different levels of simplifications, and gained detailed insight into the global and local operation characteristics of relevant direct current (DC) discharges. These results will serve as input to further models of dust charging and dust dynamics and also provide valuable information to the fundamental understanding of nonequilibrium and non-linear processes at multiple time scales in DC discharges.
In the next introductory sections a brief review of the PK-4 experiment and related results is provided, followed by summaries of the earlier studies on dust particle chains relevant to our system. We introduce our new numerical simulation in section 2, and our ground based PK-4 replica experiments in section 3. This is followed by the presentation and discussion of the results in section 4. A brief summary is given in section 5.

The PK-4 experiment
The PK-4 system is the latest generation in the line of micro-gravity dusty plasma experiments currently in operation on-board the International Space Station. The experiment utilizes a DC discharge in a glass tube with inner diameter of 30 mm and a length of ∼ 800 mm of which the central 400 mm is accessible for observations. The experiment is equipped with both neon and argon gases. The DC driving voltage can be modulated with arbitrary wave-forms and includes polarity switching up to 1 kHz repetition rate (with rise times of about 200 μs between the opposite-polarity DC voltages), an important tool for the prevention of a constant drift of the dust cloud due to electrophoresis. In addition to the DC drive two radio frequency (RF) coils, a heating ring, and a laser system are available for manipulation of the dust cloud in the positive column of the discharge. Various stable gas flows can be established and melamine-formaldehyde micro-particles with six different sizes between 1.31 and 10.41 μm diameters can be injected into the plasma. A detailed description of the setup, its capabilities and some results of the commissioning experiments can be found in [1].
Although the launch and commissioning of the apparatus took place in 2014 and 2015, the scientific program started years earlier utilizing different development versions of the device both on the ground and in micro-gravity via several parabolic flight campaigns since 2003 [2]. Langmuir probe characterization of the plasma is presented in [3], mapping the electron density as function of the gas pressure and discharge current in neon, resulting in values in the positive column in the range of 2V cm −1 < E < 3V cm −1 for the electric field, 2 × 10 14 m −3 < n e < 1 × 10 15 m −3 for the electron density, and 4 eV < ε < 6 eV for the mean electron energy. The capabilities of the 3D dust particle detection system are demonstrated for a liquid-like ensemble in [4].
The effect of the dust particle cloud on the emission spectra of the neon discharge in the PK-4 experiment is presented in [5,6]. As shown, spectral lines having 3p-levels shown an increase in their intensity in the presence of the dust cloud. A new kind of dusty plasma instability was reported in [7], where in the elongated cylindrical dust cloud the particles vibrated synchronously in the radial direction perpendicular to the tube axis attended by discharge glow fluctuations. The presence of dust density waves and their dependence on the DC driving voltage is presented in [8]. Bifurcation and remerging of new and old wave crests triggered by the polarity switching was identified. The radial velocity distribution of the dust particles was studied in [9] and a combination of Maxwellian + kappa distribution was identified. It has been reported in [10] that when applying polarity switching the dust particles show a strong tendency to form long chains oriented along the discharge axis. The reason for this phenomena is not yet understood, as all known processes that can induce large enough anisotropy of the electric field around a charged dust grains are driven by large, directed ion flow around the grains. In the following, a brief review of these processes and related findings are introduced.

Dust particle chains
In complex plasma, equally charged particles repel each other by the electrostatic Coulomb force. At large enough coupling (i.e., low temperature and high density) they can form ordered structures, where the distance between the particles is maximized, as in BCC or triangular lattices in three and two dimensions, respectively [11]. Since the early years of laboratory dusty plasma experiments, vertical alignment of dust particle pairs have been observed in the sheath above the lower electrode of an RF plasma chamber, apparently contradicting this simple principle. This phenomena is now attributed to the formation of an ion-wake field downstream of the dust particles [12].
Studies on the charging and coupling of a vertically aligned particle pair in the plasma sheath are presented in [13,14]. To enhance and control this phenomena, the horizontal confinement and local plasma environment can be modified by placing an open glass box onto the RF electrode, making the formation of long (on the order of 10 particles) vertical dust particle chains possible [15]. Numerical simulations of the charging, wake field formation and dust particle chain formation in plasma environments with subsonic and near supersonic ion flows is presented in [16,17]. Sophisticated theoretical results for the wake field potential are derived in [18]. Based on these investigations we can conclude that below the ion Mach number M = v ion drift /v ion sound < 0.2 no significant wake field formation is expected, which would still require v ion drift > 1000 m s −1 (for Ne + ions and 5 eV mean electron energy). This might be a realistic scenario in the sheath region of a capacitively coupled RF discharge or in the cathode fall of a DC discharge, but certainly not in the positive column of a DC discharge.
The first experimental observation and theoretical description of the electrorheological effect in dusty plasmas is given in [19]. This effect provides the possibility of tuning the inter-grain interaction with alternating external electric fields. Based on the early experiments obtained with the PK-4 setup, Ivlev et al derived a formula for the electric potential around a dust grain in plasma in an external electric field, where the anisotropy term is governed by the ion drift through the ion thermal Mach number [20]. This formula expanded to the second order (quadrupole-like) term takes on the form: where (r, θ) coordinates are the distance and polar angle (angle between the alternating electric field and the vector between the dust particles) and λ D is the Debye-length. We note that since the derivation of (1) several studies incorporating it into molecular dynamics simulations have failed [21,22]. However, more recently good agreement with PK-4 experimental results on string-fluid structures was demonstrated in [23] concluding that keeping only the quadrupolelike term and dropping the dipole-like term in (1) provides acceptable accuracy but only when M th was set to a much larger value than one would expect in dust free discharge. These results emphasize the non-linearity (in M th ) of the electrorheological effect. This has significant consequences as will be explained later. Recent experimental results for the dependence of dust particle chain formation on discharge parameters are given in [24]. Specifically, the dependence of the shape of the dust structures on discharge current and pressure in neon glow DC discharge at 77 K and 295 K has been studied in [25].
To summarize, both the ion wake effect and the electrorheological effect require large ion drift velocities as quantified by the sonic and thermal Mach numbers, respectively. However, the electric field in the steady positive column of a DC discharge is small (few V cm −1 ) resulting in ion drift velocities on the order of 100 m s −1 [26] and Mach numbers well below one. We aim to provide insight into the details of the plasma environment on multiple spatial and temporal scales, relevant to the PK-4 experiment, that will lead to a better understanding of the observed dust particle chain formation.

The PIC/MCC simulation
The primary tool of this numerical investigation is our newly implemented axi-symmetric two-dimensional electrostatic particle-in-cell with Monte Carlo collisions (2d3v PIC/MCC) discharge simulation. The principle of operation is analogous to that of other widely used 1d3v (one dimensional in real space and three-dimensional in velocity space) codes [27][28][29][30][31][32][33][34][35][36], but due to the higher spatial dimension it is able to include the effect of the glass cylinder that is confining the discharge. Unfortunately, at the same time, higher dimensionality results in the need for higher particle numbers and more sophisticated algorithms for solving Poisson's equation, demanding significantly more computations. Recent improvements in computing capabilities, e.g. those provided by graphics processing units (GPUs) lead to drastic shortening of simulation times and make time-consuming problems tractable [37][38][39]. To benefit from these new architectures we have implemented our PIC/MCC code on GPU architecture utilizing the CUDA programming language extension.
The geometry employed in the simulation is defined as an L = 400 mm long dielectric tube with 2R = 30 mm inner diameter and two flat metallic disk electrodes at the ends. A numeric grid in the (r, x) cylindrical coordinate space is defined with a resolution of N r = 64 grid points along the r-axis and N x = 2048 grid points along the x-axis providing an equidistant spatial resolution of Δr = 238 μm and Δx = 195 μm, resolving the Debye screening length in the plasma. This numerical grid is used for the solution of Poisson's equation by means of the parallel successive over-relaxation method and the interpolation of the electric field to the actual positions of every particle at each time-step. Due to the singularity of the cylindrical grid at r = 0 we also define an equi-volume mesh (Δr(r) = r 2 + R 2 /(N r − 1) − r) for the charge-to-mesh assignment step of the PIC cycle. The charge density distribution on this equi-volume mesh is then interpolated onto the equidistant mesh for further processing.
In the PIC/MCC simulation the large quantity of real particles is represented by a reduced number of super-particles, each of which represent a given number of real particles defined with the w super-particle weight parameter. In our simulations we have set w = 50 000, which resulted in superparticle counts on the order of 10 7 for each charged species in the converged solution.
The trajectories of the charged particles are integrated in time using the leap-frog method with a time step of Δt = 24.6 ps chosen to resolve the electron plasma oscillations, assuming only the electrostatic electric field, computed on the numerical grid from the global charge distribution, as the source of acceleration in the frame of classical Newtonian mechanics. Gas phase collisions between charged particles (electrons and Ne + ions) and the neutral background gas atoms are modeled using the standard Monte Carlo method as discussed in great detail in [32] using the scattering cross sections shown in figure 1. The gas temperature is set to T g = 300 K, which results in a neutral gas particle density of n g = 1.45 × 10 22 m −3 at p = 60 Pa gas pressure. Gas phase and surface processes that are included in the model together with the governing parameters are listed in table 1, while details are described below.
All individual binary collisions assume an isotropic scattering model in the center of mass system of reference. The initial velocity of the target background gas atoms is sampled randomly from a Maxwell-Boltzmann velocity distribution with temperature T g . In case of electron induced ionization, the energy of the impacting electron is reduced by the ionization energy (ε ion = 21.56 eV in this case), and the remaining kinetic energy is shared between the scattered impacting and the new emitted electrons with a uniformly chosen random ratio.
All particles that are traced in the discharge volume at some point hit one of the surfaces. In case of electrons arriving at one of the electrodes, there is a chance (defined by the electron reflection parameter η) that the electron gets reflected without energy loss, otherwise it gets absorbed and contributes to the electric current measured at the electrode surfaces. Ions reaching an electrode lose their charge via recombination at the surface and in addition can induce electron emission defined by the ion induced electron emission yield γ m . Similarly γ g defines the ion induced electron emission from the glass cylin-  [43] Surface: Electron emission γ g = 0.4 - Figure 1. Impact energy dependent scattering cross sections of gas phase processes for electron (a) [40] and Ne + ion (b) [41] collisions with neutral Ne atoms.
der. If charged particles reach the glass cylinder, they deposit their charge at the surface element of the x position of the impact, contributing to ρ s (x), the charge state of the surface element.
The electric field (potential) distribution is computed at every time-step by solving Poisson's equation on the equidistant grid with mixed boundary conditions defined in the following manner. One of the electrodes (at x = L) is at fixed ground potential. The driven electrode (at x = 0) has a uni-form potential Φ 0 , which is set by Φ 0 = V DC ± IR 0 , where I is the electric current measured by counting the charged particles crossing the driven electrode surface taking into account all absorption and emission processes and the direction of the current as expressed by the '±' symbol. R 0 = 20 kΩ is a series resistor typically used experimentally to enable stable operation of the DC discharge. V DC is the fixed voltage applied by the power supply. Neumann boundary conditions are used for the glass cylinder instead of the potential. The surface charge and with this the surface normal electric field is available as Metastable Ne m atoms, being charge neutral, are treated as a continuous constituent and their transport is incorporated by solving the standard cylindrical diffusion equation with the diffusion coefficient D m and the source term determined from individual electron induced excitation processes. The gas phase loss channel is taken into account through pooling ionization with the rate coefficient k pi . Although pooling ionization is a loss for Ne m atoms, at the same time it serves as an additional source for electrons and Ne + ions, which is coupled back to the charged particle tracking PIC/MCC model. The boundary conditions for the Ne m diffusion model assumes zero density (100% deexcitation probability) at all surfaces. Note that although the Ne m is included in our model, the results show that at the discharge conditions studied here, the contribution of the pooling ionization to the total charge production is only about 1%. From the perspective of the charge balance, this small value allows us to leave further discussion of metastables. However, previous experimental [44] and numerical [45] studies have shown that in the presence of dust the metastable atoms can contribute by heating of the grains and potentially inducing electron emission from their surface. As a result, metastables are efficiently quenched into the ground state decreasing their density in the dust cloud. Further, in the presence of molecular impurities, which are unavoidable in dusty plasmas due to evaporation of the grain material, Penning-ionization of the impurities can significantly affect the charged particle composition of the plasma.

The experiments
In order to validate our numerical computations we have used two experimental systems to obtain data that can be compared with the simulations and at the same time represent conditions relevant to the PK-4 experiment. First is the 'PK-4BU' experiment, which is a ground based replica of the original PK-4 flight model operated at Baylor University in Waco, Texas. The second is a simplified straight DC discharge operated in Budapest, Hungary and referred to as 'BUD-DC' in the following. Some relevant details of both systems follow.

PK-4BU experiment
A detailed introduction of the PK-4BU experiment together with its first characterization (Langmuir probe and spectroscopic) results can be found in a recent publication by Schmidt and Hyde [46]. Here we just point out that it is close to an exact replica of the PK-4 flight model in terms of the discharge tube geometry ('U' shaped glass tube with a total length of about 800 mm with a straight center part of about 400 mm length) and electric drive (hivolt.de HA51U programmable high voltage power supply).
The diagnostic relevant to the present work is a Photron FASTCAM Mini UX50 high-speed CCD camera that is used to record image sequences of the discharge light emission with frame rates up to 50 000 frames per second (fps) in the visible spectral range. The field of view of the camera was restricted to a ∼ 200 mm long section of the central straight part of the discharge.
Both the flight model and the PK-4BU system lack well defined electrode surfaces; the power supply is connected to the steel KF vacuum flanges that seal the glass tube interior from ambient environment. From the modeling aspect both the curved corners of the glass tube and the ill-defined electrode configuration are problematic, therefore we have introduced the BUD-DC experiment, that matches exactly the conditions of our PIC/MCC simulation (geometry and electrode configuration) while remaining as close as possible to the PK-4 system (diameter, gas, pressure range, etc).

BUD-DC experiment
In the BUD-DC experiment 6.0 purity neon gas is introduced by means of an Aalborg GFC17A mass flow controller with a constant gas flow of about 1 sccm into a glass cylinder with 30 mm inner diameter. In the glass tube two flat stainless steel electrodes are placed at a distance of 400 mm and connected to the high-vacuum flanges at the tube ends. DC voltage is applied using a PS325 Stanford Research DC power supply. Two-stage pumping is realised by a Leybold TURBOVAC-151 turbomolecular pump backed by an Adixen 205SDMHEM rotary vane pump. Base pressures down to 10 −7 mbar can be reached in this system.
Without the availability of high-speed cameras on this system, only images with long exposure times of the discharge were recorded utilizing an Allied Vision Prosilica GX1050 CCD camera. On the other hand, due to the simplified geometry, the whole discharge volume could be imaged, making comparison with simulations possible for more than just the positive column region. Besides imaging, fast single channel light intensity measurement was realized using a detection systems composed of a Hamamatsu H7732P-11 photomultiplier tube (PMT) based photon counting module fed with an optical fiber assembly collecting light from a cross section with ∼1 mm width from the discharge. The PMT signal was amplified using an FEMTO HCA-400M-5K-C current amplifier and the individual photo-current pulses were counted with an Ortec MCS-PCI multichannel scaler card using 2 μs dwell time and 65 536 time channels.

Experimental validation
It is the freedom of the model constructor to decide which processes to include and which to leave out as well as what approximations and simplifications to incorporate. In most cases a numerical simulation provides a converged solution which may properly or poorly represent the real physical system. The simulation on its own can not validate itself; independent, preferably experimental, data is needed to provide some kind of quality control of the numerical results. To provide experimental validation of our simulation we compare the global discharge structure recorded by imaging the whole discharge of the BUD-DC experiment at 60 Pa gas pressure at different discharge currents. Figure 2 shows the comparison of the experimental light intensity distributions and the lineof-sight integrated (Abel transformed) total electron induced excitation source function. Strictly speaking the comparison of these two quantities is not exact, but it becomes reasonable once we accept two approximations: (i) the life-time of the excited electronic states is short with respect to characteristic atom diffusion times, and (ii) the discharge plasma is optically thin, meaning that radiation trapping is negligible.
From the above, it can be concluded that the simulation properly captures the global discharge structure and its discharge current dependence. Starting from the cathode one can identify the cathode dark space (0 < x < 10 mm), the negative glow (10 < x < 60 mm), the Faraday dark space (60 < x < 90 mm), and the positive column (x > 90 mm). Both in the experiment and the simulations the positive column starts with two well defined bright features, followed by increasingly homogeneous light emission. Keeping in mind that in the PK-4 and the PK-4BU experiments the observation volume is deep in the positive column, these long exposure/acquisition time data support the observation of a stable homogeneous (along the x-axis) plasma environment. In the following we will focus our attention on this range of the positive column.
One of our key findings of this work is introduced in figure 3(a), where the computed excitation rate distribution is now shown for a short data acquisition time t acq = 1 μs. One can clearly identify high intensity striations, which persist throughout the whole positive column without losing contrast along the x-axis. The comparison between panels (a) and (b) of figure 3 illustrates the time evolution of the instantaneous structure with a 10 μs time difference between the plots. We observe the rapid motion of the bright blobs in the direction toward the cathode with different velocities. Panels (c) and (d) show the computed electron density distribution for the same times as in panels (a) and (b), and one can conclude that the bright features coincide with regions of high plasma density. Additionally the plasma density shows high intensity spatial variation along the x-axis reaching density ratios on the order of 10 between density maxima and minima.
Before drawing any conclusions we again have to address the question regarding the accuracy of our PIC/MCC simulations. To validate the reliability of the simulation data with respect to the moving striation we performed experiments on the PK-4BU system utilizing the high-speed camera setup. Comparison of the measured and computed space-time (x − t distributions) of the on-axis excitation source/emitted light are shown in figure 4. For this comparison the discharge parameters had to be adjusted to a higher pressure (p = 133 Pa) and to the maximum available current (3 mA) in order to achieve high enough light intensity needed for the short exposure time of the high speed imaging at 50 000 fps. The simulation parameters were adjusted accordingly.
We find that the experiment (panel (c)) fully supports the numerical finding (panel (b)) of the presence of the rapidly moving striations. The average spatial distance between the bright features is about 15%-20% larger in the experiment, and the propagation velocity is slightly lower compared to the simulation, but otherwise all prominent features (bending, brightness variation, bifurcation) appear very similar in both the numerical and experimental data. These x − t plots also provide a direct way to measure the propagation velocity of the moving striations, which is found to be in the range of 300 m s −1 and 600 m s −1 , preferring the lower value near the anode region and the higher value closer to the cathode.

Striations and ionization waves
In the following we will use the terms striations and ionization waves for axial inhomogeneities in the positive column, which are quasiperiodic and static, and quasiperiodic and moving phenomena, respectively. Ionization waves are often referred to as moving striations in earlier publications. In all cases the variations of the local plasma parameters are linked to variations of the luminosity of the plasma column.
DC discharges, being one of the most frequently studied systems of gaseous electronics [47], are fairly well understood with a large number of details available in the literature. Nevertheless, these systems are of ongoing interest still today [48][49][50][51][52][53]. It is known that in DC discharges self-excited oscillations may occur over a wide range of frequencies and discharge conditions [54]. Striations are most prominent in molecular gases or in contaminated rare gases, but can also coexist with ionization waves, observed as early as 1841 [55,56].
According to [57] striations are manifestations of ionization oscillations and waves caused by alternation of regions of predominant production and loss of electrons. The wavelength of these longitudinal instabilities is bounded from below by the tube radius R, since perturbations in n e are actively destroyed by radial ambipolar diffusion, and bounded from above by the electron temperature relaxation length, which is governed by the electron drift velocity and the momentum transfer collision frequency. Both theoretical results and fluid simulations support these properties [58][59][60][61][62][63][64]. A discharge parameter diagram for neon distinguishing regions of stability and several oscillatory modes has been constructed based on experimental data in [65]. According to this the PK-4 discharge operates in the 'diffusive striated' regime, in agreement with our observations. A theoretical model based on non-equilibrium electron transport in a simplified model discharge was derived in [66] concluding that the only phenomenon which can lead to the amplification of the instability is the spatial phase shift of the electron temperature curve with respect to the course of the electric field. This is caused by the non-local nature of the electron transport, meaning that the energy of the electrons at a certain location is given not only by the local electric field, but also depends partly on the potential difference the electrons have passed at earlier times along their trajectories.
Based on a recent review on the types and properties of ionization waves referred to as kinetic resonances in [50] we can identify two types of ionization waves in our system in figure 4(a). Based on the kinetic computation in [67] the principal wavelength is related to E x the axial electric field and 1 the first atom excitation level by L S ≈ 1 /(eE x ). This approximation is valid at low pressures, where the energy relaxation due to elastic electron-atom collisions is inefficient. In our case with E x ≈ 2.5 V cm −1 and 1 ≈ 16.6 eV we obtain L S ≈ 6.6 cm. The wave modes, referred to as integer resonances in [50] are characterized by wavelength L = L S /k, where k = 1, 2, 3, . . . integer values. Observations report the appearance of S-striations (k = 1) and P-striations (k = 2) in various discharges. In our case, the fast waves at the cathode side of the positive column have a wavelength around L fast ≈ 7 cm, while in the case of the slow waves near the anode L slow ≈ L fast /2. Figure 5 shows the axial distribution of the electric potential and the oscillations of the excitation source at an arbitrarily chosen time instance for the p = 133 Pa and I = 3 mA case. The potential difference between the plateaus is on the order of ∼ 20 V for striations on the cathode side and ∼ 10 V on the anode side, which are similar to those reported in [50] for S-, and P-striations in neon.
To more closely match the discharge parameters used in the dusty plasma experiments in PK-4 we continue our analysis the p = 60 Pa and I = 2 mA base case. Figure 6 depicts the computed space and time resolved excitation rates, which when In many of the earlier studies the wave modes in neon and argon appeared at distinct discharge conditions, typically P-striations at lower current and S-striations at higher currents [50]. However, bi-modal mode structures and transitions between them were reported in [68,69] for pulsed discharges and externally excited waves. To verify the bi-modal mode structure predicted by our simulations we have performed an experiment on the BUD-DC system utilising the fast single channel light detection system. The detector was placed at six different locations with 5 cm steps and time-series of the photon count rates were recorded and Fourier transformed. The series of measured power spectral densities are shown in figure 7. At the cathode side we find a strongest spectral feature with a frequency of f S = 24.3 kHz and a weaker peak at f P = 18.8 kHz. Approaching the anode side the bi-modal structure transforms gradually by changing the intensity ratios in favor of the f P peak. Very close to the anode at x ≈ 38 cm the f S peak is almost completely depressed and the governing feature in the spectra is a broad asymmetric peak with maximum close to f P and extending down to ∼ 10 kHz, confirming the lower degree of regularity of oscillations at the anode side already seen in figure 6 from the simulations.
Simulations with different discharge currents in the range between 1 mA and 4 mA under otherwise identical conditions at p = 60 Pa confirm that the axial position of the transition region between the wave modes moves away from the cathode with increasing current. As a result, if investigating the oscillations in a range around a fixed position one would observe P-striations at low currents and S-striations at high currents. Based on the above we can conclude that our simulation results are compatible with theoretical predictions and experimental observations of S-, and P-striations and we observe both wave modes and the transition between them spatially separated in our discharge.
Recent experiments on ionization waves [69][70][71][72][73][74] provide experimental spatiotemporal patterns with light-emission fluctuations similar to that shown in figure 4, providing independent experimental validation of our numerical results. The impact of the dust particles on the stabilization of the striations has been studied in [75] in the PK-4 system as well as by means of fluid type numerical discharge modeling.
The working hypothesis of this study is that the apparently stable PK-4 positive column includes fast, optically unresolved oscillatory features (ionization waves), which may be enhanced during the polarity switching procedure in a fast transient process. To support this hypothesis first we discuss the static global discharge parameters in detail, followed by an analysis of the time-evolution of the local plasma parameters at the position of the dust particles in the PK-4 experiment, namely the geometrical center of the discharge tube. Following the discussion of the stationary state, the effect of the polarity switching on the local plasma environment will be examined.  Figure 8 shows a set of plasma parameter distributions across the discharge volume. The simulation parameters employed are p = 60 Pa neon gas pressure and I = 2 mA discharge current. For this case the DC voltage on the powered electrode was Φ 0 ≈ −770 V. The data acquisition time is t acq = 1 μs, and all panels show the same time instance as already introduced in figures 3(b) and (d) for the excitation source and the electron density distributions, respectively.

Global structure
Comparing panel (a) to figure 3(d) we see that the ion density is closely correlated to the electron density with the important difference that the ion density decays to a higher value while approaching the glass surface. This is a signature of the presence of a sheath region in the vicinity of the surface and the negative charge state of the surface with respect to the plasma potential facing the surface. Additional detail of the radial electric field structure will be given later.
As shown, the ionization source distribution in panel (b) also follows the striated structure but shows a well defined shift of the maxima positions of the striation heads of approximately 5 mm toward the cathode with respect to the charged particle densities. The fact that the enhanced ionization precedes the density maxima during the wave propagation is an indication that the instability is driven by an inhomogeneous ionization profile.
The axial electric field profile in panel (c) shows that the regions of enhanced ionization are preceded by regions of high negative fields accelerating electrons toward the anode and that the electric field crosses zero at the position of the density maxima, followed by a low amplitude field reversal. These details suggest the striation features are sustained by a mechanism very similar to that found in the cathode region, where electrons emitted from the cathode are accelerated by the cathode fall, releasing their energy in the field free negative glow before being further decelerated by a weak field reversal. Also there seems to be a direct connection between the magnitude of the acceleration field and the amplitude of the density variation.
The radial electric field in panel (d), on the other hand, always acts to confine the electrons, and the strength of the cylindrical sheath electric field follows the striated structure. The region of strongest radial field coincides with the region of strongest axial electric fields and is weakest where the axial electric field shows a reversed (positive) sign.
The electron mean energy shown in panel (e) closely follows the distribution of the electric field incorporating both axial and radial contributions causing the strongly curved shapes appearing in the figure.
The axial electron drift velocity distribution in panel (f ) shows an unidirectional flow of electrons in the positive column despite the presence of the field reversal regions supporting the concept, arising from current conservation principles, that field reversals serve only deceleration (cooling) of the electrons but are not strong enough to reverse net motion. The axial velocity maxima coincide with the maxima on the axial electric field and the electron mean energy.
The radial electron velocity in panel (g) shows a weak and fluctuating radial motion of electrons near the center and a steady out-flux near the glass cylinder controlled by the radial ambipolar diffusion. At regions following the density maxima, where the axial field is reversed, the radial electron diffusion becomes reversed as well in the bulk region of the positive column.
From the aspect of the dust particle potential (wake fields and electrorheological effect) and inter-dust interactions the behavior of the ion species in such striated discharges is of key importance. The characteristic time scales in these ionization waves (on the order of 10 kHz) are slow enough for ions to react to the instantaneous local plasma environment in contrast to the case of RF (MHz) driven discharges. The Ne + mean energy distribution in panel (h) shows a strong global ion heating (radial acceleration) toward the glass cylinder by the sheath electric field. Within the bulk, in regions of large negative axial electric field the ions are accelerated but lose their energy quickly in regions of low electric field due to efficient momentum transfer between ions and the background gas.
The axial ion velocity profile shown in panel (i) illuminates an interesting feature of the ion transport. In contrast to the electrons, the ion transport is very close to equilibrium transport and the axial ion mean velocity changes direction in regions of field reversals. The places of zero axial velocity coincide with positions of the density maxima and the position of the axial electric field zero crossing from negative to positive. A consequence is that the global ion current in the positive column is carried by the motion of the striation in bunches of ions, which is further supported by the fact that ion drift velocity values match the velocity of the striation heads.
Another aspect of the dust particle structure in dusty plasmas is the radial confinement. Due to high dust particle inertia the dust component cannot directly respond to the several kHz field oscillations; it is only the time averaged field that plays a primary role. Time averaged radial density and radial electric potential/field profiles are depicted in figure 9 at the center cross-section at x = 200 mm.
Both the density profile in panel (a) peaking at the value n e ≈ n i ≈ 1.18 × 10 15 m −3 and the electric field distribution in panel (b) shows the presence of a quasi-neutral central bulk at r < 10 mm and an electron depleted sheath region at r > 10 mm. Thus, the electric field in the center of the bulk (r < 5 mm) can be well approximated using a harmonic trap with a linear electric field profile and time averaged slope of ∂E r /∂r = 8.94 V cm −2 .

Time dependent local plasma environment
Here we discuss the local plasma parameters that a dust particle would experience if placed in the geometrical center of our discharge tube (r = 0, x = 200 mm) in a DC discharge at p = 60 Pa and I = 2 mA. The purpose of this study is (i) to provide input data for future dust charging and dust-dust interaction models related to PK-4 discharges and (ii) to provide insight into the fundamental differences of a heavily striated positive column plasma environment with respect to a homogeneous situation. We emphasize that our simulation does not contain the dust component and the self-consistent mutual interaction between the dust and the discharge plasma. For this reason we expect our results to be relevant to cases where Z d n d n e , the dust charge density is low with respect to the electron density. The effect of high dust charge density on the striation structure has been previously discussed in [75]. Figure 10 shows the time evolution of the plasma density (a), the axial electric field (b) and the ionization rate (c) at the center of the discharge tube. The immediate conclusion is that we observe a large amplitude, non-linear oscillation of all quantities with a governing frequency of about f = 24 kHz. The electron and ion densities are not distinguishable at the resolution of the graph, even at the zoomed right section magnifying the crossing of a single ionization wave front. The plasma density oscillates between n min ≈ 0.4 × 10 15 m −3 and n max ≈ 2.7 × 10 15 m −3 representing a contrast of almost 7:1. The axial electric field oscillation is even stronger, alternating time intervals of low positive values with peaks with large negative amplitudes, representing the crossing of the field reversal. The ionization rate at the discharge center shows a pulsed-like waveform with a zero baseline and sharp peaks. Comparing the order of appearance of these wave fronts as shown on the magnified right side of the graphs one can clearly identify a sequence of causal connections. First the high electric field pulse causes strong electron heating, wherein high energy electrons overcome the ionization threshold energy and contribute to enhanced charge production. These charges accumulate and cause a rapid increase of the plasma density. The ions are then carried away by drift and ambipolar diffusion on a longer time scale on the order of several tens of microseconds.
In the theoretical models of both the ion wake field and the electrorheological effect the ion drift motion plays a central role. The key quantity for the characterization of the streaming plasma environment is the local velocity distribution function (EVDF for electrons and IVDF for ions) at the position of the dust particle (which would be at the center of the discharge in our case). To quantify the time evolution of the EVDF and the IVDF we compute the first four moments, the mean (v, representing the drift velocity), the standard deviation (σ, representing the kinetic temperature), the skewness (S), and the kurtosis (K) based on the following particle sums: where N is the accumulated number of particles in the central numerical grid cell sampled in every time step for the duration of the short data acquisition (t acq = 1 μs in our case) period. The sums run over these particles found only in the discharge center. Here and in the following the bar notation e.g. inv represents the ensemble average collected over a short data acquisition time t acq , while the long time average will be denoted by brackets below. Here the velocity appears as a scalar quantity, which represent the axial, radial, and tangential vector components individually. For an ideal Maxwell-Boltzmann velocity distribution the drift velocity v drift =v, the kinetic temperature k B T = mσ 2 , the skewness S ≡ 0 and the kurtosis K ≡ 3. Any deviation of the last two moments from these constant values is a direct measure of the deviation of the whole distribution function from thermal equilibrium. Figure 11 shows the time evolution of the first four moments of both the EVDF and the IVDF at the center of the discharge operated as before at p = 60 Pa and I = 2 mA. Focusing first on the electrons we observe that all four moments of the axial velocity distribution closely follow the high amplitude oscillations already discussed in figure 10. The odd numbered moments of the radial and the tangential components show only statistical fluctuations around zero, while the even numbered moments coincide with the axial component. This parity rule simply means that the distributions of the offaxis components are symmetric around v rad/tan = 0 but otherwise share the non-equilibrium nature of the axial component. The mean axial velocity in panel (a) fluctuates between v min = 2 × 10 4 m s −1 andv max = 16 × 10 4 m s −1 with a long time average v = 6.5 × 10 4 m s −1 . These values are put into perspective when compared to the standard deviation shown in panel (c), which fluctuates in-phase with the mean velocity but between limits that are about an order of magnitude higher. This property is in agreement with EVDF results in other discharge plasmas, like in [76], where the randomized isotropic component of the EVDF dominates and the drift component can be viewed as an anisotropic perturbation. The oscillation of the second moment corresponds to an oscillation of the isotropic electron temperature between 2.6 eV and 6.6 eV. The fact that the radial and tangential components copy the behavior of the axial component means that the isotropization and coupling between velocity components is effective on the submicrosecond time scale. Despite the fast isotropization of the electron velocity, the thermalization with the background gas appears to be a much slower process, which is clearly shown by the non-equilibrium values of the third and fourth moments as shown in panels (e) and (g). The evolution of the skewness in panel (e) supports the ionization driven nature of the wave development because it shows that during axial acceleration the high energy tail becomes enhanced (positive skewness), while during the fast decay of the axial mean velocity the high energy tail becomes depopulated by inelastic collisions (negative skewness), such as excitation and ionization, the reaction channels of which are characterized by threshold energies, therefore acting only on high energy electrons. This variation of the EVDF was recently investigated experimentally in a stratified argon DC discharge in agreement with our results [53].
The IVDF properties significantly differ from the EVDF. The mean velocity in panel (b) shows not only an oscillation of the magnitude but also changes in direction between the large amplitude peaks. The average drift velocity | v | = 120 m s −1 is somewhat below the value one would expect in case of equilibrium transport in a constant electric field that corresponds to the mean value of the fluctuating axial field | E x | = 2.18 V cm −1 (or in terms of the reduced electric field of E /n g = 15 Td at 60 Pa pressure) based on [26], which is about v drift (15 Td) = 160 m s −1 . This difference illuminates an important consequence: a high amplitude oscillation together with non-linear transport properties (v drift (E/n)) can cause substantial deviation of the time averaged quantities compared to equilibrium values in homogeneous and constant fields. The average velocity value is far below the thermal velocity depicted in panel (d) that is on average v th = 455 m s −1 , representing a thermal Mach number of M th = 0.26, which is too low to form stable particle chains according to [20]. The amplitude of the mean velocity oscillation, on the other hand reaches maximum values around |v max | = 700 m s −1 , which is almost equal to the thermal velocity at the same time instance (v th,max = 760 m s −1 ) resulting in instantaneous thermal Mach numbers during the phase of the oscillation maxima around unity. This represents a large value causing potentially strong dust particle chain formation according to [20]. The fact that the third and fourth moments of the IVDF show only statistical fluctuations around the equilibrium values indicate that the ion transport is close to equilibrium. The thermalization of the ions with the background gas occur on the microsecond time scales. This supports our implicit assumption that the standard deviation simply relates to the thermal velocity as v th = σ = k B T/m. The base value of σ between the peaks is σ ≈ 400 m s −1 , which is about 12% higher than the values one would expect in the case of perfect thermalization with a 300 K background gas, a consequence of the driven-dissipative nature of the ion dynamics.
Here we emphasize that the subject of the present study is the DC neon glow discharge of the PK-4 experiment, and the effects of the dust component on the properties of the gas discharge are not included in the model. This condition limits the validity of the quantitative result given in the discussion above to cases with low dust densities. Earlier experiments have shown that the presence of high dust densities cause charge depletion [77] by enhancing surface recombination and by substituting negatively charged dust grains for part of the electrons in the charge balance. In this case the discharge current has to be carried by fewer electrons, as the dust grains are practically immobile, resulting in the increase of the longitudinal electric field in both DC [78,79] and RF discharges [80].

Polarity switching
In the PK-4 experiment, in order to prevent the dust particle cloud from permanently drifting, the polarity of the applied DC voltage can be switched between positive and negative values with up to a 1 kHz repetition rate. The power supply installed in the experiment provides typical voltage rise times of about 200 μs. For this reason we took the simulation at p = 60 Pa and I = 2 mA discussed so far and introduced a linear voltage ramp between V DC = −820 V to V DC = 820 V with a transit time of 200 μs. Figure 12(a) shows the time evolution of Φ 0 , the electric potential at the driven electrode and the global discharge current. The current characteristic verifies that a voltage switching time of 200 μs is long enough for the discharge to completely die so that it has to re-ignite during every polarity switching cycle. Nevertheless it is interesting to observe the asymmetry of the off-on cycle in the current profile. The discharge decays slowly to zero charge density (in the positive column), as shown in panel (c), but requires reaching the breakdown voltage to ignite and another ∼ 100 μs time to reach stationary current. During the ignition phase the axial electric field in panel (b) shows an extra pre-ignition peak at t ≈ 180 μs in coincidence with a small current peak. This electric field peak is followed by a small density peak, but more importantly from the aspect of dust particle dynamics is the development of a strong radial electric field gradient as shown in panel (d), which is about two times stronger compared to the periodic peaks observed during the stationary ionization waves. This strong radial confinement pulse contributes significantly to the stabilization of the dust particle chain and originates from the interaction of the first population of charges created during preignition and the charges on the glass surface that survived the polarity switch cycle.

Summary
To date, theoretical investigations of observed dust particle chain formation in the PK-4 system have assumed a time-independent stable plasma environment. In this work, a numerical PIC/MCC gas discharge simulation of a DC neon discharge with parameters relevant to the PK-4 micro-gravity dusty plasma experiment has shown that the apparently quiet positive column assumed in the PK-4 is actually a region of rapidly moving high amplitude ionization waves. As a result, all local plasma parameters relevant to dust particle charging, confinement and chain formation fluctuate at frequencies on the order of 10 kHz with amplitudes up to 10 times the time-averaged values.
This highly oscillatory environment directly impacts the electrorheological effect, currently the most prominent candidate mechanism for chain formation. As introduced in equation (1) and discussed above, the potential contribution due to a streaming plasma environment depends on the ion drift via the ion thermal Mach number as ϕ stream ∼ M 2 th . Although dust particle dynamics are largely insensitive to fluctuations on this time-scale leaving expectations that only timeaveraged quantities matter, in the case of such a non-linear phenomena, these results can vary significantly depending on where the time-averaging is performed. In this particular case (p = 60 Pa and I = 2 mA neon PK-4 discharge) averaging the plasma parameters first yields M 2 th = ( v drift / v th ) 2 = 0.069, while computing the time-dependent thermal Mach number and computing its time-average yields M th 2 = ( v drift /v th ) 2 = 0.097.
The difference between these represents a 40% stronger anisotropy of the inter-particle interaction potential. The anisotropy term becomes even larger if averaging is performed after the square: M 2 th = v 2 drift /v 2 th = 0.146, yielding an enhancement factor of 2.1 over the initial result. The exact anisotropy value sufficient to reproduce the experimentally observed dust particle chain formation will need to be determined from future detailed dust dynamics simulations. However, the radial confinement (i.e., the trap frequency) will certainly be a key ingredient in those calculations, given this is the point where the polarity switching can contribute. In this paper we have shown that during the 200 μs long polarity transition the discharge must completely re-ignite allowing a buildup of a strongly enhanced transient radial confinement field during the pre-ignition phase. When repeated at a kHz rate, this may be able to stabilize dust particle chains forming in the stationary period of the discharge operation.