plus1sp

[1]Gaurav Arya

# End-to-end Optimization of Metasurfaces for Imaging with Compressed Sensing

William F. Li  Research Laboratory of Electronics, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA  Charles Roques-Carmes  Ginzton
Laboratory, Stanford University, Stanford, CA 94305, USA  Marin Soljačić
Research Laboratory of Electronics and Department of Physics, Massachusetts
Institute of Technology, Cambridge, MA 02139, USA  Steven G. Johnson
Department of Mathematics, Massachusetts Institute of Technology, Cambridge,
MA 02139, USA  Zin Lin  Bradley Department of Electrical and Computer
Engineering, Virginia Tech, Blacksburg, VA 24061, USA

###### Abstract

We present a framework for the end-to-end optimization of metasurface imaging
systems that reconstruct targets using compressed sensing, a technique for
solving underdetermined imaging problems when the target object exhibits
sparsity (i.e. the object can be described by a small number of non-zero
values, but the positions of these values are unknown). We nest an iterative,
unapproximated compressed sensing reconstruction algorithm into our end-to-end
optimization pipeline, resulting in an interpretable, data-efficient method
for maximally leveraging metaoptics to exploit object sparsity. We apply our
framework to super-resolution imaging and high-resolution depth imaging with a
phase-change material. In both situations, our end-to-end framework discovers
optimal metasurface structures for compressed sensing recovery, automatically
balancing a number of complicated design considerations to select the best
imaging measurement matrix from a complex, physically-constrained manifold
with millions of dimensions. The optimized metasurface imaging systems are
robust to noise, significantly improving over random scattering surfaces and
approaching the ideal compressed sensing performance of a Gaussian matrix,
showing how a _physical_ metasurface system can demonstrably approach the
_mathematical_ limits of compressed sensing.

##  1 Introduction

In this paper, we combine a subwavelength flat-optics (``metasurface'')
physical platform [27, 30, 17, 10] with a compressed sensing (CS)
reconstruction algorithm [13] to end-to-end optimize the performance of
underdetermined imaging systems. Underdetermined imaging systems, where the
imaging sensor has fewer pixels than the target object, are becoming
increasingly important as researchers attempt to extract more information from
objects while keeping imaging systems compact, inexpensive, and fast. CS is a
powerful, principled way of solving underdetermined imaging problems using a
_sparsity_ prior, where the object (or some transformation of it) is known to
have a large number of nearly-zero values [19, 61, 11]. However, in CS
problems, the parameters of the physical platform (which determine the
``measurement matrix'' of the system) have traditionally been assumed to be
fixed.

Fig. 1: The bilevel optimization problem for a metasurface imaging system in
an underdetermined and potentially multi-channel imaging setting. A compressed
sensing minimization problem is _nested_ within the optimization pipeline. We
find the optimal measurement matrix that is physically realizable by our
metasurface optics, along with the best corresponding reconstruction
hyperparameters, through end-to-end optimization of the reconstruction error.

In this work, we show how to ``end-to-end'' optimize the metasurface imaging
system to find the _physically realizable_ measurement matrix that minimizes
the CS reconstruction error. Thus, we propose a framework (Fig. 1, Section
2.1) for answering the following question: how can we maximally leverage
metaoptics to exploit object sparsity? We embed an unapproximated, iterative
$\ell_{1}$-regularized CS algorithm into our end-to-end pipeline, performing
an adjoint sensitivity analysis of the Karush-Kuhn-Tucker (KKT) optimality
conditions of the reconstructed CS solution for efficient backpropagation
(Section 2.2). Like any classical CS algorithm, the reconstruction process is
interpretable, data-efficient, and generalizable to data outside the training
distribution: the _only_ assumption made is sparsity, and no additional data
is needed to formulate and solve the optimization problem. This distingishes
our approach from previous end-to-end metasurface [35, 59, 36, 17] and
diffractive-optics [51, 5, 40, 54, 60, 14] optimization that considered
overdetermined problems where CS was not necessary and/or took a data-
intensive deep-learning approach to image reconstruction.

Our end-to-end optimization forms composite lens-like metasurface
configurations that dramatically improve upon random scattering surfaces. For
example, for pixel-limited imaging of sparse objects with a small sensor, we
discover computational metasurface designs that approach the theoretically
optimal performance of physically unrealizable Gaussian measurement matrices
(Fig. 3). Importantly, metasurfaces exploit full-wave electromagnetic
physics—such as multiple scattering, strong wavelength/polarization
dispersion, and resonances [38, 20]. Even greater inference capabilities can
potentially be unlocked by phase-change materials (PCMs) [50, 48, 24], which
allow one to collect more data with a single device. We show that our end-to-
end framework can maximize complex electromagnetic interactions in PCMs to
tackle severely underdetermined multi-channel imaging problems (Section 3.2)
such as high-resolution 3D imaging. Our work represents a new paradigm for
enabling _maximal symbiosis_ between metaoptics and CS. Our framework is also
immediately applicable to diffractive optics, which are simply a special
limiting case of metasurfaces (when the surface is locally uniform on the
wavelength scale). These results open up new possibilities for high-resolution
depth, spectral, and polarization imaging and beyond (Section 4). For example,
in subsequent work we have applied our end-to-end CS framework to design 3D
volumetric non-paraxial metasurfaces for detecting angles and wavelengths of
sparse incoming beams, yielding order-of-magnitude improvements in resolution
$\times$ size [32].

Our metasurface physics platform is an ultra-thin nanostructured interface
composed of an aperiodic array of subwavelength nanopillars, similar to
structures previously fabricated for various metaoptics applications [27, 42,
35, 28, 33, 16]. In contrast to diffractive optical elements, which can be
modeled by scalar diffraction theory [12], metasurfaces consist of nanoscale
elements which require the simulation of ``full-wave'' Maxwell's equations.
This results in enhanced dispersion, which can be used to realize multi-
functional devices [16, 45, 44]. Furthermore, metasurfaces can be fabricated
over large scales with chip-scale fabrication techniques in various mature
(e.g., titania [18] and silicon [4]) and emerging (e.g., phase-change
materials (PCM) [50, 48, 24]) material platforms. The large number of degrees
of freedom (up to $\gtrsim 10^{6}$ in this paper) and the rich physics at the
nanoscale have spurred an emerging field of large-scale metasurface inverse
design [34, 41, 43]. However, most metaoptics designs have been confined to
optics-only applications, such as lensing, beam-steering, and holography [15].
In those applications, the desired image is prescribed _a priori_. In
contrast, recent work has introduced _end-to-end_ metaoptics inverse design
[36, 35, 59], whereby metasurfaces are directly optimized to minimize the
reconstruction error in overdetermined imaging settings, computationally
discovering the best sensor image for subsequent reconstruction. This paper
extends end-to-end metasurface optimization to the underdetermined regime,
where sensor data is limited, for both single-channel (2D) imaging (Fig. 3)
with a small sensor and multi-channel (3D) high-resolution depth imaging with
a phase-change material (Section 3.2).

Compressed sensing (CS) has become increasingly prevalent in computational
imaging, whereby optical components are paired with a CS reconstruction
algorithm, which in some sense finds the ``sparsest'' object that produces the
observed image. A key contribution of our work is to perform end-to-end
optimization with a _nested_ CS minimization problem that is solved in a fully
iterative manner (not approximated, unrolled, or truncated), as opposed to
optimizing proxies for good CS performance [62, 47] or employing the
approximation of a differentiable unrolled network (a fixed number of
iterations) [23, 63, 40]. We solve the full bilevel optimization problem by
performing a sensitivity analysis of the KKT conditions [9, 2] associated with
the CS convex-optimization problem, using an adjoint formulation to
efficiently backpropagate the gradient (Section 2.2). Crucially, we extend
previous work that performed hyperparameter optimization of CS [7] by also
optimizing with respect to the metasurface parameters underlying the
measurement matrix itself.

##  2 Method

###  2.1 The bilevel optimization problem

Our pipeline consists of a physics platform, which simulates light propagation
through the metasurface to a grayscale sensor, and a compressed-sensing
reconstruction algorithm, which solves a convex optimization problem (Fig. 1).
We optimize this pipeline end-to-end by minimizing the expected relative mean-
square error (MSE) of the reconstructions over a distribution of training
objects. Given a training object $u$ (flattened into a vector), a raw noisy
image $y$ is formed by:

$y=G(p)u+\eta,\leavevmode\nobreak\ \eta\sim\mathcal{N}(0,\sigma),$ (1)

where $G$ is the measurement matrix of the imaging system, dependent on the
metasurface geometry $p$ (see below), and $\eta$ is an additive Gaussian noise
[51]. We now write a bilevel optimization problem for our end-to-end pipeline:

$\displaystyle\min_{p,\alpha,\beta}\mathcal{L}=\Bigg{\langle}{|u-u_{\text{est}}|_{2}^{2}\over|u|_{2}^{2}}\Bigg{\rangle}_{u,\eta},$
(2) $\displaystyle
u_{\text{est}}=\operatorname{arg\,min}_{x}|Gx-y|_{2}^{2}+\alpha|\Psi
x|_{1}+\beta|x|_{2}^{2}.$ (3)

The reconstructed object $u_{\text{est}}$ is formed by solving a generalized
Lasso problem [56] with an additional ElasticNet $\ell_{2}$ term [65] (which
accelerates convergence as explained below). Here,
$\langle\leavevmode\nobreak\ \rangle_{u,\eta}$ denotes averaging over multiple
training objects and noise realizations, and $\Psi$ is a _sparsifying_
transformation under which we expect the object to have a large number of
coefficients that are nearly zero. The $\ell_{1}$ regularization term $|\Psi
x|_{1}=\sum_{j}|(\Psi x)_{j}|$ reflects this prior by encouraging sparsity of
the vector $\Psi x$ [13]. In the imaging settings we treat in this paper,
$\Psi$ is the identity operator ($\ell_{1}$ regularization), but we also
develop our sensitivity analysis in Section 2.2 for the 2D/3D anisotropic
gradient operator with periodic boundary conditions—total variation (TV)
regularization [46]—to show the generality of our framework for arbitrary
linear transformations $\Psi$. We solve Eq. 3 using the matrix-free fast
iterative soft-thresholding algorithm (FISTA) [6, 26], which allows us to
exploit the convolutional structure of the large matrix $G$ [3, 35].

We work in the paraxial regime, where the response of the physics platform is
characterized by a set of point spread functions (PSFs)—the 2D image of a
point source on the optical axis—one for each depth and/or wavelength channel
[22]. The computation of the PSFs consists of spherical wave propagation,
followed by a simulation of light scattering at the metasurface, followed by a
near-to-far field transformation to compute the electric field in the sensor
plane. To model scattering at the metasurface, we use the locally periodic
approximation [43], in which the large-area metasurface is locally
approximated by sub-wavelength _periodic_ unit cells, each containing a
nanopillar. We use rigorous coupled-wave analysis (RCWA) [37] to simulate the
complex-valued transmission coefficients of a library of unit cells, from
which we construct Chebyshev polynomials that interpolate the transmission
coefficients as a function of the nanopillar geometrical parameters [43]. This
differentiable Chebyshev surrogate model allows for rapid simulation of
metasurface scattering during optimization.

From these PSFs, we assemble the measurement matrix $G$. The action of $G$ is
to convolve each channel of the object with its corresponding PSF, sum up the
resulting contributions to the electric field intensity in the sensor plane,
and then crop the field onto the sensor's area. Here, the number of sensor
pixels is smaller than the number of pixels in the (potentially multi-channel)
object, and thus $G$ has fewer rows than columns. We add independent Gaussian
noise to the sensor readings, which approximates Poisson shot noise in the
limit of large photon flux. The standard deviation $\sigma$ of the Gaussian
noise $\eta$ is chosen to be a few percent of the mean image intensity, i.e.
proportional to $|Gu|_{1}$.

In our bilevel formulation, the choice of ElasticNet reconstruction
($\beta>0$) instead of pure Lasso ($\beta=0$) serves a crucial purpose
_during_ the optimization. The initial metasurface geometry leads to a
measurement matrix $G$ that is highly unsuitable for compressed sensing; with
no $\ell_{2}$ term, the number of FISTA iterations required for convergence
can be orders of magnitude larger for unoptimized metasurface geometries than
for optimized geometries. Intuitively, for poorly performing $G$, nearby
pixels are difficult to distinguish, so FISTA spends a large number of
iterations making an arbitrary, badly conditioned decision about which non-
zero entries to include. With the $\ell_{2}$ term, the algorithm instead
``gives up'' on this futile task and spreads the values across many entries.
In practice, we find that a large initial $\beta$ makes the number of nested
iterations roughly _independent_ of $G$, fixing the slow-start issue. As we
approach the optimized matrix, the need for the $\ell_{2}$ term diminishes and
the optimization chooses to shrink $\beta$ to $\approx 0$, recovering the pure
Lasso (see Fig. 3).

###  2.2 Sensitivity Analysis of Compressed Sensing Reconstruction

For large-scale optimization ($\sim 10^{6}$ parameters) to be feasible, we
must compute the sensitivity of the error $\mathcal{L}$ with respect to the
parameters $p,\alpha$, and $\beta$. An adjoint method allows us to
backpropagate gradients through our pipeline, at a computational cost similar
to that of $\mathcal{L}$. For the majority of the pipeline, this can be done
effortlessly using automatic differentiation tools [25].

Fig. 2: The formation of the projection operator $P_{S}$, given the Lasso
solution found in the forward pass. For simplicity we depict $6\times 6$ 2D
objects, where different colors represent different intensity values, in the
cases of both $\ell_{1}$ and TV regularization. The support of $\Psi
u_{\text{est}}$ is denoted in red (boxes around non-zero pixels for $\ell_{1}$
regularization, arrows representing non-zero differences for TV
regularization). The operator $P_{S}$ is then used to set up the adjoint
system, which is solved in the backward pass.

However, differentiating through Eq. 3 requires special handling. Since we
compute the reconstruction $u_{\text{est}}$ (a vector of length $\approx
10^{5}$ in our optimizations) through an iterative solver, directly
backpropagating through the solver's iterations would require a prohibitive
amount of memory. Furthermore, due to the non-differentiability of the
$\ell_{1}$ term $|\Psi u|_{1}$, the reconstruction $u_{\text{est}}$ is itself
not strictly differentiable with respect to $G$, $\alpha$, and $\beta$.

To surmount these hurdles, our starting point is the KKT conditions for the
generalized Lasso [57]. We expect that $\Psi u_{\text{est}}$ is sparse; let
$S$ be the support of $\Psi u_{\text{est}}$ (the set of indices of non-zero
values). We can then define $P_{S}$, the orthogonal projection operator onto
the subspace of objects $x$ that satisfy $\operatorname{support}\Psi
x\subseteq S$. When $\Psi$ is the identity ($\ell_{1}$ regularization),
$P_{S}$ restricts an object $x$ to the indices in $S$, zeroing out the other
values. When $\Psi$ is the anisotropic gradient operator (TV regularization),
$P_{S}$ sums an object $x$ over each of the piecewise constant regions induced
by the sparse gradient $\Psi u_{\text{est}}$, divides each region's sum by the
square root of the size of the region (the division ensures $P_{S}$ is an
orthogonal projection, i.e. $P_{S}=P_{S}^{T}$), and sets all elements in that
region equal to this summed and normalized value. These two cases are
illustrated in Fig. 2. The KKT conditions then imply that

$\underbrace{(P_{S}G^{T}GP_{S}^{T}+\beta
I)}_{A}\underbrace{u_{\text{est}}}_{x}=\underbrace{P_{S}G^{T}y-\frac{1}{2}\alpha
P_{S}\Psi^{T}\operatorname{sgn}{\Psi u_{\text{est}}}}_{b},$ (4)

which is a square linear system with dimension equal to the rank of $P_{S}$.
Intuitively, the hard work performed in solving Eq. 3 by convex optimization
methods is to discover the support $S$. Once $S$ (and the signs
$\operatorname{sgn}{\Psi u_{\text{est}}}$) of the elements of $\Psi
u_{\text{est}}$) have been found in the forward pass, we know a subgradient of
the non-differentiable $\ell_{1}$ term $|\Psi u|_{1}$ at $u_{\text{est}}$.
This information allows us to form an efficient backward pass.

Specifically, we can use Eq. 4 to compute the sensitivity of $u_{\text{est}}$.
Given the gradient $\frac{\partial\mathcal{L}}{\partial u_{est}}$, we write
$\frac{\partial\mathcal{L}}{\partial
x}=P_{S}\frac{\partial\mathcal{L}}{\partial u_{\text{est}}}$. Given
$\frac{\partial\mathcal{L}}{\partial x}$, we now apply the standard adjoint
method for a linear system $Ax=b$ [53] to find
$\frac{\partial\mathcal{L}}{\partial A}$ and
$\frac{\partial\mathcal{L}}{\partial b}$ by solving a single linear system
with left operator $A^{T}$ (the adjoint system). Since
$A=P_{S}G^{T}GP_{S}^{T}$ is symmetric, we can solve this system using the
matrix-free conjugate-gradient method [49]. By backpropagating through the
formation of $A$ and $b$ from $G,\alpha$ and $\beta$, we then find the desired
sensitivities $\frac{\partial\mathcal{L}}{\partial G}$,
$\frac{\partial\mathcal{L}}{\partial\alpha}$, and
$\frac{\partial\mathcal{L}}{\partial\beta}$. In this manner, we can
backpropagate the gradient through Eq. 3 and thus the entire pipeline. An
important aspect of our algorithm is that it is ``matrix-free'' [31]: in
metasurface optics, the measurement matrix is an enormous convolution with a
Green's function [22] that would be impractical to construct explicitly, but
our algorithm is constructed to only employ the matrix implicitly as a linear
operator (exploiting FFT-accelerated convolutions [52]).

We remark that the solution to the generalized Lasso is unique almost surely
for the cases of $\ell_{1}$ and TV regularization [55, 1], and hence $A$ is
non-singular. In our optimizations, we find the adjoint system to be well-
conditioned (it is strictly easier than the convex-optimization forward
problem). It is important to note that $\frac{\partial\mathcal{L}}{\partial
G}$, $\frac{\partial\mathcal{L}}{\partial\alpha}$, and
$\frac{\partial\mathcal{L}}{\partial\beta}$ are only _subgradients_ , as the
support of $\Psi u_{\text{est}}$ can change. However, the support of $\Psi
u_{\text{est}}$ is locally constant almost everywhere [57]; in practice, we
observe excellent results when plugging the subgradient
$\frac{\partial\mathcal{L}}{\partial p}$ into standard gradient-based
optimization algorithms. We employ the stochastic-gradient algorithm Adam [29]
because our training objects $u$ are drawn from a random distribution of
sparse objects.

##  3 Results

###  3.1 Single-channel imaging with small sensor

Fig. 3: Metasurface 2D imager with compressed sensing. (a) The optimized point
spread function, a multi-focal metalens with 27 foci (8x binned to emphasize
lens positions). The relative size of the small sensor is shown; the multi-
focal design spans beyond it so that the edges of the larger object are
captured on the sensor. (An animation of the optimization is provided as
supplementary material.) (b) Ground truth $256\times 256$ sparse object (7%
sparsity) with wavelength $550\mathrm{nm}$. (c) Image formed from ground truth
object by the optimized design, with additive Gaussian noise (2%). (d) CS
reconstruction of the ground truth object from the noisy image. (e)
Reconstruction error on _unphysical objects_ whose non-zero values are drawn
from a distribution with mean 0 as a function of object sparsity. Shown for
the optimized design, a random design, and the ideal Gaussian baseline under
the same noise setting (5%). (f) Evolution of $\alpha$ during optimization.
(g) Evolution of $\beta$ during optimization. (h) Convergence of the
reconstruction error during training on random objects with $5\%$ sparsity.
The objects have 3250 positive values uniformly sampled from the volume, with
intensity values uniformly drawn from the interval $[0.8,1.2]$. We also plot
the error of a random metasurface. (i) The $1.12\mathrm{mm}\times
1.12\mathrm{mm}$ TiO2 metasurface (refractive index $n\approx 2.4$) on a
silica substrate, consisting of a grid of $2048\times 2048$ sub-wavelength
square unit cells, where each unit cell has width $\approx 470\mathrm{nm}$.
The square pillar in each cell has a height of $600\mathrm{nm}$, and width
ranging from $\approx 60\mathrm{nm}$ to $\approx 410\mathrm{nm}$. The
metasurface is depicted as a heatmap of the pillar widths, and we also zoom in
on a $0.06\times 0.06\mathrm{mm}$ section of the metasurface.

We first apply our optimization technique to single-channel imaging of sparse
objects in the visible regime. We consider the case of a _small sensor_ that
has fewer pixels ($128\times 128$) than the object we wish to reconstruct
($256\times 256$). This is an underdetermined situation in which a traditional
lens would perform very poorly because the object would not fit on the sensor.
However, working under the assumption that our object is sparse in the spatial
domain (thus we set $\Psi$ to the identity), one may hope for a
``multiplexing'' design that compresses the object onto the small sensor in
such a way that the image can be deconvolved effectively through CS. It is
unclear _a priori_ what metasurface geometry would best realize this goal.

As shown in Fig. 3, the end-to-end optimization _spontaneously configures_ an
arrangement of lens-like focusing elements, leading to a number of sharp peaks
in the far field (Fig. 3). Over 250 iterations of optimization on training
objects with 5% sparsity (the proportion of non-zero values), the
metasurface—initialized as an array of identical nanonpillars—evolves into the
structure shown in Fig. 3. The relative RMSE (the square root of the relative
MSE $\mathcal{L}$) converges to $7.8\%$ (Fig. 3). This is an order of
magnitude improvement over a randomly structured metasurface (96% relative
RMSE).

Effectively, the optimized metasurface system takes many copies of a sparse
object (Fig. 3 shows a test object with 7% sparsity) and ``interlaces'' them
on the smaller sensor (Fig. 3). The object is faithfully recovered (18%
relative RMSE, Fig. 3) even though the small sensor image is further corrupted
by noise (2%). Our optimization also automates the choice of the
reconstruction-algorithm hyperparameters $\alpha$ and $\beta$. In particular,
$\alpha$ varies significantly on a logarithmic scale (Fig. 3) and settles at
approximately $10^{5}$, while $\beta$ shrinks to a negligible value that is
orders of magnitude smaller than $\alpha$ (Fig. 3).

We now discuss how the measurement matrix of our optimized system compares to
the ``gold standard'' of compressed sensing, a Gaussian random matrix. Such a
matrix is _unphysical_ for our optical system, but by comparing to it we can
characterize the performance of our system relative to _all_ possible
matrices, whether physical or not. Representing Gaussian matrices would
require a prohibitive amount of memory, so we use partial Gaussian circulant
matrices as a computationally tractable equivalent [58]. As shown in Fig. 3, a
random metasurface significantly underperforms a Gaussian matrix. Remarkably,
however, the optimized metasurface matches the performance of the Gaussian
matrix throughout the entire error-sparsity ``phase transition'' [8] (the
upswing segments of the error-sparsity relations in Fig. 3).

We qualify this apparent optimality result by noting that the comparison in
Fig. 3 is made with _unphysical_ objects. (In the following, "unphysical"
objects have mean zero and can have negative values; "physical" objects have
only non-negative values.) On physical objects, optical systems are less
robust to noise. This is due to the non-negativity of the measurement matrix
of the metasurface system, which is a fundamental optical limitation.
Intuitively, in the case of physical objects and measurement matrices, light
can only be added to the sensor and not subtracted, making it harder to form
sharp noise-robust features. This limitation is not faced by the Gaussian
measurement matrices used as benchmarks in Fig. 3, which have negative values
and are therefore unphysical (not realizable with optics).

Mathematically, recalling our noise model (Section 2.1), the worse performance
on physical objects manifests as a relative increase in the mean intensity
magnitude of the image relative to unphysical objects (which include negative
``intensities''), quantified by the following ratio:

$\text{Image Mean Gap of
}G=\frac{\langle|Gu|_{1}\rangle_{\text{phys.}}}{\langle|Gu|_{1}\rangle_{\text{unphys.}}}\cdot\frac{\langle|Xu|_{1}\rangle_{\text{unphys.}}}{\langle|Xu|_{1}\rangle_{\text{phys.}}},$
(5)

where $G$ is the measurement matrix of the metasurface system, $X$ is a
Gaussian matrix, and $\langle\cdot\rangle$ denotes averaging over particular
distributions of physical and unphysical objects. The Image Mean Gap
characterizes the ability of the metasurface system to form sharp noise-robust
image features on _physical_ non-negative objects. From another perspective,
the Image Mean Gap measures the ``noisiness'' of the PSF, i.e. its deviation
from a sharply focused lens $G=I$.

We can use the Image Mean Gap to link the results of Fig. 3 to the physical
case. We consider our training object distribution, where objects have $5\%$
sparsity and their non-zero values are drawn uniformly from the interval
$[0.8,1.2]$. In this setting, the Image Mean Gap is $\approx 5.4$ for a random
metasurface and $\approx 2.1$ for the optimized structure. This poorer Image
Mean Gap, combined with the random metasurface's poorer performance on
unphysical objects in Fig. 3, explains the random metasurface's very poor
performance on physical objects (96% relative RMSE). The Image Mean Gap of
$\approx 2.1$ ($2\times$ improvement) for the optimized structure is
significantly better: it implies our metasurface system performs equivalently
to an _unphysical_ Gaussian matrix system that has $\approx 2.1\times$ more
noise. It is unclear if the Image Mean Gap can be reduced even further while
maintaining the design's optimality on unphysical objects. One might hope that
a reduced number of foci in the design, leading to less overlap of the object
``copies'' on the sensor, would lead to sharper image features on non-negative
objects and thus a better Image Mean Gap. But this could come at the expense
of optimality on unphysical objects due to reduced multiplexing of the object
on the sensor, illustrating a tradeoff between the two factors. We leave more
refined lower bound analyses for future work.

For physical imaging situations, we emphasize the importance of sparsity _of
the measurement matrix_ , which corresponds to lens-like PSFs in the
convolutional model. Indeed, under a non-negativity constraint, one can show
that _sparse_ random matrices can perform well [21], and furthermore that
convolution with certain random sparse patterns can lead to excellent
performance [39]. Without any human input, the end-to-end optimization forms a
_sparse_ PSF that fits this characterization with optimal choices of the
number and positions of lenses. In doing so, the optimization forms a
measurement matrix with a reduced Image Mean Gap ($\approx 2.1$) and a near-
exact match with a Gaussian matrix on unphysical objects, therefore finding a
_physical_ metasurface geometry that significantly improves over random
scattering patterns.

###  3.2 High-resolution 3D imaging with PCM

Fig. 4: Design of GSS4T1 3D imager with compressed sensing. (a) Optimized
point spread functions, which exhibit high-frequency features, depth
sensitivity, and phase change contrast (4x binned to emphasize lens
positions). (b) The $0.77\mathrm{mm}\times 0.77\mathrm{mm}$ GSS4T1 metasurface
on a CaF2 substrate, consisting of a grid of $256\times 256$ sub-wavelength
square unit cells,. Each unit cell has width $\approx 3\mathrm{\mu m}$ and
index phase contrast $\Delta n=1.40$. The square pillar in each cell has a
height of $1.1\mu\mathrm{m}$, and width ranging from $\approx
0.29\mu\mathrm{m}$ to $\approx 2.9\mu\mathrm{m}$. The metasurface is depicted
as a heatmap of the pillar widths. (c) Zoom-in of PSFs at adjacent depths,
illustrating the depth sensitivity and phase change contrast.

To show the flexibility of our framework, we now consider a multi-channel
64-depth imaging problem in the mid-infrared regime. In particular, we
endeavor to image $\ell_{1}$-sparse, volumetric distributions of point sources
at very fine depth resolutions. This leads to a severely underdetermined
problem where an object with discretization $32\times 32\times 64\textrm{
depths}$ is captured on a $96\times 96$ sensor. To alleviate this problem, we
incorporate a _reconfigurable_ material (PCM), GSS4T1, which can undergo a
rapid phase change from amorphous to crystalline states via electrical
switching [64]. Thus, a GSS4T1 metasurface allows us to obtain two images from
a single structure design with two material phases of refractive index
contrast $\Delta n=1.40$ (Fig. 5). This presents a fascinating problem for
end-to-end optimization: the fineness of the required depth resolution
necessitates high depth sensitivity, and the metasurface design must also
ensure contrast between the images corresponding to the two material phases
(but the same geometry) so that the second shot provides useful additional
information. Our two-shot scheme using a PCM is significantly faster than
multi-shot schemes that use digital micromirror devices or spatial light
modulators [61].

In this setting, the optimization produces the metasurface shown in Fig. 4.
The PSFs of the optimized metasurface exhibit high-frequency lens-like
features (Fig. 4), necessary for good lateral resolution. Moreover, the PSFs
demonstrate spatial dispersion and phase change contrast (Fig. 4) despite
significant physical challenges. For instance, the thickness of the GSS4T1
metasurface is a tenth of the wavelength in the mid-infrared regime, limiting
the amount of wave scattering and dispersion available at the interface. Over
the depth range, we notice dramatic evolution in the point-spread function, as
various lenses of the multi-focal design come in and out of focus. However,
physical limitations prevent the design from achieving such dramatic shifts
between every pair of adjacent depths. Instead, the optimization tends to
judiciously shift lens positions between adjacent depths by a few pixels,
thereby shifting the images of adjacent channels on the sensor sufficiently to
achieve the desired $z$-resolution. Crucially, the optimized design achieves
these goals for the images of both configurations of the GSS4T1 metasurface.
Furthermore, it achieves high contrast between configurations for the same
depth channel, thus collecting _complementary_ information from the two shots.

Fig. 5: Analysis of GSS4T1 3D imager with compressed sensing. (a) A 64-depth
sparse ground truth test object ($\approx 1.5\%$ sparsity, or 1000 non-zero
values) with wavelength $3.2\mu\mathrm{m}$, and the depth channels inversely
spaced from $1.65\mathrm{mm}$ to $6.6\mathrm{mm}$. Patterns of varying
structure are ``hidden'' in three of the 64 channels, plotted in grayscale.
(b) The two images formed of the ground truth object by the amorphous and
crystalline metasurface states of the optimized two-shot structure, with
additive Gaussian noise (1%). (c) CS reconstruction of the ground truth object
for the optimized two-shot system (18% relative RMSE). (d) CS reconstruction
of the ground truth object for a two-shot system with a random metasurface
(59% relative RMSE). (e) CS reconstruction of the ground truth object for an
optimized single-shot system (74% relative RMSE). (f) Reconstruction error on
_unphysical objects_ with mean 0 as a function of object sparsity for the
optimized design, a random design, an ideal Gaussian baseline, and the
corresponding single-shot designs under the same noise setting (5%). (The
single-shot end-to-end design is directly optimized for the single shot case,
and the single-shot Gaussian baseline has half as many rows as the two-shot
matrix.) (g) GSS4T1 metasurface pillar in amorphous ($n=3.25$) and crystalline
($n=4.65$) states.

While these qualitative considerations give us a clue of the 3D imaging
system's resolving power, we now analyze the system's CS performance on
complex sparse scenes. In Fig. 5, we plot the reconstruction error as a
function of object sparsity on _unphysical objects_ for the end-to-end
optimized design, a random metasurface, and a partial Gaussian circulant
baseline. With 5% noise, we observe that the optimized structure closely (but
not perfectly) matches the performance of an ideal Gaussian matrix, while
substantially outperforming random metasurfaces (e.g. $\approx 3\times$
smaller RMSE at 3% sparsity). Now, as in Fig. 3, we study the performance gap
due to the non-negativity of the imaging system. At 3% object sparsity, the
Image Mean Gap (Eq. 5) between unphysical objects and non-negative training
objects with intensity values in the range $[0.8,1.2]$ is 46 for a random
structure and 18 for an optimized structure. Once again, the Image Mean Gap is
worse for a random structure ($\times 2.5$), showing how the end-to-end
optimization (which is trained on physical non-negative objects) favors PSFs
that are more lens-like and lead to sharper image features. The lens-like
features of the optimized PSFs (Figs. 4 and 4) are not perfectly focused,
however. This is reflected in the Image Mean Gap of 18 for the optimized
structure, which is still relatively large; we expect that incorporating more
complex nanophotonic structures than a single-layer metasurface will alleviate
this issue (Section 4).

We also perform a visual test of the system using a physical ground truth
object with $\approx 1000$ non-zero values ($\approx 1.5\%$ sparsity) where
patterns of varying structure have been hidden at 3 of the 64 depths. This
ground truth object differs substantially from the objects in the training
distribution which have uniformly distributed non-zero values across all
depths. Nevertheless, the optimized imaging system accurately recovers the
support and forms a reconstruction with $18\%$ relative RMSE (Fig. 5), despite
the image being corrupted with $2\%$ noise, showing how the system generalizes
to out-of-distribution sparse data. In comparison to the optimized system, a
random metasurface system achieves a much less accurate reconstruction of the
same object (59% relative RMSE, Fig. 5), and fails to accurately recover the
support. We also note the significant improvement over _single-shot_ designs
where the GSS4T1 metasurface is kept in only its crystalline state. On the
same ground truth object with $3\%$ sparsity (Fig. 5), an end-to-end optimized
single-shot design obtains a significantly worse reconstruction (74% relative
RMSE, Fig. 5). This is due in large part to the delayed sparsity-error phase
transition for the two-shot design (the rightward shift and slower incline of
the sparsity-error relations in Fig. 5), emphasizing the crucial role of the
material index contrast and the ability of our end-to-end optimization to
exploit it.

##  4 Conclusion and Outlook

Our results show how end-to-end optimization can navigate complicated,
physically-constrained optimization spaces to discover measurement matrices
with improved compressed-sensing performance. In certain cases where there are
enough physical degrees of freedom for a given task (such as 2D imaging with a
small sensor), the optimization discovers a measurement matrix that approaches
the performance of a Gaussian matrix for _noise-robust_ compressed sensing.
Even under severe physical constraints as in the case of a deeply sub-
wavelength PCM interface, end-to-end optimization can balance a complicated
set of considerations to obtain measurement matrices which significantly
outperform random two-shot designs and optimized single-shot designs for high-
resolution 3D imaging. We expect that our technique—which, in generality,
allows one to pair any measurement system with an $\ell_{1}$-regularized CS
minimization algorithm and optimize the system _end-to-end_ —will find use in
a wide range of CS recovery systems whose design involves a large number of
degrees of freedom.

In the design of such systems, it is crucial to take into account the type of
prior imposed on the object. We note that our end-to-end framework supports
any linear sparsifying transformation $\Psi$; we need only be equipped with
the projection operator $P_{S}$ to perform the sensitivity analysis. Future
work may build upon our method by performing end-to-end optimization of
physical systems that incorporate more complicated non-linear priors that
cannot be straightforwardly put in the form of the generalized Lasso.

Looking forward, as one leverages richer physics to extract more and more
information from a scene with limited sensor data, we anticipate a growing
need to cope with severely underdetermined problems. Our end-to-end framework,
which supports arbitrary multi-channel (and multi-configuration) settings,
paves the way for underdetermined nanophotonics imaging systems that can
extract hyperfine depth, spectral, and/or polarization information from a
single shot. The ultra-compact all-dielectric metasurface designs we have
presented are directly amenable to large-scale fabrication platforms [33] and
may enable the realization of metasurface cameras for imaging sparse scenes
(e.g., locating stars in the night sky). We are particularly excited about
future applications of our technique to thermal imaging, where one may
endeavour to construct a full spectral imager that can reconstruct both
emissivity and temperature profiles from a single shot. Furthermore, by
incorporating multi-layered metasurfaces and photonic crystal slabs, one may
exploit even richer _non-local_ physics and develop imaging systems whose
measurement matrix differs substantially in structure from a simple
convolutional model. In this case, there are few conventional principles to
guide the design; it is our hope that end-to-end optimization with compressed
sensing will allow one to exploit the full power of nanophotonics in
challenging underdetermined imaging problems.

Research Funding: This work was supported in part by the Simons Foundation and
by the U.S. Army Research Office through the Institute for Soldier
Nanotechnologies under Award No. W911NF-18-2-0048.

Conflict of interest statement: The authors declare no conflicts of interest
regarding this article.

Data availability: Data underlying the results presented in this paper may be
obtained from the authors upon reasonable request. The code used is publicly
available at https://github.com/gaurav-arya/ImagingOpt.jl.

## Bibliography

  * Ali, Tibshirani [2019] Ali Alnur, Tibshirani Ryan J.  The Generalized Lasso Problem and Uniqueness  // Electronic Journal of Statistics. 2019. 13, 2. 2307 – 2347.
  * Amos, Kolter [2017] Amos Brandon, Kolter J Zico.  Optnet: Differentiable optimization as a layer in neural networks  // International Conference on Machine Learning. 2017. 136–145.
  * Antipa et al. [2018] Antipa Nick, Kuo Grace, Heckel Reinhard, Mildenhall Ben, Bostan Emrah, Ng Ren, Waller Laura.  DiffuserCam: Lensless single-exposure 3D imaging  // Optica. Jan 2018. 5, 1. 1–9.
  * Arbabi et al. [2015] Arbabi Amir, Horie Yu, Bagheri Mahmood, Faraon Andrei.  Dielectric metasurfaces for complete control of phase and polarization with subwavelength spatial resolution and high transmission  // Nature nanotechnology. 2015. 10, 11. 937–943.
  * Baek et al. [2021] Baek Seung-Hwan, Ikoma Hayato, Jeon Daniel S, Li Yuqi, Heidrich Wolfgang, Wetzstein Gordon, Kim Min H. Single-shot Hyperspectral-Depth Imaging with Learned Diffractive Optics  // Proc. IEEE International Conference on Computer Vision (ICCV). October 2021.
  * Beck, Teboulle [2009] Beck Amir, Teboulle Marc.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems  // SIAM Journal on Imaging Sciences. 2009. 2, 1. 183–202.
  * Bertrand et al. [2020] Bertrand Q., Klopfenstein Q., Blondel M., Vaiter S., Gramfort A., Salmon J. Implicit differentiation of Lasso-type models for hyperparameter optimization  // Proceedings of the 37th International Conference on Machine Learning. 2020.
  * Blanchard [2013] Blanchard Jeffrey D.  Toward deterministic compressed sensing  // Proceedings of the National Academy of Sciences. 110, 4. 2013. 1146–1147.
  * Blondel et al. [2021] Blondel Mathieu, Berthet Quentin, Cuturi Marco, Frostig Roy, Hoyer Stephan, Llinares-López Felipe, Pedregosa Fabian, Vert Jean-Philippe.  Efficient and modular implicit differentiation  // arXiv preprint arXiv:2105.15183. 2021.
  * Bomzon et al. [2001] Bomzon Ze’ev, Kleiner Vladimir, Hasman Erez.  Pancharatnam–Berry phase in space-variant polarization-state manipulations with subwavelength gratings  // Optics letters. 2001. 26, 18. 1424–1426.
  * Boominathan et al. [2022] Boominathan Vivek, Robinson Jacob T., Waller Laura, Veeraraghavan Ashok.  Recent advances in lensless imaging  // Optica. Jan 2022. 9, 1. 1–16.
  * Born, Wolf [2013] Born Max, Wolf Emil.  Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of light. 2013.
  * Candes, Wakin [2008] Candes Emmanuel J., Wakin Michael B. An Introduction To Compressive Sampling  // IEEE Signal Processing Magazine. 2008. 25, 2. 21–30.
  * Chang, Wetzstein [2019] Chang Julie, Wetzstein Gordon.  Deep optics for monocular depth estimation and 3d object detection  // Proceedings of the IEEE/CVF International Conference on Computer Vision. 2019. 10193–10202.
  * Chen et al. [2016] Chen Hou-Tong, Taylor Antoinette J, Yu Nanfang.  A review of metasurfaces: physics and applications  // Reports on Progress in Physics. 2016. 79, 7. 076401.
  * Chen et al. [2018] Chen Wei Ting, Zhu Alexander Y, Sanjeev Vyshakh, Khorasaninejad Mohammadreza, Shi Zhujun, Lee Eric, Capasso Federico.  A broadband achromatic metalens for focusing and imaging in the visible  // Nature Nanotechnology. 2018. 13, 3. 220–226.
  * Del Hougne et al. [2020] Del Hougne Philipp, Imani Mohammadreza F, Diebold Aaron V, Horstmeyer Roarke, Smith David R.  Learned integrated sensing pipeline: reconfigurable metasurface transceivers as trainable physical layer in an artificial neural network  // Advanced Science. 2020. 7, 3. 1901913.
  * Devlin et al. [2016] Devlin Robert C, Khorasaninejad Mohammadreza, Chen Wei Ting, Oh Jaewon, Capasso Federico.  Broadband high-efficiency dielectric metasurfaces for the visible spectrum  // Proceedings of the National Academy of Sciences. 2016. 113, 38. 10473–10478.
  * Duarte et al. [2008] Duarte Marco F., Davenport Mark A., Takhar Dharmpal, Laska Jason N., Sun Ting, Kelly Kevin F., Baraniuk Richard G. Single-pixel imaging via compressive sampling  // IEEE Signal Processing Magazine. 2008. 25, 2. 83–91.
  * Engelberg, Levy [2020] Engelberg Jacob, Levy Uriel.  The advantages of metalenses over diffractive lenses  // Nature Communications. Apr 2020. 11, 1. 1991.
  * Gilbert, Indyk [2010] Gilbert Anna, Indyk Piotr.  Sparse Recovery Using Sparse Random Matrices  // Proceedings of the IEEE. 98, 6. 04 2010. 937–947.
  * Goodman [2005] Goodman Joseph W.  Introduction to Fourier optics, 3rd ed.  // CO: Roberts & Co. Publishers. 2005. 1.
  * Gregor, LeCun [2010] Gregor Karol, LeCun Yann.  Learning Fast Approximations of Sparse Coding  // Proceedings of the 27th on International Conference on Machine Learning. Madison, WI, USA: Omnipress, 2010. 399–406.  (ICML'10).
  * He et al. [2019] He Qiong, Sun Shulin, Zhou Lei.  Tunable/Reconfigurable Metasurfaces: Physics and Applications  // Research. 2019. 2019.
  * Innes et al. [2019] Innes Mike, Edelman Alan, Fischer Keno, Rackauckas Chris, Saba Elliot, Shah Viral B, Tebbutt Will.  A Differentiable Programming System to Bridge Machine Learning and Scientific Computing. 2019.
  * Kamilov [2017] Kamilov Ulugbek S. A Parallel Proximal Algorithm for Anisotropic Total Variation Minimization  // IEEE Transactions on Image Processing. 2017\. 26, 2. 539–548.
  * Khorasaninejad et al. [2016a] Khorasaninejad Mohammadreza, Chen Wei Ting, Devlin Robert C, Oh Jaewon, Zhu Alexander Y, Capasso Federico.  Metalenses at visible wavelengths: Diffraction-limited focusing and subwavelength resolution imaging  // Science. 2016a. 352, 6290. 1190–1194.
  * Khorasaninejad et al. [2016b] Khorasaninejad Mohammadreza, Zhu Alexander Yutong, Roques-Carmes Charles, Chen Wei Ting, Oh Jaewon, Mishra Ishan, Devlin Robert C., Capasso Federico.  Polarization-insensitive metalenses at visible wavelengths  // Nano letters. 2016b. 16, 11. 7229–7234.
  * Kingma, Ba [2014] Kingma Diederik P, Ba Jimmy.  Adam: A method for stochastic optimization  // arXiv preprint arXiv:1412.6980. 2014.
  * Lalanne et al. [1998] Lalanne Philippe, Astilean Simion, Chavel Pierre, Cambril Edmond, Launois Huguette.  Blazed binary subwavelength gratings with efficiencies larger than those of conventional échelette gratings  // Optics letters. 1998. 23, 14. 1081–1083.
  * Langville, Meyer [2011] Langville Amy N., Meyer Carl D. Google's PageRank and Beyond. Princeton: Princeton University Press, 2011\.
  * Li et al. [In Press.] Li William F, Arya Gaurav, Roques-Carmes Charles, Lin Zin, Johnson Steven G, Soljačić Marin.  Transcending shift-invariance in the paraxial regime via end-to-end inverse design of freeform nanophotonics  // Optics Express. In Press.
  * Li et al. [2021] Li Zhaoyi, Pestourie Raphaël, Park Joon-Suh, Huang Yao-Wei, Johnson Steven G., Capasso Federico.  Inverse design enables large-scale high-performance meta-optics reshaping virtual reality  // arXiv preprint arXiv:2104.09702. 2021.
  * Lin et al. [2018] Lin Zin, Groever Benedikt, Capasso Federico, Rodriguez Alejandro W, Lončar Marko.  Topology-optimized multilayered metaoptics  // Physical Review Applied. 2018. 9, 4. 044030.
  * Lin et al. [2021a] Lin Zin, Pestourie Raphaël, Roques-Carmes Charles, Li Zhaoyi, Capasso Federico, Soljačić Marin, Johnson Steven G.  End-to-end metasurface inverse design for single-shot multi-channel imaging  // arXiv preprint arXiv:2111.01071. 2021a.
  * Lin et al. [2021b] Lin Zin, Roques-Carmes Charles, Pestourie Raphaël, Soljačić Marin, Majumdar Arka, Johnson Steven G.  End-to-end nanophotonic inverse design for imaging and polarimetry  // Nanophotonics. 2021b. 10, 3. 1177–1187.
  * Liu, Fan [2012] Liu Victor, Fan Shanhui.  S4: A free electromagnetic solver for layered periodic structures  // Computer Physics Communications. 2012. 183, 10. 2233–2244.
  * Mait et al. [2020] Mait Joseph N., Athale Ravindra A., Gracht Joseph van der, Euliss Gary W. Potential Applications of Metamaterials to Computational Imaging  // Frontiers in Optics / Laser Science. 2020. FTu8B.1.
  * Marcos et al. [2016] Marcos Diego, Lasser Theo, López Antonio, Bourquard Aurélien.  Compressed imaging by sparse random convolution  // Opt. Express. Jan 2016. 24, 2. 1269–1290.
  * Markley et al. [2021] Markley Eric, Liu Fanglin Linda, Kellman Michael, Antipa Nick, Waller Laura.  Physics-Based Learned Diffuser for Single-shot 3D Imaging  // NeurIPS 2021 Workshop on Deep Learning and Inverse Problems. 2021.
  * Molesky et al. [2018] Molesky Sean, Lin Zin, Piggott Alexander Y, Jin Weiliang, Vucković Jelena, Rodriguez Alejandro W.  Inverse design in nanophotonics  // Nature Photonics. 2018. 12, 11. 659–670.
  * Mueller et al. [2017] Mueller JP Balthasar, Rubin Noah A, Devlin Robert C, Groever Benedikt, Capasso Federico.  Metasurface polarization optics: independent phase control of arbitrary orthogonal states of polarization  // Physical Review Letters. 2017. 118, 11. 113901.
  * Pestourie et al. [2018] Pestourie Raphaël, Pérez-Arancibia Carlos, Lin Zin, Shin Wonseok, Capasso Federico, Johnson Steven G.  Inverse design of large-area metasurfaces  // Opt. Express. Dec 2018. 26, 26. 33732–33747.
  * Roques-Carmes et al. [2022] Roques-Carmes Charles, Lin Zin, Christiansen Rasmus E, Salamin Yannick, Kooi Steven E, Joannopoulos John D, Johnson Steven G, Soljačić Marin.  Toward 3D-Printed Inverse-Designed Metaoptics  // ACS Photonics. January 2022. 9, 1. 43–51.
  * Rubin et al. [2019] Rubin Noah A, D’Aversa Gabriele, Chevalier Paul, Shi Zhujun, Chen Wei Ting, Capasso Federico.  Matrix Fourier optics enables a compact full-Stokes polarization camera  // Science. 2019. 365, 6448.
  * Rudin et al. [1992] Rudin Leonid I., Osher Stanley, Fatemi Emad.  Nonlinear total variation based noise removal algorithms  // Physica D: Nonlinear Phenomena. 1992. 60, 1. 259–268.
  * Satat et al. [2017] Satat Guy, Tancik Matthew, Raskar Ramesh.  Lensless Imaging With Compressive Ultrafast Sensing  // IEEE Transactions on Computational Imaging. 2017. 3, 3\. 398–407.
  * Shalaginov et al. [2021] Shalaginov Mikhail, An Sensong, Zhang Yifei, Yang Fan, Su Peter, Liberman Vladimir, Chou Jeffrey, Roberts Christopher, Kang Myungkoo, Rios Ocampo Carlos, Du Qingyang, Fowler Clayton, Agarwal Anuradha, Richardson Kathleen, Baleine Clara, Zhang Hualiang, Hu Juejun, Gu Tian.  Reconfigurable all-dielectric metalens with diffraction-limited performance  // Nature Communications. 2021. 12, 1. 1–8.
  * Shewchuk [1994] Shewchuk Jonathan R.  An Introduction to the Conjugate Gradient Method Without the Agonizing Pain. USA: Carnegie Mellon University, 1994.
  * Simpson, Cao [2020] Simpson Robert E, Cao Tun.  Phase Change Material Photonics  // arXiv preprint arXiv:2006.14441. 2020.
  * Sitzmann et al. [2018] Sitzmann Vincent, Diamond Steven, Peng Yifan, Dun Xiong, Boyd Stephen, Heidrich Wolfgang, Heide Felix, Wetzstein Gordon.  End-to-end optimization of optics and image processing for achromatic extended depth of field and super-resolution imaging  // ACM Transactions on Graphics (TOG). 2018. 37, 4. 1–13.
  * Smith [1997] Smith Steven W.  The scientist and engineer's guide to digital signal processing. 1997\.
  * Strang [2019] Strang Gilbert.  Computational Science and Engineering. 2019.
  * Sun et al. [2020] Sun Qilin, Zhang Jian, Dun Xiong, Ghanem Bernard, Peng Yifan, Heidrich Wolfgang.  End-to-end learned, optically coded super-resolution spad camera  // ACM Transactions on Graphics (TOG). 2020. 39, 2. 1–14.
  * Tibshirani [2013] Tibshirani Ryan J. The lasso problem and uniqueness  // Electronic Journal of Statistics. 2013. 7, none. 1456 – 1490.
  * Tibshirani, Taylor [2011] Tibshirani Ryan J., Taylor Jonathan.  The solution path of the generalized lasso  // The Annals of Statistics. 2011. 39, 3. 1335 – 1371.
  * Tibshirani, Taylor [2012] Tibshirani Ryan J, Taylor Jonathan.  Degrees of freedom in lasso problems  // The Annals of Statistics. 2012. 40, 2. 1198 – 1232.
  * Tropp et al. [2006] Tropp J.A., Wakin M.B., Duarte M.F., Baron D., Baraniuk R.G. Random Filters for Compressive Sampling and Reconstruction  // 2006 IEEE International Conference on Acoustics Speech and Signal Processing Proceedings. 3. 2006. III–III.
  * Tseng et al. [2021a] Tseng Ethan, Colburn Shane, Whitehead James, Huang Luocheng, Baek Seung-Hwan, Majumdar Arka, Heide Felix.  Neural nano-optics for high-quality thin lens imaging  // Nature Communications. 2021a. 12, 1. 6493\.
  * Tseng et al. [2021b] Tseng Ethan, Mosleh Ali, Mannan Fahim, St-Arnaud Karl, Sharma Avinash, Peng Yifan, Braun Alexander, Nowrouzezahrai Derek, Lalonde Jean-Francois, Heide Felix.  Differentiable compound optics and processing pipeline optimization for end-to-end camera design  // ACM Transactions on Graphics (TOG). 2021b. 40, 2. 1–19.
  * Willett et al. [2011] Willett Rebecca M, Marcia Roummel F, Nichols Jonathan M.  Compressed sensing for practical optical imaging systems: a tutorial  // Optical Engineering. 2011. 50, 7. 072601.
  * Yanny et al. [2020] Yanny Kyrollos, Antipa Nick, Liberti William, Dehaeck Sam, Monakhova Kristina, Liu Fanglin Linda, Shen Konlin, Ng Ren, Waller Laura.  Miniscope3D: optimized single-shot miniature 3D fluorescence microscopy  // Light: Science & Applications. Oct 2020. 9, 1. 171.
  * Zhang, Ghanem [2018] Zhang Jian, Ghanem Bernard.  ISTA-Net: Interpretable Optimization-Inspired Deep Network for Image Compressive Sensing  // 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition. 2018. 1828–1837.
  * Zhang et al. [2019] Zhang Yifei, Chou Jeffrey B, Li Junying, others .  Broadband transparent optical phase change materials for high-performance nonvolatile photonics  // Nature Communications. 2019. 10, 1. 1–9.
  * Zou, Hastie [2005] Zou Hui, Hastie Trevor.  Regularization and Variable Selection via the Elastic Net  // Journal of the Royal Statistical Society. Series B (Statistical Methodology). 2005. 67, 2. 301–320.


11institutetext:  LERMA & UMR8112 du CNRS, Observatoire de Paris, PSL
University, Sorbonne Universités, CNRS, F-75014 Paris, France
11email<EMAIL_ADDRESS>22institutetext: Green Bank Observatory, P.O.
Box 2, Green Bank, WV 24944, USA 33institutetext: Lycée Jean-Baptiste Corot,
F-91605 Savigny-sur-Orge, France 44institutetext: Institut de Radioastronomie
Millimétrique (IRAM), 300 rue de la Piscine, 38400 Saint-Martin d’Heères,
France

# Radio telescope total power mode: improving observation efficiency††thanks:
Based on observations carried out with the Green Bank Telescope and the IRAM
30-m. The Green Bank Observatory is a facility of the National Science
Foundation operated under cooperative agreement by Associated Universities,
Inc. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain).

Laurent Pagani 11 David Frayer 22 Bruno Pagani 33 Charlène Lefèvre 44

(Received July 20, 2020; accepted September 19, 2020)

###### Abstract

Aims. Radio observing efficiency can be improved by calibrating and reducing
the observations in total power mode rather than in frequency, beam, or
position-switching modes.

Methods. We selected a sample of spectra obtained from the Institut de Radio-
Astronomie Millimétrique (IRAM) 30-m telescope and the Green Bank Telescope
(GBT) to test the feasibility of the method. Given that modern front-end
amplifiers for the GBT and direct Local Oscillator injection for the 30 m
telescope provide smooth pass bands that are a few tens of megahertz in width,
the spectra from standard observations can be cleaned (baseline removal)
separately and then co-added directly when the lines are narrow enough (a few
km s-1), instead of performing the traditional ON minus OFF data reduction.
This technique works for frequency-switched observations as well as for
position- and beam-switched observations when the ON and OFF data are saved
separately.

Results.  The method works best when the lines are narrow enough and not too
numerous so that a secure baseline removal can be achieved. A signal-to-noise
ratio improvement of a factor of $\sqrt{2}$ is found in most cases, consistent
with theoretical expectations.

Conclusions. By keeping the traditional observing mode, the fallback solution
of the standard reduction technique is still available in cases of suboptimal
baseline behavior, sky instability, or wide lines, and to confirm the line
intensities. These techniques of total-power-mode reduction can be applied to
any radio telescope with stable baselines as long as they record and deliver
the ONs and OFFs separately, as is the case for the GBT.

###### Key Words.:

Telescopes – Techniques: spectroscopic – Methods: observational – Radio lines:
general

††offprints: L.Pagani

##  1 Introduction

For a long time, heterodyne receivers in the millimeter (mm) domain were made
of front-end mixers with optical injection of the Local Oscillator followed by
intermediate-frequency amplifiers and production of spectra via filter banks
for which the homogeneity of the gain was not secured from channel to channel.
Combined to on-axis telescopes, these suffer from standing waves between the
secondary mirror and the receiver. Local Oscillator injections via Martin-
Pupplet optical diplexer devices with relatively narrow bandpass were adding
to the frequency-dependent gain variations of the whole system (a few
telescope descriptions can be found in, e.g., Castets et al. 1988, Booth et
al. 1989, Schuster et al. 2004). Combining the effects from the frequency
response of the mixer, the noisy receivers, the inhomogeneous backend, and the
intermediate-frequency amplifiers with their imperfect impedance match to the
mixer output resulted in complex bandpass profiles with high amplitudes,
making the weak astronomical lines difficult to detect without subtracting a
reference spectrum with a similar profile to the observations. This
subtraction has to be performed relatively quickly to protect against gain
variations and atmospheric absorption fluctuations. Various strategies are
employed to get a reference spectrum that is as close as possible to the
scientific spectrum. These include the well-known position-switching (the
telescope primary mirror is shifted to an emission-free region in the sky),
beam-switching (the beam is deviated from the telescope pointing direction by
a mirror or simply by wobbling the secondary mirror within a few arcminutes
away from the pointed direction), and frequency-switching observing modes
(where there is no mechanical movement, but the Local Oscillator reference
frequency is slightly changed so that the subtraction does not cancel the line
but only the baseline offset).

For the Rayleigh-Jeans temperature scale, the noise fluctuation (root mean
square deviation, $\sigma$ or rms hereafter) of a radio spectrum is given by
the radiometer equation:

$\sigma(T)=\frac{\eta T_{\mathrm{sys}}}{\sqrt{\delta\nu\times\tau}},$ (1)

where $T_{\mathrm{sys}}$ is the total system noise temperature111For the GBT,
Tsys is measured in the Ta scale, while for the IRAM 30-m telescope, it is
measured in the T${}_{\mathrm{a}}^{*}$ scale, which includes atmospheric
attenuation and forward beam efficiency corrections: T${}_{\mathrm{a}}^{*}$ =
T${}_{\mathrm{a}}\times\frac{e{{}^{tau}}}{\eta_{\mathrm{ffs}}}$ (Kutner &
Ulich 1981), including the receiver noise and all noise from the sky and
ground spillover, $\delta\nu$ is the spectral resolution, and $\tau$ the
integration time. Here, $\eta$ is sometimes introduced to take into account
other losses such as one- or two-bit autocorrelator conversion losses. When
using beam-switching or position-switching modes, a second spectrum of similar
characteristics is obtained to be subtracted from the first one. Because the
noise fluctuations of both spectra are uncorrelated, the subtraction increases
the noise fluctuation temperature by $\sqrt{2}$. If we consider that $\tau$
represents the total (ON+OFF) time, then the noise fluctuation temperature
increases by another $\sqrt{2}$, as only half of the time was spent on each
individual spectrum, effectively doubling the final noise fluctuation of the
spectrum (and more time is lost in overheads, like the mechanical displacement
of mirrors, or the rotation of the whole telescope). If the frequency-switch
mode has been used instead, then the OFF spectrum becomes a second ON spectrum
and the penalty is only $\sqrt{2}$ but at the expense of spoiling the baseline
because of the frequency dependence of the receiver (and the standing wave
patterns) that prevents an exact superposition of the baseline structure from
the ON1 and ON2 spectra. Improved strategies consist of integrating over
longer periods of time for the OFF spectrum and subtracting the same OFF
spectrum from several ONs, minimizing the time spent OFF source and the noise
fluctuations added in the subtraction. These strategies require very stable
receivers and sky conditions to be efficient, but can marginally beat the
performance of the frequency-switch mode in terms of noise fluctuations and
provide much flatter baselines; however, they become less efficient in the end
if spatial smoothing is applied, because the OFF is identical for the adjacent
ON-OFF pairs such that their average does not diminish the noise fluctuations
efficiently, and then frequency-switching becomes preferable.

With the advent of high-frequency amplifiers, drop of LO optical diplexer
injection, Fourier-transform spectrometers, and the benefit of an off-axis
telescope such as the Green Bank Observatory 100-m Telescope (GBT), the
quality and stability of the bandpass has considerably improved, meaning that
total power (staring in a fixed direction at fixed frequency) mode
observations can now be reconsidered. This would avoid the quadratic addition
of noise fluctuations and therefore save the $\sqrt{2}$ factor discussed
above. However, there are caveats in operating in pure total power mode, and
we propose that observations be run in the usual way, but with the data
reduction performed on individual ON and OFF spectra to exploit the benefit of
total power mode when possible while not losing the fallback possibility of
standard data reduction. We present the method in Sect. 2 and discuss a few
cases in Sect. 3.

##  2 Method

The calibration of mm-wave radiotelescopes has been discussed in Kutner &
Ulich (1981) and a technical report from the Institut de Radio-Astronomie
Millimétrique (IRAM) 30-m presents the details of the calibration
procedure222http://www.iram.es/IRAMES/mainWiki/CalibrationPapers?action=AttachFile&do=view&target=kramer_1997_cali_rep.pdf
for that telescope which is representative of mm-wave radiotelescope
calibrations. Nevertheless, there are differences between the GBT and the IRAM
30-m telescope. For example, the GBT has amplifiers before the heterodyne
mixer and has no image sideband.

Briefly, the observations of two internal loads at different temperatures give
a linear fit between backend voltages or detector counts and temperatures (in
the Rayleigh-Jeans approximation333The physical temperature of the loads is
usually considered to be equal to their Rayleigh-Jeans radiation temperature.
This is only true for $h\nu<<kT$, with $h$ and $k$ being the Heisenberg and
Boltzmann constants, $\nu$ the frequency, and $T$ the physical temperature.
See footnote 2, and Sect. 1.3 of that paper for further considerations on the
validity of the approximation.). The slope of the fit is referred to as the
gain and is defined by

$g=\frac{T_{\mathrm{hot}}-T_{\mathrm{cold}}}{V_{\mathrm{hot}}-V_{\mathrm{cold}}},$
(2)

where $T_{\mathrm{hot}}$ (resp. cold) represents the temperature of the hot
(resp. cold) load and $V_{\mathrm{hot}}$ (resp. cold) represents the voltage
(or counts) of the hot (resp. cold) load measured at the receiver backend. If
$T_{\mathrm{rec}}$ is the receiver equivalent noise temperature (its value is
derived from the gain measurement), any signal from the sky delivers a voltage
$V_{\mathrm{sky}}$ which can be converted to a Rayleigh-Jeans temperature with

$T_{\mathrm{sky}}=g\times V_{\mathrm{sky}}-T_{\mathrm{rec}}.$ (3)

The sky signal is a composite of ground-emission spillover, atmospheric
emission, and cosmic signal attenuated by the atmospheric absorption and is
beyond the scope of this discussion (see Kutner & Ulich 1981 for more
information). If two measurements of the sky are performed along one of the
switching procedures described in Sect. 1, we can compute their difference to
cancel all unwanted signals and retrieve only the cosmic signal of interest
($T_{\mathrm{sky}}$ is hereafter referred to with the more commonly used
$T_{\mathrm{a}}$ or antenna temperature, and similarly $V_{\mathrm{sky}}$
becomes $V_{\mathrm{a}}$):

$\Delta
T_{\mathrm{a}}=T_{\mathrm{a,ON}}-T_{\mathrm{a,OFF}}=g\times(V_{\mathrm{a,ON}}-V_{\mathrm{a,OFF}}),$
(4)

where $\Delta T_{\mathrm{a}}$ represents the cosmic signal we want to observe
(yet uncorrected for various losses). If we define

$T_{\mathrm{sys}}=T_{\mathrm{a}}+T_{\mathrm{rec}},$ (5)

we get

$T_{\mathrm{sys}}=g\times V_{\mathrm{a}}$ (6)

and eq. 4 becomes the more familiar form

$\Delta
T_{\mathrm{a}}=T_{\mathrm{sys}}\times\frac{(V_{\mathrm{a,ON}}-V_{\mathrm{a,OFF}})}{V_{\mathrm{a,OFF}}}.$
(7)

Here, $\Delta T_{\mathrm{a}}$ suffers from various sources of noise which we
can analyse:

$\begin{aligned} \sigma^{2}(\Delta
T_{\mathrm{a}})=\sigma^{2}(g)\times(V_{\mathrm{a,ON}}-V_{\mathrm{a,OFF}})^{2}\\\
+g^{2}\times(\sigma^{2}(V_{\mathrm{a,ON}})+\sigma^{2}(V_{\mathrm{a,OFF}}))\end{aligned},$
(8)

where $\sigma$ is the root mean square error expressed in a similar form as in
eq. 1, that is, depending on the inverse square root of the integration time
and frequency sampling of the radiometer. Where the ON and OFF signals do not
cancel (i.e. inside the observed line), the equation can be expressed in the
familiar form,

$\frac{\sigma^{2}(\Delta T_{\mathrm{a}})}{\Delta
T_{\mathrm{a}}^{2}}=\frac{\sigma^{2}(g)}{g^{2}}+\frac{\sigma^{2}(V_{\mathrm{a,ON}})+\sigma^{2}(V_{\mathrm{a,OFF}})}{(V_{\mathrm{a,ON}}-V_{\mathrm{a,OFF}})^{2}}.$
(9)

In the rest of the spectrum (which is usually named the baseline), the noise
is simply444for a long time, mm heterodyne receivers were much noisier than
the strength of the line and the total noise inside the line was dominated by
the receiver noise. This is no longer the case for the strongest lines with
modern receivers and this should be taken into account; though generally, the
total noise becomes negligible in such cases.

$\sigma^{2}(\Delta
T_{\mathrm{a}})=g^{2}\times(\sigma^{2}(V_{\mathrm{a,ON}})+\sigma^{2}(V_{\mathrm{a,OFF}})).$
(10)

If we suppose that the hot and cold load temperatures are known with enough
accuracy to have a negligible contribution to the noise, the gain noise can be
expressed as

$\frac{\sigma^{2}(g)}{g^{2}}=({T_{\mathrm{hot}}-T_{\mathrm{cold}}})\frac{\sigma^{2}(V_{\mathrm{hot}})+\sigma^{2}(V_{\mathrm{cold}})}{(V_{\mathrm{hot}}-V_{\mathrm{cold}})^{2}}.$
(11)

If instead of observing with a switching mode (SM) one observes in total power
mode (TP), the signal is not extracted from the receiver plus sky noise during
the observations; its expression is simply eq. 3 and the error budget is
(supposing $T_{\mathrm{rec}}$ is known with enough precision to neglect its
contribution)

$\frac{\sigma^{2}(T_{\mathrm{a}})}{T_{\mathrm{a}}^{2}}=\frac{\sigma^{2}(g)}{g^{2}}+\frac{\sigma^{2}(V_{\mathrm{a}})}{V_{\mathrm{a}}^{2}}.$
(12)

Equations 9 and 12 can be compared only if we have manually subtracted the sky
plus receiver noise contribution from the TP spectrum (this will be discussed
in the following section), in which case $T_{\mathrm{a}}(\mathrm{TP})=\Delta
T_{\mathrm{a}}(\mathrm{SM})$. and
$V_{\mathrm{a}}(\mathrm{TP})=(V_{\mathrm{a,ON}}-V_{\mathrm{a,OFF}})(\mathrm{SM})$.
When the receiver and the sky transparency are constant inside the bandpass of
an observation, the gain and receiver noise temperature can be calculated by
integrating over all the channels of the spectrometer, largely decreasing —
suppressing in practice— their contribution to the final noise. The noise then
simply depends on $\sigma(V_{\mathrm{a}})$, which is inversely proportional to
the square root of the integration time, everything else being equal. Using
the same total integration time $\tau$ for the TP observations and for the sum
of the two phases of the SM observations, we get

$\frac{\sigma^{2}(\Delta
T_{\mathrm{a}})}{\sigma^{2}(T_{\mathrm{a}})}=\frac{\sigma^{2}(V_{\mathrm{a,ON}})+\sigma^{2}(V_{\mathrm{a,OFF}})}{\sigma^{2}(V_{\mathrm{a}})}=\frac{{\frac{2}{\tau}}+{\frac{2}{\tau}}}{{\frac{1}{\tau}}}=4.$
(13)

The SM observations are therefore twice as noisy as the TP mode ones. This is
true for the mechanical (position or beam) SMs because only one phase is
exposed to the signal, but in the case of the frequency SM, the subsequent
folding of the spectrum reduces the noise further by $\sqrt{2}$ by averaging
two independent realizations of the measurement and the TP mode advantage is
only $\sqrt{2}$.

It is not always possible to use a constant gain and $T_{\mathrm{rec}}$
throughout the spectrum and the balance between the gain noise and the
observation noise contributions must be addressed when the calibration is
performed channel-wise555At the IRAM 30-m the calibration scheme uses an
intermediate configuration where the variations of the receiver temperature
and sky opacity across the band are averaged by chunks of 20 MHz with their
new calibration program MRTCAL. https://www.iram-
institute.org/medias/uploads/mrtcal-check.pdf. We have seen (Eqs. 9 and 11)
that both sources of noise depend primarily on the voltage or count
measurement noise which depends upon integration time (eq. 1). Since the
calibrators are usually observed on a short timescale (1 – 5 seconds
typically), while the sky observations can last 60 seconds or more666we do not
discuss on-the-fly observing mode here, the rms of the calibration phases is
the highest. In TP observations, the denominator is high
($V_{\mathrm{a}}=V_{\mathrm{rec}}$+$V_{\mathrm{atmosphere}}$+$V_{\mathrm{source}}$)
and comparable to the load measurements
($V_{\mathrm{load}}=V_{\mathrm{rec}}+V_{\mathrm{cold}}$ or
$V_{\mathrm{hot}}$). Therefore, the dominant error term comes from the
calibration part itself. In SM observations, the observation denominator
($V_{\mathrm{a,ON}}$-$V_{\mathrm{a,OFF}}$) is much smaller than
$V_{\mathrm{load}}$ and despite a lower rms, this term dominates the gain and
calibration errors, which is preferable.

Figure 1: 12CO (_J_ :1 – 0) frequency SM observation obtained at the IRAM 30-m
telescope with the EMIR receiver and FTS backend. Top: Standard folded
spectrum after subtraction and calibration. Middle: Original two phases
calibrated independently. Bottom: Two phases directly averaged after
realignment without using noise-averaged gain or subtracting a baseline from
the raw data. The baseline shows no ripples but the noise is three times
higher than in the upper panel. The vertical axis has the same amplitude for
all three boxes (16 K). The color of the rms figure corresponds to the color
of the spectra to which it pertains.

Figure 1 shows an example of the problem encountered when separately
calibrating the two phases on a channel-wise mode. The final spectrum noise is
dominated by the calibration noise and is therefore noisier, which confers no
advantage. To benefit from the $\sqrt{2}$ improvement, two possibilities are
available: (1) use a noise-averaged gain (either constant or smoothed over an
intermediate width, as is the case for the IRAM 30-m telescope with MRTCAL),
or (2) subtract a baseline from the raw data (backend voltage or counts $V$)
before applying the calibration to minimize the contribution of the gain
(calibration) noise $\sigma(g)$.

##  3 Examples

###  3.1 The IRAM 30-m telescope

Though the IRAM 30-m telescope is still equipped with front-end mixers and the
antenna is on-axis, the replacement of the optical LO injection by line
injection in the new EMIR receivers (for Eight MIxer Receivers, Carter et al.
2012) has improved the baseline performance of the observations, and even
frequency-switched observations are only subject to relatively limited
baseline ripples. We therefore study a few cases with two different backends:
the fast Fourier Transform Spectrometer (FTS) in its narrow mode (8 $\times$
37275 channels, 48.83 kHz each) and the VErsatile SPectrometer Array (VESPA),
which is an autocorrelator used in a narrow-window, high-resolution mode (20
MHz, 10 kHz respectively). Presently, IRAM does not deliver the individual
phases (ON/OFF or Freq 1/Freq 2) of the observations but only the final
calibrated spectrum (ON - OFF or Freq 1 - Freq 2 with subsequent folding). We
modified the calibration code (Millimeter RadioTelescope CALibration, MRTCAL)
to obtain access to the individual raw phases for this work777This version is
incomplete and archaic and not meant to be distributed..

####  3.1.1 Fourier Transform Spectrometer data

Figure 2: Same observations as in Fig. 1. Top: Sky observation with the full
band IRAM 30-m FTS backend (one setup, one polarization). Bottom: Resulting
folded frequency SM observation after one-minute integration. The strong
positive peak with two negative half-intensity ones is the 12CO line displayed
in the upper panel of Fig. 1.

Figure 2 clearly shows that the subtraction of a baseline from the single
phase observation on the full bandwidth (upper panel) will never manage to
compete with the differential observation done in frequency SM (lower panel).
There is therefore no hope to use this total power technique with wide lines
and one must focus on a part of the backend small enough to be able to recover
a flat baseline around the line of interest.

Figure 3: Same observations as in Fig. 1. Top: Raw autocorrelator counts of
the two phases re-aligned on top of each other and averaged together. A
polynomial of 25th order is fitted (in red). Red vertical lines delimit the
weak telluric and strong astrophysics 12CO lines. Middle: Baseline-subtracted
raw average. Bottom: Comparison between the standard data reduction of the
observations (black) and TP mode reduction of the data (red). Table 1:
Parameter optimization for TP data reduction.

window width | polynomial | rms | line
---|---|---|---
MHz | channels | km s-1 | degree | ratioaa$a$Values close to the expected $\sqrt{2}$ improvement are boldfaced. | ratiobb$b$line ratios within 1 $\sigma$ (in terms of the integrated intensity, i.e. 0.12 Kkm s-1, which is 0.986$<$R$<$1.014) are boldfaced.
53.53 | 1096 | 139 | 40 | 0.871 | 0.998
– | – | – | 60 | 1.240 | 0.979
– | – | – | 65 | 1.275 | 1.012
– | – | – | 70 | 1.288 | 1.020
– | – | – | 75 | 1.341 | 0.998
– | – | – | 100 | 1.419 | 1.020
49.53 | 1014 | 129 | 55 | 1.242 | 0.980
– | – | – | 60 | 1.274 | 1.018
– | – | – | 65 | 1.294 | 1.008
– | – | – | 70 | 1.341 | 0.988
– | – | – | 90 | 1.409 | 0.995
– | – | – | 95 | 1.411 | 0.983
45.53 | 932 | 118 | 35 | 0.929 | 0.983
– | – | – | 55 | 1.252 | 1.010
– | – | – | 60 | 1.267 | 1.017
– | – | – | 65 | 1.329 | 1.000
– | – | – | 80 | 1.389 | 0.994
– | – | – | 95 | 1.424 | 0.988
41.53 | 851 | 108 | 30 | 0.829 | 1.020
– | – | – | 50 | 1.251 | 1.017
– | – | – | 55 | 1.277 | 1.006
– | – | – | 85 | 1.408 | 1.004
37.53 | 769 | 98 | 45 | 1.236 | 1.006
– | – | – | 50 | 1.256 | 1.016
– | – | – | 55 | 1.332 | 0.992
– | – | – | 75 | 1.410 | 1.016
33.53 | 687 | 87 | 25 | 0.823 | 1.021
– | – | – | 45 | 1.259 | 1.013
– | – | – | 60 | 1.379 | 1.000
29.53 | 605 | 77 | 35 | 1.218 | 1.010
– | – | – | 40 | 1.259 | 1.013
– | – | – | 60 | 1.409 | 1.013
25.53 | 523 | 66 | 35 | 1.316 | 0.996
– | – | – | 45 | 1.412 | 1.019
– | – | – | 50 | 1.433 | 1.020
– | – | – | 60 | 1.492 | 0.985
21.53 | 441 | 56 | 25 | 1.177 | 1.007
– | – | – | 30 | 1.276 | 0.981
– | – | – | 40 | 1.369 | 0.984
– | – | – | 50 | 1.424 | 1.011
– | – | – | 65 | 1.513 | 0.987
13.53 | 277 | 35 | 15 | 1.209 | 0.992
– | – | – | 25 | 1.421 | 1.012
– | – | – | 30 | 1.473 | 1.009
– | – | – | 45 | 1.597 | 1.023
– | – | – | 50 | 1.619 | 1.008
9.53 | 195 | 25 | 10 | 1.287 | 0.994
– | – | – | 20 | 1.536 | 1.011
5.53 | 113 | 14 | 5 | 1.459 | 0.987
– | – | – | 10 | 1.667 | 0.993
– | – | – | 15 | 1.784 | 0.977
888

Figure 3 shows a small section of the bandpass (35 km s-1) centered on the
12CO (_J_ :1 – 0) line. The upper panel shows the average of the two phases
expressed in counts (analog to digital units, ADU), the cosmic and telluric
line emissions are masked away by two windows, and a polynomial baseline of
25th order has been subtracted from it. The flattened spectrum is shown in the
middle window. The gain to convert the spectrum to the antenna temperature
scale is then applied, and the final spectrum is compared to the original
spectrum obtained by folding the two phases together. The spectra are
identical and their noise is in the $\sqrt{2}$ ratio
($\frac{0.214}{0.151}=1.42$). The polynomial degree of the baseline necessary
to achieve this result is a function of the width of the spectrum onto which
the baseline is fitted. Though this might be specific to each set of
observations, we have explored the combination of frequency window size and
baseline polynomial degree to check their interdependence and evaluate the
minimum polynomial degree needed to retrieve a consistent line-integrated
intensity and a baseline flat enough to recover the theoretical $\sqrt{2}$
improvement on the noise (Table 1). It can be noted that for polynomial
degrees that are too high compared to the channel numbers, the noise drops
even lower than expected. This indicates that the high polynomial degree
induces the fit of the noise itself and removes part of it. Consequently, the
TP noise data reduction cannot be lead blindly but must be compared to the
standard SM data reduction to check that both the noise and the line-
integrated intensity achieve their expected values. In the case of no signal,
the same procedure should also apply to get the correct noise improvement.

####  3.1.2 VESPA data

Figure 4: IRAM 30-m telescope observations of the 13CO (_J_ :1 – 0) line with
the VESPA backend, full window. The original bandwidth is $\sim$20 MHz and
drops to $\sim$12 MHz because of the $\pm$3.9 MHz frequency SM. The three
panels display the same type of information as in Fig. 3.

VESPA can be split into many small windows of 10 to 40 MHz with high frequency
resolution, i.e., 3.3 to 40 kHz. There are many other modes in addition to
those discussed here, but these are beyond the scope of this study. On such
narrow windows, the baseline is relatively smooth but still needs polynomial
fitting of high order (18-20) to remove small ripples and retrieve the correct
line intensity (Fig. 4). However, the noise gain is less than the expected
$\sqrt{2}$. We have identified the origin of this discrepancy to the greater-
than $\sqrt{2}$ noise diminution when folding the spectrum. The measured
unfolded spectrum noise in the present case is 0.37 K, and we therefore expect
a folded spectrum noise of 0.26 K instead of the measured 0.23 K (Fig. 4, with
an unfolded noise of 0.33 K, expected noise of 0.23 K, and measured noise of
0.20 K for the other polarization). Though adding or subtracting spectra
should have the same impact on noise, the noise we obtain in the TP mode data
reduction is close to expected (0.37 K $\rightarrow$ 0.185 K, measured 0.181
K). The standard folding method returns a spectrum with a noise lower than
expected because the frequency displacement could be a fractional number of
channels, and when the folding is performed the folded channels are split into
two-subchannels to be added to the two adjacent channels. This introduces a
noise correlation between the channels which can artificially lower the noise.
We note that when half-channel intensities are added, their noise is not
changed, and therefore this is similar to smoothing the spectroscopic
resolution. The diminution is maximal when the channel displacement is exactly
a whole number plus one-half of a channel and can reach $\sqrt{2}$ of
supplemental noise diminution. This is not specific to VESPA data and can
happen to any spectrometer when using the frequency SM999Observers should be
aware of this artifact and should select their frequency throw to adjust to a
whole number of channels though this might prove difficult when several
backends with different resolutions are used in parallel.. However this noise
is correlated and a subsequent smoothing in frequency would not reduce the
noise as much as expected.

###  3.2 The GBT 100-m

The GBT presently uses the VErsatile GBT Astronomical Spectrometer (VEGAS) as
a backend, which is a somewhat similar autocorrelator to VESPA in its agility,
multi-windowing, and variable resolution
capabilities101010https://www.cv.nrao.edu/ aroshi/VEGAS/skyfretobb.pdf.
Contrary to those from the IRAM 30-m telescope, the data are provided raw and
uncalibrated to the user. Any user can therefore treat the two phases of the
observations separately in order to perform TP mode data reduction. Another
autocorrelator, known as the GBT spectrometer, was in use before VEGAS, and
the results are not sensitive to the backend itself for similar bandwidth and
resolution. We visit the K-band and W-band cases here.

####  3.2.1 NH3 (K-band) observations

  * •

The GBT spectrometer

Figure 5: NH3 observations with the GBT 100-m in position SM. The full 50 MHz
sub-window of the GBT spectrometer is represented. Only the ON phase is shown.

Figure 5 shows the complete 50 MHz spectrum on source of a NH3 line for an
astrophysical source. Except at the edges where the bandpass drops rapidly to
zero, the band is relatively smooth and flat as for VESPA, and the lines are
readily visible in the band.

Figure 6: Same observations as in Fig. 5 but dropping the ends, with the same
three types of panels as in Fig. 3

A polynomial fit of 13th degree is sufficient here to remove the baseline
cleanly as seen in Fig. 6 which proves that in such a case, the TP mode can
provide baselines as flat as those obtained naturally with mechanical SMs (for
narrow lines). The noise is improved by the expected $\sqrt{2}$ value, and
this is also the case for the frequency SM observations of NH3 with the same
spectrometer (Fig. 7.)

Figure 7: Comparison between standard reduction (position and frequency SM,
black) with TP mode reduction (red) for NH3 observations with the old GBT
spectrometer in two different sources. The spectrum shown at the top is
identical to the one displayed in Fig. 6 but zoomed in on the lines. Noise
improvement is close to $\sqrt{2}$. The difference between both reduction
modes is displayed in blue and shifted by -0.5K.

  * •

VEGAS

Figure 8: GBT NH3 observations with the VEGAS backend and the KFPA receiver,
central pixel. Only one phase of the frequency SM is shown. The full 50 MHz
sub-window of VEGAS is displayed.

VEGAS observations of the NH3 line are quite similar to those provided by the
GBT spectrometer, as expected, but the baseline is even flatter, due to the
new seven-pixel K-Band Focal Plane Array receiver, (KFPA, Fig. 8) allowing the
use of a polynomial of very low degree (4 or 5) to zero the baseline (Fig. 9).
The baseline is also flat enough so that the channel-wise gain is hardly
variable and can be replaced by a constant value, therefore suppressing the
calibration noise from eq. 12. In this figure, the phases have indeed been
calibrated separately before being flattened. The use of a constant gain to
calibrate the data allows the user to calibrate individual phases without
subtracting any offset, but conversely the shoulders in the frequency SM data
are hidden and their contribution to the noise lowers the average while it
should increase it if attenuation is properly taken into account. If we
measure the noise in this spectrum to the edges of the plot, which remain
inside the original shoulders as displayed in Fig. 8, the standard deviation
is 0.165 K instead of 0.173, with constant gain applied to the whole spectrum,
and 0.179 K if calibrated channel-wise.

Figure 9: GBT NH3 observations with the VEGAS backend showing a comparison
between TP mode and frequency SM. The top window shows the individual spectra
at two different frequencies. The frequency axis is narrowed to mask the
shoulders at both ends (seen in Fig. 8). The middle box displays the two
spectra realigned and averaged (TP mode). Because of the frequency shift, both
shoulders are brought back in the frequency range. A fifth degree polynomial
is fitted in between them (red line). The bottom box shows the comparison
between the standard frequency SM folded spectrum (black) and the TP mode
spectrum (red, shoulders have been erased).

####  3.2.2 W-band observations

Figure 10: N2D+ (_J_ :1 – 0) raw spectrum observed with the GBT spectrometer
in frequency SM. The full 50 MHz bandwidth is displayed. Only one phase is
shown. Figure 11: N2D+ (_J_ :1 – 0) data reduction in TP mode compared to the
original frequency SM reduction. The residual in the central part of the
spectrum is shown in blue with a negative offset of -0.5 K. The ends have been
erased because the negative lines from the frequency SM folding have no
counterpart in the TP mode data reduction. Red vertical lines mark the masks
used for the TP mode data reduction to compute the baseline, while the black
vertical lines mark those used for the frequency SM data reduction. Part of
the negative lines are outside the window.

The W-band observations we present here were taken with the old GBT
spectrometer. Though the baseline is smooth in appearance (Fig. 10), the
baseline subtraction from the ON+OFF spectra realigned on top of each other
leaves small ripples which can be reduced by narrowing the spectral window but
cannot be completely erased. However, the line residual is null on average,
and the noise improvement is still $\sqrt{2}$ for a baseline polynomial degree
of 16th order (Fig. 11).

###  3.3 Other methods

When the line is too wide for this kind of baseline fitting, there are still
various advantages to treating the ON and the OFF observations separately. For
example, several adjacent OFF spectra in space and time (e.g., offsets of a
map) can be averaged together to be subtracted from each individual ON
spectrum, reconstructing a posteriori the ‘one long OFF–several short ONs’
observation method (which is not available on all telescopes). The OFF can
also be smoothed as strongly as possible, as long as narrow features from the
baseline itself are not smoothed out because that would replace Gaussian noise
with systematic error.

Figure 12: Same NH3 data from GBT + VEGAS spectrometer as in Fig. 9. The top
window shows the first frequency-shifted spectrum with the right-hand side
smoothed by a median of 0.4 MHz in width. The middle window shows the second
frequency-shifted spectrum smoothed on the left-hand side. The bottom window
shows a comparison between the original folded spectrum and the half-smoothed
folded spectrum. The baseline is computed between the two vertical dashed
lines as in Fig. 9.

For frequency SM spectra, a similar method can be applied except that the two
sides of the spectra have to be treated separately because both contain
signal. We applied this method to the VEGAS GBT NH3 data with success (Fig.
12)

##  4 Conclusions

Total power data reduction with GBT data can easily be conducted with the
present data-delivery and data-reduction tools (GBTIDL, a GBT customization of
the interactive data language$\texttrademark$) offered at the GBO. The
advantage of the present procedure is a gain of $\sqrt{2}$ in sensitivity at
no cost whilst retaining the ability to check the result by comparing it to
the normal differential data reduction or even to fall back on the standard
mode reduction in case of problems (especially if the line is too large to
support a high-degree polynomial fitting without damage). This method is
particularly suited for narrow lines in cold clouds that are usually weak and
few in number per GHz. In principle, this method could also be used with IRAM
30-m data or any other radiotelescope data provided the baseline is flat
enough and single phase observations are retrievable.

Though not illustrated in this paper because of a lack of suitable data, in
case of wide lines, all OFF observations of different positions can be
averaged together (and smoothed if necessary) to reach similar results in the
same manner as on-the-fly position SM is proceeding, and frequency SM
observations can be half averaged and/or smoothed to reach sensitivities
comparable to those of TP mode observations.

Finally, OTF data can be advantageously treated in this manner: after the
window width and baseline degree have been adjusted to optimize the data
reduction, the treatment can be applied to all spectra automatically. In
particular, in the case of position SM OTF data, the noise contribution of the
OFF spectrum does not decrease when convolving the data to produce the final
map because it is the same OFF for all consecutive ONs. Getting rid of this
noise contribution could have a noticeable impact.

###### Acknowledgements.

We thank the anonymous referee and the editor, Thierry Forveille, for their
remarks which helped to improve the clarity of the paper. We thank Sébastien
Bardeau for his help to modify MRTCAL and explanations on data resampling.

## References

  * Booth et al. (1989) Booth, R. S., Delgado, G., Hagstrom, M., et al. 1989, A&A, 216, 315
  * Carter et al. (2012) Carter, M., Lazareff, B., Maier, D., et al. 2012, A&A, 538, A89
  * Castets et al. (1988) Castets, A., Lucas, R., Lazareff, B., et al. 1988, A&A, 194, 340
  * Kutner & Ulich (1981) Kutner, M. L. & Ulich, B. L. 1981, ApJ, 250, 341
  * Schuster et al. (2004) Schuster, K. F., Boucher, C., Brunswig, W., et al. 2004, A&A, 423, 1171