—As Multiple Inputs Multiple Outputs (MIMO) is becoming one of the enable techniques in modern wireless communication like 5G/6G and beyond, it is important to design efﬁcient controls for the MIMO antenna ar- rays to realize critical functions such as high-throughput communication and beamforming. Efﬁcient and rigorous Computational Electromagnetics (CEM) algorithms are key for such control design problem, especially for massive MIMO antenna arrays and beyond. Here we present such universal Fast Fourier Transform (FFT) based iterative CEM algorithm for efﬁcient control design of massive MIMO antenna array and beyond, for example Recon- ﬁgurable Intelligent Surface and Large Intelligent Surface (RIS/LIS) MIMO. Our FFT-based iterative CEM algorithm is universal and works for both discrete MIMO and quasi- *Shaolin Liao is the corresponding author and contributes equally as the ﬁrst author. Shaolin Liao is partially supported by start-up grant of SYSU’s “100 Top Talents” Program under its “Academic Leaders” Category. Lu Ou is partially supported by the National Natural Science Foundation of China under Grant No.62002112, China Postdoctoral Science Foundation No. 2020M672488, Science and Technology Projects of Changsha City (kq2004027).continuous RIS/LIS MIMO under various realistic ampli- tude and phase control constraints. It makes use of the translation invariant property of the dyadic function of the MIMO antenna arrays and compute the convolution operator in the spectral domain. Then it updates the control parameters iteratively under the control amplitude and control phase constraints, with the incident ﬁeld and the target output ﬁeld as inputs. Due to the use of FFT, the computational effort of the algorithm scales as NlogN forN=NxNyMIMO antennas, compared to N3of the direct solving method of the linear equation. Numerical simulation for 6G beamforming of both discrete MIMO and RIS/LIS MIMO under various control design constraints has been carried out to show the efﬁciency of the algorithm. Index Terms —Computational Electromagnetics (CEM), massive MIMO, Reconﬁgurable Intelligent Surface (RIS), Large Intelligent Surface (LIS), 5G, millimeter wave/THz 6G. I. I NTRODUCTION Millimeter Wave (mmW: 30 GHz -300 GHz) and Terahertz (THz: 300 GHz -3 THz) boasts of a broad 1 spectrum and has been widely considered for one of the next-generation 6G wireless communication techniques to resolve the ever increasing demand on data capacity and speed. Due to its smaller wavelength, mmW/THz 6G wireless communication allows for easy deployment of massive antenna arrays for Multiple Input Multiple Output (MIMO) wireless communication [1], [2], for example, Reconﬁgurable Intelligent Surface [3]–[8] and Large Intelligent Surface [9]–[11] (RIS/LIS) MIMO and even holographic MIMO [12], to realize many important functions such as beamforming [13], massive MIMO (mMIMO) wireless communication [14] and even Com- pressed Sensing (CS [15]–[19]) that makes use of the spatial diversity for higher data rate. Thus MIMO control design becomes an important task to achieve the optimal performance of mmW/THz 6G wireless communication: to design the optimal control pa- rameters for the MIMO antenna arrays, with the incident wave and the target output ﬁelds as the inputs. However, as the unknown control parameters scale as the MIMO antenna array size N, the direct solving method for the linear equation system requires the inverse of a NNma- trix, which means a computation effort of N3. Such cubic scaling of the computation effort is unbearable due to the ever-increasing MIMO antenna array size N: from the discrete MIMO, to massive MIMO and quasi-continuous sub-wavelength RIS/LIS MIMO. Another issue that arises with the MIMO antenna arrays is the realistic constraints on control amplitude and phase of the MIMO antenna arrays, for which usually analytical optimized solution is not possible and numerical algorithm is required. All of the above call for efﬁcient and scalable Computational Electromagnetics (CEM) design algorithm. In this paper, we present such efﬁcient CEM control design algorithm based on the Fast Fourier Transform (FFT) to accelerate the optimal control design of the MIMO antenna arrays. The FFT makes full use of the translation invariant property of the Green’s function and solve the convolution equation in the spectral domain to realize anNlogNcomputational time, as compared to theN3required by the direct solving method. II. MIMO D ESIGN : THEORETICAL FORMULATION Fig. 1a) shows the problem formulation for the MIMO design: the rectangular MIMO array consists ofN=NxNyantennas that has length periods ofLx;Lyin^x;^ydirection. The unknown parameters are the distributed surface current value J(x;y) = fJx(x;y);Jy(x;y);Jz(x;y)gof theN=NxNyan- tennas. Assuming linear excitation of the MIMO antenna unit cell, the distributed surface current value J(x;y)is proportional to only the excitation ﬁeld for all (x;y)lo- cation. So the unknown parameter J(x;y)can be reduced to the surface current of some anchor point of the MIMO antenna unit cell, for example, the geometric center. Let’sdenote such unknown surface current values at the anchor points of the MIMO array in the vector format as JN. Similarly, the incident ﬁeld upon the MIMO antenna array is denoted as vector Ei Nand the target output ﬁeld is denoted as vector Eo N. Without loss of generality, we assume that the con- trolled surface current values at the anchor points of the controlled MIMO array is given by JN=CNEi N; (1) whereCi Nis the unknown complex control parameter vector that remains to be solved. The output ﬁeld Eo Ncan be calculated from the surface currentINas follows Eo N=L~JN; (2) where Lis the dyadic function and ~denotes the spatial convolution. Thus, the surface current values at the anchor points of the MIMO antenna array JNcan be solved as follows, JN=L1 ~Eo N; (3) where L1 is the inverse dyadic function of Lthat satisﬁes the following property L1 ~L=I; (4) where Iis the unit matrix. Substituting Eq. (3) into Eq. (1), the linear equation of the collocation method is obtained JN=L1 ~Eo N=CNEi N; (5) which is the mathematical formulation of the MIMO control design problem. A. Discrete MIMO Antenna Arrays For discrete MIMO antenna arrays consisting of unit cell of ﬁnite size, the dyadic function Lis actually the radiation pattern of each individual antenna. For example, the dyadic function in the far ﬁeld for plane wave incidence upon a square metallic patch antenna in the far ﬁeld under the Physical Optics (PO) approximation is given as follows [20], Lx/sin cossinc (X)sinc(Y) ; (6) L x/cos cossinc (X)sinc(Y); (7) where the far ﬁeld is given for x-polarized incident ﬁeld in the spherical coordinates (; ). Similar expression can be obtained for y-polarized incident ﬁeld. 2 Fig. 1. MIMO control design problem formulation of an N=NxNyMIMO antenna array with a unit cell area of A=LxLy: a) control design for MIMO beamforming; b) the FFT-based iterative CEM algorithm. B. Quasi-continuous RIS/LIS MIMO For continuous MIMO antenna arrays of inﬁnitesimal small unit cell, the dyadic function Lreduces to the dyadic Green’s function Ge given as follows [21], Ge = 1 +1 k2rr g(r;r0)I; (8) whereg(r;r0)is the scalar Green’s function at the observation point rdue to the source point r0. III. MIMO C ONTROL DESIGN :THE FFT- BASED CEM ALGORITHM One issue to solve the surface current values at the anchor points of the MIMO antenna array JNof Eq. (5) is that, it requires the inverse of the dyadic function L, which means N3computation effort. However, due to the spatial translation invariant of the convolution, FFT and Inverse FFT (IFFT) can be used to accelerate the computation as follow [21], [22] JN=IFFT L1 Eo N ; (9) whereLandCNare the Fourier transformed spectra of the dyadic function L, and the control parameter at the anchor points of the MIMO antenna array respectively. A. Exact Complex Control For MIMO antenna array that can be controlled exactly in both amplitude and phase, the control parameter CN of Eq. (5) can be readily obtained as follows CN=IFFT L1 Eo N Ei N: (10)B. Digital Control Another issue remains to be solved is that, due to real- istic constraints of technique and cost, MIMO amplitude and phase control might be limited. For example, only digital amplitude and digital phase are available. Under such situation, an analytical solution of Eq. (5) is not possible and calls for iterative algorithm that updates the control amplitude and phase at each iteration [22]–[25]. Here we study such iterative CEM algorithm based on FFT for efﬁcient control design of MIMO antenna arrays under various amplitude and phase constraints. As shown in Fig. 1b), the FFT-based iterative CEM algorithm follows the procedure given below, 1) Compute the surface current JNwith the target output ﬁeld Eo Naccording to Eq. (9); 2) Update the control parameter CNaccording to Eq. (1) with the surface current JNand the incident ﬁeldEi N; 3) Apply digital amplitude and phase constraints on the control parameter CN ; 4) Update the constrained surface current JN with the constrained control parameter CN ; 5) Compute output ﬁeld Eo Naccording to Eq. (2) with the constrained surface current JN ; 6) Update the target output ﬁeld Eo Nby keeping its phase while replacing its amplitude pattern by the designed amplitude pattern; 7) Go to Step 1) and repeat till the preset convergence criterion is satisﬁed. C. Phase-only Control When only phase control is available, the iterative CEM algorithm is similar to that for digital control in Section III-B, except applying uniform amplitude constraint. 3 Fig. 2. Sparse beamforming MIMO control design: a) Beam pattern at 3m for exact complex control; b) Beam pattern at 3m for digital control; c) Beam pattern at 3m for digital phase-only control; d) control amplitude of the exact complex control; e) control amplitude of the digital control; f) control phase of the exact complex control; g) digital control phase of the digital control and the phase-only control. IV. MIMO C ONTROL DESIGN : EXAMPLES The quasi-continuous sub-wavelength RIS/LIS MIMO array has attracted increasing attention, especially for the upcoming millimeter-wave/THz 6G wireless com- munication. Here two RIS/LIS MIMO control design examples are given to show how Eq. (10) can be used to efﬁciently solve the unknown control parameter CN of the MIMO antenna array: 1) sparse beamforming; 2) OAM beamforming. In all simulation, a digital amplitude of 0.1 incremental scale or resolution and digital phase of 45incremental scale or resolution have been used. Also, square RIS/LIS MIMO antennas with negligible gaps are used, with a unit cell area of 0:10:1inside which the surface current can be approximately considered as uniform. Finally, the incident ﬁeld is a Gaussian beam. A. Sparse beamforming In the dynamic wireless communication environment, mobile users are moving around. Thus steering commu- nication energy to mobile users can avoid wasting energy and increase the Signal to Noise Ratio (SNR) of the wireless communication, which means higher data rate according to the Shannon’s theorem. For such reason,it is desirable to design optimal sparse beamforming RIS/LIS MIMO control under its amplitude and phase constraints. Here we design such sparse beamforming pattern with 5 smartphone communication locations with equal communication energy. Finally, the Gaussian width of the incident ﬁeld is 25, while that of the target sparse output ﬁeld is 5. 1) Exact Complex Control: Fig. 2a) shows the ob- tained sparse MIMO beamforming pattern at a distance of 3m for exact complex control and Fig. 2d), Fig. 2f) show the corresponding control amplitude and phase respectively. It can be seen that, with exact complex control without any amplitude and phase constraint, the designed beamforming pattern behaves well and all ﬁve beam spots are well identiﬁed with small side lobes are other locations. 2) Digital Control: Similarly, Fig. 2b) shows the ob- tained sparse MIMO beamforming pattern at a distance of 3m for digital control and Fig. 2e), Fig. 2g) show the corresponding control amplitude and phase respectively. It can be seen that, with digital amplitude and digital phase control, the designed beamforming pattern contains bigger side lobes at other locations. 4 Fig. 3. OAM beamforming MIMO control design: a) Beam pattern at 3m for exact complex control; b) Beam pattern at 3m for digital control; c) Beam pattern at 3m for digital phase-only control; d) control amplitude of the exact complex control; e) control amplitude of the digital control; f) control phase of the exact complex control; g) digital control phase of the digital control and the phase-only control. 3) Phase-only Control: At last, Fig. 2c) shows the obtained sparse MIMO beamforming pattern at a distance of 3m for digital phase-only control, from which it can be seen that more and bigger side lobes appear at a much larger area. However, it still roughly forms 5 communication spot locations. B. OAM beamforming OAM is a special type of MIMO communication that has its unique property such as coherent spatial phase that can be used to increase the SNR and thus the wireless communication data rate. Here we design such OAM beamforming pattern with 4 smartphone communication locations with equal communication energy. Finally, the Gaussian width of the incident ﬁeld is 100, while that of the target OAM output ﬁeld is 50. 1) Exact Complex Control: Fig. 3a) shows the ob- tained OAM MIMO beamforming pattern at a distance of 100m for exact complex control and Fig. 3d), Fig. 3f) show the corresponding control amplitude and phase respectively. It can be seen that four OAM communication spot locations clearly show up, as expected.2) Digital Control: Similarly, Fig. 3b) shows the ob- tained OAM MIMO beamforming pattern at a distance of 100m for digital control and Fig. 3e), Fig. 3g) show the corresponding control amplitude and phase respectively. It can be seen that satisﬁed OAM beamforming pattern still can be obtained, though with bigger side lobes spreading into a larger area. 3) Phase-only Control: At last, Fig. 3c) shows the obtained OAM MIMO beamforming pattern at a distance of 100m for digital phase-only control, from which it can be seen that the four OAM spot location overlaps with each other and bigger side lobes spread into a even larger area. V. C ONCLUSION An FFT-based iterative CEM algorithm has been devel- oped for universal control of massive MIMO and beyond. The algorithm is universal and works for both discrete MIMO antenna arrays and quasi-continuous RIS/LIS MIMO antenna arrays. Numerical simulation has been performed for both sparse beamforming and OAM beam- forming, under three control constraints: exact complex control, digital control and phase-only control. The results 5 show that the FFT-based iterative CEM algorithm is efﬁcient in design of arbitrary control amplitude and phase for massive MIMO and beyond, with a computation effort scales as NlogNfor anN=NxNyMIMO antenna array.