1 Introduction
The sheath formation in the interface of a plasma with a wall is an important research area in plasma physics. The electrons have a thermal speed much higher than the ions and thus their flux received on the wall is larger than that of the ions. Then, the wall has negative potential with respect to the main body of the plasma. Due to the negative potential, the electrons are repelled from the wall while positive ions are attracted towards it, so a layer with net positive charge is formed near the wall, which is called the plasma sheath. In most cases, the plasma contains dust nanoparticles in addition to electrons and ions (Daugherty & Graves Reference Daugherty and Graves1995; Samarian et al. Reference Samarian, James, Vladimirov and Cramer2001; Samarian & Vladimirov Reference Samarian and Vladimirov2003; Samarian, Vladimirov & James Reference Samarian, Vladimirov and James2005; Adhikary et al. Reference Adhikary, Bailung, Pal, Chutia and Nakamura2007). The problem of sheath formation in a plasma in the absence (Godyak & Sternberg Reference Godyak and Sternberg1990; Riemann Reference Riemann1990, Reference Riemann1991; Sheridan & Goree Reference Sheridan and Goree1991; Sheridan Reference Sheridan2001; Franklin Reference Franklin2003; Sternberg & Godyak Reference Sternberg and Godyak2003; Franklin Reference Franklin2005; Riemann et al. Reference Riemann, Seebacher, Tskhakaya and Kuhn2005; Tskhakaya et al. Reference Tskhakaya, Shukla, Eliasson and Kuhn2005; Lee et al. Reference Lee, Severn, Oksuz and Hershkowitz2006) and presence of dust nanoparticles (Barnes et al. Reference Barnes, Keller, Forster, ONeill and Coultas1992; Ma, Liu & Yu Reference Ma, Liu and Yu1997; Tsytovich, Vladimirov & Benkadda Reference Tsytovich, Vladimirov and Benkadda1999; Resendes, Sorasio & Shukla Reference Resendes, Sorasio and Shukla2002; Davoudabadi & Mashayek Reference Davoudabadi and Mashayek2005; Pandey, Samarian & Vladimirov Reference Pandey, Samarian and Vladimirov2007; Foroutan, Mehdipour & Zahed Reference Foroutan, Mehdipour and Zahed2009; Foroutan Reference Foroutan2010; Foroutan & Akhoundi Reference Foroutan and Akhoundi2012a,Reference Foroutan and Akhoundib; Jalilpour & Foroutan Reference Jalilpour and Foroutan2017, Reference Jalilpour and Foroutan2018) has been extensively studied under different conditions. It was found that the presence of dust nanoparticles in the plasma modifies strongly the structure of the plasma sheath (Barnes et al. Reference Barnes, Keller, Forster, ONeill and Coultas1992; Ma et al. Reference Ma, Liu and Yu1997; Tsytovich et al. Reference Tsytovich, Vladimirov and Benkadda1999; Resendes et al. Reference Resendes, Sorasio and Shukla2002; Davoudabadi & Mashayek Reference Davoudabadi and Mashayek2005; Pandey et al. Reference Pandey, Samarian and Vladimirov2007; Foroutan et al. Reference Foroutan, Mehdipour and Zahed2009; Foroutan Reference Foroutan2010; Foroutan & Akhoundi Reference Foroutan and Akhoundi2012a,Reference Foroutan and Akhoundib; Jalilpour & Foroutan Reference Jalilpour and Foroutan2017, Reference Jalilpour and Foroutan2018).
In gas discharges the energy distribution of the electrons is very important in determining the plasma properties. It is well known that, when a system of particles is in complete thermal equilibrium, the distribution of electrons can be represented by the Maxwellian distribution. The electrons in a discharge are not in equilibrium and the electron energy distribution function (EEDF) in non-stationary plasmas shows usually a non-Maxwellian form, and it is necessary to use a suitable function to formulate these EEDFs. (Lieberman & Lichtenberg Reference Lieberman and Lichtenberg1994). Druyvesteyn and Penning have shown that, in most gas plasmas, the energy distribution has a Druyvesteyn type (Druyvesteyn & Penning Reference Druyvesteyn and Penning1940). During the discharge the energetic electrons diffuse rapidly to the wall, so the electron distribution function displays a depleted high-energy tail and in many cases, the EEDF, 
$f(\unicode[STIX]{x1D700})$ has a form with the medium-energy range being enhanced while the low- and high-energy ranges are depleted (Amemiya Reference Amemiya2012). In this case, the EEDF may be represented by a Druyvesteyn-like distribution with a parameter 
$x$ controlling the shape of the EEDF. As 
$x$ varies from 
$1$ to 
$2$, the electron energy distribution function transforms from a Maxwellian to a Druyvesteyn-like distribution.
Due to the importance of the electron distribution function in gas discharges, many authors have included its study in their research on different gases, for example, in nitrogen (Nighan Reference Nighan1970; Loureiro & Ferreira Reference Loureiro and Ferreira1986, Reference Loureiro and Ferreira1989), hydrogen (Hassouni, Grotjohn & Gicquel Reference Hassouni, Grotjohn and Gicquel1999; Diomede, Longo & Capitelli Reference Diomede, Longo and Capitelli2005; Longo & Diomede Reference Longo and Diomede2009), oxygen (Gudmundsson et al. Reference Gudmundsson, Marakhtanov, Patel, Gopinath and Leiberman2000) and argon/oxygen mixtures (Gudmundsson, Kimura & Leiberman Reference Gudmundsson, Kimura and Leiberman1999). Behringer and Fantz used the Druyvesteyn distribution family to fit measurements of the electron distribution function in gas discharges for a nitrogen discharge (Behringer & Fantz Reference Behringer and Fantz1994). A sum of two Druyvesteyn distributions has been used to match both high- and low-energy electrons with different parameters at different conditions. The criterion for the formation of a stable ion sheath has been considered for an arbitrary electron energy distribution function 
$f(\unicode[STIX]{x1D700})$ to estimate the ion flux for the cases of Maxwellian, Druyvesteyn and other energy distributions (Amemiya Reference Amemiya1997, Reference Amemiya1998). A global (volume-averaged) model was applied to investigate the effect of the electron energy distribution function on the plasma parameters of a high-density, low-pressure, argon discharge (Gudmundsson Reference Gudmundsson2001). In Zhao et al. (Reference Zhao, Liao, Lin, Chang and Fu2011) a modified EEDF was introduced into the fluid model and its effects on the high power microwave breakdown in air and argon have been investigated. The Druyvesteyn EEDF was used for a global model in the study of argon–methane–hydrogen plasmas, applied for growth of vertically aligned carbon nanostructures (Denysenko et al. Reference Denysenko, Xu, Long, Rutkevych, Azarenkov and Ostrikov2004a). Analytical expressions describing the variation of the EEDF, including Maxwellian and Druyvesteyn distributions in the afterglow of a plasma, were obtained (Denysenko, Kersten & Azarenkov Reference Denysenko, Kersten and Azarenkov2016). A complex low-pressure argon discharge containing dust grains was studied and it was shown that, for sufficiently high grain density and/or size, in a low-pressure argon glow discharge, the Druyvesteyn-type electron distribution in pristine plasmas can become nearly Maxwellian (Denysenko et al. Reference Denysenko, Yu, Ostrikov, Azarenkov and Stenflo2004b).
However, related to the problem of the sheath formation in non-stationary plasmas, some authors attempted to include the effects of a non-Maxwellian distribution function by two temperature Maxwellian distributions (Foroutan & Akhoundi Reference Foroutan and Akhoundi2012b; Ou et al. Reference Ou, Xiang, Gan and Yang2013; Driouch, Chatei & Bojaddaini Reference Driouch, Chatei and Bojaddaini2015) and a 
$q$ non-extensive electron distribution (Safa, Ghomi & Niknam Reference Safa, Ghomi and Niknam2014; Hatami Reference Hatami2015; Driouch & Chatei Reference Driouch and Chatei2017; Borgohain & Saharia Reference Borgohain and Saharia2018a,Reference Borgohain and Sahariab; Ghani, Driouch & Chatei Reference Ghani, Driouch and Chatei2019), where 
$q$ is the parameter which represents the degree of nonextensivity. Since a Druyvesteyn-type function for the electrons is the actual distribution in a plasma discharge, in the present work we take into account the Druyvesteyn-type distribution for electrons and study the plasma sheath structure.
This paper is organized as follows. In § 2 we introduce the EEDF and evaluate the number density and current of the electrons for a given 
$f(\unicode[STIX]{x1D700})$. Then, § 3 presents the sheath model and sets up the basic equations to be solved numerically. The results of the numerical simulations are presented and discussed in § 4. Finally, the main findings are summarized and enumerated in § 5.
2 Electron energy distribution function 
$f(\unicode[STIX]{x1D700})$
In comparison with a Maxwellian distribution, the Druyvesteyn distribution is characterized by depletion in the high-energy region and a shift of the maximum towards higher energies. The general form for the electron energy distribution function can be expressed as (Amemiya Reference Amemiya1998; Gudmundsson Reference Gudmundsson2001)
$$\begin{eqnarray}f(\unicode[STIX]{x1D700})=\frac{c_{1}\sqrt{\unicode[STIX]{x1D700}}}{\unicode[STIX]{x1D700}_{m}^{3/2}}\exp \left(-c_{2}\left(\frac{\unicode[STIX]{x1D700}}{\unicode[STIX]{x1D700}_{m}}\right)^{x}\right).\end{eqnarray}$$ Here, 
$x=1$ and 
$x=2$ correspond to the Maxwellian and Druyvesteyn electron energy distributions, respectively.
In (2.1) 
$\unicode[STIX]{x1D700}$ is the electron energy and 
$\unicode[STIX]{x1D700}_{m}=\frac{3}{2}T_{e}$ is the average energy of the electrons, where 
$T_{e}$ is the electron temperature. The parameters 
$c_{1}$ and 
$c_{2}$ are defined as
$$\begin{eqnarray}\displaystyle & \displaystyle c_{1}=x\frac{[\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D709}_{2})]^{3/2}}{[\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D709}_{1})]^{5/2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle c_{2}=\left[\frac{\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D709}_{2})}{\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D709}_{1})}\right]^{x}. & \displaystyle\end{eqnarray}$$Here
$$\begin{eqnarray}\unicode[STIX]{x1D6E4}(a)=\int _{0}^{\infty }t^{a-1}\text{e}^{-t}\,\text{d}t\end{eqnarray}$$ is the gamma function and 
$\unicode[STIX]{x1D709}_{1}=3/2x$ and 
$\unicode[STIX]{x1D709}_{2}=5/2x$.
Profiles of the EEDF versus 
$\unicode[STIX]{x1D700}$ with 
$T_{e}=1$ eV for different values of 
$x$: 
$x=1.0$ (solid), 
$x=1.5$ (dashed) and 
$x=2.0$ (dotted).
It is noted that (2.1) is not the most general EEDF expression. The actual EEDF can be calculated using freely available codes like BOLSIG+, based on the two-term approximation in spherical harmonics of the angular dependence of the velocity distribution (Rockwood Reference Rockwood1973; Hagelaar Reference Hagelaar2016).
The factor 
$x$ controls the behaviour of the Druyvesteyn-type distribution. Figure 1 shows the electron energy distribution function as a function of the electron energy 
$\unicode[STIX]{x1D700}$ for 
$T_{e}=1~\text{eV}$ for three different values of 
$x$. It is seen that, as 
$x$ increases, 
$f(\unicode[STIX]{x1D700})$ decreases in the low- and high-energy regions while at intermediate energies, its value is enhanced and the maximum of the function moves towards the higher-energy regime. At high energies the reduction in electron number is small, therefore, one finds that, with increasing 
$x$, the number of electrons with a higher mean energy will be increased.
According to the Bohm criterion, for the electrostatic sheath formation, ions must have a supersonic flow speed as they leave the quasi-neutral plasma and enter the sheath. The generalized Bohm criterion can be evaluated from the electron energy distribution function 
$f(\unicode[STIX]{x1D700})$ as (Lieberman & Lichtenberg Reference Lieberman and Lichtenberg1994; Amemiya Reference Amemiya1997)
$$\begin{eqnarray}\frac{1}{m_{i}}\langle v_{i}^{-2}\rangle =\int _{0}^{\infty }\frac{f(\unicode[STIX]{x1D700})}{2\unicode[STIX]{x1D700}}\,\text{d}\unicode[STIX]{x1D700},\end{eqnarray}$$ where 
$m_{i}$ and 
$v_{i}$ are the ion mass and ion velocity, respectively. For a Druyvesteyn-type plasma the Bohm speed is a function of 
$x$ and is given by Lieberman & Lichtenberg (Reference Lieberman and Lichtenberg1994) and Takahashi et al. (Reference Takahashi, Charles, Boswell and Fujiwara2011) as
$$\begin{eqnarray}v_{i}(x)=\sqrt{\frac{T_{e}}{m_{i}}}\frac{\sqrt{3}\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D709}_{1})}{[\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D709}_{2})\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D709}_{3})]^{1/2}},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D709}_{3}=1/2x$ and for a Maxwellian plasma 
$v_{i}^{x=1}=\sqrt{T_{e}/m_{i}}$ and for the Druyvesteyn this yields 
$v_{i}^{x=2}=1.17v_{i}^{x=1}$.
2.1 Electron number density
Knowing the electron energy distribution function allows us to evaluate the electron number density and the incident electron current on the dust. For a non-zero potential 
$\unicode[STIX]{x1D719}$, the electron number density 
$n_{e}$ is given by (Amemiya Reference Amemiya1997)
$$\begin{eqnarray}n_{e}=n_{e0}\int _{e\unicode[STIX]{x1D719}}^{\infty }\sqrt{1-\frac{e\unicode[STIX]{x1D719}}{\unicode[STIX]{x1D700}}}f(\unicode[STIX]{x1D700})\,\text{d}\unicode[STIX]{x1D700}.\end{eqnarray}$$ Substitution of (2.1) into (2.7) and using 
$\unicode[STIX]{x1D700}/\unicode[STIX]{x1D700}_{m}=w$ we have
$$\begin{eqnarray}n_{e}=n_{e0}c_{1}\int _{e\unicode[STIX]{x1D719}/\unicode[STIX]{x1D700}_{m}}^{\infty }\sqrt{w-\frac{e\unicode[STIX]{x1D719}}{\unicode[STIX]{x1D700}_{m}}}\exp (-c_{2}w^{x})\,\text{d}w.\end{eqnarray}$$ For a Maxwellian distribution, 
$x=1$ and 
$c_{1}=2(1.5)^{1.5}/\sqrt{\unicode[STIX]{x03C0}}$, 
$c_{2}=3/2$, the electron energy distribution function is given as:
$$\begin{eqnarray}f(\unicode[STIX]{x1D700})=\frac{2}{\sqrt{\unicode[STIX]{x03C0}}T_{e}^{3/2}}\sqrt{\unicode[STIX]{x1D700}}\exp (-\unicode[STIX]{x1D700}/T_{e})\end{eqnarray}$$and we find the well-known relation for the electron number density for a Maxwellian distribution
$$\begin{eqnarray}n_{e}=n_{e0}\exp (e\unicode[STIX]{x1D719}/T_{e}).\end{eqnarray}$$ For other values of 
$x$ the integration of (2.8) will be evaluated numerically.
2.2 The dust charging currents
To evaluate the electron current on the dust, consider a finite size dust immersed in a plasma. The radius of a dust grain is 
$r_{d}$ and the potential of the dust relative to the plasma potential is 
$\unicode[STIX]{x1D719}_{d}$ and it is proportional to the dust charge 
$q_{d}$ via 
$\unicode[STIX]{x1D719}_{d}=q_{d}/r_{d}$. The electrons are collected by the dust and the cross-section of charging collisions is given by (Sutherland Reference Sutherland1893)
$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{e}=\unicode[STIX]{x03C0}r_{d}^{2}\left(1+\frac{e\unicode[STIX]{x1D719}_{d}}{\unicode[STIX]{x1D700}}\right).\end{eqnarray}$$ Using the orbital motion limited (OML) theory, the relation between the electron energy distribution function 
$f(\unicode[STIX]{x1D700})$ and the current of electrons 
$I_{e}$ on a dust particle is given by Amemiya (Reference Amemiya1997), Shukla & Mamun (Reference Shukla and Mamun2002) and Denysenko et al. (Reference Denysenko, Kersten and Azarenkov2016) as
$$\begin{eqnarray}I_{e}=-en_{e}\sqrt{\frac{2}{m_{e}}}\int _{-e\unicode[STIX]{x1D719}_{d}}^{\infty }\unicode[STIX]{x1D70E}_{e}\sqrt{\unicode[STIX]{x1D700}}f(\unicode[STIX]{x1D700})\,\text{d}\unicode[STIX]{x1D700}\end{eqnarray}$$or
$$\begin{eqnarray}I_{e}=-en_{e}\unicode[STIX]{x03C0}r_{d}^{2}\sqrt{\frac{2}{m_{e}}}\int _{-e\unicode[STIX]{x1D719}_{d}}^{\infty }\sqrt{\unicode[STIX]{x1D700}}\left(1+\frac{e\unicode[STIX]{x1D719}_{d}}{\unicode[STIX]{x1D700}}\right)f(\unicode[STIX]{x1D700})\,\text{d}\unicode[STIX]{x1D700}.\end{eqnarray}$$ Using (2.1) and considering 
$w=\unicode[STIX]{x1D700}/T_{e}$, 
$\unicode[STIX]{x1D6F7}_{d}=e\unicode[STIX]{x1D719}_{d}/T_{e}$, 
$\acute{c} _{1}=c_{1}/(1.5)^{1.5}$ and 
$\acute{c} _{2}=c_{2}/(1.5)^{x}$ we have the electron current on the dust for 
$e\unicode[STIX]{x1D719}_{d}<0$
$$\begin{eqnarray}I_{e}=-\unicode[STIX]{x03C0}r_{d}^{2}n_{e}e\sqrt{\frac{2}{m_{e}}}\acute{c} _{1}T_{e}^{1/2}\left[\int _{-\unicode[STIX]{x1D6F7}_{d}}^{\infty }w\exp (-\acute{c} _{2}w^{x})\,\text{d}w+\unicode[STIX]{x1D6F7}_{d}\int _{-\unicode[STIX]{x1D6F7}_{d}}^{\infty }\exp (-\acute{c} _{2}w^{x})\,\text{d}w\right]\end{eqnarray}$$and finally,
$$\begin{eqnarray}I_{e}=-\unicode[STIX]{x03C0}r_{d}^{2}n_{e}e\sqrt{\frac{2}{m_{e}}}\acute{c} _{1}T_{e}^{1/2}\left[\frac{1}{x\acute{c} _{2}^{2/x}}\unicode[STIX]{x1D6E4}(2/x,\acute{c} _{2}\unicode[STIX]{x1D6F7}_{d}^{x})+\frac{\unicode[STIX]{x1D6F7}_{d}}{x\acute{c} _{2}^{1/x}}\unicode[STIX]{x1D6E4}(1/x,\acute{c} _{2}\unicode[STIX]{x1D6F7}_{d}^{x})\right].\end{eqnarray}$$Here,
$$\begin{eqnarray}\unicode[STIX]{x1D6E4}(a,z)=\int _{a}^{\infty }t^{z-1}\text{e}^{-t}\,\text{d}t\end{eqnarray}$$ is the incomplete gamma function, and for 
$\unicode[STIX]{x1D6F7}_{d}>0$
$$\begin{eqnarray}I_{e}=-en_{e}\unicode[STIX]{x03C0}r_{d}^{2}\sqrt{\frac{2}{m_{e}}}\acute{c} _{1}T_{e}^{1/2}\left[\frac{1}{x\acute{c} _{2}^{2/x}}\unicode[STIX]{x1D6E4}(2/x)+\frac{\unicode[STIX]{x1D6F7}_{d}}{x\acute{c} _{2}^{1/x}}\unicode[STIX]{x1D6E4}(1/x)\right].\end{eqnarray}$$ For a Maxwellian EEDF, 
$\acute{c} _{1}=2/\sqrt{\unicode[STIX]{x03C0}}$ and 
$\acute{c} _{2}=1$ respectively, and so the well-known electron currents in a Maxwellian plasma can be obtained (Shukla & Mamun Reference Shukla and Mamun2002)
$$\begin{eqnarray}I_{e}=-en_{e}\unicode[STIX]{x03C0}r_{d}^{2}\sqrt{\frac{8T_{e}}{\unicode[STIX]{x03C0}m_{e}}}\exp \left(\frac{eq_{d}}{r_{d}T_{e}}\right)\quad q_{d}<0\end{eqnarray}$$and
$$\begin{eqnarray}I_{e}=-en_{e}\unicode[STIX]{x03C0}r_{d}^{2}\sqrt{\frac{8T_{e}}{\unicode[STIX]{x03C0}m_{e}}}\left(1+\frac{eq_{d}}{r_{d}T_{e}}\right)\quad q_{d}>0.\end{eqnarray}$$ In the Druyvesteyn limit 
$\acute{c} _{1}=1/\sqrt{\unicode[STIX]{x03C0}}$ and 
$\acute{c} _{2}=0.243$ respectively, and we have (Khrapak Reference Khrapak2010)
$$\begin{eqnarray}\displaystyle I_{e} & = & \displaystyle -en_{e}\unicode[STIX]{x03C0}r_{d}^{2}\sqrt{\frac{8T_{e}}{\unicode[STIX]{x03C0}m_{e}}}\left(\frac{3}{2\sqrt{2}}\right)^{1/2}\nonumber\\ \displaystyle & & \displaystyle \quad \times \left[\exp \left(-0.243\left(\frac{eq_{d}}{r_{d}T_{e}}\right)^{2}\right)\right.\nonumber\\ \displaystyle & & \displaystyle \quad -\left.\sqrt{0.243\unicode[STIX]{x03C0}}\left(\frac{eq_{d}}{r_{d}T_{e}}\right)\text{erfc}\left(\sqrt{0.243}\left(\frac{eq_{d}}{r_{d}T_{e}}\right)\right)\right]\quad q_{d}<0\end{eqnarray}$$and
$$\begin{eqnarray}I_{e}=-en_{e}\unicode[STIX]{x03C0}r_{d}^{2}\sqrt{\frac{8T_{e}}{\unicode[STIX]{x03C0}m_{e}}}\left(\frac{3}{2\sqrt{2}}\right)^{1/2}\left(1+\sqrt{0.243\unicode[STIX]{x03C0}}\frac{eq_{d}}{r_{d}T_{e}}\right)\quad q_{d}>0.\end{eqnarray}$$3 Plasma sheath model
In this section we first present the multi-fluid sheath model which contains the electrons, ions and nanosized charged dust particles. Our study is based on the hydrodynamic model, but it should be mentioned that the sheath region can also be studied using a kinetic approach via particle in cell simulations with Monte Carlo collisions (Birdsall Reference Birdsall1991). It is assumed that the plasma sheath boundary is located at 
$z=0$. The influences of the electric, gravitational, ion drag and natural forces are included in the kinetics of the dust particles. The number density of electrons is given by (2.8) and for other particles the governing equations are the continuity
$$\begin{eqnarray}\frac{\text{d}}{\text{d}z}(n_{d}v_{d})=0,\end{eqnarray}$$and
$$\begin{eqnarray}\frac{\text{d}}{\text{d}z}(n_{i}v_{i})=\unicode[STIX]{x1D708}_{iz}n_{e}-\unicode[STIX]{x1D708}_{id}n_{d},\end{eqnarray}$$and momentum equations
$$\begin{eqnarray}\displaystyle & \displaystyle m_{d}v_{d}\frac{\text{d}v_{d}}{\text{d}z}=-q_{d}\frac{\text{d}\unicode[STIX]{x1D719}}{\text{d}z}+m_{d}g+F_{id}+F_{nd}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle m_{i}v_{i}\frac{\text{d}v_{i}}{\text{d}z}=-\text{e}\frac{\text{d}\unicode[STIX]{x1D719}}{\text{d}z}-\frac{T_{i}}{n_{i}}\frac{\text{d}n_{i}}{\text{d}z}-m_{i}v_{i}\unicode[STIX]{x1D708}_{in}, & \displaystyle\end{eqnarray}$$ for the dust grains and ions. Here, 
$m_{j}$ is the mass, 
$n_{j}$ is the number density, 
$v_{j}$ is the fluid velocity of the particle species 
$j$ and 
$\unicode[STIX]{x1D719}$, 
$T_{i}$, 
$\unicode[STIX]{x1D708}_{iz}$, 
$\unicode[STIX]{x1D708}_{in}$ and 
$\unicode[STIX]{x1D708}_{id}$ are the sheath electric potential, the ion temperature, the ionization frequency and the collision frequency of the ions with neutrals and and dust grains, respectively. The Poisson equation for the electrostatic potential is given by
$$\begin{eqnarray}\displaystyle \frac{\text{d}^{2}\unicode[STIX]{x1D719}}{\text{d}z^{2}}=4\unicode[STIX]{x03C0}(e(n_{e}-n_{i})-n_{d}q_{d}) & & \displaystyle\end{eqnarray}$$ and the quasi-neutrality condition applies at 
$z=0$
$$\begin{eqnarray}e(n_{i0}-n_{e0})+n_{d0}q_{d0}=0,\end{eqnarray}$$ where the index 
$0$ refers to the values of the plasma and dust parameters at the sheath edge, 
$z=0$. The ion drag force consists of two parts, the collection drag and Coulomb drag (Barnes et al. Reference Barnes, Keller, Forster, ONeill and Coultas1992; Shukla & Mamun Reference Shukla and Mamun2002)
$$\begin{eqnarray}\displaystyle & \displaystyle F_{id}^{\text{coll}}=\unicode[STIX]{x03C0}r_{d}^{2}n_{i}m_{i}v_{si}\tilde{v}_{i}\left(1-\frac{2eq_{d}}{r_{d}m_{i}v_{si}^{2}}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle F_{id}^{\text{coul}}=2\unicode[STIX]{x03C0}b_{0}^{2}n_{i}m_{i}v_{si}\tilde{v}_{i}\ln \left(\frac{b_{0}^{2}+\unicode[STIX]{x1D706}_{De}^{2}}{b_{0}^{2}+b_{c}^{2}}\right), & \displaystyle\end{eqnarray}$$ where 
$v_{si}=(\tilde{v}_{i}^{2}+8T_{i}/\unicode[STIX]{x03C0}m_{i})^{1/2}$ is the mean ion speed, 
$\tilde{v}_{i}=v_{i}-v_{d}$ is the directed ion speed relative to a dust particle with speed 
$v_{d}$, 
$b_{0}=eq_{d}/m_{i}v_{i}^{2}$ is the impact radius corresponding to a 
$90^{\circ }$ deflection, 
$\unicode[STIX]{x1D706}_{De}=\sqrt{T_{e}/4\unicode[STIX]{x03C0}n_{e0}e^{2}}$ is the electron Debye length and 
$b_{c}=r_{d}(1-2b_{0}/r_{d})^{1/2}$ is the direct collision impact parameter. The neutral drag force 
$F_{nd}$ is given by the Epstein expression
$$\begin{eqnarray}F_{nd}=-{\textstyle \frac{3}{8}}\sqrt{2\unicode[STIX]{x03C0}}r_{d}^{2}m_{n}n_{n}v_{thn}v_{d},\end{eqnarray}$$ where 
$v_{d}\ll v_{thn}$ and 
$v_{thn}=\sqrt{T_{n}/m_{n}}$ is the thermal speed of the neutrals.
To obtain the dust charge 
$q_{d}$, it is assumed that the charging process is a local phenomenon, i.e. the dust charging time is much smaller than the characteristic time for the dust motion. Therefore, the dust charge is determined from the condition that, at equilibrium, the net current flowing onto the grains is equal to zero
$$\begin{eqnarray}I_{e}+I_{i}=0.\end{eqnarray}$$ Using (2.15) and (2.17) and replacing 
$\unicode[STIX]{x1D6F7}_{d}$ with 
$eq_{d}/r_{d}$, the electron and ion currents 
$I_{e}$ and 
$I_{i}$ are given by the OML theory as
$$\begin{eqnarray}\displaystyle I_{e} & = & \displaystyle -\unicode[STIX]{x03C0}r_{d}^{2}n_{e}e\left(\frac{2T_{e}}{m_{e}}\right)^{1/2}\acute{c} _{1}\nonumber\\ \displaystyle & & \displaystyle \times \left[\frac{1}{x\acute{c} _{2}^{2/x}}\unicode[STIX]{x1D6E4}(2/x,\acute{c} _{2}(eq_{d}/r_{d})^{x})+\frac{(eq_{d}/r_{d})}{x\acute{c} _{2}^{1/x}}\unicode[STIX]{x1D6E4}(1/x,\acute{c} _{2}(eq_{d}/r_{d})^{x})\right]\end{eqnarray}$$
$$\begin{eqnarray}I_{i}=\unicode[STIX]{x03C0}r_{d}^{2}ev_{i}n_{i}\left(1-\frac{2eq_{d}}{r_{d}m_{i}v_{i}^{2}}\right)\end{eqnarray}$$ for 
$q_{d}<0$ and
$$\begin{eqnarray}\displaystyle & \displaystyle I_{e}=-\unicode[STIX]{x03C0}r_{d}^{2}n_{e}e\left(\frac{2T_{e}}{m_{e}}\right)^{1/2}\acute{c} _{1}\left[\frac{1}{x\acute{c} _{2}^{2/x}}\unicode[STIX]{x1D6E4}(2/x)+\frac{(eq_{d}/r_{d})}{x\acute{c} _{2}^{1/x}}\unicode[STIX]{x1D6E4}(1/x)\right] & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle I_{i}=\unicode[STIX]{x03C0}r_{d}^{2}ev_{i}n_{i}\exp \left(\frac{2eq_{d}}{r_{d}m_{i}v_{i}^{2}}\right), & \displaystyle\end{eqnarray}$$ for 
$q_{d}<0$. By using (3.10), the normalized dust charge 
$Q=eq_{d}/r_{d}T_{e}$ is obtained for 
$Q<0$ from
$$\begin{eqnarray}\displaystyle & & \displaystyle n_{e}\left(\frac{T_{e}}{m_{e}}\right)^{1/2}\acute{c} _{1}\left[\frac{1}{x\acute{c} _{2}^{2/x}}\unicode[STIX]{x1D6E4}(2/x,\acute{c} _{2}Q^{x})+\frac{Q}{x\acute{c} _{2}^{1/x}}\unicode[STIX]{x1D6E4}(1/x,\acute{c} _{2}Q^{x})\right]\nonumber\\ \displaystyle & & \displaystyle \quad =n_{i}\left(\frac{4T_{i}}{\unicode[STIX]{x03C0}m_{i}}\right)^{1/2}\tilde{v}_{si}\left[1-\frac{T_{e}}{T_{i}}\frac{\unicode[STIX]{x03C0}Q}{4\tilde{v}_{si}^{2}}\right]\end{eqnarray}$$ and for 
$Q>0$
$$\begin{eqnarray}\left(\frac{T_{e}}{m_{e}}\right)^{1/2}n_{e}\acute{c} _{1}\left[\frac{1}{x\acute{c} _{2}^{2/x}}\unicode[STIX]{x1D6E4}(2/x)+\frac{Q}{x\acute{c} _{2}^{1/x}}\unicode[STIX]{x1D6E4}(1/x)\right]=n_{i}\left(\frac{4T_{i}}{\unicode[STIX]{x03C0}m_{i}}\right)^{1/2}\tilde{v}_{si}\exp \left[-\frac{T_{e}}{T_{i}}\frac{\unicode[STIX]{x03C0}Q}{4\tilde{v}_{si}^{2}}\right],\end{eqnarray}$$ where 
$\tilde{v}_{si}=(1+v_{i}^{2}/v_{thi}^{2})^{1/2}$ and 
$v_{thi}=\sqrt{8T_{i}/\unicode[STIX]{x03C0}m_{i}}$ is the ion mean thermal speed. For numerical simulation, we normalize equations (3.1)–(3.6) using the following set of dimensionless variables:
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D6F7}=e\unicode[STIX]{x1D719}/T_{e},\quad Q=eq_{d}/r_{d}T_{e},\quad \unicode[STIX]{x1D709}=z/\unicode[STIX]{x1D706}_{De},\\ N_{e,i,d}=n_{e,i,d}/n_{e0},\quad u_{i}=v_{i}/c_{i},\quad u_{d}=v_{d}/c_{d},\\ \bar{\unicode[STIX]{x1D708}}_{iz,in,id}=\unicode[STIX]{x1D708}_{iz,in,id}/\unicode[STIX]{x1D714}_{pi},\end{array}\right\}\end{eqnarray}$$ where 
$c_{i}=\sqrt{T_{e}/m_{i}}$ and 
$\unicode[STIX]{x1D714}_{pi}=\sqrt{4\unicode[STIX]{x03C0}n_{i0}e^{2}/m_{i}}$ are the ion acoustic speed and ion plasma frequency, and 
$c_{d}=\sqrt{ZT_{e}/m_{d}}$ is the dust-acoustic speed. The parameter 
$Z=r_{d}T_{e}/e^{2}$ is the dust charge number at 
$Q=1$, i.e. when the floating potential on the grain surface is 
$T_{e}/e$. Using the above normalization, the following set of coupled differential equations can be obtained from (2.8) and (3.1)–(3.6)
$$\begin{eqnarray}\displaystyle & \displaystyle N_{e}=c_{1}\int _{2/3\unicode[STIX]{x1D6F7}}^{\infty }\sqrt{w-\frac{2}{3}~\unicode[STIX]{x1D6F7}}\exp (-c_{2}w^{x})\,\text{d}w, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}N_{d}}{\text{d}\unicode[STIX]{x1D709}}=\frac{N_{d}}{u_{d}^{2}}\left[Q\frac{\text{d}\unicode[STIX]{x1D6F7}}{\text{d}\unicode[STIX]{x1D709}}-\unicode[STIX]{x1D707}(m_{d}\text{g}+F_{id}+F_{nd})\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}u_{d}}{\text{d}\unicode[STIX]{x1D709}}=-\frac{1}{u_{d}}\left[Q\frac{\text{d}\unicode[STIX]{x1D6F7}}{\text{d}\unicode[STIX]{x1D709}}-\unicode[STIX]{x1D707}(m_{d}\text{g}+F_{id}+F_{nd})\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}N_{i}}{\text{d}\unicode[STIX]{x1D709}}=\frac{T_{e}}{T_{e}u_{i}^{2}-T_{i}}\left[N_{i}\frac{\text{d}\unicode[STIX]{x1D6F7}}{\text{d}\unicode[STIX]{x1D709}}+u_{i}\unicode[STIX]{x1D6FF}(N_{e}\bar{\unicode[STIX]{x1D708}}_{iz}+N_{i}\bar{\unicode[STIX]{x1D708}}_{in}-N_{d}\bar{\unicode[STIX]{x1D708}}_{id})\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}u_{i}}{\text{d}\unicode[STIX]{x1D709}}=\frac{T_{e}u_{i}}{T_{i}-T_{e}u_{i}^{2}}\left[\frac{\text{d}\unicode[STIX]{x1D6F7}}{\text{d}\unicode[STIX]{x1D709}}+u_{i}\unicode[STIX]{x1D6FF}\bar{\unicode[STIX]{x1D708}}_{in}+\frac{\unicode[STIX]{x1D6FF}T_{i}}{u_{i}N_{i}T_{e}}(N_{e}\bar{\unicode[STIX]{x1D708}}_{iz}-N_{d}\bar{\unicode[STIX]{x1D708}}_{id})\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}^{2}\unicode[STIX]{x1D6F7}}{\text{d}\unicode[STIX]{x1D709}}=N_{e}-N_{i}-ZN_{d}Q, & \displaystyle\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FF}=\unicode[STIX]{x1D706}_{De}\unicode[STIX]{x1D714}_{pi}/c_{i}=\sqrt{N_{i0}}$, 
$\unicode[STIX]{x1D707}=\unicode[STIX]{x1D706}_{De}/m_{d}c_{d}^{2}$. Equations (3.18)–(3.23) with boundary conditions are solved using a fourth-order Runge–Kutta method. Simultaneously, equations (3.15) and (3.16) are solved to obtain the dust charge at every position self-consistently.
4 Numerical results
This section presents the results of numerical simulations of the multi-fluid equations describing an electrostatic plasma sheath including the electrons, ions and charged dust nanoparticles. The used parameters are listed in table 1 and are chosen as default parameters of the simulations. In addition, at 
$\unicode[STIX]{x1D709}=0$ we have 
$\unicode[STIX]{x1D6F7}=0$, 
$\text{d}\unicode[STIX]{x1D6F7}/\text{d}\unicode[STIX]{x1D709}=\bar{\unicode[STIX]{x1D708}}_{in}$ and the wall potential is 
$\unicode[STIX]{x1D6F7}_{w}=-10$.
The parameters used for simulations.
As mentioned above, the change in 
$x$ affects the shape of the electron distribution function and consequently sheath parameters will be influenced by 
$x$. Here, we investigate the effects of 
$x$ on the plasma sheath dynamics. Figures 2(a)–2(d) display the profiles of the sheath potential, number densities of the electrons and ions and the ion speed functions of 
$\unicode[STIX]{x1D709}$ for three different values of, 
$x=1$, 
$1.5$ and 
$x=2$. It is found that the sheath thickness increases with increasing 
$x$. This is due to the fact that, when 
$x$ increases, as discussed above, the mean electron energy increases and, as a result, their flux increases towards the wall leading to an increased number density of electrons in the sheath (figure 2b) and then the sheath is broadened. Figures 2(c) and 2(d) show that the ion density increases with 
$x$ but the ion speed decreases.
Figures 3(a)–3(c) show the charge, speed and number density of the dust nanoparticles for three different values of 
$x$ as in figure 2. It is seen that the absolute dust charge is reduced but the dust number density is increased with an increase in 
$x$. The enhancement in the dust number density can be explained on the basis of the quasi-neutrality condition at the sheath edge. With a reduction of the absolute dust charge with 
$x$, the quasi-neutrality condition at the sheath edge, for given electron and ion number densities, leads to a higher dust density.
Profiles of (a) the sheath potential, (b) the electron number density, (c) ion number density and (d) ion speed, as functions of 
$\unicode[STIX]{x1D709}$ for three different values of the parameter 
$x$ as in figure 1.
Spatial distributions of (a) the dust charge, (b) dust speed and (c) dust number density for the same values of 
$x$ as in figure 1.
From figure 3(b) we find that, initially, the dust speed decreases and then increases with increasing 
$x$. To explore this behaviour we look again at figure 3(a). Near the sheath edge the absolute dust charge first decreases but at larger 
$\unicode[STIX]{x1D709}$ it slightly increases. With decreasing dust charge both the decelerating electric force and the accelerating ion drag force are reduced. The dependences of these two forces on 
$q_{d}$ are different. The electric force is directly proportional to 
$q_{d}$, whereas the Coulomb component of the ion drag is proportional to 
$q_{d}^{2}$. With an increase in 
$x$ and thus decrease in 
$q_{d}$, the electric force overcomes the ion drag and decelerates the dust particles and the speed of the dust particles decreases at small 
$\unicode[STIX]{x1D709}$ but for larger 
$\unicode[STIX]{x1D709}$ the inverse is true and, at larger sheath widths, with increasing 
$q_{d}$, the decelerating electric force is weakened and the accelerating ion drag will be dominant, which results in an enhancement of the dust speed.
Now, we investigate the dependence of the sheath thickness 
$\unicode[STIX]{x1D709}_{s}$, on the parameter of the electron energy distribution function 
$x$. Figure 4 represents 
$\unicode[STIX]{x1D709}_{s}$ as a function of 
$x$. It is seen that the sheath thickness increases with 
$x$, which is in agreement with the results of figures 2 and 3.
The sheath width 
$\unicode[STIX]{x1D709}_{s}$ as a function of the 
$x$.
Next, we study the effects of the dust parameters i.e. the dust radius and dust number density, on the spatial distribution of the plasma and dust particles in a plasma with a Druyvesteyn electron distribution (
$x=2$). Figure 5 displays the profiles of the electron and ion number densities and the charge and number density of the nanoparticles for three different dust radii. First, we see that the sheath thickness and the number densities of the electrons, ions and dust particles decrease monotonically with the dust radius. As the dust radius increases, more electrons and ions are collected by the dust grains. As a result, the number densities of electron and ions decrease in the sheath, as seen from figures 5(a) and 5(b). Also, due to the high mobility of the electrons the number of electrons collected by the dust is higher than that of ions, and consequently, the ion to electron ratio increases in the sheath and the sheath thickness is reduced.
Spatial profiles of the electron number density (a), ion number density (b), dust charge (c) and dust number density (d) with 
$x=2.0$ for three different values of 
$r_{d}=5.0\times 10^{-6}~\text{cm}$ (solid), 
$r_{d}=1.0\times 10^{-5}~\text{cm}$ (dashed) and 
$r_{d}=2.0\times 10^{-5}~\text{cm}$ (dotted).
To investigate the effects of the dust number density, figure 6 represents the number densities of the electrons, ions and dust nanoparticles, as well as the dust charge for three different values of the ion to electron number density ratio 
$\unicode[STIX]{x1D6FF}^{2}$ at the sheath edge for 
$x=2$. Because of the quasi-neutrality condition at this point, 
$\unicode[STIX]{x1D6FF}^{2}$ is proportional to the dust number density at the sheath edge. With increasing 
$\unicode[STIX]{x1D6FF}^{2}$ (or equivalently the dust number density), more electrons attach to the dust particles, so the ion to electron density ratio rises inside the sheath (figure 6b), and the sheath thickness is reduced. Figure 6(c) shows that the absolute dust charge decreases with the dust number density.
Spatial profiles of the electron number density (a), ion number density (b), dust charge (c) and dust number density (d) with 
$x=2.0$ for three different values of 
$\unicode[STIX]{x1D6FF}^{2}=3.5$ (solid), 
$\unicode[STIX]{x1D6FF}^{2}=5.5$ (dashed) and 
$\unicode[STIX]{x1D6FF}^{2}=6.5$ (dotted).
In the above simulations we have assumed that the ions enter the sheath with a constant velocity. But, as mentioned in § 2, the Bohm velocity is a function of 
$x$. Now we explore the effects of the change in the ion velocity at the sheath edge with 
$x$. Figure 7 represents the number densities of the electrons and ions, ion velocity and dust charge, velocity and dust number density for three different values of 
$x$. The values of 
$u_{i0}$ according to equation (2.5) are 1.0, 1.11 and 1.17 for 
$x=1.0$, 1.5 and 2.0, respectively. The results are noticeably deferent from the pervious results. When we include the effects of a change in 
$u_{i0}$, it is seen that the sheath thickness decreases with increasing 
$u_{i0}$. This is reasonable because, when 
$u_{i0}$ increases at the sheath edge, the ion density inside the sheath increases and thus the sheath thickness decreases.
Spatial profiles of the electron number density (a), ion number density (b), ion velocity (c), dust charge (d), dust velocity (e) and dust number density (f) for three different values of 
$x$, with 
$x=1$ and 
$u_{i0}=1.0$ (solid), 
$x=1.5$ and 
$u_{i0}=1.11$ (dashed) and 
$x=2$ and 
$u_{i0}=1.17$ (dotted).
5 Discussion
Using numerical simulations of the multi-fluid equations, we investigated the structure of a dusty plasma sheath that consists of electrons, ions and charged nanoparticles. It was assumed that the electron distribution function is a Druyvesteyn-type distribution. The spatial profiles of the plasma and dust nanoparticles, the ion speed and the charge and speed of the nanoparticles were investigated for different values of 
$x$, which characterizes the behaviour of the electron distribution function. The main results can be summarized as follows:
For a Druyvesteyn distribution, 
$f(\unicode[STIX]{x1D700})$ has a form where the medium-energy range is enhanced while the low- and high-energy ranges are depleted. This means that the electrons in the Druyvesteyn electron distribution have a higher mean energy than the Maxwellian one. The behaviour of the plasma sheath depends on 
$x$ and the shape of the electron energy distribution function. As 
$x$ increases, it was shown that the sheath thickness increases due to the enhancement of the electron number density in the sheath. At the sheath edge, the absolute dust charge decreases and, as a result, according to the quasi-neutrality condition, the dust number density increases. Also, the effects of dust parameters on the sheath dynamics were investigated and it was found that, with an increase in the dust radius, the sheath width declines. This is due to the fact that, as 
$r_{d}$ increases, more plasma particles are collected by the nanoparticles and due to the higher mobility of the electrons, the ion to electron density ratio increases at the sheath and consequently the thickness of the sheath decreases. It was shown that the sheath width also decreases with the dust number density at the sheath edge.
Acknowledgement
This work supported by Payame Noor University.
