Qualitative properties for elliptic problems

with CKN operators

Huyuan<EMAIL_ADDRESS>Yishan<EMAIL_ADDRESS>

Department of Mathematics, Jiangxi Normal University,
Nanchang, Jiangxi 330022, PR China

###### Abstract

The purpose of this paper is to study basic property of the operator

${\mathcal{L}}_{\mu_{1},\mu_{2}}u=-\Delta+\frac{\mu_{1}}{|x|^{2}}x\cdot\nabla+\frac{\mu_{2}}{|x|^{2}},$

which generates at the origin due to the critical gradient and the Hardy term,
where $\mu_{1},\mu_{2}$ are free parameters. This operator arises from the
critical Caffarelli-Kohn-Nirenberg inequality. We analyze the fundamental
solutions in a weighted distributional identity and obtain the Liouville
theorem for the Lane-Emden equation with that operator, by using the
classification of isolated singular solutions of the related Poisson problem
in a bounded domain $\Omega\subset\mathbb{R}^{N}$ ($N\geq 2$) containing the
origin.

Keywords:  Degenerate operator, fundamental solution, Liouville theorem, Lane-
Emden problem.

MSC2010: 35J75; 35B44; 35B53.

##  1 Introduction

The well-known Caffarelli-Kohn-Nirenberg inequality (CKN inequality for short)
proposed in [7] states as following

$\Big{(}\int_{\mathbb{R}^{N}}|x|^{-b(p+1)}|u|^{p+1}dx\Big{)}^{\frac{2}{p+1}}\leq
C_{a,b,N}\int_{\mathbb{R}^{N}}|x|^{-2a}|\nabla u|^{2}dx,$

where $N\geq 2$,

$-\infty<a<\frac{N-2}{2},\quad a\leq b\leq a+1\quad{\rm and}\quad
p=\frac{N+2(1+a-b)}{N-2(1+a-b)}.$

A critical CKN inequality with $b=a+1$ (see [8, 1, 29]) is

$\int_{\mathbb{R}^{N}}|x|^{-2a}|\nabla
u|^{2}dx\geq\Big{(}\frac{N-2-2a}{2}\Big{)}^{2}\int_{\mathbb{R}^{N}}|x|^{-2(a+1)}|u|^{2}dx,$
(1.1)

which, for $a=0$ and $N\geq 3$, reduces to the classical Hardy inequality

$\int_{\mathbb{R}^{N}}|\nabla
u|^{2}dx\geq\frac{(N-2)^{2}}{4}\int_{\mathbb{R}^{N}}\frac{|u|^{2}}{|x|^{2}}dx.$
(1.2)

It is known that the related elliptic operator arising form (1.1)

$\displaystyle{\mathcal{L}}_{a}u$ $\displaystyle=$ $\displaystyle-{\rm
div}(|x|^{-2a}\nabla
u)-\frac{\big{(}N-2-2a\big{)}^{2}}{4}\frac{u}{|x|^{2(a+1)}}$ $\displaystyle=$
$\displaystyle|x|^{-2a}\Big{(}-\Delta u+\frac{2a}{|x|^{2}}x\cdot\nabla
u-\frac{\big{(}N-2-2a\big{)}^{2}}{4}\frac{u}{|x|^{2}}\Big{)},\quad\forall\,u\in
C^{2}_{c}(\mathbb{R}^{N}),$

where $a<\frac{N-2}{2}$. If we take $\mu_{1},\mu_{2}$ to replace $2a$ and
$-\frac{1}{4}\big{(}N-2-2a\big{)}^{2}$ respectively in ${\mathcal{L}}_{a}$, we
propose a degenerate operator

${{\mathcal{L}}}_{\mu_{1},\mu_{2}}=-\Delta+\frac{\mu_{1}}{|x|^{2}}x\cdot\nabla+\frac{\mu_{2}}{|x|^{2}}$
(1.3)

with two free parameters $\mu_{1},\mu_{2}$. Note that the operator
${{\mathcal{L}}}_{\mu_{1},\mu_{2}}$ degenerates at the origin both for the
gradient term and the critical Hardy term, and here we call it the CKN
operator. Our aim is to consider the qualitative properties of the solutions
of the elliptic equations with that operator.

When $\mu_{2}=0$,
${{\mathcal{L}}}_{\mu_{1},0}=-\Delta+\frac{\mu_{1}}{|x|^{2}}x\cdot\nabla$ is a
type of degenerate elliptic operator, which together with its divergence form,
plays an important role in the harmonic analysis, see [27], it attracts lots
of attentions, the basic regularities of related equations [18], qualitative
properties for equation with more general degenerate operators in divergence
form [28].

When $\mu_{1}=0$, ${{\mathcal{L}}}_{0,\mu_{2}}$ reduces to the Hardy-Leray
operator, here we can write
${{\mathcal{L}}}_{\mu}=-\Delta+\frac{\mu}{|x|^{2}}$, which is the prototype of
the degenerate operators. The equations with Hardy-Leray operators has been
studied extensively in the last decades. The authors in [24] initiated the
analysis of isolated singular solutions of semilinear problems
$u\mapsto{\mathcal{L}}_{\mu}u+g(u)$ under the condition
$\mu\geq-\frac{(N-2)^{2}}{4}$, where $\frac{(N-2)^{2}}{4}$ is the best
constant of Hardy inequality (1.2), more related Hardy inequalities refer to
[5, 20, 13]. Normally, the distributional solution of the Hardy problem
$\mathcal{L}_{\mu}u+g=0\ {\rm in}\ \Omega$ would be proposed as

$\int_{\Omega}u\mathcal{L}_{\mu}\xi\,dx+\int_{\Omega}g\xi\,dx=0,\quad\forall\,\xi\in
C^{\infty}_{c}(\Omega),$

where $g$ is a nonlinearity of $x$ and $u$. In this distributional sense, [4,
17, 19] show the existence for particular nonlinearity of $u$ under some
restriction that $N\geq 3$ and $\mu\in[\mu_{0},0)$. Later on, Cîrstea at el in
[14], Cîrstea in [15] classified the isolated singular classical solution of
$\mathcal{L}_{\mu}u+b(x)h(u)=0$ in $\Omega\setminus\\{0\\}$, where both $b$
and $h$ consist of regularly varying and slowly varying parts (see their
definitions in §1.2.2 of [15]). There a solution is considered as a
$C^{1}(\Omega\setminus\\{0\\})$-solution in the sense of distributions in
$\Omega\setminus\\{0\\}$, that is,

$\int_{\Omega}\nabla u\nabla\varphi
dx-\int_{\Omega}\frac{\lambda}{|x|^{2}}u\varphi
dx+\int_{\Omega}b(x)h(u)\varphi dx=0,\quad\forall\,\varphi\in
C^{1}_{c}(\Omega\setminus\\{0\\}).$

Thanks to a notion of weak solutions of ${\mathcal{L}}_{\mu}u=0$ combined with
a dual formulation of the equation introduced in [9] the equation

${\mathcal{L}}_{\mu}u+g(u)=\nu\quad{\rm in}\,\,\Omega\qquad u=0\quad{\rm
on}\,\partial\Omega,$

where $\Omega$ is a bounded smooth domain, $g$ is a continuous nondecreasing
function and $\nu$ is a Radon measure. When the pole of the Leray-Hardy
potential is addressed on the boundary of the domain $\Omega$, [11, 12] extend
the approach to classify the boundary isolated singular solutions of Poisson
problem

${\mathcal{L}}_{\mu}u+g(u)=\nu_{1}\quad{\rm in}\ \,\Omega,\qquad
u=\nu_{2}\quad{\rm on}\ \partial\Omega,$

where $\nu_{1},\,\nu_{2}$ are bounded Radon measures respectively on $\Omega$
and $\partial\Omega$. Recently, the singularities of the Hardy problems have
been also studied extensively [2, 3, 26, 21, 22, 23] and the references
therein.

When $\mu_{1},\mu_{2}\not=0$, the CKN operator
${{\mathcal{L}}}_{\mu_{1},\mu_{2}}$ defined in (1.3) degenerates at the origin
thanks to both the gradient term and the Hardy potential. It is worth noting
that the gradient term has the singularity $|x|^{-1}$ indeed, which is also
critical at the origin, compared with the Laplacian operator. The CNK operator
in the diverging form ${\mathcal{L}}_{a}$ has been discussed in [18, 25]. So
our aim is to characterize the roles of the two critical terms in the related
fundamental solution and the related Poisson problems.

To consider the fundamental solution of ${\mathcal{L}}_{\mu_{1},\mu_{2}}$, we
provide the following setting of $\mu_{1},\mu_{2}$ in this article is the
following:

$\mu_{1}\in\mathbb{R},\quad\ \ \mu_{2}\geq-\frac{(2-N+\mu_{1})^{2}}{4}.$

In this setting, direct computation shows that the homogeneous problem

${\mathcal{L}}_{\mu_{1},\mu_{2}}u=0\quad{\rm in}\
\,\mathbb{R}^{N}\setminus\\{0\\}$ (1.4)

has two the radially symmetric solutions

$\Phi_{\mu_{1},\mu_{2}}(x)=\left\\{\begin{array}[]{lll}|x|^{\tau_{-}(\mu_{1},\mu_{2})}&{\rm
if}\ \,\mu_{2}>-\frac{(2-N+\mu_{1})^{2}}{4},\\\\[5.69054pt]
-|x|^{\tau_{-}(\mu_{1},\mu_{2})}{\ln|x|}&{\rm if}\
\,\mu_{2}=-\frac{(2-N+\mu_{1})^{2}}{4}\end{array}\right.$ (1.5)

and

$\Gamma_{\mu_{1},\mu_{2}}(x)=|x|^{\tau_{+}(\mu_{1},\mu_{2})},$ (1.6)

where

$\tau_{-}(\mu_{1},\mu_{2})=\frac{(2-N+\mu_{1})-\sqrt{(2-N+\mu_{1})^{2}+4\mu_{2}}}{2}$

and

$\tau_{+}(\mu_{1},\mu_{2})=\frac{(2-N+\mu_{1})+\sqrt{(2-N+\mu_{1})^{2}+4\mu_{2}}}{2}.$

When $\mu_{2}=-\frac{(2-N+\mu_{1})^{2}}{4}$, we observe that

$\tau_{-}(\mu_{1},\mu_{2})=\tau_{+}(\mu_{1},\mu_{2})=\frac{2-N+\mu_{1}}{2}:=\tau_{0}(\mu_{1})$

and

$\tau_{0}(\mu_{1})<0\ \iff\ \mu_{1}<N-2.$

In the following, we use the notation of $\tau_{\pm}$ replaced by
$\tau_{\pm}(\mu_{1},\mu_{2})$, $\tau_{0}$ done by $\tau_{0}(\mu_{1})$ for
simplicity if there is no confusion.

###### Theorem 1.1

Assume that $N\geq 2$, $\mu_{1}\in\mathbb{R}$ and

$\mu_{2}\geq-\frac{(2-N+\mu_{1})^{2}}{4}.$

Let $d\gamma_{\mu_{1},\mu_{2}}:=|x|^{\tau_{+}-\mu_{1}}dx$ and

${{\mathcal{L}}}_{\mu_{1},\mu_{2}}^{*}=-\Delta+(-2\tau_{+}+\mu_{1})\frac{x}{|x|^{2}}\cdot\nabla,$

then we have

$\int_{{R}^{N}}\Phi_{\mu_{1},\mu_{2}}\,{\mathcal{L}}^{*}_{\mu_{1},\mu_{2}}(\xi)d\gamma_{\mu_{1},\mu_{2}}=c_{\mu_{1},\mu_{2}}\xi(0),\
\ \ \forall\,\xi\in C^{2}_{c}(\mathbb{R}^{N}),$ (1.7)

where

$c_{\mu_{1},\mu_{2}}=\left\\{\begin{array}[]{lll}\sqrt{(2-N+\mu_{1})^{2}+4\mu_{2}}\,|\mathbb{S}^{N-1}|&{\rm
if}\ \,\mu_{2}>-\frac{(2-N+\mu_{1})^{2}}{4},\\\\[5.69054pt]
|\mathbb{S}^{N-1}|&{\rm if}\
\,\mu_{2}=-\frac{(2-N+\mu_{1})^{2}}{4}.\end{array}\right.$

###### Remark 1.1

$(i)$ The identity (1.7) means that

${\mathcal{L}}_{\mu_{1},\mu_{2}}\Phi_{\mu_{1},\mu_{2}}=c_{\mu_{1},\mu_{2}}\delta_{0}$

in the weighted distributional sense (1.7), where $\delta_{0}$ is Dirac mass
concentrates at the origin.

$(ii)$ When $\mu_{1}=0$, the fundamental solution expressed by Dirac mass in
the weighted distributional is derived in [9] and singular solution of Lane-
Emden equations with Hardy operators is classified in [10] in this framework.

$(iii)$ When $\mu_{2}=0$, $\mu_{1}\leq N-2$, $\Gamma_{\mu_{1},\mu_{2}}=1$ and
$\Phi_{\mu_{1},\mu_{2}}(x)=-|x|^{\tau_{-}}\ln|x|$, which is the fundamental
solution of ${\mathcal{L}}_{\mu_{1},0}$ in the weighted distributional
identity (1.7), which coincides the normal distributional identity and
${\mathcal{L}}^{*}_{\mu_{1},0}={\mathcal{L}}_{\mu_{1},0}$.

Our second concern is to consider the nonexistence of positive solution for
Lane-Emden equation with CKN operators

$\left\\{\begin{array}[]{lll}{\mathcal{L}}_{\mu_{1},\mu_{2}}u\geq Qu^{p}&{\rm
in}\ \ \Omega\setminus\\{0\\},\\\\[5.69054pt]
\phantom{{\mathcal{L}}_{\mu_{1},\mu_{2}}}u\geq 0&{\rm in}\ \
\,\partial\Omega,\end{array}\right.$ (1.8)

where $p>0$, $\Omega$ is a bounded domain containing the origin and the
potential $Q\in C^{\beta}_{loc}(\mathbb{R}^{N}\setminus\\{0\\})$ with
$\beta\in(0,1)$ is a positive function such that for some $\theta>-2$,

$\liminf_{|x|\to 0^{+}}Q(x)|x|^{-\theta}>0.$ (1.9)

Before stating our main results, let us involve two critical exponents

$p^{\\#}_{\mu_{1},\mu_{2},\theta}=1+\frac{2+\theta}{-\tau_{+}}\quad{\rm
for}\;\ -\frac{(N-2-\mu_{1})^{2}}{4}\leq\mu_{2}<0.$ (1.10)

Here $p^{\\#}_{\mu_{1},\mu_{2},\theta}$ is a particular exponent appearing
only for $\tau_{+}<0$, which is the essence for the nonexistence of positive
solutions to (1.8). Note that any positive solution of (1.8) blows up at least
like $\Gamma_{\mu}$ at the origin and this singularity would be improved by
the interact of the nonlinearity $Qu^{p}$, which may lead to an unadmissible
singularity in some weighted $L^{1}$ space. Inspired by this observation, [6,
17] showed the nonexistence of positive solution of (1.8) for $p\geq
p^{*}_{0,\mu_{2},0}$ when $\mu_{1}=0$, $\mu_{2}\in[-\frac{(N-2)^{2}}{4},0)$
and $Q\equiv 1$.

Our interest in this paper is to obtain the nonexistence of positive solutions
of (1.8) with the parameters of $\mu_{1},\mu_{2}$ in some suitable range. Here
a function $(u)$ is a positive solution of (1.8), if $u$ satisfies the
inequalities $\mathcal{L}_{\mu_{1},\mu_{2}}u(x)\geq Q(x)u^{p}(x)\quad{\rm for\
any\ }x\in\Omega\setminus\\{0\\}.$

###### Theorem 1.2

Let $N\geq 3$,

$\mu_{1}<N-2,\quad\quad-\frac{(N-2-\mu_{1})^{2}}{4}\leq\mu_{2}<0,$

potential $Q$ verify (1.9) with $\theta>-2$ and
$p^{\\#}_{\mu_{1},\mu_{2},\theta}$ be defined in (1.10).

Then for $p\geq p^{\\#}_{\mu_{1},\mu_{2},\theta}$, problem (1.8) has no
positive solution.

This type of nonexistence is based on the classification of Poisson problem
with the CKN operator

${\mathcal{L}}_{\mu_{1},\mu_{2}}u=g\quad{\rm in}\
\,\Omega\setminus\\{0\\},\qquad u=0\quad{\rm on}\ \,\partial\Omega,$

where $\Omega$ is a bounded domain containing the origin. We provide sharp
conditions of $g$ for the existence and nonexistence of a positive solution to
the Poisson problem. The precise results see Section 3 below.

The rest of this paper is organized as follows. In Section 2, we build a
weighted distributional identity for $\Phi_{\mu_{1},\mu_{2}}$. Section 3 is
devoted to classify the isolated singular solution of the related Poisson
problem and some important estimates. Finally, we deal with the nonexistence
of super solutions for semilinear problem (1.8) in Section 4.

##  2  Fundamental solutions

Here we first remark that the derivation of
$\tau_{\pm}=\frac{(2-N+\mu_{1})\pm\sqrt{(2-N+\mu_{1})^{2}+4\mu_{2}}}{2}$ is
based on the calculation

${\mathcal{L}}_{\mu_{1},\mu_{2}}|x|^{\tau}=\Big{(}-\tau^{2}+(2-N+\mu_{1})\tau+\mu_{2}\Big{)}{|x|^{\tau-2}}$

and $\tau_{\pm}$ are zero points of $\tau(N-2-\mu_{1}+\tau)-\mu_{2}=0$ and
$\tau_{+}+\tau_{-}=2-N+\mu_{1}$. Furthermore, we observe that

$\displaystyle{\mathcal{L}}_{\mu_{1},\mu_{2}}\left(|x|^{\tau}(-\ln|x|)\right)$
$\displaystyle=$
$\displaystyle\Big{(}-\tau^{2}+(2-N+\mu_{1})\tau+\mu_{2}\Big{)}{|x|^{\tau-2}(-\ln|x|)}$
$\displaystyle-\Big{(}N-2+\mu_{1}+2\tau\Big{)}|x|^{\tau-2},$

which implies that for $\mu_{2}=-\frac{(2-N+\mu_{1})^{2}}{4}$ and
$\tau_{0}=\frac{2-N+\mu_{1}}{2}$,

${\mathcal{L}}_{\mu_{1},\mu_{2}}\left(|x|^{\tau}(-\ln|x|)\right)=0\quad{\rm
in}\ \,\mathbb{R}^{N}\setminus\\{0\\}.$

From the definition ${\mathcal{L}}^{*}_{\mu_{1},\mu_{2}}$, direct computation
shows that

$|{\mathcal{L}}^{*}_{\mu_{1},\mu_{2}}\xi(x)|\leq
c\Big{(}\|\xi\|_{C^{2}}+\frac{1}{|x|}\|\xi\|_{C^{1}}\Big{)},\quad\forall\,\xi\in
C^{2}_{c}(\mathbb{R}^{N})$ (2.1)

for some $c>0$ independent of $\xi$.

Proof of Theorem 1.1. Case 1: $\mu_{2}>-\frac{(2-N+\mu_{1})^{2}}{4}$. Let
$u_{0}=\Phi_{\mu_{1},\mu_{2}}$, then we have that $\xi\in
C^{2}(\mathbb{R}^{N})$

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}\big{(}{|x|^{\tau_{+}-\mu_{1}}}{\xi}){\mathcal{L}}_{\mu_{1}\mu_{2}}u_{0}dx$
$\displaystyle=$ $\displaystyle\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}(-\Delta)u_{0}\,(|x|^{\tau_{+}-\mu_{1}}{\xi)}dx+\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}\frac{\mu_{1}}{|x|^{2}}{x}\cdot\nabla
u_{0}(|x|^{\tau_{\tau+}-\mu_{1}}{\xi})dx$
$\displaystyle+\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}\frac{\mu_{2}}{|x|^{2}}|x|^{\tau_{+}-\mu_{1}}\xi u_{0}dx,$

where $B_{r}(z)$ is the ball with radius $r>0$ centered at $z$, particularly,
$B_{r}=B_{r}(0)$.

We note that

$\displaystyle\mu_{1}\int_{\mathbb{R}^{N}\setminus
B_{\varepsilon}}\frac{1}{|x|^{2}}{x}\cdot\nabla
u_{0}\,|x|^{\tau_{+}-\mu_{1}}{\xi}dx$ $\displaystyle=$
$\displaystyle-\mu_{1}[N+(\tau_{+}-\mu_{1}-2)]\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}|x|^{-N}\xi dx$
$\displaystyle-\mu_{1}\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}|x|^{-N}x\cdot\nabla\xi dx-\mu_{1}\int_{\partial
B_{\varepsilon}}|x|^{1-N}\xi d\omega$

and

$\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}\frac{\mu_{2}}{|x|^{2}}|x|^{\tau_{+}-\mu_{1}}\xi
u_{0}dx=\mu_{2}\int_{{\mathbb{R}}^{N}\setminus B_{\varepsilon}}|x|^{-N}\xi
dx.$

Direct computation shows

$\displaystyle\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}(-\Delta)u_{0}\,(|x|^{\tau_{+}-\mu_{1}}{\xi})dx$
$\displaystyle=$ $\displaystyle\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}\nabla{u_{0}}\nabla(|x|^{\tau_{+}-\mu_{1}}\xi)dx+\int_{\partial
B_{\varepsilon}}\frac{x\cdot\nabla{u_{0}}}{|x|}|x|^{\tau_{+}-\mu_{1}}\xi
d\omega$ $\displaystyle=$ $\displaystyle\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}{u_{0}}(-\Delta)(|x|^{\tau_{+}-\mu_{1}}\xi)dx+\int_{\partial
B_{\varepsilon}}\frac{x\cdot\nabla{u_{0}}}{|x|}|x|^{\tau_{+}-\mu_{1}}\xi
d\omega$ $\displaystyle-\int_{\partial
B_{\varepsilon}}\frac{x\cdot\nabla(|x|^{\tau_{+}-\mu_{1}}\xi)}{|x|}u_{0}d\omega,$

where $\nu=-\frac{x}{|x|}$ be the unit outside normal vector of
${\mathbb{R}}^{N}\setminus\ B_{\varepsilon}$ and $d\omega$ be the Hausdorff
measure of $\mathbb{S}^{N-1}$,

$\displaystyle\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}{u_{0}}(-\Delta)(|x|^{\tau_{+}-\mu_{1}}\xi)dx$
$\displaystyle=$ $\displaystyle\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}u_{0}(-\Delta)\xi
d\gamma_{\mu_{1},\mu_{2}}-2(\tau_{+}-\mu_{1})\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}u_{0}\frac{x}{|x|^{2}}\cdot\nabla\xi
d\gamma_{\mu_{1},\mu_{2}}$
$\displaystyle-(\tau_{+}-\mu_{1})(N-2+\tau_{+}-\mu_{1})\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}|x|^{-N}\xi dx,$ $\displaystyle\int_{\partial
B_{\varepsilon}}\frac{x\cdot\nabla{u_{0}}}{|x|}|x|^{\tau_{+}-\mu_{1}}\xi
d\omega$ $\displaystyle=$ $\displaystyle\tau_{-}\int_{\partial
B_{\varepsilon}}|x|^{1-N}\xi d\omega=\tau_{-}\int_{\partial
B_{\varepsilon}}|x|^{1-N}\big{(}\xi(0)+\nabla\xi(0)\cdot x\big{)}d\omega$
$\displaystyle\to$ $\displaystyle\tau_{-}|\mathbb{S}^{N-1}|\xi(0)\quad{\rm
as}\ \ \varepsilon\to 0^{+}$

and

$\displaystyle\int_{\partial
B_{\varepsilon}}\frac{x\cdot\nabla(|x|^{\tau_{+}-\mu_{1}}\xi)}{|x|}u_{0}d\omega$
$\displaystyle=$ $\displaystyle(\tau_{+}-\mu_{1})\int_{\partial
B_{\varepsilon}}|x|^{1-N}\xi d\omega+\int_{\partial
B_{\varepsilon}}|x|^{1-N}x\cdot\nabla\xi d\omega$ $\displaystyle=$
$\displaystyle(\tau_{+}-\mu_{1})\int_{\partial
B_{\varepsilon}}|x|^{1-N}\big{(}\xi(0)+\nabla\xi(0)\cdot
x\big{)}d\omega+\int_{\partial B_{\varepsilon}}|x|^{1-N}x\cdot\nabla\xi
d\omega$ $\displaystyle\to$
$\displaystyle(\tau_{+}-\mu_{1})|\mathbb{S}^{N-1}|\xi(0)\quad{\rm as}\ \
\varepsilon\to 0^{+},$

where $|\nabla\zeta|$ is bounded.

Now we conclude that

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}(-\Delta)u_{0}\,(|x|^{\tau_{+}-\mu_{1}}{\xi)}dx+\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}\frac{\mu_{1}}{|x|^{2}}{x}\cdot\nabla
u_{0}(|x|^{\tau_{+}-\mu_{1}}{\xi})dx$ $\displaystyle\
+\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}\frac{\mu_{2}}{|x|^{2}}|x|^{\tau_{+}-\mu_{1}}\xi u_{0}dx$
$\displaystyle=$ $\displaystyle\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}u_{0}(-\Delta)\xi
d\gamma_{\mu_{1},\mu_{2}}+(-2\tau_{+}+\mu_{1})\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}u_{0}\frac{x}{|x|^{2}}\cdot\nabla\xi
d\gamma_{\mu_{1},\mu_{2}}$
$\displaystyle-\Big{(}\mu_{1}\big{(}N+(\tau_{+}-\mu_{1}-2)\big{)}+(\tau_{+}-\mu_{1})(N-2+\tau_{+}-\mu_{1})-\mu_{2}\Big{)}\int_{\mathbb{R}^{N}\setminus
B_{\varepsilon}}|x|^{-N}\xi dx$
$\displaystyle-(\tau_{+}-\tau_{-})\int_{\partial
B_{\varepsilon}}|x|^{1-N}\big{(}\xi(0)+O(|x|)\big{)}d\omega-\int_{\partial
B_{\varepsilon}}|x|^{1-N}x\cdot\nabla\xi d\omega$ $\displaystyle=$
$\displaystyle\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}u_{0}{\mathcal{L}}^{*}_{\mu_{1},\mu_{2}}\xi
d\gamma_{\mu_{1},\mu_{2}}-(\tau_{+}-\tau_{-})\int_{\partial
B_{\varepsilon}}|x|^{1-N}\big{(}\xi(0)+O(|x|)\big{)}d\omega$
$\displaystyle-\Big{(}\tau_{+}(N-2+\tau_{+}-\mu_{1})-\mu_{2}\Big{)}\int_{\mathbb{R}^{N}\setminus
B_{\varepsilon}}|x|^{-N}\xi dx-\int_{\partial
B_{\varepsilon}}|x|^{1-N}x\cdot\nabla\xi d\omega$ $\displaystyle\to$
$\displaystyle\int_{{\mathbb{R}}^{N}}u_{0}{\mathcal{L}}^{*}_{\mu_{1},\mu_{2}}\xi
d\gamma_{\mu_{1},\mu_{2}}-(\tau_{+}-\tau_{-})|\mathbb{S}^{N-1}|\xi(0)\quad{\rm
as}\ \ \varepsilon\to 0^{+},$

where $\tau_{+}(N-2+\tau_{+}-\mu_{1})-\mu_{2}=0$,

$\tau_{+}-\tau_{-}=\sqrt{(2-N+\mu_{1})^{2}+4\mu_{2}}$

and thanks to (2.1)

$\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}u_{0}{\mathcal{L}}^{*}_{\mu_{1},\mu_{2}}\xi
d\gamma_{\mu_{1},\mu_{2}}\to\int_{{\mathbb{R}}^{N}}u_{0}{\mathcal{L}}^{*}_{\mu_{1},\mu_{2}}\xi
d\gamma_{\mu_{1},\mu_{2}}\quad{\rm as}\ \,\epsilon\to 0^{+}.$

Here we note that

$\sqrt{(2-N+\mu_{1})^{2}+4\mu_{2}}>0\ \iff\
\mu_{2}>-\frac{(2-N+\mu_{1})^{2}}{4}.$

As a consequence, we obtain that

$\int_{{\mathbb{R}}^{N}}u_{0}{\mathcal{L}}^{*}_{\mu_{1},\mu_{2}}\xi
d\gamma_{\mu_{1},\mu_{2}}=c_{\mu_{1},\mu_{2}}\xi(0),$ (2.2)

where

$c_{\mu_{1},\mu_{2}}=\Big{(}\sqrt{(2-N+\mu_{1})^{2}+4\mu_{2}}\Big{)}{|\mathbb{S}^{N-1}|}.$

Case 2: $\mu_{2}=-\frac{(2-N+\mu_{1})^{2}}{4}$. When

$\mu_{1}\in\mathbb{R}\ \ \ \ {\rm and}\
\,\mu_{2}=-\frac{(N-2)^{2}-2(N-2)\mu_{1}}{4},$

the solution $u_{0}=\Phi_{\mu_{1},\mu_{2}}$ has different form and then

$\displaystyle\int_{\partial
B_{\varepsilon}}\frac{x\cdot\nabla{u_{0}}}{|x|}|x|^{\tau_{+}-\mu_{1}}\xi
d\omega$ $\displaystyle=$ $\displaystyle\tau_{-}\int_{\partial
B_{\varepsilon}}|x|^{1-N}(-\ln|x|)\xi d\omega-\int_{\partial
B_{\varepsilon}}|x|^{1-N}\xi d\omega$

and

$\displaystyle\int_{\partial
B_{\varepsilon}}\frac{x\cdot\nabla(|x|^{\tau_{+}-\mu_{1}}\xi)}{|x|}u_{0}d\omega$
$\displaystyle=$ $\displaystyle(\tau_{+}-\mu_{1})\int_{\partial
B_{\varepsilon}}|x|^{1-N}(-\ln|x|)\xi d\omega$ $\displaystyle+\int_{\partial
B_{\varepsilon}}|x|^{1-N}(-\ln|x|)x\cdot\nabla\xi d\omega,$

then

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}(-\Delta)u_{0}\,(|x|^{\tau_{+}-\mu_{1}}{\xi)}dx+\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}\frac{\mu_{1}}{|x|^{2}}{x}\cdot\nabla
u_{0}(|x|^{\tau_{+}-\mu_{1}}{\xi})dx$ $\displaystyle\
+\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}\frac{\mu_{2}}{|x|^{2}}|x|^{\tau_{+}-\mu_{1}}\xi u_{0}dx$
$\displaystyle=$ $\displaystyle\int_{{\mathbb{R}}^{N}\setminus
B_{\varepsilon}}u_{0}{\mathcal{L}}^{*}_{\mu_{1},\mu_{2}}\xi
d\gamma_{\mu_{1},\mu_{2}}-(\tau_{+}-\tau_{-})\int_{\partial
B_{\varepsilon}}|x|^{1-N}(-ln|x|)\xi d\omega-\int_{\partial
B_{\varepsilon}}|x|^{1-N}\big{(}\xi(0)+O(|x|)\big{)}d\omega$
$\displaystyle-\int_{\partial
B_{\varepsilon}}|x|^{1-N}(-\ln|x|)x\cdot\nabla\xi d\omega$ $\displaystyle\to$
$\displaystyle\int_{{\mathbb{R}}^{N}}u_{0}{\mathcal{L}}^{*}_{\mu_{1},\mu_{2}}\xi
d\gamma_{\mu_{1},\mu_{2}}-|\mathbb{S}^{N-1}|\xi(0)\quad{\rm as}\ \
\varepsilon\to 0^{+},$

where

$c_{\mu_{1},\mu_{2}}={|\mathbb{S}^{N-1}|}\quad\ {\rm if}\ \
\mu_{2}=\frac{\mu_{1}^{2}}{4}-\frac{(2-N+\mu_{1})^{2}}{4}.$

We complete the proof. $\Box$

##  3  Basic tools

An important tool is the comparison principle.

###### Lemma 3.1

Assume that $\mu_{1}\in\mathbb{R},$ $\mu_{2}\geq-\frac{(N-2-\mu_{1})^{2}}{4}$,
$O$ is a bounded $C^{2}$ domain containing the origin and $u_{i}$ with $i=1,2$
are classical solutions of

$\left\\{\begin{array}[]{lll}\displaystyle\mathcal{L}_{\mu_{1},\mu_{2}}u_{i}=f_{i}&{\rm
in}\quad O\setminus\\{0\\},\\\\[4.2679pt] \phantom{L_{\mu}}\displaystyle
u_{i}=0&{\rm on}\ \ \partial O\end{array}\right.$ (3.1)

and

$\limsup_{|x|\to
0^{+}}u_{1}(x)\Phi_{\mu_{1},\mu_{2}}^{-1}(x)\leq\liminf_{|x|\to
0^{+}}u_{2}(x)\Phi_{\mu_{1},\mu_{2}}^{-1}(x).$

If $f_{1}\leq f_{2}$ in $O\setminus\\{0\\}$, then

$u_{1}\leq u_{2}\quad{\rm in}\quad O\setminus\\{0\\}.$

Proof. Case 1: $\mu_{2}>-\frac{(N-2-\mu_{1})^{2}}{4}$.  Let $u=u_{1}-u_{2}$
satisfy that

$\mathcal{L}_{\mu_{1},\mu_{2}}u\leq 0\quad{\rm in}\
\,O\setminus\\{0\\}\qquad{\rm and}\quad\limsup_{|x|\to
0^{+}}u(x)\Phi_{\mu_{1},\mu_{2}}^{-1}(x)\leq 0,$

then for any $\epsilon>0$, there exists $r_{\epsilon}>0$ converging to zero as
$\epsilon\to 0$ such that

$u\leq\epsilon\Phi_{\mu_{1},\mu_{2}}\quad{\rm
in}\quad\overline{B_{r_{\epsilon}}}\setminus\\{0\\}.$

We see that

$u=0<\epsilon\Phi_{\mu_{1},\mu_{2}}\quad{\rm on}\ \ \partial O,$

then by standard comparison principle for $\mu\geq\mu_{0}$, see CKN
inequality, we have that $u\leq\epsilon\Phi_{\mu_{1},\mu_{2}}$ in
$O\setminus\\{0\\}.$ By the arbitrary of $\epsilon$, we have that $u\leq 0$ in
$O\setminus\\{0\\}$.

Case 2: $\mu_{2}=-\frac{(N-2-\mu_{1})^{2}}{4}$.  Replace
$\Phi_{\mu_{1},\mu_{2}}$ in case 1 by

$w_{t_{0}}:=\Phi_{\mu_{1},\mu_{2}}+t_{0}\Gamma_{\mu_{1},\mu_{2}},$

where $t_{0}\geq 0$ be such that

$w_{t_{0}}\geq 0\quad{on}\ \partial O$

thanks to the boundedness of $O$.

Repeat the above argument of the case 1 replacing $\Phi_{\mu_{1},\mu_{2}}$ by
$w_{t_{0}}$. $\Box$

###  3.1  Poisson problems

Our second purpose in this article is to classify the isolated singular
solutions of Poisson problem with CKN operator

$\left\\{\begin{array}[]{lll}{{\mathcal{L}}}_{\mu_{1},\mu_{2}}u=f\quad{\rm
in}\ \,\,\Omega\setminus\\{0\\},\\\\[5.69054pt] \qquad\ \ u=0\quad{\rm on}\
\,\partial\Omega,\end{array}\right.$ (3.2)

where $f:\Omega\mapsto\mathbb{R}$ is a measurable function and $\Omega$ is a
bounded smooth domain contains the origin. For Poisson problem (3.2), we have
the following results.

###### Theorem 3.1

Let $\mu_{1}\in\mathbb{R}$ and

$\mu_{2}\geq-\frac{(N-2-\mu_{1})^{2}}{4}$

and $f$ be a function in $C^{\gamma}_{loc}(\overline{\Omega}\setminus\\{0\\})$
for some $\gamma\in(0,1)$.

$(i)$ Assume that

$f\in L^{1}(\Omega,d\gamma_{\mu_{1},\mu_{2}})\qquad{\rm
i.\,e.}\;\;\int_{\Omega}|f|\,d\gamma_{\mu_{1},\mu_{2}}<+\infty,$ (3.3)

then for any $k\in\mathbb{R}$ problem (3.2) has a classical solution $u_{k}$
with the asymptotic behavior

$\lim_{x\to 0}u_{k}(x)|x|^{-\tau_{-}}=k.$ (3.4)

$(ii)$ Assume that $f\geq 0$ and

$\lim_{r\to 0^{+}}\int_{\Omega\setminus
B_{r}(0)}f\,d\gamma_{\mu_{1},\mu_{2}}=+\infty,$ (3.5)

then problem (3.2) has no nonnegative solutions.

###### Theorem 3.2

Assume that

$\mu_{1}\in\mathbb{R},\quad\mu_{2}<-\frac{(N-2-\mu_{1})^{2}}{4}$

and $f$ be a nonnegative nonzero function, then problem (3.2) has no
nonnegative solutions.

Our method is to transform the CKN Poisson problem into a Hardy Poisson
problem and the existence and nonexistence could be derived by the results in
[9].

We recall the Poisson problem with Hardy-Leray operator

$\left\\{\begin{array}[]{lll}\displaystyle\mathcal{L}_{\mu}u=f\qquad{\rm
in}\quad{\Omega}\setminus\\{0\\},\\\\[4.2679pt] \phantom{L_{\mu}}\displaystyle
u=0\qquad{\rm on}\quad\partial{\Omega},\end{array}\right.$ (3.6)

where $f:\bar{\Omega}\setminus\\{0\\}\mapsto\mathbb{R}$ is a Hölder continuous
locally in $\bar{\Omega}\setminus\\{0\\}$ and

$\mathcal{L}_{\mu}=-\Delta+\frac{\mu}{|x|^{2}}$

is the Hardy-Leray operator with $\mu\geq-\frac{(N-2)^{2}}{4}$. [9] classifies
the isolated singularities of solutions to (3.6) by building the connection
with weak solutions of

$\left\\{\begin{array}[]{lll}\displaystyle\mathcal{L}_{\mu}u=f+c_{\mu}k\delta_{0}&{\rm
in}\quad{\Omega},\\\\[4.2679pt] \phantom{L_{\mu}}\displaystyle u=0&{\rm
on}\quad\partial{\Omega}\end{array}\right.$ (3.7)

in the $|x|^{\tau_{+}(\mu)}dx$-distributional sense, that is, $u\in
L^{1}(\Omega,|x|^{\tau_{+}(\mu)}dx)$ and satisfying

$\int_{\Omega}u\mathcal{L}_{\mu}^{*}(\xi)\,|x|^{\tau_{+}(\mu)}dx=\int_{\Omega}f\xi\,|x|^{\tau_{+}(\mu)}dx+c_{\mu}k\xi(0),\quad\forall\,\xi\in
C^{1.1}_{0}(\Omega),$ (3.8)

where $k\in\mathbb{R}$ and $\mathcal{L}_{\mu}^{*}u=-\Delta
u-2\tau_{+}(\mu)\frac{x}{|x|^{2}}\cdot\nabla u$ and

$\tau_{\pm}(\mu)=\frac{2-N}{2}\pm\sqrt{\frac{(N-2)^{2}}{4}+\mu}.$

It is worth noting that $\tau_{\pm}(\mu)=\tau_{\pm}(0,\mu)$ with $\mu_{1}=0$
and $\mu_{2}=\mu$.

The classification of isolated singular solutions of (3.6) states as
following:

###### Theorem 3.3

[9, Theorem 1.3] Let $f$ be a function in
$C^{\gamma}_{loc}(\overline{\Omega}\setminus\\{0\\})$ for some
$\gamma\in(0,1)$.

$(i)$ Assume that

$f\in L^{1}(\Omega,|x|^{\tau_{+}(\mu)}dx)\qquad{\rm
i.\,e.}\;\;\int_{\Omega}|f(x)|\,|x|^{\tau_{+}(\mu)}dx<+\infty,$ (3.9)

then for any $k\in\mathbb{R}$, problem (3.7) admits a unique weak solution
$u_{k}$, which is a classical solution of problem (3.6). Furthermore, if
assume more that

$\lim_{|x|\to 0^{+}}f(x)|x|^{2-\tau_{-}(\mu)}=0,$ (3.10)

then we have the asymptotic behavior

$\lim_{|x|\to 0^{+}}u_{k}(x)|x|^{-\tau_{-}(\mu)}(x)=k.$ (3.11)

$(ii)$ Assume that $f$ verifies (3.9) and $u$ is a nonnegative solution of
(3.6), then $u$ is a weak solution of (3.7) for some $k\geq 0$.

$(iii)$ Assume that $f\geq 0$ and

$\lim_{r\to 0^{+}}\int_{\Omega\setminus
B_{r}(0)}f(x)\,|x|^{\tau_{+}(\mu)}dx=+\infty,$ (3.12)

then problem (3.6) has no nonnegative solutions.

Now we are in a position to show the isolated singular solution of (3.2).

Proof of Theorem 3.1.  Let $\tau=\frac{\mu_{1}}{2}$ and $u\in
C^{2}(\Omega\setminus\\{0\\})$ be a solution of (3.2) and

$u(x)=v(x)|x|^{\tau}\quad{\rm in}\ \,\Omega\setminus\\{0\\},$

then we have that

$\displaystyle{\mathcal{L}}_{\mu_{1},\mu_{2}}(|x|^{\tau}v(x))$
$\displaystyle=$ $\displaystyle|x|^{\tau}\Big{(}-\Delta
v+(\mu_{1}-2\tau)\frac{x}{|x|^{2}}x\cdot\nabla
v+\big{(}\mu_{2}+\tau\mu_{1}-\tau(N-2+\tau)\big{)}\frac{v}{|x|^{2}}\Big{)}$
$\displaystyle=$ $\displaystyle|x|^{\tau}{{\mathcal{L}}}_{\tilde{\mu}}v,$

where

${{\mathcal{L}}}_{\tilde{\mu}}v=-\Delta v+\tilde{\mu}\frac{v}{|x|^{2}}$

with

$\tilde{\mu}=\mu_{2}+\frac{\mu_{1}^{2}}{4}-\frac{\mu_{1}}{2}(N-2).$

Hence $v$ is a classical solution of

$\left\\{\begin{array}[]{lll}{{\mathcal{L}}}_{\tilde{\mu}}v=|x|^{-\frac{\mu_{1}}{2}}f&{\rm
in}\ \,\,\Omega\setminus\\{0\\},\\\\[5.69054pt] \quad\,u=0&{\rm on}\
\,\partial\Omega,\end{array}\right.$ (3.13)

where we note that $\tilde{\mu}\geq-\frac{(N-2)^{2}}{4}$ is equivalent to

$\mu_{2}\geq-\frac{(N-2-\mu_{1})^{2}}{4}.$

It is remarkable that

$\displaystyle\tau_{\pm}(\tilde{\mu})$ $\displaystyle=$
$\displaystyle-\frac{N-2}{2}\pm\sqrt{\frac{(N-2)^{2}}{4}+\mu_{2}+\frac{\mu_{1}^{2}}{4}-\frac{\mu_{1}}{2}(N-2)}$
$\displaystyle=$
$\displaystyle\frac{2-N\pm\sqrt{(2-N+\mu_{1})^{2}+4\mu_{2}}}{2}$
$\displaystyle=$ $\displaystyle\tau_{\pm}(\mu_{1},\mu_{2})+\frac{\mu_{1}}{2},$

which implies that

$f\in L^{1}(\Omega,d\gamma_{\mu_{1},\mu_{2}})\ \iff\
|x|^{-\frac{\mu_{1}}{2}}f\in L^{1}(\Omega,|x|^{\tau_{+}(\tilde{\mu})}dx)$

and for $k\in\mathbb{R}$

$\lim_{|x|\to 0^{+}}v(x)|x|^{-\tau_{-}(\tilde{\mu})}(x)=k\ \iff\ \lim_{|x|\to
0^{+}}u(x)|x|^{-\tau_{-}(\tilde{\mu})-\frac{\mu_{1}}{2}}=\lim_{|x|\to
0^{+}}u(x)|x|^{-\tau_{-}(\mu_{1},\mu_{2})-\mu_{1}}=k.$

Therefore, the existence and nonexistence of problem (3.2) follows Theorem 3.3
part $(i)$ and $(iii)$. We complete the proof. $\Box$

Proof of Theorem 3.2.  Observe that

$\mu_{2}<-\frac{(N-2-\mu_{1})^{2}}{4}\ \iff\ \tilde{\mu}<-\frac{(N-2)^{2}}{4}$

and as shown in the proof Theorem 3.1, the nonexistence of problem (3.2)
follows by the nonexistence of positive solution of (3.13) with
$\tilde{\mu}<-\frac{(N-2)^{2}}{4}$ from [9, Proposition 5.2]. $\Box$

###  3.2 Some estimates

In order to improve the blowing up rate at the origin or decay at infinity, we
need the following estimates. Let

${\bf c}(\tau)=-\tau(N-2-\mu_{1}+\tau)+\mu_{2}>0,$ (3.14)

the map: $\tau\mapsto c(\tau)$ is concave and $\tau_{+}\geq\tau_{-}$ are the
two zero points of ${\bf c}(\tau)=0.$ There holds

${\mathcal{L}}_{\mu_{1},\mu_{2}}|\cdot|^{\tau}={\bf
c}(\tau)|\cdot|^{\tau-2}\qquad{\rm in}\ \ \mathbb{R}^{N}\setminus\\{0\\}.$
(3.15)

###### Lemma 3.2

Assume that $\mu_{1}\in\mathbb{R},$ $\mu_{2}\geq-\frac{(N-2-\mu_{1})^{2}}{4}$,
the nonnegative function $g\in C^{\beta}_{loc}(B_{r}\setminus\\{0\\})$ for
some $\beta\in(0,1)$ and there exists $\tau\in(\tau_{-},\tau_{+})$, $c>0$ and
$r_{0}>0$ such that

$g(x)\geq c|x|^{\tau-2}\quad{\rm in}\ \,B_{r}\setminus\\{0\\}.$

Let $u_{g}$ be a positive solution of problem

${\mathcal{L}}_{\mu_{1},\mu_{2}}u\geq g\quad{\rm in}\ \
B_{r}\setminus\\{0\\},\qquad u\geq 0\quad{\rm on}\ \ \partial B_{r},$ (3.16)

then there exists $c_{1}>0$ such that

$u_{g}(x)\geq c_{1}|x|^{\tau}\quad{\rm in}\ \ B_{\frac{r}{2}}.$

Proof.  For $\tau\in(\tau_{-},\tau_{+})$, we have that

${\mathcal{L}}_{\mu_{1},\mu_{2}}|x|^{\tau}={\bf c}(\tau)|x|^{\tau-2}\quad{\rm
in}\ \ \mathbb{R}^{N}\setminus\\{0\\}.$

For $r>0$, let

$w(x)=|x|^{\tau}-r^{\tau-\tau_{+}}|x|^{\tau_{+}}\quad{\rm in}\ \
\mathbb{R}^{N}\setminus\\{0\\},$

which verifies that

${\mathcal{L}}_{\mu_{1},\mu_{2}}w={\bf c}(\tau)|x|^{\tau-2}\quad{\rm in}\
B_{r}\setminus\\{0\\},\qquad w=0\ \ {\rm on}\ \ \partial B_{r},$

where ${\bf c}(\tau)>0$. Then our argument follows by Lemma 3.1. $\Box$

##  4  Nonexistence in a punctured domain

We prove the nonexistence of positive solution of (1.8) in a punctured domain
$\Omega\setminus\\{0\\}$ by contradiction, i.e. (1.8) is assumed to have a
positive solution $u_{0}$ and we will obtain a contradiction from Theorem 3.1
and Theorem 3.2. In this section, we can assume that

$B_{2}\subset\Omega$

and

$Q(x)\geq q_{0}|x|^{\theta}\quad{\rm in}\ \,B_{1}\setminus\\{0\\}.$

Proof of Theorem 1.1.  By contradiction, we assume that $u_{0}$ is a positive
super solution of (1.8) in $\Omega\setminus\\{0\\}$. Let

$v_{0}(x)=|x|^{\tau_{+}}-1\quad{\rm for}\ \,B_{1}\setminus\\{0\\},$

then

${\mathcal{L}}_{\mu_{1},\mu_{2}}v_{0}(x)=\frac{\mu_{2}}{|x|^{2}}<0$

for $\mu_{2}<0$. Here we recall that $\tau_{+}<0$ for
$-\frac{(N-2-\mu_{1})^{2}}{4}\leq\mu_{2}<0$.

Since $B_{2}\subset\Omega$ and there exists $t_{1}>0$ such that

$u_{0}\geq t_{1}\quad{\rm for}\ \ |x|=1.$

Therefore, for any $\epsilon>0$, there holds

$\liminf_{x\to
0}\big{(}u_{0}-t_{1}v_{0}+\epsilon\Phi_{\mu_{1},\mu_{2}}\big{)}(x)\Phi_{\mu_{1},\mu_{2}}^{-1}\geq\epsilon,$
$\big{(}u_{0}-t_{1}v_{0}+\epsilon\Phi_{\mu_{1},\mu_{2}}\big{)}(x)\geq
0\quad{\rm for}\ \,|x|=1$

and

${\mathcal{L}}_{\mu_{1},\mu_{2}}(u_{0}-t_{1}v_{0}+\epsilon\Phi_{\mu_{1},\mu_{2}})\geq
0\quad{\rm in}\ B_{1}\setminus\\{0\\}$

It follows by Lemma 3.1 that

$u_{0}-t_{1}v_{0}+\epsilon\Phi_{\mu_{1},\mu_{2}}\geq 0\quad{\rm in}\
B_{1}\setminus\\{0\\}$

and by the arbitrary of $\epsilon>0$, we have that

$u_{0}\geq t_{1}v_{0}\quad{\rm in}\ B_{1}\setminus\\{0\\}.$

Let

$q^{\\#}_{\mu_{1},\mu_{2},\theta}=\frac{N+\theta}{-\tau_{+}}-1.$

Then

$q^{\\#}_{\mu_{1},\mu_{2},\theta}=p^{\\#}_{\mu_{1},\mu_{2},\theta}\quad{\rm
for}\ \ \mu_{2}=-\frac{(N-2-\mu_{1})^{2}}{4}$

and

$q^{\\#}_{\mu_{1},\mu_{2},\theta}>p^{\\#}_{\mu_{1},\mu_{2},\theta}\quad{\rm
for}\ \ -\frac{(N-2-\mu_{1})^{2}}{4}<\mu_{2}<0.$

Set $g(x)=Q(x)u_{0}(x)^{p}$, then for some $d_{0}>0$

$u_{0}(x)\geq d_{0}|x|^{\tau_{+}}\quad{\rm in}\ \
B_{\frac{1}{2}}\setminus\\{0\\}.$

Part 1: Nonexistence for $-\frac{(N-2-\mu_{1})^{2}}{4}\leq\mu_{2}<0$ and
$p\geq q^{\\#}_{\mu_{1},\mu_{2},\theta}$.  Note that

$\displaystyle{\mathcal{L}}_{\mu_{1},\mu_{2}}u_{0}(x)\geq g(x)\geq
d_{0}^{p}|x|^{\theta+\tau_{+}p}\quad{\rm in}\ \
B_{\frac{1}{2}}\setminus\\{0\\},$

where $\theta+\tau_{+}p+\tau_{+}\leq-N$ and

$\lim_{r\to 0^{+}}\int_{B_{r_{0}}(0)\setminus
B_{r}(0)}Q(x)|x|^{\tau_{+}p}\Gamma_{\mu_{1},\mu_{2}}dx=+\infty$

by the fact that $p\geq q^{\\#}_{\mu_{1},\mu_{2},\theta}$. As a consequence,
we see that $u_{0}$ is a solution of

${\mathcal{L}}_{\mu_{1},\mu_{2}}u_{0}\geq g\quad{\rm in}\
\,\Omega\setminus\\{0\\},\qquad u_{0}\geq 0\quad{\rm in}\ \
\mathbb{R}^{N}\setminus\\{0\\}.$

Then there is a contradiction from Theorem 3.1 part $(ii)$.

Part 2: Nonexistence for $-\frac{(N-2-\mu_{1})^{2}}{4}<\mu_{2}<0$ and
$p\in\big{(}p^{\\#}_{\mu_{1},\mu_{2},\theta},\,q^{\\#}_{\mu_{1},\mu_{2},\theta}\big{)}$.
Let $\tau_{0}=\tau_{+}<0$, then

${\mathcal{L}}_{\mu_{1},\mu_{2}}u_{0}(x)\geq
q_{0}d_{0}^{p}|x|^{p\tau_{0}+\theta}=q_{0}d_{0}^{p}|x|^{\tau_{1}-2}\quad{\rm
in}\ \,B_{\frac{1}{2}}\setminus\\{0\\},$

where

$\tau_{1}:=p\tau_{0}+\theta+2.$

By Lemma 3.2, we have that

$u_{0}(x)\geq d_{1}|x|^{\tau_{1}}\quad{\rm in}\
\,B_{\frac{1}{2}}\setminus\\{0\\}.$

Iteratively, we recall that

$\tau_{j}:=p\tau_{j-1}+\theta+2,\quad j=1,2,\cdots.$

Note that

$\tau_{1}-\tau_{0}=(p-1)\tau_{0}+\theta+2<0$

for $p\in(p^{\\#}_{\theta,\mu},q^{\\#}_{\theta,\mu})$.

If $\tau_{1}p+\theta+2\leq\tau_{-},$ then $\tau_{1}p+\tau_{+}+\theta<-N$ and

$\displaystyle{\mathcal{L}}_{\mu_{1},\mu_{2}}u_{0}(x)\geq g(x)\geq
q_{0}d_{1}^{p}|x|^{\theta+\tau_{1}p}\quad{\rm in}\ \
B_{\frac{1}{2}}\setminus\\{0\\}$

and a contradiction comes from Theorem 3.1 part $(ii)$.

If not, we iterate above procedure. If

$\tau_{j+1}:=\tau_{j}p+\theta+2\in(\tau_{-},\tau_{+}),$

it following by Lemma 3.2 that

$u_{0}(x)\geq d_{j+1}|x|^{\tau_{j+1}},$

where

$\tau_{j+1}=p\tau_{j}+2+\theta<\tau_{j}.$

If $\tau_{j+1}p+\tau_{+}+\theta\leq-N,$ then

$\displaystyle{\mathcal{L}}_{\mu_{1},\mu_{2}}u_{0}(x)\geq
q_{0}d_{j+1}^{p}|x|^{\theta+\tau_{j+1}p}\quad{\rm in}\ \
B_{\frac{1}{2}}\setminus\\{0\\}$

then we have a contradiction and we are done.

Furthermore, this iteration must stop by finite times since

$\tau_{j}-\tau_{j-1}=p(\tau_{j-1}-\tau_{j-2})=p^{j-1}(\tau_{1}-\tau_{0})\to-\infty\
\ {\rm as}\ \,j\to+\infty.$

As a consequence, we obtain the nonexistence for the case
$p\in\big{(}p^{\\#}_{\mu_{1},\mu_{2},\theta},\,q^{\\#}_{\mu_{1},\mu_{2},\theta}\big{)}$.

Part 3: Nonexistence for $-\frac{(N-2-\mu_{1})^{2}}{4}<\mu_{2}<0$ and
$p=p^{\\#}_{\mu_{1},\mu_{2},\theta}>1$.  Note that

$\displaystyle{\mathcal{L}}_{\mu_{1},\mu_{2}}u_{0}(x)\geq
q_{0}d_{0}^{p}|x|^{\theta+\tau_{+}p}\quad{\rm in}\ \
B_{\frac{1}{2}}\setminus\\{0\\},$

where in this case

$\theta+\tau_{+}p+2=\tau_{+}.$

Here we see that for some $\sigma_{0}>0$

$\frac{1}{2}Qu_{0}^{p-1}(x)\geq\sigma_{0}|x|^{-2}\quad{\rm in}\ \
B_{\frac{1}{2}}\setminus\\{0\\}.$

If $\mu_{2}-\sigma_{0}<-\frac{(N-2-\mu_{1})^{2}}{4}$, then we have that

${\mathcal{L}}_{\mu_{1},\mu_{2}-\sigma_{0}}u_{0}(x)=f\geq 0,$

which has no positive solution by Theorem 3.2.

If $\mu_{2}-\sigma_{0}\geq-\frac{(N-2-\mu_{1})^{2}}{4}$, we can write problem
(1.8) as following

${\mathcal{L}}_{\mu_{1},\mu_{2}-\sigma_{0}}u_{0}\geq\frac{1}{2}Q(x)u_{0}^{p}\quad{\rm
in}\ \ B_{\frac{1}{2}}\setminus\\{0\\},$ (4.1)

which the critical exponent

$p^{\\#}_{\mu_{1},\mu_{2}-\sigma_{0},\theta}=1+\frac{2+\theta}{-\tau_{+}(\mu_{1},\mu_{2}-\sigma_{0})}<1+\frac{2+\theta}{-\tau_{+}(\mu_{1},\mu_{2})}=p^{\\#}_{\mu_{1},\mu_{2},\theta},$

since $t\mapsto\tau_{+}(\mu_{1},t)$ is decreasing. Thus it reduces to the
case: part 2 for (4.1) for $p=p^{\\#}_{\mu_{1},\mu_{2},\theta}$ is
supercritical and we obtain the nonexistence. $\Box$

Acknowledgements: This work is supported by the Natural Science Foundation of
China, No. 12071189, by Jiangxi Province Science Fund for Distinguished Young
Scholars, No. 20212ACB211005.

## References

  * [1] B. Abdellaoui, E. Colorado, I. Peral, Some improved Caffarelli-Kohn-Nirenberg inequalities, Calc. Var. PDE 23, 327–345 (2005)
  * [2] C. Bandle, M. Marcus, V. Moroz, Boundary singularities of solutions of semilinear elliptic equations in the half-space with a Hardy potential, Israel J. Math. 222(1), 487–514 (2017).
  * [3] M. Bhakta, D. Mukherjee, P. Nguyen, Multiplicity and uniqueness for Lane-Emden equations and systems with Hardy potential and measure data, J. Diff. Eq. 304, 29–72 (2021).
  * [4] L. Boccardo, L. Orsina, I. Peral, A remark on existence and optimal summability of solutions of elliptic problems involving Hardy potential, Discrete Contin. Dyn. Syst., 16, (2006), 513-523.
  * [5] H. Brézis, M. Marcus, Hardy’s inequalities revisited, Ann. Sc. Norm. Super. Pisa Cl. Sci. 25, 217–237 (1997).
  * [6] H. Brézis, L. Dupaigne, and A. Tesei, On a semilinear elliptic equation with inverse-square potential, Selecta Math. (N.S.) 11, 1–7 (2005).
  * [7] L. Caffarelli, R. Kohn, L. Nirenberg, First order interpolation inequalities with weights. Compos. Math. 53, 259–275 (1984).
  * [8] F. Catrina, Z.Q. Wang, On the Caffarelli-Kohn-Nirenberg inequalities: sharp constants, existence (and nonexistence), and symmetry of extremal functions, Comm. Pure Appl. Math. 54(2), 229–258 (2001).
  * [9] H. Chen, A. Quaas and F. Zhou, On nonhomogeneous elliptic equations with the Hardy-Leray potentials, J. d’Anal. Math. 144(1), 305–334 (2021).
  * [10] H. Chen, F. Zhou, Isolated singularities for elliptic equations with Hardy operator and source nonlinearity, Disc. Cont. Dyn. Syst. A 38, 2983–3002 (2018).
  * [11] H. Chen, L. Véron. Schrödinger operators with Leray-Hardy potential singular on the boundary. J. Diff. Equ. 269, 2091–2131(2020).
  * [12] H. Chen, L. Véron. Boundary singularities of semilinear elliptic equations with Leray-Hardy potential, Comm. Contemp. Math. doi:10.1142/S0219199721500516 (2021).
  * [13] J.L. Vázquez, E. Zuazua, The Hardy inequality and the asymptotic behaviour of the heat equation with an inverse-square potential, J. Funct. Anal. 173(1), 103–153 (2000).
  * [14] N. Chaudhuri, F. Cîrstea, On trichotomy of positive singular solutions associated with the Hardy-Sobolev operator, C. R. Math. Acad. Sci. Paris 347, 153–158 (2009).
  * [15] F. Cîrstea, A complete classification of the isolated singularities for nonlinear elliptic equations with inverse square potentials, Mem. Amer. Math. Soc. 227, no. 1068 (2014).
  * [16] J. Dolbeault, M. Esteban, M. Loss, G. Tarantello, On the symmetry of extremals for the Caffarelli-Kohn-Nirenberg inequalities Adv. Nonlin. Stud.  713–726 (2009)
  * [17] L. Dupaigne, A nonlinear elliptic PDE with the inverse square potential, J. Anal. Math. 86, 359–398(2002).
  * [18] E.B. Fabes, C.E. Kenig, R.P. Serapioni, The local regularity of solutions of degenerate elliptic equations, Comm. Part. Diff. Eq. 7(1), 77–116 (1982).
  * [19] M. Fall, R. Musina, Sharp nonexistence results for a linear elliptic inequality involving Hardy and Leray potentials, J. Inequal. Appl., 21 pp (2011).
  * [20] N. Ghoussoub, F. Robert, The Hardy-Schrödinger operator with interior singularity: the remaining cases, Calc. Var. Part. Diff. Eq. 56(5), Art. 149, 54 pp (2017).
  * [21] K. Gkikas, P. Nguyen, Elliptic equations with Hardy potential and gradient-dependent nonlinearity. Adv. Nonlinear Stud. 20(2), 399–435 (2020).
  * [22] K. Gkikas, P. Nguyen, Semilinear elliptic equations with Hardy potential and gradient nonlinearity, Rev. Mat. Iberoam. 36(4) 1207–1256 (2020).
  * [23] K. Gkikas, L. Véron, Boundary singularities of solutions of semilinear elliptic equations with critical Hardy potentials, Nonlinear Anal. 121, 469–540 (2015).
  * [24] B. Guerch, L. Véron. Local properties of stationary solutions of some nonlinear singular Schrödinger equations, Rev. Mat. Iberoamericana 7, 65–114 (1991).
  * [25] J. Heinonen, T. Kilpeläinen, O. Martio, Nonlinear Potential Theory of Degenerate Elliptic Equations, Dover Publications, Inc., Mineola, NY, (2006).
  * [26] M. Marcus, P. Nguyen, Moderate solutions of semilinear elliptic equations with Hardy potential, Ann. Inst. H. Poincaré - Anal. Non Linéaire 34, 69–88 (2017).
  * [27] B. Muckenhoupt, Weighted norm inequalities for the Hardy maximal function, Trans. Am. Math. Soc. 165, 207–226 (1972).
  * [28] Y. Sire, S. Terracini, S. Vita, Liouville type theorems and regularity of solutions to degenerate or singular problems part I: even solutions, Comm. Part. Diff. Eq.  Doi: 10.1080/03605302.2020.1840586 (2020).
  * [29] Z.Q. Wang, M. Willem, Caffarelli-Kohn-Nirenberg inequalities with remainder terms, J. Funct. Anal. 203(2), 550–568 (2003).


# V606 Cen: a newly formed massive contact binary in a hierarchical triple
system

Li F.-X  Yunnan Observatories, Chinese Academy of Sciences(CAS), P.O. Box 110,
650216 Kunming, P.R. China  University of Chinese Academy of Sciences, Yuquan
Road 19#, Sijingshang Block, 100049 Beijing, China  Key Laboratory of the
Structure and Evolution of Celestial Objects, CAS, Kunming 650216, P.R. China
Liao W.-P  Yunnan Observatories, Chinese Academy of Sciences(CAS), P.O. Box
110, 650216 Kunming, P.R. China  Key Laboratory of the Structure and Evolution
of Celestial Objects, CAS, Kunming 650216, P.R. China  Qian S.-B  Yunnan
Observatories, Chinese Academy of Sciences(CAS), P.O. Box 110, 650216 Kunming,
P.R. China  University of Chinese Academy of Sciences, Yuquan Road 19#,
Sijingshang Block, 100049 Beijing, China  Key Laboratory of the Structure and
Evolution of Celestial Objects, CAS, Kunming 650216, P.R. China  Fernández
Lajús E  Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional
de La Plata, Paseo del Bosque s/n,
1900 La Plata, Buenos Aires, Argentina  Instituto de Astrofísica de La Plata
(CCT La Plata $-$ CONICET/UNLP), 1900 La Plata, Argentina  Zhang J  Yunnan
Observatories, Chinese Academy of Sciences(CAS), P.O. Box 110, 650216 Kunming,
P.R. China  Key Laboratory of the Structure and Evolution of Celestial
Objects, CAS, Kunming 650216, P.R. China  Zhao E.-G  Yunnan Observatories,
Chinese Academy of Sciences(CAS), P.O. Box 110, 650216 Kunming, P.R. China
Key Laboratory of the Structure and Evolution of Celestial Objects, CAS,
Kunming 650216, P.R. China

###### Abstract

V606 Centauri (V606 Cen) is an early B-type close binary with an orbital
period of 1.4950935 d and the complete light curves are very difficult to be
observed on the ground. By analyzing the continuous light curve obtained by
TESS, we found that it is a marginal contact binary with a very low fill-out
factor of about 2%. The O-C diagram of V606 Cen is constructed for the first
time based on the 118.8-years of eclipse times. The O-C diagram has been found
to show a downward parabolic change together with a cyclic oscillation with a
semi-amplitude of 0.0545 d and a period of 88.3 yr. The downward parabolic
variation reveals a linear period decrease at a rate of
$dP/dt=-2.08\times{10^{-7}}d\cdot yr^{-1}$ that can be explained by the mass
transfer from the more massive component to the less massive one. Both the
marginal contact configuration and the continuous period decrease suggest that
V606 Cen is a newly formed contact binary via Case A mass transfer. The cyclic
change in the O-C diagram can be explained by the Light-Travel Time Effect via
the presence of a third body. The lowest mass of the tertiary companion is
determined as M3 = 4.51($\pm 0.43$)M⊙ and the tertiary is orbiting around the
central eclipsing binary in a nearly circular orbit (e=0.33). All the results
indicate that V606 Cen is a newly formed massive contact binary in a
hierarchical triple system.

binary(including multiple):close - stars:binaries:eclipsing - stars:evolution
\- stars:individual (V606 Cen)

##  1 Introduction

Massive stars play an important role in the properties of galaxies (Massey
2003). Among the early OB-type stars, binary stars are key in exploring
evolutionary process. More significantly, there is a very high fraction of
binaries in OB-type stars. For example, 56% of O-type objects are probable
close binaries, and over 70% of all massive stars will exchange mass with a
companion (Sana et al. 2012). The presence of a nearby component will alter
the evolution (de Mink et al. 2013). Therefore, we can study evolutionary
models by analyzing the evolution stages in many early-type binary stars
(Lorenzo et al. 2014; Qian et al. 2014; Polushina 2004). The massive binaries
may evolve into a contact configuration via a fast/slow Case A or Case B mass
transfer (Sybesma 1985, 1986; Pols 1994), and many of these contact binaries
will merge (Langer 2012) by common envelope evolution or form to be compact
objects and core-collapse supernovas and so on. Eclipsing binaries are
important in the determination of the fundamental, and it can determine the
fundamental stellar parameters from photometry and spectroscopy.

V606 Cen is an early B-type massive eclipsing binary, and it was found to be a
variable star by Swope (1939). It was classified as a $\beta$ Lyr-type light
curve with the period 1.495108 d, and the first ephemeris was determined
(Hertzsprung 1950). The spectral type was first obtained in 1971 (Stephenson &
Sanduleak), and the outcome was $OB^{-}$. Four years later, its spectral type
was revised as B1-2 Ib/II by Houk (1975). It was first noted as a contact
binary by Eggen (1978). The light curve was presented by Mayer et al. (1992)
who observed three minima. V606 Cen was studied in detail for the first time
by Lorenz et al. (1999), and no authors have studied it again after that. They
presented high-resolution CCD spectroscopy and the photometric light curve in
the Johnson $UBV$ bands. On the one hand, the radial velocity curve and the
spectral mass ratio $q_{spec}$ were determined. On the other hand, the
photometric solution was accomplished based on the MORO code, and the contact
configuration for V606 Cen was confirmed. The absolute parameters of the two
components were also derived. When considering the effect of radiation
pressure on early spectral type over-contact binaries, Bauer (2005) modeled
V606 Cen with radiation pressure. The results were in good agreement with that
of Lorenz et al. (1999) calculation and corresponded well to the predictions
based on detached binaries. We now have the opportunity to further study the
evolution state of V606 Cen with better data and better methods, such as the
high-precision continuous photometry data of the Transiting Exoplanet Survey
Satellite (TESS, Ricker et al. 2015) and the Digital Access to a Sky Century
at Harvard (DASCH, Grindlay et al. 2012a, 2012b) database that provides about
100 years of a photometric times series. More information about the derived
minima is described in Section 2.

For the evolution stage of V606 Cen, Lorenz et al. (1999) made the initial
discussion without consideration of the orbital. They proposed that V606 Cen
has probably evolved into a contact binary during the slow phase of case A
mass transfer. We collected minima spanning over 118.8-years, and the O-C
diagram was constructed and analyzed for the first time. Firstly, the
existence of the additional body was confirmed by the Light-Travel Time
Effect. Secondly, we can estimate the parameters of the additional body. The
formation of early massive contact binary stars is still unsolved. At present
there are relatively few known OB-type massive contact binaries (Table 1).
Moreover, the results and evolution stage of V606 Cen are different from from
these cases (e.g., evolutionary state, mass transfer, geometrical structure,
etc.). The explanations of the columns in Table 1: Column 1 is stars name.
Column 2 is period of binaries. Column 3 is mass ratio
$\mathnormal{q}=M_{less}/M_{more}$. Column 4 is the mass of more massive
components. Column 5 is the mass of less massive components. Column 6 is
degree of contact factor. Column 7 is structure type of O-C. Column 8 is
period of third body. Column 9 is lowest mass of third body. Column 10 is
rates of periodic change of binaries. Column 11 is rates of the mass transfer.
Column 12 is spectral types. Column 13 is references. In this paper, we
present a new analysis of the light curve that demonstrates that V606 Cen is
actually a triple system observed at a special point in its evolutionary
history.

Table 1: The samples of the early OB-type contact binaries

Name | Period | q | $M_{1}$ | $M_{2}$ | f | O$-$C | $P_{3}$ | $M_{3}$ | dP$/$dt | dM$/$dt | Spec. | Ref.
---|---|---|---|---|---|---|---|---|---|---|---|---
| (d) |  $M_{2}$/$M_{1}$ | ($M_{\odot}$) | ($M_{\odot}$) | % | Type | (yr) | min($M_{\odot}$) | $\times{10}^{-7}$ | $\times{10}^{-7}$ |  |
CT Tau | 0.667 | 0.983 | 14.25 | 14.01 | 99 | b | 64.6 | 1.29 |  $\cdots$ |  $\cdots$ | B1V+B1V | 1
V593 Cen | 0.755 | 0.952 | 6.2 | 5.9 | 45.4 | b | 50.9 | 4.3 |  $\cdots$ |  $\cdots$ | B5+B5 | 2
V701 Sco | 0.762 | 0.996 | 9.78 | 9.74 | 55.4 | a | 41.2 | 1.86 | +0.00118 |  $\cdots$ | B1+B1.5 | 1,3,4,5
BH Cen | 0.792 | 0.89 | 9.4 | 7.9 | 46.4 | a | 50.3 | 2.2 | +1.26 | 28 | B3+B3 | 3,6
BR Mus | 0.798 | 0.96 | $\cdots$ | $\cdots$ | 11 |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ | B3 | 7
GU Mon | 0.897 | 0.976 | 8.79 | 8.58 | 72.4 | c | 34.2 | 1.89 | -5.09 |  $\cdots$ | B1V+B1V | 1,8
MY Cam | 1.175 | 0.84 | 37.7 | 31.6 | $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ | O5.5V+O6.5V | 9
TU Mus | 1.387 | 0.652 | 16.7 | 10.4 | $\cdots$ | a | 47.73 | 1.55 | +4.0 | 4.2 | O7V$+$O8V | 10,11,12
V745 Cas | 1.410 | 0.572 | 18.31 | 10.47 | $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ | B0V+B1-2V | 13
SV Cen | 1.659 | 0.707 | 8.56* | 6.05* | 79 |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ | B1V+B6.5III | 14,15
V599 Aql | 1.849 | 0.672 | 6.1 | 4.1 | $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ | B3+B5 | 16
V382 Cyg | 1.886 | 0.742 | 26 | 19.3 | 10 | a | 47.7 | 2.56 | +4.4 | 4.3 | O7.3V$+$O7.7V | 10,17,18
LY Aur | 4.002 | 0.528 | 30 | 18.6 | $\cdots$ | a | 12.5 | 3.4 | +7.2 |  $\cdots$ | O9V+O9.5V | 3,19
V729 Cyg | 6.598 | 0.303 | 30.4 | 9.2 | $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ | O7+O8 | 20
OGLE SMC | 0.883 | 0.845 | 16.9 | 14.7 | 70 |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ | O9+O9.5-B0 | 18,21
SC10 108086 |  |  |  |  |  |  |  |  |  |  |  |
VFTS 352 | 1.120 | 0.98 | 25.6 | 25.1 | 29 |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ | O5.5V+O6.5V | 18,22
VFTS 066 | 1.140 | 0.508 | 13 | 6.6 | $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ | O9V+B0.2V | 22
OGLE-SMC | 2.206 | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | b | 12.4 |  $\cdots$ |  $\cdots$ |  $\cdots$ | O5+O6 | 23
ECL-4690 |  |  |  |  |  |  |  |  |  |  |  |

Note: *M$sin^{3}(i)$; (a) Upward parabola trend superimposed cyclic variation;
(b) Cyclic variation; (c) Downward parabola trend superimposed cyclic
variation;
Reference: (1) Yang et al. (2019); (2) Zhao et al. (2019); (3) Zhao et al.
(2014); (4) Qian et al. (2006); (5) Bell & Malcolm (1987); (6) Zhao et al.
(2018); (7) Lapasset et al. (1987); (8) Li et al. (2021); (9) Lorenzo et al.
(2014); (10) Qian et al. (2007); (11) Penny et al. (2008); (12) Menon et al.
(2021); (13) Çakırlı et al. (2014); (14) Ruciński et al. (1992); (15) Zola et
al. (2013); (16) Hill & Khalesseh (1991); (17) Harries et al. (1997); (18)
Abdul-Masih et al. (2021); (19) Stickland et al. (1994); (20) Antokhina et al.
(2016); (21) Hilditch et al. (2005); (22) Mahy et al. (2020); (23) Zasche et
al. (2017).

##  2 Data acquisition and new photometric observation

The photometric data used in this paper for V606 Cen are based upon our
observations and open available light curve databases. Photometric
observations were carried out on 2019 May 1, 4 and 2008 April 7 from the
Complejo Astronomico El Leoncito (CASLEO), San Juan, Argentina, using a 0.60 m
Helen Sawyer Hogg (HSH) telescope with a Johnson $V$-band filter. There are
three eclipse minima shown in Fig. 1. Another neglected minimum was found from
the light curves in the Optical Monitoring Camera (OMC, Alfonso-Garzón et al.
2012), as shown in the lower right panel of Fig. 1. The OMC is a small
refractor telescope with a 5 cm diameter aperture and Johnson V-filter. We
also used the online open databases from the Transiting Exoplanet Survey
Satellite (TESS, Ricker et al. 2015), the All Sky Automated Survey (ASAS,
Pojmanski 2002), the Bochum Survey of the Southern Galactic Disk (GDS,
Hackstein et al. 2015), and the Digital Access to a Sky Century at Harvard
(DASCH, Grindlay et al. 2012a, 2012b).

TESS is a space telescope project. It is equipped with four cameras with a
large field of view (24$\times$24 degree) and used to detect brighter stars.
Each strip was monitored for about 27 days. The orbital parameters were
initially derived using the TESS data, which is a continuous and unbroken
time-series photometry data. Since the period of V606 Cen is close to 1.5 d,
it is challenging to observe a complete phased light curve on the earth (e.g.,
Mayer et al. (2010) took more than two years to obtain the ground light curve
from a total of 22 observation nights). However, the space TESS project can
achieve a continuous light curve. For V606 Cen, we downloaded the images from
TESSCut (Brasseur et al. 2019). There are two segments, tess-s0011-2-2 and
tess-s0038-2-1, and their exposure time is 30 and 10 minutes, respectively. We
extracted the light curve by applying the aperture photometry method to the
candidate object and background (Dotson et al. 2019; Barentsen et al. 2019;
Lightkurve Collaboration et al. 2018; Liu et al. 2021), which has proved to be
reliable and trustworthy for eclipsing binaries. In the extraction, we set the
percentile threshold of background and aperture to 15% and 85% of total pixels
with a flux order, respectively, and the aperture radius was three pixels (Liu
et al. 2021). Both extracted phase curves are equally good, and their mean
errors are 0.0002 and 0.0003 mag. The results are shown in Fig.2.

We collected many eclipse minima from all the available databases and our
observations. Two methods were used: parabola fitting directly for continuous
observation data, and parabola fitting by reconstructing a phase with
dispersed data of more than one cycle. The photographic plates of V606 Cen
were digitized by the DASCH project, and they cover the years between 1889 and
1989. These data show large scatter and lower time resolution with a mean
error of 0.32 mag. However, they are beneficial for orbital period analysis
because of the long time span. Many authors have also used the second method
(Shi et al. 2021; Zasche et al. 2017; Qian et al. 2021), and its reliability
and accuracy have been demonstrated. There are 34 eclipse times from DASCH,
eight eclipse times from ASAS and one eclipse time from GDS, that were
determined by the method. In the ASAS database, only the data of the quality
flag A and B were adopted. The samples from ASAS, DASCH and GDS are shown in
Fig. 3. All the eclipse times are listed in Table 2.

Figure 1: The eclipse times by the CCD photometric observations from CASLEO and OMC. Figure 2: The light curve and the phased curve of TESS data. Table 2:  All the times of light minimum for V606 Cen Eclipse-times | Errors | Eclipse | Method | Source | Epoch | $O-C$ | Weight | Reference
---|---|---|---|---|---|---|---|---
HJD-2400000 |  ($\pm d$) | P/S |  |  |  | (d) |  |
15959.27684 | 0.00763 | S | PG | DASCH | -8021.5 | -0.18465 | 1 | 1
16270.24973 | 0.00821 | S | PG | DASCH | -7813.5 | -0.19120 | 1 | 1
16544.61053 | 0.00756 | P | PG | DASCH | -7630 | -0.18007 | 1 | 1
16963.23875 | 0.00839 | P | PG | DASCH | -7350 | -0.17803 | 1 | 1
17461.10923 | 0.00992 | P | PG | DASCH | -7017 | -0.17368 | 1 | 1
18015.79964 | 0.00965 | P | PG | DASCH | -6646 | -0.16296 | 1 | 1
18645.23932 | 0.00852 | P | PG | DASCH | -6225 | -0.15764 | 1 | 1
18926.31490 | 0.00992 | P | PG | DASCH | -6037 | -0.15964 | 1 | 1
19548.28716 | 0.00852 | P | PG | DASCH | -5621 | -0.14628 | 1 | 1
20826.62548 | 0.00916 | P | PG | DASCH | -4766 | -0.11290 | 1 | 1
20976.12821 | 0.00835 | P | PG | DASCH | -4666 | -0.11952 | 1 | 1
21491.94110 | 0.00685 | P | PG | DASCH | -4321 | -0.11388 | 1 | 1
22271.66405 | 0.00818 | S | PG | DASCH | -3799.5 | -0.08219 | 1 | 1
23344.40151 | 0.00843 | P | PG | DASCH | -3082 | -0.07433 | 1 | 1
24990.51767 | 0.01088 | P | PG | DASCH | -1981 | -0.05611 | 1 | 1
27165.92788 | 0.01360 | P | PG | DASCH | -526 | -0.00694 | 1 | 1
27952.35400 | $\cdots$ | P | PG | GCVS | 0 | 0.00000 | 1 | 4
28079.44022 | 0.01002 | P | PG | DASCH | 85 | 0.00327 | 1 | 1
28988.46130 | 0.00973 | P | PG | DASCH | 693 | 0.00751 | 1 | 1
29525.21128 | 0.00999 | P | PG | DASCH | 1052 | 0.01892 | 1 | 1
29857.12932 | 0.01195 | P | PG | DASCH | 1274 | 0.02621 | 1 | 1
30266.78011 | 0.01024 | P | PG | DASCH | 1548 | 0.02137 | 1 | 1
30796.04216 | 0.01459 | P | PG | DASCH | 1902 | 0.02032 | 1 | 1
31216.92583 | 0.00970 | S | PG | DASCH | 2183.5 | 0.03518 | 1 | 1
31702.09793 | 0.00963 | P | PG | DASCH | 2508 | 0.04943 | 1 | 1
32792.02629 | 0.00941 | P | PG | DASCH | 3237 | 0.05463 | 1 | 1
42912.28966 | 0.00974 | P | PG | DASCH | 10006 | 0.03010 | 1 | 1
43996.23569 | 0.00873 | P | PG | DASCH | 10731 | 0.03334 | 1 | 1
44159.18444 | 0.00807 | P | PG | DASCH | 10840 | 0.01690 | 1 | 1
44643.60631 | 0.01281 | P | PG | DASCH | 11164 | 0.02848 | 1 | 1
45213.22513 | 0.01516 | P | PG | DASCH | 11545 | 0.01667 | 1 | 1
45415.06169 | 0.00718 | P | PG | DASCH | 11680 | 0.01561 | 1 | 1
46020.55414 | 0.00817 | P | PG | DASCH | 12085 | -0.00481 | 1 | 1
46132.67760 | 0.00941 | P | PG | DASCH | 12160 | -0.01336 | 1 | 1
46341.99088 | 0.00840 | P | PG | DASCH | 12300 | -0.01317 | 1 | 1
48684.81310 | 0.00040 | P | CCD | ESO 0.5m | 13867 | -0.00246 | 10 | 2
48687.80230 | 0.00030 | P | CCD | ESO 0.5m | 13869 | -0.00345 | 10 | 2
48690.79300 | 0.00030 | P | CCD | ESO 0.5m | 13871 | -0.00294 | 10 | 2
49016.72400 | 0.00030 | P | CCD | ESO 0.5m | 14089 | -0.00232 | 10 | 2
49019.71420 | 0.00040 | P | CCD | ESO 0.5m | 14091 | -0.00231 | 10 | 2
49153.52610 | 0.00030 | S | CCD | ESO 0.5m | 14180.5 | -0.00128 | 10 | 2
49457.77940 | 0.00020 | P | CCD | ESO 0.5m | 14384 | 0.00050 | 10 | 2
51968.04600 | $\cdots$ | P | CCD | ASAS | 16063 | 0.00511 | 10 | 4
52415.08235 | 0.00147 | P | CCD | ASAS | 16362 | 0.00850 | 1 | 1
52715.59627 | 0.00203 | P | CCD | ASAS | 16563 | 0.00863 | 1 | 1
53124.5083 | 0.0003 | S | PE | SAAO | 16836.5 | 0.01259 | 10 | 3
53204.49226 | 0.00192 | P | CCD | ASAS | 16890 | 0.00905 | 1 | 1
53351.01126 | 0.00365 | P | CCD | ASAS | 16988 | 0.00888 | 1 | 1
53419.7860 | $\cdots$ | P | CCD | GCVS | 17034 | 0.00932 | 1 | 4
53575.27621 | 0.00223 | P | CCD | ASAS | 17138 | 0.00981 | 1 | 1
54153.87958 | 0.00398 | P | CCD | ASAS | 17525 | 0.01200 | 1 | 1
54331.79796 | 0.00401 | P | CCD | ASAS | 17644 | 0.01424 | 1 | 1
54563.53619 | 0.00053 | P | CCD | CASLEO 0.6m | 17799 | 0.01264 | 10 | 1
54695.11153 | 0.00118 | P | CCD | OMC | 17887 | 0.02009 | 10 | 1
54911.89675 | 0.00428 | P | CCD | ASAS | 18032 | 0.01675 | 1 | 1
56600.61013 | 0.00228 | S | CCD | GDS | 19161.5 | 0.02203 | 1 | 1
58597.28547 | 0.00248 | P | CCD | TESS | 20497 | 0.03173 | 1 | 1
58598.03302 | 0.00239 | S | CCD | TESS | 20497.5 | 0.03318 | 1 | 1
58598.78056 | 0.00267 | P | CCD | TESS | 20498 | 0.03263 | 1 | 1
58599.52811 | 0.00297 | S | CCD | TESS | 20498.5 | 0.03308 | 1 | 1
58600.27566 | 0.00321 | P | CCD | TESS | 20499 | 0.03154 | 1 | 1
58601.02320 | 0.00269 | S | CCD | TESS | 20499.5 | 0.03199 | 1 | 1
58601.77075 | 0.00223 | P | CCD | TESS | 20500 | 0.03044 | 1 | 1
58602.51830 | 0.00334 | S | CCD | TESS | 20500.5 | 0.03289 | 1 | 1
58603.26584 | 0.00273 | P | CCD | TESS | 20501 | 0.03335 | 1 | 1
58604.01339 | 0.00358 | S | CCD | TESS | 20501.5 | 0.03280 | 1 | 1
58604.76094 | 0.00268 | P | CCD | TESS | 20502 | 0.03225 | 1 | 1
58604.79310 | 0.00035 | P | CCD | CASLEO 0.6m | 20502 | 0.03216 | 10 | 1
58605.50848 | 0.00425 | S | CCD | TESS | 20502.5 | 0.03370 | 1 | 1
58606.25603 | 0.00213 | P | CCD | TESS | 20503 | 0.02915 | 1 | 1
58607.00358 | 0.00293 | S | CCD | TESS | 20503.5 | 0.03261 | 1 | 1
58607.75112 | 0.00221 | P | CCD | TESS | 20504 | 0.03206 | 1 | 1
58607.78405 | 0.00053 | P | CCD | CASLEO 0.6m | 20504 | 0.03293 | 10 | 1
58608.49867 | 0.00361 | S | CCD | TESS | 20504.5 | 0.03351 | 1 | 1
58609.24622 | 0.00267 | P | CCD | TESS | 20505 | 0.03096 | 1 | 1
58611.48886 | 0.00310 | S | CCD | TESS | 20506.5 | 0.03232 | 1 | 1
58612.23640 | 0.00219 | P | CCD | TESS | 20507 | 0.03077 | 1 | 1
58612.98395 | 0.00286 | S | CCD | TESS | 20507.5 | 0.03223 | 1 | 1
58613.73150 | 0.00196 | P | CCD | TESS | 20508 | 0.03168 | 1 | 1
58614.47904 | 0.00307 | S | CCD | TESS | 20508.5 | 0.03313 | 1 | 1
58615.22659 | 0.00237 | P | CCD | TESS | 20509 | 0.03359 | 1 | 1
58615.97414 | 0.00225 | S | CCD | TESS | 20509.5 | 0.03404 | 1 | 1
58616.72168 | 0.00163 | P | CCD | TESS | 20510 | 0.03149 | 1 | 1
58617.46923 | 0.00321 | S | CCD | TESS | 20510.5 | 0.03294 | 1 | 1
58618.21678 | 0.00226 | P | CCD | TESS | 20511 | 0.03040 | 1 | 1
58618.96433 | 0.00290 | S | CCD | TESS | 20511.5 | 0.03185 | 1 | 1
58619.71187 | 0.00202 | P | CCD | TESS | 20512 | 0.03230 | 1 | 1
58620.45942 | 0.00300 | S | CCD | TESS | 20512.5 | 0.03375 | 1 | 1
58621.20697 | 0.00230 | P | CCD | TESS | 20513 | 0.03321 | 1 | 1
58621.95451 | 0.00233 | S | CCD | TESS | 20513.5 | 0.03366 | 1 | 1
58622.70206 | 0.00274 | P | CCD | TESS | 20514 | 0.03211 | 1 | 1
58623.44961 | 0.00568 | S | CCD | TESS | 20514.5 | 0.03356 | 1 | 1
59334.40320 | 0.00041 | P | CCD | TESS | 20990 | 0.03664 | 1 | 1
59335.15108 | 0.00046 | S | CCD | TESS | 20990.5 | 0.03697 | 1 | 1
59335.89814 | 0.00043 | P | CCD | TESS | 20991 | 0.03648 | 1 | 1
59336.64635 | 0.00043 | S | CCD | TESS | 20991.5 | 0.03714 | 1 | 1
59337.39309 | 0.00046 | P | CCD | TESS | 20992 | 0.03634 | 1 | 1
59338.14125 | 0.00045 | S | CCD | TESS | 20992.5 | 0.03695 | 1 | 1
59338.88859 | 0.00041 | P | CCD | TESS | 20993 | 0.03674 | 1 | 1
59339.63652 | 0.00043 | S | CCD | TESS | 20993.5 | 0.03713 | 1 | 1
59340.38350 | 0.00042 | P | CCD | TESS | 20994 | 0.03656 | 1 | 1
59341.13175 | 0.00045 | S | CCD | TESS | 20994.5 | 0.03726 | 1 | 1
59341.87866 | 0.00041 | P | CCD | TESS | 20995 | 0.03663 | 1 | 1
59342.62668 | 0.00042 | S | CCD | TESS | 20995.5 | 0.03710 | 1 | 1
59343.37323 | 0.00049 | P | CCD | TESS | 20996 | 0.03610 | 1 | 1
59344.12154 | 0.00047 | S | CCD | TESS | 20996.5 | 0.03687 | 1 | 1
59344.86860 | 0.00042 | P | CCD | TESS | 20997 | 0.03638 | 1 | 1
59345.61705 | 0.00043 | S | CCD | TESS | 20997.5 | 0.03728 | 1 | 1
59347.85854 | 0.00046 | P | CCD | TESS | 20999 | 0.03613 | 1 | 1
59348.60686 | 0.00044 | S | CCD | TESS | 20999.5 | 0.03691 | 1 | 1
59349.35402 | 0.00041 | P | CCD | TESS | 21000 | 0.03652 | 1 | 1
59350.10176 | 0.00047 | S | CCD | TESS | 21000.5 | 0.03671 | 1 | 1
59350.84882 | 0.00043 | P | CCD | TESS | 21001 | 0.03623 | 1 | 1
59351.59714 | 0.00043 | S | CCD | TESS | 21001.5 | 0.03700 | 1 | 1
59352.34396 | 0.00041 | P | CCD | TESS | 21002 | 0.03627 | 1 | 1
59353.09200 | 0.00043 | S | CCD | TESS | 21002.5 | 0.03677 | 1 | 1
59354.58671 | 0.00046 | S | CCD | TESS | 21003.5 | 0.03638 | 1 | 1
59355.33409 | 0.00041 | P | CCD | TESS | 21004 | 0.03622 | 1 | 1
59356.08237 | 0.00046 | S | CCD | TESS | 21004.5 | 0.03695 | 1 | 1
59356.82921 | 0.00043 | P | CCD | TESS | 21005 | 0.03624 | 1 | 1
59357.57722 | 0.00045 | S | CCD | TESS | 21005.5 | 0.03671 | 1 | 1
59358.32386 | 0.00046 | P | CCD | TESS | 21006 | 0.03580 | 1 | 1
59359.07213 | 0.00047 | S | CCD | TESS | 21006.5 | 0.03652 | 1 | 1
59359.81934 | 0.00041 | P | CCD | TESS | 21007 | 0.03619 | 1 | 1
Table 2: (Continued)

References: (1) This paper; (2) Lorenz et al. (1999); (3) Mayer et al. (2006);
(4) O-C gateway 111http://var2.astro.cz/ocgate/

Figure 3: Eclipse times from three databases (left for ASAS, upper right for
DASCH, lower right for GDS). In each case, the upper panel show the times
series and the lower panel shows the orbital phase light curve with the fitted
parabola around the time of minimum.

##  3 INVESTIGATION ON PERIOD VARIATIONS

The analysis of periodic variations in the eclipse timings is a common way to
study close binaries. It may help us to detect visible or invisible additional
objects by the Light-Travel Time Effect and understand the mass transfer
between their components by the parabolic trend. V606 Cen is an early B-type
contact binary, whose eclipse timings have not been studied before this work.
We have analyzed its O-C diagram for the first time, and the data adopted span
over 118.8-years. So far, there are few studies of the early OB-type contact
binaries, which makes it difficult to confirm the evolution theory. We
collected data on the known OB-type contact binaries in Table 1. The sample
includes galactic (upper part) and extragalactic (lower part) binaries. The
O-C curve of V606 Cen with respect to the linear ephemeris was given from the
O-C gateway (see earlier footnote):

$Min.I=HJD\,2427952.354+1.4950935\times{E}.$ (1)

In the analysis of the O-C diagram, the eccentric orbit of the third celestial
body rotating around the central binary was considered first. It was
calculated based on equation (1), and the resulting O-C estimates are shown in
the upper panel of Fig. 4. The O-C variations were fit based on the Light-
Travel Time Effect by the following equations (Irwin 1952).

$O-C=\Delta{T_{0}}+\Delta{P_{0}}\times{E}+{\beta}E^{2}+{\tau},$ (2)
${\tau}=K\frac{1}{(1-e^{2}\cos^{2}\omega)^{1/2}}[\frac{1-e^{2}}{1+e\cos\nu}\sin(\nu+\omega)+e\sin\omega].$
(3)

Where ${\tau}$, ${\beta}$ are the cyclic change and the rate of the linear
period decrease. $\Delta{T_{0}}$ and $\Delta{P_{0}}$ in the formula are the
revised epoch and period. K is the semi-amplitude of O-C, e, $\nu$, $\omega$
are the eccentricity, the true anomaly and the longitude of the periastron
from the ascending node for the additional body. The parameters in equation
(2) and equation (3) were estimated using a weighted, least-squares fit with
weights of 10 and 1, the former was assigned to photoelectric (PE) and CCD
data from literature and observation except for TESS, because many eclipse-
times of TESS were obtained in a very short time. The latter was applied to
the photographic (PG), TESS and the reconstructing fitting minima.

The (O-C)1 diagram shows a noticeable cyclic oscillation with a semi-amplitude
of 0.0545 d and a period of 88.3 yr. We tried to analyze the O-C diagram when
the additional body with a circular orbit. Its (O-C)1 also shows a clear
periodic change with a small semi-amplitude in Fig.4. The residuals of both
methods are shown in the lower panel of Fig. 4. The corresponding results of
these two cases are listed in Table 3. From these two cases, their O-C
diagrams appear as a downward parabola superimposed with a periodic change.
The $\chi^{2}_{\nu}$ value of these two cases are 1.708 and 2.212, and the
$\chi^{2}_{\nu}$ value of the eccentric orbit case is smaller. In order to
detect the actual orbit parameters of the third body, more high precise minima
data need to be observed, which may be obtained in the next cycle of about 20
years in order to record more than the 1.5 cycles of the Light-Travel Time
Effect shown here.

Figure 4: The O-C diagrams of V606 Cen fitted with an eccentric orbit and a circular orbit, the data sources are listed in upper panel. The red dash lines are parabolic change. The black solid lines refer to downward parabolic plus cyclic variation in upper panel, and the cyclic variation alone is shown in the middel panel. Table 3: Orbital parameters of the third body in the V606 Cen system. Parameters | Eccentric orbit case | Circular orbit case
---|---|---
Revised epoch, $\Delta{T_{0}}$(d) | $-0.02347\,(\pm 0.00256)$ | $-0.02601\,(\pm 0.0016)$
Revised period, $\Delta{P_{0}}$(d) | $1.11\,(\pm 0.01)\times{10^{-5}}$ | $1.13\,(\pm 0.01)\times{10^{-5}}$
Semiamplitude, K(d) | $0.0545\,(\pm 0.0046)$ | $0.0441\,(\pm 0.0008)$
Orbital period, $P_{3}$ (yr) | $88.3\,(\pm 1.9)$ | $85.9\,(\pm 1.3)$
Rate of the period change, $\dot{P}({d\,yr^{-1}})$ | $-2.08\,(\pm 0.07)\times 10^{-7}$ | $-2.22\,(\pm 0.06)\times 10^{-7}$
Longitude of the periastron passage, $\omega$(deg) | 14.9 ($\pm 10.0$) | $\cdots$
Eccentricity, $e$ | $0.33\,(\pm 0.07)$ | $\cdots$
Mass function, $f(m)$ ($M_{\odot}$) | $0.126\,(\pm 0.033)$ | $0.060\,(\pm 0.016)$
Projected semi-major axis, $a_{12}\sin{i^{\prime}}$ (au) | $9.95\,(\pm 0.85)$ | $7.64\,(\pm 0.15)$
Projected tertiary mass, $M_{3}\sin{i^{\prime}}$($M_{\odot}$) | $4.51\,(\pm 0.43)$ | $3.42\,(\pm 0.08)$
$\chi^{2}_{\nu}$ | $1.708\,(\pm 0.018)$ | $2.212\,(\pm 0.034)$

##  4 Photometric Investigations with W-D Method

The Wilson-Devinney (W-D) code used to model binary light curves has been
constantly updated and improved since it came out in 1971 (Wilson & Devinney
1971; Wilson 1990, 2012; Wilson & Wyithe 2003). It has gradually become the
most commonly used tool to analyze the light curves of binary stars. Up to
now, only Lorenz et al. (1999) have analyzed the light curve of V606 Cen. They
studied the light curve data from Mayer et al. (2010). The MORO code was used
based on the W-D approach. Now we use the latest W-D version to analyze the
TESS curve with higher accuracy. It is a continuous light curve, as shown in
Fig. 2. We selected the data of tess-s0011-2-2 to study, because the
observations are high quality and cover a continuous span of 27 days.

Table 4: Photometric solutions for V606 Cen Parameters | Solution H | Complete data | TESS data
---|---|---|---
$\mathnormal{q}$(M2/M1) | 0.541 ($\pm 0.001$) | 0.54845 ($\pm 0.00074$) | 0.5263 ($\pm 0.0025$)
T1(K) | 29200 (fixed) | 29000 (fixed) | 29000 (fixed)
T2(K) | 21770 ($\pm 20$) | 22073 ($\pm 25$) | 22113 ($\pm 16$)
$\mathnormal{i}$(∘) | 87.3 ($\pm 0.1$) | 88.30 ($\pm 0.09$) | 86.44 ($\pm 0.04$)
${\Delta}$T(K) | 7430 ($\pm 20$) | 6926 ($\pm 25$) | 6886 ($\pm 16$)
R2/R1 | $\cdots$ | 0.75787 ($\pm 0.00059$) | 0.7434 ($\pm 0.0002$)
L1/(L1+L2)($\mathnormal{U}$) | 0.776 ($\pm 0.002$) | 0.77911 ($\pm 0.00095$) |  $\cdots$
L2/(L1+L2)($\mathnormal{U}$) | 0.224 ($\pm 0.002$) | 0.22089 ($\pm 0.00029$) |  $\cdots$
L3/(L1+L2+L3)($\mathnormal{U}$) | 0.1% ($\pm 0.05$) | 0.00% ($\pm 0.14$) |  $\cdots$
L1/(L1+L2)($\mathnormal{B}$) | 0.758 ($\pm 0.003$) | 0.75275 ($\pm 0.00042$) |  $\cdots$
L2/(L1+L2)($\mathnormal{B}$) | 0.242 ($\pm 0.003$) | 0.24725 ($\pm 0.00042$) |  $\cdots$
L3/(L1+L2+L3)($\mathnormal{B}$) | 1.0% ($\pm 0.8$) | 0.34% ($\pm 0.18$) |  $\cdots$
L1/(L1+L2)($\mathnormal{V}$) | 0.742 ($\pm 0.002$) | 0.74532 ($\pm 0.00038$) |  $\cdots$
L2/(L1+L2)($\mathnormal{V}$) | 0.257 ($\pm 0.002$) | 0.25468 ($\pm 0.00038$) |  $\cdots$
L3/(L1+L2+L3)($\mathnormal{V}$) | 0.2% ($\pm 0.6$) | 0.00% ($\pm 0.16$) |  $\cdots$
L1/(L1+L2)($\mathnormal{TESS}$) | $\cdots$ |  $\cdots$ | 0.74156 ($\pm 0.00012$)
L2/(L1+L2)($\mathnormal{TESS}$) | $\cdots$ |  $\cdots$ | 0.25844 ($\pm 0.00012$)
L3/(L1+L2+L3)($\mathnormal{TESS}$) | $\cdots$ |  $\cdots$ | 1.775% ($\pm 0.046$)
${\Omega}_{1}$=${\Omega}_{2}$ | 2.942 ($\pm 0.002$) | 2.96249 ($\pm 0.00067$) | 2.92263 ($\pm 0.00028$)
r1(pole) | 0.410 ($\pm 0.001$) | 0.40750 ($\pm 0.00094$) | 0.41057 ($\pm 0.00031$)
r1(side) | 0.435 ($\pm 0.001$) | 0.43219 ($\pm 0.00029$) | 0.43571 ($\pm 0.00041$)
r1(back) | 0.464 ($\pm 0.001$) | 0.46106 ($\pm 0.00020$) | 0.46419 ($\pm 0.00060$)
r2(pole) | 0.308 ($\pm 0.001$) | 0.30811 ($\pm 0.00029$) | 0.30440 ($\pm 0.00012$)
r2(side) | 0.322 ($\pm 0.001$) | 0.32193 ($\pm 0.00036$) | 0.31791 ($\pm 0.00014$)
r2(back) | 0.356 ($\pm 0.001$) | 0.35490 ($\pm 0.00061$) | 0.35069 ($\pm 0.00024$)
the degree of contact($f$) | 0.04 | 0.0165 ($\pm 0.0021$) | 0.01095 ($\pm 0.00089$)
Residual | $\cdots$ | 0.00063 | 0.00052

Although the spectral type of V606 Cen was initially considered as B1-2 IB/II
(Houk & Cowley 1975), Lorenz et al. thought that it should be B0-0.5 V. Based
on V=9.887 mag for the primary star alone and E(B-V)=0.51 (Lorenz et al. 1999)
and the distance obtained from the Gaia EDR3 (Gaia Collaboration 2020; Zari et
al. 2021), we can estimate the absolute magnitude using the equation
$M_{V}=m_{V}-5log(1000/\pi)+5-A_{V}$, where $\pi$ is the parallax in units of
mas. We found an absolute magnitude for the primary of $M_{V}$ = -3.47 mag.
This is just as expected for a B0.5 V star as listed on the online
table222http$://$www.pas.rochester.edu/$\sim$emamajek/EEM$\\_$dwarf$\\_$UBVIJHK$\\_$colors$\\_$Teff.txt
by Mamajek (based mainly on Table 5 from Pecaut & Mamajek 2013). Therefore,
the temperature of the primary star is set to be T = 29000 K, which is similar
to the temperature of solution H from Lorenz et al. (1999). The parameters for
their solution H are listed in column 2 of Table 4. Next we performed our own
fit of the photometry from Mayer et al. (2010) that was analyzed by Lorenz et
al. However, instead of rebinning the data into averages of 11 to 22
individual measurements (as done by Lorenz et al.), we used all of the
individual points in the fit. The results of this fit (column 3 of Table 4)
are the same within errors as the solution H values. This indcates that time
averaging (similar to the TESS cadence with 30 minute sampling) has no
influence on the fit.

The space TESS project data can determine more reliable fundamental parameters
by photometric solution. The specific settings in W-D calculation are as
follows: the gravity-darkening coefficients g1 = g2 = 1.0 (Lucy 1979) and the
bolometric albedos A1 = A2 = 1.0 (Ruciński 1969) were adopted. The specific
intensities for the TESS wavelength band (dependent on local temperature,
gravity, and viewing angle) were kindly provided to us from Professor Van
Hamme. The contact model (model 3) was selected. The free parameters includes
orbital inclination $\mathnormal{i}$, modified dimensionless surface potential
of star 1 ${\Omega}_{1}$ (${\Omega}_{1}$ = ${\Omega}_{2}$ for mode 3),
effective surface temperature of star 2 (T2), bandpass luminosity of star 1
(L1) and the mass ratio $\mathnormal{q}$. T1 = 29000 K was fixed, and the
initial estimate of the mass ratio, 0.527, was adopted from the radial
velocity measurements of Lorenz et al. (1999). First, we derived a convergent
solution. Then, we tried to find solutions with other model configurations
(W-D models 2, 4, and 5), but these did not converge to an acceptable
solution. Therefore, we adopted the contact binary geometry for V606 Cen. The
next step was to add the light of the third body, because of evidence of
motion about a tertiary presented in Section 3. The flux contribution of the
tertiary star $l_{3}$ was an adjustable parameter used in the final results.
The corresponding photometric parameters are listed in Table 4. The
theoretical light curves (solid line) are shown in Fig. 5 (TESS) and Fig. 6
(UBV from Mayer et al. 2010). The geometrical appearance of V606 Cen at
orbital phases 0.0 and 0.25 is displayed in Fig. 7.

The second column of Table 4 is the solution H of Lorenz et al. (1999), the
third column is the result by using the complete unbinned data of Mayer et al.
(2010) with W-D code, and the fourth column gives the results from the TESS
data. The radii are given as fractions of the semimajor axis, $R_{1}$ and
$R_{2}$ are the equivalent volume radii, and the degree of contact is $f$ =
($\Omega_{in}$-$\Omega_{star}$)/($\Omega_{in}-\Omega_{out}$), where
$\Omega_{star}$ is the modified dimensionless potential of the star surface,
$\Omega_{in}$ and $\Omega_{out}$ are the dimensionless potential of the inner
and outer Roche lobe, respectively. All the solutions suggest that V606 Cen is
a total-eclipsing shallow-contact binary ($f$ is less than 2$\%$). It is an
A-subtype (the larger star is the hotter one), and the third light
$\mathnormal{l}_{3}$ is less than 2$\%$ of the monochromatic flux. The results
also reveal that the ${\Delta}$T is lower about 500 K than the solution H, but
the more massive component is about 7000 K hotter than its less massive one.
The mass ratio M2/M1 = 0.5263 of the TESS solution is closer to the mass ratio
of $\mathnormal{q}_{spec}$ from the spectroscopic data.

Figure 5: The orange line and black open circles represent the theoretical and
the observational light curve obtained by TESS, and Cyan open circles show the
corresponding residuals Figure 6: The solid lines and open symbols represent
the theoretical and observational light curves obtained by Mayer et al., and
below are their residuals (offset for $U$, $B$ and $V$ from top to bottom).
Figure 7: The geometrical structure of V606 Cen for two orbital phases.

##  5 Discussions and conclusions

We have analyzed the periodic variation of V606 Cen for the first time based
on our observations and a variety of online photometric databases. Its O-C
diagram shows an obvious periodic oscillation, which is probably caused by the
Light-Travel Time Effect due to the presence of the additional body (Liao &
Qian 2009; Qian et al. 2013a, 2013b; Zhao et al. 2019). We constructed and
analyzed O-Cs in two cases: one is the eccentric orbit of the third body
rotating around the central binary, a variation with a semi-amplitude of
0.0545 d and a period of 88.3 yr was detected. The other is a circular orbit
solution with a semi-amplitude of 0.0441 d and a period of 85.9 yr. The value
of $\chi^{2}_{\nu}$ for these two cases are 1.708 and 2.212. We can calculate
the mass of third body using the equation (4):

${f}(m)=\frac{(M_{3}sini)^{3}}{(M_{1}+M_{2}+M_{3})^{2}}=\frac{4\pi^{2}}{GP_{3}^{2}}\times(a_{12}sini)^{3}.$
(4)

If we adopt binary masses of $M_{1}=14.7\,M_{\odot}$ (Lorenz et al. 1999),
$M_{2}=7.74\,M_{\odot}$ (obtained by $M_{2}$/$M_{1}$), then the minimum mass
of the third body is $M_{3}=4.51(\pm 0.43)\,M_{\odot}$ and $3.42(\pm
0.08)\,M_{\odot}$ respectively for the eccentric and circular cases. If the
orbital inclination of the third body is $\mathnormal{i}=90^{\circ}$, the
distance between the tertiary companion and binary system is $59.8(\pm 6.4)$
and $58.5(\pm 1.5)$ AU, respectively.

The photometric solutions imply that V606 Cen is a contact binary, and this
binary just evolved to the contact stage with a fill-out factor of contact
close to 0. The effective temperature of the more massive component is around
7000 K hotter than that of its component. A high orbital inclination
($\mathnormal{i}$ = 86.44∘) was determined. Both the orbital inclination and
the phased light curves indicate that this system is a totally eclipsing
binary. Combining with the long-term orbital change shown in O-C diagram, the
two cases show a downward parabolic change revealing a period decrease at a
rate of $dP/dt=-2.08\times{10^{-7}}d\cdot yr^{-1}$ and
$dP/dt=-2.22\times{10^{-7}}d\cdot yr^{-1}$ respectively for the eccentric and
circular solutions. Considering the conservation of mass and angular momentum,
the mean mass transfer rates can be yielded by the equation (5),

$\dot{P}/P=3\dot{M}(1/M_{1}-1/M_{2}).$ (5)

$\dot{M}=+7.58\times{10^{-7}}M_{\odot}\cdot{yr^{-1}}$ and
$+8.09\times{10^{-7}}M_{\odot}\cdot{yr^{-1}}$, respectively. The period
decrease in V606 Cen can be explained by the mass transfer from the more
massive component to the less massive one. This fact shows that the direction
of mass transfer (primary to secondary) is the opposite of what Lorenz et al.
expected for the evolutionary state of V606 Cen. For V606 Cen, its progenitor
is a semi-detached binary. The less massive component continues to gain
material from its companion, and a new contact binary star is born with both
stars filling their Roche lobes.

The absolute magnitude of the secondary is $M_{V}$ = -2.53 mag by using the
secondary to primary $V$-band flux ratio (Lorenz et al. 1999). Based on the
table provided by Mamajek, we estimate the absolute magnitude of the secondary
star from $T_{2}$ and $M_{2}$, $M_{V}$ is -2.11 and -1.93 mag. This indicates
that the secondary component is overluminous for its mass and temperature. The
reason may be that the higher temperature material was transferred from the
primary star to the secondary component, and the secondary star obtains part
of the energy, which makes the luminosity of the secondary component higher
than a normal main sequence.

We can estimate the absolute magnitude of the tertiary from the minimum mass
and the main sequence relations from Mamajek, $-0.71$ and $-0.03$ mag for the
eccentric and circular orbits, respectively. This implies that the tertiary
has a flux contribution of 5.2% and 2.9% in V-band flux for the eccentric and
circular cases, respectively. Those are all larger than the third light (less
than 2%) estimated from the light curve. Thus, the light curve analysis
appears to rule out any main sequence tertiary this bright. Is the tertiary
body a compact body or a black hole? Further observations (e.g., the maximum
angular separation, X-ray, etc.) might help establish whether or not the
tertiary is a black hole.

Table 1 lists most of the known early OB-type contact binaries that have been
investigated in detail. There are 18 contact binaries, 14 of which are located
within the Milky Way and four beyond the Galaxy. There are nine objects whose
orbital periodic changes have been studied, including one in an extragalactic
target. However, only one of all the periodic variations of these binaries, GU
Mon, is similar to that of V606 Cen. Both have the downward parabolic trend
superimposed on a periodic variation. Interestingly, these two binaries are
completely different in terms of their mass ratio and fill-out factor, i.e.,
the mass ratio is about 1 and 0.5, and the degree of contact is about 100% and
0%, respectively. Meanwhile, the contact degree of V606 Cen is the smallest in
these binaries. We also discovered that the difference in effective
temperature between the two components of V606 Cen is the largest except for
SV Cen from the spectral type. All of the evidence shows that V606 Cen is in a
special evolution stage, which provides a valuable test of the evolutionary
models of the early OB-type binaries.

According to Kepler’s third law and the parameters of photometric solution,
the absolute parameters, such as the semi-major axes $a=15.53\,R_{\odot}$,
$R_{1}=6.76\,R_{\odot}$, $R_{2}=5.05\,R_{\odot}$, are determined for V606 Cen.
It is obvious that the component stars are very close to each other. This
raises the question of how the binary could have avoided contact during the
pre-main sequence stage when their radii were much larger. The additional body
may be the key to this mystery. Due to the existence of the third body, the
two stars with a relatively distant were brought together and evolved into an
interacting binary (Naoz 2016; Qian et al. 2013c, 2014; Zhao et al. 2021;
Eggleton 2012; Wang et al. 2021).

The formation and evolution of massive contact binaries is still poorly
understood. There is some theoretical support (Sybesma 1985, 1986; Pols 1994;
Qian et al. 2013a) that these early contact binaries may be formed through
Case A mass transfer and the contact state is only a short-lived stage during
the mass transfer. Some observations are in agreement with the prediction of
the theories, such as BH Cen, V382 Cyg, TU Mus listed in Table 1. They are on
the slow phase of Case A mass transfer with increasing periods. V606 Cen is in
a unique stage of evolution different from that of these systems, and it is a
newly formed massive contact binary created by mass transfer from the primary
to the secondary. If its orbital period keeps decreasing, it will form a
merger by common envelope evolution (Langer 2012), perhaps creating the kind
of rapidly rotating object that may be a progenitor of gamma-ray burst during
its ultimate collapse.

The authors thank the teams of DASCH database, ASAS and GDS survey for
providing open data. This work is partly supported by the Chinese Natural
Science Foundation (Nos. 11933008, 11873017 and 11903076) and the Yunnan
Natural Science Foundation (No. 202001AT070051), and it is also partly
supported by Science popularization department of Guizhou Association for
Science and Technology and by Chinese Research Entity of Astronomical
Technology and Education (CREATE). New photometric CCD observations were
obtained at the Complejo Astronomico El Leoncito (CASLEO), San Juan, Argentina
with a 0.60 m telescope. This paper includes data collected by the TESS
mission. Funding for the TESS mission is provided by the NASA’s Science
Mission Directorate. Finally, we thank the referee for providing so many
helpful suggestions.

## References

  * Abdul-Masih et al. (2021) Abdul-Masih M., Sana H., Hawcroft C., Almeida L. A., Brands S. A., de Mink S. E., Justham S., et al., 2021, A&A, 651, A96
  * Alfonso-Garzón et al. (2012) Alfonso-Garzón J., Domingo A., Mas-Hesse J. M., Giménez A., 2012, A&A, 548, A79
  * Antokhina et al. (2016) Antokhina E. A., Kumsiashvili M. I., Chargeishvili K. B., 2016, Ap, 59, 68
  * Barentsen et al. (2019) Barentsen G., Hedges C. L., De Miranda Cardoso J. V., Saunders N., Gully-Santiago M. A., Cody A. M., Dotson J., et al., 2019, AAS
  * Bauer (2005) Bauer M., 2005, Ap&SS, 296, 255.
  * Bell & Malcolm (1987) Bell S. A., Malcolm G. J., 1987, MNRAS, 226, 899
  * Brasseur et al. (2019) Brasseur, C. E., Phillip, C., Fleming, S. W., et al., 2019, Astrophysics Source Code Library. ascl:1905.007
  * Çakırlı et al. (2014) Çakırlı Ö., Ibanoglu C., Sipahi E., 2014, MNRAS, 442, 1560
  * de Mink et al. (2013) de Mink S. E., Langer N., Izzard R. G., Sana H., de Koter A., 2013, ApJ, 764, 166
  * Dotson et al. (2019) Dotson J., Barentsen G., Hedges C. L., Gully-Santiago M. A., Cody A. M., Saunders N., et al., 2019, AAS
  * Eggen (1978) Eggen O. J., 1978, AJ, 83, 288
  * Eggleton (2012) Eggleton P. P., 2012, JASS, 29, 145
  * Gaia Collaboration (2020) Gaia Collaboration, 2020, yCat, I/350
  * Grindlay & Griffin (2012) Grindlay J. E., Griffin R. E., 2012a, IAUS, 285, 243
  * Grindlay et al. (2012) Grindlay J., Tang S., Los E., Servillat M., 2012b, IAUS, 285, 29
  * Hackstein et al. (2015) Hackstein M., Fein C., Haas M., Ramolla M., Pozo Nuñez F., Barr Domínguez A., Kaderhandt L., et al., 2015, AN, 336, 590
  * Harries et al. (1997) Harries T. J., Hilditch R. W., Hill G., 1997, MNRAS, 285, 277
  * Hertzsprung (1950) Hertzsprung E., 1950, AnLei, 20, 119
  * Hilditch et al. (2005) Hilditch R. W., Howarth I. D., Harries T. J., 2005, MNRAS, 357, 304
  * Hill & Khalesseh (1991) Hill G., Khalesseh B., 1991, A&A, 245, 517
  * Houk & Cowley (1975) Houk N., Cowley A. P., 1975, mcts.book
  * Irwin (1952) Irwin J. B., 1952, ApJ, 116, 211
  * Langer (2012) Langer N., 2012, ARA&A, 50, 107
  * Lapasset et al. (1987) Lapasset E., Gomez M. N., Claria J. J., 1987, RMxAA, 14, 402
  * Li et al. (2021) Li K., Xia Q.-Q., Kim C.-H., Gao X., Hu S.-M., Guo D.-F., Gao D.-Y., et al., 2021, AJ, 162, 13
  * Liao & Qian (2009) Liao W.-P., Qian S.-B., 2009, PASJ, 61, 777
  * Lightkurve Collaboration et al. (2018) Lightkurve Collaboration, Cardoso, J. V. de M., Hedges, C., et al., 2018, Astrophysics Source Code Library. ascl:1812.013
  * Liu et al. (2018) Liu, N.-P., Qian, S.-B., Huang, Y., Yuan, Z.-L., 2021, PASP (in preparation)
  * Lorenz et al. (1999) Lorenz R., Mayer P., Drechsel H., 1999, A&A, 345, 531
  * Lorenzo et al. (2014) Lorenzo J., Negueruela I., Baker A. K. F. V., García M., Simón-Díaz S., Pastor P., Méndez M. M., 2014, A&A, 572, A110
  * Lucy (1967) Lucy L. B., 1967, ZA, 65, 89
  * Mahy et al. (2020) Mahy L., Almeida L. A., Sana H., Clark J. S., de Koter A., de Mink S. E., Evans C. J., et al., 2020, A&A, 634, A119
  * Massey (2003) Massey P., 2003, ARA&A, 41, 15
  * Mayer et al. (2010) Mayer P., Božić H., Lorenz R., Drechsel H., 2010, AN, 331, 274
  * Mayer et al. (1992) Mayer P., Drechsel H., Lorenz R., 1992, IBVS, 3805, 1
  * Mayer et al. (2006) Mayer P., Wolf M., Niarchos P. G., Gazeas K. D., Manimanis V. N., Chochol D., 2006, Ap&SS, 304, 39
  * Menon et al. (2021) Menon A., Langer N., de Mink S. E., Justham S., Sen K., Szécsi D., de Koter A., et al., 2021, MNRAS, 507, 5013
  * Naoz (2016) Naoz S., 2016, ARA&A, 54, 441
  * Pecaut & Mamajek (2013) Pecaut M. J., Mamajek E. E., 2013, ApJS, 208, 9
  * Penny et al. (2008) Penny L. R., Ouzts C., Gies D. R., 2008, ApJ, 681, 554
  * Pojmanski (2002) Pojmanski G., 2002, AcA, 52, 397
  * Pols (1994) Pols O. R., 1994, A&A, 290, 119
  * Polushina (2004) Polushina T. S., 2004, A&AT, 23, 213
  * Qian et al. (2006) Qian S.-B., Liu L., Kreiner J. M., 2006, NewA, 12, 117
  * Qian et al. (2013) Qian S.-B., Liu N.-P., Li K., He J.-J., Zhu L.-Y., Zhao E. G., Wang J.-J., et al., 2013a, ApJS, 209, 13
  * Qian et al. (2013) Qian S.-B., Liu N.-P., Liao W.-P., He J.-J., Liu L., Zhu L.-Y., Wang J.-J., et al., 2013b, AJ, 146, 38
  * Qian et al. (2021) Qian, S.-B., Shi, X.-D., Liao, W.-P., He, J.-J., Liu N.-P., Zhang, J., Zhao, E.-G., et al., 2021, ApJL (in preparation)
  * Qian et al. (2014) Qian S.-B., Wang J.-J., Zhu L.-Y., Snoonthornthum B., Wang L.-Z., Zhao E. G., Zhou X., et al., 2014, ApJS, 212, 4
  * Qian et al. (2007) Qian S.-B., Yuan J.-Z., Liu L., He J.-J., Fernández Lajús E., Kreiner J. Z., 2007, MNRAS, 380, 1599
  * Qian et al. (2013) Qian S.-B., Zhang J., Wang J.-J., Zhu L.-Y., Liu L., Zhao E. G., Li L.-J., et al., 2013c, ApJS, 207, 22
  * Ricker et al. (2015) Ricker G. R., Winn J. N., Vanderspek R., Latham D. W., Bakos G. Á., Bean J. L., Berta-Thompson Z. K., et al., 2015, JATIS, 1, 014003
  * Ruciński (1969) Ruciński S. M., 1969, AcA, 19, 245
  * Ruciński et al. (1992) Ruciński S. M., Baade D., Lu W. X., Udalski A., 1992, AJ, 103, 573
  * Sana et al. (2012) Sana H., de Mink S. E., de Koter A., Langer N., Evans C. J., Gieles M., Gosset E., et al., 2012, Sci, 337, 444
  * Shi et al. (2021) Shi X.-D., Qian S.-B., Li L.-J., Liao W.-P., 2021, MNRAS, 505, 6166
  * Stephenson & Sanduleak (1971) Stephenson C. B., Sanduleak N., 1971, PW&SO, 1, 1
  * Stickland et al. (1994) Stickland D. J., Koch R. H., Pachoulakis I., Pfeiffer R. J., 1994, Obs, 114, 107
  * Swope (1939) Swope H., 1939, Harvard Annals 90, 108
  * Sybesma (1985) Sybesma C. H. B., 1985, A&A, 142, 171
  * Sybesma (1986) Sybesma C. H. B., 1986, A&A, 159, 108
  * Wang & Zhu (2021) Wang Z.-H., Zhu L.-Y., 2021, MNRAS, 507, 2804.
  * Wilson (1990) Wilson R. E., 1990, ApJ, 356, 613
  * Wilson (2012) Wilson R. E., 2012, AJ, 144, 73
  * Wilson & Devinney (1971) Wilson R. E., Devinney E. J., 1971, ApJ, 166, 605
  * Wilson & Wyithe (2003) Wilson R. E., Wyithe S. B., 2003, ASPC, 298, 313
  * Yang et al. (2019) Yang Y.-G., Yuan H.-Y., Dai H.-F., 2019, AJ, 157, 111
  * Zari et al. (2021) Zari E., Rix H.-W., Frankel N., Xiang M., Poggio E., Drimmel R., Tkachenko A., 2021, A&A, 650, A112
  * Zasche et al. (2017) Zasche P., Wolf M., Vraštil J., 2017, MNRAS, 469, 2952
  * Zhao et al. (2014) Zhao E.-G., Qian S.-B., Li L.-J., He J.-J, Liu L., Wang J.-J., Zhang J., 2014, NewA, 26, 112
  * Zhao et al. (2019) Zhao E.-G., Qian S.-B., Soonthornthum B., Poshyachinda S., Zhou X., Sarotsakulchai T., Zhang J., et al., 2019, ApJL, 871, L10
  * Zhao et al. (2018) Zhao E.-G., Qian S.-B., Zejda M., Zhang B., Zhang J., 2018, RAA, 18, 059
  * Zhao et al. (2021) Zhao E.-G., Qian S.-B., Zhou X., Li L.-J., Zhu L.-Y., Zhang L., 2021, MNRAS, 504, 5155
  * Zola et al. (2013) Zola S., Zejda M., Siwak M., de Villiers S. N., 2013, CEAB, 37, 79