\section{Introduction}
\label{sec:introduction}

We report the results of an all-sky search for continuous, nearly monochromatic gravitational waves from
rapidly rotating isolated neutron stars using the first six months of data from the third observing run (O3) of the
Advanced Laser Interferometer Gravitational wave Observatory (Advanced LIGO \cite{bib:LIGO,bib:DetectorPaper}). This first search of the
early O3 data uses the {\footnotesize\textsc{PowerFlux}}\ hierarchical search program~\cite{bib:S4allsky,bib:S5allsky1,bib:S5allsky2,bib:S6allsky,bib:O1allsky1,bib:O1allsky2}
with loose-coherence follow-up~\cite{bib:loosecoherence1,bib:loosecoherence2} of outliers and covers a broad band of frequencies
(20-2000 Hz) and frequency derivatives ($\sci{[-1.0, +0.1]}{-8}$~Hz/s).
Although more than $\sci{1.4}{5}$ search outliers are found in the initial stage of analysis, successive follow-up stages failed to
confirm an astrophysical signal. Hence the primary results from this analysis are the
upper limits on strain amplitude presented in section~\ref{sec:upperlimits}.

All-sky searches for continuous gravitational waves from isolated neutron stars
have been carried out in Advanced LIGO and Virgo data
previously~\cite{bib:O1allsky1,bib:O1allsky2,bib:O1EatH,bib:O2allsky,bib:O2allskyRome,bib:O2EatH,bib:O2Falcon1,bib:O2Falcon2,bib:O2Falcon3,bib:deepsearch171-172Hz}.
The results presented here are the most sensitive to date in strain amplitude for a broadband, all-sky search with high allowed
spin-down magnitudes, with improvement factors ranging from $\sim$1.2 at signal frequencies of
$\sim$100 Hz to $\sim$2 at the higher frequencies.
The improvement with respect to similar, previous analyses of O2 data~\cite{bib:O2allsky,bib:O2allskyRome}
is greater at higher frequencies, in part because of the improved detector noise ($\sim$3db) achieved with
quantum squeezing~\cite{bib:O3squeezing} and in part because
longer Fourier transform coherence times are used here (for frequencies up to 1475 Hz) than in the O2 analyses.

Our primary targets in this analysis are fast-spinning, non-axisymmetric neutron stars in the Milky Way.
Given the immense pressure on its nuclear matter, one expects a neutron star to assume a
highly spherical shape in the limit of no rotation and, with rotation, to form
an axisymmetric oblate spheroid. A number of physical processes can disrupt the symmetry,
however, to produce quadrupolar gravitational waves from the stellar rotation.
Those processes include crustal distortions from cooling or accretion,
buried magnetic field energy and excitation of r-modes. Comprehensive reviews of
continuous gravitational wave emission mechanisms from neutron stars can be found in \cite{bib:Lasky_review,bib:GandG_review}

This article is organized as follows: Section~\ref{sec:data} describes the data set used, including steps taken to
mitigate extremely loud and relatively frequent instrumental glitches seen in the O3 LIGO data,
a phenomenon not seen in previous LIGO observing runs. Section~\ref{sec:method} briefly describes the
{\footnotesize\textsc{PowerFlux}}\ and loose-coherence algorithm used in this and previous searches. Section~\ref{sec:results} presents
the results of the analysis. Section~\ref{sec:conclusions} concludes with a discussion of the results and
prospects for future searches.

\section{Data sets used}
\label{sec:data}

Advanced LIGO consists of two detectors, one in Hanford, Washington (designated H1), and the other in Livingston, Louisiana (designated L1), separated by a $\sim$3000-km baseline~\cite{bib:LIGO}.
Each site hosts one, 4-km-long interferometer inside a vacuum envelope with the primary interferometer optics suspended by a cascaded, quadruple suspension system in order to isolate them from external disturbances.
The interferometer mirrors act as test masses, and the passage of a gravitational wave induces a differential-arm length change that is proportional to the gravitational-wave strain amplitude.

Advanced LIGO's first observing run (O1) occurred between September 12, 2015 and
January 19, 2016 and led to the discovery on September 14, 2015 of gravitational waves
from binary black hole (BBH) coalescences~\cite{bib:GW150914}.
The O2 observing run began November 30, 2016 and ended August 25, 2017 and included the first detection
of a binary neutron star (BNS) merger~\cite{bib:GW170817}.
The O3 run began April 1, 2019 and ended March 27, 2020, for which the first six months (April 1, 2019 to October 1, 2019),
prior to a 1-month commissioning break, is designated as the O3a epoch.
The analysis presented here is based primarily on the O3a data set, with data from the remainder of the run (O3b epoch)
used only for following up on promising signal candidates.
From the O1, O2 and O3a data sets, the LIGO Scientific Collaboration
and Virgo Collaboration have reported a total of 50 compact binary coalescences~\cite{bib:O1catalog,bib:O2catalog,bib:O3acatalog},
primarily BBH events, with two binary neutron star (BNS) detections~\cite{bib:GW170817,bib:GW190425}
and one potential binary neutron star / black hole detection~\cite{bib:GW190814}.

The Virgo interferometer~\cite{bib:Virgo} observed during August 2017 near the end of the O2 run
and throughout the O3 run. The Virgo data has not
been used in this analysis, however, because of an unattractive tradeoff in computational cost for sensitivity gain, given
the interferometer's higher noise level during the O3 run.

Prior to searching the O3a data for continuous wave (CW) signals, the quality of the data was assessed
and steps taken to mitigate the effects of instrumental artifacts. As in previous Advanced LIGO observing
runs~\cite{bib:linespaper}, instrumental ``lines'' (sharp peaks in fine-resolution run-averaged H1 and L1 spectra)
are marked, and where
possible, their instrumental or environmental sources identified~\cite{bib:O3aLinelist}. The resulting database of artifacts proved helpful in
eliminating spurious signal candidates emerging from the search. In general, the line multiplicity for H1 O3a data was
similar to that observed in the O2 run, while the line multiplicity for L1 O3a data was substantially reduced.

Another type of artifact observed in the O3a data for both H1 and L1 were relatively frequent and loud
``glitches'' (short, high-amplitude instrumental transients) with most of their spectral power lying below $\sim$500 Hz.
Although loud glitches have been observed in previous runs, their frequency in O3a was dramatically higher.
At present, investigations of the source of these glitches remain inconclusive. An effort to identify and mitigate
them during the October 2019 commissioning break (between the O3a and O3b epochs) did not succeed.
Unlike in previous LIGO data runs, the sheer spectral power in the glitches increased the effective broadband noise
floor of the data below $\sim$500 Hz quite substantially, as seen in run-averaged spectra computed from 1800s and 7200s
discrete Fourier transforms known as ``SFTs'' (for ``short Fourier transforms''). Most CW searches based on
summing strain spectral power from SFTs, including the {\footnotesize\textsc{PowerFlux}}\ program used here,
apply weightings that disfavor SFTs with high noise. Because the average time interval between loud glitches
in O3a data is comparable to or smaller than the coherence time of the SFTs, however, inverse noise weighting of SFTs
proved much less effective than in previous runs, especially for the 7200s SFTs.

To mitigate the effects of these glitches on O3a CW searches for signals below 475 Hz, a simple glitch-gating
algorithm was applied ~\cite{bib:gatingdoc,bib:detcharpaper} to excise the transients from the data.
For each 1/16-second for which a whitened version of the H1 and L1 strain data channels
had excess RMS power in the 25-50 Hz or 70-110 Hz bands, the strain channel was set to zero. To reduce artifacts
from discontinuous data, a 1/4-second half-Hann-window ramping from unity to zero was multiplied against the data stream
prior to each zeroed interval, and a 1/4-second half-Hann window ramping from zero back to unity was multiplied at the end
of zeroed intervals. This gating can be considered to be inverse-Tukey-windowed. Although the Tukey windowing
mitigates severe spectral artifacts, the resulting spectra still suffer visible spectral leakage
very near loud instrumental lines, such as from 60-Hz power mains and ``violin modes'' (near 500 Hz) due to
ambient vibrations of the silica fibers from which LIGO mirrors are suspended.
All gated intervals longer than 30 seconds are excluded from analysis, as are 7200s SFTs
containing total gate durations longer than 120 seconds. Altogether, the applied gating leads in this analysis to losses
of about 1\%\ and 11\%\ of the H1 and L1 observation times during the O3a epoch for search frequencies below 475 Hz.
Details of the gating, including validation that
low-frequency CW ``hardware injections'' (simulated signals imposed on interferometer mirrors during data taking, see
Section~\ref{sec:hwinjections}) are recovered with higher signal-to-noise ratio (SNR) in gated data than in the
original, ungated data can be found in a technical report~\cite{bib:gatingdoc}.

As discussed in section~\ref{sec:method}, the {\footnotesize\textsc{PowerFlux}}\ program searches many narrow frequency bands, where the SFT coherence time chosen depends on frequency. In general, longer coherence times give
improved sensitivity, but incur larger computational costs from the need to search more finely in parameter space.
In addition, efficiency loss from spectral leakage due to Doppler modulations from the Earth's motion increases
at higher frequency, especially for sources near the ecliptic plane. Given these tradeoffs, the same coherence-time
choices made for the O1 analysis~\cite{bib:O1allsky2} are chosen here, as shown in Table~\ref{tab:sft_info}.
The SFTs are created from the C01 calibrated strain data~\cite{bib:O3calibpaper}, using Hann windowing and 50\%\ overlap.
Upper limits on the calibration uncertainties over the 20-2000 Hz band in the O3a epoch are estimated
to be $<$ 7\%\ in magnitude and $<$ 4 deg in phase (68\% confidence interval).

\begin{table}[htb]
\begin{center}
  \begin{tabular}{lccc}\hline
    \T\B                                &  20-475 Hz      & 475-1475 Hz       & 1475-2000 Hz \\
    \hline\hline
    \T\B Data set\strut                       & Gated C01       & Ungated C01       & Ungated C01  \\
    \hline
    \T\B Coherence time\strut             & 7200s           & 3600s             & 1800s        \\
    \hline
    \multirow{2}*{Number of SFTs}  & H1:\T 2287        & H1: 5535          & H1: 11735     \\
                                   & L1:\B 2247        & L1: 5973          & L1: 12574     \\
\hline
\end{tabular}
  \caption{Information on the Short Fourier Transforms (SFTs) used in the search. All SFTs are
    Hann-windowed with 50\%\ overlap. The numbers of gated SFTs used below 475 Hz are affected
    by livetime loss from avoiding data stretches with high rates of gated gltiches, with a
    larger effect seen for L1 data.}
\label{tab:sft_info}
\end{center}
\end{table}

\section{Analysis Method}
\label{sec:method}

This search uses the {\footnotesize\textsc{PowerFlux}}\ program~\cite{bib:S4allsky,bib:O1allsky1,bib:O1allsky2,bib:S5allsky2} with loose-coherence follow-up
of outliers~\cite{bib:loosecoherence1,bib:loosecoherence2}. In brief, strain power is summed over many SFTs after
correcting for Doppler modulations, for a large bank of templates based on sky location, frequency, frequency
derivative and stellar orientation. The maximum strain powers detected over the entire sky and for all
frequency derivatives in each narrow frequency sub-band (see Table~\ref{tab:tuning_info}
for the frequency-dependent widths) define strict frequentist upper limits.
Initial outliers, defined by high values of SNR, are followed up in a
multi-stage search with increasing levels of effective coherence at each stage
to increase the expected SNR for true signals. Candidates passing all stages are
examined individually to determine their credibility as gravitational wave sources,
with further warranted follow-up carried out in the full O3 data set.

\subsection{Signal model and parameter space searched}
\label{sec:model}

The signal templates assume a classical model of a spinning neutron star with a time-varying quadrupole moment that produces
circularly polarized gravitational radiation along the rotation axis, linearly polarized radiation
in the directions perpendicular to the rotation axis and elliptical polarization for the general case.
The star's orientation, which determines the polarization, is parametrized by the inclination angle $\iota$ of its
spin axis relative to the detector line-of-sight and by the angle $\psi$ of the axis projection on the plane of the sky.
The linear polarization case ($\iota=\pi/2$) is the most unfavorable
because the gravitational wave flux impinging on the detectors
is smallest for an intrinsic strain amplitude $h_0$,
possessing eight times less incident strain power than
for circularly polarized waves ($\iota = 0,\>\pi$).

The strain signal model $h(t)$ for a periodic source is assumed to be the following function of time $t$:
\begin{eqnarray}
  h(t) & = & h_0\bigl(F_+(t, \alpha_0, \delta_0, \psi)\frac{1+\cos^2(\iota)}{2}\cos(\Phi(t))\vphantom{\frac{1+\cos^2(\iota)}{2}}\nonumber\\
  & & \quad\>+F_\times(t, \alpha_0, \delta_0, \psi)\cos(\iota)\sin(\Phi(t))\bigr)\textrm{~,}
\end{eqnarray}

\noindent where $h_0$ is the intrinsic strain amplitude, $\Phi(t)$ is the signal phase, $F_+$ and $F_\times$ characterize the detector responses to signals with ``$+$'' and ``$\times$''
quadrupolar polarizations \cite{bib:S4allsky}, and the sky location is described by right ascension $\alpha_0$ and
declination $\delta_0$.

In a rotating triaxial ellipsoid model for a
star at distance $r$ spinning at frequency $f_{\rm rot}$ about its (approximate) symmetry axis ($z$),
the amplitude $h_0$ can be expressed as
\begin{eqnarray}
  \label{eqn:hexpected}
h_0 & = & {4\,\pi^2G\epsilon\Izzf_{\rm GW}^2\over c^4r}  \\
    & = & [1.1\times10^{-24}]\Bigl[{\epsilon\over10^{-6}}\Bigr]\Bigl[{I_{\rm zz}\over I_0}\Bigr]\Bigl[{f_{\rm GW}\over1\>{\rm kHz}}\Bigr]^2
\Bigl[{1\>{\rm kpc}\over r}\Bigr],
\end{eqnarray}
where $I_0=10^{38}$ kg$\cdot$m$^2$ (10$^{45}$ g$\cdot$cm$^2$) is a nominal neutron star moment of inertia $I_{\rm zz}$ about $z$,
and the gravitational radiation is emitted at frequency $f_{\rm GW}=2\,f_{\rm rot}$.
The equatorial ellipticity $\epsilon$ is a convenient, dimensionless measure of stellar non-axisymmetry:
\begin{equation}
\label{eqn:ellipticity}
\epsilon \quad \equiv \quad {I_{xx}-I_{yy}\over I_{zz}}.
\end{equation}

The phase evolution of the signal is given in the reference frame of the Solar System barycenter (SSB) by the second-order approximation:
\begin{equation}
\label{eqn:phase_evolution}
\Phi(t)=2\pi\left(f_\textrm{source}\cdot (t-t_0)+f^{(1)}\cdot (t-t_0)^2/2\right)+\phi\textrm{~,}
\end{equation}
where $f_\textrm{source}$ is the SSB source frequency, $f^{(1)}$ is the first frequency derivative (which, when negative,
is termed the spin-down), $t$ is the SSB time, and the initial phase $\phi$ is
computed relative to reference time $t_0$. When expressed as a function of the local time of ground-based detectors,
Eq.~\ref{eqn:phase_evolution} acquires sky-position-dependent Doppler shift terms.

We search a frequency band 20-2000 Hz and a frequency derivative $\dot{f}$ range
of $[-\sci{1}{-8},+\sci{1}{-9}]~\mathrm{Hz}~\mathrm{s}^{-1}$.
Figure~\ref{fig:parameterspace} shows this parameter space coverage together with
those of previous all-sky searches of the LIGO O2 data using different methods.
Most natural sources are expected to have a negative first frequency derivative, as the energy lost in gravitational or
electromagnetic waves would make the source spin more slowly.
A small number of isolated pulsars in globular clusters exhibit slight apparent spin-up,
believed to arise from acceleration in the Earth's direction;  known apparent spin-up values have
magnitudes too small to prevent source detection with the zero-spin-down templates
used in this search, given a strong enough signal.
The frequency derivative can also be positive when the source
is affected by a strong slowly-varying Doppler shift, such as due to a long-period orbit with a companion.
A more exotic source of spin-up is gravitational wave superradiance from a boson cloud
in the vicinity of an isolated black hole~\cite{bib:bosoncloud}.

All known isolated pulsars spin down more slowly
than the maximum value of $|\dot{f}|_\textrm{\mbox{\tiny{max}}}$ used here, and as seen in the results
section, the equatorial ellipticity required for higher $|\dot{f}|$ is improbably high
for a source losing rotational energy primarily via gravitational radiation at low
frequencies. More plausible is a source with spin-down dominated by electromagnetic
radiation energy loss, but for which detectable gravitational radiation is also emitted.

One measure of gravitational wave detectability of known pulsars is the ``spin-down strain limit"
defined by equating inferred rotational kinetic energy loss to gravitational wave energy emission.
For the frequency band searched here, such limits range from as high as O($10^{-24}$) for energetic,
low-frequency young stars, such as the Crab and Vela pulsars, down to below $10^{-27}$ for high-frequency
millisecond pulsars~\cite{bib:O2targetedpaper}, where the highest such limit at high frequencies
is O($10^{-26}$).

\subsection{Methodology}
\label{sec:methodology}

The {\footnotesize\textsc{PowerFlux}}\ pipeline has a hierarchical structure
that permits systematic follow-up of loud outliers from the
initial stage. The later stages improve intrinsic strain sensitivity by
increasing effective coherence time while dramatically
reducing the parameter space volume over which the follow-up
is pursued. The pipeline uses loose coherence~\cite{bib:loosecoherence1}
with stages of improving refinement via steadily increasing effective
coherence times. Any outliers that survive
all stages of the search pipeline are examined manually for
contamination from known instrumental artifacts and for evidence of
contamination from a previously unknown single-interferometer artifact.
Those for which no artifacts are found are subjected to further
follow-up described below.

In the pipeline's initial stage, the main {\footnotesize\textsc{PowerFlux}}\
algorithm \cite{bib:S4allsky,bib:O1allsky1,bib:O1allsky2,bib:S5allsky2,bib:S5allsky1,bib:S6allsky}
establishes upper limits and produces lists of outliers.
The program sets strict frequentist upper limits
on detected strain power in circular and linear polarizations that
apply everywhere on the sky except for small regions near the ecliptic
poles, where signals with small Doppler modulations can be masked by
stationary instrumental spectral lines.
The procedure defining these excluded regions is described in \cite{bib:S5allsky2} and applies to less than $0.2$\% of the sky over the entire run,
where the precise shapes of the regions near the poles depend on assumed signal frequency and spin-down.
Initial outliers are defined by a
joint H1-L1 signal-to-noise ratio (SNR) greater than a threshold of 7, with consistency among
corresponding H1, L1 and joint H1-L1
outliers (criteria described in section~\ref{sec:outlier_followup}). These outliers are then followed up with
a loose-coherence detection pipeline \cite{bib:loosecoherence1,bib:loosecoherence2,bib:S5allsky2}, which is
used to reject or confirm the outliers.

The power calculation of the data can be expressed as a bilinear form of the input matrix $\left\{a_{t,f}\right\}$
constructed from the SFT coefficients with indices representing time and frequency:

\begin{equation}
P[f] = \sum_{t_1, t_2,D_i,D_j} a_{t_1, f+\Delta f(t_1)}^{(D_i)} a_{t_2, f+\Delta f(t_2)}^{(D_j)*} K_{t_1, t_2, f}^{D_iD_j}\textrm{~.}
\end{equation}

\noindent In this expression $\Delta f(t)$ is the detector-frame frequency drift due to the effects from both Doppler shifts and the first
frequency derivative. The sum is taken over all times $t$ corresponding to the midpoints of the SFT time intervals.
The kernel $K_{t_1, t_2, f}^{(D_iD_j)}$ includes the contribution of time-dependent SFT noise weights, antenna response, signal
polarization parameters, and relative phase terms~\cite{bib:loosecoherence1,bib:loosecoherence2} for detectors $D_{i,j}$ (= H1, L1).
Separate power sums are computed for H1, L1 and combined H1-L1 data.

The fast first-stage (stage 0) {\footnotesize\textsc{PowerFlux}}\ algorithm uses a kernel with diagonal terms only (including separate single-detector
contributions $D_i=D_j$). The second stage (stage 1)
increases effective coherence time while still allowing for controlled deviation in
phase~\cite{bib:loosecoherence1} via kernels that increase effective coherence length by inclusion of
limited single-detector, off-diagonal terms. The third stage (stage 2) maintains the stage-1 effective coherence time, but
adds SFT coefficients from H1 and L1 data coherently ($D_i\ne D_j$) to improve SNR and parameter resolution.

The effective coherence length is captured in a parameter $\delta$~\cite{bib:loosecoherence1}, which describes the degree of phase drift
allowed between SFTs. A value of $\delta=0$ corresponds to a fully coherent case,
and $\delta=\pi$ corresponds to incoherent power sums.

Depending on the terms used, the data from different interferometers can be combined incoherently (such as in stages 0
and 1, see Table \ref{tab:followup_parameters}) or coherently (as used in stage 2). The coherent
combination is more computationally expensive but improves parameter estimation.

\subsection{Upper limits determination}
\label{sec:upperlimit_determination}

The 95\% confidence-level (CL) upper limits presented in section~\ref{sec:results}
are reported in terms of the worst-case value of $h_0$ (linear polarization) and for
the most sensitive case of circular polarization.

These upper limits,
produced in stage 0, are based
on the overall noise level and largest outlier in strain found for every template in each narrow sub-band in
the first stage of the pipeline. Sub-bands are analyzed by separate instances of {\footnotesize\textsc{PowerFlux}}~\cite{bib:S5allsky2}

To allow robust analysis of the entire spectrum, including regions with severe spectral artifacts, a
{\em Universal statistic} algorithm~\cite{bib:universal_statistic,bib:O1allsky1} is used for establishing upper limits.
The algorithm is derived from the Markov inequality and shares its independence from the underlying noise distribution.
It produces upper limits less than $5$\% above optimal in case of Gaussian noise. In non-Gaussian bands,
it can report values larger than what would be obtained if the true underlying
distribution were known, but the upper limits are always
at least 95\% valid. Figure~\ref{fig:ul_vs_strain} shows results of a high-statistics
``software injections'' simulation run performed as described in~\cite{bib:S5allsky2}.
Correctly established upper limits lie above the dashed diagonal lines (defining equality between upper limit obtained and
true injection strain) in each panel, corresponding to
four selected sub-bands [SFT coherence times]: 20-60 Hz [7200s], 60-475 Hz [7200s], 475-1475 Hz [3600s] and 1475-2000 Hz [1800s].
Performance for the 7200s-SFT 20-60 Hz and 60-475 Hz bands are shown separately because of the proliferation
of spectral line artifacts below 60 Hz, primarily in the H1 data. The breakpoint frequencies of 475 Hz and 1475 Hz for
decreasing SFT coherence time (Table~\ref{tab:sft_info}) are those used in the O1 PowerFlux
search~\cite{bib:O1allsky1,bib:O1allsky2}, marking the starts of bands disturbed by 1st and 3rd violin mode harmonics.
Additional band-specific parameters for the initial stage of the search are listed in Table~\ref{tab:tuning_info}.

\begin{table*}[htb]
\begin{center}
  \begin{tabular}{lcccc}\hline
    \T\B                                         &  20-475 Hz      & 475-1475 Hz       & 1475-2000 Hz \\
    \hline\hline
    \T\B Frequency bin width                        &  0.139 mHz       & 0.278 mHz          & 0.556 mHz  \\
    \hline
    \T\B Number of spin-down templates\strut        &   66            &  34        &  18     \\
    \hline
    \T\B Spin-down step size (Hz/s)\strut        &  $\sci{1.692}{-10}$ & $\sci{3.333}{-10}$  & $\sci{6.667}{-10}$ \strut        \\
    \hline
     \T\B H1/L1 frequency mismatch tolerance (mHz)  &   2.5          &   2.5             & 2.5         \\
    \hline
     \T\B H1/L1 spin-down mismatch tolerance (Hz/s)   &  $\sci{3.0}{-10}$ &   $\sci{3.0}{-10}$    & $\sci{3.0}{-10}$        \\
    \hline
\end{tabular}
  \caption{Information on the band-specific parameters for the initial search (stage 0) and outlier follow-up stages.
    The number of spin-down steps and the spin-down step size refer to the templates used in the stage-0 search.
    The tolerances refer to the consistency requirements between H1 and L1 outliers used to define candidates selected
    for stage-1 follow-up. Other search parameters relevant to all search bands are described in the text.}
\label{tab:tuning_info}
\end{center}
\end{table*}

\begin{table*}[htbp]
\begin{center}
\begin{tabular}{llccccc}\hline
\T\B Stage & Instrument sum & {Phase coherence $\delta$} & \multicolumn{1}{c}{Spin-down step} & \multicolumn{1}{c}{Sky refinement} & \multicolumn{1}{c}{Frequency refinement} & \multicolumn{1}{c}{SNR increase} \\
 & & \multicolumn{1}{c}{rad} & \multicolumn{1}{c}{Hz/s} & $\alpha\times\delta_0$  &  & \multicolumn{1}{c}{\%}\\
\hline \hline
\multicolumn{7}{c}{\T\B 20-60\,Hz frequency range, 7200\,s SFTs, 0.0625\,Hz frequency bands} \\
  0 & Initial/upper limit incoherent & NA      & $\sci{1.692}{-10}$ & $1$   & $1/2$ & -- \\
  1 & incoherent                        & $\pi/2$ & $\sci{1.0}{-10}$   & $1/4\times1/4$ & $1/8$ & 20 \\
  2 & coherent                          & $\pi/2$ & $\sci{5.0}{-11}$   & $1/4\times1/4$ & $1/8$ & 10  \\
 \hline
\multicolumn{7}{c}{\T\B 60-475\,Hz frequency range, 7200\,s SFTs, 0.0625\,Hz frequency bands} \\
  0 & Initial/upper limit incoherent & NA      & $\sci{1.692}{-10}$ & $1$   & $1/2$ & -- \\
  1 & incoherent                        & $\pi/2$ & $\sci{1.692}{-10}$ & $1/4\times1/4$ & $1/8$ & 20 \\
  2 & coherent                          & $\pi/2$ & $\sci{5.0}{-11}$   & $1/4\times1/4$ & $1/8$ & 10  \\
 \hline
  \multicolumn{7}{c}{\T\B 475-1475\,Hz frequency range, 3600\,s SFTs, 0.125\,Hz frequency bands} \\
  0 & Initial/upper limit incoherent & NA      & $\sci{3.33}{-10}$  & $1$   & $1/2$ & -- \\
  1 & incoherent                        & $\pi/2$ & $\sci{3.33}{-10}$  & $1/4\times1/4$ & $1/8$ & 15 \\
  2 & coherent                          & $\pi/2$ & $\sci{5.0}{-11}$   & $1/4\times1/4$ & $1/8$ & 10  \\
 \hline
 \multicolumn{7}{c}{\T\B 1475-2000\,Hz frequency range, 1800\,s SFTs, 0.25\,Hz frequency bands} \\
  0 & Initial/upper limit incoherent & NA      & $\sci{6.67}{-10}$  & $1$   & $1/2$ & -- \\
  1 & incoherent                        & $\pi/2$ & $\sci{3.33}{-10}$  & $1/4\times1/4$ & $1/8$ & 15 \\
  2 & coherent                          & $\pi/2$ & $\sci{5.0}{-11}$   & $1/4\times1/4$ & $1/8$ & 10  \\
 \hline
\end{tabular}
\caption[Outlier follow-up parameters]{{\footnotesize\textsc{PowerFlux}}\ outlier follow-up parameters. Stage 1 and higher use a loose-coherence algorithm for demodulation. The sky (both right ascension $\alpha$ and declination $\delta_0$) and frequency refinement parameters are relative to values used in the stage 0 {\footnotesize\textsc{PowerFlux}}\ search}
\label{tab:followup_parameters}
\end{center}
\end{table*}

\begin{figure*}[htbp]
  \begin{center}
    \ifshowfigs
    \includegraphics[width=7.5in]{paramspace3.png}
    \fi
    \caption[Comparison of parameter space coverage]{Comparison in frequency and spin-down for this search with those of previous
      all-sky searches of LIGO O2 data. The shaded rectangle with vertical bars shows the 20--2000 Hz and
      $-10^{-8}$--$10^{-9}$ Hz/s range for this O3a search. The slightly larger rectangle with horizontal bars shows the
      region searched in the O2 data with the Frequency Hough method~\cite{bib:O2allsky,bib:O2allskyRome}. The smaller rectangle
      with crossed diagonal bars shows the region searched by the distributed-computing project Einstein@Home~\cite{bib:O2EatH}.
      The solid line at zero spin-down depicts the specialized O2 search for low-ellipticity millisecond pulsars using the
      Falcon method~\cite{bib:O2Falcon1,bib:O2Falcon2,bib:O2Falcon3} (the thickness of the line overstates the coverage in spin-down range).
      The dotted curves indicate contours of constant equatorial ellipticity $\epsilon$ = ($10^{-8}$, $10^{-7}$, $10^{-6}$, $10^{-5}$ and $10^{-4}$) for a star with stellar spin-down dominated by gravitational wave emission.}
    \label{fig:parameterspace}
  \end{center}
\end{figure*}

\begin{figure*}[htbp]
  \begin{center}
    \ifshowfigs
    \includegraphics[width=3.3in]{ULvalidation20-60.png}
    \includegraphics[width=3.3in]{ULvalidation60-475.png}
    \includegraphics[width=3.3in]{ULvalidation475-1475.png}
    \includegraphics[width=3.3in]{ULvalidation1475-2000.png}
    \fi
    \caption[Upper limit versus injected strain]{Illustration of {\footnotesize\textsc{PowerFlux}}\ upper limit validation for very low, low, mid and high frequencies.
      Each point in the upper left, upper right, lower left and lower right panels represents a separate
      injection in the 20-60\,Hz, 60-475\,Hz, 475-1475\,Hz and 1475-2000\,Hz frequency ranges, respectively.
      Each established upper limit (vertical axis) is compared against the
      injected strain value defined by the horizontal axis (diagonal line defines equality of upper limit obtained and
      true injected strain). These injections are for random polarizations. Corresponding scatterplots for linearly
      polarized injections confirm the desired 95\%\ coverage, regardless
      of polarization.}
    \label{fig:ul_vs_strain}
  \end{center}
\end{figure*}

\subsection{Outlier follow-up}
\label{sec:outlier_followup}

A follow-up search for detection is carried out for high-SNR outliers found in stage 0.
The outliers are subject to an initial coincidence test.  For each
outlier with $\textrm{SNR}>7$ in the combined H1 and L1 data, we require there
to be outliers in the individual detector data of the same small sky patch,
approximately square with side length $\sim$30 mrad $\times$ (100 Hz / frequency) that have $\textrm{SNR}>5$ and
match the parameters of the combined-detector outlier within 2.5~mHz in frequency and
$\sci{3}{-10}$\,Hz/s in spin-down.  The combined-detector SNR is additionally required to be above both
single-detector SNRs, in order to suppress single-detector instrumental artifacts, except for unusually
loud outliers (combined, H1 and L1 SNRs all greater than 20).

The identified outliers using combined data are then passed to the follow-up stage using a loose-coherence
algorithm~\cite{bib:loosecoherence1} with progressively
improved determination of frequency, spin-down, and sky location.

As the initial stage 0 sums only powers, it does not use the relative phase between interferometers,
which results in some degeneracy among sky position, frequency, and spin-down. The first loose-coherence
follow-up stage (1) demands greater temporal
coherence (smaller $\delta$) within each interferometer, which should
boost the SNR of viable outliers, but combines H1 and L1 power sums incoherently,
The subsequent stage (2) uses combined H1 and L1 data coherently, providing tighter bounds on outlier location.

Testing of the stages 0--2 pipeline is performed for frequency bands searched via software injections
using the same follow-up procedure. The recovery criteria also require that an outlier close to the true injection location (within 2.5~mHz in frequency
$f$, $\sci{3}{-10}$~Hz/s in spin-down and 28.5~rad$\cdot$Hz$/f$ in sky location) be found and successfully pass
through all stages of the detection pipeline.

Injection recovery efficiencies from simulations covering the major sub-bands (20-475 Hz, 475-1475 Hz, 1475-2000 Hz) are
shown in Fig.~\ref{fig:injection_recovery} for stages 0, 1 and 2, which confirm that 95\%\ signal
recovery is comparable to the 95\%\ upper limit in all bands, as desired, except for the region
below 60 Hz for which spectral line artifacts heavily contaminate the H1 data set.
Injections in vetoed frequency bands (see section~\ref{sec:outliers} and Table~\ref{tab:outlier_badbands}) are not included in these graphs.

As in previous {\footnotesize\textsc{PowerFlux}}\ analyses with loose-coherence follow-up~\cite{bib:O1allsky1,bib:O1allsky2},
only a mild influence from parameter mismatch is expected, as the parameters are chosen to accommodate the worst
few percent of injections. The follow-up procedure establishes very wide margins for outlier follow-up.
For example, when transitioning from the semi-coherent stage 0 to the loose-coherence stage 1 below 475 Hz, the effective
coherence length increases by a factor of 4. The average true signal SNR should then increase by more than $40$\%.
But the threshold used in follow-up is only 15--20\%, depending on frequency, which accommodates unfavorable noise conditions,
template mismatch, and detector artifacts.

Although a similar prior analysis~\cite{bib:O1allsky1,bib:O1allsky2} of O1 data used additional stages (3, 4) of
loose coherence, we choose after stage 2 to inspect candidates manually. The small number not found to be contaminated by
obvious instrumental artifacts are followed up in the full O3 data set by
a search similar to the {\footnotesize\textsc{PowerFlux}}\ O3a stage 1, but using much finer spin-down stepping of $\sci{1}{-11}$ Hz/s, and refinement factors of 1/4
for both sky (right ascension and declination each) and frequency stepping, to exploit the improved SNR and parameter
resolution possible in the nearly doubled observation span.
Any outlier surviving this full-O3 follow-up is explored via
a more sensitive method, implemented in {\footnotesize\textsc{PyFstat}}~\cite{bib:PyFstatMethod1,bib:PyFstatZenodo} which
uses Markov Chain Monte Carlo (MCMC) techniques~\cite{bib:AshtonPrix,bib:PyFstatMethod2}
to explore small regions of parameter space in successive stages of increasing coherence times.

The full-O3 {\footnotesize\textsc{PowerFlux}}\ follow-up search uses nearly 11 months of
data (compared to the 6-month O3a period), leading to an expected increase in SNR by a factor of approximately $\sqrt{11/6} = 1.35$,
with an additional $\sim$40\%\ increase from the use of stage-1 loose coherence, leading to a nominal expectation of SNR gain of
a factor of $\sim$1.9. Software simulations of signals with amplitudes as weak as 1/3 of the upper limit in a given band,
however, reveal large fluctuations above and below this expectation for all search bands.
These simulations suggest applying a conservative (negligible efficiency loss)
threshold on the SNR increase of only 20\%, which is nonetheless effective in suppressing most noise artifacts.

The final follow-up method applied to any survivors of the full-O3 loose-coherence check
uses the Python-based {\footnotesize\textsc{PyFstat}}~\cite{bib:PyFstatMethod1,bib:PyFstatZenodo} software
infrastructure to combine a
MCMC approach~\cite{bib:AshtonPrix,bib:PyFstatMethod2} with semi-coherent summing of the well known $\mathcal{F}$-statistic\ detection
statistic~\cite{bib:JKS}. In this approach, the parameter space near to those values from a stage-2 outlier is sampled
randomly according to a certain probability density function determined by the $\mathcal{F}$-statistic\ likelihood function.
We closely follow an implementation~\cite{bib:PyFstatMethod2} applied to a recent analysis of O2 data to follow up
on outliers emerging from several prior O2 continuous wave searches.

Briefly, the O3a observation time is divided into $N_{\textrm{seg}}$\ segments, for each of which the $\mathcal{F}$-statistic\ is computed
over a coherence time approximately equal to the observation time divided by $N_{\textrm{seg}}$. For each point sampled
in parameter space, the sum of the $\mathcal{F}$-statistic\ values is computed to form a total detection statistic.
This procedure is repeated, decreasing the number of segments (increasing the coherence time),
using the resulting MCMC-maximized  $\mathcal{F}$-statistic\ sum as the seed for the next stage, with a
consequent reduction in parameter space volume searched.
For this O3a analysis we choose five successive stages of follow-up with decreasing values of $N_{\textrm{seg}}$\ = 500 (coherence
time of 0.36 day), 250, 55, 5 and 1. A random sample of 600 off-source sky locations having the same declination as the putative
signal direction, but separated by more than 90 deg from that direction, is used to determine a non-signal expectation for
the background distributions in the same frequency band~\cite{bib:PyFstatMethod2}.
A Bayes factor $B_{\rm SN}$ is computed from the change in
$\mathcal{F}$-statistic\ values for a nominal signal compared to the  empirical background  distribution in the last stage ($N_{\textrm{seg}}$=5 to $N_{\textrm{seg}}$=1).

\begin{figure*}[htbp]
  \begin{center}
    \ifshowfigs
    \includegraphics[width=3.3in]{EfficVsRelstrain_20_60Hz.jpg}
    \includegraphics[width=3.3in]{EfficVsRelstrain_60_475Hz.jpg}
    \includegraphics[width=3.3in]{EfficVsRelstrain_475_1475Hz.jpg}
    \includegraphics[width=3.3in]{EfficVsRelstrain_1475_2000Hz.jpg}
    \fi
    \caption[Injection recovery]{Injection (software simulations) recovery efficiencies
      in the 20-60\,Hz, 60-475\,Hz, 475-1475\,Hz and 1475-2000\,Hz frequency bands are shown
      in the upper left, upper right, lower left and lower right panels, respectively, for stages 0, 1 and 2 of the search.
      The injected strain divided by the 95\%\ CL upper limit in its band (without injection)
      is shown on the horizontal axis.
      The percentage of surviving injections is shown on the vertical axis,
      with a horizontal dashed line drawn at the $95$\% level. The vertical dashed line
      marks a relative strain of unity. Ideally, the recovery efficiencies should lie
      to the left of the vertical line or above the horizonal line. One observes, however,
      that the ideal recovery efficiency is significantly degraded below 60 Hz, where large instrumental
      artifacts in the H1 data and only modest expected Doppler modulations make clean signal recovery more challenging.
      Most of the degradation stems from the band below 27 Hz for which H1 line contamination is severe and for
      which the H1 noise floor is substantially higher than the L1 noise floor.}
   \label{fig:injection_recovery}
  \end{center}
\end{figure*}

\subsection{Hardware injections}
\label{sec:hwinjections}

During the O3 run 18 hardware injections were used to simulate
particular CW signals, as part of detector response validation, including long-term
phase fidelity. The injections were imposed via radiation pressure from auxiliary lasers~\cite{bib:O1hwinjections}.
For reference, Table~\ref{tab:hwinjections} lists the key source parameters for the
14 injections relevant to this analysis (labeled Inj0-Inj12 and Inj14), namely those that simulate isolated
neutron stars with nominal frequencies between 20 and 2000 Hz.
In general, the injection amplitudes used in the O3 run are substantially lower than those used in
previous observing runs. For this reason, some injections with identical source parameters, except for stronger signal amplitudes,
were detected with high SNR in the O2 run, but are not detected in this analysis,
although they can be found with targeted matched-filter analyses.

For each injection in Table~\ref{tab:hwinjections}, the column labeled ``UL sig bin''
gives the 95\%\ UL for the corresponding band. Ideally,
that value should exceed the true injected amplitude $h_0$ for at least 95\%\ of the injections.
In this case, that statement holds for all 14 injections. The column labeled ``UL ctrl bins'' shows the average
95\%\ UL for the six nearest neighboring frequency subbands as a rough guide to the expected value in
the absence of an injection (or true signal).

The last column states whether or not the injection survives all stages of loose-coherence follow-up.
As hoped, all five injections with a true $h_0$ amplitude above the expected background estimation of the
UL do survive these follow-up stages, as do two additional injections with somewhat smaller amplitudes.

\begin{table*}[htbp]
\begin{center}
\begin{tabular}{lD{.}{.}{6}rD{.}{.}{5}D{.}{.}{5}rrrc}
\hline
Label & \multicolumn{1}{c}{Frequency} & \multicolumn{1}{c}{Spindown} & \multicolumn{1}{c}{$\textrm{RA}_{\textrm J2000}$} & \multicolumn{1}{c}{$\textrm{DEC}_{\textrm J2000}$} & \multicolumn{1}{c}{$h_0$} & \multicolumn{1}{c}{UL} & \multicolumn{1}{c}{UL} & \multicolumn{1}{c}{Detected?}\\
 & \multicolumn{1}{c}{Hz} & \multicolumn{1}{c}{nHz/s} & \multicolumn{1}{c}{degrees} & \multicolumn{1}{c}{degrees} & \multicolumn{1}{c}{true} & \multicolumn{1}{c}{sig bin}  & \multicolumn{1}{c}{ctrl bins} & \\
\hline \hline
Inj0   &  265.575086  & $\sci{-4.15}{-3}$   &   71.55193     &  -56.21749 & $\sci{6.12}{-26}$ & $\sci{1.8}{-25}$ & $\sci{1.8}{-25}$ & {\color{blue}Yes} \\
Inj1   &  848.937350  & $\sci{-3.00}{-1}$   &   37.39385     &  -29.45246 & $\sci{5.47}{-25}$ & $\sci{8.1}{-25}$ & $\sci{3.6}{-25}$ & {\color{blue}Yes}\\
Inj2   &  575.163506  & $\sci{-1.37}{-4}$   &  215.25617     &    3.44399 & $\sci{7.59}{-26}$ & $\sci{2.7}{-25}$ & $\sci{2.7}{-25}$ & No\\
Inj3   &  108.857159  & $\sci{-1.46}{-8}$   &  178.37257     &  -33.4366  & $\sci{1.30}{-25}$ & $\sci{1.9}{-25}$ & $\sci{1.9}{-25}$ & No\\
Inj4   & 1390.806586  & $\sci{-2.54}{+1}$   &  279.98768     &  -12.4666  & $\sci{1.07}{-24}$ & $\sci{4.9}{-25}$ & $\sci{4.9}{-25}$ & {\color{blue}Yes\footnote{True spin-down value outside of nominal search range. Injection is recovered when search range is extended for this band.}}\\
Inj5   & 52.808324  & $\sci{-4.03}{-9}$     &  302.62664     &  -83.83914 & $\sci{3.99}{-25}$ & $\sci{5.7}{-25}$ & $\sci{3.2}{-25}$ & {\color{blue}Yes}\\
Inj6   &  145.444975 & $\sci{-6.73}{0\mathrm{\;\;\;}}$     &  358.75095     &  -65.42262 & $\sci{3.84}{-25}$ & $\sci{4.3}{-25}$ & $\sci{2.2}{-25}$ & {\color{blue}Yes}\\
Inj7   & 1220.434717 & $\sci{-1.12}{0\mathrm{\;\;\;}}$     &  223.42562     &  -20.45063 & $\sci{1.65}{-25}$ & $\sci{4.5}{-25}$ & $\sci{4.4}{-25}$ & No\\
Inj8   &  190.100215 & $\sci{-8.65}{0\mathrm{\;\;\;}}$     &  351.38958     &  -33.41852 & $\sci{1.30}{-25}$ & $\sci{1.8}{-25}$ & $\sci{1.7}{-25}$ & {\color{blue}Yes} \\
Inj9   &  763.847316 & $\sci{-1.45}{-8}$    &  198.88558     &   75.68959 & $\sci{1.30}{-25}$ & $\sci{3.1}{-25}$ & $\sci{3.0}{-25}$ & No\\
Inj10  &   26.332768 & $\sci{-8.50}{-2}$    &  221.55565     &   42.87730 & $\sci{6.26}{-25}$ & $\sci{1.9}{-24}$ & $\sci{1.2}{-24}$ & No\\
Inj11  &   31.424704 & $\sci{-5.07}{-4}$    &  285.09733     &  -58.27209 & $\sci{3.17}{-25}$ & $\sci{9.4}{-25}$ & $\sci{7.8}{-25}$ & No\\
Inj12  &   37.805210 & $\sci{-6.25}{0\mathrm{\;\;\;}}$     &  331.85267     &  -16.97288 & $\sci{2.63}{-25}$ & $\sci{4.8}{-25}$ & $\sci{4.9}{-25}$ & No\\
Inj14  & 1991.092292 & $\sci{-1.00}{-3}$     & 300.80284     &  -14.32394 & $\sci{1.83}{-24}$ & $\sci{2.1}{-24}$ & $\sci{8.1}{-25}$ & {\color{blue}Yes}\\
\hline
\end{tabular}\\
\caption[Parameters of hardware injections]{Key parameters of the hardware-injected simulated isolated-source continuous wave
  signals during the O3 data run (epoch GPS $1238166018$) with nominal frequencies
  between 20 and 2000 Hz, along with upper limits in the nominal signal bins and in the six nearest
  control bins. The last column states whether or not the injection is detected in outlier follow-up.
  Most of the injection amplitudes are significantly weaker than those used in prior observing runs,
  leading to fewer detections than in past analyses. The amplitudes quoted here are the {\it intended} values;
  uncertainties in the approximate actuation functions used in the online injection system may be as large
  as a few percent, depending on frequency. Also, a now-known error in the actuation function used for L1 injections
  led to an unintended time delay of approximately 122 $\mu$s, which corresponds to a phase lag in the highest-frequency
  injection (2991 Hz) of $\sim$130$^\circ$ and degrades coherent H1-L1 injection recovery.}
\label{tab:hwinjections}
\end{center}
\end{table*}

\section{Search results}
\label{sec:results}

Carrying out the stage-0 analysis described above leads to a set of all-sky upper limits
(95\%\ CL) on strain amplitude for worst-case, linear polarization (relative stellar spin
orientation) and for best-case, circular polarization. That analysis also leads to an
initial set of outliers for follow-up with later analysis stages.
Whether or not a particular narrow frequency band contains an outlier,
the upper limits obtained remain valid.
As shown in Table~\ref{tab:hwinjections}, hardware injections were reliably recovered when
their injected strain was at least 0.7 times the recorded upper limit (and in one case as low as 0.34 times the upper limit).
The upper limits obtained from stage-0 analysis
are presented in section~\ref{sec:upperlimits}, along with corresponding astrophysical
sensitivities. Section~\ref{sec:outliers} presents the results of outlier follow-up.

\subsection{Upper limits and sensitivity}
\label{sec:upperlimits}

Figure~\ref{fig:upper_limits} shows the upper limits obtained in this search as three curves, along with results from
other all-sky searches in the O2 data set~\cite{bib:O2allsky,bib:O2allskyRome,bib:O2EatH,bib:O2Falcon1,bib:O2Falcon2,bib:O2Falcon3}.
The upper (blue) curve shows the upper limits for a worst-case (linear) polarization. The lowest (red) curve shows
upper limits for an optimally oriented source (circular polarization). Both curves represent strict frequentist upper
limits on these worst-case and best-case orientations with respect to location on the sky, spin-down and frequency
within each narrow search band. For sensitivity comparisons with the previous O2 results,
the intermediate (green) curve represents approximate population-averaged upper limits
(over random sky locations and polarizations) and is
derived from the circular polarization curve with a simple scale factor (2.3), based on injection studies in test bands.

Each linear-polarization or circular-polarization point in Fig.~\ref{fig:upper_limits} represents a maximum over the sky,
except for a small excluded portion of the sky near the ecliptic poles,
which is highly susceptible to detector artifacts due to stationary frequency evolution produced by the combination of
frequency derivative and Doppler shifts~\cite{bib:S5allsky2}.

The O3a results presented here improve upon the previous O2 results in strain
sensitivity with factors ranging from $\sim$1.2 at signal frequencies of
$\sim$100 Hz to $\sim$2.0 at higher frequencies below the 1475-Hz breakpoint. The improvements at high frequencies come in part from
the improved detector noise achieved with quantum squeezing~\cite{bib:O3squeezing} and in part from using
longer Fourier transform coherence times than were used in the similar O2 analyses.
At frequencies below $\sim$40 Hz, the O2 data was badly contaminated by
instrumental lines~\cite{bib:linespaper}, leading to poorer upper limits, and in some bands, precluding reliable upper limits
using the method of~\cite{bib:O2allskyRome}. The dramatically improved sensitivies obtained for low frequencies in this analysis come in part from
mitigation of many spectral lines between the O2 and O3 runs, especially in the L1 interferometer, and in part from the
use here of the ``universal statistic'' described in section~\ref{sec:upperlimit_determination}.

\begin{figure*}[htbp]
  \begin{center}
    \ifshowfigs
    \includegraphics[width=7.5in]{Limits_O3avsO2_020_2000_Hz_semlgy_withfalcon.png}
    \fi
    \caption{Upper limits on gravitational strain amplitude for this O3a analysis and for previous O2 analyses.
      The dark blue curve shows worst-case (linearly polarized)
      $95$\%~CL upper limits in analyzed ($62.5$\,mHz, $125$\,mHz, $250$\,mHz) sub-bands for the
      broad (20-475 Hz, 475-1475 Hz, 1475-2000 Hz) search bands. The lowest, purple curve shows upper limits assuming a
      circularly polarized source. The middle (green) curve shows approximate
      population-averaged all-sky upper limits inferred from the circularly polarized limits.
      For comparison, the uppermost yellow-green diamonds show previous population-averaged 95\%\ CL all-sky upper limits
      derived from LIGO O2 data using the Frequency Hough method~\cite{bib:O2allsky,bib:O2allskyRome}.
      The orange curve shows 90\%\ CL upper limits from the O2 Einstein@Home search~\cite{bib:O2EatH}
      over the band 20-585.15 Hz and over a spin-down range that is 26\%\ of the range used here.
      The light blue curve shows 95\%\ CL upper limits from the O2 Falcon search~\cite{bib:O2Falcon1,bib:O2Falcon2,bib:O2Falcon3}
      over the 20-2000 Hz band and a severely restricted spin-down range.
      See Fig.~\ref{fig:parameterspace} for a comparison
      of parameter space coverage from these various searchers.}
    \label{fig:upper_limits}
  \end{center}
\end{figure*}

One can recast these upper limits on source amplitude as lower limits on the range at which
neutron stars with assumed equatorial ellipticity can reside.
Figure~\ref{fig:range} shows the ranges (kpc) of this search (circular-polarization and population-averaged)
{\it vs.}\ frequency for assumed equatorial ellipticity values
$\epsilon$ = $10^{-4}$, $10^{-6}$ and $10^{-8}$ and under the assumption that the {\it total} stellar spin-down magnitude does
not exceed the $10^{-8}$ Hz/s maximum used in this search.
The population-averaged ranges from prior O2 searches of comparable parameter space
coverage~\cite{bib:O2allsky,bib:O2allskyRome} are also shown for
these ellipticities. One sees a gain in range with respect to O2 analyses over the full frequency band,
with the greatest improvements seen at highest frequencies, as expected. The number of accessible neutron
stars is expected to rise at least quadratically as the range extends out into the plane of the galaxy and
approaches the denser regions near its center, at a distance of $\sim$8.5 kpc, enhancing the detection
probabilities from these range gains.

For reference, Fig.~\ref{fig:range} also shows this search's implied maximum ranges {\it vs.}\ frequency for a
``gravitar'' (star with spin-down dominated by gravitational radiation) for
maximum spin-down magnitudes of $10^{-8}$, $10^{-10}$ and $10^{-12}$ Hz/s.
If a neutron star's spin-down has major contributions from
electromagnetic radiation, the range to which this search is sensitive for fixed $\epsilon$ may be
reduced by the choice of maximum spin-down permitted.
Conversely, searching for spin-down magnitudes much higher than expected for a gravitar of given $\epsilon$
could allow detection of a star with more realistic emission contributions.

\begin{figure*}[htbp]
  \ifshowfigs
  \includegraphics[width=7.5in]{Range_O3avsO2_020_2000_Hz.png}
  \fi
  \caption[Search range]{
    Ranges (kpc) of the {\footnotesize\textsc{PowerFlux}}\ search for neutron stars
    spinning down solely due to gravitational radiation (``gravitars'') under different assumptions.
    The three sets of three curves (purple, green dotted, yellow-green diamonds)
    that generally rise with frequency are the ranges for which the O3a circular-polarization,
    O3a population-averaged and O2 population-averaged Frequency
    Hough~\cite{bib:O2allsky,bib:O2allskyRome} upper limits apply for three
    different assumed equatorial ellipticities of $\epsilon$ = $10^{-4}$, $10^{-6}$ and $10^{-8}$.
    In each case, the curve stops at the frequency at which the implied spindown magnitude
    at the corresponding range exceeds the maximum range allowed by the maximum spindown magnitude ($10^{-8}$ Hz/s)
    used in these analyses.
    The three black curves (solid, dashed, dot-dashed) that peak below 100 Hz represent the maximum possible ranges
    for which the O3a circular-polarization limits apply under the assumption that the
    maximum spindown magnitude is ($10^{-8}$, $10^{-10}$, $10^{-12}$) Hz/s.
    Searching for high spindown magnitudes simultaneously probes higher
    ellipticities and larger regions of the galaxy, reaching well beyond the galactic center ($\sim$ 8.5 kpc)
    at high ellipticities and low frequencies. More realistically, a neutron star's spin-down may be dominated
    by electromagnetic radiation, in which case searching for high spindown magnitude achieves a shorter search
    range, but enables detection of stars that would not be detected under the gravitar assumption.
    (color online)}
  \label{fig:range}
\end{figure*}

\subsection{Outliers}
\label{sec:outliers}

Outliers seen in the stage-0 search are followed up with loose-coherence stages of
increasing effective coherence time, as described in section~\ref{sec:outlier_followup},
with the SNR expected to increase for true signals. Table~\ref{tab:outlier_counts}
shows the counts of outliers seen at each stage for major sub-bands.
In the presence of Gaussian noise and no signal, one expects the outlier counts to
decrease monotonically with increasing stage as SNR increases are required for each advancement.
For particular sub-bands, however, one can see count increases from at least two contributions,
both associated with the finer sampling of parameter space in successive stages.
1) Hardware injections (section~\ref{sec:hwinjections}) from
a simulated signal naturally satisfy the SNR increase requirements demanded of a signal;
and 2) stationary line artifacts of finite bandwidth can be compatible with
signal templates having limited Doppler modulation or having partial cancellation
between seasonal modulation and assumed frequency derivative
for a certain region of the sky~\cite{bib:S4allsky}.

In the extreme, the number of outliers produced by a particular instrumental artifact
can be so large as to make systematic follow-up impracticable and pointless.
Table~\ref{tab:outlier_badbands} shows particular frequency bands (about 13\%\ of the original search band) for which no
outlier follow-up is attempted for the initial or later stages because of instrumental artifacts
that lead to such an unmanageable flood of initial or later-stage outliers. The table includes the stage for which the
vetos are applied in order to reduce artifact-induced outlier counts.
Artifacts include loud mechanical resonances, such as higher harmonics of violin modes, and especially
loud hardware injections. The widest bands excluded lie in the regions of test-mass violin modes and their higher harmonics.
Figure~\ref{fig:vetobands} shows these vetoed regions over the full search band, with a magnification of an example
100-Hz band (1900-2000 Hz) that includes a ``forest'' of 4th harmonics of violin modes
from both interferometers.
For the narrow search bands affected by loud hardware injections, only outliers within 1.0 radian
of the injection's sky location are excluded, except for the top 20 (highest SNR),
which are followed up to verify successful injection recovery.

\def\mc#1{\multicolumn{1}{c}{#1}}
\begin{table*}[htb]
\begin{center}
  \begin{tabular}{lrrrr}\hline
    \T\B                                       &    \mc{20-60  Hz}        &  \mc{60-475 Hz}         & \mc{475-1475 Hz}           & \mc{1475-2000 Hz} \\
    \hline\hline
    \T\B Stage 0                                     &  {\color{magenta}\arabic{VlowBandcountzero}} &  {\color{magenta}\arabic{LowBandcountzero}}  & {\color{magenta}\arabic{MediumBandcountzero}}   & {\color{magenta}\arabic{HighBandcountzero}} \\
    \hline
    \T\B Stage 1                                     &  {\color{magenta}\arabic{VlowBandcountone}}  &  {\color{magenta}\arabic{LowBandcountone}}   & {\color{magenta}\arabic{MediumBandcountone}}    & {\color{magenta}\arabic{HighBandcountone}}  \\
    \hline
    \T\B Stage 2                                     &  {\color{magenta}\arabic{VlowBandcounttwo}}  &  {\color{magenta}\arabic{LowBandcounttwo}}   & {\color{magenta}\arabic{MediumBandcounttwo}}    & {\color{magenta}\arabic{HighBandcounttwo}}   \\
    \hline\hline
    \T\B HW injections                               &  {\color{magenta}\arabic{VlowBandcountHWinj}}  &  {\color{magenta}\arabic{LowBandcountHWinj}}   & {\color{magenta}\arabic{MediumBandcountHWinj}}    & {\color{magenta}\arabic{HighBandcountHWinj}}   \\
    \T\B Visible artifacts                           &  {\color{magenta}\arabic{VlowBandcountStrainhist}}  &  {\color{magenta}\arabic{LowBandcountStrainhist}}   & {\color{magenta}\arabic{MediumBandcountStrainhist}}    & {\color{magenta}\arabic{HighBandcountStrainhist}}   \\
    \T\B Full-O3 follow-up                           &  {\color{magenta}\arabic{VlowBandcountFullO3}}  &  {\color{magenta}\arabic{LowBandcountFullO3}}   & {\color{magenta}\arabic{MediumBandcountFullO3}}    & {\color{magenta}\arabic{HighBandcountFullO3}}   \\
    \hline\hline
    \T\B Full-O3 survivors                           &  0 & 0 & 0 & 21 \\
    \hline
  \end{tabular}
  \caption{Counts of outliers surviving different stages of the hierarchical
    follow-up for four frequency bands.
    Survivors from stage-2 follow-up are broken down into the categories of hardware injections,
    visible artifacts identified via strain histograms and candidates for loose-coherence follow-up in the
    full-O3 data. Only one cluster of O3a outliers (near 1663.6 Hz) survives the full-O3 fine-grained follow-up, leading to
    a cluster of 21 full-O3 outliers.
    Additional follow-up of these weak candidates, which appear in a visibly disturbed band of L1 data,
    was carried out using the MCMC {\footnotesize\textsc{PyFstat}}\ follow-up procedure described in section~\ref{sec:outlier_followup}.
    The {\footnotesize\textsc{PyFstat}}\ follow-up confirmed these 21 outliers to be uninteresting.}
\label{tab:outlier_counts}
\end{center}
\end{table*}

\begin{table*}[htb]
  \begin{center}
    \input{band_linevetos_expanded_v4.tex}
    \caption{Frequency bands vetoed from outlier follow-up because of excessive
      outlier counts from instrumental artifacts. The stage at
      which outlier followup is aborted is also shown for each band.
      Altogether, about 257 Hz out of the 1980-Hz original search band was vetoed, or about 13\%.
      Figure~\ref{fig:vetobands} shows these vetoed bands graphically.}
\label{tab:outlier_badbands}
\end{center}
\end{table*}

Nearly all outliers that survived all stages of the loose-coherence
follow-up correspond to hardware injections (see Table~\ref{tab:hwinjections})
or lie in highly disturbed bands, for which contamination of the putative signal by an instrumental spectral line is
apparent. To identify these contaminations, we construct ``strain histograms'' in which the summed power over the
observation period from a simulation of the nominal
signal candidate is superposed on a background estimate of the noise estimated via interpolation between neighboring
frequency bands.
Except for signal templates with high-magnitude spin-downs, the histograms typically display at least one ``horn''
(narrow peak) from
an epoch during the 6-month O3a period when the orbitally modulated frequency is relatively stationary. We discard
outliers for which the signal template's shape aligns with a spectral artifact known to be instrumental or appearing loudly
in one detector but not the other. Before visual inspection, outliers are clustered in frequency, spin-down and sky location,
with the 10 highest-SNR members receiving scrutiny. For reference, Table~\ref{tab:histogramvetoes} lists the parameters of
the single loudest outlier in each cluster discarded after inspecting strain histograms. Figure~\ref{fig:strainhistogram}
shows an example of a histogram for an outlier contaminated by an L1 spectral line at 1740.39 Hz.

Any remaining survivors (see Table~\ref{tab:followups})
are followed up in the full O3 data set, as described in section~\ref{sec:outlier_followup},
and required to yield an SNR increase of at least 20\%\ relative to the original stage-0 values.
All but one cluster of outliers fail to satisfy this (conservative) requirement
from the fine-grained full-O3 stage-1 follow-up.

The one cluster with (21) full-O3 outliers passing the requirement occurs near 1663.6 Hz (see Table~\ref{tab:followups}) for
which the maximum SNR increase observed is 29\%.
Examination of the corresponding strain histogram shows evidence of a mild non-Gaussian artifact in L1
data not visible in H1 data, consistent with an instrumental source, but not at a strong enough level to justify visual discarding.

Hence this outlier is subjected
to the MCMC {\footnotesize\textsc{PyFstat}}\ follow-up procedure described in section~\ref{sec:outlier_followup}.
The resulting (log) Bayes factor comparing signal to empirical noise expectation
for the maximum $\mathcal{F}$-statistic\ in the last stage
is found to be log$_{10}$($B_{\rm SN}$) = 1.1. Based on the analysis in~\cite{bib:PyFstatMethod2},
this value is lower by about 50 than that expected for a standard CW signal of the
outlier's putative strain amplitude in the O3a data set.

We conclude that there is no significant evidence for a continuous wave signal from this search.

\begin{table}[htb]
  \begin{center}
    \input{strainhist_vetos.tex}
    \caption{Parameters of the outliers surviving stage-2 follow-up, but discarded after visual inspection of a strain histogram confirmed instrumental
      contamination. Each outlier listed is the loudest after clustering in frequency, spin-down and sky location. Figure~\ref{fig:strainhistogram} shows
      a strain histogram example.}
\label{tab:histogramvetoes}
\end{center}
\end{table}

\begin{table}[htb]
  \begin{center}
    \input{followups.tex}
    \caption{Parameters of the outliers surviving stage-2 follow-up and followed up in the full-O3 data set.
      Each outlier listed is the loudest after clustering in frequency, spin-down and sky location.}
\label{tab:followups}
\end{center}
\end{table}

\begin{figure*}[htbp]
  \ifshowfigs
  \includegraphics[width=6.5in]{O3a_H1_L1_spectra_20_000_2000_000_Hz_withvetoes_v3.png}
  \includegraphics[width=6.5in]{O3a_H1_L1_spectra_1900_000_2000_000_Hz_withvetoes_v3.png}
  \fi
  \caption[Vetoed bands]{
    Shaded regions (yellow) overlaid on the run-averaged spectra (inverse-noise weighted) from the
    H1 (magenta) and L1 (blue) O3a data sets
    indicate vetoed bands for which outlier follow-up is impeded by artifacts. The widest affected
    bands correspond to test-mass ``violin modes'' (ambient vibrations of the silica fibers from which LIGO mirrors are suspended)
    near multiples of 500 Hz.
    {\it Top panel:} Full 20-2000 Hz band.
    {\it Bottom panel:} Magnification of example 1900-2000 Hz band, for which line artifacts
    are dominated by the 4th harmonics of violin modes.
    (color online)}
  \label{fig:vetobands}
\end{figure*}

\begin{figure*}[htbp]
  \ifshowfigs
  \includegraphics[width=6.5in]{SampleStrainHistogram.png}
  \fi
  \caption[Strain Histogram]{
    Example of a ``strain histogram'' graph used in vetoing outliers surviving stage-2 follow-up
    for which instrumental contamination is apparent. The solid curves show the O3a-run-averaged
    H1 (red) and L1 (blue) amplitude spectral densities in a narrow band containing an artifact at 1740.39 Hz.
    The dotted curves show histograms of expected strain excess from H1 (black) and L1 (magenta)
    signal templates added to
    smooth backgrounds interpolated from neighboring frequency bands.

    In this depiction, the strain amplitude of the signal template has been
    magnified to make its structure clear. The large excess power in the L1 data, not seen in the H1 data, despite
    comparable strain sensitivities and comparable sidereal-averaged antenna pattern sensitivities to any point in the sky,
    excludes an astrophysical source for the L1 artifact. The fact that the artifact aligns in frequency
    with the putative signal's template peak in power confirms contamination of the outlier from an
    instrumental source.
    (color online)}
  \label{fig:strainhistogram}
\end{figure*}

\section{\label{sec:conclusions}Conclusions}

We have performed the most sensitive all-sky search to date for continuous gravitational waves
in the range 20-2000\,Hz while probing spin-down magnitudes as high as $10^{-8}$ Hz/s, using the {\footnotesize\textsc{PowerFlux}}\ search program with loose coherence.
The overall improvements in strain sensitivity come primarily
from the improved noise floors of the Advanced LIGO interferometers over previous LIGO data sets.
Improvements in strain sensitivity over our previous O2 results
range from a factor of 1.2 at $\sim$100 Hz to $\sim$2 at $\sim$2000 Hz, with still larger
improvements at the lowest frequencies for which spectral line artifacts degraded O2 results.

No credible gravitational wave signals are observed, allowing upper limits to be placed on
possible source signal amplitudes.
Fig.~\ref{fig:upper_limits} shows the strain amplitude upper limits obtained.
The lowest upper limits on worst-case (linearly polarized) strain amplitude $h_0$ are  $\sim\sci{1.7}{-25}$ near 200\,Hz.
For a circularly polarized source (most favorable orientation), the smallest
upper limits are  $\sim\sci{6.3}{-26}$.
These upper limits refer to all sky locations and the entire range of frequency derivative values.
For a population-averaged ensemble of sky locations and stellar orientations,
the lowest 95\%\ CL upper limits on the strain amplitude are $\sim\sci{1.4}{-25}$.

At the highest frequencies ($\sim$2000 Hz) we are sensitive to neutron stars with an equatorial
ellipticity $\epsilon$ as small as $\sci{4}{-7}$ and as far away as 6 kpc for
favorable spin orientations (see Fig.~\ref{fig:range}).
For a higher ellipticity $\epsilon$ = $10^{-6}$ and favorable spin orientations,
we are sensitive to neutron stars beyond the galactic center.
The maximum ellipticity that a conventional neutron star can theoretically support
is at least $\sci{1}{-5}$ according to \cite{bib:crust_limit1,bib:crust_limit2}. Our results
are sensitive to such maximally deformed pulsars above $100$\,Hz pulsar rotation
frequency ($200$\,Hz gravitational-wave frequency) within $3$\,kpc.
Outliers from initial stages of each search method are followed up
systematically, but no candidates from any search survived scrutiny.

A recent, similar all-sky search of the LIGO O2 data~\cite{bib:O2Falcon1,bib:O2Falcon2,bib:O2Falcon3},
over the same frequency range (20-2000 Hz),
but severely restricted in spin-down range ($|\dot f|<\sci{3}{-12}$ Hz/s), achieved better strain sensitivity
(improvements ranging from $\sim$30\%\ lower at 500 Hz to $\sim$15\%\ lower at 1500 Hz) by using loose coherence
with an effective coherence time of 12 hours in its initial search stage, at a substantially higher computational cost
than that of this analysis. A similar approach~\cite{bib:deepsearch171-172Hz},
but further restricting the search frequency range to 171--172 Hz,
also achieved a higher sensitivity $h_0 \sim 1.07{\times}10^{-25}$ at 95\%\ confidence.
Figure~\ref{fig:range} illustrates the potential gain in ellipticity and
range sensitivity from searching for high-spin-down magnitudes. A recent deep search of the O2 data using
the distributed-computing project Einstein@Home~\cite{bib:O2EatH} over the frequency band 20-585.15 Hz and
spindown magnitude up to $\sci{2.6}{-9}$ Hz/s achieved similar sensitivity to this search over the common
search parameter space, quoting 90\%\ CL upper limits slightly lower than the 95\%\ CL values presented here.
Figure~\ref{fig:upper_limits} shows the strain upper limits achieved in this O3a search and in the previous O2 searches,
and Fig.~\ref{fig:parameterspace} shows a comparison of their parameter space coverages.

As the LIGO, Virgo and KAGRA gravitational wave detectors improve their strain sensitivities in the coming
decade~\cite{bib:scenariospaper}, searches will probe still smaller neutron star deformations and
explore further out into the galaxy.

\section{Acknowledgments}
This material is based upon work supported by NSF’s LIGO Laboratory which is a major facility
fully funded by the National Science Foundation.
The authors also gratefully acknowledge the support of
the Science and Technology Facilities Council (STFC) of the
United Kingdom, the Max-Planck-Society (MPS), and the State of
Niedersachsen/Germany for support of the construction of Advanced LIGO
and construction and operation of the GEO600 detector.
Additional support for Advanced LIGO was provided by the Australian Research Council.
The authors gratefully acknowledge the Italian Istituto Nazionale di Fisica Nucleare (INFN),
the French Centre National de la Recherche Scientifique (CNRS) and
the Netherlands Organization for Scientific Research,
for the construction and operation of the Virgo detector
and the creation and support  of the EGO consortium.
The authors also gratefully acknowledge research support from these agencies as well as by
the Council of Scientific and Industrial Research of India,
the Department of Science and Technology, India,
the Science \& Engineering Research Board (SERB), India,
the Ministry of Human Resource Development, India,
the Spanish Agencia Estatal de Investigaci\'on,
the Vicepresid\`encia i Conselleria d'Innovaci\'o, Recerca i Turisme and the Conselleria d'Educaci\'o i Universitat del Govern de les Illes Balears,
the Conselleria d'Innovaci\'o, Universitats, Ci\`encia i Societat Digital de la Generalitat Valenciana and
the CERCA Programme Generalitat de Catalunya, Spain,
the National Science Centre of Poland and the Foundation for Polish Science (FNP),
the Swiss National Science Foundation (SNSF),
the Russian Foundation for Basic Research,
the Russian Science Foundation,
the European Commission,
the European Regional Development Funds (ERDF),
the Royal Society,
the Scottish Funding Council,
the Scottish Universities Physics Alliance,
the Hungarian Scientific Research Fund (OTKA),
the French Lyon Institute of Origins (LIO),
the Belgian Fonds de la Recherche Scientifique (FRS-FNRS),
Actions de Recherche Concertées (ARC) and
Fonds Wetenschappelijk Onderzoek – Vlaanderen (FWO), Belgium,
the Paris \^{I}le-de-France Region,
the National Research, Development and Innovation Office Hungary (NKFIH),
the National Research Foundation of Korea,
the Natural Science and Engineering Research Council Canada,
Canadian Foundation for Innovation (CFI),
the Brazilian Ministry of Science, Technology, and Innovations,
the International Center for Theoretical Physics South American Institute for Fundamental Research (ICTP-SAIFR),
the Research Grants Council of Hong Kong,
the National Natural Science Foundation of China (NSFC),
the Leverhulme Trust,
the Research Corporation,
the Ministry of Science and Technology (MOST), Taiwan,
the United States Department of Energy,
and
the Kavli Foundation.
The authors gratefully acknowledge the support of the NSF, STFC, INFN and CNRS for provision of computational resources.
This work was supported by MEXT, JSPS Leading-edge Research Infrastructure Program, JSPS Grant-in-Aid for Specially Promoted Research 26000005, JSPS Grant-in-Aid for Scientific Research on Innovative Areas 2905: JP17H06358, JP17H06361 and JP17H06364, JSPS Core-to-Core Program A. Advanced Research Networks, JSPS Grant-in-Aid for Scientific Research (S) 17H06133, the joint research program of the Institute for Cosmic Ray Research, University of Tokyo, National Research Foundation (NRF) and Computing Infrastructure Project of KISTI-GSDC in Korea, Academia Sinica (AS), AS Grid Center (ASGC) and the Ministry of Science and Technology (MoST) in Taiwan under grants including AS-CDA-105-M06, Advanced Technology Center (ATC) of NAOJ, and Mechanical Engineering Center of KEK.

This document has been assigned LIGO Laboratory document number \texttt{LIGO-P2000334-v9}.

\newpage


\section{Introduction}

Superheavy elements with atomic numbers $104 \le Z \le 113$ were discovered in cold fusion reactions in which closed-shell $^{208}$Pb or $^{209}$Bi target nuclei were bombarded with projectiles ranging from Ti to Zn~\cite{Hofmann,Riken}.

The production cross section for superheavy nuclei (SHN) can be considered as the product of three factors: the cross section for the projectile to overcome the entrance channel barrier (capture cross section), the probability that the resulting system fuses and reaches the compound nucleus configuration, and the probability, that the excited compound nucleus survives fission during deexcitation.

The optimal conditions for obtaining a given superheavy nucleus result from various factors. The increase in the symmetry between reaction partners requires higher bombarding energies to overcome the entrance channel barrier and enhances the contribution of fast non-equilibrium deep-inelastic (DIC) and quasi fission (QF) processes preventing fusion. The fusion probability rapidly drops when the product of projectile and target nuclei atomic numbers $Z_{1}\times Z_{2} \geq 1600 $ ~\cite{Schmidt-91}. Therefore, the compound nucleus formation cross section represents only part of the capture cross section.

In cold fusion reactions, merging the strongly bound target and projectile nuclei leads to a weakly bound compound nucleus. Typically, SHN have higher thresholds for neutron emission~\cite{Jach2021} than the heights of the fission barrier~\cite{Jach2017}, making fission the dominant deexcitation process. At each step of the deexcitation cascade, neutron evaporation competes with fission, which additionally reduces the final evaporation residue cross section.

The cross sections for the production of nuclei with $Z\ge102$ in cold fusion reactions drop approximately seven orders of magnitude as the projectile atomic number changes from 20 (Ca) to 30 (Zn). The question what is the mechanism preventing the synthesis of SHN is still under discussion. A low survival probability is not enough to explain the extremely low production cross sections. One way of thinking about the additional hindrance mechanism is the concept of an internal barrier holding back (counteracting) the fusion process. Overcoming this barrier by a diffusion process and thermal fluctuations could help to reach the state of a compound nucleus.

Recently, the probability of compound nucleus formation $P_{CN}$, at energies around and above the interaction barrier $B_0$ was measured for $^{48}$Ca, $^{50}$Ti, and $^{54}$Cr projectile incident on a $^{208}$Pb target~\cite{HindePRL}. The authors of Ref.~\cite{HindePRL} concluded that ``the energy dependence of $P_{CN}$ indicates that cold fusion reactions (involving $^{208}$Pb) are not driven by a diffusion process''. This letter aims to investigate the fusion probabilities for these reactions using the diffusion approach. Calculations are performed within the $l$-dependent fusion-by-diffusion model (FbD)~\cite{FBD-11} using the new nuclear data tables for SHN by Jachimowicz et al.~\cite{Jach2021} as input.

\section{FbD model}~\label{SECT:THEORY}

The fusion-by-diffusion model in its first form was a simple tool to calculate cross sections and optimum bombarding energies for a class of $1n$ cold fusion reactions~\cite{FBD-Acta,FBD-05}. A significant development of this model was the incorporation of the angular momentum dependence, that is, the contributions from successive partial waves to the reaction cross section~\cite{FBD-11}.

Due to the different time scales of the particular reaction stages, the partial evaporation residue cross section, $\sigma_{ER}(l)$, can be factorized as the product of the partial capture cross section $\sigma_{cap}(l) = \pi {\mkern0.75mu\mathchar '26\mkern -9.75mu\lambda}^2(2l + 1)T(l)$, the fusion probability $P_{fus}(l)$, and the survival probability $P_{surv}(l)$. Thus, the total evaporation residue cross section for the production of a given superheavy nucleus in its ground state is
\begin{equation}
\label{factorize}
 \sigma_{ER} = \pi {\mkern0.75mu\mathchar '26\mkern -9.75mu\lambda}^2 \sum_{l = 0}^{\infty}(2l+1)T(l)\times P_{fus}(l)\times P_{surv}(l),
\end{equation}
where ${\mkern0.75mu\mathchar '26\mkern -9.75mu\lambda}$ is the wavelength, and ${\mkern0.75mu\mathchar '26\mkern -9.75mu\lambda}^2=\hbar^2/2\mu E_{c.m.}$. Here $\mu$ is the reduced mass of the colliding system, and $E_{c.m.}$ is the center-of-mass energy at which the reaction takes place.

The method of calculating the capture cross section is described in the next section. The fusion probability is described in detail in section~\ref{SOBSECT:FUS}.

The last factor in Eq.~\ref{factorize}, the survival probability, is calculated by applying classical transition state theory using nuclear data from Ref.~\cite{Jach2021}. Details regarding this reaction stage for $1n$ cold fusion reactions can be found in Ref. \cite{FBD-11}.

\subsection{Capture cross section}~\label{SOBSECT:CAP}
The capture transmission coefficients $T(l)$ in Eq. \ref{factorize} are calculated in a simple sharp cut off approximation, where the upper limit $l_{max}$ of full transmission, $T(l)=1$, is determined from the empirical systematics of the capture cross sections for heavy nuclear systems.

Following the experimental results, the entrance channel barrier is not described by a single value but by a distribution that can be approximated by a Gaussian shape described by two parameters, the mean barrier $B_{0}$ and the distribution width $\omega$~\cite{KSW04}. Folding the Gaussian barrier distribution with the classical expression for the fusion cross section leads to the formula for the capture cross section
\begin{eqnarray}
\label{capture}
\sigma_{cap} &=& \pi R^{2} \frac{\omega}{E_{c.m.} \sqrt {2\pi}} \Big[X\sqrt{\pi} (1 + \textrm{erf}(X))\nonumber\\&+&\exp(-X^{2})\Big] = \pi {\mkern0.75mu\mathchar '26\mkern -9.75mu\lambda}^2 (2l_{max}+1)^2,
\end{eqnarray}
where $X = \frac{E_{c.m.}-B_{0}}{\omega \sqrt{2}}.$
The empirical systematics of $B_0$, $\omega$, and the normalization factor $R$ were obtained from analyzing precisely measured fusion or capture excitation functions for about 50 heavy nuclear systems for which the fusion probability is equal or close to unity~\cite{KSW04}. In this paper we use the parametrizations of $B_0$, $w$, and $R$ of Ref.~\cite{FBD-11}.

\subsection{Fusion probability}~\label{SOBSECT:FUS}

The second factor in Eq.~\ref{factorize}, $P_{fus}(l)$, is the probability that after reaching the capture configuration, the colliding system will eventually overcome the fusion saddle point and merge, avoiding reseparation. It is assumed in the FbD model that after sticking, a neck between the target and projectile nuclei rapidly grows at an approximately fixed mass asymmetry and elongation \cite{FBD-Acta, FBD-05} bringing the system to the ``injection point'' somewhere along the bottom of the asymmetric fission valley. Let us denote the elengation of the system at the ``injection point'' by $L_{inj}$. The localization of this point with respect to the macroscopic conditional saddle (at the elongation $L_{sd}$) is crucial for the fusion process. If $L_{inj} > L_{sd}$ the system is still ``outside" the barrier separating the ``injection point'' from the compound nucleus configuration and must climb uphill to overcome the saddle. If $L_{inj} < L_{sd}$ the ``injection point'' configuration is more compact than the saddle configuration, and the system is already ``inside" (behind the barrier). In this case, the barrier guards the system against reseparation by reducing the outgoing flux of particles.

In the diffusion approach, transition over the barrier happens by thermal fluctuations in the shape degrees of freedom. The fusion probability, $P_{fus}(l)$, may be derived by solving the Smoluchowski diffusion equation. With the assumption that the internal barrier has height $H(l)$ and is of inverted parabola form one gets
\cite{FBD-Acta}
\begin{align}\label{pfus}
  P_{fus}(l) = \frac{1}{2} \left\{
\begin{array}{lr}
    1 + { \rm erf}\sqrt{ H(l)/T} & : L_{inj} < L_{sd} \\
    1 - { \rm erf}\sqrt{ H(l)/T} & : L_{inj} \ge L_{sd} \\
  \end{array},
\right.
\end{align}
where $T$ is the average temperature of the fusing system (see~\cite{FBD-11} for details).

The energy threshold in Eq.~\ref{pfus} is taken as the difference between the energy of the fusion saddle point $E_{sd}$ and the energy of the combined system at the ``injection point'' $E_{inj}$, corrected by the appropriate rotational energies,
\begin{equation}\label{H}
H(l) = (E_{sd} + E_{sd}^{rot}) - (E_{inj} + E_{inj}^{rot}).
\end{equation}

The energies $E_{sd}$ and $E_{inj}$ are calculated using simple algebraic expressions that approximate the potential energy surface \cite{FBD-11}. The shape parametrization used to describe the interacting system is that of two spheres joined smoothly by a third quadratic surface. The corresponding values of the rotational energies at the injection point $E_{inj}^{rot}$ and the saddle point $E_{sd}^{rot}$ are calculated assuming the rigid-body moments of inertia for the respective shapes~\cite{FBD-11}.

The distance between the nuclear surfaces of two colliding nuclei at the injection point, $s_{inj}$, is
the only adjustable parameter of the model. It defines the onset of the diffusion process, thus, the moment when the available kinetic energy that remains after passing the entrance barrier is already transformed into internal degrees of freedom in the over-damped regime.

\begin{figure}[t!!]
 \center{\includegraphics[width=0.8\columnwidth,angle=-90]{sinj_new.pdf}}
\caption{\label{Fig_sinj} The injection point systematics obtained for the set of $1n$ cold fusion reactions \cite{HOFMANN200493,Hofmann_1998,PhysRevC.78.024605,ISI:000172568800008,PhysRevC.78.034604,PhysRevC.79.027602,PhysRevC.78.024606,ISI:A1989U697600008,PhysRevLett.100.022501,PhysRevC.73.014611,PhysRevC.79.011602,ISI:A1997XX26900004,PhysRevC.79.027605,PhysRevC.67.064609,PhysRevLett.93.212702,ISI:000223717300009,MORITA2004101,ISI:000167664500002,doi:10.1143/JPSJ.73.1738,doi:10.1143/JPSJ.76.043201,doi:10.1143/JPSJ.73.2593,New113} using new nuclear data tables~\cite{Jach2021}. If not indicated otherwise, targets were $^{208}$Pb or $^{209}$Bi. The color of the points indicates the laboratory where the reaction was studied: LBNL (red), GSI (black), RIKEN (blue). See text for details. }
\end{figure}

In this paper, we redefine the systematics of this crucial parameter based on a new set of ground state and saddle point properties of SHN \cite{Jach2021}. The new parametrization of the ``injection point distance'' is shown in Fig.~\ref{Fig_sinj} as a function of the excess of the center-of-mass energy $E_{c.m.}$ over the mean barrier $B_{0}$. Each point represents the value of the $s_{inj}$ distance obtained by fitting Eq.~\ref{factorize} to the experimentally measured $1n$ evaporation residue cross sections for 27 cold fusion reactions (see Ref.~\cite{FBD-11} for fitting protocol details).

It can be seen from Fig.~\ref{Fig_sinj} that for energies up to a few MeV above $B_{0}$ the $s_{inj}$ distance can be well approximated by a straight line given by
\begin{equation}\label{sinj}
s_{inj} = 0.878~\textrm{fm} - 0.294\times(E_{c.m.} - B_{0})~\textrm{fm/MeV}.
\end{equation}
A similar linear trend of $s_{inj}$ as a function of $E_{c.m.} - B_{0}$ was reported in Ref.~\cite{Boilley} by solving Langevin type equations.

The shaded area in Fig.~\ref{Fig_sinj} represents an error corridor of $\pm 1$ fm, which allows the uncertainty of the calculated fusion probabilities to be determined. The parametrization given by Eq.~\ref{sinj} should be used for interpolation rather than extrapolation far beyond the explored range of $ E_{c.m.} - B_{0}$ values, especially if the extrapolation leads below the physically acceptable limit of the touching configuration $(s_{inj} \approx 0)$. Negative values of this parameter would correspond to a large overlap of the density distributions at the sticking stage, an effect that is impossible in nuclear collisions at low kinetic energies. Therefore, in collisions at energies higher than a few MeV above $B_{0}$, we assume $s_{inj}=0$ (allowing a deviation in the range of 1 fm).

\section{Results and discussion}~\label{SECT:RES}

\begin{figure}[t!!]
\center{\includegraphics[width=8cm,height=7cm,angle=0.0]{3drob.jpg}}
\caption{\label{3d} Fusion probability $P_{fus}(l)$ as a function of temperature $T$ and the barrier height  opposing fusion $H(l)$ for angular momenta $l = 0 \hbar, 20 \hbar, 40 \hbar$, and $60 \hbar$. Calculations for $^{48}$Ca +$^{208}$Pb (green lines), $^{50}$Ti +$^{208}$Pb (blue lines) and $^{54}$Ca +$^{208}$Pb (red lines) fusion reactions. The color of the surface marks the temperature gradient of the synthesized system. Labels ``outside" and ``inside" refer to the ``injection point'' position to the saddle (see Eq.~\ref{pfus} and its discussion).}
\end{figure}

The FbD model with the new ``injecton point distance'' parameterization and input data taken from Ref.~\cite{Jach2021} was used to analyze fusion probabilities for $^{48}$Ca, $^{50}$Ti, and $^{54}$Cr reactions incident on a $^{208}$Pb target.

The fusion probability given by Eq.~\ref{pfus} depends on the ratio of the height of the barrier opposing fusion $H(l)$ and the average temperature $T$ of the system during the merging process. The fusion probabilities, $P_{fus}(l)$, for $^{48}$Ca +$^{208}$Pb  (green lines),  $^{50}$Ti +$^{208}$Pb (blue lines), and $^{54}$Ca +$^{208}$Pb  (red lines) reactions as a functions of $H(l)$ and $T$ are shown in Fig.~\ref{3d} for a few selected $l$ values.

Let us start the discussion by analyzing the $l=0$ case in which the height of the barrier is simply the difference between the asymmetric conditional saddle point energy and the energy of the combined system of the projectile and target nuclei separated by the distance $s_{inj}$ (see Eq.~\ref{H}). As the available energy increases, the ``injection point distance'' decreases (see Fig.~\ref{Fig_sinj}), leading to a lowering of the barrier height, and thus to the rapid growth of the fusion probability (see  Fig. 2). When the separation distance reaches zero at the energy corresponding to $T \approx 0.6 -0.8$ MeV, the fusion probability reaches the maximum in all three analyzed reactions. Further energy increase does not change the height of the barrier ($s_{inj}$ remains equal to zero) but heats the system up and thus affects the fusion probability. For $^{48}$Ca +$^{208}$Pb reaction, the touching configuration is behind the asymmetric saddle point (as seen from the entrance channel perspective; ``inside" regime in Fig.~\ref{3d}). In this case, the rising temperature increases the flux of particles escaping through the asymmetric saddle point and thus slightly reduces the fusion probability (see $L_{inj} < L_{sd}$ case in Eq.~\ref{pfus}). For $^{50}$Ti and $^{54}$Cr projectiles, the touching configuration is ``outside" the barrier and the fusion probability slowly increases with the increase of the incindent energy ($L_{inj} \ge L_{sd}$ case in the Eq.~\ref{pfus}).

The inclusion of the higher partial waves affects the entire potential energy surface topology and influences the competition between the existing symmetric and asymmetric saddle points. In particular, the symmetric saddle, being more compact and having a lower moment of inertia, is more sensitive to the increase of the angular momenta. Above a certain $l$-value, the symmetric saddle begins to dominate and becomes the main point to overcome in the fusion process for all studied systems. In this case, the solution of the Smoluchowski diffusion equation can also be applied to calculate the fusion probability. However, the barrier height $H(l)$ should be calculated with respect to the symmetric saddle point. The barrier height increases with the increase of the $l$-value due to the difference in the rotational energies in the symmetric saddle and ``injection'' points. Therefore, the contribution of higher partial waves to the fusion cross section is suppressed. The systematic decrease of the fusion probability with the increase of the $l$-value observed in Fig.~\ref{3d} for all three reactions might be viewed as a manifestation of the well-known effect of the critical angular momentum.

As one can see in Fig.~\ref{3d}, the dominant contribution to the analyzed cold fusion reactions comes from near-central collisions. The more peripheral collisions are less favorable and lead to the re-separation of the system at the beginning of the nuclear reaction, rather than merging of target and projectile nuclei.
\begin{figure*}[t!!]
\center{\includegraphics[width=1.1\columnwidth,angle=-90]{multiplot_2.pdf}}
\caption{\label{multiplot} Capture cross section $\sigma_{cap}$ (top panels) and averaged  fusion probability $P_{fus}$ (bottom panels) for the reactions: $^{48}$Ca+$^{208}$Pb - panels (a), (b), $^{50}$Ti +$^{208}$Pb -  panels (c), (d), $^{54}$Cr +$^{208}$Pb -  panels (e), (f). Solid lines show the FbD model calculations of $\sigma_{cap}$ and $P_{fus}$. Dashed lines in the top panels show calculated $\sigma_{cap}$ scaled by the indicated suppression factors. The arrows indicate the value of the mean entrance channel barrier, $B_{0}$, for each reaction. The error corridors resulting from the $s_{inj}$ systematics uncertainty are shown as shaded areas in the bottom panels. Points represent relevant experimental data taken from Refs.~\cite{HindePRL,BOCK,PROKH,Pacheco,CLERC,Naik}. If not shown, error bars are smaller than the symbol sizes.}
\end{figure*}

To study the effective fusion probability for a given reaction, one can define the quantity
\begin{equation}
\label{fuseq}
 P_{fus} =\frac {1}{(2l_{max}+1)^2}\sum_{l = 0}^{l_{max}}(2l+1)\times P_{fus}(l),
\end{equation}
which is the fusion probability ``averaged'' over all angular momenta contributing to the fusion cross section.

In Fig.~\ref{multiplot} we present a comparison of the FbD model predictions with the experimental data. The top panels show the capture cross sections for each of the reactions (i.e. cross sections for overcoming the entrance channel barrier) calculated using Eq.~\ref{capture}. Model calculations are compared with experimental data taken from Ref.~\cite{HindePRL,BOCK,PROKH,Pacheco,CLERC,Naik}. The arrows in panels (a), (c), and (e) indicate the values of the mean entrance channel barriers $B_0$, calculated using the empirical parametrization~\cite{FBD-11}, $173.0$, $191.2$, and $208.3$ MeV for $^{48}$Ca, $^{50}$Ti, and $^{54}$Cr respectively.

The experimentally measured fission-like cross sections shown in the top panels of Fig.~\ref{multiplot} lie below our calculations (solid lines). The deviation increases with increasing projectile atomic number. As proposed in \cite{HindePRL}, we estimated scaling factors, $S$, for our calculations to reproduce the experimental results in the energy range above $B_0$. These factors are 0.96 for $^{48}$Ca+$^{208}$Pb, 0.55 for $^{50}$Ti+$^{208}$Pb and 0.24 for $^{54}$Cr+$^{208}$Pb (in \cite{HindePRL} the respective factors are 0.75, 0.48 and 0.22). Scaled capture cross sections are shown as dashed lines in Fig.~\ref{multiplot}.

Our scaling factors are in reasonable agreement with the results presented in Ref.~\cite{HindePRL}, where they were estimated as a deviation from the CCFULL model based on the coupled channels formalism \cite{Hagino}. The capture cross section suppression might be associated with mass-asymmetric fast non-equilibrium processes, such as QF or DIC, appearing just after the interacting system passes the entrance channel barrier. It should be emphasized that both in this work and Ref.~\cite{HindePRL}, the obtained scaling factors are model-dependent.

Calculated average fusion probabilities (see Eq.~\ref{fuseq}) for $^{48}$Ca, $^{50}$Ti, and $^{54}$Cr reactions on a $^{208}$Pb target are shown in the lower panels of Fig.~\ref{multiplot}. Full points in panels (b), (d), and (f) represent upper limits on the compound nucleus formation probabilities $P_{sym}$ taken from Ref.~\cite{HindePRL}. $P_{sym}$ is derived as the ratio of the measured symmetric-peaked fission cross section $\sigma_{sym}$ to the capture cross section taken as the measured total fission-like cross section $\sigma_{fis}$ divided by the appropriate scaling factor $S$ ($P_{sym} = \frac{\sigma_{sym}}{\sigma_{fis}/S}$, see Eq. 1 in the supplementary material for~\cite{HindePRL}).

For the $^{50}$Ti+$^{208}$Pb reaction additional experimental points (open triangles in Fig.~\ref{multiplot}(d)) taken from \cite{Naik} are shown. These data were derived by measuring the angular distribution of mass-symmetric fission.

The calculated average fusion probabilities (Eq.~\ref{fuseq}, solid lines in panels (b),(d), and (f)) are in good agreement with the experimental data for all studied reactions. A rapid decrease of the fusion probability in the energy region below $B_0$ reported in~\cite{Naik} for the $^{50}$Ti+$^{208}$Pb reaction is reproduced in our calculations (see Fig.~\ref{multiplot}(d)). Unfortunately, the data for other reactions are limited in this energy region.

For each reaction, the maximum value of $P_{fus}$ is reached for an energy a few MeV above $B_{0}$ (when $s_{inj} \approx 0$). Thus, the steady decrease of $P_{fus}$ with increasing energy is due to the dependence of $H(l)$ on the angular momentum only (see Fig.~\ref{3d}).

Finally, in Fig.~\ref{PCN}, we show the FbD model calculations of compound nucleus formation cross sections, defined as
\begin{equation}
\label{PCNeq}
 \sigma_{fus} = \pi {\mkern0.75mu\mathchar '26\mkern -9.75mu\lambda}^2 \sum_{l = 0}^{l_{max}}(2l+1)
 T(l)P_{fus}(l) = \sigma_{cap}\times P_{fus}.
\end{equation}
The model calculations (solid lines) are compared with the symmetric-peaked fission cross sections measured in~\cite{HindePRL}, but not given in that paper. Therefore, we deduced $\sigma_{sym}$ values from the data using the relation $\sigma_{sym} = P_{sym} \times (\sigma_{fis}/S)$ (see Fig. 4 in Ref.~\cite{HindePRL} and Fig. 5 in the corresponding supplementary material). Although the comparison of $\sigma_{sym}$ with $\sigma_{fus}$ is not entirely unequivocal, it seems adequate. Qualitatively, one can see very similar behavior of the calculated and experimental cross sections as a function of the center-of-mass energy.
Some additional subtle effects related to the deformation of the $^{54}$Cr projectile might be expected. However, as we have checked (by analyzing collisions at the tip to tip and equatorial configurations; dashed lines in Fig.~\ref{PCN}), the results for extreme orientations in the entrance channel are within the error corridor resulting from the systematics of $s_{inj}$.
\begin{figure}[t!!]
 \center{\includegraphics[width=1.0\columnwidth,angle=-90]{fis_all2.pdf}}
\caption{\label{PCN} Calculated compound nucleus formation cross sections, $\sigma_{fus}$, for $^{48}$Ca +$^{208}$Pb, $^{50}$Ti +$^{208}$Pb, and $^{54}$Ca +$^{208}$Pb  fusion reactions. Points are derived from the experimental data presented in Ref.~\cite{HindePRL}. Dashed lines show calculations for two extreme orientations of target and $^{54}$Cr projectile in the entrance channel: tip to tip (bottom line) and body to body (upper line). The arrows indicate the value of the mean entrance channel barrier, $B_{0}$, for each reaction.
See text for details.}
\end{figure}

\section{Conclusions}~\label{SECT:CONCS}

The presented results show that the compound nucleus formation cross sections and related fusion probabilities for $^{48}$Ca, $^{50}$Ti, and  $^{54}$Cr incident on a $^{208}$Pb target can be well reproduced within the framework of the FbD model.

The experimentally observed dependence of the fusion probability on the energy can be reproduced using the diffusion approach. In the energy range below $B_0$, the fusion probability growth comes from the reduction in the height of the internal barrier opposing fusion with increasing bombarding energy. The fusion probability saturation above $B_0$ results from suppression of the contributions from higher partial waves and can be linked to the critical angular momentum. The difference between rotational energies in the fusion saddle and the contact (sticking) configuration at the beginning of the fusion process plays a major role in compound nucleus formation at energies above $B_0$.

\section*{ACKNOWLEDGEMENTS}
M.K. was co-financed by the National Science Centre under Contract No. UMO-2013/08/M/ST2/00257  (LEA COPIGAL).