Using Loaded N-port Structures
to Achieve the Continuous-Space
Electromagnetic Channel Capacity Bound

We firstly constrain the radiated power for the derivation of channel capacity. This optimization problem can be formulated as

By decomposing arbitrary radiation pattern $\mathbf{e}$ onto the set of orthonormal radiation pattern basis in $\mathbf{B}_{T}$, we have

with

so that each entry in $\boldsymbol{{\beta}}$ refers to the magnitude allocated to each orthonormal pattern basis. The radiated power constraint (12) can be transformed to

where $\mathbf{R}_{\beta}=\mathrm{E}\left\{\boldsymbol{{\beta}}\boldsymbol{{\beta}}^{H}\right\}$ is the covariance matrix of orthonormal pattern basis magnitude. Accordingly, the optimization problem (11) and (12) is re-formulated as

It can be noticed that the problem (16) and (17) is the same as the capacity maximization problem of ideal $M\times N$ MIMO. However, due to some extremely small eigenvalues in $\Lambda$, the norm of current $\mathbf{i}$ could be extremely large and not implementable. In other words, the maximum current is unconstrained.

Next, we consider imposing a constraint on the current norm $\mathrm{Tr}\left(\mathbf{R}_{\mathbf{i}}\right)$. The optimization problem of channel capacity can then be formulated as

By decomposing the current $\mathbf{i}$ onto the set of current basis $\mathbf{q}_{i}$, $i=1,2,...,N$, we can obtain the coefficient for the decomposition as $\gamma_{i}=\left\langle\mathbf{q}_{i},\mathbf{i}\right\rangle=\mathbf{q}_{i}^{H}\mathbf{i}$ which refers to the magnitude allocated to the $i$th current basis. We collect the coefficient $\gamma_{i}$ into a vector $\boldsymbol{{\gamma}}=\left[\gamma_{1},\gamma_{2},...,\gamma_{N}\right]^{T}$ which can be related to $\mathbf{i}$ by

$\boldsymbol{{\gamma}}$ has the same norm as $\mathbf{i}$ so that the current constraint (19) can be transformed to

where $\mathbf{R}_{\gamma}=\mathrm{E}\left\{\boldsymbol{{\gamma}}\boldsymbol{{\gamma}}^{H}\right\}$ is the covariance matrix of current basis magnitude. Then the problem (18) and (19) can be transformed to

The optimal solution for the currents can be found by performing equal power (EP) or water-filling (WF) allocation directly on $\mathbf{R}_{\gamma}$.

It can be observed that due to the existence of $\Lambda^{\frac{1}{2}}$ in (22), the radiated power of each orthogonal basis in $\mathbf{B}_{T}$ is not equal and proportional to the corresponding eigenvalues in $\Lambda$. Therefore, in this formulation, the WF method tends to allocate more power to those basis with larger eigenvalues since it increases radiated power and resulting system capacity. In other words, the radiated power is unconstrained.

Finally, we consider dual constraints imposed by the current norm and radiated power simultaneously. That is, in practical setups the radiated power is limited and the norm of current $\mathbf{i}$ must also be restricted. This optimization problem is formulated as

We follow the transformation in (14) as well as the formulation in (16) and (17), so the problem (24) to (26) can be transformed to

It can be observed that the capacity formulation is only related to $\Lambda$, which is a diagonal matrix with diagonal entries being the eigenvalues of $\mathbf{K}_{T}$.

To solve the problem (27) to (29), we firstly consider EP allocation method. Although the system provides $N$ orthogonal radiation pattern basis, some eigenvalues in $\Lambda$ are too small so that large currents are required to radiate their corresponding basis (5) which cannot be used for practical MIMO antenna design. To avoid generating large currents, i.e. satisfying the current constraint (29), $P_{\mathrm{rad}}$ should be allocated to those basis with largest radiation resistance (i.e. larger eigenvalues in $\Lambda$). Therefore, we equally allocate power $P_{\mathrm{rad}}$ to the first $N_{\mathrm{eff}}$ $\left(N_{\mathrm{eff}}\leqslant N\right)$ basis and drop the remaining $N-N_{\mathrm{eff}}$ basis. These $N_{\mathrm{eff}}$ basis can be effectively used with maximum current norm $I_{\mathrm{in}}^{2}$ and radiated power $P_{\mathrm{rad}}$, and thus can be regarded as the effective aerial degrees-of-freedom (EADoF) of the transmitter [28].

The diagonal entries in $\mathbf{R}_{\beta}$ are then given by

The constraint (29) can then be written as

We define $\epsilon=\frac{I_{\mathrm{in}}^{2}}{P_{\mathrm{rad}}}$ in (31) and it can be interpreted loosely as a conductance which has the unit of $\Omega^{-1}$. In essence for a given $P_{\mathrm{rad}}$, a larger $I_{\mathrm{in}}$ ($\epsilon$ high) indicates that the antenna will have a smaller input resistance overall and if it is too small the antenna will not be implementable. Alternatively if $I_{\mathrm{in}}$ is too low ($\epsilon$ low), the input impedance will be required to be too high and also not implementable. Therefore we should set it to be in a range centered around 1/50 $\Omega^{-1}$ so that it is in line with the antennas desired input impedance.

Then the EADoF $N_{\mathrm{eff}}$ is related to $\epsilon$ for a fixed SNR as

where the left term refers to the average of the reciprocal of eigenvalues for the first $N_{\mathrm{eff}}$ basis. It can be observed from (32) that a larger $\epsilon$ indicates that the number of available orthogonal basis, i.e. $N_{\mathrm{eff}}$, and the consequent capacity bound can be larger. The maximum $N_{\mathrm{eff}}$ can be obtained by solving (32) for a fixed $\epsilon$. That is EADoF $N_{\mathrm{eff}}$ is dependent on $\epsilon$ (the maximum allowable current norm given $P_{\mathrm{rad}}$), indicating the number of basis that can be used with the constraint of current norm and radiation power.

Next we consider using the WF method to maximize capacity in problem (27) to (29). We start with the channel gain of $N$ sub-channels in channel matrix $\mathbf{H}_{\mathrm{iid}}$. By performing EVD on $\mathbf{H}_{\mathrm{iid}}\mathbf{H}_{\mathrm{iid}}^{H}$, we obtain the channel gain matrix $\mathbf{S}=\mathrm{diag}\left(s_{1},s_{2},...,s_{N}\right)$ which is a real diagonal matrix with $s_{i}$, $i=1,2,...,N$, being the channel gain of the $i$th sub-channel. Then we use the Lagrangian method whose function is given by

where $\mu_{\mathrm{rad}}$ and $\mu_{\mathrm{in}}$ are Lagrangian multipliers. Taking the partial derivative of $L$ with respect to $\mathbf{R}_{\beta}$, we have

Combining the solution of $\mathbf{R}_{\beta}$ in (34), two constraints (28), (29) can be written as [21]

The optimal multipliers $\mu_{\mathrm{rad}}^{\star}$ and $\mu_{\mathrm{in}}^{\star}$ in (36) and (35) can be found via binary search [21] and the optimal $\mathbf{R}_{\beta}$ can be solved in (34).

It should be noted that the above closed-form expressions for system capacity are formulated without considering the physical realization and practical excitation for MIMO antennas. Therefore, the MIMO antenna structure should be carefully designed and optimized to achieve the capacity bounds with dual constraints, which will be described in the next section.

To begin to construct the MIMO antenna, we divide the $N$-ports in Fig. 1(b) into $Q$ active feeding ports and $N-Q$ loaded ports. The advantage of this approach is that the feeds for the final potential $Q$-port MIMO antenna are incorporated at the beginning of the analysis through the $N$-port structure. The equivalent circuit model for the $N$-port network with $Q$ active feeding ports (numbered 1 to $Q$) and $N-Q$ parasitic loaded ports (numbered $Q+1$ to $N$) is shown in Fig. 2. The loaded ports are either lossless with inductive and capacitive reactance or open. The active feeding ports are connected to a matching network to achieve the necessary 50 $\Omega$ input impedance. By optimizing the $N-Q$ load reactances, we wish to generate $Q$ orthogonal radiation patterns from the designated $Q$ feeding ports so that the capacity (8) can be maximized. In essence we will attempt to find loads that create the required current distributions for approaching the capacity bounds. The design procedure for selecting which ports are fed and which are loaded to achieve the optimum channel capacity using (10) is described later.

In Fig. 2, the $q$th feeding port is excited by a voltage source $v_{q}$, $q=1,2,...,Q$, and the $\left(Q+p\right)$th parasitic port is loaded with reactance $x_{Q+p}^{L}$, $p=1,2,...,N-Q$. We group the voltage source excitation $v_{q}$, $q=1,2,...,Q$, into a vector as $\mathbf{v}_{\mathrm{A}}=\left[v_{1},v_{2},...,v_{Q}\right]{}^{T}\in\mathbb{C}^{Q\times 1}$ and define $\mathbf{v}_{\mathrm{P}}=\left[0,0,...,0\right]{}^{T}\in\mathbb{C}^{\left(N-Q\right)\times 1}$ as a $\left(N-Q\right)$-dimension zero vector due to there being no excitation at the parasitic ports. We also group current at the feeding and parasitic ports into vectors as $\mathbf{i}_{\mathrm{A}}=\left[i_{1},i_{2},...,i_{Q}\right]{}^{T}\in\mathbb{C}^{Q\times 1}$ and $\mathbf{i}_{\mathrm{P}}=\left[i_{Q+1},i_{Q+2},...,i_{N}\right]{}^{T}\in\mathbb{C}^{\left(N-Q\right)\times 1}$, respectively. The voltage and current in the $N$-port network are then related by

where $\mathrm{\mathbf{X}}^{L}=\mathrm{diag}\mathnormal{\left(x_{Q+\text{1}}^{L},x_{Q+\text{2}}^{L},\cdots,x_{N}^{L}\right)}$ represents the load reactance connected to each parasitic port. $\mathbf{Z}_{\mathrm{A}}\in\mathbb{C}^{Q\times Q}$ and $\mathbf{Z}_{\mathrm{P}}\in\mathbb{C}^{\left(N-Q\right)\times\left(N-Q\right)}$ are the impedance sub-matrix of the $Q$ feeding ports and the $N-Q$ parasitic ports, respectively. $\mathbf{Z}_{\mathrm{AP}}\in\mathbb{C}^{Q\times\left(N-Q\right)}$ and $\mathbf{Z}_{\mathrm{PA}}\in\mathbb{C}^{\left(N-Q\right)\times Q}$ are the sub-matrix referring to the mutual impedance between $Q$ feeding and $N-Q$ parasitic ports with $\mathbf{Z}_{\mathrm{AP}}=\mathbf{Z}_{\mathrm{PA}}^{T}$. From (37) we can obtain the relationship between the current on the feeding and parasitic ports as

Our interest is obtaining $Q$ orthogonal far-field radiation patterns to form the beamspace MIMO system. To that end, we collect the radiation pattern of each feeding port $\mathrm{\mathbf{e}}_{T,q}\left(\Omega\right)$, $q=1,2,...,Q$, in the antenna system into matrix form $\mathrm{\mathbf{E}}_{T}=\left[\mathrm{\mathbf{e}}_{T,1}\left(\Omega\right),\mathrm{\mathbf{e}}_{T,2}\left(\Omega\right),...,\mathrm{\mathbf{e}}_{T,Q}\left(\Omega\right)\right]$, which can be given as

where $\Omega=\left(\theta,\phi\right)$ denotes the spatial angle with $\theta$ and $\phi$ representing the elevation and azimuth angles in spherical coordinates, respectively. $\mathbf{E}_{\mathrm{A}}\left(\Omega\right)=\left[\mathbf{e}_{\mathrm{A},1}\left(\Omega\right),\ldots,\mathbf{e}_{\mathrm{A},Q}\left(\Omega\right)\right]$ collects $Q$ open-circuit radiation patterns of the feeding ports with $\mathbf{e}_{\mathrm{A},q}\left(\Omega\right)$, $q=1,2,...,Q$ being the open-circuit radiation pattern of the $q$th feeding port excited by a unit current when all the other feeding and parasitic ports are open-circuit. Similarly, $\mathbf{E}_{\mathrm{P}}\left(\Omega\right)=\left[\mathbf{e}_{\mathrm{P},Q+1}\left(\Omega\right),\ldots,\mathbf{e}_{\mathrm{P},N}\left(\Omega\right)\right]$ collects $N-Q$ open-circuit radiation patterns of the parasitic ports with $\mathbf{e}_{\mathrm{P},Q+p}\left(\Omega\right)$, $p=1,2,...,N-Q$, being the open-circuit radiation pattern of the $\left(Q+p\right)$th parasitic port excited by a unit current when all the other feeding and parasitic ports are open-circuited. It can be observed that $\mathrm{\mathbf{E}}_{T}$ consists of the original open-circuit radiation pattern of the $Q$ feeding ports $\mathbf{E}_{\mathrm{A}}$ and a perturbation term affected by the load reactance $\mathrm{\mathbf{X}}^{L}$ across the parasitic ports.

We use a full electromagnetic solver, CST studio suite [34], to simulate the open-circuit radiation patterns of all $N$ ports of the loaded $N$-port structure. It should be noted that the simulation only needs to be performed once because any radiation pattern excited by any current distribution can then be found by using (39) [35]. This reduces the computation complexity enormously as full electromagnetic simulation is not needed during the load optimization and capacity simulation stages.

To obtain the optimal design of the MIMO antenna, we wish to select $Q$ active feeding ports and optimize the load reactances at the $N-Q$ parasitic ports. The objective is to generate a set of orthogonal radiation patterns in $\mathrm{\mathbf{E}}_{T}$ at $Q$ feeding ports so that the resulting channel capacity can be maximized.

To achieve an optimization result that also meets practical design constraints, we need to note that not all ports would be suitable for use as feeding. For example, ports positioned on the edges of the structure are more likely to be accessible for feeds than those completely surrounded by ports. For this reason, we partition the $N$ ports into two sets for optimization. Assume $S$ $\left(S>Q\right)$ out of $N$ ports are feasible feeding ports and whose indices in the set are given by $\mathcal{S}=\left\{1,2,...,S\right\}$. The objective is therefore to find the optimum $Q$ ports from the $S$ feasible ports where the set of indices for the $Q$ feeding ports is given by $\mathcal{P}=\left\{p_{1},p_{2},...,p_{Q}\right\}\subset\mathcal{S}$ . In addition, we must find the optimum loads for the remaining $N-S$ ports so that together with the $S-Q$ unselected feeding ports, the load reactance matrix $\mathbf{X}^{L}$ can meet the requirements of orthogonality for the $Q$ selected feeding ports.

We use the correlation coefficient between the radiation patterns of two feeding ports, i.e. $\text{$\mathrm{\mathbf{e}}_{T,j}$}\left(\Omega,\mathrm{\mathbf{X}}^{L}\right)$ of the $j$th feeding port and $\text{$\mathrm{\mathbf{e}}_{T,k}$}\left(\Omega,\mathrm{\mathbf{X}}^{L}\right)$ of the $k$th feeding port, to evaluate their similarity, and this is defined as

where $\rho_{jk}$ satisfies $0\leq\left|\rho_{jk}\right|\leq 1$. When $\left|\rho_{jk}\right|=0$, $\mathrm{\mathbf{e}}_{T,j}$ and $\mathrm{\mathbf{e}}_{T,k}$ are orthogonal to each other. Leveraging the correlation coefficient (40), we can formulate the optimization problem as

It should be noted that the optimization variables $\mathcal{P}$ and $\mathbf{X}^{L}$ are highly coupled with each other in this problem. This is because altering feeding port indices greatly changes the impedance matrix $\mathbf{Z}_{\mathrm{A}}$, $\mathbf{Z}_{\mathrm{AP}}$ and $\mathbf{Z}_{\mathrm{P}}$ in (37), making the problem (41) complicated to solve. To overcome this difficulty and meet these constraints, we propose the use of the alternating optimization method to iteratively optimize the load reactances and selection of the feeding ports [36], [37]. For simplicity, the $S-Q$ unselected feeding ports are left open in the entire optimization process to ease the computational burden. However they could also be included in the load optimization process if deemed necessary.

To start with, we randomly select $Q$ feeding ports, out of $S$, as the initial guess at iteration 0. The set consisting of the indices for $Q$ initial random feeding ports is given by $\mathcal{P}^{\left(0\right)}=\left\{p_{1}^{\left(0\right)},p_{2}^{\left(0\right)},...,p_{Q}^{\text{$\left(0\right)$}}\right\}\subset\mathcal{S}$ with $p_{q}^{\left(0\right)}$ being the index of the $q$th initial feeding port. Presetting the $S-Q$ unused feeding ports to be open, we then need to optimize $N-S$ load reactances at the parasitic ports where the initial guess for the load reactance matrix is given as $\mathbf{X}^{L\left(0\right)}=\mathrm{diag}\mathnormal{\left(\infty,...,\infty,\text{0},...,\text{0}\right)}$ with the first $S-Q$ diagonal entries being fixed as infinity (open).

At the $i$th iteration, we wish to find the load reactances $\mathbf{X}^{L\left(i\right)}$ to produce near orthogonal radiation patterns with feeding ports indices set $\mathcal{P}^{\left(i-1\right)}$. This can be performed by solving the optimization problem

This makes the correlation coefficient norm between any two patterns in $\mathrm{\mathbf{E}}_{T}$ as close to zero as possible. To solve the unconstrained optimization problem (42), we can use the quasi-Newton method [38] which guarantees convergence to a stationary point of the problem (42). The optimal load reactances at the $i$th iteration are denoted by $\mathbf{X}^{L\left(i\right)}$.

We then use $\mathbf{X}^{L\left(i\right)}$ from (42) to optimize the indices of the feeding ports $\mathcal{P}^{\left(i\right)}$, that is,

which is an NP-hard optimization problem. To solve this problem, we propose a low-complexity approach which is based on sequentially optimizing each entry in $\mathcal{P}^{\left(i\right)}$, i.e. $p_{1}^{\left(i\right)}$, $p_{2}^{\left(i\right)}$,…, $p_{Q}^{\left(i\right)}$, one-by-one. Specifically, we define $\mathcal{P}_{q}^{\left(i\right)}=\left\{p_{1}^{\left(i\right)},...,p_{q}^{\text{$\left(i\right)$}},p_{q+1}^{\left(i-1\right)},p_{q+2}^{\left(i-1\right)},...,p_{Q}^{\left(i-1\right)}\right\}\subset\mathcal{S}$ as the intermediate state of feeding port indices at the $i$th iteration where the indices of the first $q$ feeds have been optimized and the remaining $Q-q$ feeds not optimized. When optimizing index $p_{q}^{\left(i\right)}$ for the $q$th feeding port, we fix the other $Q-1$ indices and select the optimal $p_{q}^{\left(i\right)}$ from the remaining $S-Q+1$ indices in $\mathcal{S}$ to minimize the correlation coefficient (40) among all $Q$ ports in $\mathcal{P}_{q}^{\left(i\right)}$, which can be formulated as

with $q$ sequentially taken as $q=1,2,...,Q$. It should be noted that $\mathcal{P}_{q-1}^{\left(i\right)}=\mathcal{P}^{\left(i-1\right)}$ when $q=1$.

By iteratively optimizing the load reactances and selecting feeding ports, the objective function (42) and (44) can converge to a local optimal solution [37]. Therefore, we obtain the optimal load reactances $\mathbf{X}^{L\star}$ of the loaded $N$-port structures with $Q$ feeds with indices of $\mathcal{P}^{\star}$. Algorithm 1 summarizes the overall algorithm for optimizing the load reactances and feeding ports. The performance of the algorithm in finding optimum MIMO antenna configuration will be shown in the next section.

In the first set of simulation results we wish to verify our approach by comparing with previous results [11]. The previous results consider a continuous source current on a 2D square planar surface with one square wavelength size at both the transmitter and receiver. The 2D surface is discretized into an $N$-port structure as shown in Fig 3 where discretizations with 4 and 180 ports are shown. This corresponds to dividing the one wavelength square surface into 4 and 100 sub-elements respectively similarly to the discretizations of previous work [11].

To compare our method to previous work [11], we only need to find the required currents at all $N$ ports and do not need to consider finding $N-Q$ loads or $Q$ feeding ports. We can therefore use the results in (22) and (23) directly without considering the loads and feeds to find the port currents.

The same channel setup utilized previously [11] is also invoked and includes single polarization in a rich scattering environment with Rayleigh fading and 2D uniform power angular spectrum (PAS) on the azimuth plane [31]. In addition, we use the current norm constraint so that the capacity optimization problem is formulated as (22) and (23), with WF to find the currents. This is similar to equation (21) in [11].

Simulation results of the system capacity are shown in Fig. 4 for various levels of discretizations. In the figure the number of discretizations per dimension is used to align with previous work [11]. That is 10 sub-elements per wavelength correspond to 100 sub-elements in total. It can be observed that by using our method, the system capacity approaches that using the method in [11] when the number of sub-elements per wavelength is more than four. The capacity gap (for 4 or less sub-elements per square wavelength) arises due to the current on the discrete ports being an approximation to the continuous current on the surface, which is related to the position of the ports in the $N$-port structure. Therefore, we can conclude that our method can obtain the continuous-space electromagnetic channel capacity bound if the dimension of the sub-element and the separation distance between ports is less than a wavelength (from Fig. 3 greater than 5 sub-elements per wavelength).

Next, we use the technique in Section IV to design a practical MIMO antenna that approaches the capacity bounds. The frequency considered in the following is 2.4 GHz so that the wavelength is 125 mm.

The first step is to propose a structure that is implementable with $Q$ feeds but general enough to allow the formation of arbitrary current distributions for capacity maximization. As such our proposed $N$-port structure is shown in Fig. 5 where a 2D square copper planar surface with one square wavelength size is again utilized. The copper surface is mounted on a substrate where the copper has electric conductivity of $5.8\times 10^{7}$ S/m while the substrate is made with Rogers 5880C which has permittivity of 2.2 and loss tangent of 0.0009. We discretize the 2D copper surface into a $11\times 11$ sub-element array where the size of the sub-element is $10\times 10\>\mathrm{mm}^{2}$. To implement the feeds, an additional copper plane is placed underneath as a ground plane so that a conventional feed can be connected between the ground and the center of a sub-element as shown in the elevation view in Fig. 5. The $N$ ports of the structure are defined as the ports across each pair of adjacent sub-elements on the surface and also between the ground and the center of each sub-element. That is there are 220 ports on the surface and 121 ports between the ground and sub-elements making up a total of $N=341$ ports. The $Q$ feeding ports need to be selected from the 121 ports between the ground and sub-elements.

While the ultimate aim is to select $Q$ feeding ports from the ground to sub-element ports, we firstly need to find the required currents on the $N$ ports to optimize capacity and approach the capacity bounds similarly as in Section V.A. This allows us to first determine the capacity performance of the structure with ideal port currents. The final load reactances and feeds arrangements of the antenna that can approximate these currents are found in the next subsection.

In the capacity simulations, we obtain them as a function of $P$ ($P=1,2,3$, etc.) basis with the highest eigenvalues (out of $N=341)$ for transmission. This allows us to determine the number of EADoF provided by the structure. In the simulations, we use $P$ ideally isolated antennas at the receiver side. We also assume a rich scattering environment so that entries in $\mathbf{H}_{v}$ satisfy $\left[\mathbf{H}_{v}\right]_{ij}\sim\mathcal{CN}\left(0,1\right)$ for $i,j=1,2,\ldots,K$. In addition, dual polarization in a rich scattering environment with Rayleigh fading and 3D uniform PAS over full sphere is utilized. For these reasons the capacity results obtained for this channel will be much greater than that obtained in Section V.A which used only one polarization, a 2D PAS and a receiver antenna that was the same as the transmitter.

To provide a useful antenna, we need to select feeding ports and the loads for all the ports in the $N$-port loaded structure. The feeding ports and loads must be selected to achieve a close match to the optimum currents found in the previous section for a feasible $P$. This then would provide the optimum antenna configuration that can achieve the channel capacity bound with dual constraints.

As discussed previously for an ideal lossless antenna, the input impedance is $Z_{\mathrm{in}}=50\ \Omega$ and therefore $\epsilon\geqslant I_{\mathrm{in}}^{2}/P_{\mathrm{rad}}=1/Z_{\mathrm{in}}=0.02$ so that $N_{\mathrm{eff}}$ is at least 18. Expecting to have higher currents internally on the antenna, we aim for 20 ports in our design. That is we will need to re-create the currents obtained when $P=20$ from the previous section.

To reduce the computational effort in finding the loads, we preset the sub-elements in the center row and center column of the 2D copper surface as being connected together (shorted) so there is a single large cross element centered on the surface. That is ports {56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 6, 17, 28, 39, 50, 72, 83, 94, 105, 116} in Fig. 5 are shorted together. The cross element separates the other sub-elements into four sub-square arrays. The cross is also isolated from the other sub-elements so that the ports between the cross and the other adjacent sub-elements are open without loads. In effect, we are presetting some of the ports to be shorted or open on the surface. This reduces the computational load of the optimization by reducing the number of load reactances to be found from 220 to 160. In addition, to reduce the search space for the $Q=20$ feeds we also restrict it to the $S=64$ ports between the ground plane and the sub-elements on the edge of each sub-square array (instead of all 121 possible locations). For those ports underneath the surface that were not selected as feeds, we set them to be open.

To find the optimum reactances across the 160 ports on the surface and the indices of the 20 feeding ports underneath, we use the method as described in Section IV so that the excited radiation patterns of the 20 feeding ports are closest to being orthogonal (42). The final indices of the 20 feeding ports after the alternating optimization are given by

and their physical location can be found from Fig. 5. The 20 feeding ports are then each fed by a feeding probe across the common ground and the center of the selected sub-elements as shown in the elevation view in Fig. 5. The impedance matching for the 20 feeding ports is performed separately with a conventional T-type matching network to achieve matching to 50 $\Omega$. The return loss (S11) results of the 20 ports are provided in Fig. 9 where the reflected power at 2.4 GHz is around -15 dB and the bandwidth is around 40 MHz.

By using the 20 excited radiation patterns from the resultant loaded structure with 20 feeds, and a beamspace channel model (7), the capacity of the $20\times 20$ MIMO system is found and the simulated results are shown in Fig. 10. These are also benchmarked with the capacity bound of the loaded structure without feeds obtained by WF method in Section III.C and the capacity of ideal MIMO. It can be observed that the MIMO antenna using the loaded structure with 20 feeds achieves performance close to ideal MIMO and the capacity bound predicted by Fig. 8 for loaded structure without feeds. The capacity gap is due to the small correlations among the radiation patterns of the feeding ports and power loss from mutual coupling. However, the key advantage of the proposed method is that the optimum antenna structure for the MIMO transmitter with size constraint can be found with only minor compromise on system capacity.