Simulating Imaging of Large Scale Radio Arrays on the Lunar Surface

Summary

A simulation framework for testing the imaging capabilities of large-scale radio arrays on the lunar surface is presented. Major noise components are discussed, and a software pipeline is walked through with details on how to customize it for novel scientific uses.

Abstract

In recent years there has been a renewed interest in returning to the Moon for reasons both scientific and exploratory in nature. The Moon provides the perfect training ground for building large scale bases that one may apply to other planets like Mars. The existence of a radio quiet zone on the lunar far side has promise for early universe studies and exoplanet searches, while the near side provides a stable base that may be used to observe low frequency emissions from Earth’s magnetosphere that may help gauge its response to incoming space weather. The construction of a large-scale radio array would provide large scientific returns as well as acting as a test of humanity’s ability to build structures on other planets. This work focuses on simulating the response of small to large-scale radio arrays on the Moon consisting of hundreds or thousands of antennas. The response of the array is dependent on the structure of the emission along with the configuration and sensitivity of the array. A set of locations are selected for the simulated radio receivers, using Digital Elevation Models from the Lunar Orbiter Laser Altimeter instrument on Lunar Reconnaissance Orbiter to characterize the elevation of the receiver locations. A custom Common Astronomy Software Applications code is described and used to process the data from the simulated receivers, aligning the lunar and sky coordinate frames using SPICE to ensure the proper projections are used for imaging. This simulation framework is useful for iterating array design for imaging any given scientific target in a small field of view. This framework does not currently support all sky imaging.

Introduction

The field of radio astronomy began in 1932 with the accidental detection of galactic radio emission by Karl G. Jansky¹ at 20 MHz, in a range now commonly called the low frequency radio. Ever since then, radio astronomy has grown rapidly, catching up with higher frequency optical observations that have been going on for centuries longer. Another breakthrough was the utilization of radio interferometry, where groups of antenna separated by large distances are used to create a synthetic aperture, providing a way to scale up the sensitivity and resolution of radio observations^2,3. This can intuitively be thought of as an extension of the regular resolution formula for optical observations:

For an observing dish of size D meters, and an observing wavelength of λ meters, Θ_HPBW is the angular size in radians of the Half Power Beam Width (HPBW), defining the resolution on the sky. This process of synthesizing a fraction of a large full dish with only scattered points across a mostly empty area is also called aperture synthesis. In the realm of radio interferometry, the resolution of an array is determined by the furthest distance between any two receivers in the array, and this distance is used as D in Equation 1.

The mathematics behind interferometry has been well documented in classic texts like Thompson’s Interferometry and Synthesis in Radio Astronomy³. The basic insight can be communicated informally as “(for planar arrays observing a small field of view) the cross correlation of signals between any 2 receivers (a visibility) will yield information about a 2D Fourier coefficient of the sky brightness pattern.” What Fourier mode is sampled depends on the separation of the receivers (the baseline), normalized by the observing wavelength. Receivers that are further apart (in the standard UVW coordinate system oriented towards the imaging target) sample higher spatial frequency features, yielding higher resolution details at smaller scales. Conversely, receivers that are close together in the same UVW frame sample lower spatial frequencies, giving information of larger scale structures at a lower resolution.

For the lowest radio frequencies, free electrons in Earth’s ionosphere prevent radio waves below 10 MHz from travelling from space to the ground, and vice versa. This so-called “ionospheric cutoff” has long prevented ground-based observations of the sky for this frequency range. The obvious answer to this limitation is to put radio receivers into space where they can record data free of the influence of Earth’s atmosphere and free electrons in its ionosphere. This has been done before with single antennas on spacecraft like Wind⁴ and STEREO⁵, which have revealed many astrophysical processes that produce emissions in this low frequency radio range. This includes emissions from the interactions of electrons with the Earth’s magnetosphere, electron acceleration from solar eruptions, and from the galaxy itself. Single antenna observations can measure the total flux density of such events, but cannot pinpoint where the emission is coming from. In order to localize this low frequency emission and make images in this frequency regime for the first time, many antennas will have to be sent to space and have their data combined to make a synthetic aperture.

Doing this would open a new window through which humankind can observe the universe, enabling a number of scientific measurements that require images of the sky in these lowest frequencies. The Moon is one possible site for a synthetic aperture in space, and it comes with pros and cons when compared to free flying orbiting arrays. The lunar far side has a unique radio quiet zone that blocks all of the usual interference coming from man-made signals, while the near side provides a static place for Earth observing arrays, and if constructed at the lunar sub-Earth point, the Earth will always be at the zenith of the sky. With a static array, it is easier to obtain short baselines to measure large scale emissions, as they are in no danger of colliding, unlike free flying arrays. The drawbacks of a lunar array are chiefly difficulties in cost and power. A large-scale array on the Moon would require a substantial amount of infrastructure and money, while smaller orbiting arrays would require far fewer resources. There is also the issue of power; most places on the Moon are exposed to sufficient sunlight for solar power generation for 1/3 of each lunar day. Surviving the large swings in temperature from lunar day to night is also an engineering concern. Putting aside these difficulties, there is still the problem of making sure that the proposed array design is suitable for its specified science target(s). The response of any given array is dependent on the structure of the emission being observed along with the configuration and sensitivity of the array.

Several conceptual arrays to go on the lunar surface have been drawn up over the decades. Early designs were not the most detailed, but still recognized the scientific advances that could be attained by such arrays^6,7,8,9,10. More arrays have also been put forth in recent years, some of which, like FARSIDE¹¹, DEX¹², and DALI¹³ seek to measure the absorption troughs of the redshifted neutral hydrogen 21-cm signal in the 10-40 MHz range to probe the so called “Dark Ages” and constrain cosmological models of the early universe. Others like ROLSS¹⁴ call out tracking bright solar type II radio bursts far into the heliosphere to identify the site of solar energetic particle acceleration within coronal mass ejections as their compelling science case. Smaller scale arrays have also been described like the 2-element interferometer RIF¹⁵, which would use a single lander and a moving rover to sample many baselines as it moves outward from the lander. RIF focuses on the ability to make a sky map of these low frequencies for the first time, and calculates the UV coverage and synthesized beam for integrated observations.

Space-based radio arrays could also enable low frequency imaging of distant radio galaxies to determine magnetic fields and astrometric measurements¹⁶. Low frequency images of these bodies would provide a more complete picture of the physics governing these systems, in particular yielding synchrotron emission data for the lower end of the electron energy distribution. There are also a range of various magnetospheric emissions that occur at these low frequencies, providing both global (constant synchrotron emission) and local (bursts, auroral kilometric radiation) signatures of electron dynamics that are not detectable from the ground¹⁷. The brightest recorded emissions of these types have come from Earth and Jupiter, as these are the nearest planets with strong magnetospheres. However, arrays with sufficient sensitivity and resolution could observe magnetospheric emission from other outer planets, or even extrasolar planets¹⁸. This topic in particular was called out as an area of interest at the recent Planetary Sciences Vision 2050 workshop.

This work focuses on simulating the response of radio arrays on the Moon consisting of anywhere from just a few antennas, to hundreds or thousands of antennas. This simulation framework is useful for iterating array design for imaging any given scientific target in a small field of view (a few square degrees) but does not currently support all sky imaging. Accurate estimates of the predicted brightness maps along with realistic noise profiles must be used to ensure that a given array size/configuration is sufficient to observe the target to a certain noise level or resolution. The geometry of the array must also be known to a high degree so that the baselines are computed accurately to enable correct imaging of the data. Currently, the best maps of the Lunar surface are Digital Elevation Models (DEMs) from Lunar Reconnaissance Orbiter’s (LRO’s)¹⁹ Lunar Orbiter Laser Altimeter (LOLA)²⁰. The simulation pipeline accepts longitude latitude coordinates for each receiver and interpolates the elevation at these points from existing DEMs to calculate the full 3D position.

From these coordinates the baselines are computed and inserted into a Common Astronomy Software Applications (CASA)²¹ Measurement Set (MS) file. The MS format can be used with many existing analysis and imaging algorithms, and holds information about the array configuration, visibility data, and alignment with the sky. However, many of these software routines are hard coded to work with arrays that rotate with the Earth’s surface, and do not work for orbiting or Lunar arrays. To circumvent this, this pipeline manually calculates the baselines and visibilities for a given array and imaging target, and inserts the data into the MS format. The SPICE²² library is used to correctly align the lunar and sky coordinate systems and track the motions of the Moon, Earth, and Sun.

The simulation framework described here follows Hegedus et al.¹⁷, and the software is archived by the University of Michigan library in the Deep Blue archive²³, stored at https://deepblue.lib.umich.edu/data/concern/data_sets/bg257f178?locale=en. Any patches or updates to this archived software can be found at https://github.com/alexhege/LunarSynchrotronArray. The following section will describe the requirements for this software, and walk through the process of forming an array, setting the appropriate noise levels, feeding the array a simulated truth image of the targeted emission, and simulating the array’s noiseless and noisy reconstructions of the emission using a CASA script.

Protocol

1. Software setup

1. First, go to https://deepblue.lib.umich.edu/data/concern/data_sets/bg257f178?locale=en and download the software package. This software has only been tested in a UNIX environment, and may not fully function in other environments. The README in this package will help guide one through the rest of the software needed and its uses.
2. Make sure python 2.7 or greater is installed. A link is provided in the README. Several common python libraries are also needed including numpy, matplotlib, pylab, scipy, subprocess, ephem, and datetime.
3. Make sure CASA 4.7.1 or greater is installed. A link in provided in the README.
4. Make sure gcc 4.8.5 or greater is installed. A link is provided in the README.
5. Make sure the C toolkit for SPICE is installed. This software is used to align different astronomical reference frames and track the relative positions of planets, moons, and satellites. A link to download this software is also included in the README.
   1. Download several kernels that contain information on astronomical and lunar reference frames, as well as the orbital dynamics of the Moon, Earth, and Sun. The specific kernels needed are listed in the README alongside a link of where to download them.
6. Obtain the final prerequisite data needed: Digital Elevation Models (DEMs) of the lunar surface created from LRO LOLA measurements. The specific file needed is listed and linked in the README.

2. Creating the array configuration

1. Customize the createArrayConfig.py script.
   1. Choose the configuration of the array by providing a list of Longitude and Latitude coordinates for each antenna.
      NOTE: The script is currently formatted for a 10 km diameter array with 1024 elements, 32 arms with 32 log space spaced antenna each, using a constant factor to convert between meters and degrees of longitude/latitude near 0 degrees latitude. The site of the array, (-1.04°, -0.43°), was chosen because it is the center of the 10x10 km patch with the lowest elevation variation (σ = 5.6 m) close to the sub-Earth point (0°, 0°) in the Moon ME frame.
   2. Change the lunarPath variable in the script to reflect the new download location of the Digital Elevation Model containing the elevation data of the lunar surface.
2. Run the createArrayConfig.py script with “python createArrayConfig.py”. This will use the lunar Digital Elevation Model to solve for the elevation at each longitude and latitude for each antenna. Save the longitude, latitude, and elevation to files and print to the screen for easy copying and pasting into the next script. Make figures showing the array configuration on top of the local lunar topography (Figure 1).

3. Using SPICE to align coordinates

1. Customize the eqArrOverTimeEarth.c script.
   1. Take the output from the previous script, the Longitude, Latitude, and elevation of each antenna and copy them over into the corresponding lists in the script, also updating the variable ‘numsc’ with the number of receivers and corresponding coordinates.
      NOTE: Since C does not have dynamic array allocation, there was no easy way to flexibly read in the data automatically, so manual copying must be done.
   2. Update the lunar_furnsh.txt included in the package with the new path names for the required frame and ephemeris files.
   3. Specify what set of dates to observe on. This will inform the ephemerides within SPICE to accurately track where the Earth and Sun are in relation to the defined array for those dates. In the script currently 48 dates occurring roughly weekly over the year 2025 are selected.
   4. Specify the targeted area of the sky for the array to track and image. Currently the script saves the RA Dec of the Earth as seen from the lunar surface, but one may easily just put in static RA Dec coordinates instead.
2. Compile the eqArrOverTime.c script
   1. Compile the script using the gcc command in the comment at the top of the script. It will be something like “gcc eqArrOverTimeEarth.c -o eqArrOverTimeEarth -I/home/alexhege/SPICE/cspice/include /home/alexhege/SPICE/cspice/lib/cspice.a -lm -std=c99”. Change the paths to reflect where the cspice libraries are located.
3. Run the eqArrOverTime executable with “./eqArrOverTime”. This should result in a number of files each with a set of variables in them. Most important are the XYZ position of each antenna in J2000 coordinates, and the Right Ascension and Declination (RA and Dec) coordinates of the targeted area in the sky (currently those of the Earth from the Moon’s perspective). The output variables are saved to .txt files containing the data for all the requested dates.

4. Using CASA to simulate array response

1. Customize the LunarEarthPicFreqIntegration.py script.
   1. Specify the observing frequency for the array to make an image at. This is currently set to 0.75 MHz.
   2. Specify a CASA compatible truth image (or create from a .fits image file) with Jansky/pixel values for the array to reconstruct (e.g., Figure 2). Constants (res, res1, width, arcMinDiv) in the code will need to be changed to reflect the size and resolution of the input truth image.
      NOTE: If using the SPICE method to provide the RA Dec coordinates, one may comment out the ‘import ephem’ statement in this script. This library requires the use of casa-pip from the casa-python package to install, but
      allows for tracking of other astronomical objects within python.
2. Run the LunarEarthPic.py script. Commented at the top of the script are examples on how to run the script. The following command is one example on how to run the script from the command line:
   “nohup casa --nologger --nologfile --nogui --agg -c LunarEarthPicFreqIntegration.py -outDir . -correlate True -numSC 1024 | tee earth.out &”
   The -numSC flag is used to inform the code how many antenna/receivers are being used, and helps unpack the data from the .txt files containing the receiver coordinates.
   NOTE: The antenna baseline vector, measured in units of the observing wavelength (λ), has length D_λ and components (υ, ν, w) = (∆x,∆y,∆z)/λ . The pipeline then calculates the visibilities, or observed cross correlated voltages for each pair of antennas. Here the small field of view approximation is used to calculate the visibilities, following the standard formula from Thompson et al.² for an infinitesimal bandwidth at frequency ν.
   
   The sky coordinates of the target the array is imaging is deemed the phase center, to which the z, or w, axis of the frame is pointed. (l, m, n) are the direction cosines from the (U, V, W) coordinate system. The sky brightness pattern around the source under observation is I_ν(l, m). Spectral flux density is often presented in the derived unit 1 Jansky (Jy) = 10−26 W/m²/Hz. Spectral brightness is simply Jy/steradian to represent the amount of flux coming from a particular area in the sky. A_ν(l, m) is the normalized antenna primary beam pattern, or how sensitive it is to radiation coming from that point in the sky.
   This script calculates the antenna separations in the appropriately projected reference frame from the coordinates output from the previous script. It then uses Equation 2 to calculate the visibility data for each pair of antenna. The resulting visibilities are stored alongside the baselines in a CASA Measurement Set file (.ms). This MS file is the primary output of this script.

5. Imaging the data – noiseless and noisy

1. Customize the noiseCopies.py script.
   1. Set the System Equivalent Flux Density (SEFD), referred to as avNoise in the script. The SEFD is a convenient way to talk about the total noise of a radio antenna since it ties in both the system temperature and the effective area, and provides a way to directly compare the signal and the noise. It is currently set to 1.38e7 Jansky, which is an optimistic noise level for 0.75 MHz.
      NOTE: For the low frequency radio regime, there are three main sources on constant noise: amplifier noise, quasithermal noise from free electrons (estimated by Meyer-Vernet et al.²⁴ to be 6.69e4 Jy at 0.75 MHz, using a electrically short dipole approximation), and Galactic background radiation from the Milky Way (estimated by Novacco & Brown²⁵ to be 4.18e6 Jansky at 0.75 MHz for the full sky, of which a lunar array will only see some portion). This optimal noise level of 1.38e7 Jy assumes that amplifier noise dominates the other terms. See Hegedus et al. for a more detailed discussion.
   2. Set the bandwidth being integrated over in variable ‘noise’ line 200. Set to 500 kHz.
   3. Set the integration time in variable ‘noise’ line 200.
2. Run the noiseCopies.py script with “nohup casa --nologger --nologfile --nogui --agg -c noiseCopies.py | tee noise.out &”.
   1. The script will first create an image from the noiseless visibility data, calling standard radio astronomy algorithm CLEAN²⁶ to create an image like Figure 3.
   2. The script will then create copies of the MS and add the appropriate noise level to the complex visibility data and image it using CLEAN. The script currently makes images for a range of integration times up to 24 hours and over several robust weighting scheme values. Depending on the configuration of the array, image quality may vary with the choice of data weighting schemes. These noisy images will look something like Figure 4, which used an integration time of 4 hours.
      NOTE: The noise is added with standard Signal to Noise formulas. From Taylor² the interferometric noise for a single polarization is
      
      Here, η_s is the system efficiency or correlator efficiency, which has been set to a conservative value of 0.8. N_ant is the number of antennas in the array (N_ant = 2 for each individual visibility), ∆ν is the bandwidth being integrated over in Hz, and ∆t is the integration time in seconds.

Discussion

Each step of the simulation pipeline is necessary and feeds into the next, taking an array configuration on the lunar surface, aligning the reference frame correctly to orient the array to the target area in the sky, calculating the visibility data, adding the appropriate noise levels, and running imaging algorithms on the resulting data.

For each step, customizations may be made. In step 2, the user defined array configuration may be any list of longitudes and latitudes. This then feeds into the SPICE script in step 3, where one may choose the exact time of the planned measurements, as well as where in the sky the array should be focused. In step 3, one can specify the simulated truth emission that the array is attempting to image by providing a suitable CASA .truth file. Then in step 4 one may change the expected level of noise depending on the observing frequency and expected hardware capabilities. This set of codes constitutes a flexible simulation framework that may be used to iterate array design for any number of uses, depending on the targeted science. These codes can all be run on an average laptop or workstation, though computation time increases with number of antennas. The slowest parts of the process are predicting the visibilities, followed by imaging. For small arrays, the entire process can be done in minutes, while for larger arrays of a few hundred or thousand receivers, hours or days may be needed.

Some next steps that could be taken with this pipeline to increase its realism include adding a channel-dependent foreground removal system. This requires building up a global sky model, dominated at low frequencies by galactic synchrotron emission and a few bright sources like Cas A, tracking which part of the sky is visible to the receivers, and convolving that brightness pattern with the primary beam, with the phase center of the array aligned towards the imaging target. For longer integration times, tracking the apparent motion of the sky is also an issue. Another improvement that could be added is a transient event/radio frequency interference (RFI) flagging system that can remove flagged channels from normal imaging, and send them to a specialized pipeline that images and characterizes the flagged data. This transient event pipeline could then use special algorithms like uvmodelfit that can take advantage of the high signal to noise ratio of these events to characterize them better than the normal resolution of the array²⁷.

There are also additional effects that need to be taken into account for a full array calibration, one of which is mutual coupling. As discussed in Ellingson²⁸, this can lead to a decrease in sensitivity in arrays if they have receivers that are within a few wavelengths of each other. This is seen in a decrease in sensitivity for the array, or equivalently, an increase in the SEFD. This is especially true for beams greater than 10 degrees away from zenith. The example array in this work targets Earth, which is always near zenith by design, so mutual coupling should not affect this particular imaging target, but studies of the SEFD over the full range of elevation angles and frequencies will need to be done in commissioning for any real array to unlock its full potential. Another shortcoming of this array simulation pipeline lies in the imperfect lunar surface maps used. DEMs from LRO LOLA measurements have at best a resolution of 60x60 meters/pixel in the 512 pixels/degree maps. One can interpolate these data for simulated arrays, but for real arrays there will need to be a commissioning/calibration period where sources with a known position will be used to determine the relative separations between all the antenna to high precision. Possible calibration sources include Cas A, periodic low frequency emission from Jupiter or Earth, or potentially the Lunar Gateway²⁹.

There is also the response of the lunar surface to consider. There is a layer of lunar topsoil called the regolith that acts like a lossy dielectric that can reflect incoming emission with some efficiency, above the lunar bedrock which can also reflect incoming emission with some better efficiency^30,31. This response is dependent on the ambient temperature and incoming frequency, as well as the chemical composition of the regolith. Studies^30,31 have found that at lower temperatures below 100 K, the regolith is nearly transparent to radio emission, and reflection occurs at the bedrock level with a reflection coefficient of around 0.5-0.6. At higher temperatures 150-200 K, the regolith can absorb emission and reflect incoming radiation at the surface with a reflection coefficient of around 0.2-0.3. At temperatures above 200 K, it is found that the dielectric properties of the regolith are diminished, and variation from reflection can be ignored. These effects can reduce the effective area of the array, reducing sensitivity and requiring longer integration times. This effect can be modeled with electromagnetic simulation software packages such as NEC4.2³² given models of relative permittivity/dielectric constant as a function of lunar depth. This will output the SEFD of a receiver for a given frequency, which can be given to the array simulation pipeline to calculate the correct noise to add to the simulated signal. Adding a grounding grid between the receiver and the lunar surface may help diminish the effect of reflected waves, but adds its own set of complications in the form of deployment.

Many of the hypothetical or fuzzy details around the implementation of a radio receiver on the lunar surface will finally solidify into reality with recent funding of single low frequency antenna projects like Radio wave Observations on the Lunar Surface of the photoElectron Sheath (ROLSES) and the Lunar Surface Electromagnetics Experiment (LuSEE)³³. LuSEE was recently funded by NASA as part of the Commercial Lunar Payload Services program. Both suites of antenna will consist mainly of flight spares for past instruments like STEREO/WAVES or PSP FIELDS and are planned for a 2021 delivery. Measurements from these receivers will finally solidify the level of quasithermal noise from the photoelectron sheath from ionized dust on the lunar surface and how it changes over the course of a lunar day. These measurements will also characterize the level of reflection and absorption from the lunar surface, and quantify how it changes the SEFD of the receiver. They will also provide statistics on the number of transient events or RFI that are received on the lunar surface. These missions will pave the way for arrays of antennas that will finally be able to make a multitude of novel scientific observations such as low frequency emission from solar radio bursts, far away galaxies, and planetary magnetospheres. The simulation pipeline described in this work provides a flexible way to iterate the design of these future arrays for a variety of scientific targets.

Materials

Name	Company	Catalog Number	Comments
No physical materials are needed, this is a purely computational work.