\section{Introduction}
\label{sec: intro}
On October 19, 2017, the first interstellar object in the Solar System, `Oumuamua (1I/2017 U1)  was discovered by the PAN-STARRS1 survey.
It has a highly hyperbolic trajectory (with eccentricity $e=1.1956 \pm 0.0006$) and pre-entry velocity of $v_{\infty}\approx 26$ km s$^{-1}$ \citep{Meech2017}.
Based on the survey properties and the single detection, \citet{do2018} estimated the interstellar density of objects like `Oumuamua or larger to be $n \approx 2 \times 10^{15}$ pc$^{-3}$, 2-8
orders of magnitude larger than expected by previous theoretical models \citep{Moro-Martin2009}.
The large variations in its apparent magnitude and the non-trivial periodicity of the lightcurve, suggest that `Oumuamua is rotating in an excited spin state (tumbling motion), and has an extreme aspect ratio of at least $5:1$   \citep{Fraser2018, Drahus2018}, an unprecedented value for previously known asteroids and comets in the Solar System.
\citet{Belton2018} have shown that if `Oumuamua rotates in its highest rotational energy state, it should be extremely oblate (pancake-like).

Recently, \citet{Micheli2018} reported the detection of non-gravitational acceleration in the motion of `Oumuamua, at a statistical significance of 30$\sigma$. Their best-fit to the data is obtained for a model with a non-constant excess acceleration which scales with distance from the Sun, $r$, as $\Delta a \propto r^{-2}$, but other power-law index values are also possible.
They concluded that the observed acceleration is most likely the result of a cometary activity.
Yet, despite its close Solar approach of $r = 0.25$ AU, `Oumuamua shows no signs of a any cometary activity, no cometary tail, nor gas emission/absorption lines were observed \citep{Meech2017, Knight2017, Jewitt2017, Ye2017, Fitzsimmons2017}.
From  a theoretical point of view, \citet{Rafikov2018} has shown that if outgassing was responsible for the acceleration \citep[as originally proposed by][]{Micheli2018}, then the associated outgassing torques would have driven a rapid evolution in `Oumuamua's spin, incompatible with observations.

If not cometary activity, what can drive the non-gravitational acceleration observed?
In this {\it Letter} we explore the possibility of `Oumuamua being a thin object accelerated by Solar radiation pressure, which would naturally result in an excess acceleration $\Delta a\propto r^{-2}$.
\footnote{Interestingly, a similar approach was adopted historically regarding the anomalous orbit of Phobos \citep{Shklovskii1962}.}
However, for radiation pressure to be effective, the mass-to-area ratio must be very small.
In \S \ref{sec: rad pressure} we derive the required mass-to-area ratio and find $(m/A)\approx 0.1$ g cm$^{-2}$, corresponding to an effective thin sheet of thickness $w \approx 0.3-0.9$ mm.
We explore the ability of such an unusually thin object to survive interstellar travel, considering collisions with interstellar dust and gas (\S \ref{sec: travel}), as well as to withstand the tensile stresses caused by rotation and tidal forces (\S \ref{sec: strength}).
Finally, in \S \ref{sec: discussion} we discuss the possible implications of the unusual requirements on the shape of `Oumuamua.

\section{Acceleration by Radiation Pressure}
\label{sec: rad pressure}
\citet{Micheli2018} had shown that `Oumuamua
experiences an excess radial acceleration, with
their best fit model
\begin{align}
\label{eq: a_obs}
    \Delta a &= a_0 \left( \frac{r}{\rm AU} \right)^{n} \\ \nonumber
    {\rm with} \ \ n &= -2 \ , \\
     a_0 & = (4.92 \pm 0.16) \times 10^{-4} \ {\rm cm \ s^{-2}} \ .
\end{align}
The value for $a_0$ is averaged over timescales much longer than `Oumuamua's rotation period.

An acceleration of this form is naturally produced by radiation pressure,
\begin{equation}
\label{eq: rad_press}
P = C_R \frac{L_{\odot}}{4 \pi r^2 c} \ ,
\end{equation}
where $L_{\odot}$ is the Solar luminosity, $c$ is the speed of  light, and $C_R$ is a coefficient of order unity which depends on the object's composition and geometry.
For a sheet perpendicular to the Sun-object vector $C_R= 1+\epsilon$ where $\epsilon$ is the reflectivity. For a perfect reflector $\epsilon=1$, and $C_R=2$, whereas for a perfect absorber $\epsilon = 0$ and $C_R = 1$.

For an object of mass $m$ and area $A$
the acceleration would be
\begin{align}
\label{eq: a_rad_press}
a &=  \frac{PA}{m}
= \left( \frac{L_{\odot}}{4 \pi r^2 c} \right)
 \left(\frac{A}{m}  \right)
  C_R \\ \nonumber
& = 4.6 \times 10^{-5}
\left( \frac{r}{\rm AU}     \right)^{-2}
\left( \frac{m/A}{\rm g \ cm^{-2} } \right)^{-1}
C_R \
{\rm cm \ s^{-2} }
\ .
\end{align}
Comparing Eq.~(\ref{eq: a_obs}) and (\ref{eq: a_rad_press}) we find that the requirement on the mass-to-area ratio is
\begin{equation}
\label{eq: constraint  m_A}
     \left( \frac{m}{A} \right)

    = (9.3 \pm 0.3) \times 10^{-2} \ C_R \ {\rm g \ cm^{-2}} \ .
\end{equation}
For a planar body with mass density $\rho$, this translates into a requirement on the body's thickness
\begin{equation}
w \ = \ \frac{m}{A \rho} = (9.3 \pm 0.3) \times 10^{-2} \rho_0^{-1} \ C_R \ {\rm cm} \ ,
\end{equation}
where $\rho_0=\rho/(10^0 \ {\rm g \ cm^{-3}})$. Typically, $\rho_0 \approx 1-3$, giving a thin sheet of 0.3--0.9 mm thickness.
Other geometries are also possible, and are discussed in \S \ref{sec: discussion}.
The force exerted by the Solar wind on a solid surface is negligible compared to that of the Solar radiation field and is neglected hereafter.

The observed magnitude of `Oumuamua constrains its area to be
$A \approx 8 \times 10^6 \alpha^{-1}$ cm$^2$, where $\alpha$ is the albedo \citep{Jewitt2017}.
This corresponds to an effective radius, $R_{\rm eff} \equiv \sqrt{A}/\pi = 16 \alpha^{-1/2}$ meters.
For our estimation of the mass-to-area ratio, this area translates into a mass of $m \approx 740  (C_R/\alpha)$ kg.

\section{Maximum  Distance for Interstellar Travel}
\label{sec: travel}
\begin{figure}[t]
	\centering
	\includegraphics[width=0.5\textwidth]{L_max}
	\caption{
	 The maximum allowed travel distance through the interstellar medium (ISM), as a function of $(m/A)$. The blue and red lines are limitations obtained by slow-down due to gas accumulation, and vaporization by dust-collisions, respectively.
The plotted results are for a mean ISM proton density of $\langle n \rangle \sim 1$ cm$^{-3}$. All lines scale as $1/\langle n \rangle$.
The dashed magenta line is our constraint on `Oumuamua based on its excess acceleration.
The Solar Galactrocentric distance is also indicated.
		}
		\label{fig: L_max}
\end{figure}

Next we explore the implications of impacts with interstellar dust-grains and gas particles, in terms of momentum and energy transfer.
We obtain general requirements for the object's mass-to-area ratio, or alternatively, for the maximum interstellar distance that can be traveled before it encounters appreciable slow-down or evaporation.

\subsection{Momentum Transfer - Slow Down}
An object with a cross sectional area $A$ traveling a distance $L$ through the ISM would accumulate an ISM gas mass of,
\begin{equation}
\label{eq: M ISM gas swept}
M_{\rm gas} = A \Sigma_{\rm gas} = 1.4 m_{p} \langle n \rangle L A \ ,
\end{equation}
where $\Sigma_{\rm gas}$ is the accumulated mass column density of interstellar gas,  $m_p$ is the proton mass, $\langle n \rangle$ is the mean proton number density averaged along the object's trajectory, and the factor 1.4 accounts for the contribution of helium to the mass density of the ISM.
For trajectories that span Galactic distances the contribution of the Solar System to the accumulated column is negligible.
The contribution of accumulated dust to the momentum transfer is also negligible, because the typical dust-to-gas mass ratio in the Galaxy is $\approx 1/100$.

Once $M_{\rm ISM}$ approaches the object's mass, $m$, the momentum of the traveling object will decrease by a significant amount.
The requirement $M_{\rm ISM}/m \ll 1$ translates into a maximum allowed value on the object's mass-to-area ratio, giving
\begin{align}
\label{eq: m_A momentum}
\left( \frac{m}{A} \right)_{\rm min, p} &= \Sigma_{\rm gas} = 1.4 m_p \langle n \rangle L \\ \nonumber
& = 7.2 \times 10^{-3} \langle n \rangle_0 L_{0} \ {\rm g \ cm^{-2}} \ .
\end{align}
In the last equality we normalized to the typical values $\langle n \rangle_0 = \langle n \rangle/(10^0 \ {\rm cm^{-3}})$ and  $L_0 = L/(10^0 \ {\rm kpc})$.

Given a mass-to-area ratio, the maximum travel distance is
\begin{align}
\label{eq: L_max p}
L_{\rm max, p} &=  \frac{1}{1.4 m_p \langle n \rangle} \left( \frac{m}{A} \right)\\ \nonumber
& = 14 \langle n \rangle_0^{-1}
\left( \frac{m}{A} \right)_{-1}
 \ {\rm kpc} \ .
\end{align}
In the second equation we denoted
$(m/A)_{-1} \equiv (m/A)/(10^{-1} {\rm \ g \ cm^{-2}})$.
Figure~\ref{fig: L_max} shows the results from Eq.~(\ref{eq: L_max p}) as a function of ($m/A$).
The dashed vertical line indicates our constraint on ($m/A$) for `Oumuamua (Eq.~\ref{eq: constraint  m_A}).
Evidently, `Oumuamua can travel Galactic distances before encountering appreciable slow-down.

\subsection{Energy Transfer - Collisions with Dust-Grains}

Collisions with dust grains at high velocities will induce crater formation by melting and evaporation of the target material.
Since the typical time between dust collisions is long compared to the solidification time, any molten material will solidify before the next collision occurs, and thus will only cause a deformation of the object's surface material, not reduction in mass.
On the other hand, atoms vaporized through collisions can escape and thus cause a mass ablation.

We would like to estimate the minimum mass-to-area ratio required for the object to not lose significant fraction of its mass upon dust-grain collisions.
Let $m_d$ be the colliding dust-grain mass, and $\phi$ the fraction of the kinetic energy that is converted into vaporization of the object's body. The total number of vaporized atoms per collision is then,
\begin{equation}
\label{eq: dust melting one grain}
    N_v =  \frac{m_v}{\bar{m}} = \phi \frac{ m_d v^2}{2U_v}\ ,
\end{equation}
where $m_v$ is the mass of vaporized material in the object, $\bar{m}$ is the mean atomic mass of the object, and $U_v$ is the vaporization energy.
Although highly simplistic, this analysis captures the results of the detailed  theoretical model of \citet{Tielens1994} and the  empirical data from   \citet{Okeefe1977}.
A good match to the numerical results is obtained for $\phi = 0.2$.

Over a distance $L$ there will be many collisions.
Adopting the conservative assumption that for each collision all the vaporized material escapes to the ISM, we can account for all the collisions along a path-length, $L$, by replacing $m_d$ in Eq.~(\ref{eq: dust melting one grain}) with the total accumulated dust mass,
\begin{equation}
\label{eq: total dust}
    M_{d} = \Sigma_{\rm dust} A = \Sigma_{\rm gas} \varphi_{dg} A = 1.4 m_p \langle n \rangle L \varphi_{dg} A \ ,
\end{equation}
where $\Sigma_{\rm dust}$ is the dust column density and $\phi_{dg}$ is the dust-to-gas mass ratio. This gives,
\begin{equation}
\label{eq: dust melting all dust}
    \frac{m_v}{\bar{m}} = \phi \frac{ 1.4 m_p \langle n \rangle L \varphi_{dg} A v^2}{2U_v} \ .
\end{equation}
Requiring that the total vaporized mass not exceed half of the object's mass, we obtain a constraint on the minimum mass-to-area ratio of the object,
\begin{align}
\label{eq: m_A dust melt}
\left( \frac{m}{A} \right)_{\rm min, m}
&= 1.4 m_p \langle n \rangle L
\left( \frac{\phi  \varphi_{dg} \bar{m} v^2}{U_v} \right)  \\ \nonumber
& = 3.1 \times 10^{-4} \langle n \rangle_0 L_{0}
\left( \frac{\varphi_{-2}\bar{m}_{12} v^2_{26}}{U_{4}} \right) \ {\rm g \ cm^{-2}} \ .
\end{align}
Here we defined the normalized parameters, $\varphi_{-2}=\varphi_{dg}/10^{-2}$, $\bar{m}_{12}=\bar{m}/(12 m_p)$, as appropriate for carbon based materials (e.g., graphite or diamond); $U_{4}=U/(4 \ {\rm eV})$, as appropriate for typical vaporization energies (e.g., for graphite, $U_v=4.2$ eV); and $v_{26}=v/(26 \ {\rm km \ s^{-1}})$, the velocity at infinity of `Oumuamua.
For a given mass-to-area ratio, the maximum allowed distance before significant evaporation is
\begin{align}
\label{eq: L_max melt}
    L_{\rm max, d} &= \frac{1}{1.4 m_p \langle n \rangle} \left( \frac{U_m}{\phi  \varphi_{dg} \bar{m} v^2} \right) \left( \frac{m}{A} \right) \\ \nonumber
    & = 330 \langle n \rangle_0^{-1} \left( \frac{U_{4}}{\varphi_{-2}\bar{m}_{12} v^2_{26}} \right) \left(
\frac{m}{A}
\right)_{-1} \ {\rm kpc} \ .
\end{align}
For our constrained value for the mass-to-area ratio, `Oumuamua can travel through the entire galaxy before a significant fraction of its mass is evaporated.
Evaporation becomes important  at higher speeds.
Comparing Eqs.~(\ref{eq: m_A momentum}) and (\ref{eq: m_A dust melt}) we find that only for speeds above
\begin{align}
    v_{\rm crit} &= \sqrt{\frac{2U_m}{\bar{m}\varphi \phi}} \\ \nonumber
    & = 130 \left( \frac{U_{4}}{\varphi_{-2}\bar{m}_{12} } \right)^{1/2} \ {\rm km \ s^{-1}} \ ,
\end{align}
vaporization dominates over slow-down.

\subsection{Energy Transfer - Collisions with Gas Particles}

When an object travels at a high speed, collisions with atoms in the ISM can potentially transfer sufficient energy to produce sputtering.
This process was studied in the context of dust grains  in hot shocks \citep{Tielens1994}.
For an object traveling at a velocity $v$, over a distance $L$ through the ISM, the total number of sputtered particles is,
\begin{equation}
    N_s  = \frac{m_s}{\bar{m}} =  2 A  L \langle n \rangle Y_{\rm tot} \ ,
\end{equation}
where $Y_{\rm tot}=\sum_i Y_i x_i$ is the total sputtering yield (which depends on the kinetic energy), summed over collisions with different species (i.e., H, He, and metals), and $x_i$ is the abundance of the colliding species  relative to hydrogen.

The minimal mass-to-area ratio below which half of the object's mass will be sputtered is,
\begin{align}
\label{eq: m_A_sputtering}
\left( \frac{m}{A} \right)_{\rm min, s} &= \bar{m} \langle n \rangle L  Y_{\rm tot}  \\ \nonumber
&= 6.2 \times 10^{-5} \bar{m}_{12}  \langle n \rangle_0 L_0 Y_{-3} \ {\rm g \ cm^{-2}} \ .
\end{align}
In the second equality we normalized to $Y_{\rm tot}=10^{-3}$, corresponding to kinetic energies $E \approx 30 - 100$ eV,  corresponding to $v \approx 40 - 70$ km s$^{-1}$.
For lower speeds, as that of `Oumuamua, the yield is even lower, further decreasing the value of $(m/A)_{\rm min, s}$.
At higher speeds, the yield increases but typically remains below $0.01$ \citep{Tielens1994}, thus at any velocity, vaporization and slow-down remain the dominating processes limiting the allowed distance an object can travel through the ISM.

Cosmic-rays are expected to cause even less damage.
Although their energy density is comparable to that of the ISM gas, they deposit only a very small fraction of their energy as they penetrate through the thin object.

\section{Tensile Stresses}
\label{sec: strength}

A thin object can be torn apart by
centrifugal forces or tidal forces
if its tensile strength is not sufficiently strong.
Typical values for the tensile strengths of various materials are shown in Table \ref{table: tensile strength}.
Next, we calculate whether centrifugal or tidal forces can destroy `Oumuamua.

\begin{table}[]
\caption{Tensile Strengths}
\centering
\begin{tabular}{l l}
\hline\hline
material  & tensile strength (${\rm dyne \  cm^{-2}}$)
  \\ [0.5ex]
  \hline
  67P/Churyumov-Gerasimenko$^1$ & 10-50
\\
Meteorites$^2$ & $(1 - 5) \times 10^7$
\\
Iron & $3 \times 10^{7}$ \\
  Diamond & $2 \times 10^{10}$ \\
  Silicon (monocrystalline) & $7 \times 10^{10}$ \\
  \hline
    \hline
$^1$ \citet{Attree2018} & \\
$^2$ \citet{Petrovic2001} & \\
\end{tabular}
\label{table: tensile strength}
\end{table}

\subsection{Rotation}
\label{sub: rot}

Oumuamua's lightcurve shows periodic modulations on an order of 6-8 hours. Ignoring the tumbling motion, let us estimate the tensile stress originating from the centrifugal force.
The largest stress is produced if the object is elongated such that the longest dimension is perpendicular to the rotation axis. We denote this dimension as $d$.
Considering the object as made of two halves, each located with a center of mass at a distance $d/4$ from the rotation axis, and ignoring self gravity, a radial force of magnitude,
\begin{equation}
F = \frac{1}{2} m \Omega^2 \frac{d}{4} \ ,
\end{equation}
will be exerted on each half.
The associated tensile stress is,
\begin{align}
P_{\rm rot} &= \frac{1}{4} \rho d^2  \Omega^2\\ \nonumber \nonumber
&= 0.25 \ \rho_0 \ d_{4}^2 \ \Omega_{-4}^2 \  {\rm dyne \ cm^{-2}} \ ,
\end{align}
where $d_4 \equiv d/(10^4 \ {\rm cm})$, $\Omega_{-4} \equiv \Omega/(10^{-4} \ {\rm s^{-1}})$.
This is much smaller than typical tensile strengths of normal materials, and even of that of the comet 67P/Churyumov-Gerasimenko (see Table \ref{table: tensile strength}).
Thus, even when self-gravity is ignored, `Oumuamua can easily withstand its centrifugal force.

\subsection{Tidal Forces}
The tidal force will be maximal if the long dimension of the object is parallel to the Sun-object vector.
Again, modeling the object as consisting of two halves as in \S \ref{sub: rot},
the difference in the gravitational force experienced by the far and near ends of the object is,
\begin{equation}
\delta F \approx \frac{1}{4} d \frac{GM_{\odot}m}{r^3} \ ,
\end{equation}
where $r$ is the distance of the center of mass from the Sun.
The associated tensile stress,
\begin{align}
P_{\rm tid} &\approx \frac{1}{4} \rho d^2 \frac{GM_{\odot}}{r^3} \\ \nonumber
& = 9.9 \times 10^{-7} \ \rho_0 \ d_4^2 \left(\frac{r}{{\rm AU}} \right)^{-3} \ {\rm dyne \ cm^{-2}} \ .
\end{align}
Even at perihelion ($r=0.25$ AU), the tensile stress is negligible.

The critical distance below which tidal forces dominate over centrifugal is
\begin{align}
R_{\rm tid} = \left(\frac{G M}{\Omega^2} \right)^{1/3}  = 3.4  \Omega_{-4}^{-2/3} \left(  \frac{M}{M_{\odot}}\right)^{1/3} \ R_{\odot} \ .
\end{align}
Thus, unless `Oumuamua encountered an extremely close approach to a star in its past, it is unlikely that tidal forces played any significant role.

\section{Summary and Discussion}
\label{sec: discussion}
We have shown that the observed non-gravitational acceleration of `Oumuamua, may be explained by Solar radiation pressure.
This requires a small mass-to-area ratio for `Oumuamua of $(m/A) \approx 0.1$ g cm$^{-2}$.
For a planar geometry and typical mass densities of 1--3 g cm$^{-2}$  this gives an effective thickness of only 0.9--0.3 mm, respectively. For a material with lower mass density, the inferred effective thickness is proportionally larger.
We find that although very thin, such an object can travel over galactic distances, maintaining its momentum and withstanding collisional destruction by dust-grains and gas,
as well as centrifugal and tidal forces.
For `Oumuamua, the limiting factor is the slow-down by accumulated ISM mass, which limits its maximal travel distance to $\sim 10$  kpc (for a mean ISM particle density of $\sim 1$ cm$^{-3}$).
Our inferred thin geometry is consistent with studies of its tumbling motion. In particular, \citet{Belton2018} inferred that `Oumuamua is likely to be an extremely oblate spheroid (pancake) assuming that it is  excited by external torques to its highest energy state.

While our scenario may naturally explains the peculiar acceleration of `Oumuamua,
it opens up the question what kind of object might have such a small mass-to-area ratio?
The observations are not sufficiently sensitive to provide a resolved image of `Oumuamua, and one can only speculate on its possible geometry and nature.
Although periodic variations in the apparent magnitude are observed, there are still too many degrees of freedom (e.g., observing angle, non-uniform reflectively, etc.) to definitely constrain the geometry.
The geometry should not necessarily be that of a planar sheet, but may acquire other shapes, e.g., involving a curved sheet, a hollow cone or ellipsoidal, etc.
Depending on the geometry our estimated value for the mass-to-area ratio will change (through $C_R$ in Eq.~\ref{eq: constraint  m_A}), but the correction is typically of order unity.

Known Solar System objects, like asteroids and comets have mass-to-area ratios orders of magnitude larger than our estimate for `Oumuamua.
If radiation pressure is the accelerating force, then `Oumuamua represents a new class of thin interstellar material, either  produced naturally,through a yet unknown process in
the ISM or in proto-planetary disks, or of an artificial origin.

Considering an artificial origin, one possibility is that `Oumuamua is a lightsail,   floating in interstellar space
as a debris from an advanced technological equipment  \citep{Loeb2018}.
Lightsails with similar dimensions have been designed and constructed by our own civilization, including the IKAROS project and the Starshot Initiative\footnote{A list of books and papers on lightsails is provided in \href{http://breakthroughinitiatives.org/research/3}{http://breakthroughinitiatives.org/research/3}. The IKAROS project is discussed in  \href{http://global.jaxa.jp/projects/sat/ikaros/}{http://global.jaxa.jp/projects/sat/ikaros/}.}.
The lightsail technology might be abundantly used for transportation of cargo between planets \citep{Guillochon2015} or between stars \citep{Lingam2017}. In the former case, dynamical ejection from a planetary System could result in space debris of equipment that is not operational any more
\footnote{Note that `Oumuamua was found not to show any radio emission down to a fraction of the power of a cell phone transmission \citep{Enriquez2018a, Tingay2018, Harp2018}.}
\citep{Loeb2018}, and is floating at the characteristic speed of stars relative to each other in the Solar neighborhood.
This would account for the various anomalies of `Oumuamua, such as the
unusual geometry
 inferred from its lightcurve \citep{Meech2017, Fraser2018, Drahus2018, Belton2018},
its low thermal emission, suggesting high reflectivity \citep{Trilling2018},
and its deviation from a Keplerian orbit \citep{Micheli2018} without any sign of a cometary tail \citep{Meech2017, Knight2017, Jewitt2017, Ye2017, Fitzsimmons2017} or spin-up torques \citep{Rafikov2018}.
 Although `Oumuamua has a red surface color, similar to organic-rich surfaces of Solar-System comets and D-type asteroids \citep{Meech2017}, this does not contradict the artificial scenario, since irrespective of the object's composition, as it travels through the ISM its surface will be covered by a layer of interstellar dust, which is itself composed of organic-rich materials \citep{Draine2003}.

Alternatively, a more exotic scenario is that
`Oumuamua may be a fully operational probe sent {\it intentionally} to Earth vicinity by an alien civilization.
Based on the PAN-STARRS1 survey characteristics, and assuming natural origins following {\it random} trajectories,
  \citet[][]{do2018} derived that the interstellar number density of `Oumuamua-like objects should be extremely high, $\sim 2 \times 10^{15}$ pc$^{-3}$, equivalent to $\sim 10^{15}$  ejected planetisimals  per star, and a factor of 100 to 10$^8$ larger than predicted by
  theoretical models \citep{Moro-Martin2009}.
This discrepancy is readily solved if `Oumuamua does not follow a random trajectory but is rather a targeted probe.
Interestingly, `Oumuamua's entry velocity is found to be extremely close to the velocity of the Local Standard of Rest, in a kinematic region that is occupied by less than 1 to 500 stars \citep{Mamajek2017}.

Since it is too late to image `Oumuamua with existing telescopes or chase it with chemical propulsion rockets \citep[][but see \citealt{Hein2017}]{Seligman2018}, its likely origin and mechanical properties could only be deciphered by searching for other objects of its type in the future. In addition to the vast unbound population,
thousands of interstellar `Oumuamua-like space-debris
 are expected to be trapped at any given time in the Solar System through gravitational interaction with Jupiter and the Sun \citep{Lingam2018}.
Deep wide-area surveys of the type expected with the Large Synoptic Survey Telescope (LSST)\footnote{See  https://www.lsst.org/} will be particularly powerful in searching for additional members of `Oumuamua's population of objects.

A survey for lightsails as technosignatures in the Solar System is warranted, irrespective of whether `Oumuamua is one of them.

\acknowledgements
We thank Manasvi Lingam, Paul Duffell, Quanzhi Ye, and an anonymous referee for helpful comments.
This work was supported in part by a grant from the Breakthrough Prize Foundation.


\section{Introduction}
Multi-messenger astrophysics produced two foundational discoveries in 2017: the detection of a binary-neutron star merger through gravitational waves (GWs) and electromagnetic emission \cite{2017ApJ...848L..12A}, and the observation of a blazar through high-energy neutrinos and electromagnetic emission \cite{ic1709022mm}. The multimessenger science reach of the GW detectors had been enabled by decades of effort preceding the discovery~\cite{LIGOG060660,2006ivoa.spec.1101S,2008CQGra..25k4039A,2008CQGra..25k4051A,2009IJMPD..18.1655V,2010JPhCS.243a2001M,2011APh....35....1B,2011CQGra..28k4013M,2011GReGr..43..437C,2011ivoa.spec.0711S,2011PhRvL.107y1101B,2012JPhCS.363a2022B,
2012PhRvD..85j3004B,2012PhRvD..86h3007B,2013APh....45...56S,2013CQGra..30l3001B,2013JCAP...06..008A,2013PhRvL.110x1101B,2013RvMP...85.1401A,2014PhRvD..90j1301B,2014PhRvD..90j2002A,2015PhRvL.115w1101B,2016PhRvD..93l2010A,2017ApJ...848L..12A,2017ApJ...850L..35A,2017PhRvD..96b2005A,2017PhRvD..96b3003B}.

The third leg of multi-messenger astrophysics will be the discovery of GWs and high-energy neutrinos from a common source \cite{2013RvMP...85.1401A,2013CQGra..30l3001B}. Such a detection could shed light to, e.g., how newly formed compact objects accelerate particles to extreme energies. In addition, some high-energy neutrinos are identified rapidly with localization accuracies much better than that available with GW detectors, which can guide observatories in their search of the electromagnetic counterparts of GW sources.

Several source candidates are considered to generate GWs and high-energy neutrinos, including core-collapse supernovae~\cite{2012PhRvD..86h3007B, kohta-sn-2018}, gamma-ray bursts (GRBs)~(See e.g.~\cite{kohta2006,meszaros13}), BNS mergers~\cite{2018arXiv180511613K}, neutron star-black hole  mergers~\cite{shigeo2017}, soft gamma repeaters~\cite{ioka2005,murphy2013}, and microquasars~\cite{2011GReGr..43..437C}.
Besides these candidate sources, searches might reveal unknown source populations or production mechanisms.
Detecting even one joint source of GWs and high-energy neutrinos will significantly increase our understanding of the underlying mechanisms that create them~\cite{2013RvMP...85.1401A,2013CQGra..30l3001B}.

Searching for joint GW+high-energy neutrino (hereafter GW+neutrino) sources has only become viable in recent years with the advent of large-scale detectors, in particular the Advanced LIGO \cite{aligo2015} and Advanced Virgo \cite{avirgo2015} observatories on the GW side, and the IceCube \cite{icecube2017}, ANTARES \cite{2011NIMPA.656...11A} and Pierre Auger \cite{2015arXiv150201323T} observatories on the neutrino side. Both sides will experience significant upgrades in the coming years. Advanced LIGO and Advanced Virgo are set to reach their design sensitivities within the next few years \cite{lrr2018}. IceCube started an upgrade towards a second generation detector, IceCube-Gen2, with several times improved sensitivity \cite{2014arXiv1412.5106I}. Another neutrino detector, KM3NeT, is being constructed in the Mediterranean \cite{2016JPhG...43h4001A}. Due to these advances, our reach to GW and neutrino sources is set to rapidly increase in the near future and beyond.

While no joint GW+neutrino discovery has been confirmed to date, there has been significant effort to search for such events. Following the first observational constraints on common sources in 2011 \cite{2011PhRvL.107y1101B}, independent searches were carried out using Initial LIGO/Virgo and the partially completed ANTARES and IceCube detectors \cite{2013JCAP...06..008A,2014PhRvD..90j2002A}. With the completion of Advanced LIGO, several searches were carried out to find the neutrino counterpart of GW discoveries \cite{2016PhRvD..93l2010A,2017PhRvD..96b2005A,2017ApJ...850L..35A}. A separate search was carried out to find joint events for which neither the GW nor the neutrino signal could be independently confirmed to be astrophysical \cite{Albert_2019}.

Most of these searches were based on the analysis method developed by Baret et al. \cite{2012PhRvD..85j3004B}. This method combines GW amplitude, neutrino reconstructed energy, temporal coincidence and directional coincidence to separate astrophysical events from chance coincidences. The method aims to be emission model agnostic and does not impose constraints on the source properties except by assuming that higher neutrino energy is more likely to indicate an astrophysical signal.

Following the success of the search method by Baret et al. \cite{2012PhRvD..85j3004B} spanning over a decade, it is time to upgrade it to enhance its sensitivity and aid newly relevant real-time searches. Two particular motivations for the upgrade are to facilitate the incorporation of astrophysical information and detector characteristics in the search. Regarding astrophysical information, while it is beneficial to keep the search largely model independent, in many cases signal constraints can be specified that do not depend strongly on particular model. Regarding detector characteristics, a more complex detector model will improve sensitivity and accuracy, but requires the incorporation of prior information on these characteristics to the search.

In this paper, we present a new search algorithm for common sources of GWs and high-energy neutrinos based on Bayesian hypothesis testing. A Bayesian framework is a natural choice to incorporate prior astrophysical and detector information. Bayesian solutions are becoming more common in GW \cite{2010PhRvD..81f2003V,2005PhRvD..72j2002D,2015CQGra..32m5012C,2016PhRvD..93b4013S} and more recently multi-messenger data analysis \cite{2008ApJ...679..301B,2018ApJ...860....6A,1997scma.conf.....B,2013ApJS..209...30N,2014ApJ...795...43F,2017PhRvL.119r1102F}.

The paper is organized as follows. The general idea for this analysis is described in Sec~\ref{sec:bayes}, following by probabilities describing signal hypothesis in Sec~\ref{sec:sig}, null hypothesis in Sec~\ref{sec:null} and chance coincidence hypothesis in Sec~\ref{sec:coincidence}. We define the use of odds ratios in Section \ref{sec:odds}. We conclude in Section \ref{sec:conclusion}.

\section{Multi-messenger Search Method}
\label{sec:bayes}
To determine whether a multi-messenger coincident signal is a real event or a random coincidence, we formulate the problem in the context of Bayesian hypothesis testing. We further incorporate detector and background characteristics as well as astrophysical information of the messenger particle and its source.

We will compare multiple hypotheses. Our signal hypothesis, $H_{\rm s}$, is that all considered messengers originated from the same astrophysical source. Our null hypothesis, $H_{\rm 0}$, is that triggers in all messengers arose from the background. Additionally, we will consider a chance coincidence hypothesis, $H_{\rm c}$, that one type of the messengers has an astrophysical trigger, but the other type of messenger only has triggers from the background. We will neglect the possibility that different messengers from distinct astrophysical signals coincide as this is highly unlikely given our low signal rate.

For GWs we use the following observational information for the search: (i) detection time $t_{\rm gw}$; (ii) reconstructed sky location probability density $\mathcal{P}_{\rm gw}=\mathcal{P}_{\rm gw}({\bf \Omega})$, called the {\it skymap}, where ${\bf \Omega}$ is the source sky location; (iii) the GW data analysis pipeline specific signal-to-noise ratio (SNR) of the GW event in the GW detectors network $\rho_{\rm gw}$, which is the individual SNRs of the signals at each detector summed in quadrature, and (iv) reconstructed distance distribution $\mathcal{D}_{\rm gw}=\mathcal{D}_{\rm gw}(r)$ where \(r\) is the distance of the event to Earth \cite{PhysRevLett.116.061102}. We define a vector containing the measured properties of a GW trigger as
\begin{equation}
{\bf x}_{\rm gw} = \{t_{\rm gw}, \mathcal{P}_{\rm gw}, \rho_{\rm gw}, \mathcal{D}_{\rm gw}\}.
\end{equation}
For multiple source types, an additional variable could be the source-dependent gravitational waveform. We omit this as a factor in the following description.

For high-energy neutrinos, the used observational information for them includes (i) their detection times $t_{\nu}$; (ii) their reconstructed sky location probability densities $\mathcal{P}_{\nu}=\mathcal{P}_{\nu}({\bf \Omega})$; and (iii) their reconstructed neutrino energies $\epsilon_{\nu}$. As high-energy neutrinos are not directly observed, the observed energies of the leptons produced in the neutrino interactions are taken as $\epsilon_{\nu}$.  Generally the reconstructed neutrino sky location, $\mathcal{P}_{\nu}$, can be described as a Gaussian distribution centered at reconstructed neutrino direction ${\bf \Omega}_{\nu}$, with reconstructed uncertainty $\sigma_{\nu}$ \cite{2008APh....29..299B,2016JInst..1111009I}. We define a matrix containing the measured properties of all neutrino triggers as
\begin{equation}
    {\bf X}_{\rm \nu}=\begin{bmatrix} {\bf x}_{\rm \nu 1} \\ {\bf x}_{\rm \nu 2} \\ ...\end{bmatrix}
\end{equation}
with rows
\begin{equation}
{\bf x}_{\nu i} = \{t_{\rm \nu i}, {\bf \Omega}_{\rm \nu i}, \sigma_{\rm \nu i}, \epsilon_{\rm \nu i}\}.
\end{equation}
Throughout the paper we will assume we have \(N\) neutrino triggers.
We define a vector containing our model parameters for the signal hypothesis as
\begin{equation}
{\boldsymbol \theta}=\{t_{\rm s}, r, {\bf \Omega}, E_{\rm gw}, E_{\rm \nu}\},
\end{equation}
where $t_{\rm s}$ is the reference time, $r$ is the luminosity distance, ${\bf \Omega}$ is the sky location, $E_{\rm gw}$ is the isotropic-equivalent total GW energy, and $E_{\nu}$ is the isotropic-equivalent total high-energy neutrino energy emitted from the astrophysical event. The reference time can be thought of as the time of a relevant astrophysical event to which we compare the other times of arrival, delayed by the travel time of information to Earth at the speed of light.
The neutrino energies considered here render the neutrino travel time practically the same as travel time at the speed of light.

At the end of our analysis we will compute a Bayes factor for our signal hypothesis given the observational data as
\begin{equation}
\label{eq:or}
\mathcal{O_{\rm gw+\nu}} = \frac{P(H_{\rm s}|{\bf x}_{\rm gw}, {\bf X}_{\rm \nu})}{P(H_{\rm 0}|{\bf x}_{\rm gw}, {\bf X}_{\rm \nu}) + P(H_{\rm c}|{\bf x}_{\rm gw}, {\bf X}_{\rm \nu})}.
\end{equation}

\section{Signal hypothesis}
\label{sec:sig}

We first introduce our signal hypothesis $H_{\rm s}$. This hypothesis considers having at least one coincident signal neutrino with the gravitational wave which is also signal. Therefore we split this hypothesis into sub-hypotheses for different number of coincident signal neutrinos and denote them by $H_{\rm s}^{n}$ where n is the number of coincident neutrinos. What this means is for example for n=1 we would have one neutrino which comes from the same source of gravitational wave and other neutrinos belong to background, or to the null hypothesis. In order to label the signal neutrinos and the background neutrinos separately we will use the notation ${\bf X}_{\rm \nu}^{\rm i}={\bf X}_{\rm \nu}\setminus{\bf x}_{\nu i}$ to refer to the ${\bf X}_{\rm \nu}$ matrix without the $\mbox{i}^{\rm th}$ row $={\bf x}_{\nu i}$. Given the observational data, the probability of the signal hypothesis being true can be written as
\begin{equation}P(H_{\rm s}|{\bf x}_{\rm gw}, {\bf X}_{\rm \nu})=\sum_{n=1}^{N} P(H_{\rm s}^{n}|{\bf x}_{\rm gw}, {\bf X}_{\rm \nu})
\end{equation}
We use Bayes' rule to this expression
\begin{equation}
\sum_{n=1}^{N} P(H_{\rm s}^{n}|{\bf x}_{\rm gw}, {\bf X}_{\rm \nu})=\sum_{n=1}^{N} \frac{P({\bf x}_{\rm gw},{\bf X}_{\rm \nu}|H_{\rm s}^{n})P(H_{\rm s}^{n})}{P({\bf x}_{\rm gw},{\bf X}_{\rm \nu})}
\label{eq:1}
\end{equation}
We are interested in the ratio of such probabilities for different hypotheses, hence the denominator above will cancel out. We therefore omit its computation.
Then we further expand the first term by specifying the signal neutrinos
\begin{widetext}
\begin{equation}
P({\bf x}_{\rm gw},{\bf X}_{\rm \nu}|H_{\rm s}^{n})=\sum_{\{i,j,...\}}P({\bf x}_{\rm gw},{\bf X}_{\rm \nu}|H_{\rm s}^{n},s=\{i,j,...\})P(s=\{i,j,...\}|H_{\rm s}^{n})
\label{eq:1.5}
\end{equation}
\end{widetext}
The sum in Eq. \ref{eq:1.5} is over all \({N \choose n}\) subsets of the set of integers from 1 to \(N\) with \(n\) elements, which are denoted with the set \(s\) in the sum which stands for the signal set
The second term on the right hand side corresponds to probability of each combination which is
\begin{equation}
P(s=\{i,j,...\}|H_{\rm s}^{n})={N \choose n}^{-1}
\label{eq:1.55}
\end{equation}
We further decompose the first term in Eq. \ref{eq:1.5} by separating the signal and background neutrino terms via using their independence with a memoryless detector assumption as
\begin{equation}
\begin{split}
P({\bf x}_{\rm gw},{\bf X}_{\rm \nu}|H_{\rm s}^{n},s=\{i,j,...\})=\\P({\bf x}_{\rm gw},{\bf x}_{\nu i},{\bf x}_{\nu j},...|H_{\rm s}^{n})P({\bf X}_{\rm \nu}^{i,j,...}|H_{\rm 0})
\end{split}
\label{eq:1.6}
\end{equation}
In Eq. \ref{eq:1.6} we dropped the \(s\) set from the conditions of the first term on the right side of the equation since there are already only \(n\) neutrinos in that probability, for convenience. Second term corresponds to the probability for all other than those \(n\) neutrinos belonging to background. Next, to obtain the first term in the right side of Eq. \ref{eq:1.6} we marginalize over the parameters
\begin{equation}
\begin{split}
P({\bf x}_{\rm gw}, {\bf x}_{\nu i},{\bf x}_{\nu j}...|H_{\rm s}^{n}) = \\ \int P({\bf x}_{\rm gw}, {\bf x}_{\nu i},{\bf x}_{\nu j}...|{\boldsymbol \theta}, H_{\rm s}^n)P({\boldsymbol \theta}|H_{\rm s}^n)d{\boldsymbol \theta}
\label{eq:2}
\end{split}
\end{equation}
Since ${\bf x}_{\rm gw}$ and all ${\bf x}_{\nu i}$ which belong to signal hypothesis are dependent on ${\boldsymbol \theta}$ but are otherwise can be considered independent of each other, we can separate the GW and the high-energy parts from each other such as
\begin{equation}
\begin{split}
P({\bf x}_{\rm gw}, {\bf x}_{\nu i},{\bf x}_{\nu j}...|{\boldsymbol \theta}, H_{\rm s}^n) = \\ P({\bf x}_{\rm gw}|{\boldsymbol \theta}, H_{\rm s}^n)P({\bf x}_{\nu i}|{\boldsymbol \theta}, H_{\rm s}^n)P({\bf x}_{\nu j}|{\boldsymbol \theta}, H_{\rm s}^n)...
\label{eq:3}
\end{split}
\end{equation}
We now specify the independent elements of Eqs. \ref{eq:1}, \ref{eq:2} and \ref{eq:3} in the context of our astrophysical and detection models. \( P({\bf X}_{\rm \nu}^{i,j,...}|H_{\rm 0})\) in Eq. \ref{eq:1.6} will be explained in section \ref{sec:null}.
\subsection{Parameter and hypothesis priors ($H_{\rm s}$)}
\label{sec:priorsHs}
There are two types of prior probabilities that we need to compute in our signal hypothesis: $P({\boldsymbol \theta}|H_{\rm s}^{n})$ and $P(H_{\rm s}^{n})$. In order to find $P({\boldsymbol \theta}|H_{\rm s}^{n})$ we again use Bayes' rule
\begin{equation}
P({\boldsymbol \theta}|H_{\rm s}^{n})=\frac{P(H_{\rm s}^n|{\boldsymbol \theta})P({\boldsymbol \theta})}{P(H_{\rm s}^n)}
\label{eq:4}
\end{equation}
\(P(H_{\rm s}^n)\) in the denominator and in Eq \ref{eq:1} cancel. So actually we need to have \(P({\boldsymbol \theta})\) and \(P(H_{\rm s}^n|{\boldsymbol \theta})\).
We first discuss the prior probability distribution of the parameters, $P({\boldsymbol \theta})$. Here we review the role of each source parameter.

\begin{enumerate}
\item {\bf Time ($t_{\rm s}$):} We assume that a signal is equally likely to occur at any time during an observation period. We further assume that no other parameter depends on the time of observation, therefore we can treat this probability independently. Taking the livetime duration of the joint observation period to be $T_{\rm obs}$, the resulting prior probability distribution is
\begin{equation}
P(t_{\rm s}) = \frac{1}{T_{\rm obs}}.
\label{eq:time}
\end{equation}
\item {\bf Source distance ($r$):} We assume a uniform distribution of sources in the universe \textcolor{black}{such that \(P(r)=3r^2/r_{max}^3\). \(r_{max}\) is the maximum value of \(r\) and its divergence is not important as it gets cancelled in the analysis.} We further assume that a GW signal can be detected if its root-sum-squared GW strain $h_{\rm rss}$ is above a detection threshold $h_{\rm rss,0}$ \cite{2010CQGra..27q3001A}. The probability density that an observed GW+neutrino event occurred at distance $r$ is dependent on $r$ since the volume in space in the distance range $[r,r+dr]$ is $\propto r^2 dr$; but the probability of detecting \(n\) neutrinos from the source falls according to Poisson probabilities for \(n\) detections whose means are proportional to $r^{-2}$. This dependency is valid up to the GW distance range $r_0 f_{\rm A}({\bf \Omega},t_{\rm gw})$ beyond which sources are not detected. Here, $r_0$ is the GW detection range for optimal source direction, and $f_{\rm A}({\bf \Omega},t_{\rm gw})$ is the antenna pattern of the GW detector network. \textcolor{black}{The latter is the square root of the quadrature sum of the antenna responses of each detector for the two polarizations of GWs as
\begin{equation}
f_A({\bf \Omega},t_{\rm gw})=\sqrt{\sum_k {F_{k,+}(\psi,{\bf \Omega},t_{\rm gw})}^2+{F_{k,\times}(\psi,{\bf \Omega},t_{\rm gw})}^2}
\end{equation} where sum over \(k\) is for summing over different detectors, \(F\) functions are the antenna responses for each polarization and \(\psi\) angle is the GW emission inclination which vanishes after the quadrature sum. }Its maximum value is 1 which corresponds to the optimal source position. The range $r_0$ satisfies $r_0(E_{\rm gw})\propto E_{\rm gw}^{1/2}$ and is defined as the distance at which an event which emits \(E_{\rm gw}\) amount of energy creates the least acceptable SNR \(\rho_{\rm min}\) or the strain \(h_{\rm rss,0}\). The relationship between \(E_{\rm gw},\rho_{\rm gw}\) and \(r_0\) is explained in Eq. \ref{eq:delta} with \(\rho_{\rm gw}=\rho_{\rm min}\)  and \(r=r_0\) such that \(r_0=\kappa_0 E_{\rm gw}^{1/2}/\rho_{\rm min}\).

\item {\bf Energy ($E_{\rm gw}$ and $E_{\rm \nu}$):} We need to specify our dependency to energies. A naive choice can be independent log-uniform distributions over the energy ranges \([E_{\rm gw}^{-}, E_{\rm gw}^{+}], [E_{\rm \nu}^{-}, E_{\rm \nu}^{+}]\) with probability density
\begin{equation}
\begin{split}
P(E_{\rm gw},E_{\rm \nu})=P(E_{\rm gw})P(E_{\nu})\\=(E_{\rm gw}E_{\rm \nu}\log(\frac{E_{\rm gw}^{+}}{E_{\rm gw}^{-}})\log(\frac{E_{\rm \nu}^{+}}{E_{\rm \nu}^{-}}))^{-1}
\end{split}
\end{equation}
Throughout the paper instead of the expression for the specific log-uniform model, \(P(E_{\rm gw},E_{\nu})\) will be used to express the universality of the method.
\item {\bf Sky position (\({\bf \Omega}\)):} We assume uniform prior distribution in the sky as \(P({\bf \Omega})=1/4\pi\)
\end{enumerate}
Overall, we find

\begin{equation}
P({\boldsymbol \theta}) =  \frac{P(E_{\rm gw},E_{\rm \nu})r^2}{4\pi T_{\rm obs}N_1}
\label{eq:parameterprior}
\end{equation}

with suitable normalization constant $N_1=r^3_{max}/3$ for maximum assumed possible distance \(r_{max}\). Its divergence is not important and all divergences in the analysis cancel.

Next we consider the term \(P(H_{\rm s}^n|{\boldsymbol \theta})\) which depends on the expected detection count of multi-messenger events.  A useful quantity for it is the expected number of detected neutrinos for given emission energy, sky location and distance per event
\begin{equation}
\langle n_{\nu}(E_{\rm \nu},r,{\bf \Omega})\rangle = n_{\nu,51,100}({\bf \Omega}) \left(\frac{E_{\nu}}{10^{51}\,\mbox{erg}}\right)\left(\frac{r}{100\,\mbox{Mpc}}\right)^{-2}
\end{equation}
Here, $n_{\nu,51,100}({\bf \Omega})$ is a detector specific parameter. Sky averaged $n_{\nu,51,100}\approx1.1$ for IceCube \cite{2018arXiv181011467B} and  $n_{\nu,51,100}({\bf \Omega})$ only depends on the declination but not to right ascension due to axial symmetry of the detector whose symmetry axis coincides with Earth's rotation axis due to its location at the South Pole.

Probability of detection of multi-messenger events given the source parameters will be
\begin{multline}
P_{det}^n({\boldsymbol \theta})=\mbox{Poiss}(n,\langle n_{\nu}(E_{\nu},r,{\bf \Omega})\rangle)\\ \begin{cases} 1 &  r\leq r_0(E_{\rm gw})\Bar{f}_A({\bf \Omega},t_s) \\ 0 & \text{otherwise}
\end{cases}
\end{multline}
with the time averaged antenna pattern \(\Bar{f}_A({\bf \Omega},t_s)\) between \([t_s+t_{\rm gw}^-,t_s+t_{\rm gw}^+]\). Parameters \(t_{\rm gw}^-\) and \(t_{\rm gw}^+\) will be explained in section \ref{sec:GWHs}. Poiss$(k,\lambda)$ is the Poisson probability density function with $\lambda$ mean and $k$ observed events. Then the expected count is
\begin{equation}
C^n_{det}({\boldsymbol \theta})=\dot{n}_{\rm gw+\nu}T_{obs}P^n_{det}({\boldsymbol \theta})
\end{equation}
where $\dot{n}_{\rm gw+\nu}$ is the total multi-messenger event rate in the whole universe, which is bounded by the distance \(r_{max}\). Its divergence cancels the divergence of \(N_1\). Then the prior probability will be
\begin{equation}
P(H_{\rm s}^{n}|{\boldsymbol \theta}) = \frac{1}{N_{2}} C^n_{det}({\boldsymbol \theta})
\label{eq:Hs}
\end{equation}
with suitable normalization constant $N_{2}$. This constant, which is the sum of the expected event counts of all hypotheses, will be the present for our other hypotheses too and will be cancelled out.

\subsection{Gravitational waves ($H_{\rm s}$)}
\label{sec:GWHs}

We now consider the probability $P({\bf x}_{\rm gw}|{\boldsymbol \theta}, H_{\rm s}^{n})$
We have:
\begin{multline}
P({\bf x}_{\rm gw}|{\boldsymbol \theta}, H_{\rm s}^{n}) = P(t_{\rm gw}|{\boldsymbol \theta},H_{\rm s}^{n})P(\rho_{\rm gw}|t_{\rm gw},{\boldsymbol \theta},H_{\rm s}^{n})\\ \times P(\mathcal{P}_{\rm gw}|t_{\rm gw},\rho_{\rm gw},{\boldsymbol \theta},H_{\rm s}^{n}) P(\mathcal{D}_{\rm gw}|\mathcal{P}_{\rm gw},t_{\rm gw},\rho_{\rm gw},{\boldsymbol \theta},H_{\rm s}^{n})
\label{eq:234243432
\end{multline}

The term $P(t_{\rm gw}|{\boldsymbol \theta},H_{\rm s}^{n})$ should only depend on the difference $t_{\rm gw}-t_{\rm s}$. We adopt the assumption that the probability $P(t_{\rm gw}|t_{\rm s},H_{\rm s}^{n})$ is uniform within a time window $t_{\rm gw}-t_{\rm s} \in [t_{\rm gw}^{-},t_{\rm gw}^{+}]$ for suitable parameters $t_{\rm gw}^{-}$ and
$t_{\rm gw}^{+}$ and is zero elsewhere:
\begin{equation}
 P(t_{\rm gw}|t_{\rm s},H_{\rm s}^{n})=
 \begin{cases} (t_{\rm gw}^{+}-t_{\rm gw}^{-})^{-1} & \text{if } t_{\rm gw}-t_{\rm s}\in [t_{\rm gw}^{-},t_{\rm gw}^{+}]  \\
                      0                                    & \mbox{otherwise}
 \end{cases}
\end{equation}
For example, previous GW+neutrino searches used parameters $t_{\rm gw}^{+}=-t_{\rm gw}^{-}=250$\,s \cite{2011APh....35....1B,2012PhRvD..85j3004B,2014PhRvD..90j2002A,2016PhRvD..93l2010A,2017PhRvD..96b2005A,2017ApJ...850L..35A}. We assume that the other source parameters are independent of $t_{\rm gw}$.

To understand the second term on the right hand side of Eq. \ref{eq:234243432}, we make use of the fact that $\rho_{\rm gw}$ on average is proportional to the GW signal's amplitude at Earth, characterized by the root-sum-squared GW strain $h_{\rm rss}$. Assuming here for simplicity that all gravitational waveforms are similar, the GW strain is fully determined by $r$, $E_{\rm gw}$, ${\bf \Omega}$ and $t_{\rm gw}$. The time dependence comes in due to Earth's rotation if we measure sky location in equatorial coordinates. Assuming that $\rho_{\rm gw}$ precisely describes $h_{\rm rss}$, this term represents a constraint on the source parameters, which need to be such that they produce at Earth the $h_{\rm rss}$ value that corresponds to the measured $\rho_{\rm gw}$ value. This means that only the combination $E_{\rm gw}^{1/2}r^{-1}f_{\rm A}({\bf \Omega},t_{\rm gw})$ is constrained, where $f_{\rm A}({\bf \Omega},t_{\rm gw})$ is the direction-dependent antenna pattern of the GW detector network. This combination is proportional to the measured GW strain amplitude.

We therefore write the probability as a constraint
\begin{equation}
\label{eq:delta}
P(\rho_{\rm gw}|t_{\rm gw},{\boldsymbol \theta},H_{\rm s}^{n}) = \delta\left[\rho_{\rm gw} - \kappa_0 E_{\rm gw}^{1/2}r^{-1}f_{\rm A}({\bf \Omega},t_{\rm gw})\right]
\end{equation}
Where $\delta$ is the Dirac-delta and $\kappa_0$ is an appropriate constant that depends on the noise spectrum in the GW detector at the time of detection and on the GW search algorithm. The delta distribution is an approximation for the accurate measurement of the SNR although in practice there is always an uncertainty in the measured SNR. We do not consider that uncertainty in our analysis.

Next, we look at the term $P(\mathcal{P}_{\rm gw}|t_{\rm gw},\rho_{\rm gw},{\boldsymbol \theta},H_{\rm s}^{n})$. Using Bayes' rule, we write
\begin{multline}
P(\mathcal{P}_{\rm gw}|t_{\rm gw},\rho_{\rm gw},{\boldsymbol \theta},H_{\rm s}^{n}) \\ = \frac{P({\boldsymbol \theta}|\mathcal{P}_{\rm gw},t_{\rm gw},\rho_{\rm gw},H_{\rm s}^{n})P(\mathcal{P}_{\rm gw}|t_{\rm gw},\rho_{\rm gw},H_{\rm s}^{n})}{P({\boldsymbol \theta}|t_{\rm gw},\rho_{\rm gw},H_{\rm s}^{n})}
\label{eq:45645465}
\end{multline}
Regarding $P(\mathcal{P}_{\rm gw}|t_{\rm gw},\rho_{\rm gw},H_{\rm s}^{n})$, we assume that the distribution of $\mathcal{P}_{\rm gw}$ is independent of the underlying hypothesis, i.e. $P(\mathcal{P}_{\rm gw}|t_{\rm gw},\rho_{\rm gw},H_{\rm s}^{n}) = P(\mathcal{P}_{\rm gw}|t_{\rm gw},\rho_{\rm gw})$. This term appears for the alternative hypothesis as well, it cancels out and therefore we can ignore it here. For the remaining terms assuming that our reconstructed GW skymap is accurate, we have
\begin{equation}
\frac{P({\boldsymbol \theta}|\mathcal{P}_{\rm gw},t_{\rm gw},\rho_{\rm gw},H_{\rm s}^{n})}{P({\boldsymbol \theta}|t_{\rm gw},\rho_{\rm gw},H_{\rm s}^{n})} = \mathcal{P}_{\rm gw}({\bf \Omega})\frac{N_{\bf \Omega}}{f_A({\bf \Omega},t_{\rm gw})}
\end{equation}
The numerator of the left side gives the skymap since the skymap determines the probability of the signal coming from a given sky location ${\bf \Omega}$. The denominator of the left side is proportional to the antenna pattern at the time of detection, normalized by the factor \(N_{\bf \Omega}\)

Now we look at the term $P(\mathcal{D}_{\rm gw}|\mathcal{P}_{\rm gw},t_{\rm gw},\rho_{\rm gw},{\boldsymbol \theta},H_{\rm s}^{n})$. The handling of this term is identical to $P(\mathcal{P}_{\rm gw}|t_{\rm gw},\rho_{\rm gw},{\boldsymbol \theta},H_{\rm s}^{n})$. Again we assume the independence of \(\mathcal{D}_{\rm gw}\) and the hypothesis, and write
\begin{multline}
P(\mathcal{D}_{\rm gw}|\mathcal{P}_{\rm gw},t_{\rm gw},\rho_{\rm gw},{\boldsymbol \theta},H_{\rm s}^{n}) \\ =   \frac{P({\boldsymbol \theta}|\mathcal{P}_{\rm gw},t_{\rm gw},\rho_{\rm gw},\mathcal{D}_{\rm gw},H_{\rm s}^{n})P(\mathcal{D}_{\rm gw}|\mathcal{P}_{\rm gw},t_{\rm gw},\rho_{\rm gw},H_{\rm s}^{n})}{P({\boldsymbol \theta}|\mathcal{P}_{\rm gw},t_{\rm gw},\rho_{\rm gw},H_{\rm s}^{n})}
\label{eq:dist1}
\end{multline}
We again ignore the second term in the denominator as it cancels with the corresponding terms in other hypotheses. Assuming that our reconstructed GW distance distribution is accurate, we have
\begin{equation}
\frac{P({\boldsymbol \theta}|\mathcal{D}_{\rm gw},\mathcal{P}_{\rm gw},t_{\rm gw},\rho_{\rm gw},H_{\rm s}^{n})}{P({\boldsymbol \theta}|\mathcal{P}_{\rm gw}, t_{\rm gw},\rho_{\rm gw},H_{\rm s}^{n})}=\mathcal{D}_{\rm gw}(r)\frac{N_r(\rho_{\rm gw})}{r^2\times r^{-1}}
\end{equation}
The first term in the right side comes from the numerator of the left side since the distance distribution determines the probability of the signal coming from a source distance $r$. The second term in the right side represents the denominator of the left side and is obtained by multiplying the prior \(r^2\) distribution and the \(r^{-1}\) distribution of the likelihood of \(r\) for known \(\rho_{\rm gw}\) which is proportional to the distribution of \(E_{\rm gw}^{1/2}\) which is \(\propto E_{\rm gw}^{-1/2}\) for a log uniformly assumed \(P(E_{\rm gw})\). Here \(N_r(\rho_{\rm gw})\) is the normalization for the \(r\) distribution between \([\sqrt{E^-_{\rm gw}},\sqrt{E^+_{\rm gw}}]\rho_{\rm gw}/\kappa_0\).

Putting everything together, we have for the GW term without the cancelling terms
\begin{widetext}
\begin{equation}
P({\bf x}_{\rm gw}|{\boldsymbol \theta}, H_{\rm s}^{n}) = N_{\bf \Omega} N_r(\rho_{\rm gw}) \delta\left[\rho_{\rm gw} - \kappa_0 E_{\rm gw}^{1/2}r^{-1}f_{\rm A}({\bf \Omega},t_{\rm gw})\right] \frac{\mathcal{P}_{\rm gw}({\bf \Omega})}{f_A({\bf \Omega},t_{\rm gw})}\frac{\mathcal{D}_{\rm gw}(r)}{r}
\begin{cases} (t_{\rm gw}^{+}-t_{\rm gw}^{-})^{-1} & \text{if } t_{\rm gw}-t_{\rm s}\in [t_{\rm gw}^{-},t_{\rm gw}^{+}]  \\
                      0                                    & \mbox{otherwise}
\end{cases}
\label{eq:gwsignal}
\end{equation}
\end{widetext}

\subsection{High-energy neutrinos ($H_{\rm s}$)}

We now turn our attention to the high-energy neutrino term $P({\bf x}_{\nu i},{\bf x}_{\nu j} ...|{\boldsymbol \theta}, H_{\rm s}^{n})$ as in Eq. \ref{eq:3}. We assume no dependency of the observables of different neutrinos on other neutrinos' observables, except the dependency through \({\boldsymbol \theta}\). Therefore we can separate each neutrino terms as
\begin{equation}
P({\bf x}_{\nu i},{\bf x}_{\nu j} ...|{\boldsymbol \theta}, H_{\rm s}^{n})=P({\bf x}_{\nu i}|{\boldsymbol \theta},H_{\rm s}^{n})P({\bf x}_{\nu j}|{\boldsymbol \theta},H_{\rm s}^{n})...
\end{equation}
We treat the temporal term similarly to the GW case. We assume that the time difference $t_{\rm \nu i}-t_{\rm s}$ is the only relevant temporal value. We further use a uniform probability density within the time interval $[t_{\nu}^{-}$,$t_{\nu}^{+}]$, and 0 outside:
\begin{equation}
 P(t_{\nu i}|t_{\rm s},H_{\rm s}^{n})=
 \begin{cases} (t_{\nu}^{+}-t_{\nu}^{-})^{-1} & \text{if } t_{\rm \nu i}-t_{\rm s}\in [t_{\nu}^{-},t_{\nu}^{+}]  \\
                      0                                    & \mbox{otherwise}
 \end{cases}
.
\end{equation}
Previous GW+neutrino searches used parameters $t_{\nu}^{+}=-t_{\nu}^{-}=250$\,s \cite{2011APh....35....1B,2012PhRvD..85j3004B,2014PhRvD..90j2002A,2016PhRvD..93l2010A,2017PhRvD..96b2005A,2017ApJ...850L..35A}.

The remaining neutrino observables, ${\bf \Omega}_{\rm \nu i}$, $\sigma_{\rm \nu i}$ and $\epsilon_{\rm \nu i}$, are not independent. The sensitivity of neutrino detectors varies with both energy and sky location, and localization accuracy depends on source direction and energy.

Let us take the remaining neutrino term \(P({\bf \Omega}_{\rm \nu i}, \sigma_{\rm \nu i}, \epsilon_{\rm \nu i} | r, {\bf \Omega}, E_{\nu}, H_{\rm s}^{n})\). We assume that the signal distribution of $\epsilon_{\nu}$ follows a power law, therefore the neutrino spectrum is independent of the source distance. Such a power-law distribution is typical in neutrino emission models \cite{1997PhRvL..78.2292W}. Consequently, parameters $r$ and $E_{\nu}$ do not affect the probability here. We further assume that the directional uncertainty variable $\sigma_{\nu}$ and the reconstructed sky position of neutrino \({\bf \Omega}_{\nu}\) don't depend on \(r\) and \(E_{\nu}\).

We use the chain rule to write
\begin{widetext}
\begin{equation}
P({\bf \Omega}_{\rm \nu i}, \epsilon_{\rm \nu i }, \sigma_{\rm \nu i}| {\bf \Omega}, H_{\rm s}^{n}) = P({\bf \Omega}_{\rm \nu i} | \sigma_{\rm \nu i}, \epsilon_{\rm \nu i}, {\bf \Omega}, H_{\rm s}^{n})P(\sigma_{\rm \nu i}|\epsilon_{\rm \nu i},{\bf \Omega}, H_{\rm s}^{n})P(\epsilon_{\rm \nu i}|{\bf \Omega}, H_{\rm s}^{n})
\label{eq:76567}
\end{equation}
\end{widetext}
Given the source direction as a parameter, the probability of reconstructing $\epsilon_\nu$ for a detected neutrino depends on the energy- and direction-dependent effective area $A_{\rm eff}(\epsilon_{\nu},{\bf \Omega})$ of the neutrino detector, as well as the source power spectral density. Here we ignore the difference between true and reconstructed energy when calculating the effective area as this should not significantly change its value. We take the neutrino spectral density to be $dN_{\nu}/d\epsilon_{\nu}\propto \epsilon_{\nu}^{-2}$, which is the standard spectrum expected from Fermi processes \cite{1997PhRvL..78.2292W}. With these dependencies, we write
\begin{equation}
P(\epsilon_{\rm \nu i} | {\bf \Omega}, H_{\rm s}^{n}) = \frac{1}{N_{\epsilon}}A_{\rm eff}(\epsilon_{\rm \nu i},{\bf \Omega})\epsilon_{\rm \nu i}^{-2},
\end{equation}
where
\begin{equation}
N_{\epsilon}=\int d{\bf \Omega} \int_{\epsilon_{\rm min}}^{\epsilon_{\rm max}} \epsilon_{\nu}^{-2}A_{\rm eff}(\epsilon_{\nu},{\bf \Omega})d\epsilon_{\nu}
\end{equation}
where the \(d {\bf \Omega}\) integral is over all sky and \(\epsilon_{\rm min}\) and \(\epsilon_{\rm max}\) are the minimum and maximum reconstructable energies.

We assume \(P(\sigma_{\rm \nu i}|\epsilon_{\rm \nu i})\) does not depend on the hypothesis of consideration or its parameters; therefore will be cancelled with the same term in the other hypotheses when the ratio of probabilities is taken at the very end. Hence we do not consider the actual value of  \(P(\sigma_{\rm \nu i}|\epsilon_{\rm \nu i},{\bf \Omega},H_{\rm s}^{n})\).

For the first term on the right-hand side in Eq. \ref{eq:76567}, we adopt the Normal distribution \cite{2008APh....29..299B}
\begin{equation}
P({\bf \Omega}_{\rm \nu i} | \sigma_{\rm \nu i}, \epsilon_{\rm \nu i}, {\bf \Omega}, H_{\rm s}^{n}) = \frac{1}{2\pi\sigma_{\rm \nu i}^{2}}\mbox{e}^{-\frac{|{\bf \Omega}_{\rm \nu i} - {\bf \Omega}|^2}{2\sigma_{\rm \nu i}^{2}}}
\end{equation}
by assuming no further dependence for \({\bf \Omega}_{\rm \nu i}\) on \(\epsilon_{\rm \nu i}\) except the one through \({\bf \Omega}\). Putting everything together, we have for the neutrino term
\begin{widetext}
\begin{equation}
P({\bf x}_{\nu i}|{\boldsymbol \theta}, H_{\rm s}^{n}) = \frac{1}{N_{\epsilon}}A_{\rm eff}(\epsilon_{\rm \nu i},{\bf \Omega})\epsilon_{\rm \nu i}^{-2} \frac{1}{2\pi\sigma_{\rm \nu i}^{2}}\mbox{e}^{-\frac{|{\bf \Omega}_{\rm \nu i} - {\bf \Omega}|^2}{2\sigma_{\rm \nu i}^{2}}}
\begin{cases} (t_{\nu}^{+}-t_{\nu}^{-})^{-1} & \text{if } t_{\rm \nu i}-t_{\rm s}\in [t_{\nu}^{-},t_{\nu}^{+}]  \\
                      0                                    & \mbox{otherwise}
\end{cases}
\label{eq:neutrinosignal}
\end{equation}
\end{widetext}

\subsection{Combination of probabilities ($H_{\rm s}$)}

We can combine the above results to obtain the probability of the joint event being a signal by taking Eqs. \ref{eq:gwsignal} and \ref{eq:neutrinosignal} and substituting them into Eq. \ref{eq:3}. We then substitute Eq. \ref{eq:3} into Eq. \ref{eq:2}.

To solve Eq. \ref{eq:2}, we further substitute $P({\boldsymbol \theta}|H_{\rm s}^{n})$ from Eq. \ref{eq:4} for which we use Eqs. \ref{eq:parameterprior} and \ref{eq:Hs}. Then, we can substitute Eq. \ref{eq:2} into Eq. \ref{eq:1.6}. Finally we substitute Eq.  \ref{eq:1.6} into Eq. \ref{eq:1.5} with Eq. \ref{eq:1.55} which is the required term for Eq. \ref{eq:1}, where we obtain $P(H_{\rm s}^{n}|{\bf x}_{\rm gw}, {\bf X}_{\rm \nu})$, except the factor $P({\bf x}_{\rm gw},{\bf X}_{\rm \nu})$ which will cancel out in comparison to the alternative hypothesis. Computation of \(P({\bf X}_{\rm \nu}^{i,j,...}|H_{\rm 0})\) will be explained in section \ref{sec:null}.

\section{Null hypothesis}
\label{sec:null}

We now move to our null hypothesis $H_{\rm 0}$. Given the observational data, the probability of the null hypothesis being true can be written as $P(H_{\rm 0}|{\bf x}_{\rm gw}, {\bf X}_{\rm \nu})$. We apply Bayes' rule to express this probability as
\begin{equation}
P(H_{\rm 0}|{\bf x}_{\rm gw}, {\bf X}_{\rm \nu}) = \frac{P({\bf x}_{\rm gw}, {\bf X}_{\rm \nu}|H_{\rm 0})P(H_{\rm 0})}{P({\bf x}_{\rm gw},{\bf X}_{\rm \nu})},
\label{eq:1b}
\end{equation}
The denominator here will cancel out with the same denominator in the signal hypothesis, therefore we do not need to further consider it.
Since the background events for GW and neutrino observations are independent, we can write
\begin{equation}
P({\bf x}_{\rm gw}, {\bf X}_{\rm \nu}|H_{\rm 0}) = P({\bf x}_{\rm gw}|H_{\rm 0})P({\bf X}_{\rm \nu}|H_{\rm 0})
\label{eq:3b}
\end{equation}
We will now specify the independent elements of Eqs. \ref{eq:1b} and \ref{eq:3b} in the context of our background model. We can perform the calculations in this section without a need for additional parameters to marginalize over due to the fact that every measured parameter is assumed to be independent of each other for background.

\subsection{Hypothesis prior ($H_{\rm 0}$)}

There is one prior probability that we need to compute in our null hypothesis: $P(H_{\rm 0})$. This probability is again proportional to the expected detection count of background events. Given the observation period \(T_{obs}\) and the background GW and rate \(R_{\rm gw,bg}\), we have \(T_{obs}R_{\rm gw,bg}\) background GW events. Hence
\begin{equation}
P(H_{\rm 0}) = \frac{1}{N_{2}}R_{\rm gw,bg}T_{obs}
\label{eq:pH0}
\end{equation}
The normalization factor $N_{2}$ will cancel out with the same factor in the signal hypothesis, see Eq. \ref{eq:Hs}.

\subsection{Gravitational waves ($H_{\rm 0}$)}

We now consider the GW component $P({\bf x}_{\rm gw}|H_{\rm 0})$. We assume \(t_{\rm gw}\) and \(\rho_{\rm gw}\) are independent. We can then define the probabilities of measuring each parameter independently:
\begin{multline}
P({\bf x}_{\rm gw}| H_{\rm 0}) =
P(t_{\rm gw}| H_{\rm 0}) P(\rho_{\rm gw}|H_{\rm 0})\\ \times P(\mathcal{P}_{\rm gw}|t_{\rm gw},\rho_{\rm gw},H_{\rm 0})P(\mathcal{D}_{\rm gw}|t_{\rm gw},\rho_{\rm gw},\mathcal{P}_{\rm gw},H_{\rm 0})
\end{multline}

Starting with the first term on the right-hand side, we expect the probability distribution of detection time for a background event to be independent of time, therefore we adopt a uniform distribution within the observation time. We therefore have
\begin{equation}
 P(t_{\rm gw}|H_{\rm 0})=\frac{1}{T_{obs}}
 \label{eq:H0gwtime}
\end{equation}

The distribution of $\rho_{\rm gw}$ depends on the detector properties as well as the properties of the reconstruction algorithm. We therefore estimate this distribution empirically. \textcolor{black}{We use the $\rho_{\rm gw}$ background distribution obtained from the GW search pipelines which produce it by time shifting data between multiple GW observatories and carrying out the full analysis algorithm over this time shifted data. We do not extrapolate beyond the maximum $\rho_{\rm gw}$ from the background distribution and conservatively choose it as the highest value for the background triggers.} We denote the empirically established distribution of $\rho_{\rm gw}$ with $P_{\rm emp}(\rho_{\rm gw}|H_{\rm 0})$. The \(\rho_{\rm gw}\) obtained by time shifts will be acquired from the GW data analysis pipelines.

Considering the terms $P(\mathcal{P}_{\rm gw}|t_{\rm gw},\rho_{\rm gw},H_{\rm 0})$ and $P(\mathcal{D}_{\rm gw}|t_{\rm gw},\rho_{\rm gw},\mathcal{P}_{\rm gw},H_{\rm 0})$, we do not have any prior information on $P(\mathcal{P}_{\rm gw}|H_{\rm 0})$ and $P(\mathcal{D}_{\rm gw}|H_{\rm 0})$, therefore we assume that it is independent of $\mathcal{P}_{\rm gw}$ and $\mathcal{D}_{\rm gw}$. Since there are similar terms in our signal hypothesis, these cancel out.  We therefore ignore these terms in the following.

Putting everything together, we have for the background GW term
\begin{equation}
P({\bf x}_{\rm gw}| H_{\rm 0}) = P_{\rm emp}(\rho_{\rm gw}|H_{\rm 0})\frac{1}{T_{obs}}
\label{eq:gwsignalH0}
\end{equation}

\subsection{High-energy neutrinos ($H_{\rm 0}$)}

Next, we examine $P({\bf X}_{\rm \nu}|H_{\rm 0})$ in Eq. \ref{eq:3b} and also \(P({\bf X}_{\rm \nu}^{i,j,...}|H_{\rm 0})\) in Eq. \ref{eq:1.6}. Given the background neutrino rate \(R_{\rm \nu, bg}\), the probability of having N background neutrinos in the observation period is \({\rm Poiss}(N,R_{\rm \nu, bg}T_{obs})\), which allows us to write
\begin{equation}
 P({\bf X}_{\rm \nu}|H_{\rm 0})=P({\bf X}_{\rm \nu}|H_{\rm 0},\#_{bg}=N)P(\#_{bg}=N|H_{\rm 0})
 \label{eq:sumbg}
\end{equation}
with \(P(\#_{bg}=N|H_{\rm 0}) ={\rm Poiss}(N,R_{\rm \nu, bg}T_{obs})\) where \(\#_{bg}\) is the number of background neutrinos. For short notation we will use \(H_{\rm 0}^N\) instead of both \(H_{\rm 0}\) and \(\#_{bg}=N\). Now we can decompose the first term into single neutrino terms as we did previously in section \ref{sec:sig}. Then we first separate the temporal term which we assume to be independent of the other parameters. We assume that the time of arrival of a background neutrino signal is time-independent, and can be anytime during the observation period. We therefore have for each neutrino
\begin{equation}
 P(t_{\nu i}|H_{\rm 0}^N)=\frac{1}{{T_{obs}}}
 \label{eq:H0nutime}
\end{equation}

The remaining measured parameters will not be independent of each other. In particular, the reconstructed neutrino direction and energy are interconnected. As we explained before the term for directional uncertainty parameter $ \sigma_{\nu i}$ will be cancelled when we decompose the remaining terms as
\begin{multline}
P({\bf \Omega}_{\rm \nu i}, \epsilon_{\rm \nu i }, \sigma_{\rm \nu i}| H_{\rm 0}^N) = \\ P({\bf \Omega}_{\rm \nu i} | \sigma_{\rm \nu i}, \epsilon_{\rm \nu i}, H_{\rm 0}^N)P(\sigma_{\rm \nu i}|\epsilon_{\rm \nu i},H_{\rm 0}^N)P(\epsilon_{\rm \nu i}|H_{\rm 0}^N)
\end{multline}
In addition we assume \(P({\bf \Omega}_{\rm \nu i} | \sigma_{\rm \nu i}, \epsilon_{\rm \nu i}, H_{\rm 0}^N)\) term doesn't have any \(\sigma_{\nu i}\) dependence. We therefore effectively need to examine the probability $P({\bf \Omega}_{\rm \nu i} | \epsilon_{\rm \nu i}, H_{\rm 0}^N)P(\epsilon_{\rm \nu i}| H_{\rm 0}^N)= P({\bf \Omega}_{\nu},\epsilon_{\nu}|H_{\rm 0}^N)$. Given a sufficient number of observed background events, this probability can be estimated empirically using observed data. Let $\{{\bf \Omega}_{\rm \nu,j},\epsilon_{\rm \nu,j}\}, j\in N_{\rm \nu,obs}$ be the reconstructed parameters observed set of $N_{\rm \nu,obs}$ neutrino candidates. For the direction we only care about the declination angle in the equatorial coordinate system \(\delta_{\nu i}({\bf \Omega}_{\nu i})\) primarily because of the axial symmetry for IceCube, which is described in section \ref{sec:priorsHs} when commenting on \(n_{\rm \nu,51,100}({\bf \Omega})\). For non-coaxial detectors with Earth's rotation axis, full \({\bf \Omega}\) should be considered. We then have the empirical estimate with the kernel density estimation
\begin{widetext}
\begin{equation}
P_{\rm emp}({\bf \Omega}_{\nu i},\epsilon_{\nu i}|H_{\rm 0}^N) =\frac{\sum_{j\in N_{\rm \nu,obs}} \left[ |\delta_{\nu i} - \delta_{\rm \nu,j}| < \Delta_\delta \,\,\,\,\&\,\,\,\, |\epsilon_{\rm \nu i} - \epsilon_{\rm \nu,j}| < \Delta_\epsilon(\epsilon_{\rm \nu i})  \right]}{4\pi N_{\rm \nu,obs} |cos(\delta_{\nu i}+\Delta_\delta)-cos(\delta_{\nu i}-\Delta_\delta)|\Delta_\epsilon(\epsilon_{\rm \nu i})}
\end{equation}
\end{widetext}
where we use here the bracket notation such that $[P]$ is $1$ if $P$ is true and $0$ if $P$ is false which corresponds to the top hat kernel. We further introduced constants $\Delta_\delta$ and $\Delta_\epsilon(\epsilon_{\rm \nu})$, which should be selected such that the uncertainty on the probability estimate is minimal.

Putting everything together, we have for the background neutrino term
\begin{equation}
P({\bf x}_{\nu i}| H_{\rm 0}^N) = P_{\rm emp}({\bf \Omega}_{\nu i},\epsilon_{\nu i}|H_{\rm 0}^N)\frac{1}{T_{obs}}
\label{eq:neutrinoH0full}
\end{equation}

\subsection{Combination of probabilities ($H_{\rm 0}$)}

We can combine the above results to obtain the probability of the joint event being from the background by taking Eq. \ref{eq:neutrinoH0full} for each neutrino and substituting them in Eq. \ref{eq:sumbg}. Then Eq. \ref{eq:gwsignalH0} and Eq. \ref{eq:sumbg} can be substitued into Eq. \ref{eq:3b}. Finally Eq. \ref{eq:3b} along with Eq. \ref{eq:pH0} can be substituted into Eq. \ref{eq:1b}. Eq. \ref{eq:1b} will miss a normalization factor both from Eq. \ref{eq:pH0} and from the denominator on the right side, both of which cancel out upon calculating the Bayes factor. For the background terms in other hypotheses such as \(P({\bf X}_{\rm \nu}^{i,j,...}|H_{\rm 0})\) in Eq. \ref{eq:1.6}, Eq. \ref{eq:sumbg} can be used similarly for \(N-n\) number of background neutrinos instead of \(N\).

\section{Chance coincidence hypothesis}
\label{sec:coincidence}
We finally calculate the probability for the chance coincidence hypothesis $H_{\rm c}$. Given the observational data the probability of the chance coincidence hypothesis being true can be written as $P(H_{\rm c}|{\bf x}_{\rm gw}, {\bf X}_{\rm \nu})$. $H_{\rm c}$ can be separated into two parts, one of which considers a background neutrino event and a foreground gravitational wave event denoted by $H_{\rm c, gw}$; and the other one considers a background gravitational wave event and a foreground neutrino event denoted by $H_{\rm c, \nu}$. Since these two cases are mutually exclusive and complementary to each other for the chance coincidence hypothesis we can write \(P(H_{\rm c}|{\bf x}_{\rm gw}, {\bf X}_{\rm \nu}) = P(H_{\rm c, gw}|{\bf x}_{\rm gw}, {\bf X}_{\rm \nu}) + P(H_{\rm c, \nu}|{\bf x}_{\rm gw}, {\bf X}_{\rm \nu})\). \(H_{\rm c, gw}\) is a special case for the signal hypothesis \(H_{\rm s}^n\) with \(n=0\), so all the explanations in section \ref{sec:sig} apply for it. Due to absence of related events, we have a simpler case. We apply Bayes' rule again.
\begin{equation}
P(H_{\rm c, gw}|{\bf x}_{\rm gw}, {\bf X}_{\rm \nu})= \frac{P({\bf x}_{\rm gw}, {\bf X}_{\rm \nu}|H_{\rm c, gw})P(H_{\rm c, gw})}{P({\bf x}_{\rm gw}, {\bf X}_{\rm \nu})}
\end{equation}
\(P({\bf x}_{\rm gw}, {\bf X}_{\rm \nu})\) is omitted as it is through the paper. Then we separate the first term in the numerator due to independent detections as
\begin{equation}
P({\bf x}_{\rm gw}, {\bf X}_{\rm \nu}|H_{\rm c, gw})=P({\bf x}_{\rm gw}|H_{\rm c, gw})P({\bf X}_{\rm \nu}|H_{\rm c, gw})
\end{equation}
Then we obtain the GW part by marginalizing over parameters \({\boldsymbol \theta}\)
\begin{multline}
P({\bf x}_{\rm gw}|H_{\rm c, gw})=\\ \int P({\bf x}_{\rm gw}|{\boldsymbol \theta}, H_{\rm c, gw})P({\boldsymbol \theta}|H_{\rm c, gw})d{\boldsymbol \theta}
\end{multline}
with
\begin{equation}
    {\boldsymbol \theta}=\{t_s, r,{\bf \Omega},E_{\rm gw},E_{\nu}\}
\end{equation}
which are the same parameters defined in Section \ref{sec:sig}.
Now we move on analyzing \(P(H_{\rm c, \nu}|{\bf x}_{\rm gw}, {\bf X}_{\rm \nu})\). We first decompose it in sub-hypotheses for different number of signal neutrinos, denoted as \(H_{\rm c, \nu}^n\) for \(n\) signal neutrinos.
\begin{equation}
P(H_{\rm c, \nu}|{\bf x}_{\rm gw}, {\bf X}_{\rm \nu})=\sum_{n=1}^N P(H_{\rm c, \nu}^n|{\bf x}_{\rm gw}, {\bf X}_{\rm \nu})
\end{equation}
We again apply Bayes' rule to each term
\begin{equation}
\label{eq:1c}
P(H_{\rm c, \nu}^n|{\bf x}_{\rm gw}, {\bf X}_{\rm \nu})= \frac{P({\bf x}_{\rm gw}, {\bf X}_{\rm \nu}|H_{\rm c, \nu}^n)P(H_{\rm c, \nu}^n)}{P({\bf x}_{\rm gw}, {\bf X}_{\rm \nu})}
\end{equation}
Here the denominator again cancels out with the same denominator in the signal and null hypotheses. We again separate the first term in numerator due to independent detections as
\begin{equation}
P({\bf x}_{\rm gw}, {\bf X}_{\rm \nu}|H_{\rm c, \nu}^n)=P({\bf x}_{\rm gw}|H_{\rm c, \nu}^n)P({\bf X}_{\rm \nu}|H_{\rm c, \nu}^n)
\end{equation}
We further decompose the neutrino term by identifying the signal neutrinos with the set \(s\) which has \(n\) elements as in Section \ref{sec:sig}
\begin{widetext}
\begin{equation}
P({\bf X}_{\rm \nu}|H_{\rm c, \nu}^n)=\sum_{\{i,j,...\}}P({\bf X}_{\rm \nu}|H_{\rm c, \nu}^n,s=\{i,j,...\})P(s=\{i,j,...\}|H_{\rm c, \nu}^n)
\end{equation}
\end{widetext}
with
\begin{equation}
P(s=\{i,j,...\}|H_{\rm c, \nu}^n)={N \choose n}^{-1}
\end{equation}
We separate the signal and background neutrinos due to independence of different detections and drop the set \(s\) in notation as we did in Eq \ref{eq:1.6}
\begin{multline}
P({\bf X}_{\rm \nu}|H_{\rm c, \nu}^n,s=\{i,j,...\})=\\ P({\bf x}_{\nu i},{\bf x}_{\nu j},...|H_{\rm c, \nu}^n)P({\bf X}_{\rm \nu}^{i,j,...}|H_{\rm 0})
\end{multline}
The second term was explained in Section \ref{sec:null}. In order to obtain the first term we marginalize over parameters \({\boldsymbol \theta}\)
\begin{equation}
\label{eq:2c}
\begin{split}
P({\bf x}_{\nu i},{\bf x}_{\nu j}...|H_{\rm c, \nu}^n)= \\
\int P({\bf x}_{\nu i},{\bf x}_{\nu j},...|{\boldsymbol \theta}, H_{\rm c, \nu}^n)P({\boldsymbol \theta}|H_{\rm c, \nu}^n)d{\boldsymbol \theta}
\end{split}
\end{equation}
We again split each neutrino as we did in Section \ref{sec:sig}
\begin{equation}
P({\bf x}_{\nu i},{\bf x}_{\nu j},...|{\boldsymbol \theta}, H_{\rm c, \nu}^n)=P({\bf x}_{\nu i}|{\boldsymbol \theta}, H_{\rm c, \nu}^n)P({\bf x}_{\nu j}|{\boldsymbol \theta}, H_{\rm c, \nu}^n)...
\end{equation}
\subsection{Parameter and hypothesis priors (${\boldsymbol {H_{\rm c, gw}}}$)}
For \(P({\boldsymbol \theta}|H_{\rm c, gw})\) we use Bayes' rule.
\begin{equation}
P({\boldsymbol \theta}|H_{\rm c, gw})=\frac{P(H_{\rm c, gw}|{\boldsymbol \theta})P({\boldsymbol \theta})}{P(H_{\rm c, gw})}
\end{equation}
with \textcolor{black}{\(P({\boldsymbol \theta})\) in Eq. \ref{eq:parameterprior}.}

We write the probability density for detecting a GW event but not a neutrino as
\begin{multline}
P_{det}^{\rm c, gw}({\boldsymbol \theta}, \alpha)=(\alpha \mbox{Poiss}(0,\langle n_{\nu}(E_{\nu},r,{\bf \Omega})\rangle)+(1-\alpha))\\ \begin{cases} 1 &  r\leq r_0(E_{\rm gw})\Bar{f}_A({\bf \Omega},t_s) \\ 0 & \text{otherwise}
\end{cases}
\end{multline}
with \(\alpha\) being the ratio of total multi-messenger event rate to total astrophysical GW event rate. Then
\begin{equation}
P(H_{\rm c, gw}|{\boldsymbol \theta})=\frac{1}{N_2}T_{obs}\dot{n}_{\rm gw}P_{det}^{\rm c,gw}({\boldsymbol \theta},\alpha)
\end{equation}
with \(\dot{n}_{\rm gw}\) being the total astrophysical GW event rate in the whole universe.
\subsection{Parameter and hypothesis priors (${\boldsymbol {H_{\rm c, \nu}}}$)}
For \(P({\boldsymbol \theta}|H_{\rm c, \nu}^n)\) we use Bayes' rule.
\begin{equation}
P({\boldsymbol \theta}|H_{\rm c, \nu}^n)=\frac{P(H_{\rm c, \nu}^n|{\boldsymbol \theta})P({\boldsymbol \theta})}{P(H_{\rm c, \nu}^n)}
\end{equation}
with \textcolor{black}{\(P({\boldsymbol \theta})\) in Eq. \ref{eq:parameterprior}.}

We write the probability density for detecting \(n\) neutrinos but not a GW event as
\begin{widetext}
\begin{equation}
P_{det}^{\rm c, \nu,n}({\boldsymbol \theta}, \beta)= \begin{cases} (1-\beta) &  r\leq r_0(E_{\rm gw})\bar{f}_A({\bf \Omega},t_s) \\ \beta \mbox{Poiss}(n,\langle n_{\nu}(E_{\nu},r,{\bf \Omega})\rangle)+(1-\beta) & \text{otherwise}
\end{cases}
\end{equation}
\end{widetext}
with \(\beta\) being the ratio of total multi-messenger event rate to total astrophysical neutrino event rate. Then
\textcolor{black}{
\begin{equation}
P(H_{\rm c, \nu}^n|{\boldsymbol \theta})=\frac{1}{N_2}T_{obs}^2R_{\rm bg,gw}\dot{n}_{\rm \nu}P_{det}^{\rm c,\nu, n}({\boldsymbol \theta},\beta)
\end{equation}}
with \(\dot{n}_{\nu}\) being the total astrophysical neutrino event rate in the universe.
\textcolor{black}{
\subsection{Remaining terms}
\begin{itemize}
\item The term $P({\bf x}_{\rm gw}|{\boldsymbol \theta}, H_{\rm c}^{\rm gw})$ is equal to the same term for our signal hypothesis, i.e. $P({\bf x}_{\rm gw}|{\boldsymbol \theta}, H_{\rm s}^n)$ (see Eq. \ref{eq:gwsignal}), since in both cases there is a detected astrophysical gravitational wave signal.
\item The term $P({\bf x}_{\nu i}|{\boldsymbol \theta}, H_{\rm c, gw})$ is equal to the same term for our null hypothesis, i.e. $P({\bf x}_{\nu i}|H_{\rm 0})$ (see Eq. \ref{eq:neutrinoH0full}), since in both cases there is a background neutrino event, and neither term depends on the GW signal.
\item The term $P({\bf x}_{\rm gw}|H_{\rm c, \nu}^n)$ is equal to the same term for our background hypothesis, i.e. $P({\bf x}_{\rm gw}|H_{\rm 0})$ (see Eq. \ref{eq:gwsignalH0}), since in both cases there is a GW false detection from the background, and neither term depends on the neutrino signal.
\item The term $P({\bf x}_{\nu i}|{\boldsymbol \theta}, H_{\rm c, \nu}^n)$ is equal to the same term for our signal hypothesis, i.e. $P({\bf x}_{\nu i}|{\boldsymbol \theta}, H_{\rm s}^n)$ (see Eq. \ref{eq:neutrinosignal}), since in both cases there is a detected astrophysical neutrino, and neither term depends on the GW signal.
\end{itemize}
Combination of probabilities are done similarly to the signal and null hypothesis cases.}

\section{Odds ratio}
\label{sec:odds}
We test our signal hypothesis using odds ratios. We compare our signal hypothesis against both null and coincident hypotheses as in Eq. \ref{eq:or}.

It should be noted that this end result doesn't depend on \(T_{obs}\) which is a quantity fixed by humans' decisions and expected not to affect the significance of any astrophysical event. In addition there is no terms with explicit \(N\) dependency as expected since it could be arbitrarily large due to linear dependence on \(T_{obs}\). However, overall there is still a dependency due to the maximum possible signal neutrino count. In other words, \(N\) dependency of all terms got cancelled; but the number of terms depend on it.

This comparison will be applicable both for (i) GW and neutrino candidates that are not independently established detections, and for (ii) detections that are already confirmed through one channel. For the prior case, the first term in the denominator will be relevant, while in the latter case it will be the second term.

Although the odds ratio can be converted to a Bayesian probability for having a signal given the observations, it will be dependent on the parameter priors and the event rate densities which can be very uncertain. Therefore the odds ratio can be used as a test statistic. We empirically characterize the required threshold values based on background data and simulations for frequentist significances, similarly to \cite{Albert_2019}. For the searches with confirmed GW detections, the simulations consist of randomly paired simulated GWs and background neutrinos from previous detections. The number of background neutrinos in the time window around a GW is determined by Poisson distribution whose mean is the actual background neutrino rate times the duration of the time window. For searches with non-established GW detections, namely subthreshold searches, besides the GW and neutrino pairs for the previous case, there are pairs of time shifted background GW detections, which are acquired from GW data analysis pipelines, and background neutrinos, and pairs of background GW detections and signal neutrinos as well. All pairs are mixed in proportion to their estimated rates. These background comparisons allow us to determine a false alarm probability, namely p-value or significance, for the given events, which can be reported to initiate electromagnetic follow-up observations.

\textcolor{black}{During the O3 public alert search for coincident GW and HENs the following parameter values have been used: \(t_{gw}^+=-t_{gw}^-=t_{\nu}^+=-t_{\nu}^-=250s,\ \dot{n}_{\rm gw+\nu}=\dot{n}_{\rm gw}=(4\pi r_{max}^3/3)1500\ Gpc^{-3}year^{-1},\ \alpha=1,\ R_{\rm \nu,bg}=6.4\times 10^{-3} Hz,\ E_{\nu}^+=10^{51}\ ergs,\ E_{\nu}^-=10^{46}\ ergs,\ \Delta_{\delta}=2.5^\circ,\ \Delta_{\epsilon}(\epsilon_{\nu})=0.3\times \epsilon_{\nu}\). Furthermore since it is a public alert search for verified GW detections. we set \(P(H_{\rm 0})=P(H_{\rm c, \nu})=0\).}

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

We presented a search algorithm for common sources of GWs and high-energy neutrinos based on Bayesian hypothesis testing. This algorithm upgrades the method of Baret et al. \cite{2012PhRvD..85j3004B} that was used in most prior joint searches. The main advantages of the new method are that (i) it incorporates astrophysical priors about the source that help differentiate between signal and background, while being largely independent of the specific astrophysical model in consideration; (ii) it incorporates a more realistic model of the detector background, for example by taking into account the direction dependent background rate and energy distribution. These detector properties are straightforward to establish empirically, and the method presents a straightforward way to incorporate them as priors.

In the presentation of the method, we made simplifications that make the algorithm easier to implement and can make the computation faster. As an example, we assumed that all GW and neutrino sources emit the same energy. It will be useful to study how these simplifications affect the sensitivity of the search, and how much model dependence they introduce. This will be carried out in a future work.

\section*{Acknowledgments}
The authors are grateful for the useful feedback of Marek Szczepanczyk, Thomas Dent, Xilong Fan, and the IceCube Collaboration. The article has been approved for publication by the LIGO Scientific Collaboration under document number LIGO-P1800303. The authors thank the University of Florida and Columbia University in the City of New York for their generous support.
The Columbia Experimental Gravity group is grateful for the generous support of the National Science Foundation under grant PHY-1708028. \textcolor{black}{DV is grateful to the Ph.D. grant of the Fulbright foreign student program.}