# Eclipse Timing Modeling of Three Post-Common Envelope Binaries: Hybrid
Solutions

Xinyu Mai,1 Robert L. Mutel1
1Department of Physics and Astronomy, University of Iowa, Iowa City, IA 52242,
USA
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We report 90 new observations of three post common envelope binaries at
primary eclipse spanning between December 2018 to February 2022. We combine
recent primary eclipse timing observations with previously published values to
search for substellar circumbinary components consistent with timing
variations from a linear ephemeris. We used a least-squares minimization
fitting algorithm weighted by a Hill orbit stability function, followed by
Bayesian inference, to determine best-fit orbital parameters and associated
uncertainties. For HS2231+2441, we find that the timing data are consistent
with a constant period and that there is no evidence to suggest orbiting
components. For HS0705+6700, we find both one and two-component solutions
stable for at least 10 Myr. For HW Vir, we find three and four-component
solutions that fit the timing data reasonably well, but are unstable on short
timescales, and therefore highly improbable. Conversely, solutions calculated
using a Bayesian orbit stability prior result in a poor fit. The stable
solutions significantly deviate from the ensemble timing data in both systems.
We speculate that the observed timing variations for these systems, and very
possibly other sdB binaries, may result from a combination of substellar
component perturbations and an Applegate-Lanza mechanism.

###### keywords:

binaries: close - binaries: eclipsing - stars: individual: HW Vir,
HS0705+6700, HS2231+2441 - planetary systems

††pubyear: 2021††pagerange: Eclipse Timing Modeling of Three Post-Common
Envelope Binaries: Hybrid Solutions–A

##  1 Introduction

Since the discovery of the first extra-solar planets nearly 30 years ago
(Wolszczan & Frail, 1992; Mayor & Queloz, 1995), more than 4,000 exoplanetary
systems have been discovered, mostly orbiting late-type single stars (Martin,
2018). However, more than 100 planet-hosting multiple-star systems have also
been found (Schwarz et al., 2016), perhaps not surprisingly since a
substantial fraction of all stars are in multiple systems (Tokovinin, 2014).

An intriguing subset of these planet-hosting binaries are evolved short-period
binaries. In these systems, the more massive star has evolved off the main
sequence and lost nearly all of its hydrogen envelope during the red-giant
phase, filling its Roche lobe and generating unstable mass transfer to the
secondary star (Icko & Livio, 1993). Since the timescale for mass transfer is
much shorter than the thermal timescale of the secondary, a common envelope
forms surrounding both components. The orbiting stars experience drag forces,
transferring angular momentum to the envelope and causing it to expand and
eventually be expelled (Webbink, 1984; Zorotovic et al., 2010; Ivanova et al.,
2013). This leaves a hot B-star or white dwarf with a low mass main-sequence
companion separated by a few solar radii, and consequent orbital period of a
few hours (Nebot Gómez-Morán et al., 2011; Heber, 2009). These systems are
referred to as post common-envelope binaries (PCEB).

Eclipsing PCEBs provide a sensitive probe of potential substellar companions,
since the eclipses are frequent and deep, and can be timed with high accuracy.
This allows one to infer the presence of orbiting substellar companions by
modeling the observed eclipse timing variations (ETV’s) caused by the reflex
motion of the binary’s barycenter due to the substellar companion[s].

To date, more than a dozen eclipsing PCEBs have been observed to exhibit
large, occasionally erratic changes in their orbital period (Heber, 2016).
Suppose circumbinary substellar components cause these timing variations. In
that case, they must have either survived the common envelope phase or formed
from the common envelope debris and therefore be much younger than the parent
stars. This ‘second generation’ scenario is considered more likely for PCEBs
by some authors (Zorotovic & Schreiber, 2013; Schleicher & Dreizler, 2014;
Schleicher et al., 2015), but this hypothesis is still controversial. For
example, Bear & Soker (2014) argue that too large a fraction of the common
envelope mass would have to form into substellar companions for this to be
viable.

Detailed modeling of eclipse timing variations (ETV) has resulted in the
claimed detection of multiple circumbinary companions, either brown dwarfs or
massive planets, in at least ten PCEBs (Marsh, 2018). However, these claims
are not universally accepted: Several of the proposed planetary orbits have
been shown to be dynamically unstable over timescales much shorter than would
be consistent with even the ‘second generation’ planetary formation hypothesis
(Horner et al., 2012, 2014; Wittenmyer et al., 2013).

Alternative explanations for observed eclipse timing variation in binary stars
have been considered, including apsidal precession (Shakura, 1985; Parsons et
al., 2013) and changes in the mass quadruple moment of the secondary star
caused by magnetic activity, a.k.a. the Applegate mechanism (Applegate &
Patterson, 1987; Applegate, 1992; Lanza et al., 1998; Lanza, 2006; Beuermann
et al., 2012b; Lanza, 2020). Apsidal precession, which results from a small
non-zero eccentricity of the binary orbit, appears unlikely, since it would
require eccentricities of order $e\sim 0.01$ to explain observed PCEB eclipse
timing variations ($\sim$100 sec). These could be ruled out by observing
primary to secondary eclipse timing asymmetries, which would be of the same
order as the observed ETVs, and would have a periodic time shift of the
secondary eclipse with respect to phase 0.5. In one case (HW Vir), Baran et
al. (2018) reported no significant deviation from phase 0.5 for secondary
eclipse times observed over 70 days. Furthermore, ETV’s from PCEB binaries are
typically aperiodic, whereas apsidal precession would result in smooth timing
variations over an orbital period.

The Applegate mechanism has been considered to explain PCEB ETV’s by several
authors (e.g. Qian et al. (2008); Lee et al. (2009); Wittenmyer et al. (2012);
Nasiroglu et al. (2017). However, detailed stellar models of the angular
momentum exchange between a finite shell and the core of the secondary star
show that the required energy to produce the quadrupole moment variations is
too high in most PCEB systems, including all three of those described here
(Völschow, M. et al., 2016; Navarrete, F. H. et al., 2018; Navarrete et al.,
2019). However, recently (Lanza, 2020) has proposed a new model based on the
angular momentum exchange between the spin of the magnetically active
secondary and the binary orbital motion. This mechanism may overcome the
energy limitations of the original Applegate mechanism, and could contribute
to the observed eclipse timing variations.

In this paper, we report new primary eclipse timing observations of three PCEB
systems, combining them with previously published observations to solve for
possible substellar companions that can account for the observed timing
variations. We also analyze the orbital stability of the proposed systems
using numerical simulation of orbital dynamics and testing for chaotic
behavior over long timescales.

This procedure could be considered fraught and even unproductive, since it is
typically degenerate (fitted model parameters are typically highly correlated)
and almost all previously proposed planetary models have failed to correctly
predict subsequent ETV observations. Indeed, as Marsh (2018) recently pointed
out, there is no independent evidence for any substellar companion deduced
from ETV variations in PCEBs. Nevertheless, the ETV variations are clearly
real and all alternate explanations have so far failed with compelling
physical arguments. Hence, both Occam’s razor and the scientific method
encourage us to assume the most plausible hypothesis, find physical models
that are in agreement with the observations, and most importantly, provide
predictions that future observations can test.

###  1.1 HW Virginis

HW Vir (V = 10.7, P = 2.8 hr) is the prototype of the eponymous HW Vir class
of eclipsing subdwarf B stars with low-mass dwarf companions. Kilkenny et al.
(1994) first discovered variations in the orbital period and reported a period
decrease over nine years. Several authors have carried out subsequent ETV
studies of HW Vir (Çakirli & Devlen, 1999; Kiss et al., 2000; Kilkenny et al.,
2000; Kilkenny et al., 2003; İbanoǧlu et al., 2004; Lee et al., 2009; Qian et
al., 2010; Beuermann et al., 2012b; Skelton, 2012; Lohr et al., 2014; Baran et
al., 2018; Esmer et al., 2021; Brown-Sevilla et al., 2021). These papers
invoked one or more substellar companions as the primary cause of the observed
ETVs. However, as new data became available, these models inevitably needed to
be modified.

After ruling out mass transfer, magnetic braking, and apsidal motion, Çakirli
& Devlen (1999) concluded that a third body with mass 23 MJup and an orbital
period of 19 yrs is the most probable cause. İbanoǧlu et al. (2004) also
proposed a brown dwarf companion, with a minimum mass of 23 MJup and orbital
period of 18.8 yrs. Lee et al. (2009) proposed two companions with periods
15.8 and 9.1 yrs and masses of 19.2 MJup and 8.5 MJup respectively. Beuermann
et al. (2012a) added more ETV data, and showed that they deviated
significantly from the solution proposed by Lee et al. (2009). They also
showed that the solution proposed by Lee et al. is dynamically unstable. They
proposed a new solution consisting of two brown dwarfs (14 MJup \+ 65 MJup )
that was stable over more than 10 Myr. However, Baran et al. (2018) showed
that the solution of Beuermann et al. does not fit more recent O-C data.

Horner et al. (2012) independently made a stability analysis of the Lee et al.
solution, and also found it highly unstable. Indeed, they searched a
significant volume of substellar component parameter space, and could not find
any stable solutions, concluding "If the [HW Vir ] system does host
exoplanets, they must move on orbits differing greatly from those previously
proposed". This failure to find a stable solution is reinforced by the recent
studies of Esmer et al. (2021) and Brown-Sevilla et al. (2021) using all
previously published O-C observations along with new times of minima up to
2019.4. They tried a number of multi-planet models but also failed to find a
stable solution.

###  1.2 HS0705+6700

Drechsel et al. (2001) first confirmed HS0705+6700 (V470 Cam, V= 14.7, P= 2.3
hr) as a post common envelope sdB binary. They fit their eclipsing timing
observations with a linear ephemeris. Qian et al. (2009) reported a cyclic
change in orbital period and proposed a substellar companion with orbital
period of 7.15 years and mass of 40 MJup is responsible for the observed ETV.
Beuermann et al. (2012b) observed an additional 28 epochs, with a similar
solution: A single brown dwarf with mass 31 MJup and orbital period of 8.4
years.

In 2013, Qian et al. (2013) fit a single brown dwarf model with a period 8.87
years and mass 32 MJup , along with a positive quadratic term
$dP/dt=+9.8\times 10^{-9}\textrm{\ days/year}$. They ruled out mass transfer
and angular momentum loss due to gravitational radiation and magnetic braking
as explanations of this period increase. Pulley et al. (2015) showed that the
more recent O-C timing data began to diverge from the Qian model. Later in
2018, they confirmed this departure by adding 65 more times of minima (Pulley
et al., 2018).

Bogensberger et al. (2017) proposed a new one planet model with a longer
orbital 11.8 year period, eccentricity 0.38, and no quadratic term. Recently,
Sale et al. (2020) added 20 additional ETV observations through mid-2018. They
found a two component solution that fit all published O-C observations, but
unfortunately it was highly unstable over a short timescale (1,000 yr) and
hence not physically plausible.

###  1.3 HS2231+2441

HS2231+2441 (V = 14.2, P = 2.7 hr) is an eclipsing white dwarf + brown dwarf
binary (Almeida et al., 2017) first reported by Østensen et al. (2007). They
fit two years of primary eclipse timing observations with a linear (constant
period) ephemeris. Qian et al. (2010), using more extensive timing data,
reported that the orbital period shows both a cyclic change and a period
derivative. They proposed a 14 MJup component to explain the cyclic variation
and magnetic braking to explain the period decrease. However, as more timing
observations became available, Lohr et al. (2014) found that a linear function
might provide a better fit to the O-C curve. Almeida et al. (2017) reported
eight more eclipse times and updated the linear ephemeris. Most recently,
Pulley et al. (2018) extended the timing database with 26 more timing minima.
They compared the quadratic ephemeris by Qian et al. (2010) and linear
ephemeris proposed by Lohr et al. (2014) concluded there is no statistically
significant difference between the two.

Table 1: Observing Summary Target | V | $N_{obs}$ | t(sec) | Filter | Epochs
---|---|---|---|---|---
HS2231+2441 | 14.2 | 26 | 40 | Luminance | 2018.9 – 2020.8
HS0705+6700 | 14.7 | 40 | 40 | Luminance | 2019.1 – 2022.2
HW Vir | 10.7 | 24 | 15 | Sloan r′ | 2018.9 – 2022.1
Table 2: Mid-Eclipse Times (BJD)

[b]  HS2231+2441 HS0705+6700 HW Vir 2458471.69752(10) 2458500.71590(3)
2459143.84436(5) 2458472.00981(2) 2458473.57744(9) 2458504.82867(24)
2459144.89662(6) 2458473.99404(4) 2458617.89415(8) 2458505.59376(7)
2459145.85304(4) 2458607.75446(2) 2458650.84983(9) 2458512.67178(8)
2459146.80949(6) 2458612.77339(2) 2458991.90258(10) 2458565.66003(3)
2459147.86162(4) 2458621.76088(2) 2458992.89767(7) 2458589.66736(5)
2459149.87025(6) 2458622.69461(3) 2458998.86963(8) 2458619.70026(7)
2459151.87879(7) 2458641.019582(41)* 2459019.88131(6) 2458847.33965(5) *
2459152.83526(7) 2458657.71049(3) 2459020.87669(7) 2458985.74066(17)
2459153.79168(6) 2458981.72376(3) 2459021.87188(6) 2458987.74896(6)
2459156.85240(5) 2458981.84055(4) 2459022.86704(9) 2458994.73099(8)
2459157.80886(8) 2458990.71110(2) 2459023.86254(8) 2459122.89784(5)
2459158.86111(5) 2458992.69525(3) 2459024.85779(8) 2459126.91491(6)
2459266.65475(11) 2458996.78043(2) 2459026.84836(6) 2459127.87140(6)
2459619.78262(3) 2458997.71417(2) 2459138.65276(7) 2459128.92344(7)
2458998.76469(2) 2459139.64791(7) 2459129.88002(6) 2459003.66701(2)
2459140.64330(7) 2459131.88864(5) 2459004.71749(2) 2459141.63860(4)
2459132.84512(5) 2459020.70807(3) 2459142.63383(5) 2459133.89723(4)
2459022.69231(4) 2459143.62915(4) 2459134.85363(6) 2459024.67650(3)
2459144.62452(6) 2459136.86217(5) 2459263.83482(3) 2459145.61985(7)
2459138.87083(7) 2459265.81903(3) 2459146.61505(9) 2459139.82732(8)
2459269.78746(3) 2459147.61035(7) 2459140.87940(6) 2459583.99637(1)
2459151.59148(8) 2459141.83589(5) 2459152.58676(6) 2459142.88799(6)

  * •

NOTE - $*$ Times of minima with asterisk are unpublished observations from
David Pulley, John Mallet and the Altair Group. $a.$ HS0705+6700 - Dec 29th,
2019. $b.$ HW Vir - June 6th, 2019.

Table 3: Binary star physical properties

[b]  Binary Type Mpri Msec $R_{pri}$ $R_{sec}$ $L_{sec}$ $P_{orb}$ [d] Refs
HS2231+2441 WD+BD 0.288 0.046 0.165 0.086 0.0009 0.1106 Almeida et al. (2017)
HS0705+6700 sdB+M4-5 0.483 0.134 0.230 0.186 0.0022 0.0956 Drechsel et al.
(2001) HW Virginis sdB+M4-5 0.485 0.142 0.183 0.175 0.0025 0.1167 Lee et al.
(2009)

  * •

NOTE - $a.$ All masses, radii, and luminosities are solar units. $b$. Almeida
et al. (2017) presents two mass solutions for HS2231: We have listed solution
2. $c$. Baran et al. (2018) derived a primary mass 0.26 M⊙ for HW Vir, but
argue that it is probably an underestimate.

##  2 Observations

Figure 1: Sample HW Vir light curve taken at the Iowa Robotic Observatory on
21 May 2020 from 04:06 — 06:01 BJD, with fitted Gaussian profile to determine
time of minima [red dashed line].

Our study is an independent analysis of previously published eclipse times and
new timing data obtained with the 0.5m telescope at the Iowa Robotic
observatory111http://astro.physics.uiowa.edu/iro, located at Winer
Observatory222https://winer.org/ in southeastern Arizona. We used a SBIG 6303e
CCD camera (2048x3072 x 12 micron pixels) with an field of view 28x18 arcmin.
We report 90 new times of minimum for three PCEB systems at primary eclipse
between December 2018 and February 2022 (Table 1). A luminance filter was used
for HS2231+2441 and HS0705+6700, and an exposure time of 40 sec was adopted. A
Sloan r′ filter with an exposure time of 15 sec was used for the brighter
object HW Vir.

After the normal CCD image calibration including dark frame subtraction and
flat fielding, we created light curves using airmass-corrected differential
photometry. Both HS 2231+2441 and HS0705+6700 are typical of HW Vir-like
eclipsing systems, their light curves show a strong reflection effect with a
very deep primary and a shallow secondary minimum. We fit a Gaussian profile
to the primary eclipse light curve using a nonlinear least-squares (Levenberg-
Marquardt) algorithm and determine the observed mid-eclipse times. A sample
light curve for HW Vir with associated Gaussian profile fit is shown in Fig.
1. We evaluated the uncertainty in time of eclipse minimum by a Monte Carlo
method consisting of a series synthetic eclipse light curves with randomly-
chosen times of minima and different exposure times. This procedure is
described in detail in Appendix 1. Table 2 lists the eclipse timing
observations and uncertainties for all three PCEBs.

##  3 Data Analysis

###  3.1 Data selection

Our eclipse timing dataset compromised 90 new time of primary eclipse minima
(Table 2) together with previously published times from the literature. Before
use, we made a few modifications.

  1. 1.

For HW Vir and HS2231+2441, we excluded a portion of archival SuperWASP times
of minima that were rejected by Lohr et al. (2014)’s outlier methodology and
by the size of their uncertainties. In particular, we did not include Lohr et
al. (2014)’s listed extra times of minima due to their precision, and
discarded timing data for which uncertainties are between 8 sec $\sim$40 sec.
We also omitted the individual superWASP times of minima for HS2231+2441,
since they had large ($\sigma>20s$) uncertainties.

  2. 2.

We also removed O-C eclipse timings obtained from the Mercator observation in
Baran et al. (2018) for HW Vir since they deviate from the overall O-C trend
line by 50 seconds. We also excluded outliers that had O-C values deviating by
more than 400 seconds from the overall O-C trend.

###  3.2 Light Travel Time Model

For a binary system with orbiting substellar components, the predicted primary
mid-eclipse time as a function of cycle number $E$ can be written as

$T(E)=T_{0}+P_{0}E+\alpha E^{2}+\sum_{i}\tau_{i}(E)$ (1)

where $T_{0}$ is the reference epoch, $P_{0}$ is the linear eclipse period,
where $\tau_{i}$ are the time-dependent contributions from each perturbing
body. and

$\alpha=\frac{1}{2}\frac{dP}{dt}P_{0},$ (2)

is a quadratic term can would result from non-planetary effects, such as
apsidal motion, mass exchange, magnetic braking, or gravitational radiation.
We used the usual formulation of Irwin (1952, 1959) to calculate the time
signature of perturbing bodies on eclipse timing as a function of each body’s
orbital properties.

###  3.3 Orbit stability

A physically credible model must not only fit all the current and previously
published ETV datasets, but it is also subject to the prior constraint that
the substellar components must be stable to orbital instabilities over
timescales at least comparable to the lifetime of the PCEB phase ($\sim$ 100
Myr) (Han et al., 2002; Han et al., 2003; Zorotovic et al., 2011). We used the
Hill stability criterion (Gladman, 1993; Chambers et al., 1996) to predict
whether any two components in orbit about a massive central body will avoid
mutual close encounters over many orbits.

For circular orbits, stability occurs if the semi-major axis difference
between pairs of components exceeds a critical distance $\Delta$ given by
(Chambers et al., 1996),

$\Delta\simeq 2\sqrt{3}\cdot\ R_{H}$ (3)

where $R_{H}$ is the Hill radius,

$R_{H}=\frac{a_{1}+a_{2}}{2}\left[\frac{(m_{1}+m_{2})}{3M}\right]^{\frac{1}{3}}$
(4)

where $m_{1},m_{2}$ and $a_{1},a_{2}$ are the masses and semi-major axes of
the components respectively, and $M$ is the mass of the central body.

For non-circular orbits, a more general stability criterion is (Hasegawa &
Nakazawa, 1990)

$\Delta^{2}=\left[12+\frac{4}{3}K^{2}e^{2}\right]\ R_{H}^{2}$ (5)

where

$K=\left[\frac{(m_{1}+m_{2})}{3M}\right]^{\frac{1}{3}},$ (6)

and $e$ is the larger eccentricity. For low eccentricity orbits and $K\ll 1$,
we can expand the second term to obtain

$\Delta=2\sqrt{3}\ R_{H}\ \cdot\left(1+\frac{K^{2}e^{2}}{18}\right).$ (7)

The second term in parentheses is much less than 1 for the components
considered here ($K\ll 1$) so the simpler expression equation 3 was used in
our calculations.

We can define a dimensionless parameter $\beta$ that characterizes each planet
pair spacing in terms of the normalized mutual Hill radii as (cf. Smullen et
al., 2016)

$\beta=\frac{\left|a_{2}-a_{1}\right|}{\Delta}$ (8)

Stable component pairs correspond to $\beta>1$ while unstable pairs have
$\beta<1$.

For systems with more than two substellar components, instability can occur
even if all components lie outside the Hill radius of their nearest companion,
although the timescale for close encounters increases as the initial component
separations become larger (Chambers et al., 1996). Hence, we used the Hill
parameter as an optimization prior but ran numerical orbit simulations to test
orbit stability after optimization, as described below.

###  3.4 Model Fitting Scheme

The model fitting scheme comprised four steps (Fig. 3):

  1. 1.

Create an initial model of between one and four sub-stellar components, with
multi-component models spaced at least one Hill radius apart (i.e., $\beta>1$
for all pairs),

  2. 2.

Search for the minimum of an objective function proportional to the product of
the summed $\chi^{2}$ of the model vs. O-C time series and an ansatz ’penalty’
function that sharply increases when the smallest Hill parameter $\beta$
becomes less than 1.0 for any component pair,

  3. 3.

Use Bayesian inference to solve for the final fitted parameters using a Markov
chain Monte Carlo (MCMC) sampling, a weighted least-squares objective
function, and an orbit stability prior ($\beta>1$),

  4. 4.

Run orbit stability tests using both the N-body orbital integration package
REBOUND (Rein & Liu, 2012; Rein & Spiegel, 2015) and the chaotic indicator
program MEGNO (Cincotta & Simó, 2000; Hinse et al., 2010).

We first chose the number of substellar companions to model. Examination of
the O-C plot for both HS0705+6700 and HW Virginis indicates that a large
number of components, although certainly possible, is not justified by the
data, since the model parameters are already highly degenerate for few
component models, as we will see. Hence we restrict our model selection to one
through four component models. Each component is characterized by five
parameters: mass, semi-major axis, eccentricity, longitude of periapsis
$\tilde{\omega}$, and time of periastron passage . We restricted the model to
co-planar orbits for the sake of computational simplicity. We also added a
period derivative parameter to account for the possible contribution of non-
orbiting components, such as apsidal motion or gravitational radiation.

We started model fitting using a downhill simplex (Nelder-Mead) minimization
(Python library LMFIT333https://lmfit.github.io/lmfit-py to fit the O-C curve.
However, instead of minimizing the sum-squared differences between the O-C
data and the model, we incorporate the orbit stability constraint by
minimizing the product,

$\zeta\cdot\chi^{2}$ (9)

Where we define an empirical Hill stability function

$\zeta=0.8\ \rm{tanh}\left(2\beta_{min}-0.2\right)^{-\alpha}$ (10)

where the Hill parameter $\beta_{min}$ is the smallest normalized Hill
parameter (equation 8) of all component pairs. This stability function
asymptotes to unity for $\beta>1$ but increases sharply for $\beta<1$ (see
Fig.2).

Figure 2: Hill penalty function used in prior function (see text). (Eqn. 10).
Figure 3: Flow chart of data analysis.

The resulting model was used as a starting point for Bayesian inference
modeling. We used the EMCEE (Foreman-Mackey et al., 2013) Python
implementation444https://github.com/dfm/emcee with 4096 walkers and 4000
iterations. The likelihood function was summed-squared difference between the
model and O-C data weighted by the data uncertainties. We used three priors:
(i) stability parameter $\beta>1.05$, (ii) all masses and semi-major axes
positive-definite, and (iii) eccentricities restricted to the range $0<e<1$.

Using a model fitting scheme with a Bayesian methodology as the final step has
two important advantages over the more commonly-used frequentist least-squares
minimization approach. First, the orbit stability constraint (Hill parameter)
can be incorporated naturally as a Bayesian prior without the ad hoc ansatz
function we used in the LMFIT minimization. Second, the Bayesian posterior
distribution provides a more comprehensive description of the parameter
uncertainties and correlations.

For each model, we tested dynamical stability by numerically integrating the
component orbits using WHFAST, a Wisdom-Holman symplectic integrator (Rein &
Tamayo, 2015). We used a time step of 36.5 days and integrated for 10 Myr for
each model. In addition, we tested for chaotic behavior using the MEGNO
algorithm (Mean Exponential Growth factor of Nearby Orbits, Cincotta & Simó
(2000)) also with a time interval of 10 Myr. Both WHFAST and MEGNO were called
from the Python package REBOUND (Rein & Liu, 2012) which also provided
convenient tools for plotting orbits and time histories of orbital elements.

##  4 Results

###  4.1 HS2231+2441

After adding 26 additional primary eclipse timing observations (listed in
Table 2) to previously published observations, we fit the linear ephemeris of
Almeida et al. (2017) with a slightly changed orbital period,

${\rm BJD}_{min}=2455428.76185(10)+0.11058786(10)\cdot E$ (11)

The O-C diagram with all of the available data is shown in Figure 4.

We find that there is no statistically significant period variation for
HS22315+2441 spanning 16 years, from 2004 to 2020, but we cannot rule out the
variations smaller than 15 seconds. A linear fit without a orbiting planet is
consistent with the full data set.

Figure 4: HS2231+2441 O-C times of primary minima computed using the linear
ephemeris given in Equation 11. Observations shown as blue data points are
times of minima from Almeida et al. (2017), while green points are from Pulley
et al. (2018), and red points are new times of minima(Table 2).

###  4.2 HS0705+6700

Figure 5: Upper panel: HS0705+6700 O-C plot with one component model (blue
line) and one component model (green line) of Bogensberger et al. 2017
overlaid. The black data points are from the literature, while the red data
points are from Table 2. Lower panels: O-C differences between observed O-C
and model predictions. Figure 6: HS0705+6700 O-C plot of one component model.
The gray shaded area is the $\pm 1\sigma$ spread of the EMCEE posterior
samples distribution. The black data points are from the literature, while the
red data points are from Table 2. Figure 7: Joint posterior model parameter
probability distributions for 1-component model for HS0705+6700, as listed in
Table 4 respectively. Our model parameters converge to the median values,
shown in red, and contours 16%, 50%, and 84% quantiles. Note the time of
periastron passage is in shown in seconds.

Figure 5 Shows the O-C timing residuals versus epoch computed using linear
ephemeris,

${\rm BJD_{min}}=2451822.760509+0.095646609\cdot E$ (12)

along with our one component model (solid blue line) and the one-component
model of Bogensberger et al. 2017 (solid green line). Table 4 lists the model
parameters of our model.

We also found a best-fit two component model with masses 21.8 MJup and 1 MJup
and semi-major axes 4.9 AU and 37.3 AU, respectively. The two component model
was stable for at least $10^{7}$ years as verified by REBOUND and MEGNO
simulations, but the fit to the timing data was no better than the one
component model. Therefore, we do not discuss this model further.

It is clear that the previously published model of Bogensberger et al.
deviates strongly from the observed timing data starting shortly after the
publication of the model. As previously noted, this discouraging trend has
been seen for several other eclipsing sdB binary exoplanet detection. Also, we
did not plot the recent two component solution proposed by Sale et al. (2020)
since we could not determine their ephemeris or interpret orbital parameters
of their model from their paper.

An expanded view of the O-C fits for our model is shown in Fig. 6, along with
the $\pm 1\sigma$ posterior spread from the EMCEE solutions. Fig. 7 shows the
two-dimensional projections of the sample covariances for our one component
model. They demonstrate that our model have well-defined parameters with small
fractional uncertainties.

Table 4: HS0705+6700 One Component Model Parameter | Fitted Values | Unit
---|---|---
Inner Binary
T0 | 2451822.760509 |  | BJD
P0 | 0.095646609 |  | day
$\dot{P}$ | ${22.71_{-0.03}^{+0.02}}$ |  |  $10^{-12}$ s/s
Substellar Component
$P$ | $13.63_{-0.06}^{+0.05}$ |  | yr
$K$ | $74.86_{-0.20}^{+0.22}$ |  | sec
$asin(i)$ | $4.91_{-0.01}^{+0.01}$ |  | AU
$e$ | $0.65_{-0.01}^{+0.01}$ |  |
$\tilde{\omega}$ | $1.08_{-0.01}^{+0.01}$ |  | radian
Min. Mass | $20.45_{-0.10}^{+0.10}$ |  | MJup
$T$ | $2451081.52413_{-17.52}^{+24.26}$ |  | day

###  4.3 HW Virginis

Figure 8: HW Vir O-C times of primary minima computed using the linear
ephemeris given in Equation 13. Black points are archival times of minima from
the literature, red points are new times of minima from Table 2, and blue
lines are model predictions. (a) Best-fit 3-component model subject to the
constraint that the model is stable i.e., all pairwise component normalized
Hill parameters exceed 1.05. (b) Eccentricity and semi-major axes of all three
components over 10 Myr. This solution is stable but is a very poor fit
($\chi^{2}_{r}\approx 93$) to the timing data. (c) Best-fit 4-component model.
(d) Eccentricity and semi-major axes of all four components over 100 yr. The
four component model provides a much improved fit ($\chi^{2}_{r}\approx 19$)
but is highly unstable. We conjecture that the observed ETV’s may be result
from a combination of gravitational perturbations by substellar companions and
a Applegate-Lanza mechanism.

Figure 8 Shows the O-C timing residuals versus epoch using the linear
ephemeris

${\rm BJD_{min}}=2450280.28596+0.116719519\cdot E$ (13)

On the right of each panel are plots of eccentricity and semi-major axes of
all model components vs. epoch as computed by numerical integration using
REBOUND with time steps of 0.1 yr. The plots also show two model fits:

  * •

The upper panel shows a three component model computed using the minimization
scheme described in section 3, with a Bayesian prior $\beta>1.05$ for all
component pairs, where $\beta$ is the normalized Hill constraint. This
solution is stable for at least 10 Myr, but is a very poor fit to the timing
data ($\chi^{2}\approx 93$).

  * •

The lower panel shows a 4-component model, also computed using the
minimization scheme in section 3, but without the Hill stability prior. This
provides a much better fit ($\chi^{2}\approx 19$) but is highly unstable over
short time internals (one hundred years).

We do not list the orbital elements for either model since they are
implausible: They are either very poor fits to the timing data (3-component
stable model) or are highly unstable (four-component model). This unfortunate
scenario (acceptable O-C model fits are invariably unstable) was also found by
(Brown-Sevilla et al., 2021) who analyzed a series of best-fit six multi-
component models for HW Vir. In the next section we consider possible
alternatives to pure substellar component models.

##  5 Alternative explanations for eclipse timing variations

As mentioned previously, eclipse timing variations can be caused by mechanisms
other than gravitational perturbation by orbiting circumbinary masses. These
include interbinary mass transfer, apsidal motion, and changes in the moment
of inertia caused by magnetic activity on one of the components (Applegate
mechanism).

For sdB binaries, mass transfer is ruled out since they are detached systems.
Apsidal precession can result from (general relativistic) gravitational
radiation and by tidal effects caused by an asymmetric mass distribution in
one (or both) components. Lee et al. (2009) has shown that gravitational
radiation for the sdB binary HW Vir, and by extension other similar systems
including HS0705+6700, would result in time variations at about 2 dex smaller
than observed.

Ruling out apsidal effects is more problematic since they depend on the tidal
distortion of the internal mass distribution (via the internal structure
constant $k_{2}$, Sirotkin & Kim 2009) which is not well-characterized for
these systems. Almeida et al. (2020) recently considered tidal and rotation
apsidal effects on eclipse timing variations for the PCEB GK Vir, They found
that apsidal is less likely than substellar component model, but could not
rule it out as a contributing factor. In addition, apsidal effects will
produce strictly periodic signatures in O-C plots, which is not observed for
the binaries considered in this paper. Hence, we do not consider apsidal
effects further, but examine the possible contribution of Applegate-type
mechanisms.

Applegate (Applegate & Patterson, 1987; Applegate, 1992) proposed that orbital
period changes in close binaries may result from a variable quadrupole moment
produced by magnetic activity in the outer convection zone of one of the
stars. For sdB binaries, this is likely to be the low-mass secondary, which is
fully convective. Although we are not aware of any direct evidence for
magnetic fields in any sdB secondary star, there is ample evidence for strong
(kiloGauss) surface fields in many late-type (M0-M5) main sequence single
stars similar to the secondaries in sdB binaries (e.g., Kochukhov, 2020, and
references therein). Also, tidal locking will ensure that the rotation of the
secondary is very rapid, enhancing the prospects for a robust magnetic dynamo.

The original Applegate mechanism assumed a time-dependent gravitational
quadrupole moment modulated by the magnetic activity cycle in an outer thin
shell of the active star. An increasing quadrupole moment results in a
stronger gravitational field, driving a decreasing orbital radius and
consequent reduction of the binary period. However, it was subsequently found
that energy required to drive quadrupole moment variations of the secondary
star that could account for the ETV’s was inconsistent with observed
luminosity changes in many close binaries (Völschow, M. et al., 2016).

Several recent papers have analyzed modifications to the original Applegate
model using more complex 2-shell models, generalizing the magnetohydrodynamic
effects of internal magnetic fields, and considering different types of dynamo
models (e.g., Lanza et al., 1998; Lanza, 2006; Navarrete, F. H. et al., 2018;
Völschow, M. et al., 2018; Navarrete et al., 2019). Here we consider the
recent model of Lanza (2020) which is based on angular momentum exchange
between the spin of the active component and the orbital motion. Spin–orbit
coupling is produced by a non-axisymmetric component of the gravitational
quadrupole moment of the active star caused by an assumed persistent non-
axisymmetric internal magnetic field.

###  5.1 Hybrid models with substellar components and Applegate-type
modulation

As we have seen, the best-fit stable orbit model fits for both HS0705 and HW
Vir are not fully in agreement with the observed O-C data, with residual
average differences of order 10 sec and 20 sec over 10-20 year timescales for
HS0705 and HW Vir respectively (Figs. 5,8). We now consider the question of
whether the ETV variations may result from a combination of substellar
component perturbation and Applegate-Lanza effects.

The change in fractional angular velocity $\Omega$ of the secondary is related
to the observed fractional period changes by (Lanza, 2020, eqn. 56)

$\frac{\Delta\Omega}{\Omega}=-\frac{ma^{2}}{3I_{s}}\frac{\Delta P}{P}$ (14)

Where $m$ in the reduced mass of the binary, $a$ is the semi-major axis of the
binary, $I_{s}$ is the moment of inertia about the active (secondary) star
spin axis, and $\Delta P/P=\Delta t/T$ is the fractional change in the period,
or equivalently the O-C deviations $t$ over the timescale $T$ of the O-C data
. The consequent change in rotational energy is

$\Delta E_{rot}=I_{s}\ \Omega\ \Delta\Omega$ (15)

Combining these equations, the dependence on $I_{p}$ cancels, giving,

$\Delta E_{rot}=-\left[mr_{0}^{2}\Omega^{2}\right]\frac{\Delta t}{3T}$ (16)

For both HS0705+6700 and HW Vir, the average timing residuals of the
substellar companion model fits to the observed O-C data are both of order 10
sec over 10 yr ($\Delta t/3T\sim 10^{-8}$), while the bracketed term evaluates
to $3\cdot 10^{40}$ J, so the change in rotational energy is

$\Delta E_{rot}\approx\ -3\cdot 10^{32}\ \rm{J}$ (17)

We can compare this with the maximum energy available from luminosity changes
over the timescale for O-C variability in the secondary star. From Table 3,

$E_{max}=L\cdot t\sim 0.002\ L_{sun}\cdot 10\rm{\ yr}\sim 3\cdot 10^{32}\ J$
(18)

Remarkably, this is the same as the energy loss required to account for
residual changes in the O-C model fits, but would be only $\sim 10\%$ of the
entire O-C deviations from a constant period were due to an Applegate-Lanza
mechanism. This suggests that the observed O-C variations could be explained
by a combination of substellar companions accounting for the bulk of the
ETV’s, with a Lanza-Applegate mechanism contributing about 10% of the total.

###  5.2 Tidal Synchronization Timescale

Another relevant dynamical timescale is the tidal synchronization time i.e.,
the timescale in which tides on the secondary dissipate the kinetic energy of
asynchronous rotation (Lanza, 2020),

$t_{sync}=\frac{2Q}{9}\frac{M_{T}}{m_{p}^{2}}\frac{I_{s}}{\sqrt{GM_{T}R}}{\left(\frac{a}{R}\right)}^{9/2}$
(19)

where Q is a dimensionless tidal quality factor that characterizes the
efficiency of tidal energy dissipation in the active (secondary) component,
$M_{T}$ is the total mass of the binary system, $m_{p}$ is the primary mass,
$a$ is the binary semi-major axis, and $R$ is the active star radius.
Following Lanza (2020), we assume $Q\sim 10^{5}-10^{6}$ corresponding to
strong tidal coupling due to nearly synchronous rotation. For the moment of
inertia $I_{s}$, we use the numerical expression of Criss & Hofmeister (2015)
for a fully convective secondary star with a polytropic index n = 1.5,

$I_{s}=0.2\ m_{sec}R_{sec}^{2}$ (20)

Using stellar parameters list in Table 3, we find tidal synchronization
timescales $t_{s}$ = 4 – 40 yr and 10 – 100 yr for HS0705 and HW Vir
respectively. These times are comparable to the observed O-C variations, so
tidal synchronization may dictate the timescale for spin-orbit coupling rather
than magnetic activity timescales.

##  6 Summary and Conclusions

In this paper, we have reported new eclipse timing observations of three post-
common envelope binary systems. We combined these observations with previously
published timing data to investigate whether the eclipse timing variations
from a linear ephemeris could be explained by the gravitational influence of
orbiting substellar companions. We used models comprising multiple substellar
companions in coplanar orbits parallel to the observer’s line of sight,
adjusting their orbital parameters and masses to best fit the observed timing
variations.

For one of the three binaries (HS2231+2441), nearly all times of minima are
consistent with a linear ephemeris. Hence, there is no evidence for substellar
companions in this system, although they cannot be ruled out if the masses are
small. However for the remaining two binaries, there are large timing
variations (ETV) from the linear trend, indicating the presence of one or more
gravitationally perturbing effects. We focused on finding orbital parameters
and masses of companions that could account for the observed ETV’s by varying
the masses and orbital elements to best-fit the observed ETV time signature
using a least-squares minimization.

In order to ensure that the substellar components had stable orbits, the
minimization function included a Hill stability ‘penalty’ ansatz that
increases sharply when the minimum separation between adjacent components was
smaller than one Hill radius. We also used a Bayesian inference algorithm to
better characterize the uncertainties of the derived orbital parameters and
masses. In addition, we used the Hill stability criterion $H>1$ as a prior to
restrict orbital solutions only to those that were plausibly long-term stable
i.e., at least as long as the assumed time scale for common-envelope evolution
(Schreiber & Gänsicke, 2003; Ivanova et al., 2013).

Finally, we checked stability timescales for each model using a numerical
N-body symplectic time-step integrator (WHFAST) and an algorithm to detect the
onset of chaotic behavior (MEGNO). We ran both the numerical integrations and
chaos indicator time evolutions for $10^{7}$ years.

For HS0705+6700, we found one and two-component models that reasonably
reproduced the observed ETV time history and were stable on long time scales
($\approx 10^{7}$ yr). However, the fits were not statistically convincing,
with reduced chi-squares of the fitted models in the range 10–19, although
this may be partly a result of underestimated uncertainties in previously
published timing measurements (cf. Hinse et al., 2014).

For HW Vir, we found no stable solutions that could fit the ETV time
signature, although it was possible to find adequate four component fits,
albeit with highly unstable orbits. This has also been the conclusion of
several previous studies of HW Vir (e.g., Horner et al., 2012; Brown-Sevilla
et al., 2021) and although not definitive, this strongly suggests that
substellar companions cannot account for all of the ETV time history observed
in these (and likely other) PCEB binaries.

Therefore, we investigated whether some fraction of the observed variations
could be caused by an Applegate mechanism in which the secondary star’s
gravitational quadrupole moment is altered by magnetic activity. Although the
original Applegate mechanism was found to be energetically insufficient to
account for the observed timing variations in sdB binaries (e.g., Völschow, M.
et al. 2016), various revisions of the mechanism have been suggested that may
account for at least a fraction of the observed variations.

We applied the model of Lanza (2020), in which spin–orbit coupling is produced
by an active (secondary) star’s non-axisymmetric internal magnetic field
creates a non-axisymmetric gravitational quadrupole moment. Using this model,
we find that the predicted timing variations could account for perhaps 10% of
the observed O-C variability, although this result is uncertain, since it is
dependent on poorly constrained binary system parameters such as the tidal
energy dissipation and the internal mass distribution of the stars. Hence, we
posit that the observed O-C variations could result from a hybrid model in
which substellar companions account for a largest fraction of the timing
variations, but that an Appleget-Lanza mechanism contributes the remainder.

## Acknowledgements

This research made use of the LMFIT nonlinear least-square minimization Python
library, and Astropy, a community developed core Python package for Astronomy.
We also made use of emcee, an MIT licensed pure-Python implementation of the
affine invariant Markov chain Monte Carlo (MCMC) ensemble sampler (Goodman &
Weare, 2010). We also used the excellent Python module corner.py to visualize
multidimensional results of Markov Chain Monte Carlo simulations (Foreman-
Mackey, 2016). Numerical orbit integrations in this paper made use of the
REBOUND N-body code (Rein & Liu, 2012). The simulations were integrated using
WHFAST, a Wisdom-Holman symplectic integrator (Rein & Tamayo, 2015).

We are grateful to David Pulley, John Mallett and the Altair Group for their
contribution of times of minima for HS0705+6700 on Dec 29th, 2019, and times
of minima for HW Vir on June 6th, 2019. We are also grateful to Prof. Antonino
Lanza for several very insightful communications. The Iowa Robotic Observatory
is funded by the College of Liberal Arts at the University of Iowa.

## Data Availability

We have provided Table 2 as our new mid-eclipse timing data of three post-
common envelope binaries. The light curves available on request from the
corresponding author.

## References

  * Almeida et al. (2017) Almeida L. A., Damineli A., Rodrigues C. V., Pereira M. G., Jablonski F., 2017, MNRAS, 472, 3093
  * Almeida et al. (2020) Almeida L. A., Pereira E. S., Borges G. M., Damineli A., Michtchenko T. A., Viswanathan G. M., 2020, Monthly Notices of the Royal Astronomical Society, 497, 4022
  * Applegate (1992) Applegate J. H., 1992, ApJ, 385, 621
  * Applegate & Patterson (1987) Applegate J. H., Patterson J., 1987, ApJ, 322, L99
  * Baran et al. (2018) Baran A. S., et al., 2018, MNRAS, 481, 2721
  * Bear & Soker (2014) Bear E., Soker N., 2014, Monthly Notices of the Royal Astronomical Society, 444, 1698
  * Beuermann et al. (2012a) Beuermann K., et al., 2012a, A&A, 540, A8
  * Beuermann et al. (2012b) Beuermann K., Dreizler S., Hessman F. V., Deller J., 2012b, A&A, 543, A138
  * Bogensberger et al. (2017) Bogensberger D., Clarke F., Lynas-Gray A. E., 2017, Open Astronomy, 26, 134
  * Brown-Sevilla et al. (2021) Brown-Sevilla S. B., et al., 2021, Monthly Notices of the Royal Astronomical Society, 506, 2122
  * Carter et al. (2008) Carter J. A., Yee J. C., Eastman J., Gaudi B. S., Winn J. N., 2008, The Astrophysical Journal, 689, 499
  * Chambers et al. (1996) Chambers J., Wetherill G., Boss A., 1996, Icarus, 119, 261
  * Cincotta & Simó (2000) Cincotta P. M., Simó C., 2000, Astron. Astrophys. Suppl. Ser., 147, 205
  * Criss & Hofmeister (2015) Criss R. E., Hofmeister A. M., 2015, New Astron., 36, 26
  * Deeg (2021) Deeg H. J., 2021, Galaxies, 9
  * Deeg, H. J. (2015) Deeg, H. J. 2015, A&A, 578, A17
  * Drechsel et al. (2001) Drechsel H., et al., 2001, A&A, 379, 893
  * Esmer et al. (2021) Esmer E. M., Baştürk Ö., Hinse T. C., Selam S. O., Correia A. C. M., 2021, A&A, 648, A85
  * Foreman-Mackey (2016) Foreman-Mackey D., 2016, Journal of Open Source Software, 1, 24
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Gladman (1993) Gladman B., 1993, Icarus, 106, 247
  * Goodman & Weare (2010) Goodman J., Weare J., 2010, Communications in Applied Mathematics and Computational Science, 5, 65
  * Han et al. (2002) Han Z., Podsiadlowski P., Maxted P. F. L., Marsh T. R., Ivanova N., 2002, MNRAS, 336, 449
  * Han et al. (2003) Han Z., Podsiadlowski P., Maxted P. F. L., Marsh T. R., 2003, Monthly Notices of the Royal Astronomical Society, 341, 669
  * Hasegawa & Nakazawa (1990) Hasegawa M., Nakazawa K., 1990, A&A, 227, 619
  * Heber (2009) Heber U., 2009, Annual Review of Astronomy and Astrophysics, 47, 211
  * Heber (2016) Heber U., 2016, Publications of the Astronomical Society of the Pacific, 128, 082001
  * Hinse et al. (2010) Hinse T. C., Christou A. A., Alvarellos J. L. A., Goździewski K., 2010, MNRAS, 404, 837
  * Hinse et al. (2014) Hinse T. C., Horner J., Lee J. W., Wittenmyer R. A., Lee C.-U., Park J.-H., Marshall J. P., 2014, A&A, 565, A104
  * Horner et al. (2012) Horner J., Hinse T. C., Wittenmyer R. A., Marshall J. P., Tinney C. G., 2012, MNRAS, 427, 2812
  * Horner et al. (2014) Horner J., Wittenmyer R., Hinse T., Marshall J., Mustill A., 2014, Wobbling Ancient Binaries - Here Be Planets? (arXiv:1401.6742)
  * İbanoǧlu et al. (2004) İbanoǧlu C., Çakırlı Ö., Tas G., Evren S., 2004, A&A, 414, 1043
  * Icko & Livio (1993) Icko J. I., Livio M., 1993, Publications of the Astronomical Society of the Pacific, 105, 1373
  * Irwin (1952) Irwin J. B., 1952, ApJ, 116, 211
  * Irwin (1959) Irwin J. B., 1959, AJ, 64, 149
  * Ivanova et al. (2013) Ivanova N., et al., 2013, A&ARv, 21, 59
  * Kilkenny et al. (1994) Kilkenny D., Marang F., Menzies J. W., 1994, MNRAS, 267, 535
  * Kilkenny et al. (2000) Kilkenny D., Keuris S., Marang F., Roberts G., van Wyk F., Ogloza W., 2000, The Observatory, 120, 48
  * Kilkenny et al. (2003) Kilkenny D., van Wyk F., Marang F., 2003, The Observatory, 123, 31
  * Kiss et al. (2000) Kiss L. L., Csák B., Szatmáry K., Furész G., Sziládi K., 2000, A&A, 364, 199
  * Kochukhov (2020) Kochukhov O., 2020, The Astronomy and Astrophysics Review, 29, 1
  * Kwee & van Woerden (1956) Kwee K. K., van Woerden H., 1956, Bull. Astron. Inst. Netherlands, 12, 327
  * Lanza (2006) Lanza A. F., 2006, MNRAS, 369, 1773
  * Lanza (2020) Lanza A. F., 2020, MNRAS, 491, 1820
  * Lanza et al. (1998) Lanza A. F., Rodono M., Rosner R., 1998, MNRAS, 296, 893
  * Lee et al. (2009) Lee J. W., Kim S.-L., Kim C.-H., Koch R. H., Lee C.-U., Kim H.-I., Park J.-H., 2009, The Astronomical Journal, 137, 3181
  * Lohr et al. (2014) Lohr M. E., et al., 2014, A&A, 566, A128
  * Marsh (2018) Marsh T. R., 2018, Circumbinary Planets Around Evolved Stars. Springer International Publishing, Cham, pp 2731–2747, doi:10.1007/978-3-319-55333-7_96, https://doi.org/10.1007/978-3-319-55333-7_96
  * Martin (2018) Martin D. V., 2018, Populations of Planets in Multiple Star Systems. Springer International Publishing, Cham, pp 2035–2060, doi:10.1007/978-3-319-55333-7_156, https://doi.org/10.1007/978-3-319-55333-7_156
  * Mayor & Queloz (1995) Mayor M., Queloz D., 1995, Nature, 378, 355
  * Nasiroglu et al. (2017) Nasiroglu I., et al., 2017, AJ, 153, 137
  * Navarrete, F. H. et al. (2018) Navarrete, F. H. Schleicher, D. R. G. Zamponi Fuentealba, J. Völschow, M. 2018, A&A, 615, A81
  * Navarrete et al. (2019) Navarrete F. H., Schleicher D. R. G., KÃ¤pylÃ¤ P. J., Schober J., VÃ¶lschow M., Mennickent R. E., 2019, Monthly Notices of the Royal Astronomical Society, 491, 1043
  * Nebot Gómez-Morán et al. (2011) Nebot Gómez-Morán A., Gänsicke B. T., Schreiber M. R., Schwope A. D., 2011, in Schmidtobreick L., Schreiber M. R., Tappert C., eds, Astronomical Society of the Pacific Conference Series Vol. 447, Evolution of Compact Binaries. p. 187
  * Østensen et al. (2007) Østensen R., Oreiro R., Drechsel H., Heber U., Baran A., Pigulski A., 2007, in Napiwotzki R., Burleigh M. R., eds, Astronomical Society of the Pacific Conference Series Vol. 372, 15th European Workshop on White Dwarfs. p. 483
  * Parsons et al. (2013) Parsons S. G., et al., 2013, Monthly Notices of the Royal Astronomical Society: Letters, 438, L91
  * Pulley et al. (2015) Pulley D., Faillace G., Smith D., Watkins A., Owen C., 2015, Journal of the British Astronomical Association, 125, 284
  * Pulley et al. (2018) Pulley D., Faillace G., Smith D., Watkins A., von Harrach S., 2018, A&A, 611, A48
  * Qian et al. (2008) Qian S.-B., Dai Z., Zhu L.-Y., Liu L., He J.-J., Liao W.-P., Li a., 2008, The Astrophysical Journal Letters, 689, L49
  * Qian et al. (2009) Qian S.-B., et al., 2009, The Astrophysical Journal, 695, L163
  * Qian et al. (2010) Qian S. B., et al., 2010, Ap&SS, 329, 113
  * Qian et al. (2013) Qian S. B., et al., 2013, MNRAS, 436, 1408
  * Rein & Liu (2012) Rein H., Liu S. F., 2012, A&A, 537, A128
  * Rein & Spiegel (2015) Rein H., Spiegel D. S., 2015, MNRAS, 446, 1424
  * Rein & Tamayo (2015) Rein H., Tamayo D., 2015, MNRAS, 452, 376
  * Sale et al. (2020) Sale O., Bogensberger D., Clarke F., Lynas-Gray A. E., 2020, Monthly Notices of the Royal Astronomical Society, 499, 3071
  * Schleicher & Dreizler (2014) Schleicher D. R. G., Dreizler S., 2014, A&A, 563, A61
  * Schleicher et al. (2015) Schleicher D. R. G., Dreizler S., VÃ¶lschow M., Banerjee R., Hessman F. V., 2015, Astronomische Nachrichten, 336, 458
  * Schreiber & Gänsicke (2003) Schreiber M. R., Gänsicke B. T., 2003, A&A, 406, 305
  * Schwarz et al. (2016) Schwarz R., Funk B., Zechner R., BazsÃ³ Ã., 2016, Monthly Notices of the Royal Astronomical Society, 460, 3598
  * Shakura (1985) Shakura N. I., 1985, Soviet Astronomy Letters, 11, 224
  * Sirotkin & Kim (2009) Sirotkin F. V., Kim W.-T., 2009, The Astrophysical Journal, 698, 715
  * Skelton (2012) Skelton P. L., 2012, African Skies, 16, 132
  * Smullen et al. (2016) Smullen R. A., Kratter K. M., Shannon A., 2016, MNRAS, 461, 1288
  * Tokovinin (2014) Tokovinin A., 2014, The Astronomical Journal, 147, 87
  * Völschow, M. et al. (2016) Völschow, M. Schleicher, D. R. G. Perdelwitz, V. Banerjee, R. 2016, A&A, 587, A34
  * Völschow, M. et al. (2018) Völschow, M. Schleicher, D. R. G. Banerjee, R. Schmitt, J. H. M. M. 2018, A&A, 620, A42
  * Webbink (1984) Webbink R. F., 1984, ApJ, 277, 355
  * Wittenmyer et al. (2012) Wittenmyer R. A., Horner J., Marshall J. P., Butters O. W., Tinney C. G., 2012, Monthly Notices of the Royal Astronomical Society, 419, 3258
  * Wittenmyer et al. (2013) Wittenmyer R. A., Horner J., Marshall J. P., 2013, Monthly Notices of the Royal Astronomical Society, 431, 2150
  * Wolszczan & Frail (1992) Wolszczan A., Frail D. A., 1992, Nature, 355, 145
  * Zorotovic & Schreiber (2013) Zorotovic M., Schreiber M. R., 2013, A&A, 549, A95
  * Zorotovic et al. (2010) Zorotovic M., Schreiber, M. R. Gänsicke, B. T. Nebot Gómez-Morán, A. 2010, A&A, 520, A86
  * Zorotovic et al. (2011) Zorotovic M., Schreiber M. R., Gänsicke B. T., 2011, Astronomy & Astrophysics, 536, A42
  * Çakirli & Devlen (1999) Çakirli Ö., Devlen A., 1999, A&AS, 136, 27

##  Appendix A Time of minimum Uncertainty estimates

In this appendix, we evaluate the uncertainty in the time of minimum
determination of an eclipsing binary light curve constructed from a series of
images with relatively long exposure times and image download times. These
pertain especially to observations using relatively small-aperture telescopes,
including those reported in this paper.

A standard technique used by many previous authors (more than 1,100 citation!)
is the folded-symmetry algorithm described by Kwee & van Woerden (1956). This
technique determines the time of minimum by adjusting a trial midpoint time
and minimizing the summed squared magnitude difference between points equally
time-displaced on either side of the midpoint. The associated uncertainty is
determined by fitting a quadratic function and applying standard propagation
of errors to fitted minimum value. The technique relies on equally-space data
samples and does not account for interpolation errors when this is not the
case. Other more sophisticated methods (e.g., Carter et al., 2008; Deeg, H.
J., 2015; Deeg, 2021) have addressed eclipse timing uncertainty, but they
generally have been restricted to trapezoidal light curves, as observed in
exoplanet transits. This is not applicable to short-period binary systems,
whose eclipse light curves are much closer to Gaussian profiles (e.g., Baran
et al., 2018).

Figure 9: (a) HS0705+6700 primary eclipse 16 April 2019 with LSQ fit Gaussian
model.(b) Simulated light curve.

We evaluated time of minimum uncertainty using a Monte Carlo approach,
simulating observed eclipse light curves with a Gaussian profile whose width,
depth, and scatter closely resembles the observed light curve. Fig. 9 shows a
light curve of HS0705+6700 observed on with the Iowa Robotic telescope, along
with the corresponding simulated light curve. The synthetic light curve was
made by sampling a Gaussian function with 0.25 magnitude eclipse depth, and
300 sec half-width. We added Gaussian random noise to each sampled point, with
a standard deviation $\sigma=0.025\times 2.5^{\Delta m}$ where $\Delta m$ is
the differential magnitude difference to the out-of-eclipse level. The
synthetic light curve was sampled every 10 msec. We then binned the simulated
date into 40 second bins with 10 second time gap between bins, corresponding
the the observed image exposure time and download times for HS0705+6700. The
magnitude assigned to each image was the computed mean value of 10 msec
samples within each bin. The resulting synthetic light curve was fit with a
Gaussian profile to determine the time of minimum. This process was repeated
10,00 times with different times of minima varying from -100 sec to +100 sec
relative to the nominal zero-point time of minimum. The difference between the
true time of minimum and the fitted time of minimum was recorded.

A histogram of the time difference between the best-fit time of minimum and
the true minimum for 1,000 realizations is shown in Fig.10. The $1\sigma$
width is 4 seconds.

Figure 10: Histogram of time difference between true time of minimum and
10,000 Monte Carlo simulations.


# Denoising Distantly Supervised Named Entity Recognition via a Hypergeometric
Probabilistic Model

Wenkai Zhang1,3,*, Hongyu Lin1,*, Xianpei Han1,2,$\dagger$, Le
Sun1,2,$\dagger$, Huidan Liu1,
Zhicheng Wei4, Nicholas Jing Yuan4

###### Abstract

Denoising is the essential step for distant supervision based named entity
recognition. Previous denoising methods are mostly based on instance-level
confidence statistics, which ignore the variety of the underlying noise
distribution on different datasets and entity types. This makes them difficult
to be adapted to high noise rate settings. In this paper, we propose
_Hypergeometric Learning (HGL)_ , a denoising algorithm for distantly
supervised NER that takes both noise distribution and instance-level
confidence into consideration. Specifically, during neural network training,
we naturally model the noise samples in each batch following a hypergeometric
distribution parameterized by the noise-rate. Then each instance in the batch
is regarded as either correct or noisy one according to its label confidence
derived from previous training step, as well as the noise distribution in this
sampled batch. Experiments show that HGL can effectively denoise the weakly-
labeled data retrieved from distant supervision, and therefore results in
significant improvements on the trained models. ††$\dagger$: Corresponding
authors. ††*: Equal Contribution. H. Lin designed the method and wrote the
paper. W. Zhang conducted the experiments and detailed analysis.

## Introduction

Named Entity Recognition (NER), aiming to identify text spans pertaining to
specific semantic types such as person, organization and location, is a
foundational NLP task. In recent years, supervised neural network-based
approaches (Lample et al. 2016; Chiu and Nichols 2016; Ma and Hovy 2016) ,
which can automatically extract underlying features from annotated texts and
conduct NER recognition, have achieved promising results in almost all NER
benchmarks.

(a) Same setting, differnet types (b) Same type, different settings

Figure 1: Precision-Recall curves of different denoising algorithms on ACE2005
datasets. We can see that HGL significantly outperforms other baselines,
especially on high noise rate ACE_ALL setting.

Even with great success, supervised learning methods heavily rely on fully-
annotated training data. However, annotated training data is too expensive to
obtain, which restricts the real-world application of current NER models in
various entity types. To tackle this problem, distant supervision based NER
(DS-NER) (Ren et al. 2017; Yang et al. 2018; Shang et al. 2018; Peng et al.
2019; Nooralahzadeh, Lønning, and Øvrelid 2019) methods are introduced.
Generally, DS-NER uses easily-available resources (commonly dictionaries in
NER) to automatically label plain texts to produce large-scale data, which is
generally referred as _weakly annotated data_. And then the generated training
data will be used to learn neural network models. For example, by matching
person name “Washington” in dictionary with text “Washington is the first U.S.
president”, DS-NER will produce a training instances with entity type PERSON.
As unlabeled plain texts and entity mention dictionaries are very easy to
obtain, distant supervision based methods can significantly reduce the
annotation efforts and produce much larger scale of annotation data for model
learning.

Figure 2: The overall architecture of HGL. At each step, a batch of instances
is sent into a neural network-based denoiser to obtain its correct/noisy
confidence. Then the confidence of each instance is incorporated with the
overall hypergeometric noise distribution of this batch to obtain the training
loss of current batch.

Unfortunately, although DS-NER boosts the size of training data, it still
faces critical high noise challenge. Due to the ambiguity of natural language,
distant supervision will inevitably introduce remarkable noise into weakly
labeled data. For example, “Washington” can refer to either a person name or a
location name and DS will generate wrong PERSON instances such as “Washington,
capital of the United States”. As a result, denoising is critical to achieve
promising DS-NER performance. Currently, most DS-NER denoising methods are
based on instance-level confidence statistics, which predict whether one
instance is noise according to its feature similarity to the majority of other
instances. These instance-level denoising algorithms, unfortunately, ignore
the underlying noise structure of different datasets. For example, the
dictionary qualities and dataset properties vary significantly in different DS
settings, which will seriously affect the denoising decisions. To illustrate
this, Figure 1 shows the noise rate of DS-NER on ACE2005 datasets. We can see
that the noise rates on weakly-labeled training data vary significantly from
different entity types and different settings (e.g., named mention only v.s.
all mentions). In some settings, even a great majority of weakly-labeled
instances are noise. Consequently, conventional instance-level methods cannot
achieve robust performance on different noise DS-NER settings adaptively, and
therefore it is critical to take the noise distribution into consideration.

To this end, this paper proposes _Hypergeometric Learning (HGL)_ , which can
effectively and adaptively denoise DS-NER by further modeling and leveraging
the underlying noise structure of a dataset. HGL learns to denoise distantly
supervised labeled data by reshaping conventional training procedure according
to a noise rate-aware hypergeometric probabilistic model. Figure 2 shows the
overall framework of HGL. Specifically, to model the noise distribution during
neural network training, we formulate each mini-batch of size $B$ as a subset
being drawn from a weakly labeled data of size $N$ with noise rate $p$. Then
the noise sample size $S$ in this mini-batch will naturally follow a
hypergeometric distribution, i.e., $S\sim H(N,N\times p,B)$. Therefore, we
regard each instance as either correct or noisy one according to its label
confidence derived from previous training step, as well as the noise
distribution in the sampled batch. For example, if one training batch is
expected to have 5 instances to be correctly labeled while the others to be
noisy, then the instances with top-5 confidence will be regarded as correct
instances, while others are treated as noisy ones.

Based on the above principles, we propose a neural network-based denoiser for
DS-NER. The training loss of the denoiser is marginalized over the
hypergeometric distribution $S$ to consider all potential noise distribution
of the current training batch. For different DS-NER settings with different
noise rates, the hypergeometric distribution can be easily adapted by only
adjusting the noise rate according to the reality. Compared with traditional
instance confidence-based denoising methods, the proposed HGL explicitly
exploits the noise distribution during the denoising procedure. This enables
more robust performance among different NER settings, especially when the
noise rate is relatively high.

Generally, the main contributions of this paper are:

1) We propose to naturally model the noise distribution of DS-NER during
training using a hypergeometric probabilistic model. The model exploits the
underlying noise structures, and therefore can be more robustly adapted to
different DS-NER Settings. To the best of our knowledge, this is the first
work which tries to consider this critical information in DS-NER denoising.

2) Based on the noise model, we propose Hypergeometric Learning (HGL), a
denoising algorithm for DS-NER.  By reshaping the training procedure according
to the noise distribution and instance-level confidence, HGL can robustly
denoise the training data and therefore leads to better NER models. Because
denoising is a critical step for all distant supervision settings, we believe
our idea and method can potentially benefit many other DS tasks, such as
relation extraction.

3) We conducted experiments under various distant supervision settings on
ACE2005 and CoNLL03 NER datasets. Experimental results show that HGL can
effectively denoise DS-NER data. Furthermore, using the datasets denoised by
HGL to train a novel NER model can also significantly boost the system
performance, compared with using only weakly-labeled data and other denoising
algorithms.

## Related Work

Neural network-based models have achieved very promising results in NER field
(Lample et al. 2016; Chiu and Nichols 2016; Ma and Hovy 2016; Lin et al.
2020). These methods commonly rely on large-scale training data to learn
effective entity mention recognizer. However, fully-annotated training data is
too expensive to obtain, which limits the application of these methods to more
NER scenario.

Distant supervision (Mintz et al. 2009), aiming to efficiently generate large-
scale training data using easily available resources, has attracted much
attention (Zeng et al. 2015; Lin et al. 2016). For NER, one common practice is
to use an entity mention dictionary to automatically annotate the plain texts
to generate large-scale training data (Ren et al. 2017; Yang et al. 2018;
Shang et al. 2018; Peng et al. 2019; Nooralahzadeh, Lønning, and Øvrelid
2019). Unfortunately, DS-NER will introduce remarkable noise into the training
data, and therefore undermines the NER performance.

To tackle this problem, one solution is to only use a high precision
dictionary to avoid introducing noise, and then denoising algorithms are
applied to handle false negative instances (Yang et al. 2018; Shang et al.
2018; Peng et al. 2019; Nooralahzadeh, Lønning, and Øvrelid 2019). These
approaches have shown promising results in some specific domains such as
medical NER. However, due to the ambiguity of natural language expressions, it
is often impossible to obtain a high precision dictionary for many domains.
For example, there are many location names (e.g., Washington) which could also
be used as person names. This limits the application of these methods to more
general NER tasks.

Instead of assuming using a high precision dictionary, this paper proposes
Hypergeometric Learning, which can be directly applied to various distant
supervision settings. Furthermore, we also proposed a _Mention Blocking_
approach to extend our framework into the circumstance where extreme low false
negative rate exists. There are some previous studies along this direction,
but the majority of them are based on instance-level confidence derived from
an Expectation-Maximization style procedure (Moon 1996). These methods mainly
reduce the impact of noise on the training procedure by rescaling the
instances according to their confidence. However, they did not explicitly
model the noise distribution during the neural network training, and therefore
unable to achieve robust performance when high noise rate exists.

## Hypergeometric Learning for Denoising Distantly Supervised NER

As we mentioned above, the quality of corpus and dictionaries varies
significantly from different DS-NER settings. As a result, denoising methods
merely based on instance-level confidence to identify noise can not achieve
promising results, because in many circumstance even the majority of the
datasets would be mislabeled instances. Instead, this paper proposes to
address this issue by jointly considering the underlying noise distribution
and the instance-level confidence. The essential of our method is to learn a
neural network-based denoiser using Hypergeometric Learning (HGL) algorithm.
In the following, we will first illustrate the structure of our denoiser, then
introduce our Hypergeometric Learning algorithm.

### Denoiser

Our denoiser is a neural network model that takes the weakly annotated
instances from distant supervision as input and predicts whether these
instances are correctly labeled entity mentions. Specifically, given an input
sentence $\\{x_{1},x_{2},...,x_{n}\\}$ and a weakly labeled candidate mention
$\\{x_{i},x_{i+1},...,x_{j}\\}$ in it, we follow Baldini Soares et al. (2019)
to represent an instance $x$ by adding a start symbol [BEG] before $x_{i}$ and
an end symbol [END] after $x_{j}$:

$\mathbf{x}=[x_{0},...,[BEG],x_{i},...,x_{j},[END],...,x_{n}].$ (1)

Then $\mathbf{x}$ is sent into a Transformer (Vaswani et al. 2017) pretrained
using BERT (Devlin et al. 2019). We use $\mathbf{h_{i}}$ to denote the encoder
output of token $x_{i}$, and $\mathbf{h_{BEG}}$ and $\mathbf{h_{END}}$ to
denote the representations at token [BEG] and [END] respectively. Then we
apply an attention-based mechanism to obtain the representation of the
candidate mention $\mathbf{r}$:

$\mathbf{r}=\sum_{k=[BEG]}^{[END]}{\alpha_{k}\mathbf{h}_{k}},$ (2)

where $\alpha_{k}$ can be regarded as an important score of token $x_{k}$ to
the entire instance, which is computed by

$\alpha_{k}=\frac{\exp(\mathbf{Wh_{i}}+\mathbf{b})}{\sum_{k=[BEG]}^{[END]}\exp(\mathbf{Wh_{k}}+\mathbf{b})}.$
(3)

The representation $\mathbf{r}$ is then sent into a multi-layer perceptron
followed by a sigmoid layer to predict the probability $f(x)$ of this instance
being a correctly labeled mention:

$f(x)=\sigma(MLP(\mathbf{r})),$ (4)

where $\sigma$ is the sigmoid function and MLP is multi-layer perceptron
without active function at the top. The denoiser is expected to output a
larger $f(x)$ (i.e., more close to 1) if the instance is more likely to be a
correctly labeled instance (i.e., an correct entity mention), while output a
smaller $f(x)$ (i.e., more close to 0) if this instance is a noisy one.

### Hypergeometric Learning for Denoiser Training

Commonly, neural network-based models are trained in a supervised learning
paradigm using batch-wise learning algorithm. However, due to the lack of
correct/noisy annotation data, it is impractical to directly train our
denoiser using conventional neural network learning criteria.

#### Hypergeometric Probabilistic Model.

In this paper, instead of using annotated data, we propose to tackle this
problem by exploiting the noise distribution during denoiser training.
Specifically, if there are $N$ weakly-labeled instances in total and the
accuracy of this dataset is expected to be $p$ (i.e., the noise rate is
$1-p$). During neural network training, a common practice is to sample $B$
instances from the overall $N$ population. Then the number of correctly
labeled entity mentions $S$ in each batch will naturally follow the
hypergeometric distribution:

$S\sim H(N,K,B),$ (5)

where $K=N\times p$ is the expected size of correctly labeled instances and
$H$ is the hypergeometric distribution, which can be computed as:

$P(S=k)=\frac{\binom{K}{k}\binom{N-K}{B-k}}{\binom{N}{B}},$ (6)

where $\binom{A}{B}=\frac{A!}{B!(A-B)!}$. For clarity, we will use $Q_{k}$ to
denote $P(S=k)$ in the rest of this paper.

#### Hypergeometric Learning for DS-NER.

Once we know the distribution of correct/noise sample size in each batch, we
will evaluate the potential of each instance being a correct or noisy one
according to its label confidence derived from previous training step as well
as the noise distribution in the current batch. For example, if in this batch
there are only 5 instances expected to be the correct instances according to
the hypergeometric distribution, we will regard instances with top-5
confidence (i.e., $f(x;\theta)$ in Equation 4) as the entity mentions, while
other instances in this batch will be regarded as noisy instances. The loss
function during training will be marginalized over the hypergeometric
distribution to consider all potential noise distribution.

Formally, let {$\mathbf{x^{(1)}},\mathbf{x^{(2)}},...,\mathbf{x^{(B)}}$}
denote the instances in the current sampled batch in a descending order
according to their correct confidence in the previous step. Our Hypergeometric
Learning will optimize the denoiser in the current training set according to
the following loss function:

$L(\theta)=\sum_{i=1}^{B}\omega_{i}\log
f(x^{(i)})+(1-\omega_{i})\log(1-f(x^{(i)})),$ (7)

where $f(x)$ and $1-f(x)$ are the denoiser output of the probability of
instance $x$ being a correct or noisy instance respectively. $\omega_{i}$ is
the weight that indicates instance $x^{(i)}$ should be regarded as a correct
entity mention or a noisy instance, which is derived from the hypergeometric
distribution by:

$\omega_{i}=\sum_{k=i}^{B}Q_{k}.$ (8)

This equation indicates that if an instance ranks at $k^{th}$ according to the
confidence, it will be regarded as a correct instance only if at least $k$
instances in this batch are expected to be correct instances. Therefore, HGL
is able to exploit both the instance-level confidence and the overall noise
distribution in a unified training procedure, and thus can achieve robust
performance on various settings only by adaptively adjusting the expected
noise rate according to the underlying noise distributions.

### Dealing with Extremely Low False Negative Rate with Mention Blocking

As described above, the Hypergeometric Learning mainly focuses on the
circumstance where the dictionary is with high recall rate. However, in
practice the applied dictionary commonly cannot cover all mentions in the
dataset, and therefore requires the model to also deal with false negative
instances. Considering the fact that entity mentions of a specific type only
cover a very minority of all noun phrases (commonly less than 1%), it is quite
difficult to distinguish false negative instances from the true negative ones
using standard denoising techniques. To tackle this problem, we further
propose a _Mention Blocking_ approach, which first screens possible entity
mentions of a specific type by leveraging the mention compositional knowledge
learned from the dictionary, then the potential mentions will be further
introduced into Hypergeometric Learning procedure, so that false negatives can
also be denoised.

Specifically, we regard all noun phrases in the entire corpus (except mentions
matched by the distant supervision dictionary) as the potential mentions. Then
inspired by Lin et al. (2019b), we learn a Bert-based classifier from given
dictionaries, which takes a context-free noun phrase (mention utterance) as
the input, and outputs a score indicating how possible the input phrase could
be a valid mention of such entity type. For example, assigning “Washington” a
high PER score and “New York” a low PER score. After that, we collect the noun
phrases with high classification scores as the possible entity mention block.
And only the instances within the block will be regarded as potential false
negative mentions. Then Hypergeometric Learning will be applied to denoise
these potential mentions, which tries to further distinguish false negative
instances (i.e., the noises in true negative instances) from all noun phrases
screened out. The learning procedure regards the false negative instances as
the noise out of all blocked noun phrases. And the denoiser is learned jointly
from denoising both positive and negative instances, which means that we will
learn one denoiser with the loss from both these blocked noun phrases and the
instances from distant supervision.

## Experiments with Synthetic Dictionaries

In this section, we conducted experiments to verify the efficacy of the
proposed Hypergeometric Learning algorithm. Specifically, we evaluate our
method to answer the following three questions:

  * •

Can Hypergeometric Learning-based denoiser effectively identify the noisy
instances in the weakly-labeled training data?

  * •

Can the proposed Mention Blocking approach effectively detect false negative
instances and improves the denoising performance?

  * •

Can using the dataset denoised by HGL significantly improve the performance of
conventional NER models such as Bert-CRF tagger?

### Experimental Settings

#### Dataset.

We conduct experiments on ACE2005 (Walker et al. 2006) benchmark, which
contains 7 entity categories. We use the same experimental setup as (Ju, Miwa,
and Ananiadou 2018; Katiyar and Cardie 2018; Wang and Lu 2018; Straková,
Straka, and Hajič 2019; Lin et al. 2019a), where the entire dataset is split
into 8:1:1 for training, developing and testing respectively. We keep only the
outmost entity mentions and therefore ignore the overlapping mentions of the
same type. To verify the adaptiveness of the proposed methods, we conducted
experiments based on two settings: 1) Vanilla ACE2005, which considers all
named, nominal and pronominal mentions; 2) ACE2005 NAM111For ACE2005 NAM, we
only conduct experiments on PER, ORG and GPE because other entity types do not
have sufficient instances in both training and test set., which considers only
named mentions in the original dataset, while ignores nominal and pronominal
mentions. We refer to these two settings as ACE_ALL and ACE_NAM respectively.
Generally speaking, ACE_ALL is with relatively large noise ratio compared with
ACE_NAM due to the stronger ambiguity of nominal and pronominal utterances.

#### Distant Supervision.

To effectively verify the adaptiveness of Hypergeometric Learning under
different dictionary quality settings, in this section we conducted
experiments using synthetic dictionaries extracted from the training data. And
in the next section we will conduct an experiment on the real-world
dictionaries to further verify the effectiveness of HGL.

Specifically, in this section we used all mentions extracted from the original
training data as our distant supervision source dictionary. We then used the
dictionary to label all documents in the original training set to obtain
weakly-labeled training instance for our denoiser. This guarantees that the
dictionary used in our experiment is with high recall rate. The macro-averaged
noise rate on ACE_ALL and ACE_NAM are 0.79 and 0.20, respectively. This
indicates the necessity of introducing denoising algorithm to distant
supervision NER, even when a golden dictionary is applied.

#### Noise Rate Estimation.

For HGL, we estimate the noise ratio of each entity type from the development
set. Because we may not estimate noises rate exactly accurately in reality, we
approximate noise rates with its nearest 5% percentile. For example, if the
noise rate is 34.1% in the development set, we will use 35% as the noise rate
during HGL learning.

#### Implementation Detail.

We use Allennlp (Gardner et al. 2018), an open-source NLP research library to
implement our method and conduct experiments. For HGL, we set batch size to
150 and use Adam (Kingma and Ba 2015) as optimizer, the learning rate is
$1\times 10^{-5}$. We used _cased_L-24_H-1024_A-16_ as pretrained encoder. We
openly release our source code at github.com/zwkatgithub/HGL.

### Baselines

We compared HGL with the following baselines:

  * •

Supervised Learning (Supervised), which used fully-annotated data to directly
train the NER model, which can serve as the upper bound of all denoising
methods.

  * •

Naive Distant Supervision (Naive), which directly used the dictionary to label
the plain texts and does not conduct any denoising on the annotation results.

  * •

Instance-level Expectation Maximization (Instance) (Moon 1996), which
identifies noisy training instances in a self-supervision manner according to
the instance-level confidence. Specifically, the loss of each instance is
reshaped by:

$\displaystyle L(x)=$ $\displaystyle f(x;\theta_{t})\log f(x)$ (9)
$\displaystyle+[1-f(x;\theta_{t})]\log[1-f(x)].$

This can be regarded as a baseline that only considers instance-level
confidence but ignores the overall noise distribution during the training
procedure.

  * •

XR-Loss (Ben Noach and Goldberg 2019) (XR), which considers the proportional
mapping between distant supervision label mentions to correctly labeled and
noisy instances. In other words, this method takes only the overall noise rate
into consideration but ignores the instance-level confidence.

  * •

Positive-Unlabeled Learning (Peng et al. 2019) (PU-Learning), which tackles
the mislabeling in distant supervision by regarding it as a positive-unlabeled
learning problem and reshapes the training procedure with a regularizer
considers the mislabeling ratio. Note that this method is to learn a NER model
directly on DS-NER datasets, which is unable to identify the exact noise
instance during learning.

Dataset | ACE_ALL | ACE_NAM
---|---|---
Method | PER | ORG | GPE | LOC | FAC | VEH | WEA | PER | ORG | GPE
Instance | 63.97 | 72.88 | 66.15 | 83.75 | 77.61 | 79.23 | 80.92 | 71.65 | 50.07 | 52.41
XR | 56.51 | 77.42 | 67.35 | 82.03 | 78.62 | 84.63 | 74.47 | 66.97 | 55.30 | 59.46
HGL | 80.35 | 82.28 | 86.50 | 93.86 | 90.56 | 90.03 | 94.00 | 71.68 | 59.22 | 64.35
Table 1: The AUC scores of denoising performance of different algorithms on
ACE_ALL and ACE_NAM. We can see that HGL significantly outperfoms other
methods by a large margin.

### Denoising Performance

The first group of our experiments was conducted to evaluate how well the
proposed HGL and baseline methods can denoise the weakly labeled data. We
evaluate the denoising performance using ranking based metrics to see how well
the algorithms can identify correct instances out of noisy instances.

(a) ACE_ALL

(b) ACE_NAM

Figure 3: Precision-Recall curves of different denoising algorithms on ACE2005
datasets. We can see that HGL significantly outperforms other baselines,
especially on high noise rate ACE_ALL setting.

Figure 3 shows the macro-averaged precision-recall curves of HGL and other
baselines on ACE_ALL and ACE_NAM respectively, while Table 1 shows the overall
AUC scores of all systems. We can see that the proposed HGL algorithm
significantly outperforms other denoising algorithms. For AUC scores, HGL
outperforms other baselines by a large margin on almost all settings, which
demonstrates the effectiveness of our method. Furthermore, we can see that
from the Precision-Recall curves, the performance of other baselines
dramatically drops as the recall increases, which indicates that HGL performs
significantly better, especially on correct instances that are hard to be
distinguished from noise. Besides, from both the AUC scores and the Precision-
Recall curves, the improvement of HGL, compared with other baselines, is more
remarkable on ACE_ALL than on ACE_NAM. This is because, as we mentioned above,
the overall noise rate of ACE_ALL is much larger than ACE_NAM, which makes
previous methods that do not consider both instance confidence and noise
distribution can not achieve promising performance in this circumstance.
Generally speaking, owing to the effective exploit both instance-level
confidence and the overall noise distribution, HGL can achieve much better and
robust denoising performance for distant supervision NER.

| ACE_ALL | ACE_NAM
---|---|---
Model | PER | ORG | GPE | LOC | FAC | VEH | WEA | PER | ORG | GPE
Supervised | 85.51 | 66.45 | 76.92 | 65.96 | 58.97 | 52.29 | 73.12 | 84.80 | 76.55 | 89.89
Naive | 60.78 | 32.70 | 33.28 | 7.91 | 15.22 | 8.12 | 7.47 | 81.03 | 70.36 | 75.63
Instance | 54.52 | 32.53 | 34.43 | 30.99 | 28.07 | 11.32 | 18.75 | 82.84 | 51.83 | 76.53
XR | 50.60 | 38.10 | 28.38 | 22.54 | 27.43 | 35.21 | 14.46 | 80.12 | 57.92 | 71.40
PU-Learning | 56.54 | 23.91 | 24.63 | 5.32 | 11.96 | 7.13 | 6.81 | 76.36 | 53.86 | 75.79
HGL | 69.65 | 47.09 | 59.94 | 60.87 | 51.56 | 42.18 | 59.77 | 84.45 | 67.51 | 77.56
Table 2: The F1 scores of Bert-CRF models with different denoising algorithms
on ACE2005 test sets. We can see that the proposed HGL significantly
outperforms other baselines, especially on those high noise-ratio settings.

### Learning NER models with Denoised Data

In this section, we investigate whether denoised data can effectively improve
the training of a commonly-used NER model. Specifically, for each denoising
algorithm, we regarded $N\times p$ weakly-labeled instances in the training
set with the highest confidence derived from the denoiser as the golden entity
mentions. Then the denoised data is used to train widely-used BERT-CRF NER
taggers to evaluate the performance. The performance is evaluated using the
final F1-score achieved in the original test set.

Table 2 shows the results. We can see that:

1) Denoising is critical for distant supervision NER. Compared with the naive
distant supervision setting, using denoised data to train Bert-CRF tagger
achieves significant improvements.

2) HGL outperforms all other denoising baselines on all settings by a large
margin, especially on those settings with larger noise rates. This result is
not surprising because, as we have shown before, the denoising performance of
HGL is significantly better than other baselines. This again demonstrates the
effectiveness of HGL in modeling noise distribution and instance-level
confidence.

3) There is still a gap between the supervised learning and distant
supervision. This is because even remarkable progress has been made by HGL,
there still exists some errors in denoising. This indicates the potential
improvement in denoising distant supervision NER in the future.

Method | PER | ORG | GPE
---|---|---|---
Instance | 74.79 | 71.70 | 69.85
XR | 79.07 | 72.93 | 73.13
HGL (False Positive Only) | 78.88 | 72.69 | 76.06
HGL + Mention Blocking | 79.80 | 73.97 | 76.13
Table 3: The AUC scores of denoising performance on ACE2005_NAM when false
negative instances exist. Figure 4: The coverage of false negative instances
w.r.t. noun phrase screening rate. We can see that Mention Blocking can cover
more than 50% false negative instances by considering less than 10% noun
phrases.

### Efficacy of Mention Blocking

This section conducted experiments to verify the effectiveness of the proposed
Mention Blocking approach and how well HGL can work when false negative
exists. For this, we randomly dropped 50% mentions in the dictionary derived
from the ACE2005 training data. Therefore, the weakly-labelled data will
contain a number of both false positive and false negative instances. We
compared HGL with or without Mention Blocking to verify the effectiveness of
processing false negative instances. We keep noun phrases that with top 10%
confidence in Mention Blocking as the potential false negative mentions, and
introduce them into the HGL procedure.

Table 3 shows the AUC scores of denoising performance under this setting
222The AUC scores does not consider the noun phrases that are not identified
by either distant supervision or Mention Blocking, because all methods will
treat them as negative instances and they will not influence the evaluation
results.. We can see that HGL with Mention Blocking significantly outperforms
the one that only focus on false positive instances, as well as other
baselines. This shows the effectiveness of the proposed Mention Blocking
approach. To further demonstrate this, Figure 4 shows the coverage of the
false negative instances with respect to the total noun phrases screened out.
From this figure, we can see that Mention Blocking is very effective to find
the potential entity mentions: top 10% noun phrases screened out can cover
more than 50% false negative mentions. This verify the efficacy of Mention
Blocking method.

## Experiments with Real-world Dictionaries

To further verify the effectiveness of the proposed hypergeometric learning
method on real-world DS-NER settings. We follow Peng et al. (2019) to further
conduct an experiments on CoNLL2003 dataset, using a real-world dictionary as
the source of distant supervision.

#### Distant Supervision.

We followed Peng et al. (2019) and used their extracted dictionary to conduct
the evaluation. The dictionary was extracted from government websites
333http://www.ons.gov.uk/ons/index.html and Wikipedia pages
444https://en.wikipedia.org/wiki/List_of_countries_by_natio
nal_capital_largest_and_secondlargest_cities. We then used the dictionaries to
obtain weakly-labeled training instances using forward maximum matching.
According to Peng et al. (2019), this dictionary is with very high precision
but relatively low coverage rate on CoNLL03, and the denoising algorithm
should be conducted to distinguish false negative instances. Therefore, we
directly apply HGL on the negative instances.

#### Overall Performance.

| PER | ORG | LOC | MISC | Overall
---|---|---|---|---|---
XR | 91.97 | 64.46 | 71.00 | 55.11 | 73.65
PU-Learning | 93.93 | 67.34 | 65.58 | 50.94 | 72.34
HGL | 92.78 | 68.30 | 71.62 | 56.69 | 74.87

Table 4: DS-NER performance on CoNLL03 test set evaluated using span-based micro-F1 scores. We can see that HGL significantly outperform other denoising algorithms. Recall | Model | PER | ORG | LOC | MISC
---|---|---|---|---|---
25% | XR | 99.00 | 71.41 | 85.38 | 65.32
PU-Learning | 95.73 | 87.22 | 86.45 | 90.03
HGL | 96.73 | 85.12 | 96.33 | 90.82
50% | XR | 98.32 | 60.20 | 82.25 | 62.26
PU-Learning | 92.76 | 77.03 | 75.90 | 80.72
HGL | 97.96 | 78.22 | 90.29 | 66.68
75% | XR | 96.70 | 48.88 | 69.13 | 45.32
PU-Learning | 86.89 | 55.88 | 41.78 | 20.36
HGL | 97.93 | 70.40 | 70.98 | 43.04
Table 5: The token-level precisions under different recall rates on the
denoised distantly-supervised CoNLL03 training set. We can see that HGL
achieved the best performance on the majority of settings.

Table 4 shows the overall experimental results measured using standard span-
based micro-F1555From the source code provided by Peng et al. (2019) on
Github, we find that they used a token-based evaluation criteria, which is
different from the most widely-used span-based criteria and therefore over-
estimated NER performance. In this paper, we reproduced their results using
the provided source code and compared it with HGL/XR using the standard span-
based micrio-F1 metric.. We can see that HGL significantly outperforms other
baselines. Furthermore, Table 5 shows the change of precision w.r.t. different
target recall rates. We can see that under the majority of experiment groups,
HGL can achieve better precision than other previous methods. This demonstrate
that HGL can achieve steadily better performance under ranking measurements,
i.e., instances obtained high confidence from HGL are more likely to be true
entity mentions. All these experimental results consist with results in the
previous section, which further demonstrate the effectiveness of HGL.

## Conclusions

This paper proposes a new denoising algorithm – Hypergeometric Learning, which
can model and exploit both the instance-level evidence and the underlying
noise distribution for denoising DS-NER. Experimental results on both high-
precision and high-recall DS-NER settings show that, by further taking the
corpus-level noise distribution into consideration, HGL can effectively
denoise the DS-labeled data, learn an effective NER model, and achieve robust
performance on different entity types and different dataset settings.

## Acknowledgements

This research work is supported by National Key Research and Development
Program of China under Grant No.2019YFC1521200, the National Natural Science
Foundation of China under Grants no. U1936207 and 61772505, Beijing Academy of
Artificial Intelligence (BAAI2019QN0502), scientific research projects of the
State Language Commission (YW135-78), and in part by the Youth Innovation
Promotion Association CAS(2018141). Moreover, we thank all reviewers for their
valuable comments and suggestions.

## References

  * Baldini Soares et al. (2019) Baldini Soares, L.; FitzGerald, N.; Ling, J.; and Kwiatkowski, T. 2019.  Matching the Blanks: Distributional Similarity for Relation Learning.  In _ACL_ , 2895–2905.
  * Ben Noach and Goldberg (2019) Ben Noach, M.; and Goldberg, Y. 2019.  Transfer Learning Between Related Tasks Using Expected Label Proportions.  In _EMNLP_ , 31–42.
  * Chiu and Nichols (2016) Chiu, J. P.; and Nichols, E. 2016.  Named Entity Recognition with Bidirectional LSTM-CNNs.  _TACL_ 4: 357–370.
  * Devlin et al. (2019) Devlin, J.; Chang, M.-W.; Lee, K.; and Toutanova, K. 2019.  BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding.  In _NAACL_ , 4171–4186.
  * Gardner et al. (2018) Gardner, M.; Grus, J.; Neumann, M.; Tafjord, O.; Dasigi, P.; Liu, N. F.; Peters, M.; Schmitz, M.; and Zettlemoyer, L. 2018.  AllenNLP: A Deep Semantic Natural Language Processing Platform.  In _Proceedings of Workshop for NLP Open Source Software (NLP-OSS)_ , 1–6.
  * Ju, Miwa, and Ananiadou (2018) Ju, M.; Miwa, M.; and Ananiadou, S. 2018.  A neural layered model for nested named entity recognition.  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long Papers)_ , 1446–1459.
  * Katiyar and Cardie (2018) Katiyar, A.; and Cardie, C. 2018.  Nested named entity recognition revisited.  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long Papers)_ , 861–871.
  * Kingma and Ba (2015) Kingma, D. P.; and Ba, J. 2015.  Adam: A Method for Stochastic Optimization.  In _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_.
  * Lample et al. (2016) Lample, G.; Ballesteros, M.; Subramanian, S.; Kawakami, K.; and Dyer, C. 2016.  Neural Architectures for Named Entity Recognition.  In _ACL_ , 260–270.
  * Lin et al. (2019a) Lin, H.; Lu, Y.; Han, X.; and Sun, L. 2019a.  Sequence-to-Nuggets: Nested Entity Mention Detection via Anchor-Region Networks.  In _ACL_ , 5182–5192.
  * Lin et al. (2019b) Lin, H.; Lu, Y.; Han, X.; Sun, L.; Dong, B.; and Jiang, S. 2019b.  Gazetteer-Enhanced Attentive Neural Networks for Named Entity Recognition.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , 6232–6237. Hong Kong, China.  doi:10.18653/v1/D19-1646.
  * Lin et al. (2020) Lin, H.; Lu, Y.; Tang, J.; Han, X.; Sun, L.; Wei, Z.; and Yuan, N. J. 2020.  A Rigorous Study on Named Entity Recognition: Can Fine-tuning Pretrained Model Lead to the Promised Land?  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , 7291–7300. Online.  doi:10.18653/v1/2020.emnlp-main.592.
  * Lin et al. (2016) Lin, Y.; Shen, S.; Liu, Z.; Luan, H.; and Sun, M. 2016.  Neural Relation Extraction with Selective Attention over Instances.  In _ACL_ , 2124–2133.
  * Ma and Hovy (2016) Ma, X.; and Hovy, E. 2016.  End-to-end Sequence Labeling via Bi-directional LSTM-CNNs-CRF.  In _ACL_ , 1064–1074.
  * Mintz et al. (2009) Mintz, M.; Bills, S.; Snow, R.; and Jurafsky, D. 2009.  Distant supervision for relation extraction without labeled data.  In _ACL_ , 1003–1011.
  * Moon (1996) Moon, T. K. 1996.  The expectation-maximization algorithm.  _IEEE Signal processing magazine_ 13(6): 47–60.
  * Nooralahzadeh, Lønning, and Øvrelid (2019) Nooralahzadeh, F.; Lønning, J. T.; and Øvrelid, L. 2019.  Reinforcement-based denoising of distantly supervised NER with partial annotation.  In _Proceedings of the 2nd Workshop on Deep Learning Approaches for Low-Resource NLP (DeepLo 2019)_ , 225–233.
  * Peng et al. (2019) Peng, M.; Xing, X.; Zhang, Q.; Fu, J.; and Huang, X. 2019.  Distantly Supervised Named Entity Recognition using Positive-Unlabeled Learning.  In _ACL_ , 2409–2419.
  * Ren et al. (2017) Ren, X.; Wu, Z.; He, W.; Qu, M.; Voss, C. R.; Ji, H.; Abdelzaher, T. F.; and Han, J. 2017.  CoType: Joint Extraction of Typed Entities and Relations with Knowledge Bases.  In _WWW_ , 1015–1024.
  * Shang et al. (2018) Shang, J.; Liu, L.; Gu, X.; Ren, X.; Ren, T.; and Han, J. 2018.  Learning Named Entity Tagger using Domain-Specific Dictionary.  In _EMNLP_ , 2054–2064.
  * Straková, Straka, and Hajič (2019) Straková, J.; Straka, M.; and Hajič, J. 2019.  Neural architectures for nested NER through linearization.  _arXiv preprint arXiv:1908.06926_ .
  * Vaswani et al. (2017) Vaswani, A.; Shazeer, N.; Parmar, N.; Uszkoreit, J.; Jones, L.; Gomez, A. N.; Kaiser, L. u.; and Polosukhin, I. 2017.  Attention is All you Need.  In _Advances in Neural Information Processing Systems 30_ , 5998–6008.
  * Walker et al. (2006) Walker, C.; Strassel, S.; Medero, J.; and Maeda, K. 2006.  ACE 2005 multilingual training corpus.  _Linguistic Data Consortium, Philadelphia_ 57\.
  * Wang and Lu (2018) Wang, B.; and Lu, W. 2018.  Neural segmental hypergraphs for overlapping mention recognition.  _arXiv preprint arXiv:1810.01817_ .
  * Yang et al. (2018) Yang, Y.; Chen, W.; Li, Z.; He, Z.; and Zhang, M. 2018.  Distantly Supervised NER with Partial Annotation Learning and Reinforcement Learning.  In _COLING_ , 2159–2169.
  * Zeng et al. (2015) Zeng, D.; Liu, K.; Chen, Y.; and Zhao, J. 2015.  Distant Supervision for Relation Extraction via Piecewise Convolutional Neural Networks.  In _EMNLP_ , 1753–1762. Lisbon, Portugal.