Imaging sensitivity of a linear interferometer array on lunar orbit

At present, full sky maps in the frequency range of 1 - 30 MHz are not available. We therefore use a sky map generated by extrapolation from higher frequencies. A number of extrapolated sky models have been developed, such as the Global Sky Model (GSM) (de Oliveira-Costa et al., 2008; Zheng et al., 2017), and the Cosmology in the Radio Band (CORA) (Shaw et al., 2014, 2015). In the present work, we make use of the high-resolution Self-consistent whole Sky Model (SSM) (Huang et al., 2019). Note that the SSM model used by the current work has not accounted for the absorption effect that would become very significant at $\nu<10\,{\rm MHz}$. A model which incorporates the free-free absorption effect has been developed by Cong et al. (2021). For the purpose of evaluating the map-making capability, however, the original SSM model is fairly adequate, and maybe even more preferable as it gives a familiar sky map that is easier to compare at different wavelengths. This model may yield a slightly higher sky brightness temperature than the actual one, as it ignores the interstellar absorption.

In the full-sky simulation below, we pixelize the map with a resolution of about $1^{\circ}$, corresponding to HEALPIX pixels with $N_{\rm side}=64$. The input sky map at 10 MHz is shown in Fig. 1.

We shall consider a nearly circular orbit for the array of satellites. In choosing the orbit parameters, we keep in mind three aspects of the problem: (1) the stability of the orbit; (2) the good observation time during which the Earth (and the Sun) is shielded. (3) the nodal precession of the orbit which generates a three dimensional distribution of baselines, which is important for breaking the mirror symmetry in image synthesis (Huang et al., 2018). Some of these may conflict with each other, for example, a lower orbit will allow a larger part of the sky being blocked, thus increasing the good observation time, but the lower orbit is also less stable, which would require more orbit maintenance to ensure the safety of the array. We have to weigh these factors and reach a compromise in our choice.

In order to evaluate the merit of each orbit choice, we carry out orbit simulations. We select the GL1500E lunar gravity field model (Kahan, 2013), in which the order of the model is $100\times 100$, and the JPL DE405 Planetary ephemeris²²2https://ssd.jpl.nasa.gov/?planet_eph_export is used. In Fig. 2, we plot the perigee height of three orbits of different average height over a period of 5 years. We see that as the orbit precesses, there is a variation in its perigee height. For the orbit of 300 km average height, the minimum height of the perigee is 180 km during the the procession cycle, while for the 200 km average, the the minimum height of the perigee is 38 km. On the other hand, the mean half angle extended by the Moon is given by

where $R=1737.4$km is the mean radius of the Moon, so $\theta=58.51^{\circ}$ and $63.74^{\circ}$ for the 300 km and 200 km respectively. We have tentatively chosen a 300 km as the fiducial value of the orbit height for the project.

The Moon revolves around the Earth in an orbit which is only slightly inclined ($\sim 5^{\circ}$) with respect to the ecliptic plane, and the lunar equator is only $1.543^{\circ}$ inclined w.r.t. the ecliptic plane, so we can approximate both with the ecliptic. Due to orbital dynamics, the orbital plane of the lunar satellite array will precess. This is very important for the imaging synthesis, because for the linear array the baselines are all distributed on the orbit plane, without precession the orientation of this plane will be fixed and there is a mirror symmetry with respect to this plane. The precession allows the plane to rotate and orient toward different directions, resulting in a three dimensional distribution of baselines, which is essential to break the mirror symmetry. Shorter orbit precession cycle allows faster filling of the three dimensional $uvw$ space. In Fig. 3, we plot the precession period (reading number from the left axis), for different orbital inclination angles for the 300 km orbit. The precession rate decreases with the inclination angle, the cycle is longer than 2 years for orbits with inclination angle $>55^{\circ}$, and unsafe for $60^{\circ}$. On the other hand, although the precession is faster at lower inclination angles, the resulting baseline distribution would be compressed into a disc in the $uvw$ space if the inclination angle is too small, limiting the angular resolution in the direction perpendicular to the ecliptic plane.

Another factor to be considered is the percentage of the good observation time, i.e. the time when the Earth and/or the Sun is shielded by the Moon from the view of the satellites. Here we take a geometric optics definition of shielding, though due to diffraction some RFI may still be observed at the edge of the shielded region. The time when the Sun and the Earth are both shielded is the intersection of them. This fraction is maximized at zero inclination, and decreases with increasing declination. This is because at zero inclination angle, the Earth, the Moon and the satellites are all on the same plane, so there is always time when the satellites are shielded when they orbit the Moon, but for large inclination angle the satellites will spend some time outside the plane of the Moon’s orbit, during which they are more easily exposed. In Fig. 3, with the right axis, we plot the fraction of good observation time for a 300 km orbit with different inclination angles during five years, and we also plot the “double good” time when both the Sun and Earth are shielded.

Based on these considerations, we choose a 300 km circular orbit with an inclination angle of $30^{\circ}$ as our fiducial orbit. For this orbit, the precession cycle period is about 1.3 year, the good observation time that the Earth is shaded is about $31\%$, and the “doubly good” observation time both Earth and Sun are shielded is about $10\%$.

The thermal noise on the measured visibility in temperature units is

where $\Delta\nu$ is the bandwidth of the channel, and $t_{\text{int}}$ is the integration time. Based on the preliminary design of the DSL, the basic channel width is set at $\Delta\nu=8\,{\rm kHz}$, but note that multiple channels may be combined to increase the sensitivity.

The system temperature $T_{\mathrm{sys}}$ receives contributions from both sky radiation and receiver noise,

Here $T_{\rm{rcv}}$ is the effective temperature of receiver noise, and $T_{\rm sky}$ is the mean brightness temperature of the sky, as average over the beam of the antenna. In the Rayleigh-Jeans limit, the total energy flux received by a dipole antenna is given by

where $A(\theta,\phi)$ is the antenna beam pattern, $G$ is the gain which equals to 1.5 for short dipoles, and the factor 1/2 in first equation is due to the fact that the antenna responds to only one polarization. The design of the receiver for the DSL low frequency interferometry system will be presented elsewhere, but based on preliminary design, each polarization of the DSL low frequency band system is made up of a pair of monopole antenna. The signal from each antenna is individually amplified before combined and digitized. The receiver noise on each side is given in terms of equivalent voltage fluctuation at its input stage, with $\approx 3\,{\rm nV}/\sqrt{\rm Hz}$ over the observation band. As the noise on the two sides are uncorrelated, the combined noise is

This is to be compared with the voltage induced in the antenna. The equivalent energy flux of the receiver noise, that is to be combined with the flux from the sky (Eq.(4)), is given by

where $Z_{0}=377\,\Omega$ is the impedance of the vacuum, and $l_{\rm eff}$ is the effective length of the dipoles. Then the effective receiver temperature is

In this work, we assume the length of the dipoles is 5 m, corresponding to $l_{\rm eff}\approx 2.5\rm m$. Thus, the receiver noise is equivalent to

or

These are plotted in Fig. 4. This receiver temperature increases as $\lambda^{2}$, while the Galactic synchrotron radiation, in the optical thin regime, scales as $\lambda^{2.6}$, until at very low frequency when absorption becomes important. Hence at the low frequency band considered here, the sky radiation dominates the system temperature even for the relatively short antennas of DSL.

The satellites will orbit the Moon with a period of $\sim 2.3$ hour. The longer baselines are moving fairly fast, and the integration time for a visibility measurement is limited to the baseline crossing time. For a given observing frequency and a specific baseline, the change in the projected baseline length should be less than $\sim 1/10$ of the wavelength within the integration time, in order to avoid the smearing in the $(u,v,w)$ sampling. Therefore, the integration time of any visibility measurement is restricted to

where $\omega$ is the angular velocity of the satellites, and $D$ is baseline length. For the DSL orbit, the maximum integration time is 0.013 seconds for the 100 km baseline at 30 MHz.

The sky temperature and noise temperature at different frequencies are shown by the dotted and solid lines respectively in Fig. 4. Note that although the longer baselines have much higher noise due to the shorter integration time, they have much more measured visibility data points.

The visibility data $V_{ij}$ measured by an interferometer array elements pair $i,j$ is related to the sky brightness $I(\bm{\hat{k}})$ by

where $A_{ij}(\bm{\hat{k}})$ is the combined antenna beam pattern that depends on the antenna primary beam response, $I(\bm{\hat{k}})$ is the sky brightness in the direction $\bm{\hat{k}}$, $\textbf{{k}}=\frac{2\pi}{\lambda}\bm{\hat{k}}$ is the wave vector, and $\textbf{{r}}_{ij}=\textbf{{r}}_{i}-\textbf{{r}}_{j}$ is the baseline vector formed by antennas $i$ and $j$. The visibility function can be re-written in the $uvw$ coordinates given in wavelength units $\textbf{{u}}=(u,v,w)\equiv\textbf{{r}}/\lambda$,

where $(l,m,n)$ are the direction cosines with respect to the three coordinate axes, and $l^{2}+m^{2}+n^{2}=1$.

For ground-based arrays with small FOV, the $w$-term in the exponential can often be neglected, such that the integral in Eq.(12) can be reduced to a 2D Fourier transform, and with the beam $A(l,m)$ known, the sky intensity can be easily recovered by inverse Fourier transform. For wider FOV, the image can still be recovered by a 3D Fourier transform, and also a number of algorithms were developed for such wide-field imaging (Thompson et al., 2017). However, for the case we are considering here, these algorithms are not suitable for several reasons. First, there is no ground to shield half of the sky, and by using electrically short dipole antennas, nearly the whole sky is within the FOV. For data taken during one orbit, where the baselines are distributed on a circular ring on the orbital plane, there is a mirror symmetry problem: a source located on one side of the plane produce exactly the same visibility as a source located at the mirror image position on the other side of the plane, the two can not be distinguished (Huang et al., 2018). Second, the precession of the orbital plane produces a three-dimensional distribution of baselines. While this breaks the mirror symmetry, the $w$-term could never be neglected. Third, the Moon will block a fraction of the sky, but at each location on the orbit, the blocked part is different. Strictly speaking, the visibility data at each time corresponds to a sky with a different blockage, so even the 3D Fourier transform method is not applicable.

However, if we neglect the reflection and radiation of the Moon, the visibility is still linearly related to the sky intensity, so in principle the sky intensity can be recovered by solving this time-dependent linear equation (Huang et al., 2018), which can work in spite of all the issues listed above. The Moon blockage can be introduced as a shade function into the visibility, i.e.

where $S_{ij}$ is the shade function that describes the time-dependent positions in the sky blocked by the Moon. We pixelize the sky with the HEALPix scheme (Górski et al., 2005), and the integral over sky angles is discretized as a sum over sky pixels,

where $\alpha$ is the pixel index, $N_{\rm pix}$ is the total number of pixels, $\Delta\Omega=4\pi/N_{\rm pix}$ is the solid angle corresponding to one pixel, and $B(\alpha,t)$ is the complex response of the array, which is an effective “beam” taking into account of both antenna beam response and the Moon blockage, i.e.

Converting the sky brightness $I$ to the brightness temperature $T$ according to the Rayleigh-Jeans approximation, and considering the thermal noise, the observed visibilities can be written in a matrix form as

With a total of $N_{\mathrm{t}}$ observation time points and $N_{\mathrm{bl}}$ baselines, $\mathbf{B}$ is a $(N_{\mathrm{t}}\cdot N_{\mathrm{bl}})\times N_{\mathrm{pix}}$ matrix, $\mathbf{T}$ is a $N_{\mathrm{pix}}$-element vector with each element corresponding to one pixel in the sky, and the visibility $\mathbf{V}$ and noise $\mathbf{n}$ have the dimension of $(N_{\mathrm{t}}\cdot N_{\mathrm{bl}})$. In the following, we adopt the convention $\sum_{\alpha=1}^{N_{\rm pix}}A_{ij}(\alpha)\,\Delta\Omega=1$, so that the visibility has the same units as the sky brightness temperature.

We model the noise as a Gaussian random variable with zero mean, i.e. $<n>=0$, variance $\sigma_{n}^{2}$ (given by Eq. (2)), and the noise covariance matrix $\mathbf{N}\equiv\langle\mathbf{nn}^{\mathrm{H}}\rangle$. The minimum variance estimator (for Gaussian noise, also the maximum likelihood estimator) of $\mathbf{T}$ is then

Note that $\mathbf{B^{\mathrm{H}}N^{-1}B}$ is a $N_{\mathrm{pix}}\times N_{\mathrm{pix}}$ matrix. However, if $\mathbf{B^{\mathrm{H}}N^{-1}B}$ is not invertible, which is often the case due to limited sampling in the $uvw$ plane, we can still compute some psuedo-inverse matrix, e.g. Moore-Penrose pseudo-inverse (Huang et al., 2018):

More generally other estimators of $\mathbf{T}$ can be formed with different choices of $\mathbf{D}$ (Dillon et al., 2015). The imaging quality can be quantified by the point spread function (PSF) of the reconstruction, and the covariance of the estimator. The expectation of the above estimator is

where

is the matrix-valued point spread function. If the matrix $\mathbf{B^{\mathrm{H}}N^{-1}B}$ is invertible we would have the idealized PSF $\mathbf{P}=\mathbf{I}$. Note that $\mathbf{B^{\mathrm{H}}N^{-1}B}$ fully characterizes the observation, and may be regarded as the dirty beam of the system. If we simply set $\mathbf{D}$ as a normalization matrix in Eqs. (19)-(20), then we obtain the so called “dirty image” without deconvolution. This PSF of the dirty beam shows the basic characteristics of the specific array configuration.

The covariance of the estimator is

In the case of sufficient $uvw$ sampling with an invertible $\mathbf{B^{\mathrm{H}}N^{-1}B}$ matrix, the covariance $\mathbf{C\to(B^{\mathrm{H}}N^{-1}B)^{-1}}$. In our Moore-Penrose pseudo-inverse approach, the matrix $\mathbf{B^{\mathrm{H}}N}^{-1}\mathbf{B}$ is approximately the inverse of the covariance, which measures the information content of the reconstructed map. In Section 4.2 we will calculate and analyze the PSF and eigenvalue spectra of the $\mathbf{B^{\mathrm{H}}N}^{-1}\mathbf{B}$ matrix for our different array configurations.

In the above, we have assumed the satellites are at rest with each other. Indeed, for an array of satellites moving on the same circular orbit and located close to each other, the velocity differences are small and this is a good approximation. However, we will make sky map by synthesizing the data collected over a period of several years, during this time, the orbital motion of the satellites, Moon, and Earth induce varying velocities with respect to the inertial frame, so here the effects induced by the motion need to be considered.

In the presence of velocity, up to the first order, $\nu^{\prime}=\nu(1-\bm{\beta}\cdot\bm{\hat{k}})$, and $\textbf{{k}}^{\prime}=\textbf{{k}}/(1+\bm{\beta}\cdot\bm{\hat{k}})$, where $\bm{\beta}=\textbf{{v}}/c$. The visibility is given by

Simple estimates show that the circular motion around the Sun has a magnitude of $v\sim 30\,{\rm km}s^{-1}$, while the lunar motion around the Earth ($v\sim 1\,{\rm km}s^{-1}$) and the satellite motion around the Moon ($v\sim 1.7\,{\rm km}s^{-1}$) are smaller. For the heliocentric motion which has the largest effect, $\beta\sim 10^{-4}$. For the continuum radiation, this induces a correction which is much smaller than the current precision of the measurement, so at first approximation, we can take $I(\nu^{\prime},\textbf{{k}}^{\prime})\approx I(\nu,\textbf{{k}})$. However, the correction on the phase is not negligible, so we have

The correction $(\bm{\beta}\cdot\bm{\hat{k}})(\textbf{{k}}\cdot\textbf{{r}}_{\text{ab}})\sim\mathcal{O}(1)$ for the longest baseline of 100 km at wavelength of 10 m. Note that in ground-based observation, this effect also exists, but it is usually corrected in real time when making the observation as the field of view is usually small and this phase factor can be directly compensated for a given direction. In our case, the visibility receives contribution from the whole sky, so the correction can not be directly applied to the data. However, this phase is also deterministic, and one can include this phase factor in the mapping matrix $\mathbf{B}$, and correct for the velocity effect in the map reconstruction procedure.

Another possible problem is that if there is a spectral line from an object in a certain direction, and if we make imaging observation using narrow channels, i.e. the channel width is smaller than or comparable with the line width, then the line may be shifted outside the channel. However, if the source of such a line emission is known, one can make the observation by making real-time phase correction on the data for the given direction, just as is used in ground-based observations. In principle, one can also apply such a correction for every sky direction to obtain the whole sky map for the narrow spectral line, though that requires much more computation.

With the electrically short antennas which are omnidirectional, and without a ground shielding one-half of the sky, there is a mirror-symmetry problem (Huang et al., 2018). However, because of the nodal precession of the orbit, we can naturally acquire a three dimensional (3-D) distribution of baselines to break this mirror-symmetry. As an example, the baselines generated by 5 satellites are shown in Fig. 5, which has an onion-like layered structure.

The number of the daughter satellites considered here is in the range of $5\sim 8$, to be determined by the mass of the payload and the capability of the launching rocket. After the number of satellites is decided, there are also different options for the array configuration, i.e. the relative positions and baseline lengths of the satellites. As a start, we set the minimum distance between any of the daughter satellites to be fixed as 1 km, which provides an ample safety margin against satellite collision. The maximum distance is set to be 100 km, within which precise angular position measurement can be achieved straightforwardly with small star sensors, and at the same time it also provides an angular resolution which approaches the static imaging resolution limit imposed by interplanetary and interstellar scintillation (Jester & Falcke, 2009).

We consider four configurations of the linear array as follows, and investigate the imaging quality in each case.

a. Even spacing. In this case the daughter satellites are equally spaced, except for the spacing between the first and second satellites which forms the shortest baseline of 1 km. The distance between the $i$-th satellite and the first satellite is

where the $N$ is total satellite number.

b. Random spacing. In this case, after fixing the first two satellites and the last one to ensure the minimum and the maximum baselines, we randomly insert the other satellites between the second and the last satellite with the minimum distance between any two satellites being larger than 1 km. One such realization is used, which we refer to as the random spacing³³3The position of the 8 satellites are set at 0 km, 1.0 km, 2.89 km, 10.17 km, 24.23 km, 43.23 km, 69.75 km, and 100.00 km with respect to the first satellite, respectively. , though strictly speaking the word “random” should refer to a distribution of baselines, not a single case.

c. Logarithmic spacing. We also test the logarithmic spacing between the satellites, in which the position of the $i$-th daughter satellite is set at

The basic idea here is to produce baselines on all different scales.

d. Optimized spacing. The optimized spacing configuration is motivated by the consideration of achieving an $uvw$-coverage on all different scales as in the logarithmic case, but with an emphasize on having more long baselines. We propose a modified logarithmic spacing configuration, in which the the $i$-th daughter satellite ($2<i<N$) is placed at a distance $L_{i}$ away from the first one, with

and $L_{2}=1\,{\rm km}$, $L_{N}=100\,{\rm km}$.

The length distribution of the 28 baselines formed by 8 satellites are shown in Fig. 6 for the above four configurations. While in all four cases the lengths are distributed across the whole range from 1 km to 100 km, the “optimized” (d) and “random” (b) spacing cases appear to spread out more smoothly. The even spacing (a) configuration has a high degree of redundancy, and the logarithmic spacing (c) strategy produce more short baselines.

In this work we shall model the motion of the satellites as simple circular motion on a slowly rotating (precessing) orbital plane, and the satellites are fixed on the assigned place. The actual orbital motion of the satellites may be more complicated, as the orbit may not be a perfect circle, and is subject to the perturbation of the anomalous lunar gravitational field, and each satellite may have slightly different orbital parameters. However, these would not change the overall imaging results, as long as the relative positions between the daughter satellites are measured with sufficient precision. The slightly different orbital parameters may actually allow the satellites to sample a slightly large part of $uvw$ space. Indeed, one may intentionally generate relative motions between the satellite, so that they periodically approach and then back away from each other. Such a “breathing mode" would increase the $uvw$ space sampled over the operation period, but here we shall not consider these more complicated cases.

Given the array configuration, the coordinates of each satellite during the orbital motion is specified, and the varying baseline vector $\textbf{{r}}_{ij}(t)$ formed by each pair of the satellites can be calculated. We shall compare the imaging results produced with the different array configurations in section 4.

The computational cost of the imaging simulation is dominated by the multiplication and inversion of the matrices. For a single frequency, the beam matrix $\mathbf{B}$ in Eq. (15) has a dimension of $(\mathrm{N_{bl}}\cdot\mathrm{N_{t}},\mathrm{N_{pix}})$. For long-term observations, $\rm N_{\rm bl}N_{\rm t}\gg\mathrm{N_{pix}}$, and it is beyond the capability of most current computers to compute the pseudo-inverse of matrix $\mathbf{B}$. Therefore, we choose to reconstruct the sky map through Eq.(17), by computing $\mathbf{(B^{\mathrm{H}}N^{-1}B)^{-1}}$ instead of $\mathbf{B^{-1}}$ as in Huang et al. (2018) by the Moore-Penrose pseudo-inverse, then the cost of the matrix inversion scales as $\mathrm{N_{pix}}^{3}$. However, the computation of the matrix $\mathbf{(B^{\mathrm{H}}N^{-1}B)}$ scales as $(\mathrm{N_{pix}}\mathrm{N_{bl}}\mathrm{N_{t}}+\mathrm{N_{pix}}^{2}\mathrm{N_{bl}}\mathrm{N_{t}})$ (where $\mathrm{N_{t}}N_{\rm bl}\gg\mathrm{N_{pix}}$), then the matrix multiplication dominates the computational cost in the entire algorithm.

In order to make the computation feasible, we divide the entire observation period into multiple time steps in simulating the imaging. For example, imaging at 1.5 MHz can be performed every 15 days, then totally 30 time-steps in one precession cycle (1.3 years) are required.

Even with such a step-by-step imaging approach, we can produce the full-sky map only at a resolution of $1^{\circ}$ (corresponding to $N_{\mathrm{side}}=64$ in Healpix), which is cruder than the resolution of the array, e.g. 6.5 arcmins at 1.5 MHz for the maximum baseline of 100 km, because the computational cost for the inversion of the matrix $\mathbf{(B^{\mathrm{H}}N^{-1}B)}$ scales as $\mathrm{N_{pix}}^{3}$, which is quite demanding beyond this resolution. We list the full-sky imaging simulation requirement on computing resources in Table 2. One will need more computing resources to achieve higher resolutions.

In order to overcome the difficulty in achieving the resolution corresponding to the maximum baseline of the DSL array, we may apply the map-making algorithm to a small patch of the sky. This is detailed in the next subsection, and we can then estimate the point source sensitivity of the array.

For the maximum baseline of 100 km, the resolution is 10.31 arcmins at 1.0 MHz (corresponding to $N_{\rm side}=512$ and $3\times 10^{6}$ pixels in the whole sky), and 0.344 arcmins at 30 MHz (corresponding to $N_{\rm side}=16384$ and $3\times 10^{9}$ pixels in the whole sky). However, with the huge number of pixels in the full-sky map, the computation needed for reconstructing the full-sky map is quite formidable. The computation can be greatly reduced by restricting the simulation and reconstruction to a small patch of the sky. Now in the real observation, the visibility will also receive contribution from the sky outside the patch. In the reconstruction, this may produce some error in the reconstructed map. We shall call this the “sidelobe effect". It is not due to the side lobe of the antenna, whose main lobe is much wider than the region of concern, rather it simply refers to the fact that in the part-sky reconstruction the radiation outside the patch may also affect the reconstruction result. We assess this side-lobe effect by checking the consistency between the imaging results for a small region of the sky, with the corresponding sky region from the full-sky imaging.

For a small patch of the sky, we set the pixels outside this patch of sky to be 0, i.e.

For a total of $N_{\mathrm{t}}$ observation time points, $\mathbf{V_{p}}$ is still a vector of $(N_{\mathrm{t}}\cdot N_{\mathrm{bl}})$, but the computational cost of $\mathbf{B_{p}}$, $\mathbf{(B_{p}^{\mathrm{H}}N_{p}^{-1}B_{p})}$, and its inverse is largely reduced, as $N_{\mathrm{pix}}$ is now given by the pixel number of this patch of sky. The expected value of the estimator is:

Fig. 9 compares the results of the full-sky reconstruction from the entire orbit precession period (1.3 year) at $1^{\degree}$ resolution (panel (a)) with the part-sky reconstruction for the region in the red rectangle (right ascension (RA): [34^∘, 74^∘], declination (DEC): [-38^∘, 2^∘]) under the same observation conditions (panel (b)). In order to eliminate the sidelobe effect, we conservatively use the central quarter of the simulated sky area to estimate the array sensitivity in the part-sky simulation. The relative errors are plotted in panel (c). The results obtained by the two approaches are quite similar, with a relative difference less than 5%. This shows that the sidelobe effect is not significant for this area of the sky, and the part-sky reconstruction approach is feasible.

We shall start with a simple, idealized case, and then add practical issues one by one to see how each factor affects the result.

(a) The Ideal Case. This is the case that is free of noise and Moon blockage, and the antennas are assumed to have isotropic response. These conditions are similar to those adopted in Huang et al. (2018), though there the visibility data were randomly sampled in the orbit-allowed region of the $uvw$ space, whereas here they are generated from the orbital motion with a given array configuration.

(b) Noise included. Based on the ideal case, we first add a Gaussian random noise, with the standard deviation of $\sigma_{\rm n}$ given by Eq. (2), to the mock visibility data, so that the sensitivity will be limited by the noise level.

(c) Effect of Moon blockage. Inevitable for a lunar orbit array, the Moon blockage makes the visibility dependent on the position of the satellites on the orbit. We have noted earlier that this effect has to be taken into account in the reconstruction. Furthermore, it also reduces the effective observation time. Here we take into account the time-varying Moon blockage.

The reconstructed sky maps for these cases are shown in the upper panels of Fig. 10 for the 10 MHz simulations. We also show in the bottom panels in each figure the corresponding relative error maps, defined as $\varepsilon=(T-\hat{T})/T,$ where $T$ and $\hat{T}$ are the brightness temperatures of a pixel in the input sky map and the corresponding pixel in the reconstructed map, respectively.

In all cases, the visual impression is that the sky is reasonably reconstructed, the error is only prominent at some pixels. As expected, in the ideal case (a) the reconstructed map has very low errors. The presence of noise (case (b)) visibly increases the relative error in the map. The error increases further with the introduction of the Moon blockage (case (c)), which is expected because at least it reduces the observation time for each pixel.

To evaluate the global quality of the maps more quantitatively, we introduce a few commonly used statistics, such as the Mean Squared Error (MSE), the R Squared value (RS), and the Structural SIMilarity (SSIM). The MSE is a representation of the absolute error, while the RS value and the SSIM are normalized and quantify relative errors. These are widely used for image comparison. The MSE is defined as

where $I$ and $R$ denotes the values of the input and reconstructed maps respectively. The MSE is the variance in the case of an unbiased estimator. Note that the value of MSE depends on the absolute magnitude of the measured quantity. The RS value is related to the MSE by

Its value ranges from 0 to 1, and is independent of the absolute magnitude of the measured quantity The structural similarity is defined as

where $\bar{I}$ and $\bar{R}$ are the mean values of the input and the reconstructed maps, $\sigma(I)$ and $\sigma(R)$ are the standard deviations of $I$ and $R$ respectively, and $\rm{Cov}(I,R)$ is the cross-covariance between the input and the output. $c_{1}=(k_{1}L)^{2}$ and $c_{2}=(k_{1}L)^{2}$ are very small constants used to maintain stability of the computation (Hengl et al., 2004), in which L is the range of the sky map data. In our analysis, we adopt $k_{1}=0.01$,$k_{2}=0.02$.

The MSE, RS and SSIM values for the cases described above are listed in Table 3. As we depart from the ideal case (a), with more practical issues taken into account, the MSE value increases while RS and SSIM values decrease, indicating that the imaging quality degrades with the inclusion of more practical issues.

In Fig.11 we show the MSE, RS and SSIM values of the reconstructed maps at 1.5 MHz as functions of total observation time, for different numbers of satellites, and different array configurations. Here, we also consider observation time much shorter than the full precession period. Although the precession is incomplete, as shown by the curves, within the span of a few months it is already possible to produce a good whole-sky map. The error level quantified by the MSE drops with the time as $1/\sqrt{t_{\rm obs}}$ as expected, and the quality of the map would be improved continuously with the accumulation of more data with the passing of time. Comparison of the different curves also shows that while the sky map could be produced even with 3 satellites, the noise as quantified by MSE will be reduced as the inverse of the number of satellites, and the image quality as quantified by the RS and SSIM will also be improved with the addition of more satellites.

We also compare the even spacing configuration, random spacing configuration, and the optimized configuration for the 8 satellites case in Fig.11. The result of 8 satellites with the optimized configuration is obviously better than the other two configurations and that of uniform distribution of satellites is the worst.

Due to the limitation of computer memory, in our study of the global image quality we have adopted a relatively low resolution ($1^{\circ}$). We now turn our attention to imaging a small part of the sky and analyze the resolution and imaging capabilities of the array, using the PSF and eigenvalue spectra (Dillon & Parsons, 2016). Our imaging algorithm can also deal with a part of the sky. As an example, we selected two particular patches of sky with a higher angular resolution of 6.7 arcmin (HEALPix $N_{\mathrm{side}}=512$), shown in Fig. 12, i.e. position A at RA $[68^{\circ},76^{\circ}]$, DEC $[-32.8^{\circ},-24.8^{\circ}$], and position B at RA $[248^{\circ},256^{\circ}]$, DEC $[-31.0^{\circ},-23.0^{\circ}$]. Position A is located on the satellite orbit plane, while position B is located away from the orbital plane. Here we will show the PSF at position A to illustrate the influence of orbital precession on map-making, and the PSF at position B to illustrate the optimization of the satellite array configuration for map-making.

In this work, we plot the PSF for the dirty beam, $\mathbf{D(B^{\mathrm{H}}N^{-1}B)}$, where $\mathbf{D}$ is the normalization matrix. In Fig. 13, we show the results of PSF at position A after different exposure times (i.e. 1 orbit (8248.7 seconds), $1/4\ T$, $1/2\ T$, and $1\ T$, from left to right respectively, where $T$ is the precession period) and the projected $vw$ coverage of the corresponding baseline distribution. When the satellite completes a single orbit, the baseline is almost on a plane, and the PSF in the direction on the orbital plane is the linear stripes shown in the first lower panel of Fig. 13. When the orbital precession completes 1/4 period, the linear striped side lobes are obviously weakened, and an incomplete ring-shaped structure appears at the same time. After 1/2 period observation, the ring-shaped stripe tend to be closed and the linear stripe structure is further weakened. After completing the whole cycle of precession, the PSF evolves to what is shown in the last panel in Fig. 13. The side lobes at the border almost disappeared, and an X-shaped structure appeared at the center. This arise from the triangular gap at the upper and lower ends of the $vw$ plane shown in the upper right panel of Fig. 13.

We compute the PSF for the $8^{\circ}\times 8^{\circ}$ facet at position B for 1, 100, and 1000 orbits’ observation with different numbers of satellites, which is shown in Fig. 14. For simplicity, here we do not consider the shielding of Sun and Earth, but simply take all orbital time as the effective observation time. The results of 3 satellites (3 baselines), 5 satellites (10 baselines), and 7 satellites (21 baselines) are shown from top to bottom rows in Fig. 14. Note that the PSF matrix contains all the relevant information about the telescope’s ability to imaging the sky. It is clearly seen that the side lobes are more suppressed with increasing number of satellites. Some of the ripples (or stripes) in the PSF is caused by the under-sampling of baselines in the torus envelope of the precessing orbital plane, and such ripples would be suppressed with longer observation time, as clearly shown in the first row of Fig. 14, but there are also stripes originated from the inadequate sampling in the upper and lower cone regions of the $uvw$ space(see the upper panels in Fig. 13), as the $vw$ (or $uw$) cross-section is blank at the upper and lower ends. This structure will not disappear because it is caused by the peculiarity of the baseline configuration generated by the orbital inclination of 30^∘ and the orbital precession.

In Fig. 15, we plot the PSF for different array configurations. From top to bottom, the four rows show the PSF formed by 8 satellites with even spacing, random spacing, logarithmic spacing, and optimized spacing configurations, respectively. The left and right columns show the PSF after 1 orbit and 1000 orbits observation respectively. It is seen that in the even spacing configuration, the structure of the PSF is composed of multiple concentric rings, and the side lobes are more prominent, showing the effect of having more redundant baselines. The PSF of logarithmic configuration also has prominent side lobes around the center, because it has more short baselines as shown in Fig. 6. Among the four configurations, the optimized spacing has the sharpest PSF, thanks to the largest number of long baselines. We also find that the side lobes in the PSF is reduced considerably by prolonging the observation time.

We also analyze the eigenvalue spectrum of the matrix $\mathbf{B^{\mathrm{H}}N^{-1}B}$, which is approximately the inverse of noise covariance matrix of the maximum-likelihood estimate. The matrix $\mathbf{B^{H}N^{-1}B}$ contains all information about the instrument correlation in our map, but we have an implicit assumption here that the noise is uncorrelated between pixels, so the matrix $N^{-1}$ is diagonal in our simulation. We anticipate that the smallest covariance will be associated with the greatest inverse covariance.

The eigenvalue spectrum of the matrix $\mathbf{B^{\mathrm{H}}N^{-1}B}$ are plotted for different array configurations and observation times in Fig. 16. We plot the eigenvalue spectra for 3, 5, and 7 satellites, labeled as (3s), (5s), (7s) in the legend respectively, after 1, 100, and 1000 orbits in Fig. 16 (a). It can be seen from the eigenvalue spectra that for longer observation times and more baselines one would obtain more independent modes with larger eigenvalues. In Fig. 16 (b), we show the eigenvalue spectra for different array configurations and the different observation times. We find that the logarithmic spacing configuration (c) and the optimized spacing configuration (d) have more modes with high eigenvalues as compared with the even spacing (a) and random spacing (b) configurations. In addition, after a period of observation, the optimized spacing configuration results in the highest eigenvalue spectrum among the four configurations. This suggests that sampling the $uvw$ space for the whole range of scales uniformly is desirable.

The results of PSF and the eigenvalue spectrum are consistent and complementary to each other, both quantifying the imaging ability of the array with different designs and observation times. From the above analysis of the PSF and the eigenvalue spectra, we conclude that the imaging ability of a lunar orbit array can be improved by increasing the number of satellites, extending the observation time, and optimizing the satellite arrangement. The proposed optimized spacing configuration is the preferred option for the DSL mission, as it results in the sharpest PSF and the lowest noise covariance.

In the above, we have assessed the global map-making capability of a DSL-like array by simulating the full-sky imaging with $1^{\circ}$ resolution (HEALPix $N_{\mathrm{side}}=64)$. We now study the sensitivity of the array for point sources by focusing on a small part of sky with higher resolution.

Theoretically, the surface brightness sensitivity of an interferometer array is

where $L_{\rm max}$ is the maximum baseline, $A_{\rm eff}$ is the effective area of the antenna, $N$ is the number of interferometric elements, and $t_{\rm obs}$ is the total observation time. The practical sensitivity of DSL may deviate from the theoretical value due to the practical issues involved. Below, we examine the point source detection from the reconstructed sky map.

We select a dark sky region away from the Galactic plane, for example the region A in Fig. 12, and then apply our synthesis method to obtain the reconstructed sky map with the full resolution of the DSL array for a given frequency. We calculate the local brightness temperature fluctuation $T_{\rm RMS}$ in the selected area, which may vary from region to region. Setting the detection threshold at $5\sigma$, the flux density limit for source detection is

where $\Omega\approx 1.13\,\theta_{\rm o}^{2}$ is the beam solid angle, in which the prefactor of 1.13 is adopted in analogy to a Gaussian beam. However, we do not assume any shape for the beam here, and the beam size in estimating the point source sensitivity is $\theta_{\mathrm{o}}\approx 1.6\times\lambda/L_{\rm max}$, which is determined by the full width of half magnitude of the PSF from our simulation, where $L_{\rm max}=100\,{\rm km}$ is the longest baseline for the DSL mission.

We run simulations for a linear array of 5 satellites in optimized spacing configuration at frequencies of 1.5 MHz, 3 MHz, 6 MHz, 12 MHz, and 24 MHz, with the pixel size of 6.87 arcmin, 3.44 arcmin, 1.72 arcmin, 0.859 arcmin, and 0.429 arcmin, respectively, which are a bit finer than the theoretical resolution of the 100 km baseline. We obtained $T_{\rm RMS}$ and the detection threshold of $S_{\rm min}$ according to the simulated beam size $\theta_{\rm o}$, at different frequencies. As shown in Fig. 17, we find that the simulated $T_{\rm RMS}$ generally agrees with the theoretical expectation of Eq.(35). The results for the flux density limit are plotted in Fig. 18, showing the increase of sensitivity with the increasing observation time. The observing band width is assumed to be 8 kHz according to the current DSL mission concept, which is fairly narrow, but note that the sensitivity may be improved if a larger bandwidth could be achieved with the inter-satellite communication link. The point source sensitivities for different frequencies and different observation times are listed in Table 4.

The expected number of detectable sources depends on the luminosity function and the spectral energy distribution of radio sources at ultra-long wavelengths, both of which are completely unknown below 30 MHz. Tentative radio source counts have been obtained at 41.7 MHz and 61 MHz by the LOFAR array, but it is still not sure if the derived shallow spectral index is due to residual calibration issues or the intrinsic property of sources at these frequencies (Shulevski et al., 2021). Therefore, here we compare three models extrapolated from observations at higher frequencies.

1. VLSS model. Jester & Falcke (2009) extrapolated the faint-end source number count derived from the VLA Low-Frequency Sky Survey (VLSS) at 74 MHz towards lower frequencies, assuming a spectral index of $\alpha=0.7$ that is typical for optically thin synchrotron sources. The number density of sources is given by

where S is the flux density limit, measured in mJy, and $\beta=-1.3$. In this model, the source number distribution function $n(S)$ follows a power-law distribution above a certain threshold.

2. SSM model. The all-sky map predicted by the SSM model is used as the input sky map in our simulation. By combining the data from the NRAO VLA Sky Survey (NVSS) at 1.4 GHz (Condon et al., 1998) and the Sydney University Molonglo Sky Survey (SUMSS) at 843 MHz (Bock et al., 1999; Mauch et al., 2003), the SSM model also includes a point source catalogue covering the entire sky, but given at higher frequencies. The source distribution function at 1.4 GHz is given by (Huang et al., 2019)

The total number of detectable sources in the whole sky is then

Here we extrapolate the source count at 1.4 GHz down to below 30 MHz using a spectral index of $\alpha=0.8157$, which is fitted from sources in the overlap region of SUMSS and NVSS (Huang et al., 2019).

3. GLEAM model. We may also use a more up-to-date source count derived from the GaLactic and Extragalactic All-sky MWA survey (GLEAM) at 200 MHz, 154 MHz, 118 MHz and 88 MHz, covering 24,831 $\mathrm{deg}^{2}$ (Franzen et al., 2019). We apply a weighted least squares fit to the GLEAM source counts at 154 MHz, with a $5^{\mathrm{th}}$ order polynomial given by

The fitted parameters are $a_{0}=3.5230$, $a_{1}=0.2979$, $a_{2}=-0.3986$, $a_{3}=-0.0204$, and $a_{4}=0.0376$ at 154 MHz (Franzen et al., 2019). The total number of detectable sources also follows Eq.(39). Then the source counts are extrapolated to 1$\sim$30 MHz with a spectral index of $0.8$.

In Fig. 19 we plot the number of detectable point sources at several frequencies for a linear array of 5 satellites in the “optimized spacing configuration". The histograms show results for 1 month and 1.3 years of total observation time, without considering the radio frequency interferences from the Sun and the Earth. The predictions for the VLSS model, SSM model and GLEAM model are plotted in the left, middle and right panels respectively. The number is quite dependent on the source count model, which is currently very uncertain. For example, according to the VLSS model, the number of $5\sigma$ detectable sources are 180 (or 590) for one month’s observation (or one precession cycle) at 12 MHz, and drops to 140 (or 466) for one month’s observation (or one precession cycle) at 24 MHz. The GLEAM model predicts more detectable sources at low frequencies, but fewer detectable sources at the high-frequencies, as compared with the VLSS model. According to the GLEAM model, the number of detectable source are 180 (595) sources for one month (one precession cycle) at 12 MHz, but only 44 (280) sources for one month (one precession cycle) at 24 MHz. The SSM model predicts much higher numbers: 10610 (18740) for one month (one precession cycle) at 12 MHz, and 5260 (10700) for one month (one precession cycle) at 24 MHz.

If we only use the data obtained when the Earth is shielded, or when both the Earth and Sun are shielded, the available integration time is further reduced to1/3 and 1/10 of the total orbital time. The total number of detectable sources would be further reduced. The corresponding results for the VLSS, SSM, and GLEAM models are plotted with dots and stars in the figure. In Fig. 20, we show the number of $5\sigma$ detectable point sources for different numbers of satellites (5 – 8). The number of detectable sources increases with the number of satellites as expected.

Note that in the above estimates, we have made several simplifications in our simulations. First, the sky model ignores the absorption effect, which would be important below $\sim 3$ MHz, leading to higher sky temperature and higher flux density limit, underestimating the detectable number of point sources. Secondly, all the above three source models ignore the possible low-frequency turnover in the source spectrum due to synchrotron self-absorption, and the number of sources could therefore be overestimated. Furthermore, here we have not considered the measurement errors on baselines, the time synchronization errors, the calibration errors on both the amplitude and phase of visibilities, etc. More practical considerations will further decrease the sensitivity and reduce the expected number of detectable sources. However, here we have noted the significant variation between the source models, due to the completely unknown nature of sources at the ultra-long wavelengths. A more realistic input sky model and more sophisticated simulations will provide more reliable results of the sensitivity of the DSL array, but the expected number of sources is still quite uncertain.

Above we have discussed the sensitivity limit of point source observation, but whether a source can be detected as a distinct one also depends on the number density of the sources and the angular resolution of the telescope. If multiple sources above the sensitivity limit are located within the same beam, they would appear to be overlapped and indistinguishable, this imposes a confusion limit for point source detection, which usually limit the mean surface density of the point sources to be less than $1/30$ per beam (Väisänen et al., 2001; Condon et al., 2012). In the most extreme case of the SSM model, if the observation is performed with 8 satellites at 1.5 MHz, there are less than 30,000 detectable sources after one precession cycle using only the Earth-shielded observing time, while there are $3\times 10^{6}$ pixels for the whole sky at the observational resolution, it amounts to less than one source in 100 pixels. Note that the results at 1.5 MHz may be severely affected by the free-free absorption effect on both the sky temperature and the source luminosity. At the more reliable frequencies of 12 MHz and 24 MHz, the SSM model predicts about 15,700 and 12,000 detectable sources respectively, after one precession cycle using 8 satellites and the Earth-shielded observing time, while there are $5\times 10^{7}$ and $2\times 10^{8}$ pixels in the whole sky for 12 MHz and 24 MHz, respectively. Note this is the sensitivity limit for a single channel with 8 kHz bandwidth. The sensitivity could be significantly improved for the continuum sources by using multi-channel data. Combining the 30 channels to achieve a bandwidth of $30\times 8$ kHz, the number of $5\sigma$ sources in the SSM model are then about 58,300 and 45,200 at 12 MHz and 24 MHz respectively, after a complete precession cycle using 8 satellites and the Earth-shielded observing time. Thus the confusion limit is not a serious concern in our case.