2.1 The data and imaging/calibration pre-processing

We use data obtained by the MWA for the GLEAM-X survey, already calibrated and imaged for the purpose of searching for transient radio sources. The GLEAM-X survey observations were structured as $\sim$ 100 drift scans, whereby the MWA beam was pointed at the local meridian and a fixed declination and over the course of a night the sky drifted through the beam. During a night, the MWA’s 30.72 MHz contiguous bandwidth was switched between five centre frequencies so that five drift scans were produced, covering a frequency range of 72–231 MHz (centre frequencies of 87.675, 118.395, 154.235, 184.955, and 215.675 MHz, hereafter referred to 88, 118, 154, 185, and 216 MHz). The data were collected in 120 s blocks, each identifed as an MWA Observation ID (ObsID). An ObsID contains data at only a single frequency band.

For each ObsID, data were delivered to us already calibrated and imaged using WSClean, (Offringa et al. Reference Offringa2014) providing images for each MWA polarisation (XX and YY) that were averaged over the full 30.72 MHz bandwidth and at a cadence of 4 s. Further, without correction for the primary beam, the XX and YY images had been added together to produce pseudo Stokes I images. The images had been cleaned and the output was stored in an HDF5 format that contains header information and metadata, clean models, and residual images. Each HDF5 container has a size of approximately 300 MB. Our starting points for the analysis in this paper are the HDF5 files. We use data collected during the night of 2020-10-10 (ObsIDs ranging from 1286363152 to 1286398432) for the $\delta=1.6^{\circ}$ drift scan, which covers the Equatorial region and thus the maximum number of GEO objects. The data consist of 254 HDF5 files and approximately 74 GB data volume. Table 1 contains relevant imaging parameters at each frequency. A more detailed description of the data processing steps leading to the HDF5 formatted data products is given in Ross et al. (Reference Ross2024). For a description of the HDF5 format itself, see Appendix 1 of Morgan et al. (Reference Morgan2018).

Table 1. Imaging parameters.

All original data are discoverable at the MWA ASVO archive portal,Footnote ^b where they can be downloaded in raw form (original visibilities, not the processed form of data used here).

2.2 Data processing and image stacking

Using recent measurements of unintended emissions from LEO satellites as a guide (Grigg et al. Reference Grigg, Tingay, Sokolowski, Wayth, Indermuehle and Prabu2023), in the tens to hundreds of Jy range, comparable emitted powers from objects in GEO would be in the mJy regime. Signals of this level would not be included in the clean models in the HDF5 containers. Also, at these levels, RFI flagging applied during initial calibration steps will not excise these signals. So, we take the approach of utilising the residual images, with image stacking to decrease the noise and increase the sensitivity of the stacked images to the point that mJy-level signal detection is plausible. Implicit in this approach is the assumption that any unintended emission fully occupies the 30.72 MHz bandwidth of the observation, as well as the full temporal period of the images being stacked. As unintended emissions from objects in LEO have been observed to be broadband in nature and constant in time, this is a reasonable assumption. In the 4 s images, the motion of the sky across the field of view (manifested as the drift of a GEO satellite across the image) corresponds to approximately 1 arcmin, approximately equal to one pixel in the images. Thus, minimal loss of signal due to smearing will be experienced.

All steps described henceforth were implemented in a Python script.Footnote ^c In order to stack images over time, centred on the location of an individual satellite, information on the predicted location of the satellite is required. We obtained this information via an API query to the space-track.org database of space situational awareness information. The query requested return of Two Line Elements (TLEs) for all tracked objects that have a semi-major axis within 5% of the geostationary distance and an ellipticity $\lt$ 0.01, with TLE generation dates contained in a one week window centred on 2020-10-10. TLEs for 1 396 unique objects were returned.

From the HDF5 containers for ObsIDs at a given frequency, the residual images (FITS format) were extracted, along with their time stamps. From the time stamps, predicted Right Ascension and Declination values for a satellite, as observed from the MWA location, were generated using the Skyfield python module (Rhodes Reference Rhodes2019).Footnote ^d Using the World Coordinate System (WCS) extracted from the residual image FITS header information, if the satellite was found to be within the imaged MWA field of view, a 100 $\times$ 100 pixel cutout of the image, centred on the predicted satellite coordinates, was generated (positions needed to be more than 50 pixels from an edge of the image).

Before being added to a stack, each individual cutout image was corrected for the response of the MWA primary beam.Footnote ^e As the residual images are psuedo Stokes I, simply the addition of the orthogonal XX and YY polarisation images, the beam model used to correct the cutouts was constructed from the average of the XX and YY beam models, generated for the appropriate frequency (centre frequency of each 30.72 MHz band) and pointing direction. The reciprocal of the beam response in the direction of the object was used to multiply all pixels in the cutout, before adding the cutout to the stack for that object. After all residual images at all times and in all relevant ObsIDs were processed, the stack was divided by the number of cutouts that contributed to the stack.

It should be noted that the stacks are produced in pixel space and the individual cutouts are not rotated and re-sampled, according to the instantaneous direction of motion of the satellite under analysis, before adding to a stack. When tracking a satellite from horizon to horizon, such rotations and re-sampling are required (Prabu et al. Reference Prabu, Hancock, Zhang, Tingay, Hodgson, Crosse and Johnston-Hollitt2022). For geostationary and geosynchronous satellites over a limited field of view near the celestial equator, however, the rotations and smearing have less impact. We found that for all 162 satellites ultimately under consideration, all signals were smeared over areas less than two pixels in width, leading to a small maximum ( $\lt$ 10%) reduction in peak intensity due to smearing. This has been disregarded in the subsequent analysis, in favour of not introducing additional complexity into the data processing for minimal gains.

In order to verify that the translation from Right Ascension and Declination to pixel coordinates via the WCS, and the stacking of cutout images, was working correctly, we performed the same set of steps but for two known strong extragalactic radio sources, 3C33 and 3C118. In this case, all processing steps were identical, but as strong radio sources their signals were captured in the clean models, so we used the combined CLEAN model plus residual image rather than just the residual image. Although these astronomical sources do not move in celestial coordinates as objects in GEO do, the drift of the sources through the images partially mimics the behaviour of objects in GEO (a geostationary object drifts through celestial coordinates at the sidereal rate). In both cases, and at all frequencies, the stacked cutouts of 3C33 and 3C118 produced the expected stacked signal at the expected centre pixel.

As the field of view of the MWA is frequency dependent, a larger number of objects was found to be present in the images at the lowest frequency, where the field of view is largest. We found that 162 objects were present for at least some of the time in the field of view at 88 MHz, with 119 at 118 MHz, 93 at 154 MHz, 76 at 185 MHz, and 63 at 216 MHz. Thus, a total of 513 stacks were generated across all frequencies.

2.3 Identification of candidate detections

The RMS of the 100 $\times$ 100 pixel image stacks was measured, to determine detection thresholds. Further, since the TLE information for satellites in GEO can have an up to 40 km position error (Racelis & Joerger Reference Racelis and Joerger2018), corresponding to approximately 4 arcmin, candidate detections were thus considered in a region of this size, centred on the predicted position. As the pixel scale is frequency dependent, this region corresponded to $\pm$ 4 pixels in x and y directions at 88 MHz, $\pm$ 5 pixels at 118 MHz, $\pm$ 7 pixels at 154 MHz, $\pm$ 8 pixels at 185 MHz, and $\pm$ 9 pixels at 216 MHz.

Assuming Gaussian statistics, indicative detection thresholds can be explored. First, however, it should be recognised that neighbouring pixels in the stacks are not independent (as extensively discussed by Hancock, Trott, & Hurley-Walker (Reference Hancock, Trott and Hurley-Walker2018) in the context of fitting source models to interferometric images), with the number of independent samples in a stack approximately equal to the number of pixels divided by the number of pixels per synthesied beam, which is approximately 10 for our data. For example, the error region at 88 MHz of 8 $\times$ 8 pixels yields approximately six independent samples. Adopting a 3 $\sigma$ detection threshold would predict approximately one detection per stack due to thermal noise fluctuations, with approximately 0.008 detections per error region per stack at this frequency. Across the 162 stacks for all objects, one would thus expect approximately one false positive due to noise fluctuations. The same calculation at our highest frequency of 216 MHz (18 $\times$ 18 pixels and 62 stacks) predicts approximately three false positives. Adopting 4 $\sigma$ as the detection threshold reduces the expected false positive rates to below unity at all frequencies. Thus, a 4 $\sigma$ detection threshold may be reasonable for the stacks at each frequency. We also consider detections above a 3 $\sigma$ threshold, looking for detections at the same location across multiple frequencies at this lower level of significance. However, it also needs to be noted that, for these two tests, the data are not likely to be perfectly Gaussian distributed in reality, making this type of analysis an approximation.

Therefore, while we undertake threshold testing of possible detections at 3 $\sigma$ (across multiple frequencies) and 4 $\sigma$ (single frequencies), we also examine other detection metrics. We use a simple metric that asks if there is an excess in the incidence of the brightest pixel in a stack being in the TLE error region. For example, for a 100 $\times$ 100 pixel stack with an error region of 8 $\times$ 8 pixels, the chance that random fluctuations produce the brightest pixel in the error region is approximately 0.006. Across 162 stacks, approximately one detection in the error region may be expected. An excess beyond a given number of detections in the error region can be ascribed a probability via binomial statistics. In python we use scipy.stats.binom.sf(x,n,p) to assess the probability of the number of detections greater than some value, x, given a number of trials, n, and the probability per trial, p. For example, for the $n=162$ trials at 88 MHz, with $p=0.006$ , a value of $x\geq3$ has a probability of 0.04. Thus, detection rates of three or above per 162 stacks can be ruled out as due to random fluctuations with 95% confidence. The test does not indicate, however, which individual detections may be identified with satellites in the case of an excess. The test is robust to the correlation between pixels noted above, but a caveat on the application of this test to this particular dataset, is that the stacks may not, strictly speaking, be fully independent. For two satellites in our field of view that have similar trajectories across the sky, image information will be common to both stacks and therefore not be fully independent. Also, it is worth noting that it is possible that real but weak astronomical radio sources appear in the stacks, at the higher sensitivity levels the stacks achieve.

Table 2. Summary of statistical test results. The 4 $\sigma$ detection column gives the NORAD ID number, detection signal-to-noise, peak intensity of detection (mJy/beam), mean distance (km), and Equivalent Isotropic Radiated Power (EIRP; W) for each detection. $N_{ER}$ is the number of maxima that occur in the error region and $P_{ER}$ is the probability of that number occurring by chance, as described in the text (with the $+$ and $-$ indicating the stacks and the negated stacks, respectively.

Figure 1. MWA imaged field of view in azimuth and elevation coordinates (blue line), with the MWA beam response as a colour greyscale. Red points are the coordinates of individual satellites at individual times. Frequencies are 88–216 MHz, from top left to bottom.

Figure 2. The detection stacks for NORAD 43432 at 88 MHz (top left), 15236 at 118 MHz (top right), and 27683 at 154 MHz (bottom). Shown in red are the 4 arcmin error region adopted in the search process and the specific error regions estimated for the individual satellites (shown using the larger of the RA and DEC region values in Table 3 as a circular error region.).

We also run these analyses after multiplying the pixel values in the stacks by -1. For detections generated by random fluctuations, there should be as many negative fluctuations as positive fluctuations. Thus, the statistical results in all three tests should be consistent if detections are due to random fluctuations.

For completeness, we also ran full processing of the satellite list stacking the CLEAN model plus residuals, rather than just using the residual images. The results were barely discernible from each other, which indicates that the temporal occupancy of CLEAN model components in our cutout regions is very low. This is understandable as the cutouts are $\sim$ 10 arcmin on a side, which represents less than one minute in time at the Equator. We thus consider the residuals-only processing for further discussion.

3.1 Individual candidates ( $\gt4\sigma$ events within error regions)

At 88 MHz, 1/162 objects exceeded the 4 $\sigma$ threshold, NORAD ID 43432 (Table 2). This object is listed as COSMOS 2526, launched by the Commonwealth of Independent States from the Tyuratam Missile and Space Complex on April 18, 2018; it is in a geostationary orbit. While definitive information is difficult to obtain, the following source claims the satellite to be a communications satellite for the Russian Ministry of Defence (https://weebau.com/satcosmos/2/2526.htm). At this frequency, the test for excess detections in the error region gave a result that is borderline significant, with a probability due to noise fluctuations of 0.05.

At 118 MHz, 1/119 objects exceeded the 4 $\sigma$ threshold, NORAD ID 15236 (Table 2). This object is listed as LEASAT 2, launched by the United States from the Air Force Eastern Test Range on August 30, 1984 via the Space Shuttle; it was one of a series of communications satellites owned by Hughes Aircraft Company and leased to the US military. A detailed description of the LEASAT satellites can be found here: https://en.wikipedia.org/wiki/Syncom.

At 154 MHz, 1/93 objects exceeded the 4 $\sigma$ threshold, NORAD ID 27683 (Table 2). This object is INTELSAT 907, one of the series of Intelsat satellites deployed to support communications over North and South America, Western Europe, and Africa; it was launched from French Guiana on February 15, 2003 (https://www.n2yo.com/satellite/?s=27683).

At 185 MHz, 0/76 objects exceeded the 4 $\sigma$ threshold.

At 216 MHz, 0/63 objects exceeded the 4 $\sigma$ threshold. At this frequency, our test for excess detections in the TLE error region yielded four detections out of the 63 stacks, which is borderline significant, with a probability due to noise fluctuations of 0.05.

Figure 2 shows the stacks for NORAD numbers 43432, 15236, and 27683.

While we can adopt a statistical error region for the purposes of our detection algorithms, it is possible to look in detail at the three candidates these algorithms produce, using error region estimates specific to these objects. If we retrieve TLE data from space-track.org for the three candidate objects, for a two week period centred on 2020 October 10, we can assess the spread in positional predictions from these TLEs, for a particular point in time (2020-10-10:00:00:00 UTC). This allows us to estimate specific error regions for each object, that we can then compare to the images in Figure 2. Table 3 describes the resulting estimates.

Table 3 shows that the error region estimates for the three candidate objects are far smaller than the 4 arcmin $\times$ 4 arcmin regions assumed above, based on a statistical 40 km positional error. In particular, when compared to the images in Figure 2, it can be seen that the 4 $\sigma$ signals lie outside these smaller error regions. This is strong evidence that the signals represent statistical fluctuations rather than any origin from these satellites.

The spreadfor NORAD 15236 is significantly larger than for the other two objects, as 15236 is an active operational satellite that manoeuvres often, such that any given TLE will be incorrect. The other two objects are in graveyard orbits, such that the objects gradually drift, meaning TLEs will be accurate for longer.

Figure 3. EIRP upper limits ( $4\sigma$ ) at 88–216 MHz, from top left to bottom right. The horizontal axis is labelled ‘Order’, which is simply the order in which TLEs for satellites were returned from the space-track.org query.

3.2 Further statistical tests

Across all five frequency bands and 513 stacks, three candidate detections above our 4 $\sigma$ threshold are apparent, one each at 88, 118, and 154 MHz. At each frequency a different satellite was identified. This result is inconsistent with the expected false positive rates calculated above, which are set to be below unity. However, we note the caveats we have made when it comes to Gaussian statistics. Thus, we explore further tests.

We note that the smaller field-of-view at higher frequencies reduces the number of opportunities for detection across multiple frequencies. Looking at the three detections listed above in individual frequency bands, COSMOS 2526 (NORAD 43432) was detected at 88 MHz and was in the field-of-view at 118 MHz (peak intensity of 17 mJy/beam and signal-to-noise of 2.8), but was not in the field-of-view at 154 MHz, 185 MHz, or 216 MHz.

LEASAT 2 (NORAD 15236) was detected at 118 MHz and was present in the field-of-view at 88 MHz (peak intensity of 57 mJy/beam and signal-to-noise of 1.6), 154 MHz (peak intensity of 18 mJy/beam and signal-to-noise of 2.9), and 185 MHz (peak intensity of 40 mJy/beam and signal-to-noise of 2.6). LEASAT 2 was not in the field-of-view at 216 MHz.

INTELSAT 907 (NORAD 27683) was detected at 154 MHz and was in the field-of-view at 88 MHz (peak intensity of 87 mJy/beam and signal-to-noise of 3.5), 118 MHz (peak intensity of 10 mJy/beam and signal-to-noise of 2.4), 185 MHz (peak intensity of 6 mJy/beam and signal-to-noise of 2.1), and 216 MHz (peak intensity of 53 mJy/beam and signal-to-noise of 1.7). INTELSAT 907 was therefore in the field-of-view at all frequencies, but only detected at one frequency.

When we lower our thresholds and test for instances of 3 $\sigma$ detections in multiple frequency bands at the same location, we find five instances of $\gt3\sigma$ detections within the TLE error region common to two frequency bands, but no instances where these detections are coincident in position within the error region.

When running all tests after negating the pixel values in the stacks we find 4 $\sigma$ detections at 88 MHz (two), 185 MHz (one), and 216 MHz (one), similar to the non-negated results in terms of detections. The instances of maxima in the error regions are listed in Table 2 and are also consistent with the non-negated results. Finally, in the negated images, we find six instances of $\gt3\sigma$ detections in the TLE error region with no coincident positions, again consistent with the result above for the non-negated images.

Table 3. The Right Ascension and Declination spreads (RMS) calculated from TLEs available in a two week period for the three candidate detections at a reference time of 2020-10-10:00:00:00 UTC.

Considering all of the above, the caveats on the use of Gaussian statistics and the consistency of the statistical results between the stacks and the negated stacks, there is no significant evidence for detection of emission from any of these objects. While we show the three 4 $\sigma$ detections in Figure 2, this is mainly to give the reader an idea of the quality of the stacked images.

We calculate upper limits on the Equivalent Isotropic Radiated Power (EIRP) based on our observed intensity upper limits, for all objects at all frequencies, and show them in Figure 3. From these results we see that, in general, we reach 4 $\sigma$ upper limits on EIRP of well below 1 mW, and below 10 $\unicode{x03BC}$ W in the best cases (see details for each satellite in the Appendix).

We conclude that RFI due to unintended emissions from satellites in geostationary and geosynchronous orbits is currently likely to be undetectable at typical MWA sensitivities, based on the assumptions that we have made (emission of constant strength filling both the observing bandwidth and observing time). There is currently little indication that this type of RFI will be an issue in the same way that unintended emissions from objects in LEO is becoming an issue.

For comparison to previous work, detections of unintended emissions from objects in LEO by Grigg et al. (Reference Grigg, Tingay, Sokolowski, Wayth, Indermuehle and Prabu2023) at levels of $\sim$ 200 Jy, in an approximately 1 MHz bandwidth, correspond to EIRPs of a few $\unicode{x03BC}$ W. In the most sensitive of our data, we are approaching 4 $\sigma$ EIRP limits at comparable levels for objects in GEO, a factor of $\sim$ 70 further away. In terms of comparison of our data to industry Electromagnetic Compatibility (EMC) standards, Di Vruno et al. (Reference Di Vruno, Winkel, Bassa, Józsa, Brentjens, Jessner and Garrington2023) notes that ‘currently no international agency or space law … requires a spacecraft to comply to a certain EMC standard’. A range of standards can be used for comparison to our data, to give a feel for the data in an EMC sense. We choose CISPR 32 Class A, for which the emission limit is 40 dB(uV/m) at a distance of 10 m, across the frequency range of 30–230 MHz, and measured in a 120 kHz band. From this limit, we calculate the EIRP we would expect for the bandwidth used for our observations (assuming that the power density in frequency is the same for the 120 kHz band and our full observing bandwidth). We calculate an equivalent CISPR 32 Class A EIRP of $1.5\times10^{-2}$ W, for our bandwidth, which is well above all the 4 $\sigma$ upper limits presented in Figure 3. The measurements we make indicate that any UEMR from GEO satellites we have observed is well below CISPR 32 Class A limits.

This is a welcome overall conclusion and the current study provides a useful baseline for comparison to future detailed monitoring of the RFI environment of the site of the SKA-Low, particularly as the SKA-Low sensitivity reaches and then significantly exceeds the current MWA sensitivity. An additional factor of 50 sensitivity in the final SKA-Low would bring EIRP sensitivity to GEO objects into at least the $\sim$ 100 nW range. Given the emerging issues with objects in LEO, the situation for objects in GEO and other orbital domains should be monitored regularly into the future.

While we assume continuous, broadband signals for the purposes of deriving our limits, limits are also possible in the case that any signals are narrow band in nature. In this case, a narrow band signal would be diluted by the full observed bandwidth such that our upper limits would be higher by the factor $\frac{32}{\Delta \nu}$ , where $\Delta \nu$ is the width of the narrow band signal in MHz. The GLEAM-X visibility data would be amenable to a possible future search for narrow band RFI from GEO satellites, as the inherent frequency resolution of the data is 10 kHz. An analysis that leads to 10 kHz images, rather than 30.72 MHz images, is possible, but would require a factor of approximately 2 000 more compute resource.