# Radio Frequency Interference from Radio Navigation Satellite Systems: simulations and comparison to MeerKAT single-dish data

Brandon Engelbrecht, Mario G Santos, José Fonseca, Yichao Li, Jingying Wang, Melis O Irfan, Stuart E Harper, Keith Grainge, Philip Bull, Isabella P Carucci, et al.

► To cite this version:

Brandon Engelbrecht, Mario G Santos, José Fonseca, Yichao Li, Jingying Wang, et al. Radio Frequency Interference from Radio Navigation Satellite Systems: simulations and comparison to MeerKAT single-dish data. Mon.Not.Roy.Astron.Soc., 2024, 536 (1), pp.1035-1055. 10.1093/mn-ras/stae2649. hal-04584573

HAL Id: hal-04584573

https://hal.science/hal-04584573v1

Submitted on 18 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Radio frequency interference from radio navigation satellite systems: simulations and comparison to MeerKAT single-dish data

Brandon N. Engelbrecht,<fnref n="1" />* Mario G. Santos<sup>1,2*</sup> José Fonseca<sup>3,4,1</sup> Yichao Li,<fnref n="5" />,<fnref n="1" /> Jingying Wang<sup>1,6,1</sup>

Melis O. Irfan,<fnref n="7" />,<fnref n="1" /> Stuart E. Harper<sup>1,8</sup> Keith Grainge,<fnref n="8" /> Philip Bull,<fnref n="8" />,<fnref n="1" /> Isabella P. Carucci<sup>1,9,10</sup>

Steven Cunnington<sup>1,8</sup> Alkistis Pourtsidou<sup>1,11,12,1</sup> Marta Spinelli<sup>1,13,1</sup> and Laura Wolz<fnref n="8" />

Accepted 2024 November 20. Received 2024 November 19; in original form 2024 May 7

## ABSTRACT

Radio frequency interference (RFI) is emitted from various sources, terrestrial or orbital, and creates a nuisance for ground-based 21-cm experiments. In particular, single-dish observations will be highly susceptible to RFI due to their wide primary beam and sensitivity. This work aimed to simulate the contamination effects from the Radio Navigational Satellite System (RNSS) within the 1100–1350 (MHz) frequency band. The simulation can be divided into two parts: (i) satellite positioning, emission power, and the beam response on the telescope, and (ii) calibration of the satellite signals to data to improve the original model. We utilize previously observed single-dish L-band data from the Meer-Karoo Array Telescope (MeerKAT), which requires special calibration to account for regions contaminated by satellite-based RFI. We find that we can recreate the satellite contamination with high accuracy around its peak frequencies provided the satellite is not too close to the telescope’s pointing direction. The simulation can predict satellite movements and signals for past and future observations, aiding in RFI avoidance and testing novel cleaning methods. The predicted signal sits below the noise in the target cosmology window in the L band (970–1015 MHz) making it difficult to confirm any out-of-band emission from satellites. However, in our simulations, this contamination still overwhelmed the 21-cm auto-power spectrum. Nevertheless, it is possible to detect the signal in cross-correlations after mild foreground cleaning. Whether such out of band contamination does exist will require further characterization of the satellite signals far away from their peak frequencies.

**Key words:** instrumentation: miscellaneous – methods: observational – telescopes – software: simulations – cosmology: observations – diffuse radiation.

## 1 INTRODUCTION

The integrated emission of the 21-centimetre (cm) line of neutral hydrogen (H1) can be used to map the Large Scale Structure (LSS) of the universe. The method, known as intensity mapping (IM), allows us to scan large areas of the sky without the need to resolve individual galaxies (Bharadwaj, Nath & Sethi 2001; Chang et al. 2008, 2010; Masui et al. 2013; Bull et al. 2015). Detecting this faint cosmological signal can be hindered by radio

frequency interference (RFI) emitters which, due to the technological evolution of communication devices (Bentum, Boonstra & Baan 2010), have been increasingly contaminating the radio bands. These sources consist of terrestrial RFI, such as frequency modulation (FM) broadcasting towers, and orbital RFI, which emanates from satellites (Baan 2010). The main techniques applied to mitigate the effects of RFI have been avoidance, filtering, and flagging (Baan 2019). Radio telescopes are planned and operate within a designated radio quiet region, which aids in restricting contamination. However, as much as we locate radio facilities away from terrestrial RFI, satellites in orbit around the Earth will always contaminate the signal and become one

<footnote n="1">¹Department of Physics & Astronomy, University of the Western Cape, Cape Town 7535, South Africa</footnote>

<footnote n="2">²South African Radio Astronomy Observatory (SARAO), 2 Fir Street, Cape Town 7925, South Africa</footnote>

<footnote n="3">³Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto CAUP, 4150-762 Porto, Portugal</footnote>

<footnote n="4">⁴Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre 687, PT4169-007 Porto, Portugal</footnote>

<footnote n="5">⁵Department of Physics, College of Sciences, Northeastern University, Wenhua Road, Shenyang, 11089, China</footnote>

<footnote n="6">⁶Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, China</footnote>

<footnote n="7">⁷Institute of Astronomy, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="8">⁸Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, The University of Manchester, Manchester M13 9PL, UK</footnote>

<footnote n="9">⁹INAF - Osservatorio Astronomico di Trieste, Via G.B. Tiepolo 11, 34131 Trieste, Italy</footnote>

<footnote n="10">¹⁰IFPU - Institute for Fundamental Physics of the Universe, Via Beirut 2, 34151 Trieste, Italy</footnote>

<footnote n="11">¹¹Institute for Astronomy, The University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="12">¹²Higgs Centre for Theoretical Physics, School of Physics and Astronomy, The University of Edinburgh, Edinburgh EH9 3FD, UK</footnote>

<footnote n="13">¹³Observatoire de la Côte d'Azur Laboratoire Lagrange, Bd de l'Observatoire, CS 34229, 06304 Nice cedex 4, France</footnote>

<footnote>*E-mail: engelbrechtbn@gmail.com (MNE); mgrsantos@uwc.ac.za (MGS)</footnote>
of its largest contaminants. This is of particular concern as satellites transmitting in the L-band frequency range between 1 and 2 gigahertz (GHz), contaminate regions of the spectrum important for H1 IM. The situation will only worsen due to the increasing number of satellites being launched annually (ITU 2021). However, real-time filtering schemes can prevent the leakage of strong RFI signals into the frequencies of interest before correlation (Buch et al. 2023, for example). Specifically, the most common approach is to flag (i.e. remove) post-correlation data contaminated by RFI (Offringa et al. 2012).

As radio telescopes become more sensitive and the sky becomes more contaminated with satellites, it becomes crucial to fully understand the radio signals emitted by them (Di Vruno et al. 2023). Such a study serves three purposes (i) to find parts of the sky and times when observations are less contaminated, (ii) to understand how much the RFI contaminates the target frequencies for our science, (iii) to test whether there are ways to further clean the signal from RFI on those target frequencies (away from the prominent peaks). There are several examples of follow-ups on satellite emissions alongside their impact on observations. The Green Bank Telescope (GBT)<fnref n="1" /> has been characterizing RFI sources at different frequency bands (NRAO 2024). Researchers have used the Five Hundred Metre Telescope (FAST)<fnref n="2" /> to estimate the power of various satellites; to achieve this, an extra antenna was used to monitor nearby satellites (Wang et al. 2021a). Additionally, observations using the Parkes<fnref n="3" /> ultra-wide bandwidth receiver included navigational and telecommunication satellites (Hobbs et al. 2020). A detailed study of the MeerKAT RFI environment as a function of direction, frequency, time of day and baseline has also been done Goedhart (2024).

Satellite emissions can be particularly adverse for single-dish (SD) H1 IM surveys, such as those planned for Baryon Acoustic Oscillations from Integrated Neutral Gas Observations (BINGO) (Battye et al. 2013) and FAST (Bigot-Sazy et al. 2016), as well as the MeerKAT (Santos et al. 2017) and Square Kilometre Array Observatory (SKAO)-MID (Santos et al. 2015) arrays when used as a collection of single dishes. Although interferometers are also affected, the signal is further reduced due to time and frequency smearing (decorrelation) away from the delay centre. A more direct mitigation strategy is to consider frequencies far away from the allocated satellite transmission band. These allocations are set by the International Telecommunications Union (ITU). However, out-of-band transmission remains an issue, as the ITU only provides guidelines for the maximum out-of-band power that can leak from these allocations (ITU 2006, 2015). Given that the H1 cosmological signal is less than a milli-kelvin (mK), it is essential to understand how satellite RFI propagates outside the allocated bands and gauge whether their contamination is negligible.

It is becoming increasingly important to study the impact of RFI emissions due to the cosmological surveys planned for MeerKAT and the SKAO site in South Africa. The situation should be better for surveys using MeerKAT's ultra-high frequency (UHF) band (580–1015 MHz) and SKAO-MID band 1 (350–1050 megahertz [MHz]), but it will be more problematic for MeerKAT's L band (900–1670 MHz) and SKAO-MID band 2 (950–1760 MHz). Harper & Dickinson (2018) (referred to as HD18 hereafter) estimated the impact of global navigational satellite systems (GNSS) on a future

SKAO-MID band 2 H1 IM survey, using models for the total power and spectral structure of GNSS signals convolved with a model SKA beam. They conclude that for frequencies >950 MHz, the emission from GNSS satellites will exceed the expected H1 signal for all angular scales. However, information on the technical specifications of satellites is lacking, namely the signal's amplitude which introduces uncertainties in the modelling. Moreover, the impact on the final H1 IM power spectrum depends on the scanning strategy used. Additionally, one could use cleaning methods to reduce the contamination level or even subtract it, as has been done in Finlay et al. (2023).

In this work, we will expand on HD18's work using a complete catalogue that incorporates radio navigational satellite systems (RNSS) extending beyond the focus of GNSS to include regional and augmented systems. We develop new simulations that are compared and calibrated to calibrated MeerKAT data from Wang et al. (2021b) (W21 hereafter). These data were used to make a recent detection of the cross-correlation power spectrum between H1 and galaxy surveys (Cunnington et al. 2023a). In this study, we focus on satellites with allocated emission bands below 2 GHZ, which should have the largest impact on the target cosmology window (970–1015 MHz). In total, we consider 73 satellites operating in MeerKAT's field-of-view (FoV) during the scanning period of an observation commencing on 2019-02-25 02:40:11 South African Standard Time (SAST). The new set of simulations can accurately estimate satellite contamination over a wide frequency range once properly calibrated with the data. In this work, we will determine in which regimes we can calibrate and characterize the satellite signals. The flowchart in Fig. 1 summarizes the steps taken to simulate the calibrated observed temperature of the satellites' constellations. One can also find the reference to the section in the paper describing them. Additionally, we will use this simulation to understand the satellites' impact on the H1 IM 3D power spectrum and test cleaning methods. These results will help in preparing observations with future surveys with MeerKAT (see e.g. MeerKLASS Collaboration et al. 2024) and the SKAO and set the stage to understand the ultimate level of contamination on the science that will be extracted from such data.

The remainder of the paper is organized into several sections hereafter. In Section 2, we describe how we construct the navigational satellite signals and positions. In Section 3, we address the calibration of the observational data to include the previously flagged RFI region. In Section 4, we discuss the effectiveness of the model. In Section 5, we explore the applications of the model. Finally, in Section 6, we provide a conclusion to the study.

## 2 NAVIGATIONAL SATELLITE SYSTEMS SIMULATION

Artificial satellites serve a variety of purposes, including communication, Earth observation, navigation, and more. In this paper, we focus exclusively on radio navigation satellite systems, including subcategories such as global, regional, and augmented systems. The different systems within these sub-categories are usually referred to as constellations. These constellations cater to and have slight differences between civilian and military uses for any given country. Currently, there are four global navigational satellite systems in orbit which provide positional information to users across the globe. These are the United States of America's (USA) -NAVSTAR (hereafter, GPS); the European Galileo system (GAL); the Russian Global'naya Navigatsionnaya Sputnikovaya Sistema (known as GLONASS, but

<footnote n="1">¹https://greenbankobservatory.org/</footnote>

<footnote n="2">²https://fast.bao.ac.cn/</footnote>

<footnote n="3">³https://www.parkes.atnf.csiro.au/</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Flowchart showing the steps taken to simulate the position, emitted and observed temperature of the constellations. The referred sections in the flowchart describe in detail the particular step.</figcaption></figure>

herein only GLO), and the Chinese Compass navigation system Beidou (BEI).

On the other hand, regional navigation systems operate across a specific country or fixed geography. Here, we consider the Indian Regional Navigation Satellite System (NavIC, but herein referred to as IRNSS) and the Japanese Quasi-Zenith Satellite System (QZS). Additionally, we incorporated the Satellite Based Augmentation System (SBAS) into this work, as this geostationary satellite system benefits the current GNSS by improving positioning accuracy. These constellations follow operational standards, such as allocated frequency bands, which follow the ITU recommendations (ITU 2006, 2015).

Table B2, lists all the navigation satellite systems considered for this work. This catalogue contains information on the aforementioned constellations, including the relevant parameters required to model their spectral signal emissions (see subsequent subsections). We have utilized several sources in the literature to gather this information. However, some details are unavailable or remain uncertain, such as the transmitted power of satellite signals in certain constellations. This has prompted us to allow some of these parameters to vary, which were later calibrated with data (see Section 2.4).

In order to build a simulation of the observed satellite signal, the following information is required:

(i) The frequency range of the observation $\Delta v_j$ [MHz];

(ii) The time-stamps $t$ [sec] of the observation;

(iii) Location of the telescope on the Earth [Latitude (Lat) & Longitude (Long)];

(iv) Observational pointing per timestamp [Right ascension (RA) & Declination (Dec.)];

(v) An all-sky telescope beam model;

(vi) Positional information of the satellites; and

(vii) The power emitted by the satellites as a function of frequency.

The location information, telescope pointing, and beam model are flexible and can be user-specific, allowing the simulation to be applied to various SD radio telescope experiments around the globe. The following subsections provide a detailed explanation of the steps taken to simulate the observed satellite signal.

## 2.1 Satellite position

Each navigational constellation or system comprises several individual satellites. These individual satellites follow a predetermined course at a specific orbital distance. Although celestial mechanics could naively predict the positions of the satellites, several factors complicate calculations over long time-scales, necessitating up-to-date information. The positional vector information for each satellite is contained within a two-line element (TLE) data set on the CelesTrak<fnref n="4" /> database. The data are saved as snapshots for the various satellite constellations. This catalogue is updated every 24 h and does not store previous TLE information. Historical TLE information from months to years back can be accessed via the Wayback Machine.<fnref n="5" /> The Wayback Machine stores snapshots from various days but not concurrently. However, when comparing historical and present snapshots, we found that a TLE can reliably describe a satellite position within three months. Therefore, when selecting historical snapshots, we use the one closest to the observation dates.

Once the necessary TLE information has been acquired, we apply the PYTHON package SKYFIELD<fnref n="6" /> (Rhodes 2019) to load the TLE vectors and retrieve the orbital positions of the individual satellites, such as Lat and Long; RA and Dec. or azimuth and altitude (elevation) for a reference position on Earth. In our case, we took the position to be the central location of the MeerKAT array; Long: 21° 26′ 38″ E and Lat: 30° 42′ 47.41″ S. Given that the diameter of the MeerKAT array is much smaller than the orbital height of the navigational satellites, each antenna will observe the satellites at approximately the same angular distance from the beam centre.

The final ingredient is the angular distance, $\theta$, between the satellite and the pointing direction of the telescope, which can be calculated using the formula from (Duffett-Smith & Zwart 2017):

$$ \cos(\theta) = \sin(al_p) \sin(al_s) + \cos(al_p) \cos(al_s) \cos(az_p - az_s), \quad (1) $$

<footnote n="4">⁴www.celestrak.com</footnote>

<footnote n="5">⁵https://archive.org/web/</footnote>

<footnote n="6">⁶https://rhodesmill.org/skyfield/</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The angular distance to the telescope pointing direction (equation 1) for various satellites in the GAL constellation (solid lines represent). The black dashed line represents an angular distance of 5°.</figcaption></figure>

where az (al) refers to the azimuth (altitude) and the subscripts p and s represent the telescope pointing and satellite position, respectively. The telescope pointing can be provided from a scanning strategy or taken directly from the observation data files. The satellite position and telescope pointing are calculated for every time-stamp in the scan (two second resolution in this work).

In Fig. 2, we present the angular distances of various satellites within the GAL constellation for one of the observational scans derived from actual MeerKAT data. The seesaw structure is related to the fast scanning strategy adopted for this observation. As individual satellites move below the telescope horizon, they are disregarded. The appearance and disappearance of satellite information within the observation window illustrates this. We observe that a single satellite from the GAL constellation passes directly through the pointing centre of the telescope.

## 2.2 Beam response

To effectively simulate the observed signal from the satellite emissions present in the data, we require an accurate description of the MeerKAT beam response. An accurate beam model will inform us not only about the telescope’s main lobe but, more importantly, about the sidelobe structure and its frequency dependence. In Asad et al. (2021), the MeerKAT beam model was measured and fitted up to a radius of 5°. However, satellite emission beyond this radius can still affect our observations (HD18). Therefore, we decided to use a MeerKAT L-band beam model based on detailed electromagnetic simulations provided EMSS Antennas<fnref n="7" />, which offers the response of the MeerKAT antenna in both the HH (horizontal) and VV (vertical) polarizations. These simulations cover a frequency range of 900–1670 MHz and an angular range of 0°–100°. Beyond this, we set the beam response to zero.

Note that, we only consider the spherically averaged beam, so the dependence is solely on the angular distance from the beam centre (the polar angle, ‘θ’). In reality, beams are asymmetrical. The MeerKAT beam was measured to deviate from this symmetry by less than 1 per cent within the main lobe, although the deviation can become much larger at the second sidelobe and beyond (see Asad et al. 2021; de Villiers 2023). The provided EM simulations have an angular resolution of 0.1° and a frequency resolution of 2 MHz therefore interpolation was performed in both angle and frequency. The code also provides the holographic beam model from Asad et al. (2021), which we refer to as ‘Eidos’, and the analytical ‘Cosine beam’ which is a reasonable fit to the MeerKAT beam (Mauch et al. 2020). The top and middle panels of Fig. 3 shows a comparison between these models. We observe that beyond 10° the Cosine model starts deviating from the simulation.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** The different beam models available in the simulation. In the 1st and 2nd panels, we show the beam models: EMSS (blue thin dashed), Eidos (orange solid), and Cosine (green thick dashed) at 1150 MHz, with the 2nd panel up to higher angular separation (θ). The 3rd panel shows the EMSS beam at different frequencies: 1100 MHz (blue thin dashed), 1150 MHz (orange thick solid), 1250 MHz (thin green solid), and 1305 MHz (red thick dashed).</figcaption></figure>

In this paper, we only consider total intensity measurements. Although satellite emission can be polarized, modelling this is challenging, so we combine both polarizations in the simulation and the observed data (commonly referred to as horizontal, [HH], and vertical [VV], polarizations). The beam model is also normalized to have a maximum value of one. The bottom panel of Fig. 3 highlights the EMSS beam model as a function of the angular separation for different frequencies.

## 2.3 Power-spectrum density

The power-spectrum density (PSD) is a satellite’s transmitted signal spectral response. The signal along frequency can be described as a combination of sine and cosine waves (the modulation of the signal). The simplest of these modulations is a constant pulse known as binary phase shift keying (BPSK), which, when Fourier transformed, is described by:

$$P_{\text{BPSK}(\nu_c)}(\nu) = \frac{\text{sinc}(\pi(\nu - \nu_c)/[n_c f_0])}{\sqrt{n_c f_0}}, \quad (2)$$

where $f_0 = 1.023$ MHz is the reference frequency, $n_c$ is the chip rate, which refers to the rate at which a phase change occurs, and $\nu_c$ is the central frequency of the signal.

The signal’s PSD is then given by the autocorrelation (Teunissen & Montenbruck 2017a):

$$S_{xx}(\nu) = |P_{xx}(\nu)|^2, \quad (3)$$

where $xx$ is the modulation type and $S_{xx}(\nu)$ is normalized to integrate to 1 over frequency. Fig. 4 shows other, more complex signal modulation types used in the simulation (see Appendix A). Our current model can simulate BPSK; binary offset carrier (BOC); alternative BOC (altBOC), multiplexed BOC (MBOC), and time-MBOC (TMBOC). We list all the modulations used in this paper and the variables required for them in Table B2. These signals will be present at all frequencies, although out-of-band emissions will be smaller. It is possible that there is additional dampening of the signal away from the main peak, but such information is not available, so we assume these modulations are a good representation of the signal at

<footnote n="7">⁷https://www.emssantennas.com/</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Normalized PSD of transmission patterns. Solid lines represent odd values of *n* and dashed lines are even ones. Starting from the top panel, BPSK; sine BOC; cosine BOC; altBOC; MBOC and TMBOC. In Table B2, we highlight the various modulation types employed by the GNSS.</figcaption></figure>

all frequencies. Estimating whether the level of out-of-band emission is correct is one of the goals of this paper and will be addressed in a later section.

## 2.4 Signal temperature

Each satellite can emit a combination of different signals with specific modulations as described in the previous section. For each of these signals, *i*, the observed brightness temperature from one satellite, is

calculated as:

$$T_{\text{sat},i}(t, \nu) = B(\theta, \nu) S_{xx_i}(\nu) \frac{P_i}{r^2} \frac{c^2}{4\pi \nu^2 k_b}, \quad (4)$$

where *ν* is the observed frequency, *k*<sub>b</sub> is the Boltzmann constant, and *c* the speed of light in the vacuum. The beam response is given by *B*(*θ*, *ν*), where again, *θ*(*t*) represents the angular separation between the satellite position and the pointing centre of the telescope and is a function of time. The radial distance between satellites and the telescope on Earth *r*(*t*) is also a function of time and is calculated
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Raw visibility maps of antenna m000 for HH (right) and VV (left) polarizations for the observation scanning period as a function of frequency and time. The bright spots correspond to saturation in the receiver due to strong RFI, such as when satellites pass close to the pointing (see Fig. 2).</figcaption></figure>

with the SKYFIELD package (Rhodes 2019). Since the channel width for MeerKAT is very small, we neglect variations in $S_{xx}$ and simply use the value at the channel's central frequency, $\nu$.

$P_i$ represents the total emitting power for the signal $i$, which is defined using the transmitter antenna gain $G_{ti}$ (dBi) and the transmitted power $P_{ti}$ (dBW):

$$P_i = \frac{G_{ti} P_{ti}}{4\pi}. \qquad (5)$$

The spectral energy distribution (SED) describes the units of power per frequency. The SED for each signal from a satellite is given by the product of $P_i$ and $S_{xx}$. In some instances, the signal power values are not known to the public, and therefore, we selected a value of $P_i = 10$ dBW and $G_t = 10$ dBi where required.

Note that, we neglect the satellite beam in the modelling, for example, we assume the satellite is always pointing towards the telescope. This beam is expected to be quite wide and so any directional effect should be negligible, especially for satellites near the telescope zenith. The situation becomes more uncertain for satellites farther from the telescope's position, as they will be pointing downwards rather than towards the telescope.

The values used for each signal can be found in Table B2. The first column corresponds to the index $i$, used to identify the signal. Not all signals listed in the table are used in our simulation, as their peak emissions are well outside our target frequency band, and their contributions are expected to be negligible compared to the others.

For simplicity, we will assume that every satellite of a given constellation emits all the signals listed for that constellation (highlighted in grey in Table B2). This may not necessarily be true, as we know that some of the listed signals belong to specific generations of a constellation. As an example, younger satellites could have additional signals on board compared to their predecessors. However, such information is generally not publicly available. Therefore, in the simulation, we sum all satellite signals for a given constellation and constellations, as described at the beginning of Section 4. The terms $B(\theta, \nu)$ and $r$ are then calculated for each satellite individually. Note again that if the satellite is below the horizon at the telescope site or is more than 100° away from the pointing centre, we remove it from the calculation.

## 3 OBSERVED SATELLITE SIGNAL

To build and test our GNSS satellite simulation, we used the MeerKAT observational data described in W21. These data are a SD H1M pilot survey that mapped a region of around 200 deg² centred

<table><caption>**Table 1.** The five observational blocks used in the analysis. The block ID, antenna number, and date and time of the observation are listed. The results, unless stated otherwise, are shown for the highlighted observation.</caption><thead><tr><th>Observation ID</th><th>Antenna no.</th><th>Date/Time (UCT)</th></tr></thead><tbody><tr><td>1551055211</td><td>m000</td><td>2019-02-25 00:40:11</td></tr><tr><td>1 553 966 342</td><td>m000</td><td>2019-03-30 17:19:02</td></tr><tr><td>1 554 156 377</td><td>m000</td><td>2019-04-01 22:06:17</td></tr><tr><td>1 556 138 397</td><td>m000</td><td>2019-04-24 20:39:57</td></tr><tr><td>1 562 857 793</td><td>m004</td><td>2019-07-11 15:09:53</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** The time-average of Fig. 5, with the HH (blue curve) and VV (orange curve) polarizations shown for the m000 antenna. The light-grey (920–960 MHz) shaded region represents the low-end RFI, which stems from GSM and aeroplane transponders (Sihlangu et al. 2020). The green-shaded region (1140–1310 MHz) represents our frequency range of interest dominated by the navigational satellite RFI. The dark grey (1520–1630 MHz) shaded region relates to the higher frequency RFI in the L band, which contains signals from a few navigational services but is dominated by telecommunication services.</figcaption></figure>

at RA ~ 11 h and Dec. ~ 0° over the L-band, which ranges between ≈980–1700 MHz over a span of 4096 channels and a frequency resolution of 0.2 MHz (see Fig. 5). The observations are split into blocks of about 1.5 h, each covering the same area with all dishes pointing in the same direction. The data is taken in both linear polarizations, HH and VV. A standard SD flux/bandpass calibrator is observed before and after each block, and noise diodes are fired every 20 s in order to calibrate gain fluctuations in time. The observation blocks used in this analysis are listed in Table 1. In the next sections, we show the results for block 1551055211, dish m000 only, unless stated otherwise.

Fig. 6 shows the different frequency regions contaminated by RFI for one of the blocks. The light-grey shaded region (920–960 MHz) corresponds to the Global System for Mobiles (GSM), aeroplane transponders, amongst others (Sihlangu, Oozeer & Bassett 2020). The dark grey shaded region at higher frequencies (1520–1630 MHz)
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Functional form of the gain as a function of time (left panels) and frequency (right panels) for m000 antenna W21. The top panels are for HH polarization and the bottom ones for VV. RFI flagging results in the missing information in frequency.</figcaption></figure>

corresponds to telecommunication services and a few navigational satellites. The region shaded in green corresponds to our frequencies of interest (1140–1310 MHz) which navigational satellites’ RFI dominates. For reference, the target cosmology window used to make a cross-correlation detection of the power spectrum with the same data (Cunnington et al. 2023a) is at 973–1015 MHz. As mentioned before, we ignore the navigational satellites present in the high L-band frequency range, since their contribution should be subdominant at 1140–1310 MHz compared to the satellites emitting there (and even less in the target cosmology band).

## 3.1 Calibration of the RFI-contaminated region

In principle, it would be ideal to use the gain solutions derived in W21. From these, we could calibrate the time-ordered data (TOD) for each block, dish and polarization using the following approach:

$$T_{\text{cal}}(t, \nu) = T_{\text{raw}}(t, \nu)/g(t, \nu), \qquad (6)$$

where $T_{\text{cal}}$ is the calibrated temperature in kelvin, $T_{\text{raw}}(t, \nu)$ the raw visibility in correlator units, and $g(t, \nu)$ is the gain. Unfortunately, the calibration pipeline used in W21 immediately flagged the 1100–1350 MHz region for all times due to the strong RFI. This poses a problem for us since we need the data to be calibrated within that frequency interval in order to compare it to our simulations. In Fig. 7, we show the derived gain from W21 for one block and one dish as well as the HH and VV polarizations (top and bottom, respectively). The left panels show the gain as a function of time (averaged over frequency), while the right panels show the time average as a function of frequency (only non-flagged times and frequencies are included in these averages).

Using a point source to calibrate the flagged region is quite challenging due to the very strong RFI signal. Noise diodes could be better, however, their frequency dependence is not very well known from the lab and again, calibrating the noise diodes beforehand with a point source would have the same contamination issues in the target RFI contaminated region. Instead, we took a different approach. We first assume that the gain can be broken into the temporal ($g_t$) and frequency ($g_\nu$) components, for example, $g(t, \nu) \approx g_t(t)g_\nu(\nu)$. This essentially means that the bandpass shape is assumed constant throughout the observation, with only an overall amplitude variation over time. Such an assumption is well supported by the analysis in W21, with frequency-dependent variations measured to be less than

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Estimate of the base raw visibility without satellites after the 4th iterative clipping (dashed line) overlaid on antenna m000 minimum raw visibility (solid line) for the HH (top) and VV (bottom) polarizations.</figcaption></figure>

1 per cent. The gain can then be written as:

$$\tilde{g}(t, \nu) = \frac{\langle g(t, \nu) \rangle_{\nu} \langle g(t, \nu) \rangle_{t}}{\langle g(t, \nu) \rangle_{t, \nu}}, \qquad (7)$$

where $\tilde{g}$ is used to represent the new derived gain. The averages, $\langle \cdot \rangle$, are over the non-flagged times and frequencies. This immediately provides us with the time variation of the gain, as shown in Fig. 7.

Obtaining $g_\nu$ is more complicated due to the missing information. A simple solution would be to interpolate the frequency curve in Fig. 7, but the gap is too wide. To obtain more points for the interpolation, we assumed that the frequency dependence of the gain follows the frequency dependence of the raw visibilities without the satellites. This implies that the input temperature should be flat across the frequency gap, except for the satellites. In reality, the ground pickup, receiver temperature, and galactic synchrotron have a slope. However, this effect will be small within the interval considered and can be corrected later when matching to the gain solution outside the RFI band. Below, we describe in detail the steps we followed to complete the missing frequency information in the gain maps:

(i) We start with the raw visibility of both HH and VV polarizations before RFI flagging. We then select the minimum visibility values across time for each frequency channel, which we take as an estimate of the channel’s floor level when satellites are not present. We do this in the frequency range spanning 1000–1500 MHz, which is wider than the RFI-contaminated band, but it allows us to obtain the tail-end information for the frequency region we wish to calibrate;

(ii) The next step is to remove all the peaks in the data caused by the satellite contamination, to get a smooth curve of the floor level. We start by performing a spline fit to the minimum raw visibility, which provides a smooth curve in frequency. We then calculate the difference between the spline and minimum raw visibility curves. From this residual, a sigma clip (Astropy Collaboration et al. 2013) is applied to remove the peaks that stem from the satellite emission. This procedure is interactively applied to the clipped curve until no further clipping is needed (after five iterations). Finally, we obtain a smooth overlay spanning the affected frequency channels as seen in Fig. 8;

(iii) The resulting smooth curve represents the underlying base’s structure of the raw visibility and we assume that the gain will be proportional to this. We then normalize the smooth curve by matching it
<figure><img src="image_10.png" /><figcaption>**Figure 9.** The bandpass of the gain for the m000 antenna in the HH (top panel) and VV (bottom panel) polarization. The blue curve represents the available gain information from W21; the orange dashed curve represents the best estimate of the gain across the frequency region of interest.</figcaption></figure>

to the known gain from W21 outside the satellite's frequency region. We did this between 1138 and 1307 MHz for the HH polarization and between 1100 and 1314 MHz for the VV polarization. Finally, we create a smooth function (spline) that becomes the final frequency component of the gain map, as seen in Fig. 9.

The final temperature map from the data will be the average of the calibrated temperatures, $T_{cal}(t, v)$, from each polarization.

## 3.2 Background model

In order to compare the satellite simulation to the data, we also need to account for the other temperature contributions, which we call the 'background model'. This background is much smaller than the satellite emission and reasonably flat across frequency, forming the baseline level we observe in Fig. 6. We follow W21 to model this background temperature ($T_{BG}$) as

$$T_{BG}(t, v) = T_{\text{rec}}(t, v) + T_{\text{el}}(t, v) + T_{\text{gal}}(t, v) + T_{\text{CMB}}, \quad (8)$$

where $T_{rec}$ is the receiver temperature, $T_{el}$ depends on the elevation and includes contributions from the ground spill and atmosphere, $T_{gal}$ is the synchrotron emission from the centre of our Galaxy and $T_{CMB}$ is cosmic microwave background (CMB) temperature.

$T_{CMB}$ is constant, and since the observations were done at constant elevation, $T_{el}$ should be constant in time. The same is expected from $T_{rec}$ since the instrument is quite stable over the course of the 1.5 h observation. $T_{gal}$ will change as the telescope pointing moves through the sky, but these fluctuations will be at ~0.1 K level. In principle, we could simply remove most of this background contribution by subtracting a constant in time. We would then have to do the same to the satellite simulation for a proper comparison. This would mean we would not be able to fit for the satellites DC level which can be useful for other situations. So in this paper, we did not remove any time offset and instead included $T_{BG}$ in the simulation. For this, we use the same models as in W21. However, since the receiver temperature was calculated as part of the calibration process in W21, the frequencies corresponding to satellite emission are also flagged (Fig. 10).

To restore the masked information, we will follow the same approach as the one used to estimate the gain in the RFI-contaminated band. We assume that the receiver temperature is a separable function in time and frequency, which is then proportional to the frequency and time average of $T_{rec}(t, v)$, respectively. These averages do not use any of the flagged data, and this is done for each block, dish and polarization. Again, the function in time is already continuous. For the frequency dependence, we interpolate the known data points (solid blue

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Receiver temperature for the m000 antenna in HH (left panel) and VV (right panel) polarization. RFI cleaning results in missing channel frequency information from 1075 to 1305 MHz.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** A time-averaged representation of the receiver temperature along frequency for the m000 antenna with HH (left) and VV (right) polarization. The blue curve is the receiver temperature from the model, whilst the orange dash curve is the RBF interpolation overlaid. In the lower and higher frequencies, the RBF follows the data and makes the best attempt to connect the high and low-frequency values.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Interpolated receiver temperature for the m000 antenna in the HH (left panel) and VV (right panel) polarization.</figcaption></figure>

in Fig. 11) using a radial basis function (RBF) Virtanen et al. (2020). This results in a curve (orange line) covering the whole frequency range (Fig. 11). In this instance, we do not have any data points in the RFI-flagged region, which makes the interpolation more uncertain. However, looking at the fluctuations in Fig. 11, one should expect any deviations to be less than 0.5 K. Since this is an additive quantity to the far stronger satellite signal, such errors should have little impact in the analysis. We plot the final receiver temperature used in Fig. 12.

## 4 CALIBRATING THE SATELLITE SIMULATION WITH OBSERVATIONS

In this section, we aim to determine how well the simulation performs in recreating the navigational satellite signals we see in the observational data of the MeerKAT telescope. We compare the simulated satellite signal $T_{sat}(t, v)$ to the calibrated observations from each block. In this regard, we have to combine the background model of Section 3.2 to that of $T_{sat}(t, v)$. This allows us to obtain the same temperature base component between the observation and simulation. Our results, which are centred around 1100–1350 MHz, extend toward the tail of the spectrum, where the background model dominates.

The signal modelling presents a few uncertainties, namely, in the overall amplitude of the temperature of a signal from a constellation. The emitting power and gain of the satellite are not constant with
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Comparison for the no-mask scenario. The 1st row shows the waterfall images for the observation (left) and the simulation: case 1 (centre, hereafter $C_1$) and case 2 (right, hereafter $C_2$). The temperature colour bar was set at the maximum and minimum values corresponding to the observation. The 2nd row shows the time-averaged frequency plane of the waterfall plots, with the left panel being $C_1$ and the right being $C_2$. The solid blue curve represents the observation, and the dashed orange curve represents the adjusted simulation for the two cases.</figcaption></figure>

time, and we expect it to be reduced as the satellite ages. Moreover, in some cases, the values tabulated in the literature appear as upper limits or are not available. Uncertainties in the telescope's primary beam or calibration errors can also affect the results. Because of these issues, we introduce a set of $\alpha_i$ parameters associated with the overall amplitudes of each signal within a given constellation:

$$T_{S,i}(t, \nu) = \alpha_{sat,i} T_{sat,i}(t, \nu), \quad (9)$$

where $T_{S,i}$ is the simulated temperature of the i-th signal and $T_{sat,i}$ is given by equation (4). Each signal in Table B2 has an associated $\alpha_{sat,i}$. These just mimic the overall amplitudes of a given spectral energy distribution. While they are assumed to be fixed in frequency, their value is usually only valid for short time intervals due to changes in the satellites that are up in the sky at any given time. These amplitudes should be seen as global calibration parameters and not intrinsic to the satellite. A full characterization is beyond the scope of this paper and requires a more elaborate model with priors on it. This will be investigated in a follow-up work. If there is a perfect match between the satellite models and the observations, these $\alpha$s should equal one. The number of parameters introduced depends on the number of relevant signals that fall within our frequency range of interest and which constellations are visible to the telescope. The total simulated signal is then:

$$T_S(t, \nu) = \sum_{\text{constellation}} \sum_{j} \alpha_{sat,i} T_{sat,i}^j(t, \nu) + T_{BG}(t, \nu), \quad (10)$$

where we sum over all constellations, i runs over the signals present in each constellation and j runs over the satellites in that constellation.

## 4.1 Estimating the $\alpha_i$ amplitudes

In order to find the value of $\alpha_i$ that best approximates the simulation to a selection of the data, we minimize a cost function (CF) which we define as

$$CF = \frac{1}{N} \sum \frac{(T_O - T_S)^2}{\sigma_D^2}, \quad (11)$$

where $T_O$ is the adjusted calibrated observational data and $T_S$ is the simulated temperature response. N is the total number of points used

and is only used for normalization when we compare cases with a different number of data points. The sum is over all times and frequencies of the TOD. $\sigma_D^2$ can be seen as a weight for each data point, and we consider two distinct cases:

$$\sigma_D = \begin{cases} T_O [\text{case 1} (C_1)] \\ \text{or} \\ 1K [\text{case 2} (C_2)] \end{cases}$$

Case 1, ($C_1$ hereafter) uses the temperature itself, which is proportional to the expected thermal noise as in the radiometer equation. Essentially, this would correspond to the standard $\chi^2$ up to a proportionality constant. This can also be seen as minimizing the percentage difference between the model and the data. The problem with this cost function is that the uncertainties in the model should dominate the expected thermal noise, since at 2 s and 0.2 MHz resolution, the thermal noise gives an error of about 0.1 per cent for each data point. This case down-weights the data in the region where the satellite's signal is the brightest so that the simulation calibration will be dominated by the (many) small temperature values in the region where the background model dominates, making it harder to use the rather important bright satellite peaks. For this reason, we also consider case 2 (denoted as $C_2$ below), where all the data points are equally weighted, for example, $\sigma_D = 1$.

The minimization of the cost function is accomplished with the SCIPY.OPTIMIZE.MINIMIZE package (Virtanen et al. 2020) using the *Powell* method. The initial conditions for the $\alpha$ parameters are set to zero for the first iteration; a lower bound of zero and an unconstrained upper bound are placed on the parameters as boundary conditions. A tolerance threshold (tol) of $10^{-6}$ is applied to the minimization as a convergence limit. If a signal is not used in a given observation (because the corresponding satellites are too far away from the pointing direction) then the $\alpha$ will remain at zero.

The results for the calibration of the simulation using block 1551 055 211 are shown in Fig. 13. Although some of the main features are present, we observe that the amplitude of the simulated signal is smaller and misses some of the peaks. This is more evident in the bottom panel, where a time-averaged waterfall plot is shown. This is true for both cost functions, although $C_2$ performs better.
The solid blue represents the observational data, while the dashed orange line represents the simulation calibrated result for $C_1$ (left) and $C_2$ (right), with a minimum value of the cost function of 0.0240 and 802.91, respectively. The reason for this is that the observation saturates when the satellites get too close to the telescope, which makes our model break down. In this block, two satellites passed through the telescope's pointing. The first belongs to the GLO constellation and breached the $\theta \le 1^\circ$ between 1622.39–1672.37 s. The second satellite belongs to the GAL constellation (see Fig. 2), and was within the $\theta \le 1^\circ$ between 3913.426–4195.30 s. This saturation is not clear in the bottom panel of Fig. 13 because we averaged over time and only a small percentage of the time is saturated. This is better seen in Fig. 14. The top panel shows the signal across frequency for the snapshot at the time when the signal is stronger (satellites close to the pointing). Note that, the saturation/cut-off level is not constant across frequency. The two middle panels show the total signal along time for the peak frequencies belonging to the GLO constellation (left panel) at 1248.3 MHz, and the GAL constellation (right panel) at 1278.75 MHz. Although other signals (and satellites) may contribute to these frequencies, they are the dominant ones. The fluctuations in time should depend on the angular distance to the pointing direction and follow the primary beam shape. The solid orange line represents the result for the calibrated simulation with $C_2$, while the dots represent the observational data, with different colours corresponding to different angular distances. Below $2^\circ$ the telescope system completely saturates at these frequencies.

The bottom panel of Fig. 14 shows the calibrated temperature from the observed block at the frequency of 1278.6 MHz against the expected movement of the GAL satellite GSAT 02029-(PRN-0E09) (EUSPA 2024). The position is calculated exactly for the observation date and time-stamps and the angular distance is derived with respect to the dish pointing centre, taking into account the movement of the dishes as well. Although other signals will ‘contaminate’ this frequency, we can clearly see the shape of the MeerKAT primary beam as expected. Aside from contamination by other satellites (at other angular distances), there are other reasons why the comparison to the primary beam is not perfect, in particular, the saturation of the signal and small errors in the satellite position. These errors have a greater impact as the satellite approaches the dish pointing. Finally, we note that the saturation not only imposes a cutoff on the observed signal but also makes the gain non-linear. This non-linearity can affect most frequencies when the signal becomes strong, even outside the peak frequency, and can show up before saturation. However, our model was built assuming the telescope response was linear, even in the presence of navigation satellite RFI. As we can see, this is not always true and will affect our calibration with the data. This raises the question of addressing the saturation issue by masking pixels/timestamps where the telescope response becomes non-linear. In other words, can we improve the model that recreates the satellite’s emissions? These possibilities are explored in the next section.

## 4.2 Testing limits of validity

From now on, we will only consider $C_2$ to test the limits of validity in our approach. In all masking scenarios, we have compared $C_1$ with $C_2$ and concluded that $C_2$ gives better overall results.

### 4.2.1 Angular masking

To address the saturation/non-linearity of the data which affects our determination of the signal's renormalization amplitudes, we

<figure><img src="image_15.png" /><figcaption>**Figure 14.** A snapshot of the data along different axes: frequency (1st panel), time (2nd and 3rd panels) and along the angular separation of a GAL satellite: GSAT0209-(PRN-E09) (4th panel). In the 1st panel, the data (blue), $C_1$ (orange), and $C_2$ (green) at two different timestamps correspond to incursions from the GLO (left panel) and GAL (right panel) constellations. The snapshot is when $\theta \le 1^\circ$ for both constellations; it also reveals the peak frequency of GLO at ~1248 MHz. However, in the case of the GAL constellation, the peak is hidden due to the saturation of the receiver. In the 2nd and 3rd panels, based on the information in the 1st panel, we select the peak frequencies from Table B2. The 2nd panel is at 1248.3 MHz (GLO), and the 3rd panel at 1278.6 MHz (GAL). The observational data are shown with blue dots, the $C_2$ is shown in orange, and the limits of incursion by the satellites are shown for the 1°, 2°, 5° in red, green, and magenta, respectively. In the 4th panel, we show the observed temperature versus the angular separation values for the GAL satellites at the frequency selected around which the satellite operates.</figcaption></figure>

start by trying to remove any points that correspond to satellites close to the observation direction. Using a simulation with the same observation specifications, we calculate the positions of all the satellites as shown in Fig. 2, and remove time-stamps corresponding to an angular distance below a certain threshold. This is done on both the observation data and simulations. We used an angular cut of $1^\circ$ and $5^\circ$. The results can be seen in Fig. 15. A satellite with $\theta \le 1^\circ$ would be found inside the main lobe of the MeerKAT beam and saturate the signal, while $\theta \le 5^\circ$ should also remove the contribution of satellites that can cause non-linearity effects on the receivers without necessarily saturating the signal. We can see that indeed the minimization of the CF performs much better for the cut with
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Comparison for the *angular masking scenario*. The waterfall images (1st and 2nd rows) for the various angular are displayed: the observation (left) and the simulation (right). The 1st row shows the 1°, and the 2nd row shows the 5° angular masking cut, respectively. In the 3rd and 4th rows are the time-averaged results as a function of frequency for the two angular cuts. The observation is in solid blue curve and the simulation in dashed orange.</figcaption></figure>

5°. On the other hand, the 5° cut means that any satellite with a weaker signal will probably be below the error and would not be accounted for. Moreover, this cut is quite aggressive, removing up to 48 per cent of the data while the 1° will remove about 6 per cent. Note that, this implies a loss of the main peaks in the data. The effect at lower frequencies where we do our cosmological analysis is very hard to see and does not impact the results.

### 4.2.2 Thermal masking

This approach tries to address the issue of saturation directly, by making cuts when the temperature goes above a certain threshold. This is applied to the observation data both across time and frequency and the mask is then transferred to the simulation. The threshold values for the thermal cut are 25, 50, and 100 K, corresponding to 19 per cent, 2 per cent, and 1 per cent of data loss, respectively. The results are shown in Fig. 16. The fitting improves as more of the saturated pixels are removed. However, this is still not as good as the 5° scenario, even for the 25 K cut, as can be seen by the frequency plots or the minimum values of the cost function. This is more obvious for the emission above 1250 MHz, which is a contribution from the GAL E6 and BEI B1-2 and B3 signals. Additionally, on the lower end of the frequency scale, around 1150 MHz, the simulation does not account for the smaller peaks, in particular the second side lobe emanating from the central peak at ~1175 MHz. It is also interesting to note that the thermal cut of 25 K limits the peak

<figure><img src="image_17.png" /><figcaption>**Figure 16.** Comparison for the *thermal masking scenario*. The waterfall images (1st, 2nd, and 3rd rows) with the observation (left), and the simulation (right). The time-averaged results below in the 4th, 5th, and 6th rows, with the observation (solid blue curve) and simulation (dashed orange curve). The various thermal thresholds given are: 100 K (1st rows), 50 K (2nd rows), and 25 K (3rd rows).</figcaption></figure>

contributions from satellite signals, removing more signal there than the 5° cut.

### 4.2.3 Full thermal masking

In the previous scenario, only the 'pixels' in time and frequency that exceeded a given threshold were flagged. We now consider a case where all the points across frequency are flagged for a given time as long as one of the values exceeds that cut-off. The reasoning is that once the receiver system saturates, the gain will become non-linear, affecting all frequencies, not just those where the signal hits the limit. Since the saturation limit seems to change between blocks, we decided to use a percentage of the maximum pixel value from each block instead of the threshold. We applied masking thresholds of: $T_{max}/2$, $T_{max}/5$, and $T_{max}/7$ (results are shown in Fig. 17). For this block (ending 211) and antenna m000, the maximum temperature $T_{max}$ is ≈ 400 K, which correlates to temperature thresholds of ≈ 200, ≈ 80, and ≈ 60 K, respectively. As seen from the time average plots in the bottom panels, one can capture the signals in the data quite well, even for the higher-level thresholds. In particular, they are clearly better than the equivalent
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Comparison for the *full thermal masking scenario*. The waterfall images (1st, 2nd and 3rd rows) with the observational (left), and the simulation (right). The time-averaged results below in the 4th, 5th, and 6th rows, with the observation (solid blue curve) and simulation (dashed orange curve). The various pixel thermal limits are set at: 50.00 per cent (1st and 4th rows), 20.00 per cent (2nd and 5th rows), and 14.29 per cent (3rd and 6th rows) of the maximum temperature.</figcaption></figure>

temperature cut in the previous scenario. This shows that there is indeed unmodelled contamination away from the peak frequencies, probably due to non-linearities, which are removed by flagging the full-time dump. At the same time, the level of data loss is smaller than, say, the 5°, with 7 per cent, 14 per cent, and 20 per cent for each threshold.

Table 2 summarizes the results from the different data cuts. There is a clear trend of improvement as we remove more of the strong peaks, which should remove the saturated data points in the system. Other things could contribute to this improvement, such as reducing the out-of-band signal which is hard to model. In the next sections, we discuss possible issues in our model that could affect this calibration, besides the expected saturation/non-linearity.

### 4.2.4 Time variation of the α parameters

The α parameters are used to renormalize the satellite total power amplitude. In that sense, they should have little to no time evolution since the satellite signal should be stable in time. However, three situations can change this. First, uncertainties in the telescope's

<table><caption>**Table 2.** The minimum cost function ($C_2$) values for various masking scenarios and the % of data that is flagged per scenario.</caption><thead><tr><th>Masking Type</th><th>C<sub>2</sub></th><th>% flagged data</th></tr></thead><tbody><tr><td></td><td><strong>Angular</strong></td><td></td></tr><tr><td>1°</td><td>261.24</td><td>5.4</td></tr><tr><td>5°</td><td>1.21</td><td>47.55</td></tr><tr><td></td><td><strong>Thermal</strong></td><td></td></tr><tr><td>100 K</td><td>24.47</td><td>1.39</td></tr><tr><td>50 K</td><td>7.99</td><td>2.25</td></tr><tr><td>25 K</td><td>1.41</td><td>9.23</td></tr><tr><td></td><td><strong>Full thermal</strong></td><td></td></tr><tr><td>T<sub>max</sub>/2</td><td>5.20</td><td>7.49</td></tr><tr><td>T<sub>max</sub>/5</td><td>2.56</td><td>14.25</td></tr><tr><td>T<sub>max</sub>/7</td><td>1.80</td><td>20.1</td></tr></tbody></table>

primary beam, such as asymmetries, could be absorbed into the αs, which can then change depending on the position of the satellite with respect to the telescope pointing. Secondly, as the observed signal is always a mix of different satellites, it would not be a surprise if different values of the same αs provide a better fit depending on what satellites appear in the FoV. However, such variations should be small when one single signal dominates the data. A third reason for time evolution (and difficulties in the simulation calibration) is that individual satellites from a given constellation might not be emitting the exact same signals as we assume. The α amplitudes will try to correct this and might change as these satellites from the same constellation go in and out of the FoV. For these reasons, a value of αᵢ is only valid for a small time interval. Furthermore, since each pixel is 2s×0.2 MHz, the uncertainties in the model should dominate over the thermal noise, rendering any tentative error bars on overall αᵢs uninformative. The errors in the fitted αᵢs are dominated by uncertainties in the satellite model and calibration errors in the data, not thermal noise. For a proper parameter estimation and characterization of the GNSS, a better model is required. The development of such a model, and fit with prior information, or jointly with instrumental parameters (Finlay et al. 2023), is beyond the scope of this paper and will be left for follow-up work.

In Fig. 18, we show the observation and calibrated simulation for three different time windows, corresponding to 775–1000 s, 2200–2400 s, and 5500–6200 s. These time periods were picked so that no satellites entered the 5° zone to avoid saturation/non-linearities. By looking at the minimum values of $C_2$, we can again confirm that variations in the fit across time are not due to the telescope thermal noise fluctuations, as these would lead to sub-per cent fluctuations in $C_2$, given the large number of data points being used. Fig. 19 shows the values of the amplitudes we found α parameters. The figure shows a non-negligible variation in some of the αs. In particular, $i = 5, 7, 9, 12, 14$, which is probably due to satellites from GAL and GLO entering the FoV. We note again that since the starting values for the αs is set to zero, an α moving from zero simply means that the satellites with the corresponding signal have moved into the FoV. However, once an α becomes non-zero, it should keep that value, if not for the abovementioned reasons.

A better way to see how well the calibration of the simulation approximates the data is through the residuals, as in Fig. 20. A reason for the observed difference could be that the frequency shapes used in the modelling (taken from literature, see Table B2) do not accurately follow the real satellite signal. When taking the time average of the residuals, the narrow spikes around ≈1225 MHz
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 18.** Comparison for different time windows. The waterfall images (1st, 2nd, and 3rd rows) with the observation (left), and C₂ (right). The time-averaged results below in the 4th, 5th, and 6th rows, with the observation (solid blue curve) and simulation (dashed orange curve). The various temporal periods are set at 775–1000 s (1st and 4th rows), 2200–2400 s (2nd and 5th rows), and 5500–6200 s (3rd and 6th rows).</figcaption></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 20.** The residual 2d and time-averaged residual 1d for the temporal masking scenario (Fig. 18). The various temporal periods are set at 775–1000 s (1st and 4th rows), 2200–2400 s (2nd and 5th rows), and 5500–6200 s (3rd and 6th rows).</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 19.** The α amplitudes for each temporal window (Fig. 18). The different periods are 775–1000 (blue ●), 2200–2400 (orange ★) and 5500–6200 (green x). The central frequency of each α signal is shown on the top y-axis. For the constellations: GPS-violet (α₁ – α₄), GLO-blue (α₅ – α₉), GAL-dark red (α₁₀ – α₁₃), BEI-teal (α₁₄ – α₁₈), IRNSS-black (α₁₉ – α₂₀) and SBAS-dark grey (α₂₁).</figcaption></figure>

show the width of the peak is not captured perfectly in the model; this peak is associated with the GPS-L2 signal. Similarly, peaks at 1246 MHz and 1248 MHz, belonging to Glonass-L2 signals, are not accurately accounted for by the simulation throughout the different time periods that are seen in the observational data. Another example is the dip at 1305 MHz, indicating that the model has a peak that is not seen within the observational data. A possible explanation is that the out-of-band transmissions are being dampened; such dampening is not accounted for in our current modelling of the satellite signals.

### 4.2.5 Positional errors

One last issue we will address with the simulation calibration is the possibility of positional errors. As the dishes move back and forth, an error in the satellite position will translate to an amplitude difference that changes with time. To test this, we isolate the movement of a single GAL satellite that approaches the telescope’s pointing position and focus on 1278.75 MHz, which corresponds to one of the peak frequency emissions from the GAL constellation. The immediate neighbouring channels to the peak frequency are averaged to increase the signal-to-noise.
<figure><img src="image_24.png" /><figcaption>**Figure 21.** Ratio between the observational data and the model for the temporal masking scenario. Isolating a single GAL satellite operating at a frequency range (averaged) of 1278.6–1278.8 MHz as it moves with respect to the MeerKAT pointing. The 1st, 2nd and 3rd rows represent the periods of 775–1000, 2200–2400 and 5500–6200 s, respectively. The left column is when no time average is applied, and the right is for a 10-s time average. The selected frequency is displayed with the CF for each panel.</figcaption></figure>

For this test, we use the data mentioned above with temporal windows that do not include any satellites below 5°. Including data below the 5° cutoff would make the results more obvious; however, it would be affected by saturation and non-linearities. In Fig. 21, we show the signal ratio between the observation and the simulation. As the dishes and the satellite move, the same angular position of the satellite will be observed at different times. Therefore, we stack any time corresponding to the same angular position of the target satellite to increase the signal-to-noise. The panels on the left reveal a ringing effect that oscillates around one, essentially following the shape of the telescope’s primary beam. The obvious implication is that there is an error in the angular position which translates to an error in the beam value we should use. There could also be smaller contributions from other satellites since we cannot fully isolate the signal. Those would further mix up the ratio as they will be at different positions for the same times considered here.

Even in the far sidelobes, the beam can change significantly in relative terms, with ratios of up to 100 between the minimum and maximum values as the beam oscillates, with a period of about 1 degree. The dishes are moving in azimuth at speeds of 5 arcmin s⁻¹, for example they would move 10 arcmin within the time resolution of 2-s, which can lead to a mismatch even if the satellite is not moving. We therefore consider smoothing both the data and simulation over 10-seconds. The idea is that any errors in the position as the dishes move back and forth could be averaged out. The results are shown in Fig. 21. The oscillating structure seems to disappear, but part of this is simply due to the fact that we now have fewer points. Overall, the ratio does get closer to 1, which improves the fitting, though not substantially.

In summary, the satellite central frequencies are well represented in the observational data. The frequency structure of the different global signals in particular *GPS*, *GAL*, and *BEI* were well described even with the overlapping contribution from *GAL* and *BEI* at 1268 MHz. One concern arises from the *GLO* whose emitted power and gain are not well documented in literature. Additionally, *GLO* at 1248 MHz had several signals with similar frequency structures which can cause difficulties in finding the amplitude's values given the uncertainties

in the tabulated specs. Apart from the *GLO* issue, the models appear to have a good representation of the data, as long as saturation is kept to a minimum. The other challenge is whether satellite signals weaken faster the further we are from the central frequency.

## 5 IMPACT ON THE COSMOLOGICAL BAND

Thus far, we have considered the impact of GNSS satellites around frequencies where they dominate. This allows for higher signal-to-noise when calibrating the satellite simulations. These simulations can, in principle, be used to plan for observations that avoid strong contamination and also provide ways to flag the signal better. For cosmology, we have focused on a frequency band that we hope is entirely free from satellite contamination, ranging from 973 to 1015 MHz, where Cunnington et al. (2023a) presented a measurement of the cross-correlation power spectrum. To further help clean the data, the flagging algorithms will usually completely remove time stamps where the signal saturates, even if the contamination is low at the target frequencies. Fig. 22 shows a comparison between the observation data and the simulation for two frequency channels within this 'cosmology band'. The same flags were applied to data and simulation, taken from W21. The αs in the simulation were first found using data from the frequencies where we clearly see satellite peaks (1100–1350 MHz).

Furthermore, to remove any background contamination, we have subtracted the time average for each channel, so that the signal will fluctuate around zero. There is no obvious pattern in the data, which seems noise-dominated, and most of the satellite signal predicted by the simulations is within that noise. This makes it hard to try to fit for any out-of-band contamination from the satellites. However, it is possible that when we combine multiple observations and calculate the H1 power spectrum, some out-of-band signals from the satellites may still contaminate the cosmological signal. Indeed, in Cunnington et al. (2023a), we needed 30 principle component analysis (PCA) modes to clean the signal and achieve a detection. It is unclear what are the main reasons for such contamination. A significant portion is likely due to issues in the bandpass calibration or other sources of RFI. However, there may still be some residual satellite
<figure><img src="image_25.png" /><figcaption>**Figure 22.** A temperature comparison across the observational scanning period between the satellite model (blue) and the residual calibrated observational data (W21) in the HH (orange) and VV (green) polarization for the m000 antenna. Masks from W21 were included. The time average per channel is subtracted from each frequency channel. The frequency channels, 973.0 MHz (left panel) and 988.7 MHz (right panel), are selected as they reside in the cosmological band (973–1015 MHz) of the observation.</figcaption></figure>

contamination, although we do not see any obvious signature in the data. In the next section, we use the calibrated simulations to predict how much satellite contamination is expected in the cosmology band. We note that there is a fair amount of uncertainty in the out-of-band emission due to the lack of information in the literature.

## 5.1 Extracting the H I power spectrum in the presence of satellites

Our goal in this section is to test the impact of satellites on MeerKAT's ability to measure the H I power spectrum. To achieve this, we will combine the satellite simulations with simulations of the H I signal and foregrounds. This will then be processed in the same way as we would for the data, with foreground cleaning followed by a power spectrum estimator. The procedure is similar to what is described in Irfan et al. (2024), where we go from TOD to data cubes and then power spectrum. To make the conclusions clearer, we will neglect instrumental noise, as this would tend to dilute any residual systematics.

We begin by applying the same methodology as shown above to four additional observational blocks. The inclusion of additional blocks (see Table 1) aims to improve the sky coverage and include a good mix of satellites that are expected to show up on different days. We calibrated the simulation to each block and tried different masking scenarios. In Fig. 23, we showcase three blocks (two new additional blocks) using the best performing full thermal masking scenario with a thermal threshold of 33 per cent of the maximum temperature to remove any saturated points. We decided to use these models for further processing. The full simulated data blocks were then flagged. At this stage, we opted to remove any data points with satellites closer than 5°. However, one of the blocks contains a GNSS satellite in geostationary orbit that is often below the 5° and would cause a large percentage of data to be masked. For this block, we applied a 1° flag instead. Since satellite movement differs across observations, the masks of the blocks do not overlap, and there are no gaps in the combined maps.

For the H I simulation, we used Gaussian density fields generated from a CAMB power spectrum (Lewis & Bridle 2002) and produced H I maps at the sky positions we are observing. The Galactic emission

<figure><img src="image_26.png" /><figcaption>**Figure 23.** Results for the full thermal masking scenario with a threshold at a third of the maximum temperature for three observation blocks.</figcaption></figure>

component was generated under the assumption that the dominant contribution is from synchrotron emission, and that synchrotron emission follows a power-law parametrization:

$$T_{\text{sync}}(v, p) = T_{\text{amp}}(v, p) \left(\frac{v}{408}\right)^{\beta_{\text{sy}}(p)}, \quad (12)$$

where the amplitude template [$T_{\text{amp}}(v, p)$] at 408 MHz is the reproduced, destriped Haslam 408 MHz map (Remazeilles et al. 2015) with the combined CMB monopole and unresolved point source background constant offset removed (Wehus et al. 2017). The spectral index template [$\beta_{\text{sy}}(p)$] is the model for the spectral index map from Miville-Deschênes et al. (2008).

Both the H I and Galactic foreground model required smoothing by the beam; whilst the EMSS beam is the optimal choice for modelling the satellite contributions, which are still of considerable magnitude up to several degrees away from the main lobe, a simple Gaussian beam model will suffice for smoothing the Galactic and H I contributions as the aim of these simulations was to assess the impact of the satellite residuals on top of a smoothed sky signal (Matchawule
<figure><img src="image_27.png" /><figcaption>**Figure 24.** A comparison between the synchrotron emission only (magenta solid line) and the combined satellite and synchrotron emission (cyan dashed line) for the observational block of 155 105 211 after a 5° angular cut is applied.).</figcaption></figure>

et al. 2020). The frequency-dependant Gaussian beam FWHM used was:

$$ \theta_{\text{FWHM}}(\nu) = 1.2^{\circ} \times \frac{1280 \text{ MHz}}{\nu}. \qquad (13) $$

Finally, following once again Irfan et al. (2024), we use the RA and Dec. information from each data block to extract the corresponding foreground and H1 signal from the HEALPIX (Górski et al. 2005) maps using tools from the HEALPY library. This way, we generate equivalent TODs for the Galactic synchrotron and H1.

For the five observational blocks, the residual satellite contributions are not visible by eye above the Galactic contributions within the TOD. Fig. 24 illustrates this by showing the frequency-averaged TOD for block 1551055211, both with and without the satellite contribution included in the simulation.

We now have five observation blocks (TOD), each including satellite signal, H1 and galactic synchrotron, convolved by the appropriate beam. As mentioned above, we do not include thermal noise so that we can test the results in the optimal case of high signal-to-noise, when small systematic residuals will become visible. A few ingredients are still missing: the receiver temperature, the elevation-dependant temperature (ground pickup and atmosphere) and the CMB monopole. However, these are constant or fairly constant over the 90-min blocks. Since we actually remove a smooth function of time in the processing of the real data, we can neglect these contributions here. This means that we have also removed the time average of the combined signal at each frequency for each of the five blocks' TOD.

The five observational blocks are then combined at each frequency to create a 2D map in RA and Dec. The TOD samples (Y) are averaged into map pixels (X) using a pointing matrix (A):

$$ X = (A'_{1}A)^{-1}A'_{1}Y, \qquad (14) $$

which is assembled using the ASTROPY coordinate library and setting the pixel spacing to be 0.4° (roughly two pixels per beam). Then, combining all the frequency maps, we have a 3D data cube. Fig. 25 shows the total system temperature map at frequency 983.5 MHz.

As the Galactic emission contribution dwarfs the H1 signal, the maps require foreground cleaning before power spectra can be made. We aim to determine how residual satellite contamination corrupts a measurement of the H1 signal in addition to standard problems of foreground cleaning, for example residual foreground contamination or conversely over-cleaning and removal of H1 signal (Cunnington et al. 2023b). Therefore, we apply the same PCA foreground cleaning technique used in Cunningham et al. (2023a). We first calculated

<figure><img src="image_28.png" /><figcaption>**Figure 25.** The combined system temperature map at $\nu = 983.5$ MHz for the five observation blocks. A pixel spacing of 0.4° was applied.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Figure 26.** The percentage difference between our estimator and the H1 power spectrum. The estimator consists of the PCA-cleaned cross-correlation power spectrum with the H1 signal itself and is a proxy for cross-correlations with galaxy surveys. The maroon curve represents H1 with synchrotron emission as the only contaminant and one PCA mode removed. The remaining curves contain residual satellite contamination with one (fuchsia), two (gold), 3 (teal), 6 (pastel blue), and 10 (orange) modes removed. The bottom panel highlights a zoomed-in version of the top panel with a larger number of modes removed.</figcaption></figure>

the frequency-frequency covariance matrix of our data cube, then removed the largest eigenvalues under the assumption that the signal with the strongest frequency correlation is the Galactic signal. We increased the number of modes removed until it was clear we had strong signal loss.

In order to better understand the residual contamination, we used an estimator based on the cross-correlation between the PCA-cleaned total emission maps and the input H1 maps (as convolved with our Gaussian beam model). The cross-correlation power spectra were calculated using functionality from the H1 simulations package FASTBOX<fnref n="8" />, which creates a 3D Cartesian data cube. In reality, future data will require a more robust regridding of the maps (Cunnington & Wolz 2024). The power was averaged between nine bins in $k$-space spanning from $k = 0.02$ to $k = 0.7$. Twenty realizations of the H1 signal were created so that we can calculate the mean and standard deviation of the mean for the estimators.

Fig. 26 shows the percentage deviation of the cross-correlation power spectra with respect to the H1 signal, calculated as:

$$ \frac{\Delta P(k)}{P_{H1}(k)} = \frac{\text{cross}(\text{cleaned})}{\text{auto}(H1)} - 1 \qquad (15) $$

We considered two scenarios: the first where the only contaminant to the H1 signal is from diffuse Galactic synchrotron emission

<footnote n="8">⁸https://github.com/philbull/FastBox</footnote>
<figure><img src="image_30.png" /><figcaption>**Figure 27.** The autocorrelation for the H1 power spectrum with different modes cleaned with PCA. The fuchsia curve represents H1 with synchrotron emission as the only contaminant and two modes removed. Meanwhile, the remaining curves contain residual satellite contamination with 13 (teal), 15 (gold), and 20 (navy blue) modes removed. The grey straight line is the ratio for the H1 only.</figcaption></figure>

and the second where the contaminants are the combined Galactic and residual satellite temperatures. It can be seen that the removal of a single mode of frequency-correlated contaminants is sufficient for the scenario where only H1 and synchrotron emission are present in the maps. For the inclusion of residual satellite emission, Fig. 26 shows the effect of different numbers of modes removed, choosing to zoom in on the most significant number of modes removed (bottom panel) due to the scales of the y-axis.

We know that the satellite contamination is not smooth in frequency, so we do not necessarily expect PCA cleaning to work. Still, the power is reduced as we increase (slightly) the number of modes, showing that the PCA is able to distinguish the satellite signal from H1 up to a certain point. At the largest *k*-scales, the residual satellite emission always adds power regardless of the level of cleaning, although increasing the number of modes decreases power. Conversely, at the smaller *k*-scales, the H1 emission is dominant, and so all the power spectra start to converge to zero. We highlight the line of $\Delta P(k)/P_{H1}(k) = 0$ across all *k*-scales in grey as a visual aide.

Between the largest and smallest scales, the residual satellite contamination can be seen to cause fluctuations in power, which seems to imply a level of mode mixing between the residual satellite contamination, residual foreground contamination and the H1 signal. Looking at the cross-power spectra for one mode removed compared to 10 modes removed, it can be seen that the level of mode mixing is inversely proportional to the number of modes removed, i.e. the more residual satellite contamination cleaned, the smoother the power spectrum ratio.

Although removing 10 modes from the total emission maps results in the smoothest cross-correlation, we have over-cleaned to such an extent that the H1 cannot be detected until $k > 0.3 \text{ Mpc}^{-1}$. Therefore, we would argue that the optimum number of modes removed is three as this preserves H1, albeit slightly contaminated, across the largest number of *k* scales. We note that these results are based on the assumption that our simulations can be extended to these frequencies since no clear satellite signal is seen in this 'cosmology' frequency band. Moreover, due to the intermittent nature of the satellite contamination, we expect this RFI contamination to average down as more blocks are combined.

Finally, we also consider looking at the H1 in the autocorrelation (see Fig. 27), which required a larger number of modes to be removed. With 20 modes removed, we could detect H1 in the auto at the two

largest *k* values but with strong signal loss across larger scales. Note as well that any robust conclusion regarding the number of modes removed requires an estimation of the signal loss, which is beyond the scope of this work.

## 6 CONCLUSION

The ability to determine the effects of RFI from satellite emissions has been and will continue to be a growing concern as more satellites enter Earth's orbit and the sensitivity of radio telescopes increase. Until now, flagging has been radio astronomers' primary tool in dealing with unwanted RFI from both terrestrial and orbital sources. However, even when removing contaminated channels, we fear emissions from these sources would progress into channels that we assumed to be clean. To investigate these effects, this paper put forth a set of detailed simulations capable of predicting the signal of GNSS satellites across the sky and calibrated these to real data from MeerKAT SD observations. Note that, while we have not gone into detail concerning other telecommunication satellites, particularly those in a geostationary orbit, we have looked at the possible impact these satellites have on our data. In Appendix C, we assessed the movement of two satellites that operate in the *L* band. However, we saw no changes in the data that would correspond to the nearby movement of the satellites. This reaffirms ITU regulations, that this *L*-sub-band is restricted to RNSS only.

We constructed a method to calibrate the usually flagged satellite emission within the 1140–1310 MHz region and applied it to data from the MeerKAT SD H1 Intensity Mapping project (W21). We allowed the power of the different satellite signals to vary in the simulation and calibrated their amplitude with the data. This was done for all frequencies and timestamps. We found that our calibration of the simulation begins to break down when the observed signal becomes too strong, which we attributed to the system going non-linear (possibly from the Analogue-to-Digital Converter). This usually occurs when a satellite operating in the frequency range of interest breaches the 5° angular separation distance from the telescope's pointing. To remedy this, we applied various masking scenarios ('angular', 'thermal', and 'full-thermal'). A 5° cut seems to perform well but leads to a high percentage of flagging. We found that masking all the timestamps when the signal goes above a certain threshold performs as well or better while keeping more data. We also noted some temporal variation on the overall amplitudes of the signals. This was likely not due to a true time dependence of the power but because of other uncertainties that were absorbed by the overall amplitudes. Of note is the fact that satellites in a constellation may have different signals while we assume they all share the same signals. Nevertheless, our calibrated simulations were able to capture most of the structures seen in the data. In this paper, we focused on building a realistic simulation while calibrating it with overall, constellation-dependent, amplitudes. While it sufficed for the purposes of testing the simulation regimes of validity, a better model, accounting for all uncertainties in the forward modelling, is required if one wants a Bayesian approach to the problem. This will also allow to provide a simulation valid at all times without requiring time dependent amplitudes. To achieve this, it will be necessary to increase the complexity of the parametrization and produce a model that fully characterizes each satellite and include possible priors. Fully sampling the posterior distribution of such model will be computationally quite expensive and is left for future work.

Finally, we looked at the effects of the satellite emission propagating to the lower frequency channels, mainly the 973–1015 MHz
‘cosmology’ frequency range. We showed that the simulated satellite signal was well within the observational data uncertainties for any given time dump and channel. The next step was to examine how satellite emission would affect the auto and cross-correlation of H1 power spectrum at this frequency range. To check this, we ran end-to-end simulations for a few observational blocks that included H1 and foregrounds. We found that the satellite signal strongly contaminates the H1 power spectrum. By applying PCA cleaning, we were able to recover the power spectrum reasonably well in cross-correlation, when removing more than three modes. The auto-power spectrum, however, was strongly contaminated, requiring up to 20 modes and leading to strong signal loss. This indicates that satellite contamination could be the dominant contaminant in our current L-band analysis, although such a conclusion requires some extrapolation in our knowledge of the satellite out-of-band emission.

Satellite communication signals, and RFI in general, are contaminants to radio observations. If all satellite details were publicly available or known, one could just subtract the satellite signals from the cosmology band. Such total knowledge would allow us to subtract their signal from the data instead of flagging the data. Presumably, some flagging would still be necessary where the satellite signal is very strong and would make the receivers’ response non-linear. But one can use such knowledge even further and use the satellites’ emission as calibrators in the sky, using them to fit for the gain and other telescope parameters (for example the beam model). In reality, we do not have a perfect understanding of the GNSS and a substantial amount of information is either unknown or classified. But one can still develop models to jointly fit such satellite unknowns and telescope parameters (see for example Finlay et al. 2023, in the case of an interferometer). Using the work developed here to model the GNSS emission, together with the approach of Wang et al. (2021b) to find the gain solution in single dish mode in a Bayesian framework is a path that one ought to explore in the future.

## ACKNOWLEDGEMENTS

We thank Chris Finlay for useful discussions throughout the work and Wenkai Hu for comments on the draft. BNE and MGS acknowledge support from the South African Radio Astronomy Observatory and National Research Foundation (grant no. 84156). We acknowledge using the ILIFU<fnref n="9" /> cloud computing facility through the Inter-University Institute for Data Intensive Astronomy. The MeerKAT telescope is operated by the South African Radio Astronomy Observatory, which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation. JF thanks the support of Fundação para a Ciência e a Tecnologia (FCT) through the research grants UIDB/04434/2020 and UIDP/04434/2020 and through the Investigador FCT Contract No. 2020.02633.CEECIND/CP1631/CT0002. This result is part of a project that has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement no. 948764). SEH acknowledges funding from the Science and Technology Facilities Council (STFC) Consolidated grant (ST/P000649/1) and the UKSA grant (ST/Y005945/1) funding LiteBIRD foreground activities. SC is

supported by a UK Research and Innovation Future Leaders Fellowship grant [MR/V026437/1]. IPC is supported by the European Union within the Next Generation EU programme [PNRR-4-2-1.2 project No. SOE_0000136]. AP is a UK Research and Innovation Future Leaders Fellow [grant MR/X005399/1]. LW is a UK Research and Innovation Future Leaders Fellow [grant MR/V026437/1].

## DATA AVAILABILITY

Both simulations and data products are available on request. Access to the raw data used in the analysis is public (for access information please contact archive@ska.ac.za).

## REFERENCES

* Asad K. M. B. et al., 2021, MNRAS, 502, 2970

* Astropy Collaboration et al., 2013, A&A, 558, A33

* Baan W., 2010, Proc. Sci., Proceedings of RFI Mitigation Workshop, Vol. 107, SISSA, Trieste, PoS#001

* Baan W. A., 2019, J. Astron. Instru., 08, 1940010

* Battye R. A., Browne I. W. A., Dickinson C., Heron G., Maffei B., Pourtsidou A., 2013, MNRAS, 434, 1239

* Bentum M., Boonstra A. J., Baan W., 2010, Proc. Sci., Proceedings of RFI Mitigation Workshop, Vol. 107, SISSA, Trieste, PoS#009

* Bharadwaj S., Nath B. B., Sethi S. K., 2001, JA&A, 22, 21

* Bigot-Sazy M. A. et al., 2016, in Qain L., Li D., eds, ASP Conf. Ser. Vol.502, Frontiers in Radio Astronomy and FAST Early Sciences Symposium 2015. Astron. Soc. Pac., San Francisco, p. 41

* Buch K. D., Kale R., Muley M., Kudale S., Ajithkumar B., 2023, JA&A, 44, 37

* Bull P., Ferreira P. G., Patel P., Santos M. G., 2015, ApJ, 803, 21

* Chang T.-C., Pen U.-L., Peterson J. B., McDonald P., 2008, Phys. Rev. Lett., 100, 091303

* Chang T.-C., Pen U.-L., Bandura K., Peterson J. B., 2010, Nature, 466, 463

* Cunningham S., Wolz L., 2024, MNRAS, 528, 5586

* Cunningham S. et al., 2023a, MNRAS, 518, 6262

* Cunningham S. et al., 2023b, MNRAS, 523, 2453

* de Villiers M. S., 2023, AJ, 165, 78

* Di Vruno F., Winkel B., Bassa C. G., Józsa G. I. G., Brentjens M. A., Jessner A., Garrington S., 2023, A&A, 676, A75

* Duffett-Smith P., Zwart J., 2017, Coordinate Systems. Cambridge Univ. Press, Cambridge, p. 66

* EUSPA, 2024, Galileo, Available at: https://www.euspa.europa.eu/eu-space-programme/galileo

* Finlay C., Bassett B. A., Kunz M., Oozeer N., 2023, MNRAS, 524, 3231

* Goedhart S., 2024, Radio Frequency Interference (RFI), Available at: https://skaafrica.atlassian.net/wiki/spaces/ESDKB/pages/305332225/Radio+Frequency+Interference+RFI

* Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759

* Harper S. E., Dickinson C., 2018, MNRAS, 479, 2024

* Hobbs G. et al., 2020, Publ. Astron. Soc. Aust., 37, e012

* ITU, 2006, Protection of the Radio Astronomy Aervice from Transmitters Operating in Adjacent Bands, Recommendation ITU-R RA.517-4, 2

* ITU, 2015, Unwanted Emissions in the Out-of-Band Domain, Recommendation ITU-R SM.1541-6, 23

* ITU, 2021, Characteristics of Radio Quiet Zones, Report ITU-R RA.2259-1, 25

* Irfan M. O., Li Y., Santos M. G., Bull P., Gu J., Cunningham S., Grainge K., Wang J., 2024, MNRAS, 527, 4717

* Lewis A., Bridle S., 2002, Phys. Rev. D, 66, 103511

* Masui K. W. et al., 2013, ApJ, 763, L20

<footnote n="9">⁹www.ilifu.ac.za</footnote>
Matshawule S. D., Spinelli M., Santos M. G., Ngobese S., 2020, MNRAS, 506, 5075

Mauch T. et al., 2020, ApJ, 888, 61

MeerKLASS Collaboration et al., 2024, preprint (arXiv:2407.21626)

Miville-Deschênes M.-A., Ysard N., Lavabre A., Ponthieu N., Macías-Pérez J. F., Aumont J., Bernard J. P., 2008, A&A, 490, 1093

NRAO, 2024, Green Bank Interference Protection Group, Available at: https://www.gb.nrao.edu/IPG/

Offringa A. R. et al., 2012, A&A, 549, A11

Remazeilles M., Dickinson C., Banday A. J., Bigot-Sazy M. A., Ghosh T., 2015, MNRAS, 451, 4311

Rhodes B., 2019, Skyfield: High precision research-grade positions for planets and Earth satellites generator, record ascl:1907.024

Rodríguez J. A. A., Hein G.W., 2007, GNSS Signals and Spectra, Available at: https://www.unoosa.org/documents/pdf/icg/activities/2007/icg2/presentations/24.pdf

Rodríguez J. A. A., 2008, Doctoral thesis, University faf Munich

Santos M. et al., 2015, Proc. Sci., Advancing Astrophysics with the Square Kilometre Array, Vol. 215, SISSA, Trieste, Pos#019

Santos M. et al., 2017, Proc. Sci., MeerKAT Science: On the Pathway to the SKA, Vol. 277, SISSA, Trieste, PoS#032

Sihlangu I., Oozeer N., Bassett B. A., 2020, preprint (arXiv:2008.08877)

Teunissen P., Montenbruck O., eds, 2017a, Springer Handbook of Global Navigation Satellite Systems, 4.3 Modulation Schemes, Springer International Publishing, Switzerland, p. 107

Teunissen P., Montenbruck O., eds, 2017b, Annex B: GNSS Parameters Springer Handbook of Global Navigation Satellite Systems, Springer International Publishing, Switzerland, p. 1235

Virtanen P. et al., 2020, Nat. Methods, 17, 261

Wang Y. et al., 2021a, Res. Astron. Astrophys., 21, 018

Wang J. et al., 2021b, MNRAS, 505, 3698

Wehus I. K. et al., 2017, A&A, 597, A131

## APPENDIX A: GNSS POWER SPECTRUM DENSITY

Below are various binary offset carriers (BOC) used in this paper. We have incorporated them into our GNSS simulation and examples of their spectral shape we shown in Section 2.3. These are a combination of sine and cosine sub-carriers. Here, the sub-carrier rate is denoted as $n_s$. A value $n$ is defined as $n = 2n_s/n_c$, which defines the number of half periods of the sub-carrier in one chip interval duration, where $n_c$ is denoted as the chip rate. We denote the sine and cosine sub-carrier for the BOC as $BOC_{sin}$ and $BOC_{cos}$, respectively. For simplicity we will set $w_c = \pi/n_c f_0$ and $w_s = \pi/n_s f_0$. $f_0$ is the reference frequency of the signal and is set at 1.023 MHz to obtain the spectrum $S_{xx} = |P_{xx}|^2$ where $xx$ denotes the specific modulation in operation.

### (i) BPSK:

$$P_{\text{BPSK}_{n_c}}(\nu) = \frac{\sinh(\pi\nu/[n_c f_0])}{\sqrt{n_c f_0}}. \qquad (A1)$$

### (ii) BOC$_{sin}$ and $n$ even:

$$P_{\text{BOC}(n_s, n_c)}(\nu) = \sqrt{n_c f_0} \frac{\sin(w_c \nu)}{\pi \nu} \tan\left(\frac{w_s \nu}{2}\right) \qquad (A2)$$

### (iii) BOC$_{sin}$ and $n$ odd:

$$P_{\text{BOC}(n_s, n_c)}(\nu) = \sqrt{n_c f_0} \frac{\cos(w_c \nu)}{\pi \nu} \tan\left(\frac{w_s \nu}{2}\right) \qquad (A3)$$

### (iv) BOC$_{cos}$ and $n$ even:

$$P_{\text{BOC}(n_s, n_c)}(\nu) = \sqrt{n_c f_0} \frac{\sin(w_c \nu)}{\pi \nu} \frac{1 - \cos\left(\frac{w_s \nu}{2}\right)}{\cos\left(\frac{w_s \nu}{2}\right)} \qquad (A4)$$

### (v) BOC$_{cos}$ and $n$ odd:

$$P_{\text{BOC}(n_s, n_c)}(\nu) = \sqrt{n_c f_0} \frac{\cos(w_c \nu)}{\pi \nu} \frac{1 - \cos\left(\frac{w_s \nu}{2}\right)}{\cos\left(\frac{w_s \nu}{2}\right)} \qquad (A5)$$

### (vi) altBOC and $n$ even (Rodríguez 2008):

$$S_{\text{alt}(n_s, n_c)}(\nu) = \frac{4n_c f_0}{\pi^2} \frac{\sin^2(w_c \nu)}{\cos^2\left(\frac{w_s \nu}{2}\right)} \phi_{n_s, n_c}(\nu) \qquad (A6)$$

### (vii) altBOC and $n$ odd:

$$S_{\text{alt}(n_s, n_c)}(\nu) = \frac{4n_c f_0}{\pi^2} \frac{\cos^2(w_c \nu)}{\cos^2\left(\frac{w_s \nu}{2}\right)} \phi_{n_s, n_c}(\nu) \qquad (A7)$$

with

$$\phi_{n_s, n_c}(\nu) = \left[ \begin{aligned} & \cos^2\left(\frac{w_s \nu}{2}\right) - \cos\left(\frac{w_s \nu}{2}\right) \\ & -2\cos\left(\frac{w_s \nu}{2}\right)\cos\left(\frac{w_s \nu}{4}\right) + 2 \end{aligned} \right]$$

### (viii) TMBOC

Note that, Composite BOC (CBOC) Rodríguez & Hein (2007) has the same modulation structure as TMBOC.

$$S_{\text{TMBOC}(6,1,1/11)}(\nu) = \frac{10}{11} S_{\text{BOC}_{\sin(1,1)}(\nu)} + \frac{1}{11} S_{\text{BOC}_{\sin(6,1)}(\nu)} \qquad (A8)$$

$$= \frac{f_0}{11\pi^2\nu^2} \sin^2\left(\frac{\pi\nu}{f_0}\right) \left[ 10\tan^2\left(\frac{\pi\nu}{2f_0}\right) + \tan^2\left(\frac{\pi\nu}{1f_0}\right) \right] \qquad (A9)$$

## APPENDIX B: GNSS SIGNAL INFORMATION

Table B1 shows the number of satellites above the horizon for the telescope during the observation. In Table B2, we summarise the details we used to model the observed emission of the satellite constellation's signals.

## APPENDIX C: GEOSTATIONARY SATELLITES

Besides the navigational satellites, many other satellites reside in a geostationary orbit around the Earth. These satellites generally

<table><caption>Table B1. Number of satellites per constellation included, based on MeerKAT's FoV during the scanning period, for the observation on 2019/02/25–00:40:11 GMT.</caption><thead><tr><th>Const.</th><th>GPS</th><th>GLO</th><th>GAL</th><th>BEI</th><th>IRNSS</th><th>SBAS</th><th>QZS</th></tr></thead><tbody><tr><td>#</td><td>15</td><td>13</td><td>13</td><td>18</td><td>5</td><td>9</td><td>0</td></tr></tbody></table>
<table><caption>**Table B2.** List of satellite signals used in the simulation. Only signals with an index were potentially present and therefore included (first column). 'Sys': constellation system; 'Band': allocation of signal transmission; 'Signal': name of radio transmission; 'Frequency': peak frequency of transmission; 'Modulation': frequency structure of signal; 'Rate': $n_c \times f_0$; '$P_t$': transmitted power of transmission; '$G_t$': antenna gain of transmission (Teunissen & Montenbruck 2017b; Harper & Dickinson 2018). The highlighted grey regions indicate the signals we considered for each of the constellations. All the satellites from a given constellation are assumed to have the same signals.</caption><thead><tr><th>#</th><th>Sys</th><th>Band</th><th>Signal</th><th>Frequency [MHz]</th><th>Modulation</th><th>Rate [MHz]</th><th>P<sub>t</sub> [dBW]</th><th>G<sub>t</sub> [dBi]</th></tr></thead><tbody><tr><td>–</td><td>GPS</td><td>L1</td><td>P(Y)</td><td>1575.420</td><td>BPSK(10)</td><td>10.2300</td><td>13.5</td><td>13.5</td></tr><tr><td>–</td><td>GPS</td><td>L1</td><td>C/A</td><td>1575.420</td><td>BPSK(1)</td><td>1.0230</td><td>16.5</td><td>13.5</td></tr><tr><td>–</td><td>GPS</td><td>L1</td><td>L1C-D</td><td>1575.420</td><td>TMBOC(6,1,4/33)</td><td>1.0230</td><td>10.0</td><td>10.0</td></tr><tr><td>–</td><td>GPS</td><td>L1</td><td>M-D</td><td>1575.420</td><td>BOC<sub>sin</sub>(10,5)</td><td>5.1150</td><td>18.2</td><td>13.5</td></tr><tr><td>1</td><td>GPS</td><td>L2</td><td>P(Y)</td><td>1227.600</td><td>BPSK(10)</td><td>10.2300</td><td>10.0</td><td>10.0</td></tr><tr><td>2</td><td>GPS</td><td>L2</td><td>L2CM</td><td>1227.600</td><td>BPSK(1)</td><td>0.5115</td><td>10.0</td><td>10.0</td></tr><tr><td>3</td><td>GPS</td><td>L2</td><td>M-D</td><td>1227.600</td><td>BOC(10,5)<sub>sin</sub></td><td>5.1150</td><td>16.0</td><td>13.5</td></tr><tr><td>4</td><td>GPS</td><td>L5</td><td>L5I</td><td>1176.450</td><td>BPSK(10)</td><td>10.2300</td><td>18.0</td><td>18.0</td></tr><tr><td></td><td>GLO</td><td>L1</td><td>L1SF(P)</td><td>1602.000</td><td>BPSK(5)</td><td>5.1100</td><td>10.0</td><td>10.0</td></tr><tr><td>5</td><td>GLO</td><td>L2</td><td>L2SF(P)</td><td>1245.100</td><td>BPSK(5)</td><td>5.1100</td><td>10.0</td><td>10.0</td></tr><tr><td>6</td><td>GLO</td><td>L2</td><td>L2OF(C/A)</td><td>1245.100</td><td>BPSK(0.5)</td><td>0.5110</td><td>10.0</td><td>10.0</td></tr><tr><td>7</td><td>GLO</td><td>L3</td><td>L3OC-D</td><td>1202.025</td><td>BPSK(10)</td><td>10.2300</td><td>10.0</td><td>10.0</td></tr><tr><td>8</td><td>GLO</td><td>L2</td><td>L2OC-D</td><td>1248.300</td><td>BPSK(1)</td><td>1.0230</td><td>13.0</td><td>12.0</td></tr><tr><td>9</td><td>GLO</td><td>L2</td><td>L2OC-P</td><td>1248.300</td><td>BOC(1,1)</td><td>0.5115</td><td>5.0</td><td>5.0</td></tr><tr><td>–</td><td>GAL</td><td>E1</td><td>OS-D(B)</td><td>1575.420</td><td>CBOC(6,1,1/11)</td><td>1.0230</td><td>10.0</td><td>10.0</td></tr><tr><td>10</td><td>GAL</td><td>E6</td><td>CS-P(C)</td><td>1278.750</td><td>BPSK(5)</td><td>5.1150</td><td>16.0</td><td>15.0</td></tr><tr><td>11</td><td>GAL</td><td>E6</td><td>PRS(A)</td><td>1278.750</td><td>BOC<sub>cos</sub>(10,5)</td><td>5.1150</td><td>18.0</td><td>16.0</td></tr><tr><td>12</td><td>GAL</td><td>E5ab</td><td>PRS(A)</td><td>1191.795</td><td>AltBOC(15,10)</td><td>10.2300</td><td>10.0</td><td>10.0</td></tr><tr><td>13</td><td>GAL</td><td>E5a</td><td>E5a-D</td><td>1176.450</td><td>AltBOC(15,10)</td><td>10.2300</td><td>6.0</td><td>6.0</td></tr><tr><td>–</td><td>BDS-2</td><td>B1-2</td><td>RS</td><td>1561.098</td><td>BPSK(2)</td><td>2.0460</td><td>10.0</td><td>10.0</td></tr><tr><td>14</td><td>BDS-2</td><td>B3</td><td>RS</td><td>1268.520</td><td>BPSK(10)</td><td>10.2300</td><td>16.0</td><td>16.0</td></tr><tr><td>15</td><td>BDS-2</td><td>B2b</td><td>OS</td><td>1207.140</td><td>BPSK(2)</td><td>2.0460</td><td>14.0</td><td>12.0</td></tr><tr><td>16</td><td>BDS-2</td><td>B2b</td><td>RS</td><td>1207.140</td><td>BPSK(10)</td><td>10.2300</td><td>18.0</td><td>18.0</td></tr><tr><td>–</td><td>BDS-3</td><td>B1-2</td><td>OS</td><td>1561.098</td><td>BPSK(2)</td><td>2.0460</td><td>10.0</td><td>10.0</td></tr><tr><td>–</td><td>BDS-3</td><td>B1</td><td>B1C-D1</td><td>1575.420</td><td>TMBOC(6,1,4/33)</td><td>1.0230</td><td>10.0</td><td>10.0</td></tr><tr><td>17</td><td>BDS-3</td><td>B3</td><td>B3C-Dm</td><td>1268.520</td><td>BPSK(10)</td><td>10.2300</td><td>15.0</td><td>13.5</td></tr><tr><td>18</td><td>BDS-3</td><td>B3</td><td>B3A-Dm</td><td>1268.520</td><td>BOC(15,2.5)</td><td>2.5575</td><td>10.0</td><td>10.0</td></tr><tr><td>–</td><td>QZS-1</td><td>L1</td><td>C/A</td><td>1575.420</td><td>BPSK(1)</td><td>1.0230</td><td>10.0</td><td>10.0</td></tr><tr><td>–</td><td>QZS-1</td><td>L1</td><td>L1C-D</td><td>1575.420</td><td>BOC(1,1)</td><td>1.0230</td><td>10.0</td><td>10.0</td></tr><tr><td>–</td><td>QZS-1</td><td>L1</td><td>L1C-D</td><td>1575.420</td><td>TMBOC(6,1,4/33)</td><td>1.0230</td><td>10.0</td><td>10.0</td></tr><tr><td>–</td><td>QZS-1</td><td>L1</td><td>SAIF</td><td>1575.420</td><td>BPSK(1)</td><td>1.0230</td><td>10.0</td><td>10.0</td></tr><tr><td>–</td><td>QZS-1</td><td>L2</td><td>L2CL</td><td>1227.600</td><td>BPSK(1)</td><td>0.5115</td><td>10.0</td><td>10.0</td></tr><tr><td>–</td><td>QZS-1</td><td>L6</td><td>L61(LEX)n</td><td>1278.750</td><td>BPSK(5)</td><td>5.1150</td><td>10.0</td><td>10.0</td></tr><tr><td>–</td><td>QZS-1</td><td>L6</td><td>L62o</td><td>1278.750</td><td>BPSK(5)</td><td>5.1150</td><td>10.0</td><td>10.0</td></tr><tr><td>–</td><td>QZS-1</td><td>L5</td><td>L5I</td><td>1176.450</td><td>BPSK(10)</td><td>10.2300</td><td>16.0</td><td>16.0</td></tr><tr><td>–</td><td>QZS-1</td><td>L5</td><td>L5Q</td><td>1176.450</td><td>BPSK(10)</td><td>10.2300</td><td>16.0</td><td>16.0</td></tr><tr><td>19</td><td>IRNSS</td><td>L5</td><td>SPS</td><td>1176.450</td><td>BPSK(1)</td><td>1.0230</td><td>16.0</td><td>14.0</td></tr><tr><td>20</td><td>IRNSS</td><td>L5</td><td>RS-D</td><td>1176.450</td><td>BOC(5,2)</td><td>2.0460</td><td>18.0</td><td>16.0</td></tr><tr><td>–</td><td>SBAS</td><td>L1</td><td>C/A</td><td>1575.420</td><td>BPSK(1)</td><td>1.0230</td><td>13.0</td><td>13.5</td></tr><tr><td>21</td><td>SBAS</td><td>L5</td><td>L5I</td><td>1176.450</td><td>BPSK(10)</td><td>10.2300</td><td>18.0</td><td>16.0</td></tr></tbody></table>

operate at higher frequencies (K, KU, Ka-band) than their navigational counterparts and belong to not only government agencies but private companies such as Starlink<fnref n="10" />, OneWeb<fnref n="11" />, among others for communication purposes.

In our constructed satellite catalogue, we excluded the telecommunication geostationary satellites due to the many unknown variables associated with them, in particular, their modulation and signal behaviour. However, we can still track and locate satellites that enter the primary beam's 5° boundary. In our analysis, we found 17 satellites passed into the 5°, of which two were actively transmitting

in the L-band. We show their $\theta$ movement concerning MeerKAT's pointing in Fig. C1.

With the additional positional information, we created a new mask targeting the timestamps $\le 5°$ and combined this mask with the already in use 5° mask of the GAL and GLO satellites. The additional masking removed ~83 per cent of the data, which is less overall than the temporal masking but ~56 per cent more than that of the 5°. In the case of the FoM, we obtained a $C_1$ value of 30 855.53 and a $C_2$ value of 0.54. This shows a factor of 3 improvements when compared to only the 5° case (see Table 2).

However, we suspect that these two geostationary satellites do not emit within our frequency range of interest (1100–1350 MHz) for two reasons. (1) According to regulations the International Telecommunication Union sets, our frequency range is reserved for navigational satellite emission. (2) When comparing Fig. C1 with

<footnote n="10">¹⁰https://www.starlink.com/</footnote>

<footnote n="11">¹¹https://oneweb.net/</footnote>
<figure><img src="image_31.png" /><figcaption>**Figure C1.** The angular distance (equation 1 of two Russian geostationary communication satellites that transmit in the L band and cross below the 5° boundary).</figcaption></figure>

the bottom panel of Figs 15 and 18, it seems that the model can account for all the peaks of emission and overlay and accurate fit. This implies that there is no need to include extra satellite emission in the model.

This paper has been typeset from a LaTeX/LATEX file prepared by the author.