{'fraction_non_alphanumeric': 0.04056785812735459, 'fraction_numerical': 0.014288339612239272, 'mean_word_length': 4.175106990014266, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Tb/s Indoor MIMO Optical Wireless Backhaul System Using VCSEL Arrays

In this paper, the design of a multiple-input multiple-output (MIMO) optical wireless communication (OWC) link based on vertical cavity surface emitting laser (VCSEL) arrays is systematically carried out with the aim to support data rates in excess of 1 Tb/s for the backhaul of sixth generation (6G) indoor wireless networks. The proposed design combines direct current optical orthogonal frequency division multiplexing (DCO-OFDM) and a spatial multiplexing MIMO architecture. For such an ultra-high-speed line-of-sight (LOS) OWC link with low divergence laser beams, maintaining alignment is of high importance. In this paper, two types of misalignment error between the transmitter and receiver are distinguished, namely, radial displacement error and orientation angle error, and they are thoroughly modeled in a unified analytical framework assuming Gaussian laser beams, resulting in a generalized misalignment model (GMM). The derived GMM is then extended to MIMO arrays and the performance of the MIMO-OFDM OWC system is analyzed in terms of the aggregate data rate. Novel insights are provided into the system performance based on computer simulations by studying various influential factors such as beam waist, array configuration and different misalignment errors, which can be used as guidelines for designing short range Tb/s MIMO OWC systems.


I. INTRODUCTION
The proliferation of Internet-enabled premium services such as 4K and 8K ultra-high-definition (UHD) video streaming, immersion into virtual reality (VR) or augmented reality (AR) with three dimensional (3D) stereoscopic vision, holographic telepresence and multi-access edge computing will extremely push wireless connectivity limits in years to come [2]. These technologies will require an unprecedented system capacity above 1 Tb/s for real-time operation, which is one of the key performance indicators of the future sixth generation (6G) wireless systems [3].
The achievability of ultra-high transmission rates of Tb/s has been addressed in the literature for both wired and wireless systems [4]- [6]. Targeting digital subscriber line (DSL) applications, in [4], Shrestha et al. have used a two-wire copper cable as a multi-mode waveguide for multipleinput multiple-output (MIMO) transmission and experimentally measured the received power for signals with 200 GHz bandwidth. They have predicted that aggregate data rates of several Tb/s over a twisted wire pair are feasible at short distances of ≤ 10 m by using discrete multitone (DMT) modulation and vector coding. In [5], the authors have reported successful implementation of a 1 Tb/s super channel over a 400 km optical single mode fiber (SMF) link based on quadrature amplitude modulation (QAM) with probabilistic constellation shaping. In [6], Petrov et al. have elaborated on a roadmap to actualize last meter indoor broadband wireless access in the terahertz band, i.e. 0.1-10 THz, in order to enable Tb/s connectivity between the wired backbone infrastructure and personal wireless devices.
The feasibility of Tb/s data rates has been actively studied for outdoor point-to-point free space optical (FSO) communications [7]- [9]. In [7], Ciaramella

Beam-steered infrared (IR) light communication is an optical wireless communication (OWC)
technology primarily tailored to provide multi-Gb/s data rates per user by directing narrow IR laser beams at mobile devices. In [10], Koonen  Indoor laser-based optical wireless access networks can generate aggregate data rates beyond 1 Tb/s [12]. Such ultra-high-speed indoor access networks impose a substantial overhead on the backhaul capacity, and a cost-effective backhaul solution is a major challenge. In this paper, a high-capacity wireless backhaul system is designed based on laser-based OWC to support aggregate data rates of at least 1 Tb/s for backhaul connectivity in next generation Tb/s indoor networks. While FSO systems suffer from outdoor channel impairments such as weather-dependent absorption loss and atmospheric turbulence, short range laser-based OWC under stable and acclimatized conditions of indoor environments potentially enhances the signal quality. This way, the need for bulky FSO transceivers equipped with expensive subsystems to counteract outdoor effects is eliminated. Moreover, the aforementioned FSO systems use dense wavelength division multiplexing (DWDM) to deliver Tb/s data rates, which significantly increases the cost and complexity of the front-end system. Different from WDM FSO systems, in this paper, a single wavelength is used to achieve a data rate of ≥ 1 Tb/s by means of vertical cavity surface emitting lasers (VCSELs). The choice of VCSELs for the optical wireless system design is motivated by the fact that, among various types of laser diodes, VCSELs are one of the strongest contenders to fulfil this role due to several important features of them, including [13], [14]: 1) a high modulation bandwidth of ≥ 10 GHz; 2) a high power conversion efficiency of > 50%; 3) cost-efficient fabrication by virtue of their compatibility with large scale integration processes; 4) possibility for multiple devices to be densely packed and precisely arranged as 2D arrays. These attributes make VCSELs appealing to many applications such as optical networks, highly parallel optical interconnects and laser printers, to name a few [15]. Single mode VCSELs, which are the focus of this paper, generate an output optical field in the fundamental transverse electromagnetic mode (TEM) (i.e. TEM 00 mode), resulting in a Gaussian profile on the transverse plane, in that the optical power is maximum at the center of the beam spot and it decays exponentially with the squared radial distance from the center [16].
For line-of-sight (LOS) OWC links, accurate alignment between the transmitter and receiver is a determining factor of the system performance and reliability. In principle, two types of misalignment may occur in the link: radial displacement between the transmitter and receiver, and orientation angle error at the transmitter or receiver side. Modeling of the Gaussian beam misalignment has been addressed in the context of terrestrial FSO systems such as the works of Farid and Hranilovic, for single input single output (SISO) [17] and MIMO [18] links. The FSO transceiver equipment is commonly installed on the rooftops of high-rise buildings and hence random building sways due to wind loads and thermal expansions cause a pointing error in the transmitter orientation angle with independent and identical random components in elevation and horizontal directions [19]. The works in [17]- [19] implicitly base their modeling methodology on the assumption of treating the effect of this angle deviation at the transmitter (with a typical value of 1 mrad) as a radial displacement of the beam spot position at the receiver (typically located at 1 km distance from the transmitter). By contrast, in [20], Huang and Safari, through applying a small angle approximation, have modeled the receiver-induced angle-of-arrival (AOA) misalignment again as a radial displacement of the optical field pattern on the photodetector (PD) plane. In [21], [22] Although the problem of misalignment constitutes an important practical challenge for designing Tb/s indoor optical wireless links as a key enabler for ultra-high-speed networks of the future, it has not been widely studied yet. For short range indoor OWC systems with compact PDs, to minimize the geometric loss, the beam spot size is required to be relatively small, comparable to the size of a PD, in which case angular misalignment can significantly influence the link performance, independent of the radial displacement error. In a previous work [1], the authors have presented preliminary results to study the effect of only displacement error on the performance of indoor Tb/s MIMO OWC. However, such high performance indoor OWC links are essentially prone to any type of misalignment errors. To the best of the authors' knowledge, there is a lack of a comprehensive and analytically tractable model of the link misalignment for laser-based OWC systems inclusive of orientation angle errors at the transmitter and receiver sides as well as the radial displacement error. This paper puts forward the modeling and design of a spatial multiplexing MIMO OWC system based on direct current-biased optical orthogonal frequency division multiplexing (DCO-OFDM) and VCSEL arrays to unlock Tb/s data rates with single mode VCSELs. The contributions of this paper are concisely given as follows: • An in-depth analytical modeling of the misalignment for SISO optical wireless channels with Gaussian beams is presented. Thereupon, a generalized misalignment model (GMM) is derived where radial displacement and orientation angle errors at both the transmitter and receiver sides are all taken into consideration in a unified manner. This model is verified by computer simulations using a commercial optical design software, Zemax OpticStudio.
• The GMM derivation is extended to MIMO OWC systems with arbitrary configurations for transmitter and receiver arrays. The geometric modeling of the VCSEL and PD arrays is explicated by highlighting critical design parameters such as array size, spacing between array elements and the size of PDs.
• The spatial multiplexing MIMO-OFDM transceiver under consideration is elucidated and the received signal-to-interference-plus-noise ratio (SINR) and aggregate data rate are analyzed.
It is shown that treating an angular pointing error of the transmitter as a radial displacement is a special case of the GMM and a tight analytical approximation of the MIMO channel direct current (DC) gains is derived for this case.
The remainder of the paper is organized as follows. In Section II, the SISO channel model for a perfectly aligned VCSEL-based OWC system is described. In Section III, the detailed analytical modeling framework for the generalized misalignment of the SISO channel is established. In Section IV, the design and analysis of the MIMO-OFDM OWC system using VCSEL and PD arrays is presented, including the incorporation of the GMM in the MIMO channel model. In Section V, numerical results are provided. In Section VI, concluding remarks are drawn and a II. VCSEL-BASED OPTICAL WIRELESS CHANNEL A MIMO OWC system can be realized by means of an array of VCSELs. Fig. 1 illustrates a 4 × 4 MIMO system configuration, comprising a 2 × 2 VCSEL array and a 2 × 2 PD array, which are perfectly aligned to each other. There is a directed LOS link from every VCSEL to its corresponding PD as represented by an exclusive color. This section deals with the channel model for a VCSEL-based SISO optical wireless link as a building block for the MIMO system.

A. Gaussian Beam Propagation
The wavefront of the Gaussian beam is initially planar at the beam waist and then expanding in the direction of propagation. The wavefront radius of curvature at distance z from the transmitter is characterized by [16]: where w 0 is the waist radius; and λ is the laser wavelength. Also, the radius of the beam spot, which is measured at the 1 e 2 normalized intensity contour on the transverse plane, takes the following form [16]: where z R is the Rayleigh range. It is defined as the distance at which the beam radius is extended by a factor of √ 2, i.e. w(z R ) = √ 2w 0 . In this case, the beam spot has an area twice that of the beam waist. The Rayleigh range is related to w 0 and λ via [16]: From (2), w(z) for z ≫ z R approaches the asymptotic value: thus varying linearly with z. Therefore, the circular beam spot in far field is the base of a cone whose vertex lies at the center of the beam waist with a divergence half-angle: The spatial distribution of a Gaussian beam along its propagation axis is described by the intensity profile on the transverse plane. By using Cartesian coordinates, the intensity distribution at distance z from the transmitter at the point (x, y) is given by [16]: where P t is the transmitted optical power; and ρ(x, y) is the Euclidean distance of the point (x, y) from the center of the beam spot.

B. Channel DC Gain
As observed from (6), the fundamental property of the Gaussian beam propagation is that the intensity distribution is given on the transverse plane as a circularly symmetric function centered about the beam axis. For a general link configuration, the transmitter and receiver are oriented toward arbitrary directions in the 3D space with n t and n r denoting their normal vectors, respectively. In this case, the received optical power P r is obtained by integrating (6) over the effective area of the PD [24]. To this end, the PD area is projected onto the transverse plane by using the cosine of the normal vector of the PD plane with respect to the beam propagation axis, which is equal to n t · n r . The DC gain of an intensity modulation and direct detection (IM-DD) channel is defined as the ratio of the average optical power of the received signal to that of the transmitted signal. Therefore, the DC gain of the channel is calculated as follows: z' x' y' where A denotes the PD area; and the factor n t ·n r accounts for Lambert's cosine law. Throughout the paper, a circular PD of radius r PD is assumed for which A = {(x, y) ∈ R 2 | x 2 + y 2 ≤ r 2 PD }. Fig. 2 illustrates a SISO OWC system in a directed LOS configuration with perfect alignment.
In this case, the beam waist plane is parallel to the PD plane so that n t · n r = 1 and the center of the beam spot is exactly located at the center of the PD. Hence, ρ(x, y) in (6) is equal to r = x 2 + y 2 on the PD plane. From (7), for a link distance of z = L, the DC gain of the channel becomes: where

III. GENERALIZED MISALIGNMENT MODELING
This section establishes a mathematical framework for the analytical modeling of misalignment errors for the SISO optical wireless channel discussed in Section II. In the following, two cases of displacement error and orientation angle error are presented.

A. Displacement Error
A displacement error between the transmitter and receiver causes the center of the beam spot at the PD plane to deviate radially, relative to the center of the PD, which is equivalent to the radial displacement in [17]. In this case, n t · n r = 1. The magnitude of the error is represented , where x DE and y DE correspond to the error components along the x ′ and y ′ axes, as shown in Fig. 3. It can be observed that the intensity value depends on the axial distance z between the beam waist and the PD plane where z = L, and the Euclidean distance ρ from the center of the beam spot to the coordinates (x, y). It follows that: Substituting (9) in (7), the DC gain of the channel turns into:

B. Orientation Angle Error
An orientation error occurs when the transmitter or receiver has a non-zero tilt angle with respect to the alignment axis. Orientation angles of the transmitter and receiver, denoted by φ and ψ, respectively, entail arbitrary and independent directions in the 3D space. Note that the transmitter and receiver orientation errors are jointly modeled, though they are separately depicted in Fig. 4 to avoid intricate geometry. Furthermore, the angles φ and ψ are decomposed into azimuth and elevation components in the 3D space using the rotation convention shown in   beam axis, with their centers lying on the beam axis. Among them, the one with a circumference intersecting the point (x, y) on the PD plane is the basis for analysis. This particular disk is referred to as principal disk hereinafter, which is drawn as a yellow disk in Fig. 4. The variables to be characterized are the axial distance z between the beam waist and the center of the principal disk and the Euclidean distance ρ of the point (x, y) to the beam axis, i.e. the radius of the principal disk.
The PD's x ′′ y ′′ z ′′ coordinate system is rotated with respect to the reference coordinate system x ′ y ′ z ′ as shown in Fig. 4(b). Based on Euler angles with clockwise rotations, the x ′ y ′ z ′ system is transformed into the x ′′ y ′′ z ′′ system by rotating first about the y ′ axis through an angle −ψ a , then about the x ′′ axis through an angle −ψ e , using the following rotation matrices: for β = −ψ a and α = −ψ e . The desired point [x y 0] ⊤ is given in the x ′′ y ′′ z ′′ system. The projected coordinates of this point in the x ′ y ′ z ′ system is obtained as: The axial distance z consists of two segments including the projection of L onto the beam axis and the additive length ℓ: From Fig. 4(a), there are two parallel planes indicated by P 1 and P 2 so that the additive length ℓ can be found as the distance between P 1 and P 2 , which is subsequently derived in (16). These planes are perpendicular to the beam axis such that P 1 passes through the point [u v w] ⊤ in the x ′ y ′ z ′ system and P 2 crosses the origin. The normal vector of both planes is: where n x ′ , n y ′ and n z ′ represent unit vectors for x ′ , y ′ and z ′ axes, respectively. Let n t = [a b c] ⊤ to simplify notation, where n t = √ a 2 + b 2 + c 2 = 1. It follows that: The additive length ℓ in (13) is derived by finding the distance from the origin to P 1 , resulting in: Combining (16) with (12) and (14), and using trigonometric identities, yields: The squared radius of the principal disk illustrated in Fig. 4(a) is given by: where: Substituting u, v and w from (12) into (19), and simplifying, leads to: The last piece required to complete the analysis of the channel gain based on (7) is the calculation of the inner product of n t and n r . From Fig. 4(b), the normal vector to the PD surface is: By using (14) and (21), the cosine of the planar angle between the surface normal and the beam axis is obtained as follows: By combining (2), (13), (17), (18) and (20), the DC gain of the channel, denoted by H 0 (φ a , φ e , ψ a , ψ e ), can be evaluated based on (7) and (22) when using: x cos(φ e ) sin(φ a − ψ a ) + y(sin(ψ e ) cos(φ e ) cos(φ a − ψ a ) + cos(ψ e ) sin(φ e )) 2 .

C. Unified Misalignment Model
In order to unify displacement and orientation errors, after the transmitter is rotated, it is Referring to the parallel planes P 1 and P 2 in (15), Consequently, the squared radius of the principal disk is determined by using (18) in conjunction with (13) and is computed based on (7) and (22) through the use of: x'

A. Structure of Arrays and MIMO Channel
The array structure shown in Fig. 1 is extended to a K ×K square, forming an N t ×N r MIMO OWC system where 1 N t = N r = K 2 . Fig. 5 depicts a K × K VCSEL array and a K × K PD array on the x ′ y ′ plane. The gap between adjacent elements of the PD array is controlled by δ > 0, which is referred to as inter-element spacing hereinafter. For those PDs that are close to the edges of the array, there is a margin of δ 2 with respect to the edges. The center-to-center distance for neighboring PDs along rows or columns of the array is: The side length for each array is W = Kd PD , leading to array dimensions of W × W .
The MIMO channel is identified by an N r × N t matrix of DC gains for all transmission paths between the transmitter and receiver arrays: where the entry H ij corresponds to the link from VCSEL j to PD i . For the array structure shown in Fig. 5, the elements are labeled by using a single index according to their row and column indices. This way, for an N r × N t = K 2 × K 2 array, the VCSEL (resp. PD) situated at the the ith element on each array. From Fig. 5, it is straightforward to show that: where m = ⌈ i K ⌉ and n = i − ⌈ i K ⌉ − 1 K, with ⌈q⌉ denoting the smallest integer that satisfies ⌈q⌉ ≥ q. In this case, evaluating H ij based on (7) leads to:

B. Generalized Misalignment of Arrays
Under the generalized misalignment, the whole transmitter and receiver arrays are affected by both displacement and orientation errors. To characterize the MIMO channel with the generalized misalignment, the entries of the MIMO channel matrix can be derived by applying the procedure y' x' Fig. 6: Orientation angle errors of the transmitter and receiver in a MIMO system configuration.
The parameters of the LOS channel between VCSEL j and PD i are highlighted in blue.
described in Section III-C for the SISO channel 2 . Fig. 6 depicts the MIMO link configuration where the transmitter and receiver arrays are rotated according to their respective orientation angles. Note that Fig. 6 does not include radial displacement for simplicity. The corresponding parameters of the SISO link between VCSEL j and PD i are highlighted in blue, as shown in Fig. 6.
Considering the generalized misalignment, the VCSEL array is first rotated by an angle φ and then its center is radially displaced relative to the center of the receiver array. The coordinates 2 In consideration of the receiver optics, a non-imaging optical design based on ideal compound parabolic concentrators (CPCs) can be used. Such a receiver structure would not alter the spatial distribution of the beam spots at the receiver array, since an ideal lossless CPC transfers all the optical power collected at the entrance aperture to the exit aperture [24]. Therefore, the analytical framework proposed for the generalized misalignment modeling in Section III would still be valid for computing the optical power incident on each receiver element due to the fact that all the integrations can be equally evaluated over the circular area of the entrance aperture of CPCs. The optimal design of non-imaging receivers is, however, beyond the scope of this paper and it is subject to a separate study.
of VCSEL j in the reference coordinate system x ′ y ′ z ′ are:  where R y ′ (−φ a ) and R x ′′ (φ e ) are given by (11) for β = −φ a and α = φ e . Also, after the receiver array undergoes a rotation by an angle ψ, the coordinates of PD i in the x ′ y ′ z ′ system are: where R y ′ (−ψ a ) and R x ′′ (ψ e ) are given by (11) for β = −ψ a and α = −ψ e . To calculate the channel gain between VCSEL j and PD i , denoted by H ij (x DE , y DE , φ a , φ e , ψ a , ψ e ), using (7), the parameters w(z) and ρ(x, y) are evaluated based on (25) and (26) by substituting the link distance L ′ =ž j −ẑ i for L, and the displacement components x ′ DE =x j −x i and y ′ DE =y j −ŷ i for x DE and y DE , respectively. This exact procedure is referred to as the MIMO GMM for brevity.

C. Approximation of the MIMO GMM
The computation of the MIMO GMM as described above entails numerical integrations. In the following, approximate analytical expressions of the MIMO channel gain H ij are derived for two special cases of radial displacement and orientation error at the transmitter. Then, the relation between them for a small angle error is elaborated. The area of a circular PD of radius r PD is approximated by an equivalent square of side length a PD = √ πr PD with the same area.

1) Radial Displacement:
In this case, Therefore, z = L, n t · n r = 1 and: From (7), H ij is then written as: which can be derived as follows: where erf(t) = 2 √ π t 0 e −s 2 ds is the error function. 2) Orientation Error of the Transmitter: For the case of orientation angle error at the transmitter, the use of (31) and (32) leads tox j −x i = x j cos(φ a ) − y j sin(φ e ) sin(φ a ) − x i ,y j −ŷ i = y j cos(φ e )−y i andž j −ẑ i = L+x j sin(φ a )+y j sin(φ e ) cos(φ a ). After simplifying, the parameters w 2 (z) and ρ 2 (x, y) are obtained as: Considering that x + x i ≪ L and y + y i ≪ L hold, for sufficiently small values of φ a and φ e , (x + x i ) cos(φ e ) sin(φ a ) + (y + y i ) sin(φ e ) ≪ L cos(φ e ) cos(φ a ), which gives rise to: This approximation means the axial distance variation of the slightly tilted beam spot over the PD surface is ignored due to its small size. Hence, (36) is simplified to w 2 (z) ≈ w 2 (L cos(φ e ) cos(φ a )).
Besides, when φ a and φ e are small enough, in the right hand side of (37), the last term is deemed negligible compared to the first two terms from the factor sin(φ e ) sin(φ a ) ≪ 1. Consequently, using n t · n r = cos(φ e ) cos(φ a ), the integration in (7) is approximated by: A closed form solution of (39) is readily derived as follows:  and L sin(φ e ) cos(φ a ), respectively. The approximations in (35) and (40) are verified in Appendix and it is shown that they are highly accurate as long as the beam spot size at the receiver is in the order of or larger than the PD size.

D. Spatial Multiplexing MIMO-OFDM
1) Transceiver Architecture: In this paper, a spatial multiplexing MIMO-OFDM system is used to maximize the transmission rate by sending independent data streams over the MIMO optical wireless channel. To ensure a high spectral efficiency for transmitting each data stream, the standard DCO-OFDM technique along with adaptive QAM is used [25]. In this case, the spatial overlapping between the beam spots of the VCSELs at the PD array causes crosstalk in the MIMO channel. Fig. 7 depicts a simplified block diagram of the MIMO-OFDM transceiver architecture. The singular value decomposition (SVD) precoding and decoding stages shown in Fig. 7 can be applied provided the MIMO channel state information (CSI) is known at both the transmitter and receiver. First, it is assumed that this is not the case to avoid the overhead associated with the CSI estimation and feedback, and the transceiver system is described without the use of SVD.
At the transmitter, the N t input binary data streams are individually mapped to a sequence of complex QAM symbols. With a digital realization of the OFDM modulation and demodulation by way of N FFT -point inverse fast Fourier transform (IFFT) and fast Fourier transform (FFT), respectively, the resulting sequences are buffered into blocks of size N t ×N QAM . They are loaded onto the N QAM data-carrying subcarriers of the N t OFDM frames in positive frequencies, where For baseband OFDM transmission in IM-DD systems, it is necessary for the time domain signal to be real-valued. To this end, for each OFDM frame, the number of symbols is extended to N FFT according to a Hermitian symmetry and the DC and Nyquist frequency subcarriers are zero-padded before the IFFT operation. Also, in order to comply with the nonnegativity constraint of IM-DD channels, a proper DC level is added in the time domain to obtain a positive signal [25]. Let x(t) = [x 1 (t) x 2 (t) · · · x Nt (t)] ⊤ be the vector of instantaneous optical powers emitted by the VCSELs at time sample t for t = 0, 1, . . . , N FFT − 1. It is given by: where P elec is the average electrical power of each OFDM symbol; is the vector of the normalized discrete time OFDM samples; x DC = P t is the DC bias with P t representing the average optical power per VCSEL; and 1 ∈ R Nt×1 is an all-ones vector.
The finite dynamic range of the VCSELs determines the available peak-to-peak swing for their modulating OFDM signal. The envelope of the unbiased OFDM signal follows a zero mean real Gaussian distribution for N FFT ≥ 64 [26]. The choice of P t = 3 √ P elec guarantees that 99.7% of the signal variations remains undistorted, thereby effectively discarding the clipping noise [27].
Thus, the average power of the OFDM signal assigned to each VCSEL is P elec = 1 9 P 2 t . At the receiver array, after filtering out the DC component and perfect sampling, the vector of received photocurrents is obtained. LetX k ∈ R Nt×1 be the vector of symbols modulated on the kth subcarrier in the frequency domain. After the FFT operation, the received symbols are extracted from the data-carrying subcarriers and then they are demodulated using maximum likelihood detection. The vector of received signals on the kth subcarrier for k = 0, 1, . . . , N FFT − 1 is written in the form: where R PD is the PD responsivity; H k is the frequency response of the MIMO channel; and Z k is the additive white Gaussian noise (AWGN) vector. Note that without SVD processing, N r = N t holds, in which case Y k ∈ R Nt×1 and H k ∈ R Nt×Nt . Considering strong LOS components when using laser beams with low divergence, the channel is nearly flat for which H k = H 0 ∀k, where H 0 refers to (28). The ith element of the noise vector comprises thermal noise and shot noise of the ith branch of the receiver and the relative intensity noise (RIN) caused by all the VCSELs which depends on the average received optical power [28]. The total noise variance is given by: where κ is the Boltzmann constant; T is temperature in Kelvin; R L is the load resistance; B is the single-sided bandwidth of the system; F n is the noise figure of the transimpedance amplifier (TIA); q is the elementary charge; and RIN is defined as the mean square of instantaneous power fluctuations divided by the squared average power of the laser source [28]. Based on (42), the received SINR per subcarrier for the ith link is derived as follows: 2) SVD Processing: When the CSI is available at the transmitter and receiver, the MIMO channel can be transformed into a set of parallel independent subchannels by means of SVD of the channel matrix in the frequency domain. The use of SVD leads to the capacity achieving architecture for spatial multiplexing MIMO systems [29]. The SVD of H k ∈ R Nr×Nt , with where U k ∈ R Nr×Nr and V k ∈ R Nt×Nt are unitary matrices; * denotes conjugate transpose; and Λ k ∈ R Nr×Nt is a rectangular diagonal matrix of the ordered singular values, i.e. λ 1 ≥ λ 2 ≥ · · · ≥ λ Nt > 0 [29]. Note that H k = H 0 ∀k as discussed so the subscript k can be dropped from the singular values. After SVD decoding at the receiver, the N t -dimensional vector of received symbols on the kth subcarrier for k = 0, 1, . . . , N FFT − 1 becomes:Ỹ whereZ k = U * k Z k . Note that the statistics of the noise vector is preserved under a unitary transformation. Therefore, the ith elements ofZ k and Z k have the same variance of σ 2 i . The received signal-to-noise ratio (SNR) per subcarrier for the ith link is derived as follows: (46) 21

E. Aggregate Rate Analysis
Based on adaptive QAM with a given bit error ratio (BER) performance in an AWGN channel, a tight upper bound for the number of bits per symbol transmitted by VCSEL i , for BER ≤ 10 −2 and 0 ≤ γ i ≤ 30 dB, is given by [30]: where: represents the SINR gap due to the required BER. With a symbol rate of 2B N FFT for DCO-OFDM, the achievable rate per subcarrier is 2B N FFT η i bit/s. According to N FFT 2 −1 data-carrying subcarriers, the transmission rate for VCSEL i is then: where ξ = N FFT −2 N FFT . Hence, the aggregate data rate of the MIMO-OFDM system is expressed as: 22

V. NUMERICAL RESULTS AND DISCUSSIONS
The performance of the VCSEL-based MIMO OWC system is evaluated by using computer simulations and the parameters listed in Table I, where the VCSEL and noise parameters are adopted from [15], [31]. Numerical results are presented for the effective radius w 0 of the beam waist over the range 10 µm ≤ w 0 ≤ 100 µm with the assumption that there is a convex lens next to each bare VCSEL to widen its output beam waist in order to reduce the far field beam divergence. For a link distance of L = 2 m, the beam spot radius and divergence angle vary For a given wavelength, the MPE value reduces with an increase in the beam waist radius [33], so the most restrictive case is when w 0 is at a maximum. The case of w 0 = 100 µm and λ = 850 nm leads to a subtense angle of α < 1.5 mrad. The subtense angle is measured as the plane angle subtended by the apparent source to the observer's pupil at MHP. For α < 1.5 mrad, the laser source is classified as a point source for which η = 1 and MPE = 50.8 W/m 2 [33]. For a circular aperture of diameter 7 mm, A pupil ≈ 38.5 mm 2 , and it follows that P t ≤ 1.95 mW for each VCSEL. Hence, P t = 1 mW is considered eye-safe.
A. Perfect Alignment 1) Spatial Distribution of SINR: Fig. 8 illustrates the spatial distribution of the received SINR on the transverse plane at the receiver for a 5 × 5 PD array, for w 0 = 50 µm and w 0 = 100 µm, representing two cases for the beam spot radius including w(L) = 10.8 mm and w(L) = 5.4 mm, respectively. For w 0 = 50 µm with a larger beam spot size at the receiver, as shown in Fig. 8(a), the SINR ranges from −6 dB to +12 dB. A dissimilar distribution of the SINR over 2) Rate vs. Beam Waist: Fig. 9 demonstrates the aggregate data rate achieved by the proposed MIMO system under perfect alignment conditions when w 0 is varied from 10 µm to 100 µm.
For all MIMO realizations under consideration, the aggregate rate monotonically increases for larger values of w 0 . At the beginning for w 0 = 10 µm, the beam spot size at the receiver is very large, i.e. w(L) = 54 mm. This renders the signal power collected by each PD from direct links very low. Besides, there is substantial crosstalk among the incident beams, which severely degrades the performance. The use of SVD yields the upper bound performance for the MIMO system. When w 0 increases, the data rate grows, and so does the gap between the performance of the MIMO system without SVD and the upper bound. The maximum difference between the performance of the two systems occurs at about w 0 = 40 µm. After this point, by increasing w 0 , the aforementioned gap is rapidly reduced and the data rate for the MIMO system without SVD asymptotically approaches that with SVD. The right tail of the curves in Fig. 9 indicates Tb/s, 9 × 9, 16 × 16 and 25 × 25 systems, respectively, require the beam waist radii of 98 µm, 60 µm and 50 µm. This target is not achievable by a 4 × 4 system for w 0 ≤ 100 µm.

B. GMM Verification
The GMM of the SISO channel developed in Section III is the underlying foundation for the MIMO misalignment modeling presented in Section IV. Therefore, its accuracy needs to be verified with a dependable benchmark. A powerful commercial optical design software by Zemax, known as OpticStudio [34], is used for this purpose. Empirical data is collected by running extensive simulations based on non-sequential ray tracing in OpticStudio. Fig. 10  This corroborates the high reliability of the proposed analytical modeling framework.

C. Impact of Misalignment
The impact of misalignment is studied for an N t × N r MIMO system with N t = 25 (i.e. 5 × 5 VCSEL array) and w 0 = 100 µm, as this is the case where the maximum data rate can be achieved, though such very low divergence beams are prone to misalignment errors. To make use of the SVD architecture, three configurations of the PD array with different number of elements as shown in Fig. 11 are examined, where N r = 25, 41, 81 in Configs. I, II and III, respectively.
For the PD array structure described in Section IV-A, the fill factor (FF) can be calculated as The corresponding FF value for Configs. I, II and III is FF = 20%, 32%, 64%. 1) Radial Displacement Error: Fig. 12 plots the aggregate data rate of the system with respect to the radial displacement error. 'Exact' and 'Approx.' labels refer to the exact MIMO GMM procedure from Section IV-B and the approximate formula in (35) used for computing the MIMO channel gains for the case without SVD. The results of the two approaches perfectly match with each other for both scenarios shown in Figs. 12a and 12b. In Fig. 12(a), the radial displacement takes place along one of the two lateral dimensions of the PD array relative to the VCSEL array, i.e. horizontally. Note that vertical displacement is similar to this case due to the symmetry in the square lattice of the arrays. It is observed from Fig. 12(a) that without SVD, the data rate rapidly drops to zero after r DE = 6 mm. When using SVD for Config. I, the rate performance exhibits a decaying oscillation with exactly five peaks, each of which corresponds to an event that a given number of columns of the VCSEL array are aligned to those of the PD array. The use of Config. II slightly improves the valleys because of the second lattice of PDs placed on the vertices of the first lattice. This keeps the aggregate data rate above 1 Tb/s for r DE ≤ 17 Config. II is more pronounced. Also, Config. III consistently retains the highest performance.
2) Transmitter Orientation Error: Fig. 13 presents the aggregate data rate of the system as a function of the orientation angle error at the transmitter. For the system without SVD, the results evaluated by using the approximate expression in (40) perfectly match with those obtained based on the MIMO GMM from Section IV-B for both cases shown in Figs. 12a and 12b. It is evident how sensitive the system performance is with respect to the transmitter orientation error such that a small error of about 1.7 • is enough to make the data rate zero. This is because the transmitter is 2 m away from the receiver, and hence small deviations in its orientation angle are translated into large displacements of the beam spots on the other end. In Fig. 13(a), the orientation error happens only in the azimuth angle φ a by assuming φ e = 0. The results in Fig. 13(a) have similar trends as those in Fig. 12(a), except for their different scales in the horizontal axis. In fact, an azimuth angle error of 1.7 • is equivalent to a horizontal displacement error of 60 mm.
Consequently, φ a = 1.7 • causes the beam spots of the VCSELs to completely miss the receiver array. Likewise, Fig. 13(b) can be paired with Fig. 12(b). Fig. 13(b) represents the case where the  orientation error comes about in both azimuth and elevation angles equally, which produces the same effect as the diagonal displacement error as shown in Fig. 12(b). Therefore, the transmitter orientation error can be viewed as an equivalent radial displacement error, if the beam spot size at the receiver array is sufficiently small, as formally established in Section IV-C. Fig. 14 shows the aggregate data rate when the orientation angle error at the receiver is variable. It can be clearly seen that the MIMO system is significantly more tolerant against the receiver misalignment as compared to the transmitter misalignment in terms of the orientation angle error. In Fig. 14(a), the azimuth angle ψ a is varied between 0 and 90 • while the elevation angle is fixed at ψ e = 0. It is observed that even without SVD, the rate is above 1 Tb/s over a wide range of ψ a , i.e. for 0 ≤ ψ a ≤ 46 • . The use of SVD gives an almost equal performance for Configs. I and II, providing a noticeable improvement with respect to the case without SVD by maintaining the rate above 1 Tb/s for 0 ≤ ψ a ≤ 65 • . Since the size of the PDs is smaller than the size of the beam spots, small rotations of the PD array about its axes have a marginal effect on the system performance, unless they are sufficiently large to alter the distribution of the received optical power on the PD array. Also, the performance of Config. III is slightly better than Configs. I and II. In Fig. 14(b), where ψ e = ψ a , the performance without  [35]. Thus, the computational complexity for Configs. I, II and III, respectively, is O(4 × 10 5 ), O(5 × 10 5 ) and O(10 6 ). Therefore, Config. III requires 2.5 times higher computational power than Config. I. To elaborate on the performance gain, in the case of radial displacement as shown in Fig. 12, the first crossover of the aggregate rate from the 1 Tb/s threshold takes place at horizontal displacements r DE = 5.3, 17.2, 38.8 mm for Configs. I, II and III, respectively. Compared with Config. I, Config. III provides more than 7 times wider room for misalignment tolerance while incurring 2.5 times higher complexity. This underlines the point that it is worthwhile to use Config. III for the system implementation.

VI. CONCLUSIONS
A VCSEL-based MIMO OWC system using DCO-OFDM and spatial multiplexing techniques was designed and elaborated. The fundamental problem of the link misalignment was supported by extensive analytical modeling and the generalized model of misalignment errors was derived for both SISO and MIMO link configurations. Under perfect alignment conditions, data rates of ≥ 1.016 Tb/s are achievable with a 25 × 25 MIMO system without SVD over a link distance of 2 m for a beam waist radius of w 0 ≥ 50 µm, equivalent to a beam divergence angle of θ ≤ 0.3 • , while fulfilling the eye safety constraint. The same system setup attains ≥ 1.264 Tb/s data rates when SVD is applied. The use of w 0 = 100 µm (i.e. θ = 0.16 • ) renders beam spots on the receiver array almost nonoverlapping and the aforementioned 25 × 25 system delivers 2.835 Tb/s with or without SVD. In fact, when the link is perfectly aligned, the use of SVD is essentially effective in the crosstalk-limited regime. The derived GMM was used to study the effect of different misalignment errors on the system performance. Under radial displacement error or orientation angle error at the transmitter, the performance of MIMO systems with SVD shows a declining oscillation behavior with increase in the error value. For a 25×25 system using w 0 = 100 µm, the aggregate rate stays above the 1 Tb/s level for horizontal displacements of up to r DE = 17 mm (0.28 relative to the array side length). The performance remains over 1 Tb/s for an orientation error of φ ≤ 0.8 • in the azimuth angle of the transmitter. In the presence of the receiver orientation error, the aggregate rate is maintained above 1 Tb/s over a wide range of ψ ≤ 65 • for the azimuth angle of the receiver. The results indicate that the orientation angle error at the transmitter is the most impactful type of misalignment. They also confirm that the impact of misalignment is alleviated by using a receiver array with densely packed PD elements, improving the system tolerance against misalignment errors. This is especially pronounced for the radial displacement and orientation angle error at the transmitter. Future research involves extended system modeling and performance evaluation under practical design limitations including multimode output profile of VCSELs, frequency-selective modulation response of VCSELs and PDs, and receiver optics. An interesting application of the proposed Tb/s OWC backhaul system is wireless inter-rack communications in high-speed data center networks, which brings an avenue for future research. APPENDIX In this appendix, the aim is to verify the accuracy of the approximations applied to the MIMO GMM in Section IV-C to derive (35) and (40). For clarity, the SISO channel gain between VCSEL j and PD i is evaluated for x i = y i = 0. Due to the circular symmetry of the beam, it is assumed that the radial displacement is located along the x ′ axis on the receiver plane such that    ratio w(L) r PD . It can be seen that the proposed approximations excellently conform with the exact values especially when w(L) r PD > 1. In the case where the beam spot size is equal to the PD size, the approximate expressions slightly underestimate the exact values in a logarithmic scale. As a measure of accuracy, normalized mean square error (NMSE) between the exact and approximate expressions are calculated in Table II. By increasing w(L) r PD , NMSE rapidly decreases for both cases and NMSE < 10 −3 always holds. Hence, the underlying approximations are highly accurate for w(L) r PD ≥ 1. Note that this conclusion is independent of the link distance as the results are presented in terms of the beam spot size at the receiver. The smallest beam spot size considered in Section V is w(L) = 5.4 mm for w 0 = 100 µm, in which case w(L) r PD = 1.8 and NMSE = 7.9055 × 10 −5 .