Performance Analysis of MIMO‑mQAM System with Pointing Errors and Beam Spreading in Underwater Málaga Turbulence Channel

: Both the long‑term beam spreading caused by ocean turbulence and the pointing errors induced by the jitter of transmitters and receivers degrade the performance of underwater wireless optical communication (UWOC) links. To effectively alleviate their effects, an in‑depth study was carried out over the M á laga turbulence channel with pointing errors and beam spreading in multiple‑ input and multiple‑output (MIMO) UWOC. First, we analyzed the long‑term beam spreading and the received light power for the finite receiving aperture in the presence of pointing error displace‑ ments. Based on this, the relationship between beam spreading, pointing errors, and signal power was established. Second, the approximate expressions of the average bit error rate (BER) and the communication outage probability were derived theoretically for this MIMO system using maximal‑ ratio combining (MRC) diversity. Third, the effects of the pointing errors on the coding and the diversity gains were explored for the MIMO links. Finally, using the observed ocean data from the Global Ocean Argo gridded dataset, we numerically verified the combined effects of ocean turbu‑ lence strength, beam spreading, and pointing errors on the average BER and outage probability of this system. These results also proved that adjusting the size of the receiving aperture or the order of the multiple quadrature amplitude modulation (mQAM) could effectively mitigate their effects.


Introduction
Underwater wireless optical communication (UWOC) has many advantages, such as ultrahigh bandwidth, ultralow delay, and high security, and it can be combined with underwater acoustic communication to implement high-performance underwater communication and networking.However, the absorption, scattering, and turbulence of seawater can give rise to severe attenuation or fading of the received optical signal.Therefore, these factors can seriously degrade the performance of UWOC links [1,2].
Light beam propagation in oceanic turbulence results in beam wander induced by large-scale eddies and beam spreading caused by small-scale eddies.Their combined effects are called long-term beam spreading [3].Therefore, long-term beam spreading in oceanic turbulence has a more serious impact on UWOC systems.Moreover, due to the effects of ocean currents and other external factors on the platforms equipped with optical transceivers and repeaters in the ocean, it is for pointing errors of different levels of severity to emerge between the optical transceivers.However, most researches have investigated beam wander and beam spreading based on the assumption that transceivers are perfectly aligned in the UWOC system [4,5].Hence, it is essential to explore the misalignment of UWOC links caused by long-term beam spreading.
Over the past couple of years, a great number of in-depth studies of free space optical (FSO) links operating over the Málaga turbulent channels with and without pointing errors have been conducted [6,7].Nevertheless, the spectra of the refractive-index variations caused by temperature or pressure inhomogeneities in the atmosphere are very different from the refractive-index spectra for temperature or salinity in seawater.Thus, the UWOC physical channel is very different from the physical channel of FSO communication.
Recently, numerous studies have been performed on misaligned UWOC links.Closedform expressions for the average symbol error probability (SEP) and the average channel capacity in a system employing the M-ary pulse position modulation (PPM) over Málaga oceanic turbulence with pointing errors are proposed in [8].The average BER in a system employing rectangular quadrature amplitude modulation (QAM) and aperture receiving technology for gamma-gamma strong oceanic turbulence fading with pointing errors was analyzed in [9], and it was proved that the contributions of four oceanic turbulence parameters must be considered in practical UWOC applications.The unified expression of the multi-layer UWOC system performance, considering generalized gamma (GG), exponential GG, exponentiated Weibull, and gamma-gamma oceanic turbulence channels with pointing errors, was analyzed in [10].However, we found that the authors did not include the relevant assumptions for multiple-input multiple-output (MIMO) technology.As is known, the MIMO technology in UWOC can alleviate signal fading, reduce the system bit error rate (BER), and improve communication reliability by using the spatial diversity gain [11][12][13].Previous studies have shown that multiple beams can effectively reduce turbulence effects [14,15].Pointing errors between the misaligned transceivers were analyzed in beam array systems for UWOC in [16], but long-term beam spreading was not involved.To the best of our knowledge, nobody has researched the performance of an MIMO system with pointing errors and long-term beam spreading suitable for the arbitrary intensity of ocean turbulence so far.Thus, further research is required on the combined effects of oceanic turbulence and pointing errors in the presence of long-term beam spreading in order to improve the achievable UWOC system performance.
In our previous work [17], we assumed the salinity was 35 ppt and the temperature was 20 • C in the ocean water.Under these conditions, we analyzed the effects of the numbers of receivers and transmitters on the average BER and the outage probability in a turbulence channel with pointing errors based on its analytical solution.In this study, we considered practical cases in which the diffusivities of temperature and salinity were based on the observed data for the Global Ocean Argo gridded dataset.Then, the optical power received by the finite receiving aperture in the presence of pointing error displacements and long-term beam spreading was analyzed.In order to obtain further insights, the approximate expressions of the average BER and the outage probability for the MIMO system employing multiple quadrature amplitude modulation (mQAM) technology over the Málaga turbulence with pointing errors and long-term beam spreading were derived theoretically.We also analyzed the relationships between the long-term beam spreading and the turbulence strength and between the diversity gain and pointing errors.Our results confirmed that both the BER and the outage performance in UWOC systems are drastically influenced by the severity of pointing errors and ocean turbulence strength and that increasing the system receiving aperture diameter to mitigate the effects of long-term beam spreading and ocean turbulence can improve UWOC reliability and efficiency better than adjusting the QAM order.
The rest of the paper is organized as follows.In Sections 2 and 3, we discuss the system model and channel model for the MIMO-UWOC.In Section 4, the performance of the MIMO-UWOC system in terms of the average BER and the outage probability is analyzed.In Section 5, we present numerical results and discussions.Finally, the conclusion is given in Section 6.

System Model
A UWOC-MIMO system with M laser diodes (LDs) at the transmitter side and N photodetectors at the receiver side was employed with rectangular mQAM and Gray coding, as shown in Figure 1.At the transmitter side, the input data bits are first modulated by an electrical mQAM modulator, which drives the LDs.The light wave propagates through the oceanic turbulence channel with pointing errors to the receiver side.At the receiver side, the avalanche photodiode (APD) converts an optical signal into an electrical signal.Finally, the output signal is obtained after the maximal-ratio combining (MRC) and the mQAM demodulator.For the sake of generality, we employed the following assumptions.
(1) The UWOC turbulence channel is memoryless, stationary, and ergodic, as well as independent and identically distributed.Its channel state information (CSI) is known by the transmitters and receivers; (2) Since the optical communication link length is much longer than the optical wavelength,thesystemignoresspatialcorrelations.Moreover,thenoise at received branches is independent and uncorrelated.

System Model
A UWOC-MIMO system with M laser diodes (LDs) at the transmitter side and N photodetectors at the receiver side was employed with rectangular mQAM and Gray coding, as shown in Figure 1.At the transmitter side, the input data bits are first modulated by an electrical mQAM modulator, which drives the LDs.The light wave propagates through the oceanic turbulence channel with pointing errors to the receiver side.At the receiver side, the avalanche photodiode (APD) converts an optical signal into an electrical signal.Finally, the output signal is obtained after the maximal-ratio combining (MRC) and the mQAM demodulator.For the sake of generality, we employed the following assumptions.
(1) The UWOC turbulence channel is memoryless, stationary, and ergodic, as well as independent and identically distributed.Its channel state information (CSI) is known by the transmitters and receivers; (2) Since the optical communication link length is much longer than the optical wavelength, the system ignores spatial correlations.Moreover, the noise at received branches is independent and uncorrelated.

Transmitte r signal QAM modulation
Underwat er Turbul ence Channel It was assumed that h(t) is a stationary random process for the system's existing pointing errors and optical signal fading caused by ocean turbulence over the period entailed by each QAM symbol.Thus, the output signal of the APD receiver can be obtained by [18]: where ℜ is the APD's responsivity, g is its average gain, PS is the transmitted optical power per symbol, and m0 is the modulation index indicating the number of bits per transmitted QAM symbol.s(t) is the output signal of the electrical QAM modulator, hc is the channel attenuation, hl is the geometric loss, and n(t) is the receiver noise.In the UWOC system, the effect of the optical signal scintillation on APD shot noise leads to uncertainty in the shot noise variance, and the thermal noise in the high-sensitivity APD is very low.The ocean turbulence links with pointing errors can be modeled as a slow-fading channel.It is thus possible to consider the samples for the channel fading process at a given time h = h(t = t0).As a result, the APD shot noise present can be largely taken into account and modeled as zero-mean Gaussian white noise with the variance conditioned on h; i.e., σ ℜ [18], where q denotes the electron charge, FA is the excess noise factor of the APD, and B indicates the symbol's effective noise bandwidth.Therefore, the received instantaneous signal-to-noise ratio (SNR) contributed by the sub-channel between the m-th light beam and n-th photodetector in this UWOC-MIMO system can be expressed as It was assumed that h(t) is a stationary random process for the system's existing pointing errors and optical signal fading caused by ocean turbulence over the period entailed by each QAM symbol.Thus, the output signal of the APD receiver can be obtained by [18]: where ℜ is the APD's responsivity, g is its average gain, P S is the transmitted optical power per symbol, and m 0 is the modulation index indicating the number of bits per transmitted QAM symbol.s(t) is the output signal of the electrical QAM modulator, h c is the channel attenuation, h l is the geometric loss, and n(t) is the receiver noise.In the UWOC system, the effect of the optical signal scintillation on APD shot noise leads to uncertainty in the shot noise variance, and the thermal noise in the high-sensitivity APD is very low.The ocean turbulence links with pointing errors can be modeled as a slowfading channel.It is thus possible to consider the samples for the channel fading process at a given time h = h(t = t 0 ).As a result, the APD shot noise present can be largely taken into account and modeled as zero-mean Gaussian white noise with the variance conditioned on h; i.e., σ 2 sh = 2qg 2 F A ℜh c h l hP s m 0 B [18], where q denotes the electron charge, F A is the excess noise factor of the APD, and B indicates the symbol's effective noise bandwidth.Therefore, the received instantaneous signal-to-noise ratio (SNR) contributed by the subchannel between the m-th light beam and n-th photodetector in this UWOC-MIMO system can be expressed as where γ = ℜP s m 0 h c h l /(2MNqF A B) is the average SNR and h mn represents the instantaneous fading value from the m-th light beam to n-th photodetector.The output SNR of the UWOC-MIMO system with MRC can be indicated as follows [19]:

Channel Model
The UWOC channel losses include those caused by absorption, scattering, turbulence, and pointing errors.The channel attenuation caused by absorption and scattering can be given as h c = exp[−c(λ)z] [20], where c(λ) is the extinction coefficient, including the absorption coefficient and the scattering coefficient, and z is the communication link length.The geometric loss can be expressed as h l = D 2 rx /(zθ) 2 , where D rx symbolizes the received aperture diameter and θ indicates the transmitter divergence angle [20].From this, the channel attenuation for h c and the geometric loss for h l can be considered constant and invariant during the symbol transmission duration interval.Therefore, the UWOC channel fading can be expressed as h(t) = h t (t)•h p (t), where h t (t) and h p (t) denote the attenuation caused by ocean turbulence and pointing errors, respectively.

Turbulence Model
The Málaga distribution model is suitable for describing the ocean turbulence channel with arbitrary levels of strength from weak to strong, with α being a positive parameter related to the effective number of large-scale cells in the scattering process that determines the turbulence strength.β represents the value of the fading parameter and is a natural number.The probability density function (PDF) of h t (t) can be represented as [21]: where , and β k is the binomial coefficient.K v (•) is a modified Bessel function of the second kind with the order of v. g = 2b 0 (1 − δ) denotes the average power of the energy scattered to the receiver by off-axis eddies.2b 0 indicates the average power of the total scatter components.The parameter δ is the factor expressing the amount of scattering power coupled to the lineof-sight (LOS) component, 0 ≤ δ ≤ 1.
symbol- izes the average power from the coherent contribution, whereas Ω represents the average power of the LOS term.The parameters φ A and φ B are the deterministic phases of the LOS scattering term and the coupled-to-LOS scatter term, respectively.
The parameters α and β can be expressed as follows for a plane wave [3]: where σ 2 R is the Rytov variance for a plane wave.The scintillation index can be described by the Rytov variance for a plane wave.Thus, σ 2 R = σ 2 I for a plane wave.In general, σ 2 I describes the turbulence strength; e.g., σ 2 I < 1 for weak turbulence, σ 2 I ≈ 1 for moderate turbulence, and σ 2 I > 1 for strong turbulence.σ 2 I can be calculated as follows [3,13]: where k = 2π/λ denotes the optical wave number at the wavelength of λ.Θ = 1 and 0 for plane and spherical waves, respectively.κ is the magnitude of the spatial frequency and ζ is the normalized distance variable, 0 ≤ ζ ≤ 1. Φ n (κ) represents the power spectrum in isotropic and heterogeneous ocean conditions [22]; i.e., where ε is the rate of dissipation of kinetic energy per unit mass of fluid, as given by ε = ν 3 η −4 , where ν denotes the kinematic viscosity.η is the Kolmogorov microscale, 35, χ T , χ S , and χ TS represent the dissipation rate of the mean-squared temperature, the dissipation rate of the mean-squared salinity, and the correlation between χ T and χ S , respectively.α t is the thermal expansion coefficient.ω is a unitless parameter known as the relative strength of temperature and salinity fluctuations.d r denotes the ratio of saline eddy diffusivity to thermal eddy diffusivity in unstable stratification ocean turbulence.Related formulas and calculations for the above parameters are reported in [22,23].
Substituting Equation ( 7) into Equation ( 6), σ 2 I can be calculated as follows for a plane wave [22]: In ocean turbulence, the long-term beam spreading of W 2 LT , which can be used to describe the effective received spot size, is calculated as follows [3,24]: where is the beam radius at the variable distance z from the transmitter.W 0 is the beam radius at the transmitter.Θ 0 and Λ 0 refer to a pair of non-dimensional quantities used as transmitter beam parameters.Θ 0 is called the curvature parameter and Θ 0 = 1-z/F 0 , where F 0 is the radius of curvature.Λ 0 is the Fresnel ratio at the transmitter plane and Λ 0 = 2z/kW 2 0 .The curvature parameter Θ 0 = 1 for the collimated beam.In Equation ( 9), T indicates the change in the on-axis mean irradiance at the receiver plane, which is caused by turbulence and can be expressed as [3,25]: where Λ = 2z/kW 2 is the Fresnel ratio at the receiver plane.

Pointing Error Model
Since the size of the received light spot is always small, the probability of pointing errors is very high.In the UWOC link, p = (p x , p y ) denotes the receiver plane spatial coordinates, and the radial displacement vector can be expressed as r = r x , r y T , where r x and r y denote the displacements along the horizontal direction (p x -axis) and the vertical direction (p y -axis) of the beam in the detector plane, respectively.Assuming that the displacements follow independent Gaussian distributions along the p x -axis and p y -axis, we can express r x and r y as r x ∼ N µ x , σ 2 x and r y ∼ N µ y , σ 2 y , respectively.Then, the radial displacement r = |r| = r 2 x + r 2 y follows a Rayleigh distribution with a mean of zero and identical variance; that is, Thus, the probability density function (PDF) of h p (t) can be expressed as [26,27]: where being the equivalent beam width at the receiver and σ p the pointing error displacement standard deviation at the receiver. 2is the fraction of the collected power when the radial distance is zero, erf(•) is the error function, and v = √ π/2/(W z /a).W 2 z is the beam waist at the distance z from the transmitter, and the detection aperture radius a = D rx /2.The channel fading due to geometric spreading and pointing errors can be approximated by h p ≈ A 0 exp −2r 2 /W 2 zeq , the approximate condition of which is W z /a > 6 [28].
In terms of the analysis in Section 3.1, for the ocean turbulence channel, the beam radius of W 2 z at the variable distance of z from the optical transmitter is characterized as the long-term beam spreading W 2  LT , considering the beam spreading and the beam wander; i.e., W 2 z = W 2 LT .Moreover, ξ can be rewritten in a new form as ξ = W zeq /a 2σ p /a , where W zeq /a can be expressed as a function of W z /a based on (07.20.03.0051.01) in [29]; i.e., where 1 F 1 (•) is the confluent hypergeometric function [29], and W zeq /a increases with the increase in W z /a according to Equation (12).The parameter for W z /a is defined as the normalized received beam waist coefficient, which represents the long-term beam spreading.The parameter σ p /a is the normalized pointing error deviation, which represents the jitter strength.Its value range is between 1 and 10.The larger the σ p /a is, the more serious the pointing error is [27].Without loss of generality, it can be assumed that the radial distance caused by the pointing error is located along the p x -axis.Within the unit symbol transmission interval, the average intensity at the finite-sized received aperture can be obtained with the help of the extended Huygens-Fresnel principle, as follows [3]: where <•> stands for a statistical average operation.The average power collected by the finite aperture at the receiver plane is expressed as [30]: Substituting Equation (13) into Equation ( 14), for the integral domain in the integration calculation, the circular plane with the received aperture radius of a is equivalent to the rectangle plane with a side length of √ πa.Then, the analytic formula for the average power can be derived as follows: where 2 , and the other parameters' meanings are the same as above.

MIMO Joint Fading Channel Model
According to the above description of channel fading, the distribution of h(t) is the same as the joint distribution of h t (t)•h p (t).It is assumed that the effects of ocean turbulence and pointing errors are approximately independent of each other.Therefore, considering the combined effects of the Málaga turbulence and Rayleigh pointing errors, the PDF can be expressed as [7]: where f h|h p h h p is the condition PDF for a given channel fading h p (t) caused by pointing errors.Based on (07.34.21.0085.01)and (07.34.16.0001.01) in [29], the PDF of the received instantaneous channel fading in the MIMO system can be expressed as [6,8,31]: where G a,b c,d (•) is the Meijer-G function [29], and the other parameters' meanings are the same as above.

System Average BER Analysis
Assuming that each symbol of the mQAM signal is sent with an equal probability, each QAM symbol will be included by the M I dimensions in-phase and M Q dimensions quadrature-phase symbols, respectively, with m = M I × M Q .The system instantaneous BER is the sum of the error probability in the quadrature component and the error probability in the in-phase component.The quadrature-to-in-phase decision distance ratio r 0 = d Q /d I .d Q and d I denote the decision distance in the quadrature and in-phase components, respectively.Without loss of generality, assuming that r 0 = 1, the UWOC-MIMO system instantaneous BER can be expressed as [7,32]: where m 0 is the modulation index, and m 0 = log 2 (M I × M Q ).γ is the SNR received by the MIMO system, which is the SNR at the output of the MRC combiner, as shown in Equation ( 4).Q(•) is the Gaussian Q function.C i , C j , D i , and D j are represented as follows [7,32]: , D j = (2j + 1) .
For each optical link of the MIMO system, assuming that the channel fading caused by the ocean turbulence and pointing errors is independent and identically distributed, the joint PDF of the channel fading vector h = (h 11 , h 12 , . . ., h MN ) T , and Therefore, the system average BER, which represents the mean of the system BER for the joint fading channel, can be expressed as: Inserting Equation (3) into Equation ( 18), using the approximation of Q(x), Q(x) ≈ exp −x 2 /2 /12 + exp −2x 2 /3 /4, and then employing (07.34.21.0011.01)from [29] and letting A z = A 0 h c h l , the average BER in the UWOC-MIMO system with MRC can be deduced as: In Equation ( 21), the calculation cost is large for Meijer's G-function and its power operation in a numerical calculation with an increasing number of transceivers and increasing QAM index.Considering the parameters of (α, β) that relate to the factual oceanic conditions, it will be found that α is always greater than β.Furthermore, there is a trend that α > 3 and β→1 with increasing turbulence strength.Thus ， we can obtain the corresponding closed-form solutions for the average BER using the series expansion corresponding to Meijer's G-function ([29], (07.34.06.0006.01)),as can be seen in Equations ( 22) and ( 23): As is well-known, the smaller ξ 2 is, the stronger the impact of the pointing error it implies.ξ 2 > >1 as the pointing error becomes smaller, and ξ 2 < 1 as the pointing error becomes stronger.Depending on the value range for k, we can describe various serious pointing errors using Equation ( 22) with a wide variety of oceanic turbulence strengths.It is straightforward to show that the average BER behaves asymptotically as (G c γ) −G d , where G d and G c denote the diversity order and the coding gain, respectively [19,33].Thus, the corresponding closed-form solutions for the average BER at a high SNR can be written as: The diversity order of G d and the coding gain of G c are written as: This further shows that the average BER is determined using the optical links (M, N), the QAM order, the oceanic turbulence (α, β), and the pointing errors ξ 2 .Obviously, the effect of the QAM order on the average BER can be weakened by the optical links (M, N) and the pointing errors ξ 2 , and the effect of the oceanic turbulence (α, β) on it can only positively correlate with the pointing errors ξ 2 .The diversity order is only determined by using the optical links (M, N) and the pointing errors ξ 2 .

Link Outage Probability Analysis
For the independent and identically distributed channel in the MIMO system with MRC, its outage probability can be expressed using the numerical calculation method as follows [6]: where

indicates taking the real part, and the chosen values for
A E , N E , and Q can be found in [6].E(A E , N E , Q) is the overall error term.M γ MRC (s) is the moment- generating function (MGF) of γ MRC at the output of the MRC combiner.γ th denotes the outage SNR threshold.However, due to the long calculation time required for Equation ( 26), we adopt the solution method for the BER in this section and use the series expansion corresponding to Meijer's G-function (07.34.06.0006.01)from [29].First, according to Equations ( 3) and ( 4), we can obtain the relationship between the PDF of γ mn and the PDF of h mn , and f γ (γ mn ) = f h (γ mn /γ)/γ.Then, the relationship for the cumulative distribution function (CDF) for γ MRC at the output of the MRC combiner and the CDF for γ mn can be obtained as [19]: Next, using the Laplace transform, F γ MRC (s) can be deduced as: The outage probability of the MIMO system with MRC can be obtained as [17,30,34]: where f γ MRC (γ) is the PDF of γ MRC at the output of the MRC combiner.It can be obtained by applying the Laplace inverse transform to Equation (28).Moreover, for the same reason as for Equation ( 22), the outage probability can be deduced as:

Simulation Parameters
The system parameters set in our simulation are shown in Table 1, and we employed observed ocean data, including the pressure, temperature, and salinity of seawater, from the Global Ocean Argo gridded dataset [35].In order to identify the reliability of the simulation results, we selected four Argo nodes from ocean areas with different longitudes and latitudes in our simulation, as shown in Table 2.All the results presented in Sections 5.2 and 5.3 were verified via Monte-Carlo simulations.First, we present the numerical results for the scintillation indexes, effective received spot size, and the Málaga parameters (α, β) for the four Argo nodes in Table 2.It can be seen from Figure 2 that the received spot size W LT and the scintillation index σ 2 I synchronously increased with increasing z and χ T .In Figure 3, the change trend for the parameters α and β with increases in z and χ T is shown.It can be clearly seen that the value of α was always larger than that of β under a wide variety of turbulence conditions, and β could move closer to 1 for strong ocean turbulence.Based on Equations ( 5)-( 8), the scintillation indexes and Málaga model parameters (α, β) were calculated when the link lengths were 26 m and 30 m, respectively, as shown in Table 3.They were used in the following simulation.
Based on Equation ( 15) and the ocean turbulence parameters from Argo node 4, the effects of the long-term beam spreading coefficient W z /a, and the received aperture diameter D rx on the received optical power with the pointing error displacements are shown in Figure 4.It can be seen that the received optical power gradually increased with the decrease in W z /a and the increase in D rx .According to Figure 4, W z /a is larger in strong turbulence than in weak turbulence when D rx is fixed.The received optical power at the finite receiving aperture is smaller in strong turbulence than in weak turbulence.Therefore, when determining the optimal received power for the optical signal, increasing D rx can mitigate the effects of pointing errors.
α and β with increases in z and χΤ is shown.It can be clearly seen that the value of α was always larger than that of β under a wide variety of turbulence conditions, and β could move closer to 1 for strong ocean turbulence.Based on Equations ( 5)-( 8), the scintillation indexes and Málaga model parameters (α, β) were calculated when the link lengths were 26 m and 30m, respectively, as shown in Table 3.They were used in the following simulation.α and β with increases in z and χΤ is shown.It can be clearly seen that the value of α was always larger than that of β under a wide variety of turbulence conditions, and β could move closer to 1 for strong ocean turbulence.Based on Equations ( 5)-( 8), the scintillation indexes and Málaga model parameters (α, β) were calculated when the link lengths were 26 m and 30m, respectively, as shown in Table 3.They were used in the following simulation.(1 × 10 −6 , 10 eter Drx on the received optical power with the pointing error displacements are shown in Figure 4.It can be seen that the received optical power gradually increased with the decrease in Wz/a and the increase in Drx.According to Figure 4, Wz/a is larger in strong turbulence than in weak turbulence when Drx is fixed.The received optical power at the finite receiving aperture is smaller in strong turbulence than in weak turbulence.Therefore, when determining the optimal received power for the optical signal, increasing Drx can mitigate the effects of pointing errors.).

System Average BER
First, we illustrate the effect of the turbulence strength on the system average BER with changes in the normalized pointing error coefficient σp/a or the average received SNR.These three ocean turbulence scenarios were set from Table 3.In Figure 5, the selected ARGO node numbers are one, two, and three, respectively.As the link length is 30 m, their scintillation indexes are 0.770, 1.523, and 2.643, respectively.It can be seen in Figure 5a that the system average BER shows an upward trend with increasing σp/a or ocean turbulence strength.The numerical results are consistent with the analytical results when σp/a > 2. In Figure 5b, when the link length is 26 m, the scintillation indexes for the three nodes are 0.553, 1.102, and 1.860, respectively.The numerical results are consistent with the analytical results for a high SNR when ξ 2 < 1.

System Average BER
First, we illustrate the effect of the turbulence strength on the system average BER with changes in the normalized pointing error coefficient σ p /a or the average received SNR.These three ocean turbulence scenarios were set from Table 3.In Figure 5, the selected ARGO node numbers are one, two, and three, respectively.As the link length is 30 m, their scintillation indexes are 0.770, 1.523, and 2.643, respectively.It can be seen in Figure 5a that the system average BER shows an upward trend with increasing σ p /a or ocean turbulence strength.The numerical results are consistent with the analytical results when σ p /a > 2. In Figure 5b, when the link length is 26 m, the scintillation indexes for the three nodes are 0.553, 1.102, and 1.860, respectively.The numerical results are consistent with the analytical results for a high SNR when ξ 2 < 1.
Secondly, to examine a QAM order m of 8, 16, or 32 with an increase in the average received SNR, we selected Argo node 4 and determined the effects of the normalized pointing error deviations σ p /a and the normalized pointing error coefficient W z /a on the system average BER, respectively.As shown in Figure 6, the system average BER performance could be weakened with increases in σ p /a, W z /a, or m.According to the numerical results shown in Figure 4, long-term beam spreading increases with increasing turbulence strength.Thus, the above results are consistent with the conclusion that the system average BER increases with increasing turbulence strength or decreasing received aperture diameter.
Secondly, to examine a QAM order m of 8, 16, or 32 with an increase in the average received SNR, we selected Argo node 4 and determined the effects of the normalized pointing error deviations σp/a and the normalized pointing error coefficient Wz/a on the system average BER, respectively.As shown in Figure 6, the system average BER performance could be weakened with increases in σp/a, Wz/a, or m.According to the numerical results shown in Figure 4, long-term beam spreading increases with increasing turbulence strength.Thus, the above results are consistent with the conclusion that the system average BER increases with increasing turbulence strength or decreasing received aperture diameter.Thirdly, to examine a QAM order m of 8, 16, or 32 with an increase in the received aperture diameter Drx, we selected node 4 and verified the effects of σp/a and Wz/a on the system average BER, respectively.The numerical results are shown in Figure 7.In Figure 7, when the link length is 26 m, the scintillation index of node 4 is 1.670.In Figure 7a, m = 32 and, when σp/a is 2, 3, or 4, the received aperture diameter Drx that the system designs for the average BER is less than the forward error correction (FEC) threshold (3.8 × 10 −3 ), Thirdly, to examine a QAM order m of 8, 16, or 32 with an increase in the received aperture diameter D rx , we selected node 4 and verified the effects of σ p /a and W z /a on the system average BER, respectively.The numerical results are shown in Figure 7.In Figure 7, when the link length is 26 m, the scintillation index of node 4 is 1.670.In Figure 7a, m = 32 and, when σ p /a is 2, 3, or 4, the received aperture diameter D rx that the system designs for the average BER is less than the forward error correction (FEC) threshold (3.8 × 10 −3 ), which is 6.1 mm, 8.1 mm, or 9 mm.Similarly, in Figure 7b, W z /a = 8 when the QAM order m is 8, 16, or 32, and D rx is chosen as 5.7 mm, 6.6 mm, or 6.9 mm.Obviously, the effect of σ p /a on the system average BER is greater than that of W z /a. which is 6.1 mm, 8.1 mm, or 9 mm.Similarly, in Figure 7b, Wz/a = 8 when the QAM order m is 8, 16, or 32, and Drx is chosen as 5.7 mm, 6.6 mm, or 6.9 mm.Obviously, the effect of σp/a on the system average BER is greater than that of Wz/a.

System Outage Performance
In this section, due to the long simulation time involved, only the Monte-Carlo simulation results up to Pout = 10 −8 are included.The scintillation indexes for the three selected Argo nodes are 0.553, 1.102, and 1.670, respectively, when the link length is 26 m.For a QAM order m of 8, 16, or 32, respectively, and an increase in the received aperture diameter Drx, we separately verified the effects of the normalized pointing error coefficient σp/a and the turbulence strength on the system outage performance.The numerical results are shown in Figure 9.It can be seen from Figure 9 that the outage probability increases gradually with increases in σp/a, the QAM order, and the turbulence strength.Moreover, the outage probability decreases with the increase in Drx.In brief, the effect of turbulence strength on the outage probability is greater than that of σp/a.For example, based on the numerical results shown in Figure 9, if the system is designed to satisfy Pout = 10 −3 , then the Drx that can be reliably transmitted over this channel is 8 mm. ).

System Outage Performance
In this section, due to the long simulation time involved, only the Monte-Carlo simulation results up to P out = 10 −8 are included.The scintillation indexes for the three selected Argo nodes are 0.553, 1.102, and 1.670, respectively, when the link length is 26 m.For a QAM order m of 8, 16, or 32, respectively, and an increase in the received aperture diameter D rx , we separately verified the effects of the normalized pointing error coefficient σ p /a and the turbulence strength on the system outage performance.The numerical results are shown in Figure 9.It can be seen from Figure 9 that the outage probability increases gradually with increases in σ p /a, the QAM order, and the turbulence strength.Moreover, the outage probability decreases with the increase in D rx .In brief, the effect of turbulence strength on the outage probability is greater than that of σ p /a.For example, based on the numerical results shown in Figure 9, if the system is designed to satisfy P out = 10 −3 , then the D rx that can be reliably transmitted over this channel is 8 mm.

System Outage Performance
In this section, due to the long simulation time involved, only the Monte-Carlo simulation results up to Pout = 10 −8 are included.The scintillation indexes for the three selected Argo nodes are 0.553, 1.102, and 1.670, respectively, when the link length is 26 m.For a QAM order m of 8, 16, or 32, respectively, and an increase in the received aperture diameter Drx, we separately verified the effects of the normalized pointing error coefficient σp/a and the turbulence strength on the system outage performance.The numerical results are shown in Figure 9.It can be seen from Figure 9 that the outage probability increases gradually with increases in σp/a, the QAM order, and the turbulence strength.Moreover, the outage probability decreases with the increase in Drx.In brief, the effect of turbulence strength on the outage probability is greater than that of σp/a.For example, based on the numerical results shown in Figure 9, if the system is designed to satisfy Pout = 10 −3 , then the Drx that can be reliably transmitted over this channel is 8 mm.

Conclusions
In this paper, we explored the system average BER and the outage probability for the Málaga turbulence channel with pointing errors and beam spreading in MIMO-UWOC.We also analyzed the relationships between the long-term beam spreading and the turbulence strength and between the diversity gain and the pointing errors.In the numerical simulation, we employed four Argo nodes' observed ocean data, including the pressure, temperature, and salinity of seawater.We illustrated the effects of the pointing errors and the ocean turbulence on the average BER, considering pointing errors with long-term beam spreading, the average received SNR, and the received aperture diameter.In addition, we analyzed their effects on the system outage performance, taking into account the received aperture diameter.The numerical results showed that the system average BER and the outage probability increased with increasing turbulence strength and pointing errors.Moreover, in the mQAM-MIMO system with MRC, the pointing errors had a larger effect on the average BER than the turbulence strength, which had a larger effect on the outage performance than the pointing errors.Therefore, increasing the system receiving aperture diameter to mitigate long-term beam spreading and ocean turbulence can better improve UWOC reliability and efficiency than adjusting the QAM order.

Figure 5 .
Figure 5. System average BER for three different Argo nodes: (a) changing the normalized jitter, σp/a; (b) changing the average received SNR.Figure 5. System average BER for three different Argo nodes: (a) changing the normalized jitter, σ p /a; (b) changing the average received SNR.

Figure 5 .
Figure 5. System average BER for three different Argo nodes: (a) changing the normalized jitter, σp/a; (b) changing the average received SNR.Figure 5. System average BER for three different Argo nodes: (a) changing the normalized jitter, σ p /a; (b) changing the average received SNR.

Figure 5 .
System average BER for three different Argo nodes: (a) changing the normalized jitter, σp/a; (b) changing the average received SNR.

Figure 7 .
Figure 7. System average BER with different received aperture diameters for Argo node 4: (a) QAM order m = 8, 16, and 32 and σp/a = 2, 3, and 4; (b) QAM order m = 8, 16, and 32 and Wz/a = 8, 9, and 10.Finally, we verified the coding gain with the different link lengths in three Argo nodes.The numerical results are shown in Figure8.Based on the results shown in Figure2, the turbulence strength increases with the increase in link length.Obviously, the coding gain is weakened due to the pointing errors and the turbulence strength.

Table 1 .
System parameters in the simulation.

Table 2 .
Selected Argo Node information.