2. Basic equation

In this paper, we consider a plasma consisting of electrons, singly charged ions, and negatively charged dust particles of the same size in a uniform magnetic field $\boldsymbol {B}_{0} = B_{0} \boldsymbol {k}$, where $\boldsymbol {k}$ is the unit vector in the $z$-direction in the Cartesian coordinate system. We consider a case in which the propagation direction of the wave is in the $x$-direction, perpendicular to the direction of the magnetic field. Furthermore, for simplicity, we assume that the charged dust stream is in the same direction with speed ${u}_{d0}$. In an equilibrium state, both electrons and ions satisfy the Maxwell distribution, whereas the dust particles satisfy a drifting Maxwell distribution.

We now use the magnetohydrodynamic equations to study our system. The basic equations are as follows:

(2.1)\begin{gather} \frac{\partial n_{d}}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot}(n_{d} \boldsymbol {u}_{d})=0 \end{gather}

(2.2)\begin{gather}n_{d} m_{d}\left(\frac{\partial}{\partial t} + \boldsymbol{u}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla}\right)\boldsymbol{u}_{d} ={-}\boldsymbol{\nabla} p_{d} + n_{d} Q_{d} \boldsymbol{E} \end{gather}

(2.3)\begin{gather}\frac{\partial n_{i}}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot}(n_{i} \boldsymbol{u}_{i})=0 \end{gather}

(2.4)\begin{gather}n_{i} m_{i}\left(\frac{\partial}{\partial t} + \boldsymbol{u}_{i} \boldsymbol{\cdot} \boldsymbol{\nabla}\right)\boldsymbol{u}_{i} ={-}\boldsymbol{\nabla} p_{i} + n_{i} Q_{i} \left(\boldsymbol{E}+\frac{\boldsymbol{u}_{i}\times \boldsymbol{B}}{c}\right) \end{gather}

(2.5)\begin{gather}\frac{\partial n_{e}}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (n_{e} \boldsymbol{u}_{e})=0\end{gather}

(2.6)\begin{gather}n_{e} m_{e}\left(\frac{\partial}{\partial t} + \boldsymbol{u}_{e} \boldsymbol{\cdot} \boldsymbol{\nabla}\right)\boldsymbol{u}_{e} ={-}\boldsymbol{\nabla} p_{e} + n_{e} Q_{e} \left(\boldsymbol{E}+\frac{\boldsymbol{u}_{e}\times \boldsymbol{B}}{c}\right) \end{gather}

(2.7)\begin{gather}\nabla^{2}\phi ={-}\frac{1}{\varepsilon_{0}}(n_{d}Q_{d}+n_{i}Q_{i}+n_{e}Q_{e}), \end{gather}

where $n_{\alpha }$, $\boldsymbol u_{\alpha }$, $m_{\alpha }$ and $Q_{\alpha }$ are the number density, velocity, mass and charge of the species $\alpha$ ($\alpha =d,i,e$, denote the dust particles, ions and electrons), respectively. Here $p_{\alpha }$ is the pressure of species $\alpha$ ($\alpha =d,i,e$), and ${\boldsymbol E}$, ${\boldsymbol {B}}$, $c$ and $\phi$ are the electric field, magnetic field, light velocity and electric potential, respectively. We use $\varepsilon _{0}$ to denote the permittivity of the vacuum. The charges on the dust particles are assumed to be constant due to the charging time of the dust particles being far less than the hydrodynamic time of the dust particles, i.e. $Q_{d}={\rm const}$.

3. Dispersion relation

In order to obtain the dispersion relation of the system, we assume that all the physical quantities vary in the following forms: $f=f_{0}+f_{1}\exp ({{\rm i}(\boldsymbol {k}\boldsymbol {\cdot } \boldsymbol r-\omega t)})$, where $f_{0}$ is the value at its equilibrium state, $f_{1}$ is its small perturbation, $\boldsymbol {k}$ and $\omega$ are the wavenumber and frequency of the perturbation wave, respectively. We assume that $\boldsymbol {k}$ is in the $xoz$ plane, i.e. a two-dimensional system, and the angle with the $z$-axis is $\theta$, i.e. $\boldsymbol {k} =k_{\bot }\hat {x}+k_{\|}\hat {z}=k \sin \theta \hat {x}+k \cos \theta \hat {z}$, where $k=\mid {\boldsymbol {k}}\mid$. Then, we have from (2.1)–(2.7) the following dispersion relation:

(3.1)\begin{align} & 1-\frac{\omega_{{\rm pe}}^{2}\left(1+ \dfrac{\varOmega_{{\rm ce}}^{2}\sin^{2}\theta}{ \omega^{2}-\varOmega_{{\rm ce}}^{2}}\right)}{\omega^{2}-k^{2}\upsilon_{e}^{2} \left(1+\dfrac{\varOmega_{{\rm ce}}^{2}\sin^{2}\theta} {\omega^{2}-\varOmega_{{\rm ce}}^{2}}\right)}- \frac{\omega_{{\rm pi}}^{2}\left(1+ \dfrac{\varOmega_{{\rm ci}}^{2}\sin^{2}\theta}{\omega^{2}-\varOmega_{{\rm ci}}^{2}}\right)} {\omega^{2}-k^{2}\upsilon_{i}^{2}\left(1+\dfrac{\varOmega_{{\rm ci}}^{2}\sin^{2}\theta} {\omega^{2}- \varOmega_{{\rm ci}}^{2}}\right)}\nonumber\\ & \quad -\frac{\omega_{{\rm pd}}^{2}}{(\omega-k \sin\theta u_{d0})^{2}-k^{2}\upsilon_{d}^{2}}=0, \end{align}

where $\omega _{p\alpha }$, $\varOmega _{c\alpha }$ and $\upsilon _{\alpha }$ are the oscillation frequency, cyclotron frequency and the thermal velocity of the species $\alpha$ ($\alpha =d, i, e$), respectively, $\omega _{p\alpha }=({n_{\alpha }e^{2}}/{\varepsilon _{0}m_{\alpha }})^{1/2}$, $\varOmega _{c\alpha }=({Q_{\alpha }B_{0}}/{m_{\alpha }})^{1/2}$ and $\upsilon _{\alpha }=({\gamma T_{\alpha }}/{m_{\alpha }})^{1/2}$. Here $T_{\alpha }$ is the temperature of the species $\alpha$ ($\alpha =d, i, e$).

Equation (3.1) is the dispersion relation of the linear wave for our system. In order to further understand it, we consider three special cases in the following because the masses of electrons, ions and dust particles in a dusty plasma are different widely, thus their oscillation frequencies are much different.

3.1. High-frequency case

First, we consider the high-frequency case: $\omega ^{2}>\omega _{{\rm pe}}^{2}$, $\varOmega _{{\rm ce}}^{2}\gg \omega _{{\rm pi}}^{2}>\varOmega _{{\rm ci}}^{2}\gg \omega _{{\rm pd}}^{2}$, then the dispersion relation becomes as follows:

(3.2)\begin{equation} \omega^{2}=(\omega_{{\rm pe}}^{2}+k^{2}\upsilon_{e}^{2}) \left(1+\frac{\varOmega_{{\rm ce}}^{2}\sin^{2}\theta}{\omega^{2}-\varOmega_{{\rm ce}}^{2}}\right). \end{equation}

If the waves propagate in the direction parallel to the magnetic field, we have the dispersion relation of the Langmuir waves (LWs)

(3.3)\begin{equation} \omega^{2}=\omega_{{\rm LM}}^{2}=\omega_{{\rm pe}}^{2}+k^{2}\upsilon_{e}^{2}. \end{equation}

If the waves propagate in the direction perpendicular to the magnetic field, we have the dispersion relation of the upper hybrid waves (UHWs)

(3.4)\begin{equation} \omega^{2}=\omega_{{\rm UH}}^{2}=\omega_{{\rm pe}}^{2}+\varOmega_{{\rm ce}}^{2}+k^{2}\upsilon_{e}^{2}. \end{equation}

3.2. Medium-frequency case

Second, we consider the medium-frequency case: $\omega _{{\rm pe}}^{2}$, $\varOmega _{{\rm ce}}^{2}\gg \omega ^{2}>\omega _{{\rm pi}}^{2}>\varOmega _{{\rm ci}}^{2}\gg \omega _{{\rm pd}}^{2}$, then the dispersion relation becomes as follows:

(3.5)\begin{equation} 1-\frac{\omega_{{\rm pe}}^{2}}{k^{2}\upsilon_{e}^{2}}- \frac{\omega_{{\rm pi}}^{2}\left(1+ \dfrac{\varOmega_{{\rm ci}}^{2}\sin^{2}\theta}{\omega^{2}-\varOmega_{{\rm ci}}^{2}}\right)} {\omega^{2}-k^{2}\upsilon_{i}^{2}\left(1+\dfrac{\varOmega_{{\rm ci}}^{2}\sin^{2}\theta} {\omega^{2}-\varOmega_{{\rm ci}}^{2}}\right)}=0. \end{equation}

If the waves propagate in the direction parallel to the magnetic field, we have the dispersion relation of the dust ion acoustic waves (DIAWs)

(3.6)\begin{equation} \omega^{2}=\omega_{{\rm DIA}}^{2}=\frac{\omega_{{\rm pi}}^{2}k^{2}\upsilon_{e}^{2}} {k^{2}\upsilon_{e}^{2}+\omega_{{\rm pe}}^{2}}+k^{2}\upsilon_{i}^{2}. \end{equation}

If the waves propagate in the direction perpendicular to the magnetic field, we have the dispersion relation of the dust ion cyclotron waves (DICWs)

(3.7)\begin{equation} \omega^{2}=\omega_{{\rm DIC}}^{2}=\frac{\omega_{{\rm pi}}^{2}k^{2}\upsilon_{e}^{2}} {k^{2}\upsilon_{e}^{2}+\omega_{{\rm pe}}^{2}}+k^{2}\upsilon_{i}^{2}+\varOmega_{{\rm ci}}^{2}. \end{equation}

3.3. Low-frequency case

Finally, we consider the low-frequency case: $\omega _{{\rm pe}}^{2}$, $\varOmega _{{\rm ce}}^{2}\gg \omega _{{\rm pi}}^{2}>\varOmega _{{\rm ci}}^{2}\gg \omega _{{\rm pd}}^{2}>\omega ^{2}$, the dispersion relation becomes as follows:

(3.8)\begin{equation} 1-\frac{\omega_{{\rm pe}}^{2}}{k^{2}\upsilon_{e}^{2}}- \frac{\omega_{{\rm pi}}^{2}}{k^{2}\upsilon_{i}^{2}}- \frac{\omega_{{\rm pd}}^{2}}{(\omega-k \sin\theta u_{d0})^{2}-k^{2}\upsilon_{d}^{2}}=0. \end{equation}

If the waves propagate in the direction parallel to the magnetic field, we have the following dispersion relation of the dust acoustic waves (DAWs)

(3.9)\begin{equation} \omega^{2}=\omega_{{\rm DA}}^{2}=\frac{k^{2}\omega_{{\rm pd}}^{2}}{k^{2}+K_{{\rm De}}^{2}+K_{{\rm Di}}^{2}}+k^{2}\upsilon_{d}^{2}. \end{equation}

If the waves propagate in the direction perpendicular to the magnetic field, we have

(3.10)\begin{equation} \omega=\omega_{{\rm DAU}}=\left(\frac{k^{2}\omega_{{\rm pd}}^{2}}{k^{2}+K_{{\rm De}}^{2}+K_{{\rm Di}}^{2}} +k^{2}\upsilon_{d}^{2}\right)^{{1}/{2}}+ku_{d0}, \end{equation}

where $K_{{\rm De}}=1/\lambda _{{\rm De}}$, $K_{{\rm Di}}=1/\lambda _{{\rm Di}}$, $\lambda _{{\rm De}}$ and $\lambda _{{\rm Di}}$ are the Debye length of the electrons and ions, respectively.

Figure 1 shows the dispersion relations given by (3.1) for our system. The dispersion relations on three different time scales are also shown in figure 1. Note that they are in agreement for small-wavenumber cases.

Figure 1. The dispersion relations of the waves excited by the relative streaming of the dust fluid to the background plasma, where the blue asterisks represent theoretical solutions of (3.1), the red solid lines represent the dispersion relations of the waves on three different time scales, and the green dots in (d) and (e) represent the particle-in-cell (PIC) simulation results on electron and ion vibration time scales, respectively.

4. Particle-in-cell simulation

We now investigate waves and instabilities driven by a charged dust beam in the ionosphere by using the particle-in-cell (PIC) simulation method. The initial condition is the equilibrium state of the system. The periodic boundary condition is used. We only study the motion of the dusty plasma in the $xoz$ plane, i.e. a two-dimensional system. The size of the simulation area is $1902\ {\rm mm} \times 951\ {\rm mm}$, and $256 \times 128$ grid nodes are chosen. The equations of motion of the super particles (SPs) are the Newton equations as follows:

(4.1)\begin{gather} m_{\alpha}\frac{{\rm d} \boldsymbol{\upsilon}_{\alpha}}{{\rm d}t} =q_{\alpha}(\boldsymbol{E}+ \boldsymbol{\upsilon}_{\alpha}\times \boldsymbol{B}), \end{gather}

(4.2)\begin{gather}\frac{{\rm d}\boldsymbol{r}_{\alpha}}{{\rm d}t}=\boldsymbol{\upsilon}_{\alpha}, \end{gather}

where $m_{\alpha }$, $q_{\alpha }$, $\boldsymbol {\upsilon }_{\alpha }$ and $\boldsymbol {r}_{\alpha }$ are the mass, charge, velocity and position of the species $\alpha$ ($\alpha =d,i,e$), respectively. Here $\boldsymbol {E}$ is the electric field and $\boldsymbol {B}$ is the external magnetic field. In the PIC simulation, the simulation area is divided into grid cells. At each time step, the velocities and the positions of SPs are weighted to all the grids to calculate the charge density $\rho _{g}$. Once $\rho _{g}$ is obtained, the Poisson equation (electrostatic model) will be used to solve the potential of each grid, and the value of $\boldsymbol {E}$ is further derived. Then, each SP will be driven by the electric field according to (4.1) and (4.2), which will be solved numerically via the leap-frog algorithm. Finally, the new positions and velocities are obtained, the procedure repeats until the simulation is completed. As the masses of electrons, ions and dust particles vary widely, we choose different time steps and time scales for the simulations. The simulation is carried out for collisional and collisionless plasma, respectively.

First, on the electron vibration time scale, electrons, ions and dust particles in the plasma are treated as SPs. As the ion mass is much larger than the electron mass ($m_{i}/m_{e}\geqslant 1836$), as a first-order approximation, we assume that ions do not move, and only form a uniform background in space. The dust particles also do not move. We only study the oscillation of electrons. The parameters in the simulation are selected as follows: the space step is ${{\rm d} x}={\rm d}z=\lambda _{{\rm Di}}$, $\lambda _{{\rm Di}}\sim 7.43\times 10^{-3}\ {\rm m}$, the time step is ${\rm d}t=2.0\times 10^{-10}\ {\rm s}$.

Second, on the ion vibration time scale, electrons, ions and dust particles in the plasma are treated as SPs. The parameters in the simulation are selected as follows: the space step is ${{\rm d} x}={\rm d}z=\lambda _{{\rm Di}}$, $\lambda _{{\rm Di}}\sim 7.43\times 10^{-3}\ {\rm m}$, the time step is ${\rm d}t=2.0\times 10^{-8}\ {\rm s}$.

Third, on the dust vibration time scale, dust particles are regarded as SPs with macroscopic drift velocity, whereas electrons and ions are modelled as a Boltzmann distributed background. The parameters in the simulation are selected as follows: the space step is ${{\rm d} x}={\rm d}z=\lambda _{{\rm Di}}$, $\lambda _{{\rm Di}}\sim 7.43\times 10^{-3}\ {\rm m}$, the time step is ${\rm d}t=2.0\times 10^{-6}\ {\rm s}$.

In addition, we have considered the collision effects between particles on the waves in the system. We use the Monte Carlo method to study the collision between ions and neutrals, dust particles and neutrals. We neglect the collision between electrons and neutrals because the collisional frequency is so small compared with the electron oscillation frequency that can be ignored.

5. Numerical simulation results

We consider the F-region of the ionosphere at $240$ km altitude. The primary ion in the F-region is $O^{+}$. The parameters chosen in the simulation are as follows: $B \sim 5 \times 10^{-5}\ {\rm T}$, $n_{e} \sim 0.994 \times 10^{5}\ {\rm cm}^{-3}$, $n_{i} \sim 1.0 \times 10^{5}\ {\rm cm}^{-3}$, $n_{d} \sim 0.006n_{i}$, $T_{e} \sim 0.2\ {\rm eV}$, $T_{d} = T_{i} \sim 0.1\ {\rm eV}$, $Z_{d}\sim 50$ and $m_{d} \sim 4\times 10^{9}m_p$, where $m_p$ is the proton mass (Mendis & Rosenberg Reference Mendis and Rosenberg1994). The dust particles are injected into the background plasma with a speed $u_{d0} \sim 2\ {\rm km}\ {\rm s}^{-1}$ (Huba et al. Reference Huba, Joyce and Fedder2000; Rosenberg & Sorasio Reference Rosenberg and Sorasio2006). Then we obtain $\omega _{{\rm pe}}=1.50\times 10^{7}\ {\rm rad}\ {\rm s}^{-1}$, $\omega _{{\rm pi}}=0.93\times 10^{5} \ {\rm rad}\ {\rm s}^{-1}$, $\omega _{{\rm pd}}=25.50\ {\rm rad}\ {\rm s}^{-1}$, $\varOmega _{{\rm ce}}=8.79\times 10^{6} \ {\rm rad}\ {\rm s}^{-1}$, $\varOmega _{{\rm ci}}=0.24\times 10^{3}\ {\rm rad}\ {\rm s}^{-1}$ and $\varOmega _{{\rm cd}}=5.99 \times 10^{-5}\ {\rm rad}\ {\rm s}^{-1}$.

Figure 2 shows the numerical results of the waves excited by the charged dust beam on the electron vibration time scale. The dependence of the electrostatic potential on time at the fixed point for collisionless dusty plasma is shown in figure 2(a), and its corresponding frequency spectrum analysis is given by figure 2(b). It seems from figures 2(a) and 2(b) that there is only one mode of the electron vibration and its frequency is about $\omega _e \sim 1.46\times 10^{7}\ {\rm rad}\ {\rm s}^{-1}$. In order to know what kind of wave is excited by this electron vibration, the variation of the electrostatic potential with the position is shown in figure 2(c), and its corresponding wavenumber spectrum analysis is given by figure 2(d). Note from figures 2(c) and 2(d) that the wavenumber is about $k_e \sim 3.28979\ {\rm m}^{-1}$. It is inferred from our numerical results of both wave frequency and wavenumber that the UHW is excited in the system by comparing our numerical results with the theoretical results in figure 1.

Figure 2. Numerical results of the electron vibration for (a)–(d) collisionless dusty plasma and (e)–(h) collisional dusty plasma. The dependence of the electrostatic potential on time is shown in (a) and (e) and the corresponding frequency spectrum is shown in (b) and ( f). The variation of the electrostatic potential with respect to the position is shown in (c) and (g) and the corresponding frequency spectrum is shown in (d) and (h).

For collisional dusty plasma, the dependence of the electrostatic potential on time at the fixed point is shown in figure 2(e), and its corresponding spectrum analysis is given by figure 2( f). The variation of the electrostatic potential with the position is shown in figure 2(g), and its corresponding wavenumber spectrum analysis is given by figure 2(h). Numerical results show that the effects of the collision between electrons and neutrals can be neglected because the collision frequency between electrons and neutrals is much smaller than the vibration frequency of the electrons.

In order to investigate whether this wave is stable, the trajectory and phase diagram of a single SP (electron) are given in figure 3. It seems that the wave is quasi-periodic. Therefore, it is stable.

Figure 3. (a) The trajectory of one of the SPs in the $xoz$ plane. (b) The phase diagram of the SP in phase space of $(x,\upsilon _{x})$.

Figure 4 shows the numerical results of the waves excited by the charged dust beam on the ion vibration time scale. The dependence of the electrostatic potential on time at the fixed point for collisionless dusty plasma is shown in figure 4(a), and its corresponding spectrum analysis is given by figure 4(b). It seems from figures 4(a) and 4(b) that there is only one mode of the ion vibration, and its frequency is about $\omega _{i} \sim 4.17\times 10^{4}\ {\rm rad}\ {\rm s}^{-1}$. In order to know what kind of wave is excited by this ion vibration, the variation of the electrostatic potential with the position is shown in figure 4(c), and its corresponding spectrum analysis is given by figure 4(d). Note from figures 4(c) and 4(d) that the wavenumber is about $k\sim 9.87\ {\rm m}^{-1}$. It is inferred from our numerical results of both wave frequency and wavenumber that the DICW is excited in the system by comparing our numerical results with the theoretical results in figure 1.

Figure 4. Numerical results of the ion vibration for (a)–(d) collisionless dusty plasma and (e)–(h) collisional dusty plasma. The dependence of the electrostatic potential on time is shown in (a) and (e) and the corresponding frequency spectrum is shown in (b) and ( f). The variation of the electrostatic potential with respect to the position is shown in (c) and (g) and the corresponding frequency spectrum is shown in (d) and (h).

For collisional dusty plasma, the dependence of the electrostatic potential on time at the fixed point is shown in figure 4(e), and its corresponding spectrum analysis is given by figure 4( f). The variation of the electrostatic potential with the position is shown in figure 4(g), and its corresponding spectrum analysis is given by figure 4(h). It seems from the numerical results that the collisional effect on the ion vibration time scale is so small that it is negligible.

In order to know whether the ion vibration is stable or not, we give the Lyapunov index in table 1. It is noted that the Lyapunov indexes are all negative. Therefore, we conclude that the wave on the ion vibration time scale may be stable no matter the collisions between particles are considered or not.

Table 1. The Lyapunov index on the ion vibration time scale.

Figure 5 shows the numerical results of the waves excited by the charged dust beam on the dust vibration time scale. The dependence of the electrostatic potential on time at the fixed point for collisionless dusty plasma is shown in figure 5(a), and its corresponding spectrum analysis is given by figure 5(b). For the collisional dusty plasma, the variation of the electrostatic potential with the time at the fixed point is shown in figure 5(c), and its corresponding spectrum analysis is given by figure 5(d).

Figure 5. Numerical results of the vibration of the dust particles for (a) and (b) collisionless dusty plasma and (c) and (d) collisional dusty plasma. The dependence of the electrostatic potential on time is shown in (a) and (c), and the corresponding frequency spectrum is shown in (b) and (d).

Note from figures 5(a)–5(d) that the modes of the dust vibration are more complex. There is no certain frequency of the vibration of dust particles for both collisional and collisionless dusty plasma. Therefore, we now focus on whether the waves excited by the dust fluid beam on the dust particles vibration time scale are stable or not. For this purpose, the corresponding Lyapunov index of the dust vibration is listed in table 2. It is found that the Lyapunov indexes of the dust vibration always have positive values no matter the collisions between particles are considered or not. We then conclude that the vibrations on the dust particles time scale are unstable.

Table 2. The Lyapunov index on the dust particles vibration time scale.