## 3 Average BER Estimation  
### On-Off Keying Modulation  
According to the On-Off keying (OOK) modulation format, the instantaneous BER for each hop can be obtained as described in [25]:  
\[P_{e,n}=p_{n}(1)p_{n}(e|1)+p_{n}(0)p_{n}(e|0), \tag{18}\]  
where \(p_{n}\) (0) and \(p_{n}\) (1) denote the probabilities of transmitting the bits \"0\" and \"1\", respectively, while \(p_{n}\) (\(e\) 10) and \(p_{n}\) (\(e\) 1) are their corresponding conditional bit-error probabilities [2]. By assuming that \(p_{n}(0)=p_{n}(1)=1/2\) and \(p_{n}(e|0)=p_{n}(e|1)\) we obtain [57]:  
\[P_{e,n}=p_{n}(e|1)=p_{n}(e|0)=\frac{1}{2}\text{erfc}\!\left(\frac{P_{l}h_{n}} {\sqrt{2}\sigma_{n}}\right), \tag{19}\]  
where \(\text{erfc}(.)\) stands for the complementary error function [70], Equation (8.250.4).  
The average BER for the \(n\)-th hop is defined as  
\[P_{e,av,n}=\int_{0}^{\infty}f_{h_{n}}(h_{n})P_{e,av,n}dh_{n}, \tag{20}\]  
Consequently, by substituting (17) and (19) by (20), we obtain  
\[\begin{array}{ll}P_{e,av,n}\!=\frac{\alpha_{n}g_{n}\psi_{n}^{2}}{2A_{0,n}g_ {n}h_{n,n}\Gamma(a_{n})\Gamma(\beta_{n})}\\ \\ \times\int_{0}^{\infty}\text{erfc}\!\left(\frac{P_{l}h_{n}}{\sqrt{2}\sigma_{n} }\right)G_{1,3}^{3,0}\!\left(\frac{a_{n}g_{n}h_{n}}{A_{0,n}g_{n}h_{n,n}}\left| \begin{array}{c}\psi_{n}^{2}\\ \psi_{n}^{2}-1,\alpha_{n}-1,\ \beta_{n}-1\end{array}\right.\right)dh_{n}\end{array} \tag{21}\]  
In order to solve the integral of (20), we initially transform the above complementary error function in terms of Meijer G-function by utilizing [71], (03.04.26.0009.01). Next, by using [72], Equation (21), and after performing some mathematical manipulations, we obtain:  
\[\begin{array}{ll}P_{e,av,n}\!=\frac{2^{\alpha_{n}+\beta_{n}-4}\psi_{n}^{2}}{ \sqrt{\pi^{3}\Gamma(a_{n})\Gamma(\beta_{n})}}\\ \\ \times G_{7,4}^{2,6}\!\left(\frac{8p_{1}^{2}A_{0,n}^{2}g_{1,n}^{2}}{\sqrt{\pi^{ 3}\Gamma(a_{n})}^{2}g_{1,n}^{2}}\left|\begin{array}{c}1,\frac{1-\psi_{n}^{2} }{2},\frac{2-\psi_{n}^{2}}{2},\frac{1-\alpha_{n}}{2},\frac{2-\alpha_{n}}{2},\frac{1-\beta_{n}}{2},\frac{2-\beta_{n}}{2}\\ 0,\frac{1}{2},-\frac{\psi_{n}^{2}}{2},\frac{1-\psi_{n}^{2}}{2}\end{array} \right.\right),\end{array} \tag{22}\]  
Next, by following the analysis performed above and correspondingly to Equation (23), we obtain the ABER for each \(L\)-PAM THz hop as  
\[\begin{array}{ll}P_{e,a\
u,n,L-PAM}&=\frac{2^{\pi n+\beta n-3}\phi_{n}^{2}}{ \sqrt{\pi}^{3}\Gamma(\alpha_{R})\Gamma(\beta_{n})}\frac{(L-1)}{(L\log_{2}L)}\\ &\times G_{6,3}^{2.5}\left(\frac{16(\log_{2}L)\Gamma_{H_{I}J}^{2}h_{I}^{2}h_{I }^{2}a_{I}^{2}a_{I}^{2}a_{I}^{2}}{N_{0,\mu}A_{\
u=2\pi}^{2}a_{I}^{2}a_{I}^{2} \theta_{n}^{2}}\right|&1,\frac{2-\psi_{n}^{2}}{2},\frac{1-\alpha_{n}}{2},\frac {2-\alpha_{n}}{2},\frac{1-\beta_{n}}{2},\frac{2-\beta_{n}}{2}\\ &0,\frac{1}{2},-\frac{\psi_{n}^{2}}{2}\end{array} \tag{29}\]  
Consequently, correspondingly to Equation (27), we eventually obtain the ABER for the total \(L\)-PAM THz multi-hop system as  
\[\begin{array}{ll}P_{e,a\
u,tot,L-PAM}=\sum\limits_{n=1}^{N}\Phi_{n}G_{6,3}^{ 2.5}\left(\frac{16(\log_{2}L)\Gamma_{H_{I}J}^{2}h_{I}^{2}a_{I}^{2}a_{I}^{2}a_ {I}^{2}}{N_{0,\mu}A_{\
u=2\pi}^{2}a_{I}^{2}a_{I}^{2}a_{I}^{2}}\right|&1,\frac{2 -\psi_{n}^{2}}{2},\frac{1-\alpha_{n}}{2},\frac{2-\alpha_{n}}{2},\frac{1-\beta _{n}}{2},\frac{2-\beta_{n}}{2}\\ &0,\frac{1}{2},-\frac{\psi_{n}^{2}}{2}\end{array} \tag{30}\]  
where we have set \(\Phi_{n}=\frac{2^{\pi n+\beta-3}\phi_{n}^{2}}{\sqrt{\pi}^{3}\Gamma(\alpha_{R}) \Gamma(\beta_{n})}\frac{(L-1)}{(L\log_{2}L)}\). It is worth noting that when \(L=2\), Equation (29) reduces into Equation (23), and Equation (30) reduces into Equation (27). The latter indicates that 2-PAM is a special case equivalent to OOK.