Broadband RF Phased Array Design with MEEP: Comparisons to Array Theory in Two and Three Dimensions

Abstract

Phased array radar systems have a wide variety of applications in engineering and physics research. Phased array design usually requires numerical modeling with expensive commercial computational packages. Using the open-source MIT Electrogmagnetic Equation Propagation (MEEP) package, a set of phased array designs is presented. Specifically, one and two-dimensional arrays of Yagi-Uda and horn antennas were modeled in the bandwidth [0.1–5] GHz, and compared to theoretical expectations in the far-field. Precise matches between MEEP simulation and radiation pattern predictions at different frequencies and beam angles are demonstrated. Given that the computations match the theory, the effect of embedding a phased array within a medium of varying index of refraction is then computed. Understanding the effect of varying index on phased arrays is critical for proposed ultra-high energy neutrino observatories which rely on phased array detectors embedded in natural ice. Future work will develop the phased array concepts with parallel MEEP, in order to increase the detail, complexity, and speed of the computations.

1. Introduction

Radio-frequency phased array antenna systems with design frequencies of order 0.1–10 GHz have applications in 5G mobile telecommunications, ground penetrating radar (GPR) systems, and scientific instrumentation [1,2,3,4]. In the one-dimensional case, a series of three-dimensional antenna elements are arranged in a line with fixed spacing [5]. Common antenna designs like loops and dipoles can be used to limit the elements to two dimensions. In this special case, phased array radiation may be modeled in two spatial dimensions plus time. In the two-dimensional case, a series of three-dimensional antenna elements are arranged in a two-dimensional pattern, often a grid with fixed element spacing in both dimensions. The elements may be strictly two-dimensional, but there is still an increase in computational complexity and the radiation is calculated in three dimensions plus time.

Proprietary RF modeling packages like XFDTD and HFSS are often used to model the response of elements within phased arrays and the behavior of arrays [6,7,8,9]. The XFDTD package, for example, relies on the finite difference time domain (FDTD) method. The FDTD approach is a computational electromagnetics (CEM) technique in which spacetime and Maxwell’s equations are broken into discrete form. One variant of the FDTD method is the conformal FDTD method (CFDTD), recently used to study phased array concepts on a large scale [9]. The NEC2 and NEC4 family of codes relies on the method-of-moments (MoM) approach [10]. Aside from the cost, a drawback of proprietary modeling software can be a lack of fine control over each individual object in the simulation. Because Maxwell’s equations are scale-invariant, in principle open-source FDTD codes designed for optical regimes could be re-purposd for RF design workflows. One such open-source package is the MIT Electromagnetic Equation Propagation (MEEP) package [11].

A recent review [12] covered how open radio design software like openEMS [13], gprMax [14], and the NEC2 family of codes [10] facilitate design workflows. In this work, the radiation patterns of one-dimensional and two-dimensional phased array designs are simulated with the MEEP package. MEEP takes advantage of the scale-invariance of Maxwell’s equations. Common MEEP applications are found in optical wavelength $\mathsf{\mu }$m-scale designs, but scale-invariance allows the user to treat designs as cm-scale RF elements (see Appendix A for details). Although MEEP has been used to optimize antenna designs [15], this work appears to be the first to model entire phased arrays in MEEP with a variable index of refraction. Two classes of phased array element are considered: Yagi-Uda and horn antennas. The former is applied to single-frequency designs, while the latter is applied to broadband design. Each element class is treated in both the one-dimensional and two-dimensional cases. The phase-steering properties and radiation patterns of all designs are shown to match theoretical predictions. The appropriate array theory is shown in Section 2, based on Chapter 1 of Reference [16]. Section 3 contains comparisons between theory and simulation for one-dimensional cases, and Section 4 contains the corresponding two-dimensional comparisons. In Section 5, the varying index of refraction is introduced. Results and future work are summarized in Section 6.

A workflow using MEEP for phased array design is outlined in Figure 1 based on Figure 1 from Reference [12]. Examples of decisions within the specifications category are: single-frequency or broadband, desired directivity and beamwidth, side-lobe tolerance, and number of antenna elements. These decisions lead to the choice of element type which must be implemented in MEEP. Simple shapes like dipoles can be modeled with built-in MEEP objects. Complex shapes like horns and dishes can be assembled from groups of objects. Radiation sources and current functions must be defined. For these studies, pure sinusoidal currents are passed to radiators which in turn radiate sinusoidal fields. The dielectric constant and boundary conditions of the simulation volume and objects within the volume are defined in the next step. The information is loaded into a simulation object and run for a number of time-steps. Once complete, near-to-far field routines are called to produce the power at a set of angles. The power versus angle is converted to normalized E and H-plane array radiation patterns and compared to theoretical models. Given a match, the frequency is updated and the process is repeated. If there is not a match, element separation and other array parameters are adjusted.

The workflow in Figure 1 represents a non-parallelized approach. Much development has gone into enhancing the speed, accuracy, and utility of the FDTD method. First, MEEP itself may be run in parallel mode, providing a speed enhancement. In a high-performance computing (HPC) environment, where each node has allocated memory (implying local RAM is not the limiting factor) running MEEP in parallel would speed up results. There has also been CEM research devoted to enhancing the FDTD approach itself. Decreasing memory usage and avoiding repetitive computations in favor of a more subtle approach is presented in [7]. A three-dimensional implementation of FDTD algorithms on GPUs via CUDA has also been explored [17]. The results of this work were obtained using the simplest version of MEEP: non-parallel with the python3 interface run in Jupyter notebooks on a laptop. Therefore the results shown in Section 3 and Section 4 could benefit from speed and memory enhancements in future studies.

2. Phased Array Antenna Theory

The basic structure of a one-dimensional phased array of RF radiating elements is shown in Figure 2. Two important numerical constants that determine the beam angle $\mathrm{\Delta }\varphi$ of the array are the inter-element spacing $d_y$ and the phase shift per antenna $\mathrm{\Delta }\mathrm{\Phi }$. Letting the subscript i label each of the N elements, the one-dimensional inter-element spacing in Figure 2 is $d_y\widehat{j}={\overrightarrow{r}}_{i+1}-{\overrightarrow{r}}_i$, where ${\overrightarrow{r}}_i$ records the position of element i. The phase shift per antenna is $\mathrm{\Delta }\mathrm{\Phi }={\mathrm{\Phi }}_{i+1}-{\mathrm{\Phi }}_i$. The relationship between $d_y$, $\mathrm{\Delta }\mathrm{\Phi }$, and $\mathrm{\Delta }\varphi$ is derived in Section 2.1. The radiation pattern for a given $\mathrm{\Delta }\varphi$ is derived in Section 2.2. For all coordinate systems, the azimuthal angle in the xy-plane is $\varphi$, and the polar angle from the z-axis is $\theta$.

2.1. Phase Steering and Beam Angle

The beam angle $\mathrm{\Delta }\varphi$ of the array given $\mathrm{\Delta }\mathrm{\Phi }$ and $d_y$ will now be derived for the coordinate system in Figure 2. First, the relevant far-field approximation will be described. Second, it will be assumed that the elements all radiate at the same frequency $\omega$ and have the same vector radiation pattern $\overrightarrow{f}(\theta ,\varphi )$ that accounts for co-polarized and cross-polarized radiated power. Third, the $\overrightarrow{E}$-field at point P will be treated as a sum of the ${\overrightarrow{E}}_i$ radiated from each element. Fourth, the calculations will be restricted to the xy-plane and the relationship between the beam angle $\mathrm{\Delta }\varphi$ and array parameters will be obtained for a one-dimensional array.

According to Figure 2, the position of P can be written

$$\overrightarrow{R}={\overrightarrow{r}}_i+{\overrightarrow{R}}_i$$

(1)

Rearranging, the displacement between the i-th element and P is

$${\overrightarrow{R}}_i=\overrightarrow{R}-{\overrightarrow{r}}_i$$

(2)

The magnitude of the displacement is

$$R_i=\sqrt{(\overrightarrow{R}-{\overrightarrow{r}}_i)·(\overrightarrow{R}-{\overrightarrow{r}}_i)}={\left(R^2-2\overrightarrow{R}·{\overrightarrow{r}}_i+r_i^2\right)}^{1/2}$$

(3)

Factoring an $R^2$, and neglecting the third term because it is small compared to the others,

$$R_i\approx R{\left(1-\frac{2\overrightarrow{R}·{\overrightarrow{r}}_i}{R^2}\right)}^{1/2}$$

(4)

Expanding in a Taylor series to first order in $2\overrightarrow{R}·{\overrightarrow{r}}_i/R^2$, with $\widehat{r}=\overrightarrow{R}/R$, yields

$$R_i\approx R\left(1-\frac{\widehat{r}·{\overrightarrow{r}}_i}{R}\right)$$

(5)

Distributing the R gives the approximation:

$$R_i\approx R-\widehat{r}·{\overrightarrow{r}}_i$$

(6)

The electric field at P due to the i-th element with individual radiation pattern ${\overrightarrow{f}}_i(\theta ,\varphi )$ is

$${\overrightarrow{E}}_i(R,\theta ,\varphi )=\frac{{\overrightarrow{f}}_i(\theta ,\varphi )exp(-jk{R}_i)}{R_i}$$

(7)

Substituting Equations (6) into (7):

$${\overrightarrow{E}}_i(R,\theta ,\varphi )={\overrightarrow{f}}_i(\theta ,\varphi )\frac{exp(-jkR)}{R}exp(jk{\overrightarrow{r}}_i·\widehat{r})$$

(8)

The element positions are written in Cartesian coordinates, while P is written in spherical coordinates using $u=sin\theta cos\varphi$ and $v=sin\theta sin\varphi$:

$$\begin{array}{cc}\hfill {\overrightarrow{r}}_i& =\widehat{x}{x}_i+\widehat{y}{y}_i+\widehat{z}{z}_i\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \widehat{r}& =\widehat{x}u+\widehat{y}v+\widehat{z}cos\theta \hfill \end{array}$$

(10)

The total field $\overrightarrow{E}$ at P requires summing over elements. The current delivered to the i-th element could have a potentially complex amplitude $a_i$. The details of how the currents $a_i$ are converted to radiated $\overrightarrow{E}$-field are taken to be part of $\overrightarrow{f}(\theta ,\varphi )$. The summation for $\overrightarrow{E}$ over elements is

$$\overrightarrow{E}(R,\theta ,\varphi )=\frac{exp(-jkR)}{R}\sum _i{a}_i{\overrightarrow{f}}_i(\theta ,\varphi )exp(jk{\overrightarrow{r}}_i·\widehat{r})$$

(11)

For identical radiating elements: ${\overrightarrow{f}}_i=\overrightarrow{f}$:

$$\overrightarrow{E}(R,\theta ,\varphi )=\overrightarrow{f}(\theta ,\varphi )\frac{exp(-jkR)}{R}\sum _i{a}_{i}exp(jk{\overrightarrow{r}}_i·\widehat{r})$$

(12)

Define the array-factor $F(\theta ,\varphi )={\sum }_i{a}_{i}exp(jk{\overrightarrow{r}}_i·\widehat{r})$:

$$\overrightarrow{E}(R,\theta ,\varphi )=\overrightarrow{f}(\theta ,\varphi )\frac{exp(-jkR)}{R}F(\theta ,\varphi )$$

(13)

Thus, if $F=1$, then the $\overrightarrow{E}$-field is a plane wave, modified only by the elemental radiation pattern. Complex amplitudes $a_i$ that cause a plane wave with wavevector pointed to $({\theta }_0,{\varphi }_0)$ are

$$a_i=\left|a_i\right|exp(-jk{\overrightarrow{r}}_i·{\widehat{r}}_0)$$

(14)

The notation for beam angle $\mathrm{\Delta }\varphi =\varphi -{\varphi }_0$ will be introduced shortly. For ${\widehat{r}}_0$, $u_0$ and $v_0$ take the corresponding ${\theta }_0$ and ${\varphi }_0$ for the angles: ${\widehat{r}}_0=\widehat{x}{u}_0+\widehat{y}{v}_0+\widehat{z}cos{\theta }_0$. The angles $({\theta }_0,{\varphi }_0)$ correspond to the plane wave because the phases in the array factor in Equation (13) are cancelled by those in Equation (14), and the summation is over just the magnitudes $|a_i|$. For a linear array in one-dimension, oriented along the y-axis as shown in Figure 2, ${\theta }_0=\pi /2$ and ${\overrightarrow{r}}_i=i{d}_y\widehat{y}$:

$$\overrightarrow{E}(R,\theta ,\varphi )=\overrightarrow{f}(\theta ,\varphi )\frac{exp(-jkR)}{R}\sum _i{a}_{i}exp\left(jk\left(i{d}_{y}v\right)\right)$$

(15)

The summation is $F(\pi /2,{\varphi }_0)$, $v=sin\left(\varphi \right)$ and $v_0=sin\left({\varphi }_0\right)$. The weights $a_i$ may be arranged to produce a plane wave at ${\varphi }_0$:

$$a_i=\left|a_i\right|exp\left(-jki{d}_y{v}_0\right)$$

(16)

The i-th phase of $\overrightarrow{E}$ in the array factor is

$${\mathrm{\Phi }}_i=ki{d}_y(sin\varphi -sin{\varphi }_0)$$

(17)

The difference $\mathrm{\Delta }\mathrm{\Phi }={\mathrm{\Phi }}_{i+1}-{\mathrm{\Phi }}_i$ for angles not far from the x-axis, $\left|\varphi \right|<1$ and $|{\varphi }_0|<1$, is

$$\mathrm{\Delta }\mathrm{\Phi }\approx d_{y}k(\varphi -{\varphi }_0)=2\pi (d_y/\lambda )(\varphi -{\varphi }_0)=2\pi (d_y/\lambda )\mathrm{\Delta }\varphi$$

(18)

The beam angle is $\mathrm{\Delta }\varphi =\varphi -{\varphi }_0$, the angular distance between a reference angle and the angle at which all contributions to $\overrightarrow{E}$ are in phase. Equation (18) reveals that the relationship between $\mathrm{\Delta }\varphi$ and $\mathrm{\Delta }\mathrm{\Phi }$ is linear, with slope $\lambda /\left(2\pi d_y\right)$. In Section 3, the relationship between $\mathrm{\Delta }\mathrm{\Phi }$ and $\mathrm{\Delta }\varphi$ obtained from FDTD calculations via MEEP are shown to match precisely the theoretical prediction. For two-dimensional grid arrays, the relationship “factors,” in that phase shift per element row and phase shift per element column govern $\mathrm{\Delta }\varphi$ and $\mathrm{\Delta }\theta$ independently. This theoretical prediction is matched precisely by the FDTD calculations shown in Section 4 as well.

2.2. Radiation Patterns and Beam Width

The radiation pattern, or relative power P emitted versus beam angle, is obtained from the array factor $F(\pi /2,\varphi )$ summation. Summation over the phased array with identical elements causes the vector element pattern $\overrightarrow{f}(\theta ,\varphi )$ and the common phase and amplitude factors $exp\left(jkR\right)/R$ to cancel upon normalization. The parameters that characterize the radiation patterns of arrays are N, the number of elements, and $d_y/\lambda$. The magnitude of the complex current to each element is assumed to be the same, $|a_i|=a$. Recall the array factor from Equation (13), with $\theta =\pi /2$ and $a_i=a$:

$$F(\varphi ,{\varphi }_0)=a\sum _{i}exp\left(jki{d}_y(v-v_0)\right)$$

(19)

Let $\chi =k{d}_y(v-v_0)$ so that $z=exp\left(j\chi \right)$. The sum is a geometric series from $i=1$ to $i=N$, the number of elements:

$$F\left(z\right)=a\sum _{i=1}^N{z}^i=a\left(\frac{1-z^N}{1-z}\right)$$

(20)

Using the Euler formula for $sin\left(\chi \right)$, the array factor simplifies to

$$F\left(\chi \right)=-aexp(j(N-1)\chi /2)\left(\frac{sin(N\chi /2)}{sin(\chi /2)}\right)$$

(21)

The radiation pattern is proportional to power, so it is prudent to take the magnitude of $F\left(\varphi \right)$:

$$\left|F\left(\chi \right)\right|=a\left(\frac{sin(N\chi /2)}{sin(\chi /2)}\right)$$

(22)

The normalized radiation pattern will be ${(F/F_{max})}^2$, so it is necessary to find $F_{max}$:

$$\underset{\chi \to 0}{lim}\left|F\left(\chi \right)\right|=a\underset{\chi \to 0}{lim}\left(\frac{sin(N\chi /2)}{sin(\chi /2)}\right)=aN$$

(23)

So $\left|F\left(\chi \right)\right|/F_{max}$ is

$$\frac{F\left(\chi \right)}{F_{max}}=\frac{sin(N\chi /2)}{Nsin(\chi /2)}$$

(24)

Finally, with $\chi =k{d}_y(v-v_0)$, $v=sin\left(\varphi \right)$, and $v_0=sin\left({\varphi }_0\right)$, the radiation pattern P is ${\left|F\left(\chi \right)/F_{max}\right|}^2$:

$$P\left(\varphi \right)={\left(\frac{sin\left(\pi N(d_y/\lambda )(sin\left(\varphi \right)-sin\left({\varphi }_0\right))\right)}{Nsin\left(\pi (d_y/\lambda )(sin\left(\varphi \right)-sin\left({\varphi }_0\right))\right)}\right)}^2$$

(25)

The −3 dB beamwidth is $0.886\lambda /L$, where $L=(N-1)d_y$. In fact, Equation (19) is a function of $v-v_0$, so altering the $\mathrm{\Delta }\mathrm{\Phi }$ in the $a_i$ only rotates $P\left(\varphi \right)$ in $\varphi$-space, corresponding to a translation in v-space. The radiation pattern in Equation (25) is shown to match precisely the main beam of FDTD calculations via MEEP for one-dimensional arrays in Section 3. For two-dimensional grid arrays, the E and H plane radiation patterns “factor,” in that $P(\theta ,\varphi )=P\left(\theta \right)P\left(\varphi \right)$. In Section 4, precision matches for two-dimensional grid arrays are shown.

2.3. Regarding Array Radiation Patterns

Because one-dimensional and two-dimensional arrays are considered, some notes about radiation patterns are necessary. First, all one-dimensional array radiation patterns correspond to the E-plane (the xy-plane). The arrays are specified using elements situated in the xy-plane, and the array extends along the y-axis. Radiators are linearly polarized such that the E-plane at some radius r is $(rcos(\varphi ),rsin(\varphi ),0)$. The H-plane at r would be $(rsin(\theta ),0,rcos(\theta \left)\right)$, but this data is not relevant for a one-dimensional array. Second, the MEEP python routine get_farfield is evaluated at a radius $r\gg L$, the length of the array, to obtain the far-fields $\overrightarrow{E}$ and $\overrightarrow{H}$. Notice that not all open-source FDTD codes offer near-field to far-field transition modeling [12].

All two-dimensional phased array elements presented in Section 4 are arrayed in the yz-plane, and the E and H-planes have the same definitions as the one-dimensional case. However, the H-plane results have been shifted so that the main beam occurs at $\theta =0$ degrees, rather than the expected 90 degrees. This is purely for visual comparison to Equation (25) cast as $P\left(\theta \right)$ with ${\theta }_0=0$, and does not mean the phased array is radiating orthogonally to broadside. Equation (25) is matched to E and H-plane two-dimensional patterns, and both are normalized to 0 dB at peak power.

3. Phased Array Designs in One Dimension: Two-Dimensional Fields

Two antenna designs were considered in modeling one-dimensional phased arrays: Yagi-Uda and horn, corresponding to narrowband and wideband applications, respectively. The two designs are depicted in Figure 3 with associated parameters described in Section 3.1 below. The Yagi-Uda antennas have 6 elements with the same radius, oriented in the xy-plane: one reflector, one radiatior, three directors and a connecting boom. The current $a_i$ is connected only to the radiator. The horn antennas have three structures: the box containing the linearly polarized radiator, the radiator which is connected to $a_i$, and the curves of the horn. An exponential function $y=f\left(x\right)=k_{1}exp\left(k_2(x-a)\right)$ describes the curves (see Figure 3), and the origin is taken to be at the center of the back edge of the box. The constants are $k_1=a/2$ and $k_2=(1/c)ln(2d/a)$. The curves are built from n slices where $n=c/dx$. All objects comprising the antenna elements have the same metallic conductivity, and the surrounding volume has an index of refraction $n=1$. At the edge of the space is a layer 1 $\mathrm{\Delta }x$ unit thick called the perfectly matched layer (PML) which cancels reflections.

3.1. Phase Steering, Beam Angle, and Beamwidth

As described in Section 2, the beam angle is controlled by the phase shift per antenna. Simulation results were run with the parameters in

Table 1. Yagi-Uda: The first and second columns contain the geometric parameters describing the antenna elements for the Yagi-Uda array. Horn: The third and fourth columns contain those for the horn array. Scan-loss: The fifth through eighth columns contain scan loss ($S{L}_{dB}$) data, reported for different frequencies and different $d_y/\lambda$ values for the $N=16$ horn array.

Yagi-Uda		Horn		SLdb
Parameter	Value	Parameter	Value	f (GHz)	$\mathrm{\Delta }\mathrm{\Phi }$ (degrees)	${\mathit{d}}_{\mathit{y}}/\mathit{\lambda }$	${\mathit{SL}}_{\mathit{dB}}$
N	8.16	N	8.16	0.5	80	0.125	−11.6
L	7.20	a	0.95	1.0	80	0.25	−1.2
$d_1$	1.80	c	15.0	2.0	80	0.5	−1.0
r	3.75	d	3.8	4.0	80	1.0	−0.9
y	2.81	$dx$	0.1
z	1.24	$n=c/dx$	150
$d_y$	3.92	$d_y$	2d
resolution	6	resolution	6

for the one-dimensional arrays. The main results are shown in Figure 4. A discussion about scan loss below references data from Table 1.

The phase-steering results are shown in Figure 4. The y-axes of Figure 4 (top left) and (top right) are the beam angles of the Yagi-Uda arrays, divided by the beam widths. The x-axes for these graphs are the phase shifts per element. The top left plot and top right plots corresond to $N=8$ and $N=16$, respectively. For the $N=16$ horn case (bottom left and right), the value of $d_y/\lambda =f{d}_y/c$ varies because the elements can radiate from ≈0.3–5.0 GHz. The black solid lines in the top left and top right graphs of Figure 4 are linear fits to the Yagi-Uda data. The gray lines represent the function $f\left(x\right)=bx$, with $b=\lambda /\left(2\pi d_y\right)$. For these models, $d_y=3.92$ cm and $\lambda =6$ cm (Table 1). The slopes match almost exactly, with slight errors arising from radiation pattern distortion at high beam angle $\mathrm{\Delta }\varphi$. At such large $\mathrm{\Delta }\mathrm{\Phi }$ values, side lobes can shift the location of the main beam by $\mathcal{O}\left(1\right)$ degree by merging with the main beam. In the $N=8$ case the fitted slope is slightly higher, and in the $N=16$ case the fitted slope is slightly lower. In each model, the observed beamwidths are within 1% of the value predicted by Equation (18).

For the broadband horn case in the bottom left of Figure 4, three frequency cases are shown: 0.3, 1.5, and 3.0 GHz. The intercepts are all consistent with zero and the slopes scale correctly: dividing the frequency by a factor of 2 increases the slope by a factor of 2, and dividing by a factor of 10 increases it by a factor of 10. Graphs like the top left and top right of Figure 4 would be misleading for horn antennas since the beamwidth depends on frequency (bottom right). The fit parameters for beam width were $a=12.0\pm 0.1$ degree GHz, and $b=1.1\pm 0.2$ degrees. For these two-dimensional antennas, the $1/f$-dependence is a good description of the beamwidth across the [0.3–5 GHz] bandwidth. The constant term b is only necessary since the array has finite length L. The beamwidth ($BW$) scales inversely with array length L: $BW\approx 0.886\lambda /L$, from Equation (25).

A discussion of scan loss is merited when analyzing normalized radiation patterns, which are shown below in Section 3.2. Scan loss may be quantified as the peak power at the given beam angle divided by the peak power at a beam angle of zero degrees. In the form of an equation in decibels, scan loss becomes a subtraction:

$$S{L}_{dB}=P_{\varphi }-P_{\varphi =0}$$

(26)

The scan loss $S{L}_{dB}$ is shown for the $N=16$ one-dimensional horn array in Table 1 (right), as it varies with frequency and $d_y/\lambda$. The conservative value $\mathrm{\Delta }\mathrm{\Phi }=80$ degrees was chosen because it is associated with the largest beam angles that do not generate side lobes larger than −15 dB. Given the beam width of the $N=16$ design (5.04 degrees), this corresponds to a scan range of $\pm 20.16$ degrees. The largest beam angles tend to have the largest scan losses, so the numbers in Table 1 should be the most conservative. Scan losses of less than 1 dB are observed at high frequencies, but at $d_y/\lambda$ values that begin to admit large side lobes (Section 3.2).

3.2. Radiation Patterns

Radiation patterns in the E-plane from $N=16$ one-dimensional Yagi and horn arrays are shown in Figure 5 and Figure 6, respectively. As described above, the x-direction ($\mathrm{\Delta }\varphi =0$) corresponds to no phase shift per element ($\mathrm{\Delta }\mathrm{\Phi }=0$). The radiation patterns are normalized to the power at the beam angle, and are shown in blue. The red curves represent Equation (25) with the correct N-value and $d_y/\lambda$-value. Equation (25) is symmetric, with identical forward and backward lobes. The front-to-back or FB ratio would be 1.0 or 0 dB for a row of ideal point sources. Although there is no backplane in either simulated one-dimensional array, the FB ratios of ≤$-15$ dB are observed.

The Yagi-Uda results are shown for 2.5 and 5.0 GHz frequencies in Figure 5, with $\mathrm{\Delta }\mathrm{\Phi }=0,20,40,$ and 60 degrees. Though the radiating elements are 6 cm long, good agreement between simulation and Equation (25) is observed at both 2.5 GHz and 5.0 GHz, including side-lobes. The beamwidth is proportionally larger at 2.5 GHz relative to 5.0 GHz, and at 5.0 GHz, the theoretical −3 dB beam width of 5.0 degrees is achieved. The amplitudes of all side-lobes are limited to ≈$-15$ dB, except at the highest beam angles where scan losses are experienced. Finally, the effect of frequency on beam steering is evident. The same $\mathrm{\Delta }\mathrm{\Phi }$ does not generate as large a $\mathrm{\Delta }\varphi$ at higher frequencies because the slope implied by Equation (18) is proportional to $\lambda$.

The horn results are shown in Figure 6 for 0.5 GHz and 5.0 GHz frequencies, corresponding to the lower and upper end of the bandwidth. The phase shifts per element are $\mathrm{\Delta }\mathrm{\Phi }=0,10,20,$ and 30 degrees. The angular range of $\mathrm{\Delta }\mathrm{\Phi }$ is restricted relative to the Yagi-Uda case. Wideband systems experience a natural trade-off in angular range versus bandwidth. A $d_y/\lambda$ value that is acceptably smaller than one at low frequencies can grow larger with increasing frequency, leading to interference patterns. At $5.0$ GHz, the horns radiate at $\pm 45$ degrees from $\mathrm{\Delta }\varphi =0$. The prediction from Equation (25) is that these side-lobes, or grating lobes, are equal in relative power to the main beam. The actual array limits them to $-15$ dB, but only if $\left|\mathrm{\Delta }\mathrm{\Phi }\right|<35$ degrees. For larger phase shifts per element, the opposite side-lobe grows above $-15$ dB. If the beam is steered too far in the $-\widehat{\varphi }$-direction, the side-lobe on the $\widehat{\varphi }$ side grows, and vice versa.

The general features of the radiation pattern compare well to the theoretical prediction. The $1/f$-dependence of the main beamwidth is evident in Figure 6. Like the Yagi-Uda array, the minimum theoretical beamwidth is reached at the highest frequencies (Figure 4 bottom right). The mini-lobes that are partially merged with the main beam widen the beam, however, the beamwidth is calculated at angles corresponding to −3 dB relative power. Since the mini-lobes are below −3 dB, the beamwidth calculation is unaffected. The simulation also matches the location and width of side-lobes to the theoretical prediction across the bandwidth. The six grating lobes at 5 GHz are a result of the pattern multiplication theorem, which states that the normalized radiation pattern is a product of the horn pattern and the pattern of an array of point sources. At 5 GHz, this multiplication suppresses the horn element pattern in the multiplication.

The field magnitude $|\overrightarrow{E}(x,y,t)|$ for the $N=16$ horn array is shown in Figure 7 for $t=0.5,1.0,1.5,2.0$ ns at a frequency of $2.5$ GHz, along with the radiation pattern (far right). The original amplitude of the radiation source within the horns is $\pm 1$ units, and the color scale for radiated $|\overrightarrow{E}(x,y,t)|$ is $\pm 0.002$. The $\mathrm{\Delta }\varphi$ is 9 degrees, with $\mathrm{\Delta }\mathrm{\Phi }=-35$ degrees, and the beamwidth is $5.5\pm 0.5$ degrees. The area is $480\times 900$ pixels describing $80\times 150$ cm${}^2$ with a resolution of 6 pixels per $\mathrm{\Delta }x$. The dimensions of the box (Table 1: $a=0.95$ cm) are smaller than $\lambda =12$ cm. As the radiation escapes to free space, the wavefront forms several $\lambda$ in front of the horns. Higher-frequency modes with $f\gg 2.5$ GHz are observed at the wavefront that correspond to start-up effects at $t=0$ ns. With 72 pixels per wavelength, minute features of $|\overrightarrow{E}(x,y,t)|$ can be interpreted as physical rather than numerical. These features can be eliminated with amplitude smoothing near $t=0$ ns, a feature available in the MEEP CustomSource class. Amplitude smoothing, however, makes the location of the wavefront less precise.

The radiation pattern in Figure 7 matches the theoretical prediction for the main lobe and first few side lobes. The origin of the side lobes is apparent from the $|\overrightarrow{E}(x,y,t)|$ images, where diffraction patterns at the edges of the array are visible. At 0.5 ns, radiation in an element is confined to the horn. By 1.0 ns, that radiation joins the waves from horns on either side. However, horns at the end of the array have no partner on one side, and some radiation leaks outside the main lobe. The side lobes can change if the total run time is not sufficient. The get_farfield routine in MEEP requires Near2FarRegion surfaces that form the near-field box that collects flux information at the radiation frequency. The parameters of the near-field box are set by the add_near2far routine. The get_farfield routine performs a near-to-far field projection to the given radius ($r=1000$ cm) where the field is computed for the radiation pattern. The propagation code must be run for sufficient units of $\mathrm{\Delta }t$ so that enough radiation can cross the near-field box. The code is run for 6.67 ns to generate the radiation pattern in Figure 7. Thus, the side lobes are averaged over many radiation periods.

4. Phased Array Designs in Two Dimensions: Three-Dimensional Fields

For two-dimensional grids of radiating elements, the array-factor $F(u,v)$ factors:

$$F(\theta ,\varphi )=F(u-u_0)F(v-v_0)$$

(27)

The radiation pattern in Equation (25) applies to the E and H plane separately. The two-dimensional arrays modeled below are square $N\times N$ arrays, so beamwidths implied by Equation (25) are equal for the E and H planes. The complex phasing of Equation (27) also indicates that $\mathrm{\Delta }{\varphi }_E\propto \mathrm{\Delta }{\mathrm{\Phi }}_E$, and $\mathrm{\Delta }{\varphi }_H\propto \mathrm{\Delta }{\mathrm{\Phi }}_H$, as shown in Equation (18) for the one-dimensional case. For the designs presented, the H-plane corresponds to the xz-plane, and to varying the phase in the z-direction (by array row). The E-plane corresponds to the xy-plane, and to varying the phase in the y-direction (by array column). In Figure 8, the basic shape of the two-dimensional array is shown in the yz-plane with $\mathrm{\Delta }{\mathrm{\Phi }}_E=15$ degrees, and $\mathrm{\Delta }{\mathrm{\Phi }}_H=15$ degrees. Section 4.1 contains results along the lines of Section 3.1 but for two-dimensional Yagi and horn arrays, and Section 4.2 contains results along the lines of Section 3.2 but for two-dimensional Yagi and horn arrays. As before, Table 1 contains the typical run parameters, with a few important exceptions.

The first exception is that for the two-dimensional case N becomes $N\times N$. However, squaring the number of antennas raises the memory requirements. In order to stay within a 16 GB memory limit, the two-dimensional horn array results had to be restricted to $N\times N=8\times 8$. The 2D horn array still has over $\mathcal{O}\left({10}^4\right)$ metal objects, compared to the $\mathcal{O}\left({10}^3\right)$ objects for the $N\times N=16\times 16$ 2D Yagi-Uda. The typical memory consumption is listed in

Table 2. The parameters for the $N\times N=8\times 8$ horn array, modified from Table 1. The Yagi-Uda $N\times N=16\times 16$ array did not require modification. The number of CPU cores was 4 in hardware, but was effectively 8 with hyperthreading. The most memory-intensive simulation was the $8\times 8$ horn array, which consumed 11.7 GB of memory out of 15.5 GB free. The code was written with the Python3 interface to MEEP, installed with the conda package manager, and run in Jupyter notebooks.

Horn
Parameter	Value	System Information
$N\times N$	$8\times 8$	Memory Consumption
a	2.0	11.7 GB out of 15.5 GB
c	15.0
d	8.0	CPU cores
$dx$	0.5	Intel i7 1.80 GHz (8)
$n=c/dx$	30
$d_y$	16	MEEP installation
resolution	4	Python3 interface (conda)
backplane location	$-2a$
backplane thickness	0.5
backplane dim.	142 × 142

, along with modified run paramters. Further, the resolution parameter was restricted to 4.0 for the horns. Restricting to 4 pixels per $\mathrm{\Delta }x$ unit limits memory consumption, but then the box containing the radiator has too few pixels. Enlarging the box allows the proper sized radiator to be fully contained. A final object was added to reduce the FB ratio: a back-plane with parameters listed in Table 2. One interesting modification is the doubling of the ratio of the box size (a) and the final horn width (d). This had the effect of limiting the maximum frequency to ≈1 GHz. At 1 GHz, $d/\lambda \approx 0.5$. A full optimization study on the horn parameters is warranted, though outside the present scope.

The horn elements radiate linearly polarized radiation in the y-direction, so the width in the z-direction does not follow the exponential functions but remains fixed at a. Initial runs were performed with horn elements that simultaneously widened according to the exponential function defined in Section 3. That design allowed reflections internal to the horns to distort the initial wavefront. Holding the horn-width constant in the z-direction produces radiation patterns that match Equation (25) because it follows the one-dimensional example of Section 3. To obtain z-polarized wavefronts, all that is necessary is to rotate the array. Practically, there are already examples of dually polarized RF band horns used in particle astrophysics [18,19], meaning that if this design were created with such elements, no rotation would be necessary.

4.1. Phase Steering, Beam Angle, and Beamwidth

The $\mathrm{\Delta }\varphi$ vs. $\mathrm{\Delta }\mathrm{\Phi }$ results for the two-dimensional $N\times N=16\times 16$ Yagi-Uda array are shown in Figure 9. Figure 9 (top left) contains $\mathrm{\Delta }{\varphi }_E$ versus $\mathrm{\Delta }{\mathrm{\Phi }}_E$ data at 5 GHz. The data match the theoretical linear slope $\lambda /\left(2\pi d_y\right)$ and $\lambda /\left(2\pi d_z\right)$, with $d_y=d_z$. The phase shift per antenna is varied over $[0,75]$ degrees in 15 degree increments independently by row and column. The circles and squares correspond to $\mathrm{\Delta }{\mathrm{\Phi }}_H=0$ and 45 degrees, respectively. Both the circles and squares follow the same line, implying the correct phase independence: when $\mathrm{\Delta }{\mathrm{\Phi }}_H$ is held at either constant, $\mathrm{\Delta }{\varphi }_E$ still varies with $\mathrm{\Delta }{\mathrm{\Phi }}_E$ correctly. Figure 9 (bottom left) contains $\mathrm{\Delta }{\varphi }_H$ versus $\mathrm{\Delta }{\mathrm{\Phi }}_H$ data at 5 GHz. The circles and squares correspond to $\mathrm{\Delta }{\mathrm{\Phi }}_E=0$ and 45 degrees, respectively. Beyond $\mathrm{\Delta }{\mathrm{\Phi }}_E=75$ degrees or $\mathrm{\Delta }{\mathrm{\Phi }}_H=75$ degrees, side lobes appear ($>-15$ dB). Figure 9 (right) contains the beam angle results after steering the beam to 36 of the possible $11\times 11$ positions in the E and H plane using increments of $\mathrm{\Delta }{\mathrm{\Phi }}_{E/H}=15$ degrees. The xy-errorbars correspond to the beamwidths.

The $\mathrm{\Delta }\varphi$ vs. $\mathrm{\Delta }\mathrm{\Phi }$ results for the two-dimensional $N\times N=8\times 8$ horn array are shown in Figure 10 (left). Figure 10 (top left) contains $\mathrm{\Delta }{\varphi }_E$ versus $\mathrm{\Delta }{\mathrm{\Phi }}_E$ data at 1 GHz. The larger horn size relative to those in Section 3 means the upper frequency is ≈1.2 GHz. The data match the theoretical slopes just as in Figure 9. The phase shift per antenna is varied in the same pattern as in Figure 9. The circles and squares correspond to $\mathrm{\Delta }{\mathrm{\Phi }}_H=0$ and 45 degrees, respectively, and both data sets follow the theory. Figure 10 (bottom left) contains $\mathrm{\Delta }{\varphi }_H$ versus $\mathrm{\Delta }{\mathrm{\Phi }}_H$ data at 1 GHz. The circles and squares correspond to $\mathrm{\Delta }{\mathrm{\Phi }}_E=0$ and 45 degrees, respectively, and both data sets follow the theory. The beamwidth as a function of frequency across the bandwidth for the design is shown in Figure 10 (right). The fit parameter mean values and standard errors are: $a=10.6\pm 0.2$ degree GHz, $b=2.8\pm 0.3$ degrees, $c=8.7\pm 0.4$ degree GHz, and $d=5.0\pm 0.9$ degrees. The width of the mouth of the horns is 16 cm in the E-plane direction and 2 cm in the H-plane direction, so some small difference in beamwidth is not surprising.

4.2. Radiation Patterns

The radiation patterns in the E and H plane for the two-dimensional Yagi-Uda array are shown in Figure 11. The phase combinations $(\mathrm{\Delta }{\mathrm{\Phi }}_E,\mathrm{\Delta }{\mathrm{\Phi }}_H)=(0,0),(30,60),(60,30)$ degrees are shown for E and H planes at 3 and 4 GHz. As in Section 3.2, Equation (25) is shown in red, and the simulation results are shown in blue. The main beam and first several side lobes are modeled correctly in each case, and the FB ratio is $\le -15$ dB. The side lobes are also at the $\approx -15$ dB level. Following Figure 10 (right), the main beam is narrower at 4 GHz than at 3 GHz. Though not generally designed to be broadband elements, the Yagi-Uda elements do display some flexibility in frequency. The log-periodic dipole array (LPDA) is a broadband example constructed from dipoles as the Yagi is [20].

Producing the radiation patterns in Figure 11 requires only $\mathcal{O}\left(10\right)$ seconds to run near-field calculations, and only another ≈5 min each to run the get_farfield routine over the E and H-planes. Modeling arrays constructed from dipole elements is orders of magnitude faster than for the array of horns, due to the two-dimensional nature of the dipole elements. Producing the radiation patterns of Figure 12 for the two-dimensional horn array requires ≈60 min combined for the E and H-plane patterns, per frequency. Unlike the Yagi case, the vast majority of time is not dedicated to the get_farfield routine, but to the near-field calculations. The near-field calculations require “sub-pixel smoothing” for the many edges of the blocks that comprise the horn structure.

The radiation patterns in the E and H plane for the two-dimensional horn array are shown in Figure 12. The phase combinations $(\mathrm{\Delta }{\mathrm{\Phi }}_E,{\mathrm{\Phi }}_H)=(0,0),(30,60),(60,30)$ degrees are shown for E and H planes at 0.5 and 1.0 GHz. As in Section 3.2, Equation (25) is shown in red, and the simulation results are shown in blue. The main beam and first several side lobes are modeled correctly in each case, and the FB ratio is ≤$-15$ dB. The side lobes are also at the ≈$-15$ dB level. Due to the higher bandwidth, a wider range of beamwidths is available (see Figure 10). The main beam is narrower at 1 GHz than at 0.5 GHz. The horns produce the correct pattern for $(\mathrm{\Delta }{\mathrm{\Phi }}_E,\mathrm{\Delta }{\mathrm{\Phi }}_H)=(0,0)$ degrees from 0.1 to 1.2 GHz. However, grating lobes above $-15$ dB are a known problem that occur when attempting to steer phased arrays built from broadband horns to wide angles (see Chapter 9 of Reference [16]). The addition of the backplane limits diffraction of the radiation around the edges of the array and therefore limits the FB ratio, but grating lobes appear at $\pm 45$ degrees from the main beam. There is occasionally a back lobe, which can be attributed to the diffraction of fields around the edge of the backplane. This effect is more pronounced when the main beam is steered to a wide angle and occurs in the hemisphere opposite to the main beam.

5. Variation of the Index of Refraction

The behavior of a one-dimensional phased array embedded within a dielectric medium with spatially-dependent index of refraction $n\left(z\right)$ is interesting to the ultra-high energy (UHE) neutrino community [3,5]. Phased arrays represent an opportunity to lower the RF detection threshold for RF pulses generated by UHE neutrinos via the Askaryan effect. Antarctic ice is the most convenient and natural medium for Askaryan pulse detection, due to the RF transparency and large pristine volumes located in Antarctic and Greenlandic ice sheets and shelves [21,22,23]. The index of refraction varies within the ice because of the transition between surface snow ($\rho \approx 0.4$ g/cm${}^3$) and the solid ice below ($\rho =0.917$ g/cm${}^3$). Most recent and intricate studies of phased array beam behavior still assume a uniform medium [24,25,26]. Embedded phased arrays with varying $n\left(z\right)$ emit signals that curve in the direction of increasing $n\left(z\right)$.

The shadow zone is the volume of ice from which RF signals do not reach a receiver due to the excess curvature of the ray trace [27]. While there is evidence that RF signals can propagate horizontally through Antarctic ice [28], data from Greenland suggests the relative strength of the effect is small compared to the curved radiation [29]. Using the tools developed in this work, it is possible to map out the shadow zone for an embedded phased array radiating sinusoidal signals at fixed frequency. Intriguingly, when the phased array radiates, the grating lobe power reflect downward from the snow-air interface, and radiates into the shadow zone. Grating lobe power also refracts into the air above the interface. Grating lobe power leaves the array at a different angle than the main beam, so their presence in the shadow zone does not represent forbidden RF propagation.

A two-parameter fit to the $n\left(z\right)$ data versus depth z below the surface is given by [28]

$$n\left(z\right)=\left\{\begin{array}{cc}1\hfill & z>0\hfill \\ n_{\mathrm{ice}}-\mathrm{\Delta }nexp(z/z_0)\hfill & z\le 0\hfill \end{array}\right.$$

(28)

The fit parameters in Equation (28) come from Reference [28]: $\mathrm{\Delta }n=0.423\pm 0.004$ and $z_0=77\pm 2$ m, with $n_{\mathrm{ice}}=1.78$ for RF frequencies. These values are derived from the SPICE core data taken in 2015 near the South Pole, and are in statistical agreement with fits from data obtained by the RICE experiment (see also Reference [28]). Equation (28) was implemented in the one-dimensional horn array case, but the horn structure surrounding each radiating element was removed. The array is therefore a one-dimensional dipole array. Further, the length scale was reinterpreted to be meters rather than centimeters, which is an ability conferred by the scale invariant FDTD algorithms. In this medium, there is no fixed in-situ value of $\lambda$ so the $\lambda /4$ dipoles were spaced by $\lambda /2$ according to their free space $\lambda$ value. At the selected frequency of 200 MHz, the dipole length is 0.375 m, and the spacing is 0.75 m.

Figure 13 and Figure 14 contain the results of a $N=8$ one-dimensional dipole array embedded in a medium with the index profile in Equation (28). Figure 13 shows the schematic of the calculation, and Figure 14 shows the magnitude of the z-component of the z-polarized array. Figure 14 represents the same physical dimensions as Figure 13. Equation (28) was sampled 100 times vertically, and with a resolution parameter of 10, the effective $\mathrm{\Delta }z$ is 0.1 m. The units in Figure 13 are meters, and the unit-less frequency in MEEP was scaled accordingly, to correspond to 200 MHz. The distances between the air-snow interface and the first phased array element is 15 m (top) and 35 m (bottom).

The color scale in Figure 14 is $\pm 0.05$ with the signal amplitude of the elements $\pm 1.0$ at 200 MHz. The amplitude scale is less important than observing where the radiation has penetrated the ice after 200 time steps. In Figure 14 (top), the main beam has curved downwards in the direction of increasing $n\left(z\right)$, while grating lobes have both diffracted to the air and reflected into the shadow zone. The rate of curvature of the main beam is controlled by the fit parameter $z_0$ in Equation (28). In Figure 14 (bottom), the physics is the same as Figure 14 (top), but the effect of $n\left(z\right)$ curvature is weakened. The beam travels farther horizontally because the gradient of $n\left(z\right)$ is smaller at the larger depth. The geometry of the larger depth is such that the reflected grating lobe power is interfering with grating lobe power that was curved downwards without reflection. This can be seen just above the C marker in Figure 14 (bottom).

6. Summary and Future Analysis

Four phased array designs have been modeled with the MIT Electromagnetics Equation Propagation (MEEP) package in non-parallel mode. Two types of individual radiating element were explored: the narrow-band Yagi-Uda and broadband horn antennas. Two phased array geometries were explored: one-dimensional and two-dimensional. The one and two-dimensional Yagi-Uda phased arrays were designed for ≤5 GHz, however, scale-invariance makes this design scalable to a variety of frequencies. The one-dimensional horn array performs in the range [0.3–5] GHz, using two-dimensional versions of the horn elements. The two-dimensional horn array had to be modified due to memory constraints. The result was an array that performed in the range [0.1–1.2] GHz. In all cases, comparisons to array theory were shown.

The one-dimensional array of Yagi-Uda antennas was analyzed in Section 3. The array demonstrated the correct linear relationship between $\mathrm{\Delta }\varphi$ and $\mathrm{\Delta }\mathrm{\Phi }$ (Figure 4) (top left and right). Although any row of point sources would obey the relationship in Equation (18), a row of point sources has two main beam solutions by symmetry. Thus Figure 4 could not be interpreted correctly were it not for the proper functioning of the Yagi elements. The radiation patterns produced with the one-dimensional Yagi array were compared to Equation (25) in Figure 5. The radiation pattern in the E-plane is shown to agree with Equation (25) in both the main beam and the first several grating lobes. The calculation takes place in two-dimensions, so an H-plane comparison is not relevant.

The one-dimensional array of horn antennas was analyzed in Section 3. The array demonstrated the correct linear relationship between $\mathrm{\Delta }\varphi$ and $\mathrm{\Delta }\mathrm{\Phi }$ (Figure 4). In that case, the slope of $\mathrm{\Delta }\varphi$ vs. $\mathrm{\Delta }\mathrm{\Phi }$ was increased by a factor of 2 and then 10 by decreasing the frequency by a factor of 2 and then 10. The bandwidth of the two-dimensional versions of the horns allows the variation of scan range. The scan range is smaller at high frequencies, as indicated in Figure 4 (bottom left). However, the beamwidth is also smaller at high frequencies, as indicated in Figure 4 (bottom right). The design trade-off is between small beamwidth and large scan range. In Figure 6 the one-dimensional horn array radiation pattern is shown to match Equation (25) at both 0.5 and 5.0 GHz. There are 2–4 side lobes at 0.5 GHz to match per pattern, and the simulation results match them as well as the wide main beam. At 5.0 GHz, the main beam is accompanied by two prominent grating lobes at $\pm 45$ degrees that should be as powerful as the main beam. The simulation finds them at the −15 dB level. The grating lobes are being suppressed by the the pattern null from the horn element pattern [16]. At lower frequencies, however, scan loss takes a toll on radiated power (Table 1).

The two-dimensional, $N\times N=16\times 16$ Yagi-Uda array was analyzed in Section 4. The array demonstrated the correct linear relationships between $\mathrm{\Delta }{\varphi }_E$ and $\mathrm{\Delta }{\mathrm{\Phi }}_E$, and $\mathrm{\Delta }{\varphi }_H$ and $\mathrm{\Delta }{\mathrm{\Phi }}_H$ (Figure 9) (left). Given the narrow beamwidth, the array design can be scanned $\pm 5$ beamwidths in the E-plane and $\pm 5$ beamwidths in the H-plane before side lobes become too large. One fourth of these scan positions are shown in Figure 9 (right). The radiation pattern of the full two-dimensional array was displayed in Figure 11 at 3 and 4 GHz, for several scan angles. In each case, the pattern matched Equation (25) in the E and H-planes in the main beam and dominant side lobes. The addition of a metal back plane helps to suppress back lobes. Peculiarly, the two-dimensional array did not match the theoretical prediction at 5 GHz as well as the one-dimensional case when scanned.

The two-dimensional, $N\times N=8\times 8$ horn array was analyzed in Section 4. The array demonstrated the correct linear relationships between $\mathrm{\Delta }{\varphi }_E$ and $\mathrm{\Delta }{\mathrm{\Phi }}_E$, and $\mathrm{\Delta }{\varphi }_H$ and $\mathrm{\Delta }{\mathrm{\Phi }}_H$ (Figure 10) (left). The beamwidth is again inversely proportional to frequency (Figure 10) (right). It is not surprising that the fits differ slightly in the E and H-planes, since the horn width changes in the E-plane but does not in the H-plane. The quality of the fits to $1/f+const$ are excellent. The additive constants in these fits are only necessary because the array cannot be infinitely long. Technically, Equation (25) implies that the beamwidth would go to zero as $N\to \infty$. The radiation patterns of the two-dimensional horn array are displayed in Figure 12 at $0.5$ and $1.0$ GHz, for the same sampling of scan angles as in Figure 11. The high-frequency beam is narrower and accompanied by grating lobes at $\pm 45$ degrees. The patterns agree with theoretical expectations, with the exception of the H-plane lobes at $\pm 90$ degrees. At low frequency, the beam is wider and is accompanied by grating lobes at $\pm 45$ degrees from the main beam. The results match in the main lobe, but the simulation does not match the theoretical grating lobes. This is pronounced when the beam is moved far from broadside in the H-plane.

Finally, a simplified version of the $N=8$ one-dimensional case of dipoles was embedded in a medium with varying index of refraction, $n\left(z\right)$. The model for $n\left(z\right)$ was a simple fit to the profile of the ice at the South Pole, which is a location of interest for planned phased array detectors designed to record Askaryan signals from UHE neutrinos passing through ice. Though the studies in this work are restricted to phased-arrays as transmitters, and not receivers, the shadow zone of the array was mapped at 200 MHz under realistic conditions. An interesting side effect of the phased array being the radiating system was that the grating lobes managed to propagate into the shadow zone.

Future work would include several enhancements to the simulations. Calculations of S-parameters for individual elements should be added, and optimization studies on horn and Yagi geometric parameters are warranted. However, other RF element types should also be studied. Due to the relevance of one-dimensional phased array receivers for UHE neutrino physics, one interesting choice is the wide-radius dipole used by the Radio Neutrino Observatory Greenland (RNO-G) collaboration [30]. Such elements already have low VSWR measurements in the relevant bandwidth. Finally, upgrading the simulation code to utilize parallel MEEP capabilities will increase the potential speed and complexity. Additional complexity will come in the form of more accurate antenna structure modeling, thereby improving the precision across a wide bandwidth.