accretion, accretion disks — black hole physics — radiative transfer – X-rays:
general

# What Determines Unique Spectra of Super-Eddington Accretors?: Origin of
Optically Thick and Low Temperature Coronae in Super-Eddington Accretion Flows

Norita Kawanaka11affiliation: Department of Astronomy, Kyoto University,
Kitashirakawa-Oiwake-cho,
Sakyo-ku, Kyoto 606-8502, Japan 22affiliation: Hakubi Center, Kyoto
University, Yoshida-Honmachi,
Sakyo-ku, Kyoto 606-8501, Japan and Shin Mineshige11affiliationmark:
<EMAIL_ADDRESS>
###### Abstract

Existence of relatively cool ($k_{B}T\lesssim 10~{}{\rm keV}$) and optically
thick ($\tau\gtrsim 3$) coronae are inferred above super-Eddington accretion
flow such as ultraluminous X-ray sources (ULXs), GRS 1915+105, and narrow-line
Seyfert 1 galaxies (NLS1), which contrasts the cases in sub-Eddington
accretion flows, which are associated with coronae with $k_{B}T\sim 100~{}{\rm
keV}$ and $\tau\sim 1$. To understand their physical origin, we investigate
the emission properties of the corona which is formed by the gas blown off the
super-Eddington inner disk by radiation pressure. We assume that the corona is
heated by the reconnection of magnetic loops emerged from the underlying disk.
We show that this radiation pressure driven wind can act as an optically thick
corona which upscatters thermal soft photons from the underlying disk, and
that with a reasonable parameter set we can theoretically reproduce the
coronal optical depth and temperature which are inferred by spectral fittings
of observational data. By contrast, the coronal optical depth cannot be so
high in sub-Eddington cases, since the coronal material is supplied from the
disk via evaporation and there is a maximum limit on the evaporation rate. We
support that the low temperature, optically thick Comptonization should be a
key signature of super-Eddington accretion flow.

##  1 Introduction

Recently, the existence of super-Eddington (or supercritical) accretion flows
onto black holes have been investigated from both observational and
theoretical points of view (see, e.g., [Kato et al. (2008)]). For example,
ultraluminous X-ray sources (ULXs) are very bright X-ray sources whose X-ray
luminosity is much larger than the Eddington luminosity of a stellar-mass
black holes $L_{X}\gtrsim 10^{39}{\rm erg}~{}{\rm s}^{-1}$, and two
possibilities have been presented to account for ULXs: sub-Eddington accretion
onto an intermediate mass black hole (IMBHs; [Colbert & Mushotzky (1999),
Makishima et al. (2000), Miller et al. (2003), Miller et al. (2004), Cropper
et al. (2004)]) and super-Eddington accretion onto a stellar-mass black hole
([Watarai et al. (2001), King et al. (2001), Poutanen et al. (2007)]). In
fact, some studies show that their spectral features can be well explained by
the model with super-Eddington accretion flows ([Ebisawa et al. (2003),
Okajima et al. (2006), Tsunoda et al. (2006), Vierdayanti et al. (2006),
Vierdayanti et al. (2008), Gladstone et al. (2009), Kawashima et al. (2012),
Sutton et al. (2013), Liu et al. (2013)]). In addition, Narrow-Line Seyfert 1
galaxies (NLS1s) are considered to harbor a black hole with a super-Eddington
accretion flow (e.g. [Zhou & Wang (1999), Mineshige et al. (2000), Collin &
Kawaguchi (2004), Jin et al. (2017)]). As for Galactic X-ray sources, some
microquasars (e.g. GRS 1915+105; [Done et al. (2004), Done et al. (2007),
Vierdayanti et al. (2010b)]) are the candidates of super-Eddington accreting
black holes. The detailed observational data of these sources would help us
interpret comprehensively the behavior of super-Eddington accretion flows.

Despite extensive spectral analyses of potential super-Eddington sources the
presence of super-Eddington accretion flow is still a controversial issue,
since there is no good theory to predict its firm observational signatures. We
definitely need information of black hole masses. One exceptional case is GRS
1915+105, for which the black hole mass is known to be $14\pm 4M_{\odot}$
(Greiner et al., 2001). This is the best target in our Galaxy to study high
luminosity states of accretion (Done et al., 2007). The other exceptional
cases are ULX pulsars, whose central objects are accreting neutron stars and
thus less massive than $3M_{\odot}$. To explain their luminosities ($\gtrsim
10^{40}~{}{\rm erg}~{}{\rm s}^{-1}$) super-Eddington accretion is required,
and eight ULX pulsars have been discovered so far (Bachetti et al., 2014;
Fürst et al., 2016; Israel et al., 2017, 2017; Carpano et al., 2018; Wilson-
Hodge et al., 2018; Rodríguez Castillo et al., 2020; Sathyaprakash et al.,
2019; Chandra et al., 2020).

As for the theoretical approach, simple analytical or semi-analytical modeling
have been attempted under one-zone approximation (e.g. Shakura & Sunyaev
(1973); Abramowicz et al. (1988); Honma et al. (1991); Poutanen et al.
(2007)). To study the gas dynamics of super-Eddington accretion flow, global,
multi-dimensional, radiation-hydrodynamic (RHD) simulations are indispensable.
Such simulations have been attempted by several authors (Eggum et al., 1988;
Okuda, 2002; Okuda et al., 2005; Ohsuga et al., 2005). Especially Ohsuga et
al. (2005) performed large simulation box, long-term RHD simulations of super-
Eddington accretion flows for the first time, and have demonstrated that
super-Eddington accretion onto a black hole in a quasi-steady fashion is
feasible. Moreover, they have shown that such super-Eddington disks naturally
give rise to high-velocity ($v\sim 0.1-0.2c$) outflows driven by radiation
pressure111Note that the possibility of such super-Eddington outflows was
already suggested by Shakura & Sunyaev (1973).. Importantly, there is no upper
limit to the mass accretion rate to a black hole nor to the mass outflow rate;
the mass accretion rate and outflow rate are roughly in proportion to the mass
supply rate (Ohsuga, 2007). Note, however, these studies adopted the
$\alpha$-prescription for the viscosity in the accretion flows. This
shortcoming was improved by Ohsuga et al. (2009), who have presented for the
first time global two-dimensional radiative magnetohydrodynamic simulations of
black hole accretion flows following the angular momentum transports via the
magnetic torque, and these simulations have also shown the existence of the
high-velocity disk outflows driven by radiation-pressure force in super-
Eddington accretion flows (see also Ohsuga & Mineshige (2011)).

The outflow associated with a super-Eddington accretion flow may be optically-
thick with respect of Thomson scattering for sufficiently high $\dot{M}>L_{\rm
Edd}/c^{2}$, where $L_{\rm Edd}$ is the Eddington luminosity, and electrons in
it would Compton-upscatter soft photons from the underlying accretion flow.
The plasma surrounding an accretion disk and upscattering soft photons from
the disk is called a corona and is considered to exist universally in
accretion disks associated with AGNs and X-ray binaries that are emitting hard
X-rays (Bisnovatyi-Kogan & Blinnikov, 1976, 1977; Sunyaev & Truemper, 1979),
and its temperature is supposed to be $k_{\rm B}T\sim 100~{}{\rm keV}$ as a
result of energy dissipation of magnetic fields amplified in the accretion
disk (Haardt & Maraschi, 1991). Kawashima et al. (2009) showed by the RHD
simulation that the Compton $y$ parameter of the radiation pressure-driven
outflow from a super-Eddington accretion flow is around unity. That is,
occurrence of significant Comptonization was expected. Later, Kawashima et al.
(2012) calculated the radiation spectra of super-Eddington accretion flows
using a Monte Carlo method and showed that the radiation pressure-driven
outflow acts as a corona above the disk and should Compton upscatters soft
photons from the underlying disk, thereby modifying the X-ray spectra from
super-Eddington accretion flows (see Kitaki et al. (2017) for a more recent
work). On the other hand, Gladstone et al. (2009) did the observational
studies about ULXs and found that most of their X-ray spectra can be fitted
well with the thermal component being Compton upscattered by relatively low-
temperature ($k_{\rm B}T\lesssim 10{\rm keV}$) and optically-thick ($\tau_{\rm
sc}\gtrsim 3$) plasma. Similar spectral fitting is also performed for other
ULXs (Yoshida et al., 2013; Shidatsu et al., 2017) and for GRS 1915+105
(Vierdayanti et al., 2010b), whose mass accretion rate is considered to be
super- or near-Eddington as stated above, and for NLS1s (Idogaki et al.,
2018). Therefore, these facts support the supercritical accretion flow model
for ULXs. However, the physical interpretation of those low temperature and
large optically depth of that corona has not been described so far.

In this study, we investigate the basic properties of outflowing coronae in
supercritical accretion flows with a simple, semi-analytical model. Especially
we aim at understanding the physical reason for the low temperature and large
optical depth of coronae associated with some ULXs, microquasars, and NLS1s.
It is obvious that the scattering optical depth can be higher for larger mass
outflow rates, though it remains a question how and why electron temperatures
should be low. In the present study we raise the following questions; why
electron temperatures tend to be low in the outflowing coronae from
supercritical accretion flows; why $\tau$ of the outflowing corona cannot be
so high in subcritical flows. We discuss these issues from the viewpoint of
energy balance in the outflowing corona. In our calculation we take into
account the heating due to magnetic reconnections and cooling due to inverse
Compton scatterings in radiation pressure-driven outflows in a simple and
self-consistent manner.

##  2 Model

In the present paper, we apply the basic concept of the magnetic reconnection
heated corona model by Haardt & Maraschi (1991) to the case of super-Eddington
accretors. The big distinction lies in that disk coronae are formed not via
evaporation of disk gas but by radiation-pressure driven outflow.

###  2.1 Coronal Density

We consider a super-Eddington accretor with $\dot{M}>L_{\rm Edd}/(\eta
c^{2})$, where $\eta\sim 0.1$ is the radiation efficiency. Emergence of a
strong disk wind is expected from the innermost region of the supercritical
disk, being driven by radiation-pressure force by the disk (Eggum et al.,
1988). Due to this mass loss process, the mass accretion rate in the disk is
effectively suppressed, and generally depends on the radius. It is often
assumed that the local accretion rate in an outflowing accretion disk varies
as a power law in radius,

$\displaystyle\dot{M}(r)=\left\\{\begin{array}[]{ll}\dot{M}_{0}&r>r_{\rm
crit}\\\ \dot{M}_{0}\left(\frac{r}{r_{\rm crit}}\right)^{s}&r\leq r_{\rm
crit}.\end{array}\right.$ (3)

Here $s$ is a constant less than unity and the critical radius $r_{\rm crit}$
is determined as (Ohsuga et al., 2005)

$\displaystyle r_{\rm crit}=\frac{3\dot{M}_{0}c^{2}}{4L_{\rm Edd}}r_{\rm S},$
(4)

where $r_{\rm S}=2GM_{\rm BH}/c^{2}$ is the Schwarzschild radius, and
$\dot{M}_{0}$ is the mass accretion rate at the outer edge. Assuming that the
wind velocity is on the same order of the escape velocity $v_{\rm
esc}=(2GM_{\rm BH}/r)^{1/2}$, the coronal density distribution $n_{\rm
cor}(r)$ should satisfy

$\displaystyle m_{\rm p}n_{\rm cor}v_{\rm esc}4\pi rdr=d\dot{M}(r),$ (5)

which gives

$\displaystyle n_{\rm cor}(r)=\frac{s}{4\pi m_{\rm p}rv_{\rm
esc}}\frac{\dot{M}(r)}{r},$ (6)

where $m_{\rm p}$ is the proton mass. The scattering optical depth of an
outflowing corona can be obtained as

$\displaystyle\tau_{\rm cor}(r)=n_{\rm cor}(r)\sigma_{T}\ell_{\rm cor},$ (7)

where $\ell_{\rm cor}$ is the scale height of the outflowing corona, which
will be given in the following section.

###  2.2 Optical Depth and Temperature of a Corona

In this subsection, we present the heating process for the outflowing corona
above a supercritical accretion disk. In the context of a subcritical
accretion flow, Haardt & Maraschi (1991) proposed a corona model in which the
corona is heated by the reconnection of magnetic loops emerging from the
rotating disk, and here we discuss the energy balance of an outflowing corona
following their idea. As for the strength of magnetic field, here we assume
that its energy density in the disk is the constant fraction of the radiation
energy density in the disk. Then the magnetic field in the disk can be
determined by

$\displaystyle\frac{B^{2}}{8\pi}\approx\eta_{B}\cdot{\rm max}\left(aT_{\rm
disk}^{4},\frac{\dot{M}\Omega}{2\pi r\alpha}\right),$ (8)

where $\eta_{B}<1$ is a parameter which should be given from the results of
radiative magnetohydrodynamic simulations (e.g. Ohsuga & Mineshige (2011)) in
the following discussion. Here $aT_{\rm eff}^{4}$ and
$\frac{\dot{M}\Omega}{2\pi r\alpha}$ are the radiation pressure and gas
pressure (Kawanaka et al. (2013)). Here $T_{\rm disk}$ is the temperature at
the equatorial plane, which is given from the standard model of a radiation
pressure-dominated disk (Shakura & Sunyaev, 1973) as

$\displaystyle T_{\rm disk}=2.1\times 10^{7}{\rm K}\left(\alpha_{-1}M_{{\rm
BH},1}\right)^{-1/4}r_{1}^{-3/8},$ (9)

where $M_{{\rm BH},1}=M_{\rm BH}/10M_{\odot}$, $\alpha=10^{-1}\alpha_{-1}$ is
the viscous paramter, and $r=10r_{\rm S}r_{1}$ is the radial distance from a
black hole. We note that this temperature is independent of the mass accretion
rate. Then we can estimate the magnetic field strength in the case with a
radiation-pressure dominated disk:

$\displaystyle B\approx 1.9\times 10^{7}~{}{\rm
Gauss}~{}\eta_{B,-2}^{1/2}\alpha_{-1}^{-1/2}M_{{\rm
BH},1}^{-1/2}r_{1}^{-3/4},$ (10)

where $\eta_{B,-2}=\eta_{B}/10^{-2}$.

Following Liu et al. (2002b), we assume that the energy equilibrium of the
outflowing corona is established between the magnetic reconnection heating and
the inverse Compton cooling. The vertical scale height of an outflowing corona
is determined by the loop length of the magnetic field emerged from the
underlying disk, $\ell_{\rm loop}$. Here we assume that $\ell_{\rm loop}$ has
the same order as the radial position of the spot where the magnetic loop
emerged: $\ell_{\rm loop}\simeq r$. In the innermost region, however, the loop
should remain long enough due to the existence of a funnel structure in the
outflow (Kitaki et al., 2017), and so we set the floor value of the loop
length as $\ell_{\rm loop}={\rm max}(r,10r_{\rm S})$. In addition, we should
take into account the fact that the corona above the super-Eddington flow is
not static but is upwardly moving. Photons can escape the outflowing corona
when their diffusion timescale, $\ell_{\rm cor}\tau_{\rm cor}/c$, is shorter
than the escape timescale for the corona, $\ell_{\rm cor}/v_{\rm esc}$, which
gives the lower limit for the scale height of the corona as

$\displaystyle\ell_{\rm cor,lim}\approx\frac{c}{v_{\rm esc}n_{\rm
cor}\sigma_{T}}.$ (11)

Hereafter we estimate the coronal scale height as $\ell_{\rm cor}\approx{\rm
min}(\ell_{\rm cor,lim},\ell_{\rm loop})$.

Then the coronal heating rate per unit surface area due to magnetic
reconnection can be described as

$\displaystyle Q^{+}=\frac{B^{2}}{4\pi}V_{R}\cdot\frac{\ell_{\rm
cor}}{\ell_{\rm loop}},$ (12)

where $V_{R}$ is the speed of energy dissipation via magnetic reconnection.
Here we simply assume that $V_{R}={\rm min}(V_{A},c)$ where $V_{A}=B/(4\pi
n_{\rm cor}m_{\rm p})^{1/2}$ is the Alfvén velocity. On the other hand, the
coronal cooling is dominated by inverse Compton scattering of seed photon
field supplied from the underlying disk:

$\displaystyle Q^{-}=\frac{4k_{B}(T_{\rm cor}-T_{\rm rad})}{m_{e}c^{2}}cU_{\rm
rad}\cdot{\rm max}(\tau_{\rm cor},\tau_{\rm cor}^{2}),$ (13)

where $U_{\rm rad}$ is the energy density of the soft photon field to be
Comptonized in the corona (Rybicki & Lightman, 1979). Here we take into
consideration that the temperature of coronal electrons $T_{\rm cor}$ may not
be much higher than that of soft photon field, $T_{\rm rad}=(U_{\rm
rad}/a)^{1/4}$, and that we cannot neglect the term $4k_{B}T_{\rm rad}$ in the
numerator of the right-hand side. Since the radiation energy density from a
super-Eddington accretion flow is limited by the local Eddington flux, we can
assume the soft photon field emitted intrinsically from the underlying disk as

$\displaystyle U_{\rm rad,intr}$ $\displaystyle=$ $\displaystyle\frac{L_{\rm
Edd}}{4\pi r^{2}c}\simeq 3.7\times 10^{12}~{}{\rm erg}~{}{\rm cm}^{-3}M_{{\rm
BH},1}^{-1}r_{1}^{-2}.$ (14)

As the soft photon field, we should also take into account the reprocessed
radiation, $U_{\rm rad,re}$, that originates from the backward Compton
scattering in the corona,

$\displaystyle U_{\rm rad,re}=F_{\rm cor,irr}/c,$ (15)

where $F_{\rm c,irr}$ is the irradiation flux from the corona onto the
underlying disk. Here the disk albedo is assumed to be unity. The soft photon
field is reevaluated as $U_{\rm rad}=U_{\rm rad,intr}+U_{\rm rad,re}$, which
should have the Planck distribution. By solving $Q^{+}=Q^{-}$ at each radius,
we can evaluate the radial distribution of coronal temperature. We should note
that $F_{\rm c,irr}$ can be obtained only after solving the radiation transfer
inside the corona by Monte Carlo simulations (see Sec.3.1) for a calculated
temporal temperature distribution of the corona. Therefore, to obtain the
coronal temperature distribution consistently, we should repeat performing
Monte Carlo simulations and updating the coronal temperature until it
converges.

Let us estimate the typical optical depth of the wind corona in our model.
Here we adopt $s=0.15$, which is implied by recent radiative hydrodynamic
simulations of super-Eddington accretion flows Kitaki et al. (2018). The
coronal density as a function of radius, black hole mass, and mass accretion
rate can be then evaluated as

$\displaystyle n_{\rm cor}(r)\simeq 1.2\times 10^{17}~{}{\rm
cm}^{-3}r_{1}^{-1.35}M_{{\rm BH},1}^{-1}\dot{M}_{1}^{0.85},$ (16)

where $\dot{M}_{1}=\dot{M}_{0}/(10\dot{M}_{\rm Edd})$, and $\dot{M}_{\rm
Edd}=L_{\rm Edd}/\eta c^{2}$ is the Eddington accretion rate, where $\eta=0.1$
is the radiation efficiency. The optical depth can be then described as

$\displaystyle\tau_{\rm cor}\simeq
2.4r_{1}^{-0.35}\dot{M}_{1}^{0.85}\left(\frac{\ell_{\rm cor}}{r}\right),$ (17)

which can be larger than unity when the mass accretion rate is well above the
Eddington value. This is consistent with the observations of ULXs and other
super-Eddington accreting black holes.

##  3 Results

###  3.1 Coronal Structure

In our model, soft thermal photons are emitted from an optically-thick disk,
and some of them are Compton upscattered by hot electrons in the outflowing
corona lying above the disk. As a result, the radiation spectrum emerged from
such a system would have a hard component in addition to the thermal component
originating from a disk.

We use Monte Carlo simulations to calculate the radiation spectrum expected
from our model of a wind-fed disk corona, based on the prescription presented
in Pozdnyakov et al. (1977). As mentioned in the previous subsection, the
coronal temperature distribution and radiation spectrum are obtained
iteratively. First, we assume at each radius the seed photon field, $U_{\rm
rad}$ as blackbody radiation whose energy density is given by Eq.(14). Second,
we compute the coronal temperature and optical depth at that radius as
described in the previous section. Third, we solve the radiation transfer in a
corona by Monte Carlo simulations. Here we adopt the vertically one-zone
approximation for the corona at each radius. Especially, we take into account
the bulk motion of a corona, whose upward velocity is given by the escape
velocity at that radius. Fourth, we obtain the downward flux from the corona
(i.e., the irradiation flux onto the underlying disk, $F_{\rm cor,irr}$), and
evaluate the reprocessed photon field by Eq.(15). Fifth, we compare the
reprocessed photon field, $U_{\rm rad,re}$, and the photon field used in the
simulations, $U_{\rm rad}$. If the former is comparable to or larger than the
latter, we should reevaluate the soft photon field as $U_{\rm rad}=U_{\rm
rad,intr}+U_{\rm rad,re}$, and go back to the second step to perform the Monte
Carlo simulations with renewed coronal temperature. If the former is much
smaller than the latter at each radius, we regard that the consistent solution
is obtained, and compute the radiation luminosity by integrating the energy
spectra of photons emerged from the upper boundary of coronae for all radii.

Fig. 1 depicts the temperature of the outflowing corona as a function of
radius. We can see that the coronal temperatures is roughly several tens of
${\rm keV}$, and we see a tendency that the higher the mass accretion rate is,
the cooler becomes the coronal temperature. Note that the outflow (corona)
solution exists only in the inner regions, and that the higher the accretion
rate is, the larger becomes the cutoff radius.

We see several breaks in the slope of each curve (see the arrows with numbers
attached to the curve of $\dot{M}=10\dot{M}_{\rm Edd}$). In the innermost
region (i.e., inside the point (1) for $\dot{M}=10\dot{M}_{\rm Edd}$), the
coronal scale height is determined by the diffusion length of a photon within
the dynamical timescale, $\ell_{\rm cor,lim}$ (see Eq. 9). At the break point
(1), $\ell_{\rm cor,lim}$ is equal to the loop length, $\ell_{\rm
loop}=10r_{\rm S}$, and the coronal scale height is determined as $\ell_{\rm
cor}=\ell_{\rm loop}=10r_{\rm S}$ in the region between the points (1) and
(2). The magnetic loop length is switched from $10r_{\rm S}$ to $r$ at the
break point (2), and the coronal scale height is determined as $\ell_{\rm
cor}=\ell_{\rm loop}=r$ in the region between the points (2) and (3). Finally,
the break point (3) corresponds to the radius where the radiation pressure is
equal to the gas pressure: outside this point the radiation pressure would
dominate over the gas pressure.

Fig. 2 depicts the scattering optical depth of the outflowing corona as a
function of radius. In the innermost region, the coronal optical depth
increases with radius as $\propto r^{1/2}$. This dependence can be explained
in the following way: in the innermost region, as shown above, the scale
height of the outflowing corona is determined as $\ell_{\rm cor,lim}$ in Eq.
9. Therefore, the optical depth is evaluated as $n_{\rm
cor}\sigma_{T}\ell_{\rm cor,lim}\approx c/v_{\rm esc}\propto r^{1/2}$. Note
that this value is independent of the mass accretion rate, which can also be
seen in Fig. 1. The optical depth peaks at the point (2) where the loop length
$\ell_{\rm loop}$ is equal to $\ell_{\rm cor,lim}$. In the region between (1)
and (2), except the case of $\dot{M}=50\dot{M}_{\rm Edd}$, $\tau_{\rm cor}$
decreases with radius steeply, $\propto r^{s-3/2}$, because the coronal scale
height is evaluated as $\ell_{\rm cor}=\ell_{\rm loop}=10r_{\rm S}$ in this
region. In the outermost region, since $\ell_{\rm cor}=r$, $\tau_{\rm cor}$
decreases with radius shallowly, $\propto r^{s-1/2}$. When the mass accretion
rate is large enough, $\ell_{\rm cor,lim}$ is always smaller than the floor
value of the scale height, $10r_{\rm S}$, the middle region disappears. In
contrast with the case of sub-Eddinton accretion, where the coronal optical
depth is around unity, the peak value of the optical depth of outflowing
corona is larger than unity, $1\lesssim\tau_{\rm cor}\lesssim 10$.

Fig. 3 depicts the Compton $y$-parameter of the outflowing corona as a
function of radius. We can see that in the inner region where the outflowing
corona is sufficiently optically-thick, the $y$-parameter is the order of
unity, which is similar to the case of sub-Eddington accretion flows.

We should note that the energy conservation is not violated in our model. Fig.
4 depicts the fraction of disk energy, $3GM_{\rm BH}\dot{M}/(8\pi r^{3})$,
which is dissipated in the corona via magnetic reconnection for various
accretion rates, assuming $\eta_{B}=7.5\times 10^{-3}$. We can see that the
ratio never exceeds unity at any accretion rate and radius, which means that
the energy conservation is not violated in the calculations. Note that this
fraction should be around unity throughout the corona in the case of sub-
Eddinton accretion (Haardt & Maraschi, 1991; Liu et al., 2002b).

Figure 1: Coronal temperatures as functions of radius for different accretion
rates. The black hole mass is fixed as $M_{\rm BH}=10M_{\odot}$. The breaks in
the plot of $\dot{M}=10\dot{M}_{\rm Edd}$ are indicated with arrows: the
breaks (1), (2), and (3) correspond to the radii where the loop length
$\ell_{\rm loop}=10r_{\rm S}$ is equal to $\ell_{\rm cor,lim}$, where the loop
length $\ell_{\rm loop}=10r_{\rm S}$ is equal to $r$ (i.e., $r=10r_{\rm S}$),
and where the gas pressure $\dot{M}\Omega/(2\pi r\alpha)$ is equal to the
radiation pressure $aT_{\rm disk}^{4}$ in the disk, respectively. Figure 2:
Scattering optical depths of the corona as functions of radius for different
accretion rates. The black hole mass is fixed as $M_{\rm BH}=10M_{\odot}$. The
explanations of the breaks (1) and (2) in the plot of $\dot{M}=10\dot{M}_{\rm
Edd}$ are the same as those in Fig. 1. Figure 3: Compton-$y$ parameters of the
corona as functions of radius for different accretion rates. The black hole
mass is fixed as $M_{\rm BH}=10M_{\odot}$. Figure 4: The fraction of disk
energy which is dissipated in the corona via magnetic reconnection as
functions of radius for various accretion rates. All the ratios do not exceed
unity.

###  3.2 Radiation Spectra

Figure 5 depicts the results of the Monte Carlo simulations for the radiation
spectra expected from our corona model for different mass accretion rates.
First, we can see that each of the computed spectra consists of a soft thermal
component whose peak frequency is $\sim 10^{17.4}~{}{\rm Hz}$ ($\sim 1~{}{\rm
keV}$), and a hard Comptonized component whose frequency ranges up to $\sim
10^{18.5-19.5}~{}{\rm Hz}$ ($\sim 10-100~{}{\rm keV}$). We can also see that
the peak energy and cutoff energy of the Comptonized component get lower with
higher accretion rate. When the mass accretion rate is $\gtrsim 5\dot{M}_{\rm
Edd}$, especially, the Comptonized component merges with the soft thermal
component, and the spectrum has a broadened peak. Such a feature is consistent
with what is seen in the so-called broadened disk state (Kaaret et al., 2017).

Figure 5: Radiation spectra from super-Eddington accretion flows with wind-fed
coronae for different mass accretion rates. The parameters are $(M_{\rm
BH},\alpha,\eta_{B})=(10M_{\odot},0.1,0.015)$

##  4 Discussion

###  4.1 Origin of an optically-thick corona

According to the spectral analyses for BH objects suspected to harbor
supercritical accretion flows done in the past studies (Gladstone et al.,
2009; Vierdayanti et al., 2010a; Shidatsu et al., 2017; Idogaki et al., 2018),
their X-ray spectra can be fitted by a thermal component and a hard X-ray
component that are supposed to come from an optically-thick accretion disk and
a hot corona surrounding the disk, respectively. However, the corona assumed
in these studies is cool ($k_{B}T\lesssim 10~{}{\rm keV}$) and optically thick
($\tau_{e}\gtrsim 3$) compared to the corona assumed in fitting the spectra of
BH objects with subcritical accretion flows.

Although the formation process of corona is still unclear even in the case
with subcritical accretion flows, it is often considered that the corona is
build up by the evaporation of the underlying disk due to the vertical heat
conduction from a hot corona to a cool disk (Meyer & Meyer-Hofmeister, 1994;
Meyer et al., 2000; Różańska & Czerny, 2000; Liu et al., 2002a, b; Mayer &
Pringle, 2007). In the model proposed by Liu et al. (2002b), they discuss the
structure of a corona that is heated by the reconnection of magnetic fields
emerged from the underlying disk. By analogy with the solar corona (Yokoyama &
Shibata, 2001), they consider the chromospheric evaporation at the interface
of a cool disk and a corona. When the conduction flux from a corona heats up
some of the chromospheric plasma into the magnetic tube, the density of a
corona is determined by the energy balance at the interface, described as the
equation of

$\displaystyle\frac{k_{0}T_{c}^{7/2}}{l}\approx\frac{\gamma}{\gamma-1}nk_{B}T_{c}\left(\frac{k_{B}T_{c}}{m_{H}}\right)^{1/2},$
(18)

where $k_{0}\approx 10^{-6}~{}{\rm erg}~{}{\rm cm}^{-1}~{}{\rm s}^{-1}~{}{\rm
K}^{-7/2}$, $\gamma=5/3$, and $l$ is the length of the magnetic tube.
Combining this equation with the energy balance in the corona (i.e., $Q^{+}$
in Eq.12 is equal to $Q^{-}$ in Eq.13), one can derive the density and
temperature of a corona as functions of radius, black hole mass, and mass
accretion rate. Especially, the scattering optical depth of a corona derived
in their model is limited below unity (see their Fig. 1) on the contrary to
our outflowing corona model. This can be explained in the following way.

In the evaporation scenario, the plasma supply rate to the corona would be
higher when the coronal temperature is higher. However, when the corona
becomes too dense and too optically thick with respect to electron scattering,
due to the inverse Compton cooling, its temperature cannot be too high, which
would regulate the evaporation rate at the interface between the corona and
disk. In the outflowing corona model, by contrast, the coronal plasma is
supplied from the disk wind, whose efficiency is mainly determined by the
radiation pressure inside the disk and the effect of the thermal conduction is
subdominant. Therefore, the scattering optical depth of the corona in our
model is not limited to below unity as in the evaporation scenario.

Since the escape velocity in the innermost region is up to $\sim$ a few ten %
of the speed of light, the outflowing corona there is also assumed to have
relativistic speed. In such a case, the scattering optical depth of the
outflowing corona $\tau_{\rm cor}$ should be modified by a factor of
$\gamma(1-\beta\cos\theta)$, where $\beta=v_{\rm esc}/c$,
$\gamma=(1-\beta^{2})^{-1/2}$, and $\theta$ is the angle between the velocity
of the outflowing corona and the line of sight (Abramowicz et al., 1991).
However, in the realistic case, the outflow is accelerated gradually and
reaches the relativistic speed only when it goes up to reach the height,
$z\sim r$, which is on the order of the scale height of the corona. Therefore,
the modification to the coronal optical depth would be rather subtle, and so
our main conclusions would not be altered.

To fully check the energy balance one should evaluate the radiation flux and
kinetic energy flux of the outflowing corona originating from the locally
liberated energy in the underlying disk per unit surface area. However, since
photons generated inside the accretion flow would be significantly trapped and
advected inward in the case with super-Eddington accretion, it is not
straightforward to check if the results are consistent energetically only from
the distributions of radiation flux. In our calculations, moreover, the
radiation flux from the underlying disk is dominated by reprocesses photons,
whose energy originally comes from the magnetic energy dissipated to heat the
corona. As for the outflow, the amount of ejected material is much smaller
than that of accreted material because we assume $\dot{M}\propto r^{s}$ where
$s=0.15$, which is inferred from the latest radiation hydrodynamical
simulations(Kitaki et al., 2018). Therefore, the kinetic energy flux of the
outflow is negligible compared to the energy liberated in the disk.

###  4.2 Radiation spectra

In each radiation spectrum shown in Figure 5, the hard component (i.e.,
Compton upscattered component) is less energetic than the thermal component.
Such spectra are seen in some ULXs (e.g., NGC 5204 X-1, NGC 4559 X-1, NGC 5408
X-1 etc.; Gladstone et al. (2009)). However, some ULXs (e.g., IC 342 X-1;
Shidatsu et al. (2017)) and other BHs (e.g., GRS 1915+105; Vierdayanti et al.
(2010b)) with supercritical accretion flows show radiation spectra that are
dominated by a hard component, and such spectra cannot be reproduced in our
model as it is, because the corona is not so hot in the inner region, where
most of seed photons are emerged from the disk. However, if we take into
account the radial motion of the outflowing corona, we would have different
results. Actually, the outflowing corona would suffer from the centrifugal
force and be expected to move radially outward. In such a case, the whole
corona would be more optically thick and as a result, the corona would be
cooled down to lower temperature via efficient inverse Compton scattering.

Another effect that may produce a hard component in the spectrum is the
formation of the funnel wall. The radiation-hydrodynamic simulations of
supercritical accretion flows done in Kitaki et al. (2017) show that the
collision between the inflowing material and the funnel wall (i.e., the low-
density region near the polar axis) produces overheated regions whose
temperature is up to $\sim 10^{8}~{}{\rm K}$ in the vicinity of a black hole,
and that significant amount of hard photons are produced there. This is beyond
the scope of the simple analysis made in this paper.

##  5 Summary

We construct a simple model of an outflowing corona formed in a black hole
accretion flow whose mass accretion rate is comparable to or above the
Eddington value. In this model, the coronal plasma is fed by the radiation
pressure-driven wind from an underlying disk, and heated by the reconnection
of magnetic loops emerged from the disk. By solving the energy equation of
this wind-fed corona taking into account the cooling due to the Compton
scattering of soft photons emitted from a disk, we found that the coronal
temperature is relatively low ($k_{B}T_{c}\sim 10~{}{\rm keV}$), and that the
coronal optical depth is relatively large ($\tau_{c}\gtrsim 3$) compared to
the corona expected in subcritical accretion flows, which is independent of
the black hole mass. The existence of such a cool and optically-thick corona
has been inferred from the spectral fitting of GRS 1915+105, some ULXs, and
some NLS1s. We calculate the radiation spectra emerged from this corona by
Monte Carlo simulations, and show that the spectral features of these sources
are qualitatively reproduced in our corona model. Since such an optically-
thick corona is not likely to be formed in a subcritical accretion flow, where
the coronal plasma is fed by the chromospheric evaporation driven by the heat
conduction from a hot corona to a cool disk, because when the corona becomes
dense and optically-thick, it would be efficiently cooled via inverse Compton
scattering, which would suppress the chromospheric evaporation. Since the
corona in our model is fed by the disk wind that is driven by the radiation
pressure in the disk, it can avoid such a regulation process and can be
sufficiently optically-thick. Our results support that the cool and optically-
thick coronae whose existence is inferred from the X-ray spectra observed in
some ULXs, microquasars, and NLS1 are the key signature of supercritical
accretion flows.

We appreciate the anonymous referee for his/her helpful comments and
suggestions. This work is supported in part by JSPS KAKENHI Grant Number,
20K04026 (S. M.). N.K. is supported by the Hakubi project at Kyoto University.

## References

  * Abramowicz et al. (1988) Abramowicz, M. A., Czerny, B., Lasota, J. P. & Szuszkiewicz, E. 1988, ApJ, 332, 646
  * Abramowicz et al. (1991) Abramowicz, M. A., Novikov, I. D., & Paczynski, B. 1991, ApJ, 369, 175. doi:10.1086/169748
  * Bachetti et al. (2014) Bachetti, M., Harrison, F. A., Walton, D. J., et al. 2014, Nature, 514, 202. doi:10.1038/nature13791
  * Bisnovatyi-Kogan & Blinnikov (1976) Bisnovatyi-Kogan, G. S. & Blinnikov, S. I. 1976, Soviet Astronomy Letters, 2, 191
  * Bisnovatyi-Kogan & Blinnikov (1977) Bisnovatyi-Kogan, G. S. & Blinnikov, S. I. 1977, A&A, 59, 111
  * Carpano et al. (2018) Carpano, S., Haberl, F., Maitra, C., et al. 2018, MNRAS, 476, L45. doi:10.1093/mnrasl/sly030
  * Chandra et al. (2020) Chandra, A. D., Roy, J., Agrawal, P. C., et al. 2020, MNRAS, 495, 2664. doi:10.1093/mnras/staa1041
  * Colbert & Mushotzky (1999) Colbert, E. J. M. & Mushotzky, R. F. 1999, ApJ, 519, 89
  * Collin & Kawaguchi (2004) Collin, S. & Kawaguchi, T. 2004, A& A, 426, 797
  * Cropper et al. (2004) Cropper, M., Soria, R., Mushotzky, R. F., Wu, K., Markwardt, C. B. & Pakull, M. 2004, MNRAS, 349, 39
  * Done et al. (2004) Done, C., Wardziński, G. & Gierliński, M. 2004, MNRAS, 349, 393
  * Done et al. (2007) Done, C., Gierlinski, M. & Kubota, A. 2007, A&ARv 15, 1
  * Ebisawa et al. (2003) Ebisawa, K., Zycki, P., Kubota, A., Mizuno, T. & Watarai, K. 2003, ApJ, 597, 780
  * Eggum et al. (1988) Eggum, G. E., Coroniti, F. V. & Katz, J. I. 1988, ApJ, 330, 142
  * Fürst et al. (2016) Fürst, F., Walton, D. J., Harrison, F. A., et al. 2016, ApJ, 831, L14. doi:10.3847/2041-8205/831/2/L14
  * Gladstone et al. (2009) Gladstone, J., Roberts, T. P. & Done, C. 2009, MNRAS, 397, 1836
  * Greiner et al. (2001) Greiner, J., Cuby, J. G. & McCaughrean, M. J. 2001, Nature, 414, 522
  * Haardt & Maraschi (1991) Haardt, F. & Maraschi, L. 1991, ApJ, 380, L51
  * Honma et al. (1991) Honma, F., Kato, S. & Matsumoto, R. 1991, PASJ, 43, 147
  * Idogaki et al. (2018) Idogaki, H. H., Ohmura, S., Takahashi, H., et al. 2018, arXiv e-prints, arXiv:1805.11314
  * Israel et al. (2017) Israel, G. L., Belfiore, A., Stella, L., et al. 2017, Science, 355, 817. doi:10.1126/science.aai8635
  * Israel et al. (2017) Israel, G. L., Papitto, A., Esposito, P., et al. 2017, MNRAS, 466, L48. doi:10.1093/mnrasl/slw218
  * Jin et al. (2017) Jin, C., Done, C., Ward, M., et al. 2017, MNRAS, 471, 706
  * Kaaret et al. (2017) Kaaret, P., Feng, H., & Roberts, T. P. 2017, ARA&A, 55, 303
  * Kato et al. (2008) Kato, S., Fukue, J. & Mineshige, S. 2008, Black-Hole Accretion Disks – Towards a New Paradigm (Kyoto: Kyoto Univ. Press), Chap. 10
  * Kawanaka et al. (2013) Kawanaka, N., Piran, T., & Krolik, J. H. 2013, ApJ, 766, 31
  * Kawashima et al. (2009) Kawashima, T., Ohsuga, K., Mineshige, S., Heizeller, D., Takabe, H. & Matsumoto, R. 2009, PASJ, 61, 769
  * Kawashima et al. (2012) Kawashima, T., Ohsuga, K., Mineshige, S., Yoshida, T., Heinzeller, D. & Matsumoto, R. 2012, ApJ, 752, 18
  * King et al. (2001) King, A. R., Davies, M. B., Ward, M. J., Fabbiano, G. & Elvis, M. 2001, ApJ, 552, L109
  * Kitaki et al. (2017) Kitaki, T., Mineshige, S., Ohsuga, K. & Kawashima, T. 2017, PASJ, 69, 92
  * Kitaki et al. (2018) Kitaki, T., Mineshige, S., Ohsuga, K., et al. 2018, PASJ, 70, 108
  * Liu et al. (2002a) Liu, B. F., Mineshige, S., Meyer, F., et al. 2002a, ApJ, 575, 117
  * Liu et al. (2002b) Liu, B. F., Mineshige, S. & Shibata, K. 2002b, ApJ, 572, L173
  * Liu et al. (2003) Liu, B. F., Mineshige, S. & Ohsuga, K. 2003, ApJ, 587, 571
  * Liu et al. (2013) Liu, J.-F., Bregman, J. N., Bai, Y., et al. 2013, Nature, 503, 500
  * Makishima et al. (2000) Makishima, K. et al. 2000, ApJ, 535, 632
  * Mayer & Pringle (2007) Mayer, M., & Pringle, J. E. 2007, MNRAS, 376, 435
  * Meyer et al. (2000) Meyer, F., Liu, B. F., & Meyer-Hofmeister, E. 2000, A&A, 361, 175
  * Meyer & Meyer-Hofmeister (1994) Meyer, F. & Meyer-Hofmeister, E. 1994, A&A, 288, 175
  * Miller et al. (2003) Miller, J. M., Fabbiano, G., Miller, M. C. & Fabian, A. C. 2003, ApJ, 585, L37
  * Miller et al. (2004) Miller, J. M., Fabian, A. C. & Miller, M. C. 2004, ApJ, 614, L117
  * Mineshige et al. (2000) Mineshige, S., Kawaguchi, T., Takeuchi, M. & Hayashida, K. 2000, PASJ, 52, 499
  * Ohsuga et al. (2005) Ohsuga, K., Mori, M., Nakamoto, T. & Mineshige, S. 2005, ApJ, 628, 368
  * Ohsuga (2007) Ohsuga, K. 2007, ApJ, 659, 205
  * Ohsuga et al. (2009) Ohsuga, K., Mineshige, S., Mori, M. & Kato, Y. 2009, PASJ, 61, L7
  * Ohsuga & Mineshige (2011) Ohsuga, K. & Mineshige, S. 2011, ApJ, 736, 2
  * Okajima et al. (2006) Okajima, T., Ebisawa, K. & Kawaguchi, T. 2006, ApJ, 652, L105
  * Okuda (2002) Okuda, T. 2002, PASJ, 54, 253
  * Okuda et al. (2005) Okuda, T., Teresi, V., Toscano, E. & Molteni, D. 2005, MNRAS, 357, 295
  * Poutanen et al. (2007) Poutanen, J., Lipunova, G., Fabrika, S., Butkevich, A. G. & Abolmasov, P. 2007, MNRAS, 377, 1187
  * Pozdnyakov et al. (1977) Pozdnyakov, L. A., Sobol, I. M., & Sunyaev, R. A. 1977, Soviet Astron., 21, 708
  * Rodríguez Castillo et al. (2020) Rodríguez Castillo, G. A., Israel, G. L., Belfiore, A., et al. 2020, ApJ, 895, 60. doi:10.3847/1538-4357/ab8a44
  * Różańska & Czerny (2000) Różańska, A., & Czerny, B. 2000, MNRAS, 316, 473
  * Rybicki & Lightman (1979) Rybicki, G. B. & Lightman, A. P. 1979, A Wiley-Interscience Publication, New York: Wiley, 1979
  * Sathyaprakash et al. (2019) Sathyaprakash, R., Roberts, T. P., Walton, D. J., et al. 2019, MNRAS, 488, L35. doi:10.1093/mnrasl/slz086
  * Shakura & Sunyaev (1973) Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337
  * Shidatsu et al. (2017) Shidatsu, M., Ueda, Y., & Fabrika, S. 2017, ApJ, 839, 46
  * Sunyaev & Truemper (1979) Sunyaev, R. A. & Truemper, J., 1979, Nature, 279, 506
  * Sutton et al. (2013) Sutton, A. D., Roberts, T. P., & Middleton, M. J. 2013, MNRAS, 435, 1758
  * Tsunoda et al. (2006) Tsunoda, N., Kubota, A., Namiki, M., Sugiho, M., Kawabata, K. & Makishima, K. 2006, PASJ, 58, 1081
  * Vierdayanti et al. (2010a) Vierdayanti, K., Done, C., Roberts, T. P. & Mineshige, S. 2010a, MNRAS, 403, 1206
  * Vierdayanti et al. (2006) Vierdayanti, K., Mineshige, S., Ebisawa, K. & Kawaguchi, T. 2006, PASJ, 58, 915
  * Vierdayanti et al. (2010b) Vierdayanti, K., Mineshige, S. & Ueda, Y. 2010b, PASJ, 62, 239
  * Vierdayanti et al. (2008) Vierdayanti, K., Watarai, K. & Mineshige, S. 2008, PASJ, 60, 653
  * Watarai et al. (2001) Watarai, K., Mizuno, T. & Mineshige, S. 2001, ApJ, 549, L77
  * Wilson-Hodge et al. (2018) Wilson-Hodge, C. A., Malacaria, C., Jenke, P. A., et al. 2018, ApJ, 863, 9. doi:10.3847/1538-4357/aace60
  * Yokoyama & Shibata (2001) Yokoyama, T., & Shibata, K. 2001, ApJ, 549, 1160
  * Yoshida et al. (2013) Yoshida, T., Isobe, N., Mineshige, S., et al. 2013, PASJ, 65, 48
  * Zhou & Wang (1999)


# Satellite Constellation Avoidance with the Rubin Observatory Legacy Survey
of Space and Time

Jinghan Alina Hu Harvey Mudd College, Claremont, CA, USA  Meredith L. Rawls
Department of Astronomy / DiRAC / Vera C. Rubin Observatory, University of
Washington, Seattle, WA, USA  Peter Yoachim Department of Astronomy / DiRAC /
Vera C. Rubin Observatory, University of Washington, Seattle, WA, USA  Željko
Ivezić Department of Astronomy / DiRAC / Vera C. Rubin Observatory, University
of Washington, Seattle, WA, USA

###### Abstract

We investigate a novel satellite avoidance strategy to mitigate the impact of
large commercial satellite constellations in low-Earth orbit on the Vera C.
Rubin Observatory Legacy Survey of Space and Time (LSST). We simulate the
orbits of currently planned Starlink and OneWeb constellations ($\sim$40,000
satellites) to test how effectively an upgraded Rubin scheduler algorithm can
avoid them, and assess how the overall survey is affected. Given a reasonably
accurate satellite orbit forecast, we find it is possible to adjust the
scheduler algorithm to effectively avoid some satellites. Overall, sacrificing
10% of LSST observing time to avoid satellites reduces the fraction of LSST
visits with streaks by a factor of two. Whether such a mitigation will be
required depends on the overall impact of streaks on science, which is not yet
well quantified. This is due to a lack of adequate information about satellite
brightness distributions as well as the impact of glints and low surface
brightness residuals on alert purity and systematic errors in cosmological
parameter estimation. A significant increase in the number of satellites or
their brightness during Rubin Operations may make implementing this satellite
avoidance strategy worthwhile.

Ground-based astronomy, Light pollution, Sky surveys, Artificial Satellites

††software: Astropy (Astropy Collaboration et al., 2013, 2018, 2022), Healpy
and HEALpix666http://healpix.sourceforge.net (Górski et al., 2005; Zonca et
al., 2019), Matplotlib (Hunter, 2007), Numpy (Harris et al., 2020), rubin_sim
(Yoachim et al., 2022b), Scipy (Virtanen et al., 2020), Skyfield (Rhodes,
2019), Shapely (Gillies et al., 2007–)

##  1 Introduction

Rubin Observatory’s Legacy Survey of Space and Time (LSST) is a ten-year
astronomical imaging survey that will begin in 2024 from a new telescope under
construction in Chile. Instead of soliciting individual requests for what the
telescope should observe, the LSST will uniformly survey the sky every few
nights using six color filters to create a decade-long high-resolution survey
of the entire southern sky, and share massive quantities of data products with
the astronomy community (Ivezić et al., 2019). To accomplish this, the LSST
will employ a scheduling algorithm that uses a modified Markov Decision
Process which can generate lists of desirable observations in real time
(Naghib et al., 2019). The LSST scheduler balances the desire to minimize slew
time, optimize signal to noise in individual images, and to maintain survey
footprint uniformity.

One challenge for the LSST is that increasing numbers of bright low-Earth
orbit (LEO) satellites (e.g., Starlink) are being launched, which may leave
streaks in astronomical pointings. LEO satellites are visible from Earth
because they reflect sunlight, especially during twilight. As the Sun-
illuminated satellites move across the field of view of an astronomical
pointing, they leave a streak in the image. While the flux from satellite
streaks can in many cases be identified and removed, the resulting pixels have
much lower signal-to-noise. For a thorough discussion of the scientific
utility of residual light after masking satellite trails, see Hasan et al.
(2022). Over the last three years, many astronomers have raised concerns about
the impact of the proliferation of commercial satellites on the LEO ecosystem
and astronomical surveys (Lawrence et al., 2022; Tyson et al., 2020). In
addition, astronomers have come together with satellite operators and other
stakeholders to create recommendations and strategies to mitigate impacts to
observational astronomy and beyond (Walker et al., 2020a; Hall et al., 2021;
Walker et al., 2020b, 2021).

Tyson et al. (2020) used a very simple algorithm to see if the LSST could
avoid imaging satellite streaks. They concluded that attempting to naively
dodge of order 48,000 LEO satellites is useless because it is operationally
inefficient. In this paper, rather than try to avoid all satellite streaks, we
incorporate satellite avoidance as a component of the LSST scheduler’s Markov
Decision Process. This allows us to avoid a significant fraction of satellite
streaks and investigate what level of avoidance might be acceptable because it
does not drastically impact the overall performance of the LSST.

There are other efforts underway to mitigate the impact of satellite streaks
in astronomical images. For example, satellite companies like SpaceX have
worked on darkening the exterior of satellites so they will be less
visible111https://api.starlink.com/public-
files/BrightnessMitigationBestPracticesSatelliteOperators.pdf. However, even
with the most effective darkening mitigations to date, satellites still appear
bright to the LSST Camera, and are likely to cause effects like non-linear
crosstalk or glints that are challenging to correct with the LSST Science
Pipelines software and may introduce systematic biases or spurious detections.
This is discussed in more detail in Tyson et al. (2020) and on the Rubin
Observatory LSST Project website222https://ls.st/satcon. Astronomers have also
developed algorithms for masking satellite trails in images, but covering the
outer wings of the trails without losing extra pixels remains a challenge
(Hasan et al., 2022). The rapid increase in population of LEO satellites
threatens to compromise the quality and scientific value of LSST images and
also requires extra human and computer resources to effectively mask trails.
Thus, we explore an additional option: incorporating the orbits of known
commercial satellites into the LSST scheduler so the worst of them may be
avoided.

In this paper, we first create realistic simulated forecasts of satellite
orbits in Section 2. We then build a tool that uses that data to create new
scheduler constraints, and test the impact of the modified scheduler algorithm
on LSST observing programs in Section 3. Finally, we discuss the resulting
trade-off between number of streaks and reduced survey depth that results from
avoiding satellites in Section 4, and lay out possibilities for the future as
the satellite population changes during Rubin Operations. We also make
available a GitHub repository with the data and software necessary to
reproduce the paper’s figures333https://github.com/lsst-sims/satellite-
dodging-ApJL.

##  2 Methods

We begin by creating realistic forecasts of three commercial satellite
constellations, which are illustrated in Figure 1. These are Starlink Gen1
(4,408 satellites, altitude $540-570$ km), OneWeb (6,372 satellites, altitude
1200 km), and Starlink Gen2 (29,988 satellites, altitude $340-614$ km). Each
constellation uses orbital inclinations and number of satellite planes
matching current plans444https://ls.st/x1o. To date, OneWeb has launched and
deployed several hundred satellites, while the number of Starlink satellites
is in the thousands.

Figure 1: Three simulated satellite constellations, one per column. Starlink
Gen1 is 4,408 satellites, Starlink Gen2 is 29,988 satellites, and OneWeb is
6,372 satellites, for a grand total of 40,768. The top row shows the 3D
distribution of each constellation around Earth. The middle row shows an
instantaneous Hammer projection of the altitude and azimuth positions of each
constellation as seen from Rubin Observatory on October 1, 2023 during
twilight (Sun altitude $-18$ degrees). Blue points are satellites illuminated
by the Sun at this time, red points are satellites not illuminated by the Sun,
and black points are satellites that are both illuminated and above the Rubin
20 degree altitude pointing limit. The bottom row is the same Hammer
projections six hours later in the middle of the night (Sun altitude $-50$
degrees). Because Starlink satellites orbit at 550 km, none are illuminated in
the middle of the night at this time of year. The OneWeb constellation at 1200
km has only a single illuminated satellite above the Rubin altitude limit at
this particular time.

To simulate LSST observations, we start with the baseline observing strategy
in Yoachim et al. (2022a). LSST observations are scheduled in visits, where a
$u$ visit is one 30s exposure and visits in all other filters ($grizy$) are
back-to-back 15s exposures. The baseline strategy attempts to take most
observations in mixed filter pairs (e.g., an $r$ visit followed by an $i$
visit 33 minutes later), and completes 215,000 visits in the first year.

The baseline LSST observing strategy uses three primary basis functions which
reward (1) minimizing slewtime, (2) maximizing the depth of images (e.g., by
avoiding the Moon and high airmass), and (3) maintaining a uniform survey
footprint. To this, we add a fourth basis function which penalizes observing
areas of the sky which will have high concentrations of illuminated
satellites. Figure 2 shows an example of this new basis function for the three
simulated satellite constellations. The positions of the illuminated
satellites are computed in 10-second intervals and marked on the sky. These
maps are then summed over 90-minute blocks to generate the basis function
maps. Thus our modified scheduler with a satellite avoidance strategy does not
try to avoid individual satellite streaks, but rather has a parameterized
method for avoiding regions of the sky where satellite streaks are more
likely. This has the additional benefit of not requiring high precision
satellite orbit forecasts.

Figure 2: Satellite avoidance maps constructed for the Rubin scheduler for
each simulated constellation. Each is for a twilight observation period of 90
minutes (beginning after sunset with a Sun altitude of $-17.1$ degrees). The
map projections are rotated so zenith is in the center of the image. Darker
regions have more illuminated satellites and therefore more negative
weighting. By varying the dodging weight placed on these maps, the scheduler
will more actively avoid regions of the sky where satellites could streak
images.

We show three example satellite avoidance maps ready for use by the LSST
scheduler in Figure 2. It is apparent from Figure 2 that the simulated OneWeb
constellation has more negative area — regions that should be avoided due to
large numbers of illuminated satellites — than the other two constellations.
Although OneWeb has fewer satellites than Starlink Gen2, the OneWeb satellites
orbit at a higher altitude (1200 km compared to 340-614 km for Starlink),
meaning that they will be illuminated for a longer portion of the night, and
also have a larger impact close to twilight. This is why one of the
recommendations from Walker et al. (2020a) is to keep LEO satellites below 600
km altitude.

To investigate whether the scheduler behaves how we expect with the new
satellite avoidance strategy, we create a testing function that measures the
length of satellite streaks in the simulated fields of view. To ensure
efficiency, only satellites that are above the altitude limit and illuminated
by the Sun are considered. Satellites below the altitude limit (indicated in
the gray region in Figure 2) cannot be reached by the Simonyi Survey Telescope
and are therefore not included. For each satellite, we first determine whether
it is in the field of view for a given pointing by calculating their distance
from the center of the field of view. If this distance is less than the radius
of the field of view, the satellite has crossed through the pointing. To
quantify the impact of the satellite on the pointing, we then project both the
satellite location and the pointing to a 2D x,y plane. In this plane, the
field of view is roughly circular and the start and end locations of the
satellite crossing are two points on the plane, and we can calculate the total
intersection length. Therefore, given a simulated satellite constellation and
a schedule of observations, we are able to record the number of satellites in
each pointing and measure the total streak length. We assume that the impact
of streaks on science is proportional to their total length, which allows us
to quantify the efficiency of the satellite avoidance strategy.

To estimate the fraction of pixels affected by streaks, following Hasan et al.
(2022) we adopt a fiducial width of 300 pixels (equivalent to 1 arcminute
given the plate scale of 0.2 arcseconds per pixel). The 3.5 degree diameter
Rubin focal plane is populated with 189 4kx4k CCDs. Assuming a length of 15
CCDs (with a CCD side of 4096 pixels), a single streak corresponds to 0.6% of
all the pixels in the focal plane. On the other hand, if a streak is so bright
that entire CCDs are rendered scientifically useless, a single streak would
wipe out 8% of all pixels in the focal plane.

We simulate observations for only the first year of the planned ten-year LSST,
as the survey strategy does not significantly change in later years. We
acknowledge this does not account for the likely satellite population increase
beyond the three simulated constellations; however, our results should scale
linearly to larger future constellations in similar orbital distributions. We
do not consider the effects of satellites launching or de-orbiting, and for
simplicity we assume each satellite’s orbital parameters are constant. This
should be an acceptable approximation as long as actual satellite orbital
parameters are available $\sim 1$ day in advance so our avoidance basis
functions can be constructed. If there is no timely information publicly
available on LEO satellite constellation orbits, or it is highly inaccurate,
the satellite avoidance strategy would be impossible to implement. We estimate
that satellite orbital solutions correct to within about a degree in space and
to within a few minutes in time would be sufficient to effectively avoid some
regions of the sky with more satellites in large constellations. Observing
current Starlink satellites, Halferty et al. (2022) find they can construct
TLEs which predict positions with sub-degree spatial precision and sub-second
temporal precision which is more than adequate for our proposed satellite
avoidance methods.

##  3 Results

We find that higher dodging weights reduces pixels lost to satellite streaks,
and that the satellite avoidance algorithm is able to effectively avoid
satellite streaks in simulated pointings. This is shown in the top two panels
of Figure 3. We also find that smaller constellations at lower orbital
altitudes (Starlink Gen1, for example) inherently cause less pixel loss per
pointing, nearly independent of the dodging weight.

Figure 3: Illustration of changes in the mean streak length in all visits
including those with no streak (top left), the fraction of visits with streaks
(top right), the number of acquired visits in year 1 (bottom left) and coadded
depth in the $g$ band (bottom right) as a function of dodging basis function
weight (starting with essentially no dodging on the left). Note that to reduce
the fraction of visits with streaks by about a factor of two, satellite
avoidance will require 10% of total observing time.

Next, we investigate the relationship between the number of exposures the
scheduler is able to complete as a function of the dodging weight. As shown in
the bottom two panels of Figure 3, higher dodging weight results in fewer
visits, most likely due to longer slew times. With a higher dodging weight,
the telescope may be prompted to slew to a location other than the most
desirable nearby pointings, resulting in fewer overall exposures. We also find
that a larger constellation (Starlink Gen2) tends to decrease the number of
exposures slightly more than a smaller constellation (Starlink Gen1), which is
expected. More satellites or satellites at higher orbital altitudes result in
larger areas of avoidance on the sky, which leads to more slewing required to
avoid the affected areas, which subsequently reduces the total number of
visits. In addition to forcing longer slews, avoiding satellite dense areas
pushes the scheduler to observe pointings with lower signal to noise (e.g.,
higher airmass, brighter sky background areas) than it normally would.

One important LSST survey goal is to collect a large number of exposures of
the whole southern sky so these may be co-added to reveal faint structures
that are not visible in individual visits. As a result, survey depth is
crucial to LSST science, and the trade-off between pixel loss from satellite
streaks versus survey depth reduction from fewer total visits must be
evaluated. With the satellite avoidance algorithm, the scheduler is prompted
to avoid regions with illuminated satellites, which sometimes results in
longer slew times or less desirable pointing conditions and contributes to a
loss in survey depth. Therefore, Figure 4 explores the trade-off between
survey depth loss and satellite avoidance. As evident from the figure, the
fraction of LSST visits with streaks can be decreased by a factor of two with
an investment of 10% of LSST observing time, corresponding to a loss of
coadded depth of 0.05 mag.

Figure 4: The trade-off between the mean streak length and final median co-
added depth in $g$ band for one year of the LSST (controlled by dodging basis
function weight). A negative change in co-added depth indicates the survey is
shallower.

So far, we have primarily considered impacts on the overall LSST. However,
another important LSST science goal involves using twilight images to search
for Near Earth Objects (NEOs) and solar system objects interior to the orbits
of Earth and Venus. These observations must be taken in the direction of the
rising or setting sun at high airmass. Because a small potential area is
targeted, our proposed satellite avoidance scheme is ineffective. Figure 5
shows how regular survey observations and twilight NEO observations would be
affected by satellite constellations using the LSST scheduler with no
satellite avoidance. While the majority of twilight NEO observations could
include a satellite streak with a Starlink Gen2 size constellation, we find
this would only result in a $\sim$0.36% loss of science pixels (for a 1
arcminute wide mask).

Figure 5: Impacts of satellite streaks on simulated LSST observations without
any satellite avoidance. Compared to standard LSST observations (top),
twilight NEO observations (bottom) cannot easily be shifted to avoid
satellites. The left panels show the altitude and azimuth distribution of
observations on the sky (zenith at the center of the plots), and the right
panels show how many streaks would result from the three simulated satellite
constellations as a function of how high above the horizon the telescope is
pointing (observation altitude). Note that most twilight NEO observations
would contain a satellite streak.

##  4 Discussion

We have demonstrated that adding a weighted term in the scheduler algorithm
for illuminated satellites can effectively reduce the amount of satellite
streaks in observations, and subsequently reduce mean pixel loss per visit and
other impacts on science, as illustrated in Figure 3. However, with the new
added priority on satellite avoidance, the telescope can be pushed to take an
observation path that does not optimize slew time, which subsequently reduces
the number of exposures and overall survey depth. The trade-off essentially
comes down to the relationship between total streak length reduction and
survey depth reduction. The final decision will depend on the overall impact
of streaks on science which is not well quantified yet due to a lack of
adequate information about the satellite brightness distribution and the
impact of glints and low surface brightness residuals on alert purity and
systematic errors in cosmological parameter estimation. In other words, when
evaluating whether to implement a weighted satellite avoidance strategy to
effectively reduce satellite streak density, it is necessary to evaluate the
trade-off between pixel loss together with non-linear crosstalk, time-domain
glint effects, and any other relevant systematics versus loss in observing
time.

We note that earlier publications (Lawrence et al., 2022; Tyson et al., 2020)
stated that the majority of LSST images are likely to contain a satellite
streak. They also included some higher-altitude Starlink orbits that are no
longer planned. Our study finds that about 10% of all LSST images will have a
streak from the three simulated constellations (Starlink Gen1, Starlink Gen2,
and OneWeb, totaling 40,768 satellites as currently planned). It is true,
however, that twilight observing campaigns at high airmass — like those
necessary to perform NEO searches — will have streaks in the majority of
images.

There is a concerning possibility of sharp increases in satellite population
in the next $5–10$ years, overlapping the LSST operations period
($2024–2035$). With a dramatic increase in satellite population, the ability
to avoid satellites might become more relevant. It is possible to linearly
extrapolate our results to consider a possible LEO satellite population in the
hundreds of thousands circa 2030, assuming the orbital distribution in LEO is
similar to that of Starlink and OneWeb. A future with 400,000 LEO satellites
rather than 40,000 — the stated goal of various companies intending to launch
very large constellations given the present
filings555https://planet4589.org/space/stats/conlist.html — could render the
trade-off of 10% of LSST observing time in order to cut the number of visits
with streaks in half worthwhile. Tyson et al. (2020) find that satellites with
AB magnitudes of $g\sim 3.2$ to $y\sim 1.5$ would saturate LSST images.
Satellite streaks from Starlink and OneWeb as presently designed are not
expected to saturate the LSST Camera’s CCD detectors as they have magnitudes
of $\sim 5$ in their final orbits (Halferty et al., 2022). While a star of
$g=10$ would be very saturated in a LSST images, LEO satellites are moving
fast enough that their effective exposure time is much lower than astronomical
targets. Satellites typically only leave streaks in images when they are both
illuminated by the Sun and visible from the observatory, and LEO satellites
spend most of the night in Earth’s shadow. Satellites from other operators may
be significantly brighter than present-day Starlink and OneWeb satellites, and
may saturate the LSST Camera’s detectors or cause overwhelming levels of non-
linear crosstalk. In particular, the Blue Walker 3 satellite is the first of a
proposed 100 satellite constellation which have a $V$ magnitude between 1-0.
Such a bright satellite would saturate LSST images, potentially causing much
higher pixel losses than satellites which have been launched to date.

Therefore, one future work direction involves adding brightness weighting to
the satellite avoidance algorithm. The idea is to only avoid satellites
brighter than a certain brightness threshold. This could potentially reduce
the region of avoidance, therefore reducing the loss in observing time and co-
added depth. It may be possible to compute optimal starting locations for a
series of observations based on satellite forecasts to further optimize
satellite avoidance. Finally, since faint trail detection and masking is not
perfect, no satellite avoidance strategy will effectively mitigate faint
glints and the resulting bogus alerts.

We thank the anonymous referee who provided helpful prompt feedback that
improved the clarity of the paper. JH acknowledges support from the Computing
Research Association-Widening Participation (CRA-WP) Distributed Research
Experiences for Undergraduates (DREU) program. JH and MLR are grateful for
LSST Corporation travel support for JH to attend the 2022 Rubin Project and
Community Workshop and present a poster. MLR acknowledges support from Bob
Blum, Leanne Guy, and all of Rubin Operations to spend a fraction of her time
on satellite mitigation work; this study would not have been possible without
formal recognition that the proliferation of bright commercial LEO satellites
poses a threat to LSST science. The authors all wish to thank Tony Tyson for
valuable discussions and feedback that helped place this work in context. This
work was facilitated through the use of advanced computational, storage, and
networking infrastructure provided by the Hyak supercomputer system at the
University of Washington.

## References

  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * 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
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Gillies et al. (2007–) Gillies, S., et al. 2007–, Shapely: manipulation and analysis of geometric objects.  https://github.com/Toblerity/Shapely
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759, doi: 10.1086/427976
  * Halferty et al. (2022) Halferty, G., Reddy, V., Campbell, T., Battle, A., & Furfaro, R. 2022, MNRAS, 516, 1502, doi: 10.1093/mnras/stac2080
  * Hall et al. (2021) Hall, J., Walker, C., Rawls, M., et al. 2021, in Bulletin of the American Astronomical Society, 2.205, doi: 10.3847/25c2cfeb.4554c01f
  * 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
  * Hasan et al. (2022) Hasan, I., Tyson, J. A., Saunders, C., & Xin, B. 2022, Astronomy and Computing, 39, 100584, doi: 10.1016/j.ascom.2022.100584
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Lawrence et al. (2022) Lawrence, A., Rawls, M. L., Jah, M., et al. 2022, Nature Astronomy, 6, 428, doi: 10.1038/s41550-022-01655-6
  * Naghib et al. (2019) Naghib, E., Yoachim, P., Vanderbei, R. J., Connolly, A. J., & Jones, R. L. 2019, AJ, 157, 151, doi: 10.3847/1538-3881/aafece
  * Rhodes (2019) Rhodes, B. 2019, Skyfield: High precision research-grade positions for planets and Earth satellites generator, Astrophysics Source Code Library, record ascl:1907.024.  http://ascl.net/1907.024
  * Tyson et al. (2020) Tyson, J. A., Ivezić, Ž., Bradshaw, A., et al. 2020, AJ, 160, 226, doi: 10.3847/1538-3881/abba3e
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Walker et al. (2020a) Walker, C., Hall, J., Allen, L., et al. 2020a, in Bulletin of the American Astronomical Society, Vol. 52, 0206, doi: 10.3847/25c2cfeb.346793b8
  * Walker et al. (2020b) Walker, C., Di Pippo, S., Aubé, M., et al. 2020b, Dark & Quiet Skies I (2020), Dark & Quiet Skies I (2020), Report of the conference held 5-9 October, 2020., doi: 10.5281/zenodo.5898785
  * Walker et al. (2021) —. 2021, Dark & Quiet Skies II (2021), Dark & Quiet Skies II (2021), Report of the conference held 3-7 October, 2021., doi: 10.5281/zenodo.5874725
  * Yoachim et al. (2022a) Yoachim, P., Jones, R. L., Eric H. Neilsen, J., & Ribeiro, T. 2022a, lsst-sims/sims_featureScheduler_runs2.2: Initial release, v1.0.0, Zenodo, doi: 10.5281/zenodo.7150784
  * Yoachim et al. (2022b) Yoachim, P., Jones, L., Eric H. Neilsen, J., et al. 2022b, lsst/rubin_sim: 0.12.1, 0.12.1, Zenodo, doi: 10.5281/zenodo.7087823
  * Zonca et al. (2019) Zonca, A., Singer, L., Lenz, D., et al. 2019, Journal of Open Source Software, 4, 1298, doi: 10.21105/joss.01298


# PPoL: A Periodic Channel Hopping Sequence with Nearly Full Rendezvous
Diversity

Yi-Jheng Lin and Cheng-Shang Chang, , Y.-J. Lin and C.-S. Chang are with the
Institute of Communications Engineering, National Tsing Hua University,
Hsinchu 300, Taiwan, R.O.C. email<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We propose a periodic channel hopping (CH) sequence, called PPoL (Packing the
Pencil of Lines in a finite projective plane), for the multichannel rendezvous
problem. When $N-1$ is a prime power, its period is $N^{2}-N+1$, and the
number of distinct rendezvous channels of PPoL is at least $N-2$ for any
nonzero clock drift. By channel remapping, we construct CH sequences with the
maximum time-to-rendezvous (MTTR) bounded by $N^{2}+3N+3$ if the number of
commonly available channels is at least two. This achieves a roughly 50%
reduction of the state-of-the-art MTTR bound in the literature.

###### Index Terms:

multichannel rendezvous, worst case analysis, finite projective planes.

##  I Introduction

The multichannel rendezvous problem that asks two users to meet each other by
hopping over their available channels is a fundamental problem in many IoT
applications, and it has received a lot of attention lately (see, e.g., the
excellent book [1] and references therein). In the multichannel rendezvous
problem, a channel is called a rendezvous channel of a periodic channel
hopping CH sequence if the two users (with any arbitrary clock drift between
them) rendezvous on that channel within the period of the sequence. The
degree-of-rendezvous (DoR) is the number of distinct rendezvous channels
within the period of the sequence. A periodic CH sequence is said to achieve
full rendezvous diversity (or maximum rendezvous diversity) for the
multichannel rendezvous problem with $N$ channels if its DoR is $N$.

For the multichannel rendezvous problem with $N$ channels, it was shown in
Theorem 1 of [2] that there do not exist deterministic periodic CH sequences
that can achieve full rendezvous diversity with periods less than or equal to
$N^{2}$. The lower bound was further strengthened in Theorem 3 of [3]:

$\displaystyle p\geq\left\\{\begin{array}[]{ll}N^{2}+N&\mbox{if $N\leq 2$}\\\
N^{2}+N+1&\mbox{if $N\geq 3$ and $N$ is a prime power}\\\
N^{2}+2N&\mbox{otherwise}\end{array}\right.,$

where $p$ is the period of a CH sequence with full rendezvous diversity.

One open conjecture in the multichannel rendezvous problem is whether it is
possible to construct a periodic CH sequence that achieves the lower bound. In
the literature, there are various periodic CH sequences that can achieve full
rendezvous diversity, see, e.g., CRSEQ [4], JS [5], DRDS [3], T-CH [6], DSCR
[7], and IDEAL-CH [8]. The asymptotic approximation ratio, defined as the
ratio of the period to the lower bound $N^{2}$ when $N\to\infty$, is still 2
for IDEAL-CH, 2.5 for T-CH and DSCR, and 3 for CRSEQ and DRDS. Clearly, there
is still a significant theoretical gap between the lower bound and the state-
of-the-art CH sequences.

In this letter, we tackle such a conjecture from another perspective. Instead
of focusing on the constructions of periodic CH sequences with full rendezvous
diversity, we consider periodic CH sequences with “nearly” full rendezvous
diversity. To formally define the notion of “nearly” full rendezvous diversity
for a periodic CH sequence $\\{c(t),t=0,1,\ldots,p-1\\}$ on $N$ channels, we
define the quantity $DoR(d)$ as the number of distinct rendezvous channels
between $\\{c(t),t=0,1,\ldots,p-1\\}$ and $\\{c(t+d),t=0,1,\ldots,p-1\\}$.
Clearly, a periodic CH sequence achieves full rendezvous diversity if
$DoR(d)=N$ for all $d\neq 0$. As such, we call a periodic CH sequence achieves
“nearly” full rendezvous diversity if $DoR(d)$ is very close to $N$ for all
$d\neq 0$. For this, we propose a periodic CH sequence, called PPoL, from
Packing the Pencil of Lines in a finite projective plane. When the number of
channels is $N$ with $N-1$ being a prime power, the period $p$ of PPoL is
$N^{2}-N+1$, which is smaller than the lower bound in Theorem 3 of [3] for a
periodic CH sequence with full rendezvous diversity. We show that the $DoR(d)$
of PPoL is at least $N-2$ for all $d$.

Another theoretical contribution of the PPoL CH sequence is the reduction of
the maximum time-to-rendezvous (MTTR). In this letter, we consider the
symmetric, asynchronous, and heterogeneous setting, where (i) the two users
are indistinguishable (and follow the same algorithm to generate their CH
sequences), (ii) their clocks may not be synchronized, and (iii) their
available channels may be different. However, we assume that the $N$ channels
are commonly labeled for these two users. By channel remapping (for the
channels not in the available channel set), we show that the MTTR can be
bounded by $N^{2}+3N+3$ (when $N+1$ is a prime power) if the number of
commonly available channels is at least two. On the other hand, the MTTR
bounds of the state-of-the-art algorithms, including ORTHO-CH [8], SRR [9],
and FRCH [10], are at least $2N^{2}$ when the number of commonly available
channels is at least one. This shows that one can achieve a 50% reduction of
the MTTR bound (for a large $N$) if we relax the assumption of the number of
commonly available channels from one to two.

##  II The PPoL CH sequence

###  II-A Difference sets and finite projective planes

To construct CH sequences with full rendezvous diversity, it is known in [11,
3] that it is equivalent to the problem of packing disjoint difference sets.
Though there are efficient algorithms (see, e.g., [12, 13]) that can find
disjoint perfect difference sets in a periodic sequence, there is no lower
bound on the number of disjoint perfect difference sets that can be found. In
other words, given a specified period $p$, in general we do not have a lower
bound for DoR.

Instead of packing disjoint difference sets, the key idea of our construction
of the PPoL CH sequence is to pack the pencil of lines in a finite projective
plane. For the letter to be self-contained, we first briefly review the
notions of difference sets and finite projective planes (even though they have
been widely used in the literature, see, e.g., [11, 3, 14, 15, 8]).

###### Definition 1

(Difference sets) Let $Z_{p}=\\{0,1,\ldots,p-1\\}$. A set
$D=\\{a_{0},a_{1},\ldots,a_{k-1}\\}\subset Z_{p}$ is called a
$(p,k,\lambda)$-difference set if for every $(\ell\;\mbox{mod}\;p)\neq 0$,
there exist at least $\lambda$ ordered pairs $(a_{i},a_{j})$ such that
$a_{i}-a_{j}=(\ell\;{\rm mod}\;p)$, where $a_{i},a_{j}\in D$. A
$(p,k,1)$-difference set is said to be perfect if there exists exactly one
ordered pair $(a_{i},a_{j})$ such that $a_{i}-a_{j}=(\ell\;{\rm mod}\;p)$ for
every $(\ell\;\mbox{mod}\;p)\neq 0$.

###### Definition 2

(Finite projective planes) A finite projective plane of order $m$ is a
collection of $m^{2}+m+1$ lines and $m^{2}+m+1$ points such that

(i)

every line contains $m+1$ points,

(ii)

every point is on $m+1$ lines,

(iii)

any two distinct lines intersect at exactly one point, and

(iv)

any two distinct points lie on exactly one line.

In the following, we state some well-known facts for the perfect difference
sets and the finite projective planes. Suppose that
$D=\\{a_{0},a_{1},\ldots,a_{m}\\}$ is an $(m^{2}+m+1,m+1,1)$-perfect
difference set. Let $p=m^{2}+m+1$.

(P1): If for some indices $i,j,k,\ell$,

$((a_{i}-a_{j})\;{\rm mod}\;p)=((a_{k}-a_{\ell})\;{\rm mod}\;p)\neq 0,$

then $a_{i}=a_{k}$ and $a_{j}=a_{\ell}$. This is due to the unique difference
representation property of a perfect difference set.

(P2): A time-shifted version of $D$ is still a perfect difference set.
Specifically, let

$D_{\ell}=\\{(a_{0}+\ell)\;\mbox{mod}\;p,(a_{1}+\ell)\;{\rm
mod}\;p,\ldots,(a_{m}+\ell)\;{\rm mod}\;p\\},$ (2)

$\ell=0,1,2,\ldots,p-1$. Then $D_{\ell}$ is an $(m^{2}+m+1,m+1,1)$-perfect
difference set (from the unique difference representation property).

(P3): There exists exactly one common element in $D$ and $D_{\ell}$ for
$\ell\neq 0$ such that $a_{i}=(a_{j}+\ell)\;{\rm mod}\;p$ (from the unique
difference representation property). Also, as there is exactly one ordered
pair $(a_{i},a_{j})$ such that $a_{i}-a_{j}=1$, without loss of generality we
can assume that $a_{0}=0$ and $a_{1}=1$ and order all the elements in $D$ in
the increasing order as follows:

$a_{0}=0<a_{1}=1<a_{2}<\ldots,<a_{m}<p.$ (3)

(P4): Singer [16] established an important connection between an
$(m^{2}+m+1,m+1,1)$-perfect difference set and a finite projective plane of
order $m$.

(i)

Let $\\{0,1,\ldots,m^{2}+m\\}$ be the $m^{2}+m+1$ points.

(ii)

Let $D_{\ell}$, $\ell=0,1,2,\ldots,p-1$ be the $m^{2}+m+1$ lines.

Then these $m^{2}+m+1$ points and $m^{2}+m+1$ lines form a finite projective
plane of order $m$.

(P5): The $m+1$ lines in the corresponding finite projective plane that
contain point $0$ are $D_{0},D_{p-a_{1}},D_{p-a_{2}},\ldots,D_{p-a_{m}}$.
These $m+1$ lines are called the pencil of lines [16] that contain point 0 (as
the pencil point). As the only intersection of the $m+1$ lines is point 0,
these $m+1$ lines, excluding point 0, are disjoint, and thus can be packed
into $Z_{p}$, i.e., $\\{D_{0},D_{p-a_{1}}^{0},\ldots,D_{p-a_{m}}^{0}\\}$ forms
a partition of $Z_{p}$, where

$D_{p-a_{i}}^{0}=D_{p-a_{i}}\backslash\\{0\\},\quad i=1,2,\ldots,m.$ (4)

###  II-B Construction of the CH sequence

In this section, we propose the PPoL algorithm for constructing CH sequences
with nearly full rendezvous diversity. For this, one first constructs an
$(m^{2}+m+1,m+1,1)$-perfect difference set, $D=\\{a_{0},a_{1},\ldots,a_{m}\\}$
with

$a_{0}=0<a_{1}=1<a_{2}<\ldots,<a_{m}<p,$ (5)

where $p=m^{2}+m+1$. This is feasible when $m$ is a prime power [16]. Then one
assigns each channel to a line packed in $Z_{p}$. Specifically, the PPoL CH
sequence $\\{c(t),t=0,1,\ldots,p-1\\}$ is constrcuted by assigning channel 0
to the time slots in $D_{0}$ and channel $i$, $i=1,2,\ldots,m$, to the time
slots in $D_{p-a_{i}}^{0}$, i.e.,

$\displaystyle c(t)=\left\\{\begin{array}[]{ll}0&\mbox{if $t\in D_{0}$}\\\
i&\mbox{if $t\in D_{p-a_{i}}^{0}$ for some $i\neq 0$}\end{array}\right..$ (8)

Algorithm 1 The PPoL CH sequence

Input A set of $m+1$ channels $\\{0,1,2,\ldots,m\\}$ with $m$ being a prime
power.

Output A CH sequence $\\{c(t),t=0,1,\ldots,p-1\\}$, where $p=m^{2}+m+1$.

1: Let $p=m^{2}+m+1$ and construct a perfect difference set
$D=\\{a_{0},a_{1},\ldots,a_{m}\\}$ in $Z_{p}$.

2: For $\ell=0,1,\ldots,p-1$, let

$D_{\ell}=\\{(a_{0}+\ell)\;{\rm mod}\;p,(a_{1}+\ell)\;{\rm
mod}\;p,\ldots,(a_{m}+\ell)\;{\rm mod}\;p\\}.$

3: Let $D_{p-a_{i}}^{0}=D_{p-a_{i}}\backslash\\{0\\}$, $i=1,2,\ldots,m$.

4: Construct the CH sequence $\\{c(t),t=0,1,\ldots,p-1\\}$ by assigning
channel 0 to the time slots in $D_{0}$ and channel $i$, $i=1,2,\ldots,m$, to
the time slots in $D_{p-a_{i}}^{0}$

TABLE I: An illustrating example for the PPoL CH sequence with the perfect difference set $D=D_{0}=\\{0,1,4,6\\}$ in $Z_{13}$. $d=a_{k}-a_{j}$ | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12
---|---|---|---|---|---|---|---|---|---|---|---|---|---
$(k,j)$ |  | (1,0) | (3,2) | (2,1) | (2,0) | (3,1) | (3,0) | (0,3) | (1,3) | (0,2) | (1,2) | (2,3) | (0,1)
$D_{0}$ shifts $d$ | 0,1,4,6 |  1,2,5,7 |  2,3,6,8 |  3,4,7,9 |  4,5,8,10 |  5,6,9,11 |  6,7,10,12 |  7,8,11,0 |  8,9,12,1 |  9,10,0,2 |  10,11,1,3 |  11,12,2,4 |  12,0,3,5
$D_{12}^{0}$ shifts $d$ | 12,3,5 | 0,4,6 |  1,5,7 | 2,6,8 |  3,7,9 | 4,8,10 |  5,9,11 |  6,10,12 | 7,11,0 |  8,12,1 | 9,0,2 |  10,1,3 | 11,2,4
$D_{9}^{0}$ shifts $d$ | 9,10,2 |  10,11,3 | 11,12,4 | 12,0,5 | 0,1,6 |  1,2,7 |  2,3,8 |  3,4,9 |  4,5,10 | 5,6,11 | 6,7,12 | 7,8,0 |  8,9,1
$D_{7}^{0}$ shifts $d$ | 7,8,11 |  8,9,12 | 9,10,0 |  10,11,1 |  11,12,2 | 12,0,3 | 0,1,4 | 1,2,5 | 2,3,6 |  3,4,7 |  4,5,8 | 5,6,9 |  6,7,10

In Table I, we provide an illustrating example for the PPoL CH sequence with
the perfect difference set $D=D_{0}=\\{0,1,4,6\\}$ in $Z_{13}$. The other
three lines that contain point 0 are $D_{13-1}=D_{12}=\\{12,0,3,5\\}$,
$D_{13-4}=D_{9}=\\{9,10,0,2\\}$ and $D_{13-6}=D_{7}=\\{7,8,11,0\\}$. Thus, we
have $D_{12}^{0}=\\{12,3,5\\}$, $D_{9}^{0}=\\{9,10,2\\}$ and
$D_{7}^{0}=\\{7,8,11\\}$. According to Algorithm 1, we assign channel 0 for
$t=0,1,4,6$, channel 1 for $t=3,5,12$, channel 2 for $t=2,9,10$ and channel 3
for $t=7,8,11$. This leads to the CH sequence

$\\{c(t),t=0,1,\ldots,12\\}=\\{0,0,2,1,0,1,0,3,3,2,2,3,1\\}.$

As $D_{0}$ is a perfect difference set, channel 0 is a rendezvous channel. In
particular, for $d=1,2,\ldots,12$, let

$t_{0}(d)=1,6,4,4,6,6,0,1,0,1,4,0$

be the time slots marked in red in the third row of Table I. Then

$c(t_{0}(d))=c(t_{0}(d)+d)=0.$

However, for $d=1,3,5,8,10,12$, channel 1 is not a rendezvous channel.
Similarly, for $d=2,3,4,9,10,11$, channel 2 is not a rendezvous channel. Also,
for $d=2,5,6,7,8,11$, channel 3 is not a rendezvous channel. As shown in Table
I, for any $d$, there are at most two channels that are not rendezvous
channels (see the cells marked in grey).

In the following theorem, we prove the main result of this letter.

###### Theorem 3

For any $d\neq 0$, the $DoR(d)$ of the PPoL CH sequence in Algorithm 1 is at
least $m-1$ for a system with $m+1$ channels, where $m$ is a prime power.

The proof of Theorem 3 requires the following lemma.

###### Lemma 4

If

$d=(a_{k}-a_{j})\;{\rm mod}\;p$

for some $j\neq i$ and $k\neq i$, then there exists $0\leq t\leq p-1$ such
that

$c(t)=c(t+d)=i,$ (9)

i.e., channel $i$ is a rendezvous channel for such a clock drift $d$.

Proof. Let $t=(a_{j}-a_{i})\mod p$. Since $j\neq i$, $t\in D_{p-a_{i}}^{0}$.
Thus, $c(t)=i$ from (8). On the other hand,

$\displaystyle(t+d)\mod p$ $\displaystyle=(a_{j}-a_{i}+(a_{k}-a_{j}))\mod p$
$\displaystyle=(a_{k}-a_{i})\mod p.$

Since $k\neq i$, $((t+d)\mod p)\in D_{p-a_{i}}^{0}$. Thus, $c(t+d)=i$ from
(8).

Proof. (Theorem 3) Since $D$ is a perfect difference set, for any $0<d<p$,
there exists a unique ordered pair $(a_{j},a_{k})$ such that
$d=(a_{j}-a_{k})\;{\rm mod}\;p$. It then follows from Lemma 4 that except
channels $j$ and $k$, every channel is a rendezvous channel when
$d=(a_{j}-a_{k})\;{\rm mod}\;p$. Thus, the $DoR(d)$ of the PPoL CH sequence in
Algorithm 1 is at least $m-1$ for a system with $m+1$ channels for any clock
drift $d$.

##  III The remapped PPoL CH sequences

In this section, we show how one remaps the PPoL CH sequence for the
multichannel rendezvous problem in the symmetric, asynchronous, and
heterogeneous setting. One popular remapping method is known as random
remapping that randomly re-assigns a channel not in the channel available set
to a channel in the channel available set. As the $DoR(d)$ of the PPoL CH
sequence is at least $N-2$ from Theorem 3, the two users following the PPoL CH
sequences with random remapping are guaranteed to rendezvous within
$m^{2}+m+1$ time slots for any prime power $m\geq N-1$ if the number of
commonly available channels between these two users is not smaller than three.

Now we show that one can further reduce the requirement for the minimum number
of commonly available channels from three to two. Specifically, for the
multiple channel rendezvous problem with $N$ channels, we first construct the
PPoL sequence from an $(m^{2}+m+1,m+1,1)$-perfect difference set, where $m$ is
the smallest prime power not smaller than $N+1$ (instead of $N-1$ in Algorithm
1). For a user, let $n$ be the number of its available channels, and
$n^{\prime}=m+1-n$ be the number of channels not in its available channels
among the $m+1$ channels in Algorithm 1. Consider the following two cases:

Case 1. $n>(N+2)/2$:

In this case, we simply use random remapping.

Case 2. $n\leq(N+2)/2$:

Since $m\geq N+1$, we know that

$n^{\prime}=m+1-n\geq\frac{N+2}{2}\geq n.$

Thus, we can remap the first $n$ channels not in the available channels to the
$n$ available channels through a one-to-one deterministic function. For the
rest $n^{\prime}-n$ channels not in the available channels, we simply use
random remapping.

One key insight of the deterministic remapping in Case 2 is that every
available channel of a user is assigned to two lines of the $m+1$ lines in the
PPoL CH sequence. The detailed steps are shown in Algorithm 2.

Algorithm 2 The remapped PPoL CH sequence

Input: A set of available channels ${\bf c}=\\{c_{0},c_{1},\ldots,c_{n-1}\\}$
that is a subset of the $N$ channels $\\{0,1,\ldots,N-1\\}$.

Output: A CH sequence $\\{c(t),t=0,1,\ldots\\}$ with $c(t)\in{\bf c}$.

1: Let $m$ be the smallest prime power such that $m\geq N+1$ and
$p=m^{2}+m+1$.

2: Use Algorithm 1 to generate the PPoL CH sequence
$\\{c(t),t=0,1,\ldots,p-1\\}$.

3: Let ${\bf c}^{c}=Z_{m+1}\backslash{\bf
c}=\\{c_{0}^{\prime},c_{1}^{\prime},\ldots,c_{m+2-n}^{\prime}\\}$ be the set
of channels not in the available channel set.

4: Case 1. $n>(N+2)/2$: Use random remapping, i.e., if $c(t)=c_{j}^{\prime}$
for some channel $c_{j}^{\prime}$ in ${\bf c}^{c}$, remap $c(t)$ randomly to a
channel in ${\bf c}$.

5: Case 2. $n\leq(N+2)/2$: Suppose $c(t)=c_{j}^{\prime}$ for some channel
$c_{j}^{\prime}$ in ${\bf c}^{c}$. If $j<n$, remap $c(t)=c_{j}$. Otherwise
remap $c(t)$ randomly to a channel in ${\bf c}$.

###### Theorem 5

Consider two users with the available channel sets ${\bf c}_{1}$ and ${\bf
c}_{2}$. Suppose that there are at least two commonly available channels
between these two users, i.e., $|{\bf c}_{1}\cap{\bf c}_{2}|\geq 2$. If the
two users use the remapped PPoL CH sequences in Algorithm 2 to generate their
CH sequences, then these two users are guaranteed to rendezvous within
$m^{2}+m+1$ time slots, where $m$ is the smallest prime power not smaller than
$N+1$.

As a direct consequence of Theorem 5, the MTTR of the remapped PPoL CH
sequences in Algorithm 2 is $N^{2}+3N+3$ if $N+1$ is a prime power. Such an
MTTR bound is substantially smaller than those from the state-of-the-art
algorithms, including $(2N+1)N$ in ORTHO-CH [8], $(2N+2)N$ in SRR [9], and
$(2N+1)N$ in FRCH [10]. The reduction of the MTTR bound is due to the
assumption that the number of commonly available channels is at least two.

For the proof of Theorem 5, we need the following extension of Lemma 4. The
result in Lemma 6 shows that if a channel of one user appears in the time
slots of two different lines, one from the original assignment and the other
from the remapping of a channel not in the available channel set, then that
channel is a rendezvous channel except for one half of the set of clock drifts
in Lemma 4.

###### Lemma 6

Consider two sequences $\\{c_{1}(t),0\leq t\leq p-1\\}$ and $\\{c_{2}(t),0\leq
t\leq p-1\\}$. Suppose that $c_{1}(t)=i$ for $t\in D_{p-a_{i}}^{0}$ and that
$c_{2}(t)=i$ for $t\in D_{p-a_{i}}^{0}\cup D_{p-a_{i_{2}}}^{0}$ for some
$i_{2}\neq i$. If

$d=(a_{k}-a_{j})\;{\rm mod}\;p,$ (10)

for some $k\neq i$, there exists $0\leq t\leq p-1$ such that

$c_{1}(t)=c_{2}((t+d)\mod p)=i.$ (11)

Thus, channel $i$ is a rendezvous channel for such a clock drift $d$.

Proof. We consider the following two cases.

Case 1. $k\neq i$ and $j\neq i$:

If $k\neq i$ and $j\neq i$, then by Lemma 4, there exists $0\leq t\leq p-1$
such that both $t$ and $((t+d)\mod p)$ are in $D_{p-a_{i}}^{0}$. Thus
$c_{1}(t)=c_{2}(t+d)=i$ from (9).

Case 2. $k\neq i$ and $j=i$:

First, we argue that $(2a_{i}-a_{i_{2}}-a_{k})\mod p\neq 0$ when $k\neq i$. If
$(2a_{i}-a_{i_{2}}-a_{k})\mod p=0$, then

$(a_{i}-a_{i_{2}})\mod p=(a_{k}-a_{i})\mod p.$ (12)

Since $k\neq i$, $(a_{k}-a_{i})\mod p\neq 0$. From the unique difference
representation property for the perfect difference set $D$, we have from (12)
that $a_{i}=a_{k}$ and $a_{i_{2}}=a_{i}$. This contracts to the fact that
$k\neq i$.

Since $(2a_{i}-a_{i_{2}}-a_{k})\mod p\neq 0$ for $k\neq i$, there exists a
unique pair of $a_{\ell_{1}}$ and $a_{\ell_{2}}$ in $D$ such that

$(a_{\ell_{1}}-a_{\ell_{2}})\mod p=(2a_{i}-a_{i_{2}}-a_{k})\mod p.$ (13)

We argue that $\ell_{1}\neq i$ and $\ell_{2}\neq i_{2}$. If $\ell_{1}=i$, then
it follows from (13) that

$(a_{i_{2}}-a_{\ell_{2}})\mod p=(a_{i}-a_{k})\mod p.$ (14)

From the unique difference representation property for the perfect difference
set $D$, we have from (14) that $i_{2}=i$ and $\ell_{2}=k$. This contracts to
the fact that $i_{2}\neq i$. On the other hand, if $\ell_{2}=i_{2}$, then it
follows from (13) that

$(a_{\ell_{1}}-a_{i})\mod p=(a_{i}-a_{k})\mod p.$ (15)

From the unique difference representation property for the perfect difference
set $D$, we have from (15) that $\ell_{1}=i$ and $i=k$. This contracts to the
fact that $i\neq k$.

Select

$t=(a_{\ell_{1}}-a_{i})\mod p.$ (16)

Since $\ell_{1}\neq i$, we know that $t\in D_{p-a_{i}}^{0}$ and thus

$c_{1}(t)=i.$ (17)

On the other hand, we have from (10) and (13) that

$\displaystyle t+d$ $\displaystyle=$
$\displaystyle(a_{\ell_{1}}-a_{i})+(a_{k}-a_{i})\mod p$ (18) $\displaystyle=$
$\displaystyle(a_{\ell_{2}}-a_{i_{2}})\mod p.$

Since $\ell_{2}\neq i_{2}$, we know that $t+d\in D_{p-a_{i_{2}}}^{0}$ and thus

$c_{2}(t+d)=i.$ (19)

Proof. (Proof of Theorem 5) Let $n_{1}=|{\bf c}_{1}|$ (resp. $n_{2}=|{\bf
c}_{2}|$) be the number of available channels of user 1 (resp. user 2). Also,
let $n_{1,2}=|{\bf c}_{1}\cap{\bf c}_{2}|$ be the number of common channels
between these two users. We consider the following two cases:

Case 1. $n_{1}+n_{2}\geq N+3$:

Since ${\bf c}_{1}$ and ${\bf c}_{2}$ are subsets of the $N$ channels,

$N\geq|{\bf c}_{1}\cup{\bf c}_{2}|=n_{1}+n_{2}-n_{1,2}.$

Since we assume that $n_{1}+n_{2}\geq N+3$ in this case, we have $n_{1,2}\geq
3$. As a result of Theorem 3, these two users are guaranteed to rendezvous
within $p$ time slots (even without the need of remapping).

Case 2. $n_{1}+n_{2}\leq N+2$:

Without loss of generality, we assume that $n_{2}\leq n_{1}$. In this case, we
have that $n_{2}\leq(N+2)/2$. According to Case 2 of Algorithm 2, the first
$n_{2}$ channels in ${\bf c}^{c}_{2}$ are remapped to the $n_{2}$ channels in
${\bf c}_{2}$ through a one-to-one deterministic function. In other words, for
each available channel $i$ in ${\bf c}^{c}_{2}$, we have $c_{2}(t)=i$ for
$t\in D_{p-a_{i}}^{0}\cup D_{p-a_{i_{2}}}^{0}$ for some $i_{2}\neq i$.

Since we assume that $|{\bf c}_{1}\cap{\bf c}_{2}|\geq 2$, there exist two
distinct channels $\alpha$ and $\beta$ in ${\bf c}_{1}\cap{\bf c}_{2}$. For
the CH sequences of these two users $\\{c_{1}(t),0\leq t\leq p-1\\}$ and
$\\{c_{2}(t),0\leq t\leq p-1\\}$, we know that $c_{1}(t)=\alpha$ for $t\in
D_{p-a_{\alpha}}^{0}$ and that $c_{2}(t)=\alpha$ for $t\in
D_{p-a_{\alpha}}^{0}\cup D_{p-a_{\alpha_{2}}}^{0}$ for some
$\alpha_{2}\neq\alpha$. Similarly, we also know that $c_{1}(t)=\beta$ for
$t\in D_{p-a_{\beta}}^{0}$ and that $c_{2}(t)=\beta$ for $t\in
D_{p-a_{\beta}}^{0}\cup D_{p-a_{\beta_{2}}}^{0}$ for some
$\beta_{2}\neq\beta$.

Now we prove by contradiction that at least one of the two commonly available
channels is a rendezvous channel. Suppose that both channels are not
rendezvous channels. Then we have from Lemma 6 that the clock drift $d$ must
satisfy

$d=(a_{\alpha}-a_{j_{1}})\mod p=(a_{\beta}-a_{j_{2}})\mod p$

for some $j_{1}$ and $j_{2}$. From the unique difference representation
property for the perfect difference set $D$, we have $\alpha=\beta$ and
$j_{1}=j_{2}$. This contradicts to the assumption that $\alpha\neq\beta$.

In Theorem 5, we have shown the MTTR bound for the remapped PPoL CH sequence
if the number of commonly available channels between the two users is at least
two. To further reduce the number of commonly available channels to 1, we need
to ensure that each available channel of a user is mapped to two lines. This
is stated in the following corollary.

###### Corollary 7

Consider two users with the available channel sets ${\bf c}_{1}$ and ${\bf
c}_{2}$. Let $m$ be the smallest prime power such that $m\geq N+1$. If the two
users use the remapped PPoL CH sequence in Algorithm 2 to generate their CH
sequences and

$\max[n_{1},n_{2}]\leq\frac{N+2}{2},$ (20)

then these two users are guaranteed to rendezvous on every commonly available
channel within $m^{2}+m+1$ time slots.

Proof. According to Case 2 of Algorithm 2 and the assumption in (20), each
available channel of each user is mapped to two lines, i.e., for each
available channel $i$ in ${\bf c}_{1}\cap{\bf c}_{2}$, we have $c_{1}(t)=i$
for $t\in D_{p-a_{i}}^{0}\cup D_{p-a_{i_{1}}}^{0}$ for some $i_{1}\neq i$, and
$c_{2}(t)=i$ for $t\in D_{p-a_{i}}^{0}\cup D_{p-a_{i_{2}}}^{0}$ for some
$i_{2}\neq i$.

Now we argue that each channel $i$ in ${\bf c}_{1}\cap{\bf c}_{2}$ is a
rendezvous channel. This holds trivially for the clock drift $d=0$. Suppose
that channel $i$ in ${\bf c}_{1}\cap{\bf c}_{2}$ is not rendezvous channels
for some $d\neq 0$. Then we have from Lemma 6 that

$d=(a_{i}-a_{j_{1}})\mod p$ (21)

for some $j_{1}$. By interchanging $c_{1}(t)$ and $c_{2}(t)$ in Lemma 6, we
also know that

$d=(a_{j_{2}}-a_{i})\mod p$ (22)

for some $j_{2}$. From the unique difference representation property for the
perfect difference set $D$, we must have $i=j_{2}$ and $j_{1}=i$. This
contradicts to the assumption that $d\neq 0$.

##  IV Conclusion

In this letter, we proposed in Algorithm 1 the PPoL CH sequence with
$DoR(d)\geq N-2$ for the multichannel rendezvous problem with $N$ channels.
Such a CH sequence guarantees the rendezvous of two users if the number of
commonly available channels is at least three. The remapped PPoL CH sequence
in Algorithm 2 further relaxes the number of commonly available channels from
three to two. These new results have a roughly 50% reduction of the state-of-
the-art MTTR bound in the literature.

Through extensive simulations, we also observed that the expected time-to-
rendezvous (ETTR) of the PPoL CH sequence (with random remapping) is almost
the same as those of the simple random algorithm and several existing
algorithms, including ORTHO-CH [8], CRSEQ [4], DRDS [3], and T-CH [6]. Due to
space limitations, these numerical results for ETTR’s are omitted.

## References

  * [1] Z. Gu, Y. Wang, Q.-S. Hua, and F. C. M. Lau, Rendezvous in Distributed Systems: Theory, Algorithms and Applications. Springer, 2017.
  * [2] K. Bian and J.-M. Park, “Maximizing rendezvous diversity in rendezvous protocols for decentralized cognitive radio networks,” IEEE Transactions on Mobile Computing, vol. 12, no. 7, pp. 1294-1307, 2013.
  * [3] Z. Gu, Q.-S. Hua, Y. Wang, and F. C. M. Lau, “Nearly optimal asynchronous blind rendezvous algorithm for cognitive radio networks,” in Proc. IEEE SECON, 2013.
  * [4] J. Shin, D. Yang, and C. Kim, “A channel rendezvous scheme for cognitive radio networks,” IEEE Communications Letter, vol. 14, no. 10, pp. 954-956, 2010.
  * [5] Z. Lin, H. Liu, X. Chu, and Y.-W. Leung, “Jump-stay based channel-hopping algorithm with guaranteed rendezvous for cognitive radio networks,” in Proc. IEEE INFOCOM 2011.
  * [6] G. Y. Chang, J. F. Huang, and Y. S. Wang, “Matrix-based channel hopping algorithms for cognitive radio networks,” IEEE Transactions on Wireless Communications, vol. 14, no. 5, pp. 2755–2768, 2015.
  * [7] B. Yang, M. Zheng, and W. Liang, “A time-efficient rendezvous algorithm with a full rendezvous degree for heterogeneous cognitive radio networks,” in Proc. IEEE INFOCOM, pp. 1–9, 2016.
  * [8] C.-S. Chang, J.-P. Sheu, and Y.-J.Lin, “On the theoretical gap of channel hopping sequences with maximum rendezvous diversity in the multichannel rendezvous problem, IEEE/ACM Transactions on Networking, early access, 2021.
  * [9] Y. Fu, Y. Wang, Z. Gu, X. Zheng, T. Wei, Z. Cao, H. Cui, and F. C. M. Lau, “How local information improves rendezvous in cognitive radio networks,” In Proc. 15th annual IEEE international conference on sensing, communication, and networking (SECON), pp. 1–9, 2018.
  * [10] G.-Y. Chang and J.-F. Huang, “A fast rendezvous channel-hopping algorithm for cognitive radio networks,” IEEE Communications Letters, vol. 17, no. 7, pp. 1475–1478, 2013.
  * [11] F. Hou, L. X. Cai, X. Shen, and J. Huang, “Asynchronous multichannel MAC design with difference-set-based hopping sequences,” IEEE Transactions on Vehicular Technology, vol. 60, pp. 1728–1739, 2011.
  * [12] K. Wu, F. Han, F. Han, and D. Kong, “Rendezvous sequence construction in cognitive radio ad-hoc networks based on difference sets,” in Proc. IEEE PIMRC, pp. 1840–1845, 2013.
  * [13] X. J. Tan, C. Zhou, and J. Chen, “Symmetric channel hopping for blind rendezvous in cognitive radio networks based on union of disjoint difference sets,” IEEE Transactions on Vehicular Technology, vol. 66, no. 11, pp. 10233-10248, 2017.
  * [14] C.-S. Chang, W. Liao and C.-M. Lien, “On the multichannel rendezvous problem: fundamental limits, optimal hopping sequences, and bounded time-to-rendezvous,” Mathematics of Operations Research, vol. 40, no. 1, pp. 1-23, 2015.
  * [15] C.-S. Chang, D.-S, Lee, and C. Wang, “Asynchronous grant-free uplink transmissions in multichannel wireless networks with heterogeneous QoS guarantees,” IEEE/ACM Transactions on Networking, vol. 27, no. 4, pp. 1584-1597, 2019.
  * [16] J. Singer, “A theorem in finite projective geometry and some applications to number theory,” Transactions of the American Mathematical Society, vol. 43, no. 3, pp. 377–385, 1938.