\section{Introduction}\label{sec:introduction}
Optical wireless communication (OWC) is a potential technology that transmits data through an unguided atmospheric channel  in the unlicensed optical spectrum \cite{Khalighi2014, Bloom2003, Kedar2004}.  However, signal transmission at a small wavelength encounters  different channel impairments such as  atmospheric turbulence,  pointing errors, and fog. The atmospheric turbulence  is caused by  the scintillation effect of light propagation whereas the  pointing errors (i.e.,  misalignment between the transmitter and receiver) happens due to the  dynamic wind loads, weak earthquakes, and thermal expansion \cite{Farid2007, Vavoulas2012}.  The impact of foggy conditions on OWC systems depends on the intensity of fog ranging between light, medium, and dense \cite{Kruse1962,Kim2001}.  Although turbulence and fog may not co-exist since both are  inversely correlated with each other \cite{Bushuev2006, Farid2007},  the effect of turbulence  can not be ignored in light foggy conditions. The combined effect of atmospheric turbulence, pointing errors, and fog has a  detrimental  effect on  the signal quality and presents a major challenge  in the  OWC deployment in outdoor environments.


The use of relaying  has  been extensively studied to improve the performance of OWC systems under the effect of turbulence and/or pointing errors \cite{Safari2008, Aghajanzadeh2011, multi_hop_turb2015, Chatzidiamantis2013,parallel_fso2015,  multi_fso2013, parallel_multi_fso2016, Karimi2011, Kazemlou2011, Bayaki2012, Kashani2012, Trinh2015,Yang2014_relay,  dual_hop_turb2017, Dabiri2018, Huang2018}. In the aforementioned and related research, the statistical  effect of foggy channels combined with pointing errors and atmospheric turbulence has not been considered. Recent studies confirm that the signal attenuation  in the fog is not deterministic but follows a probabilistic model  \cite{Khan2009, Esmail2016_Photonics,Esmail2017_Access}.  In \cite{Esmail2016_Photonics}, the authors developed Johnson SB based probability distribution function (PDF) as a model for the random fog channel. They studied numerically the system performance in terms of average bit error rate (BER) and channel capacity. Considering the intractability of the Johnson SB for performance analysis,  the authors in \cite{Esmail2017_Access} proposed Gamma distribution for the signal attenuation in foggy weather and evaluated average signal-to-noise ratio (SNR), ergodic rate, outage probability, and BER. In our recent paper \cite {Rahman2020}, we  proposed a multi-aperture OWC system to mitigate the effect of fog. The authors in \cite{Schimmel2018} presented ultrashort high-intensity
laser filaments for high-bit-rate transmissions over the fog.  These studies show that  the OWC performance is significantly limited in dense fog but can provide  acceptable performance in the light fog over shorter links.  However, combining the effect of pointing errors with fog shows high degradation in performance even in light foggy conditions\cite{Esmail2017_Photonics, rahman2020cl}. The authors in \cite{Esmail2017_Photonics} have considered three techniques to  mitigate this effect: multi-hop relay systems using decode-and-forward (DF), active laser selection, and parallel radio frequency/free space optical (RF/FSO)  link.  Laser selection and hybrid transmission techniques require feedback from the receiver to the transmitter, thereby increasing the overhead.  On the other hand,  Multi-hop relaying requires channel state information (CSI) at each relay to decode the signal, which can be hard in practice.  In \cite{Esmail2017_Photonics}, the outage probability of a multi-hop FSO system under the combined effect of fog and pointing errors has been analyzed. Analyzing the system performance using different performance metrics such as average SNR, ergodic rate, and average BER is desirable for an efficient deployment of OWC systems under the combined effect of fog and pointing errors. It is noted that there are no analyses of average SNR, ergodic rate, BER, and outage probability for the performance improvement of a relay-assisted OWC system under the statistical effect of random fog, pointing errors, and atmospheric turbulence.

 
In this paper,  we  analyze the end-to-end performance of a relay-assisted OWC system under the combined effect of  fog, pointing errors, and atmospheric turbulence by considering a single  DF relaying with  no direct transmission to the destination. The major contributions of the paper are as follows:\vspace{-0.1mm}
\begin{itemize}[leftmargin=*]
\item We consider a generalized  model for the end-to-end channel of the relay-assisted system with independent and  not identically distributed (i.ni.d.) pointing errors, random fog with Gamma distributed attenuation coefficient,  exponentiated Weibull turbulence, and asymmetrical distance between the source and destination. We assume independent and identically distributed (i.i.d.)  atmospheric turbulence and random fog since these channel impairments may effect both  links independently generated from the same probabilistic model.
 
\item	We derive a novel asymptotic PDF of the combined channel,  and  simplify the existing cumulative distribution function (CDF) of the channel with random fog and pointing errors in terms of a single incomplete gamma function and  without involving the exponential integral.
	\item We derive analytical expressions of the outage probability,  average SNR, ergodic rate, and average BER of the relay-assisted OWC system.   The derived expressions are expressed in well-known mathematical functions and consist of  a few summation terms  over the shape parameter  of the foggy channel.  We show that the dual-hop relaying can mitigate fog,  pointing errors,  and turbulence-induced fading for high-speed OWC links.
	
\item We show  that there is a significant gap in the performance using the existing visibility range based path-loss computation as compared to the probabilistic modeling of fog and  the effect of fading due to the atmospheric turbulence is  negligible in the presence of fog.

\end{itemize}

\subsection{Related Research}

Traditionally, signal attenuation due to the fog was assumed to be deterministic and quantified using a visibility range, for example, less attenuation in light fog and more in the dense fog \cite{Kruse1962,Kim2001,Kedar2003, Naboulsi2004, Awan2008, Esmail2016_ICC,Esmail2016_Photonics, Berenguer2018}. Kruse model in \cite{Kruse1962} is based on experimental data, whereas Kim \cite{Kim2001} used Mie scattering theory to predict the signal attenuation. The authors in \cite{Naboulsi2004} \cite{Awan2008} updated the earlier models considering  modified Gamma distribution for the particle size of fog. The authors in \cite{Esmail2016_ICC} developed  a power delay model for the attenuation coefficient based on extensive measurement data. On the other hand, there are  quite a few  statistical models for the atmospheric turbulence, for example, log normal \cite{Zhu2002}, Gamma-Gamma  \cite{Andrews20015},  Mal\'aga \cite{malaga2011}, and  ${\mathcal{F}}$-distribution \cite{Peppas2020}.  Recently, Barrios and Dios \cite{Barrios2012}  proposed  the exponentiated Weibull (EW) distribution model for  atmospheric turbulence. This model provides a good fit between simulation and experimental data under the moderate-to-strong turbulence for aperture averaging conditions, as well as for point-like apertures. A distinguishing  feature of the EW model is its simple closed-form	expression of the PDF \cite{Wardhan2020}. It should be noted that the model of pointing errors \cite{Farid2007} is used widely,  assuming independent identical distributed Gaussian  for the elevation and the horizontal displacement.

Relay-assisted communication is a potential technique to deal with the channel fading in wireless systems \cite{Nosratinia2004, Li2012}. Here, a single relay  or many intermediate nodes   can assist data transmission between a single source and destination. In particular, for OWC systems,  there is a vast literature on the relaying  using amplify-and-forward (AF) and DF protocols in  \cite{Safari2008, Aghajanzadeh2011, multi_hop_turb2015, Chatzidiamantis2013,parallel_fso2015,  multi_fso2013, parallel_multi_fso2016}, all-optical relaying in  \cite{Karimi2011, Kazemlou2011, Bayaki2012, Kashani2012, Trinh2015,Yang2014_relay,  dual_hop_turb2017, Dabiri2018, Huang2018}, and relaying for hybrid  RF/FSO  systems in  \cite{Lee2011, Ansari2013, Samimi2013, assym_rf_fso2015, series_hybrid_m_channel2015, dual_hop_rf_fso_turb2016, Bag2018,Zhang2020}. In the seminal work \cite{Safari2008},  multi-hop and cooperative relaying  using AF and DF protocols  have been considered for an aggregated channel model which takes into account both path-loss and turbulence-induced log-normal fading. The authors in \cite{Aghajanzadeh2011} investigated a multi-hop relaying to mitigate the effect of fading in  FSO systems over log-normal atmospheric turbulence channels.  The end-to-end performance in terms of outage probability, the average BER, and the ergodic capacity for a multi-hop relaying with AF and DF protocols under the combined effect of Gamma-Gamma turbulence and pointing errors  have been investigated in \cite{multi_hop_turb2015}. Although the complex multi-hop relaying can provide a better performance, a dual-hop relay system (with no direct link between the source and the destination) that selects a single relay opportunistically is considered in \cite{Chatzidiamantis2013}. The authors in \cite{parallel_fso2015} studied the information-theoretic performance of parallel relaying for FSO communications over Gamma-Gamma fading channels  with a single relay  but with a line-of-sight link between the source and the destination.  An optimal relay placement
scheme for serial and parallel relaying along with a diversity gain analysis has been considered  in \cite{multi_fso2013}.  The authors in \cite{parallel_multi_fso2016} have considered the inter-relay cooperation on the outage probability and diversity order performance of the DF cooperative FSO communication systems. 

  
An all-optical relaying scheme is  efficient  since signals are processed in the optical domain  without requiring  optical-to-electrical and electrical-to-optical conversions \cite{Karimi2011, Kazemlou2011, Bayaki2012, Kashani2012, Trinh2015,Yang2014_relay,  dual_hop_turb2017, Dabiri2018, Huang2018}. The references \cite{Yang2014_relay, dual_hop_turb2017} provided analytical expressions for the outage probability, average BER, and ergodic capacity over strong atmospheric turbulence  channels with misalignment-induced pointing errors by considering a single optical AF  relay with fixed and variable gain.  Hybrid RF/FSO systems, where relays act as an interface between RF and optical links, have been studied  \cite{Lee2011, Ansari2013, Samimi2013, assym_rf_fso2015, series_hybrid_m_channel2015, dual_hop_rf_fso_turb2016, Bag2018,Zhang2020}.   A dual-hop relay system over the asymmetric links has been considered for both RF and FSO environments in \cite{Lee2011} and derived exact expression for outage probability.  The BER performance and the capacity analysis of an AF-based dual-hop mixed RF–FSO is presented in \cite{Ansari2013}, where the RF links are Rayleigh distributed and FSO links are characterized by the Gamma–Gamma distributed turbulence and pointing errors.    Considering a dual-hop transmission with a single-relay and ignoring the direct transmission, the authors in  \cite{assym_rf_fso2015, series_hybrid_m_channel2015, dual_hop_rf_fso_turb2016} have analyzed the OWC performance under the  turbulence and pointing errors.    The authors in \cite{Bag2018} have considered an FSO/RF-FSO link adaptation scheme for hybrid FSO systems and analyzed  different performance metrics like outage probability, average BER, and ergodic rate. Recently, a hybrid dual-hop relaying with mmWave  and FSO scheme  is studied in \cite {Zhang2020}. 
 
\subsection{Notations and Organization}
We list the main notations  in Table I. This paper is organized as follows: In Section II, we discuss the relay assisted OWC system model. In Section III, we analyze the OWC system performance by deriving closed-form expressions for the outage probability, average SNR, ergodic rate, and average BER. In Section IV, the simulation results of the proposed system are presented. Finally, in Section V, we provide conclusions. 

\begin{table}[tp]

	\caption{The List of Main  Notations}
	\label{table:notation_parameters}
	\centering
	\begin{tabular}{|c|p{5cm}|c|}
		\hline
		 $(\cdot)_1$ & Notation for the first link \\ \hline  $(\cdot)_2$ & Notation for the second link\\ 
		 \hline
		$y$& Received signal \\ \hline	$h$ & Random channel state\\
		\hline
		$R$ & Detector responsivity\\	
		
		\hline
		$x$ & Transmit signal intensity \\ \hline 		$w$ & AWGN  \\
		\hline
		$P_t$ & Transmit power \\ \hline 
		$d_1$ & Distance: transmitter to relay \\ \hline 
		$d_2$ & Distance: relay and destination\\ \hline 
		$d$ & $d_1+d_2$ \\ \hline $k$ & Shape parameter of foggy channel\\
		\hline
		$\beta$ & Scale parameter of foggy channel \\ \hline 
		$P_{\rm out}$ & Outage probability\\ \hline
		$\gamma_0$ &  SNR without fading\\\hline
		$\gamma$  & SNR with fading \\ \hline 	
		
			$\bar{\gamma}$ & Average SNR\\
		\hline
		
				$\bar{\eta}$ & Ergodic rate \\ \hline	
				$\bar{P_e}$ & Average BER \\ \hline	
							$\Gamma(a)$ & $\int\limits_{0}^{\infty}t^{a-1} e^{-t}dt$ \\ \hline
				$\Gamma(a,t)$ & $\int_{t}^{\infty}s^{a-1}e^{-s}ds$\\ \hline 	$Q(\gamma)$ & $\frac{1}{\sqrt{2 \pi}}\int_{\gamma}^{\infty}e^{-\frac{u^2}{2}}du$\\ \hline
 $\rm{erf}(\gamma)$ &$\frac{2}{\sqrt{\pi}}\int_{0}^{\gamma}e^{-u^2}du$\\ \hline 
${}_1F_{1}(a;b;z)$ & $\sum_{k=0}^\infty \frac{\Gamma(a+k)\Gamma(b)}{\Gamma(a)\Gamma(b+k)}\frac{z^k}{k!}$ \\ \hline
		
	\end{tabular}
	
\end{table}
\begin{figure}[tp]	
	\centering
	\includegraphics[scale=0.3]{system_model_revised_ver1}
	\vspace{-2cm}
	\caption{Relay assisted OWC system.}
	\label{system_model}
	
\end{figure}
\section{System Model}\label{sec:system_model}
We consider an OWC system using intensity modulation/direct detection (IM/DD). It consists of a single-aperture transceiver system, with a negligible line-of-sight (LOS) direct link under fog, pointing errors, and atmospheric turbulence,  as shown in Fig~\ref{system_model}.
The signal $y$ at a receiver aperture is given as 
\begin{eqnarray}
	y = h_{f} h_{p} h_tRx + w,
	\label{eq:received_signal}
\end{eqnarray} 
where $x$ is the transmitted signal, $R$ represents the detector responsivity (in amperes per watt), and $w$ represents an additive white Gaussian noise (AWGN) with variance $\sigma^2_{w}$. The terms $h_{f}$, $h_{p}$, and $h_{t}$ are the random states of the foggy channel, pointing errors, and atmospheric turbulence induced fading, respectively.

We use the PDF of fog as given in \cite{Esmail2017_Access}:
\begin{equation}
f_{h_{f}}(h_{f}) = \frac{z^k}{\Gamma(k)}\left(\ln\frac{1}{h_{f}}\right)^{k-1} h_{f}^{z-1},
\label{eq:fh}
\end{equation}
where $0<h_{f}\leq 1$, $z=4.343/\beta_{} d$, $k>0$ is the shape parameter and $\beta_{}>0$ is the scale parameter. It is noted that different pairs of $k, \beta$  determine the severity of the foggy channel such as  $\{k=2.32, \beta=13.12\}$,  $\{k=5.49, \beta=12.06\}$, $\{k=6.0, \beta=23.06\}$ for light, moderate and thick foggy conditions, respectively. The PDF of pointing errors fading $h_{p}$ is given in \cite{Farid2007}:
\begin{equation}
\begin{aligned}
f_{h_p}(h_p) &= \frac{\rho^2}{A_{0}^{\rho^2}}h_{p}^{\rho^{2}-1},0 \leq h_p \leq A_0,
\end{aligned}
\label{eq:pdf_hp}
\end{equation}
where $A_0=\mbox{erf}(\upsilon)^2$ with $\upsilon=\sqrt{\pi/2}\ a/\omega_z$ and $\omega_z$ is the beam width,	and $\rho = {\frac{\omega_{z_{\rm eq}}}{2 \sigma_{s}}}$ with  $\omega_{z_{\rm eq}}$ as the equivalent beam width at the receiver and $\sigma_{s}$ as the variance of pointing errors displacement characterized by the horizontal sway and elevation \cite{Farid2007}. Finally, the random atmospheric turbulence induced fading can be modeled over EW distribution with PDF \cite{Barrios2012}:	
\begin{equation}
\begin{aligned} f_{h_t}\left(h_t\right)= &\frac{\alpha_t \beta_t}{\eta_t}\left(\frac{h_t}{\eta_t}\right)^{\beta_t-1} \exp \left[-\left(\frac{h_t}{\eta_t}\right)^{\beta_t}\right] \\ & \times\left\{1-\exp \left[-\left(\frac{h_t}{\eta_t}\right)^{\beta_t}\right]\right\}^{\alpha_t-1}, h_t \geq 0 
\end{aligned}
\label{eq:pdf_ha_ew}
\end{equation}
where $\beta_t > 0$ is the shape parameter of the scintillation index (SI), $\eta_t > 0$ is a scale parameter of the mean value of the irradiance and $\alpha_t> 0$ is an extra shape parameter that is strongly dependent on the receiver aperture size. The specific values of the parameters $\alpha$, $\beta$ and $\eta$ as well as some expressions for evaluating these parameters is given in \cite{Barrios2012}. 

We denote by $d_1$ the distance between the source and relay and $d_2$ the distance between the relay and  destination.  For the case of relayed transmission using the DF protocol, the expressions for signals received at the relay and destination when   $x$ is the transmitted signal:
\begin{eqnarray}
	y_r = h_{f1} h_{p1}h_{t1}R_1x + w_1
	\label{eq:received_signal_at_relay}
\end{eqnarray}
\begin{eqnarray}
y_d = h_{f2} h_{p2}h_{t2}R_2x + w_2
\label{eq:received_signal_at_destination}
\end{eqnarray}
where $h_{f1}$, $h_{p1}$, $h_{t1}$ and $h_{f2}$, $h_{p2}$, $h_{t2}$ are random fog, pointing errors, and atmospheric turbulence channel states between source-relay and relay-destination, respectively, each having $w_1$ and $w_2$ as AWGNs. Note that $h_{1} = h_{f1}h_{p1} h_{t1}$ is the combined channel between the source and the relay, and $h_{2} = h_{f2}h_{p2} h_{t2}$ is the combined channel between the relay and destination.

We use a general relaying scenario  $d_2> d_1$  with different values of pointing errors parameters for both  links giving an i.ni.d fading model for pointing errors.  However, we assume i.i.d. atmospheric turbulence and random fog since these channel impairments may affect both the links independently generated from the same probabilistic model. As a special case, we also consider that  the relay is situated at the midway  between the source and destination (i.e., $d_1=d_2$) and that the channel parameters are same for both links (i.e., i.i.d. condition). 

 For comparison, we also consider a channel-assisted AF such that the received signal at the destination is: 
\begin{eqnarray}
y_d = h_{f2} h_{p2}h_{t2}Gy_r + w_2
\label{eq:received_signal_at_destination_af}
\end{eqnarray}
where $G = \sqrt{\frac{P_t}{h_{1}^2P_t^2R_1^2 + \sigma^2_{w_1}}}$ is the relaying gain.
The end-to-end SNR with the AF relaying is given by
\begin{eqnarray}
		\gamma_{}=\frac{\gamma_1\gamma_2}{\gamma_1+\gamma_2+1}
		\label{eq:end_to_end_snr}
\end{eqnarray}
It should be noted that the performance analysis using the AF relaying in closed form is intractable due to the underlying distribution function of the combined fog and pointing errors, as shown in Appendix D.
\section{Performance Analysis}\label{sec:performance_analysis}
In this section, we analyze the performance of a relay-assisted system. First, we present an asymptotic expression for the PDF of SNR under the combined effect fog, pointing errors, and atmospheric turbulence, and an exact 	CDF expression under the  combined effect fog, pointing errors.  Next, we use the derived PDF’s to analyze the outage probability, average SNR, ergodic rate, and average BER performance of the OWC system. The derived expressions show the system behavior in a relay-assisted environment under the combined effect of channel impairments that can help the network operator for an efficient design of relay-assisted OWC links.

\subsection{Probability Distribution Functions}
We define $\gamma=\gamma_0 |h|^2$ as the SNR, where $h = h_fh_p h_t$, $\gamma_0=2P^2_tR^2/\sigma^2_w $,  and $P_t$ is the average optical transmitted power. The probability distribution of the SNR $\gamma$ is given in the following Proposition:
\begin{my_proposition}\label{snr_combine}
The PDF of SNR $\gamma$ under the combined effect of random fog, pointing errors, and  atmospheric turbulence for the single OWC link is given as
	\begin{eqnarray}
		&{f_\gamma}(\gamma)= \frac{C^{(1)}}{\sqrt{\gamma\gamma_0}}\left(\sqrt{\frac{\gamma}{\gamma_0}}\right)^{\phi^2-1}-\frac{C^{(2)}}{\sqrt{\gamma\gamma_0}}\left(\sqrt{\frac{\gamma}{\gamma_0}}\right)^{\phi^2-1}\times\nonumber \\&\Gamma(k,m \ln(a/\sqrt{\gamma/\gamma_0})),~~\gamma\leq a^2 \gamma_0
	\label{eq:snr_pdf_1}
	\end{eqnarray}
	where $	{f_\gamma}(\gamma)$ is the  exact PDF under the combined effect of the random fog and pointing errors with negligible turbulence with  parameters $a=A_0$, $\phi=\rho$, and $	{f_\gamma}(\gamma)$ is the asymptotic  PDF for  the combined channel effect of fog, pointing errors, and atmospheric turbulence with the  parameters $a=\frac{\eta_t A_0}{\big(M_{r^2}\big(\frac{2\alpha_t\beta_t}{\omega^{2}_{z_{eq}}}\big)\big)^{\frac{1}{\alpha_t \beta_t}}}$, $\phi=\sqrt{\alpha_t \beta_t}$, $C^{(1)}=\frac{z^k \phi^2}{2 m^k a^{\phi^2}}$, and  $C^{(2)}=\frac{C^{(1)}}{\Gamma(k)}$.
	
\end{my_proposition}
\begin{IEEEproof}
	An exact PDF for the combined effect of pointing errors and fog (i.e., with $a=A_0$ and $\phi=\rho$ in \eqref{eq:snr_pdf_1}) is already derived in \cite{Esmail2017_Photonics}. To include the effect of atmospheric turbulence, we use the asymptotic PDF of  $h_{pt}= h_p h_t$ (i.e., combined effect of atmospheric turbulence and pointing errors) as given in \cite{Boluda-Ruiz2017}:
	\begin{eqnarray}
	f_{h_{pt}}(h_{pt})= \frac{\alpha_t \beta_t M_{r^2} \Big( \frac{2 \alpha_t \beta_t}{\omega_{z_{eq}}^2} \Big)  }{(\eta A_{0}) ^ {\alpha_t \beta_t} } h_{pt}^{\alpha_t \beta_t -1}, 0\leq h_{pt}\leq a
	\label{eq:pdf_h_atm_point_asymp}
	\end{eqnarray}
	where $a= \frac{\eta_t A_0}{\big(M_{r^2}\big(\frac{2\alpha_t\beta_t}{\omega^{2}_{z_{eq}}}\big)\big)^{\frac{1}{\alpha_t \beta_t}}}$,
	and $M_{r^2}$ is  the squared Beckmann distribution given as
	
	\begin{equation}
	\begin{aligned}
	M_{r^{2}}(t)=\frac{\exp \left(\frac{\mu_{x}^{2} t}{1-2 t \sigma_{x}^{2}}+\frac{\mu_{y}^{2} t}{1-2 t \sigma_{y}^{2}}\right)}{\sqrt{\left(1-2 t \sigma_{x}^{2}\right)\left(1-2 t \sigma_{y}^{2}\right)}}
	\end{aligned}
	\label{eq:M_r2}
	\end{equation}
		Here,  $\sigma_{x}$ and $\sigma_{y}$ represent different jitters for the horizontal displacement $x$ and the elevation $y$, and $\mu_{x}$ and $\mu_{y}$ represent different boresight errors in each axis of the receiver plane, i.e.	$x \sim N(\mu_{x}$, $\sigma_{x}$) and $y \sim N(\mu_{y}$ , $\sigma_{y}$).			
The PDF of the  combined channel  $h= h_f h_{pt}$ can be expressed as  \cite{Papoulis2001}
	\begin{align}
	\begin{split}
	f_h(h)=\int_{}^{}f_{h_{pt}|h_{f}}({h_{pt}|h_{f}}) {f_h}_{f}(h_{f})dh_{f},
	\label{pe_5}
	\end{split}
	\end{align}
Using \eqref{eq:fh} and \eqref{eq:pdf_h_atm_point_asymp}  in \eqref{pe_5}, we get
	\begin{align}
	f_h(h) &= \frac{\alpha_t\beta_t z^k h^{\alpha_t\beta_t-1}}{a^{\alpha_t\beta_t}\Gamma(k)}\int_{h/a}^{1}\left[\ln \left(\frac{1}{h_{f}}\right)\right]^{k-1} h_{f}^{m-1}dh_{f}
	\label{eq_int}
	\end{align}
	where $m=z - \alpha_t\beta_t$. Substituting $t=\ln(1/h_{f})$ and using the identity [[7]-pp. 348], \eqref{eq_int} yields 	
	\begin{align}
	f_h(h) &= \frac{\alpha_t\beta_t z^k}{a^{\alpha_t\beta_t}\Gamma(k)m^k}h^{\alpha_t\beta_t-1}\big[\Gamma(k)-\Gamma(k,m\ln(a/h))\big]
	\label{eq:f(h)}
	\end{align}
Substituting $\phi = \sqrt{\alpha_t \beta_t}$ and $h=\sqrt{\frac{\gamma}{\gamma_0}}$  in \eqref{eq:f(h)},  we  get the asymptotic  PDF of  SNR for the combined effect of fog, pointing errors, and atmospheric turbulence in \eqref{eq:snr_pdf_1}.
\end{IEEEproof}

It is interesting to note that the form of PDF  is similar in both  cases facilitating a unified performance analysis with and without atmospheric turbulence. Further, we show  in Section IV using numerical analysis that the asymptotic performance with atmospheric turbulence is close to the exact results.

\begin{my_proposition}
\label{prop_CDF}
The CDF of the SNR for a single link under foggy channel with pointing errors and atmospheric turbulence is given as
\begin{eqnarray}
&{F_\gamma}(\gamma)=D^{(1)}\left(\frac{a}{\sqrt{\gamma/\gamma_0}}\right)^{-\phi^2}-D^{(2)}(k-1)!\nonumber\\&\times\sum_{i=0}^{k-1}\frac{m^i z^{-i-1}}{i!}\Gamma \left(i+1,z \ln \left(a/\sqrt{\gamma/\gamma_0}\right)\right)
\label{eq:snr_cdf_3}
\end{eqnarray}                                                                     
where $D^{(1)}=\frac{z^k}{m^k}$ and $D^{(2)}=\frac{D^{(1)}\phi^2}{\Gamma(k)}$ are constants. 
\end{my_proposition}

\begin{IEEEproof}
Substituting $u=\ln \left(\frac{a}{\sqrt{\gamma/\gamma_0}}\right)$ in the second term of \eqref{eq:snr_pdf_1}, the CDF ${F_\gamma}(\gamma)= \int_{0}^\gamma f(\gamma) d \gamma$ is given as:
\begin{eqnarray}
{F_\gamma}(\gamma)=D^{(1)}\left(\frac{a}{\sqrt{\gamma/\gamma_0}}\right)^{-\phi^2}- D^{(2)}\int\limits_{u}^{\infty}e^{-\phi^2 u}\Gamma[k,mu]du
\label{eq:snr_cdf_2}
\end{eqnarray} 
Using $\Gamma(a,t) \triangleq (a-1)! e^{-t} \sum_{m=0}^{a-1} \frac{t^m}{m!}$ in \eqref{eq:snr_cdf_2}, and applying the definition of incomplete Gamma function, we get \eqref{eq:snr_cdf_3}.
	
\end{IEEEproof}
It should be noted that the derived CDF in \eqref{eq:snr_cdf_3} is different from the CDF of \cite{Esmail2017_Photonics}: it consists of a single incomplete gamma function without the exponential integral, and can useful for performance analysis in  a closed-form using integer $k$.  

Now, we discuss the distribution functions with DF relaying. We assume equal transmit power at the source and relay (i.e., $P_t = P_r$) to get the instantaneous SNRs of signals received at the relay and receiver as $\gamma_1 = \gamma_{0}h_{1}^2$ and $\gamma_2 = \gamma_{0}h_{2}^2$, respectively. Assuming that $\gamma_1$ and $\gamma_2$ are independent for analytical tractability, the expression of end-to-end SNR for the DF relaying  is given as:
\begin{eqnarray}
\gamma_{}=\min(\gamma_1, \gamma_2)
\label{eq:half_harmonic_mean}
\end{eqnarray}
It is true that $k$ and $\beta$ parameters of the random foggy channel will be identical in both the hops. However,  channel realizations at two distinct points separated over several hundred meters might not be the same. The i.i.d. assumption on the foggy channel has been considered in \cite{Esmail2017_Photonics}.

 
The CDF and PDF of end-to-end SNR for DF relaying scheme can be given as \cite{Yang2014_relay}:
\begin{eqnarray}
F(\gamma)= F_1(\gamma)+F_2(\gamma)-F_1(\gamma)F_2(\gamma)
\label{eq:aprox_snr_cdf}
\end{eqnarray}
\begin{eqnarray}
f(\gamma)= f_1(\gamma)+f_2(\gamma)-f_1(\gamma)F_2(\gamma)-f_2(\gamma)F_1(\gamma)
\label{eq:aprox_snr_pdf}
\end{eqnarray}
where $f_1(\gamma)$, $f_2(\gamma)$ are the PDF of the first link and second link, respectively. Similarly,  $F_1(\gamma)$ and $F_2(\gamma)$ are the CDF of the first link and the second link, respectively.  

 It is seen that distribution functions in  \eqref{eq:snr_pdf_1} and \eqref{eq:snr_cdf_3} involves incomplete gamma functions with logarithmic argument requiring novel approaches to performance analysis.  The use of simple  approximation of incomplete Gamma function $\Gamma[k,m\ln u] \approx u^{-m} (m\ln u)^{k-1}$ can simplify the analysis, but the derived approximate expressions grossly overestimate/underestimate the exact performance.

\subsection{Outage Probability}	\label{subsec:outage_probability}
Outage probability is a performance measure to demonstrate the effect of the fading channel. It is defined as the probability of
failing to reach  an SNR threshold value $\gamma_{th}$, i.e., $P_{\rm out}=P(\gamma < \gamma_{th})$.
 
	\begin{my_proposition}
	\label{prop_outage_probability}
		If $k$ and $\beta$ are  the parameters of foggy channel,  $a_1$, $a_2$ and $\phi_1$, $\phi_2$ are the channel parameters (of pointing errors and atmospheric turbulence), then  an exact closed-form expression for the outage probability is given as
	\begin{eqnarray}
P_{\rm out}=P_{\rm out_1}+P_{\rm out_2}-P_{\rm out_1}P_{\rm out_2},
\label{eq:outage_probability}
\end{eqnarray}
where outage probabilities of source-to-relay ($P_{\rm out_1}$) and relay-to-source ($P_{\rm out_2}$) are 
	\begin{eqnarray}
	&P_{\rm out_1}= D_1^{(1)}\left(\frac{a_1^2\gamma_0}{\gamma_{\rm th}}\right)^{-\frac{\phi_1^2}{2}}-D_1^{(2)}(k-1)!\left(\frac{a_1^2\gamma_0}{\gamma_{\rm th}}\right)^{-\frac{z_1}{2}}\nonumber\\&\sum_{i=0}^{k-1}\sum_{j=0}^{i}\frac{m_1^i z_1^{j-i-1}\left(\ln\left(\frac{{a_1}\sqrt{\gamma_0}}{\sqrt{\gamma_{th}}}\right)\right)^{j}}{j!}
	\label{eq:outage_exact_1}
	\end{eqnarray}

\end{my_proposition}

\begin{eqnarray}
&P_{\rm out_2}= D_2^{(1)}\left(\frac{a_2^2\gamma_0}{\gamma_{\rm th}}\right)^{-\frac{\phi_2^2}{2}}-D_2^{(2)}(k-1)!\left(\frac{a_2^2\gamma_0}{\gamma_{\rm th}}\right)^{-\frac{z_2}{2}}\nonumber\\&\sum_{i=0}^{k-1}\sum_{j=0}^{i}\frac{m_2^i z_2^{j-i-1}\left(\ln\left(\frac{{a_2}\sqrt{\gamma_0}}{\sqrt{\gamma_{th}}}\right)\right)^{j}}{j!}
\label{eq:outage_exact_2}
\end{eqnarray}
where $z_1=4.343/\beta d_1$, $z_2=4.343/\beta d_2$ with $d_2\geq d_1$ as the transmission link lengths for a relay-assisted  OWC system.
\begin{IEEEproof}
Using $\gamma={\gamma_{th}}$, $m_1=z_1-\phi_1^2$, $m_2=z_2-\phi_2^2$ and series expansion of incomplete Gamma function $\Gamma(a,t) \triangleq (a-1)! e^{-t} \sum_{m=0}^{a-1} \frac{t^m}{m!}$\cite{Zwillinger2014} in \eqref{eq:snr_cdf_3}, we get the term in \eqref{eq:outage_exact_1} and \eqref{eq:outage_exact_2}, thus  \eqref{eq:outage_probability}. 	
\end{IEEEproof}

It is well known that  the diversity order  of the FSO channel with atmospheric turbulence (Gamma-Gamma model with parameters $\alpha_{\rm gamma}, \beta_{\rm gamma}$) and  pointing errors is given as $M= \frac{1}{2}\min\{\alpha_{\rm gamma}^2, \beta_{\rm gamma}^2, \phi^2\}$   \cite{Song2013, Yang2014}. Using \eqref{eq:outage_probability} and following the similar steps of  \cite{rahman2020cl}, the diversity order can be derived as $M=  \frac{1}{2}\min\{\phi_1^2, \phi_2^2, z_1, z_2\}$. The expression shows various possibilities of mitigating the impact of pointing errors and fog. For example, under certain conditions  $M=\frac{z_1}{2}=\frac{2.171}{\beta d_1}$, and thus the performance depends only on the fog parameter $\beta$ and source-to-relay distance $d_1$.
\subsection{Average SNR and Ergodic Rate}\label{subsec:average_snr_rate}
The average SNR and ergodic rate performance is computed using the PDF of SNR, which is given in \eqref{eq:aprox_snr_pdf}.   Thus,  expressions for the average SNR and ergodic rate for IM/DD \cite{Aghajanzadeh2011}:
\begin{align}
& \bar{\gamma}=\int\limits_{0}^{\min (a_1^2\gamma_0, a_2^2\gamma_0)}  \gamma f(\gamma)d\gamma
\label{eq:avg_snr_approx}
\end{align}
\begin{eqnarray}
&\bar{\eta}=\int\limits_{0}^{\min (a_1^2\gamma_0, a_2^2\gamma_0)}  \log_2(1+\frac{e}{2\pi}\gamma) f(\gamma)d\gamma
\label{eq:avg_rate_approx}
\end{eqnarray}
It is noted that we have assumed that relay requires negligible time to relay the data  while computing the ergodic rate.  Further, considering the shorter symbol duration (in the range of few nanoseconds), the  OWC  channel can be considered   a slow-fading channel \cite{Esmail2017_Access}. We analyze the ergodic rate performance  to provide an estimate on the throughput of the system. It should be noted that there is  a vast literature  on the ergodic rate performance  on the slow fading FSO channels (See \cite{Yang2014_relay, multi_hop_turb2015, dual_hop_turb2017}, and references therein).

\begin{my_theorem}
	\label{theorem:snr_asymt}
	If $k$ and $\beta$ are the parameters of the foggy channel,  $a_1$, $a_2$, $\phi_1$, and $\phi_2$ are the channel parameters (of pointing errors and atmospheric turbulence), and $z_1=4.343/\beta d_1$, $z_2=4.343/\beta d_2$ with $d_2\geq d_1$, then a closed-form expression of average SNR and a lower bound on ergodic rate are given by
	\begin{eqnarray*}
	&\bar{\gamma}={\cal{F}}_{\gamma}(a_1,a_2,\phi_1,\phi_2,z_1,z_2,k)
\\&\bar{\eta}={\cal{F}}_{\eta}(a_1,a_2,\phi_1,\phi_2,z_1,z_2,k)
	\end{eqnarray*}
	where  ${\cal{F}}_{\gamma}(a_1,a_2,\phi_1,\phi_2,z_1,z_2,k)$ and ${\cal{F}}_{\eta}(a_1,a_2,\phi_1,\phi_2,z_1,z_2,k)$ are given in \eqref{eq:asymt_snr_general_k} and \eqref{eq:asymt_rate_general_k} (in Appendix A),  respectively.
\end{my_theorem}

\begin{IEEEproof}	
We use the upper limit of the integration as $a_2^2 \gamma_0$ since $d_2\geq d_1 \rightarrow a_2\leq a_1$.  We substitute \eqref{eq:snr_pdf_1} and \eqref{eq:snr_cdf_3} in  \eqref{eq:aprox_snr_pdf} to get $f(\gamma)$ in  terms of system parameters. Since each $f_1(\gamma)$,  $f_2(\gamma)$,  $F_1(\gamma)$, and $F_2(\gamma)$ consist of two terms resulting into twelve terms of integration in  \eqref{eq:avg_snr_approx}.   
 We use $m_1=(z_1-\phi_1^2)$, $m_2=(z_2-\phi_2^2)$ and  the series expansion $\Gamma(a,t) \triangleq (a-1)! e^{-t} \sum_{m=0}^{a-1} \frac{t^m}{m!}$, and substitute $u=\frac{a_1}{\sqrt{\gamma/\gamma_0}}$, and $v=\frac{a_2}{\sqrt{\gamma/\gamma_0}}$ to simplify the integrations into  algebraic functions. Apart from simple integrations, we also encounter the following integration terms:
\begin{align}
\begin{split}
\small
\int\limits_{1}^{\infty} u^{-n}\left(\ln(u)\right)^{p} du, \hspace{-1mm} \int\limits_{a}^{\infty} u^{-n-3}\left(\ln(u)\right)^{p} du, \hspace{-1mm} \int\limits_{1}^{\infty} u^{-n} \Gamma[k,n\ln u] du
\label{eq:identity_1_asmyt2}
\end{split}
\end{align}
We solve the integrations of  \eqref{eq:identity_1_asmyt2} in closed-forms, as represented in  Appendix B.  Using \eqref{eq:identity_1_asmyt}, \eqref{eq:identity_2_asmyt} and \eqref{eq:identity_3_asmyt} of Appendix B in \eqref{eq:avg_snr_approx} with some algebraic simplifications, we get \eqref{eq:asymt_snr_general_k} of Theorem \ref{theorem:snr_asymt}. Similarly, to obtain the ergodic rate, we use the inequality $\log_2(1+\gamma)\geq\log_2(\gamma)$  in \eqref{eq:avg_rate_approx} and follow the same procedure used in deriving the average SNR expression. Here, in addition to terms in \eqref{eq:identity_1_asmyt2}, we also need to integrate the following:
\begin{align}
\begin{split}
\int\limits_{1}^{\infty} u^{-n}\left(\ln(au)\right)^{p} \left(\ln(u)\right)^{t} du
\label{eq:identity_4_asmyt2}
\end{split}
\end{align}
A closed-form expression of \eqref{eq:identity_4_asmyt2} is given in Appendix B. Thus, using \eqref{eq:identity_1_asmyt}, \eqref{eq:identity_2_asmyt},  \eqref{eq:identity_3_asmyt}, and  \eqref{eq:identity_4_asmyt} of Appendix B in \eqref{eq:avg_rate_approx}, we get \eqref{eq:asymt_rate_general_k} of Theorem 1.
\end{IEEEproof}	
In the following Lemma 1,  we simplify the derived expression by considering an i.i.d  model and when the  relay is located in the middle of the source and the destination.
\begin{my_lemma}
	\label{lemma:snr_symt}
	If $k$ and $\beta$ are  the parameters of the foggy channel,  $a$ and $\phi$ are the parameters of pointing errors, and a relay is at the mid-point with $z=8.686/\beta d$,  then a  closed-form expression of of average SNR and a lower bound on the ergodic rate are:
	
	\begin{eqnarray}\label{func:symt_snr}
&\bar{\gamma}={\cal{F}}_{\gamma}(a,\phi,z,k)\\&\bar{\eta}={\cal{F}}_{\eta}(a,\phi,z,k)
\label{func:symt_rate}
\end{eqnarray}	
	where ${\cal{F}}_{\gamma}(a,\phi,z,k)$ and ${\cal{F}}_{\eta}(a,\phi,z,k)$ are given in \eqref{eq:close_form_snr_symt_1} and \eqref{eq:close_form_rate_symt_1} respectively.			
	\end{my_lemma}


\begin{IEEEproof}
Using \eqref{eq:aprox_snr_pdf} in  \eqref{eq:avg_snr_approx} and \eqref{eq:avg_rate_approx}, and noting that $a_1=a_2=a$, and  $\phi_1=\phi_2=\phi$ under the symmetric i.i.d fading model,  we get:
\begin{eqnarray}
\label{eq:snr_smyt_proof}
&\bar{\gamma}= 2 \int\limits_{0}^{a^2\gamma_0}  \gamma \left[f_\gamma(\gamma)-f_\gamma(\gamma)F_\gamma(\gamma) \right]d\gamma\\&\bar{\eta}= 2 \int\limits_{0}^{a^2\gamma_0}  \log_2(1+\frac{e}{2\pi}\gamma) \left[f_\gamma(\gamma)-f_\gamma(\gamma)F_\gamma(\gamma) \right]d\gamma
\label{eq:rate_smyt_proof}
\end{eqnarray}		
Substituting $u=\frac{a}{\sqrt{{\gamma}/{\gamma_0}}}$, $m=(z_r-\phi^2)$ and using series expansion $\Gamma(a,t) \triangleq (a-1)! e^{-t} \sum_{m=0}^{a-1} \frac{t^m}{m!}$ in \eqref{eq:snr_smyt_proof} for the average SNR and \eqref{eq:rate_smyt_proof} with the inequality $\log_2(1+\gamma)\geq\log_2(\gamma)$ for the ergodic rate, we encounter some simple integration terms along with the first and third integration terms of \eqref{eq:identity_1_asmyt2}.  Using \eqref{eq:identity_1_asmyt} and \eqref{eq:identity_3_asmyt} of Appendix B  with some algebraic simplifications, we get \eqref{func:symt_snr} and  \eqref{func:symt_rate} of Lemma \ref{lemma:snr_symt}. 
\end{IEEEproof} 
\begin{figure*}[tp]
	\begin{eqnarray}
	&{\cal{F}}_{\gamma}(a,\phi,z,k)= 2\bigg(\frac{2C^{(1)}a^{2+\phi^2}\gamma_0}{(2+\phi^2)}+2C^{(2)}a^{2+\phi^2}\gamma_0\frac{\left(-1+\left(\frac{z-\phi^2}{2+z}\right)^k\right)\Gamma[k]}{(2+\phi^2)}-\frac{C^{(1)}D^{(1)}a^{2+\phi^2}\gamma_0}{(1+\phi^2)}+C^{(1)}D^{(2)}(k-1)!\nonumber\\&\times\sum_{i=0}^{k-1}\frac{m^iz^{-i-1}}{i!}\frac{\Gamma (i+1) \left(\phi^2+z+2\right)^{-i-1} \left(\left(\phi^2+z+2\right)^{i+1}-z^{i+1}\right)}{\phi^2+2}-C^{(2)}D^{(1)}a^{2+\phi^2}\gamma_0\frac{\left(-1+\left(\frac{m}{2+m+2\phi^2}\right)^k\right)\Gamma[k]}{(1+\phi^2)}\nonumber\\&-C^{(2)}D^{(2)}a^{\phi^2+2}i!((k-1)!)^2\sum_{i=0}^{k-1}\sum_{j=0}^{k-1}\sum_{t=0}^{i}\frac{m^{i+j}z^{t-i-1}\Gamma (j+t+1)}{i!j!t!\left(\phi^2+2\right)^{j+t+1}} \bigg)
	\label{eq:close_form_snr_symt_1}
	\end{eqnarray}
\end{figure*}

\begin{figure*}[tp]
	\begin{eqnarray}
	&{\cal{F}}_{\eta}(a,\phi,z,k)= \log_2(\frac{e}{2\pi })+2.8854\bigg(\frac{2C^{(1)}a^{\phi^2}\left(-2+\phi^2\ln (a^2\gamma_0)\right)}{\phi^4}-2C^{(2)}a^{\phi^2}(k-1)!\sum_{i=0}^{k-1}\frac{m^i\left(-2(1+i)+z\ln (a^2\gamma_0)\right)}{z^{2+i}}-\nonumber\\&\frac{C^{(1)}D^{(1)}a^{\phi^2}\left(-1+\phi^2\ln (a^2\gamma_0)\right)}{\phi^4}+2C^{(1)}D^{(2)}a^{\phi^2}(k-1)!\sum_{i=0}^{k-1}\sum_{j=0}^{i}\frac{m^iz^{j-i-1}\left(-2(1+j)+(\phi^2+z)\ln (a^2\gamma_0)\right)}{(\phi^2+z)^{2+j}}+\nonumber\\&2C^{(2)}D^{(1)}a^{\phi^2}(k-1)!\sum_{i=0}^{k-1}\frac{m^i\left(-2(1+i)+(\phi^2+z)\ln (a^2\gamma_0)\right)}{(\phi^2+z)^{2+i}}+\nonumber\\&C^{(2)}D^{(2)}a^{\phi^2}((k-1)!)^2\sum_{i=0}^{k-1}\sum_{j=0}^{k-1}\sum_{t=0}^{i}\frac{m^{i+j}z^{-i-j-3}(j+t)!\left(1+j+t-z\ln (a^2\gamma_0)\right)}{j!t!2^{j+t}}\bigg)	
	\label{eq:close_form_rate_symt_1}
	\end{eqnarray}
	\hrulefill
\end{figure*}


Further, we consider light foggy conditions  (i.e., $k=2$) to derive  simple analytical expressions on the average SNR and ergodic rate for the i.i.d model in the following Corollary.

\begin{my_corollary}
	\label{corollary_snr_rate_exact_k2}
	If $k=2$ and $\beta$ are  the parameters of the foggy channel,  $a$ and $\phi$ are the parameters of pointing errors, and relay is at the mid-point with $z=8.686/\beta d$,  then a  closed-form expression of average SNR and a lower bound on the ergodic rate are given as

	\begin{eqnarray}
		&\bar{\gamma}= 2 (a \phi z)^2\gamma_0\Big(\frac{1}{(2+\phi^2)(2+z)^2}-\nonumber\\&\frac{(2 (1 + z)^3 + \phi^4 (1 + 2 z) + \phi^2 (3 + 4 z (2 + z)))}{4 (1 + \phi^2) (1 + z)^3 (2 + \phi^2 + z)^2}\Big)
	\label{eq:snr_k2}
	\end{eqnarray}
	\begin{eqnarray}
	&\bar{\eta} \geq  \log_2(\frac{e}{2\pi })+ 2\Big(\big(\frac{\phi^2 z (2 \ln a + \ln \gamma_0)-2 (2 \phi^2 + z)}{\phi^2 z \ln 2}\big)\nonumber \\ &-0.36\big(\phi^2(\phi^2 + z)^{-2}-2\phi^{-2} - 5z^{-1}  - \nonumber \\&3(\phi^2 + z)^{-1}+ 4 \ln a + 2 \ln \gamma_0\big)\Big)
	\label{eq:rate_k2}
	\end{eqnarray}
\end{my_corollary}
\begin{IEEEproof}
The proof follows the similar procedure used in Lemma \ref{lemma:snr_symt} with $k=2$.
\end{IEEEproof}

In our earlier work \cite{rahman2020cl}, we have shown that the average SNR without relaying is  $\bar{\gamma}_{\rm direct}= \frac{(z a \phi)^2 \gamma_0}{(2+\phi^2)(2+z_{\rm direct})^2}$, where $z_{\rm direct}=4.343/\beta d $. Thus,  the first term in \eqref{eq:snr_k2}  corresponds to twice of  the average SNR  without relaying.  Since the negative term in \eqref{eq:snr_k2} is insignificant, we expect a higher average SNR with relaying. Similar conclusions can be made for the ergodic rate performance. These have been extensively studied through numerical analysis in Section IV.


Finally, we demonstrate the impact of randomness in the path loss due to the fog by considering that the fog causes a deterministic path gain $L = e^{-\psi d_{}}$, where $d_{}$ is the link distance (in \mbox{km}) and $\psi$ is the atmospheric attenuation factor, which depends  on the visibility range and may depend on the wavelength \cite{Kim2001}. 
\begin{my_corollary}
	\label{proposition_snr_rate_asymptotic}
	If the fog  causes a deterministic path loss  $L_r= e^{-\psi d/2} $ with relay-assisted transmission  $d_1=d_2= d/2$ for an i.i.d fading model  	under random pointing errors and atmospheric turbulence, then  asymptotic expressions of average SNR and ergodic rate are given as
		\begin{eqnarray}
	&	\bar{\gamma}= \frac{(aL_r)^2\phi^4\gamma_0}{(2+3\phi^2+\phi^4)},\nonumber\\&~\bar{\eta}=  \log_2 (\frac{e}{2\pi })+\frac{1.4427(\phi^2\ln((aL_r)^2\gamma_0)-3)}{\phi^2}
	\label{eq:snr_rate_asymp}
	\end{eqnarray}
\end{my_corollary}
\begin{IEEEproof}
Substituting $h_{pt}=\sqrt{\frac{\gamma}{\gamma_0}}$ in \eqref{eq:pdf_h_atm_point_asymp}, we get an asymptotic  PDF of the SNR for an OWC system under the combined effect of  atmospheric turbulence and pointing errors with atmospheric path gain:
	\begin{eqnarray}
	f_{\gamma}(\gamma)= \frac{\phi^2}{2\sqrt{ \gamma \gamma_{0}} (aL)^ {\phi^2} } \left(\sqrt{\frac{\gamma}{\gamma_{0}}}\right) ^{\phi^2-1},  0\leq\gamma\leq (aL)^2\gamma_0
	\label{eq:pdf_SNR_atm_point_asymp}
	\end{eqnarray}
Using \eqref{eq:pdf_SNR_atm_point_asymp} and the CDF $F_\gamma(\gamma)=\int_{0}^{\gamma}f(\gamma)d\gamma$  in \eqref{eq:snr_smyt_proof} and \eqref{eq:rate_smyt_proof}, it is straightforward to prove \eqref{eq:snr_rate_asymp}.
\end{IEEEproof}

Comparing 	\eqref{eq:close_form_snr_symt_1} and 	\eqref{eq:close_form_rate_symt_1} with the  expressions in 	\eqref{eq:snr_rate_asymp} shows that the randomness in fog significantly complicates  the system analysis. Further, the attenuation coefficient using  deterministic path loss may overestimate/underestimate the performance obtained  with the random fog distribution.

 To compare with the relayed transmissions, we also find the average SNR and ergodic for the direct link $d$ with $L= e^{-\psi d}$:
\begin{eqnarray}
&\bar{\gamma} = \frac{\phi^2 (aL)^2}{2 + \phi^2} \gamma_{o},\nonumber\\&~~
\bar{\eta} =\log_2\left(\frac{e}{2\pi }\right)+ \frac{2}{ \log(4) \phi^2} \Big( -2 + \phi^2 \log((aL)^2 \gamma_{o}) \Big)	
\end{eqnarray}

	\begin{figure*}[tp]
	
		\begin{eqnarray}
		&{\cal{F}}_{P_e}(a_1,a_2,\phi_1,\phi_2,z_1,z_2)=\frac{A D_1^{(1)}P_1^{-\frac{\phi_1^2}{2}}}{2\sqrt{\pi}}\left(\Gamma\left(\frac{1}{2}(1+\phi_1^2)\right)-\Gamma\left(\frac{1}{2}(1+\phi_1^2),P_2\right)\right)-\frac{AD_1^{(2)}P_1^{-\frac{z_1}{2}}}{2\sqrt{\pi}z_1}\nonumber\\& \times \Bigg[\left(\Gamma\left(\frac{1}{2}(1+z_1)\right)-\Gamma\left(\frac{1}{2}(1+z_1),P_2\right)\right)\left(1+\frac{m_1}{z_1}\right)+\frac{m_1}{2}\bigg(\Gamma\left(\frac{1}{2}(1+z_1),\frac{a_2^2P_1}{a_1^2}\right)\ln\frac{a_2^2}{a_1^2}+\nonumber\\&G_{3,0}^{2,3} \left(\begin{matrix} 1,1 \\0,0,\frac{1}{2}(1+z_1)\end{matrix} \bigg| \frac{a_2^2 P_1}{a_1^2}\right)+\Gamma\left(\frac{1}{2}(1+z_1)\right)\left(\ln P_1-\psi ^{(0)}\left(\frac{1+z_1}{2}\right)\right)\bigg)\Bigg]+\frac{A D_2^{(1)}P_2^{-\frac{\phi_1^2}{2}}}{2\sqrt{\pi}}\nonumber\\& \times \left(\Gamma\left(\frac{1}{2}(1+\phi_2^2)\right)-\Gamma\left(\frac{1}{2}(1+\phi_2^2),P_2\right)\right)-\frac{AD_2^{(2)}P_2^{-\frac{z_2}{2}}}{2\sqrt{\pi}z_2}\Bigg[\left(\Gamma\left(\frac{1}{2}(1+z_2)\right)-\Gamma\left(\frac{1}{2}(1+z_2),P_2\right)\right)\left(1+\frac{m_2}{z_2}\right)\nonumber\\&+\frac{m_2}{2}\bigg(G_{3,0}^{2,3} \left(\begin{matrix} 1,1 \\0,0,\frac{1}{2}(1+z_2)\end{matrix} \bigg| P_2\right)+\Gamma\left(\frac{1}{2}(1+z_2)\right)\left(\ln P_2-\psi ^{(0)}\left(\frac{1+z_2}{2}\right)\right)\bigg)\Bigg]-\nonumber\\&\frac{A D_1^{(1)}D_2^{(1)}P_1^{-\frac{\phi_1^2}{2}}P_2^{-\frac{\phi_2^2}{2}}}{2\sqrt{\pi}}\left(\Gamma\left(\frac{1}{2}(1+\phi_1^2+\phi_2^2)\right)-\Gamma\left(\frac{1}{2}(1+\phi_1^2+\phi_2^2),P_2\right)\right)+\nonumber\\&\frac{AD_1^{(1)}D_2^{(2)}}{2\sqrt{\pi}z_2}\Bigg[P_1^{-\frac{\phi_1^2}{2}}P_2^{-\frac{z_2}{2}}\left(\Gamma\left(\frac{1}{2}(1+\phi_1^2+z_2)\right)-\Gamma\left(\frac{1}{2}(1+\phi_1^2+z_2),P_2\right)\right)\left(1+\frac{m_2}{z_2}\right)+\nonumber\\&\frac{m_2P_2^{\frac{1}{2}(-\phi_1^2-z_2)}}{2\left(\frac{a_1}{a_2}\right)^{\phi_1^2}}\bigg(G_{3,0}^{2,3} \left(\begin{matrix} 1,1 \\0,0,\frac{1}{2}(1+\phi_1^2+z_2)\end{matrix} \bigg| P_2\right)+\Gamma\left(\frac{1}{2}(1+\phi_1^2+z_2)\right)\left(\ln P_2-\psi ^{(0)}\left(\frac{1+\phi_1^2+z_2}{2}\right)\right)\bigg)\Bigg]\nonumber\\&+\frac{AD_1^{(2)}D_2^{(1)}}{2\sqrt{\pi}z_1}\Bigg[P_1^{-\frac{z_1}{2}}P_2^{-\frac{\phi_2^2}{2}}\left(\Gamma\left(\frac{1}{2}(1+\phi_2^2+z_1)\right)-\Gamma\left(\frac{1}{2}(1+\phi_2^2+z_1),P_2\right)\right)\left(1+\frac{m_1}{z_1}\right)+\nonumber\\&\frac{m_1P_1^{\frac{1}{2}(-\phi_2^2-z_1)}}{2\left(\frac{a_2}{a_1}\right)^{\phi_2^2}}\bigg(\Gamma\left(\frac{1}{2}(1+\phi_2^2+z_1),\frac{a_2^2P_1}{a_1^2}\right)\ln\frac{a_2^2}{a_1^2}+G_{3,0}^{2,3} \left(\begin{matrix} 1,1 \\0,0,\frac{1}{2}(1+\phi_2^2+z_1)\end{matrix} \bigg| \frac{a_2^2P_1}{a_1^2}\right)+\nonumber\\&\Gamma\left(\frac{1}{2}(1+\phi_2^2+z_1)\right)\left(\ln P_1-\psi ^{(0)}\left(\frac{1+\phi_2^2+z_1}{2}\right)\right)\bigg)\Bigg]-\nonumber\\&\frac{AD_1^{(2)}D_2^{(2)}}{2\sqrt{\pi}z_1z_2}\Big[P_1^{\frac{-z_1}{2}}P_2^{\frac{-z_2}{2}}\left(\Gamma\left(\frac{1}{2}(1+z_1+z_2)\right)-\Gamma\left(\frac{1}{2}(1+z_1+z_2),P_2\right)\right)\left(1+\frac{m_1}{z_1}+\frac{m_2}{z_2}+\frac{m_1m_2}{z_1z_2}\right)+\nonumber\\&\frac{m_2P_2^{\frac{1}{2}(-z_1-z_2)}}{2\left(\frac{a_1}{a_2}\right)^{z_1}}\left(G_{3,0}^{2,3} \left(\begin{matrix} 1,1 \\0,0,\frac{1}{2}(1+z_1+z_2)\end{matrix} \bigg|P_2\right)+\Gamma\left(\frac{1}{2}(1+z_1+z_2)\right)\left(\ln P_2-\psi ^{(0)}\left(\frac{1+z_1+z_2}{2}\right)\right)\left(1+\frac{m_1}{z_1}\right)\right)\nonumber\\&+\frac{m_1P_1^{\frac{1}{2}(-z_1-z_2)}}{2\left(\frac{a_2}{a_1}\right)^{z_2}}\bigg(\Gamma\left(\frac{1}{2}(1+z_1+z_2),\frac{a_2P_1}{a_1}\right)\ln \frac{a_2}{a_1}+G_{3,0}^{2,3} \left(\begin{matrix} 1,1 \\0,0,\frac{1}{2}(1+z_1+z_2)\end{matrix} \bigg|\frac{a_2^2P_1}{a_1^2}\right)+\nonumber\\&\Gamma\left(\frac{1}{2}(1+z_1+z_2)\right)\left(\ln P_1-\psi ^{(0)}\left(\frac{1+z_1+z_2}{2}\right)\right)\bigg)\left(1+\frac{m_2}{z_2}\right)+\frac{m_1m_2P_2^{\frac{1}{2}(-z_1-z_2)}}{4\left(\frac{a_1}{a_2}\right)^{z_1}}\nonumber\\& \times \big(2\ln \frac{a_1}{a_2}G_{3,0}^{2,3} \left(\begin{matrix} 1,1 \\0,0,\frac{1}{2}(1+z_1+z_2)\end{matrix} \bigg|P_2\right)-2G_{4,0}^{3,4} \left(\begin{matrix} 1,1,1 \\0,0,0,\frac{1}{2}(1+z_1+z_2)\end{matrix} \bigg|P_2\right)\nonumber\\&+\Gamma\big(\frac{1}{2}(1+z_1+z_2)\big)\big(2\ln \frac{a_1}{a_2} \ln P_2+\ln ^2P_2-2\ln \frac{a_1P_2}{a_2}\psi ^{(0)}\left(\frac{1+z_1+z_2}{2}\right)+\psi ^{(0)}\big(\frac{1+z_1+z_2}{2}\big)^2+\psi ^{(1)}\big(\frac{1+z_1+z_2}{2}\big)\big)\big)\Big]
		\label{eq:asymt_ber_1}\nonumber\\& .
		\end{eqnarray}
		$\text{where}~P_1=\frac{Ba_1^2\gamma_0}{2}~\text{and}~P_1=\frac{Ba_2^2\gamma_0}{2}$.
	
			\hrule
	\end{figure*}
	\subsection{Average BER}
In this subsection, we derive  the average BER for the proposed relay-assisted scheme. Assuming IM/DD with on-off keying (OOK) modulation, the BER can be obtained as \cite{Proakis2007}:
\begin{eqnarray}
P_e(\gamma)=AQ(\sqrt{B\gamma})
\end{eqnarray}
where $A=1$ and $B=\frac{1}{2}$ for OOK modulation. However, we use the CDF in \eqref{eq:aprox_snr_cdf} to derive the average BER $\bar{P_e}$ as \cite{Boulogeorgos2020}
\begin{align}
\begin{split}
&\bar{P_e}=-\int_{0}^{{\min (a_1^2\gamma_0, a_2^2\gamma_0)}}\frac{\mbox{d} P_e(\gamma)}{\mbox{d}\gamma}F(\gamma)d\gamma 
\\&=A\sqrt{\frac{B}{8\pi}}\int_{0}^{{\min (a_1^2\gamma_0, a_2^2\gamma_0)}}\gamma^{-\frac{1}{2}}e^{-\frac{B\gamma}{2}}F(\gamma)d\gamma
\label{func:asymt_ber}
\end{split}
\end{align}
A general solution to the integration in \eqref{func:asymt_ber} is intractable due to the exponential function and incomplete Gamma function with logarithmic argument raised to the power $k$. A general expression is intractable even with the application of  Meijer's G representation. 


 We derive a closed-form expression for the average BER for a particular integer value of $k$ (i.e., $k=2$, $k=3$, and so on). Specifically, we consider light foggy condition (with a shape parameter $k = 2$) and  general pointing errors  with asymptotic atmospheric turbulence  to provide exact closed-form expression on the average BER. As noted in \cite{Esmail2017_Photonics},  the OWC performance is not sufficient for high-speed data transmission under  moderate ($k>2$) foggy conditions with  pointing errors.  It is also noted that light foggy conditions represent many OWC scenarios \cite{Awan2008}, \cite{Esmail2016}. 

\begin{my_theorem}
\label{theorem:ber_asymt}
	If $k=2$ and $\beta$ are  the parameters of the foggy channel,  $a_1$, $a_2$ and $\phi_1$, $\phi_2$ are the parameters of pointing errors and atmospheric turbulence, and $z_1=4.343/\beta d_1$, $z_2=4.343/\beta d_2$ with $d_2\geq d_1$, then closed-form expression of average BER is:
	\begin{eqnarray}
	&\bar{P_e}={\cal{F}}_{P_e}(a_1,a_2,\phi_1,\phi_2,z_1,z_2)
	\label{eq:asymt_ber}
	\end{eqnarray}
where ${\cal{F}}_{P_e}(a_1,a_2,\phi_1,\phi_2,z_1,z_2)$ is given in \eqref{eq:asymt_ber_1}.
\end{my_theorem}




\begin{IEEEproof} 
First, we substitute \eqref{eq:aprox_snr_cdf} and  \eqref{eq:snr_cdf_3} in \eqref{func:asymt_ber}.  We use $m_1=(z_1-\phi_1^2)$, $m_2=(z_2-\phi_2^2)$ and  the series expansion $\Gamma(a,t) \triangleq (a-1)! e^{-t} \sum_{m=0}^{a-1} \frac{t^m}{m!}$, and substitute $u=\frac{a_1}{\sqrt{\gamma/\gamma_0}}$, and $v=\frac{a_2}{\sqrt{\gamma/\gamma_0}}$ to simplify the integrations into  algebraic functions. Apart from simple integration terms, we also encounter the following  terms:
\begin{eqnarray}
&\int_1^{\infty } u^{-n-2} e^{-\frac{p}{u^2}} du, \int_{a}^{\infty }u^{-n-2} e^{-\frac{p}{u^2}} \ln(u)  du, \nonumber \\ &\int_1^{\infty } u^{-n-2} e^{-\frac{p}{u^2}} \ln^2(u) du
\label{eq:identity_1_asmyt_ber1}
\end{eqnarray}
 We provide  closed-form expressions of \eqref{eq:identity_1_asmyt_ber1} in Appendix C.  Using \eqref{eq:identity_1_asmyt_ber2}, \eqref{eq:identity_3_asmyt_ber2}, and \eqref{eq:identity_4_asmyt_ber2} of Appendix C with some algebraic simplifications, we get the \eqref{eq:asymt_ber} of Theorem \ref{theorem:ber_asymt}.
\end{IEEEproof} 

Similarly, the  average BER expression for other values of $k$ can be obtained. It should be noted that the use of PDF to derive the average BER results into more terms compared to the CDF based approach of \eqref{func:asymt_ber}. Although the average BER in 	\eqref{eq:asymt_ber_1} consists of many terms, it is useful in analyzing the system performance using known mathematical functions in contrast to a single-term  integral representation in \eqref{func:asymt_ber}. To simplify the underlying expressions, we consider  deterministic path loss in the following Corollary.

\begin{my_corollary}
	\label{proposition_ber_asymptotic}
	If the fog  causes a deterministic path loss  $L_r= e^{-\psi d/2} $ with relay-assisted transmission  $d_1=d_2= d/2$ for an i.i.d fading model  	under random pointing errors and atmospheric turbulence, then  an asymptotic expression of average BER is given as
\small
	\begin{align}
\begin{split}
		&\bar{P_e}= \frac{2^{-1+\frac{\phi^2}{2}}A\sqrt{\pi}}{\big((aL)^2B\gamma_0\big)^{\phi^2}}\big[2^{\frac{\phi^2}{2}}\big(\Gamma\big(\frac{1}{2}+\phi^2,\frac{(aL)^2B\gamma_0}{2}\big)-\Gamma\big(\frac{1}{2}+\phi^2\big)\big)\\&+2\big((aL)^2B\gamma_0\big)^{\frac{\phi^2}{2}}\big(\Gamma\big(\frac{1+\phi^2}{2}\big)-\Gamma\big(\frac{1+\phi^2}{2},\frac{(aL)^2B\gamma_0}{2}\big)\big)\big]
	\label{eq:ber_asymp}
		\end{split}
	\end{align}
\end{my_corollary}

\begin{IEEEproof} 
	Using \eqref{eq:aprox_snr_cdf} in  \eqref{func:asymt_ber}, and noting that $a_1=a_2=a$, and  $\phi_1=\phi_2=\phi$ under the symmetric i.i.d fading model, we get:
	\small
	\begin{align}
	\begin{split}
	\bar{P_e}=A\sqrt{\frac{B}{8\pi}}\int_{0}^{{\min (a_1^2\gamma_0, a_2^2\gamma_0)}}\gamma^{-\frac{1}{2}}e^{-\frac{B\gamma}{2}}[2F_{\gamma}(\gamma)-(F_{\gamma}(\gamma))^2]d\gamma
	\label{func:symt_ber_tp}
	\end{split}
	\end{align}
Using \eqref{eq:pdf_SNR_atm_point_asymp} and substituting $F_\gamma(\gamma)=\int_{0}^{\gamma}f(\gamma)d\gamma$  in \eqref{func:symt_ber_tp} and applying the definition of Gamma function and incomplete Gamma function, we get \eqref{eq:ber_asymp}.
\end{IEEEproof}  
Similar to the average SNR and ergodic capacity, it can be seen  from	\eqref{eq:ber_asymp} that the consideration of  deterministic path loss simplifies the  BER expression compared with the random path loss due  to the fog.


\begin{table}[h]	
	\renewcommand{\arraystretch}{01}
	\caption{Simulation Parameters}
	\label{table:simulation_parameters}
	\centering
	\begin{tabular}{|c|c|p{2.5cm}|}
		\hline 	
		Transmitted power &$P_t$ & $0$ to $40$ \mbox{dBm} \\ \hline
		Responsitivity &	$R$ & $0.5$ \mbox{A/W}\\ \hline
		
		AWGN variance &$\sigma_w^2$ & $10^{-14}~\rm {A^2/GHz}$\footnotemark{} \\ \hline	
		Link distance &$d$ & $500$ \mbox{m} and $800$ \mbox{m}\\ \hline	
		Shape parameter of fog &$k$ & \{2.32, 5.49, 6.00\} \\ \hline	
		Scale parameter of fog &$\beta$ & \{13.12, 12.06, 23.00\}\\ \hline Aperture diameter & $D=2a_r$ & $10$ \mbox{cm} \\ \hline Normalized beam-width &$w_z/a_r$ & \{15, 20, 25\}\\ \hline Normalized jitter &$\sigma_{s}/a_r$ & \{3, 5\} \\ \hline Refractive index&$C_n^2$ & $8 \times10^{-14}$\\ \hline Wavelength &$\lambda$ & $1550$ \mbox{nm} \\ \hline  Turbulence parameters &$\{\alpha, \beta, \eta\}$ & \{3.02, 2.80, 0.84\}\\ \hline
	\end{tabular}	 
\end{table} 
\footnotetext[1]{In many papers, it is wrongly written as $\rm A^2/Hz$.}

\section{Simulation and Numerical Analysis}
In this section, we use numerical analysis and Monte Carlo (MC) simulation (averaged over $10^7$ channel realizations) to demonstrate the performance of the relay-assisted OWC system  under the combined effect of fog with pointing errors. We also provide simulation results using the exact EW turbulence model  and numerical results using the derived asymptotic analysis  when the atmospheric turbulence is also considered. We consider three simulation scenarios: fog, pointing errors, and atmospheric turbulence (FPT), fog and pointing errors (FP), deterministic path loss with pointing errors, and atmospheric turbulence (PT).   Although a direct link between the source and destination may not exist, we compare the performance of direct and relay-assisted transmissions at various link distances  and pointing errors parameters. We also analyze the channel-assisted AF relaying  by  developing an integral-form expression for the PDF of the end-to-end SNR, as presented in Appendix D.  We  validate our derived analytical expressions with numerical and simulation results.  We use standard simulation parameters of the OWC system  as given in Table \ref{table:simulation_parameters}. 




\begin{figure*}[tp]
	\begin{center}
		\subfigure[Symmetrical: $d_1=d_2=400$ with $w_z/a_r=25$, $\sigma_{s}/a_r=3$ \mbox{m}.]{\includegraphics[width=\columnwidth]{snr_ltp_pf_tpf_ver1}}\quad
		\subfigure[Asymmetrical: $d_1=300$ \mbox{m}, $d_2=500$ \mbox{m}]{\includegraphics[width=\columnwidth]{SNR_revised_800m_Asym_ver2}}
		\caption{Average SNR performance of relay-assisted OWC system for a transmission link $d=800$ \mbox{m}.}
		\label{snr_d_600m_snr_d_1}
	\end{center}
\end{figure*}

First, we demonstrate the mutual effects of different channel impairments on the OWC system by plotting the average SNR versus transmit power for a link distance of $800$ \mbox{m}, as shown in     Fig.~\ref{snr_d_600m_snr_d_1}. Comparing the PT (where deterministic path loss is considered) and FPT plots in Fig.~\ref{snr_d_600m_snr_d_1}a, we can find that the fixed visibility range based path-loss computation overestimate the average SNR by $10$ \mbox{dB} with respect to the random path loss consideration. Further, the OWC system performs similarly for FP and FPT channels, as shown in  Fig.~\ref{snr_d_600m_snr_d_1}a, since the effect of fading due to the atmospheric turbulence is negligible in the presence  of fog.  Moreover, it can also be seen that there is a small difference of around $2$ \mbox{dB} in average SNR  between our analysis for FP/FPT channels (which is based on integer-valued $k=2$) with simulation results on real-valued $k=2.32$.  In Fig.~\ref{snr_d_600m_snr_d_1}b, we analyze the effect of pointing errors on the average SNR performance over  light fog under the FPT channel and moderate fog under the FP channel considering the asymmetric placement of the relay.  It can be seen from the figure that the average SNR performance decreases with an increase in normalized beam width and jitter.  It should be noted that the effect of normalized beam width $(w_z/a_r)$ on the average SNR performance is more as compared to the  standard deviation of the normalized jitter  $(\sigma_{s}/a_r)$. It can be observed that there is a small benefit in the average SNR performance with the relaying system. However, it should be emphasized that the direct link may not exist, and  relaying is necessary for the data transmission between the source and destination.

\begin{figure*}[tp]
	\begin{center}
		\subfigure[FP channel for different fog conditions at $\sigma_{s}/a_r=3$ and $w_z/a_r=25$, asymmetrical: $d_1=200$ \mbox{m}, $d_2=300$ \mbox{m}, and symmetrical: $d_1=d_2=250$ \mbox{m}.]{\includegraphics[width=\columnwidth]{rate_revised_500m_ver2}}\quad
		\subfigure[FPT channel for  different $w_z/a_r$ at  $\sigma_{s}/a_r=3$  and  asymmetrical  relay position: $d_1=300$ \mbox{m}, $d_2=500$ \mbox{m}.]{\includegraphics[width=\columnwidth]{rate_revised_800m_ver2}}
		\caption{Ergodic rate performance of relay-assisted OWC system.}
		\label{rate_perf}
	\end{center}
\end{figure*}


\begin{figure*}[tp]
	\begin{center}
		\subfigure[FP channel with $w_z/a_r=25$, $\sigma_{s}/a_r=3$, asymmetrical: $d_1=200$ \mbox{m}, $d_2=300$ \mbox{m}, and symmetrical: $d_1=d_2=250$ \mbox{m}.]{\includegraphics[width=\columnwidth]{outage_revised_500m_ver2}}\quad
		\subfigure[FPT channel with $\sigma_{s}/a_r=3$ at different $w_z/a_r$, asymmetrical: $d_1=300$ \mbox{m}, $d_2=500$ \mbox{m}, and symmetrical: $d_1=d_2=400$ \mbox{m}.]{\includegraphics[width=\columnwidth]{outage_revised_800m_ver2}}
		\caption{Outage performance of relay-assisted OWC system at a threshold SNR $\gamma_{th} = 6$ \mbox{dB}.}
		\label{out_perf}
	\end{center}
\end{figure*}


\begin{figure*}[tp]
	\begin{center}
		\subfigure[FP channel with  $\sigma_{s}/a_r=3$ at different $w_z/a_r$, asymmetrical: $d_1=200$ \mbox{m}, $d_2=300$ \mbox{m}, and symmetrical: $d_1=d_2=250$ \mbox{m}.]{\includegraphics[width=\columnwidth]{BER_revised_500m_ver2}}\quad
		\subfigure[FPT channel with  $\sigma_{s}/a_r=3$ at different $w_z/a_r$, asymmetrical: $d_1=300$ \mbox{m}, $d_2=500$ \mbox{m}, and symmetrical: $d_1=d_2=400$ \mbox{m}.]{\includegraphics[width=\columnwidth]{BER_revised_800m_ver2}}
		\caption{Average BER performance of relay-assisted OWC system.}
		\label{ber_perf}
	\end{center}
\end{figure*}


The ergodic rate performance in Fig.~\ref{rate_perf} shows  a significant benefit of the relay-assisted system under  FP and FPT channels. However, similar to the average SNR, the impact of atmospheric turbulence is found to be negligible on the ergodic rate in the presence of  fog. The relay-assisted system gives more increment  over moderate fog than the light fog as compared to the no-relay system. Moreover, the slopes  at high transmit power  show  greater improvement with the relay-assisted transmission compared with the direct transmissions.  Comparing Fig.~\ref{rate_perf}a and  Fig.~\ref{rate_perf}b, it can be seen a greater improvement in the ergodic rate performance using the asymmetric placement of the relay with a longer link length (i.e., $800$ \mbox{m}) than a shorter one (i.e., $500$ \mbox{m}). It can also be seen from  Fig.~\ref{rate_perf}b that the normalized beam width has a significant impact on the ergodic rate performance. 
 
 In Fig.~\ref{out_perf}a, we demonstrate the outage probability performance of the OWC system for three different fog conditions (i.e., light, moderate, and thick fog) at a normalized beam-width of $(w_z/a_r=25)$ and normalized standard deviation of the jitter $(\sigma_{s}/a_r=3)$.  The figure shows that the relay-assisted transmission provides a significant gain in the outage performance comparing to the direct  transmission for both symmetrical and asymmetrical relay positions.  A better outage probability can be attributed to the lower path loss in the relayed system.  Moreover, an increase in the fog density deteriorates the outage probability performance of the OWC system: transmission in moderate fog requires $15$ \mbox{dBm} more transmit power to achieve the same outage probability ($10^{-3}$ for the asymmetric condition at $25$ \mbox{dBm})  in the light fog condition.  Further, Fig.~\ref{out_perf}a shows that the symmetric relay position provides a higher gain in the outage performance comparing the asymmetric position of the relay between the source and destination. In Fig.~\ref{out_perf}b, we demonstrate the effect of pointing errors ( i.e., normalized beam width and jitter)  on the outage probability  performance of an OWC system under FPT channel for light foggy condition and atmospheric turbulence  ($C_n^2= 8\times 10^{-14}~{\rm m^{-2/3}}$).   We find a similar trend in the performance behavior of Fig.~\ref{out_perf}a,   in addition to that,  the outage probability decreases  with an increase in normalized beam width since the  fraction of  collected power $A_0$ decreases. 
 
 
 
Finally, in Fig.~\ref{ber_perf}, we demonstrate the average BER performance of the OWC system with  OOK modulation as a function of transmit power over the light foggy condition at different normalized beam-width of $(w_z/a_r=15, 25)$ and a normalized standard deviation of the jitter $(\sigma_{s}/a_r=3)$.  Fig.~\ref{ber_perf}a shows that the relay-assisted system provides a significant improvement in the average BER performance  compared with the direct transmission, even with the same fading scenarios under the FP channel for a link distance of $500$ \mbox{m}. The considered dual-hop relay assisted OWC system requires $15$  \mbox{dBm}  transmit power to achieve an  average BER ($10^{-3}$ which is considered reasonable for OWC applications) under the  FP channel with the  symmetric condition. Moreover, the symmetric relay position provides more gain than the asymmetric condition for the average BER compared with the average SNR and ergodic rate. Further, the impact of normalized beam width is more on the  average BER, as shown in Fig.~\ref{ber_perf}a. We also demonstrate the average BER performance under the FPT channel or a  $800$ \mbox{m} link system in Fig.~\ref{ber_perf}b , which shows  similar performance behavior, as observed for the FP channel. 
 

In all the above plots (Fig.~\ref{snr_d_600m_snr_d_1} to  Fig.~\ref{ber_perf}), we have also verified our derived expressions with the simulation and numerical results. It can be seen that the derived  closed-form expressions of the outage probability, average SNR, ergodic rate, and average BER have an  excellent match with the simulation results. However, there is a gap between simulation and analysis for the ergodic rate at a lower transmit power due to the use of inequality $\log_2(\gamma)\leq \log_2(1+\gamma)$. Further, we have compared the AF relaying with the DF.  It is expected that the DF relaying would  perform better than the AF since $ \gamma^{AF}=\frac{\gamma_1 \gamma_2}{\gamma_1+\gamma_2+1}\leq \min (\gamma_1, \gamma_2)= \gamma^{DF} $, but the difference is marginal under the considered channel impairments. The  performance channel-assisted AF and DF  were also shown to be similar under certain scenarios  over strong turbulence and pointing errors in \cite{Yang2014_relay}. 
 
 \begin{figure*}[tp]

	\begin{eqnarray}
	&{\cal{F}}_{\gamma}(a_1,a_2,\phi_1,\phi_2,z_1,z_2,k)=2\Bigg(\frac{C_1^{(1)} a_2^{2+\phi_1^2}\gamma_0 }{(2+\phi_1^2)}-\frac{C_1^{(2)} a_1^{2+\phi_1^2} \gamma_0(k-1)!}{(2+z_1)}\sum_{j=0}^{k-1}\frac{\left(\frac{m_1}{2+z_1}\right)^j\Gamma\left[1+j,(2+z_1)\ln\left(\frac{a_1}{a_2}\right)\right]}{j!}\nonumber\\&+\frac{C_2^{(1)}a_2^{2+\phi_2^2}\gamma_0}{(2+\phi_2^2)}-\frac{C_2^{(2)}a_2^{2+\phi_2^2}\gamma_0\Gamma[k]}{(2+\phi_2^2)}\left(1-\left(\frac{2+z_2}{m_2}\right)^{-k}\right)-\frac{C_1^{(1)} D_2^{(1)}a_2^{2+\phi_1^2} \gamma_0}{(2+\phi_1^2+\phi_2^2)}+\nonumber\\&C_1^{(1)} D_2^{(2)}a_2^{2+\phi_1^2}\gamma_0(k-1)!\sum_{i=0}^{k-1}\frac{m_2^iz_2^{-i-1}\Gamma (i+1) \left(\phi_1^2+z_2+2\right)^{-i-1} \left(\left(\phi_1^2+z_2+2\right)^{i+1}-z_2^{i+1}\right)}{i!(\phi_1^2+2)}+\nonumber\\&\frac{C_1^{(2)} D_2^{(1)}a_1^{2+\Upsilon_1^2} \left(\frac{a_1}{a_2}\right)^{\phi_2^2}\gamma_0(k-1)!}{(2+z_1+\phi_2^2)}\sum_{j=0}^{k-1}\frac{\left(\frac{m_1}{2+z_1+\phi_2^2}\right)^j\Gamma\left[1+j,(2+z_1+\phi_2^2)\ln\left(\frac{a_1}{a_2}\right)\right]}{j!}-\nonumber\\&\frac{C_1^{(2)} D_2^{(2)}a_2^{2+\phi_1^2} \gamma_0 }{\left(\frac{a_1}{a_2}\right)^{m_1}((k-1)!)^{-2}}\sum_{i=0}^{k-1}\sum_{j=0}^{k-1}\sum_{t=0}^{i}\frac{m_2^i m_1^jz_2^{t-i-1} }{j!t!\Gamma[-j]}\bigg(\Gamma[-j]\Gamma[1+j+t](2+z_1+z_2)^{-1-j-t}\nonumber\\&\times{}_1F_1\left(-j;-j-t;(2+z_1+z_2)\ln\left(\frac{a_1}{a_2}\right)\right)+\Gamma[1+t]\Gamma[-1-j-t]\left(\ln\left(\frac{a_1}{a_2}\right)\right)^{1+j+t}\nonumber\\&\times{}_1F_1\left(1+t;2+j+t;(2+z_1+z_2)\ln\left(\frac{a_1}{a_2}\right)\right)\bigg)-\frac{C_2^{(1)} D_1^{(1)}a_2^{2+\phi_2^2}\left(\frac{a_1}{a_2}\right)^{-\phi_1^2}\gamma_0}{(2+\phi_1^2+\phi_2^2)}+\nonumber\\&\frac{C_2^{(1)} D_1^{(2)}a_2^{2+\phi_2^2} \gamma_0(k-1
		)!}{(2+z_1+\phi_2^2)}\sum_{i=0}^{k-1}\frac{m_1^iz_1^{-i-1}}{i!}\left(\frac{\Gamma\left(1+i,z_1\ln \left(\frac{a_1}{a_2}\right)\right)-\left(\frac{a_1}{a_2}\right)^{2+\phi_2^2}\left(\frac{z_1}{2+\phi_2^2+z_1}\right)^{i+1}\Gamma\left(1+i,(2+\phi_2^2+z_1)\ln \left(\frac{a_1}{a_2}\right)\right)}{2+\phi_2^2}\right)\nonumber\\&+\frac{C_2^{(2)}D_1^{(1)}a_2^{2+\phi_2^2}\left(\frac{a_1}{a_2}\right)^{-\phi_1^2}\gamma_0\Gamma[k]}{(2+\phi_1^2+\phi_2^2)}\left(1-\left(\frac{2+z_2+\phi_1^2}{m_2}\right)^{-k}\right)-\nonumber\\&\frac{C_2^{(2)} D_1^{(2)}a_2^{2+\phi_1^2} \gamma_0 }{\left(\frac{a_1}{a_2}\right)^{z_1}((k-1)!)^{-2}}\sum_{i=0}^{k-1}\sum_{j=0}^{k-1}\sum_{t=0}^{i}\frac{m_1^i m_2^jz_1^{t-i-1} }{j!t!\Gamma[-t]}\bigg(\Gamma[-t]\Gamma[1+j+t](2+z_1+z_2)^{-1-j-t}\nonumber\\&\times{}_1F_1\left(-t;-j-t;(2+z_1+z_2)\ln\left(\frac{a_1}{a_2}\right)\right)+\nonumber\\&\Gamma[1+j]\Gamma[-1-j-t]\left(\ln\left(\frac{a_1}{a_2}\right)\right)^{1+j+t}{}_1F_1\left(1+j;2+j+t;(2+z_1+z_2)\ln\left(\frac{a_1}{a_2}\right)\right)\bigg)\Bigg)
	\label{eq:asymt_snr_general_k}
	\end{eqnarray}
	\hrulefill
\end{figure*}
\begin{figure*}[tp]

	\begin{eqnarray}
	&{\cal{F}}_{\eta}(a_1,a_2,\phi_1,\phi_2,z_1,z_2,k)=\log_2 (\frac{e}{2\pi })+2.8854\Bigg[\frac{C_1^{(1)}a_2^{\phi_1^2}\left(-2+\phi_1^2\ln (a_2^2\gamma_0)\right)}{\phi_1^2}+\frac{C_2^{(1)}a_2^{\phi_2^2}\left(-2+\phi_2^2\ln (a_2^2\gamma_0)\right)}{\phi_2^2 }-\nonumber\\& C_1^{(2)}a_1^{\phi_1^2}(k-1)!\sum_{i=0}^{k-1}\frac{m_1^i\left(z_1\Gamma\left(1+i,z_1\ln\left(\frac{a_1}{a_2}\right)\right)\ln (a_1^2\gamma_0)-2\Gamma\left(2+i,z_1\ln\left(\frac{a_1}{a_2}\right)\right)\right)}{i!z_1^{2+i}}-\nonumber\\&C_2^{(2)}a_2^{\phi_2^2}(k-1)!\sum_{i=0}^{k-1}\frac{m_2^i\left(-2(1+i)+z_2\ln (a_2^2\gamma_0)\right)}{z_2^{2+i}}-\frac{C_1^{(1)}D_2^{(1)}a_2^{\phi_1^2}(-2+(\phi_1^2+\phi_2^2)\ln (a_2^2\gamma_0))}{(\phi_1^2+\phi_2^2)^2}+\nonumber\\&C_1^{(1)}D_2^{(2)}a_2^{\phi_1^2}(k-1)!\sum_{i=0}^{k-1}\sum_{j=0}^{i}\frac{m_2^iz_2^{j-i-1}(-2(1+j)+(\phi_1^2+z_2)\ln (a_2^2\gamma_0))}{(\phi_1^2+z_2)^{2+j}}+\nonumber\\&C_1^{(2)}D_2^{(1)}a_1^{\phi_1^2}\left(\frac{a_1}{a_2}\right)^{\phi_2^2}(k-1)!\sum_{i=0}^{k-1}\frac{m_1^i\left(-2\Gamma\left(2+i,(\phi_2^2+z_1)\ln \left(\frac{a_1}{a_2}\right)\right)+(\phi_2^2+z_1)\Gamma\left(1+i,(\phi_2^2+z_1)\ln \left(\frac{a_1}{a_2}\right)\right)\ln (a_1^2\gamma_0)\right)}{i!(\phi_2^2+z_1)^{2+i}}\nonumber\\&- C_1^{(2)}D_2^{(2)}a_2^{\phi_1^2}\left(\frac{a_1}{a_2}\right)^{-m_1}((k-1)!)^2\sum_{k-1}^{i=0}\sum_{j=0}^{k-1}\sum_{t=0}^{i}\frac{m_1^jm_2^iz_2^{t-i-1}}{j!t!\Gamma[-j]}\nonumber\\&\times\Bigg((z_1+z_2)^{-1-j-t}\bigg(\Gamma[-j]\Gamma[1+j+t]{}_1F_1\left(-j;-j-t;(z_1+z_2)\ln\left(\frac{a_1}{a_2}\right)\right)+\nonumber\\&\Gamma[1+t]\Gamma[-1-j-t]{}_1F_1\left(1+t;2+j+t;(z_1+z_2)\ln\left(\frac{a_1}{a_2}\right)\right)\left((z_1+z_2)\ln \left(\frac{a_1}{a_2}\right)\right)^{1+j+t}\bigg)\nonumber\\&\times(2 \ln (a_2) + \ln (\gamma_0))-2\Gamma[2+t]\Gamma[-2-j-t]{}_1F_1\left(2+t;3+j+t;(z_1+z_2)\ln\left(\frac{a_1}{a_2}\right)\right)\ln ^{2+j+t}\left(\frac{a_1}{a_2}\right)\nonumber\\&-2(z_1+z_2)^{-2-j-t}\Gamma[-j]\Gamma[2+j+t]{}_1F_1\left(-j;-1-j-t;(z_1+z_2)\ln\left(\frac{a_1}{a_2}\right)\right)\Bigg)\nonumber\\&-\frac{C_2^{(1)}D_1^{(1)}a_2^{\phi_2^2}\left(\frac{a_2}{a_1}\right)^{\phi_1^2}(-2+(\phi_1^2+\phi_2^2)\ln (a_2^2\gamma_0))}{(\phi_1^2+\phi_2^2)^2}+C_2^{(1)}D_1^{(2)}a_1^{\phi_2^2}(k-1)!\nonumber\\&\times\sum_{i=0}^{k-1}\sum_{j=0}^{i}\frac{m_1^iz_1^{j-i-1}\left(-2\Gamma\left(2+j,(\phi_2^2+z_1)\ln \left(\frac{a_1}{a_2}\right)\right)+(\phi_2^2+z_1)\Gamma\left(1+j,(\phi_2^2+z_1)\ln \left(\frac{a_1}{a_2}\right)\right)\ln (a_1^2\gamma_0)\right)}{j!(\phi_2^2+z_1)^{2+j}}+\nonumber\\&C_2^{(2)}D_1^{(1)}a_2^{\phi_2^2}\left(\frac{a_1}{a_2}\right)^{-\phi_1^2}(k-1)!\sum_{i=0}^{k-1}\frac{m_2^i(-2(1+i)+(\phi_1^2+z_2)\ln (a_2^2\gamma_0))}{(\phi_1^2+z_2)^{2+i}}-\nonumber\\&C_2^{(2)}D_1^{(2)}a_2^{\phi_2^2}\left(\frac{a_1}{a_2}\right)^{-z_1}((k-1)!)^2\sum_{k-1}^{i=0}\sum_{j=0}^{k-1}\sum_{t=0}^{i}\frac{m_1^im_2^jz_1^{t-i-1}}{j!t!\Gamma[-t]}\nonumber\\&\times\Bigg((z_1+z_2)^{-1-j-t}\bigg(\Gamma[-t]\Gamma[1+j+t]{}_1F_1\left(-t;-j-t;(z_1+z_2)\ln\left(\frac{a_1}{a_2}\right)\right)+\nonumber\\&\Gamma[1+j]\Gamma[-1-j-t]{}_1F_1\left(1+j;2+j+t;(z_1+z_2)\ln\left(\frac{a_1}{a_2}\right)\right)\left((z_1+z_2)\ln \left(\frac{a_1}{a_2}\right)\right)^{1+j+t}\bigg)\times\nonumber\\&(2 \ln (a_2) + \ln (\gamma_0))-2\Gamma[2+j]\Gamma[-2-j-t]{}_1F_1\left(2+j;3+j+t;(z_1+z_2)\ln\left(\frac{a_1}{a_2}\right)\right)\ln ^{2+j+t}\left(\frac{a_1}{a_2}\right)\nonumber\\&-2(z_1+z_2)^{-2-j-t}\Gamma[-t]\Gamma[2+j+t]{}_1F_1\left(-t;-1-j-t;(z_1+z_2)\ln\left(\frac{a_1}{a_2}\right)\right)\Bigg)\Bigg]
	\label{eq:asymt_rate_general_k}
	\end{eqnarray}
	\hrulefill
\end{figure*} 
 
 
\section{Conclusions and Future Work}

In this paper, we investigated the performance of a DF  relaying based OWC system under the combined effect of  random fog,  pointing errors, and atmospheric turbulence. We have shown that the dual-hop relaying is capable of mitigating fog,  pointing errors,  and turbulence-induced fading for high-speed OWC links.
  We analyzed the distribution functions of SNR and provided a detailed analysis for the outage probability, average SNR, ergodic rate, and average BER  considering a generalized i.ni.d fading model  with asymmetrical distance for relaying between the source and destination.  We have shown that there is a significant gap in the performance using the existing visibility range based path-loss computation as compared to the  probabilistic modeling of fog. Further, the effect of fading due to  atmospheric turbulence is shown to be negligible in the presence  of fog. Moreover,  the effect of normalized beam width  on the OWC performance is more than the standard deviation of  normalized jitter.    It was also demonstrated that the relay-assisted system shows better performance than the direct-link transmissions as a benchmark,  providing a more significant benefit in the denser fog, for longer link lengths, and when the relay is located symmetrically between the  source and destination.  Numerical analysis  shows that the derived closed-form expressions excellently match with simulation results, and thus can be implemented for real-time tuning of the system parameters to  optimize the OWC performance. 
  
  The proposed work can be augmented with several research directions.   Although the asymptotic analysis using the combined effect of random  fog, pointing errors, and atmospheric turbulence is  close to the exact results, analyzing  the exact PDF of the combined channel with different atmospheric turbulence models can be useful.    We have also shown that there is a marginal difference between the performance of DF and AF relaying under the considered channel impairments. Thus, it would be interesting to investigate low complexity all-optical relaying schemes for OWC systems under the combined effect of fog, pointing errors, and atmospheric turbulence.  
   
\begin{appendices}
\section{(Closed-form expressions of Theorem 1)}
For a better presentation, we present the closed-form expressions of  average SNR and ergodic rate  as derived in Theorem \ref{theorem:snr_asymt}  in \eqref{eq:asymt_snr_general_k} and \eqref{eq:asymt_rate_general_k} respectively.

\section{(Integrals of Theorem \ref{theorem:snr_asymt})}
To solve the first and second integrals of \eqref{eq:identity_1_asmyt2}, we substitute $\ln u=t$  and apply the definition of Gamma function and incomplete gamma function to get the following closed-form expressions:
\begin{eqnarray}\label{eq:identity_1_asmyt}
&	\int\limits_{1}^{\infty} u^{-n}\left(\ln(u)\right)^{p} du=\frac{\Gamma[p+1]}{(n-1)^{p+1}}
\end{eqnarray}
\begin{eqnarray}
\int\limits_{a}^{\infty} u^{-n-3}\left(\ln(u)\right)^{p} du=\frac{\Gamma[p+1,(2+n)\ln(a)]}{(2+n)^{p+1}}
\label{eq:identity_2_asmyt}
\end{eqnarray}

Further, to solve the third integral of \eqref{eq:identity_1_asmyt2}, we substitute $n\ln u=t$ and use the well-known identity $\int_{0}^{\infty}e^{-at}\Gamma(b,t)dt=a^{-1}\Gamma(b)(1-(a+1)^{-b})$ [\cite{Zwillinger2014}--pp.657, eq. 6.451.2] to get the following:
\begin{eqnarray}
 \int\limits_{1}^{\infty} u^{-n} \Gamma[k,n\ln u] du=\frac{\left(1-n^k(2n-1)^{-k}\right)\Gamma[k]}{n-1}
\label{eq:identity_3_asmyt}
\end{eqnarray}
Finally, we solve \eqref{eq:identity_4_asmyt2} by substituting $\ln u =t$ and $\ln a+t=v$ and apply the identity $\int_{u}^{\infty} x^{v-1} (x-u)^{\mu-1} e^{-bx} dx=b^{-\frac{\mu+v}{2}}u^{\frac{\mu+v-2}{2}}\Gamma(\mu)e^{-\frac{bu}{2}}W_{\frac{v-\mu}{2},\frac{1-\mu-v}{2}}(bu)$ [\cite{Zwillinger2014}--pp.348, eq. 3.383.4] to get a closed-form expression in terms of Whittaker function, $W_{\lambda,\mu}(z)$ defined as:
\begin{align}
\begin{split}
W_{\lambda,\mu}(z)=\frac{\Gamma(-2\mu)}{\Gamma(\frac{1}{2}-\mu-\lambda)}M_{\lambda,\mu}(z)+\frac{\Gamma(2\mu)}{\Gamma(\frac{1}{2}+\mu-\lambda)}M_{\lambda,-\mu}(z)
\label{eq:whittaker_function}
\end{split}
\end{align}
where $M_{\lambda,\mu}(z)$ denotes the Whittaker $M$-function. Now,  we express the Whittaker  $M$-function in terms of confluent hypergeometric function [\cite{Zwillinger2014}--pp.1024, eq. 9.220.4]
 $M_{\lambda,\mu}=z^{\mu+\frac{1}{2}}e^{-z/2}\Phi(\mu-\lambda+\frac{1}{2},2\mu+1:z)$, $M_{\lambda,-\mu}=z^{-\mu+\frac{1}{2}}e^{-z/2}\Phi(-\mu-\lambda+\frac{1}{2},-2\mu+1:z)$, where $\Phi(\alpha,\lambda;z)={}_1F_1(\alpha;\lambda;z)$, we get the following:
\small
\begin{eqnarray}
&\int\limits_{1}^{\infty} u^{-n}\left(\ln(au)\right)^{p} \left(\ln(u)\right)^{t} du=\frac{\Gamma[1+p+t]{}_1F_1(-p;-p-t;(n-1)\ln(a))}{(n-1)^{1+p+t}}\nonumber\\&+\frac{\Gamma[-1-p-t]\Gamma[1+t]\left(\ln(a)\right)^{1+p+t}{}_1F_1(1+t;2+p+t;(n-1)\ln(a))}{\Gamma[-p]}
\label{eq:identity_4_asmyt}
\end{eqnarray}
\section{(Integrals of Theorem \ref{theorem:ber_asymt})}
To solve the first integral of \eqref{eq:identity_1_asmyt_ber1}, we use $\int_1^{\infty}f(t)dt=\int_0^{\infty}f(t)dt-\int_0^{1}f(t)dt$,   and apply the definition of Gamma function and  incomplete Gamma function. After a few simplifications,  we get the following closed-form expression:  
\small
\begin{eqnarray}
\int_1^{\infty } u^{-n-2} e^{-\frac{p}{u^2}} du=\frac{1}{2} \left(\frac{1}{p}\right)^{\frac{n+1}{2}} \left(\Gamma \left(\frac{n+1}{2}\right)-\Gamma \left(\frac{n+1}{2},p\right)\right)
\label{eq:identity_1_asmyt_ber2}
\end{eqnarray}	

For the second integral of \eqref{eq:identity_1_asmyt_ber1}, we again use $\int_1^{\infty}f(t)dt=\int_0^{\infty}f(t)dt-\int_0^{a}f(t)dt$, and  apply the integration by parts on both the terms. Using the  identity $\int_0^{\infty } t^{v-1} e^{-ut} \ln (t)  dt=u^{-v} \Gamma (v) (\psi ^{(0)}(v)-\ln (u))$ [\cite{Zwillinger2014}--pp.573, eq. 4.352.1], and the  identity of Meijer's G function [\cite{Wolfram}, eq. (07.34.21.0084.01)],  we get the following:
\begin{eqnarray}
&\int_{a}^{\infty }u^{-n-2} e^{-\frac{p}{u^2}} \ln(u)  du=\frac{1}{4} p^{\frac{1}{2} (-n-1)} \nonumber\\&\bigg(\ln \left(\frac{1}{a^2}\right) \Gamma \left(\frac{1}{2} \left(n+1\right),\frac{p}{a^2}\right)+G_{2,3}^{3,0}\left(\frac{p}{a^2}\left|
\begin{array}{c}
1,1 \\
0,0,\frac{n+1}{2} \\
\end{array}
\right.\right)\nonumber\\&+\Gamma \left(\frac{n+1}{2}\right) \left(\ln (p)-\psi ^{(0)}\left(\frac{n+1}{2}\right)\right)\bigg)
\label{eq:identity_3_asmyt_ber2}
\end{eqnarray}
Finally, to solve the third integral of \eqref{eq:identity_1_asmyt_ber1},  we apply the  integration by parts  twice, and then follow the same procedure used in deriving \eqref{eq:identity_3_asmyt_ber2} to get:
\begin{eqnarray}
&\int_1^{\infty } u^{-n-2} e^{-\frac{p}{u^2}} \ln^2(u) du=\nonumber \\&\frac{1}{8} \bigg(-2 G_{3,4}^{4,0}\left(p\left|
\begin{array}{c}
\frac{1-n}{2},\frac{1-n}{2},\frac{1-n}{2} \nonumber\\
0,\frac{1}{2} (-n-1),\frac{1}{2} (-n-1),\frac{1}{2} (-n-1) \nonumber\\
\end{array}
\right.\right)\\&-2 \left(\ln \left(\frac{1}{p}\right)+\ln (p)\right) G_{2,3}^{3,0}\left(p\left|
\begin{array}{c}
\frac{1-n}{2},\frac{1-n}{2} \\
0,\frac{1}{2} (-n-1),\frac{1}{2} (-n-1) \\
\end{array}
\right.\right) \nonumber \\&+\left(\frac{1}{p}\right)^{\frac{n+1}{2}} \bigg(\Gamma \left(\frac{n+1}{2}\right) \bigg(2 \psi ^{(0)}\left(\frac{n+1}{2}\right) \ln \left(\frac{1}{p}\right)\nonumber\\&+\psi ^{(0)}\left(\frac{n+1}{2}\right)^2+\psi ^{(1)}\left(\frac{n+1}{2}\right)+\ln ^2\left(\frac{1}{p}\right)\bigg)-\nonumber\\& \hspace{-2cm}\left(\ln \left(\frac{1}{p}\right)+\ln (p)\right)^2 \Gamma \left(\frac{n+1}{2},p\right)\bigg)\bigg) 
\label{eq:identity_4_asmyt_ber2}
\end{eqnarray}
\section{(PDF of Amplify-and-Forward Relaying)}
To derive a tractable  analytical  expression of the PDF of end-to-end SNR for the AF relaying, we ignore the $1$ (since the PDF value is small for $\gamma<1$) in the numerator of  \eqref{eq:end_to_end_snr} to get:\begin{eqnarray}
\gamma_{}<\frac{\gamma_1\gamma_2}{\gamma_1+\gamma_2}
\label{eq:end_to_end_snr2}
\end{eqnarray}

 Denoting the PDF of $\gamma_1$ and $\gamma_2$  as $f_1(\gamma)$ and $f_2(\gamma_2)$, respectively, we  evaluate the PDFs of $1/\gamma_1$ and $1/\gamma_2$ using the inverse distribution, and employ the PDF of their sum [\cite{Prob_stats} eq.5.55] to get the resulting PDF of  $\gamma_{e2e}$ as
\begin{eqnarray}
f(\gamma) = \gamma\int\limits_{t_{min}}^{t_{\rm max}} f_1\left(\frac{\gamma}{t}\right)f_2\left(\frac{\gamma}{1-t}\right)\frac{dt}{t^2(1-t)^2}
\label{eq:exact_snr_pdf}
\end{eqnarray}

Considering  the logarithmic term in \eqref{eq:snr_pdf_1} positive and applying this constraint to $f_1\left(\frac{\gamma}{t}\right)$ and $f_2\left(\frac{\gamma}{1-t}\right)$, we get  the limits of integration in \eqref{eq:exact_snr_pdf} as
$\ln\left(\frac{a^2\gamma_{0}t}{\gamma}\right) > 0 \implies t>\frac{a^2\gamma_{0}}{\gamma}= t_{\rm min}$, and 
$
\ln\left(\frac{a^2\gamma_{0}(1-t)}{\gamma}\right) > 0 \implies t<1-\frac{a^2\gamma_{0}}{\gamma}= t_{\rm max}$.  Note that using these limits,  we can get the exact PDF in \eqref{eq:exact_snr_pdf}. 

\end{appendices}
	\linespread{0.94}
\bibliographystyle{ieeetran}
\typeout{} 
