Nonlinear dynamics in dusty plasmas subjected to photo-discharging

Abstract

Experimental research into the control of particle charge in dusty plasmas conducted at Auburn University indicates that photocurrents generated by exposing dust to intense, near-ultraviolet light can provide a reliable and novel method of independently controlling dust charge without radically altering the background plasma; the experiment also showed that some particles may respond differently to this photo-discharge, with some exhibiting highly periodic responses to the discharge and others exhibiting chaotic behaviour. Since the dust particles in the experiment were a polydisperse sample of different sizes and shapes, particle geometry may play a role in explaining this difference. Simulations of particle discharge and dynamics are used in an attempt to reproduce experimental results and investigate a possible correlation between particle symmetry and dynamic periodicity.

1 Introduction

A dusty or complex plasma is a plasma system consisting of electrons, ions and neutral gas particles, in addition to an ensemble of charged microparticles (dust). These particles collect charge from the surrounding plasma environment and, consequently, become a third charged component, coupled to the rest of the plasma. The addition of this third charged species alters the existing plasma dynamics, and the comparatively larger scale and inertia of the dust leads to unique phenomena, including dust waves (Barkan, Merlino & D'Angelo 1995), dust crystals (Thomas et al. 1994) and dust voids (Goree et al. 1999).

While the dust material indirectly affects its equilibrium properties and dynamics, dusty plasmas are not limited to any particular composition, and frequently appear or form in almost any plasma regime, including particles ablated off plasma-facing components in fusion experiments (Rubel et al. 2018), in situ particle growth in processing plasma devices (Watanabe 1997), ice crystals in planetary rings (Horányi 2004) and trace atmospheres of interplanetary or orbiting bodies (Morooka et al. 2011).

The coupling between dust and plasma properties also makes dust particles useful as a minimally perturbative probe (Thomas & Morfill 1996), and a tool for investigating soft-body dynamics, statistical mechanics and imposed or driven behaviours in plasmas (Wong et al. 2017; Hall et al. 2018). These applications are thanks, in large part, to the significantly greater inertia of dust, compared with the other species in plasmas. This inertia results in dynamics that typically occur on macroscopic time and length scales. Scattered light from particles illuminated by a laser or other intense light source can be observed using inexpensive cameras or the naked eye. Laboratory experiments using moderate frame rate cameras ($\sim$100 frames per second) can capture time-resolved dust behaviour, and video analysis techniques like particle tracking velocimetry (PTV) or particle image velocimetry (PIV) (as described by Thomas 1999; Lynch, Konopka & Thomas 2016) can allow researchers to construct velocity distribution functions or other physical data about dusty plasma systems from direct visual observation instead of inferring such information from indirect instrument measurements (Fisher & Thomas 2010).

1.1 Forces

In non-magnetized dusty plasmas, the dynamics of the dust particles are primarily governed by four main interactions; two field forces and two collisional forces:

(1.1)\begin{equation} m_d\boldsymbol{a_d}=\boldsymbol{F_g}+\boldsymbol{F_e}+\boldsymbol{F_{nd}}+\boldsymbol{F_{id}}, \end{equation}

where $m_d$ and $\boldsymbol {a_d}$ are the mass and acceleration of the dust. The field forces are the gravitational force, $\boldsymbol {F_g}$, and electric force, $\boldsymbol {F_e}$, defined as usual:

(1.2)\begin{gather} \boldsymbol{F_g}=m_d\boldsymbol{g}, \end{gather}

(1.3)\begin{gather}\boldsymbol{F_e}=Q_d\boldsymbol{E}, \end{gather}

where $Q_d$ is the dust charge, and $\boldsymbol {g}$ and $\boldsymbol {E}$ are the gravitational and electric fields, respectively.

There are also collisional interactions between the dust and the background plasma species. Collisions between the dust and neutrals in a weakly ionized plasma produce a neutral drag force, $\boldsymbol {F_{nd}}$, while direct and Coulomb collisions between the dust and mobility-limited background ions create an ion drag force, $\boldsymbol {F_{id}}$. For spherical particles in a low temperature plasma, the neutral drag is generally given by the Epstein drag model:

(1.4)\begin{equation} \boldsymbol{F_{nd}} ={-}m_n n_n {\rm \pi}{r_d}^2 C_n v_{tn}\boldsymbol{v_d}, \end{equation}

while the ion drag is typically described using the Hutchinson–Khrapak model (as summarized by Hall & Thomas (2016) and derived by Khrapak et al. 2002; Hutchinson 2003, 2004, 2006):

(1.5)\begin{equation} \boldsymbol{F_{id}} = m_i n_i {\rm \pi}{r_d}^2 (H+K\ln{\varLambda}){v_{ti}}^2\boldsymbol{\hat{E}}, \end{equation}

where $m_s$, $n_s$ and $v_{ts}$ are the mass, number density and thermal velocity of each species, $r_d$ is the dust radius, $\boldsymbol {v_d}$ is the dust velocity with respect to the neutral background, $\ln {\varLambda }$ is the Coulomb logarithm, and $C_n$, $H$ and $K$ are unitless coefficients of integration. While these coefficients are introduced here to simplify the expression of the forces, it should be noted these are non-trivial functions of the dust properties and plasma conditions.

1.2 Charging

Dust particles in the plasma will collect electrons and ions from the surrounding environment in the form of charging currents. In the orbital motion-limited (OML) charging model, these currents generally take on the following form (Shukla & Mamun 2002):

(1.6)\begin{gather} I_{ed}={-}e n_e {\rm \pi}{r_d}^2 v_{te} C_e, \end{gather}

(1.7)\begin{gather}I_{id}={+}e n_i {\rm \pi}{r_d}^2 v_{ti} C_i, \end{gather}

where $e$ is the electron charge, and $C_e$ and $C_i$ are unitless coefficients. Like the coefficients introduced in (1.4) and (1.5), these are non-trivial functions of the dust properties and plasma conditions. These charging currents exist simply by virtue of the dust existing in a plasma environment, but any additional physical process which can add, remove, release or neutralize charge on or from the dust will also result in a current.

Photoelectric currents generated by exposure to energetic photons can strip electrons from the surface, reducing the net negative charge of the dust. For negatively charged dust, this photocurrent, $I_{pe}$, can be modelled as

(1.8)\begin{equation} I_{pe}={+}e J_p Y_p A_{{\rm proj}}, \end{equation}

where $J_p$ is the photon flux, $Y_p$ is the photoelectron yield and $A_{{\rm proj}}$ is the projected area of the particle along the direction of the light source.

The instantaneous charge of the dust particle is typically assumed to satisfy a floating current condition where the charging and discharging currents balance each other out:

(1.9)\begin{equation} I_{ed}+I_{id}+I_{pe}=0. \end{equation}

Since the charging currents and forces acting on the dust are all functions of the dust properties and plasma conditions, this produces a loop between the dust charge and the dust dynamics: the dust's position determines the charging currents, which change the charge; while the dust's charge determines the forces which change the position. This constant feedback loop makes the charging behaviour and dynamics of dust in a plasma highly susceptible to disruptions and changes in initial conditions.

1.3 Photo-discharging experiment

Recent experimental research at Auburn University into methods of controlling dust charging behaviour have demonstrated the viability of using high-intensity, near-ultraviolet (NUV) solid-state sources to discharge dust without radically altering the surrounding plasma environment (McKinlay & Thomas 2021).

Figure 1 shows a schematic of the Auburn Dusty Plasma Experiment (DPX) during the photo-discharging experiment. The device consists of two stainless-steel six-way crosses: The primary chamber on the left contains electrodes which generate the plasma and suspend and confine the dust, and a tray containing the dust sample; while the secondary chamber on the right contains gas and vacuum access. Windows on the sides and front of the apparatus provide access to a camera, a laser for illuminating the dust, and a Solis-365C high-intensity microscopy source, which generates NUV light at 365 nm ($\sim$3.4 eV).

Figure 1. A schematic of the Auburn DPX apparatus during the photo-discharging experiment. The primary chamber (left) contains a plasma generating anode (disk) connected to a high voltage, current-limited power supply; a confining electrode (ring) connected to a low-voltage DC bias; and a tray containing the dust sample. The secondary chamber (right) provides access to a vacuum pump and argon gas via a mass flow controller. Windows on the sides and front provide access to optical equipment.

In the experiment, particles of lanthanum hexaboride (LaB$_6$) were suspended in an argon plasma with a neutral pressure of $\sim$120 mTorr, an electron temperature of $\sim$2 eV and an ion density of ${\sim }{10}^{14}$–${10}^{15}\,{\rm m}^{-3}$. A background electric field ranging from $\sim$2 to 15 V cm$^{-1}$ suspended the LaB$_6$ in the plasma near the centre of the primary chamber, either as small dust clouds (with densities of ${\sim }{10}^9$–${10}^{10}\,{\rm m}^{-3}$) or as sets of isolated particles. LaB$_6$ was chosen for the experiment because its low work function ($\sim$2.7 eV); the material has a modest photoelectron yield at room temperature (${\sim }{10}^{-7}$–${10}^{-6}$) and is often used as a base or coating for hot cathodes (Torgasin et al. 2017). Figure 2 shows a sample of the LaB$_6$ as seen under a scanning electron microscope (SEM).

Figure 2. (a) Samples of the lanthanum hexaboride (LaB$_6$) particles used in the experiment, imaged by an SEM. The particles have a wide distribution of sizes and shapes compared with particles used in typical laboratory experiments. (b) A dust particle undergoing periodic displacement in response to pulses of NUV light during the photo-discharging experiment. The image shows the overlay of $\sim$20 seconds of video frames, or $\sim$19 complete pulses (McKinlay & Thomas 2021).

The NUV source was connected to a function generator, which output a signal with a 50 % duty cycle at a frequency of 1 Hz, exposing the interior of the chamber to periodic flashes of NUV light. Several isolated particles in the path of the light exhibited significant, periodic changes in displacement in sync with the NUV pulses, consistent with a reduction in the dust charge. Figure 2 shows an example of this behaviour (the image shows the result of overlapping frames from $\sim$19 pulses).

While the photo-discharging experiment opens up new possibilities for controlling dust equilibrium and dynamic behaviour in laboratory and practical plasmas, it has also inspired new questions about dust charging, and the effect of dust shape and composition on the behaviour seen in these photo-discharging experiments. Figure 3 shows the displacement of two different dust particles over time and their motion through phase space in response to the NUV exposure. Some of the isolated particles exhibited more consistent, periodic behaviour (like the particle emphasized in the overlaid image in figure 2), while others exhibited more erratic, chaotic motion in response to the NUV pulses (like the particle in the centre of the image).

Figure 3. (a,b) Vertical positions of two example particles from the photo-discharging experiment as a function of pulse-time (with the pulse turning on at 0.0 s and off at 0.5 s). One particle in panels (a,c) exhibits periodic motion, while another in panels (b,d) exhibits more chaotic motion. (c,d) Same behaviour plotted in the velocity versus position phase space. Each plot contains data from 19 complete pulses. The colour and brightness of the data are used to show the passage of time (as indicated by the bar at the top of each figure); becoming darker with each pulse.

Since the spatial and temporal profiles of the NUV pulses are identical, this suggests that the cause must be intrinsic to the properties or presence of the dust material. As shown in the SEM image in figure 2, the LaB$_6$ used in the experiment was a highly polydisperse sample of different sizes and shapes, which is in sharp contrast to the highly monodisperse microspheres of silica (${\rm SiO}_2$) or melamine formaldehyde (MF) used in most laboratory experiments. This difference in particle geometry suggests one possible hypothesis for the difference in dust responses to the NUV.

1.4 Hypothesized role of particle geometry

Unlike the collection of background electrons and ions by the dust, the photocurrent is anisotropic – only the illuminated region of the surface is discharged. If we assume, as in (1.8), that the current is proportional to the particle's projected area, then as a spherical dust particle tumbles about on its axes, its cross-section (and the resulting discharge) will remain constant with each pulse; a more irregularly shaped particle, meanwhile, will have a cross-section (and discharge) which changes over time (an illustration of this can be seen in figure 4). Because of the feedback loop between dust charging and dust position, this time-dependent discharge could quickly result in increasingly chaotic behaviour.

Figure 4. An illustration of how particle geometry may affect photo-discharging. The highly symmetric spherical particle in panel (a) will have the same cross-sectional area, and thus the same photocurrent, for a light pulse of a given intensity, while the ellipsoidal particle in panel (b) has a cross-section (and thus photocurrent) which may change as it rotates.

Since the instantaneous charge distribution of the dust cannot be measured, an experimental or computational analogue must be used to test the plausibility of this hypothesis. The investigation presented here attempts to answer the following questions. Can a computational model, in which the dynamics and charging behaviour of a freely rotating particle depend on its projected area, reproduce the different periodic and chaotic motions seen in the Auburn experiment? If so, does the model indicate a correlation between the periodicity of a particle's response to photo-discharging, and its shape?

2 Computational model

For computational feasibility and to simplify the analysis of any simulation results, some assumptions and approximations are introduced into the models of the particle's geometry, dynamics and charging behaviour. The models used in the simulations and the justifications for those models are presented here.

2.1 Particle geometry

While we could choose any shape for our test particles (even shapes as irregular as those seen in the SEM images in figure 2), simulating highly irregular shapes, especially with sharp edges or concavities, introduces a degree of computational complexity beyond the scope of this investigation. The simulated particles used here are modelled as uniform ellipsoids with three specified axial lengths, not unlike some models used to interpret light curve data from asteroids and other small bodies (Zappalá 1980), which provides a compromise between too much geometric complexity and too little. This choice also implies that for a parallel light source, exactly half of the surface area of any particle will be illuminated at any given instant, regardless of orientation. This will inform the model used for the discharging of the dust outlined in § 2.3.

For the purposes of describing particle shape, we will introduce a parameter, $\kappa _{ij}$, which we will here define as the natural logarithm of a ratio of axial lengths:

(2.1)\begin{equation} \kappa_{ij} = \ln{\left(\frac{r_i}{r_j}\right)}, \end{equation}

where $r_i$ and $r_j$ represent any two axial lengths of a particular ellipsoid. Figure 5 shows the simulated surfaces or ‘hulls’ of four example particles: a sphere, a prolate spheroid, an oblate spheroid and a scalene ellipsoid; all illuminated by a simulated parallel light source and with some arbitrary orientation. Table 1 shows the axial lengths and $\kappa$-values for the same particles. These four examples are included among the simulated particles, and their resulting dynamics are discussed explicitly later on.

Figure 5. Simulated surfaces or ‘hulls’ of four example particles: a sphere (a), a prolate spheroid (b), an oblate spheroid (c) and a scalene ellipsoid (d). The lighter violet colouring indicates the region of the particle surface illuminated by a simulated parallel light source pointing from the left (as viewed here). While half of the surface area of these regular, convex particles is always illuminated by the pulse, the cross-section of those particles which are more elongated or flattened will change with their orientation.

Table 1. Axial lengths (in $\mathrm {\mu }$m) and $\kappa$-values for the example particles in figure 5. A $\kappa$-value of 0 indicates symmetry about an axis, while positive and negative values indicate stretching or compression along an axis, respectively. These particles were among those simulated and are discussed later as part of the simulation results.

Transforming the axial lengths into this ‘$\kappa$-space’ eliminates any dependence on a particle's physical size while preserving information about its shape. It also greatly simplifies our description of particle symmetry, as well as our visualization of the data produced by the simulation. Since the ratios of axial lengths must satisfy

(2.2)\begin{equation} \frac{r_x}{r_y} \frac{r_y}{r_z} \frac{r_z}{r_x} = 1, \end{equation}

the $\kappa$-values must therefore satisfy

(2.3)\begin{equation} \kappa_{xy} + \kappa_{yz} + \kappa_{zx} = 0, \end{equation}

which describes a plane in $\kappa$-space. We can also define three ‘$\kappa$-axes’ within this plane:

(2.4a)\begin{equation} \hat{\kappa}_{xy}=\frac{1}{\kappa_0}(0,-1,+1),\quad \hat{\kappa}_{yz}=\frac{1}{\kappa_0}({+}1,0,-1),\quad \hat{\kappa}_{zx}=\frac{1}{\kappa_0}({-}1,+1,0), \end{equation}

where $\kappa _0 = \sqrt {2}$. While these vectors do not form an orthogonal basis, they provide a convenient coordinate system for describing particle symmetry: particles at the origin of $\kappa$-space are spheres, while those located along the axes above are spheroids, which become more elongated or prolate as one moves in the positive direction and more flattened or oblate as one moves in the negative direction. The farther a particle is from the origin, the more exaggerated its shape becomes (more splinter-like or flake-like), and the further off-axis a particle is, the more asymmetric or scalene it becomes. With a simple rotation, all of the information about the shapes and motion of the simulated particles can be visualized as a single two-dimensional heat map.

2.2 Motion and forces

As outlined in § 1.1, the field forces to which the dust is subjected are trivial to model, while the collisional forces (which depend on particle shape) are not. For the purpose of this investigation, translational motion and forces are restricted to the vertical axis. However, the particles are allowed to rotate freely about all three principal axes. The spatial profiles for the plasma potential, temperature and density are, likewise, restricted to one dimension. While this model was adopted primarily for computational simplicity, it is not without experimental justification; Langmuir probe results from the photo-discharge experiment indicate that the systematic changes in the plasma properties were predominantly oriented along the vertical axis.

The gravitational and electric forces remain as described in (1.2) and (1.3), but because of the complexity and geometry dependence of the models for collisional forces, several approximations were made for the neutral and ion drag forces shown in (1.4) and (1.5):

(2.5)\begin{gather} \boldsymbol{F_{nd}} ={-}m_n n_n A_z C_n v_{tn} \boldsymbol{v_d}, \end{gather}

(2.6)\begin{gather}\boldsymbol{F_{id}} = m_i n_i A_z(H+K \ln{\varLambda}){v_{ti}}^2\boldsymbol{\hat{E}}. \end{gather}

Here, $A_z$, the instantaneous projected cross-sectional area of the ellipsoidal particle along the vertical or $z$-axis, replaces the circular cross-section of the spherical particle. A value of $C_n \sim 10$ was used for the neutral drag coefficient (this value is of the order of the usual coefficient determined in the Epstein model for a spherical particle and proved to be adequate for producing damped motion similar to that observed in the experiment), while the functions H, K and $\ln {\varLambda }$ are evaluated using the simulated particle's instantaneous charge and treating it as a sphere of equal volume. While these approximations may slightly oversimplify the effects of the particle shape on the collisional interactions, they should at least yield forces of an equal order of magnitude to what more complex models would yield, while still taking into consideration that the particle should experience proportionally greater or smaller collisional forces as the projected cross-sectional area increases or decreases as the particle rotates.

2.3 Charging and discharging

As discussed in § 1.2, even the simplified OML charging model for spherical particles is non-trivial and is only one model of the charging behaviour of the dust. The objective of this computational research is not to exhaustively investigate the charging of ellipsoidal particles, but merely to investigate the effect (if any) that particle asymmetry has on particle dynamics in the case of photo-discharging. For this reason, the equilibrium charge and discharging model of the simulated particles is also simplified.

The instantaneous charge distribution of the dust is modelled as a dipole: two point charges, $Q_1$ and $Q_2$, located at the geometric centroids of the illuminated and unilluminated surfaces. This gives the particle an instantaneous total charge,

(2.7)\begin{equation} Q_d=Q_1+Q_2, \end{equation}

and an instantaneous dipole moment, $\boldsymbol {p}$,

(2.8)\begin{equation} \boldsymbol{p}=(Q_1-Q_2)\boldsymbol{d}, \end{equation}

where $\boldsymbol {d}$ is the centroid of the illuminated surface. The charges on the two halves of the particle are modelled as

(2.9)\begin{gather} Q_1=(1+\delta_Q)\tfrac{1}{2} Q_{d0}, \end{gather}

(2.10)\begin{gather}Q_2=\tfrac{1}{2}Q_{d0}, \end{gather}

where $Q_{d0}$ is the initial charge of the particle (chosen to give the particle a random equilibrium position near the centre of the simulated chamber) and $\delta _Q$ is the fractional discharge on the illuminated side of the particle when the simulated NUV pulse is on. For our investigation, this discharge is modelled as being proportional to the ratio of the illuminated surface's projected area to its surface area:

(2.11)\begin{equation} \delta_Q={-}I(t)\frac{A_x}{S / 2}, \end{equation}

where $A_x$ is the projected area normal to the simulated source, $S$ is the total surface area of the particle (recall that half the surface is illuminated, regardless of orientation) and $I(t)$ is a unitless function of time which represents the light source turning on (1) and off (0) at a frequency of 1 Hz with a 50 % duty cycle, and a resistance and capacitance (RC) time of 0.1 s to approximate non-instantaneous charging and discharging. In the experiment, the intensity of the source is also a function of position, though for this investigation, it was treated as uniform.

In this simplified charging model, a spherical particle would lose half the charge from its illuminated side (25 % of its net charge) with each pulse. For contrast, in this model, a disk-like particle would lose all of the charge on its illuminated side when its face is normal to the source. When the simulated source is off, the charges on the two sides equalize, the total charge returns to its initial equilibrium value and the dipole moment reduces to zero.

2.4 Rotation and torque

The final piece of the model is the rotational dynamics of the particle. The interaction between the background electric field and the dust dipole moment gives rise to a torque:

(2.12)\begin{equation} \boldsymbol{\tau_e}=\boldsymbol{p}\times\boldsymbol{E}, \end{equation}

which twists the particle, trying to align its dipole moment with the field. This is resisted by torque from neutral collisions as the particle rotates, which has the differential form:

(2.13)\begin{equation} \boldsymbol{d\tau_{nd}}=\boldsymbol{r}\times\boldsymbol{dF_{nd}}={-}m_n n_n v_{tn} (\boldsymbol{r}\times\boldsymbol{v_d})\,\mathrm{d}S, \end{equation}

which can be integrated over the dust's surface area. If a no-slip condition is applied to the flow of neutrals normal to the dust's surface at a point, $\boldsymbol {r_s}$, as it rotates with an angular velocity, $\boldsymbol {\varOmega }$,

(2.14)\begin{equation} \boldsymbol{v_d}=\boldsymbol{\varOmega \times r_s}. \end{equation}

Then, this implies

(2.15)\begin{equation} \boldsymbol{r_s}\times\boldsymbol{v_d}=\boldsymbol{r_s}\times (\boldsymbol{\varOmega \times r_s})=\boldsymbol{\varOmega}{r_s}^2 \end{equation}

and a neutral drag torque,

(2.16)\begin{equation} \boldsymbol{\tau_{nd}}={-}\int \int m_n n_n v_{tn}\boldsymbol{\varOmega}{r_s}^2\,{\rm d} S. \end{equation}

If the surface is modelled as a set of discrete points (as is the case in our simulation), then the sum can be reduced to the total surface area, $S$ times a quantity, $\varSigma {r_s}^2$, which can be calculated numerically. This reduces the expression for the torque to

(2.17)\begin{equation} \boldsymbol{\tau_{nd}}={-}\sum_i m_n n_n v_{tn} \boldsymbol{\varOmega}{r_s}^2\Delta S_i={-}(m_n n_n v_{tn}S\varSigma {r_s}^2)\boldsymbol{\varOmega}. \end{equation}

Finally, if we impose that the rotational motion of the particle is mobility-limited by the neutral drag, then this implies an instantaneous rotational velocity:

(2.18)\begin{equation} \boldsymbol{\varOmega}=\frac{\boldsymbol{p}\times\boldsymbol{E}}{m_n n_n v_{tn}S\varSigma {r_s}^2}. \end{equation}

In principle, there should also be some torque due to uneven distribution in the ion drag along the surface of the dust. In lieu of developing a new theoretical model for this effect, a small, random, background rotational velocity was given to each particle, serving as statistical noise in the results.

The models for the forces, discharge and torques laid out in this section create a feedback loop between the position, orientation and charge distribution of the particles, which mirrors realistic dust charging behaviour and dynamics: the position determines the torques (which change the orientation) and the discharge (which changes the charge distribution); the orientation determines the forces (which change the position) and the discharge (which changes the charge distribution); and the charge distribution determines the forces (which change the position) and torques (which change the orientation).

2.5 Comparing simulation with experiment

Efforts were made to match many of the parameters of the simulation with the conditions of the experiment. The same forms and coefficients used to fit the Langmuir probe data for the floating potential ($V_f$), ion density ($n_i$) and electron temperature ($T_e$) in the analysis of the photo-discharging experiment were used as the profiles of the simulated plasma. The functions listed below were constructed such that the electric field, ion density and electron temperature, and their respective spatial derivatives, are all continuous:

(2.19)\begin{gather} V(z) = \left\{\begin{array}{@{}lr@{}} c_1 (z-c_5)^3 + c_2 (z-c_5)^2 + c_3 (z-c_5) + c_4 & z \leq c_5 \\ c_2 (z-c_5)^2 + c_3 (z-c_5) + c_4 & z > c_5 \end{array}\right\}, \end{gather}

(2.20)\begin{gather}n_i(z) = \left\{\begin{array}{@{}lr@{}} c_1 (\tan^{{-}1}(c_2 (z-c_4))+c_3) & z \leq c_4 \\ c_1 (c_2 (z-c_4) + c_3) & z > c_4 \end{array}\right\}, \end{gather}

(2.21)\begin{gather}T_e(z) = c_1 z + c_2. \end{gather}

The coefficients fitted to the Langmuir probe data and used in the simulation are listed in table 2, and plots of the experimental data and the models can be seen in figure 6.

Table 2. Coefficients for the functional forms of the floating potential (in V), ion density (in $10^{15}\,{\rm m}^{-3}$) and electron temperature (in eV) shown in (2.19), (2.20) and (2.21), after being fitted to data from the Langmuir probes. All coefficients shown above assume a vertical position, $z$, measured in cm.

Figure 6. Plots of the Langmuir probe data and (a) fits of the floating potential, (b) ion density and (c) electron temperature from the photo-discharging experiment. To match the simulation as closely to the experiment as possible, the same functions fitted to the Langmuir probe data from the experiment were used as the plasma profiles for the simulation.

Steps were also taken to try and match some of the characteristics of the dust used in the experiment. Using the SEM images of the LaB$_6$ samples from the experiment, measurements of the major and minor diameters of the particles were made. The distribution of axial diameters for the LaB$_6$ particles is a roughly Maxwellian distribution with a mean axial diameter of $10.77\pm 5.42\,\mathrm {\mu }$m. A distribution with these parameters was used to generate an initial sample of random axial diameters for the simulated particles.

While this initial sample of randomly generated axial diameters closely matched the distribution of the experimental samples, it did not generate a significant number of particles with a high degree of symmetry (nearly spherical or spheroidal particles), for which our hypothesis predicts more periodic motion. Since these particles are of particular interest to this investigation, a subset of nearly spherical and spheroidal particles were also included in the simulated sample (including the particles described in table 1). In total, 2500 ellipsoidal particles were generated for the simulated sample; $\sim$75 % of which had fully random axial lengths.

Figure 7 shows the distribution of the axial diameters of the LaB$_6$ particles, based on the measurements made using the SEM images, and those of the simulated particles. The mean axial diameter of the simulated particles was $14.56\pm 8.02\,\mathrm {\mu }$m. The figure also shows the distribution of effective diameters, defined here as the geometric mean of the axial diameters for each experimental or simulated particle. This effective diameter can be conceptualized as the diameter of a spherical particle with an equivalent cross-section or volume. The mean effective diameter of the particles was $10.45\pm 4.66\,\mathrm {\mu }$m in the experiment compared with $13.41\pm 5.93\,\mathrm {\mu }$m in the simulation.

Figure 7. Distribution of dust sizes for the particles in the experiment and the simulation and the overlap between distributions. (a) Distribution of all the measured and simulated axial diameters. The mean axial diameter of dust in the experiment was $10.77\pm 5.42\,\mathrm {\mu }$m, while the mean axial diameter of the simulated particles was $14.56\pm 8.02\,\mathrm {\mu }$m. (b) Distribution of the effective diameters of the real and simulated particles; defined here as the geometric mean of axial diameters (the diameter of a spherical particle with an equivalent cross-section or volume). The mean effective diameter of the particles was $10.45\pm 4.66\,\mathrm {\mu }$m in the experiment compared with $13.41\pm 5.93\,\mathrm {\mu }$m in the simulation.

2.6 Simulation procedure

Along with the values for axial lengths, random initial positions, velocities and orientations were also generated for each of the simulated particles. The parameters for each ellipsoid were used to generate 500 randomly distributed points, which form the surface of the simulated particle and are used to numerically calculate the surface area and $\varSigma {r_s}^2$.

In the iterative part of the simulation, the positions of the surface points are used to calculate $A_x$ and $A_z$. A normal vector, $\hat {\boldsymbol {n}}$, is also calculated for each surface point and the value of the inner product, $\hat {\boldsymbol {x}}\boldsymbol {\cdot }\hat {\boldsymbol {n}}$, is calculated for each point to determine if that point faces towards or away from the light source. The centroid of the illuminated half of the surface is calculated, along with the instantaneous charge and dipole moment.

The position of the particle is used to calculate the local plasma properties, and the forces and torques acting on the particle. Velocity, angular velocity, position and orientation are updated using the implicit Euler method. These properties are fed back into the next iteration for a simulated time of 24 seconds. The first five seconds of the simulation serve as a settling time to account for any fluctuations due to the random starting conditions, after which data are recorded for a total of 19 simulated pulses.

3 Results

The simulation outputs a significant amount of data, but for this investigation, the focus is on the particle's position and velocity over the duration of the simulated pulses.

3.1 Individual results

First, to provide context for the rest of the data, we consider the results of the simulation for the four particles from figure 5: A sphere, prolate and oblate spheroids, and a scalene ellipsoid. Figure 8 shows position versus pulse-time for the four particles over the course of the 19 simulated pulses. As was the case in figure 3, the colour and shading of the data points are used to indicate the passage of time, becoming darker with each pulse.

Figure 8. Vertical position of the particles from figure 5 – the sphere (a), prolate (b) and oblate (c) spheroids, and the scalene ellipsoid (d) – over the course of each simulated pulse (with the pulse turning on at 0.0 s, and off at 0.5 s). The colour and brightness are used to indicate the passage of time, with the data for each pulse growing steadily darker. The motion of the spherical particle is virtually identical for each pulse, while the motions of the less symmetric particles diverge.

Figure 9 shows the motion of the same four particles in phase space: the spherical particle has a tight, closed orbit with no deviation. The prolate and oblate spheroids are nearly closed orbits, with minor to moderate deviation. The scalene ellipsoid exhibits major deviations and its orbits are open, although the particle still returns to the same broad region in phase space at the end of each period. Qualitatively, we can describe the sphere's behaviour as clearly periodic, the scalene ellipsoid's behaviour as clearly chaotic, and the prolate and oblate spheroids’ behaviours as falling somewhere in between.

Figure 9. Motion of the particles from figure 5 in phase space. The colour and brightness of the data are used to indicate the passage of time, with each pulse growing darker. The orbits of the sphere in panel (a) are closed and nearly identical; those of the spheroids in panels (b,c) are closed or nearly closed, but with moderate variation; and those of the scalene ellipsoid in panel (d) are slightly open and extremely varied.

Quantitatively, we can describe the periodicity of each particle's motion in terms of the deviation in its area in phase space, $A$, where, for a continuous trajectory over a closed orbit, $A$ would be given by

(3.1)\begin{equation} A = \oint \boldsymbol{v_z}\boldsymbol{\cdot}\boldsymbol{{d} z}. \end{equation}

To approximate this integral, the phase space data are divided into one-second pulses and the midpoint method is used to estimate the phase space area for each pulse. The mean and standard deviations of the area, $\mu _A$ and $\sigma _A$, respectively, are then calculated from the results for all 19 simulated pulses. We consider a parameter, $\delta _A$, defined as the standard deviation of the phase space area divided by the mean:

(3.2)\begin{equation} \delta_A=\frac{\sigma_A}{\mu_A}. \end{equation}

The smaller values of $\delta _A$ will indicate a smaller deviation in the shape of the particle's phase space orbit over the recording period and, therefore, more periodic motion, while higher values indicate more chaotic motion. Furthermore, when converted to a logarithmic scale, this effectively describes a crude signal-to-noise ratio (SNR) in the particle's motion:

(3.3)\begin{equation} {\rm SNR}={-}\ln{\delta_A}. \end{equation}

Table 3 shows the linear and logarithmic values for the motion of the four particles shown in figure 9 and the two experimental particles shown in figure 3. Quantitatively, the sphere had the smallest numerical deviation in phase space area and can be considered the most periodic behaviour of the any of the particles (simulated or not), while the scalene ellipsoid had the greatest numerical deviation in phase space area and can be considered the most chaotic behaviour.

Table 3. Values of $\delta _A$ and SNR for the four simulated ellipsoids described in figure 9, and the two particles from the experiment described in figure 3, ordered from most periodic to least.

3.2 Ensemble results

Calculations of $\delta _A$ were made for 2500 simulated particles. Figure 10 shows the resulting data for the particles in both the real-space, described by the axial lengths of the simulated ellipsoids, and the same data projected into the $\kappa$-space we defined earlier. The colour of the data points indicates the value of $\delta _A$, with lighter points indicating more periodic behaviour and darker points indicating more chaotic behaviour. Most of the periodic points are concentrated in the regions of the plots corresponding to spherical particles.

Figure 10. Simulated particles viewed in (a) the real-space determined by the axial lengths of the ellipsoids and (b) the logarithmic $\kappa$-space. The transformation to $\kappa$-space projects the particles into a plane, which preserves information about particle shape at the expense of information about the physical scale. The colour and brightness of the points indicate the value of $\delta _A$, with lighter points corresponding to more periodic motion and darker points to more chaotic motion. (c,d) Same as panels (a,b) when looking ‘down’ through the centre of the panels (along the line defining spherical particles). The three axes in the $\kappa$-space in panel (d) define prolate spheroids, with each tick indicating a factor of $\kappa _0$.

The data projected into the $\kappa$-plane were divided into a 25-by-25 grid and the average values of $\delta _A$ were calculated for each interrogation region. The final results for $\delta _A$ and SNR were rotated into the two-dimensional heat maps in figure 11. As in figure 10, the brightest points in the grid indicate more periodic motion and the darker points in the grid indicate more chaotic motion. The greatest values of SNR are concentrated near the origin of $\kappa$-space (${\rm SNR}\gtrsim 3$). Moderate SNR values are clustered along the positive $\kappa$-axes, and minor SNR values appear to fill out the negative and off-axis space. The lowest SNR values (the most chaotic motion) appear near the fringes of our data set as it appears in $\kappa$-space (the black regions indicate values of ${\rm SNR} \lesssim 0$).

Figure 11. (a) Results for $\delta _A$ and (b) the signal-to-noise ratio averaged over a grid of interrogation cells in $\kappa$-space (note that this is the same data as figure 10, just looking out away from the $\kappa$ plane instead of down towards it). In both panels, colour and shading are used to indicate periodicity, with lighter colours corresponding to more periodic behaviour and darker colours to more chaotic behaviour. Regions where $\delta _A$ exceeds unity are coloured black in both panels. The highest periodicity is located at the origin, with behaviour becoming less periodic and more chaotic as particles move away from the centre, and towards more scalene and oblate particles. The rings indicate circles of constant-$|\kappa |$ ($\kappa _0$, 2$\kappa _0$ and 3$\kappa _0$).

4 Discussion and summary

This computational investigation set out to answer two questions. The first was whether a computational model, in which the dynamics and charging behaviour of a freely rotating particle depended on the particle's area, could reproduce the different types of periodic and chaotic motion seen in the Auburn experiment. The second was whether the same model would indicate a correlation between a particle's shape and the periodicity of its behaviour in response to the photo-discharging.

4.1 Reproduction of experimental results

While the behaviours of the example particles shown in figure 9 are not identical to those of the experimental particles in figure 3, our computational model was successful in reproducing the broad types of behaviours seen in the experiment: periodic-like behaviour, where particles follow closed or nearly closed orbits through phase space, with minor deviations, and chaotic-like behaviour, where particles experience large deviations in their motion through phase space.

4.2 Symmetry–periodicity correlation

Taking the individual results from table 3 and the ensemble results from figure 11 into consideration, there appears to be a correlation between particle shape and dynamic periodicity under the constraints of our computational model, which can be divided into four, broadly defined behavioural regions, summarized in table 4: strongly periodic behaviour describes tight, closed orbits with little or no deviation in phase space area ($\delta _A < e^{-3}$); weakly periodic behaviour describes minor deviations in phase space area and are qualitatively judged to be more periodic than chaotic ($\delta _A < e^{-1.5}$); weakly chaotic behaviour describes moderate deviations in phase space area and are qualitatively judged to be more chaotic than periodic ($\delta _A < 1$); strongly chaotic behaviour describes major deviations in phase space area, where $\delta _A$ approaches or exceeds unity.

Table 4. Descriptions of behavioural regimes of simulated particle motion correlated with approximate ranges of SNR. Behaviour becomes more periodic as SNR increases and more chaotic as it decreases. Strongly periodic behaviour describes cases where ${\rm SNR} \gtrsim 3$. Strongly chaotic behaviour describes cases where SNR drops below 0 (regions where $\delta _A$ is over unity).

Figure 12 separates the results of figure 11 into the four regimes described in table 4 to more clearly illustrate the divide in dynamic behaviour: Strongly periodic behaviour is limited to spherical particles at the origin; weakly periodic behaviour appears concentrated within a radius of ${\sim }\kappa _0$ from the origin and appears to favour prolate particles; weakly chaotic behaviour appears concentrated between $\sim 1$ and $2\kappa _0$, and appears to favour oblate and scalene particles; strongly chaotic behaviour appears concentrated beyond $\sim$2$\kappa _0$, and appears to favour oblate and scalene particles. A slight asymmetry is visible along the $\hat {\kappa }_{yz}$-axis, most notably in the weakly periodic and strongly chaotic regimes, which may be the result of an oversampling of spheroidal particles oriented along that particular axis.

Figure 12. Same data as figure 11, but divided into the behavioural regimes described in table 4: strongly periodic (SP) behaviour is limited to spherical particles at the centre of the grid; weakly periodic (WP) behaviour is mostly confined to within a radius of ${\sim }\kappa _0$ from the centre and appears to indicate a slight preference towards prolate particles; weakly chaotic (WC) behaviour dominates the region between $\sim 1$ and $2\kappa _0$; strongly chaotic (SC) behaviour seems to be limited to particles beyond $\sim 2\kappa _0$. A slight asymmetry can be seen along the $\hat {\kappa }_{yz}$-axis; this asymmetry may be the result of an oversampling of spheroidal particles along this axis.

The initial results of the investigation, best summarized by figure 11 and figure 12, appear to support the hypothesis that when particle dynamics and discharge are dependent on the projected area of the particle, the feedback loop between dynamics and charging behaviour can drive more periodic or more chaotic behaviours, depending on the particle's shape.

4.3 Future work

While these initial results are promising, further computational investigation and experimental follow-up should be considered.

The computational model could be further refined and more particles could be simulated to generate more statistical data across a wider range of parameters; this could reveal more detailed structure in the distribution of the different behavioural regimes. The model could also be expanded to consider more irregular particle shapes, sharp edges, concavities, etc. More complex treatments of the charge distribution beyond the simple dipole model could be considered. The possibility of analytic or quasi-analytic approximations to parts of the computational model could be investigated as well.

While the model presented here succeeds in reproducing and explaining the observations seen in the experiments using the LaB$_6$, further experimental work will be necessary to fully test this hypothesis. Since the individual LaB$_6$ particles cannot be extracted and analysed under the SEM to directly compare the dynamics and geometry of each particle, the next best option would be to conduct an experimental survey involving the discharge of spherical particles; either microspheres coated in LaB$_6$ or, preferably, composed of another more readily available material which also reacts to NUV exposure. In the ideal experimental scenario, examples of both spherical and non-spherical particles of the same material could be subjected to photo-discharging to see if there is a clear distinction between the behaviour of one sample versus the other.