# Vector Approximate Message Passing with Arbitrary i.i.d. Noise Priors

###### Abstract

Approximate message passing (AMP) algorithms are devised under the Gaussianity
assumption of the measurement noise vector. In this work, we relax this
assumption within the vector AMP (VAMP) framework to arbitrary independent and
identically distributed (i.i.d.) noise priors. We do so by rederiving the
linear minimum mean square error (LMMSE) to accommodate both the noise and
signal estimations within the message passing steps of VAMP. Numerical results
demonstrate how our proposed algorithm handles non-Gaussian noise models as
compared to VAMP. This extension to general noise priors enables the use of
AMP algorithms in a wider range of engineering applications where non-Gaussian
noise models are more appropriate.

Index Terms—  Approximate message passing, expectation propagation, non-
Gaussian noise, inference algorithms.

##  1 Introduction

###  1.1 Background and related work

One of the fundamental problems in the compressed sensing (CS) literature is
recovering a sparse vector $\bm{x}\in\mathbb{R}^{N}$ from an observation
vector, $\bm{y}\in\mathbb{R}^{M}$, obtained from a noisy linear measurement of
the form

$\bm{y}=\bm{A}\,\bm{x}+\bm{w}.$ (1)

Here, $\bm{A}\in\mathbb{R}^{M\times N}$ is the measurement matrix with $M<N$,
and $\bm{w}$ is unstructured additive noise vector. Solving the inverse
problem in (1) is commonly achieved under the assumption that $\bm{w}$ is an
additive white Gaussian noise (AWGN) vector whose entries are assumed to be
mutually independent with mean zero and variance $\gamma_{w}$, i.e.,
$w_{i}\sim\mathcal{N}(\mathsf{w}_{i};0,\gamma_{w}^{-1})$. Since the
introduction of the approximate message passing (AMP) algorithm [1], a lot of
interest has been devoted to devising algorithmic solutions to (1) within the
CS framework. Specifically, the generalized AMP (GAMP) [2] algorithm extends
over its AMP predecessor by $i)$ accommodating statistical priors on the
sparse vector $\bm{x}$ and $ii)$ handling nonlinear output transformations.
The divergence issues of the GAMP algorithm due to ill-conditioning of the
measurement matrix $\bm{A}$ has been alleviated by the introduction of the
vector AMP (VAMP) [3] 111Note that VAMP was independently derived in [4] by
another research group under the name orthogonal AMP (OAMP) algorithm..

$\bm{x}^{+}$$\bm{x}^{-}$$\bm{w}^{-}$$\bm{w}^{+}$$\bm{y}$$p_{\bm{\mathsf{x}}^{+}}(\bm{x}^{+})$$\delta(\bm{x}^{+}-\bm{x}^{-})$$\underbrace{\delta(\bm{y}-\bm{A}\,\bm{x}^{-}-\bm{w}^{-})}_{\triangleq~{}f_{\bm{\mathsf{y}}|\bm{\mathsf{x}}^{-},\bm{\mathsf{w}}^{-}}}$$\delta(\bm{w}^{+}-\bm{w}^{-})$$p_{\bm{\mathsf{w}}^{+}}(\bm{w}^{+})$$\mu_{f\xrightarrow{\,}\bm{x}^{-}}$$\mu_{f\xrightarrow{\,}\bm{w}^{-}}$$\mu_{\delta_{\text{x}}\xrightarrow{\,}\bm{x}^{-}}$$\mu_{\delta_{\text{w}}\xrightarrow{\,}\bm{w}^{-}}$
Fig. 1: Factor graph of the VAMP algorithm for arbitrary i.i.d. noise priors.
Circles represent variable nodes and squares represent factor nodes.

When the noise vector $\bm{w}$ is Gaussian, VAMP strikes a proper balance
between the reconstruction performance and the computational complexity as
compared to traditional convex optimization-based and iterative soft-
thresholding algorithms [5]. Despite the immense popularity of Gaussian noise
models, estimators derived with a Gaussian assumption are sensitive to
outliers [6]. In practice, many engineering applications (e.g., electronic
devices, lasers, relay switching) and natural phenomena (e.g., atmospheric
noise, lightning spikes, and ice cracking) are more accurately characterized
by heavy-tailed non-Gaussian measurement noise models [7]. To sidestep the
requirement for Gaussian noise models, this paper extends the VAMP algorithm
in the presence of arbitrary i.i.d. noise priors, thereby handling a broader
class of practical applications.

###  1.2 Paper organization and notations

We structure the rest of this paper as follows. In Section 2, we present how
we incorporate the non-Gaussian noise prior to the expectation propagation
steps of the VAMP algorithm. Then, we rederive the new linear minimum mean
square error (LMMSE) estimation step of our proposed algorithm to handle
arbitrary i.i.d. noise priors. In Section 3, we discuss the simulation results
of the algorithm. Finally, we draw out some concluding remarks in Section 4.

We also mention the common notations used in this paper. We use Sans Serif
fonts (e.g., $\mathsf{x}$) for random variables and Serif fonts (e.g., $x$)
for their realizations. We use boldface lowercase letters for vectors (e.g.,
$\bm{\mathsf{x}}$ and $\bm{x}$) and we denote the $i$th component of $\bm{x}$
as $x_{i}$. The operator diag($\bm{X}$) stacks the diagonal elements of a
matrix $\bm{X}$ in a vector, $\bm{I}$ stands for the identity matrix, $\bm{1}$
denotes the all-ones vector. The notation $\bm{\mathsf{x}}\sim$
$p_{\bm{\mathsf{\bm{x}}}}(\bm{x};\bm{\theta})$ means that $\bm{\mathsf{x}}$ is
distributed according to the pdf $p_{\bm{\mathsf{x}}}(\bm{x},\bm{\theta})$
which is parameterized by a parameter vector $\bm{\theta}$. Moreover,
$\mathcal{N}(\bm{x};\widehat{\bm{x}},\bm{R})$ stands for the multivariate
Gaussian pdf with mean $\widehat{\bm{x}}$ and covariance matrix $\bm{R}$. We
also use $\mathbb{E}[\bm{\mathsf{x}}|d(\bm{x})]$ and
$\text{Cov}[\bm{\mathsf{x}}|d(\bm{x})]$ to denote the expectation and the
covariance of $\bm{{\mathsf{x}}}\sim d(\bm{x})$. Finally, $\delta(\bm{x})$
denotes the Dirac delta distribution and $\langle\bm{x}\rangle\triangleq$
$\frac{1}{N}\sum_{i=1}^{N}x_{i}$ for $\bm{x}\in\mathbb{R}^{N}$.

##  2 VAMP with arbitrary i.i.d. noise priors

Before delving into the derivation details, we first introduce the VAMP
algorithm for arbitrary i.i.d. noise priors which run iteratively according to
the algorithmic steps of Algorithm 1. There, $t$ stands for the iteration
index, and subscripts $\mathsf{p}$ and $\mathsf{e}$ are used to distinguish
“posterior” and “extrinsic” variables, respectively. As a visual reminder, we
also use the hat symbol “$~{}\,\widehat{}~{}\,$” to refer to mean values.

To derive the VAMP algorithm with arbitrary i.i.d. noise priors, we start by
factoring the joint pdf of all the observed and unobserved variables in (1) as
follows:

$\displaystyle p(\bm{y},\bm{x},\bm{w})$
$\displaystyle=p(\bm{y}|\bm{x},\bm{w})\,p(\bm{x})\,p(\bm{w})$ (2)
$\displaystyle=\delta(\bm{y}-\bm{A}\,\bm{x}-\bm{w})\,p(\bm{x})\,p(\bm{w}).$

Then, we split both $\bm{x}$ and $\bm{w}$ into two identical variables, i.e.,
$\bm{x}^{+}=\bm{x}^{-}$ and $\bm{w}^{+}=\bm{w}^{-}$, thereby transforming the
joint pdf in (2) to the equivalent factorization:

$\displaystyle\small{\\!\\!p(\bm{y},\bm{x}^{+},\bm{x}^{-},\bm{w}^{+},\bm{w}^{-})}$
(3) $\displaystyle\hskip
28.45274pt\small{~{}=~{}p(\bm{x}^{+})\,p(\bm{w}^{+})\,\delta(\bm{x}^{+}-\bm{x}^{-})\times}$
$\displaystyle\hskip
56.9055pt\small{\delta(\bm{w}^{+}-\bm{w}^{-})\,\delta(\bm{y}-\bm{A}\,\bm{x}^{-}-\bm{w}^{-}),}$

whose associated factor graph is depicted in Fig. 1.

Similarly to the derivation of standard VAMP, we follow an expectation
propagation (EP)-like approximation of the sum-product (SP) belief propagation
(BP) algorithm based on the two following rules:

  * •

EP approximation: Given a variable node $\bm{\mathsf{x}}^{-}$ receiving the
message $\mu_{f\rightarrow\bm{\mathsf{x}}^{-}}$ from the factor node
$f_{\bm{\mathsf{y}}|\bm{\mathsf{x}}^{-},\bm{\mathsf{w}}^{-}}=\delta(\bm{y}-\bm{A}\,\bm{x}^{-}-\bm{w}^{-})$,
EP approximates the SP belief $b_{\text{sp}}(\bm{x}^{-})$ with the approximate
belief
$b_{\text{app}}(\bm{x}^{-})=\mathcal{N}(\bm{x}^{-};\widehat{\bm{x}}_{\mathsf{p}}^{-},\gamma_{\bm{x}_{\mathsf{p}}^{-}}^{-1})$,
where $\widehat{\bm{x}}_{\mathsf{p}}^{-}=\mathbb{E}[\bm{x}^{-}|b_{\text{sp}}]$
and
$\gamma_{\bm{x}_{\mathsf{p}}^{-}}^{-1}=\left\langle\operatorname{diag}\left(\operatorname{Cov}\left[\bm{x}^{-}|b_{\text{sp}}\right]\right)\right\rangle$
are the mean and average variance of the SP belief
$b_{\text{sp}}(\bm{x}^{-})$.

  * •

Extrinsic belief computation: As shown in Fig 2, given a posterior Gaussian
belief
$\mathcal{N}(\bm{x}^{-};\widehat{\bm{x}}_{\mathsf{p}}^{-},\gamma_{\bm{x}_{\mathsf{p}}^{-}}^{-1}\,\bm{I})$
and an incoming extrinsic Gaussian belief
$\mathcal{N}(\bm{x}^{-};\widehat{\bm{x}}_{\mathsf{e}}^{+},\gamma_{\bm{x}_{\mathsf{e}}^{+}}^{-1}\,\bm{I})$,
the extrinsic mean and precision of the outcoming Gaussian belief
$\mathcal{N}(\bm{x}^{+};\widehat{\bm{x}}_{\mathsf{e}}^{-},\gamma_{\bm{x}_{\mathsf{e}}^{-}}^{-1}\,\bm{I})$
are computed as follows:

$\displaystyle\\!\\!\\!\\!\gamma_{\bm{x}_{\mathsf{e}}^{-}}$
$\displaystyle\,\,=\,\,\gamma_{\bm{x}_{\mathsf{p}}^{-}}-\gamma_{\bm{x}_{\mathsf{e}}^{+}}~{}\triangleq~{}\text{EXT}\left(\gamma_{\bm{x}_{\mathsf{p}}^{-}}\right),$
(4a) $\displaystyle\\!\\!\\!\\!\widehat{\bm{x}}_{\mathsf{e}}^{-}$
$\displaystyle\,\,=\,\,\gamma_{\bm{x}_{\mathsf{e}}^{+}}^{-1}\left(\gamma_{\bm{x}_{\mathsf{p}}^{-}}\,\widehat{\bm{x}}_{\mathsf{p}}^{-}-\gamma_{\bm{x}_{\mathsf{e}}^{+}}\,\widehat{\bm{x}}_{\mathsf{e}}^{+}\right)~{}\triangleq~{}\text{EXT}\left(\widehat{\bm{x}}_{\mathsf{p}}^{-}\right).$
(4b)

$p_{\bm{\mathsf{x}}}(\bm{x})$ $\bm{y}=\bm{A}\,\bm{x}+\bm{w}$
$p_{\bm{\mathsf{w}}}(\bm{w})$EXTEXTEXTEXT$\widehat{\bm{x}}_{\mathsf{p}}^{+}$$\gamma_{\bm{x}_{\mathsf{p}}^{+}}$$\widehat{\bm{x}}_{\mathsf{e}}^{+}$$\gamma_{\bm{x}_{\mathsf{e}}^{+}}$$\widehat{\bm{w}}_{\mathsf{p}}^{-}$$\gamma_{\bm{w}_{\mathsf{p}}^{-}}$$\widehat{\bm{w}}_{\mathsf{p}}^{-}$$\gamma_{\bm{w}_{\mathsf{e}}^{-}}$$\widehat{\bm{w}}_{\mathsf{p}}^{+}$$\gamma_{\bm{w}_{\mathsf{p}}^{+}}$$\widehat{\bm{w}}_{\mathsf{e}}^{+}$$\gamma_{\bm{w}_{\mathsf{e}}^{+}}$$\widehat{\bm{x}}_{\mathsf{p}}^{-}$$\gamma_{\bm{x}_{\mathsf{p}}^{-}}$$\widehat{\bm{x}}_{\mathsf{e}}^{-}$$\gamma_{\bm{x}_{\mathsf{e}}^{-}}$
Fig. 2: Block diagram of VAMP for arbitrary i.i.d. noise priors with its three
modules: two denoising MMSE modules incorporating the prior information,
$p_{\bm{\mathsf{x}}}(\cdot)$ and $p_{\bm{\mathsf{w}}}(\cdot)$, and the LMMSE
module. The three modules exchange extrinsic information/messages through the
ext blocks. The color of each module matches the color of the corresponding
line numbers in Algorithm 1.

Unlike the factor graph of standard VAMP, however, the two variable nodes
$\bm{x}^{-}$ and $\bm{w}^{-}$ in Fig. 1 are connected through the factor node
$f_{\bm{\mathsf{y}}|\bm{\mathsf{x}}^{-},\bm{\mathsf{w}}^{-}}=\delta(\bm{y}-\bm{A}\,\bm{x}^{-}-\bm{w}^{-})$
stemming from the measurement model in (1). This means that the LMMSE estimate
of $\bm{\mathsf{x}}^{-}$ must account for the estimate of
$\bm{\mathsf{x}}^{-}$ and vice versa. Note that the denoising steps
corresponding to the red and green modules in Fig. 2 are similar to the
denoising procedure in VAMP. Thus, to customize the original VAMP framework so
as to recover the signal $\bm{x}$ under an arbitrary noise prior
$p_{\bm{\mathsf{w}}}(\bm{w})$, one should derive from scratch the LMMSE step
for both $\bm{\mathsf{x}}^{-}$ and $\bm{\mathsf{w}}^{-}$.

Algorithm 1 VAMP with arbitrary i.i.d. noise priors

1:$\mathbf{Require:}$ The channel matrix $\bm{A}\in\mathbb{R}^{M\times N}$;
the received vector $\bm{y}\in\mathbb{R}^{M}$; the maximum number of
iterations $T_{\text{max}}$.

2:$\mathbf{Initialize\,:\,}~{}\widehat{\bm{x}}^{-}_{\mathsf{e},0},\gamma_{\bm{x}^{-}_{\mathsf{e}},0}$,
$\widehat{\bm{w}}^{-}_{\mathsf{e},0}$ and $\gamma_{\bm{w}^{-}_{\mathsf{e}},0}$

3:for $t=0,\dots,T_{\text{max}}-1$ do

4: $\triangleright$ COMMENTDenoising $\bm{x}$

5:
$\bm{\widehat{x}}^{+}_{\mathsf{p},t}\;=\;g_{\bm{\mathsf{x}}}\left(\widehat{\bm{x}}^{-}_{\mathsf{e},t},\gamma_{\bm{x}^{-}_{\mathsf{e}},t}\right)$

6:
$\alpha_{\bm{x}_{\mathsf{p}}^{+},t}=\frac{1}{L}\,\sum\limits_{j=1}^{L}\\!g_{\bm{\mathsf{x}}}^{\prime}\left(\widehat{x}^{-}_{j,\mathsf{e},t},\gamma_{\bm{x}^{-}_{\mathsf{e}},t}\right)$

7:
$\gamma_{\bm{x}_{\mathsf{p}}^{+},t}=\gamma_{\bm{x}_{\mathsf{e}}^{-},t}/\alpha_{\bm{x}_{\mathsf{p}}^{+},t}$

8:
$\gamma_{\bm{x}_{\mathsf{e}}^{+},t}\;=\;\gamma_{\bm{x}_{\mathsf{p}}^{+},t}-\gamma_{\bm{x}_{\mathsf{e}}^{-},t}$

9:
$\bm{\widehat{x}}^{+}_{\mathsf{e},t}\;=\;(\gamma_{\bm{x}_{\mathsf{p}}^{+},t}\bm{\widehat{x}}^{+}_{\mathsf{p},t}-\gamma_{\bm{x}^{-}_{\mathsf{e}},t}\widehat{\bm{x}}^{-}_{\mathsf{e},t})/\gamma_{\bm{x}_{\mathsf{e}}^{+},t}$

10: $\triangleright$ COMMENTDenoising $\bm{w}$

11:
$\bm{\widehat{w}}^{+}_{\mathsf{p},t}\;=\;g_{\bm{\mathsf{w}}}\left(\widehat{\bm{w}}^{-}_{\mathsf{e},t},\gamma_{\bm{w}^{-}_{\mathsf{e}},t}\right)$

12:
$\alpha_{\bm{w}_{\mathsf{p}}^{+},t}=\frac{1}{L}\,\sum\limits_{j=1}^{L}g_{\bm{\mathsf{w}}}^{\prime}\left(\widehat{w}^{-}_{j,\mathsf{e},t},\gamma_{\bm{w}^{-}_{\mathsf{e}},t}\right)$

13:
$\gamma_{\bm{w}_{\mathsf{p}}^{+},t}=\gamma_{\bm{w}_{\mathsf{e}}^{-},t}/\alpha_{\bm{w}_{\mathsf{p}}^{+},t}$

14:
$\gamma_{\bm{w}_{\mathsf{e}}^{+},t}\;=\;\gamma_{\bm{w}_{\mathsf{p}}^{+},t}-\gamma_{\bm{w}_{\mathsf{e}}^{-},t}$

15:
$\bm{\widehat{w}}^{+}_{\mathsf{e},t}\;=\;(\gamma_{\bm{w}_{\mathsf{p}}^{+},t}\,\bm{\widehat{w}}^{+}_{\mathsf{p},t}-\gamma_{\bm{w}^{-}_{\mathsf{e}},t}\,\widehat{\bm{w}}^{-}_{\mathsf{e},t})/\gamma_{\bm{w}_{\mathsf{e}}^{+},t}$

16: $\triangleright$ COMMENTLMMSE estimation of $\bm{x}$ and $\bm{w}$

17:
$\bm{\widehat{x}}^{-}_{\mathsf{p},t}=\textswab{f}\big{(}\bm{x}_{\mathsf{e},t}^{+},\gamma_{\bm{x}_{\mathsf{e},t}^{+}},\bm{w}_{\mathsf{e},t}^{+},\gamma_{\bm{w}_{\mathsf{e},t}^{+}}\big{)}$

18:
$\bm{\widehat{w}}^{-}_{\mathsf{p},t}=\textswab{g}\big{(}\bm{x}_{\mathsf{e},t}^{+},\gamma_{\bm{x}_{\mathsf{e},t}^{+}},\bm{w}_{\mathsf{e},t}^{+},\gamma_{\bm{w}_{\mathsf{e},t}^{+}}\big{)}$

19:
$\gamma_{\bm{x}_{\mathsf{p}}^{-},t}\;=\;\gamma_{\bm{x}_{\mathsf{e}}^{+},t}/\alpha_{\bm{x}_{\mathsf{p}}^{-},t}$
with
$\alpha_{\bm{x}_{\mathsf{p}}^{-},t}\;=\;\left\langle\textswab{f}^{\prime}\right\rangle$

20:
$\gamma_{\bm{w}_{\mathsf{p}}^{-},t}\;=\;\gamma_{\bm{w}_{\mathsf{e}}^{+},t}/\alpha_{\bm{w}_{\mathsf{p}}^{-},t}$
with
$\alpha_{\bm{w}_{\mathsf{p}}^{-},t}\;=\;\left\langle\textswab{g}^{\prime}\right\rangle$

21:
$\gamma_{\bm{x}_{\mathsf{e}}^{-},t+1}\;=\;\gamma_{\bm{x}_{\mathsf{p}}^{-},t}-\gamma_{\bm{x}_{\mathsf{e}}^{+},t}$

22:
$\gamma_{\bm{w}_{\mathsf{e}}^{-},t+1}\;=\gamma_{\bm{w}_{\mathsf{p}}^{-},t}-\gamma_{\bm{w}_{\mathsf{e}}^{+},t}$

23:
$\widehat{\bm{x}}^{-}_{\mathsf{e},t+1}\;=\;\left(\gamma_{\bm{x}_{\mathsf{p}}^{-},t}\,\bm{\widehat{x}}^{-}_{\mathsf{p},t}-\gamma_{\bm{x}^{+}_{\mathsf{e}},t}\,\widehat{\bm{x}}^{+}_{\mathsf{e},t}\right)\Big{/}\gamma_{\bm{x}_{\mathsf{e}}^{+},t+1}$

24:
$\widehat{\bm{w}}^{-}_{\mathsf{e},t+1}\;=\;\left(\gamma_{\bm{w}_{\mathsf{p}}^{-},t}\bm{\widehat{w}}^{-}_{\mathsf{p},t}-\gamma_{\bm{w}^{+}_{\mathsf{e}},t}\widehat{\bm{w}}^{+}_{\mathsf{e},t}\right)\Big{/}\gamma_{\bm{w}_{\mathsf{e}}^{+},t+1}$

25:end for

26:Return $\widehat{\bm{x}}^{+}_{\mathsf{p},t}$

To derive the joint LMMSE estimates $\bm{\widehat{x}}^{-}_{\mathsf{p}}$ and
$\bm{\widehat{w}}^{-}_{\mathsf{p}}$ of $\bm{\mathsf{x}}^{-}$ and
$\bm{\mathsf{w}}^{-}$, we write the joint belief of $\bm{x}^{-}$ and
$\bm{w}^{-}$ pertaining to their joint pdf in (2) as the product of the three
beliefs depicted in Fig. 1:

$\displaystyle b(\bm{x}^{-},\bm{w}^{-})$
$\displaystyle\propto\mu_{\delta_{\text{x}}\rightarrow\bm{x}^{-}}(\bm{x}^{-})\,\cdot\,f_{\bm{\mathsf{y}}|\bm{\mathsf{x}}^{-},\bm{\mathsf{w}}^{-}}\,\cdot\,\mu_{\delta_{\text{w}}\rightarrow\bm{w}^{-}}(\bm{w}^{-})$
(5)
$\displaystyle=\mathcal{N}(\bm{x}^{-};\widehat{\bm{x}}_{\mathsf{e}}^{+},\gamma_{\bm{x}_{\mathsf{e}}^{+}}^{-1}\,\mathbf{I}_{N})~{}\delta(\bm{y}-\bm{A}\,\bm{x}^{-}-\bm{w}^{-})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\hskip
28.45274pt\times\mathcal{N}(\bm{w}^{-};\widehat{\bm{w}}_{\mathsf{e}}^{+},\gamma_{\bm{w}_{\mathsf{e}}^{+}}^{-1}\mathbf{I}_{M}).$

Finding the LMMSE estimates $\bm{\widehat{x}}^{-}_{\mathsf{p}}$ and
$\bm{\widehat{w}}^{-}_{\mathsf{p}}$ boils down to evaluating the following
integrals:

$\displaystyle\bm{\widehat{x}}^{-}_{\mathsf{p}}$
$\displaystyle=\frac{\int\int\bm{x}^{-}\,b(\bm{x}^{-},\bm{w}^{-})\,\text{d}\bm{x}^{-}\,\text{d}\bm{w}^{-}}{\int\int\,b(\bm{x}^{-},\bm{w}^{-})\,\text{d}\bm{x}^{-}\,\text{d}\bm{w}^{-}}$
$\displaystyle=\left(\gamma_{\bm{x}_{\mathsf{e}}^{+}}\,\mathbf{I}_{N}+\gamma_{\bm{w}_{\mathsf{e}}^{+}}\bm{A}^{\top}\bm{A}\right)^{-1}$
$\displaystyle\hskip
22.76228pt\times\left(\gamma_{\bm{x}_{\mathsf{e}}^{+}}\,\bm{\widehat{x}}_{\mathsf{e}}^{+}+\gamma_{\bm{w}_{\mathsf{e}}^{+}}\,\bm{A}^{\top}\left(\bm{y}-\bm{\widehat{w}}_{\mathsf{e}}^{+}\right)\right)$
$\displaystyle\triangleq\textswab{f}\left(\bm{\widehat{x}}_{\mathsf{e}}^{+},\gamma_{\bm{x}_{\mathsf{e}}^{+}},\bm{\widehat{w}}_{\mathsf{e}}^{+},\gamma_{\bm{w}_{\mathsf{e}}^{+}}\right),$
(6a) $\displaystyle\bm{\widehat{w}}^{-}_{\mathsf{p}}$
$\displaystyle=\frac{\int\int\bm{w}^{-}\,b(\bm{x}^{-},\bm{w}^{-})\,\text{d}\bm{w}^{-}\,\text{d}\bm{x}^{-}}{\int\int\,b(\bm{x}^{-},\bm{w}^{-})\,\text{d}\bm{w}^{-}\text{d}\bm{x}^{-}}$
$\displaystyle=\left(\gamma_{\bm{w}_{\mathsf{e}}^{+}}\,\mathbf{I}_{M}+\gamma_{\bm{x}_{\mathsf{e}}^{+}}\bm{Q}\right)^{-1}$
$\displaystyle\hskip
22.76228pt\times\left(\gamma_{\bm{w}_{\mathsf{e}}^{+}}\,\bm{\widehat{w}}_{\mathsf{e}}^{+}+\gamma_{\bm{x}_{\mathsf{e}}^{+}}\,\bm{Q}\left(\bm{y}-\bm{A}\,\bm{\widehat{x}}_{\mathsf{e}}^{+}\right)\right)$
$\displaystyle\triangleq\textswab{g}\left(\bm{\widehat{x}}_{\mathsf{e}}^{+},\gamma_{\bm{x}_{\mathsf{e}}^{+}},\bm{\widehat{w}}_{\mathsf{e}}^{+},\gamma_{\bm{w}_{\mathsf{e}}^{+}}\right).$
(6b)

where $\bm{Q}\triangleq\left(\bm{A}\,\bm{A}^{\top}\right)^{-1}$. Finally, to
find the posterior precisions $\gamma_{\bm{x}_{\mathsf{p}}^{-}}$ and
$\gamma_{\bm{w}_{\mathsf{p}}^{-}}$ of $\bm{x}^{-}$ and $\bm{w}^{-}$, we define
the divergences of $\textswab{f}(\cdot)$ and $\textswab{g}(\cdot)$ at
$\bm{x}_{\mathsf{e}}^{+}$ and $\bm{w}_{\mathsf{e}}^{+}$ as follows [3]:

$\displaystyle\textswab{f}^{\prime}$
$\displaystyle=\text{diag}\left(\frac{\partial\,\textswab{f}}{\partial\,\widehat{\bm{x}}_{\mathsf{e}}^{+}}\right)=\gamma_{\bm{x}_{\mathsf{e}}^{+}}\,\text{diag}\small{\left(\left(\gamma_{\bm{x}_{\mathsf{e}}^{+}}\,\mathbf{I}_{N}+\gamma_{\bm{w}_{\mathsf{e}}^{+}}\bm{A}^{\top}\bm{A}\right)^{-1}\right)},$
(7a) $\displaystyle\textswab{g}^{\prime}$
$\displaystyle=\text{diag}\left(\frac{\partial\,\textswab{g}}{\partial\,\widehat{\bm{w}}_{\mathsf{e}}^{+}}\right)=\gamma_{\bm{w}_{\mathsf{e}}^{+}}\,\text{diag}\left(\left(\gamma_{\bm{w}_{\mathsf{e}}^{+}}\,\mathbf{I}_{M}+\gamma_{\bm{x}_{\mathsf{e}}^{+}}\bm{Q}\right)^{-1}\right).$
(7b)

Now, after defining
$\alpha_{\bm{x}^{-}}\triangleq\left\langle\textswab{f}^{\prime}\right\rangle$
and
$\alpha_{\bm{w}^{-}}\triangleq\left\langle\textswab{g}^{\prime}\right\rangle$,
we use the fact that [3]

$\displaystyle\left\langle\textswab{f}^{\prime}\right\rangle=\frac{\gamma_{\bm{x}_{\mathsf{e}}^{+}}}{\gamma_{\bm{x}_{\mathsf{p}}^{-}}},~{}\text{and}~{}~{}\left\langle\textswab{g}^{\prime}\right\rangle=\frac{\gamma_{\bm{z}_{\mathsf{e}}^{+}}}{\gamma_{\bm{z}_{\mathsf{p}}^{-}}},$
(8)

to deduce the posterior precisions $\gamma_{\bm{x}_{\mathsf{p}}^{-}}$ and
$\gamma_{\bm{w}_{\mathsf{p}}^{-}}$ as shown in lines 19–20 of Algorithm 1.

##  3 Simulation Results

In this section, we assess the estimation performance of the proposed VAMP
algorithm with i.i.d. priors and benchmark it against the standard VAMP
algorithm. In all simulations, we set the number of time steps to
$T_{\textrm{max}}=100$, and perform $N_{\textrm{MC}}=100$ Monte-Carlo trials
for different values of the SNR defined in dB as

$\textrm{SNR}\,=\,\textrm{10
log}_{\textrm{10}}\Bigg{(}\frac{\left\lVert\bm{A}\,\bm{x}\right\rVert_{2}^{2}}{\left\lVert\bm{w}\right\rVert_{2}^{2}}\Bigg{)}.$
(9)

Here, each element $a_{ij}$ of $\bm{A}$ is drawn from the standard Gaussian
distribution (i.e., $a_{ij}\sim\mathcal{N}(0,1)$), and the vector
$\bm{x}_{\ell}$ is drawn from a Bernoulli-Gaussian density:

$p_{\bm{\mathsf{x}}}(\bm{x})=\rho\,\delta(\bm{x})+(1-\rho)\,\mathcal{N}\left(\bm{x};\bm{0},\mathbf{I}_{N}\right),$
(10)

where $\rho$ is the sparsity rate (i.e., the percentage of the non-zero
components). We use the normalized root MSE (NRMSE) as a performance measure
which is defined as:

$\textrm{NRMSE}~{}=~{}\frac{1}{N_{\textrm{MC}}}\sum_{\ell=1}^{N_{\textrm{MC}}}\frac{\|\bm{A}\,\bm{x}_{\ell}-\bm{A}\,\widehat{\bm{x}}_{\ell}\|_{2}}{\|\bm{A}\,\bm{x}_{\ell}\|_{2}}$

where $\bm{x}_{\ell}$ is $\ell$th realization of $\bm{\mathsf{x}}$ and
$\widehat{\bm{x}}_{\ell}$ is its reconstruction during the $\ell$th Monte-
Carlo trial.

We examine the NRMSE for two i.i.d. non-Gaussian priors modelling each
component $w_{i}$ of the noise vector $\bm{w}$ as a realization draw from

  * •

the Laplace prior
$\,p_{{\mathsf{w}}}(w_{i};\mu,b)=\frac{1}{2b}\,\textrm{exp}\left(-\frac{|w_{i}-\mu|}{b}\right)$
which better models heavy-tailed non-Gaussian noise sources (e.g., in
electronic devices).

  * •

the binary prior
$\,p_{{\mathsf{w}}}(w_{i};s)=\frac{1}{2}\,\delta(w_{i}-s)+\frac{1}{2}\,\delta(w_{i}+s)$
which can represent erasing noise sources or bit-flip distortions and
corruptions (e.g., in printing process).

We start by comparing standard VAMP to our algorithm when the noise vector $w$
is drawn from the Laplace distribution $\mu=0$ while varying $b$ to meet the
SNR level defined in (9). We also set $\rho=95\%$. Fig. 3 depicts the NRMSE of
both algorithms and shows how incorporating the noise model improves the
reconstruction accuracy, especially in the low SNR region. The fact that
Laplace and Gaussian distributions are not very far (in the KL-divergence
sense) diminishes the improvement of our algorithm in the high SNR regime
since both distributions are centered around the same mean $\mu=0$.

Fig. 3: NRMSE of VAMP with arbitrary i.i.d. priors vs. standard VAMP as a
function of the SNR with the noise vector drawn from the Laplace distribution
with $\mu=0$.

Significant improvements can be observed as soon as the noise considerably
deviates from the Gaussian model. When the noise vector $\bm{w}$ is drawn from
the binary prior with $a=1$, our algorithm reveals more robust to discrete
noise sources and outperforms the standard VAMP over the entire SNR range [0
dB, 20 dB] as depicted in Fig. 4. There, unlike our observation in Fig. 3, it
is seen how the improvement of our algorithm increases as the SNR increases
because of the substantial mismatch between the binary and Gaussian
distributions.

Fig. 4: NRMSE of VAMP with arbitrary i.i.d. priors vs. standard VAMP as a
function of the SNR with the noise vector drawn from the binary distribution.

##  4 Conclusion

In this paper, we build upon the vector approximate message passing algorithm
to handle non-Gaussian measurement noise models. We did so by incorporating
the arbitrary noise prior into standard VAMP by solving a joint LMMSE
estimation problem in addition to the MMSE noise denoiser developed in the
paper. Computer simulations with Laplace and binary noise models confirmed
that our proposed algorithm exhibits significant reconstruction improvements
as a function of the mismatch between the considered noise distribution and
the Gaussian one.

## References

  * [1] David L Donoho, Arian Maleki, and Andrea Montanari,  “Message-passing algorithms for compressed sensing,”  Proceedings of the National Academy of Sciences, vol. 106, no. 45, pp. 18914–18919, 2009.
  * [2] Sundeep Rangan,  “Generalized approximate message passing for estimation with random linear mixing,”  in 2011 IEEE International Symposium on Information Theory Proceedings. IEEE, 2011, pp. 2168–2172.
  * [3] Sundeep Rangan, Philip Schniter, and Alyson K. Fletcher,  “Vector approximate message passing,”  IEEE Trans. Inf. Theory, vol. 65, no. 10, pp. 6664–6684, 2019\.
  * [4] Junjie Ma and Li Ping,  “Orthogonal amp,”  IEEE Access, vol. 5, pp. 2020–2033, 2017.
  * [5] YC Eldar and G Kutyniok,  “Compressed sensing: Theory and applications cambridge univ,” 2012.
  * [6] Saleem A Kassam,  Signal detection in non-Gaussian noise,  Springer Science & Business Media, 2012.
  * [7] Ioannis Pitas and Anastasios N Venetsanopoulos,  Nonlinear digital filters: principles and applications, vol. 84,  Springer Science & Business Media, 2013.


Further author information: (Send correspondence to B.R.B.S.)
E-mail<EMAIL_ADDRESS>
# Mechanical and Optical Design of the HIRAX Radio Telescope

Benjamin R. B. Saliwanchik  Department of Physics, Yale University, New Haven,
CT, USA  Department of Physics, Brookhaven National Laboratory, Upton, NY, USA
Aaron Ewall-Wice  NASA Jet Propulsion Laboratory, Pasadena, CA, USA
Department of Astronomy, University of California, Berkeley, CA, USA  Devin
Crichton  School of Mathematics, Statistics, & Computer Science, University of
KwaZulu-Natal, Durban, South Africa  Institute for Particle Physics and
Astrophysics, ETH Zürich, Zürich, Switzerland  Emily R. Kuhn  Department of
Physics, Yale University, New Haven, CT, USA  Deniz Ölçek  Department of
Physics, McGill University, Quebec, Canada  Kevin Bandura  Department of
Computer Science and Electrical Engineering, and Center for Gravitational
Waves and Cosmology, West Virginia University, Morgantown, WV, USA  Martin
Bucher  School of Mathematics, Statistics, & Computer Science, University of
KwaZulu-Natal, Durban, South Africa  Astroparticle and Cosmology Laboratory,
University of Paris, Paris, France  Tzu-Ching Chang  NASA Jet Propulsion
Laboratory, Pasadena, CA, USA  H. Cynthia Chiang  School of Mathematics,
Statistics, & Computer Science, University of KwaZulu-Natal, Durban, South
Africa  Department of Physics, McGill University, Quebec, Canada  Kit Gerodias
Department of Physics, McGill University, Quebec, Canada  Kabelo Kesebonye
School of Mathematics, Statistics, & Computer Science, University of KwaZulu-
Natal, Durban, South Africa  Vincent MacKay  Department of Physics, University
of Toronto, Toronto, Canada  Kavilan Moodley  School of Mathematics,
Statistics, & Computer Science, University of KwaZulu-Natal, Durban, South
Africa  Laura B. Newburgh  Department of Physics, Yale University, New Haven,
CT, USA  Viraj Nistane  Department of Theoretical Physics and Centre for
Astroparticle Physics, Université de Genève, Genève, Switzerland  Jeffrey B.
Peterson  Department of Physics, Carnegie Mellon University, Pittsburgh, PA,
USA  Elizabeth Pieters  Department of Physics, McGill University, Quebec,
Canada  Carla Pieterse  School of Mathematics, Statistics, & Computer Science,
University of KwaZulu-Natal, Durban, South Africa  Keith Vanderlinde
Department of Physics, University of Toronto, Toronto, Canada  Jonathan L.
Sievers  Department of Physics, McGill University, Quebec, Canada  School of
Chemistry and Physics, University of KwaZulu-Natal, Durban, South Africa
Amanda Weltman  Department of Mathematics, University of Cape Town, Cape Town,
South Africa  Dallas Wulf  Department of Physics, McGill University, Quebec,
Canada

###### Abstract

The Hydrogen Intensity and Real-time Analysis eXperiment (HIRAX) is a planned
interferometric radio telescope array that will ultimately consist of 1024
close packed 6 m dishes that will be deployed at the SKA South Africa site.
HIRAX will survey the majority of the southern sky to measure baryon acoustic
oscillations (BAO) using the 21 cm hyperfine transition of neutral hydrogen.
It will operate between 400-800 MHz with 391 kHz resolution, corresponding to
a redshift range of $0.8<z<2.5$ and a minimum $\Delta z/z$ of $\sim$0.003
(frequency resolution $500<R<1000$). One of the primary science goals of HIRAX
is to constrain the dark energy equation of state by measuring the BAO scale
as a function of redshift over a cosmologically significant range. Achieving
this goal places stringent requirements on the mechanical and optical design
of the HIRAX instrument which are described in this paper. This includes the
simulations used to optimize the mechanical and electromagnetic
characteristics of the instrument, including the dish focal ratio, receiver
support mechanism, and instrument cabling. As a result of these simulations,
the dish focal ratio has been reduced to 0.23 to reduce inter-dish crosstalk,
the feed support mechanism has been redesigned as a wide (35 cm diam.) central
column, and the feed design has been modified to allow the cabling for the
receiver to pass directly along the symmetry axis of the feed and dish in
order to eliminate beam asymmetries and reduce sidelobe amplitudes. The beams
from these full-instrument simulations are also used in an astrophysical
m-mode analysis pipeline which is used to evaluate cosmological constraints
and determine potential systematic contamination due to physical non-
redundancies of the array elements. This end-to-end simulation pipeline was
used to inform the dish manufacturing and assembly specifications which will
guide the production and construction of the first-stage HIRAX 256-element
array.

##  1 Introduction

The fact that the universe is in a state of accelerated expansion today is
supported by the observational evidence from various cosmological tools such
as the Cosmic Microwave Background (CMB), Type 1a Supernovae (SN1a), and
Baryon Acoustic Oscillations (BAO)[1]. The increasing rate of expansion is
attributed to an unknown cosmological component called Dark Energy, which
constitutes around 70% of the total energy density of the current universe.
Current data shows Dark Energy began to affect the dynamics of the universe
around a redshift of z $\sim$ 2 (10.5 Gyr ago), and became the dominant
component of the energy density at around a redshift of z $\sim$ 0.5 (5.2 Gyr
ago). Constraining the properties of Dark Energy, and the evolution of the
universe over cosmic timescales, is an important goal of modern cosmology.

Baryon acoustic oscillations provide a standard ruler, a structure of a fixed
co-moving scale, which can be used to measure the expansion history of the
universe over a wide range of redshifts[2]. Baryon acoustic oscillations in
the observed matter density power spectrum arise because all modes are
initially excited with the same phase on superhorizon scales, and only the
growing mode is excited. Before recombination the photons and baryons are
tightly coupled to each other owing to Thomson scattering and oscillate much
like ordinary sound, following horizon crossing. However, after recombination
the photons free stream, causing the phase oscillations of the plasma at the
moment of recombination to become frozen into the power spectrum. The
characteristic scale of the first peak of this power spectrum is the sound
horizon (comoving distance that density waves in the photon-baryon fluid could
have travelled until decoupling), which is strongly constrained by the CMB to
be 146.8 $\pm$ 1.8 Mpc [3]. Due to the variations in matter density produced
by the BAO, and structure formation favoring areas of over-density, the
distribution of galaxies subsequently displays a spatial correlation function
which corresponds to the BAO power spectrum. This large scale structure
produced by the BAO has been detected at high significance in optical galaxy
surveys of the low redshift universe, and has provided an important constraint
on cosmological parameters, including Dark Energy[4, 5, 6].

One method to measure the BAO spectrum at higher redshift is to directly
measure the neutral hydrogen density associated with galaxies by using the 21
cm emission line. This probe is particularly attractive because of the
omnipresence of neutral Hydrogen at all redshifts. Since the characteristic
scale of the BAO is very large, there is no need to detect individual galaxies
at high resolution, and instead a measurement of the power spectrum of neutral
hydrogen density on large scales can completely capture the BAO structure.
This technique is called intensity mapping [7]. Galaxies are observed
collectively via low spatial resolution measurements of redshifted 21 cm lines
of neutral Hydrogen, which also has the potential to make 3D intensity mapping
more efficient as compared to optical galaxy surveys.

The Hydrogen Intensity and Real-time Analysis eXperiment (HIRAX) will measure
the HI density field over the redshift range $0.8<z<2.5$, bracketing the
redshift at which dark energy begins to affect the dynamics of the universe.
The large daily survey area (approximately 1,000 deg2), and real-time
processing capabilities of HIRAX will also make it an effective observatory
for detecting and monitoring radio transients, such as fast radio bursts
(FRBs) and pulsars[8, 9, 10, 11, 12]. The Southern Hemisphere survey location
also provides overlap with the survey fields of cosmology surveys undertaken
by ACT[13], SPT[14], DESI[15], and the Vera Rubin Observatory[16], allowing
for significant cross correlation studies. HIRAX will also deliver a blind HI
21 cm line survey at $0.8<z<2.5$. This will complement the MeerKAT Absorption
Line Survey (MALS; $0<z<1.4$) by extending the exploration of cold atomic gas
in the circum-galactic and inter-galactic medium to higher redshifts[17].

Detecting the BAO in the presence of significant astronomical foregrounds,
most notably synchrotron radiation from the Milky Way galaxy, is a significant
challenge. The high level requirements for the array to achieve its science
goals are sensitivity, redundancy, and control of systematics. The array is
designed to meet the sensitivity requirements by a combination of low system
temperature and large collecting area. The later aspect is what drives us to a
large array, consisting of approximately one thousand 6 m dishes. The former
is achieved by developing low loss antennas and amplifiers, which are
addressed in Kuhn et al.[18]. The simulations presented in this work are part
of the effort to design and produce an array with very high levels of
redundancy between elements, and with very low systematic levels.

Tight control over redundancy is a new design driver that has emerged for the
next generation of 21 cm observatories such as HIRAX[19]. We are aiming to
solve the problem of redundancy in hardware, which requires a level of
precision that is significantly higher than existing experiments such as
CHIME[20] and HERA[21]. That is, rather than correcting for the differences
between array element structures and positions completely in software, we aim
to reduce the element differences as much as possible in hardware, before
applying fine corrections in software. This requires significantly more
stringent specifications than normal for a radio telescope, on the order of 1
part in 1000 relative to the wavelength, driven largely by the high foreground
to signal ratio for the 21 cm signal. Edge effects due to the finite size of
the array are also a source of non-redundancy, and must be taken into
consideration. Additionally, we are also concerned with minimizing crosstalk
between the elements in the array, which can result from either sky signals
bouncing between elements, or from amplifier noise being transmitted from one
element and picked up in another. The issues we want to address in this work
are, broadly, how to optimize the design of the instrument to minimize
systematics, and to determine what level of element redundancy is necessary,
and whether this level is feasible to achieve in hardware, within cost.

##  2 Telescope Design

The full HIRAX array will consist of 1024 parabolic dishes 6 m in diameter,
with a focal ratio of f/D = 0.23. The optimization of the focal ratio is
discussed in Section 3.1. A first-stage array consisting of 256 elements is
fully funded, and bidding for construction of the array is underway. The full
array will be arranged in a close-packed $32\times 32$ element configuration
to enhance sensitivity on the BAO length scales, and to increase the
redundancy of the array. The increased array redundancy in turn simplifies
calibration, reduces the number of correlations to be calculated, and reduces
the data volume to be stored. The telescope design and array layout can be
seen in Figure 1.

Figure 1: Top: Conceptual illustration of the HIRAX telescope design. The
dishes are 6 m in diameter, with a focal ratio of 0.23. The receiver is
supported by a fiberglass column which allows for axisymmetric cabling, which
is important for beam symmetry, polarization, and low sidelobe amplitudes, and
which is rigid under high winds. The receiver is a dual-polarization
cloverleaf antenna, with a metal “can” structure to provide a backplane, and
reduce spillover and crosstalk. Analog signals from the dishes are transmitted
via RF over fiber to the back end correlator. Bottom: A rendering of the full
1024 element array. The array is close-packed to provide improved sensitivity
on the BAO angular scales, and to provide highly redundant baselines, which
facilitate calibration and correlation.

The dish design includes a low mount that allows for easier access to the
feed, reduces wind force on the base, and reduces cost. The HIRAX array
operates in drift scan, so the only degree of freedom necessary in the dish
pointing is to vary the elevation between $\pm 30^{\circ}$ degrees from zenith
to encompass the intended survey field.

The receiver is supported by a central fiberglass column. The architecture of
this support structure was motivated by the optimal path of the cables from
the feed to the surface of the dish. Because the cables are metal and lie in
the optical path, the sidelobes of the primary beam depend sensitively upon
the cable placement and angle. Simulations show that sidelobe levels and
asymmetry are minimized when the cables run straight down the boresight axis
of the dish, rather than at an inclined angle, as discussed in Section 3.1. A
column was therefore chosen to provide a natural environmental enclosure
following this cable path. The support structure includes additional provision
for fully enclosing the feed and the radio-frequency over fiber (RFoF) modules
that are co-located with the feed, for weatherproofing and protection of the
full receiver system.

The HIRAX feed is a dual-polarization cloverleaf antenna consisting of copper
layers in an FR-4 (printed circuit board) composite (shown in Figure 3). This
method of manufacturing is easy to produce at scale, easy to assemble with a
minimal jig for alignment, and the FR-4 protects the metallic elements from
corrosion. The cloverleaf antenna is a proven design, having been deployed on
the Canadian Hydrogen Intensity Mapping Experiment (CHIME) [22]. While the
original CHIME feeds were passive, the versions developed for HIRAX integrate
the first stage low noise amplifier (LNA) with the antenna balun to reduce the
system noise. As a result of the lower system noise, the feed material was
changed to FR-4 (compared to teflon in CHIME), which has higher loss, but is
substantially cheaper. In an additional change from the CHIME architechture,
the signal from the individual array elements is transmitted to the correlator
by RFoF, instead of via coaxial cable. This reduces cost relative to coax for
a large array, with no loss to performance. Noise temperature measurements of
the HIRAX feeds are being conducted using a cryogenic RF chamber to perform a
differential y-factor measurement between hot (295 K) and cold (77 K) loads.
Details of those measurements are presented in Kuhn et al. [18].

Since we intend to use our EM simulations to optimize elements of the feed and
dish design, we wanted to verify that those simulations were accurate. We
conducted the first range measurements of the HIRAX feed beams at the MESA
antenna test facility at the NASA Jet Propulsion Laboratory, and further
measurements at North Carolina State University. Figure 2 shows the measured
beams of the HIRAX cloverleaf antenna and can, and compares this to the CST
simulations. These measurements confirm the accuracy of the simulations, and
that the feeds were produced according to design.

[capbesideposition=right,center,capbesidewidth=0.3 ]figure[]

Figure 2: E-plane co-polarization beam measurements of the HIRAX cloverleaf
antenna and can (orange), compared to the simulated beams (blue) in 50 MHz
increments across the HIRAX observing band. These measurements confirm the
beam width and gain are as designed and simulated. Differences in the depth of
nulls are due to low amplitude reflections in the range, which are not present
in free-space simulations. Measurements were taken at the North Carolina State
University test range.

Work is currently underway to develop a drone calibration platform, which will
allow for beam measurements of the full instrument, including the 6m dish and
amplification chains, in order to verify the simulated performance of the
instrument, fully characterize the instrument, and aid in calibration. This
platform has already been used to perform beam measurements of the Baryon
Mapping Experiment (BMX) array at Brookhaven National Laboratory[23]. Beam
measurements of the HIRAX dishes is currently being delayed due to the ongoing
Covid-19 pandemic, but will be performed when international research travel is
permitted by institutions again.

##  3 Electromagnetic Simulations

###  3.1 Design Simulations

Electromagnetic simulations of the instrument were performed to select between
design options, and to optimize the design elements. Instrument elements
explored include:

  * •

The diameter of the feed “can” structure

  * •

The nature and positioning of the power cabling for the feed

  * •

The feed mechanical support mechanism, in particular deciding between a system
with several “feed legs” and a monolithic “feed column”

  * •

Optimizing the dish focal ratio to reduce crosstalk and other array effects

Figure 3: The HIRAX feed and CST models. Top left: The HIRAX dual-polarization
cloverleaf antenna (green PCB) and can (white). Top right: CST model of the
HIRAX feed and can. Only one polarization is fed in the simulation, as
indicated by the blue bar and red arrow. Bottom left: The 6 m dish and 10 cm
diameter fiberglass feed column model. Bottom right: The dish and 35 cm
diameter feed column model.

####  3.1.1 Feed Can

The HIRAX cloverleaf antenna is backed with a metal structure referred to as
the “can” (Figure 3), a combination of a ground plane and a cylindrical
surface which helps to circularize the beam and reduce the beam FWHM, to
prevent over illuminating the dish, and reduce spillover. Spillover is
especially a concern in a close-packed array, as proximity of the neighboring
dishes increased pickup of spilled power. Without the can, the FWHM of the
cloverleaf alone is $90^{\circ}$ ($60^{\circ}$) at 400 MHz (800 MHz), while
with a nominal 330 mm diameter can the FWHM is reduced to $70^{\circ}$
($67^{\circ}$) at 400 MHz (800 MHz). Increasing the can diameter reduces
$S_{11}$ because it provides a better ground plane, but it also decreases
aperture efficiency, because it blocks more of the beam reflected from the
dish. The effect on $S_{11}$ was found to be a slow function of radius, while
the loss of aperture efficiency scales roughly with $r^{2}$. Therefore a small
diameter can was used; $330$ mm in diameter, slightly larger than the
cloverleaf itself.

####  3.1.2 Feed Support

In initial prototype dishes the feed was supported above the dish by four
cylindrical aluminum conduit “legs” 25 mm in diameter, forming the edges of a
square right pyramid, with the legs meeting at an apex behind the feed, and an
apex angle of approximately $45^{\circ}$. However, simulations showed that
scattering off the metal feed legs introduced strongly polarized components
into the beam, and increased sidelobe levels, motivating a design change to
non-conductive support structure.

Similar feed legs could have been constructed out of a non-conductive
composite material such as fiberglass, but there were additional mechanical
problems with the feed leg model. Positioning the feed could be accomplished
by a mechanism that slid over the conduit legs, and could adjust the position
of the receiver independently on each leg. However, this mechanism did not
translate easily into adjustments in a simple orthogonal coordinate system.
Furthermore, it was prone to motion in several axes, including compression
towards the dish (z-axis), and rotation around the z-axis. Weatherproofing and
feed access were also issues with the feed leg support mechanism.

A proposed solution was to instead use a fiberglass column, extending from the
vertex of the dish to the feed (Figure 3). Mechanical FEA simulations showed
this feed support mechanism simultaneously solved several problems. It was
significantly more stable in all axes, particularly the z-axis, which is
especially important for ensuring the instrument is in focus. Large diameter
cylinders are especially strong, and a 35 cm diameter cylinder, just larger
than the feed can, was found to be able to withstand the peak survival wind
speed of 44.4 m/s at the SKA South Africa site, while deflecting from the
focal point by less than 1 mm. A feed column also provided a weather proof
enclosure, easy access to the feed and RFoF module by a detachable section at
the top of the cylinder, and a stable point from which fine adjustments in
feed position could be made in an orthogonal basis system. A similar system
was employed by the CBASS collaboration, using an RF transparent receiver
support column constructed of Plastazote LD45, a closed-cell polyethylene
foam[24].

A small diameter (10 cm) feed column was also proposed, and a prototype was
fielded in the dish shown in Figure 8. However, in addition to greater
mechanical stability for fixed wall thickness, a large diameter column also
has some RF advantages. EM simulations showed that a small diameter column
reduces the main beam amplitude by up to 0.5 dB, due to the long path length
of dielectric material along the main beam line of sight (see Figure 4). Both
models slightly alter the sidelobe amplitudes, but in different directions at
different frequencies. A larger diameter column does slightly reduce the
instrument radiative efficiency compared to a smaller diameter column, for
fixed wall thickness, because of the greater net amount of dielectric in the
beam path, but this effect is also not significant. There was a $<1\%$
difference in radiative efficiency between 10 cm and 35 cm diameter columns
with 5mm wall thickness, several times thicker than what is mechanically
required.

[capbesideposition=right,center,capbesidewidth=0.35 ]figure[]

Figure 4: Comparison of 10 cm diameter receiver support column (orange) and 35
cm column (blue), with 5 mm wall thickness, at 400 MHz. A smaller diameter
support column slightly decreases main beam amplitude (by 0.5 dB). Sidelobes
can be slightly decreased or enhanced, depending on frequency.

####  3.1.3 Power Cabling

In connection with the issue of feed support, it was suspected that the cables
for powering the LNA in the feed might affect the instrument beams.
Simulations showed that a single coax cable running from the feed to the dish
edge induced polarization structures in the beam and enhanced sidelobe
structure at a level comparable to the feed support legs themselves, because
an incident plane wave polarized in the same direction as the cable will
excite currents along it. Additionally, if there is only one cable, and not a
symmetric set of feed legs, the sidelobe effect is asymmetric, which is
undesirable.

This effect can be reduced, but not eliminated, by reducing the cable
diameter. Significant improvements were seen by reducing from coax
($\mbox{$\sim$}1$ cm diam.) to 15 AWG wire ($\mbox{$\sim$}1.5$ mm diam.).
Reductions below approximately 1 mm diameter produce diminishing returns. 15
AWG wire is sufficient for the current required to power the LNA (200mA), and
since in the full array the sky signal will be transmitted on RFoF, neither
coax, nor any other metal cable elements will be required to run from the dish
to the feed, further reducing scattering.

The asymmetric sidelobes from the cabling can be further reduced by orienting
the power cables parallel to the beam, that is, running them from the feed to
the dish vertex. Simulations were run exploring a range of cable angles,
between 90 degrees (perpendicular to the beam axis) and zero degrees (parallel
to the beam axis), and with the feed side of the cable positioned at a radius
from the feed center greater than the can radius (because it mechanically had
to wrap around the can to access the feed). The sidelobe behavior was not a
monotonic function of angle, but the asymmetry was generally reduced with
increasing angle. However, even at zero degrees, the effect was not completely
eliminated, as there is still a small component of the cable perpendicular to
the beam, due to the offset from the feed center. Running the cable down the
exact symmetry axis of the feed and telescope, however, does completely remove
the asymmetric sidelobe effect (See Figure 5). This configuration also reduces
the sidelobe amplitude by up to 7 dB relative to the initial (perpendicular)
wiring case.

Figure 5: HIRAX telescope beams at 400 MHz using 15 AWG wire
($\mbox{$\sim$}1.5$ mm diameter) to power the LNA embedded in the balun of the
feed. Top: wire routed from focal point to edge of dish. Middle: wire routed
from the edge of the feed can ($\mbox{$\sim$}18$ cm offset from the axis of
symmetry of the dish and antenna) to the vertex of dish. Bottom: wire routed
though cloverleaf balun, precisely along the boresight axis of the antenna, to
the dish vertex. The completely symmetric wiring case eliminates all
asymmetric beam features, and reduces the sidelobe amplitudes by up to 7 dB
relative to the initial configuration.

We determined that the feed could easily be modified to allow for this cable
routing since the balun is already hollow, and a hole in the base board could
be added to the PCB design to allow the cable to pass through. There was a
minor concern that this would require slightly offsetting the feed point of
the antenna, which has also been shown to produce beam asymmetry. However,
simulations showed that as long as both the cable and the feed point were both
within approximately the balun diameter ($\mbox{$\sim$}1$ cm), which is
mechanically feasible, then the effects of both offsets are negligible. The
feed design has therefore been modified to allow the power cable to be exactly
on axis, and the feed point to be $<1$ cm offset from center.

####  3.1.4 Dish Focal Ratio

Reflections within and between an interferometer’s antenna elements introduce
spectral structure in the antenna gain at frequency scales corresponding to
the time delay of the reflection. If these reflections reach to fine frequency
scales, they can cause the otherwise smooth continuum foregrounds to
contaminate the spectral scales important for cosmology. Lowering the focal
ratio reduces cross coupling between feeds, reducing reflections between
elements but also has the potential to increase reflections within a single
antenna element.

We evaluated the relative delay performance of dishes with different focal
ratios, by running CST simulations of a plane-wave excitation of the dish and
feed from zenith, and use the voltages in a 50 $\Omega$ termination of the
feed to estimate the time-domain response kernel that leaks these foregrounds
using the formalism from Ewall-Wice et al.[25]. In order to evaluate the
relative contributions from inter-dish versus intra-dish reflections, we
perform our simulations with two different sets of boundary conditions. First,
we set our boundary conditions to be open, simulating a dish in isolation,
which gives us the contributions from reflections within the dish. Second, we
perform simulations with periodic rectangular boundary conditions to simulate
one of the HIRAX antennas embedded within an infinite rectangular packed
array. This response function corresponds to the spectral structure that would
appear in an autocorrelation for a sky with only a source at zenith. While
auto-correlations are not very sensitive to 21 cm fluctuations, the
simulations are useful for evaluating the relative levels of spectral
structure. The impact of intra-dish reflections on cross-correlations between
dishes is the subject of ongoing work.

In Figure 6, we compare the time-domain responses of antennas with different
focal ratios with open and periodic boundary conditions. Comparing the open
boundary curves, we see that decreasing the focal ratio worsens intra-dish
reflections by several dB. However, comparing the curves from periodic
boundary conditions, we see that inter-dish reflections dominate the time-
response of a single element, especially at large delays, which are most
sensitive to cosmology. From the perspective of reducing spectral structure,
it makes sense for us to decrease the focal ratio for a fixed dish diameter as
much as is allowed by mechanical constraints and cost. Reducing the focal
ratio also reduces the potential for line-of-sight noise coupling by reducing
$S_{31}$. We demonstrate this reduction in the right hand panel of Figure 7
where we plot the $S_{31}$ coupling parameter between two parallel polarized
feeds of two adjacent antennas.

However, reducing the focal ratio can also decrease the illumination of the
dish and degrade sensitivity. In the bottom panel of Fig 7, we see that the
aperture efficiency of the dish is moderately affected by reducing the focal
ratio. A focal ratio of 0.23 reduces $S_{31}$ by $\mbox{$\sim$}10$ dB relative
to 0.25, while only reducing aperture efficiency by $\mbox{$\sim$}5\%$. The
selected focal ratio value of 0.23 was optimal for minimizing $S_{31}$ and the
delay kernel amplitude at high delay ($\tau>50$ ns), within the constraints of
mechanical support and per-element cost.

Figure 6: The delay kernel of the HIRAX dish towards a plane-wave at zenith
with open boundary conditions (dashed lines) and periodic boundary conditions
with two meters between dish edges (solid lines). This kernel determines the
extent to which foregrounds are leaked from small line-of-sight Fourier modes
where they exist intrinsically to large line-of-sight Fourier modes. In order
to minimize foreground leakage and maximize our ability to recover
cosmological fluctuations, this kernel should be kept as narrow as possible in
delay. Different colors represent different values for the dish focal ratio
for a fixed dish diameter of six meters. Reducing the focal ratio deepens the
dish and lowers the feed below the rim. Comparing the open boundary curves, we
see that reducing the focal ratio increases intra-dish reflections, raising
the time-domain response by several dB. The curves with periodic boundaries
include reflections of the sky signal off of nearby antennas. While a higher
focal ratio is preferred for a single dish in isolation, we see that the
inter-dish reflections dominate over intra-dish reflections, especially at
higher delays. These inter-dish reflections are mitigated by lowering the
focal ratio. Thus, when considering the HIRAX dish embedded in an array, we
find that lowering the focal ratio has a net beneficial impact on the time-
domain response.

Figure 7: Comparison of nearest-neighbor crosstalk ($S_{31}$, left), and
aperture efficiencies (bottom), for a wide range of focal ratios. Decreasing
focal ratio (deeper dishes) reduces crosstalk. Decreasing focal ratio also
reduces the aperture efficiency (right), reducing overall system sensitivity.
Geometries below f/D$\sim$0.225 are also increasingly difficult to
mechanically support, requiring significantly more backing structure due to
the deep shape. The selected value of 0.23 reduces $S_{31}$ by
$\mbox{$\sim$}10$ dB relative to 0.25, while only reducing aperture efficiency
by $\mbox{$\sim$}5\%$.

###  3.2 Design Tolerance Simulations

In addition to the design selection and optimization simulations, an extensive
suite of simulations was performed exploring the results of inaccuracies in
the manufacture or assembly of the instrument which degrade the redundancy of
the array, in order to place specifications on the components and assembly.
These simulations included, among others:

  * •

Varying feed position in six axes (translations and rotations)

  * •

Varying dish diameter, holding focal ratio constant

  * •

Varying dish focal ratio, holding diameter constant

  * •

Perforations in dish conductive surface

  * •

Reduced dish surface conductivity

The resulting beams from these simulations were in turn used to mock-observe
simulated 21 cm skies, as described in Section 5, in order to examine the
effects of these mechanical errors on observations, and the recovered data
products. In addition to non-uniform manufacturing, edge effects due to the
finite extent of the array are also a source of non-redundancy. This category
of effects is not explored here, but will be the subject of future work.

####  3.2.1 Feed Position

For all simulations, the coordinate system origin is at the nominal feed
position. The positive z-axis points from the feed to the dish vertex, the
x-axis is horizontal when the dish is pointed at the horizon, and the y-axis
is vertical. For the feed position simulations an array of simulations were
performed for each axis, shifting the feed along that axis. The positions
were: from 0 mm to 10 mm from the nominal focal point in 1 mm steps along the
x and y-axes, and from -5 mm to 5 mm in 1 mm steps along the z-axis. For the
three rotational axes, the feed was rotated from $0^{\circ}$ to $5^{\circ}$ in
$0.5^{\circ}$ increments, around each axis.

The primary effect of translation along the z-axis is to throw the instrument
out of focus, which reduces gain, and alters beam sidelobe structure.
Translations in the other two dimensions primarily result in beam pointing
errors, and asymmetric sidelobes. Rotation around the z-axis does not change
the shape of the beam, but due to the non-azimuthally symmetric sidelobe does
result in relative changes in sidelobe amplitudes between pairs of dishes with
different rotations. Rotations around the x and y-axis primarily change the
pointing of the beam, and creates asymmetric sidelobes. The resulting design
tolerance was $\pm 1$ mm in all linear translation directions, $\pm 1.5$
arcmin for the azimutal rotation, and $\pm 2.5$ armin for the other rotations.
The summary of design specifications can be seen in Table LABEL:tab:specs.

####  3.2.2 Dish Diameter

The dish diameter was varied from 599 cm to 601 cm in 1 mm steps, while
holding the focal ratio constant (so as not to conflate the results with those
of the focal ratio simulations). Changing the dish diameter throws the
instrument out of focus, and changes the beam width. We determined the
accuracy, averaging over they array, must be within $\pm 3$ mm, and the
individual dish precision within $\pm 1$ mm.

####  3.2.3 Focal Ratio

The focal ratio was varied from f/D = 0.245 to 0.255 in steps of 0.001, while
holding the dish diameter fixed at 6 m. Changing the focal ratio changes the
aperture efficiency, and can alter the sidelobe levels. These simulations were
distinct from those exploring the focal ratio in Section 3.1, as they were
intended to explore fine variations around a nominal value, to constrain
manufacturing tolerances, not to explore large changes in the nominal value.
These simulations contributed to the dish shape accuracy specification in
Table LABEL:tab:specs ($\pm 3$ mm maximum deviation from the ideal paraboloid
curve).

####  3.2.4 Dish Perforations

Radio dishes are often constructed of metal mesh instead of solid metal
panels, for purposes of lightweighting, lowering wind cross section, and cost
reduction. As long as the gaps in the mesh are substantially smaller than the
wavelength, there is no reduction in the performance, but as the gap size
approaches the wavelength, the dish becomes increasingly transparent. A common
criteria given for maximum gap sizes is $\lambda/10$. We explored the effects
of hole size on the instrument beams in order to set specifications on what
types of mesh were permissible and to regulate production errors resulting in
accidental perforations or tears in the dish surface. As a first-pass method
of simulating this effect, the dish was perforated with a regular array of
circular holes, arranged in 16 equally spaced azimuthal angles ($22.5^{\circ}$
separation), and 6 radial distances from the center of the dish with 25 cm
spacing, for a total of 96 perforations in the dish. The hole diameter was
varied from 0.5 cm to 5.0 cm, in 5 mm steps.

This method was used because a realistically large number of holes in the dish
would vastly increase the complexity of the meshing, and therefore the
simulation run times. The holes were restricted to a central 1.5 m radius of
the dish to again simplify the model: at this diameter holes that are circular
in cross section and projected from a plane normal to the dish vertex are not
significantly distorted by the curvature of the dish. If holes were projected
out to the full 3 m radius of the dish, they would be extremely elongated in
the azimuthal direction, and instead of exploring one hole radius, we would be
exploring a radially and azimuthally varying range of values. Altering the
projection angle of the holes with the radius would allow the full surface
area to be perforated, and will be explored in future analysis analysis.

Perforations create additional structure within the beam, which increases with
increasing hole size, and at higher frequencies. At larger hole sizes it
decreases gain, and increases power from the ground (though this effect was
not modeled in our simulations). The resulting specification on perforation
size from this analysis (5 mm maximum gap dimension, see Table
LABEL:tab:specs) is significantly stronger than the $\lambda/10$ criteria (as
large as 7.5 cm at 400 MHz). In future this set of simulations will be
expanded to a more realistic model using a pseudo-randomly placed array of
perforations over a wider area of the dish.

####  3.2.5 Conductivity

Lastly, we simulated different conductivity values of the dish’s reflective
surface to determine any impact on the beams (if for example, the metal mesh
used in the surface was excessively fine, and the resistance per strand
increased); or if a non-uniformity in the conductivity would have an effect on
the beams (as if the mesh density varied, or if different metals were used in
different parts of the dish). Two sets of simulations were run to explore
these possibilities, one in which the conductivity of the entire dish was
varied, and one in which the dish was modeled as two halves with different
conductivities. The latter represents a limiting case, and the beam effects
calculated from a dish surface with maximally inhomogeneous conductivity can
be interpreted as an upper limit. In the first set of simulations, the
conductivity was varied from $6\times 10^{7}$ S/m (approximately the
conductivity of copper) to $1\times 10^{6}$ S/m (roughly the conductivity of
stainless steel, and the lowest conductivity of materials one might reasonably
use) in steps of $5\times 10^{6}$ S/m. In the second set of simulations, one
half of the dish was held at $6\times 10^{7}$ S/m, while the other was stepped
down from $6\times 10^{7}$ S/m in increments of $5\times 10^{6}$ S/m to
$1\times 10^{6}$ S/m. This series of simulations showed essentially no changes
in the beam over the range of explored values. This is in agreement with
experimental measurements of the radiation efficiency of dipole and meander
antennas as a function of conductivity, which show no change in efficiency
until the conductivity is below $1\times 10^{6}$ S/m [26]. The model of
conductivity variations employed in these simulations is a simplified model of
the large possible parameter space of variations in the dish surface
conductivity, and will be expanded upon in future simulations.

Table 1: HIRAX Telescope Element Specifications

Element | Specification | Notes
---|---|---
Axial symmetry of |  $\pm$1 mm |
receiver support |  |
Receiver support |  $<0.5$ dB |
RF attenuation |  |
Deviation of power |  $\pm$2 mm |
cabling from boresight |  |
Rigidity of |  $\pm 0.5$ mm | In x,y, and z dimensions
receiver support |  |
Positioning of receiver |  $\pm 0.5$ mm | In x,y, and z dimensions
relative to focal point |  |
Orientation of receiver |  $\pm 2.5$ arcmin | polar angle
relative to boresight |  $\pm 1.5$ arcmin | azimuthal angle
Dish diameter |  $\pm 3$ mm | Accuracy
|  $\pm 1$ mm | Precision
Dish shape accuracy |  $\pm 3$ mm | Deviation from ideal paraboloid
Dish electrical connectivity |  $<5$ mm | Maximum dimension of gaps
Dish surface conductivity |  $>1\times 10^{6}$ S/m |

  * Instrument specifications determined by EM simulations. This is a subset of the total system specifications.

##  4 Photogrammetry Measurements

The specifications outlined above are significantly more stringent than those
typically required for radio telescopes at these wavelengths, due to the
increased accuracy necessary to control array redundancy. Achieving these
specifications within cost is a significant engineering challenge.
Photogrammetry measurements of our first prototype dishes were performed to
ascertain whether the specifications had been meet, and to inform the design
and production of future dishes. The results show that already in the first
prototype we are close to achieving the specified goals, though improvements
still need to be made.

The depth of the dishes is designed to reduce crosstalk between the elements
in the close-packed array, but this is an uncommon design for radio telescopes
(though similarly deep dishes have been used for other close-packed radio
arrays such as HERA[21], CHIME[20], and CHORD[27]), and requires greater
support than more typical shallow profile dishes. Figure 8 shows a prototype
dish located at the Hartebeesthoek Radio Astronomy Observatory (HartRAO). The
pictured prototype dish was manufactured by MMS Technology Ltd. in South
Africa (mmstechnology.co.za). A second metal dish prototype was constructed
jointly by NVJ and Rebcon Engineering, but is not discussed here.

Figure 8: A 6 m diameter f/D = 0.25 prototype dish at the Hartebeesthoek Radio
Astronomy Observatory in South Africa. The dish is aluminum embedded in
fiberglass composite, and the receiver is supported from the dish vertex by a
10 cm diameter fiberglass column. HIRAX will operate in drift scan, so the
mount has only an elevation axis, which allows for pointings between $\pm
30^{\circ}$ degrees off zenith. PI for scale.

Photogrammetry measurements were performed on the prototype dish in Figure 8
to assess the accuracy of the dish surface and to compare gravitational
deflections to those predicted by mechanical finite element analysis (FEA)
simulations. Figure 9 shows FEA simulations111Dish FEA simulations performed
by MMS Technology (Pty) Ltd. Figure courtesy of Heinrich Bauermeister. of the
f/D = 0.25 dish, and corresponding photogrammetry results222Photogrammetry
performed by the South African Radio Astronomy Observatory (SARAO). Figure
courtesy of Mattieu de Villiers. for zenith angles of $0^{\circ}$ and
$30^{\circ}$.

The photogrammetry data points are fit to a rotationally symmetric paraboloid,
and the plots show the data residuals with respect to this fit. The dominant
deformation is the quadrupolar “potato chip” mode, which arises from the
cross-shaped dish backing structure and increases in amplitude with zenith
angle. The FEA simulations confirm the quadrupolar pattern of the distortions
under gravitational load; however, the expected amplitudes are several times
smaller than those observed in measurements. Deformations at zenith pointing
are predicted to be less than 1 mm, but are measured to have a maximum
amplitude of 5 mm. At $30^{\circ}$ zenith angle the maximum amplitude
increases to 4 mm in the model and 8 mm in the measured data.

The deformations in both the zenith and $30^{\circ}$ pointing are also most
significant towards the edge of the dish. This may be an indication that a
more substantial backing structure is necessary to support the large diameter
of the dish.

Figure 9: HIRAX prototype dish photogrammetry and FEA simulations. The top row
shows FEA simulations with predicted total deformations, while the bottom row
shows the measured deviations from the theoretical paraboloid surface,
measured normal to the surface. The left column shows zenith pointing, the
right column shows $30^{\circ}$ elevation below zenith, the lowest pointing
for the HIRAX survey. The measured deformations, while close to the design
specifications, are several times larger than those predicted by the
simulations. Dish FEA simulations were performed by MMS Technology (Pty) Ltd.,
and figures are courtesy of Heinrich Bauermeister. Photogrammetry was
performed by the South African Radio Astronomy Observatory (SARAO), and
figures are courtesy of Mattieu de Villiers.

Figure 10 shows the dish surface displacements after the best-fit paraboloid
and quadrupole are subtracted at three different angles: $0^{\circ}$ (zenith),
$15^{\circ}$, and $30^{\circ}$. The residual displacements are less than 1 mm,
showing the dish surface shape is well described by the combination of the
ideal paraboloid and a quadrupolar mode. Figure 10 also shows the differences
between pairs of angles, which shows the systematic changes in the dish
surface as the elevation angle is adjusted. These differences are larger, with
maximum amplitudes of between one and two millimeters. This is further
evidence that the dish needs addition support structure, to reduce the
variations in dish shape as it is tilted.

This data informs our future design and production methods; while close to
acceptable limits, the dish surface needs to be more precise, and the backing
structure may need to be improved as well. The dish at zenith has a measured
RMS surface deviation of 2.8 mm, and the requirement is $<1$ mm RMS deviation
from the ideal paraboloid surface. We are confident that the mechanical
improvements required to achieve the design specifications are feasible. The
derivation of the specifications is discussed below.

Figure 10: Residual displacements after subtracting best-fit paraboloid and
quadrupole. Color scales are in millimeters. Top row (left to right):
$0^{\circ}$ (zenith), $15^{\circ}$, and $30^{\circ}$. Deviations from the
model are below 1 mm, showing that the surface is well described by a
paraboloid plus a quadrupolar distortion arising from the backing structure.
Bottom row, differences between residuals at angle pairs, showing the
systematic changes in the dish surface as the dish is tilted (left to right):
$0^{\circ}$ and $15^{\circ}$, $0^{\circ}$ and $30^{\circ}$, and $15^{\circ}$
and $30^{\circ}$.

##  5 Implications for 21 cm Cosmology

In order to examine the effects of beam shape systematics related to the
mechanical specification of the dishes, we have propagated visibility
perturbations derived from the CST simulations described above through to 21
cm power spectrum errors. The method for this is outlined briefly here but
will be more comprehensively described in future work.

Firstly, we construct a summarised representation of the deviations observed
in the beam shapes in the CST parameter sweep simulations outlined in Section
3.2. This is done by applying a principle component analysis (PCA) to the
observed radial profile deviations of the beam directivities exported from the
CST simulations (currently only deviations in the co-pol. beams are
considered). From this analysis, we select the 2 modes (functions of frequency
and radial angular coordinate) that encapsulate the most variance in the beam
radial profile for a given parameter in the dish-feed system that has been
swept. The appropriately weighted linear combination of these modes then
represents an approximation of the linear-order perturbation of the beam
shape, which can be expressed as a derivative with respect to the physical
parameter under consideration. With this linear approximation of the effect of
varying a given physical parameter on the beam shape, we can efficiently
generate realisations of perturbed beams motivated from what is observed in
the CST simulation parameter sweeps.

To evaluate a potential systematic’s effect on 21 cm cosmology, we generate
random realisations of the perturbed beams based on a distribution of physical
parameters across the dishes of the array. Figure 11 shows the effect of such
perturbations in the feed position along the focal axis relative to the dish,
$\delta z$, on the recovered 21 cm power spectrum. Here we have assumed the
distribution of $\delta z$ across dishes is that of a Gaussian with standard
deviation of 10 mm. Using the perturbed beams we first generate visibility-
space perturbations in simulated data assuming a nominal sky model including
galactic foregrounds and cosmological 21 cm signal, and a nominal instrument
and survey. These perturbed visibilities are then used to construct a
foreground filtered power spectrum estimate using a modified version of the
$m$-mode formalism [28, 29]. In Figure 11 a plot of the deviations of the
recovered power spectrum with respect to the input 21 cm power spectrum due to
these systematic perturbations is shown.

Currently this effort has only considered small scale simulations of the array
(limited to core baselines and subsets of the full frequency range), and
radially symmetric perturbations in the co-pol. beam shapes. Future work will
extend this to full scale simulations using a less summarised prescription for
the beam deviations as well as incorporating full polarization information.

These results have been used to inform the telescope specifications listed in
Table LABEL:tab:specs. For this purpose we first determined our tolerance on
systematic contributions to the required power spectrum sensitivity from high-
level Fisher matrix forecasts on Dark Energy parameter constraints. This
tolerance was set such that the impact of the systematic effects under
evaluation would reduce the Fisher-forecasted Dark Energy Figure of Merit
(FoM) by no more than 20% of the fiducial forecasted level that assumes a
power spectrum noise comprised of purely statistical noise with a nominal
residual foreground contribution. Table LABEL:tab:specs represent a
mechanically achievable set of specifications that is consistent with this
tolerance as informed by the simulations described here. A more quantitative
description of this process will be presented in future work.

Figure 11: Deviations in the recovered 21 cm power spectrum due to a single
realization of systematic beam perturbations due to Gaussian distributed
positional offsets (with standard deviation=10 mm) of the feed along the focal
axis, as calibrated by CST simulations. The results shown here are scaled by
the estimated statistical noise on the power spectrum bins such that values of
magnitudes greater than one would indicate systematic contributions larger
than the estimated statistical noise level. The simulations from which these
results were generated consider only the frequency range of 600-650 MHz and
utilises only a subset of the baselines available to HIRAX, but assumes the
same level of redundancy as in the 1024 element array. Upcoming publications
will expand upon these results.

##  6 Summary and Future Work

The simulations and verification measurements performed for this work have
significantly impacted the design of the HIRAX array. The cabling and receiver
support structure were completely re-designed to improve RF performance and
mechanical stability. Routing the cabling precisely along the telescope axis
of symmetry reduces sidelobe amplitudes by 7 dB from the previous prototype
design, and eliminates asymmetries in the sidelobe structure. Selecting a
large (35 cm) diameter receiver support column over a small (10 cm) support
column increases the main beam amplitude by 0.5 dB. The focal ratio was
adjusted, resulting in a significant decrease in crosstalk. Decreasing the
focal ratio from 0.25 to 0.23 reduces inter-dish reflections by up to 5 dB, as
measured by the delay kernel, and line-of-sight noise coupling ($S_{31}$) by
up to 10 dB, while only reducing aperture efficiency by $5\%$, from
approximately $60\%$ to $55\%$. We have also conducted extensive range
measurements of the HIRAX feed to verify the accuracy of our simulations, and
are preparing for field measurements of the full 6 m dish elements with a
drone based calibration platform.

We have traced the effects of the modeled system non-redundancies to recovered
cosmological parameters, and have derived specifications for the manufacture
and assembly of the dish from those cosmological simulations, as summarized in
Table LABEL:tab:specs. This is of particular importance for the 21 cm
cosmology field: this is the first time the redundancy requirements of a 21 cm
array have been explicitly derived from the telescope mechanical and
electromagnetic performance to the cosmology simulations pipeline.

Future simulations will continue to examine the effects of intra-dish and
inter-dish reflections on array performance, will expand the simulation models
of dish perforations and conductivity beyond the first-pass models used here,
and examine the effects of large scale (quadrupolar) and small scale (RMS)
deformations of the dish surface on the instrument beams and survey
cosmological constraints. Future publications will expand on the details of
the HIRAX electromagnetic simulations, beam verification measurements, and
cosmological simulations. There is currently a tender out for production of
the dishes for the HIRAX 256-element array, based in part on the
specifications derived here, and construction of the array will begin once the
tender is complete, in 2021.

###### Acknowledgements.

This work was supported by the National Science Foundation (NSF) under Grant
No. 1751763. AEW acknowledges support from the Berkeley Center of Cosmological
Physics. ERK is supported by a NASA Space Technology Research Fellowship
(NSTRF). KM acknowledges support from the National Research Foundation (NRF)
of South Africa. The HIRAX radio array is funded by the NRF. DC acknowledges
the financial assistance of the South African Radio Astronomy Observatory
(SARAO, www.sarao.ac.za). Results presented here were produced with the
support of the computing and technical resources of the Yale Center for
Research Computing and Wright Laboratory, and with the support of personnel
and facilities at NCSU, Caltech, and NASA JPL.

## References

  * [1] Weinberg, D. H., Mortonson, M. J., Eisenstein, D. J., Hirata, C., Riess, A. G., and Rozo, E., “Observational probes of cosmic acceleration,” Physics Reports 530, 87–255 (Sept. 2013).
  * [2] Bassett, B. A. and Hlozek, R., [Dark Energy ], ch. Baryon Acoustic Oscillations, Cambridge University Press (2010).
  * [3] Kogut, A. et al., “Three-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Foreground Polarization,” Astrophys. J. 665, 355–362 (2007).
  * [4] Eisenstein, D. J. e. a., “Detection of the Baryon Acoustic Peak in the Large‐Scale Correlation Function of SDSS Luminous Red Galaxies,” The Astrophysical Journal 633, 560–574 (Nov 2005).
  * [5] Percival, W. J. et al., “Measuring the Baryon Acoustic Oscillation scale using the Sloan Digital Sky Survey and 2dF Galaxy Redshift Survey,” Monthly Notices of the Royal Astronomical Society 381, 1053–1066 (Sep 2007).
  * [6] Abbott, T. M. C. e. a., “Dark Energy Survey Year 1 results: measurement of the baryon acoustic oscillation scale in the distribution of galaxies to redshift 1,” Monthly Notices of the Royal Astronomical Society 483, 4866–4883 (Dec 2018).
  * [7] Chang, T.-C., Pen, U.-L., Peterson, J. B., and McDonald, P., “Baryon Acoustic Oscillation Intensity Mapping as a Test of Dark Energy,” Phys. Rev. Lett. 100, 091303 (2008).
  * [8] The CHIME/FRB Collaboration, “Observations of fast radio bursts at frequencies down to 400 megahertz,” Nature 566, 230–234 (Jan 2019).
  * [9] The CHIME/FRB Collaboration, “A second source of repeating fast radio bursts,” Nature 566, 235–238 (Jan 2019).
  * [10] The CHIME/FRB Collaboration, “CHIME/FRB Detection of Eight New Repeating Fast Radio Burst Sources,” The Astrophysical Journal Letters 885(1) (2019).
  * [11] Fonseca, E. et al., “Nine New Repeating Fast Radio Burst Sources from CHIME/FRB,” The Astrophysical Journal 891, L6 (Feb 2020).
  * [12] The CHIME/FRB Collaboration, “Periodic activity from a fast radio burst source,” Nature 582, 351–355 (Jun 2020).
  * [13] Swetz, D. S. et al., “The Atacama Cosmology Telescope: The Receiver and Instrumentation,” The Astrophysical Journal Supplement Series 194 (2018).
  * [14] Benson, B. A. et al., “SPT-3G: A Next-Generation Cosmic Microwave Background Polarization Experiment on the South Pole Telescope,” Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 9153 (July 2014).
  * [15] Michael Levi and the DESI collaboration, “The DESI Experiment, a whitepaper for Snowmass 2013,” (2013).
  * [16] Sarah Brough et al., “The Vera Rubin Observatory Legacy Survey of Space and Time and the Low Surface Brightness Universe,” (2020).
  * [17] Gupta, N. et al., “The MeerKAT Absorption Line Survey (MALS),” (2017).
  * [18] Kuhn et al., “Noise temperature testing for the Hydrogen Intensity and Real-time Analysis eXperiment (HIRAX),” Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 11445 (Dec. 2020).
  * [19] Newburgh, L. B. et al., “HIRAX: A Probe of Dark Energy and Radio Transients,” Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 9906 (July 2016).
  * [20] Bandura, K. et al., “Canadian hydrogen intensity mapping experiment (chime) pathfinder,” Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 9145 (July 2014).
  * [21] DeBoer, David R. et al., “Hydrogen Epoch of Reionization Array (HERA),” Publications of the Astronomical Society of the Pacific 129, 045001 (Mar 2017).
  * [22] Deng, M., Campbell-Wilson, D., and the CHIME Collaboration, “The cloverleaf antenna: A compact wide-bandwidth dual-polarization feed for chime,” International Symposium on Antenna Technology and Applied Electromagnetics (ANTEM) 16, IEEE (July 2014).
  * [23] O’Connor, P. et al., “The Baryon Mapping Experiment (BMX), a 21cm intensity mapping pathfinder,” Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 11445 (Dec. 2020).
  * [24] King, Oliver G. et al., “The C-Band All-Sky Survey: instrument design, status, and first-look data,” Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 7741 (Jul 2010).
  * [25] Ewall-Wice, A. et al., “The Hydrogen Epoch of Reionization Array Dish. II. Characterization of Spectral Structure with Electromagnetic Simulations and Its Science Implications.,” The Astrophysical Journal 831, 196 (Nov. 2016).
  * [26] Shahpari, M. and Thiel, D. V., “The impact of reduced conductivity on the performance of wire antennas,” IEEE Transactions on Antennas and Propagation 63, 4686–4692 (Nov 2015).
  * [27] Vanderlinde, K. et al., “LRP 2020 Whitepaper: The Canadian Hydrogen Observatory and Radio-transient Detector (CHORD),” (2019).
  * [28] Shaw, J. R., Sigurdson, K., Pen, U.-L., Stebbins, A., and Sitwell, M., “All-sky Interferometry with Spherical Harmonic Transit Telescopes,” Astrophysical Journal 781, 57 (Feb. 2014).
  * [29] Shaw, J. R., Sigurdson, K., Sitwell, M., Stebbins, A., and Pen, U.-L., “Coaxing cosmic 21 cm fluctuations from the polarized sky using m -mode analysis,” Phys. Rev. D 91, 083514 (Apr. 2015).