# Motif-topology and Reward-learning improved Spiking Neural Network for
Efficient Multi-sensory Integration

###### Abstract

Network architectures and learning principles are key in forming complex
functions in artificial neural networks (ANNs) and spiking neural networks
(SNNs). SNNs are considered the new-generation artificial networks by
incorporating more biological features than ANNs, including dynamic spiking
neurons, functionally specified architectures, and efficient learning
paradigms. In this paper, we propose a Motif-topology and Reward-learning
improved SNN (MR-SNN) for efficient multi-sensory integration. MR-SNN contains
13 types of 3-node Motif topologies which are first extracted from independent
single-sensory learning paradigms and then integrated for multi-sensory
classification. The experimental results showed higher accuracy and stronger
robustness of the proposed MR-SNN than other conventional SNNs without using
Motifs. Furthermore, the proposed reward learning paradigm was biologically
plausible and can better explain the cognitive McGurk effect caused by
incongruent visual and auditory sensory signals.

Index Terms—  Spiking Neural Network, Multi-sensory Integration, Motif
Topology, Reward Learning

##  1 Introduction

Spiking neural networks (SNNs) are considered as the third generation of
artificial neural network (ANNs) [1], which are biologically plausible at both
network architectures and learning paradigms. The neurons, synapses, networks,
and learning principles in SNNs are far more complex and powerful than those
used in ANNs [2].

This paper highlights two important features of SNNs, which are also the most
differences between SNNs and ANNs, including specific network architectures
and efficient learning principles. For the architectures, specific cognitive
topologies learned from evolution are highly sparse and efficient in SNNs [3],
instead of pure densely-recurrent ones in counterpart ANNs. For the learning
principles, SNNs are more tuned by biologically-plausible plasticity
principles, e.g., the spike timing-dependent plasticity (STDP) [4], short-term
plasticity (STP) [5] (which further includes facilitation and depression),
lateral inhibition, Long-Term Potentiation (LTP), Long-Term Depression (LTD),
Hebbian learning, synaptic scaling, synaptic redistribution and reward-based
plasticity [6], instead of by the pure multi-step backpropagation (BP) of
errors in ANNs. The SNNs encode spatial information by fire rate and temporal
information by spike timing, giving us hints and inspiration that SNNs are
also powerful in integrating visual and auditory sensory signals.

In this paper, we focus more on the key feature of SNNs at information
representation, integration, and classification. Hence, a Motif-network and
Reward-learning improved SNN (MR-SNN) is proposed and then will be verified
efficient on multi-sensory integration. The MR-SNN contains at least three key
advantages. First, specific Motif circuits can improve accuracy and robustness
at single-sensory and multi-sensory classification tasks. Second, MR-SNN can
reach a relatively little higher and lower computation cost than other state-
of-the-art SNNs without Motifs. Third, the reward learning paradigm can better
describe the McGurk effect [7], which describes an interesting psychological
phenomenon that a new but reasonable audio concept might generate as a
consequence of giving incongruent visual and auditory inputs. It exhibits a
biologically-like behavior by using biologically plausible learning
principles.

Fig. 1: The architecture of MR-SNN. (a) The architecture of MR-SNN on the
multi-sensory integration task. (b) The example of artificial 3-node Motifs.
(c) The spiking neuron with dynamic membrane potential.

##  2 Related works

For biological connections, the visual-haptic integration pathway [8], visual-
vestibular integration pathway [9], and visual-auditory integration pathway
[10] have been identified and played important roles in opening the black box
of cognitive functions in the biological brain [11].

For learning paradigms, besides biologically plausible principles (e.g., STDP,
STP), some efficient algorithms have been proposed, such as ANN-to-SNN
conversion (i.e., directly train ANNs with BP first and then equivalently
convert to SNNs) [12], proxy gradient learning (i.e., replacing the non-
differential membrane potential at firing threshold by an infinite gradient
value) [13], and temporal BP learning (e.g., SpikeProp) [14].

For application tasks, SNNs have shown powerful on motion planning [15],
visual pattern recognition [16, 5] and probabilistic inference [17].

##  3 Methods

###  3.1 The Motif topology

All the 13 types of 3-node Motifs in Fig. 1(b) have been used to analyze
functions in various types of biological networks [11]. Here, for simplicity,
we transform the synaptic weights in SNNs into the range of 0 to 1 by using
the Sigmoid function first and then count the Motif distribution to generate
the Motif mask, named as $M_{t}^{r,l}$ at the recurrent layer $l$.

###  3.2 The spiking neurons and recurrent architectures

The leaky integrated-and-fire (LIF) neurons are the most simple and basic
neurons for simulating biologically plausible spiking dynamics, including non-
differential membrane potential and refractory period, as shown in Fig. 1(c).

The recurrent connections between LIF neurons are usually used to simulate the
network-scale dynamics in SNNs. As shown in Fig. 1(a), we designed a four-
layer SNN architecture containing visual and auditory input encoding, multi-
sensory integration in a recurrent hidden layer, and readout layer. The
synaptic weights between neurons in hidden layers are adaptive, while
predefined Motif masks decide the connections between neurons. The membrane
potentials in hidden layers are the integration of feedforward potential and
recurrent potential, shown as follows:

$\left\\{\begin{array}[]{l}\begin{matrix}S(t)=S^{f}(t)+S^{r}(t)\end{matrix}\\\
\begin{matrix}V_{i}(t)=V_{i}^{f}(t)+V_{i}^{r}(t)\end{matrix}\\\
\begin{matrix}C\frac{dV_{i}^{f}(t)}{dt}=g(V_{i}(t)-V_{rest})(1-S(t))+\sum_{j=1}^{N}W^{f}_{i,j}X_{j}\end{matrix}\\\
\begin{matrix}C\frac{dV_{i}^{r}(t)}{dt}=\sum_{j=1}^{N}W^{r}_{i,j}S(t)\cdot
M_{t}^{r,l}\end{matrix}\\\ \end{array}\right.\textbf{,}$ (1)

where $C$ is the capacitance, $S(t)$ is the firing flag at timing $t$,
$V_{i}(t)$ is the membrane potential of neuron $i$ that incorporates feed-
forward $V_{i}^{f}(t)$ and recurrent $V_{i}^{r}(t)$, $V_{rest}$ is the resting
potential, $W_{i,j}^{f}$ is the feed-forward synaptic weight from the neuron
$i$ to the neuron $j$, $W_{i,j}^{r}$ is the recurrent synaptic weight from the
neuron $i$ to the neuron $j$. $M_{t}^{r,l}$ is the mask that incorporates
Motif topology to influence the feedforward propagation further. The
historical information is stored in the forms of recurrent membrane potential
$V_{i}^{r}(t)$, where spikes are generated after potential reaching a firing
threshold, shown as follows:

$\left\\{\begin{array}[]{l}V_{i}^{f}(t)=V_{reset},S^{f}(t)=1\quad
if(V_{i}^{f}(t)=V_{th})\\\ V_{i}^{r}(t)=V_{reset},S^{r}(t)=1\quad
if\left(V_{i}^{r}(t)=V_{th}\right)\\\ S^{f}(t)=1\quad
if\left(t-t_{s^{f}}<\tau_{ref},t\in\left(1,T_{1}\right)\right)\\\
S^{r}(t)=1\quad
if\left(t-t_{s^{r}}<\tau_{ref},t\in\left(1,T_{2}\right)\right)\end{array}\right.\textbf{,}$
(2)

where $V_{i}^{f}(t)$ is the feed-forward membrane potential, $V_{i}^{r}(t)$ is
the recurrent membrane potential, $S^{f}(t)$ and $S^{r}(t)$ are spike flags of
feed-forward and recurrent membrane potentials, respectively, $V_{reset}$ is
reset membrane potential.

###  3.3 The local principle of gradient approximation

The membrane potential at the firing time is a non-differential spike, so
local gradient approximation (pseudo-BP) [18] is usually used to make the
membrane potential differentiable by replacing the non-differential part with
a predefined number, shown as follows:

$Grad_{local}=\frac{\partial S_{i}(t)}{\partial
V_{i}(t)}=\left\\{\begin{array}[]{cc}1&if\left(\left|V_{i}(t)-V_{th}\right|<V_{win}\right)\\\
0&else\end{array}\right.\textbf{,}$ (3)

where $Grad_{local}$ is the local gradient of membrane potential at the hidden
layer, $S_{i}(t)$ is the spike flag at neuron $i$, $V_{i}(t)$ is the membrane
potential of neuron $i$, $V_{th}$ is the firing threshold. This approximation
makes the membrane potential $V_{i}(t)$ differentiable at the spiking time
between an upper bound of $V_{th}+V_{win}$ and a lower bound of
$V_{th}-V_{win}$.

###  3.4 The global principle of reward learning

The reward propagation has been proposed in our previous work [18], where the
reward signal is directly given to all hidden neurons without layer-to-layer
backpropagation, shown as follows:

$\left\\{\begin{array}[]{l}Grad_{R}=B_{rand}^{f,l}\cdot R_{t}-h^{f,l}\\\
\Delta W_{t}^{f,l}=-\eta^{f}(Grad_{R})\\\ \Delta
W_{t}^{r,l}=-\eta^{r}\left(Grad_{t+1}+Grad_{R}\right)\cdot
M_{t}^{r,l}\end{array}\right.\text{,}$ (4)

where $h^{f,l}$ is the current state of layer $l$, $R_{t}$ is the predefined
reward for current input signal. A predefined random matrix $B_{rand}^{f,l}$
is designed to generate the reward gradient $Grad_{R}$. $W_{t}^{f,l}$ is the
synaptic weight at layer $l$ in feed-forward phase, $\Delta W_{t}^{r,l}$ is
the recurrent-type synaptic modification at layer $l$ which is defined by both
$Grad_{R}$ by reward learning and $Grad_{t+1}$ by iterative membrane-potential
learning [19]. The $M_{t}^{r,l}$ is the mask that incorporates Motif topology
to further influence the propagated gradients.

##  4 Experiments

###  4.1 Visual and auditory Datasets

The MNIST dataset [20] was selected as the visual sensory dataset, containing
70,000 28$\times$28 one-channel gray images of handwritten digits from zero to
nine. Among them, 60,000 images are selected for training, while the remaining
10,000 ones are left for testing. The TIDigits dataset [21] was selected as
the auditory sensory dataset, containing 4,144 spoken digit recordings from
zero to nine, corresponding to those in the MNIST dataset. Each recording was
sampled as 20KHz for around 1 second. Some examples are shown in Fig. 1(a).

###  4.2 Experimental configurations

We built the SNN in Pytorch and trained on TITAN Xp GPU. The network
architectures for MNIST and TIDigits are the same, containing one
convolutional layer (with a kernel size of 5$\times$5), one full-connection or
integrated layer (with 200 LIF neurons), and one output layer (with ten output
neurons). The capacitance $C$ is 1$\mu$ F/cm2, conductivity $g$ is 0.2 nS,
time constant $\tau_{ref}$ is 1 ms, resting potential $V_{rest}$ is equal to
reset potential $V_{reset}$ with 0 mV. The learning rate is $1e$-$4$, the
firing threshold $V_{th}$ is 0.5 mV, the simulation time $T$ is set as 28 ms,
the gradient approximation range $V_{win}$ is 0.5 mV.

As shown in Fig. 1(a), before being given to the input layer, the raw input
signals were encoded to spike trains first by random sampling (for that in
spatial image data) or by the temporal encoding of MFCC [22] (for that in
temporal auditory data). The encoding aimed to convert the input to spike
trains by comparing each number with a random number generated from Bernoulli
sampling at each time slot of time window $T$.

Algorithm 1 The MR-SNN algorithm.

1\. Initialize the network by resetting weights and all related parameters.

2\. Encode to spike trains from the raw numbers in datasets.

3\. Train the proposed MR-SNN on single-sensory datasets.

3.1 Learn synaptic weights $w_{ij}$ and Motif masks $M_{t}^{r,l}(s)$,
$M_{t}^{r,l}(t)$ of networks with pseudo-BP and reward-learning in two
datasets, respectively.

3.2 Save the Motif masks during single-sensory classification, in which the
spatial $M^{r,l}_{t}(s)$ and temporal $M^{r,l}_{t}(t)$ were generated from
visual and auditory tasks, respectively.

4\. Synthesize the integrated masks $M^{r,l}_{t}$ from spatial and temporal
masks, where
$M^{r,l}_{t}=1/2\times\left(M^{r,l}_{t}(s)+M^{r,l}_{t}(t)\right)$. And retrain
the $w_{i,j}$ of the network with frozen Motif mask $M^{r,l}_{t}$.

5.Test the performance of SNNs using these new masks in the multi-sensory
classification tasks, and make comparison.

###  4.3 Learning Motif-topology from single-sensory tasks

The Motif distributions in visual and auditory categorizations were shown in
Fig. 2(a,c), learned from different single-sensory classification tasks. We
further set “credible frequency” by multiplying the occurrence frequency and
$1-P$, where the $P$ is the P-value of a selected Motif after comparing it to
2,000 matrices that each element is uniformly and randomly distributed,
indicating more credibility given by a lower P-value.

The 3rd Motif in MNIST and the 13th Motif in TIDigits are the most credibly
frequent topologies. The corresponding visualization of these two types of
Motif distributions is shown in Fig. 2(b,d), where the distribution on visual
task is sparser than that on auditory task, which is also corresponding with
the biological findings [23, 24]. The following sections will verify these key
Motif connections playing important roles in improving accuracy and robustness
during spatial and temporal classifications.

Fig. 2: Key Motif distributions and examples of visualization learned from the
MNIST (a, b) and TIDigits (c, d) datasets. All figures are averaged over 5
repeating experiments with different random seeds.

###  4.4 Motif-topology for stronger robust network

As shown in Table 1, the higher accuracies than SVM demonstrated the capacity
of SNNs for dealing with sensory information. Furthermore, the experimental
results showed that the SNNs using feedforward and Motif-topology (FF-Motif)
were more robust than those using the feedforward connection (FF). The
standard SNNs-Motif with dopamine learning [18] reached a higher accuracy than
other algorithms, with an improvement of around 0.03%, 0.43%, and 0.03% for
MNIST, TIDigits, and integrated sensory, respectively. After giving uniformly
distributed random noise to raw input, the accuracy improvement of SNN-Motif
was higher than not giving noise, reaching around 7.81%, 7.36%, and 0.76% for
MNIST, TIDigits, and integrated sensory, respectively. It showed that these
special Motif-topologies (the third one in MNIST and the 13th one in TIDigits)
also contributed to the robust computation of SNNs during classification.

###  4.5 Motif-topology for better multi-sensory integration

For the experiments of multi-sensory integration, the visual and auditory
signals were simultaneously given to an SNN using integrated Motif
distributions learned from single-sensory tasks. As shown in Table 1, the
learning accuracies for the SNNs using single visual (98.50%) or auditory
(98.20%) sensories were lower than that using multi sensories (99.55%). In
addition, the means of accuracies for single-sensory tasks were improved up to
0.43% for SNNs using Motif distributions than those without using them.

Table 1: The comparison of accuracy (%) after giving Motif topology with ($*$) or without ($-$) 80% additional noise to the input layers. Topology | MNIST | TIDigits | Integration
---|---|---|---
SVM[25] | 97.92 | 94.53 | –
FF($-$)[18] |  98.50$\pm$0.02 |  98.20$\pm$0.14 |  99.55$\pm$0.06
FF-Motif($-$) |  98.53$\pm$0.09 |  98.63$\pm$0.10 |  99.58$\pm$0.09
FF($*$) |  90.55$\pm$1.17 |  83.71$\pm$2.23 |  98.37$\pm$0.09
FF-Motif($*$) |  98.36$\pm$0.11 |  91.07$\pm$1.42 |  99.13$\pm$0.11

We found a notable accuracy increase for the paradigms with input noises,
reaching 98.36%, 91.07%, and 99.13% for visual, auditory, and integrated
sensory, respectively. In conclusion, the Motif-topology would improve the
accuracy of all three tasks under circumstances of whether giving additional
noise or not.

###  4.6 Motif-topology for the explainable McGurk effect

The McGurk effect [7] describes an interesting psychological phenomenon where
incongruent voice saying and face articulating will result in a new concept.
For example, we believe we hear as [d] but actually voice [b] and face [g] are
given, as shown in Fig. 3(a).

Fig. 3: McGurk effect and probability density distribution of feature
distances at integrated layers. [$tu$:] and [$\theta ri$:] are the
pronunciations of spoken two and three, [2] and [3] are the hand-digit images
of two and three.

We simulated this phenomenon with our proposed MR-SNN model. For simplicity,
we used handwritten digits (images of [2] and [3]) and spoken digits
(pronunciation of [$tu$:] and [$\theta ri$:]) to represent face articulating
and voice saying, respectively.

The two histogram figures in Fig. 3(b) were calculated to approximate the
probability density of feature distance in the integrated layer (see
experimental configurations for more details). We could only see two
distinguished normally distributed distributions representing congruent
sensory signals to SNNs using pseudo-BP, including [$\theta ri$:, 3] (orange
bar) and [$tu$:, 2] (blue bar). However, besides two old concepts, a new
concept (green) would generate for SNNs using reward learning and integrated
Motif circuits, especially after receiving incongruent sensory signals [$tu$:,
3] (green bar). This is a small step in understanding the McGurk effect from
simple algorithmic comparison. We will further discuss it in our further
works.

##  5 Conclusion

We proposed a new Motif-topology and Reward-learning improved SNN (MR-SNN),
exhibiting two important features. First, the Motif topology learned from
spatial or temporal data could improve accuracy and robustness than standard
SNNs without using Motifs. Second, with the biologically plausible reward
learning, the proposed MR-SNN could simulate the McGurk effect found in
cognitive integration of human brains, where the traditional SNNs using
pseudo-BP would be failed. The source code of the models and experiments can
be found at https://github.com/thomasaimondy/Motif-SNN.

## References

  * [1] Wolfgang Maass,  “Networks of spiking neurons: the third generation of neural network models,”  Neural Networks, vol. 10, no. 9, pp. 1659–1671, 1997.
  * [2] Demis Hassabis, Dharshan Kumaran, Christopher Summerfield, and Matthew Botvinick,  “Neuroscience-inspired artificial intelligence,”  Neuron, vol. 95, no. 2, pp. 245–258, 2017.
  * [3] Liqun Luo,  “Architectures of neuronal circuits,”  Science, vol. 373, no. 6559, pp. eabg7285, 2021.
  * [4] Tielin Zhang, Yi Zeng, Dongcheng Zhao, and Mengting Shi,  “A plasticity-centric approach to train the non-differential spiking neural networks,”  in The 32th AAAI Conference on Artificial Intelligence (AAAI-2018), 2018.
  * [5] Tielin Zhang, Yi Zeng, Dongcheng Zhao, and Bo Xu,  “Brain-inspired balanced tuning for spiking neural networks.,”  in IJCAI, 2018, pp. 1653–1659.
  * [6] Wickliffe C Abraham and Mark F Bear,  “Metaplasticity: the plasticity of synaptic plasticity,”  Trends in neurosciences, vol. 19, no. 4, pp. 126–130, 1996.
  * [7] K. Tiippana,  “What is the mcgurk effect?,”  Frontiers in Psychology, vol. 5, 2014.
  * [8] Marc O Ernst and Martin S Banks,  “Humans integrate visual and haptic information in a statistically optimal fashion,”  Nature, vol. 415, no. 6870, pp. 429–433, 2002.
  * [9] Zhixian Cheng and Yong Gu,  “Distributed representation of curvilinear self-motion in the macaque parietal cortex,”  Cell Reports, vol. 15, no. 5, pp. 1013–1023, 2016.
  * [10] Barry E Stein, M Alex Meredith, W Scott Huneycutt, and Lawrence McDade,  “Behavioral indices of multisensory integration: orientation to visual cues is affected by auditory stimuli,”  Journal of Cognitive Neuroscience, vol. 1, no. 1, pp. 12–24, 1989\.
  * [11] K. Shen, G. Bezgin, R. M. Hutchison, J. S. Gati, and A. R. Mcintosh,  “Information processing architecture of functionally defined clusters in the macaque cortex,”  Journal of Neuroscience the Official Journal of the Society for Neuroscience, vol. 32, no. 48, pp. 17465–76, 2012.
  * [12] Peter U Diehl, Daniel Neil, Jonathan Binas, Matthew Cook, Shih-Chii Liu, and Michael Pfeiffer,  “Fast-classifying, high-accuracy spiking deep networks through weight and threshold balancing,”  in The 2015 International Joint Conference on Neural Networks (IJCNN-2015). IEEE, 2015, pp. 1–8.
  * [13] Jun Haeng Lee, Tobi Delbruck, and Michael Pfeiffer,  “Training deep spiking neural networks using backpropagation,”  Frontiers in Neuroscience, vol. 10, 2016.
  * [14] Sander M Bohte, Joost N Kok, and Han La Poutre,  “Error-backpropagation in temporally encoded networks of spiking neurons,”  Neurocomputing, vol. 48, no. 1, pp. 17–37, 2002.
  * [15] Elmar Rueckert, David Kappel, Daniel Tanneberg, Dejan Pecevski, and Jan Peters,  “Recurrent spiking networks solve planning tasks,”  Scientific reports, vol. 6, pp. 21142, 2016.
  * [16] Tielin Zhang, Yi Zeng, Dongcheng Zhao, and Mengting Shi,  “A plasticity-centric approach to train the non-differential spiking neural networks,”  in Thirty-Second AAAI Conference on Artificial Intelligence, 2018\.
  * [17] Alireza Soltani and Xiao-Jing Wang,  “Synaptic computation underlying probabilistic inference,”  Nature neuroscience, vol. 13, no. 1, pp. 112–119, 2010.
  * [18] Tielin Zhang, Shuncheng Jia, Xiang Cheng, and Bo Xu,  “Tuning convolutional spiking neural network with biologically plausible reward propagation,”  IEEE Transactions on Neural Networks and Learning Systems, pp. 1–11, 2021.
  * [19] Paul J Werbos,  “Backpropagation through time: what it does and how to do it,”  Proceedings of the IEEE, vol. 78, no. 10, pp. 1550–1560, 1990.
  * [20] Yann LeCun,  “The mnist database of handwritten digits,”  http://yann. lecun. com/exdb/mnist/, 1998.
  * [21] R. Gary Leonard and George Doddington,  “Tidigits ldc93s10,”  Web Download. Philadelphia: Linguistic Data Consortium, 1993.
  * [22] Beth Logan,  “Mel frequency cepstral coefficients for music modeling,”  in In International Symposium on Music Information Retrieval. Citeseer, 2000.
  * [23] William E Vinje and Jack L Gallant,  “Sparse coding and decorrelation in primary visual cortex during natural vision,”  Science, vol. 287, no. 5456, pp. 1273–1276, 2000.
  * [24] Tomáš Hromádka, Michael R DeWeese, and Anthony M Zador,  “Sparse representation of sounds in the unanesthetized auditory cortex,”  PLoS biology, vol. 6, no. 1, pp. e16, 2008.
  * [25] John Platt,  “Sequential minimal optimization : A fast algorithm for training support vector machines,”  Microsoft Research Technical Report, 1998.


# Plaskett 1.8 metre Observations of Starlink Satellites

Aaron C. Boley Department of Physics and Astronomy, University of British
Columbia
6224 Agricultural Rd
Vancouver, BC V6T 1Z1, Canada  Ewan Wright  ISGP, University of British
Columbia
2357 Main Mall
Vancouver, BC V6T 1Z1, Canada  Samantha Lawler  Campion College and the
Department of Physics, University of Regina
3737 Wascana Parkway, Regina, SK S4S 0A2  Paul Hickson  Department of Physics
and Astronomy, University of British Columbia
6224 Agricultural Rd
Vancouver, BC V6T 1Z1, Canada  Dave Balam  Dominion Astrophysical Observatory
5071 West Saanich Road
Victoria, BC V9E 2E7, Canada

###### Abstract

We present observations of 23 Starlink satellites in the $g^{\prime}$
bandpass, obtained from the Dominion Astrophysical Observatory’s Plaskett 1.8
metre telescope. The targets include a mixture of satellites with and without
brightness mitigation measures (i.e., visors). At the time of the observations
(16 July 2021), Starlink satellites were sunlight throughout the night, and
even with strict elevation and azimuth limits, there were over 800 candidate
Starlink arcs. The satellites altogether have a median absolute brightness
(550 km) of $\overline{H}_{g}^{550}=5.3$ mag. Dividing the targets into those
without and with visors, their median absolute magnitudes are
$\overline{H}_{g}^{550}(no~{}visor)=5.1$ and
$\overline{H}^{550}_{g}(visor)=5.7$ mag, respectively. While the visor sample
is dimmer in aggregate, the absolute brightness distribution ranged from
$H_{g}^{550}=4.3$ mag to 9.4 mag, with the brightest being a visored satellite
and the dimmest a satellite with no mitigation. The intrinsic brightness
dispersion among the full sample is $\sigma_{g}=0.5$ mag.

Artificial Satellites, Astronomical Site Protection

††software: This research was made possible by the open-source projects
astrometry.net (Lang et al., 2010), astropy (Astropy Collaboration et al.,
2018), python (Van Rossum & Drake, 2009), numpy (Harris et al., 2020), and
matplotlib (Hunter, 2007; Droettboom et al., 2016).

##  1 Introduction

Large constellations of satellites will have a substantial impact on
astronomy, particularly wide-field surveys, should they be deployed in their
proposed numbers. These so-called ‘mega-constellations’ or ‘satcons’ seek to
place tens of thousands of communications satellites into low-Earth orbit
(LEO), with the total number of proposed satellites exceeding 100,000.
Starlink (FCC, 2016, 2017, 2020a, 2020b), OneWeb (FCC, 2021), Amazon/Kuiper
(FCC, 2019), and GW/StarNet (Press, 2021) alone could place approximately
65,000 satellites into orbit, although this number is being continuously
amended. While not all of the proposed satellites are expected to be deployed,
other operators are entering the market, making a plausible situation in which
LEO will have 100,000 or more satellites should no restrictions be put into
place.

There are ongoing efforts by astronomers to address the potential for NewSpace
satellite light pollution to interfere with astronomy, including engaging in
dialogue with satcon operators (e.g., SATCON2111https://noirlab.edu/satcon2/
and Dark and Quiet
Skies222https://www.iau.org/news/announcements/detail/ann21043/). Such
dialogue seeks to provide opportunities for exchanges of views and to
establish a basis for cooperation, with identifying voluntary mitigation
efforts for the short term and defining feasible and effective regulations or
guidelines for the protection of the night sky in the long-term. As part of
this effort, it is incumbent upon astronomers to routinely measure satellites
for the purpose of providing independent brightness assessments in
astronomical observing bands and to verify mitigation practices.

Several groups have reported on LEO satellite magnitudes under various
conditions (Tyson et al., 2020; Tregloan-Reed et al., 2021; Mallama, 2021a).
This includes efforts to focus on specific satellites to assess mitigation
attempts, such as ‘DarkSat’ and ‘Visorsats’, representing albedo and
reflection modifications, respectively.

We contribute to this effort by observing 23 Starlink satellites using the
Plaskett 1.8 metre telescope at the Dominion Astrophysical Observatory. The
sample includes a mixture of unmitigated Starlinks and visored Starlinks,
which we analyze together and separately.

##  2 Target Selection and Observing Method

The Two-Line Elements (TLEs) for all Starlink satellites were obtained from
Celestrak333https://celestrak.com/ on July 15, 2021 UTC. Each element was then
run through the JPL Horizons ephemeris
service444https://ssd.jpl.nasa.gov/horizons.cgi using an API to determine the
observing opportunities on July 16, 2021 UTC. An initial selection was made
for Starlink satellites that (1) were visible from DAO during the observing
night, including astronomical twilight, (2) reached a maximum elevation that
was greater than $45^{\circ}$, and (3) reached their maximum elevation in the
southern sky, with this latter being mainly a telescope slewing consideration.
This resulted in over 800 observing opportunities during the approximately
three hours of useful night.

Refining the targets further was motivated largely by the satellite observing
method. The Plaskett 1.8 metre is unable to track LEO objects. Instead, a
‘wait and catch’ methodology is used. In other words, the telescope points at
a field, and tracks it at the sidereal rate. An exposure is then made such
that the satellite of interest is expected to pass entirely through the
telescope’s field of view (FOV) during the image integration.

One of the largest uncertainties in the observing plan is the time needed to
slew between all of the different targets. A conservative approach was taken,
and altogether, we required that the spacing between observations was 12
minutes or more – carrying out this program demonstrated that a slightly
faster cadence would be possible for similar future observing plans. In the
end, the time between the planned observations varied between 12 and 15
minutes, allowing the selection of 23 targets, all of which were successfully
observed.

Among the sample, nine are unmitigated in brightness (i.e., the original
design) and the rest, while having the same overall design, are thought to
have visors. The addition of such visors is a voluntary action by SpaceX in
response to the initial outcry from the astronomical community. Early attempts
to characterize this retrofitted mitigation showed some promise in the
reduction of light pollution (e.g., Mallama, 2021a), but as we will show,
there is significant variability in satellite brightness. For ease of
discussion, we distinguish between the original Starlink design and the
visorsat addition as Starlink and Starlink-V, respectively.

The target satellites and select observing details are listed in Table 1. We,
unfortunately, cannot confirm that the satellites indicated as Starlink-Vs do
indeed have visors or that the visors are correctly deployed. Starlink-1436 is
the first visorsat (e.g., see discussion in Tregloan-Reed et al., 2021),
placed on orbit as a test case and widely discussed on social media.
Additional Starlink satellites were launched without visors until the tenth
Starlink mission, which according to SpaceX (SpaceX, 2020), included
deployable visors for all 57 Starlink satellites on board the 6 August 2020
Falcon 9 launch. Subsequent mission summaries on SpaceX’s website do not
clarify whether the satellites include visors. However, SpaceX personnel have
stated publicly that all satellites after the August launch would indeed be
equipped with visors555see, e.g., the summary at
https://directory.eoportal.org/web/eoportal/satellite-missions/s/starlink..
Public satellite information does not confirm either way, other than
referencing some public statements. Nonetheless, if we take the available
information at face value, then satellites Starlink-1522666This is the lowest
Starlink number associated with the August launch, according to the satellite
catalogue, available at celestrak.com. onward, as well as Starlink-1436, have
deployable visors.

Table 1: Target Information

Name | UTCaaUTC times represent the expected passage of the satellite through the centre of the FOV | RA(J2000) |  $\delta$(J2000) | Azimuth | Elevation | Rate | Range | STO
---|---|---|---|---|---|---|---|---
|  |  |  | deg | deg | $\arcsec/s$ | km | deg
STARLINK-2077(v)bbVisorsat Starlink satellites are denoted with a ’(v)’ | 05:45:11 | 15 28 26.31 | +08 53 52.7 | 216.00 | 45.16 | 2021 | 749.0 | 69
STARLINK-1392 | 05:57:29 | 15 52 24.66 | +13 22 50.4 | 214.71 | 50.41 | 2171 | 697.1 | 66
STARLINK-1747(v) | 06:13:47 | 18 57 30.34 | +16 29 47.3 | 145.59 | 53.97 | 2263 | 668.3 | 39
STARLINK-1728(v) | 06:28:56 | 18 45 17.69 | +27 54 54.4 | 148.44 | 66.96 | 2537 | 595.5 | 51
STARLINK-1355 | 06:41:17 | 17 42 31.87 | +40 37 24.5 | 207.97 | 81.21 | 2705 | 558.0 | 68
STARLINK-1300 | 06:55:00 | 19 35 38.38 | +17 45 39.1 | 145.94 | 55.45 | 2300 | 657.5 | 39
STARLINK-2476(v) | 07:07:22 | 19 59 51.14 | +12 51 33.2 | 144.76 | 49.83 | 2154 | 702.5 | 34
STARLINK-2565(v) | 07:11:25 | 18 35 44.49 | +48 04 10.5 | 154.91 | 89.52 | 2734 | 551.7 | 71
STARLINK-1529(v) | 07:28:02 | 18 13 43.49 | +34 24 40.8 | 209.57 | 74.29 | 2642 | 571.7 | 60
STARLINK-2530(v) | 07:43:20 | 17 26 41.52 | +08 48 58.9 | 216.05 | 45.12 | 2022 | 748.8 | 45
STARLINK-1549(v) | 07:55:38 | 17 51 22.27 | +13 20 05.0 | 214.55 | 50.47 | 2172 | 696.6 | 44
STARLINK-1012 | 08:11:54 | 20 55 40.42 | +16 21 13.5 | 145.74 | 53.91 | 2263 | 668.6 | 42
STARLINK-1009 | 08:28:28 | 20 42 08.85 | +28 54 50.0 | 148.93 | 68.18 | 2558 | 590.7 | 52
STARLINK-1498 | 08:42:05 | 21 41 52.96 | +10 30 13.8 | 143.73 | 47.02 | 2076 | 729.3 | 43
STARLINK-1561(v) | 08:55:55 | 19 24 09.88 | +26 48 08.4 | 211.13 | 65.84 | 2518 | 600.0 | 48
STARLINK-1576(v) | 09:09:34 | 20 32 40.24 | +48 16 44.6 | 192.59 | 89.83 | 2733 | 552.1 | 71
STARLINK-1037 | 09:24:58 | 21 53 23.21 | +23 00 29.2 | 147.19 | 61.48 | 2460 | 614.7 | 55
STARLINK-1060 | 09:39:47 | 22 07 49.16 | +22 57 10.8 | 147.42 | 61.47 | 2436 | 620.7 | 57
STARLINK-2063(v) | 09:55:00 | 20 23 29.08 | +26 28 37.0 | 210.80 | 65.58 | 2518 | 600.2 | 49
STARLINK-1464 | 10:10:00 | 22 54 30.14 | +16 16 48.5 | 145.60 | 53.83 | 2261 | 669.2 | 60
STARLINK-2252(v) | 10:25:17 | 21 12 44.06 | +34 10 07.6 | 208.36 | 74.25 | 2643 | 571.5 | 59
STARLINK-2249(v) | 10:41:41 | 20 56 12.20 | +21 04 56.1 | 212.77 | 59.38 | 2396 | 631.2 | 46
STARLINK-2195(v) | 10:53:10 | 22 51 39.83 | +35 54 22.5 | 149.43 | 75.93 | 2598 | 580.8 | 73

ccfootnotetext: Solar-Target-Observer (STO) angle (essentially the solar phase
angle)

##  3 Observations and Measurements

All observations were taken using a $g^{\prime}$ filter with 30 second
exposures. Each integration began 15 seconds prior to the central UTC time
listed in Table 1 to ensure that the satellite would be observed as it passed
through the Plaskett’s $22^{\prime}\times 11^{\prime}$ (edge-to-edge on the
CCD) FOV.

Images were processed in two ways to check for consistency. The first method
used the Spaceguard pipeline (Tatum et al., 1994), which includes bias
subtraction, flat-fielding, image flattening, and filtering for the removal of
image artefacts (e.g., hot pixels, cosmic rays, etc.). For each image, source
catalogues are generated of all stars in the field using Gaussian fitting. One
of a selection of astrometric catalogues are then queried (Tycho2, PPM,
Carlsberg Meridian Circle, USNO, and 2MASS) and predictions of the detector
positions of each catalogue source are calculated using the world coordinate
system (WCS) information in the image header as a first pass. The two lists
are then cross-correlated, the catalogue source positions are re-calculated
using Gaussian fitting, and the final WCS is determined. All pixels are then
re-sampled to tangential projections with constant scale ($0\arcsec.62$) in
the ICRS reference frame. All re-sampling maintained flux conservation for
each image.

The second method used astropy and python scripts to bias subtract, flat-
field, flatten, and median filter the images. The flattening was done by
subtracting from a given pixel the median value as determined from 128x128
pixel grid. The large grid was needed to avoid forming peak-and-valley
structure in the images along the streaks and the brightest stars. To speed up
this process, a subgrid consisting of 20x20 pixel bins was used, with the
median subtraction done on all subgrid pixels using the large window centred
on the subgrid. Each image was then astrometrically calibrated using the
astrometry.net software package (Lang et al., 2010).

With the astrometrically calibrated images, bright stars were selected that
sample the relevant image area (i.e., in the vicinity of the satellite
streaks) and, to the extent possible, were free of contamination from other
stars or image artefacts. Some fields were very crowded, so the selected stars
were often not the brightest stars in the field. Fortunately, the 30 second
integration times used on each image provided plenty of very high signal-to-
noise ratio (SNR) stars for calibration.

The selected stars were cross-referenced with the Gaia catalogue, and their
$g$ magnitudes were obtained by using the Gaia $G$, $G_{BP}$, and $G_{RP}$
photometry with the SDSS12-Gaia photometric transformation fits for GEDR3,
with a reported uncertainty of $\sigma_{trans}=0.075$ mag.777Available at:
https://gea.esac.esa.int/archive/documentation/GEDR3/Data_processing/chap_cu5pho/cu5pho_sec_photSystem/cu5pho_ssec_photRelations.html
The two different methods for astrometric calibration were self-consistent
with source identification and cross referencing with Gaia positions.

The transformed Gaia magnitudes, $m_{g}^{G}$, were then compared with the
corresponding instrumental magnitudes, $m_{inst}$, as determined through
aperture photometry. The zero point for each star was set according to
$m_{0}=m_{g}^{G}-m_{inst}$. The resulting distribution of zero points were
clipped to remove outliers. The $m_{0}$’s for an image, which included a
minimum of five stars after clipping, were averaged to determine the image’s
zero point. The clip-and-average results were typically within about 0.01 mag
of the median value for the zero points.

We note that we use a $g^{\prime}$ filter, but have calibrated the photometry
using $g$ magnitudes. The difference between the two bands is small ($\lesssim
0.01$; Tucker et al., 2006), much smaller than the uncertainty
$\sigma_{trans}$ in the transformation from Gaia bands to $g$. We therefore
report the magnitudes as $g$ band to reflect the calibration.

In determining the satellite magnitudes, there are several complications that
are introduced by the wait-and-catch method, all of which were anticipated:
(1) A satellite remains in the FOV for only a small fraction of the 30-sec
exposure; (2) the streaks cut across multiple stars and image artefacts; and
(3), the streaks could, in principle, vary in brightness across the image.
Fortunately, we did not see large streak brightness variations, although there
were measured variations on the level of a few tenths of a mag in a few cases.

To address the second issue, we did not use the entire streak. Rather, we
selected as long of a section as practicable that was free of substantial
contamination. In a few cases, several regions were selected and combined.
Altogether, the shortest streak length that we used was $62\arcsec$, while the
longest was about $370\arcsec$. The lowest SNR (just taken to be Poisson
dominated) is 120, which corresponded to the dimmest satellite in the sample.
The typical SNR was a few hundred. The total flux inside a rectangle region
along the selected section of the streak was used for the photometric
measurements, with adjacent regions used to determine the background for
subtraction.

The extrapolated magnitude for each satellite is then calculated using

$m_{g}=\left<m_{0}\right>+m_{inst}-2.5\log_{10}\left(\frac{R}{L}t_{EXP}\right),$
(1)

where $R$ is the sky rate of the satellite as determined from the JPL
ephemerides, $L$ is the measured streak length, and $t_{EXP}$ is the exposure
time, which again is 30 seconds for these observations.

Finally, we define an ‘absolute magnitude’, $H_{g}^{550}$, for the satellites
as the brightness the satellite would be if seen from 550 km. We note that
different papers use a different standard for this, but we use 550 km to
reflect the orbital altitude for the satellite shell in question. To be
verbose,

$H_{g}^{550}=m_{g}-5\log_{10}\left(\frac{r}{550~{}\rm{km}}\right)$ (2)

##  4 Results

The measurements are presented in Table 2. The apparent brightness
distribution varies between $m_{g}=4.7$ mag and 10.0 mag, with the brightest
being a Starlink-V and the dimmest a Starlink satellite.

Table 2 is visualized in Figure 1. The median for the apparent magnitudes is
$\overline{m}_{g}=5.7$ mag, while it is $\overline{H}_{g}^{550}=5.3$ mag for
the absolute distribution. Their respective standard deviations are 1.4 mag
and 1.3 mag. This includes the general trend of the satellites with STO, which
is discussed further below. When the trend is removed, the magnitude variation
is 0.5 mag for the apparent magnitudes. Isolating just the Starlink-V and
Starlink distributions, the apparent magnitude medians are
$\overline{m}_{g}({\rm visor})=6.1$ mag and $\overline{m}_{g}({\rm
no~{}visor})=5.2$ mag, respectively. The absolute magnitudes are likewise
$\overline{H}_{g}^{550}({\rm visor})=5.7$ mag and $\overline{H}_{g}^{550}({\rm
no~{}visor})=5.1$ mag. Thus, the Starlink-V mitigation efforts through visors
do have a measurable effect, with an overall dimming by $\Delta
H_{g}^{550}=0.6$ mag. However, the medians of both distributions are still
bright and are naked-eye visible. Moreover, the substantial spread further
suggests that mitigation needs to be part of the design process, and that
retrofits may have limited impacts as they can still be bright a large
fraction of time.

Table 2: Satellite Brightness Number | Name | Apparent Magnitude ($m_{g}$) | Absolute Magnitude ($H_{g}^{550}$)aaAs seen at 550 km, but without a STO (phase angle) correction.
---|---|---|---
1 | STARLINK-2077(v) | 6.5 | 5.8
2 | STARLINK-1392 | 5.7 | 5.2
3 | STARLINK-1747(v) | 5.3 | 4.8
4 | STARLINK-1728(v) | 5.4 | 5.2
5 | STARLINK-1355 | 5.2 | 5.1
6 | STARLINK-1300 | 4.8 | 4.4
7 | STARLINK-2476(v) | 4.9 | 4.3
8 | STARLINK-2565(v) | 6.5 | 6.5
9 | STARLINK-1529(v) | 6.4 | 6.3
10 | STARLINK-2530(v) | 7.9 | 7.2
11 | STARLINK-1549(v) | 8.2 | 7.7
12 | STARLINK-1012 | 8.5 | 8.1
13 | STARLINK-1009 | 5.2 | 5.0
14 | STARLINK-1498 | 10.0 | 9.4
15 | STARLINK-1561(v) | 5.7 | 5.5
16 | STARLINK-1576(v) | 6.8 | 6.8
17 | STARLINK-1037bbInactive | 6.2 | 5.9
18 | STARLINK-1060 | 5.2 | 4.9
19 | STARLINK-2063(v) | 5.0 | 4.8
20 | STARLINK-1464 | 5.1 | 4.7
21 | STARLINK-2252(v) | 5.4 | 5.3
22 | STARLINK-2249(v) | 4.7 | 4.4
23 | STARLINK-2195(v) | 6.9 | 6.7

The satellite distributions can further be shown against their STO, range, and
time (Fig. 2). There is a general trend of dimming with increasing STO. This
behaviour is non-trivial in that a satellite passing near the observer’s
zenith will have an STO that approaches $90^{\circ}$, but adjusted for the
Sun’s angular distance below the horizon. This is also when the satellite is
closest to the observer, which by itself would tend to make the satellite
brighter. However, a shallower STO is able to compensate for the difference in
range, up to a point, which places the shallower STOs here among the brightest
observed satellites. Of course, the STO is not unique to a given range, which
can be seen in Figure 2. For the restrictions placed on our observing
elevation and azimuth, the brightest satellites occur at a range of about 650
km. These same restrictions set the limits on the range of STOs observed.

A surprising feature in the data is the dimming of the satellites at an STO of
about $43^{\circ}$. This dimming cannot be explained by range. The time of
observation also does not show clear evidence alone for this feature given
that the dimmest two satellites are separated by one of the brightest.
Moreover, there is not a significant change in the image zero points to
suggest that there is a local sky phenomenon that is compromising some of the
measurements.

Figure 1: Magnitude distributions for the Starlink satellites observed in this
study. Each satellite was observed once only. Left: apparent $m_{g}$
magnitudes, with blue showing the contribution to each bin from Starlinks and
red the Starlink-Vs. Right: similar to the left panel, but for the absolute
magnitude, as defined by seeing the satellite at a range of 550 km. No STO
correction is applied. Because some values fall on bin edges, we use the
interval $[a,b)$ for edges $a$ and $b$ when assigning a magnitude to a bin.

Figure 2:  Select observational properties for the observed satellites,
emphasizing Starlinks (square) Starlink-Vs (circle). A slight dimming of the
satellites is seen with increasing STO. A sudden change in brightness occurs
around $43^{\circ}$ for both Starlink and Starlink-V. This feature is not
understood. The brightest satellites are seen at a range near 650 km due to
the combination of distance and favourable STO. The satellites with the
smallest possible range (550 km), due to the orbital altitude, occur near an
STO of about $70^{\circ}$. Higher STOs were not observed due to restrictions
on telescope azimuth for slew considerations.

##  5 Discussion

Characterizing the brightness distribution of satellites, their variability,
and dependence on STO (solar phase angle) is necessary for assessing the
degree to which satellites will interfere with astronomy and stargazing. It is
further required for modelling the impacts of proposed satellite systems
(e.g., Lawler et al., 2021), as well as monitoring mitigation efforts and, to
the extent practicable, identifying possible regulation and industry best
practices that would have the greatest impact on reducing light pollution
while allowing the development of orbital infrastructure.

The observations here show that the satellites have significant scatter in
brightness, which is not improved when correcting to a range of 550 km. This
is not surprising given the range of STOs among the satellites. We thus can
consider different phase models, with the simplest being the diffuse sphere
(Pradhan et al., 2019). Specifically, the model $g$ magnitude is

$m_{g}=-26.47-2.5\log_{10}\left\\{\frac{2}{3\pi^{p+1}}\zeta\left[\left(\pi-{\rm
STO}\right)\cos\left({\rm STO}\right)+\sin\left({\rm
STO}\right)\right]^{p}\right\\}+5\log_{10}\left(R_{\rm m}\right).$ (3)

where $\zeta$ is the effective area of the satellite in square metres (i.e.,
the product of the albedo and cross-sectional area), $R_{\rm m}$ is the range
in metres, and the index $p=1$.

The result is shown in Figure 3. We use $\zeta=0.7\rm~{}m^{2}$, which provides
a reasonable match to the low STO satellites, but does not dim fast enough for
higher STOs. Noticing that a faster dependency is required, we explore
possible variations of the diffuse sphere model in which $p\neq 1$.

A Monte Carlo Markov Chain routine is written in python using a Metropolis-
Hastings algorithm with a Gibbs sampler (sufficient for this application),
using uniform priors to explore possible $p$ and $\zeta$. We assume the data
have uncertainties of $0.1$ mag and exclude the four peculiarly faint
satellites at ${\rm STO}\approx 43^{\circ}$. The result is $\zeta\approx
1.1\rm~{}m^{2}$ and $p\approx 3.1$. We hesitate to place uncertainties on
these very approximate relations. The results are nonetheless shown in Figure
3.

The heuristic model magnitudes are subtracted from the data (Fig. 3),
resulting in a dispersion of 0.5 mag among the residuals (excluding again the
four dimmest satellites). We caution against reading too much into this model
at this time, as additional data are needed to test it. Moreover, we only
cover a small range of STOs.

With this in mind, an additional behaviour can be noted. First, the Starlink-
Vs appear to fall off in brightness slightly more rapidly with STO than the
Starlink satellites. This is consistent with the differences in the median
magnitudes for the two populations, but also shows that so far the
effectiveness of Starlink-Vs is highly variable and does not help with some of
the brightest orientations.

Recent work by Mallama (2021b) combines about 10,000 MMT-9 robotic
observations of Starlink-Vs using a clear filter. They report that the
magnitudes are within about 0.1 mag of the Johnson V band. The median of their
Starlink-V observations corrected to 550 km is 5.5 $V_{\rm MMT}$ mag.888The
subscript is used to denote that the observations are based in the clear MMT-9
filter and transformed to an approximate $V$ band. When we take into account
the $g-V=0.3$ colour of the Sun, assuming no wavelength dependence for the
albedo, their median $\overline{H}_{g}^{MMT}=5.8$ mag compared with our
$\overline{H}_{g}^{550}=5.7$ mag, consistent to within uncertainties.

Mallama (2021b) also find that a quadratic fit to the STO-brightness
dependence can explain their observations. We note that the quadratic fit
behaves nearly identically to a diffuse Lambertian sphere over the range shown
in Figure 3 (this work), but with the leading term set to 3 instead of their
nominal value. That correspondence means our data exhibits a slightly
different dependence on STO than the MMT-9 data. The origin of this is
unclear, and we refrain from discussing it further other than to note that
because the satellites are not actually spheres, the detailed viewing angle
will matter. Regardless, this study and Mallama (2021b) find an inherent
satellite magnitude dispersion of about 0.5 mag after subtracting their
respective STO trend. We conclude that this is an inherent dispersion due to
the shape of the satellite and that the visors do not prevent this from
occurring.

The results can be used in modelling the impacts of satellites on the night
sky, as done in Lawler et al. (2021), at least for Starlink-like satellites.
As additional photometric data are acquired in astrophysical bands, we will be
able to ensure higher fidelity between observations and models, at least in a
statistical sense.

The scatter demonstrates that analysis of mitigation efforts must be done
statistically, with the need to consider observations that span a wide range
of STOs and orientations of the spacecraft relative to the observer.

These observations only measure the brightness of Starlink satellites in a
single filter. Additional observations in different bandpasses are needed, as
well as measurements of different satellite designs (e.g., OneWeb). As part of
seeking cooperation with satcon operators toward minimizing the impact that
their satellites have on the sky, we need to be vigilant with independent
verification of their mitigation efforts. We also need to consider assessments
of single satellite magnitudes, while investigating cumulative impacts when
practicable. Astronomy is a fundamental way to explore space, test our
understanding of physical laws, and detect impact hazards, among many other
things – interference with astronomy has wide ramifications for science,
education, and safety.

The observational planning scripts, calibrated images, associated photometry,
and analysis scripts used for this paper are available at
https://github.com/norabolig/Boleyetal-Plaskett-2021. The astropy simple
reduction tools are available at https://github.com/norabolig/ABIRL.

The University of British Columbia is situated on the traditional, ancestral,
and unceded territory of the Musqueam people. The University of Regina is
located on Canadian Treaty 4 land, which is the traditional territories of the
nêhiyawak, Anihšināpēk, Dakota, Lakota, and Nakoda, and the homeland of the
Métis/Michif Nation. This work is based on observations obtained at the
Dominion Astrophysical Observatory, NRC Herzberg, Programs in Astronomy and
Astrophysics, National Research Council of Canada. It also used the facilities
of the Canadian Astronomy Data Centre operated by the National Research
Council of Canada with the support of the Canadian Space Agency. This research
was supported in part by NSERC Discovery Grants RGPIN-2020-04635 (ACB) and
RGPIN-2020-04111 (SML), as well as the Canada Research Chairs program.

Figure 3: Similar to Figure 2, but including a diffuse sphere model for the
apparent magnitude (red, downward pointing triangle) and a heuristic power law
STO (phase) model (blue, upward pointing triangle). The right panel shows the
difference between the observed magnitudes and the predictions from the
heuristic phase model.

## References

  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Droettboom et al. (2016) Droettboom, M., Hunter, J., Caswell, T. A., et al. 2016, matplotlib: matplotlib v1.5.1, doi: 10.5281/zenodo.44579
  * FCC (2016) FCC. 2016, IBFS, SAT-LOA-20161115-00118, Attachment A
  * FCC (2017) —. 2017, IBFS, SAT-LOA-20170301-00027, Attachment A
  * FCC (2019) —. 2019, IBFS, SAT-LOA-20190704-00057, Technical Appendix
  * FCC (2020a) —. 2020a, IBFS, SAT-MOD-20200417-00037, Attachment A
  * FCC (2020b) —. 2020b, IBFS, SAT-LOA-20200526-00055, Attachment A
  * FCC (2021) —. 2021, IBFS, SAT-MPL-20210112-00007, Attachment A
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Lang et al. (2010) Lang, D., Hogg, D. W., Mierle, K., Blanton, M., & Roweis, S. 2010, AJ, 139, 1782, doi: 10.1088/0004-6256/139/5/1782
  * Lawler et al. (2021) Lawler, S. M., Boley, A. C., & Rein, H. 2021, arXiv e-prints, arXiv:2109.04328.  https://arxiv.org/abs/2109.04328
  * Mallama (2021a) Mallama, A. 2021a, arXiv e-prints, arXiv:2101.00374.  https://arxiv.org/abs/2101.00374
  * Mallama (2021b) —. 2021b, arXiv e-prints, arXiv:2109.07345.  https://arxiv.org/abs/2109.07345
  * Pradhan et al. (2019) Pradhan, B., Hickson, P., & Surdej, J. 2019, Acta Astronautica, 164, 77, doi: 10.1016/j.actaastro.2019.07.008
  * Press (2021) Press, L. 2021, A New Chinese Broadband Internet Constellation.  https://circleid.com/posts/20201002-a-new-chinese-broadband-satellite-constellation
  * SpaceX (2020) SpaceX. 2020, STARLINK MISSION, August 06, 2020.  https://www.spacex.com/updates/starlink-mission-08-06-2020/index.html
  * Tatum et al. (1994) Tatum, J. B., Balam, D. D., & Aikman, G. C. L. 1994, Planet. Space Sci., 42, 611, doi: 10.1016/0032-0633(94)90036-1
  * Tregloan-Reed et al. (2021) Tregloan-Reed, J., Otarola, A., Unda-Sanzana, E., et al. 2021, A&A, 647, A54, doi: 10.1051/0004-6361/202039364
  * Tucker et al. (2006) Tucker, D. L., Kent, S., Richmond, M. W., et al. 2006, Astronomische Nachrichten, 327, 821, doi: 10.1002/asna.200610655
  * Tyson et al. (2020) Tyson, J. A., Ivezić, Ž., Bradshaw, A., et al. 2020, AJ, 160, 226, doi: 10.3847/1538-3881/abba3e
  * Van Rossum & Drake (2009) Van Rossum, G., & Drake, F. L. 2009, Python 3 Reference Manual (Scotts Valley, CA: CreateSpace)