© 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 two-satellite 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 low-Earth orbit (LEO, altitude less than 2000 km), aimed to provide cell phone service with worldwide coverage. The so-called 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 mega-constellation 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 two-way 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 low-Earth orbit.

The contrast with traditional constellations is overwhelming. The examples mentioned above, or the plans to launch new constellations of nano- and pico-satellites, imply challenges of their own and contribute to the issues of orbital crowding, space debris management, radio-electric noise, etc., but they imply numbers of satellites orders of magnitude smaller and, most often, with elements much fainter, than the mega-constellations discussed in here.

Several new generation mega-constellations 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¹. Furthermore, other satellite operators have indicated their intentions to launch mega-constellations, 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.

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 mega-constellations on ground-based optical and near-infrared 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 multi-fibre 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 near-infrared 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 near-to-re-entry 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 end-of-life 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 non-operational 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 mid-latitude 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 over-densities 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°). Over-densities 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 over-densities are also seen in the all-sky 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 all-sky maps show increasing satellite densities towards the horizon, though additional over-densities 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, over-densities 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, all-sky simulations may be iterated, including some random initialisation of parameters at each run to average-out 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 all-sky simulations, an observation-oriented approach is also needed. This requires specifying all the parameters required for the all-sky simulations and, also, selecting the observation direction (azimuth, elevation), field-of-view (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 fields-of-view. 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 field-of-view towards zenith were counted and averaged over 1000 simulations.

Fig. 3. Number of satellites visible in an exposure depending on the exposure time texp and the instrument field-of-view, specified by its circular radius Rfov = 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 texp and field-of-view (). The solid lines are predictions based on Eq. (1).

The simulations show a dependence with both field-of-view 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 field-of-view during the exposure, and has the form

(1)

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 field-of-view area of A_fov and width L_fov, and exposure time t_eff. For comparison with the simulations we use a field-of-view 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 fields-of-view, 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 fields-of-view 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². For a single satellite at orbital altitude h_sat, the probability density per unit surface area as a function of ϕ is:

(2)

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 field-of-view of an instrument from an observatory located on Earth to obtain the fraction of the sample present within the instrument field-of-view. For the remainder of the paper, we make the simplifying assumption that P(ϕ, i, h_sat) is constant over the instrument field-of-view. This assumption is generally true for the small fields-of-view under consideration in the remaining analysis. It has the advantage that it removes any dependence on the precise shape and orientation of the instrument field-of-view, and instead solely depends on the sky area covered by the field-of-view. The assumption allows us to evaluate Eq. (2) for a line-of-sight (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

(3)

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 line-of-sight 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 (field-of-view 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 line-of-sight 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 Lfov = 6′ and the exposure time is 300 s. The satellites are those from Table 1. In the black region at the south-east, all satellites are already in the shadow of the Earth. The edges running from north-east to south-west correspond to constellation shells at lower altitudes, whose south-east 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 mega-constellations 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 near-infrared (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 state-of-the-art). 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:

(4)

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 cross-section, 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 plane-parallel 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

(5)

where both d_sat and R_sat are expressed in the same units.

The photometric cross-section is the only parameter that depends on the satellite’s physical properties in this model. Hainaut & Williams (2020) adopted m² for the first-generation 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

(6)

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 cross-section 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.

2.2.2. Effective magnitude and limiting magnitude

During an exposure of duration t_exp, a satellite will leave a trail of length ω_satt_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 point-like 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:

(7)

Figure 6 shows the effective magnitude for an example. While not directly relevant for low-altitude 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 m1000 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

(8)

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 8-m 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 signal-to-noise ratio (S/N) by up to θ_sat/θ_atm and its peak intensity by up to (θ_sat/θ_atm)², 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 wide-field 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.

All-sky cameras will record only the brightest satellites and flares. Only the deepest wide-field 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 wide-field astrophotography is considered as follows: a wide-angle lens (75° ×55°) is used with a 15 mega-pixel 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 wide-angle 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 fibre-fed spectrographs, the resulting data contain no spatial information at all; for long-slit spectrographs, the spatial information is available in only one direction. Except in the case of integral-field spectrographs, the data will therefore not include a tell-tale 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 low-resolution spectrographs on small to medium telescopes and high-resolution 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 solar-type 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 science-case 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 high-altitude satellites (OneWeb, GuoWang) would remain illuminated during a large fraction or even the totality of the night; also, high-altitude 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 high-altitude 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 over-density 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 over-densities 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):

(9)

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⁻². 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⁻² range (0.3–0.7 μcd m⁻²), with peaks around 26.5 mag arcsec⁻² (2.7 μcd m⁻²) at the cusps of the constellations. The surface brightness of the dark night sky is around V = 21.7 (Patat 2008) (225 μcd m⁻²), 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 v1000 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⁻⁷ m to 5 m at altitudes above 200 km. They estimate that this excess is 16.2 μcd m⁻² (24.6 mag arcsec⁻²) and can reach 21.1 μcd m⁻² (24.3 mag arcsec⁻²) 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³ 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.

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₁ × L₂ for the field of view (with L₁ the length and L₂ the width of the FOV, L₁ > L₂) and L_fov = L₁ in the second term of Eq. (1) – this maximises the cross-section 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 non-saturated 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) multi-fibre 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 Monte-Carlo 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 4MOST-HiRes; 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. 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 field-of-view 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 sky-noise 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 trans-Neptunian 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 gamma-ray 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 wide-field cameras, the limiting magnitude for satellite trails scales inversely with the focal length of the lens, with every other parameter remaining constant. A wide-angle 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 strings-of-pearls 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 high-resolution spectrographs or échelle spectrographs, even on very large telescopes (see the examples of UVES and ESPRESSO on the VLT). However, low- to medium-resolution spectrographs on medium to large telescopes will detect all or many satellites. Furthermore, contrary to imagers, where a satellite leaves a tell-tale trail in the data, slit and fibre spectrographs do not record spatial information. While high-S/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 real-time. 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 N-band (8–13 μm) and up to 50 Jy in the M and Q-bands (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⁻⁶ 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⁻⁴ to 10⁻⁶ 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 trans-Neptunian 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, high-altitude constellations remain illuminated by the Sun much longer than low-altitude ones. The apparent magnitude of a high-altitude 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 cross-section 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⁵ 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⁶ 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 on-the-fly, 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 re-open it just after it left. As the satellites move at apparent speeds of ∼0.1 to ∼1° s⁻¹, 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 up-to-date 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 field-of-view. 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 up-to-date (as non-keplerian effects, including active orbital corrections, modify the orbit with a time-scale of up to a few days); it must be precise (the current standard Two-Line 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 real-time 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 read-out 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 un-recoverable. 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 well-quantified 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 tell-tale trail reveals the passage of a satellite, which can only be detected as an additional solar-type 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 re-acquire. Others will be lost forever (a short transient phenomenon, like the optical counterpart of a gravitational wave event).

Overall, no mitigation method will single-handedly 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.

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

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 line-of-sight 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:

(A.1)

(A.2)

(A.3)

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:

(A.4)

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 = Λκ:

(A.5)

Next, we obtain the expression for the latitudinal angular probability density due to one single orbit with one single satellite:

(A.6)

But, as shown in Eqn. A.3, cos β = cos i/cos ϕ and, thus:

(A.7)

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₂ = 2p:

(A.8)

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):

(A.9)

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)². Also, by using that cos²ϕ − cos²i = sin²i − sin²ϕ and that both expressions are differences of squares, we end up with four equivalent formulations of our probability surface density function:

(A.10)

(A.11)

(A.12)

(A.13)

(A.14)

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:

(A.15)

(A.16)

(A.17)

This last integral is reduced to a straightforward arcsine with the variable change x = sin ϕ:

(A.18)

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,

(A.19)

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, ϕ₀) and of the satellite (θ, ϕ), Eqn. A.19 becomes

(A.20)

(A.21)

(A.22)

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:

(A.23)

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.24)

A.4. Apparent angular velocity

A satellite observed in a given detection may either be on the north-bound half of its orbit, or on the south-bound 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:

(A.25)

from which we get the longitudes of the ascending nodes for both possible orbits:

(A.26)

(A.27)

The geocentric velocity vector of a satellite is obtained from:

(A.28)

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:

(A.29)

(A.30)

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:

(A.31)

(A.32)

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:

(A.33)

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 line-of-sight and the spherical shell at the satellite

• α_⊙: solar phase angle, angle Sun-Satellite-observer

• χ: 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⁻¹ s⁻¹)

• κ: 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_satb: 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, point-like 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₁: 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°). Over-densities 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, over-densities 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 texp and the instrument field-of-view, specified by its circular radius Rfov = 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 texp and field-of-view (). 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 Lfov = 6′ and the exposure time is 300 s. The satellites are those from Table 1. In the black region at the south-east, all satellites are already in the shadow of the Earth. The edges running from north-east to south-west correspond to constellation shells at lower altitudes, whose south-east 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 m1000 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 wide-angle 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 high-altitude 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 v1000 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 4MOST-HiRes; 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