## 1 Introduction
Additionally, the atmospheric turbulence effect over THz link channels has been experimentally studied in a laboratory environment in [14; 38], as well as theoretically in [39], where in the authors estimated turbulence-induced THz scintillation through the estimation of the scintillation index metric. Their results mainly indicate that although atmospheric turbulence is less significant for THz compared to FSO channels, it should not be neglected, especially for the longer propagation distances that modern THz links can support. Thus, stochastic distribution models are needed to emulate these turbulence-induced THz signal scintillations. In this respect, the authors of [19; 40] recently introduced in the THz area the well-known (from the FSOs open literature) Gamma-Gamma, Lognormal and Exponentiated Weibull turbulence distribution models, while also even more recently the authors of [41] proposed the Gamma distribution turbulence model as a simplified alternative to the Lognormal distribution for weak-to-moderate THz turbulence conditions. Furthermore, the influence of THz beams misalignment between transmitter and receiver terminals was initially considered as a part of the shadowing effect in [42], while it was later evaluated in [43; 44] through deterministic models. Consequently, these works could not incorporate the stochastic nature of pointing errors. Their stochastic impact for THz links has been first evaluated by the authors of [27; 45; 46], relying mainly on the well-known zero boresight pointing error model which was first proposed in [24] for FSO links. In this context, the authors of [19; 40] and then the authors of [41] studied the joint impact of turbulence and zero boresight stochastic pointing errors on THz link performance for different turbulence distributions. Nevertheless, in terms of the impact of the more realistic non-zero boresight (NZB) stochastic pointing errors, there is still a research gap in the THz literature. Additionally, the use of relays has been recently introduced in the THz band in [47; 48] and then in [49; 50; 51; 52; 53; 54], which demonstrate the feasibility of establishing relay-assisted THz links to address several effects that mitigate THz performance and availability. Motivated by the above and inspired also by several multi-hop FSO relaying configurations that have been reported to successfully enhance the total propagation distance [22; 23; 55; 56; 57; 58], in this work a multi-hop THz system is investigated over different channel and link characteristics that incorporate for the first time the combined impact of all the major effects described above. The evaluation of the impact of each factor as well as their joint influence on the THz system's performance is performed by means of the crucial ABER performance metric. Therefore, the first key contribution of this paper is the establishment of a THz system and channel model including atmospheric turbulence, stochastic ZB or NZB pointing errors, attenuation due to FSPL, attenuation due to molecular absorption that arises especially from water vapor and attenuation due to weather conditions and especially due to the presence of rain or fog. It is also remarkable that the proposed system and channel model include different critical design parameters, such as the selected operational THz frequency value, the transmitter and receiver antennas' gain and the number of intermediate relay nodes along with the total distance of the multi-hop THz link. Based on the proposed model and after performing an ABER analysis, novel closed-form ABER expressions are derived. To the best of the authors' knowledge, the extracted ABER expressions are the first closed-form ABER expressions that incorporate all the above major effects and design parameters in the open technical THz literature. In this context, proper analytical results are presented for different multi-hop THz link configurations that reveal the impact of each of these performance factors along with their joint impact.