## 2 Virtual Force Algorithm Based on Plasma Yukawa Potential Model  
Generally, one sensor node should include both the sensor (for collecting and exchanging information) and the movement devices (such as the GPS equipped by the UAV itself). In our paper, we consider all these hardware devices as a whole sensor node. The dust plasma can generally form stable regular hexagonal lattices in 2D plane in the physical world, which can be used for the the optimal coverage model of the wireless sensor deployment in a virtual force algorithm. Here we mainly concentrated on a theoretical algorithm, which presents a possible 3D node deployment strategy in simulations. Due to the self-crystallizing feature of dust plasma under the effect of plasma Yukawa potential, we exploit this capability to deploy sensor nodes in an innovative way. Each sensor carried by an unmanned aerial vehicle is regarded as a dust particle in our VFA-DP-3D deployment algorithm. The node position information will then be computed by the virtual force (in the dusty plasma system) of the surrounding nodes. Based on the molecular dynamics simulation, when a sensor node is regarded as a dust particle, the corresponding motion equation can be written as follows:  
\[m\frac{d^{2}\vec{r}_{i}}{dt^{2}}=-\
abla\phi_{i}-m\gamma\frac{d\vec{r}_{i}}{dt} \tag{1}\]  
or  
\[m\vec{a}_{i}=-\
abla\phi_{i}-m\gamma\vec{v}_{i}, \tag{2}\]  
where \(m\) is mass of the particle, \(\vec{r}_{i}\) is the displacement vector of particle \(i\) and \(\gamma\) is the friction coefficient. It should be noted that the gradient \(\
abla\phi_{i}\) of the electrostatic potential \(\phi_{i}\) in 3D Cartesian Coordinates can be expressed as \(\
abla\phi_{i}=(\partial\phi_{i}/\partial x,\partial\phi_{i}/\partial y, \partial\phi_{i}/\partial z)\). Therefore, the formula of \(\phi_{i}\) in a 3D space is  
\[\phi_{i} = k\vec{r}_{i}^{2}(t)+\sum_{i\
eq j}^{N}\frac{Q^{2}}{4\pi\epsilon \phi_{i}r_{ij}}e^{-\frac{r_{ij}}{\lambda_{D}}} \tag{3}\] \[=k(x_{i}^{2}+y_{i}^{2}+z_{i}^{2})+\sum_{i\
eq j}^{N}\frac{Q^{2}}{ 4\pi\epsilon_{0}}\] \[\cdot[(x_{i}-x_{j})^{2}+(y_{i}-y_{j})^{2}+(z_{i}-z_{j})^{2}]^{-1/2}\] \[\cdot e^{-\frac{[(x_{i}-x_{j})^{2}+(y_{i}-y_{j})^{2}+(z_{i}-z_{j} )^{2}]^{1/2}}{\lambda_{D}}}.\]  
where \(k\) is the spring coefficient, \(r_{i}\) is the radial distance of particle \(i\), \(Q\) is the amount of charge, \(r_{ij}\) is the distance from particle \(i\) and particle \(j\), \(\lambda_{D}\) is the Debye length of the plasma and \(N\) is the total number of all particles in the communication range. On the right side of Equation (3), the first term \(k\vec{r}_{i}^{2}\) is the external harmonic potential, which provides the radial constraint so that the spring constant \(k\) can regulate the separation between particles; the second term \((Q^{2}/(4\pi\epsilon_{0}r_{ij}))\exp(-r_{ij}/\lambda_{D})\) is a Yukawa potential. When \(r_{ij}\) is greater than \(\lambda_{D}\), the Yukawa potential quickly decreases and the interactions between the particles become very weak. This means that \(\lambda_{D}\) is the shielding length which controls the range of particle interactions in the system. \(\Sigma\) is the sum of the Yukawa potential of particle \(i\) subjected to other particles in the computational scale \(R_{c}\). In actual WSN applications, \(\lambda_{D}\) is usually chosen to be far less than the communication distance \(R_{c}\) so that the interaction between two nodes is manly obtained from the effect of Yukawa potential. In order to analyze how all nodes in wireless sensor networks are deployed from the initial random distribution to the final steady topology, the molecular dynamic equations are used for simulations in this paper. A second-order leap-frog scheme was adopted and the time discretization operation is divided into the time step \(dt\). The force on each node is constant in \(dt\) and results in a uniform motion.