I Introduction
In the foreseeable future, there would be huge amounts of lowpower devices in wireless networks to perform information forwarding, data collecting, situation sensing [9]. These devices could be the nodes in Internet of Things (IoT), wireless sensor networks (WSNs) or DevicetoDevice (D2D) systems. These lowpower devices are usually powered by batteries and deployed in a broad area [14]. In order to prolong the network lifetime and maintain the sustainability of these nodes, farfield wireless power transfer (WPT) via radio frequency (RF) has been considered as a promising technology to energize massive batterypowered devices. The reason is that RFWPT can provide a flexible and longdistance charging service, while being compatible with the existing wireless information networks [47, 23]. At present, RFWPT has been successfully applied in various wireless networks [9, 14, 47, 23].
For addressing the demand of everincreasing data transmission rates, millimeter wave (mmWave) frequencies are leveraged to fulfill the multiGigabit information transmission requirement [31]. Due to the small wavelength and blockage sensitivity of mmWave signals, strong directional antenna and small cell structure are suggested to enhance the energy and spectrum efficiency of mmWave communications [3]. As WPT also suffers from the severe power propagation loss and needs to avoid interference to existing wireless information networks, mmWave also benefits WPT. It has been proven by [21, 40] that mmWave WPT outperforms WPT with lower frequencies. In [50], lens array based mmWave WPT was proposed to charge multiple energy receivers. In [11], simultaneous wireless information and power transfer (SWIPT) was applied in the hybrid precoding mmWave system. Besides, the effect of rain attenuation on mmWave WPT was investigated by [20]. All that literature shows mmWave WPT is feasible.
To evaluate the systemlevel performance of mmWave WPT, stochastic geometry has been utilized to capture the effects of propagation loss and blockage, which dominate the received signal power [5, 15]. The energy coverage probability and average harvested energy of mmWave WPT were studied by [21], where the location of energy transmitters follows the Poisson point process (PPP). Furthermore, SWIPT was also introduced into stochastic mmWave networks [40, 38]. Both works verified that mmWave could improve the performance of SWIPT compared to lower frequency solutions. In [22], the energy coverage probability in the presence of human blockage was derived. Discretizing harvested energy into a finite number of power levels, the total coverage probability integrating information and energy transmission was derived by [51]. A beamtraining based mmWave WPT scheme was proposed in [18], where the energy transmitter steers the energy beam along the direction in which it receives the strongest training signal. Considering the nonlinear behavior of energy harvesting, the coverage probability and average harvested energy were studied by [37]. In [41], the location of mmWave powered users was modeled as the Poisson cluster process and the energy and information coverage probabilities were derived.
Most of the aforementioned works, such as [21, 40, 38, 22, 51, 18, 37, 41], employ the flattop antenna pattern for mathematical tractability. Nevertheless, it may incur great inaccuracy for evaluating systemlevel performance [13, 46]. Therefore, some literature considers more realistic antenna patterns in the performance analysis of mmWave WPT. In [16], the microstrip patch antenna and endfire antenna models were used to analyze the coverage performance of SWIPT. [30] adopted the Fejér kernel model to represent the actual array antenna gain of energy transmitters. Both [16] and [30] assumed each energy transmitter is equipped with a directional antenna, while the antenna of energy receiver was assumed to be isotropic. In [39], both Tx and Rx were assumed to be equipped with directional antennas and the average harvested energy was derived.
Most of the existing works on mmWave communications/WPT have assumed that the beam direction of the associated TxRx pair is perfectly aligned. However, due to direction estimation error and hardware imperfection, the beam alignment error (BAE) is always inevitable, and affects the received signal strength greatly
[42]. Therefore, it is necessary to investigate the performance of the mmWave transmission systems with imperfect beam alignment (IBA). By now, much effort has been dedicated to evaluating the performance of the mmWave wireless information networks in the IBA scenario, such as [8, 4, 34, 46, 32, 10]. Taking into account of the BAEs in the associated TxRx pair and the Interfering TxsRx pairs, the information coverage probability with the flattop antenna gain model was derived by [8]. With the same BAE setup, while, [4] adopted the cosine antenna and isotropic antenna patterns to model the antenna gains of the transmitters and receivers, respectively. In [34], the approximated ergodic capacity loss of the mmWave ad hoc network with both flattop and Gaussian antenna gain models was derived. In [46], the authors employed the sinc and cosine antenna patterns to derive the information coverage probability, just considering the BAEs between the typical receiver and interfering transmitters. Alternatively, using the antenna gain model suggested by the 3rd Generation Partnership Project (3GPP), [32] investigated the impact of beam misalignment from interfering transmitters, and resorted to the curvefitting method to derive the probability density function of the antenna gain with BAE. Furthermore, [10] employed both the 2D and 3D directional antenna gains and derived the information coverage probability with BAE incurred by the strongest interfering link.To the best of our knowledge, the impact of BAE on mmWave WPT system has not been investigated and well understood to date. Additionally, in this paper, we consider the comprehensive IBA scenarios and the realistic 3GPP directional antenna model, which has not been well studied from the perspective of systemlevel performance evaluation in either mmWave communications or mmWave WPT. Specifically, there are two aspects of our system model to be highlighted. First, unlike [42, 8], we employ the 3GPP Gaussian gain model, which can approximate the realistic directional antenna pattern exactly, to explore the impact of BAE. The reason is that the flattop model lacks the capability of depicting the rolloff effect of the mainlobe and thus is not suitable for accurately evaluating the impact of BAE on mmWave systems [34, 46, 32, 4]. Second, different from [34, 46, 32, 10, 4], we consider a more realistic IBA scenario, where all Txs and Rxs are equipped with directional antennas. Moreover, not only the BAE between the associated TxRx pair but also the BAEs of the nonassociated Txs to the typical Rx are simultaneously taken into consideration. Besides, the nonlinear energy harvesting model in the typical receiver is also assumed. With these models, we derive the analytic expression of the systemlevel performance of the mmWave WPT in the IBA scenario. For clarity, we summarize our main contribution as follows:

We adopt the Gaussian antenna gain model suggested by the 3GPP to represent the realistic directional antenna pattern and also assume all transmitters and receivers are equipped with the same directional antenna. Considering the truncated Gaussian and uniform BAE distributions, we derive the exact probability density functions (PDFs) of the cascaded antenna gains of a TxRx pair. Taking into account of the distribution characteristics of BAE and the strong directivity of mmWave antenna, we also provide the approximated PDFs for tractability.

With the help of the Fox’s H function and its series expansion, we provide a novel solution to derive the analytic expression of the energy coverage probability of the mmWave WPT in the presence of BAE. Besides, the closedform expression of average harvested RF energy is derived. We define the relative energy loss (REL) in order to further investigate the performance degradation incurred by BAE.

Through MonteCarlo simulations, we verify the derived energy coverage probability and average harvested RF energy. It is found that the widely used flattop antenna gain model cannot always exactly evaluate the performance of the mmWave WPT system in the presence of BAE. We also conclude that BAE indeed lowers the energy coverage probability and the average harvested RF energy.
Notation: , , and represent the sets of all integers, real numbers and complex numbers, respectively.
is the expectation of random variable (r.v.)
and is the probability of event . For r.v. , andstand for the probability density function (PDF) and cumulative distribution function (CDF) of
, respectively. For , is the modulus of . Given , means the Euclidean norm of. For a complex vector
, and stand for transpose and conjugate transpose of , respectively. means the r.v.follows Gamma distribution with PDF
, where is the Gamma function. is the Gaussian error function and is the Dirac delta function. is the natural logarithm of . For and , represents the binomial coefficient and equals to . For r.v. , is the Laplace transform of .Ii System Model
Iia Network and Node Models
In this paper, we consider a mmWave wirelesspowered network which consists of two types of nodes, i.e., Energy Transmitter (ETx) and Energy Receiver (ERx) (See the Fig.1 on the next page). The ETxs are connected to stable power sources, thus having the capability of emitting energy signals. While, the ERxs have to harvest energy to maintain their routine operation and information transmission. The location of all ETxs follows the homogeneous Poisson point process (HPPP) with intensity on the twodimensional Euclidean plane. Suppose a saturated service scenario, where each ETx is assumed to serve one dedicated ERx in a WTP block [19, 12]. We study the performance of the typical ERx located at the origin. Let the ETx associated by the typical ERx be located at and .^{1}^{1}1Herein can be seen as the maximum allowable WPT distance for an ETxERx pair. Accordingly, in this work we investigate the worst performance of the considered mmWave WPT system [19, 12]. Then, due to Slivnyak’s Theorem [19], except the ETx at , the location of other ETxs still forms a PPP with the same intensity , denoted by . Without loss of generality, the WPT duration is assumed as unit time. It means in following context, the harvested energy and the harvested power have an equivalence.
IiB Channel Model
For mmWave links, lineofsight (LOS) and nonlineofsight (NLOS) channels have sharply different propagation characteristics. Given a propagation distance , the path loss of a LOS mmWave link can be modeled as , while in the NLOS case, the path loss is . Here ( ) is the LOS (NLOS) path loss exponent and () represents the intercept of LOS (NLOS) link. In general, we have and [15, 3]. Moreover, we model the smallscale fading of each mmWave link as independent Nakagami fading with parameters and for LOS and NLOS scenarios, respectively. Thus, we denote and as the smallscale channel gains for LOS and NLOS links, respectively, and then express the LOS power gain as and NLOS power gain as .
The probability of a mmWave link experiencing LOS is modeled as a function with respect to the distance from Tx to Rx [31]. Several empirical and analytical blockage models had been reported in [31] and [3]. Considering tractability and accuracy, we employ the random shape theory model, by which the analytical LOS probability is coincident with the empirical expression of the 3GPP blockage model [3]. Hence, given the distance , the LOS and the NLOS probabilities can be separately expressed as and , where the blockage parameter is determined by the average number and average perimeter of buildings in the interested area. Similar to most literature, e.g., [Bai and Heath [5]][Comisso and Babich [10]], the correlation of LOS probabilities for different mmWave links is ignored. Consequently, the ETxs in can be divided into two independent sets, namely, LOS ETx set and NLOS ETx set. Both of which follow the inhomogeneous PPP (IPPP) with intensity and , denoted as and , respectively. As the ETx at has the distance to the typical ERx, the performance analysis in the LOS scenario is very similar to that in the NLOS scenario. Furthermore, considering the harvested energy plays a crucial role in maintaining the lifetime of wirelesspowered nodes, we assume the ETx only selects the ERx experiencing LOS as its target receiver like [46, 12, 43].
IiC Antenna Pattern
We assume that all ETxs and ERxs are equipped with the same directional antenna to perform mmWave beamforming. It is also assumed that all ETxs have the same transmit power . For the typical ERx, the received RF energy transmitted by the ETx at can be given by
(1) 
where and stand for the antenna gains of ETx and ERx, respectively. and are the orientation angles relative to the boresights of ETx and ERx, respectively, and belong to the interval . Note that we ignore the energy harvested from the noise, as it is trivial compared with the received RF energy [37].
In order to theoretically evaluate the systemlevel performance, it is inevitable to depict the PDF of the cascaded antenna gain . For the uniform linear array (ULA), the Fejér kernel based sinc and cosine antenna patterns are employed to represent the directional antenna gain [46, 2, 30, 4]. Unfortunately, it is difficult to directly derive the exact analytic expression of the PDF of using the sinc or cosine antenna gain model. In [46, 30, 4], the receiver was assumed to be equipped with omnidirectional antenna to avoid the cascaded directional antenna gain. Although [2] considered a cascaded antenna gain with the sinc antenna pattern, the nodes were not stochastically deployed.
Alternatively, the flattop antenna gain model is mostly used in stochastic geometry coverage analysis for its mathematical tractability, e.g., [21, 38, 51, 40, 42, 8, 12, 18, 41, 22, 43, 37]. In the perfect beam alignment or slight BAE scenarios, the flattop antenna gain model can provide the tractable theoretical expression of the systemlevel performance with acceptable accuracy. However, the flattop antenna model is lack of the capability of depicting the rolloff effect of the mainlobe and incurs significantly inaccuracy for evaluating the systemlevel performance in the IBA scenario [46, 34, 35, 44, 32]. It is worth mentioning that in [13] a generalized flattop model was proposed to approximate the practical antenna gain. Nevertheless, in the IBA case, it could incur relatively high complexity to derive the PDF of the cascaded antenna gain.
In [34, 35], a Gaussian antenna pattern, i.e.,
where is the maximum mainlobe gain, is the sidelobe gain and is determined by the 3dB beamwidth, was used to represent the mmWave directional antenna gain. Nevertheless, to obtain the analytic PDF of the cascaded antenna gain , the side lobe was ignored in the analysis of [34, 35]. Moreover, the loss in ergodic capacity derived by [34, 35] only involves the truncated Gaussian BAE. Differently, in the nonstochastic mmWave networks [44, 49, 32] adopted the 3GPP Gaussian antenna model, which can depict the rolloff characteristic and match the measurement well. In this paper, therefore, we employ the 3GPP Gaussian antenna model, i.e.,
(2) 
where
, is the main lobe (20dB) beamwidth, and is the halfpower (3dB) beamwidth. Since , the continuity is ensured. According to the practical measurement reported by [36, 44], when , is approximately equal to . Then, with this empirical approximation, we further obtain , , and . For convenience, we introduce the normalized antenna gain , i.e.,
(3) 
in which and . Obviously, with the empirical expression , the Gaussian antenna pattern is only determined by . We herein employ (3) to reduce the parameter number of the antenna radiation pattern. Besides, the cascaded normalized antenna gain is denoted by .
IiD Imperfect Beam Alignment Models
According to (3), if we intend to maximize the harvested energy, we can let and . In practice, however, and are not necessarily equal to zero due to the direction estimation error and hardware imperfection [42, 8, 4].
For the associated ETxERx pair, as the BAE appears after the beam aligning procedure, and
are usually modeled as independent and identically distributed truncated Gaussian random variables with zero mean and standard deviation
[34, 8, 4]. Such that, the PDF of can be expressed as(4) 
which is named as the Gaussian BAE model. The standard deviation is usually used to indicate the variability of BAE. The larger is, the stronger statistical dispersion the BAE exhibits, which means the BAE becomes more severe. Observing (4), if , gradually converges to . It is corresponding to the fact that if the beam is perfectly aligned, the antenna gain equals to with probability 1.
For the nonassociated ETxERx pairs, an ETx in or just aligns its beam with the boresight of its paired ERx. To facilitate understanding and representation, for the typical ERx, given the ETx at , , we can also treat and
as BAEs, which are usually modeled as independent uniform distribution over
[46, 8, 4], i.e.,(5) 
It is named as the uniform BAE model. Note that it is assumed that BAEs at ETxs and ERxs are independently distributed [34, 8, 42].
IiE Energy Harvesting Model
As the smallscale gains of different mmWave links are independently distributed, the harvested RF power of the typical ERx can be written by
(6) 
where is the harvested power from the ETxs in , is the harvested power from the ETxs in , and is the harvested power from the ETx located at . Then, the harvested direct current (DC) power at the typical ERx is
(7) 
Note that is the RFDC power conversion function. In practice, is a nonlinear function with respect to the input RF power [37, 30, 7]. Using the practical nonlinear energy harvesting model proposed in [7], we can write the harvested DC power as
(8) 
where is the maximum DC power that can be harvested by the ERx and and are the constants determined by the rectifier circuit [7].
Iii The PDFs of Antenna Gains
In this section, we derive the PDFs of the normalized antenna gains with the BAE following the truncated Gaussian or uniform distributions. Then, the PDFs of the cascaded antenna gains with two BAE models are derived. For tractability, we also provide the approximated PDFs of the cascaded antenna gains.
Iiia The PDFs of the Normalized Antenna Gains
Lemma 1.
If the PDF of the stochastic BAE is , the PDF of the normalized Gaussian antenna gain is given by
(9) 
where and .
Proof:
See Appendix A. ∎
Corollary 1.
If the BAE
follows the truncated Gaussian distribution with zero mean and variance
over , the PDF of can be expressed as(10) 
where and .
Thus, for the associated ETxERx pair, the PDF of the normalized antenna gain of ERx or ETx is equal to (10).
Corollary 2.
If the BAE follows the uniform distribution over , the PDF of can be written by
(11) 
where and .
For the nonassociated ETxERx pairs, the PDF of involved antenna gains is illustrated by the Corollary 2.
IiiB The PDFs of the Cascaded Antenna Gains
By (1), the cascaded antenna gain plays a crucial role in the system performance. Herein, we intend to derive the PDFs of with both BAE models respectively. To this end, we have following two theorems.
Theorem 1.
The PDF of the cascaded normalized antenna gain with truncated Gaussian BAE model can be written as (12) at the top of the next page.
(12) 
Proof:
See Appendix B ∎
Theorem 2.
The PDF of the cascaded normalized antenna gain with uniform BAE model can be written as (13) at the top of the next page.
(13) 
Proof:
The proof of Theorem 2 is similar to that of Theorem 1, therefore, we omit the proof for clarity. ∎
IiiC The Approximated PDFs of the Cascaded Antenna Gains
Although the Theorem 1 and 2 show the exact PDFs of with the Gaussian and uniform BAE models respectively, the arctan functions in (12) and (13) make further analysis less tractable. Therefore, we provide two approximated PDFs, by considering the distribution characteristics of BAE models and the strong directivity of mmWave antenna.
IiiC1 The Approximated PDF With the Gaussian BAE model
In the mmWave WPT system, the associated ETxERx pair is expected to employ the elaborately designed beam alignment algorithms, such as [48], to minimize and as much as possible. It is therefore reasonable to infer that and are less than in most cases. For instance, it is straightforwardly assumed by [28] that and lie in . With the same consideration, the authors in [34] and [35] ignored the cascaded antenna gain involving the sidelobe gain, because it has the relatively small value and happens in a very low probability in the Gaussian BAE scenario. Following these works, we also ignore sidelobe gain in the cascaded antenna gain, i.g., the component of (12) over . As a result, with the Gaussian BAE model can be approximately presented as
(14) 
When , we have .
IiiC2 The Approximated PDF With the Uniform BAE model
Recalling the expression of , the strong directional antenna means is far less than . Hence, in the uniform BAE case, the event that equals to the product of two mainlobes occurs in an extremely small probability, i.e., . Furthermore, by (40) in the Appendix B, it can be inferred that the arctan term in (13) is generated by the product of two mainlobes. Therefore, it is reasonable to neglect the arctan term of (13) in the uniform BAE case. Then, we can approximate (13) by
(15) 
IiiC3 Verification of the Approximated PDFs
To verify our approximations, we draw with Gaussian and uniform BAE models in Fig. 2 and Fig. 3, respectively. For the Gaussian BAE model, we draw the PDF of with as an example to verify our approximation.^{2}^{2}2, the approximated PDF of shows the similar accuracy with various . From the left subfigure of Fig. 2, we can see that with and , the approximated PDFs match the theoretical PDFs closely. As decreases, the PDF of tends to be a pulselike function. We can infer that it asymptotically converges to as . While, in the right subfigure of Fig. 2, for , there is a slight difference between the results of (12) and (14). It is because the probability that the sidelobe appears in the cascaded antenna gain becomes larger when grows. Thus, with , the difference generated by our approximation seems extremely apparent. Accordingly, if , (14) is an appropriate approximation of (12). ^{3}^{3}3 Given , there is always for . While, given , we have for . Consequently, the probabilities of two mainlobes cascading with and , i.e., , are about 0.9946 and 0.9111, respectively. Apparently, for , , there is . That is why we here choose as the approximation condition. Note that it is weaker than the assumption that adopted by [28]. Moreover, from Fig. 3 we can see that the curves of (15) approach those of (13) extremely closely. Summarily, it is verified that both approximated PDFs can be used to analyze the systemlevel performance instead of the exact PDFs under our considered circumstance.
Iv Energy Coverage Analysis
In this section, we focus on analyzing the energy coverage probability of the typical ERx. Energy coverage probability is defined as the probability that the harvested DC energy is larger than a predefined threshold, which is always the minimum required energy for information transmission or other operations.
In [33], the Meijer Gfunction was used to derive the analytic expression of information coverage probability of mmWave transmission. Theoretically speaking, we can also adopt this successful approach in our analysis. By [29], however, to obtain the analytic expression in the form of the Meijer Gfunction, the path loss exponents should be positive integers, which limits the application of the Meijer Gfunction based method. Alternatively, in this paper, without path loss exponent limitation, we provide an analytic expression of energy coverage probability with the help of Fox’s H function.^{4}^{4}4The Fox’s H function is a general function which can encompass almost all commonly used functions, e.g., Meijer Gfunction. Although the Fox’s H function is defined by an integral in a nonanalytic form, like the widely used Gamma function, Qfunction, Hypergeometric function, Meijer Gfunction, etc., a lookuptable (LUP) storing the values of Fox’s H function can be generated via numerical methods. A Matlab program for evaluating the Fox’s H function was provided in [26]. More details about the Fox’s H function can be found in [24].
Letting be the DC energy threshold, we can write the energy coverage probability of the typical ERx as
(16) 
where is the equivalent RF energy threshold. In addition, according to , we let to observe the behavior of . If , there is . Note that , and are independent random variables. Following the widely adopted Gamma r.v. approximation [21, 37], i.e., using , , instead of , we can rewrite as
(17) 
in which , and . Next, we derive the analytic expressions of , , and , respectively.
Iva The Analytic Expression of
IvB The Analytic Expression of
Define , . As follows the IPPP with intensity , the 2tuple forms a marked IPPP (MIPPP). Due to the probability generating functional (PGFL) of MIPPP [19], we have
(20) 
In , due to the Fubini’s theorem [6], we exchange the order of integral and expectation operations. Due to [Mathai and Haubold [25], (1.9.5)], we have
(21) 
Note that is the Fox’s H function and is defined by [24, 25]. Consequently, there is
(22) 
Before solving the expectation of with respect to , we introduce the Lemma 2 as follows.
Lemma 2.
For , there is
Proof:
See [Mathai et al. [24], (1.88)]. ∎
By the Lemma 2, we can further attain
In , we apply Binomial theorem. Hence, there is
(23) 
where , .
Proposition 1.
If the PDF of follows (15), , there is
(24) 
Proof:
By (15), it is straightforward to obtain Proposition 1. The detailed proof is omitted to save space. ∎
Proposition 2.
If , , , there is
(25) 
Proof:
By the PDF of , it is conveniently to prove Proposition 2. ∎
Therefore, substituting Proposition 1 and 2 into (23), we can obtain the analytic expression of .
IvC The Analytic Expression of
Firstly, after some manipulations, it is easy to know
(27) 
where , , which can be also written in the analytic form due to Proposition 1 and 2. Following the derivations from (21) to (23), similarly, we can obtain
(28) 
Consequently, the analytic expression of is also derived. Therefore, we can get the analytic expression of with , , and .
It is worth noting that for applying Proposition 1 and 2, we need to let . Clearly, it always holds for the piratical condition . So, we put no limitation on the path loss exponents.
V Average Harvested Energy
Although the derived can be used to evaluate the energy coverage performance of mmWave WPT, it may not provide explicit and direct insight into the effect of BAE. Thus, in this section, the average harvested energy is addressed to further investigate the effect of BAE.