###### 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.

###### keywords:

FDTD methods, MEEP, phased array antennas, antenna theory, Askaryan effect,
UHE neutrinos

xx 10 4 Received: date; Accepted: date; Published: date Broadband RF Phased
Array Design with MEEP: Comparisons to Array Theory in Two and Three
Dimensions Jordan C. Hanson Jordan C. Hanson

##  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 Syrytsin et
al. (2017); Kikuchi et al. (2017); Vieregg et al. (2016); Munekata et al.
(2014). In the one-dimensional case, a series of three-dimensional antenna
elements are arranged in a line with fixed spacing Avva et al. (2017). 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
Ansys, Inc., Canonsburg, Pennsylvania (2020); Feng et al. (2019); Zhu et al.
(2017); Ho et al. . 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 Ho et al. . The NEC2 and NEC4 family of codes relies on the method-of-
moments (MoM) approach Burke et al. (2004). 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 Oskooi et al.
(2010).

A recent review Fedeli et al. (2019) covered how open radio design software
like openEMS Liebig et al. (2013), gprMax Warren et al. (2016), and the NEC2
family of codes Burke et al. (2004) 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 $\mu$m-scale designs, but scale-invariance allows the user
to treat designs as cm-scale RF elements (see Appendix for details). Although
MEEP has been used to optimize antenna designs Richie and Ababei (2017), 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 Mailloux (2018). 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.

Figure 1:  A detailed workflow for phased array design with MEEP. See text for
details.

A workflow using MEEP for phased array design is outlined in Figure 1 based on
Figure 1 from Reference Fedeli et al. (2019). 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 Feng et al. (2019). A three-dimensional implementation of FDTD
algorithms on GPUs via CUDA has also been explored Wahl et al. (2013). 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 Sections 3 and 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 $\Delta\phi$ of the array are the inter-element spacing $d_{y}$ and
the phase shift per antenna $\Delta\Phi$. Letting the subscript $i$ label each
of the $N$ elements, the one-dimensional inter-element spacing in Figure 2 is
$d_{y}\hat{j}=\vec{r}_{i+1}-\vec{r}_{i}$, where $\vec{r}_{i}$ records the
position of element $i$. The phase shift per antenna is
$\Delta\Phi=\Phi_{i+1}-\Phi_{i}$. The relationship between $d_{y}$,
$\Delta\Phi$, and $\Delta\phi$ is derived in Section 2.1. The radiation
pattern for a given $\Delta\phi$ is derived in Section 2.2. For all coordinate
systems, the azimuthal angle in the xy-plane is $\phi$, and the polar angle
from the z-axis is $\theta$.

Figure 2:  Definitions for the coordinate system, element label $i$, position
vectors, and phase shift per antenna for a one-dimensional phased array of RF
radiating elements. An example phase shift per antenna of
$\Delta\Phi=\Phi_{2}-\Phi_{1}=\Phi_{3}-\Phi_{2}=\Phi_{i+1}-\Phi_{i}=10^{\circ}$
value is shown. Example position vectors for the 12th element are shown:
$\vec{R}=\vec{r}_{12}+\vec{R}_{12}$.

###  2.1 Phase Steering and Beam Angle

The beam angle $\Delta\phi$ of the array given $\Delta\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 $\vec{f}(\theta,\phi)$ that accounts for co-polarized and
cross-polarized radiated power. Third, the $\vec{E}$-field at point P will be
treated as a sum of the $\vec{E}_{i}$ radiated from each element. Fourth, the
calculations will be restricted to the xy-plane and the relationship between
the beam angle $\Delta\phi$ and array parameters will be obtained for a one-
dimensional array.

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

$\vec{R}=\vec{r}_{i}+\vec{R}_{i}$ (1)

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

$\vec{R}_{i}=\vec{R}-\vec{r}_{i}$ (2)

The magnitude of the displacement is

$R_{i}=\sqrt{(\vec{R}-\vec{r}_{i})\cdot(\vec{R}-\vec{r}_{i})}=\left(R^{2}-2\vec{R}\cdot\vec{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\vec{R}\cdot\vec{r}_{i}}{R^{2}}\right)^{1/2}$
(4)

Expanding in a Taylor series to first order in
$2\vec{R}\cdot\vec{r}_{i}/R^{2}$, with $\hat{r}=\vec{R}/R$, yields

$R_{i}\approx R\left(1-\frac{\hat{r}\cdot\vec{r}_{i}}{R}\right)$ (5)

Distributing the $R$ gives the approximation:

$R_{i}\approx R-\hat{r}\cdot\vec{r}_{i}$ (6)

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

$\vec{E}_{i}(R,\theta,\phi)=\frac{\vec{f}_{i}(\theta,\phi)\exp(-jkR_{i})}{R_{i}}$
(7)

Substituting Equation 6 into Equation 7:

$\vec{E}_{i}(R,\theta,\phi)=\vec{f}_{i}(\theta,\phi)\frac{\exp(-jkR)}{R}\exp(jk\vec{r}_{i}\cdot\hat{r})$
(8)

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

$\displaystyle\vec{r}_{i}$
$\displaystyle=\hat{x}x_{i}+\hat{y}y_{i}+\hat{z}z_{i}$ (9)
$\displaystyle\hat{r}$ $\displaystyle=\hat{x}u+\hat{y}v+\hat{z}\cos{\theta}$
(10)

The total field $\vec{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
$\vec{E}$-field are taken to be part of $\vec{f}(\theta,\phi)$. The summation
for $\vec{E}$ over elements is

$\vec{E}(R,\theta,\phi)=\frac{\exp(-jkR)}{R}\sum_{i}a_{i}\vec{f}_{i}(\theta,\phi)\exp(jk\vec{r}_{i}\cdot\hat{r})$
(11)

For identical radiating elements: $\vec{f}_{i}=\vec{f}$:

$\vec{E}(R,\theta,\phi)=\vec{f}(\theta,\phi)\frac{\exp(-jkR)}{R}\sum_{i}a_{i}\exp(jk\vec{r}_{i}\cdot\hat{r})$
(12)

Define the array-factor
$F(\theta,\phi)=\sum_{i}a_{i}\exp(jk\vec{r}_{i}\cdot\hat{r})$:

$\vec{E}(R,\theta,\phi)=\vec{f}(\theta,\phi)\frac{\exp(-jkR)}{R}F(\theta,\phi)$
(13)

Thus, if $F=1$, then the $\vec{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},\phi_{0})$ are

$a_{i}=|a_{i}|\exp(-jk\vec{r}_{i}\cdot\hat{r}_{0})$ (14)

The notation for beam angle $\Delta\phi=\phi-\phi_{0}$ will be introduced
shortly. For $\hat{r}_{0}$, $u_{0}$ and $v_{0}$ take the corresponding
$\theta_{0}$ and $\phi_{0}$ for the angles:
$\hat{r}_{0}=\hat{x}u_{0}+\hat{y}v_{0}+\hat{z}\cos{\theta_{0}}$. The angles
$(\theta_{0},\phi_{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 $\vec{r}_{i}=id_{y}\hat{y}$:

$\vec{E}(R,\theta,\phi)=\vec{f}(\theta,\phi)\frac{\exp(-jkR)}{R}\sum_{i}a_{i}\exp\left(jk(id_{y}v)\right)$
(15)

The summation is $F(\pi/2,\phi_{0})$, $v=\sin(\phi)$ and
$v_{0}=\sin(\phi_{0})$. The weights $a_{i}$ may be arranged to produce a plane
wave at $\phi_{0}$:

$a_{i}=|a_{i}|\exp\left(-jkid_{y}v_{0}\right)$ (16)

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

$\Phi_{i}=kid_{y}(\sin\phi-\sin\phi_{0})$ (17)

The difference $\Delta\Phi=\Phi_{i+1}-\Phi_{i}$ for angles not far from the
x-axis, $|\phi|<1$ and $|\phi_{0}|<1$, is

$\Delta\Phi\approx
d_{y}k(\phi-\phi_{0})=2\pi(d_{y}/\lambda)(\phi-\phi_{0})=2\pi(d_{y}/\lambda)\Delta\phi$
(18)

The beam angle is $\Delta\phi=\phi-\phi_{0}$, the angular distance between a
reference angle and the angle at which all contributions to $\vec{E}$ are in
phase. Equation 18 reveals that the relationship between $\Delta\phi$ and
$\Delta\Phi$ is linear, with slope $\lambda/(2\pi d_{y})$. In Section 3, the
relationship between $\Delta\Phi$ and $\Delta\phi$ 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 $\Delta\phi$
and $\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,\phi)$ summation. Summation over the
phased array with identical elements causes the vector element pattern
$\vec{f}(\theta,\phi)$ and the common phase and amplitude factors
$\exp(jkR)/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(\phi,\phi_{0})=a\sum_{i}\exp\left(jkid_{y}(v-v_{0})\right)$ (19)

Let $\chi=kd_{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(z)=a\sum_{i=1}^{N}z^{i}=a\left(\frac{1-z^{N}}{1-z}\right)$ (20)

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

$F(\chi)=-a\exp(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(\phi)$:

$|F(\chi)|=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}$:

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

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

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

Finally, with $\chi=kd_{y}(v-v_{0})$, $v=\sin(\phi)$, and
$v_{0}=\sin(\phi_{0})$, the radiation pattern $P$ is
$\left|F(\chi)/F_{max}\right|^{2}$:

$\boxed{P(\phi)=\left(\frac{\sin(\pi
N(d_{y}/\lambda)(\sin(\phi)-\sin(\phi_{0})))}{N\sin(\pi(d_{y}/\lambda)(\sin(\phi)-\sin(\phi_{0})))}\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 $\Delta\Phi$ in the
$a_{i}$ only rotates $P(\phi)$ in $\phi$-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,\phi)=P(\theta)P(\phi)$. 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 $(r\cos(\phi),r\sin(\phi),0)$. The H-plane at $r$ would be
$(r\sin(\theta),0,r\cos(\theta))$, 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
$\vec{E}$ and $\vec{H}$. Notice that not all open-source FDTD codes offer
near-field to far-field transition modeling Fedeli et al. (2019).

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(\theta)$ 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(x)=k_{1}\exp(k_{2}(x-a))$ 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 $\Delta x$ unit thick
called the perfectly matched layer (PML) which cancels reflections.

(a) The Yagi-Uda antenna, and the $N=16$ array.

(b) The horn antenna, and the $N=16$ array.

Figure 3:  The two-dimensional antenna designs used in the one-dimensional
phased array simulations. (See Section 3.1 for details). The black region has
$\epsilon=\epsilon_{0}$, the green borders are perfectly matched layers (PML),
and the blue surface is a MEEP Near2FarRegion where flux is recorded for near-
to-far projection. The white lines represent metal structures, and the red
lines represent the radiating elements.

###  3.1 Phase Steering, Beam Angle, and Beamwidth

Yagi-Uda |  | Horn |  | Scan-loss, N=16 Horn array |  |  |
---|---|---|---|---|---|---|---
Parameter | Value | Parameter | Value | Frequency (GHz) |  $\Delta\Phi$ (degrees) | $d_{y}/\lambda$ | $SL_{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 |  |  |  |
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 ($SL_{dB}$) data, reported for
different frequencies and different $d_{y}/\lambda$ values for the $N=16$ horn
array.

Figure 4:  (Top left) The beam angle $\Delta\phi$ divided by the beam width
$BW$ for the $N=8$ one-dimensional Yagi array versus $\Delta\Phi$, the phase
shift per element. (Top right) The same results for the $N=16$ array. (Bottom
left) $\Delta\phi$ versus $\Delta\Phi$ for the $N=16$ version of the one-
dimensional horn array, for several frequencies. (Bottom right) The dependence
of the beam width on frequency for the one-dimensional $N=16$ horn array.

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 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=fd_{y}/c$ varies because the elements can radiate from $\approx
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(x)=bx$, with $b=\lambda/(2\pi d_{y})$. 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 $\Delta\phi$. At such large $\Delta\Phi$ values, side lobes can shift
the location of the main beam by $\mathcal{O}(1)$ 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:

$SL_{dB}=P_{\phi}-P_{\phi=0}$ (26)

The scan loss $SL_{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 $\Delta\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 Figures 5 and 6, respectively. As described above, the
x-direction ($\Delta\phi=0$) corresponds to no phase shift per element
($\Delta\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
$\leq-15$ dB are observed.

Figure 5:  Yagi-Uda results, two-dimensional elements, one-dimensional array.
(Top row) $f=2.5$ GHz, and $\Delta\Phi=0,20$ degrees from left to right.
(Second row) $f=2.5$ GHz, and $\Delta\Phi=40,60$ degrees from left to right.
(Third row) $f=5.0$ GHz, and $\Delta\Phi=0,20$ degrees from left to right.
(Bottom row) $f=5.0$ GHz, and $\Delta\Phi=40,60$ degrees from left to right.
The radial units are dB, and the angular units are degrees.

Figure 6:  Horn results, two-dimensional elements, one-dimensional array. (Top
row) $f=0.5$ GHz, and $\Delta\Phi=0,10$ degrees from left to right. (Second
row) $f=0.5$ GHz, and $\Delta\Phi=20,30$ degrees from left to right. (Third
row) $f=5.0$ GHz, and $\Delta\Phi=0,10$ degrees from left to right. (Bottom
row) $f=5.0$ GHz, and $\Delta\Phi=20,30$ degrees from left to right.

The Yagi-Uda results are shown for 2.5 and 5.0 GHz frequencies in Figure 5,
with $\Delta\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 $\approx-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 $\Delta\Phi$ does not generate as large a
$\Delta\phi$ 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 $\Delta\Phi=0,10,20,$ and $30$ degrees. The angular range of
$\Delta\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 $\Delta\phi=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 $|\Delta\Phi|<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 $-\hat{\phi}$-direction, the side-lobe on the $\hat{\phi}$ 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.

Figure 7:  (From left to right) The $N=16$ horn one-dimensional linearly
polarized electric field $|\vec{E}(x,y,t)|$ at $t=0.5$ ns, $t=1.0$ ns, $t=1.5$
ns, and $t=2.0$ ns. The 2D area is $80\times 150$ cm2, with resolution of 6
pixels per distance unit ($480\times 900$ pixels). The frequency is 2.5 GHz,
and the beam angle is 9 degrees from broadside. (Far right) The corresponding
radiation pattern.

The field magnitude $|\vec{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
$|\vec{E}(x,y,t)|$ is $\pm 0.002$. The $\Delta\phi$ is $9$ degrees, with
$\Delta\Phi=-35$ degrees, and the beamwidth is $5.5\pm 0.5$ degrees. The area
is $480\times 900$ pixels describing $80\times 150$ cm2 with a resolution of 6
pixels per $\Delta x$. The dimensions of the box (Tab. 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 $|\vec{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 $|\vec{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 $\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,\phi)=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
$\Delta\phi_{E}\propto\Delta\Phi_{E}$, and
$\Delta\phi_{H}\propto\Delta\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 $\Delta\Phi_{E}=15$ degrees,
and $\Delta\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.

Figure 8:  The two-dimensional $N\times N=16\times 16$ Yagi-Uda/horn
y-polarized array layout. The alignment with 3D Cartesian coordinates is
depicted, along with the array spacing variables $d_{y}$ and $d_{z}$, the
wavelength $\lambda$ (to scale), and the phases for each row and column if
$\Delta\Phi_{E}=\Delta\Phi_{H}=15$ degrees.

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}(10^{4})$ metal objects, compared to the
$\mathcal{O}(10^{3})$ objects for the $N\times N=16\times 16$ 2D Yagi-Uda. The
typical memory consumption is listed in Table 2, along with modified run
paramters. Further, the resolution parameter was restricted to 4.0 for the
horns. Restricting to 4 pixels per $\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 $\approx 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.

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 $\times$ 142 |
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.

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 Gorham et al.
(2009, 2007), meaning that if this design were created with such elements, no
rotation would be necessary.

###  4.1 Phase Steering, Beam Angle, and Beamwidth

The $\Delta\phi$ vs. $\Delta\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 $\Delta\phi_{E}$ versus $\Delta\Phi_{E}$ data at 5 GHz. The data
match the theoretical linear slope $\lambda/(2\pi d_{y})$ and $\lambda/(2\pi
d_{z})$, 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 $\Delta\Phi_{H}=0$ and $45$ degrees,
respectively. Both the circles and squares follow the same line, implying the
correct phase independence: when $\Delta\Phi_{H}$ is held at either constant,
$\Delta\phi_{E}$ still varies with $\Delta\Phi_{E}$ correctly. Figure 9
(bottom left) contains $\Delta\phi_{H}$ versus $\Delta\Phi_{H}$ data at 5 GHz.
The circles and squares correspond to $\Delta\Phi_{E}=0$ and $45$ degrees,
respectively. Beyond $\Delta\Phi_{E}=75$ degrees or $\Delta\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 $\Delta\Phi_{E/H}=15$
degrees. The xy-errorbars correspond to the beamwidths.

Figure 9:  (Left) The beam angles $\Delta\phi_{E}$ and $\Delta\phi_{H}$ versus
the phase shifts for the $N\times N=16\times 16$ Yagi-Uda array at 5 GHz. The
black lines represent the theoretical prediction of a linear dependence with
slope $\lambda/(2\pi d_{y})$ or $\lambda/(2\pi d_{y})$, ($d_{y}=d_{z}$). In
each graph, the circles and squares correspond to two different $\Delta\Phi$
constant values for the other array plane. In these examples, location of zero
phase on the array is chosen to cause a negative beam angle. (Right) The data
points correspond to beam angles in the E and H-planes, with the associated
beamwidths as errorbars. These data represent one-quarter of the possible scan
positions with $\Delta\Phi_{E/H}=15$ degrees.

The $\Delta\phi$ vs. $\Delta\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 $\Delta\phi_{E}$ versus $\Delta\Phi_{E}$ data at 1 GHz. The larger
horn size relative to those in Section 3 means the upper frequency is $\approx
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 $\Delta\Phi_{H}=0$ and $45$ degrees, respectively,
and both data sets follow the theory. Figure 10 (bottom left) contains
$\Delta\phi_{H}$ versus $\Delta\Phi_{H}$ data at 1 GHz. The circles and
squares correspond to $\Delta\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.

Figure 10:  (Left) The beam angles $\Delta\phi_{E}$ and $\Delta\phi_{H}$
versus the phase shifts per column (for the E-plane) and per row (for the
H-plane) for the $N\times N=8\times 8$ horn array at 1 GHz. The black lines
represent the theoretical prediction. In each graph, the circles and squares
correspond to two different $\Delta\Phi$ constant values for the other array
plane. (Right) The beamwidth in the E and H-planes versus frequency.

###  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
$(\Delta\Phi_{E},\Delta\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 $\leq-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 Barwick et al.
(2015).

Figure 11:  Yagi-Uda results, two-dimensional array(First row) $f=3$ GHz, with
(from left to right) E-plane $(\Delta\Phi_{E},\Delta\Phi_{H})=(0,0)$ degrees,
E-plane $(\Delta\Phi_{E},\Delta\Phi_{H})=(60,30)$ degrees, E-plane
$(\Delta\Phi_{E},\Delta\Phi_{H})=(30,60)$ degrees. (Second row) $f=3$ GHz,
with (from left to right) H-plane $(\Delta\Phi_{E},\Delta\Phi_{H})=(0,0)$
degrees, H-plane $(\Delta\Phi_{E},\Delta\Phi_{H})=(60,30)$ degrees, H-plane
$(\Delta\Phi_{E},\Delta\Phi_{H})=(30,60)$ degrees. (Third row) $f=4$ GHz, with
(from left to right) E-plane $(\Delta\Phi_{E},\Delta\Phi_{H})=(0,0)$ degrees,
E-plane $(\Delta\Phi_{E},\Delta\Phi_{H})=(60,30)$ degrees, E-plane
$(\Delta\Phi_{E},\Delta\Phi_{H})=(30,60)$ degrees. (Fourth row) $f=4$ GHz,
with (from left to right) H-plane $(\Delta\Phi_{E},\Delta\Phi_{H})=(0,0)$
degrees, H-plane $(\Delta\Phi_{E},\Delta\Phi_{H})=(60,30)$ degrees, H-plane
$(\Delta\Phi_{E},\Delta\Phi_{H})=(30,60)$ degrees.

Producing the radiation patterns in Figure 11 requires only $\mathcal{O}(10)$
seconds to run near-field calculations, and only another $\approx 5$ minutes
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 $\approx 60$ minutes 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.

Figure 12:  Horn results, two-dimensional array (First row) $f=0.5$ GHz, with
(from left to right) E-plane $(\Delta\Phi_{E},\Delta\Phi_{H})=(0,0)$ degrees,
E-plane $(\Delta\Phi_{E},\Delta\Phi_{H})=(60,30)$ degrees, E-plane
$(\Delta\Phi_{E},\Delta\Phi_{H})=(30,60)$ degrees. (Second row) $f=0.5$ GHz,
with (from left to right) H-plane $(\Delta\Phi_{E},\Delta\Phi_{H})=(0,0)$
degrees, H-plane $(\Delta\Phi_{E},\Delta\Phi_{H})=(60,30)$ degrees, H-plane
$(\Delta\Phi_{E},\Delta\Phi_{H})=(30,60)$ degrees. (Third row) $f=1.0$ GHz,
with (from left to right) E-plane $(\Delta\Phi_{E},\Delta\Phi_{H})=(0,0)$
degrees, E-plane $(\Delta\Phi_{E},\Delta\Phi_{H})=(60,30)$ degrees, E-plane
$(\Delta\Phi_{E},\Delta\Phi_{H})=(30,60)$ degrees. (Fourth row) $f=1.0$ GHz,
with (from left to right) H-plane $(\Delta\Phi_{E},\Delta\Phi_{H})=(0,0)$
degrees, H-plane $(\Delta\Phi_{E},\Delta\Phi_{H})=(60,30)$ degrees, H-plane
$(\Delta\Phi_{E},\Delta\Phi_{H})=(30,60)$ degrees.

The radiation patterns in the E and H plane for the two-dimensional horn array
are shown in Figure 12. The phase combinations
$(\Delta\Phi_{E},\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 $\leq-15$
dB. The side lobes are also at the $\approx-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 $(\Delta\Phi_{E},\Delta\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 Mailloux (2018)). 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(z)$ is interesting to
the ultra-high energy (UHE) neutrino community Vieregg et al. (2016); Avva et
al. (2017). 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 Hanson (2011); Hanson et al.
(2015); Avva et al. (2014). The index of refraction varies within the ice
because of the transition between surface snow ($\rho\approx 0.4$ g/cm3) and
the solid ice below ($\rho=0.917$ g/cm3). Most recent and intricate studies of
phased array beam behavior still assume a uniform medium Adamidis et al.
(2020); Ahn et al. (2019); Gampala and Reddy (2016). Embedded phased arrays
with varying $n(z)$ emit signals that curve in the direction of increasing
$n(z)$.

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 Dookayka (2011). While
there is evidence that RF signals can propagate horizontally through Antarctic
ice Barwick et al. (2018), data from Greenland suggests the relative strength
of the effect is small compared to the curved radiation Deaconu et al. (2018).
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(z)$ data versus depth $z$ below the surface is
given by Barwick et al. (2018)

$n(z)=\begin{cases}1&z>0\\\ n_{\rm ice}-\Delta n\exp(z/z_{0})&z\leq
0\end{cases}$ (28)

The fit parameters in Equation 28 come from Reference Barwick et al. (2018):
$\Delta n=0.423\pm 0.004$ and $z_{0}=77\pm 2$ meters, with $n_{\rm 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 Barwick et al. (2018)).
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 meters, and the spacing is
0.75 meters.

Figure 13:  The simplified $N=8$ one-dimensional vertical phased array with
dimensions in meters. For this array, $d_{z}=\lambda/2.0$, and the length of
the dipole radiators is $\lambda/4.0$. The colorscale represents $n(z)$ in
Equation 28.

Figures 13 and 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 $\Delta z$ is 0.1 meters. 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 meters (top) and 35 meters
(bottom).

Figure 14:  The magnitude of the z-component of the z-polarized dipoles as
they radiate as a phased array with $\Delta\Phi=0$ degrees. The air, surface,
and ice regions are the same as Figure 13, with the same dimensions. (Top) The
array depth is -15 meters. (Bottom) The array depth is -35 meters. (A)
Radiation refracting through the surface into the air. (B) Grating lobes
reflect from the surface into the shadow zone. (C) Grating lobes propagating
through the shadow zone. (D) The main beam bent downward due to the gradient
in $n(z)$.

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(z)$, 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(z)$ curvature
is weakened. The beam travels farther horizontally because the gradient of
$n(z)$ 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 $\leq 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 $\Delta\phi$ and
$\Delta\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 $\Delta\phi$ and
$\Delta\Phi$ (Figure 4). In that case, the slope of $\Delta\phi$ vs.
$\Delta\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 Mailloux
(2018). At lower frequencies, however, scan loss takes a toll on radiated
power (Tab. 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
$\Delta\phi_{E}$ and $\Delta\Phi_{E}$, and $\Delta\phi_{H}$ and
$\Delta\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
$\Delta\phi_{E}$ and $\Delta\Phi_{E}$, and $\Delta\phi_{H}$ and
$\Delta\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\rightarrow\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(z)$. The model for
$n(z)$ 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 Aguilar et al. (2020). 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 trustworthiness across a wide bandwidth.

This research was funded by the Office of Naval Research (ONR) under the
Summer Faculty Research Program (SFRP).

Not applicable. Not applicable. The data presented in this study are available
on request from the corresponding author. Due to large file sizes and
restricted bandwidth, please contact corresponding author to set up
collaborative sharing.

###### Acknowledgements.

We would like to thank the Office of Naval Research (ONR) for helping to
support this research. In particular, we would like to thank the Naval Surface
Warfare Center Corona Division and their continued support of the ONR Summer
Faculty Research Program (SFRP). Conversations with Christopher Clark and Gary
Yeakley were especially helpful. We are also grateful to Van Nguyen, Jeffery
Benson, and Golda McWhorter. Karon Myles deserves our special thanks for
helping to coordinate the SFRP program. The authors declare no conflict of
interest. no

## Appendix A

MEEP FDTD is more often applied to $1\mu$m scale lengths than the 1 cm-scale
RF elements, so scale-invariance must be highlighted. Scale invariant units
with $c=1$ are used in MEEP when solving Maxwell’s equations with the FDTD
technique. The typical length scale in MEEP analysis is usually called the
“a-value.” Systems with dimensions of order $1\mu$m are said to have an
“a-value of $1\mu$m.” A value of $a=1$ cm is chosen for models presented in
Secs. 3 and 4. For example, if the frequency $f=7.5$ GHz, the wavelength is
$4.0$ cm $=4.0a$, or simply $4.0$ with $f=1/\lambda=0.25$. Since $c=1$,
$\lambda=f^{-1}=T$. The period is $4.0$, so a simulation run of $50T=200$ time
units corresponds to 6.67 ns. Assuming 1 pixel/a-value, simulated radiation
would therefore propagate 200 units of $\Delta x$ in a straight line before
time was up. A resolution parameter sets the number of pixels per distance
unit and is usually larger than 1.0. Selecting the right resolution is often a
subtle balance between capturing the most relevant effects while limiting the
memory usage of the simulation results.

References

yes

## References

  * Syrytsin et al. (2017) Syrytsin, I.; Zhang, S.; Pedersen, G.F.  Circularly Polarized Planar Helix Phased Antenna Array for 5G Mobile Terminals.  2017 International Conference on Electromagnetics in Advanced Applications (ICEAA), Verona, Italy, September 2017 2017, pp. 1105–1108.  doi:black10.1109/iceaa.2017.8065458.
  * Kikuchi et al. (2017) Kikuchi, K.; Mikada, H.; Takekawa, J.  Improved Imaging Capability of Phased Array Antenna in Ground Penetrating Radar Survey.  Conference Proceedings, 79th EAGE Conference and Exhibition 2017, Paris, France, June 2017 2017.  doi:black10.3997/2214-4609.201701121.
  * Vieregg et al. (2016) Vieregg, A.; Bechtol, K.; Romero-Wolf, A.  A technique for detection of PeV neutrinos using a phased radio array.  Journal of Cosmology and Astroparticle Physics 2016, 2016, 005.  doi:black10.1088/1475-7516/2016/02/005.
  * Munekata et al. (2014) Munekata, T.; Yamamoto, M.; Nojima, T.  A Wideband 16-Element Antenna Array Using Leaf-Shaped Bowtie Antenna and Series-Parallel Feed Networks.  2014 IEEE International Workshop on Electromagnetics (iWEM), Sapporo Hokkaido, Japan, August 2014 2014, pp. 80–81.  doi:black10.1109/iwem.2014.6963645.
  * Avva et al. (2017) Avva, J.; Bechtol, K.; Chesebro, T.; Cremonesi, L.; Deaconu, C.; Gupta, A.; Ludwig, A.; Messino, W.; Miki, C.; Nichol, R.; Oberla, E.; Ransom, M.; Romero-Wolf, A.; Saltzberg, D.; Schlupf, C.; Shipp, N.; Varner, G.; Vieregg, A.; Wissel, S.  Development toward a ground-based interferometric phased array for radio detection of high energy neutrinos.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 2017, 869, 46 – 55.  doi:blackhttps://doi.org/10.1016/j.nima.2017.07.009.
  * Ansys, Inc., Canonsburg, Pennsylvania (2020) Ansys, Inc., Canonsburg, Pennsylvania.  3D Electromagnetic Field Simulator for RF and Wireless Design, 2020\.
  * Feng et al. (2019) Feng, N.; Zhang, Y.; Tian, X.; Zhu, J.; Joines, W.T.; Wang, G.P.  System-Combined ADI-FDTD Method and Its Electromagnetic Applications in Microwave Circuits and Antennas.  IEEE Transactions on Microwave Theory and Techniques 2019, 67, 3260–3270.  doi:black10.1109/tmtt.2019.2919838.
  * Zhu et al. (2017) Zhu, L.; Hwang, H.S.; Ren, E.; Yang, G.  High Performance MIMO Antenna for 5G Wearable Devices.  2017 IEEE International Symposium on Antennas and Propagation & USNC/URSI National Radio Science Meeting, San Diego, CA, USA, July 2017 2017, pp. 1869–1870.  doi:black10.1109/apusncursinrsm.2017.8072977.
  * (9) Ho, T.Q.; Hunt, L.N.; Hewett, C.A.; Ready, T.G.; Mittra, R.; Yu, W.; Zolnick, D.A.; Kragalott, M.  Analysis of Electrically Large Patch Phased Arrays via CFDTD.  2006 IEEE Antennas and Propagation Society International Symposium, Albuquerque, New Mexico, July 2006.
  * Burke et al. (2004) Burke, G.J.; Miller, E.K.; Poggio, A.J.  The Numerical Electromagnetics Code (NEC) - A Brief History.  IEEE Antennas and Propagation Society Symposium, 2004 2004, 3, 2871–2874.  doi:black10.1109/aps.2004.1331976.
  * Oskooi et al. (2010) Oskooi, A.F.; Roundy, D.; Ibanescu, M.; Bermel, P.; Joannopoulos, J.; Johnson, S.G.  Meep: A flexible free-software package for electromagnetic simulations by the FDTD method.  Computer Physics Communications 2010, 181, 687–702.  doi:black10.1016/j.cpc.2009.11.008.
  * Fedeli et al. (2019) Fedeli, A.; Montecucco, C.; Gragnani, G.L.  Open-Source Software for Electromagnetic Scattering Simulation: The Case of Antenna Design.  Electronics 2019, 8, 1506.  doi:black10.3390/electronics8121506.
  * Liebig et al. (2013) Liebig, T.; Rennings, A.; Held, S.; Erni, D.  openEMS – a free and open source equivalent-circuit (EC) FDTD simulation platform supporting cylindrical coordinates suitable for the analysis of traveling wave MRI applications.  International Journal of Numerical Modelling: Electronic Networks, Devices and Fields 2013, 26, 680–696.  doi:black10.1002/jnm.1875.
  * Warren et al. (2016) Warren, C.; Giannopoulos, A.; Giannakis, I.  gprMax: Open source software to simulate electromagnetic wave propagation for Ground Penetrating Radar.  Computer Physics Communications 2016, 209, 163–170.  doi:black10.1016/j.cpc.2016.08.020.
  * Richie and Ababei (2017) Richie, J.E.; Ababei, C.  Optimization of patch antennas via multithreaded simulated annealing based design exploration.  Journal of Computational Design and Engineering 2017, 4, 249–255.  doi:black10.1016/j.jcde.2017.06.004.
  * Mailloux (2018) Mailloux, R.  The Phased Array Antenna Handbook; Artech House: Norwood, MA, 2018\.
  * Wahl et al. (2013) Wahl, P.; Gagnon, D.S.L.; Debaes, C.; Erps, J.V.; Vermeulen, N.; Miller, D.A.B.; Thienpont, H.  B-CALM: An Open-Source Multi-GPU-based 3D-FDTD with Mutli-pole dispersion for Plasmonics.  Progress In Electromagnetics Research 2013, 138, 467–478.  doi:black10.2528/pier13030606.
  * Gorham et al. (2009) Gorham, P.; Allison, P.; Barwick, S.; Beatty, J.; Besson, D.; Binns, W.; Chen, C.; Chen, P.; Clem, J.; Connolly, A.; Dowkontt, P.; DuVernois, M.; Field, R.; Goldstein, D.; Goodhue, A.; Hast, C.; Hebert, C.; Hoover, S.; Israel, M.; Kowalski, J.; Learned, J.; Liewer, K.; Link, J.; Lusczek, E.; Matsuno, S.; Mercurio, B.; Miki, C.; Miočinović, P.; Nam, J.; Naudet, C.; Nichol, R.; Palladino, K.; Reil, K.; Romero-Wolf, A.; Rosen, M.; Ruckman, L.; Saltzberg, D.; Seckel, D.; Varner, G.; Walz, D.; Wang, Y.; Williams, C.; Wu, F.  The Antarctic Impulsive Transient Antenna ultra-high energy neutrino detector: Design, performance, and sensitivity for the 2006–2007 balloon flight.  Astroparticle Physics 2009, 32, 10–41, [0812.1920].  doi:black10.1016/j.astropartphys.2009.05.003.
  * Gorham et al. (2007) Gorham, P.; Barwick, S.; Beatty, J.; Besson, D.; Binns, W.; Chen, C.; Chen, P.; Clem, J.; Connolly, A.; Dowkontt, P.  Observations of the Askaryan effect in ice.  Physical review letters 2007, 99, 171101, [hep-ex/0611008].
  * Barwick et al. (2015) Barwick, S.; Berg, E.; Besson, D.; Duffin, T.; Hanson, J.; Klein, S.; Kleinfelder, S.; Piasecki, M.; Ratzlaff, K.; Reed, C.; Roumi, M.; Stezelberger, T.; Tatar, J.; Walker, J.; Young, R.; Zou, L.  Time-domain response of the ARIANNA detector.  Astroparticle Physics 2015, 62, 139–151, [1406.0820].  doi:black10.1016/j.astropartphys.2014.09.002.
  * Hanson (2011) Hanson, J.C.  Ross Ice Shelf Thickness, Radio-frequency Attenuation and Reflectivity: Implications for the ARIANNA UHE Neutrino Detector.  In Proceedings of the 32nd International Cosmic Ray Conference, Beijing, China, August 2011 2011.  doi:black10.7529/ICRC2011/V04/0340.
  * Hanson et al. (2015) Hanson, J.C.; Barwick, S.W.; Berg, E.C.; Besson, D.Z.; Duffin, T.J.; Klein, S.R.; Kleinfelder, S.A.; Reed, C.; Roumi, M.; Stezelberger, T.; Tatar, J.; Walker, J.A.; Zou, L.  Radar absorption, basal reflection, thickness and polarization measurements from the Ross Ice Shelf, Antarctica.  Journal of Glaciology 2015, 61, 438–446.  doi:black10.3189/2015jog14j214.
  * Avva et al. (2014) Avva, J.; Kovac, J.; Miki, C.; Saltzberg, D.; Vieregg, A.  An in situ measurement of the radio-frequency attenuation in ice at Summit Station, Greenland.  Journal of Glaciology 2014, [1409.5413].  doi:black10.3189/2015jog15j057.
  * Adamidis et al. (2020) Adamidis, G.A.; Vardiambasis, I.O.; Ioannidou, M.P.; Kapetanakis, T.N.  Design and implementation of an adaptive beamformer for phased array antenna applications.  Microwave and Optical Technology Letters 2020, 62, 1780–1784.  doi:black10.1002/mop.32231.
  * Ahn et al. (2019) Ahn, B.; Hwang, I.J.; Kim, K.S.; Chae, S.C.; Yu, J.W.; Lee, H.L.  Wide-Angle Scanning Phased Array Antenna using High Gain Pattern Reconfigurable Antenna Elements.  Scientific Reports 2019, 9, 18391.  doi:black10.1038/s41598-019-54120-2.
  * Gampala and Reddy (2016) Gampala, G.; Reddy, C.J.  Advanced Computational Tools for Phased Array Antenna Applications.  2016 IEEE International Symposium on Phased Array Systems and Technology (PAST), Waltham, MA, USA, October 2016 2016, pp. 1–5.  doi:black10.1109/array.2016.7832612.
  * Dookayka (2011) Dookayka, K.  Characterizing the Search for Ultra-High Energy Neutrinos with the ARIANNA Detector.  PhD thesis, Univeristy of California at Irvine, 2011.
  * Barwick et al. (2018) Barwick, S.W.; Berg, E.C.; Besson, D.Z.; Gaswint, G.; Glaser, C.; Hallgren, A.; Hanson, J.C.; Klein, S.R.; Kleinfelder, S.; Köpke, L.; Kravchenko, I.; Lahmann, R.; Latif, U.; Nam, J.; Nelles, A.; Persichilli, C.; Sandstrom, P.; Tatar, J.; Unger, E.  Observation of classically ‘forbidden’ electromagnetic wave propagation and implications for neutrino detection.  Journal of Cosmology and Astroparticle Physics 2018, 2018, 055–055.  doi:black10.1088/1475-7516/2018/07/055.
  * Deaconu et al. (2018) Deaconu, C.; Vieregg, A.G.; Wissel, S.A.; Bowen, J.; Chipman, S.; Gupta, A.; Miki, C.; Nichol, R.J.; Saltzberg, D.  Measurements and Modeling of Near-Surface Radio Propagation in Glacial Ice and Implications for Neutrino Experiments.  Phys. Rev. D 2018, 98, 043010.  doi:black10.1103/PhysRevD.98.043010.
  * Aguilar et al. (2020) Aguilar, J.A.; Allison, P.; Beatty, J.J.; Bernhoff, H.; Besson, D.; Bingefors, N.; Botner, O.; Buitink, S.; Carter, K.; Clark, B.A.; Connolly, A.; Dasgupta, P.; de Kockere, S.; de Vries, K.D.; Deaconu, C.; DuVernois, M.A.; Feigl, N.; Garcia-Fernandez, D.; Glaser, C.; Hallgren, A.; Hallmann, S.; Hanson, J.C.; Hendricks, B.; Hokanson-Fasig, B.; Hornhuber, C.; Hughes, K.; Karle, A.; Kelley, J.L.; Klein, S.R.; Krebs, R.; Lahmann, R.; Magnuson, M.; Meures, T.; Meyers, Z.S.; Nelles, A.; Novikov, A.; Oberla, E.; Oeyen, B.; Pandya, H.; Plaisier, I.; Pyras, L.; Ryckbosch, D.; Scholten, O.; Seckel, D.; Smith, D.; Southall, D.; Torres, J.; Toscano, S.; Van Den Broeck, D.J.; van Eijndhoven, N.; Vieregg, A.G.; Welling, C.; Wissel, S.; Young, R.; Zink, A.  Design and Sensitivity of the Radio Neutrino Observatory in Greenland (RNO-G).  arXiv e-prints 2020, p. arXiv:2010.12279, [arXiv:astro-ph.IM/2010.12279].


# Web accessibility trends and implementation in dynamic web applications

Timothy W. Hostetler1, Shinyi Chen, Sergi Blanco-Cuaresma, Alberto Accomazzi,
Michael J. Kurtz, Carolyn S. Grant, Edwin Henneken, Donna M. Thompson, Roman
Chyla, Golnaz Shapurian, Matthew R. Templeton, Kelly E. Lockhart, Nemanja
Martinovic, Stephen McDonald, Felix Grezes

###### Abstract

The NASA Astrophysics Data System (ADS), a critical research service for the
astrophysics community, strives to provide the most accessible and inclusive
environment for the discovery and exploration of the astronomical literature.
Part of this goal involves creating a digital platform that can accommodate
everybody, including those with disabilities that would benefit from
alternative ways to present the information provided by the website. NASA ADS
follows the official Web Content Accessibility Guidelines (WCAG) standard for
ensuring accessibility of all its applications, striving to exceed this
standard where possible. Through the use of both internal audits and external
expert review based on these guidelines, we have identified many areas for
improving accessibility in our current web application, and have implemented a
number of updates to the UI as a result of this. We present an overview of
some current web accessibility trends, discuss our experience incorporating
these trends in our web application, and discuss the lessons learned and
recommendations for future projects.

1Center for Astrophysics | Harvard & Smithsonian, Cambridge, MA, USA<EMAIL_ADDRESS>
##  1 Introduction

Every year the WebAIM(IDRPP 2021), a nonprofit organization specializing in
accessible web solutions holds a survey of the top 1 million websites using
their evaluation tool called WAVE. The organization has found that even in
some of the most popular home pages on the web, there were on average 51.4
accessibility errors per page. Many accessibility errors are easily fixed by
conforming to web standards, using good web page structure, and following
guidelines like those found in the Web Content Accessibility
Guidelines111https://www.w3.org/WAI/standards-guidelines/wcag/ (WCAG) provided
by W3C Web Accessibility Initiative222https://www.w3.org/WAI/ (WAI). The NASA
Astrophysics Data System333https://ui.adsabs.harvard.edu (ADS; Kurtz et al.
2000) web site is a modern single-page application, which presents a unique
challenge for web accessibility due to its dynamic content. Utilizing basic
techniques and some special attributes, we can provide accessible content even
when the information on the page is changing over time.

##  2 Web Content Accessibility Guidelines

WCAG provides a set of guidelines that cover a large range of web content and
components. The current standard is WCAG 2.1, which is further broken down
into three levels of conformance: A (lowest), AA, AAA (highest). There are 4
main sections which are defined: Perceivable, Operable, Understandable and
Robust.

###  2.1 Perceivable

Web developers should ensure in their apps, if possible, that there is at
least one alternative way for users to get access to content. Non-text
elements like images should utilize alt tags, and videos should have
transcripts and/or audio descriptions. HTML standard components like lists and
tables should have proper structure and predictability, allowing users
utilizing tools to perceive content on the page to have an easier time doing
so. In general, any visual changes to the site should not cause a loss of
information.

###  2.2 Operable

All content that is accessible by mouse should also be accessible using a
keyboard. Content does not induce seizures or physical reactions of any kind
and users are given enough time to perform any actions. Make sure the page
layout and structure makes finding and getting to content as easy as possible.
All potential input modalities should be taken into account, such as touch and
voice.

###  2.3 Understandable

Text is readable and understandable, and presented in a clear and predictable
way. Users are helped, where possible, to correct and avoid mistakes.

###  2.4 Robust

Application markup should be well-structured and reliably interpreted in order
to maximise compatibility with existing and future browsers, assistive
technologies and any other type of user agent.

##  3 Web Accessibility Basics

###  3.1 Page Structure

Semantic HTML provides structural elements like <nav>, <header>, <main>, and
<footer>, these are known as page regions. Some or all of them should be used
to layout a page, since they allow screen readers and other tools to infer
some information about the content within. There are also other elements, like
<section> and <article> which should be utilized as much as possible, to
provide more context to users.

According to a screen reader survey performed by WebAIM, users have relied
less and less on page regions to navigate pages. Instead, nearly 70% of
respondents indicated that they prefer to use headings to find information on
the page. This is an important insight, since headings (or landmarks) provide
a more flexible way of laying out content on the page than the more rigid
region tags. Web apps should therefore be utilizing heading tags to provide
"jump" points for screen reader users.

ADS’s web app uses heading tags extensively to provide structure and layout.
Screen readers will follow the hierarchy of headings to jump to areas of
interest or to list the titles of search result, for example.

###  3.2 Text Alternatives

Non-text content like images and icons need to have some kind of textual
description. This can be accomplished using semantic tags like <figure> and
<figcaption> and in the case of images also using the alt attribute. Some
visually hidden inline text can be added next to the icon or image to provide
some extra information to screen readers.

Aria attributes are a set of accessibility-focused attributes which supplement
HTML to provide extra information, especially when complex structures are
involved or when creating dynamic content. The use of aria attributes like
aria-labelledby can be used to refer to a section title or some other element
for example.

##  4 Dynamic Content

One of the challenges with modern web applications is the need for handling
page transitions and information change in an accessible way. In ADS, almost
every page has some kind of dynamic content which is changing in some way over
time. In order to present changing information to screen readers, we can
utilize some special attributes. Using the role="status" or role="alert"
attribute we can inform screen readers that this element has information in it
that will change over time. It will automatically announce to the user when
something within this element changes. Another attribute aria-atomic="true"
can be provided to have the tool read the entire element text again, not just
what changed. This can be quite useful for status messages and error messages,
for example.

##  5 Tools

During development on the web, it is useful to regularly be testing that your
WCAG level of conformance is being maintained. Some static analysis tools like
WAVE and axe DevTools are browser extensions which can be used to directly
analyze the website. Screen readers should also be employed as often as
possible to get practice with use and to provide some limited real-world
experience, without the use of actual users. Some common screen readers are
JAWS, NVDA, and VoiceOver.

##  6 Accessibility Audits

ADS contracted a user experience and accessibility consultancy based at
Harvard University called Harvard Web Publishing (HWP). HWP’s in-house
accessibility experts utilized VoiceOver and NVDA screen readers and provided
us with a detailed report of several ADS web pages. Within the provided
report, we found many small areas of improvement, which were fixed.

One benefit to having an external audit like this is to provide us with
insight into how a particular user can be confounded by a lack of information
or context. By not providing appropriate semantic elements or not labelling
form elements properly, screen reader users were unable to understand the
purpose of some controls. Also, how we manage focus needed to improve, with
some areas inaccessible via keyboard.

##  7 Conclusion

Web accessibility will always be an ongoing effort, with some areas never
being able to fully conform to WAI standards. However, all sites should
conform at least to the minimal accessibility level (A). At ADS, we continue
to improve accessibility in our web application and any future enhancements of
our user interface will have the highest conformity as a core goal.

#### Acknowledgments

The NASA Astrophysics Data System is operated by the Smithsonian Astrophysical
Observatory under NASA Cooperative Agreement 80NSSC21M0056.

## References

  * IDRPP (2021) IDRPP 2021, Webaim: Screen reader user survey #9 results. URL https://webaim.org/projects/screenreadersurvey9
  * Kurtz et al. (2000) Kurtz, M. J., Eichhorn, G., Accomazzi, A., Grant, C. S., Murray, S. S., & Watson, J. M. 2000, Astronomy and Astrophysics Supplement Series, 143, 41. astro-ph/0002104