Satellite Optical Brightness

The apparent brightness of satellites is calculated as a function of satellite position as seen by a ground-based observer in darkness. Both direct illumination of the satellite by the Sun as well as indirect illumination due to reflection from the Earth are included. The reflecting properties of the satellite components and of the Earth must first be estimated (the Bidirectional Reflectance Distribution Function, or BRDF). The reflecting properties of the satellite components can be found directly using lab measurements or accurately inferred from multiple observations of a satellite at various solar angles. Integrating over all scattering surfaces leads to the angular pattern of flux from the satellite. Finally, the apparent brightness of the satellite as seen by an observer at a given location is calculated as a function of satellite position. We develop an improved model for reflection of light from Earth's surface using aircraft data. We find that indirectly reflected light from Earth's surface contributes significant increases in apparent satellite brightness. This effect is particularly strong during civil twilight. We validate our approach by comparing our calculations to multiple observations of selected Starlink satellites and show significant improvement on previous satellite brightness models. Similar methodology for predicting satellite brightness has already informed mitigation strategies for next-generation Starlink satellites. Measurements of satellite brightness over a variety of solar angles widens the effectiveness of our approach to virtually all satellites. We demonstrate that an empirical model in which reflecting functions of the chassis and the solar panels are fit to observed satellite data performs very well. This work finds application in satellite design and operations, and in planning observatory data acquisition and analysis.


Introduction
In recent years, numerous large Low Earth Orbit (LEO) satellite constellations have been proposed. There are currently more than 6000 LEO satellites in operation, a sixfold increase over just two years. This is expected to increase exponentially over the next decade. The impact on astronomy research (Tyson et al. 2020;Hu et al. 2022) and on the night sky environment (Venkatesan et al. 2020;Lawrence et al. 2022; Barentine et al. 2023) has been discussed widely. Technical mitigation involves innovation in satellite design, satellite operations, and astronomy data processing and analysis. The science pursued by groundbased wide-field sky surveys such as Rubin Observatory's Legacy Survey of Space and Time (LSST; Ivezić et al. 2019), as well as all other optical observatories, large and small, is impacted by satellite streaks.
After dusk and before dawn, LEO satellites scatter sunlight onto the Earth's surface. This sunlight is both direct and indirect (reflected from Earth). This scattered light can interfere with both casual stargazing and large ground-based observatories. The net effect depends on several variables including: satellite geometry, satellite material properties, satellite orientation, wavelength, satellite location, observatory location, satellite range, and number of satellites. In order to quantify this effect, it is necessary to predict satellite brightness. To make this prediction, we must measure the material properties of satellite surfaces, either directly in the lab or indirectly using satellite brightness observations. We must also know both the orientations and areas of the satellite's surfaces.
This paper presents techniques for calculating the brightness of satellites seen by observers on the Earth's surface. We consider two sources of light that can be scattered by a satellite. First, there is light directly incident from the Sun. Second, we include light scattered from the portion of Earth's surface illuminated by the Sun and visible to the satellite. We refer to the latter as earthshine. We treat the Sun as a planewave source and the Earth as a sphere. Because the geometry and the light sources are relatively simple, we can directly calculate fluxes incident on the satellite. Then, using a simple model for the satellite's reflectance and the position of the satellite in the sky, we calculate the overall satellite brightness. This technique offers a respectable improvement to previous diffuse sphere models of satellites and is computationally efficient. A diagram of the geometry and the light sources is shown in Figure 1.
We model a satellite as a collection of opaque surfaces. The light scattered from each surface is defined by an isotropic Bidirectional Reflectance Distribution Function (BRDF; Greynolds 2015). Physically, this means that a rotation of a surface about its normal vector does not change scattering properties. Even though individual surfaces are isotropic, in most cases the overall effective scatter from the satellite will be anisotropic. The BRDF depends on both the surface's material and surface finish. For example, a smooth metallic surface such as bare aluminum is very specular, but a rough painted surface is mostly diffuse. At this time, we do not consider shadowing between surfaces (mutual shadowing) nor anisotropic BRDFs. If there are complex components on a satellite, a ray-tracing analysis can be used to determine an "effective BRDF" for individual components. Ray tracing allows the inclusion of shadowing and scatter from multiple bounces. We find that this complexity is not required to achieve good observation correlation for Starlink satellite architecture. If lab-measured BRDFs are not available, the BRDF can also be estimated by fitting to satellite brightness observations taken over a variety of solar angles. This corresponds to an indirect measurement of the BRDF.
These calculations have been wrapped in a publicly available Python package called Lumos-Sat (Fankhauser 2023a). Tools for predicting satellite brightness are critical for both constellation operators and observatories. Constellation operators can use Lumos-Sat to include satellite brightness as a design constraint, by quantifying the brightness effects of changing satellite material, geometry, or orientation. Meanwhile, Lumos-Sat lets observatories predict and mitigate the impact of existing satellites on science.

Light Transport
First, we compute a simplified light transport equation. Because the distances between the Sun, the Earth, and the satellite are much larger than the scale of a satellite, we can make a variety of simplifying assumptions. Given a distant light source, a surface, and a distant observer, our goal is to find the flux scattered from the source onto the observer. For a given geometry, the BRDF is a function of the unit vector to the light sourceŵ i and the unit vector to the observerŵ o . The BRDF is defined as follows: The ingoing radiance is L i and the outgoing radiance is L o .
Here, f i is defined as the angle between the surface normaln and the vector to the sourceŵ i . Likewise, f o is defined as the angle between the surface normaln and the vector to the observerŵ o . 3 In this analysis, we consider a very distant point source, which can be treated as a plane wave with flux I in at the surface. The surface has an area A. This geometry is shown in Figure 2.
We can therefore rearrange Equation (1) to find: The fraction of flux scattered by the surface from the source to an observer at distance d is given: Let us recall that f r is the BRDF of the surface. The fraction of light scattered to the observer increases with surface area perpendicular to the source or observer and decreases with distance as an inverse square law. The angular distribution is determined by the BRDF. We will use the light transport equation given in Equation (3) to calculate both the flux of light scattered by the Earth's surface and the flux scattered by the satellite. The next step in our analysis is to find good BRDF models, so that our light transport equation is accurate.

BRDF Models
Our brightness calculations will only be as accurate as our BRDF models. We recommend fitting most measured data to a binomial BRDF model. The model parameters are b ik , c ik , and q: The vectors ρ and ρ 0 are the projection of the outgoing unit vector onto the surface and the projection of the specularly reflected unit vectorŝ onto the surface, respectively. Both of these vectors can be written as functions of the incoming vector w i , the outgoing vectorŵ o , and the surface normal vectorn: r r We note thatŝ is the specularly reflected unit vector. The binomial model is ideal because it enforces physical realism and has been used extensively and proven in 3 When the f ( ) cos o term is removed from the denominator of Equation (1), the resulting quantity is referred to as the Cosine Corrected BRDF (CCBRDF) or Angular Resolved Scatter (ARS). commercial optical analysis (Greynolds 2015). Constructing the BRDF from variables δ and ν is convenient because it enforces positivity, continuity, reciprocity, and isotropic scatter. However, is important to note that the binomial model is empirical. The parameters b ik , c ik , and g must be fit to measured BRDF data. Fitting binomial models does require a degree of caution. The number of coefficients should be kept as small as possible to avoid overfitting. Binomial fits should always be reviewed carefully before using them in a calculation. In general, however, any BRDF model could be used, e.g., the Harvey-Shack model (Nattinger 2020), Ross-Li model (Wanner et al. 1995), or Phong BRDF (Phong 1973. BRDFs can also be interpolated from measured data; however, data are usually taken at a small number of incident angles, so interpolation can lead to extrapolation errors or enhanced measurement noise. Additionally, directly interpolating in spherical coordinates will not work, because the specular peak of the BRDF shifts. For successful interpolation, it is necessary to first make a coordinate transform. We find the best results by interpolating the BRDF as a function of ( is the outgoing direction of scattered light and (θ s , f s ) is specularly reflected direction. BRDFs can also be interpolated as a function of δ, as given in Equation (5). When done cautiously, binomial model development will overcome the shortfalls and complexity of interpolation.
The BRDF of each surface on a satellite can be experimentally measured (Germer & Asmail 1997) or estimated from a catalog of known material BRDFs (Matusik et al. 2003). An effective BRDF for Earth's surface can be fit to data gathered from remote sensors. In particular, we use measurements from NASA's Cloud Absorption Radiometer (Gatebe & King 2016). The aircraft's radiometer has one degree angular resolution in 14 spectral bands 340-2300 nm. We use the 479 nm data averaged over hundreds of images, uncorrected for atmospheric absorption. This data is then fit to a Phong BRDF, of the following form: The parameter K d controls the magnitude of the diffuse component of the BRDF, while K s controls the magnitude of a specular lobe and n controls the width of the specular peak. The vectorsŵ r andŵ o are the specularly reflected unit vector and the outgoing unit vector, respectively. We find that the Phong model yields more reliable results when fit to aggregate data than does a binomial model.

Calculations in the Satellite-centered Frame
Our goal now is to use the equation for the fraction of scattered flux G(source → observer) given in Equation (3) to calculate the flux scattered by a satellite onto an observer on Earth's surface. We start by introducing a coordinate frame that simplifies this calculation. Next, we find the contribution of light directly scattered by the satellite from the Sun. Finally, we include the contribution of light scattered from earthshine.

The Satellite-centered Frame
To simplify calculations, we would like to use a reference frame that reduces dependent variables as much as possible. Our choice of frame is called the satellite-centered frame. 4 The z-axis points along geodetic zenith. The y-axis is in the plane defined by the center of the Earth, the Sun, and the satellite, and is perpendicular to the z-axis. The angle between the y-axis and the vector from the center of the earth to the Sun must be less than 180°. The x-axis is defined by the right-hand rule. This frame is shown in Figure 3.
In the satellite-centered frame, the flux seen by an observer depends only on the angle of the satellite past terminator α, the vector from the satellite to the observer v, the radius of Earth R E , and the geodetic height of the satellite h.

Flux from the Sun
The fraction of flux scattered by a single surface has been derived in Equation (3). We simply sum over all N s surfaces in a satellite to arrive at: The area, BRDF, and normal of satellite surface s are A s , f s , andn s , respectively. The distance from the satellite to the observer is ∥x obs − x sat ∥. The unit vector from the satellite to the Sun is  v sat sun and the unit vector from the satellite to the observer is  v sat obs .

Flux from Earthshine
The flux seen by the observer caused by earthshine can be calculated similarly. Light is scattered first by the Earth's surface, then by the satellite. The flux scattered by a single satellite surface is calculated by integrating over the portion of Earth's surface E that is illuminated by the Sun and is visible to the satellite. This is then summed over the number of satellite surfaces N s . We find the flux seen by an observer due to earthshine: s N E s observer sun 0 s ∂A is a differential area of Earth's surface. G s is the fraction of earthshine flux scattered from ∂A to the observer by a surface. ∂G is the differential fraction of Sun flux scattered from the Sun to the satellite by ∂A. We can approximate Equation (9) by discretizing the integral. This yields: We refer to each area element on the Earth's surface ΔA p as a patch. As the number of patches on the Earth's visible and illuminated surface, N p , goes to infinity, Equations (9) and (10) become equivalent. The fraction of earthshine scattered to the observer is then given by Let us recall again that the area, BRDF, and normal of satellite surface s are A s , f s , andn s , respectively. The distance from the satellite to the observer is ∥x obs − x sat ∥. The unit vector from the satellite to the patch on the Earth's surface is  v sat patch , and the unit vector from the satellite to the observer is  v sat obs . The fraction of Sun flux scattered from the Sun to the satellite by a differential area of of Earth's surface is The area, BRDF, and normal of a patch p on the Earth's surface are ΔA p , f p , andn p , respectively. The distance from the satellite to the observer is ∥x obs − x sat ∥. The unit vector from the satellite to the patch on the Earth's surface is  v sat patch , and the unit vector from the satellite to the observer is  v sat obs .

Discretization of Earth's Surface
It is now necessary to discretize the portion of Earth that is visible to the satellite and illuminated by the Sun. At first glance, using standard spherical coordinates seems like the easiest solution. Unfortunately, this results in a discretization that is heavily weighted at the poles of the Earth. We instead propose the following coordinate system: y The variables ψ and Ω represent the angle off-plane and angle on-plane, respectively, where the plane is defined byŷ andẑ . The radius of the Earth is given as R E . Using this coordinate system results in patches that have much more even spacing. Figure 4 shows 400 patches covering a quarter-sphere. We can see that using standard spherical coordinates results in "bunching up" at the poles.
In order to calculate the flux, we need to know the area ΔA p of each patch p in the mesh. This is approximated using the Jacobian determinant as follows: p E Finally, we must only include patches that are both visible to the satellite and illuminated by the Sun. Consider a patch at the point (x, y, z), measured in the satellite-centered frame. The  . We compare discretizations schemes by uniformly sampling 400 patches in both standard spherical coordinates and our custom coordinates. We see that our custom coordinate system generates a discretization that is more evenly distributed over the Earth's surface.
satellite is visible to a patch if the following is true: A patch is illuminated by the Sun if y > 0. Using Equation (13), these two constraints can be related back to the angle off-plane and the angle on-plane, (Ψ, Ω).
In practice, we find that discretization causes some noise in our brightness calculations. The amplitude of the noise decreases with number of patches and the frequency increases. We recommend applying some smoothing to results for calculations which include earthshine.

Converting Flux to AB Magnitude
It is useful to convert from incident flux to AB magnitude. This allows for comparison between satellite brightness and the brightness of celestial objects. This conversion is simply: The flux incident on the observer from the satellite is I in units of W/m 2 , and f is the frequency of the light in Hz. AB magnitude is defined for a flat spectrum, such as a very hot star. The apparent AB magnitude in a given spectral filter can be found by integrating the flux over the filter's bandpass.

Other Light Sources
Although it is beyond the scope of this paper to include in our simulations, we would also like to offer a "back-of-theenvelope" calculation for the brightness contributions of other light sources. These potentially include celestial objects such as stars, planets, and the Moon. We can modify Equation (3) to estimate the flux scattered by a surface from a source to an observer: Let us recall that f i is the angle between the surface normal and the source and f o the angle between the surface normal and the observer, A is the area of the surface, and d is the range of the satellite. To find an order-of-magnitude calculation for the worst-case scenario, , A ≈ 1m 2 , and d ≈ 250km (very low Earth orbit). We assume light from our celestial source is very specularly reflected, so that BRDF ≈ 10 4 . We can then use Equation (16)  satellites are not point sources and they move quickly across the focal plane of a telescope, it is unlikely that light scattered from stars or city lights will be detectable by observatories like LSST. For a typical LEO satellite at 500 km, the motion across the focal plane is about 0°. 5 per second. This means the effective exposure time on a 0 7 PSF footprint is just 0.5 milliseconds, independent of the camera exposure time. For LSST, a 0 7 PSF has 50 pixels. This means the peak LEO satellite trail surface brightness of a 15 AB magnitude satellite is only 1.8 electrons per pixel. This is buried in the night sky noise from a 15 second LSST exposure.
It is outside the scope of this paper, but future work in satellite brightness modeling should seek to incorporate incident light from the Moon and possibly Venus.

Calculations in the Observer Frame
The satellite-centered frame is ideal for calculations, but we need to know what brightness a ground-based observer will see. We use the horizontal coordinate system shown in Figure 5. We note thatẑ corresponds to an altitude angle of 90°. Given the position of an observer on Earth and the position of a satellite in the sky, our goal is to transform variables from the satellite-centered frame S to the observer frame O. 5 We are given a unit vector from the observer to the satellite [ˆ] v O sat and the satellite's geodetic height h. The unit vector from the observer to the satellite as well as the satellite's height can either be measured by an observer or calculated using the satellite's orbital elements. We also know the unit vector toward the Sun [ˆ] v O sun , which can be calculated using the latitude and longitude of the observer and the time of observation. The basis vectors in the observer's frame are  Figure 5. Geometry as seen in the observer's frame O.
5 A similar process to the one shown below can be used to convert to or from the satellite-centered frame. In particular, satellite operators can use these transformations to convert surface normal vectors from a satellite body frame to our satellite-centered frame. This allows brightness calculations to be generated in simulations or in real time using data from onboard satellite sensors.
Moving forward, we will drop the basis notation for vectors in the observer's frame, O. Quantities that are measured in the satellite-centered frame will be marked with an S. We first need to find the basis vectors of the satellitecentered frame, as measured in the observer's frame.
From geometric inspection, we find the following: Let us remember that all vectors are in the observer's frame. The range of the satellite is d. The angle f is the angular separation between the z-axis of the satellite-centered frame and the z-axis of the observer's frame. Angle θ is the rotation of the satellite-centered frame's z-axis as measured from the x-axis of the observer frame. Using these values, we can find an expression [ˆ] z S O , which is theẑ of the satellite-centered basis measured in the observer frame: x y z sin cos sin sin cos . 23

S O
It is defined thatŷ S is in the plane containing the vector toward the Sun andẑ S . Additionally,ŷ S is orthonormal toẑ S . This gives three constraints which fully defineŷ S : We can solve the three constraints given in Equation (24) to find a and b. Physically, a and b are the components ofŷ S in thê z S direction and thev sun direction, respectively: Finally,x S is defined by the right-hand rule: The transform from the observer reference frame to the satellite-centered frame is We can then find the quantities we need to do calculations in the satellite-centered frame. The vector from the satellite to the observer in the satellite-centered frame is Second, the angle of the satellite past the terminator is given: Let us recall that R E is the radius of the earth, h is the geodetic height of the satellite, and d is the satellite's range. Here, [ˆ] v O sat is the vector from the observer to the satellite as measured in the observer's frame, and [ˆ] v O sun is the vector from the observer to the Sun as measured in the observer's frame.
If we are interested in a particular satellite, we must know that satellite's position in the sky relative to an observer and geodetic height. These quantities can be found using a satellite's orbital elements. SpaceX and other constellation operators publish traditional TLEs (Two-Line Element Sets) on Space-Track.org. For Starlink, more accurate supplemental TLEs are published on celestrak.org. These supplemental TLEs are fit to Starlink propagated ephemerides and covariances, which are available from "Public Files" on Space-Track.org.
Once a satellite's orbital parameters are known from a TLE, its position at a past or future time can be calculated using a Simplified General Perturbations algorithm. In particular, we use the SGP4 Python package. A satellite's altitude and azimuth as seen from a given location on Earth is then found using tools provided by the Astropy software (Astropy Collaboration et al. 2022).

Model Validation
In order to validate our calculations, we create a simple brightness model for existing satellite-selected configurations of Starlink V1.5 and compare our AB magnitude calculations to observations. We show that the BRDFs for each satellite surface can either be measured in a laboratory or found by fitting to brightness measurements. We also compare our calculations and observations to the previously standard satellite brightness model-a diffuse sphere.

Starlink v1.5 Brightness Model
The two largest surfaces on Starlink v1.5 are the solar array and the chassis. The solar array has an area of 22.68 m 2 , and the chassis nadir has an area of 3.64 m 2 . SpaceX has provided BRDF data for each surface and information about the normal vectors of these surfaces in the brightness regime of interest. Using this information, we can create a brightness model for a subset of the Starlink v1.5 satellites. The Starlink v1.5 satellites have gone through a variety of design changes to reduce brightness, so we only consider the Starlink v1.5 satellites with the latest brightness configuration. These satellites have a reflective sticker on the chassis nadir and dark pigmented backsheet on the solar arrays (SpaceX 2022). In nominal operations, the chassis points directly nadir and the solar array is perpendicular to the chassis and in the direction of the Sun.
SpaceX has also provided experimentally measured BRDF data for each of these surfaces. This SpaceX contracted BRDF data was measured by Scatterworks using an SS4 scatterometer operating at a wavelength of 532 nanometers. The data is taken at an angular resolution of 1°. This BRDF data is fit to binomial models, using the methods described in Section 3. The data and the fits are shown in Figure 6.
To find a representative BRDF for Earth's surface, we fit Phong models to CARs data, a technique also described in Section 3. We use data from two missions. First, the CLASIC mission, which gathered BRDF data for vegetation over Oklahoma. Using this data to construct an "effective BRDF" for generic vegetation, we find parameters K d = 0.53, K s = 0.28, and n = 7.31. Second, the CLAMS mission, which gathered BRDF data for the ocean off the East Coast of the United States. For ocean water, we find K d = 0.48, K s = 0.08, and n = 16.45.
We can then feed our knowledge of the satellite's primary surfaces and BRDF models for the Earth's surface into our software and calculate satellite brightness.

Starlink v1.5 Model without BRDF Data
For many satellites, BRDF measurements for primary surfaces may not be available. Not only do BRDF measurements require specialized equipment, but also collaboration between some satellite operators and astronomers may prove difficult. In this case, we need to use a different approach to find an accurate satellite brightness model. We assume minimal information of the satellite-just the pointing directions of primary surfaces. For Starlink v1.5, these surfaces are the chassis, which points directly nadir, and the solar array, which is perpendicular to the chassis and toward the Sun. We then assign each of these two surfaces a Phong BRDF, which has three free parameters. The areas of each surface are set to 1 m 2 -this means that the best-fit albedo of each surface may be greater than unity. In total, our empirical satellite model has six unknown parameters. These parameters can be found by fitting to observed Starlink brightness data. As long as there are a variety of observations over many different solar angles, this makes for an accurate brightness model. This technique uses brightness observations to indirectly measure the overall effective BRDF of the satellite. The satellite's BRDF is then decomposed into the individual surface BRDFs by fitting to brightness observations. The ability to find effective BRDFs for satellite surfaces from only ground-based brightness data is extremely important to astronomers. It makes our brightness modeling approach viable for most satellites.

Diffuse Sphere Model
We can also compare our brightness model to the commonly used diffuse sphere model. For this model, the flux of light scattered by the satellite and incident on an observer is simply: Here, f is the solar phase angle-the angle between the observer, the satellite, and the Sun. The effective area and effective albedo of the satellite are A and ρ, respectively. The satellite's range is d. The flux of the Sun incident on the satellite is I Sun . There is one free parameter, the albedo-area product ρA, which must be best fit to satellite brightness observations. While this model is extremely simple, it has little basis in reality and does not correlate well with brightness measurements. As shown in Figure 7, Starlink v1.5 brightness has a complicated dependence on solar phase angle. This dependence is not fully captured by the diffuse sphere model.  At high phase angles, light is more specularly reflected by the chassis. We note that second-generation Starlinks will off-point solar arrays and use a chassis material with lower diffuse reflection. These two changes reduce brightness at low and intermediate phase angles. Figure 8. The diffuse sphere model shows little variation in brightness. Satellites in the western sky are shadowed by Earth (1). The controlling factors of brightness are satellite range and the illuminated fraction of the sphere. Satellite brightness falls off with satellite range (2). Modeling satellite brightness as a function of just these two variables is overly simplistic. Figure 9. In our model, two distinct brightness peaks are clearly visible. The first is visible earlier in the morning and is caused by light forward-scattered from the chassis nadir. This specular peak is visible just above the eastern horizon (3). The second brightness peak is caused by light backscattering from the solar array. It is most prominent in the western sky near dawn (5). There is also a transition where both peaks are visible (4). Figure 11. With BRDFs best fit to observed data, our model captures both forward scatter from the chassis (7) and backscatter from the solar array (8) to some degree. This model could be improved by using more Starlink v1.5 observations over a wider variety of solar angles and less noisy observations.  (6). This effect becomes more and more pronounced as the Sun rises.

Brightness Distributions on the Night Sky
We run our calculations for Sun altitudes ranging from −27°t o −3°. The Sun's azimuth is fixed at 90°. This corresponds to a period of time prior to sunrise for an observer on the Equator. Satellite altitude is fixed at 550 km 6 . We use our new brightness models both with and without the earthshine contributions. For calculations with earthshine, we use an earthshine discretization of 151 × 151 patches. The diffuse sphere calculation is also shown, to provide a comparison. Results are shown in Figures 8-11. It should be noted that these plots show the brightness that a Starlink v1.5 satellite would appear if it was at a given point on the sky. The center of each image corresponds to a satellite directly zenith. The outer edge is the horizon. Cardinal directions are marked: north is toward the top of the image and east is toward the left side. The grid lines show altitude increments of 10°and azimuth increments of 90°.
Our brightness modeling shows significantly different patterns over the night sky compared to the diffuse sphere model. The diffuse sphere model shows little change in brightness with respect to Sun altitude or satellite position. Additionally, the diffuse model predicts that the brightest satellites will be in the western sky, throughout dawn. On the other hand, our models shows how the location of peak satellite brightness changes from the eastern to the western sky as the Sun rises. Earlier in the morning, bright satellites in the eastern sky are caused by light forward scattering off of the specular chassis. Near dawn, bright satellites in the western sky are caused by light backscattered by the solar array. We see that earthshine creates additional brightness in the eastern sky, particularly during civil and nautical twilight. Earthshine drives satellite brightness up toAB magnitude in a region of the sky where satellites are predicted to be invisible when only illuminated directly by the Sun.

Comparison to Observed Brightness
Using SpaceX contracted data gathered at Mount Lemmon in Arizona by the Pomenis Observatory (Pearce et al. 2018), we can compare our brightness calculations to actual observations of satellite brightness. A scatter plot of these observations is shown in Figure 7. For this correlation, we once again only use observations of Starlink satellites with a pigmented solar array backsheet and reflective chassis sticker. We also exclude observations of satellites that do not have their solar array and chassis in the orientation we used for brightness modeling. The majority of these excluded observations are from satellites, which were raising or lowering their orbit at the time of observation. Starlink satellites performing orbital maneuvers use a distinct brightness mitigation orientation as much as possible. Despite removing these observations, the data set is still somewhat noisy. The Pomenis observatory has a typical measurement error of 0.1 AB magnitude. Additionally, satellite brightness can be driven by a number of small specular components that are very difficult to model. The albedo-area product in the diffuse sphere model is found using a best fit to minimize rms error between observations and the model. For these observations, we find ρA = 0.65. The parameters for our empirical model (without BRDF lab measurements) are similarly found using a best fit to observations. For the chassis, we find Phong parameters K d = 0.34, K s = 0.40, and n = 8.9. For the solar array, we find K d = 0.15, K s = 0.25, and n = 0.26. On the other hand, our model with BRDF lab measurements uses no prior knowledge of brightness observations. Rather, it is modeled using only the geometry and BRDF data provided by SpaceX. Figure 12 shows a comparison between the correlation of our models 7 with observations and the correlation of the diffuse sphere model.
We see that our models are roughly a 50% improvement over the diffuse sphere model, and both can be used to accurately predict satellite brightness. We notice that our model using laboratory measured BRDFs generally underpredicts brightness. This is likely due to brightness driven by satellite components we did not model in our analysis. Astronomers could account for the underestimation by fitting one free parameter (which accounts for light scatter from unmodeled surfaces) to observations. The underprediction is not an issue for satellite operators, as they are most interested in how different satellite designs increase or decrease brightness. It is important to note that this model does not incorporate any prior knowledge of observed brightness, while both our empirical model and the diffuse sphere model require many brightness observations of satellites over a range of solar angles. Additionally, we note that the diffuse sphere model only coincidentally fits the Starlink v1.5 data well, because the solar array backscatters light in a similar manner to the diffuse sphere model. Starlink V2 off-points solar arrays, so the dominant brightness is from forward scattering. A diffuse sphere model will not represent these satellites well. Our model with measured BRDFs is good enough to provide satellite operators with directional knowledge about how changing satellite design can reduce brightness. Additionally, both of our new models can be used by astronomers to estimate which areas of the night sky are least impacted by existing and future satellite constellations.

Discussion
Our motivation for this paper is to provide both satellite constellation operators and astronomers with tools for modeling satellite optical brightness. In particular, we developed a software package known as Lumos-Sat for satellite optical brightness predictions. In order to validate our predictions, we chose a selection of satellites for which there are time-resolved ground-based brightness observations and for which we also have BRDF data for both the satellite components and the Earth. Using SpaceX's existing Starlink v1.5 satellites, we show that our models have better predictive power than the traditional diffuse sphere calculation. For the Starlink v1.5 satellites, we find that our models are about 50% better than a calibrated diffuse sphere model. Our modeling technique can use either laboratory measured BRDFs or BRDFs that are bestfit to satellite observation data. We also note that the diffuse gray sphere model does not capture the specular nature of scattering from the satellite, in particular from the chassis deck. Despite this, the diffuse sphere model has been used in a variety of papers (Hainaut & Williams 2020;Lawler et al. 2021). This is largely because the diffuse sphere model is very 6 Typical Starlink v1.5 orbital altitude per FCC filings. 7 Correlation does not change significantly when including earthshine, because the Pomenis data set does not include many satellite observations where earthshine dominates brightness. For simplicity, our models are shown without earthshine. simple to implement and no better brightness calculation was previously known.
Improving correlation with observation is the next major challenge for this work. Moving forward, our brightness model could be improved by including mutual shadowing between satellite surfaces. For example, our current modeling does not include the effect of the chassis nadir blocking the solar array from an observer's view. This is a logical next step for our work (Cole 2020). Another shortcoming of this technique is that many satellite parameters will not be known by the public. Our model requires the normal vectors, areas, and BRDFs of the primary surfaces of a satellite. While it is known that, during nominal operation, the Starlink v1.5 chassis deck points directly nadir and that the solar array is perpendicular to the chassis deck in the brightness regime of concern, this information may not be known during orbit raising or lowering, or other off-nominal operations. However, once a BRDF model for the primary surfaces of a satellite has been found, it may be possible to use satellite brightness observations to estimate the orientations of the satellite's surfaces. This is a challenge that we leave to future research. Additionally, there is a limit to how much BRDF and geometric data will be known about satellites in the future. Generating these brightness models will require some degree of collaboration between satellite operators and astronomers. As a workaround if cooperation is not possible, BRDFs for a satellite's primary surfaces could be calibrated based on onorbit brightness observations and limited knowledge of satellite geometry made available by FCC filings. Many LEO satellites, by necessity, will have a solar array that points toward the Sun to generate power and a chassis that points toward nadir to provide internet or other communications.
One of the most difficult parts of this analysis is accurate modeling of earthshine. It is beyond the scope of this paper, but earthshine deserves further investigation. We suggest two possible avenues for improvement. First, BRDF data for the Earth's surface could be gathered from the MODIS instrument, which operates on NASA's Terra and Aqua satellites. The MODIS instrument gathers Ross-Li BRDF coefficients for a variety of terrain (Strahler et al. 1999). Another option would be to implement analytic BRDFs, such as those used in the SHARM radiative transfer software (Lyapustin 2005).
Even when there is good collaboration with satellite operators, accounting for every small component of a satellite will be very difficult. Specularly reflecting objects on the scale of centimeters can cause "glints" (changing satellite brightness on short timescales). This includes MLI thermal blankets or even small pieces of aluminum. For these specular materials, scattered flux can change by orders of magnitude over a change in surface orientation of a few degrees or less. This issue is exacerbated by additional moving parts on a satellite busparabolic dishes, laser connections, or other components may rotate quickly and cause glints that are extremely difficult to predict. These glints, if not mitigated by design, are a potentially serious contributor to bogus alerts. 8 Bright glints Figure 12. Comparison of satellite brightness models. We use Pearson correlation to measure model trend. Histograms of the residuals are also shown. For our calculations using BRDFs measured in the lab, we find a Pearson correlation of R = 0.69. This model systematically underpredicts brightness, because it only accounts for the two primary satellite surfaces and ignores smaller components on the satellite. For our empirical model, we find R = 0.68. Because this model has free parameters, it does not underpredict brightness. Both of these models are a respectable improvement over the optimal diffuse sphere model, which has a correlation of R = 0.47. We note that this model has systematic error, as the residual distribution is not Gaussian-like. This demonstrates that, if either lab-measured BRDFs or sufficient satellite brightness observations are available, our techniques are a better choice than the diffuse sphere model. such as iridium flares are so obvious that they would be rejected in the science analysis. Faint glints or flares are more problematic; they can imitate astrophysical flares or interfere with asteroid detection. The only mitigation is better satellite design, which can be driven by our software. Designs can be improved to reduce glints by eliminating or covering offending specular surfaces and using diffuse coatings or materials on moving components and complex geometries. While satellite glints have been detected and are a known issue (Krantz et al. 2021), little has been done to quantify or predict the impact of glints on time-domain astronomy.
The ability to predict satellite brightness has major impacts for constellation operators. It allows satellite operators to use brightness as a constraint during the design process. If satellite brightness is readily predictable, operators can choose to use satellite materials and configurations or satellite conops that reduce brightness. Observations presented in this paper correspond to the first-generation Starlink satellites. While these satellites did not include brightness as a design constraint early in development, SpaceX has used brightness analysis similar to that presented in this paper to inform brightness mitigating designs on the recently developed second-generation Starlinks. Using this predictive analysis ensured that the most effective and efficient brightness mitigations were implemented on second-generation Starlinks. These mitigations include the use of specular materials on the chassis that have order-ofmagnitude reductions in diffuse scatter, 9 dark paint where specular materials are not effective, and off-pointing of the solar arrays (SpaceX 2022). These improvements can be quantified using Lumos-Sat before satellites are in space.
Our calculations show that the sky near zenith can have bright satellite trails in astronomical twilight during dawn and dusk night observatory operations (Ivezić et al. 2019). In our simulations, earthshine causes notable satellite brightness on the eastern horizon. This possibly impacts the search for potentially hazardous asteroids (PHAs), which scans the sky toward the Sun at dawn and dusk to find asteroids interior to the earth's orbit. This works by connecting candidate detections in pairs of exposures during a night (tracklets) and extending the tracklets to additional nights (Schwamb et al. 2023). PHAs can be detected higher than 30°above the horizon. Unfortunately, earthshine can produce visible satellite trails in this region of the sky up to 40°. Bogus detections caused by satellites create noise in the tracklet extrapolation process, which causes PHA detections to fail. Our model makes it possible to quantify the impact of satellites on the PHA detection process and telescope operations more generally. Improved brightness analysis in this regime can be incorporated into PHA detection processes and satellite design and operation.
For the astronomy community, accurate brightness modeling is important to predict when and where satellites are most likely to interfere with observatory operations and data quality. Satellites have complex light scattering properties, so their optical brightness is highly dependent on their location in the sky. This knowledge can be used to inform better telescope scheduling algorithms to "dodge" LEO satellites. Satellite brightness predictions could even be used to create better satellite streak removal software. Using orbital parameters, as well as brightness prediction, the width and magnitude of a satellite streak in an image can be estimated ahead of time.
Additionally, astronomers can better quantify the scientific impacts of LEO satellites.
As this is one of the first attempts to accurately predict satellite brightness, not all potential applications are known. Brightness models still need to be developed for other satellites, such as Starlink's V2 and V2 Mini, Amazon's Project Kuiper, AST's SpaceMobile, and OneWeb. We hope this promising technique and the open-source software will be further developed and utilized in the future to improve satellite design, as well as astronomy operations and data analysis.