Modeling and Performance Analysis of Large-Scale Backscatter Communication Networks with Directional Antennas

Backscatter communication (BackCom) constitutes intriguing technology that enables low-power devices in transmitting signals by reflecting ambient radio frequency (RF) signals that consume ultra-low energy. Applying the BackCom technique in large-scale networks with massive low-power devices can effectively address the energy issue observed in low-power devices. Prior studies only consider large-scale BackCom networks equipped with omni-directional antennas, called Omn-BackCom Net. To improve the network’s performance, we employ directional antennas in large-scale BackCom networks, called Dir-BackCom Nets. This article establishes a theoretical model for analyzing the performance of Dir-BackCom Nets. The performance metrics include both connectivity and spatial throughput. Our model is genaralized for both Dir-BackCom Nets and Omn-BackCom Net. The accuracy of our theoretical model is verified by extensive simulations. Results indicate that Dir-BackCom Nets can improve connectivity and spatial throughput. Moreover, results show that the throughput can be maximized by choosing an optimal density of BTs. In addition, both the connectivity and spatial throughput of BackCom Nets can be improved by choosing a directional antenna with a proper beamwidth and gain of the main lobe. Our theoretical model and results can offer beneficial implications for constructing Dir-BackCom Nets.


Introduction
In the era of the Internet of things (IoT), there is an expectation that massive low-power devices will be used in large-scale IoT networks to monitor, sense, and generate information to support various IoT applications such as smart city/factory/agriculture [1][2][3][4][5]. Powering such massive devices is challenging. As a promising technique, the backscatter communication (BackCom) technique enables devices to transmit data by reflecting environmental radio frequency (RF) signals with consuming ultra-low energy. Therefore, deploying the BackCom technique into large-scale IoT networks can effectively overcome the battery issue of low-power devices [6].
In large-scale BackCom networks (BackCom Nets), backscatter transmitters (BTs) need to transmit data by reflecting RF signals from ambient energy/signal providers: carrier emitters (CEs) [7]. Consequently, the transmission of BTs is dependent on the location relationship between BTs and CEs. As a result, the communication performance of BackCom is heavily dependent on the location distribution of both CEs and BTs [6][7][8]. Therefore, it is necessary to investigate the performance of large-scale BackCom Net and analyze the relationship between the distribution of both CEs and BTs and the performance of BackCom Nets.
In the performance analysis of BackCom networks, because CEs and BTs are usually distributed with randomness in large-scale networks [6][7][8][9][10], both the energy received at BTs and interference received at a backscatter receiver (BR) are usually variables. A popular approach for analyzing the energy/interference is to utilize a mathematical tool known as receive RF signals first and then reflect signals). Therefore, applying directional antennas to a large-scale BackCom Net can potentially improve network performance. In addition, considering that directional antennas have been widely deployed in 5G communication devices and WiFi devices with 802.11ad to compensate for the high attenuation of the path loss of MMwave [19][20][21], it is worth developing a theoretical model to analyze the performance of large-scale BackCom Nets in which BTs and BRs are equipped with directional antennas, which are called Dir-BackCom Nets.
To the best of our knowledge, there is no study that analyzes the performance of large-scale Dir-BackCom Nets. The performance analysis of large-scale Dir-BackCom Nets is far more complex than large-scale BackCom Nets mounted with omni-directional antennas, called Omn-BackCom Nets. This is because the transmission/reception strength of a directional antenna varies in different directions (i.e., it is strong in some directions and weak in some directions), leading to the fact that the received energy/interference at a directional antenna in a large-scale network is spatially inhomogeneous. For example, Figure 1 illustrates a part of a Dir-BackCom Net around a transmission pair BT 0 -BR 0 , where BT 0 receives the energy (i.e., RF signals) provided from CEs and then reflects signals to BR 0 . We can observe that BT 0 receives energy from CE 1 and CE 2 by different antenna gains, thereby leading to inhomogeneous energy reception in the space around BT 0 . Then, BR 0 receives interference by different antenna gains (i.e., it receives interference from BT 1 and BT 4 by a high gain and receives interference from BT 2 and BT 3 by a low gain). In addition, the interference from BT 1 and BT 3 is transmitted by a high gain, while the interference from BT 2 and BT 4 is transmitted by a low gain. As a result, the received interference at BR 0 from all BT 1 to BT 4 is different from each other, resulting in the fact that the interference received at BR 0 is more spatially inhomogeneous. Therefore, compared with Omn-BackCom Nets where the received energy/interference is spatially homogeneous, the complexity of the performance analysis in Dir-BackCom Nets significantly increased. Therefore, the main goal of this paper is to develop a theoretical model to analyze the performance of Dir-BackCom Nets. The main contributions are summarised as follows. The rest of the article is organized as follows: Section 2 presents system models. Then, we analyze the connectivity of Dir-BackCom Nets in Section 3. Next, Section 4 derives the network throughput of Dir-BackCom Nets. The simulation and numerical results are presented in Section 5. Finally, the article is concluded in Section 6.

Network Model
We consider a large-scale Dir-BackCom network as shown in Figure 2. In the network, CEs follow a homogeneous Poisson point process (HPPP) denoted by Θ CE = {CE i , i = 1, 2, . . .} with density λ CE ; BTs follow another independent HPPP denoted by Θ BT = {BT j , j = 1, 2, . . .} with density λ BT , where λ BT λ CE [6]. Then, we consider that each BT is associated with a BR, where the location orientation of the BR is uniformly distributed in [0, 2π] around its paired BT [18]. We assume that each BT/BR is mounted with a directional antenna, where the orientation of the directional antenna of a BT/BR is directed towards its paired BR/BT [18]. As a result, both orientations of the directional antennas of BTs and BRs are uniformly distributed in [0, 2π]. In addition, we consider that each CE is mounted with an omnidirectional antenna to provide RF signals for all BTs [6]. We assume that a BT can receive RF signals from all CEs as energy [6]. When the received energy at a BT is higher than a threshold, BT can be activated from an idle state and conduct backscatter communications.
It is worth noting that our performance analyses do not consider the devices near network border, because their performance analysis is quite complex due to the fact that their performance is related to the specific shape of the network's edge [22]. Therefore, similarly to [6][7][8][9][10]15,18], this article only considers the communication performance of devices that are not affected by border effects.

Channel Model
In the large-scale Dir-BackCom Net, multiple orthogonal frequency channels can be allocated to CEs. Because available frequency channels are limited while there are massive CEs in the large-scale network, multiple CEs reuse an identical frequency channel. Assume that the number of available frequency channels is N. Then, the density of co-channel CEs is λ CE /N. In addition, we assume that BTs can randomly choose one frequency channel to receive RF signals from CEs. Consequently, the density of co-channel BTs is λ BT /N. During propagation, we consider a typical propagation model r −α hG t G r , where r is the propagation distance, α is the path loss exponent with 2 < α < 6 [18], h is the small-scale fading factor following an exponential distribution with mean 1 [6,18], and G t and G r are the antenna gains of the transmitters and receivers, respectively.

Antenna Model
The antenna gain of an omni-directional antenna, denoted by G o , is G o = 1. The radiation pattern of a realistic directional antenna is shown in Figure 3 [16,[23][24][25], which comprises one main lobe and several side/back lobes. We can observe that each lobe of the radiation pattern varies in directions. This leads to the intractability of the performance analysis in large-scale BackCom Nets. For analytical tractability, similarly to [17,[26][27][28], we simplify the realistic directional antenna to a keyhole antenna model, as shown in Figure 3. The radiation pattern of the keyhole antenna model consists of a main lobe with gain G m in beamwidth θ m and a constant side lobe with gain G s in beamwidth θ s . The antenna gain of the keyhole antenna model, denoted by G d , can be described by the following equation: where the detailed derivation process of the antenna gain is provided in Appendix A. To assure the generality of our model, we consider that directional antennas equipped by BTs and BRs can have different values of both G m and θ m . Therefore, we let θ tm and θ ts denote the beamwidth of the main lobe and the side lobe of the antenna equipped by BTs, respectively. Then, let G tm and G ts denote the corresponding gains within θ tm and θ ts , respectively. Similarly, we let θ rm and θ rs denote the beamwidth of the main lobe and the side lobe of the antenna equipped by a BR, respectively. Then, let G rm and G rs denote the corresponding gains within θ rm and θ rs , respectively. Figure 4 shows the adopted backscatter communications model [6,8,9]. RF signals generated from CEs are used for both energy harvesting and signal transmission at BTs. Specifically, on the one hand, RF signals are harvested as energy at BTs to run the circuits of backscatter communications; on the other hand, RF signals are reflected to BRs for information transmission. The reflected signals are modulated by a micro-controller: RF signals are reflected in different power levels to stand for various bits when the microcontroller switches to various impedances [29]. The ratio of the reflected power to the total received power at BTs is named as reflection coefficient and is denoted by η. Assume that the received power at a BT is denoted by P r BT , the reflected power at BTs, denoted by P t BT , and the power for energy harvesting at BTs, denoted by P h BT , can be expressed as the following two equations.

Backscatter Communication Model
RF signals used for energy harvesting can be stored at BTs by a energy harvester. Following [30][31][32], we adopt a non-linear energy harvesting model proposed in [33]. The harvested energy, denoted by P e (P h BT ), can be expressed as follows: where E max is the saturated energy power that can be stored at the energy harvester maximally; C is a constant to ensure a zero-input or zero-output response at the energy harvester, c 1 and c 2 are two constants related to the specific capacitance, resistance, and diode turn-on voltage in energy harvesting circuits. The term Ω(P h BT ) is a function of P h BT , which is given by the following.
The harvested energy of BTs (P e (P h BT )) can be used to activate backscatter communications. Let P c denote the power threshold for activating the circuits. Then, the trigger condition for backscatter communications is provided by the following.

BR CE BT
Energy harvester

Microcontroller
Variable impedance

Connectivity
The connectivity of BackCom networks is defined as follows: Definition 1. The connectivity is the probability that a BT-BR pair can establish a link successfully.
The condition of a successful link between a BT-BR pair is as follows: 1.
The BT can receive sufficient power to activate backscatter communications.

2.
The signal reflected from the BT can be successfully received by the BR.
We call the probability that a BT can be activated as the active probability. We first analyze the active probability in Section 3.1. Then, we analyze the connectivity in Section 3.2.

Active Probability
The backscatter communications at a BT can be activated when the harvested energy of the BT is higher than a threshold P c (i.e., P e ≥ P c ). Since the harvested energy at a BT is dependent on the received power at a BT, we first analyze the received power at a BT (i.e., P r BT ). Figure 5 provides an example for illustrating the energy reception at a BT. We can observe that a BT receives energy by both the main lobe within θ tm and the side lobe within θ tm . Therefore, we divide the surrounding space of a BT into two regions:

•
The region in which a BT receives signals by the main lobe (the shadowed region in Figure 5), named region A tm . • The region in which a BT receives signals by the side lobe, named region A ts .
Let P A tm BT and P A ts BT denote the received power at a BT in regions A tm and A ts , respectively. The total received power at a BT can be given by P r BT = P A tm BT + P A ts BT . We let Θ CE (A tm ) and Θ CE (A ts ) denote the HPPPs of CEs providing energy for a BT in regions A tm and A ts , respectively, where Θ CE (A tm ) and Θ CE (A ts ) are mutually independently. If a PPP is distributed in a region, the PPPs distributed in non-lapping sub-regions are mutually independent [34]; let P CE denote the transmitted power of CEs. Then, the power of received energy at a BT can be provided by the following: where r i is the distance between CE i and the studied BT, h i is the small-scale fading factor for the path between CE i and the studied BT, and G ts = Figure 5. Example to illustrate the energy reception at a BT.
Then, we provide both the cumulative distribution function (CDF) and the probability density function (PDF) of the received power at a BT (P r BT ) as the following proposition.

Proposition 1.
The CDF and the PDF of the received power at a BT (P r BT ), denoted by F P r BT (x) and f P r BT (x), respectively, are given by the following: where the following is the case. W = Proof. The proof is given in Appendix B.
Then, we provide the active probability, denoted by P a , as the following proposition.

Proposition 2.
The active probability of a BT is given by the following: and W is a function of both G tm and θ tm given by Equation (9).
Proof. The proof is provided in Appendix C. (10) is general for BTs that are equipped with directional antennas or omni-directional antennas. When we set G tm = 1 and θ tm = 2π, we can have the active probability of a BT equipped with an omni-directional antenna. We can observe from Equation (10) that the active probability is a function of the gain of the main lobe of BTs (G tm ), the beamwidth of the main lobe of BTs (θ tm ), the density of CEs (λ CE ), and the path loss exponent (α).

Connectivity
For a BT-BR transmission pair, we consider that the signal transmitted from a BT can be received successfully when the signal-to-interference-noise-ratio (SINR) at the BR is not less than a threshold δ. Given the distance between the BT-BR pair by R 0 , the connectivity of the BT-BR pair is provided by the following: where , h 0 is the small-scale fading factor for the path between the studied BT and the studied BR, and σ 2 is the power of noise.
Next, we derive the Laplace transforms L I (B). We first analyze the interference reception at a BR (i.e., I). Figure 6 provides an example to illustrate the interference received at a BR. Because the directional antenna of a BR receives interference by two antenna gains (i.e., the gain of the main lobe and the gain of the side lobe), we divide the surrounding region of the BR into two regions: • The region where a BR receives interference by the main lobe within θ rm (the shadowed region in Figure 6), named region A rm ; • The region where a BR receives interference by the side lobe within θ rm , named region A rs . Figure 6. Example for illustrating the interference reception at a BR.
Let I A rm and I A rs denote the received interference in regions A rm and A rs , respectively. The received interference at a BR can be given by I = I A rm + I A rs .
We first analyze the received interference from region A rm (i.e., I A rm ). The received interference in region A rm can be split into two parts: (1) the interference transmitted by the main lobe of BTs (such as the BT 1 in Figure 6), where the density of this type of BTs (pointing the BR by their main lobe) is θ tm /(2π) · (λ BT /N); (2) the interference transmitted by the side lobe of BTs (such as BT 2 in Figure 6), where the density of this type of BTs (pointing the BR by their side lobe) is θ ts /(2π) · (λ BT /N). Because the orientation of directional antennas of BTs is mutually independently distributed, the BTs pointing a BR by their main lobes and the BTs pointing a BR by their side lobes are two independent HPPP distributions, denoted by Θ BT (A m rm ) and Θ BT (A s rm ), respectively. Let BT 0 denote the BT paired with the BR being studied. Then, the received interference at the studied BR in region A rm can be expressed as follows: where r j is the distance between BT j and the studied BR, and h j is the small-scale fading factor for the path between BT j and the studied BR. Similarly, the received interference at the BR from region A rs can be expressed as follows: where Θ BT (A m rs ) and Θ BT (A s rs ) are HPPPs that BTs pointing the BR in region A rs by their main lobes and side lobes, respectively. The two HPPPs Θ BT (A m rs ) and Θ BT (A s rs ) are independent, and G rs = With the received interference I = I A rm + I A rs , we then provide the Laplace transform L I (B) as the following Lemma. where , and P a is given by Proposition 2.
Proof. The proof is provided in Appendix D.
Finally, the connectivity of large-scale Dir-BackCom Nets can be given by the following Theorem.
where M n = ∞ 0 (1 − Φ n )r j dr j for n = {1, 2, 3, 4}. The term Φ n (n = 1, 2, 3, 4) is given by α dγ is given by Proposition 1. (14) and (9) into Equation (11), we can obtain the connectivity given by Equation (15). (15) is general for large-scale BackCom Nets equipped with omn-directional antennas or directional antennas. Specifically, when G tm = G rm = 1 and θ tm = θ rm = 2π, we can have the connectivity of large-scale Omn-BackCom Nets. We can observe from Equation (15) that the connectivity is a function of the density of BTs (λ BT ), the antenna gain of main lobe of BTs and BRs (G tm and G rm ), the beamwidth of the main lobe of BTs and BRs (θ tm and θ rm ), and path loss exponent α. The specific relationship will be investigated in Section 5.

Network Throughput
Considering our BackCom Net is a large-scale network, we adopt the spatial throughput as the metric of the network throughput [18,35]. The spatial throughput of BackCom Nets is defined as follows.
Definition 2. The spatial throughput of BackCom Nets is the throughput of BTs within a unit area of the network space.
The spatial throughput of large-scale Dir-BackCom Nets can be expressed by the following.
Theorem 2. The spatial throughput of large-scale Dir-BackCom Nets is given by the following: where M n = ∞ 0 (1 − Φ n )r j dr j for n = {1, 2, 3, 4}, and Φ n is given by Lemma 1. The term f P r BT (x) is given by Proposition 1. Note that term B in Equation (17) needs to be replaced by B = (e t − 1)/(ηP r BT R −α 0 G tm G rm ).
Proof. The proof is provided in Appendix E.

Simulations and Numerical Results
We conduct simulations to verify our proposed analytical model. Specifically, we adopt commercial software MATLAB as our simulation tool. In simulations, fixed parameters/factors are provided in Table 1. Then, simulations are conducted in two groups: In Group 1, we make a comparison of performance among different BackCom Nets, such as Dir-BackCom Nets, Dir-Omn-BackCom Nets, and Omn-BackCom Nets; in Group 2, we make a comparison of performance among Dir-BackCom Nets with different directional antennas having various parameters. Our simulation results are obtained by the Monte Carlo simulation averaged by 3000 topological trials. In the following figures, the abbreviations 'ana' and 'sim' in legends stand for 'analysis' and 'simulation', respectively.

Parameters Values
The number of available frequency channels (N) 5 The transmitted power of CEs (P CE ) 19.5 dBm [7] The reflection coefficient (η) 0.375 [7,8] The power of noise (σ 2 ) −60 dBm [7,8] The threshold of SINR (δ) 5 dB [7] The power threshold for activating the circuits of backscatter communications (P c ) 10.6 µW [36] The saturated energy power at the energy harvester (E max ) 48.86 µW [37] The factors in the energy harvester (c 1 , c 2 ) 26515.46, −0.00002981 [37] 5.1. Comparison among Dir-BackCom Nets, Dir-Omn-BackCom Nets, and Omn-BackCom Nets Table 2 gives the difference among Dir-BackCom Nets, Dir-Omn-BackCom Nets, and Omn-BackCom Nets. We can see in Table 2 that in both Dir-BackCom Nets and Dir-Omn-BackCom Nets, BTs are equipped with directional antennas, while in Omn-BackCom Nets, BTs are equipped with omni-directional antennas. Then, we first compare the active probability in the network in which BTs are equipped with a directional antenna or an omnidirectional antenna. Figure 7 shows the active probability of a BT versus the density of CEs under different path loss exponents in the network in which BTs have different antennas. We can observe that the active probability of a BT with a directional antenna is always higher than that with an omni-directional antenna under different network environments (i.e., α = 3 and α = 4). This phenomenon indicates that compared with the omni-directional reception, directional antennas can receive more RF signals by beamforming reception with a higher gain, implying that employing directional antennas at BTs can improve the active probability. In addition, we can observe that the active probability increases with the increment of the density of CEs and decreases with the increment of the path loss exponent. This is because BTs can receive more RF signals when there are more CEs or when the path loss is low in a network environment.

Dir-BackCom Nets Dir-Omn-BackCom Nets Omn-BackCom Nets BTs Directional antennas Directional antennas Omni-directional antennas
BRs Directional antennas Omni-directional antennas Omni-directional antennas Figure 8 shows the connectivity among Dir-BackCom Nets, Dir-Omn BackCom, and Omn-BackCom Nets under different density of BTs and various path loss exponents. We can observe that the connectivity of Dir-Omn-BackCom Nets is always higher than Omn-BackCom Nets under different path loss environments (i.e., α = 3 and α = 4). From Figure 7, we know that BTs with directional antennas can have a higher active probability than BTs with omni-directional antennas. Then, the higher active probability of BTs with directional antennas can increase the quantity of transmitting BTs, thereby increasing the overall interference of the network. However, Figure 8 shows that compared with omni-directional antennas, using directional antennas at BTs can still improve the network connectivity, although with higher interferences. The reasons of this phenomenon are as follows: (1) BTs with directional antennas can have higher active probabilities; (2) the beamforming with a higher gain also can strengthen transmitted signals, thereby leading to a higher SINR. In addition, we can observe in Figure 8 that the connectivity of Dir-BackCom Nets is higher than Dir-Omn-BackCom Nets under different values of α, indicating that using directional antenna at BRs can also improve the network's connectivity. This is because the narrower beamforming with a higher gain of directional antennas can receive a higher SINR by both augmenting the signal strength in the intended direction and reducing the interference in undesired directions. Therefore, we can conclude that using directional antennas at BTs or BRs can both improve the connectivity of BackCom Nets. In addition, we can observe that the connectivity decreases with the increment of the density of BTs. This is because more BTs lead to higher interferences in BackCom Nets.    Figure 9 shows the spatial throughput among Dir-BackCom Nets, Dir-Omn Back-Com, and Omn-BackCom Nets under different path loss exponents. We can observe that the spatial throughput of Dir-BackCom Nets is significantly higher than both Omn-Dir-BackCom and Dir-BackCom Nets, indicating that equipping directional antennas at BTs and BRs instead of omni-directional antennas can greatly improve the spatial throughput of BackCom Nets. Moreover, we can observe in Figure 9 that there is an optimal density of BTs to maximize the spatial throughput. The optimal density of BTs of Dir-BackCom Nets is higher than both Dir-Omn-BackCom Nets and Omn-BackCom Nets, indicating that Dir-BackCom Nets can accommodate more BTs to obtain a maximum throughput in a unit area compared with Dir-Omn-BackCom Nets and Omn-BackCom Nets. This phenomenon also indicates that we can choose a proper density of BTs to improve the spatial throughput according to different BackCom Nets.

Comparison among Dir-BackCom Nets with Different Directional Antennas
In this group, we investigate the impact of different directional antennas on the performance of Dir-BackCom Nets. Specifically, we consider that these different directional antennas possess different beamwidths. Meanwhile, similarly to [38], the considered directional antennas with a narrower beamwidth of the main lobe have a higher gain of the main lobe. Then, we consider that for different Dir-BackCom Net, directional antennas of BTs have the following varied parameters: (θ tm , G tm ) = {(2.13 rad, 4), (1.75 rad, 6), (1.35 rad, 10)}. In addition, the antenna parameters of BRs are fixed as (θ rm , G rm ) = (1.75 rad, 6). It is worth noting that we do not investigate the impact of antenna beamwidth of BRs on performance, because BRs have a narrower beamwidth and the higher gain of the main lobe can obviously lead to a better network performance due to its high signal strength and lower interference.             Figure 12 shows the spatial throughput of Dir-BackCom Nets, where the directional antennas of BTs in different Dir-BackCom Nets have different values of antenna beamwidths and gains of the main lobe. We can observe in Figure 12 that Dir-BackCom Nets equipped with directional antennas possessing the narrowest beamwidth and the highest gain of the main lobe can obtain the highest spatial throughput when 0.01 m −2 ≤ λ BT ≤ 0.13 m −2 under α = 3 (λ CE = 0.003 m −2 ) and when 0.01 m −2 ≤ λ BT ≤ 0.15 m −2 under α = 4 (λ CE = 0.01 m −2 ). This phenomenon indicates that the spatial throughput can be improved by choosing a directional antenna with a proper antenna beamwith and gain. In addition, for each directional antenna with specific beamwidth and gain of the main lobe, there is an optimal density of BTs to maximize the spatial throughput, indicating that we can choose an optimal density to improve the spatial throughput of Dir-BackCom based on different directional antennas.

Conclusions
This article establishes a theoretic model to analyze the performance of large-scale Dir-BackCom Nets, where both BTs and BRs are equipped with directional antennas. The performance metrics include both connectivity and spatial throughput. Our theoretic model is general for BackCom Nets where BTs/BRs are equipped with directional antennas or omni-directional antennas. The accuracy of our theoretic model is verified by extensive simulations. In conclusion, this paper provides the following major findings:

•
Equipping directional antennas instead of omni-directional antenna at BTs can improve the active probability. • Employing directional antennas at either BTs or BRs can improve the connectivity and spatial throughput of BackCom Nets. • The spatial throughput can be maximized by choosing an optimal density of BTs. • Both the connectivity and spatial throughput of BackCom Nets can be improved by choosing a directional antenna with a proper antenna beamwidth and gain of the main lobe.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Deviation of the Antenna Gain of the Keyhole Antenna Model
The radiation power of the keyhole antenna model, denoted by P rad , consists of the main lobe power denoted by P m and the side lobe power denoted by P s , i.e., the radiation power of the keyhole antenna model can be given the following: where P rad = 4πU 0 , and U o is the constant radiation strength for a reference isotropic antenna in a solid angle (θ, φ). The terms P m and P s can be expressed by the following two equations, respectively.
After combining Equations (A1)-(A3), the gain of the side lobe of the keyhole antenna model, denoted by G s , can be given by the following.
After projecting the three-dimension antenna gain to two dimensions, we can have the two-dimension antenna gain of the keyhole antenna model provided in Equation (1).

Appendix B. Proof of Proposition 1
We recall two assumptions of our channel model: (1) the density of co-channel CEs is λ CE /N; (2) a BT only can choose one channel to receive RF signals, and we can see that the density of CEs that provide RF signals for a BT is λ CE /N. Thus, the density of both HPPPs Θ CE (A tm ) (CEs providing power in region A tm ) and Θ CE (A ts ) (CEs providing power in region A ts ) is λ CE /N.
Then, the CDF of the received power at a BT (P r BT ) can be derived by utilizing a Laplace transform and inverse Laplace transform. The Laplace transform of P r BT can be expressed as follows: where E[x] denotes the expectation of x, L x (s) is the Laplace transform and L x (s) = E[e −sx ]. The process (a) results from the fact that the two HPPPs of CEs distributed in regions A tm and A ts are independent with each other. As a result, variables P A tm BT and P A ts BT are independent. The Laplace transform L P A tm BT (s) can be given by the following: where W 1 = πλ CE θ tm (P CE G tm ) where v(r) is a function of r [7,39]. After employing a similar calculation with Equation (A6) and simultaneously replacing terms CE i ∈ Θ A tm CE , G tm , and θ tm by terms CE i ∈ Θ A ts CE , G ts , and θ ts , respectively, the Laplace transform L P A ts BT (s) in Equation (A5) can be given by the following: where Nα sin( 2π α ) .
Then, L P r BT (s) can be given by the following: .
Next, the CDF of P r BT can be derived by employing the inverse Laplace transform on L P r BT (s) as follows: where (d) can be obtained by applying Bromwich inversion with the modified contour [40,41].
After taking derivative of the CDF of P r BT given by Equation (A9), we can obtain the PDF of P r BT provided by Equation (9).

Appendix C. Proof of Proposition 2
A BT can be activated if the backscatter condition: P e (P h BT ) ≥ P c is satisfied. By combining Equations (3)-(6), the backscatter condition can be transferred as follows.
Let Q be equal to the right hand side of Equation (A10), i.e., Then, P a can be given by the following.
By inserting Equation (8) in Equation (A11), we can have P a provided in Equation (10).

Appendix D. Proof of Lemma 1
The Laplace transform L I (B) can be expressed as follows: where (e) results from the fact that the two HPPPs of BTs distributed in region A rm and A rs are independent with each other [34]. The terms E e −BI Arm and E e −BI Ars can be expressed as the following two equations: where ( f ) results from the fact that both point distribution processes BT j ∈ Θ BT (A m rm )\BT 0 and BT j ∈ Θ BT (A m rm ) are proven to be an identical point distribution process according to Slivnyak's Theorem [42]. Process (g) is based on the fact that the two HPPPs, Θ BT (A m rm ) and Θ BT (A s rm ), are independent. Process (h) results from the fact that the two HPPPs, Θ BT (A m rs ) and Θ BT (A s rs ), are independent. Next, we first derive Σ 1 . Then, Σ 2 , Σ 3 , and Σ 4 can be derived by using the similar calculation process with Σ 1 . Expression Σ 1 can be provided by the following: where (i) results from the fact that variable r j is independent with other variables (i.e., h j and P BT ). The process (j) is obtained by employing the PGFL property of a PPP: E[∏ r∈Θ v(r)] = exp(−λ S (1 − v(r))rdrdθ).
After inserting the expression of P r BT given by Equation (7), expression Φ 1 in Equation (A15) can be expressed as follows: where (k) results from the fact that the two HPPPs of CEs distributed in regions A tm and A ts are independent with each other. Expression φ 1 can be expressed as follows: where ξ = λ CE π(BP a ηP CE h) , and (l) is obtained by the fact that r i , h i , and h i are mutually independent. Process (m) employs the PGFL property of a PPP: E[∏ r∈Θ v(r)] = exp(−λ S (1 − v(r))rdrdθ).
The expression φ 2 can be derived by using a similar calculation process in Equation (A17) but by replacing Θ CE (A tm ) with Θ CE (A ts ), replacing G tm G tm G rm with G ts G tm G rm , and replacing θ tm wtih θ ts . Then, φ 2 can be given by the following.

Appendix E. Proof of Theorem 2
The spatial throughput can be expressed as follows: where (n) is obtained by letting e t − 1 > δ. Then, we can have t > log 2 (1 + δ). Process (o) results from P c (e t − 1) = P(SINR ≥ e t − 1, P r BT ≥ Q) according to Equation (11). Then, after inserting Equation (15) into Equation (A24), the spatial throughput can be expressed by Equation (17).