Issue 
A&A
Volume 657, January 2022



Article Number  A75  
Number of page(s)  19  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202142101  
Published online  14 January 2022 
Analytical simulations of the effect of satellite constellations on optical and nearinfrared observations
^{1}
ASTRON Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The Netherlands
email: bassa@astron.nl
^{2}
European Southern Observatory, KarlSchwarzschildStrasse 2, 85748 Garching bei München, Germany
email: ohainaut@eso.org
^{3}
Observatorio de Calar Alto, Sierra de los Filabres, 04550 Gérgal (Almería), Spain
email: dgaladi@caha.es
Received:
27
August
2021
Accepted:
14
October
2021
Context. The number of satellites in lowEarth orbit is increasing rapidly and many tens of thousands of satellites are expected to be launched in the coming years. There is a strong concern among the astronomical community about the contamination of optical and nearinfrared observations by satellite trails, what has led to several initial investigations of the impact of large satellite constellations.
Aims. We expand the impact analysis of such constellations on groundbased optical and nearinfrared astronomical observations in a more rigorous and quantitative way, using updated constellation information and considering imagers and spectrographs and their very different characteristics.
Methods. We introduce an analytical method that allows us to rapidly and accurately evaluate the effect of a very large number of satellites, accounting for their magnitudes and the effect of trailing of the satellite image during the exposure. We use this to evaluate the impact on a series of representative instruments, including imagers (traditional narrow field instruments, widefield survey cameras, and astrophotographic cameras) and spectrographs (longslit and fibrefed), taking their limiting magnitude into account.
Results. Confirming earlier findings, the effect of satellite trails is more damaging for highaltitude satellites, on widefield instruments, or essentially during the first and last hours of the night. Thanks to their brighter limiting magnitudes, low and midresolution spectrographs will be less affected, but the contamination will be at about the same level as that of the science signal, introducing additional challenges. Highresolution spectrographs will essentially be immune. We propose a series of mitigating measures, including one that uses the described simulation method to optimise the scheduling of the observations. We conclude that no single mitigation measure will solve the problem of satellite trails for all instruments and all science cases.
Key words: light pollution / site testing / space vehicles / telescopes / surveys
© ESO 2022
1. Introduction: Satellite constellations
In the most general sense, a constellation of artificial satellites is a set of spacecraft that share a common design, distributed among different orbits to provide a service and/or a land coverage that cannot be achieved by means of just one single satellite. Satellite constellations have been in use since the early times of the space age and their applications range from telecommunications (civil or military), meteorology (Meteosat), remote sensing (as the twosatellite constellations Sentinel of ESA), or global navigation satellite systems (GNSS; GPS, Glonass, Galileo or Beidou systems).
Until recently, the satellite constellation with a largest number of elements was the Iridium system, with some 70 satellites in lowEarth orbit (LEO, altitude less than 2000 km), aimed to provide cell phone service with worldwide coverage. The socalled Iridium flares caused by reflected sunlight from the flat antennae of the first generation of Iridium satellites awoke the awareness of the astronomical community about possible deleterious effects of satellite constellations on astronomical observations, both optical and in radio (James 1998). Global navigation satellites systems also require worldwide coverage that, in this case, can be achieved with smaller constellations, typically of the order of 30 satellites, placed at higher orbits (∼20 000 km).
For the purpose of this paper, we consider a megaconstellation to refer to any satellite constellation consisting of a number of satellites significantly larger than the Iridium constellation, in excess of 1000 satellites. Several such systems have been proposed during the last decade, in all cases to provide fast internet access worldwide. The technical requirements for that application (large twoway bandwidth, short delay time, large number of potential users) lead to designs of constellations made up from huge numbers of satellites, counted in the thousands, placed at lowEarth orbit.
The contrast with traditional constellations is overwhelming. The examples mentioned above, or the plans to launch new constellations of nano and picosatellites, imply challenges of their own and contribute to the issues of orbital crowding, space debris management, radioelectric noise, etc., but they imply numbers of satellites orders of magnitude smaller and, most often, with elements much fainter, than the megaconstellations discussed in here.
Several new generation megaconstellations are in the planning stages and some number of satellites of two of these constellations, namely SpaceX’s Starlink and OneWeb, have already been placed in orbit. In Table 1 we list some of the planned constellation configurations for which information is publicly available. Several more satellite operators indicated their intentions to build about a dozen of additional similar constellations. Many companies are also planning constellations of much smaller satellites (cubesats, nanosats...), which are less relevant for optical astronomy. We note that satellite operators, even those with satellites already in orbit, frequently modify the configuration of their constellation projects^{1}. Furthermore, other satellite operators have indicated their intentions to launch megaconstellations, but without so far submitting any application or publishing any details. Therefore, the configurations used in this paper should be treated as representative and their impact on astronomical observations obtained from the following simulations and their interpretation are an illustration of the situation as it could be in the late 2020’s. In first approximation, the numbers scale linearly with the total number of satellites in the constellation, so the effects be scaled.
Orbit configurations.
Of all these projects, the SpaceX Starlink constellation is clearly in the forefront. Since the launch of their first group of satellites in May 2019, the trains of very bright pearls that form the satellites illuminated by the Sun have caught the attention even of the general public and have triggered the alarms of the astronomical community (see, for example Witze 2019). These bright trains are formed only at the first stages of the deployment of each Starlink launch, with the satellites becoming fainter as they later climb to the final operational orbits and attain an operational attitude that reduces their apparent brightness.
Our aim in this paper is to assess the impact of satellite megaconstellations on groundbased optical and nearinfrared astronomy from computer simulations of two kinds. Such simulations are described in Sect. 2, where we depict the usual discrete simulations and, also, a new approach consisting on formulating statistical predictions from analytical probability density functions that describe the main properties of the constellations. The same section includes consideration on how to estimate the apparent brightness of satellites illuminated by sunlight, with special attention to the concept of effective magnitude.
The effect on observations is addressed in Sect. 3. There we provide details on the specific details of a set of simulations, their implications for several observing modes (direct imaging and spectroscopy, including multifibre instruments) and we discuss some possibilities to mitigate the impact in Sect. 4. The appendix includes the details on how the equations of our analytical model are derived.
2. Simulations
2.1. Visibility
A satellite will be visible at optical and nearinfrared wavelengths if it is above the horizon for the geographical location of the observer and if it is illuminated by the Sun, making it to reflect and diffuse sunlight back to the observer. Hence, the visibility of a given satellite will depend upon its location within its orbit, the geographic location of the observer, and the location of the Sun, all as a function of time.
While the motion of a satellite in the Earth’s gravitational field and through the higher residual layers of the atmosphere is usually computed through either numerical integration or using simplified analytical models (e.g., Montenbruck & Gill 2000; Vallado et al. 2006), these offer an accuracy that is not required for the statistical results that we want to derive from the simulations presented here. Instead, we assume that the satellites move in circular Keplerian orbits, and we neglect the effects of atmospheric drag, the lack of spherical symmetry of the Earth’s gravitational potential and perturbations by the Sun and the Moon. Furthermore, for simplicity, satellite positions are defined in an inertial geocentric reference system (Celestial Intermediate Reference System), neglecting the small effects of precession, nutation and Earth orientation parameters. These assumptions allow us to consider only the idealised motion of the satellites, the motion of the observer (uniform Earth rotation) and the apparent motion of the Sun.
Satellite constellations for worldwide coverage are generally configured as Walker constellations (Walker 1984). Here, satellites are in circular orbits of distinct orbital shells, with each shell described by its orbital altitude h_{sat} and orbital inclination i. Within each shell, equally spaced orbital planes are populated by satellites, also equally spaced within each plane. In Table 1 we use n_{sat} to denote the number of satellites in a single orbital plane, and n_{plane} to indicate the number of orbital planes. The total number of satellites within each shell is N_{sat} = n_{sat} × n_{plane}, and the total number of satellites in a constellation the sum of the number of satellites in each shell.
In this paper we consider only idealised, complete constellations. Individual satellites, as well as the trains of satellites in very low orbits right after launch and the satellites in neartoreentry orbits, are not considered. Even taking into account the large number of launches required to replenish the constellations and the hopefully equally large number of satellites on endoflife orbits, the number of such satellites will be one to two orders of magnitude smaller than the number of satellites in operation. Also, as these satellites are in lower orbits their impact is concentrated during the very beginning and end of night. The constellations that are now in the building phase have a larger proportion of satellites in nonoperational orbits, but at the final steady state their contribution to the overall situation, caused by the constellations in operation, will therefore be small.
We expand on the work by Hainaut & Williams (2020), who used a simplistic geometric approximation to estimate the density of satellites and their effects. That model had the advantage of being extremely fast and numerically lean, and to yield acceptable results for midlatitude observatories, but it had the major shortcoming not to account for latitude effects. This aspect is rigorously addressed in the following sections.
2.1.1. Discrete simulations
To obtain a realisation of a satellite constellation at a given time, the positions and velocities of all satellites in the constellation are computed using the assumptions listed earlier. Figure 1 shows such a realisation for the Starlink Generation 2 constellation on a map of the Earth near the June solstice when the Sun is at its highest northern declination. Hence, for locations at high geographical latitudes in the northern hemisphere, satellites will be visible throughout the whole night, whereas for lower latitudes and in the southern hemisphere satellites will only be visible around twilight. This figure depicts overdensities of satellites at geographical latitudes near the orbital inclination of different orbital shells.
Fig. 1.
Discrete realisation of the Starlink Generation 2 constellation overlaid on a map of Earth. The day and night sides of Earth are shown for 2022 June 21 at 23:00UTC, near the June solstice. Satellites are shown as either sunlit (yellow), or eclipsed by the Earth (black) and above the horizon for two geographical locations, Vera Rubin Observatory in Chile (latitude −30°) and London in the United Kingdom (latitude ∼50°). Overdensities of satellites in geographic latitude b are visible at latitudes b = 30°, 40° and 53°, corresponding to the inclination of the most populous orbital shells. Due to their orbital altitudes, satellites remain visible for locations where the Sun is below the horizon and will remain visible throughout all the night for locations at high geographical latitudes. 
These overdensities are also seen in the allsky maps shown in Fig. 2 for two locations, the Vera Rubin Observatory (VRO) at Cerro Pachón in Chile (at latitude −30°) and London in the United Kingdom (at latitude +50°). The increased distance to the orbital shells when looking at lower elevations above the horizon leads to larger sampling volumes, which explains that these allsky maps show increasing satellite densities towards the horizon, though additional overdensities exist near the projected location of the northern and southern limits of the most populous orbital shells.
Fig. 2.
Discrete realisation of the Starlink Generation 2 constellation from Fig. 1. The satellites are plotted on all sky maps representing the celestial sphere above Vera Rubin Observatory (left) and London (right). Satellites are shown as sunlit (yellow) or eclipsed by the Earth (black). The number of satellites increases towards the horizon because, for a given solid angle, looking at lower elevations from the horizon the distance to the shell surface increases and, hence, the field of view contains larger portions of the orbital shells. Due to the orbital inclination and orbital altitude of the orbital shells, overdensities of satellites are also present in these maps. The map for London shows this very clearly for the two orbital shells with i = 53° and h = 345 km (7178 satellites) and 499 km (4000 satellites), where more satellites are present south from the dotted line. The dotted line in the all sky map for Vera Rubin Observatory is for the orbital shell with i = 30°, h = 328 km, with a larger number of satellites north from that line. 
This discrete approach has been followed to assess the impact of satellite constellations in terms of the number of satellites that are visible during the night above a given elevation from the horizon, as shown, for instance, in McDowell (2020). To this end, a specific constellation is selected (in terms of shell number and, for each of them, the altitude, orbital inclination and number of satellites), an observatory location is specified (in general only latitude is relevant) and the illumination conditions are fixed through Sun declination and hour angle. The simulation of the Keplerian movement of each satellite of the constellation allows counting the number of satellites that are above a specific elevation and their illumination conditions (sunlit or eclipsed). Additional considerations may allow computing some photometric estimates (Sect. 2.2). This kind of discrete, allsky simulations may be iterated, including some random initialisation of parameters at each run to averageout systematic effects due to the spatial and temporal texture induced by the constellation structure.
2.1.2. Number of satellite trails in an observation
Besides the allsky simulations, an observationoriented approach is also needed. This requires specifying all the parameters required for the allsky simulations and, also, selecting the observation direction (azimuth, elevation), fieldofview (FOV) and exposure time of an observation. In these observations, the motion of the satellite during the exposure will leave a satellite trail on the images. Repeating this simulation leads to the average estimate of the number of satellite trails that would affect the observation. The geometry of the problem allows also computing the position angle of the trails and the apparent angular velocity of each satellite that crosses the FOV.
In Fig. 3 we show the results of a discrete simulation on observations with different exposure times and fieldsofview. The discrete simulation used a constellation of 10 000 satellites in a single orbital shell of 100 orbital planes with 100 satellites within each plane, at 1000 km altitude and 53° inclination. For an observer at −30° latitude, the number of satellite trails visible within an observation of an exposure time and circular fieldofview towards zenith were counted and averaged over 1000 simulations.
Fig. 3.
Number of satellites visible in an exposure depending on the exposure time t_{exp} and the instrument fieldofview, specified by its circular radius R_{fov} = L/2. The effect of a discrete constellation of 10 000 satellites in a single orbital shell of 100 orbital planes and 100 satellites per plane at 1000 km altitude and 53° inclination is simulated for an observer at −30° latitude, observing towards zenith. The points are the results of the discrete simulations for different values of exposure time t_{exp} and fieldofview (). The solid lines are predictions based on Eq. (1). 
The simulations show a dependence with both fieldofview and exposure time (solid dots). This relation can be understood as the number of satellites present at the start of the exposure, plus those traveling through the fieldofview during the exposure, and has the form
Here, ρ_{sat} is the instantaneous satellite number density, meaning the number of satellites per unit area on the sky, and ω_{sat} the angular velocity of the satellites in the direction of the exposure. The exposure itself has fieldofview area of A_{fov} and width L_{fov}, and exposure time t_{eff}. For comparison with the simulations we use a fieldofview with a circular radius R_{fov}, such that and L_{fov} = 2R_{fov}.
The simulations provide the instantaneous satellite density ρ_{sat} and the average angular velocity ω_{sat} and using these values, Eq. (1) provides the predicted number of trails for the different observation parameters, plotted as lines in Fig. 3. The predictions match the simulations, though for short exposures and/or small fieldsofview, the simulations suffer from noise due to the few satellites passing through the observation.
Equation (1) is valid for a single shell in a constellation, as the instantaneous density ρ_{sat} and angular velocity ω_{sat} depend on the properties of the shell. To obtain the total number of trails in an observation for a satellite constellation with multiple shells, Eq. (1) can be computed for each shell and the results summed.
2.1.3. Analytical simulations
The averaging over many randomly initialised parameters of a satellite constellation for the discrete simulations is computationally expensive, and may lead to noise due to insufficient satellites passing through observations with small fieldsofview and/or short exposures to obtain a valid average (see Fig. 3). Given that the averaging over randomly initialised parameters has the effect of smoothing out the satellite locations within a given orbital shell, we instead treat the satellite locations as probability density functions, which have the advantage that these are analytical expressions.
Figure 1 shows that the satellite locations are uniformly spread in geocentric longitude, but strongly peaked towards geocentric latitudes ϕ close to the values equal to the orbital inclination −i and +i of a constellation shell. For a satellite with nodal anomaly κ measured from the ascending node, the geocentric latitude is given by sin ϕ = sin κ sin i and its probability distribution follows the Arcsine probability distribution^{2}. For a single satellite at orbital altitude h_{sat}, the probability density per unit surface area as a function of ϕ is:
Here R_{⊕} is the radius of the Earth. This expression is valid for ϕ< i and zero otherwise. A detailed derivation of Eq. (2) is given in Appendix A.1.
Equation (2) can be integrated over the surface of the orbital shell spanned by the fieldofview of an instrument from an observatory located on Earth to obtain the fraction of the sample present within the instrument fieldofview. For the remainder of the paper, we make the simplifying assumption that P(ϕ, i, h_{sat}) is constant over the instrument fieldofview. This assumption is generally true for the small fieldsofview under consideration in the remaining analysis. It has the advantage that it removes any dependence on the precise shape and orientation of the instrument fieldofview, and instead solely depends on the sky area covered by the fieldofview. The assumption allows us to evaluate Eq. (2) for a lineofsight (specified by azimuth and elevation) of the observation from an observatory on Earth intersecting the orbital shell at distance d and with an impact angle α (α = 90° at zenith and α < 90° for lower elevations). The instantaneous surface density ρ_{sat} can then be obtained by scaling the probability by the surface area of the orbital shell covered an angular area A of 1 square degree, providing
where N_{sat} is the number of satellites in the orbital shell with inclination i and orbital altitude h_{sat}. Equations to derive the distance d and impact angle α given the location of the observatory and azimuth and elevation of the observation are given in Appendices A.2 and A.3.
The angular velocity of the satellites in an orbital shell towards the lineofsight can be determined using the equations provided in Appendix A.4. Due to the rotation of the Earth within the orbital shell of a satellite constellation, the velocity vectors project differently on the sphere of the sky and hence satellites have somewhat different angular velocities depending on their north or southbound trajectory. Given that for a full constellation, an equal amount of satellites are on northbound as well as southbound trajectories, we can take the average of both angular velocities to obtain ω_{sat}.
For a given satellite constellation, observatory latitude and observation parameters (fieldofview and exposure time), the analytical simulations predict the number of trails as a function of azimuth and elevation which is inherently static with time. The final step to complete the analytical simulation is taking into account the illumination of the different orbital shells by the Sun, as this modulates the visibility of a shell as a function of time of day and time of year. The impact of solar illumination is implemented by computing whether the intersection of a given lineofsight with an orbital shell is in the shadow of the Earth or not. If the intersection point is located in the shadow and if the satellites are not visible, we set N_{trail} = 0 for that shell in the sum of satellite trails over the different orbital shells.
Figure 4 shows all sky maps using the analytical simulations, providing the number of satellite trails from Eq. (1). The contribution of the different orbital shells of multiple satellite constellations is apparent, as is the impact of solar illumination.
Fig. 4.
Sky maps with an example of the resulting number of trails per exposures. The circles mark 0°, 10°, 20°, 30°, and 60° elevation. The observatory is located at +50° latitude; the sun at −18° elevation; the camera has a FOV with diameter L_{fov} = 6′ and the exposure time is 300 s. The satellites are those from Table 1. In the black region at the southeast, all satellites are already in the shadow of the Earth. The edges running from northeast to southwest correspond to constellation shells at lower altitudes, whose southeast parts are already in the shadow. The sharp features running from east to west correspond to the edges of constellation shells whose inclinations are close to the observatory latitude. 
2.2. Photometry
Providing a reliable estimation of the apparent brightness of the satellites is an obvious requirement for any model that intends to assess the impact of megaconstellations on astronomical observation. The celestial mechanics part of the problem admits an accurate solution that, even in a simplified frame (as described in Sects. 2.1.1 and 2.1.3), leads to sound predictions of the spatial parameters (satellite density and their motion). However, the photometric part of the problem faces additional difficulties, due to the complex geometrical and reflective properties of the satellite that, also, are different from one constellation to another.
In the visible and nearinfrared (NIR), the light from the satellite is reflected sunlight with a specular and a diffuse component. The specular reflection happens on flat panels: antennas, satellite bus, possibly also solar panels (although while the satellites are in operations, these are perpendicular to the Sun). A complete and accurate representation of the reflection by a satellite would require detailed knowledge of the shape and material of the satellites (see Walker et al. 2020, for a summary of the stateoftheart). A simplified model can be assembled from photometric observations covering a range of zenithal distance and solar illumination.
Empirical models (for instance, a flat panel model, Mallama 2020a) are being developed and theoretical approaches are also used to define photometric parameters of the satellites (Walker et al. 2020). However, as of today, the available observations are sufficient only for simplistic photometric models. Hopefully, dedicated observation campaigns will refine the characterisation of the satellites in the coming years.
2.2.1. Apparent magnitude
The simple photometric model we use is based on the Lambertian sphere model. The theoretical foundations of the Lambertian sphere model can be seen, for instance, in Karttunen et al. (1996). The different specific formulations, such as those by Hainaut & Williams (2020) and McDowell (2020), can be unified under a single formula:
In Eq. (4), m_{⊙} is the Sun’s apparent magnitude as seen from Earth, in the photometric band of interest. Typically, this is Johnson’s V band, with m_{⊙} = −26.75. The second term considers the object’s intrinsic photometric properties: is the photometric crosssection, with p the object’s geometric albedo and R_{sat} the radius of the (spherical) satellite. The third term includes several distances: d_{sat⊙} is the distance from the satellite to the Sun, d_{sat} represents the distance from the observer to the satellite. For our problem, d_{sat⊙} = 1 au. The fourth term is the correction for the solar phase α_{⊙}. Finally, kχ represents the extinction term, k being the extinction coefficient (in magnitudes per unit airmass), and χ the airmass (equal to 1/sin e_{sat} in the planeparallel approximation, with e_{sat} representing the satellite’s elevation above the horizon; here, as we know the orbits of the object, we use the exact χ = d_{sat}/h_{sat}). In the V band, k = 0.12 is a typical value (Patat et al. 2011).
For a Lambertian sphere, υ(α_{⊙}) = (1 + cos α_{⊙})/2. However, a large number of photometric measurements of Starlink original satellites (Mallama 2020b) and Starlink VisorSat (Mallama 2021) indicate an extremely weak dependency of the magnitude (corrected for distance and extinction) with the solar phase angle, a circumstance that has to be related to the morphology of the satellite, which is very different to a sphere. We therefore consider υ = 1, leading to
where both d_{sat} and R_{sat} are expressed in the same units.
The photometric crosssection is the only parameter that depends on the satellite’s physical properties in this model. Hainaut & Williams (2020) adopted m^{2} for the firstgeneration Starlink. In order to facilitate the comparisons between satellites, one introduces the absolute magnitude m_{1000 km}, normalised to a standard distance d_{sat} = 1000 km and without atmospheric extinction. With that value for d_{sat}, Eq. (5) becomes
with d_{sat} expressed in km. Sometimes, m_{550 km} = m_{1000 km} − 1.3 is used instead.
Table 2 lists the absolute magnitudes measured for different satellite types and the corresponding photometric crosssection and visual magnitude at zenith. The dispersion of the measurements of m_{1000 km} is around 0.7 mag; Fig. 5 shows an even wider dispersion. Consequently, the simplistic model presented above can only represent the general trend of the satellite magnitude, as illustrated by the comparison with observations (Fig. 5). As the simulation of the effect on the observations, in the following sections, is not very sensitive to that magnitude, this level of precision is sufficient for our purposes. In the simulations described below, we use m_{1000 km} = 7, equivalent to m_{550 km} = 5.7.
Fig. 5.
Magnitude of the satellites as a function of the range between the observer and the satellite. The dots are measurements of original Starlink satellites (Pomenis telescope, from Otarola et al. 2020), and the lines are obtained using a simplified Lambertian sphere model for two altitudes. The unknown solar phase angle contributes to the dispersion of the measurements. 
Representative magnitudes of the satellites.
2.2.2. Effective magnitude and limiting magnitude
During an exposure of duration t_{exp}, a satellite will leave a trail of length ω_{sat}t_{exp} (with ω_{sat} being the apparent angular speed of the satellite), typically much longer than the FOV of the instrument. The signal corresponding to the apparent magnitude is therefore spread along the length of the trail. The count level on the detector amounts to the light accumulated inside an individual resolution element (whose angular size is r_{inst}) during the time t_{eff} = r_{inst}/ω_{sat} that the satellite takes to cross that element. This leads to the concept of effective magnitude, m_{eff}, defined as the magnitude of a static pointlike object that, during the total exposure time t_{exp}, would produce the same accumulated intensity in one resolution element than the artificial satellite during a time t_{eff}:
Figure 6 shows the effective magnitude for an example. While not directly relevant for lowaltitude constellation satellites, it is worth noting that the dependency of m_{eff} with the altitude of the satellites is shallower than that of the apparent magnitude.
Fig. 6.
Visual (red) and effective magnitudes of a satellite at zenith, as a function of its altitude, for various exposure times (see legend). The satellite used is a Starlink VisorSat with m_{1000 km} = 7, considered as a trailed point source. 
In this approximation we are assuming that the satellite’s PSF has the same shape as the stellar PSF at the telescope focal surface, a crude approximation if the distance to the satellite is small enough for it to be spatially resolved, or out of focus, or both (Tyson et al. 2020; Ragazzoni 2020). The apparent angular width of the satellite trail is
where θ_{atm} is the stellar FWHM (the seeing, typically from a good site), D_{sat} the physical diameter of the satellite, D_{m} the diameter of the telescope mirror, and d_{sat} the distance to the satellite (Tyson et al. 2020). For an 8m telescope like the ESO Very Large Telescope (VLT), or the Simonyi Survey Telescope (formerly LSST), a 2 m satellite at an altitude of 300 to 550 km, θ_{sat} ∼ 6″ to 3″. The spreading of the signal from the satellite over this larger area will decrease its signaltonoise ratio (S/N) by up to θ_{sat}/θ_{atm} and its peak intensity by up to (θ_{sat}/θ_{atm})^{2}, that is from 2 to 4 magnitudes fainter than m_{eff} from Eq. (7).
For imaging, the resolution element is typically the seeing (of the order of 1″) for telescopic observations, or the pixel (a few to a few tens arcsec) for widefield astrophotography. Figure 7a displays the visual magnitude of the faintest satellite that will leave a trail as a function of the limiting magnitude of the instrument, for t_{exp} = 60 s imaging observations. The limiting magnitude for the satellites is obtained computing the effective magnitude (using Eq. (7)) equal to the limiting magnitude of the instrument, accounting for the size of the resolution element of the instrument and a typical apparent angular velocity for the satellite at the elevation considered. This shows that observations performed on telescopes with a diameter of a metre or more will be sensitive to all or most satellites.
Fig. 7.
Detection limit for the satellite apparent magnitude (coloured lines) as a function of the limiting magnitude of the instrument (both are at the same S/N). Panel a: results for 60 s exposures with an imager, while (b): prediction for 1200 s exposures with a spectrograph. The trailing of a satellite below the coloured line will make its effective magnitude too faint for detection. The shaded areas highlight the range where constellation satellites are likely to be detected. Two elevations (zenith and airmass 2) are considered, with the corresponding typical satellite visual magnitudes and angular velocities. 
Allsky cameras will record only the brightest satellites and flares. Only the deepest widefield astrophotography (with a limiting magnitude V ∼ 15 in t_{exp} = 1 min or fainter) will record the bulk of the satellites. More specifically, the case of widefield astrophotography is considered as follows: a wideangle lens (75° ×55°) is used with a 15 megapixel detector (leading to 1′ pixels), with 15 s exposure time. With these assumptions, an average exposure would contain ∼100 satellites (varying widely with the elevation). However, the typical effective magnitude of the satellites will fall in the V = 12 − 15 range, that is below the detection threshold (a typical value for a 15 s exposure is around 9–10). Furthermore, the dependency of the effective magnitude with time almost cancels out the effect on limiting magnitude, as illustrated at Fig. 8: a camera that is not sensitive to satellite at a given exposure time will essentially remain immune at other exposure times.
Fig. 8.
Limiting magnitude for an astrophotography camera with wideangle lens (red) as a function of the exposure time. Dark current noise is neglected (see text for the camera details). This magnitude is converted into the apparent magnitude of a satellite assuming the angular velocity typical of a Starlink at 30° elevation), using Eq. (7). 
The situation is slightly different for spectroscopy. In the case of fibrefed spectrographs, the resulting data contain no spatial information at all; for longslit spectrographs, the spatial information is available in only one direction. Except in the case of integralfield spectrographs, the data will therefore not include a telltale trail indicating the contamination. For an exposure time t_{exp} = 1200 s, representative of individual exposures in the visible, Fig. 7b displays the visual magnitude of the satellite that will be detected as a function of the limiting magnitude of the instrument. Spectrographs having a limiting magnitude brighter than V = 20 in t_{exp} = 1200 s will essentially be immune: the signal from a constellation satellite will be be too faint to be detected. That will be the case for lowresolution spectrographs on small to medium telescopes and highresolution spectrographs on large telescopes.
If the S/N of the contamination is much lower than that of the science target, the contamination will result in a small increase of the background noise, what can probably be neglected for most science cases. Also, cases where the S/N of the contamination is much larger than that of the science spectrum are trivial: the effect is obvious and the observation is lost. The situation for spectrographs with a limiting magnitude in the v = 20 − 23 range in t_{exp} = 20 min is more problematic: the satellite trail will have an S/N of 2–15, so that contamination caused by the satellite will be at a level comparable to that of the science signal. It is therefore plausible that the contamination will not be immediately apparent and will be discovered only at the time of the data analysis, where a solartype spectrum (reflected by the satellite) will be superimposed to that of the science target. These intermediate cases where both S/N values are similar are much more problematic and sciencecase dependent: if the science target is a distant galaxy, a solar spectrum will be identified as a contamination. However, if the target was a stellar object, a solar contamination might cause spurious conclusions.
3. Results
3.1. Time and Solar declination dependence
The analytical models for visibility and photometry allow us to compute their dependence on time of day as well as time of year. Similar results were already presented using a simplified geometric model by Hainaut & Williams (2020), or discrete simulations by McDowell (2020). Figure 9 shows the number of satellites illuminated by the Sun for the local summer and winter seasons for an observatory at latitude . During local summer, satellites remain visible above 30° elevation throughout the night. Figure 10 shows the visibility dependence as a function of solar elevation, for each shell and for the total populations from Table 1 Several hundred of the highaltitude satellites (OneWeb, GuoWang) would remain illuminated during a large fraction or even the totality of the night; also, highaltitude constellations contribute a larger number of satellites than larger constellations at a lower altitude.
Fig. 9.
Number of illuminated satellites above the horizon (top panel) and above 30° elevation, as a function of the local solar time, for Paranal (latitude ). The figure accounts for all the satellites from Table 1. Left is for the summer solstice () and right for the winter solstice (). The twilights are indicated with blue shading. The black line marks the total number of satellites above the elevation considered, the blue line those that are illuminated and the orange and red lines those brighter than magnitude 6 and 5, respectively, using the photometric model described in Sect. 2.2. This shows that during the summer months, a few hundred satellites would remain illuminated thorough the whole night, including a few tens above 30° elevation. 
Fig. 10.
Number of illuminated satellites above the horizon as a function of the Sun elevation, for the constellations listed in Table 1 seen from Paranal (latitude ; the dependency with latitude is not strong). The twilights are indicated with blue shadings and the elevation of the sun at midnight by grey shadings for the solstices and equinoxes. The thin dotted lines represent the individual shells and the thick lines the totals for each constellation. The upper thick black line is the grand total. This figure illustrates the importance of the altitude of a constellation: the highaltitude satellites (OneWeb, Kuiper) are visible much longer into the night than lower altitude ones. Also, for most of the night, the much smaller OneWeb constellation contributes to more satellites in the sky than the largest one. 
3.2. Spatial fine structure
An outstanding effect of the orbital shells of satellite constellations is the amount of spatial fine structure that arises in the quantity of satellites visible on the local celestial sphere, as illustrated in Figs. 2 and 4. Of course, first of all we find the effect of the Earth’s shadow, whose behaviour is, as expected, dominated by the diurnal rotation of the planet and by the interplay between observatory latitude and Sun declination. But the structure of satellite shells, combined with orbital mechanics, adds a far from negligible spatial fine structure in the quantity of satellites visible on the local celestial sphere. In particular, as Eq. (2) shows, each individual shell induces an unavoidable overdensity high in the sky over observatories placed at geographic latitudes ϕ whose absolute value is close to the orbital inclination i. The northern or southern boundary of each shell lies along a line on the local celestial sphere that crosses zenith if ϕ = i, what incidentally happens for Vera Rubin Observatory and some shells currently considered in several constellation designs. These overdensities may coincide with the culmination elevation of some key objects (let us say, for instance, LMC or SMC in the south, or M31 and M33 in the north). Given the inclination i of one orbital shell and the latitude ϕ of the observatory, the shell boundary cuts the local meridian at a declination δ that can be deduced from the following equation (justified in Appendix A.5):
3.3. Contribution to the sky brightness
The satellites, including those that are not directly detected, contribute to the sky background. To evaluate this effect, surface brightness maps were computed. The magnitude of a satellite from a constellation was evaluated using Eq. (6), then converted into flux, and finally used to weight the satellite density map described in Sect. 2.1.3. A total flux density map was obtained summing the contributions of all satellite shells and transformed into surface brightness in mag arcsec^{−2}. An example is presented in Fig. 11. In the illuminated part of the shells, the satellites contribute to a surface brightness in the 28–29 mag arcsec^{−2} range (0.3–0.7 μcd m^{−2}), with peaks around 26.5 mag arcsec^{−2} (2.7 μcd m^{−2}) at the cusps of the constellations. The surface brightness of the dark night sky is around V = 21.7 (Patat 2008) (225 μcd m^{−2}), which means that the satellites from Table 1 will contribute at most an additional ∼1% to the sky brightness in the worst areas of the sky. The contribution to the sky brightness is therefore small and the simulations and mitigation focus on the discrete contamination by individual satellites.
Fig. 11.
Sky brightness contribution from the satellites using the constellations from Table 1. The results are computed at astronomical twilight for latitudes +50° (a) and −25° (b). Typical sky surface brightness in the visible is 21.7 mag arcsec^{−2} or brighter. The satellites used all have v_{1000 km} = 7, resulting in visual magnitudes in the range indicated at the bottom right corner. 
Kocifaj et al. (2021) have evaluated the increase in diffuse sky brightness caused by all current space objects with sized between 5 × 10^{−7} m to 5 m at altitudes above 200 km. They estimate that this excess is 16.2 μcd m^{−2} (24.6 mag arcsec^{−2}) and can reach 21.1 μcd m^{−2} (24.3 mag arcsec^{−2}) at astronomical twilight, corresponding to about 10% of the natural sky brightness. That excess is dominated by the small objects (mm and below), what is usually called space debris. The macroscopic satellites composing the constellations discussed in this paper will not contribute much to the diffuse sky brightness provided they are not ground into microscopic debris.
3.4. Effect on observations
To evaluate in more detail the effects of satellite constellations on observations, we studied a series of representative instruments and telescopes. For each of them, we consider the field of view (size or diameter in case of 2D field, length and width for a slit, diameter of the aperture in case of a fibre), a typical exposure time, and the limiting magnitude (obtained from the Exposure Time Calculators^{3} for ESO instruments, documentation or private communications for others). We also estimate the magnitude causing heavy saturation either as 5 mag brighter (100 times) than the saturation level or from publications. Table 3 lists the parameters of the exposures and instruments.
Characteristics of instruments and exposures used for the simulations.
For each constellation shell, the instantaneous satellite density, angular velocity and apparent and effective magnitudes were estimated. The number of trails affecting an exposure was obtained using Eq. (1), using A_{fov} = L_{1} × L_{2} for the field of view (with L_{1} the length and L_{2} the width of the FOV, L_{1} > L_{2}) and L_{fov} = L_{1} in the second term of Eq. (1) – this maximises the crosssection for trails. N_{sat} was computed for each shell accounting for the effective magnitude of the satellites in that shell.
The effect on the observations is computed as follows: Those satellites with an effective magnitude fainter than the 1σ detection limit were ignored, considering that their trail would be lost in the background noise. Those between the detection limit and heavy saturation limit were counted, and each one was considered to ruin a 5″wide trail across the whole detector, a size acceptable as a good general value as can be deduced from Eq. (8). In case of a long slit, they ruin 5″ of the slit. In real observations, it is plausible that all or part of the data below a nonsaturated trail could be recovered, so this is a pessimistic limit. In the case of a fibre contaminated by a satellite, we consider that the whole spectrum is lost. For trails brighter than the heavy saturation limit, the whole exposure is considered damaged by the charge bleeding and/or electronic and/or optical ghosts. This was repeated for each shell in the constellation, and the effects were summed, resulting in maps of lost fractions. A value of, say, 50% indicates that either 50% of the individual exposures are entirely lost, or that 50% of the pixels in each exposure are lost or, more likely, a combination in between. This was then repeated for several solar elevations ranging from twilight to midnight. Figure 12 displays the resulting sky maps of trail count and fraction lost for an example, with instrument specific all sky plots provided for imagers in Fig. 13 and spectrographs in Fig. 14. In this figure, VRO’s SST Cam appears much more affected than other wide field instruments such as VST or G96. This reflex the fact that the brightest satellites will heavily saturate the SST Cam detector, causing the full exposure to be lost; similar cameras on smaller telescopes do not suffer from that effect. For the spectrographs, the 4MOST (low resolution mode) multifibre spectrograph will be crossed by a large number of trails, but the overall effect is rather small: this is because each trail will affect on average 1.3 fibre (based on MonteCarlo simulations of the field of view) on a total of 2436 fibres.
Fig. 12.
Sky maps of the number of detectable satellite trails (a) and effect on the observations (b), for all the satellites from Table 1 on a FORS2 image (6′ field of view, 5 min exposure time) on Paranal ( latitude) at equinox. The circles indicate elevations 0°, 10°, 30°, and 60°. The legend of each plot gives the Sun elevation and the average number of trails (a) and the losses they cause (b) for observations above 30° elevation All satellites are brighter than the detection limit and none is bright enough to cause heavy saturation. 
Fig. 13.
Sky maps of the number of detectable satellite trails in an exposure (left) and effect on the observations (right) for a series of imagers (see Table 3 for their characteristics). The legend of each plot also lists the average number of trails above 30° elevation. The Sun declination is 0° and its elevation −20°. 
Fig. 14.
Sky maps of the number of detectable satellite trails in an exposure (left) and effect on the observations (right) for a series of spectrographs (see Table 3 for their characteristics). The legend of each plot also lists the average number of trails above 30° elevation. The Sun declination is 0° and its elevation −20° (−18° for 4MOSTHiRes; the values for −20° are 0 trail and 0%). 
The average number of trails per exposure and the average fraction of the exposure lost were computed for the region of the sky above 30° elevation by integrating the results over that region of the sky. These averages are shown in Fig. 15.
Fig. 15.
Average number of trails per exposure (a) and average fraction of the exposure lost (b) as a function of the sun elevation, for representative exposures at elevation > 30° on the instruments listed in Table 3. Twilights are shaded in blue; inaccessible solar elevations are shaded in grey for the equinoxes and solstices for Paranal latitude (). 
Our software to predict the effect of satellites on observations is available online^{4}.
4. Discussion
As apparent from Fig. 15 and expected from Eq. (1), the number of trails in an exposure increases with the size L_{fov} of the fieldofview and with exposure time t_{exp}. The effect of the effective magnitude (Eq. (7)) is less intuitive: for the altitudes of the constellations considered in this study, m_{1000 km} and exposure times, the effective magnitudes fall in the range 13–23 and become fainter as the exposure time increases, with a linear dependency in t_{exp}. As the limiting magnitude for an exposure goes fainter with a dependency in (considering the simple skynoise dominated case) there is, for each instrument, an exposure time beyond which a satellite trail will no longer be detectable. In other words, the contribution from a satellite to the intensity in a resolution element is independent of the exposure time, while the noise increases with . Overall, the S/N of the satellite trail decreases with ; if an exposure t_{exp} is immune to satellite trails, longer exposures will also be immune.
Imagers on all but the smallest telescopes typically have a limiting magnitude fainter than the faintest satellite effective magnitude: they are, therefore, affected to some extent by all satellite constellations. For many science cases, the presence of a trail will only result in a loss of useful imaged area (of the order of 0.1 to 1% for a 5″ wide trail crossing a 1° or 8′ field of view). There will be, however, some science cases in which even a faint trail will ruin the whole exposure (such as photometry of a faint transNeptunian object overrun by a satellite), leaving no other choice than repeating the exposure, if this were possible at all (sometimes the repetition is not possible as, for instance, for the photometry of a transient gammaray burst). Furthermore, for the most sensitive cameras, some satellites have an effective magnitude brighter than the heavy saturation limit, wreaking havoc in the affected exposures, as on the LSST camera at the Vera Rubin Observatory (VRO), and resulting in much heavier losses (Tyson et al. 2020).
For astrophotography widefield cameras, the limiting magnitude for satellite trails scales inversely with the focal length of the lens, with every other parameter remaining constant. A wideangle camera with 30 mm focal length will therefore be 5 mag less sensitive than a 3 m focal length telescope with the same focal ratio. As a consequence, astrophotography will be immune to most satellites in their operational orbits. They can, however, be affected by brighter satellites such as larger satellites, or telecommunication satellites in low altitude transfer orbits (such as the bright stringsofpearls of 60 very bright satellites observed after the early Starlink launches), or specular reflections. Fortunately, these are not numerous: it is foreseen that there will be of the order of 10 trains of satellites around the Earth at any time to replenish the constellations. While potentially spectacularly damaging, these are statistically unlikely and visible only during the brightest parts of twilights.
For spectrographs, the limiting magnitude for a single exposure often falls in the range of the satellites effective magnitudes. As a consequence, those fainter than the limit are not detected and only slightly contribute to the background noise. This is the case for all satellites for highresolution spectrographs or échelle spectrographs, even on very large telescopes (see the examples of UVES and ESPRESSO on the VLT). However, low to mediumresolution spectrographs on medium to large telescopes will detect all or many satellites. Furthermore, contrary to imagers, where a satellite leaves a telltale trail in the data, slit and fibre spectrographs do not record spatial information. While highS/N contamination would be easy to notice (with the exposure level is much higher than expected and the spectral shape does not match that expected for the target), many satellites will leave a signal with a low to moderate S/N. In many cases, the contamination will be at a level comparable to or below that of the science target and therefore unlikely to be detected in realtime. Unless the contamination is flagged using other means (see below), there will be cases for which it will become apparent only at the time of the scientific analysis of the spectra. As the contamination will have a solar spectrum, some science cases will be better protected (the study of distant quasars, for example) than others (such as study of double stars, where a solar spectrum may not be surprising).
In the thermal IR domain, the overall signal is dominated by the very strong thermal emission from the sky and the telescope. The individual exposure time is therefore kept extremely short (few tens of milliseconds) and the background is registered by chopping (performing a small position offset by tilting the secondary mirror of the telescope) at about 1 Hz and nodding (another small offset by moving the whole telescope) every few seconds. In Hainaut & Williams (2020), we conservatively estimated the flux from a satellite at 2000 km at zenith up to 100 Jy in Nband (8–13 μm) and up to 50 Jy in the M and Qbands (5 and 18–20 μm, respectively). The variations between an illuminated and a shadowed satellite are negligible. These fluxes are well above the detection threshold of the thermal IR instruments in Table 3, even accounting for trailing. Because of the extremely short exposure time and small field of view, on average only 6 × 10^{−6} trails would be found in a single exposure. However, for most observations, the images are not individually recorded but averaged over a nodding cycle. While the S/N of the trail will be washed away by this average, the values plotted in Fig. 15 correspond to the duration of these averages (10 s). In the case of VISIR on the VLT, it is considered that a satellite ruins a 5″ trail across the detector and for METIS on the larger ELT, the full (smaller) image is ruined by the (broader) trail. Even with these extremely pessimistic assumptions, a negligible fraction of the thermal IR data is affected. For spectroscopy, even at low resolution, the satellite effective fluxes will be below the detection limit.
The case of stellar occultations was discussed in Hainaut & Williams (2020). The effect was found to be small: a 0.02 to 10 millimag for 10 s and 0.1 s exposure time and extremely improbable: 10^{−4} to 10^{−6} exposures would be affected. The simulations presented here do not change these estimates. As the eclipse by a satellite would affect only one measurement in a series, even if it could be measured, it would not be similar to the occultation by an exoplanet or a transNeptunian object.
The case of visual observations – either naked eye, or through binoculars or telescope – will be considered in a separate paper. In summary, 15 to 50 satellites would be visible in the sky with the naked eye when the sun elevation is between −12° and −24°. When the sun is higher than −12°, the sky is too bright and no satellite is visible. When the sun is lower than −24°, no satellite is bright enough to be visible.
First order of mitigation refers to the satellites themselves. As seen in Fig. 10, the number of illuminated satellites in sight is of course a function of the number of satellites in the constellations, but also of the altitude of the constellation. Furthermore, highaltitude constellations remain illuminated by the Sun much longer than lowaltitude ones. The apparent magnitude of a highaltitude satellite will be fainter than that of the same satellite on a lower altitude (Eq. (6)), which is advantageous for small telescopes. However, for large telescopes, the satellite appears extended (Eq. (8)), so that its surface brightness will not decrease much. Overall, a constellation at 1000 km will be more damaging than a three times larger constellation at 500 km altitude.
Obviously, keeping the brightness of all satellites below the detection limit of all telescopes would be ideal. With the increasing size of the telescopes and sensitivity of the instruments, this is not realistic. An achievable goal would be to reduce the effective crosssection of the satellite so that they always remain below the saturation threshold of the most sensitive instrument. Today and in the foreseeable future, this threshold is set by the SST camera at VRO, at V_{550 km} > 7 (Tyson et al. 2020), or V_{1000 km} > 8.3. The changes introduced by SpaceX^{5} with VisorSat and modified attitude of the solar array are very promising steps in the right direction: most of the satellites are now below the heavy saturation threshold, while still causing electronic cross talk (Tyson et al. 2020). More systematic measurements of the satellites are needed (as suggested by SATCON1 Recommendation 8, see Walker et al. 2020) and awareness of the satellite operators is a must. SATCON1 Recommendation 5 and DARK & QUIET SKIES Recommendation 15^{6} formalise the brightness limit at V_{550 km} > 7.
Once the satellites are in orbit, the next level of mitigation is at the time of preparation and scheduling of the observations. Because of the progression of the Earth shadow through the constellation shells and of the fine structure in the apparent density of satellites (Sect. 3.2), the fraction of losses can change dramatically by pointing the telescope in a slightly different direction. For a given time at a given observatory, sky maps such as those in Fig. 12 would allow an observer to pick objects in the region of the sky that are least affected by satellites. As these maps can be generated onthefly, they could be integrated, for instance, in a queue observation optimisation algorithm. Another way to consider the scheduling is to pick the best time slots to observe a given object with a specific instrument. Using the same methodology as for the sky maps, a calendar can be populated with the expected density of satellites for an object seen from an observatory. Such calendars, some examples of which are displayed in Fig. 16, could be used when allocating specific telescope time to an observation programme in traditional visitor mode. Obviously, both the sky maps and calendars will introduce additional complexity in the scheduling process, and will not resolve all issues. For instance, some observations must be performed at a given time (for instance an exoplanetary transit observation).
Fig. 16.
Two examples of calendars showing the visibility of an object from an observatory. Top: galactic Centre from the VLT. Bottom: large Magellanic Cloud from VRO. The fraction of observing time lost due to satellites is indicated by the colour scale (1 indicating that all exposures are damaged), for 300 s exposures with FORS2 and 15 s on the LSST camera, respectively. The satellites are the 60+ thousand from Table 1. The blue shading marks daytime and the blue contours indicate the twilights. The elevation of the object is indicated by the greenish contour lines and the grey shading indicates the times when the object is below 20° elevation. 
A more aggressive mitigation would be to close the shutter of the instrument just before a damaging satellite enters the field of view and reopen it just after it left. As the satellites move at apparent speeds of ∼0.1 to ∼1° s^{−1}, the interruption would be extremely short, virtually nullifying the losses of exposure time. The challenge is to send the signal to the shutter at the right time. We discuss below the two methods that can be envisioned.
The first would rely on a complete, accurate and uptodate database of orbital elements so that the position of all satellites can be computed at any time and offending ones identified in times. However, the accuracy of ephemerides, both in position and in timing, must be of the order of the fieldofview. An accuracy of a fraction of a degree (and 1 s) may be achievable, making this viable for imagers. An accuracy of ∼1″, required for spectrographs, would imply predicting the position of the satellites at 2–5 m accuracy, with a timing precision of 1/100 s. This method presents various challenges: the database must be complete; it must be uptodate (as nonkeplerian effects, including active orbital corrections, modify the orbit with a timescale of up to a few days); it must be precise (the current standard TwoLine Elements do not provide the required precision); the computations must be done with the appropriate precision, and scan the whole database for each observation. Furthermore, this method would not work for all instruments: large survey cameras tend to be fairly heavy and slow, precluding rapid and repetitive shutting and opening.
The second method would rely on an auxiliary camera mounted in parallel with the main telescope, with a field of view of 5 to 10° (a few degrees larger than the field of the main instrument). This camera would take an image of the field about every second and an analysis system would detect any transient object. An object moving towards the science field of view would trigger the closing of the shutter. The challenge here is to build a system fast enough to process the data in realtime and a camera sensitive enough to detect the satellites. For instruments sensitive to all the satellites down to the faintest (for example, imagers on large telescopes), it may not be realistically feasible. However, for spectrographs, which have a brighter limiting magnitude, a 30 cm auxiliary telescope with a fast readout detector is promising. Spectrographs would strongly benefit from this mitigation, as their exposures tend to be much longer than those of imagers (and therefore the loss of an exposure more costly) and because they lack spatial information that makes it possible that the contamination will remain unnoticed until the data are analysed.
The final stage of mitigation is an a posteriori subtraction of the satellite trail from the data. This also comes with some limitations. Saturated trails, such as those on a survey camera on a large telescope like the SST Cam at VRO, are unrecoverable. Fainter trails may be identified on images, modeled and subtracted. However, atmospheric scintillation will make these trails irregular and difficult to model. Even assuming that they can be cleanly subtracted (and that this subtraction is trusted by the scientist), they will result in an increase of the photon noise. It may be safer to mask them and filter them out by combining several exposures of the same field. This is already systematically done for various types of blemishes affecting the images, such as cosmic ray hits and gaps between chips in the detector mosaic, and would result in a wellquantified loss of total exposure time in the affected strips, as already taken into account for detector gaps.
The case of slit or fibre spectroscopy is again more difficult: no telltale trail reveals the passage of a satellite, which can only be detected as an additional solartype spectrum added to the data. A solar spectrum can be iteratively subtracted from the spectrum until the residual shows no hint of the satellite. While this would work for some science cases, it will be more difficult or even impossible in other situations (for instance a programme studying stellar abundances or binary stars).
Finally, the last resort of mitigation consists in repeating the observations that have been damaged by a satellite. In some cases, this will be immediately obvious and easily detected when controlling the quality of the data. In other cases, the effect will be more subtle. Some observations will be simple to reacquire. Others will be lost forever (a short transient phenomenon, like the optical counterpart of a gravitational wave event).
Overall, no mitigation method will singlehandedly work for all instruments and all science cases. Moreover, each of these mitigations comes with a cost that should be carefully compared to the cost of the observing time loss: in many cases, repeating the observation may be cheaper (economically and scientifically) than protecting it.
The way forward is to improve the situation at each step, starting with a collaboration with the satellite operators to make the satellite less bright, continuing with smarter scheduling of the observations, where the work presented in this paper will help thanks to the fast and accurate information it provides, shutter control where possible and accounting for the inevitable remaining trails in the data.
For example, SpaceX submitted an amendment to the orbital configurations of their Starlink Generation 2 constellation in August 2021 (https://fcc.report/IBFS/SATAMD2021081800105/12943361).
See p. 153 of the report: https://www.iau.org/static/publications/dqskiesbook291220.pdf
Acknowledgments
We thank the referee for the careful review and constructive suggestions. This work originated as part of the SATCON1 and DARK & QUIET SKIESworkshops. We thank the organisers of these workshops.
References
 Appenzeller, I., Fricke, K., Fürtig, W., et al. 1998, The Messenger, 94, 1 [NASA ADS] [Google Scholar]
 Buzzoni, B., Delabre, B., Dekker, H., et al. 1984, The Messenger, 38, 9 [NASA ADS] [Google Scholar]
 de Jong, R. S., Barden, S. C., BellidoTirado, O., et al. 2016, in Groundbased and Airborne Instrumentation for Astronomy VI, eds. C. J. Evans, L. Simard, & H. Takami, SPIE Conf. Ser., 9908, 99081O [Google Scholar]
 Dekker, H., D’Odorico, S., Kaufer, A., Delabre, B., & Kotzlowski, H. 2000, in Optical and IR Telescope Instrumentation and Detectors, eds. M. Iye, & A. F. Moorwood, SPIE Conf. Ser., 4008, 534 [Google Scholar]
 Hainaut, O. R., & Williams, A. P. 2020, A&A, 636, A121 [EDP Sciences] [Google Scholar]
 James, N. D. 1998, J. Br. Astron. Assoc., 108, 187 [Google Scholar]
 Karttunen, H., Kröger, P., Oja, H., Poutanen, M., & Donner, K. J. 1996, Fundamental Astronomy [Google Scholar]
 KisslerPatig, M., Pirard, J. F., Casali, M., et al. 2008, A&A, 491, 941 [Google Scholar]
 Kocifaj, M., Kundracik, F., Barentine, J. C., & Bará, S. 2021, MNRAS, 504, L40 [Google Scholar]
 Kuijken, K., Bender, R., Cappellaro, E., et al. 2002, The Messenger, 110, 15 [NASA ADS] [Google Scholar]
 Lagage, P. O., Pel, J. W., Authier, M., et al. 2004, The Messenger, 117, 12 [NASA ADS] [Google Scholar]
 Mallama, A. 2020a, ArXiv eprints [arXiv:2003.07805] [Google Scholar]
 Mallama, A. 2020b, ArXiv eprints [arXiv:2006.08422] [Google Scholar]
 Mallama, A. 2020c, ArXiv eprints [arXiv:2012.05100] [Google Scholar]
 Mallama, A. 2021, ArXiv eprints [arXiv:2101.00374] [Google Scholar]
 McDowell, J. C. 2020, ApJ, 892, L36 [Google Scholar]
 Montenbruck, O., & Gill, E. 2000, Satellite Orbits: Models, Methods, and Applications, Physics and astronomy online library (Berlin Heidelberg: Springer) [Google Scholar]
 Otarola, A., Allen, L., Pearce, E., et al. 2020, in United Nations Office for Outer Space Affairs, Dark and Quiet Skies for Science and Society (published by the International Astronomical Union), https://www.iau.org/static/publications/dqskiesbook291220.pdf [Google Scholar]
 Patat, F. 2008, A&A, 481, 575 [Google Scholar]
 Patat, F., Moehler, S., O’Brien, K., et al. 2011, A&A, 527, A91 [Google Scholar]
 Pepe, F., Cristiani, S., Rebolo, R., et al. 2021, A&A, 645, A96 [Google Scholar]
 Ragazzoni, R. 2020, PASP, 132, 114502 [Google Scholar]
 TregloanReed, J., Otarola, A., Ortiz, E., et al. 2020, A&A, 637, L1 [EDP Sciences] [Google Scholar]
 Tyson, J. A., Ivezić, Ž., Bradshaw, A., et al. 2020, AJ, 160, 226 [Google Scholar]
 Vallado, D., Crawford, P., Hujsak, R., Kelso, T., et al. 2006, AIAA/AAS Astrodynamics Specialist Conference and Exhibit [Google Scholar]
 Walker, J. G. 1984, J. Br. Interplanet. Soc., 37, 559 [NASA ADS] [Google Scholar]
 Walker, C., Hall, J., Allen, L., et al. 2020, Bull. Am. Astron. Soc., 52, 0206 [Google Scholar]
 Witze, A. 2019, Nature, 575, 268 [Google Scholar]
Appendix A: Derivation of key expressions
In the following paragraphs we first show the derivation of Eqn. 2, the probability density function P(ϕ, i, h_{sat}) (§A.1). We later get into some geometric details on the relation between the geocentric and topocentric positions of the satellites, impact angle of the lineofsight on an orbital shell, and the apparent and observed angular velocity of satellites (§A.2, §A.3 and §A.4). Finally, in §A.5, we provide some equations on how to locate the intersection of constellation shell boundaries on the sky with the local meridian (Eqn. 9 and some related expressions).
A.1. Probability density
The analytical probability density function P(ϕ, i, h_{sat}) displayed in Eqn. 2 is at the core of the analytical simulations included in this paper. We derive this probability density function from first principles and we also prove that its integral is equal to unity.
For one single satellite in a circular orbit with inclination i and orbital altitude h_{sat} we first obtain the angle β between the orbit and any parallel, as a function of latitude ϕ. At the nodes the parallel is the equator and we have β = i. The general situation is depicted in Fig. A.1. In that diagram Q + i = 90° and Q′ = 90° +β. Our problem is finding Q′ as a function of i and ϕ. This is solved via the law of sines:
Fig. A.1.
The spherical triangle to obtain the angle β between the orbit and the parallel of latitude ϕ. The angle κ indicates the anomaly of the satellite measured from the ascending node along its orbit. With λ we indicate the difference between the geocentric longitude of the satellite and that of the ascending node. 
The only satellite in our orbit induces a linear density of satellites per angular unit along the orbit itself that is Λ = 1/(2π). We have assumed that the satellite is uniformly distributed in time along its trajectory, which is true for a circular orbit. Now we have to transform this angular linear density (satellites per unit angle along the orbit) into an angular surface density (satellites per unit solid angle on the orbital shell).
We note that for a single orbit, even considering that its only satellite is uniformly distributed along the orbit with linear density Λ, this will not induce a uniform surface density over the shell, neither in longitude nor in latitude. In longitude, the density will be higher around the nodes and lower at the longitudes corresponding to extreme latitudinal excursions (at 90° from any node). In latitude we have more or less the opposite: the latitudinal density will be higher where the orbit is tangent to parallels (extreme latitude excursions) and lower at the equator. In any case, it is clear that the density will be zero outside the range −i < ϕ < i.
For an orbital shell where the orbital planes are uniformly distributed in longitude, with the satellites distributed uniformly within each orbital plane, the density distribution with longitude will smooth out due to the averaging over random initialisations. Hence, we conclude and impose that the probability density function is uniform in longitude.
For a band in latitude of infinitesimal angular width dϕ at latitude ϕ, as shown in Fig. A.2, we find that:
Fig. A.2.
Transformation from angular linear density to angular latitudinal density. 
We introduce the density Λ to pass from arc κ to satellites per unit arc in latitude p = Λκ:
Next, we obtain the expression for the latitudinal angular probability density due to one single orbit with one single satellite:
But, as shown in Eqn. A.3, cos β = cos i/cos ϕ and, thus:
This density takes into account only one side of the orbit, but each parallel is crossed twice by the same orbit, so that the final true latitudinal angular density is twice this value, p_{2} = 2p:
This latitudinal angular density would be measured in fractions of the sample per unit angle in latitude. We now divide by the length of a parallel at latitude ϕ, 2π cos ϕ, to transform this angular latitudinal density into true density per unit solid angle and label the resulting probability density function as c(ϕ, i):
This density per unit solid angle can be converted into the true, physical surface density, in units of fractions of the sample per surface unit on the shell, by dividing by the radius of the shell squared, (R_{⊕} + h_{sat})^{2}. Also, by using that cos^{2}ϕ − cos^{2}i = sin^{2}i − sin^{2}ϕ and that both expressions are differences of squares, we end up with four equivalent formulations of our probability surface density function:
We use the last one of the four options in our computations.
Function P(ϕ, i, h_{sat}) is a probability density distribution and its integral over the space occupied by the sample is unity, as stated in the main text of this article. This may be shown from any of the four forms of the function, but is easier from the first one, Eqn. A.11. Noting R_{sat} = R_{⊕} + h_{sat}, we have:
This last integral is reduced to a straightforward arcsine with the variable change x = sin ϕ:
A.2. Geocentric and topocentric position of the satellites
The number of trails affecting an exposure, given by Eqn. 1 in §2.1.2, relies on the apparent angular velocity of the satellites at that position in the sky. To compute that velocity, we first have to relate the geocentric longitude and latitude of the satellite (θ and ϕ) to the topocentric position vector of a satellite (pointing direction of the telescope, given by OS), which is characterised by its right ascension (more precisely, the hour angle) and declination, or equivalently by its azimuth and elevation.
Writing the vectorial relation between the centre of Earth C, the observer O and the satellite S,
in geocentric rectangular coordinates (x equatorial at the meridian of the observer, z to the pole, and y completing the reference frame) using the longitude and latitude of the observer (θ_{o}, ϕ_{0}) and of the satellite (θ, ϕ), Eqn. A.19 becomes
where R_{⊕} is the radius of the Earth, R_{sat} = R_{⊕} + h_{sat} is the radius of the satellite orbit, d = OS, and (x_{OS}, y_{OS}, z_{OS}) the unit vector, so that OS = d(x_{OS},y_{OS},z_{OS}). The coordinates are obtained from the hour angle and declination of the satellite. We eliminate θ and ϕ by summing quadratically these three equations:
Solving this quadratic equation leads to two solutions for d. The positive one is the distance to the satellite (the negative one is the distance to the satellite shell in the opposite direction, below ground).
Once d is computed, Eqn. A.22 can be solved for ϕ and then Eqn. A.20 and A.21 for θ. We have now the geocentric longitude and latitude (θ and ϕ) of the satellite that corresponds to a specific, topocentric pointing direction on the local sky.
A.3. Impact angle of the line of sight with the shell
This is the angle α between the line of sight and the normal to the shell, that is, the angle . From the cosine law in triangle CSO, we have:
A.4. Apparent angular velocity
A satellite observed in a given detection may either be on the northbound half of its orbit, or on the southbound half. Considering the right spherical triangle formed by the satellite, the ascending node of its orbit and the equator (see Fig. A.1), we obtain the longitude difference λ between the longitude of the satellite and that of the ascending nodes from:
from which we get the longitudes of the ascending nodes for both possible orbits:
The geocentric velocity vector of a satellite is obtained from:
where v is obtained from elementary celestial mechanics as , N is the unit vector normal to the orbit (built from the longitude of the ascending node and i), R = CS/CS is the unit vector pointing from the centre of the Earth to the satellite and × marks the cross product. This is repeated for both Ω_{N} and Ω_{S}.
The topocentric velocity vectors are obtained by subtracting the geocentric velocity vector of the observer from that of the satellite. The components of these vectors perpendicular to the line of sight OS are the two apparent angular speeds ω_{N} and ω_{S}.
A.5. Declination of shell boundaries
Equation 9 relates observatory latitude ϕ and shell inclination i to the declination δ at which the shell boundary cuts the local meridian. It is deduced from elementary geometry (plane sine theorem) applied to triangle COS in Fig. A.3:
Fig. A.3.
Geometry of shell boundaries at the local meridian. C: Earth centre. O: Observatory. S: Intersection of the northern shell boundary with the local meridian. 
The same elementary procedure applied to other combinations of sides and angles of triangle COS leads to two additionalrelations:
To use these two relations, it is necessary to introduce the distance d_{sat} from the observatory to the shell boundary at the meridian, a quantity that is also easily deduced from Fig A.3 through the plane cosine theorem:
Appendix B: Index of symbols
This is a list of the symbols used in this paper, together with a short definition.

α: impact angle between the lineofsight and the spherical shell at the satellite

α_{⊙}: solar phase angle, angle SunSatelliteobserver

χ: airmass

δ: declination angle

δ_{sat}: density of satellites in a field of view, in satellites per unit solid angle (n/sq.deg).

δ_{trail}: density of satellite trails crossing a the field of view, in satellite per linear degree per unit of time (deg^{−1} s^{−1})

κ: nodal anomaly of the satellite, that is, the angular longitude of the satellite measured from the ascending node along the orbit

Λ: linear density of satellites per unit angle along the orbit

λ: difference of longitudes of the satellite and the ascending node of its orbit (longitude of the satellite measured from the ascending node)

ω_{sat}: apparent angular velocity of the satellite as seen by the observer

ρ(ϕ, i, h_{sat}) = N_{sat}b: density of satellites at latitude ϕ

θ, ϕ: geocentric longitude and latitude of the satellite

P(ϕ, i, h_{sat}): probability density function of finding a satellite (with i, h_{sat}) at latitude ϕ

d, d_{sat}: distance between the observer and the satellite

d_{sat⊙}: distance from the Sun to the satellite, ∼1 AU

h_{sat}: Altitude of the satellites’ orbit, km

i: orbital inclination of the satellites’ orbit, degrees

k: atmospheric extinction coefficient, mag/airmass

L_{fov}: diameter of the (circular) field of view, degrees

l: latitude of the observer

m_{1000 km}, m_{500 km}: zenithal magnitude of the satellite normalised to a distance d_{sat} = 1000 km, 500 km, and zero airmass (no atmospheric extinction)

m_{eff}: effective magnitude of the satellite: the magnitude of a static, pointlike object that, during the exposure time considered, would produce the same accumulated intensity in a resolution element than the satellite crossing over a resolution element

m_{sat}: (visual) magnitude of the satellite

m_{⊙}: magnitude of the Sun

n_{1}: number of satellites in a single orbital plane

n_{plane}: number of orbital planes in the constellation shell

n_{sat}: number of satellites present in the field of view

n_{trail}: number of satellite trails crossing the field of view during the exposure

N_{sat}: total number of satellites in the constellation shell

p: geometric albedo of the satellite

p(S): probability of finding a satellite in region S

r_{inst}: angular size in the plane of sky of the resolution element or pixel size of a detector

R_{⊕}: radius of the spherical Earth,

R_{sat}: radius of the (spherical) satellite or, in other contexts, radius of the (circular) orbit of a satellite

t_{exp}: exposure time, in seconds

t_{eff}: effective exposure time for a satellite, the duration it takes the satellite to cross a resolution element of the detector
All Tables
All Figures
Fig. 1.
Discrete realisation of the Starlink Generation 2 constellation overlaid on a map of Earth. The day and night sides of Earth are shown for 2022 June 21 at 23:00UTC, near the June solstice. Satellites are shown as either sunlit (yellow), or eclipsed by the Earth (black) and above the horizon for two geographical locations, Vera Rubin Observatory in Chile (latitude −30°) and London in the United Kingdom (latitude ∼50°). Overdensities of satellites in geographic latitude b are visible at latitudes b = 30°, 40° and 53°, corresponding to the inclination of the most populous orbital shells. Due to their orbital altitudes, satellites remain visible for locations where the Sun is below the horizon and will remain visible throughout all the night for locations at high geographical latitudes. 

In the text 
Fig. 2.
Discrete realisation of the Starlink Generation 2 constellation from Fig. 1. The satellites are plotted on all sky maps representing the celestial sphere above Vera Rubin Observatory (left) and London (right). Satellites are shown as sunlit (yellow) or eclipsed by the Earth (black). The number of satellites increases towards the horizon because, for a given solid angle, looking at lower elevations from the horizon the distance to the shell surface increases and, hence, the field of view contains larger portions of the orbital shells. Due to the orbital inclination and orbital altitude of the orbital shells, overdensities of satellites are also present in these maps. The map for London shows this very clearly for the two orbital shells with i = 53° and h = 345 km (7178 satellites) and 499 km (4000 satellites), where more satellites are present south from the dotted line. The dotted line in the all sky map for Vera Rubin Observatory is for the orbital shell with i = 30°, h = 328 km, with a larger number of satellites north from that line. 

In the text 
Fig. 3.
Number of satellites visible in an exposure depending on the exposure time t_{exp} and the instrument fieldofview, specified by its circular radius R_{fov} = L/2. The effect of a discrete constellation of 10 000 satellites in a single orbital shell of 100 orbital planes and 100 satellites per plane at 1000 km altitude and 53° inclination is simulated for an observer at −30° latitude, observing towards zenith. The points are the results of the discrete simulations for different values of exposure time t_{exp} and fieldofview (). The solid lines are predictions based on Eq. (1). 

In the text 
Fig. 4.
Sky maps with an example of the resulting number of trails per exposures. The circles mark 0°, 10°, 20°, 30°, and 60° elevation. The observatory is located at +50° latitude; the sun at −18° elevation; the camera has a FOV with diameter L_{fov} = 6′ and the exposure time is 300 s. The satellites are those from Table 1. In the black region at the southeast, all satellites are already in the shadow of the Earth. The edges running from northeast to southwest correspond to constellation shells at lower altitudes, whose southeast parts are already in the shadow. The sharp features running from east to west correspond to the edges of constellation shells whose inclinations are close to the observatory latitude. 

In the text 
Fig. 5.
Magnitude of the satellites as a function of the range between the observer and the satellite. The dots are measurements of original Starlink satellites (Pomenis telescope, from Otarola et al. 2020), and the lines are obtained using a simplified Lambertian sphere model for two altitudes. The unknown solar phase angle contributes to the dispersion of the measurements. 

In the text 
Fig. 6.
Visual (red) and effective magnitudes of a satellite at zenith, as a function of its altitude, for various exposure times (see legend). The satellite used is a Starlink VisorSat with m_{1000 km} = 7, considered as a trailed point source. 

In the text 
Fig. 7.
Detection limit for the satellite apparent magnitude (coloured lines) as a function of the limiting magnitude of the instrument (both are at the same S/N). Panel a: results for 60 s exposures with an imager, while (b): prediction for 1200 s exposures with a spectrograph. The trailing of a satellite below the coloured line will make its effective magnitude too faint for detection. The shaded areas highlight the range where constellation satellites are likely to be detected. Two elevations (zenith and airmass 2) are considered, with the corresponding typical satellite visual magnitudes and angular velocities. 

In the text 
Fig. 8.
Limiting magnitude for an astrophotography camera with wideangle lens (red) as a function of the exposure time. Dark current noise is neglected (see text for the camera details). This magnitude is converted into the apparent magnitude of a satellite assuming the angular velocity typical of a Starlink at 30° elevation), using Eq. (7). 

In the text 
Fig. 9.
Number of illuminated satellites above the horizon (top panel) and above 30° elevation, as a function of the local solar time, for Paranal (latitude ). The figure accounts for all the satellites from Table 1. Left is for the summer solstice () and right for the winter solstice (). The twilights are indicated with blue shading. The black line marks the total number of satellites above the elevation considered, the blue line those that are illuminated and the orange and red lines those brighter than magnitude 6 and 5, respectively, using the photometric model described in Sect. 2.2. This shows that during the summer months, a few hundred satellites would remain illuminated thorough the whole night, including a few tens above 30° elevation. 

In the text 
Fig. 10.
Number of illuminated satellites above the horizon as a function of the Sun elevation, for the constellations listed in Table 1 seen from Paranal (latitude ; the dependency with latitude is not strong). The twilights are indicated with blue shadings and the elevation of the sun at midnight by grey shadings for the solstices and equinoxes. The thin dotted lines represent the individual shells and the thick lines the totals for each constellation. The upper thick black line is the grand total. This figure illustrates the importance of the altitude of a constellation: the highaltitude satellites (OneWeb, Kuiper) are visible much longer into the night than lower altitude ones. Also, for most of the night, the much smaller OneWeb constellation contributes to more satellites in the sky than the largest one. 

In the text 
Fig. 11.
Sky brightness contribution from the satellites using the constellations from Table 1. The results are computed at astronomical twilight for latitudes +50° (a) and −25° (b). Typical sky surface brightness in the visible is 21.7 mag arcsec^{−2} or brighter. The satellites used all have v_{1000 km} = 7, resulting in visual magnitudes in the range indicated at the bottom right corner. 

In the text 
Fig. 12.
Sky maps of the number of detectable satellite trails (a) and effect on the observations (b), for all the satellites from Table 1 on a FORS2 image (6′ field of view, 5 min exposure time) on Paranal ( latitude) at equinox. The circles indicate elevations 0°, 10°, 30°, and 60°. The legend of each plot gives the Sun elevation and the average number of trails (a) and the losses they cause (b) for observations above 30° elevation All satellites are brighter than the detection limit and none is bright enough to cause heavy saturation. 

In the text 
Fig. 13.
Sky maps of the number of detectable satellite trails in an exposure (left) and effect on the observations (right) for a series of imagers (see Table 3 for their characteristics). The legend of each plot also lists the average number of trails above 30° elevation. The Sun declination is 0° and its elevation −20°. 

In the text 
Fig. 14.
Sky maps of the number of detectable satellite trails in an exposure (left) and effect on the observations (right) for a series of spectrographs (see Table 3 for their characteristics). The legend of each plot also lists the average number of trails above 30° elevation. The Sun declination is 0° and its elevation −20° (−18° for 4MOSTHiRes; the values for −20° are 0 trail and 0%). 

In the text 
Fig. 15.
Average number of trails per exposure (a) and average fraction of the exposure lost (b) as a function of the sun elevation, for representative exposures at elevation > 30° on the instruments listed in Table 3. Twilights are shaded in blue; inaccessible solar elevations are shaded in grey for the equinoxes and solstices for Paranal latitude (). 

In the text 
Fig. 16.
Two examples of calendars showing the visibility of an object from an observatory. Top: galactic Centre from the VLT. Bottom: large Magellanic Cloud from VRO. The fraction of observing time lost due to satellites is indicated by the colour scale (1 indicating that all exposures are damaged), for 300 s exposures with FORS2 and 15 s on the LSST camera, respectively. The satellites are the 60+ thousand from Table 1. The blue shading marks daytime and the blue contours indicate the twilights. The elevation of the object is indicated by the greenish contour lines and the grey shading indicates the times when the object is below 20° elevation. 

In the text 
Fig. A.1.
The spherical triangle to obtain the angle β between the orbit and the parallel of latitude ϕ. The angle κ indicates the anomaly of the satellite measured from the ascending node along its orbit. With λ we indicate the difference between the geocentric longitude of the satellite and that of the ascending node. 

In the text 
Fig. A.2.
Transformation from angular linear density to angular latitudinal density. 

In the text 
Fig. A.3.
Geometry of shell boundaries at the local meridian. C: Earth centre. O: Observatory. S: Intersection of the northern shell boundary with the local meridian. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.