# A new modified Newton iteration for computing nonnegative $Z$-eigenpairs of
nonnegative tensors

Chun-Hua Guo  Department of Mathematics and Statistics, University of Regina,
Regina, SK S4S 0A2, Canada<EMAIL_ADDRESS>,  Wen-Wei Lin  Department
of Applied Mathematics, National Yang Ming Chiao Tung University, Hsinchu 300,
Taiwan<EMAIL_ADDRESS>and  Ching-Sung Liu  Department of Applied
Mathematics, National University of Kaohsiung, Kaohsiung 811, Taiwan
<EMAIL_ADDRESS>
###### Abstract.

We propose a new modification of Newton iteration for finding some nonnegative
$Z$-eigenpairs of a nonnegative tensor. The method has local quadratic
convergence to a nonnegative eigenpair of a nonnegative tensor, under the
usual assumption guaranteeing the local quadratic convergence of the original
Newton iteration.

###### Key words and phrases:

Nonnegative tensor; Nonnegative $Z$-eigenpair; Modified Newton iteration;
Quadratic convergence.

###### 1991 Mathematics Subject Classification:

Primary 65F15; Secondary 15A50

##  1\. Introduction

A real-valued $m$th-order $n$-dimensional tensor $\mathcal{A}$ has the form

$\mathcal{A}=(A_{i_{1}i_{2}\ldots i_{m}})\text{, }\quad A_{i_{1}i_{2}\ldots
i_{m}}\in\mathbb{R}\text{, }\quad 1\leq i_{1},i_{2},\ldots,i_{m}\leq n.$

The set of all such tensors is denoted by $\mathbb{R}^{[m,n]}$, and the set of
all nonnegative tensors $\mathcal{A}\in\mathbb{R}^{[m,n]}$, for which
$A_{i_{1}i_{2}\ldots i_{m}}\geq 0$ for all $i_{1},i_{2},\ldots,i_{m}$, is
denoted by $\mathbb{R}_{+}^{[m,n]}$. We use $x_{i}$ or $(\mathbf{x})_{i}$ to
represent the $i$th element of a column vector $\mathbf{x}$.

######  Definition 1 ([9, 3]).

Let $\mathcal{A}\in\mathbb{R}^{[m,n]}$. We say that
$(\mathbf{x},\lambda)\in\left(\mathbb{R}^{n}\backslash\\{0\\}\right)\times\mathbb{R}$
is a $Z$-eigenpair (eigenvector-eigenvalue) of $\mathcal{A}$ if

(1) $\mathcal{A}\mathbf{x}^{m-1}=\lambda\mathbf{x},\quad\|\mathbf{x}\|=1,$

where
$(\mathcal{A}\mathbf{x}^{m-1})_{i}=\sum_{i_{2},\ldots,i_{m}=1}^{n}A_{ii_{2}\ldots
i_{m}}x_{i_{2}}\ldots x_{i_{m}}$ for $i=1,\ldots,n$.

The vector norm in (1) can be $2$-norm or $1$-norm. A $Z$-eigenpair is called
a $Z_{2}$-eigenpair when $2$-norm is used, and a $Z_{1}$-eigenpair when
$1$-norm is used. It is shown in [3] that, for $\mathbf{x}$ with
$\|\mathbf{x}\|_{1}=1$, $(\mathbf{x},\lambda)$ is a $Z_{1}$-eigenpair if and
only if
$\left(\frac{\mathbf{x}}{\|\mathbf{x}\|_{2}},\frac{\lambda}{\|\mathbf{x}\|_{2}^{m-2}}\right)$
is a $Z_{2}$-eigenpair. In this paper, we are interested in $Z_{1}$-eigenpairs
and a $Z_{1}$-eigenpair will be referred to as a $Z$-eigenpair or simply an
eigenpair.

A modified Newton iteration (MNI) proposed in [5] can usually find some
nonnegative eigenpairs of a nonnegative tensor using different initial
vectors. A local quadratic convergence result is proved there for positive
eigenpairs. (All nonnegative $Z$-eigenpairs are positive when the tensor is
irreducible [2]). It is observed in [5] that MNI can also compute nonnegative
eigenpairs with some zero components in the eigenvectors, but the convergence
is much slower. Recently, another modified Newton iteration (called PNI) is
proposed in [1], where a projection is used at a later stage than in [5]. PNI
can significantly speed up the convergence in computing nonnegative eigenpairs
with some zero components in the eigenvectors. In this paper, we will address
some issues associated with PNI, particularly the determination of
$\lambda_{k}$ in the sequence $(\lambda_{k},\mathbf{x}_{k})$ approximating a
nonnegative eigenpair of the tensor, and present a new modification of the
Newton iteration.

##  2\. A new modified Newton iteration

To compute an eigenpair $(\mathbf{x}_{\ast},\lambda_{\ast})$ with
$\mathbf{x}_{\ast}\geq 0$ and
$\|\mathbf{x}_{\ast}\|_{1}=\mathbf{e}^{T}\mathbf{x}_{\ast}=1$, where
$\mathbf{e}=[1,\ldots,1]^{T},$ we can apply Newton’s method to
$\mathbf{f}(\mathbf{x},\lambda)=0$, where

(2)
$\mathbf{f}(\mathbf{x},\lambda)=\left[\begin{array}[]{c}\lambda\mathbf{x}-\mathcal{A}\mathbf{x}^{m-1}\\\
\mathbf{e}^{T}\mathbf{x}-1\end{array}\right].$

The Jacobian of $\mathbf{f}(x,\lambda)$ is given by

(3) $\mathbf{Jf}(\mathbf{x},\lambda)=\left[\begin{array}[]{cc}\lambda
I-T(\mathbf{x})&\mathbf{x}\\\ \mathbf{e}^{T}&0\end{array}\right],$

where the entries of $T(\mathbf{x})$ are

$T(\mathbf{x})_{ij}=\frac{\partial}{\partial
x_{j}}\left(\mathcal{A}\mathbf{x}^{m-1}\right)_{i}.$

For a given approximation $(\widehat{\mathbf{x}}_{k},\widehat{\lambda}_{k})$
to $(\mathbf{x}_{*},\lambda_{*})$, Newton’s method determines the next
approximation $(\widehat{\mathbf{x}}_{k+1},\widehat{\lambda}_{k+1})$ by

(10)
$\displaystyle\left[\begin{array}[]{cc}\widehat{\lambda}_{k}I-T({\widehat{\mathbf{x}}}_{k})&\widehat{\mathbf{x}}_{k}\\\
{\mathbf{e}}^{T}&0\end{array}\right]\left[\begin{array}[]{c}\mathbf{d}_{k}\\\
\delta_{k}\end{array}\right]$
$\displaystyle=\left[\begin{array}[]{c}\widehat{\lambda}_{k}\widehat{\mathbf{x}}_{k}-\mathcal{A}\widehat{\mathbf{x}}_{k}^{m-1}\\\
{\mathbf{e}}^{T}\widehat{\mathbf{x}}_{k}-1\end{array}\right],$ (11)
$\displaystyle\widehat{\mathbf{x}}_{k+1}$
$\displaystyle=\widehat{\mathbf{x}}_{k}\,-\mathbf{d}_{k},$ (12)
$\displaystyle\widehat{\lambda}_{k+1}$
$\displaystyle=\widehat{\lambda}_{k}-\delta_{k},$

assuming that

(13)
$\left[\begin{array}[]{cc}\widehat{\lambda}_{k}I-T({\widehat{\mathbf{x}}}_{k})&\widehat{\mathbf{x}}_{k}\\\
{\mathbf{e}}^{T}&0\end{array}\right]$

is nonsingular.

To put our new modified Newton iteration in a proper setting, we give below a
brief description of MNI in [5] and PNI in [1]. We also mention that a Newton-
based method has been studied in [6] for computing $Z$-eigenpairs of symmetric
tensors.

When ${\mathbf{e}}^{T}\widehat{\mathbf{x}}_{k}=1$ and
$\widehat{\lambda}_{k}I-T({\widehat{\mathbf{x}}}_{k})$ is nonsingular,
$(\widehat{\mathbf{x}}_{k+1},\widehat{\lambda}_{k+1})$ is also given by ([5,
1])

(14) $\displaystyle\widehat{\mathbf{x}}_{k+1}$
$\displaystyle=\frac{1}{m-1}\left((m-2)\widehat{\mathbf{x}}_{k}+\frac{1}{\mathbf{e}^{T}\widehat{\mathbf{w}}_{k}}\widehat{\mathbf{w}}_{k}\right),$
(15) $\displaystyle\widehat{\lambda}_{k+1}$
$\displaystyle=\frac{1}{m-1}\left(\widehat{\lambda}_{k}-\frac{1}{{\mathbf{e}}^{T}\widehat{\mathbf{w}}_{k}}\right),$

where

(16)
$\widehat{\mathbf{w}}_{k}=\left(\widehat{\lambda}_{k}I-T({\widehat{\mathbf{x}}}_{k})\right)^{-1}\widehat{\mathbf{x}}_{k}.$

That $\mathbf{e}^{T}\widehat{\mathbf{w}}_{k}\neq 0$ is equivalent to the
nonsingularity of (13).

In [5], a projection is applied to the vector $\widehat{\mathbf{w}}_{k}$ used
in (14) (but not to the vector $\widehat{\mathbf{w}}_{k}$ used in (15)):

(17)
${\mathbf{w}}_{k}=\left\\{\begin{array}[]{ll}\max(\widehat{\mathbf{w}}_{k},\mathbf{0}),&\text{if
}|\max\widehat{\mathbf{w}}_{k}|>|\min\widehat{\mathbf{w}}_{k}|,\\\
\min(\widehat{\mathbf{w}}_{k},\mathbf{0}),&\text{if
}|\max\widehat{\mathbf{w}}_{k}|\leq|\min\widehat{\mathbf{w}}_{k}|,\\\
\end{array}\right.$

where $\max$ and $\min$ are taken elementwise. With $\widehat{\mathbf{w}}_{k}$
in (14) replaced by ${\mathbf{w}}_{k}$, the vector
$\widehat{\mathbf{x}}_{k+1}$ from (14) is denoted by ${\mathbf{x}}_{k+1}$.
This gives a new approximation to ${\mathbf{x}}_{\ast}$:
${\mathbf{x}}_{k+1}>0$ with ${\mathbf{e}}^{T}{\mathbf{x}}_{k+1}=1$. Note that
${\mathbf{x}}_{k+1}$ can be obtained even when
$\mathbf{e}^{T}\widehat{\mathbf{w}}_{k}=0$.

For a pair of $n$-vectors $\mathbf{v}$ and $\mathbf{w}$ with $v_{i}\neq 0$ for
each $i$, we define

(18)
$\max\left(\frac{\mathbf{w}}{\mathbf{v}}\right)=\underset{i}{\max}\left(\frac{w_{i}}{v_{i}}\right),\text{
\
}\min\left(\frac{\mathbf{w}}{\mathbf{v}}\right)=\underset{i}{\min}\left(\frac{w_{i}}{v_{i}}\right).$

A new approximation $\lambda_{k+1}$ to $\lambda_{\ast}$ is suggested in [5] to
be any value in the interval
$[\underline{\lambda}_{k+1},\overline{\lambda}_{k+1}],$ where

$\underline{\lambda}_{k+1}=\min\left(\frac{\mathcal{A}{\mathbf{x}}_{k+1}^{m-1}}{{\mathbf{x}}_{k+1}}\right),\quad\overline{\lambda}_{k+1}=\max\left(\frac{\mathcal{A}{\mathbf{x}}_{k+1}^{m-1}}{{\mathbf{x}}_{k+1}}\right),$

such that ${\lambda}_{k+1}I-T({\mathbf{x}_{k+1}})$ is not singular or nearly
singular.

The following modified Newton iteration (Algorithm 1) has been presented in
[5] for finding a nonnegative eigenpair of a nonnegative tensor $\mathcal{A}$.

  1. (1)

Given $\mathbf{x}_{0}>0$ with $\|\mathbf{x}_{0}\|_{1}=1$, and ${\sf tol}>0$.

  2. (2)

Compute
$\overline{\lambda}_{0}=\max\left(\frac{\mathcal{A}\mathbf{x}_{0}^{m-1}}{\mathbf{x}_{0}}\right)$
and
$\underline{\lambda}_{0}=\min\left(\frac{\mathcal{A}\mathbf{x}_{0}^{m-1}}{\mathbf{x}_{0}}\right)$.

  3. (3)

Choose $\lambda_{0}\in[\underline{\lambda}_{0},\overline{\lambda}_{0}]$ such
that ${\lambda}_{0}I-T({\mathbf{x}_{0}})$ is nonsingular.

  4. (4)

for $k=0,1,2,\dots$ until
$\left\|\mathcal{A}\mathbf{x}_{k}^{m-1}-{\lambda}_{k}\mathbf{x}_{k}\right\|_{1}<{\sf
tol}$.

  5. (5)

Solve the linear system
$\left({\lambda}_{k}I-T(\mathbf{x}_{k})\right)\widehat{\mathbf{w}}_{k}=\mathbf{x}_{k}$.

  6. (6)

Determine the vector $\mathbf{w}_{k}$ by (17).

  7. (7)

Compute the vector
$\widetilde{\mathbf{x}}_{k+1}=(m-2)\mathbf{x}_{k}\,+\mathbf{w}_{k}/(\mathbf{e}^{T}\mathbf{w}_{k})$.

  8. (8)

Normalize the vector $\widetilde{\mathbf{x}}_{k+1}$:
$\mathbf{x}_{k+1}=\widetilde{\mathbf{x}}_{k+1}/\|\widetilde{\mathbf{x}}_{k+1}\|_{1}$.

  9. (9)

Compute
$\overline{\lambda}_{k+1}=\max\left(\frac{\mathcal{A}\mathbf{x}_{k+1}^{m-1}}{\mathbf{x}_{k+1}}\right)$
and
$\underline{\lambda}_{k+1}=\min\left(\frac{\mathcal{A}\mathbf{x}_{k+1}^{m-1}}{\mathbf{x}_{k+1}}\right)$.

  10. (10)

Choose $\lambda_{k+1}\in[\underline{\lambda}_{k+1},\overline{\lambda}_{k+1}]$
such that ${\lambda}_{k+1}I-T({\mathbf{x}_{k+1}})$ is nonsingular.

Algorithm 1 Modified Newton iteration (MNI)

The default value for $\lambda_{0}$ in line 3 of Algorithm 1 is
$\lambda_{0}=\overline{\lambda}_{0}$. When ${\bf
e}^{T}\widehat{\mathbf{w}}_{k}\neq 0$ for $\widehat{\mathbf{w}}_{k}$ in line
5, we use (15) to get

$\widehat{\lambda}_{k+1}=\frac{1}{m-1}\left({\lambda}_{k}-\frac{1}{{\mathbf{e}}^{T}\widehat{\mathbf{w}}_{k}}\right).$

The default value for $\lambda_{k+1}$ in line 10 is then given by

$\lambda_{k+1}=\left\\{\begin{array}[]{ll}\overline{\lambda}_{k+1}&\mbox{ if
}{\bf e}^{T}\widehat{\mathbf{w}}_{k}=0\mbox{ or
}\widehat{\lambda}_{k+1}>\overline{\lambda}_{k+1},\\\
\underline{\lambda}_{k+1}&\mbox{ if
}\widehat{\lambda}_{k+1}<\underline{\lambda}_{k+1},\\\
\widehat{\lambda}_{k+1}&\mbox{ if
}\widehat{\lambda}_{k+1}\in[\underline{\lambda}_{k+1},\overline{\lambda}_{k+1}].\end{array}\right.$

If the chosen $\lambda_{k}$ is such that ${\lambda}_{k}I-T({\mathbf{x}_{k}})$
is singular or nearly singular, we can always adjust it within the interval
$[\underline{\lambda}_{k},\overline{\lambda}_{k}]$. Note that we have
$\underline{\lambda}_{k}<\overline{\lambda}_{k}$ unless
$(\mathbf{x}_{k},\lambda_{k})$ is already an eigenpair.

The PNI in [1] can be described through a comparison with MNI. In PNI, line 6
of MNI is not performed. So the $\mathbf{w}_{k}$ in line 7 is still the
$\widehat{\mathbf{w}}_{k}$ in line 5. It is assumed in [1] that ${\bf
e}^{T}\widehat{\mathbf{w}}_{k}\neq 0$. Then ${\bf
e}^{T}\widetilde{\mathbf{x}}_{k+1}=m-1>0$ and $\widetilde{\mathbf{x}}_{k+1}$
has some positive components. In PNI, all negative components of
$\widetilde{\mathbf{x}}_{k+1}$ are replaced by $0$. That is, the update
$\widetilde{\mathbf{x}}_{k+1}\leftarrow\max(\widetilde{\mathbf{x}}_{k+1},\mathbf{0})$
is performed before line 8 is performed. In this way, PNI can potentially find
nonnegative eigenvectors with some zero components much more quickly. But it
has also created a new problem: in line 9 of MNI the definition in (18) must
be modified. For a pair of nonnegative vectors $\mathbf{w}$ and
$\mathbf{v}\neq\mathbf{0}$, it is defined in [1] that

$\displaystyle\max\left(\frac{\mathbf{w}}{\mathbf{v}}\right)=\left\\{\begin{array}[]{ll}\max\left\\{\underset{i\in
S_{2}}{\max}\left(\frac{w_{i}}{v_{i}}\right),\underset{i\in
S_{1}\setminus(S_{1}\cap
S_{2})}{\max}\left(w_{i}\right)\right\\}&\quad\mbox{if
}S_{1}\setminus(S_{1}\cap S_{2})\neq\emptyset,\\\ \underset{i\in
S_{2}}{\max}\left(\frac{w_{i}}{v_{i}}\right)&\quad\mbox{if
}S_{1}\setminus(S_{1}\cap S_{2})=\emptyset,\end{array}\right.$
$\displaystyle\min\left(\frac{\mathbf{w}}{\mathbf{v}}\right)=\left\\{\begin{array}[]{ll}0&\quad\mbox{if
}S_{1}\setminus(S_{1}\cap S_{2})\neq\emptyset,\\\ \underset{i\in
S_{2}}{\min}\left(\frac{w_{i}}{v_{i}}\right)&\quad\mbox{if
}S_{1}\setminus(S_{1}\cap S_{2})=\emptyset,\end{array}\right.$

where $S_{1}$ and $S_{2}$ are the index sets of all nonzero elements of
$\mathbf{w}$ and $\mathbf{v}$, respectively. Note that the above definition
reduces to definition (18) when $\mathbf{v}>\mathbf{0}$.

In PNI [1], $\lambda_{k+1}$ is determined using

(21)
$\lambda_{k+1}=\left\\{\begin{array}[]{ll}\widehat{\lambda}_{k+1}+\beta_{k+1}\left(\overline{\lambda}_{k+1}-\widehat{\lambda}_{k+1}\right)&\mbox{
if
}\widehat{\lambda}_{k+1}\leq(\underline{\lambda}_{k+1}+\overline{\lambda}_{k+1})/2,\\\
\widehat{\lambda}_{k+1}+\beta_{k+1}\left(\underline{\lambda}_{k+1}-\widehat{\lambda}_{k+1}\right)&\mbox{
if
}\widehat{\lambda}_{k+1}>(\underline{\lambda}_{k+1}+\overline{\lambda}_{k+1})/2,\end{array}\right.$

where $0\leq\beta_{k+1}\leq 1$ is a given constant such that
$\lambda_{k+1}I-T(\mathbf{x}_{k+1})$ is nonsingular. It is suggested in [1]
that $\beta_{k+1}$ be adjusted such that
$\lambda_{k+1}\in[\underline{\lambda}_{k+1},\overline{\lambda}_{k+1}]$ to get
a better approximation $\lambda_{k+1}$ to $\lambda_{\ast}$.

To prove local quadratic convergence of PNI, it is assumed in [1] that

(22)
$\beta_{k+1}\leq\min\left\\{1,\frac{\max\left\\{\|\mathbf{x}_{k+1}-\mathbf{x}_{\ast}\|_{1},|\widehat{\lambda}_{k+1}-{\lambda}_{\ast}|\right\\}}{\max\left\\{|\overline{\lambda}_{k+1}-\widehat{\lambda}_{k+1}|,|\underline{\lambda}_{k+1}-\widehat{\lambda}_{k+1}|\right\\}}\right\\}$

in (21). It is likely that the assumption holds for a small $\beta_{k+1}$ but
it is difficult to make sure, since the unknown
$(\mathbf{x}_{\ast},{\lambda}_{\ast})$ is involved. Moreover, a small
$\beta_{k+1}$ may not ensure
$\lambda_{k+1}\in[\underline{\lambda}_{k+1},\overline{\lambda}_{k+1}]$. This
has led us to re-examine the idea of using $\underline{\lambda}_{k+1}$ and
$\overline{\lambda}_{k+1}$ to choose $\lambda_{k+1}$.

We consider the following example in [5].

######  Example 1.

Consider $\mathcal{A}\in\mathbb{R}_{+}^{[4,2]}$ defined by

$A_{1111}=1.1,\ \ A_{2222}=1.2,\ \ A_{1112}=A_{1222}=0.25,\text{ and
}A_{ijkl}=0\text{ elsewhere}.$

The tensor has three nonnegative $Z$-eigenpairs, one of them is (rounded to
four digits) $\left(\left[0.1874,0.8126\right]^{T},0.7923\right)$. We take
$\mathbf{x}$ to be vectors near the eigenvector and determine the
corresponding intervals $[\underline{\lambda},\overline{\lambda}]$:

  * •

For $\mathbf{x}=[0.19,0.81]^{T}$,
$[\underline{\lambda},\overline{\lambda}]=[0.7774,0.7873]$.

  * •

For $\mathbf{x}=[0.187,0.813]^{T}$,
$[\underline{\lambda},\overline{\lambda}]=[0.7932,0.7949]$.

  * •

For $\mathbf{x}=[0.1875,0.8125]^{T}$,
$[\underline{\lambda},\overline{\lambda}]=[0.7919,0.7922]$.

Note that the eigenvalue $0.7923$ is not in any of these intervals. This means
that adjusting $\widehat{\lambda}_{k+1}$ from Newton’s method to
$\lambda_{k+1}\in[\underline{\lambda}_{k+1},\overline{\lambda}_{k+1}]$ will
not necessarily give a better approximation to $\lambda_{\ast}$.

We will assume that $\mathbf{J}\mathbf{f}(\mathbf{x}_{\ast},\lambda_{\ast})$
is nonsingular, but ${\lambda}_{\ast}I-T(\mathbf{x}_{\ast})$ may still be
singular.

######  Example 2.

Consider $\mathcal{A}\in\mathbb{R}_{+}^{[3,3]}$ defined by

$A_{2j3}=A_{3j2}=A_{3j3}=1\text{ for }j=1,2,3,\text{ and }A_{ijk}=0\text{
elsewhere}.$

One eigenpair of $\mathcal{A}$ is
$(\mathbf{x}_{\ast},\lambda_{\ast})=([1,0,0]^{T},0)$ and

$\mathbf{J}\mathbf{f}(\mathbf{x}_{\ast},\lambda_{\ast})=\left[\begin{array}[]{rrrr}0&0&0&1\\\
0&0&-1&0\\\ 0&-1&-1&0\\\ 1&1&1&0\end{array}\right]$

is nonsingular, while ${\lambda}_{\ast}I-T(\mathbf{x}_{\ast})$ is singular.

In situations like this, it is more appropriate to determine
$(\widehat{\mathbf{x}}_{k+1},\widehat{\lambda}_{k+1})$ using (10)–(12),
instead of (14)–(15).

These considerations lead to a simple modification of Newton iteration
(Algorithm 2) for finding a nonnegative eigenpair of a nonnegative tensor
$\mathcal{A}$. For a vector $\mathbf{x}$ having both positive and negative
components, a nonnegative (and nonzero) vector is obtained as in [4] and [1]:

$\textsf{proj}({\mathbf{x}})=\frac{\max({\mathbf{x}},\mathbf{0})}{\|\max({\mathbf{x}},\mathbf{0})\|_{1}}.$

  1. (1)

Given $\mathbf{x}_{0}>0$ with $\|\mathbf{x}_{0}\|_{1}=1$, and ${\sf tol}>0$.

  2. (2)

Compute
$\lambda_{0}=\max\left(\frac{\mathcal{A}\mathbf{x}_{0}^{m-1}}{\mathbf{x}_{0}}\right).$

  3. (3)

for $k=0,1,2,\dots$ until
$\left\|\mathcal{A}\mathbf{x}_{k}^{m-1}-{\lambda}_{k}\mathbf{x}_{k}\right\|_{1}<{\sf
tol}$.

  4. (4)

Compute $(\widehat{\mathbf{x}}_{k+1},\widehat{\lambda}_{k+1})$ using
(10)–(12), with
$(\widehat{\mathbf{x}}_{k},\widehat{\lambda}_{k})=({\mathbf{x}}_{k},{\lambda}_{k})$.

  5. (5)

Set ${\mathbf{x}}_{k+1}=\textsf{proj}(\widehat{\mathbf{x}}_{k+1})$.

  6. (6)

Set $\lambda_{k+1}=\max(\widehat{\lambda}_{k+1},0)$.

Algorithm 2 Modified projected Newton iteration (MPNI)

In PNI [1], effort is made to choose $\lambda_{k}$ such that
${\lambda}_{k}I-T({\mathbf{x}_{k}})$ is nonsingular. Actually, the choice of
$\lambda_{k}$ should also ensure that

(23)
$\left[\begin{array}[]{cc}{\lambda}_{k}I-T({{\mathbf{x}}}_{k})&{\mathbf{x}}_{k}\\\
{\mathbf{e}}^{T}&0\end{array}\right]$

is nonsingular. Otherwise, division by zero occurs in the computation of the
vector $\widetilde{\mathbf{x}}_{k+1}$ since ${\bf
e}^{T}\widehat{\mathbf{w}}_{k}=0$ . Our Algorithm 2 also requires the
nonsingularity of (23). As $\lambda\to\infty$, we have

$\left|\begin{array}[]{cc}{\lambda}I-T({{\mathbf{x}}}_{k})&{\mathbf{x}}_{k}\\\
{\mathbf{e}}^{T}&0\end{array}\right|=|\lambda I-T({{\mathbf{x}}}_{k})|(-{\bf
e}^{T}(\lambda I-T({{\mathbf{x}}}_{k}))^{-1}\mathbf{x}_{k})\sim-\lambda^{n-1}$

since ${\mathbf{e}}^{T}{\mathbf{x}}_{k}=1$. So the determinant is a polynomial
of degree $n-1$. Therefore, when the matrix (23) is singular, we can choose a
small $\epsilon_{k}>0$ and use the update
$\lambda_{k}\leftarrow\lambda_{k}+\epsilon_{k}$ to ensure that the new matrix
(23) is nonsingular.

We can easily prove that Algorithm 2 (MPNI) has local quadratic convergence
under the usual assumption that guarantees the local quadratic convergence of
the original Newton iteration.

######  Theorem 1.

Let $\left(\mathbf{x}_{\ast},\lambda_{\ast}\right)$ be a nonnegative eigenpair
of a nonnegative tensor $\mathcal{A}$, with
$\mathbf{J}\mathbf{f}(\mathbf{x}_{\ast},\lambda_{\ast})$ being nonsingular,
and let $\left\\{(\mathbf{x}_{k},{\lambda}_{k})\right\\}$ be generated by
Algorithm 2. Suppose that $(\mathbf{x}_{k_{0}},{\lambda}_{k_{0}})$ is
sufficiently close to $\left(\mathbf{x}_{\ast},\lambda_{\ast}\right)$ for some
$k_{0}\geq 0$. Then $(\mathbf{x}_{k},{\lambda}_{k})$ converges to
$(\mathbf{x}_{\ast},\lambda_{\ast})$ quadratically.

###### Proof.

Since $\mathbf{J}\mathbf{f}(\mathbf{x}_{\ast},\lambda_{\ast})$ is nonsingular,
$\mathbf{J}\mathbf{f}(\mathbf{x}_{k},\lambda_{k})$ is also nonsingular when
$(\mathbf{x}_{k},{\lambda}_{k})$ is sufficiently close to
$\left(\mathbf{x}_{\ast},\lambda_{\ast}\right)$. Let
$(\widehat{\mathbf{x}}_{k+1},\widehat{\lambda}_{k+1})$ be obtained by Newton’s
method as in (10)–(12), with
$\left(\widehat{\mathbf{x}}_{k},\widehat{{\lambda}}_{k}\right)=\left(\mathbf{x}_{k},{\lambda}_{k}\right)$.
By a basic result of Newton’s method (see [7, Theorem 5.1.2] for example),
there is a constant $c$ such that

$\left\|\left[\begin{array}[]{c}\widehat{\mathbf{x}}_{k+1}\\\
\widehat{\lambda}_{k+1}\end{array}\right]-\left[\begin{array}[]{c}\mathbf{x}_{\ast}\\\
\lambda_{\ast}\end{array}\right]\right\|_{1}\leq
c\left\|\left[\begin{array}[]{c}\mathbf{x}_{k}\\\
{\lambda}_{k}\end{array}\right]-\left[\begin{array}[]{c}\mathbf{x}_{\ast}\\\
\lambda_{\ast}\end{array}\right]\right\|_{1}^{2}.$

Since $\lambda_{\ast}\geq 0$, we have
$|{\lambda}_{k+1}-\lambda_{\ast}|\leq|\widehat{\lambda}_{k+1}-\lambda_{\ast}|.$
Since ${\bf e}^{T}\mathbf{x}_{k}=1$, we have ${\bf e}^{T}\mathbf{d}_{k}=0$
from (10) and ${\bf e}^{T}\widehat{\mathbf{x}}_{k+1}={\bf
e}^{T}\mathbf{x}_{k}-{\bf e}^{T}\mathbf{d}_{k}=1$ from (11). It follows from
[8, Lemma 4] that
$\|\mathbf{x}_{k+1}-\mathbf{x}_{\ast}\|_{1}\leq\|\widehat{\mathbf{x}}_{k+1}-\mathbf{x}_{\ast}\|_{1}$.
Therefore,

$\left\|\left[\begin{array}[]{c}{\mathbf{x}}_{k+1}\\\
{\lambda}_{k+1}\end{array}\right]-\left[\begin{array}[]{c}\mathbf{x}_{\ast}\\\
\lambda_{\ast}\end{array}\right]\right\|_{1}\leq
c\left\|\left[\begin{array}[]{c}\mathbf{x}_{k}\\\
{\lambda}_{k}\end{array}\right]-\left[\begin{array}[]{c}\mathbf{x}_{\ast}\\\
\lambda_{\ast}\end{array}\right]\right\|_{1}^{2}$

for all $(\mathbf{x}_{k},{\lambda}_{k})$ sufficiently close to
$\left(\mathbf{x}_{\ast},\lambda_{\ast}\right)$. ∎∎

We remark that PNI in [1] is mathematically equivalent to MPNI here under some
favorable assumptions:

  1. (1)

$\widehat{\lambda}_{k}\geq 0$ for all $k$.

  2. (2)

$\widehat{\lambda}_{k}I-T({\bf x}_{k})$ is nonsingular for all $k$ and
$\beta_{k}=0$ is used in (22).

  3. (3)

${\bf Jf}({\bf x}_{k},\lambda_{k})$ is nonsingular for all $k$.

Therefore, the good numerical performance of PNI reported in [1] carries over
to the MPNI here. But MPNI is simpler and more robust, and its local quadratic
convergence can be easily proved without any additional assumptions such as
the one in (22). As shown in the proof of Theorem 1, MPNI has a nice feature:
the approximation $({\mathbf{x}}_{k+1},\lambda_{k+1})$ from MPNI is always
better than $(\widehat{\mathbf{x}}_{k+1},\widehat{\lambda}_{k+1})$ from the
Newton iteration (starting from the same $({\mathbf{x}}_{k},{\lambda}_{k})$).

## Acknowledgments

C.-H. Guo was supported in part by an NSERC Discovery Grant, W.-W. Lin was
supported in part by the Ministry of Science and Technology, the National
Center for Theoretical Sciences, and ST Yau Center at Chiao-Da in Taiwan, and
C.-S. Liu was supported in part by the Ministry of Science and Technology in
Taiwan.

## References

  * [1] P. Bi, W. Li, D. Liu, M. Xiao, The projected Newton iteration approach for computing the nonnegative $Z$-eigenpairs of nonnegative tensors, CSIAM Trans. Appl. Math. 2 (2021) 376–394.
  * [2] K. C. Chang, K. J. Pearson, T. Zhang, Some variational principles for $Z$-eigenvalues of nonnegative tensors, Linear Algebra Appl. 438 (2013) 4166–4182.
  * [3] K. C. Chang, T. Zhang, On the uniqueness and non-uniqueness of the positive $Z$-eigenvector for transition probability tensors, J. Math. Anal. Appl. 408 (2013) 525–540.
  * [4] D. F. Gleich, L.-H. Lim, Y. Yu, Multilinear PageRank, SIAM J. Matrix Anal. Appl. 36 (2015) 1507–1541.
  * [5] C.-H. Guo, W.-W. Lin, C.-S. Liu, A modified Newton iteration for finding nonnegative $Z$-eigenpairs of a nonnegative tensor, Numer. Algorithms 80 (2019) 595–616.
  * [6] A. Jaffe, R. Weiss, B. Nadler, Newton correction methods for computing real eigenpairs of symmetric tensors, SIAM J. Matrix Anal. Appl. 39 (2018) 1071–1094.
  * [7] C. T. Kelley, Iterative Methods for Linear and Nonlinear Equations, SIAM, Philadelphia, PA, 1995.
  * [8] D. Liu, W. Li, S.-W. Vong, Relaxation methods for solving the tensor equation arising from the higher-order Markov chains, Numer. Linear Algebra Appl. 26 (2019) e2260.
  * [9] L. Qi, Eigenvalues of a real supersymmetric tensor, J. Symb. Comput. 40 (2005) 1302–1324.


# Symbolic Comparison of Geometric Quantities in GeoGebra

Zoltán Kovács The Private University College of Education of the Diocese of
Linz
Linz, AustriaThe work was partially supported by the grant
PID2020-113192GB-I00 from the Spanish MICINN<EMAIL_ADDRESS>Bolyai
Institute, University of Szeged
Szeged, Hungary Róbert Vajda Bolyai Institute, University of Szeged
Szeged, Hungary<EMAIL_ADDRESS>
###### Abstract

Comparison of geometric quantities usually means obtaining generally true
equalities of different algebraic expressions of a given geometric figure.
Today’s technical possibilities already support symbolic proofs of a
conjectured theorem, by exploiting computer algebra capabilities of some
dynamic geometry systems as well. We introduce GeoGebra’s new feature, the
Compare command, that helps the users in experiments in planar geometry. We
focus on automatically obtaining conjectures and their proofs at the same
time, including not just equalities but inequalities too. Our contribution can
already be successfully used to support teaching geometry classes at secondary
level, by getting several well-known and some previously unpublished result
within seconds on a modern personal computer.

##  1 Introduction

Planar Euclidean geometry always played an important role in history. Euclid’s
famous book Elements, until recent times, was supposedly the second-best
selling book of all time, after the Bible. In the modern era, teaching
geometry had different focus in different countries during different
educational and political systems. But recently the wide-spread of dynamic
geometry applications allowed for today’s schools to rediscover the beauty of
this topic (see [Davis95] for a summary), either in frontal teaching, or in
small groups, or even in individual work. Dynamic geometry applications also
influence contemporary research in planar geometry—they can help finding new
theorems (see e.g. [LRV2011]).

For geometric equalities there are already integrated tools in some dynamic
geometry systems (DGS). Maybe the most well-known application is GeoGebra—it
provides an automated reasoning toolset (ART) that supports symbolic check of
equality of lengths of line segments in a geometric configuration, or, the
equality of expressions. Both features are supported by the Relation tool or
command, or the low-level commands, Prove and ProveDetails, directly [RelTool-
ADG2014]. On one hand, GeoGebra’s ART successfully proves hundreds of well-
known theorems in elementary planar geometry111See https://prover-
test.geogebra.org/job/GeoGebra_Discovery-
provertest/69/artifact/fork/geogebra/test/scripts/benchmark/prover/html/all.html
for a recent benchmark of 291 test cases., and, on the other hand, there are
examples of obtaining automated proofs of more difficult, much recent
statements222See e.g. https://matek.hu/zoltan/blog-20201229.php., or of
completely new results333See e.g. [LNAI11006-isoptics]..

In this paper we describe our current work that goes one step forward. In our
work we use partly the same algebro-geometric theory that plays a crucial role
in GeoGebra’s toolset, based on the revolutionary work of Wu [Wu78], Chou
[Chou88], and improved later by Recio and Vélez [RecioVelez99] with
elimination theory. This is the required algebraic basis for gaining a
conjecture on a possibly fixed ratio of two quantities. Once the elimination
method is able to provide the exact ratio, this conjecture is in fact a
mathematically proven proposition.

On the other hand, we use partly general purpose real quantifier elimination
(RQE) methods to find the best possible geometric constants in the related
inequalities if two expressions do not have a fixed ratio. This approach
supports studying the class of general non-degenerate triangles, or a more
specific class of triangles such as isosceles triangles.

In this study we do not go into the further detail on technical difficulties,
but refer to two recent papers: [mcs2020] explains the mathematical background
on obtaining ratios in regular polygons, while [scsc-2020] summarizes the RQE
related issues, and points to a large set of benchmarks based on our tool.
Instead, our current paper is rather focusing on the practical use: how our
work can be fruitful for the student, the teacher and the researcher. We
present an experimental version of GeoGebra, namely, GeoGebra Discovery, which
tries to solve the problem with elimination first, and if this step is
unsuccessful, then the given RQE problem setting will be outsourced to an
external tool realgeom. We highlight that realgeom can work together with
Wolfram’s Mathematica [Mathematica] (in this case the end user is expected to
have a Mathematica subscription, or to own a Raspberry Pi system that offers
free access to Wolfram’s tool, or download a free developer command-line
version of Mathematica called WolframScript), or the free system Tarski
[ValeEnriquez-Brown] (which outsources some computations to the QEPCAD B
system [Brown03]).

Our experimental system is already capable of solving a large set of open
questions in planar Euclidean geometry. We believe that our contribution to
automated reasoning can be fruitful for several levels of mathematics
education. What is more, researchers can double check new conjectures and may
obtain even some non-trivial or yet unknown results.

We emphasize that the mathematical methods we use in this paper are mostly
well-known, however, a full implementation of the various techniques in a
single graphical application, being freely available for millions of potential
users, is completely new.

##  2 Comparison of expressions at school

School curriculum usually includes several relationships between expressions
in a planar construction. Every student has to learn the Pythagorean theorem
which is an equality of two expressions, namely $a^{2}+b^{2}$ and $c^{2}$
where $a$, $b$ and $c$ correspond to the lengths of sides in a right triangle.
In general, however, a simple question is not discussed: what happens if we
omit the assumption that the triangle has a right angle? By using our tool it
can be quickly shown mechanically that in general the inequality
$\displaystyle{a^{2}+b^{2}>\frac{c^{2}}{2}}$ holds.

In practical life everybody solves inequalities. Given the task that somebody
from point $A$ has to reach point $B$ as far as possible, no external point
$C\not\in AB$ will be chosen to visit that point first and just then
$B$—unless the terrain is not planar. Naturally, the shortest path is “as the
crow flies”, or, in mathematical means: the triangle inequality
$|AB|\leq|AC|+|BC|$ always holds.

Mathematics curriculum seems to have less to do with inequalities. This is
especially true for the subtopic of geometry. We rather prefer relationships
that express equality. To name a few, the intercept theorem or the geometric
mean theorem can be mentioned (the latter is a special case of the chord
theorem), or, at a higher level, Ptolemy’s theorem, the Ceva-Menelaus formulas
or Heron’s formula may appear in certain school types. For inequalities we can
think of the well-known correspondence between the angles and sides in a
triangle (namely, $a\leq b\Leftrightarrow\alpha\leq\beta$), or, as a kind of
generalization of Thales’ circle theorem, the inscribed angle theorem gives a
suggestion how to compare angles lying on different arcs that share the same
chord. From advanced topics we can mention Euler’s inequality (1765, but found
earlier by Chapple in 1746, [Chapple1746]) between the radii of the
circumcircle and incircle of a triangle, or the Erdős-Mordell inequality
(1937, [Mordell1937]).

These examples are well-known theorems. In fact, comparison of given
quantities is a more general task, and can be a very rich field of new
experiments, exercises, problems or even new theorems. On the other hand,
equality-based problems are much more usual in the literature. For example, by
having a deeper look on the problem lists of International Mathematical
Olympiad for the last 20 years, very few geometrical proofs of inequalities
were given, even if one-third of the problem settings are traditionally
related to planar Euclidean geometry. (See https://www.imo-
official.org/problems.aspx for a complete list. In particular, problems
2001/1, 2002/6, 2006/1, 2013/2 and 2020/6 are related to geometrical
inequalities.)

We need to admit that a certain type of inequalities has already been studied
quite exhaustively during the last decades. Here the main reference is,
without any doubt, [Bottema69]. In a part of that work several inequalities of
homogeneous type are considered, like

$3\cdot(ab+bc+ca)\leq(a+b+c)^{2}<4\cdot(ab+bc+ca)$ (Bottema 1.1)

where $a$, $b$ and $c$ denote the lengths of the sides of an arbitrary non-
degenerate triangle in the plane. In general we can reformulate the above
result to an open question, namely: Are there sharp constants $m,M$ such that

$m\cdot(ab+bc+ca)\underset{(=)}{<}(a+b+c)^{2}\underset{(=)}{<}M\cdot(ab+bc+ca)?$
(1)

The answer is, according to [Bottema69], yes, namely $m=3$ and $M=4$.

To mechanically give an exact answer we already pointed to the RQE approach.
But we will also focus on the “$m=M$” case in our algorithm—that is, we assume
that for the input expressions $w_{1}$ and $w_{2}$ there exists a positive
number $\mu$ such that

$\mu\cdot w_{1}=w_{2}$ (2)

always holds, for (almost) all possible positions of the free points in the
input figure. (In our example above $w_{1}=ab+bc+ca$ and $w_{2}=(a+b+c)^{2}$.)
Of course, this is often not the case, for example, in (1) not, since
$\mu\in[3,4)$, that is, $m=3$ and $M=4$.

Here we highlight that the case “$m=M$” can usually be found by utilizing
elimination theory as well. This remark can be extremely important to reduce
the computational complexity of the mechanical solution in many cases. In the
next section, after discussing an example when $m\neq M$ (Fig. 1) we will show
a non-trivial example for the case $m=M$ (Fig. 2, 3(a) and 3(b)).

##  3 The Compare command in GeoGebra Discovery

As a first example we show how the question (1) can be asked in a recent
version 2021Sep03 of GeoGebra Discovery444https://github.com/kovzol/geogebra-
discovery. The user constructs a triangle $ABC$. Implicitly the sides $a$, $b$
and $c$ will also be added by the program. Now the command
Compare($(a+b+c)^{2}$,$ab+bc+ca$) is to be typed and the program shows the
algebraic result in the bottom of the Graphics View as seen in Fig. 1.

Figure 1: A mechanical solution by GeoGebra Discovery for the problem Bottema
1.1

In fact, GeoGebra Discovery is programmed to realize that the expressions
$(a+b+c)^{2}$ and $ab+bc+ca$ both are homogeneous quadratic formulas.

We illustrate the same idea on an equality in a regular pentagon. In Fig. 2 a
regular 5-gon $ABCDE$ is shown with the midpoint $M$ of its circumcircle. We
introduce the quantities $f=AB$, $k=AC$, $l=AD$, $j=AE$, $R=AM$. After some
experiments one can find that there is an interesting relationship among the
product $f\cdot k\cdot l\cdot j$ and $R$. In this example, however, we use the
command Relation($f\cdot k\cdot l\cdot j$,$R$) to do a numerical check first
(Fig. 3(a)), and also to have a nicely looking symbolic result when clicking
on the button “More$\ldots$” (Fig. 3(b)).

Figure 2: Using the Relation command to obtain a relationship in a regular
pentagon

(a) Numerical check

(b) Symbolic check

Figure 3: Results on comparing $f\cdot k\cdot l\cdot j$ and $R$

The Relation command in GeoGebra is a well-known way to compare two objects,
from its first versions, but the comparison used to be performed only
numerically. If one wants to compare just two objects (that is, not
expressions, but, for example, segments), then a point-and-click comparison is
also possible. This is demonstrated in Fig. 4: after selecting the Relation
tool, two segments, e.g. $f$ and $k$, we learn that their ratio is
$\frac{1}{2}\cdot(\sqrt{5}-1)$ or $\frac{1}{2}\cdot(\sqrt{5}+1)$ (see Fig. 5).
Confusingly, it remains an unanswered question here whether the first or the
second one is the correct result.

Figure 4: The Relation tool in the toolbar in GeoGebra 6 Figure 5: Symbolic
relationship between the side and a diagonal of a regular pentagon

To explain this issue, we need to distinguish between two main GeoGebra
versions: 5 and 6. Version 5 has support to handle inequalities, but version 6
not at the moment. In version 6 the underlying computations in the algebro-
geometric setup cannot restrict the geometric figure to a single regular
pentagon. Instead of a single regular polygon, all possible star-regular
polygons are observed at the same time. In the case of a regular pentagon this
means that the computation for a regular pentagram will also be computed—two
results will be given concurrently. In general this means that the final
result always includes a solution according each star-regular variant.
(Fortunately, for a square or a regular hexagon there are no star-regular
variants, but for a regular heptagon there are two heptagrams—so three cases
are computed at the same time.)

In version 5 we are already capable of issuing a restriction that point $C$
takes place in the half-plane whose boundary line includes $B$ and is
perpendicular to the segment $AB$. In this case we can exclude the pentagram
variant (Fig. 6). With such considerations we can extend the equation based
approach to a much broader set of applications. Among others, it is possible
to attach a point on a segment or inside a triangle, or to explicitly define
the center of the incircle of a triangle.

Figure 6: Unambiguous result on the ratio between the side and a diagonal of a
regular pentagon (supported only in GeoGebra 5)

In the beginning of this section we mentioned the Compare command in GeoGebra
Discovery, and finally we arrived to a higher level command/tool: Relation. In
fact, Relation uses the internal routine of the Compare command, so, at the
end of the day, the same computations are performed, independently of the
launched command.

Finally we give some more description on the system layout of the above
mentioned two GeoGebra versions, 5 and 6. At the moment only GeoGebra 5 (a
native Java application) can give any result on inequalities, because all RQE
computations are outsourced to Mathematica or Tarski. Fig. 7 gives an overview
of the underlying technical hierarchy. GeoGebra internally uses the Giac
computer algebra system to perform a large set of symbolic computations
[GiacGG-RICAM2013].

Figure 7: Computation hierarchy for GeoGebra 5

By contrast, GeoGebra 6 (a web browser based variant of GeoGebra) is capable
only of obtaining equalities, because no connection has been made yet to any
systems that provide RQE computations. See Fig. 8 for the simple schema being
available for the moment.

Figure 8: Computation hierarchy for GeoGebra 6

##  4 Mathematical background

We point the reader to the first experiments on this topic to
[SturmWeispfenning96] (on proving inequalities) and [RecioVelez99] (on
obtaining equalities). By using a non-trivial but still simple example we
explain the main idea how our algorithms work.

The problem setting, as visible in Fig. 9, is to compare $f+g$ and $c$ in an
arbitrary triangle $ABC$—here $f$ and $g$ are two medians of the triangle, and
$c$ is the side of the triangle exactly as drawn in the figure. Now Fig. 10
shows how this problem is translated into an algebraic setup: each of the free
points $A$, $B$ and $C$ are expressed by two coordinates
($v_{1},v_{2},v_{3},v_{4},v_{5},v_{6}$), dependent points $D$ and $E$ are
given with another two times two coordinates ($v_{7},v_{8},v_{9},v_{10}$), and
another three variables describe the length of $c$, $g$ and $f$
($v_{11},v_{12},v_{13}$). All of these mappings are done automatically by
GeoGebra by following a mechanical algorithm, and also some equations are set
up to describe the relationships between these variables—in our case four
equations will be created.

Figure 9: Problem setting Figure 10: Algebraization of the problem setting

Four or five more technical variables are introduced:

  * •

$v_{14}$ is used to prescribe a condition on the non-degeneracy of $ABC$,
namely,

$v_{14}\cdot\left|\begin{matrix}v_{1}&v_{2}&1\\\ v_{3}&v_{4}&1\\\
v_{5}&v_{6}&1\end{matrix}\right|-1=0,$

by using the Rabinowitsch trick [Rabinowitsch1929]. Here the geometrical
meaning is an inequation, even if we use an equation for its algebraic
translation. This ensures that only algebraic equations are set up in the
first round and no algebraic inequality will be present. To use effective
elimination via Gröbner bases this step cannot be avoided.

  * •

$w_{1}$ and eventually $w_{2}$ are used to describe the two expressions with a
single variable, respectively.

  * •

$m$ plays the role of $m$, $M$ and $\mu$ as explained above.

  * •

$n$ is used to express that one of the two expressions differs from $0$. Here,
again, the Rabinowitsch trick is used.

After these settings GeoGebra uses elimination first if a suitable $\mu$
exists. To achieve this, the following Giac code is executed:

⬇

[assume(m>0),solve(eliminate(subst([2*v7-v5-v3,2*v8-v6-v4,2*v9-v5-v1,2*v10-v6-v2,-v12^2+v10^2+v9^2-2*v10*v4+v4^2-2*v9*v3+v3^2,-v11^2+v4^2+v3^2-2*v4*v2+v2^2-2*v3*v1+v1^2,-v13^2+v8^2+v7^2-2*v8*v2+v2^2-2*v7*v1+v1^2,-1-v14*v5*v4+v14*v6*v3+v14*v5*v2-v14*v3*v2-v14*v6*v1+v14*v4*v1,-w1+(v13+v12)^1,v11*m-(w1),(v11)*n-1],[v1=0,v2=0,v3=1,v4=0]),[v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,w1,n])[0],m)][1]

We note that the first free two points are positioned to $(0,0)$ and $(1,0)$,
respectively. This is allowed without loss of generality—with the assumption
that the algorithm always verifies if the expressions are homogeneous and of
the same degree. (See also [Davenport2017].)

If Giac gave a non-erroneous output here, then we managed to find the sought
$\mu$ (or, eventually a set of possible values for $\mu$). Otherwise there is
a second round needed—like in our case in this example. GeoGebra now
outsources the problem to a realgeom server, by composing an http request,
something like
http://192.168.10.20:8765/euclideansolver?lhs=w1&rhs=v11&polys=2*v7-v5-v3,2*v8-v6-v4,2*v9-v5-v1,2*v10-v6-v2,-v12This
triggers realgeom's internal computations that finally launch a request
towards Mathematica (or Tarski) like this:


#  $R$-mode Stability of GW190814’s Secondary Component as a Supermassive and
Superfast Pulsar

Xia Zhou  Xinjiang Astronomical Observatory, Chinese Academy of Sciences,
Urumqi, Xinjiang 830011, China  Ang Li  Department of Astronomy, Xiamen
University, Xiamen, Fujian 361005, China  Bao-An Li  Department of Physics and
Astronomy, Texas A&M University-Commerce, Commerce, TX 75429, USA

###### Abstract

The nature of GW190814’s secondary component $m_{2}$ of mass
$(2.50-2.67)\,\text{M}_{\odot}$ in the mass gap between the currently known
maximum mass of neutron stars and the minimum mass of black holes is currently
under hot debate. Among the many possibilities proposed in the literature, the
$m_{2}$ was suggested as a superfast pulsar while its r-mode stability against
the run-away gravitational radiation through the Chandrasekhar-Friedman-Schutz
mechanism is still unknown. Using those fulfilling all currently known
astrophysical and nuclear physics constraints among a sample of 33 unified
equation of states (EOSs) constructed previously by Fortin et al. (2016) using
the same nuclear interactions from the crust to the core consistently, we
compare the minimum frequency required for the $m_{2}$ to rotationally sustain
a mass higher than $2.50\,\text{M}_{\odot}$ with the critical frequency above
which the r-mode instability occurs. We use two extreme damping models
assuming the crust is either perfectly rigid or elastic. Using the stability
of 19 observed low-mass x-ray binaries as an indication that the rigid crust
damping of the r-mode dominates within the models studied, we find that the
$m_{2}$ is r-mode stable while rotating with a frequency higher than 870.2 Hz
(0.744 times its Kepler frequency of 1169.6 Hz) as long as its temperate is
lower than about $3.9\times 10^{7}K$, further supporting the proposal that
GW190814’s secondary component is a supermassive and superfast pulsar.

Pulsars (1306); Neutron stars (1108)

##  1 Introduction

The recent discovery by the LIGO/Virgo Collaborations that the compact binary
merger GW190814 has a secondary component $m_{2}$ with mass
$(2.50-2.67)\,\text{M}_{\odot}$ (Abbott et al., 2020) at $90\%$ credible level
has generated much interest in the astrophysics community. As the $m_{2}$ is
in the mass gap between the currently known maximum mass of neutron stars and
the minimum mass of black holes, its nature has significant ramifications on
many interesting issues in astrophysics and cosmology. Many possibilities
regarding its nature and formation path have been proposed very recently in
the literature; see, e.g., Roupas et al. (2020); Biswas et al. (2020); Bombaci
et al. (2020) among the latest reports and references therein.

Since the maximum mass $M_{\rm TOV}$ of non-rotating neutron stars is
predicted to be about $2.4\,\text{M}_{\odot}$ on the causality surface with
the equation of state (EOSs) satisfying all known constraints (e.g., Zhang &
Li, 2019; Li et al., 2020) and it is well known that the centrifugal force in
neutron stars rotating at Kepler frequencies can increase their maximum masses
by about 20% with respect to the $M_{\rm TOV}$, the mass of $m_{2}$ is well
within reach with the rotational support. Indeed, this possibility is among
the first one considered. However, the conclusions have been diverse mainly
because of the different nuclear EOSs used (see e.g., Abbott et al., 2020;
Biswas et al., 2020; Dexheimer et al., 2020; Godzieba et al., 2020; Huang et
al., 2020; Tews et al., 2020; Most et al., 2020; Zhang & Li, 2020; Zhang &
Mann, 2020). Among the conclusions supporting the $m_{2}$ as a superfast
pulsar, Most et al. (2020) found that $m_{2}$ has a rotation frequency of 1210
Hz assuming it has a typical neutron star radius of 12.5 km and a $M_{\rm
TOV}$=$2.08\,\text{M}_{\odot}$. In another study by Zhang & Li (2020), the
$m_{2}$ was found to have a rotation frequency of 971 Hz using a EOS leading
to an equatorial radius of 11.9 km for $m_{2}$ and a $M_{\rm
TOV}$=$2.39\,\text{M}_{\odot}$. More recently, Biswas et al. (2020) derived a
minimum frequency of $1143^{+194}_{-155}$ Hz and an equatorial radius
$R_{e}=15.7^{+1.0}_{-1.7}$ km for $m_{2}$ at 90% confidence level assuming a
$M_{\rm TOV}$=$2.14\,\text{M}_{\odot}$ which is the mass of MSR J0740+6620
(Cromartie et al., 2020). The minimum frequencies found for $m_{2}$ in these
studies are significantly higher than the fastest known pulsar PSR
J1748-2446ad with a frequency of 716 Hz (Hessels et al., 2006), thus making
GW190814’s secondary the most massive and fastest known pulsar if confirmed.

While the LIGO/Virgo Collaborations tightly constrained the primary spin of
GW190814, the spin of its secondary $m_{2}$ remains unconstrained (Abbott et
al., 2020). Thus, the possibility for the $m_{2}$ as a superfast pulsar is
neither confirmed nor ruled out observationally. However, an outstanding
question critical for the $m_{2}$ to be a superfast pulsar remains to be
answered. Namely, it is unknown if the $m_{2}$ with the required minimum
frequencies are stable against the well-known r-mode instability leading to
the exponential growth of gravitational radiation. As it was pointed out
already in Zhang & Li (2020); Biswas et al. (2020), for the $m_{2}$ to be a
supermassive and superfast pulsar, its rotational instabilities would have to
be suppressed, otherwise, it is more likely to be a black hole instead.

In this work, we examine the r-mode instability of $m_{2}$ using well
established EOSs satisfying all currently known constraints within the r-mode
formalisms established by Lindblom et al. (1998, 2000). While a multitude of
damping mechanisms exist, an instability can develop only if its growth is
faster than its strongest damping mechanism (Biswas et al., 2020). It is known
that a rigid crust provides the strongest r-mode damping, and it can well
explain the stability of all observed low-mass x-ray binaries (LMXBs) (e.g.,
Ho et al., 2011; Haskell et al., 2012). Using this damping mechanism, we found
that the $m_{2}$ can be r-mode stable as long as its temperature is
sufficiently low, e.g., lower than about $3.9\times 10^{7}K$ for the $m_{2}$
rotating at 1169.6 Hz (0.744 times its Kepler frequency). Since this
temperature is about an order of magnitude higher than that of some known old
neutron stars (see, e.g., Wijngaarden et al., 2019), the r-mode instability
should not be a concern for the $m_{2}$ as a supermassive and superfast pulsar
within the theoretical framework and models considered.

The r-mode instability (Andersson, 1998; Friedman & Morsink, 1998; Andersson &
Kokkotas, 1998; Ho & Lai, 2000) can trigger the exponential growth of
gravitational wave (GW) emission in rapidly rotating neutron stars through the
Chandrasekhar-Friedman-Schutz mechanism (Chandrasekhar, 1970; Friedman &
Schutz, 1978) if the GW growth rate is faster than its damping rate. Over the
past two decades, significant efforts have been devoted to better
understanding the r-mode instability and its damping mechanisms with many
interesting findings (see, e.g., Madsen, 2000; Sá, 2004; Sá & Tomé, 2005;
Andersson et al., 2010; Haskell & Andersson, 2010; Yang et al., 2011; Wen et
al., 2012; Vidaña, 2012; Haskell et al., 2012; Alford & Schwenzer, 2014a, b;
Idrisy et al., 2015; Moustakidis, 2015; Dai et al., 2016; Papazoglou &
Moustakidis, 2016; Ofengeim et al., 2019; Wang et al., 2019; Krüger &
Kokkotas, 2020). In particular, it was found that the r-mode instability
window (a region in the spin frequency versus temperature plane) in which the
r-mode is unstable depends sensitively on the EOS of neutron-rich matter,
especially its symmetry energy term which encodes the information about the
energy cost to make nuclear matter more neutron-rich. It depends even more
sensitively on the poorly known properties of neutron star crust, especially
whether it is rigid or elastic (see more discussions in recent reviews, e.g.,
Andersson & Kokkotas, 2001; Haskell, 2015; Andersson, 2017; Kokkotas &
Schwenzer, 2016). To realize the goal of this work, it is thus critical to use
unified EOSs for neutron stars not only satisfying all currently known
constraints from both nuclear physics and astrophysics but also describing the
crust, crust-core transition, and the core consistently based on the same
nuclear interactions within the same nuclear many-body theories. Moreover, it
is important to limit the uncertain r-mode damping mechanisms as much as
possible and to identify the strongest one using existing astrophysical
observations.

For the purposes outlined above, we describe in the next section how we
selected the unified EOSs satisfying the conditions mentioned above. In
Section 3, we compare 19 observed LMXBs with known frequencies and
temperatures with the r-mode instabilities boundaries calculated with the
selected unified EOSs assuming the crust is either perfectly rigid or elastic.
As expected, the stability of the 19 LMXBs prefers the perfectly rigid crust.
Assuming the same damping mechanism is at work in other neutron stars, we
infer in Section 4 the maximum temperature of $m_{2}$ by comparing the minimum
frequency it must have to sustain rotationally a mass of
$2.50\,\text{M}_{\odot}$ with the maximum frequency it can have before the
r-mode instability occurs. We then summarize our work in Sec. 5.

##  2 Unified EOSs for neutron stars

Various EOSs for neutron stars have been used in the literature. Often, they
are constructed from combining the EOSs for the core and crust calculated
using different models and/or interactions. For most purposes, these EOSs are
perfectly fine. As we discussed above, for the purpose of this work, it is
advantageous to use the EOSs with the crust, crust-core transition, and core
all calculated using the same nuclear interaction. We thus adopt the 33
unified EOSs derived by Fortin et al. (2016) within the Skyrme Hartree-Fock
and the Relativistic mean-field models for the core and the Thomas–Fermi model
for the crust using the same interactions. Because new progress has been made
in constraining the EOS since these 33 EOSs were derived, we shall first
filter them through the following three tests: (1) satisfying new constraints
on the density dependence of nuclear symmetry energy from analyzing both
nuclear experiments and astrophysical observations especially the radii and
tidal deformability of neutron stars (Newton & Crocombe, 2020; Xie & Li,
2020), (2) being stiff enough to support the presently known maximum mass of
neutron stars, i.e., the $M=2.14^{+0.10}_{-0.09}\,\text{M}_{\odot}$ (68%
confidence level) of MSP J0740+6620 (Cromartie et al., 2020), (3) be
consistent with the simultaneous mass and radius measurements by NICER for PSR
J0030+0451 (Miller et al., 2019; Raaijmakers et al., 2019; Riley et al.,
2019).

Figure 1: $K_{\rm sym}$ parameter as a function of $L$ for the 33 unified
neutron star EOSs from Fortin et al. (2016), together with two recent
constraints, to the $68\%$ confidence level, on the slope $L$ and curvature
$K_{\rm sym}$ of the symmetry energy at saturation (Newton & Crocombe, 2020;
Xie & Li, 2020). Note that the 95% confidence boundaries shown in Newton &
Crocombe (2020) are rescaled to the 68% credible ranges here for a comparison.
Only 7 EOSs: DD2, DDME2, KDE0v1, SKb, SkI6, SLy2, and SLy9 fulfill both
constraints.

The nuclear symmetry energy is well-known to have significant effects on the
structure of neutron stars (Lattimer & Prakash, 2000; Steiner et al., 2005; Li
et al., 2019), including the crust-core transition, the crust thickness, as
well as the composition of the star. For a comprehensive review on nuclear
symmetry energy and its astrophysical effects, we refer the reader to Li et
al. (2014). It is defined as

$\displaystyle E_{\rm
sym}(n)=\frac{1}{2}\frac{\partial^{2}E(n,\beta)}{\partial\beta^{2}}|_{\beta=0}\
,$ (1)

where $E(n,\beta)$ is the binding energy of neutron-rich matter as a function
of density $n$ and isospin asymmetry $\beta=(n_{n}-n_{p})/n$. The $E_{\rm
sym}(n)$ can be characterized by using its slope $L$ and curvature $K_{\rm
sym}$ at the nuclear saturation density $n_{0}$ defined as

$\displaystyle L$ $\displaystyle=$ $\displaystyle 3n_{0}[\partial E_{\rm
sym}(n)/\partial n]|_{n=n_{0}}\ ,$ (2) $\displaystyle K_{\rm sym}$
$\displaystyle=$ $\displaystyle 9n_{0}^{2}[\partial^{2}E_{\rm sym}(n)/\partial
n^{2}]|_{n=n_{0}}\ .$ (3)

Figure 2: Mass-radius relations of neutron stars using the 6 EOSs from Fig. 1
satisfying the latest constraints on $L$ and $K_{\rm sym}$. The astrophysical
constraining bands on the maximum mass (Cromartie et al., 2020) and the mass-
radius correlations for PSR J0030+0451 from NICER (Miller et al., 2019;
Raaijmakers et al., 2019; Riley et al., 2019) are shown for comparisons.

It was shown earlier that both the $L$ and especially $K_{\rm sym}$ as well as
their correlation affect the crust-core transition density and pressure
significantly (e.g., Xu et al., 2009; Vidaña et al., 2009; Zhang et al., 2018;
Providência et al., 2019; Li & Magno, 2020). They also affect the EOS and
especially the composition of the core (Lattimer et al., 1991). Fortunately,
some significant progresses have been made recently in experimentally
constraining the slope parameter $L$ (Li & Han, 2013; Baldo & Burgio, 2016;
Oertel et al., 2017), while the curvature $K_{\rm sym}$ is still much less
constrained. In particular, a recent Bayesian analysis of the radii and tidal
deformations of canonical neutron stars (Xie & Li, 2020) inferred the most
probable values of $L=66^{+12}_{-20}\,\text{MeV}$ and $K_{\rm
sym}=-120^{+80}_{-100}\,\text{MeV}$ at $68\%$ confidence level. In another
Bayesian analysis by Newton & Crocombe (2020) using combined data of neutron
skin in 48Ca, 208Pb and tin isotopes as well as the best theoretical
information about the EOS of pure neutron matter from ab initio microscopic
nuclear many-body theories, the most probable values of $L$ and $K_{\rm sym}$
were found to be $L=40^{+34}_{-26}\,\text{MeV}$ and $K_{\rm
sym}=-209^{+270}_{-182}\,\text{MeV}$, respectively, at $95\%$ confidence
level.

Shown in Fig. 1 are the two constraints on $L$ and $K_{\rm sym}$ in comparison
with predictions of the 33 unified EOSs from Fortin et al. (2016). It is seen
that only 7 of them fall into the overlapping area of the two latest
constraints (Newton & Crocombe, 2020; Xie & Li, 2020). We further test them
against the latest astrophysical observations in Fig. 2. It is seen that the
KDE0v1 EOS is further excluded by the mass measurement of MSP J0740+6620, to
the $68.3\%$ credibility interval, while the remaining 6 EOSs can support a
maximum mass of about $2.14\,\text{M}_{\odot}$ (Cromartie et al., 2020) and
satisfy the mass-radius constraints from NICER (Miller et al., 2019; Riley et
al., 2019). Thus, we will only use these 6 unified EOSs in our further studies
in the following.

Figure 3: Locations of LMXBs in the frequency-temperature plane and lower
bounds $\nu_{\rm crit}$ of the r-mode instability windows for
$1.4\,\text{M}_{\odot}$ (left) and $2.0\,\text{M}_{\odot}$ (right) neutron
stars with 6 selected unified EOSs. The upper (lower) group of lines are
calculated with rigid (elastic) crusts. The observational data and internal
redshifted temperatures ($T^{\infty}=T\sqrt{1-2GM/Rc^{2}}$) on neutron stars
in LMXBs are taken from Gusakov et al. (2014); Chugunov et al. (2017), with
the error bars due to the uncertainty in the composition of the envelope (see
detail in Gusakov et al., 2014; Chugunov et al., 2017).

##  3 R-mode instability in neutron stars

All rotating neutron stars are generically r-mode unstable due to the emission
of gravitational waves. However, there are several possible damping mechanisms
preventing the r-mode oscillation from growing exponentially. Here we adopt
two extreme damping models assuming the crust is either rigid or elastic. For
easy of discussions, in the following, we first recall briefly some of the
basic definitions and formulas most relevant for our study about the r-mode
instability of $m_{2}$. We refer the reader to the original papers for more
details.

###  3.1 Time scales of r-mode growth and damping

Whether the r-mode oscillation can grow exponentially or not depends on the
competition between the growth rate of gravitational wave emission and the
r-mode damping rate (mainly due to viscosity). Thus, the r-mode instability
boundary is determined by setting the damping $\tau_{\rm diss}$ and driving
$\tau_{\rm gw}$ timescales equal to each other (Andersson & Kokkotas, 2001),
i.e., by solving for the zeros of

$-\frac{1}{2\tilde{E}}\left(\frac{d\tilde{E}}{dt}\right)=\frac{1}{\tau_{\rm
gw}(\nu)}+\sum\frac{1}{\tau_{\rm diss}(\nu,T)}=0\ ,$ (4)

where $\tilde{E}$ is the total energy of an r-mode oscillation. It is known
that the shear viscosity dominates the energy dissipation for temperatures
around $10^{8}~{}\rm{K}$, while at higher temperatures (of the order of
$10^{10}$ K) the bulk viscosity becomes important (Andersson & Kokkotas,
2001). Since all the observed LMXBs have temperatures around $10^{8}~{}\rm{K}$
and the old neutron stars are cooler, in the present study we consider only
the shear viscosity resulting from both neutron-neutron (nn) and electron-
electron (ee) scatterings. We adopt the following parameterizations for the
two viscosities as functions of temperature $T$ (Flowers & Itoh, 1979; Cutler
& Lindblom, 1987; Shternin & Yakovlev, 2008),

$\displaystyle\eta_{nn}(r)$ $\displaystyle=$ $\displaystyle
347[\rho(r)]^{9/4}T^{-2}~{}~{}\rm{(g~{}cm^{-1}s^{-1})}\ ,$ (5)
$\displaystyle\eta_{ee}(r)$ $\displaystyle=$ $\displaystyle 6\times
10^{6}[\rho(r)]^{2}T^{-2}~{}~{}\rm{(g~{}cm^{-1}s^{-1})}\ ,$ (6)

where $\rho(r)$ is the mass density profile of the star.

The crust is expected to play an important role in determining the r-mode
stability of neutron stars. Here we consider two extreme cases. Assume the
crust is perfectly rigid (i.e., the “slippage” factor $S=1$), so the fluid
motion associated with the r-mode would be significantly damped at the crust-
core boundary, the $\tau_{\rm gw}$ is given by (Lindblom et al., 2000),

$\displaystyle\frac{1}{\tau_{\rm gw}^{\rm c}}=-\frac{32\pi
G\Omega^{2l+2}}{c^{2l+3}}\frac{(l-1)^{2l}}{[(2l+1)!!]^{2}}\left(\frac{l+2}{l+1}\right)^{2l+2}$
$\displaystyle\times\int_{0}^{R_{\rm c}}\rho(r)r^{2l+2}dr\ ,$ (7)

where $G$ is the gravitational constant, $\Omega=2\pi\nu$ is the angular
velocity of the star, $c$ is the speed of light and $R_{c}$ is the stellar
radius at the crust-core transition mass density $\rho_{c}$. While the
$\tau_{\rm diss}$ can be evaluated from (Lindblom et al., 2000),

$\tau_{\eta}^{\rm
c}=\frac{1}{2\Omega}\frac{2^{l+3/2}(l+1)!}{l(2l+1)!!\mathcal{C}_{l}}\sqrt{\frac{2\Omega
R_{c}^{2}\rho_{c}}{\eta_{c}}}\int_{0}^{R_{c}}\frac{\rho}{\rho_{c}}\left(\frac{r}{R_{c}}\right)^{2l+2}\frac{dr}{R_{c}},$
(8)

where $\eta_{c}$ is the viscosity at the crust-core transition density.
Following Lindblom et al. (1998), we only consider the lowest order
contribution of the multipole moment ($l=2$), with $\mathcal{C}_{2}=0.80411$
(Rieutord, 2001).

In the other extreme case neglecting the crustal damping, or the crust is
regarded as extremely elastic [i.e., the “slippage” factor $S\rightarrow 0$
(see e.g., Glampedakis & Andersson, 2006)], the crust will also participate in
the oscillation. In this case, the $\tau_{\rm gw}$ is given by (Lindblom et
al., 1998),

$\displaystyle\frac{1}{\tau_{\rm gw}}=-\frac{32\pi
G\Omega^{2l+2}}{c^{2l+3}}\frac{(l-1)^{2l}}{[(2l+1)!!]^{2}}\left(\frac{l+2}{l+1}\right)^{2l+2}$
$\displaystyle\times\int_{0}^{R}\rho(r)r^{2l+2}dr\ .$ (9)

While the $\tau_{\eta}$ is written as (e.g., Vidaña, 2012),

$\displaystyle\frac{1}{\tau_{\eta}}=(l-1)(2l+1)\left(\int_{0}^{R}\rho(r)r^{2l+2}dr\right)^{-1}$
$\displaystyle\times\int_{0}^{R}\eta(r)r^{2l}dr\ .$ (10)

The actual neutron star crust should have some elasticity, therefore by
calculating the r-mode boundary in the two extreme cases: Eqs. (7)-(8) (where
a rigid crust is considered) and Eqs. (9)-(10) (where no crustal damping is
included), we may estimate the uncertainty due to the crustal modelling.
Moreover, confronting the theoretical results with the observations of neutron
stars’ frequency and temperatures in LMXBs, we may get some hints about which
damping mechanism dominates.

###  3.2 Lower boundaries of r-mode instability windows

From Eq. (4), for each star at a given temperature, one can find its critical
frequency $\nu_{\rm crit}$ above which the star becomes unstable against the
running-off of gravitational radiation, by setting $1/{\tau_{\rm gw}(\nu_{\rm
crit})}+1/\tau_{\eta}(\nu_{\rm crit},T)=0$.The $\nu_{\rm crit}$ serves as an
upper limit of stable pulsars’ frequencies for a given temperature. The region
above this boundary in the frequency-temperature plane is the so-called r-mode
instability window. In the following, we examine the r-mode instability
boundaries with the 6 selected EOSs with respect to the locations of the 19
LMXBs. Since the masses of the neutron stars in these LMXBs are not measured
accurately, the calculations are done for both canonical
($M=1.4\,\text{M}_{\odot}$) and massive ($M=2.0\,\text{M}_{\odot}$) neutron
stars.

Our results are shown in Fig. 3 for both cases assuming a perfectly rigid
crust (upper groups of lines) or an elastic one (lower group of lines).
Several interesting observations can be made:

  * •

The r-mode instability boundary depends most sensitively on the crust’s
elasticity while effects of the EOS and neutron star mass are appreciable.
Most importantly, for the present study, the rigid crust provides the
strongest damping.

  * •

The damping mechanism with a perfectly rigid crust can accommodate all neutron
stars in the 19 LMXBs, while the one with an elastic crust can only
accommodate a few slowly-rotating neutron stars, such as NGC 6440 X-2, XTE
J1807-294, and XTE J0929-314.

  * •

The r-mode instability boundaries are sensitive to the EOS especially the $L$
parameter (see e.g., Wen et al., 2012), but the effects shown here are much
less than previously reported in the literature mainly because the 6 unified
EOSs used here are much more stringently constrained compared to those used
before. Nevertheless, one can still clearly see that the $L$ parameter (shown
in Fig. 1) can not be very large, e.g., smaller than $\sim 60\,\text{MeV}$ if
LMXB neutron stars are massive (e.g., $M=2.0M\,\text{M}_{\odot}$), to ensure
the rapidly-rotating neutron stars in, e.g., 4U 1608-522 and SAX J1705.8-2900,
stay r-mode stable.

  * •

Comparing the results for $1.4\,\text{M}_{\odot}$ and $2.0\,\text{M}_{\odot}$
neutron stars, we see that the r-mode instability window is generally broader
for the latter in agreement with earlier findings (e.g., Vidaña, 2012),
especially in the case with the rigid crust. Moreover, the uncertainty band
due to the EOS is relatively wider for neutron stars of mass
$2.0\,\text{M}_{\odot}$, due to the more significant diversity of the EOS at
higher densities reached.

Figure 4: Temperature dependence of the reduced critical frequency $\nu_{\rm
crit}/\nu_{\rm K}$ for $M=1.4\,\text{M}_{\odot}$ and $M=2.0\,\text{M}_{\odot}$
neutron stars with rigid (upper groups of lines) or elastic crusts (lower
groups of lines).

Since the Kepler frequency $\nu_{\rm K}$ is the absolute upper limit of spin
frequencies of all stars with a given mass, it is interesting and informative
to examine the reduced critical frequency $\nu_{\rm crit}/\nu_{\rm K}$ as
functions of the temperature $T$ for both $1.4\,\text{M}_{\odot}$ and
$2.0\,\text{M}_{\odot}$ neutron stars. Our results for $\nu_{\rm
crit}/\nu_{\rm K}$ are shown in Fig. 4. Compared to the un-scaled results
shown in Fig. 3, it is seen that the r-mode instability boundaries of the
reduced frequency are more sensitive to the neutron star mass with both rigid
and elastic crusts. This is because the $\nu_{\rm K}$ is proportional to the
square root of the average density of neutron stars (Andersson & Kokkotas,
2001). The average density depends on the mass; thus, the reduced frequency
becomes more sensitive to the neutron star mass. The lighter neutron stars
have a lower average density and a lower $\nu_{\rm K}$, making them have
higher reduced frequencies than the massive ones.

Since the Kepler frequency is EOS dependent (see, e.g., Li et al., 2016,
2017), the reduced frequencies thus obtain a more complicated EOS dependence.
With the rigid crust, as shown in Eqs. (7)-(8), the location $R_{c}$ of the
crust-core transition plays an important role. In this case, the strong EOS
dependence of $R_{c}$ makes the reduced frequency more strongly EOS dependent
than the case with an elastic crust. While in the latter case, both the
$\nu_{\rm crit}$ and $\nu_{\rm K}$ are calculated by integrating from the
center to the surface. The reduced frequency obtained appears rather EOS
insensitive due to some canceling effects in taking the ratio. Finally, the
difference in the reduced frequency between the two cases at the same
temperature are very large. More quantitatively, for a typical temperature
$T=10^{8}~{}\rm K$ relevant for the LMXBs, $\nu_{\rm crit}$ is lifted from
$\sim 0.16~{}\nu_{\rm K}$ ($\sim 0.13~{}\nu_{\rm K}$) to $\sim 0.79~{}\nu_{\rm
K}$ ($\sim 0.48~{}\nu_{\rm K}$) for $M=1.4\,\text{M}_{\odot}$
($M=2.0\,\text{M}_{\odot}$) neutron stars when assuming a rigid crust with
respect to the elastic one.

In short, the main lessons we learned from this section that are most
important for our following study about the r-mode stability of $m_{2}$ are
(1) the rigid crust provides the strongest damping, (2) the stability of the
19 neutron stars in LMXBs favors the rigid crust damping, (3) the reduced
critical frequency $\nu_{\rm crit}/\nu_{\rm K}$ with the rigid crust is
sensitive to both the neutron star mass and the EOS used.

##  4 R-mode stability of GW190814’s secondary component of mass
$2.50\,\text{M}_{\odot}$

We now turn to examine the r-mode stability of GW190814’s secondary component
$m_{2}$ as a supermassive and superfast pulsar. For this purpose, we first
find the minimum frequency $\nu_{\rm min}$ to support GW190814’s secondary
component with a minimum mass of $2.50\,\text{M}_{\odot}$ (Abbott et al.,
2020) using the well-tested ${\rm RNS}$ code (Stergioulas et al., 1995) for a
given EOS. Among the 6 EOSs used above, we found that the SLy2 EOS can not
rotationally support a $2.50\,\text{M}_{\odot}$ star even at the Kepler
frequency. Thus, only the remaining 5 EOSs are used. As discussed in the
previous section, only the damping with the rigid crust is necessary for this
discussion as it provides the strongest/quickest damping.

Table 1: Neutron star properties based on the 5 unified EOSs models (Fortin et
al., 2016). $M_{\rm TOV}$ is the maximum mass of non-rotating neutron stars
supported by the given EOS. $\nu_{\rm min}$ is the minimum frequency (in unit
of the Kepler frequency $\nu_{\rm K}$) to support GW190814’s secondary
component $m_{2}$ with a minimum mass of $2.50\,\text{M}_{\odot}$. $T_{\rm
max}$ is the maximum temperature for the $m_{2}$ to remain r-mode stable.

Model | $M_{\rm TOV}$ | $\nu_{\rm K}$ | $\nu_{\rm min}$ | $T_{\rm max}$
---|---|---|---|---
| $(\,\text{M}_{\odot})$ | (Hz) | ($\nu_{\rm K}$) | (K)
DD2 | $2.42$ | $1196.6$ | $0.755$ | $3.0\times 10^{7}$
DDME2 | $2.48$ | $1169.6$ | $0.744$ | $3.9\times 10^{7}$
SKb | $2.20$ | $1446.5$ | $0.809$ | $1.3\times 10^{7}$
SkI6 | $2.20$ | $1433.4$ | $0.827$ | $1.1\times 10^{7}$
SLy9 | $2.16$ | $1515.4$ | $0.855$ | $7.5\times 10^{6}$

Listed in Table 1 are the maximum mass of non-rotating neutron stars $M_{\rm
TOV}$, Kepler frequency $\nu_{\rm K}$ and the minimum frequency $\nu_{\rm
min}$ to support a neutron star of mass $2.50\,\text{M}_{\odot}$ for the 5
EOSs used. The $\nu_{\rm min}$ is then plotted as a horizontal line in Fig. 5
where the r-mode stability boundary for each EOS is shown in the frequency-
temperature plane. Their cross point is marked with a diamond indicating the
maximum temperatures $T_{\rm max}$ below which neutron stars remain r-mode
stable. The resulting 5 maximum temperatures are listed in the last column in
Table 1. While the allowed $\nu-T$ parameter spaces are reported with the
shaded regions in Fig. 5.

Figure 5: Temperature dependence of the reduced critical frequency $\nu_{\rm
crit}/\nu_{\rm K}$ for neutron stars with rigid crusts. For each EOS, the
$\nu-T$ space allowing the secondary component of GW190814 as a r-mode stable
pulsar is shown as the shaded region. The maximum temperatures $T_{\rm max}$
below which neutron stars remain r-mode stable are marked with diamonds on the
r-mode instability thresholds for the 5 EOSs used.

As expected, the $\nu_{\rm min}$ depends sensitively on the $M_{\rm TOV}$
determined by the stiffness of the EOS. For the 5 EOSs used, the $\nu_{\rm
min}$ and $M_{\rm TOV}$ range between 0.744 to $0.855\,\nu_{\rm K}$ and 2.16
to $2.42\,\text{M}_{\odot}$, respectively. The corresponding maximum
temperature $T_{\rm max}$ is between $7.5\times 10^{6}$ to $3.9\times 10^{7}$
K. This temperature limit is about 10 times cooler than the LMXBs but also
about 10 times hotter than some of the known old neutron stars (Wijngaarden et
al., 2019). This information may be useful for further understanding the
dynamical path of the GW190814 event and the nature of its secondary
component.

The required minimum frequencies $\nu_{\rm min}$ between 870 Hz and 1295 Hz
found here with the 5 selected EOSs satisfying all currently known constraints
from both nuclear physics and astrophysics cover the minimum frequencies found
necessary for the $m_{2}$ to rotationally sustain a minimum mass of
$2.50\,\text{M}_{\odot}$ in Most et al. (2020); Zhang & Li (2020); Biswas et
al. (2020). As discussed above, as long as the $m_{2}$ has a temperature below
about $3.9\times 10^{7}$ K, the $m_{2}$ spinning with the high frequencies
found above can remain r-mode stable. While it is not the purpose of this work
to study how the high mass and spin can be obtained, it is very interesting to
note that it was already pointed out by Safarzadeh & Loeb (2020) that they
could be supplied through the circumbinary accretion disk if the $m_{2}$ was
born as a neutron star where a significant amount of the supernova ejecta mass
from its formation remained bound to the binary due to the presence of the
massive black hole companion in the GW190814 event.

##  5 Conclusions

In conclusion, the GW190814’s secondary component $m_{2}$ can be an r-mode
stable supermassive and superfast pulsar rotating with a frequency higher than
870.2 Hz (0.744 times its Kepler frequency of 1169.6 Hz) as long as its
temperate is lower than about $3.9\times 10^{7}K$ which is still about 10
times hotter than some of the known old neutron stars. Thus, the r-mode
instability should not be a concern for the $m_{2}$. The minimum frequency and
limiting temperature found here may be useful for further understanding the
formation mechanism of $m_{2}$ and the merger dynamics of GW190814.

Our study is carried out within the minimum model of neutron stars consisting
of only nucleons and leptons. Moreover, the fast rotation is probably the
simplest mechanism to provide the additional support besides the nuclear
pressure against the strong gravity of massive neutron stars. Furthermore, we
selected the EOSs satisfying all currently known astrophysical and nuclear
physics constraints from a group of 33 unified EOSs constructed from the crust
to the core consistently using the same nuclear interactions. Thus, our
approach is probably among the most conservative ones used in the literature
in investigating the nature of GW190814’s secondary component.

A major caveat of our work is the assumption that the rigid crust of neutron
stars provides the strongest r-mode damping among a multitude of possible
damping mechanisms. While both our current work and previous theoretical
calculations by others as well as the stability of the 19 neutron stars in
LMXBs favor the dominating rigid crust damping, there may be other even
stronger damping mechanisms and/or model ingredients affecting significantly
the location of the r-mode instability windows. Nevertheless, our current
study can successfully describe all relevant observations of neutron stars. We
are thus confident that our present conclusions are physically sound and
useful for the community to finally solve the mystery regarding the nature of
GW190814’s secondary component.

We would like to thank Morgane Fortin for valuable discussions. This work is
supported in part by National SKA Program of China (No. 2020SKA0120300), the
National Natural Science Foundation of China (Grant Nos. 11873040, 12033001),
the Youth Innovation Fund of Xiamen (No. 3502Z20206061), the CAS “Light of
West China” Program (No. 2019-XBQNXZ-B-016), the Tianshan Youth Program (No.
2018Q039), the U.S. Department of Energy, Office of Science, under Award
Number DE-SC0013702, the CUSTIPEN (China-U.S. Theory Institute for Physics
with Exotic Nuclei) under the US Department of Energy Grant No. DE-SC0009971.

## References

  * Abbott et al. (2020) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020, ApJ, 896, L44. doi:10.3847/2041-8213/ab960f
  * Alford & Schwenzer (2014a) Alford, M. G. & Schwenzer, K. 2014, Nucl. Phys. A, 931, 740. doi:10.1016/j.nuclphysa.2014.08.077
  * Alford & Schwenzer (2014b) Alford, M. G. & Schwenzer, K. 2014, Phys. Rev. Lett., 113, 251102. doi:10.1103/PhysRevLett.113.251102
  * Andersson & Kokkotas (1998) Andersson, N. & Kokkotas, K. D. 1998, MNRAS, 299, 1059. doi:10.1046/j.1365-8711.1998.01840.x
  * Andersson (1998) Andersson, N. 1998, ApJ, 502, 708. doi:10.1086/305919
  * Andersson & Kokkotas (2001) Andersson, N. & Kokkotas, K. D. 2001, International Journal of Modern Physics D, 10, 381. doi:10.1142/S0218271801001062
  * Andersson et al. (2010) Andersson, N., Haskell, B., & Comer, G. L. 2010, Phys. Rev. D, 82, 023007. doi:10.1103/PhysRevD.82.023007
  * Andersson (2017) Andersson, N. 2017, Journal of Astrophysics and Astronomy, 38, 58. doi:10.1007/s12036-017-9463-8
  * Baldo & Burgio (2016) Baldo, M. & Burgio, G. F. 2016, Progress in Particle and Nuclear Physics, 91, 203. doi:10.1016/j.ppnp.2016.06.006
  * Biswas et al. (2020) Biswas, B., Nandi, R., Char, P., et al. 2020, arXiv:2010.02090
  * Bombaci et al. (2020) Bombaci, I., Drago, A., Logoteta, D., et al. 2020, arXiv:2010.01509
  * Chandrasekhar (1970) Chandrasekhar, S. 1970, Phys. Rev. Lett., 24, 611. doi:10.1103/PhysRevLett.24.611
  * Chugunov et al. (2017) Chugunov, A. I., Gusakov, M. E., & Kantor, E. M. 2017, MNRAS, 468, 291. doi:10.1093/mnras/stx391
  * Cromartie et al. (2020) Cromartie, H. T., Fonseca, E., Ransom, S. M., et al. 2020, Nature Astronomy, 4, 72. doi:10.1038/s41550-019-0880-2
  * Cutler & Lindblom (1987) Cutler, C. & Lindblom, L. 1987, ApJ, 314, 234. doi:10.1086/165052
  * Dai et al. (2016) Dai, Z. G., Wang, S. Q., Wang, J. S., et al. 2016, ApJ, 817, 132. doi:10.3847/0004-637X/817/2/132
  * Dexheimer et al. (2020) Dexheimer, V., Gomes, R. O., Klähn, T., et al. 2020, arXiv:2007.08493
  * Flowers & Itoh (1979) Flowers, E. & Itoh, N. 1979, ApJ, 230, 847. doi:10.1086/157145
  * Fortin et al. (2016) Fortin, M., Providência, C., Raduta, A. R., et al. 2016, Phys. Rev. C, 94, 035804. doi:10.1103/PhysRevC.94.035804
  * Friedman & Schutz (1978) Friedman, J. L. & Schutz, B. F. 1978, ApJ, 222, 281. doi:10.1086/156143
  * Friedman & Morsink (1998) Friedman, J. L. & Morsink, S. M. 1998, ApJ, 502, 714. doi:10.1086/305920
  * Glampedakis & Andersson (2006) Glampedakis, K. & Andersson, N. 2006, Phys. Rev. D, 74, 044040. doi:10.1103/PhysRevD.74.044040
  * Godzieba et al. (2020) Godzieba, D. A., Radice, D., & Bernuzzi, S. 2020, arXiv:2007.10999
  * Gusakov et al. (2014) Gusakov, M. E., Chugunov, A. I., & Kantor, E. M. 2014, Phys. Rev. D, 90, 063001. doi:10.1103/PhysRevD.90.063001
  * Haskell & Andersson (2010) Haskell, B. & Andersson, N. 2010, MNRAS, 408, 1897. doi:10.1111/j.1365-2966.2010.17255.x
  * Haskell et al. (2012) Haskell, B., Degenaar, N., & Ho, W. C. G. 2012, MNRAS, 424, 93. doi:10.1111/j.1365-2966.2012.21171.x
  * Haskell (2015) Haskell, B. 2015, International Journal of Modern Physics E, 24, 1541007. doi:10.1142/S0218301315410074
  * Hessels et al. (2006) Hessels, J. W. T., Ransom, S. M., Stairs, I. H., et al. 2006, Science, 311, 1901. doi:10.1126/science.1123430
  * Ho & Lai (2000) Ho, W. C. G. & Lai, D. 2000, ApJ, 543, 386. doi:10.1086/317085
  * Ho et al. (2011) Ho, W. C. G., Andersson, N., & Haskell, B. 2011, Phys. Rev. Lett., 107, 101101. doi:10.1103/PhysRevLett.107.101101
  * Huang et al. (2020) Huang, K., Hu, J., Zhang, Y., et al. 2020, arXiv:2008.04491
  * Idrisy et al. (2015) Idrisy, A., Owen, B. J., & Jones, D. I. 2015, Phys. Rev. D, 91, 024001. doi:10.1103/PhysRevD.91.024001
  * Kokkotas & Schwenzer (2016) Kokkotas, K. D. & Schwenzer, K. 2016, European Physical Journal A, 52, 38. doi:10.1140/epja/i2016-16038-9
  * Krüger & Kokkotas (2020) Krüger, C. J. & Kokkotas, K. D. 2020, Phys. Rev. Lett., 125, 111106. doi:10.1103/PhysRevLett.125.111106
  * Lattimer et al. (1991) Lattimer, J. M., Pethick, C. J., Prakash, M., et al. 1991, Phys. Rev. Lett., 66, 2701. doi:10.1103/PhysRevLett.66.2701
  * Lattimer & Prakash (2000) Lattimer, J. M. & Prakash, M. 2000, Phys. Rep., 333, 121. doi:10.1016/S0370-1573(00)00019-3
  * Li et al. (2016) Li, A., Zhang, B., Zhang, N.-B., et al. 2016, Phys. Rev. D, 94, 083010. doi:10.1103/PhysRevD.94.083010
  * Li et al. (2017) Li, A., Zhu, Z.-Y., & Zhou, X. 2017, ApJ, 844, 41. doi:10.3847/1538-4357/aa7a00
  * Li et al. (2020) Li, A., Zhu, Z.-Y., Zhou, E.-P., et al. 2020, arXiv:2007.05116
  * Li & Han (2013) Li, B.-A. & Han, X. 2013, Physics Letters B, 727, 276. doi:10.1016/j.physletb.2013.10.006
  * Li et al. (2014) Li, B.-A., Ramos, À., Verde, G., et al. 2014, European Physical Journal A, 50, 9. doi:10.1140/epja/i2014-14009-x
  * Li et al. (2019) Li, B.-A., Krastev, P. G., Wen, D.-H., et al. 2019, European Physical Journal A, 55, 117. doi:10.1140/epja/i2019-12780-8
  * Li & Magno (2020) Li, B.-A. & Magno, M. 2020, Phys. Rev. C, 102, 045807. doi:10.1103/PhysRevC.102.045807
  * Lindblom et al. (1998) Lindblom, L., Owen, B. J., & Morsink, S. M. 1998, Phys. Rev. Lett., 80, 4843. doi:10.1103/PhysRevLett.80.4843
  * Lindblom et al. (2000) Lindblom, L., Owen, B. J., & Ushomirsky, G. 2000, Phys. Rev. D, 62, 084030. doi:10.1103/PhysRevD.62.084030
  * Madsen (2000) Madsen, J. 2000, Phys. Rev. Lett., 85, 10. doi:10.1103/PhysRevLett.85.10
  * Miller et al. (2019) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJ, 887, L24. doi:10.3847/2041-8213/ab50c5
  * Most et al. (2020) Most, E. R., Papenfort, L. J., Weih, L. R., et al. 2020, MNRAS, 499, L82. doi:10.1093/mnrasl/slaa168
  * Moustakidis (2015) Moustakidis, C. C. 2015, Phys. Rev. C, 91, 035804. doi:10.1103/PhysRevC.91.035804
  * Newton & Crocombe (2020) Newton, W. G. & Crocombe, G. 2020, arXiv:2008.00042
  * Oertel et al. (2017) Oertel, M., Hempel, M., Klähn, T., et al. 2017, Reviews of Modern Physics, 89, 015007. doi:10.1103/RevModPhys.89.015007
  * Ofengeim et al. (2019) Ofengeim, D. D., Gusakov, M. E., Haensel, P., et al. 2019, Phys. Rev. D, 100, 103017. doi:10.1103/PhysRevD.100.103017
  * Papazoglou & Moustakidis (2016) Papazoglou, M. C. & Moustakidis, C. C. 2016, Ap&SS, 361, 98. doi:10.1007/s10509-016-2692-5
  * Providência et al. (2019) Providência, C., Fortin, M., Pais, H., et al. 2019, Frontiers in Astronomy and Space Sciences, 6, 13. doi:10.3389/fspas.2019.00013
  * Raaijmakers et al. (2019) Raaijmakers, G., Riley, T. E., Watts, A. L., et al. 2019, ApJ, 887, L22. doi:10.3847/2041-8213/ab451a
  * Rieutord (2001) Rieutord, M. 2001, ApJ, 550, 443. doi:10.1086/319705
  * Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJ, 887, L21. doi:10.3847/2041-8213/ab481c
  * Roupas et al. (2020) Roupas, Z., Panotopoulos, G., & Lopes, I. 2020, arXiv:2010.11020
  * Sá (2004) Sá, P. M. 2004, Phys. Rev. D, 69, 084001. doi:10.1103/PhysRevD.69.084001
  * Sá & Tomé (2005) Sá, P. M. & Tomé, B. 2005, Phys. Rev. D, 71, 044007. doi:10.1103/PhysRevD.71.044007
  * Safarzadeh & Loeb (2020) Safarzadeh, M. & Loeb, A. 2020, ApJ, 899, L15. doi:10.3847/2041-8213/aba9df
  * Shternin & Yakovlev (2008) Shternin, P. S. & Yakovlev, D. G. 2008, Phys. Rev. D, 78, 063006. doi:10.1103/PhysRevD.78.063006
  * Steiner et al. (2005) Steiner, A. W., Prakash, M., Lattimer, J. M., et al. 2005, Phys. Rep., 411, 325. doi:10.1016/j.physrep.2005.02.004
  * Stergioulas et al. (1995) Stergioulas, N. et al., 1995, RNS code, http://www.gravity.phys.uwm.edu/rns/
  * Tews et al. (2020) Tews, I., Pang, P. T. H., Dietrich, T., et al. 2020, arXiv:2007.06057
  * Vidaña et al. (2009) Vidaña, I., Providência, C., Polls, A., et al. 2009, Phys. Rev. C, 80, 045806. doi:10.1103/PhysRevC.80.045806
  * Vidaña (2012) Vidaña, I. 2012, Phys. Rev. C, 85, 045808. doi:10.1103/PhysRevC.85.045808
  * Wang et al. (2019) Wang, Y.-B., Zhou, X., Wang, N., et al. 2019, Research in Astronomy and Astrophysics, 19, 030. doi:10.1088/1674-4527/19/2/30
  * Wen et al. (2012) Wen, D.-H., Newton, W. G., & Li, B.-A. 2012, Phys. Rev. C, 85, 025801. doi:10.1103/PhysRevC.85.025801
  * Wijngaarden et al. (2019) Wijngaarden, M. J. P., Ho, W. C. G., Chang, P., et al. 2019, MNRAS, 484, 974. doi:10.1093/mnras/stz042
  * Xie & Li (2020) Xie, W.-J. & Li, B.-A. 2020, ApJ, 899, 4. doi:10.3847/1538-4357/aba271
  * Xu et al. (2009) Xu, J., Chen, L.-W., Li, B.-A., et al. 2009, ApJ, 697, 1549. doi:10.1088/0004-637X/697/2/1549
  * Yang et al. (2011) Yang, S.-H., Pi, C.-M., & Zheng, X.-P. 2011, ApJ, 735, L29. doi:10.1088/2041-8205/735/2/L29
  * Zhang et al. (2018) Zhang, N.-B., Li, B.-A., & Xu, J. 2018, ApJ, 859, 90. doi:10.3847/1538-4357/aac027
  * Zhang & Mann (2020) Zhang, C. & Mann, R. B. 2020, arXiv:2009.07182
  * Zhang & Li (2019) Zhang, N.-B. & Li, B.-A. 2019, European Physical Journal A, 55, 39. doi:10.1140/epja/i2019-12700-0
  * Zhang & Li (2020) Zhang, N.-B. & Li, B.-A. 2020, ApJ, 902, 38. doi:10.3847/1538-4357/abb470