# Vacuum and singularity formation for compressible Euler equations with time-
dependent damping ††thanks: This work is supported in part by the National
Natural Science Foundation of China (Grant No. 11671237)

Ying Sui a Weiqiang Wang b Huimin Yu c
a,c School of Mathematics and Statistics, Shandong Normal University, Jinan
250014, China.
b Academy of Mathematics and Systems Science, Chinese Academy of Sciences,
Beijing 100190, China.
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>author,
Email<EMAIL_ADDRESS>
Abstract: In this paper, vacuum and singularity formation are considered for
compressible Euler equations with time-dependent damping. For $1<\gamma\leq
3$, by constructing some new control functions ingeniously, we obtain the
lower bounds estimates on density for arbitrary classical solutions. Basing on
these lower estimates, we succeed in proving the singular formation theorem
for all $\lambda$, which was open in [1] for some cases. Moreover, the
singularity formation of the compressible Euler equations when $\gamma=3$ is
investigated, too.

Keywords: Singularity formation, compressible Euler equations, time-dependent
damping, shock wave.

##  1 Introduction

We consider the one dimensional compressible Euler equations with time-
dependent damping in Lagrangian coordinates:

$\left\\{\begin{aligned} &\tau_{t}-u_{x}=0,\\\
&u_{t}+p_{x}=-\frac{\alpha}{(1+t)^{\lambda}}u,\end{aligned}\right.$ (1.1)

with the initial data

$\tau(x,0)=\tau_{0}(x),\;u(x,0)=u_{0}(x),$ (1.2)

where $\tau(x,t)$ and $u(x,t)$ are the specific volume and velocity of the
flow respectively at location $x\in\mathbb{R}$ and time $t\in\mathbb{R}_{+}$.
For simplicity, we assume the gas is deal polytropic and the gas pressure

$p=K\tau^{-\gamma},~{}~{}\text{~{}for~{}}K>0\text{~{}and~{}}\gamma>1,$ (1.3)

where $\gamma$ is called adiabatic index. Besides, $\alpha\geq 0,\
\lambda\in\mathbb{R}$ are two constants, and the term
$-\frac{\alpha}{(1+t)^{\lambda}}u$ is the so-called damping effect on the
fluid when $\alpha>0$.

For $\alpha=0$, the system (1.1) is the classical compressible Euler equation,
which are the most basic model of hyperbolic conservation law equations and
can be used to describe many physical phenomena. It is well known that, due to
nonlinearity, classical solutions of hyperbolic conservation laws, may form
gradient blowup in finite time, even when initial data is smooth and small.
The study on breakdown of classical solutions for hyperbolic conservation laws
has a long history [2]. One can trace back to Stokes in [3] for a breakdown
example for some scalar equation. In 1964, Lax [4] studied the formation of
singularities in finite time for $2\times 2$ general systems of strictly
hyperbolic conservation laws when initial data is small and contains some
compression, then followed by John, Li and Liu in [5, 6, 7] and etc., for more
general systems of conservation laws. For $\alpha>0$, the system (1.1) can be
used to describe the flow of fluids through porous media. However, the time-
dependent damping makes some fantastic variety of the dynamic system.
Therefore, the study on this problem is full of challenges and quite
incomplete. Due to the practicability and strong relation with Euler equation,
it has attracted extensive attention of many scholars in recent years, see [8,
9, 10, 11, 12, 13, 14, 15, 16, 17] for example.

In this paper, we study the vacuum and singularity formation of compressible
Euler equations with time-dependent damping. In particular, we extend the
results and give the proofs of some unsolved cases in [1]. That is, we
consider the cases of
$\lambda\in\left(\min\\{1,\frac{\alpha(\gamma-1)}{\gamma-3}\\},\max\\{1,\frac{\alpha(\gamma-1)}{\gamma-3}\\}\right)$
when $\gamma>3$ and $\lambda<\frac{\alpha(\gamma-1)}{\gamma-3}$ when
$1<\gamma<3$, which are all open in [1]. We also investigate the case of
$\gamma=3$ and give the corresponding breakdown results. Here, we still use
the idea of constructing decoupled Riccati type equations

$y^{\prime}=a_{0}-a_{2}y^{2}$

for smooth solutions of compressible Euler equations. The difference between
the two papers is that $a_{0}$ can change sign in this paper, while $a_{0}$
must be positive or negative in [1]. Exactly speaking, for the unsolved case
in [1], there exists a $t_{0}>0$ such that $a_{0}<0$ for $t>t_{0}$, while
$a_{0}>0$ for $0<t<t_{0}$. To overcome this difficulty, we firstly show the
time-dependent lower bound on density for all $\lambda\in\mathbb{R}$ by
constructing some new functions. The lower bound on density implies that
$a_{0}$ is uniformly bounded for $0<t<t_{0}$. Then using the uniform bound on
$a_{0}$, we succeed in obtaining the singularity formation on the derivatives
of velocity or specific volume if the initial compression reaches a certain
level. Especially, our main results indicate that, for $1<\gamma\leq 3$,
vacuum never occurs at any finite time provided the initial data away from
vacuum, while the singularity must appear once the initial compression reaches
a certain level, even for the over-damping case $\lambda<0.$

##  2 The main Theorems

For convenience, let’s first introduce some variables and notations which have
been used in [1]. Denote

$\phi:=\int_{\tau}^{\infty}c(\tau)d\tau=\frac{2\sqrt{K\gamma}}{\gamma-1}\tau^{-\frac{\gamma-1}{2}}>0,$
(2.1)

where the nonlinear Lagrangian sound speed $c$ is defined as

$c:=\sqrt{-p_{\tau}}=\sqrt{K\gamma}\tau^{-\frac{\gamma+1}{2}}.$ (2.2)

It follows from (2.1)-(2.2) that

$\tau=K_{\tau}\phi^{-\frac{2}{\gamma-1}},\ \
p=K_{p}\phi^{\frac{2\gamma}{\gamma-1}},\ \
c=K_{c}\phi^{\frac{\gamma+1}{\gamma-1}},$

where $K_{\tau},K_{p}$ and $K_{c}$ are positive constants given by

$K_{\tau}:=\left(\frac{2\sqrt{K\gamma}}{\gamma-1}\right)^{\frac{2}{\gamma-1}},\
\ K_{p}:=KK^{-\gamma},\ \ \text{ and }\ \
K_{c}:=\sqrt{K\gamma}K_{\tau}^{-\frac{\gamma+1}{2}}.$

A direct calculation shows that the system (1.1) has two characteristic speeds

$\lambda_{1}=-\lambda_{2}=c.$

The forward and backward characteristic lines are described by

$\frac{dx_{+}(t)}{dt}=c\quad\text{ and }\quad\frac{dx_{-}(t)}{dt}=-c,$

and we denote the corresponding directional derivatives along them by

$\prime=\frac{\partial}{\partial t}+c\frac{\partial}{\partial x}\quad\text{
and }\quad\backprime=\frac{\partial}{\partial t}-c\frac{\partial}{\partial x}$

respectively. Furthermore, we denote the Riemann invariants by

$w:=u+\phi\quad\text{ and }\quad z:=u-\phi.$ (2.3)

By a direct calculation, we have

$w^{\prime}=-\frac{\alpha}{2(1+t)^{\lambda}}(z+w),\qquad
z^{\backprime}=-\frac{\alpha}{2(1+t)^{\lambda}}(z+w).$ (2.4)

In [1], the authors have proved the following theorems.

######  Theorem 2.1 ([1, Theorem 3.2]).

Suppose the initial data $(\tau_{0},u_{0})(x)\in C^{1}(\mathbb{R})$ and there
exists a positive constant $C_{0}$ such that

$\|u_{0}(x)\|_{C^{1}}\leq C_{0},~{}\|\tau_{0}(x)\|_{C^{1}}\leq
C_{0},~{}C_{0}^{-1}<\tau_{0}(x).$ (2.5)

Let $1<\gamma<3$, $\lambda\geq\frac{\alpha(\gamma-1)}{\gamma-3}$ and
$(\tau,u)$ is a $C^{1}$ solution of (1.1)-(1.2). Assume there exists one point
$x_{0}$ such that

$u_{x}(x_{0},0)+\phi_{x}(x_{0},0)<-\frac{\alpha(\gamma-1)}{K_{c}(3-\gamma)}(\phi(x_{0},0))^{\frac{-2}{\gamma-1}}$
(2.6)

or

$u_{x}(x_{0},0)-\phi_{x}(x_{0},0)<-\frac{\alpha(\gamma-1)}{K_{c}(3-\gamma)}(\phi(x_{0},0))^{\frac{-2}{\gamma-1}}.$
(2.7)

Then $\|(\tau_{x},u_{x})\|_{L^{\infty}}$ must blow up in finite time.

######  Theorem 2.2 ([1, Theorem 3.1]).

Suppose the initial data $(\tau_{0},u_{0})(x)\in C^{1}(\mathbb{R})$ and there
exists a positive constant $C_{0}$ such that

$|u_{0}(x)|\leq C_{0},\quad C_{0}^{-1}<\tau_{0}(x).$ (2.8)

Let $\gamma>3,\lambda<\min\\{1,\frac{\alpha(\gamma-1)}{\gamma-3}\\}$ or
$\lambda>\max\\{1,\frac{\alpha(\gamma-1)}{\gamma-3}\\}$, and $(\tau,u)$ is a
$C^{1}$ solution of (1.1)-(1.2). Assume there exists one point $x_{0}$ such
that

$u_{x}(x_{0},0)+\phi_{x}(x_{0},0)<K_{1}(\phi(x_{0},0))^{-\frac{2}{\gamma-1}}-K_{2}(\phi(x_{0},0))^{-\frac{\gamma+1}{2(\gamma-1)}}$
(2.9)

or

$u_{x}(x_{0},0)-\phi_{x}(x_{0},0)<K_{1}(\phi(x_{0},0))^{-\frac{2}{\gamma-1}}-K_{2}(\phi(x_{0},0))^{-\frac{\gamma+1}{2(\gamma-1)}},$
(2.10)

where $K_{1},K_{2}$ are positive constants which depend only on
$C_{0},\gamma,K,\lambda$ and $\alpha$. Then
$\|(\tau_{x},u_{x})\|_{L^{\infty}}$ must blow up in finite time.

For the rest cases, that is,

$\lambda<\frac{\alpha(\gamma-1)}{\gamma-3}\text{ when }1<\gamma<3,$ (2.11)

and

$\lambda\in\left(\min\\{1,\frac{\alpha(\gamma-1)}{\gamma-3}\\},\max\\{1,\frac{\alpha(\gamma-1)}{\gamma-3}\\}\right)\text{
when }\gamma>3,$ (2.12)

the authors assert that this is just a technical problem and the Riemann
invariants may blow up for some certain initial data, too. Indeed, we consider
the rest cases in this paper and have obtained the following theorems.

######  Theorem 2.3 (Singularity formation for $1<\gamma<3$ and $\lambda\neq
1$).

Suppose the initial data $(\tau_{0},u_{0})(x)\in C^{1}(\mathbb{R})$ satisfying
(2.5) and $(\tau,u)$ is a $C^{1}$ solution of (1.1)-(1.2). Then
$\|u\|_{L^{\infty}}$ will be uniformly bounded for all time $t>0$ and the
density $\rho=\tau^{-1}$ has a time-dependent lower bound,

$\displaystyle\rho(x,t)=\tau^{-1}(x,t)\geq{\left\\{\begin{matrix}C\left(1+t\right)^{\frac{4}{\gamma-3}}\quad\quad\quad\quad\quad&\quad\text{for
}\lambda\geq 0,\\\
C\left(1+e^{\frac{-\lambda(1+t)^{1-\lambda}}{1-\lambda}}t\right)^{\frac{4}{\gamma-3}}&\quad\text{for
}\lambda<0,\end{matrix}\right.}$ (2.13)

where $C$ is a positive constants depending only on $C_{0},\gamma,K,\lambda$
and $\alpha$. (2.13) indicates that the vacuum will never occur in finite
time. Moreover, assume there exists one point $x_{0}$ such that

$u_{x}(x_{0},0)+\phi_{x}(x_{0},0)<K_{3}(\phi(x_{0},0))^{-\frac{2}{\gamma-1}}-K_{4}(\phi(x_{0},0))^{-\frac{\gamma+1}{2(\gamma-1)}}$
(2.14)

or

$u_{x}(x_{0},0)-\phi_{x}(x_{0},0)<K_{3}(\phi(x_{0},0))^{-\frac{2}{\gamma-1}}-K_{4}(\phi(x_{0},0))^{-\frac{\gamma+1}{2(\gamma-1)}},$
(2.15)

where $K_{3},K_{4}$ are positive constants which depend only on
$C_{0},\gamma,K,\lambda$ and $\alpha$. Then
$\|(\tau_{x},u_{x})\|_{L^{\infty}}$ must blow up in finite time.

######  Theorem 2.4 (Singularity formation for $\gamma>3$ and $\lambda\neq
1$).

Suppose the initial data $(\tau_{0},u_{0})(x)\in C^{1}(\mathbb{R})$ satisfying
(2.8) and $(\tau,u)$ be a $C^{1}$ solution of (1.1)-(1.2). Assume there exists
one point $x_{0}$ such that

$u_{x}(x_{0},0)+\phi_{x}(x_{0},0)<K_{5}(\phi(x_{0},0))^{-\frac{2}{\gamma-1}}-K_{6}(\phi(x_{0},0))^{-\frac{\gamma+1}{2(\gamma-1)}}$
(2.16)

or

$u_{x}(x_{0},0)-\phi_{x}(x_{0},0)<K_{5}(\phi(x_{0},0))^{-\frac{2}{\gamma-1}}-K_{6}(\phi(x_{0},0))^{-\frac{\gamma+1}{2(\gamma-1)}},$
(2.17)

where $K_{5},K_{6}$ are two positive constants which depend only on
$C_{0},\gamma,K,\lambda$ and $\alpha$. Then
$\|(\tau_{x},u_{x})\|_{L^{\infty}}$ must blow up in finite time.

######  Remark 2.5.

(2.13) holds for all $\lambda\in\mathbb{R}$, and the condition $\lambda\neq 1$
in Theorem 2.3 is only used in proving the blow-up of
$\|(\tau_{x},u_{x})\|_{L^{\infty}}$. Singularity formation for $\lambda=1$
with all $\gamma>1$ can be found in [1, Theorems 4.1 and 4.2], also see the
Remark 3.1 below.

Theorem 2.3 and Theorem 2.4 remove the restrictions on $\lambda$ in [1] and
include the cases $(\ref{a2})$ and $(\ref{a1})$. In addition, the singularity
formation of the compressible Euler equations when $\gamma=3$ is given in the
following theorem.

######  Theorem 2.6 (Singularity formation for $\gamma=3$).

Suppose the initial data $(\tau_{0},u_{0})(x)\in C^{1}(\mathbb{R})$ satisfying
(2.5) and $(\tau,u)(x,t)$ is a $C^{1}$ solution of $\eqref{1.1}-\eqref{1.2}$
with $\gamma=3$. Then $\|u\|_{L^{\infty}}$ will be uniformly bounded for all
time $t>0$ and the density $\rho=\tau^{-1}$ has a time-dependent lower bound,

$\displaystyle\rho(x,t)=\tau^{-1}(x,t)\geq{\left\\{\begin{matrix}C_{0}^{-1}e^{-\hat{M}t}\quad\quad\quad\quad&\quad\quad\text{for
}\lambda\geq 0,\\\
C_{0}^{-1}e^{-\hat{M}te^{\frac{-\lambda(1+t)^{1-\lambda}}{1-\lambda}}}&\quad\quad\text{for
}\lambda<0,\end{matrix}\right.}$ (2.18)

where $\hat{M}$ is a positive constant depending only on
$C_{0},K,\gamma,\lambda$ and $\alpha$. (2.18) indicates that the vacuum will
never occur in finite time. Moreover, for any given $t_{0}>0$, there exists a
non-decreasing (w.r.t. $t_{0}$) function $\tilde{M}(t_{0})$, depending on
$C_{0},K,\gamma,\lambda$ and $\alpha$, such that if there exists one point
$x_{0}$ satisfying

$\displaystyle
u_{x}(x_{0},0)+\phi_{x}(x_{0},0)<\frac{\alpha}{2K_{c}}\phi^{-1}(x_{0},0)\ln\phi(x_{0},0)-\max\\{\frac{2}{K_{c}t_{0}},\tilde{M}(t_{0})\\},$
(2.19)

or

$\displaystyle
u_{x}(x_{0},0)-\phi_{x}(x_{0},0)<\frac{\alpha}{2K_{c}}\phi^{-1}(x_{0},0)\ln\phi(x_{0},0)-\max\\{\frac{2}{K_{c}t_{0}},\tilde{M}(t_{0})\\}.$
(2.20)

Then $\|u_{x},\tau_{x}\|_{L^{\infty}}$ must blow up before $t_{0}$.

######  Remark 2.7.

We mention that we only need $\eqref{1.10}$ to obtain the time-dependent lower
bound on density (2.13) for $1<\gamma<3$, and (2.18) for $\gamma=3$.
Therefore, (2.13) and (2.18) may be viewed as a generalization of [1, Lemma
3.3] and [17, Lemma 12], in which the authors impose some restrictions on
$\lambda$ and $\alpha$.

######  Remark 2.8.

For $1<\gamma\leq 3$ and $\lambda=0$, we can get the lower bound of density
independent of $t$. In other words, using the proof of Theorem 2.3 and Theorem
2.6, we can respectively obtain the fixed lower bound of density for
$1<\gamma\leq 3$ and $\gamma=3$ by simple calculation.

######  Remark 2.9.

The time-dependent lower bound on density plays an important role in obtaining
blow-up of $\|(\tau_{x},u_{x})\|_{L^{\infty}}$ in Theorems 2.3 and 2.6. For
$\gamma>1$, under the assumption (2.5), an optimal time-dependent lower bound
on density is obtained in [18, 19] for compressible isentropic Euler equations
and in [20] for full compressible Euler equations. However, for $\gamma\geq
3$, as stated in [17, Remark 3], it is open whether vacuum occurs in finite
time or not for damped Euler equations. Theorem 2.6 give a negative answer to
this open question for $\gamma=3$, i.e., the vacuum will never occur in finite
time.

For latter use, we recall some useful lemmas, whose proof can be found in [1].

######  Lemma 2.10 ([1, Theorem 2.1]).

Let $(\tau_{0},u_{0})(x)\in C^{1}(\mathbb{R})$ satisfying (2.8) and $(\tau,u)$
be the $C^{1}$ solution to system (1.1)-(1.2). Then in the domain where the
solution exists, $(\tau,u)(x,t)$ is uniform bounded in the following form

$\displaystyle|u(x,t)|\leq\tilde{C}_{0},\quad\tau(x,t)\geq\tilde{C}_{0}^{-1},$
(2.21)

where $\tilde{C}_{0}$ can be any constant bigger than
$\max\\{C_{0}+C_{0}^{\frac{\gamma-1}{2}},(C_{0}+C_{0}^{\frac{\gamma-1}{2}})^{\frac{2}{\gamma-1}}\\}$.

######  Lemma 2.11 ([1, Lemma 2.2]).

The $C^{1}$ solutions $(\tau,u)$ of (1.1)-(1.2) satisfy

$\displaystyle{A}^{\prime}=-\frac{\alpha}{2(1+t)^{\lambda}}({A}+{B})+K_{d}({A}{B}-{A}^{2})$

and

$\displaystyle{B}^{\backprime}=-\frac{\alpha}{2(1+t)^{\lambda}}({A}+{B})+K_{d}({A}{B}-{B}^{2}),$

where

${A}={w}_{x},\quad{B}={z}_{x},\quad
K_{d}=K_{c}\frac{\gamma+1}{2(\gamma-1)}\phi^{\frac{2}{\gamma-1}}.$

######  Lemma 2.12 ([1, Lemma 3.1]).

For the $C^{1}$ solutions $(\tau,u)$ of (1.1)-(1.2), denoting

$y:=\left(\phi^{\frac{\gamma+1}{2(\gamma-1)}}w_{x}-\frac{\alpha(\gamma-1)}{K_{c}(\gamma-3)(1+t)^{\lambda}}\phi^{\frac{\gamma-3}{2(\gamma-1)}}\right)e^{\frac{\alpha(3\gamma-1)}{2(\gamma-3)(1-\lambda)}(1+t)^{-\lambda}},$
(2.22)

and

$q:=\left(\phi^{\frac{\gamma+1}{2(\gamma-1)}}z_{x}-\frac{\alpha(\gamma-1)}{K_{c}(\gamma-3)(1+t)^{\lambda}}\phi^{\frac{\gamma-3}{2(\gamma-1)}}\right)e^{\frac{\alpha(3\gamma-1)}{2(\gamma-3)(1-\lambda)}(1+t)^{-\lambda}},$
(2.23)

then $y$ and $q$ satisfy

$y^{\prime}=a_{0}-a_{2}y^{2},\quad q^{\backprime}=a_{0}-a_{2}q^{2},$ (2.24)

where

$a_{2}=\frac{K_{c}(\gamma+1)}{2(\gamma-1)}\phi^{-\frac{\gamma-3}{2(\gamma-1)}}e^{-\frac{\alpha(3\gamma-1)}{2(\gamma-3)(1-\lambda)}(1+t)^{1-\lambda}}>0,$
(2.25)

$a_{0}=\frac{\lambda\alpha(\gamma-1)(\gamma-3)(1+t)^{\lambda-1}-\alpha^{2}(\gamma-1)^{2}}{K_{c}(\gamma-3)^{2}(1+t)^{2\lambda}}\phi^{\frac{\gamma-3}{2(\gamma-1)}}e^{\frac{\alpha(3\gamma-1)}{2(\gamma-3)(1-\lambda)}(1+t)^{1-\lambda}}.$
(2.26)

##  3 Proof of the main theorems

###  3.1 Proof of Theorem 2.3

We divide the proof into five steps. In steps 1-4, we first give the proof of
the lower bound on density, that is, the proof of (2.13). In step 5, we prove
the blow-up of $\|(\tau_{x},u_{x})\|_{L^{\infty}}$.

Step 1. It follows from (2.4) and the fact $c_{w}=-c_{z}$ that

$\displaystyle(w_{x})^{\prime}$ $\displaystyle=(w^{\prime})_{x}-c_{x}w_{x}$
(3.1)
$\displaystyle=-\frac{\alpha}{2(1+t)^{\lambda}}({z}+{w})_{x}-(c_{w}w_{x}^{2}+c_{z}z_{x}w_{x})$
$\displaystyle=-\frac{\alpha}{2(1+t)^{\lambda}}({z}+{w})_{x}-c_{w}w_{x}^{2}+c_{w}z_{x}w_{x},$

and

$\displaystyle(z_{x})^{\backprime}$
$\displaystyle=(z^{\backprime})_{x}+c_{x}z_{x}$ (3.2)
$\displaystyle=-\frac{\alpha}{2(1+t)^{\lambda}}({z}+{w})_{x}+(c_{w}w_{x}z_{x}+c_{z}z_{x}^{2})$
$\displaystyle=-\frac{\alpha}{2(1+t)^{\lambda}}({z}+{w})_{x}-c_{w}z_{x}^{2}+c_{w}z_{x}w_{x}.$

Noticing from (2.1)-(2.2), it holds

$\displaystyle\phi_{\tau}=-c,$

which together with (1.1), yields

$\displaystyle\phi^{\prime}$
$\displaystyle=\phi_{t}+c\phi_{x}=-c\tau_{t}+c\phi_{x}=-cu_{x}+c\phi_{x}=-c(u-\phi)_{x}=-cz_{x},$
$\displaystyle\phi^{\backprime}$
$\displaystyle=\phi_{t}-c\phi_{x}=-c\tau_{t}-c\phi_{x}=-cu_{x}-c\phi_{x}=-c(u+\phi)_{x}=-cw_{x}.$

Thus, we have

$z_{x}=-\frac{1}{c}\phi^{\prime}=-\frac{1}{c}\phi_{\tau}\tau^{\prime}=\tau^{\prime},\qquad
w_{x}=-\frac{1}{c}\phi^{\backprime}=-\frac{1}{c}\phi_{\tau}\tau^{\backprime}=\tau^{\backprime}.$
(3.3)

Moreover, using (3.3), a direct calculation shows that

$\displaystyle\frac{1}{2}c^{-\frac{1}{2}}c^{\prime}w_{x}$
$\displaystyle=\frac{1}{2}c^{-\frac{1}{2}}(c_{w}w_{\tau}\tau^{\prime}+c_{z}z_{\tau}\tau^{\prime})w_{x}=c^{-\frac{1}{2}}c_{w}w_{\tau}\tau^{\prime}w_{x}=-c^{\frac{1}{2}}c_{w}z_{x}w_{x},$
(3.4)

where we have used $c_{w}=-c_{z}$, $w_{\tau}=-z_{\tau}$ and
$w_{\tau}=\phi_{\tau}=-c$.

Hence it follows form (3.1)-(3.2), (3.4) that

$\displaystyle(c^{\frac{1}{2}}w_{x})^{\prime}$
$\displaystyle=c^{\frac{1}{2}}(w_{x})^{\prime}+\frac{1}{2}c^{-\frac{1}{2}}c^{\prime}w_{x}$
(3.5) $\displaystyle=-\frac{\alpha
c^{\frac{1}{2}}}{2(1+t)^{\lambda}}({z}+{w})_{x}-c^{\frac{1}{2}}c_{w}w_{x}^{2}+c^{\frac{1}{2}}c_{w}z_{x}w_{x}+\frac{1}{2}c^{-\frac{1}{2}}c^{\prime}w_{x}.$
$\displaystyle=-\frac{\alpha
c^{\frac{1}{2}}}{2(1+t)^{\lambda}}({z}+{w})_{x}-c^{\frac{1}{2}}c_{w}w_{x}^{2},$

and

$\displaystyle(c^{\frac{1}{2}}z_{x})^{\backprime}$
$\displaystyle=c^{\frac{1}{2}}(z_{x})^{\backprime}+\frac{1}{2}c^{-\frac{1}{2}}c^{\backprime}z_{x}$
(3.6) $\displaystyle=-\frac{\alpha
c^{\frac{1}{2}}}{2(1+t)^{\lambda}}({z}+{w})_{x}-c^{\frac{1}{2}}c_{w}z_{x}^{2}+c^{\frac{1}{2}}c_{w}z_{x}w_{x}+\frac{1}{2}c^{-\frac{1}{2}}c^{\backprime}z_{x}$
$\displaystyle=-\frac{\alpha
c^{\frac{1}{2}}}{2(1+t)^{\lambda}}({z}+{w})_{x}-c^{\frac{1}{2}}c_{w}z_{x}^{2}.$

Step 2. For $\lambda\geq 0$ and $\lambda\neq 1$, we introduce two new
variables

$G(t,x)=c^{1\over 2}w_{x}+{\alpha\over{(1+t)^{\lambda}}}h(\tau),$ (3.7)
$H(t,x)=c^{1\over 2}z_{x}+{\alpha\over{(1+t)^{\lambda}}}h(\tau),$ (3.8)

where

$h(\tau)=\int^{\tau}_{0}c(s)^{\frac{1}{2}}ds=\frac{4(K\gamma)^{\frac{1}{4}}}{3-\gamma}\tau^{\frac{3-\gamma}{4}}.$
(3.9)

Using (3.3), (3.5)-(3.6), we obtain

$\displaystyle G^{\prime}$
$\displaystyle=(c^{\frac{1}{2}}w_{x})^{\prime}+\frac{\alpha}{(1+t)^{\lambda}}h^{\prime}(\tau)-\frac{\lambda\alpha}{(1+t)^{\lambda+1}}h(\tau)$
(3.10) $\displaystyle=-\frac{\alpha
c^{\frac{1}{2}}}{2(1+t)^{\lambda}}({z}+{w})_{x}-c^{1\over
2}c_{w}w_{x}^{2}+\frac{\alpha}{(1+t)^{\lambda}}c^{\frac{1}{2}}z_{x}-\frac{\lambda\alpha}{(1+t)^{\lambda+1}}h(\tau)$
$\displaystyle=\frac{\alpha}{2(1+t)^{\lambda}}(H-G)-c^{1\over
2}c_{w}w_{x}^{2}-{\lambda\alpha\over{(1+t)^{\lambda+1}}}h(\tau)$
$\displaystyle:=-\frac{\alpha}{2(1+t)^{\lambda}}(G-H)+R_{1},$

and

$\displaystyle H^{\backprime}$
$\displaystyle=(c^{\frac{1}{2}}z_{x})^{\backprime}+\frac{\alpha}{(1+t)^{\lambda}}h^{\backprime}(\tau)-\frac{\lambda\alpha}{(1+t)^{\lambda+1}}h(\tau)$
(3.11) $\displaystyle=-\frac{\alpha
c^{\frac{1}{2}}}{2(1+t)^{\lambda}}({z}+{w})_{x}-c^{1\over
2}c_{w}z_{x}^{2}+\frac{\alpha}{(1+t)^{\lambda}}c^{\frac{1}{2}}w_{x}-\frac{\lambda\alpha}{(1+t)^{\lambda+1}}h(\tau)$
$\displaystyle=\frac{\alpha}{2(1+t)^{\lambda}}(G-H)-c^{1\over
2}c_{w}z_{x}^{2}-{\lambda\alpha\over{(1+t)^{\lambda+1}}}h(\tau)$
$\displaystyle:=-\frac{\alpha}{2(1+t)^{\lambda}}(H-G)+R_{2}.$

Noticing from (2.2) and (2.3) that
$c_{w}=\frac{\gamma+1}{2(\gamma-1)}K_{c}\phi^{\frac{2}{\gamma-1}}\geq 0$, then

$R_{1}:=-c^{1\over
2}c_{w}w_{x}^{2}-{\alpha\lambda\over{(1+t)^{\lambda+1}}}h(\tau)\leq 0,$
$R_{2}:=-c^{1\over
2}c_{w}z_{x}^{2}-{\alpha\lambda\over{(1+t)^{\lambda+1}}}h(\tau)\leq 0,$

which together with (3.10)-(3.11), yields

$G^{\prime}\leq-\frac{\alpha}{2(1+t)^{\lambda}}(G-H),$
$H^{\backprime}\leq-\frac{\alpha}{2(1+t)^{\lambda}}(H-G).$

It follows from (2.1), (2.3), (2.5) and (3.7)-(3.9) that there exists a
positive constant $M$ depending only on $C_{0},K,\gamma,\lambda$ and $\alpha$
such that $\sup_{x\in\mathbb{R}}|(G(0,x),H(0,x))|\leq M$. Now let

$G_{1}(t,x)=G(t,x)-M,\quad H_{1}(t,x)=H(t,x)-M,$ (3.12)

then $G_{1}(0,x)\leq 0,H_{1}(0,x)\leq 0$ and

$\displaystyle
G_{1}^{\prime}={G}^{\prime}\leq-\frac{\alpha}{2(1+t)^{\lambda}}(G-H)=-\frac{\alpha}{2(1+t)^{\lambda}}G_{1}+\frac{\alpha}{2(1+t)^{\lambda}}H_{1},$
(3.13) $\displaystyle
H_{1}^{\backprime}={H}^{\backprime}\leq-\frac{\alpha}{2(1+t)^{\lambda}}(H-G)=-\frac{\alpha}{2(1+t)^{\lambda}}H_{1}+\frac{\alpha}{2(1+t)^{\lambda}}G_{1}.$

Since $\lambda\neq 1$, we denote
$A(t)=e^{\frac{\alpha(1+t)^{1-\lambda}}{2(1-\lambda)}}$ and

$G_{2}(t,x)=A(t)G_{1}(t,x),\quad H_{2}(t,x)=A(t)H_{1}(t,x),$

then $G_{2}(0,x)\leq 0,H_{2}(0,x)\leq 0$ and

$\displaystyle
G_{2}^{\prime}=A(t)G_{1}^{\prime}+\frac{\alpha}{2(1+t)^{\lambda}}A(t)G_{1}\leq\frac{\alpha}{2(1+t)^{\lambda}}A(t)H_{1}=\frac{\alpha}{2(1+t)^{\lambda}}H_{2},$
(3.14) $\displaystyle
H_{2}^{\backprime}=A(t)H_{1}^{\backprime}+\frac{\alpha}{2(1+t)^{\lambda}}A(t)H_{1}\leq\frac{\alpha}{2(1+t)^{\lambda}}A(t)G_{1}=\frac{\alpha}{2(1+t)^{\lambda}}G_{2}.$

Using a continuity argument, it is directly from (3.14) to see

$G_{2}(x_{+}(t),t)\leq 0,\quad H_{2}(x_{-}(t),t)\leq 0,\quad\text{for all
}t\geq 0.$

Therefore,

$G_{1}(x_{+}(t),t)\leq 0,\quad H_{1}(x_{-}(t),t)\leq 0,\quad\text{for all
}t\geq 0,$

which together with (3.12) yields that

$G(t,x)\leq M,\quad H(t,x)\leq M,\quad\text{for all
}(t,x)\in\mathbb{R}_{+}\times\mathbb{R}.$ (3.15)

Hence by (3.7), (3.8) and (3.15), we obtain

$G(t,x)+H(t,x)=c^{1\over
2}(w+z)_{x}+{2\alpha\over{(1+t)^{\lambda}}}h(\tau)\leq 2M,$

which together with $\eqref{1.1}_{1}$ implies that

$c^{1\over 2}\tau_{t}+{\alpha\over{(1+t)^{\lambda}}}h(\tau)\leq M.$ (3.16)

Combing (2.2), (3.9) with (3.16), we have

$(K\gamma)^{1\over
4}\tau^{-\frac{\gamma+1}{4}}\tau_{t}+{\alpha\over{(1+t)^{\lambda}}}\frac{4(K\gamma)^{\frac{1}{4}}}{3-\gamma}\tau^{\frac{3-\gamma}{4}}\leq
M,$

which implies that

$\left(e^{\frac{\alpha(1+t)^{1-\lambda}}{1-\lambda}}\tau^{\frac{3-\gamma}{4}}\right)_{t}\leq\left(\frac{3-\gamma}{4}(K\gamma)^{-\frac{1}{4}}M\right)e^{\frac{\alpha(1+t)^{1-\lambda}}{1-\lambda}}.$
(3.17)

Integrating (3.17) from $0$ to $t$, one obtains

$\displaystyle
e^{\frac{\alpha(1+t)^{1-\lambda}}{1-\lambda}}\tau^{\frac{3-\gamma}{4}}-e^{\frac{\alpha}{1-\lambda}}\tau_{0}^{\frac{3-\gamma}{4}}$
$\displaystyle\leq\frac{3-\gamma}{4}(K\gamma)^{-\frac{1}{4}}M\int_{0}^{t}e^{\frac{\alpha(1+s)^{1-\lambda}}{1-\lambda}}ds$
(3.18)
$\displaystyle\leq\frac{3-\gamma}{4}(K\gamma)^{-\frac{1}{4}}Me^{\frac{\alpha(1+t)^{1-\lambda}}{1-\lambda}}t,$

therefore

$\displaystyle\tau^{\frac{3-\gamma}{4}}$
$\displaystyle\leq\left[e^{\frac{\alpha}{1-\lambda}}\tau_{0}^{\frac{3-\gamma}{4}}+\frac{3-\gamma}{4}(K\gamma)^{-\frac{1}{4}}Me^{\frac{\alpha(1+t)^{1-\lambda}}{1-\lambda}}t\right]e^{-\frac{\alpha(1+t)^{1-\lambda}}{1-\lambda}}$
(3.19)
$\displaystyle=e^{\frac{\alpha(1-(1+t)^{1-\lambda})}{1-\lambda}}\tau_{0}^{\frac{3-\gamma}{4}}+\frac{3-\gamma}{4}(K\gamma)^{-\frac{1}{4}}Mt,$

i.e.,

$\displaystyle\rho$
$\displaystyle\geq\left[\tau_{0}^{\frac{3-\gamma}{4}}+\frac{3-\gamma}{4}(K\gamma)^{-\frac{1}{4}}Mt\right]^{\frac{4}{\gamma-3}}$
(3.20)
$\displaystyle\geq\left[C_{0}^{\frac{3-\gamma}{4}}+\frac{3-\gamma}{4}(K\gamma)^{-\frac{1}{4}}Mt\right]^{\frac{4}{\gamma-3}}$
$\displaystyle\geq C(1+t)^{\frac{4}{\gamma-3}},$

where we have used $1<\gamma<3$, (2.5) and

$C:=\left(\max\\{C_{0}^{\frac{3-\gamma}{4}},\frac{3-\gamma}{4}(K\gamma)^{-\frac{1}{4}}M\\}\right)^{\frac{4}{\gamma-3}}.$
(3.21)

Step 3. For $\lambda=1$, we denote $A_{1}(t)=(1+t)^{\frac{\alpha}{2}}$ and

$G_{3}(t,x)=A_{1}(t)G_{1}(t,x),\quad H_{3}(t,x)=A_{1}(t)H_{1}(t,x),$

then it follows from (3.13) that $G_{3}(0,x)\leq 0,H_{3}(0,x)\leq 0$ and

$\displaystyle
G_{3}^{\prime}=A_{1}(t)G_{1}^{\prime}+\frac{\alpha}{2(1+t)}A_{1}(t)G_{1}\leq\frac{\alpha}{2(1+t)}A_{1}(t)H_{1}=\frac{\alpha}{2(1+t)}H_{3},$
(3.22) $\displaystyle
H_{3}^{\backprime}=A_{1}(t)H_{1}^{\backprime}+\frac{\alpha}{2(1+t)}A_{1}(t)H_{1}\leq\frac{\alpha}{2(1+t)}A_{1}(t)G_{1}=\frac{\alpha}{2(1+t)}G_{3}.$

Hence, by a similar argument in Step 2, we obtain

$G(t,x)\leq M,\quad H(t,x)\leq M,\quad\text{for all
}(t,x)\in\mathbb{R}_{+}\times\mathbb{R},$

which implies that

$G(t,x)+H(t,x)=c^{\frac{1}{2}}(w+z)_{x}+\frac{2\alpha}{(1+t)}h(\tau)\leq 2M.$
(3.23)

Using $\eqref{1.1}_{1}$, (2.2) and (3.9), then (3.23) becomes

$(\tau^{\frac{3-\gamma}{4}})_{t}+{\alpha\over{(1+t)}}\tau^{\frac{3-\gamma}{4}}\leq\frac{3-\gamma}{4}(K\gamma)^{-\frac{1}{4}}M,$
(3.24)

and a direct calculation shows that

$\big{(}(1+t)^{\alpha}\tau^{\frac{3-\gamma}{4}}\big{)}_{t}\leq\frac{3-\gamma}{4}(K\gamma)^{-\frac{1}{4}}M(1+t)^{\alpha}.$
(3.25)

Integrating (3.25) from $0$ to $t$, we have

$\displaystyle(1+t)^{\alpha}\tau^{\frac{3-\gamma}{4}}-\tau_{0}^{\frac{3-\gamma}{4}}$
$\displaystyle\leq\frac{3-\gamma}{4}(K\gamma)^{-\frac{1}{4}}M\int_{0}^{t}(1+s)^{\alpha}ds$
$\displaystyle\leq\frac{3-\gamma}{4}(K\gamma)^{-\frac{1}{4}}M(1+t)^{\alpha}t,$

and therefore

$\displaystyle\rho$
$\displaystyle\geq\left[(1+t)^{-\alpha}\tau_{0}^{\frac{3-\gamma}{4}}+\frac{3-\gamma}{4}(K\gamma)^{-\frac{1}{4}}Mt\right]^{\frac{4}{\gamma-3}}$
(3.26)
$\displaystyle\geq\left[\tau_{0}^{\frac{3-\gamma}{4}}+\frac{3-\gamma}{4}(K\gamma)^{-\frac{1}{4}}Mt\right]^{\frac{4}{\gamma-3}}$
$\displaystyle\geq C(1+t)^{\frac{4}{\gamma-3}},$

where we have used $1<\gamma<3$, (2.5) and (3.21).

Step 4. For $\lambda<0$, we introduce two new variables

$\hat{G}(t,x)=c^{1\over
2}w_{x}+{\alpha-\lambda\over{(1+t)^{\lambda}}}h(\tau),$ (3.27)
$\hat{H}(t,x)=c^{1\over
2}z_{x}+{\alpha-\lambda\over{(1+t)^{\lambda}}}h(\tau),$ (3.28)

where $h(\tau)$ is defined (3.9). By a similar calculation in (3.10)-(3.11),
one has

$\displaystyle\hat{G}^{\prime}$
$\displaystyle=\frac{\alpha-2\lambda}{2(1+t)^{\lambda}}\hat{H}-\frac{\alpha}{2(1+t)^{\lambda}}\hat{G}-c^{1\over
2}c_{w}w_{x}^{2}+\frac{\lambda(\alpha-\lambda)(1-(1+t)^{\lambda-1})}{(1+t)^{2\lambda}}h(\tau),$
(3.29)
$\displaystyle:=\frac{\alpha-2\lambda}{2(1+t)^{\lambda}}\hat{H}-\frac{\alpha}{2(1+t)^{\lambda}}\hat{G}+V_{1},$
$\displaystyle\hat{H}^{\backprime}$
$\displaystyle=-\frac{\alpha}{2(1+t)^{\lambda}}\hat{H}+\frac{\alpha-2\lambda}{2(1+t)^{\lambda}}\hat{G}-c^{1\over
2}c_{w}z_{x}^{2}+\frac{\lambda(\alpha-\lambda)(1-(1+t)^{\lambda-1})}{(1+t)^{2\lambda}}h(\tau)$
$\displaystyle:=-\frac{\alpha}{2(1+t)^{\lambda}}\hat{H}+\frac{\alpha-2\lambda}{2(1+t)^{\lambda}}\hat{G}+V_{2},$

where

$V_{1}:=-c^{1\over
2}c_{w}w_{x}^{2}+\frac{\lambda(\alpha-\lambda)(1-(1+t)^{\lambda-1})}{(1+t)^{2\lambda}}h(\tau)\leq
0,$ $V_{2}:=-c^{1\over
2}c_{w}z_{x}^{2}+\frac{\lambda(\alpha-\lambda)(1-(1+t)^{\lambda-1})}{(1+t)^{2\lambda}}h(\tau)\leq
0,$

for $\lambda<0$.

Let

$\hat{G}_{1}(t,x)=\hat{G}(t,x)-e^{-\frac{\lambda(1+t)^{1-\lambda}}{1-\lambda}}M,\quad\hat{H}_{1}(t,x)=\hat{H}(t,x)-e^{-\frac{\lambda(1+t)^{1-\lambda}}{1-\lambda}}M,$

where $M$, like the constant $M$ defined in (3.12), is a big enough constant
determined by the initial data, such that

$\hat{G}_{1}(0,x)=\hat{G}(0)-e^{-\frac{\lambda}{1-\lambda}}M\leq
0,\quad\hat{H}_{1}(0,x)=\hat{H}(0)-e^{-\frac{\lambda}{1-\lambda}}M\leq 0.$
(3.30)

A direct calculation by using (3.29) shows that

$\displaystyle\hat{G}_{1}^{\prime}$
$\displaystyle=\hat{G}^{\prime}+\lambda(1+t)^{-\lambda}e^{-\frac{\lambda(1+t)^{1-\lambda}}{1-\lambda}}M$
(3.31)
$\displaystyle\leq\frac{\alpha-2\lambda}{2(1+t)^{\lambda}}\left(\hat{H}_{1}+e^{-\frac{\lambda(1+t)^{1-\lambda}}{1-\lambda}}M\right)-\frac{\alpha}{2(1+t)^{\lambda}}\left(\hat{G}_{1}+e^{-\frac{\lambda(1+t)^{1-\lambda}}{1-\lambda}}M\right)+\frac{\lambda
e^{-\frac{\lambda(1+t)^{1-\lambda}}{1-\lambda}}M}{(1+t)^{\lambda}}$
$\displaystyle=-\frac{\alpha}{2(1+t)^{\lambda}}\hat{G}_{1}+\frac{\alpha-2\lambda}{2(1+t)^{\lambda}}\hat{H}_{1}+\left(\frac{\alpha-2\lambda}{2(1+t)^{\lambda}}-\frac{\alpha}{2(1+t)^{\lambda}}+\frac{\lambda}{(1+t)^{\lambda}}\right)e^{-\frac{\lambda(1+t)^{1-\lambda}}{1-\lambda}}M$
$\displaystyle=-\frac{\alpha}{2(1+t)^{\lambda}}\hat{G}_{1}+\frac{\alpha-2\lambda}{2(1+t)^{\lambda}}\hat{H}_{1},$

and similarly

$\displaystyle\hat{H}_{1}^{\backprime}$
$\displaystyle\leq-\frac{\alpha}{2(1+t)^{\lambda}}\hat{H}_{1}+\frac{\alpha-2\lambda}{2(1+t)^{\lambda}}\hat{G}_{1}.$
(3.32)

Since $\lambda<0$, denoting
$A(t)=e^{\frac{\alpha(1+t)^{1-\lambda}}{2(1-\lambda)}}$ and

$\hat{G}_{2}(t,x)=A(t)\hat{G}_{1}(t,x),\quad\hat{H}_{2}(t,x)=A(t)\hat{H}_{1}(t,x),$

then it follows from (3.29)-(3.30), (3.31)-(3.32) that

$\hat{G}_{2}(0,x)\leq 0,\quad\hat{H}_{2}(0,x)\leq 0,$

and

$\hat{G}_{2}^{\prime}\leq\frac{\alpha-2\lambda}{2(1+t)^{\lambda}}\hat{H}_{2},\quad\hat{H}_{2}^{\backprime}\leq\frac{\alpha-2\lambda}{2(1+t)^{\lambda}}\hat{G}_{2}.$

Hence by a similar argument in Step 2, we obtain

$\hat{G}(t,x)+\hat{H}(t,x)=c^{1\over
2}(w+z)_{x}+{2(\alpha-\lambda)\over{(1+t)^{\lambda}}}h(\tau)\leq
2e^{-\frac{\lambda(1+t)^{1-\lambda}}{1-\lambda}}M,$

which together with $\eqref{1.1}_{1}$, yields

$c^{1\over 2}\tau_{t}+{\alpha-\lambda\over{(1+t)^{\lambda}}}h(\tau)\leq
e^{-\frac{\lambda(1+t)^{1-\lambda}}{1-\lambda}}M.$ (3.33)

Thus using a similar calculation in (3.16)-(3.20), we have

$\displaystyle\tau^{\frac{3-\gamma}{4}}$
$\displaystyle\leq\left[e^{\frac{\alpha-\lambda}{1-\lambda}}\tau_{0}^{\frac{3-\gamma}{4}}+\frac{3-\gamma}{4}(K\gamma)^{-\frac{1}{4}}Me^{\frac{(\alpha-2\lambda)(1+t)^{1-\lambda}}{1-\lambda}}t\right]e^{-\frac{(\alpha-\lambda)(1+t)^{1-\lambda}}{1-\lambda}}$
$\displaystyle=e^{\frac{(\alpha-\lambda)(1-(1+t)^{1-\lambda})}{1-\lambda}}\tau_{0}^{\frac{3-\gamma}{4}}+\frac{3-\gamma}{4}(K\gamma)^{-\frac{1}{4}}Me^{\frac{-\lambda(1+t)^{1-\lambda}}{1-\lambda}}t,$

i.e.,

$\displaystyle\rho$
$\displaystyle\geq\left[e^{\frac{(\alpha-\lambda)(1-(1+t)^{1-\lambda})}{1-\lambda}}\tau_{0}^{\frac{3-\gamma}{4}}+\frac{3-\gamma}{4}(K\gamma)^{-\frac{1}{4}}Me^{\frac{-\lambda(1+t)^{1-\lambda}}{1-\lambda}}t\right]^{\frac{4}{\gamma-3}}$
(3.34)
$\displaystyle\geq\left[\tau_{0}^{\frac{3-\gamma}{4}}+\frac{3-\gamma}{4}(K\gamma)^{-\frac{1}{4}}Me^{\frac{-\lambda(1+t)^{1-\lambda}}{1-\lambda}}t\right]^{\frac{4}{\gamma-3}}$
$\displaystyle\geq
C\left[1+e^{\frac{-\lambda(1+t)^{1-\lambda}}{1-\lambda}}t\right]^{\frac{4}{\gamma-3}},$

where we have used $1<\gamma<3$, (2.5) and (3.21).

Step 5. Now we are going to prove the blow-up of
$\|(\tau_{x},u_{x})\|_{L^{\infty}}$. Since the case
$\lambda\geq\frac{\alpha(\gamma-1)}{\gamma-3}$ is solved in Theorem 2.1, then
we only need to prove the rest case
$\lambda<\frac{\alpha(\gamma-1)}{\gamma-3}$. Noticing that
$\lambda<\frac{\alpha(\gamma-1)}{\gamma-3}<0$ for $1<\gamma<3$, it follows
from (2.26) that there exists $t_{0}>0$ such that

$(1+t_{0})^{\lambda-1}=\frac{\alpha(\gamma-1)}{\lambda(\gamma-3)}$

and

$a_{0}\geq 0,\ \ \forall\ 0<t\leq t_{0};\ \ a_{0}<0,\ \ \forall\ t>t_{0}.$

Using (2.1), (2.25)-(2.26) and Lemma 2.10, there exist two constants
$\tilde{K}_{3}$ and $\tilde{K}_{4}$, depending only on
$C_{0},\gamma,K,\lambda$ and $\alpha$, such that

$0\leq a_{0}\leq\tilde{K}_{3},\quad a_{2}\geq\tilde{K}_{4}>0,\quad\forall\
0<t\leq t_{0}.$ (3.35)

When $t\geq t_{0}$, using (3.34) and a direct calculation shows that

$\displaystyle\int_{t_{0}}^{\infty}a_{2}(x,s)ds$ (3.36) $\displaystyle=$
$\displaystyle\frac{K_{c}(\gamma+1)}{2(\gamma-1)}\left(\frac{2\sqrt{K\gamma}}{\gamma-1}\right)^{-\frac{\gamma-3}{2(\gamma-1)}}\int_{t_{0}}^{+\infty}\rho^{\frac{3-\gamma}{4}}e^{-\frac{\alpha(3\gamma-1)}{2(\gamma-3)(1-\lambda)}(1+s)^{1-\lambda}}ds$
$\displaystyle\geq$
$\displaystyle\frac{K_{c}(\gamma+1)}{2(\gamma-1)}\left(\frac{2\sqrt{K\gamma}}{\gamma-1}\right)^{-\frac{\gamma-3}{2(\gamma-1)}}C^{\frac{3-\gamma}{4}}\int_{t_{0}}^{+\infty}\left[1+e^{\frac{-\lambda(1+s)^{1-\lambda}}{1-\lambda}}s\right]^{-1}e^{-\frac{\alpha(3\gamma-1)}{2(\gamma-3)(1-\lambda)}(1+s)^{1-\lambda}}ds$
$\displaystyle\geq$
$\displaystyle{\left\\{\begin{matrix}\frac{K_{c}(\gamma+1)}{2(\gamma-1)}\left(\frac{2\sqrt{K\gamma}}{\gamma-1}\right)^{-\frac{\gamma-3}{2(\gamma-1)}}C^{\frac{3-\gamma}{4}}\displaystyle\int_{t_{0}}^{+\infty}\left[e^{\frac{\lambda(1+s)^{1-\lambda}}{1-\lambda}}+s\right]^{-1}e^{\left(\frac{\alpha(3\gamma-1)}{2(3-\gamma)}+\lambda\right)\frac{(1+s)^{1-\lambda}}{1-\lambda}}ds\\\
\text{for
}\frac{\alpha(3\gamma-1)}{2(3-\gamma)}+\lambda<0\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\\\
\frac{K_{c}(\gamma+1)}{2(\gamma-1)}\left(\frac{2\sqrt{K\gamma}}{\gamma-1}\right)^{-\frac{\gamma-3}{2(\gamma-1)}}C^{\frac{3-\gamma}{4}}e^{\left(\frac{\alpha(3\gamma-1)}{2(3-\gamma)(1-\lambda)}+\frac{\lambda}{1-\lambda}\right)}\displaystyle\int_{t_{0}}^{+\infty}\left[e^{\frac{\lambda(1+s)^{1-\lambda}}{1-\lambda}}+s\right]^{-1}ds\\\
\text{for }\frac{\alpha(3\gamma-1)}{2(3-\gamma)}+\lambda\geq
0\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\end{matrix}\right.}$
$\displaystyle\geq$
$\displaystyle{\left\\{\begin{matrix}\tilde{N}_{1}>0,&\quad\text{for
}\frac{\alpha(3\gamma-1)}{2(3-\gamma)}+\lambda<0\\\ +\infty,&\quad\text{for
}\frac{\alpha(3\gamma-1)}{2(3-\gamma)}+\lambda\geq 0.\end{matrix}\right.}$

Assume there exists one point $x_{0}$ such that

$y(x_{0},0)<-\max\big{\\{}\frac{1}{\tilde{N}_{1}},\sqrt{\frac{\tilde{K}_{3}}{\tilde{K}_{4}}}\big{\\}}:=-N_{1},$

which means (2.14) holds with

$K_{3}:=\frac{\alpha(\gamma-1)}{K_{c}(\gamma-3)},\quad
K_{4}:=N_{1}e^{-\frac{\alpha(3\gamma-1)}{2(\gamma-3)(1-\lambda)}}.$

Noticing from (2.24) and (3.35), along the forward characteristic starting
form $x_{0}$, we have

$y^{\prime}<0\quad\text{ and }\quad y(x_{+}(t),t)<-N_{1},\quad\forall\ 0<t\leq
t_{0}.$ (3.37)

If $y(x_{+}(t),t)$ does not blow up before $t_{0}$, then

$y^{\prime}=a_{0}-a_{2}y^{2}<-a_{2}y^{2},\quad\forall\ t\geq t_{0},$

which implies that

$0>y^{-1}(x_{+}(t),t)\geq
y^{-1}(x_{+}(t_{0}),t_{0})+\int_{t_{0}}^{t}a_{2}(x_{+}(s),s)ds.$ (3.38)

Since it follows from (3.37) that
$-\frac{1}{N_{1}}<y^{-1}(x_{+}(t_{0}),t_{0})<0$, then using (3.36) and (3.38),
we obtain that $y(x_{+}(t),t)$ must blow up in finite time. This holds for $q$
if (2.15) holds.

On the other hand, it follows from (2.1) and (3.34) that
$\phi^{\frac{\gamma-3}{2(\gamma-1)}}$ will remain bounded for any finite time,
then recalling (2.22)-(2.23), we obtain that $y(x_{+}(t),t)$ or
$q(x_{-}(t),t)$ blows up in finite time means
$\|(\tau_{x},u_{x})\|_{L^{\infty}}$ blows up in finite time. Therefore the
proof of Theorem 2.3 is completed.$\hfill\square$

###  3.2 Proof of Theorem 2.4

We divide the proof into two cases: $\lambda>1$ and $\lambda<1$.

Case 1. For $\gamma>3$ and $\lambda>1$, it follows from Lemma 2.10 and
(2.25)-(2.26) that there exist two positive constants
$\tilde{K}_{5},\tilde{K}_{6}$, which depends only on $C_{0},K,\gamma,\lambda$
and $\alpha$, such that

$|a_{0}|\leq\tilde{K}_{5},\quad a_{2}\geq\tilde{K}_{6},\quad\text{for all
}t\geq 0.$

Assume there exists one point $x_{0}$ such that

$y(x_{0},0)<-\sqrt{\frac{\tilde{K}_{5}}{\tilde{K}_{6}}}:=-N_{2},$ (3.39)

which means (2.16) holds with

$K_{5}:=\frac{\alpha(\gamma-1)}{K_{c}(\gamma-3)},\quad
K_{6}:=N_{2}e^{-\frac{\alpha(3\gamma-1)}{2(\gamma-3)(1-\lambda)}}.$

Noticing from (3.39), there exists a $\varepsilon>0$ such that

$y(x_{0},0)<-(1+\varepsilon)N_{2}.$

Recalling (2.24), along the forward characteristic lines $x_{+}(t)$ starting
from $x_{0}$, we have

$y^{\prime}(x_{+}(t),t)<0\quad\text{ and }\quad
y(x_{+}(t),t)<-(1+\varepsilon)N_{2},\quad\forall\ t\geq 0,$

which together with the definition of $N_{2}$ implies that

$a_{0}-a_{2}\frac{y^{2}(x_{+}(t),t)}{(1+\varepsilon)^{2}}<0,\quad\forall\
t\geq 0.$

Hence, using (2.24) and $a_{2}>0$, we have

$y^{\prime}(x_{+}(t),t)<\left(-1+\frac{1}{(1+\varepsilon)^{2}}\right)a_{2}y^{2},\quad\forall\
t\geq 0.$

Integrating both sides of the above inequality with respect to $t$ along the
forward characteristic $x_{+}(t)$ starting form $x_{0}$, we obtain

$0>y^{-1}(x_{+}(t),t)\geq
y^{-1}(x_{0},0)+\left(1-\frac{1}{(1+\varepsilon)^{2}}\right)\int_{0}^{t}a_{2}(x_{+}(s),s)ds.$
(3.40)

Noticing that

$\int_{0}^{+\infty}a_{2}(x_{+}(s),s)ds\geq\int_{0}^{+\infty}\tilde{K}_{6}ds=+\infty,$

which shows that the right hand side of (3.40) will be positive in finite
time, then $y(x_{+}(t),t)$ must blow up in finite time. Similarly, if (2.17)
holds, we will obtain that $q(x_{-}(t),t)$ blows up in finite time. Moreover,
it follows from (2.1) and Lemma 2.10 that
$\phi^{\frac{\gamma+1}{2(\gamma-1)}}$ and
$\phi^{\frac{\gamma-3}{2(\gamma-1)}}$ remain bounded for all time when
$\gamma>3$. Hence $y(x_{+}(t),t)$ or $q(x_{-}(t),t)$ blows up means
$\|(\tau_{x},u_{x})\|_{L^{\infty}}$ blows up.

Case 2. For $\gamma>3$ and $\lambda<1$, $a_{0}$ may not be uniformly bounded
for all $t>0$ and $a_{2}$ may not have a positive lower bound. However, if
$\lambda\leq\min\\{\frac{\alpha(\gamma-1)}{\gamma-3},1\\}$ or $\alpha=0$, we
will have $a_{0}\leq 0$ for all $t>0$ and hence using (2.24), we have

$y^{\prime}=a_{0}-a_{2}y^{2}<-a_{2}y^{2}.$

This is the case in Theorem 2.2 and is solved in [1, Theorem 3.1].

If $\alpha>0$ and $0<\frac{\alpha(\gamma-1)}{\gamma-3}<\lambda<1$, there
exists $t_{0}>0$ satisfying

$(1+t_{0})^{\lambda-1}=\frac{\alpha(\gamma-1)}{\lambda(\gamma-3)}$

such that

$a_{0}\geq 0,\ \ \forall\ 0<t\leq t_{0}\ \ \text{and }\ a_{0}<0,\ \ \forall\
t>t_{0}.$

When $0<t\leq t_{0}$, using Lemma 2.10, there exists two constants
$\hat{K}_{5}$ and $\hat{K}_{6}$, which are dependent only on
$C_{0},\gamma,K,\lambda$ and $\alpha$, such that

$0\leq a_{0}\leq\hat{K}_{5},\quad a_{2}\geq\hat{K}_{6},\quad\forall\ 0<t\leq
t_{0}.$

When $t\geq t_{0}$, using Lemma 2.10, a direct calculation shows that

$\displaystyle\int_{t_{0}}^{+\infty}a_{2}(x(s),s)ds$ (3.41)
$\displaystyle=\int_{t_{0}}^{+\infty}\frac{K_{c}(\gamma+1)}{2(\gamma-1)}\phi^{-\frac{\gamma-3}{2(\gamma-1)}}e^{-\frac{\alpha(3\gamma-1)}{2(\gamma-3)(1-\lambda)}(1+s)^{1-\lambda}}ds$
$\displaystyle=\int_{t_{0}}^{+\infty}\frac{K_{c}(\gamma+1)}{2(\gamma-1)}\left(\frac{2\sqrt{K\gamma}}{\gamma-1}\tau^{-\frac{\gamma-1}{2}}\right)^{-\frac{\gamma-3}{2(\gamma-1)}}e^{-\frac{\alpha(3\gamma-1)}{2(\gamma-3)(1-\lambda)}(1+s)^{1-\lambda}}ds$
$\displaystyle\geq\frac{K_{c}(\gamma+1)}{2(\gamma-1)}\left(\frac{2\sqrt{K\gamma}}{\gamma-1}\right)^{-\frac{\gamma-3}{2(\gamma-1)}}\tilde{C}_{0}^{-\frac{\gamma-3}{4}}\int_{t_{0}}^{+\infty}e^{-\frac{\alpha(3\gamma-1)}{2(\gamma-3)(1-\lambda)}(1+s)^{1-\lambda}}ds$
$\displaystyle\geq\frac{K_{c}(\gamma+1)}{2(\gamma-1)}\left(\frac{2\sqrt{K\gamma}}{\gamma-1}\right)^{-\frac{\gamma-3}{2(\gamma-1)}}\tilde{C}_{0}^{-\frac{\gamma-3}{4}}\int_{t_{0}}^{+\infty}(1+s)^{-\lambda}e^{-\frac{\alpha(3\gamma-1)}{2(\gamma-3)(1-\lambda)}(1+s)^{1-\lambda}}ds$
$\displaystyle\geq\frac{K_{c}(\gamma+1)}{2(\gamma-1)}\left(\frac{2\sqrt{K\gamma}}{\gamma-1}\right)^{-\frac{\gamma-3}{2(\gamma-1)}}\tilde{C}_{0}^{-\frac{\gamma-3}{4}}\frac{2(\gamma-3)}{\alpha(3\gamma-1)}e^{-\frac{\alpha(3\gamma-1)}{2(\gamma-3)(1-\lambda)}}:=\tilde{N}_{3}.$

If there exists a point $x_{0}$ such that

$y(x_{0},0)<-\max\\{\frac{1}{\tilde{N}_{3}},\sqrt{\frac{\hat{K}_{5}}{\hat{K}_{6}}}\\}:=-N_{3},$

which is equivalent to (2.16) with

$K_{5}:=\frac{\alpha(\gamma-1)}{K_{c}(\gamma-3)},\quad
K_{6}:=N_{3}e^{-\frac{\alpha(3\gamma-1)}{2(\gamma-3)(1-\lambda)}}.$ (3.42)

Then following the same discussion in (3.37)-(3.38), we obtain $y(x_{+}(t),t)$
must blow up in finite time. This holds for $q$ if (2.17) holds. Hence
$\|(\tau_{x},u_{x})\|_{L^{\infty}}$ must blow up in finite time. This
completes the proof of Theorem 2.4.

$\hfill\square$

######  Remark 3.1.

When $\lambda=1$, the authors in [1] have also proved similar blow up results
for $\gamma>3,$ $\alpha\geq\frac{\gamma-3}{\gamma-1}$ and
$1<\gamma<3,\alpha\geq 0$. In fact, when $\gamma>3$ and $\lambda=1$, following
the same idea in the proofs of Theorem 2.3 and Theorem 2.4, the similar blow
up phenomenon can be proved for
$0\leq\alpha\leq\frac{2(\gamma-3)}{3\gamma-1}$. However, for the rest case
$\alpha\in\left(\frac{2(\gamma-3)}{3\gamma-1},\frac{\gamma-3}{\gamma-1}\right)$,
there is no blow up results yet now. The main reason is that we can not get
the upper bound of $|{{b_{0}}\over{b_{1}}}|$ in (4.3) and (4.4) of [1].

###  3.3 Proof of Theorem 2.6

We divide the proof into four steps. In steps 1-3, we first give the proof of
the lower bound on density, that is, the proof of (2.18). In step 4, we prove
the blow-up of $\|(\tau_{x},u_{x})\|_{L^{\infty}}$.

Step 1. For $\gamma=3$, the previous notations given in (2.1)-(2.2) becomes

$\displaystyle\phi:=\int^{\infty}_{\tau}cd\tau=\sqrt{3K}\tau^{-1}>0,$ (3.43)

and

$\displaystyle c:=\sqrt{-p_{\tau}}=\sqrt{3K}\tau^{-2}:=K_{c}\phi^{2}.$ (3.44)

where $K_{c}=1/\sqrt{3K}$.

We introduce two gradient variables

$\displaystyle y_{1}:=\phi{A}-\frac{\alpha}{2K_{c}(1+t)^{\lambda}}\ln\phi,$
(3.45)

and

$\displaystyle q_{1}:=\phi{B}-\frac{\alpha}{2K_{c}(1+t)^{\lambda}}\ln\phi,$
(3.46)

where $A$ and $B$ are given in Lemma 2.11 and satisfy

$\displaystyle{A}^{\prime}=-\frac{\alpha}{2(1+t)^{\lambda}}({A}+{B})+K_{d}({A}{B}-{A}^{2}),$
(3.47)
$\displaystyle{B}^{\backprime}=-\frac{\alpha}{2(1+t)^{\lambda}}({A}+{B})+K_{d}({A}{B}-{B}^{2}).$

For $\gamma=3$, $K_{d}=K_{c}\phi$ from Lemma 2.11.

######  Lemma 3.2.

For $C^{1}$ solutions of (1.1)-(1.2), $y_{1}$ and $q_{1}$ satisfy the
following Riccati equations:

$\displaystyle y_{1}^{\prime}=-a_{2}y_{1}^{2}-a_{1}y_{1}+a_{0},$ (3.48)
$\displaystyle q_{1}^{\backprime}=-a_{2}q_{1}^{2}-a_{1}q_{1}+a_{0},$ (3.49)

where

$\displaystyle a_{2}=K_{c}>0,$ (3.50) $\displaystyle
a_{1}=\frac{\alpha(1+2\ln\phi)}{2(1+t)^{\lambda}},$ $\displaystyle
a_{0}=\alpha\ln\phi\frac{2\lambda(1+t)^{\lambda-1}-\alpha(1+\ln\phi)}{4K_{c}(1+t)^{2\lambda}}.$

Proof of Lemma 3.2. It follows from (3.43), (1.1) and (2.3) that

$\displaystyle\phi^{\prime}$
$\displaystyle=\phi_{t}+c\phi_{x}=\left(\sqrt{3K}\tau^{-1}\right)_{t}+c\phi_{x}$
$\displaystyle=-\sqrt{3K}\tau^{-2}\tau_{t}+c\phi_{x}=-cu_{x}+c\phi_{x}$
$\displaystyle=-c({z}+\phi)_{x}+c\phi_{x}=-c{z}_{x}-c\phi_{x}+c\phi_{x}$
$\displaystyle=-c{B},$

which implies that

$\displaystyle{B}=-\frac{1}{c}\phi^{\prime}.$ (3.51)

Substituting (3.51) into (3.47), we get

$\displaystyle{A}^{\prime}=-\frac{\alpha}{2(1+t)^{\lambda}}({A}-\frac{1}{c}\phi^{\prime})+K_{d}(-\frac{1}{c}\phi^{\prime}{A}-{A}^{2}).$
(3.52)

Multiplying (3.52) by $\phi$ to see

$\displaystyle{A}^{\prime}\phi-\frac{\alpha}{2c(1+t)^{\lambda}}\phi^{\prime}\phi+\frac{K_{d}}{c}\phi^{\prime}{A}\phi=-\frac{\alpha}{2(1+t)^{\lambda}}{A}\phi-
K_{d}{A}^{2}\phi.$ (3.53)

For the left hand side of (3.53), a direct calculation shows that

$\displaystyle{A}^{\prime}\phi-\frac{\alpha}{2c(1+t)^{\lambda}}\phi^{\prime}\phi+\frac{K_{d}}{c}\phi^{\prime}{A}\phi$
(3.54) $\displaystyle=$
$\displaystyle{A}^{\prime}\phi-\frac{\alpha}{2(1+t)^{\lambda}K_{c}\phi^{2}}\phi^{\prime}\phi+\frac{K_{c}\phi}{K_{c}\phi^{2}}\phi^{\prime}{A}\phi$
$\displaystyle=$
$\displaystyle{A}^{\prime}\phi-\frac{\alpha}{2K_{c}(1+t)^{\lambda}}\phi^{\prime}\phi^{-1}+\phi^{\prime}{A}$
$\displaystyle=$
$\displaystyle\left({A}\phi-\frac{\alpha}{2K_{c}(1+t)^{\lambda}}\ln\phi\right)^{\prime}-\frac{\lambda\alpha\ln\phi}{2K_{c}}(1+t)^{-\lambda-1}.$

Recalling (3.45), we have

$\displaystyle{A}=y\phi^{-1}+\frac{\alpha\ln\phi}{2K_{c}(1+t)^{\lambda}}\phi^{-1}.$
(3.55)

Then the right hand side of (3.53) becomes

$\displaystyle-\frac{\alpha}{2(1+t)^{\lambda}}{A}\phi-K_{d}{A}^{2}\phi$ (3.56)
$\displaystyle=$
$\displaystyle-\frac{\alpha}{2(1+t)^{\lambda}}\phi\left(y_{1}\phi^{-1}+\frac{\alpha\ln\phi}{2K_{c}(1+t)^{\lambda}}\phi^{-1}\right)-K_{c}\phi^{2}\left(y_{1}\phi^{-1}+\frac{\alpha\ln\phi}{2K_{c}(1+t)^{\lambda}}\phi^{-1}\right)^{2}$
$\displaystyle=$
$\displaystyle-\frac{\alpha}{2(1+t)^{\lambda}}y_{1}-\frac{\alpha^{2}\ln\phi}{4K_{c}(1+t)^{2\lambda}}-K_{c}y_{1}^{2}-\frac{\alpha^{2}(\ln\phi)^{2}}{4K_{c}(1+t)^{2\lambda}}-\frac{\alpha\ln\phi}{(1+t)^{\lambda}}y$
$\displaystyle=$ $\displaystyle-
K_{c}y_{1}^{2}-\frac{\alpha+2\alpha\ln\phi}{2(1+t)^{\lambda}}y_{1}-\frac{\alpha^{2}\ln\phi+\alpha^{2}(\ln\phi)^{2}}{4K_{c}(1+t)^{2\lambda}}.$

Substituting (3.54) and (3.56) into (3.53), we obtain

$\displaystyle y_{1}^{\prime}=$
$\displaystyle\frac{\lambda\alpha\ln\phi}{2K_{c}}(1+t)^{-\lambda-1}-\frac{\alpha}{2(1+t)^{\lambda}}{A}\phi-
K_{d}{A}^{2}\phi$ (3.57) $\displaystyle=$ $\displaystyle-
K_{c}y_{1}^{2}-\frac{\alpha+2\alpha\ln\phi}{2(1+t)^{\lambda}}y_{1}-\frac{\alpha^{2}\ln\phi+\alpha^{2}(\ln\phi)^{2}}{4K_{c}(1+t)^{2\lambda}}+\frac{\lambda\alpha\ln\phi}{2K_{c}}(1+t)^{-\lambda-1}$
$\displaystyle=$ $\displaystyle-
K_{c}y_{1}^{2}-\frac{\alpha+2\alpha\ln\phi}{2(1+t)^{\lambda}}y_{1}+\frac{2\lambda\alpha\ln\phi(1+t)^{\lambda-1}-\alpha^{2}\ln\phi(1+\ln\phi)}{4K_{c}(1+t)^{2\lambda}}$
$\displaystyle:=$ $\displaystyle-a_{2}y_{1}^{2}-a_{1}y_{1}+a_{0},$

where $a_{0},a_{1}$ and $a_{2}$ are defined in (3.50). This completes the
proof of (3.48).

For (3.49), noticing that

$\displaystyle\phi^{\backprime}$
$\displaystyle=\phi_{t}-c\phi_{x}=\left(\sqrt{3K}\tau^{-1}\right)_{t}-c\phi_{x}$
$\displaystyle=-\sqrt{3K}\tau^{-2}\tau_{t}-c\phi_{x}=-cu_{x}-c\phi_{x}$
$\displaystyle=-c({w}-\phi)_{x}-c\phi_{x}=-c{w}_{x}+c\phi_{x}-c\phi_{x}$
$\displaystyle=-c{A},$

then using (3.47), we have

$\displaystyle{B}^{\backprime}=-\frac{\alpha}{2(1+t)^{\lambda}}({B}-\frac{1}{c}\phi^{\backprime})+K_{d}(-\frac{1}{c}\phi^{\backprime}{B}-{B}^{2}).$

Thus following a similar calculation in (3.53)-(3.57), we obtain

$\displaystyle q_{1}^{\backprime}=$
$\displaystyle\frac{\lambda\alpha\ln\phi}{2K_{c}}(1+t)^{-\lambda-1}-\frac{\alpha}{2(1+t)^{\lambda}}{B}\phi-
K_{d}{B}^{2}\phi$ (3.58) $\displaystyle=$ $\displaystyle-
K_{c}q_{1}^{2}-\frac{\alpha+2\alpha\ln\phi}{2(1+t)^{\lambda}}q_{1}-\frac{\alpha^{2}\ln\phi+\alpha^{2}(\ln\phi)^{2}}{4K_{c}(1+t)^{2\lambda}}+\frac{\lambda\alpha\ln\phi}{2K_{c}}(1+t)^{-\lambda-1}$
$\displaystyle=$ $\displaystyle-
K_{c}q_{1}^{2}-\frac{\alpha+2\alpha\ln\phi}{2(1+t)^{\lambda}}q_{1}+\frac{2\lambda\alpha\ln\phi(1+t)^{\lambda-1}-\alpha^{2}\ln\phi(1+\ln\phi)}{4K_{c}(1+t)^{2\lambda}}$
$\displaystyle:=$ $\displaystyle-a_{2}q_{1}^{2}-a_{1}q_{1}+a_{0},$

where $a_{2},$ $a_{1},$ $a_{0}$ are defined in (3.50). Hence the proof of
Lemma 3.2 is completed.

Step 2. For $\gamma=3$ and $\lambda\geq 0$, we define

$h_{1}(\tau)=\int_{\tilde{C}_{0}^{-1}}^{\tau}c(s)^{\frac{1}{2}}ds,$

where $\tilde{C}_{0}$ is the constant given in Lemma 2.10. It follows from
(3.44) that $h_{1}(\tau)=(3K)^{\frac{1}{4}}\ln(\tilde{C}_{0}\tau)$. By
replacing $h(\tau)$ with $h_{1}(\tau)$ and using the same calculations in the
Steps 2-3 in subsection 3.1, we obtain (3.16), that is,

$(3K)^{1\over
4}\tau^{-1}\tau_{t}+{\alpha\over{(1+t)^{\lambda}}}(3K)^{\frac{1}{4}}\ln(\tilde{C}_{0}\tau)\leq
M,$

which implies that

$(\ln(\tilde{C}_{0}\tau))_{t}+{\alpha\over{(1+t)^{\lambda}}}\ln(\tilde{C}_{0}\tau)\leq(3K)^{-\frac{1}{4}}M.$

Thus for $\lambda\neq 1$,

$\left(e^{\frac{\alpha(1+t)^{1-\lambda}}{1-\lambda}}\ln(\tilde{C}_{0}\tau)\right)_{t}\leq(3K)^{-\frac{1}{4}}Me^{\frac{\alpha(1+t)^{1-\lambda}}{1-\lambda}},$
(3.59)

and for $\lambda=1$,

$\left((1+t)^{\alpha}\ln(\tilde{C}_{0}\tau)\right)_{t}\leq(3K)^{-\frac{1}{4}}M(1+t)^{\alpha},$
(3.60)

where $M$ is the found given in (3.12).

Integrating (3.59)-(3.60) from $0$ to $t$ respectively to see

$\displaystyle
e^{\frac{\alpha(1+t)^{1-\lambda}}{1-\lambda}}\ln(\tilde{C}_{0}\tau)-e^{\frac{\alpha}{1-\lambda}}\ln(\tilde{C}_{0}\tau_{0})$
$\displaystyle\leq(3K)^{-\frac{1}{4}}M\int_{0}^{t}e^{\frac{\alpha(1+s)^{1-\lambda}}{1-\lambda}}ds$
$\displaystyle\leq(3K)^{-\frac{1}{4}}Me^{\frac{\alpha(1+t)^{1-\lambda}}{1-\lambda}}t,\quad\quad\quad\text{for
}\lambda\neq 1,$

and

$\displaystyle(1+t)^{\alpha}\ln(\tilde{C}_{0}\tau)-\ln(\tilde{C}_{0}\tau_{0})$
$\displaystyle\leq(3K)^{-\frac{1}{4}}M\int_{0}^{t}(1+s)^{\alpha}ds$
$\displaystyle\leq(3K)^{-\frac{1}{4}}M(1+t)^{\alpha}t,\qquad\qquad\qquad\quad\text{for
}\lambda=1.$

Therefore

$\displaystyle\ln(\tilde{C}_{0}\tau)$
$\displaystyle\leq\left[e^{\frac{\alpha}{1-\lambda}}\ln(\tilde{C}_{0}\tau_{0})+(3K)^{-\frac{1}{4}}Me^{\frac{\alpha(1+t)^{1-\lambda}}{1-\lambda}}t\right]e^{-\frac{\alpha(1+t)^{1-\lambda}}{1-\lambda}}$
$\displaystyle=e^{\frac{\alpha(1-(1+t)^{1-\lambda})}{1-\lambda}}\ln(\tilde{C}_{0}\tau_{0})+(3K)^{-\frac{1}{4}}Mt$
$\displaystyle\leq\ln(\tilde{C}_{0}\tau_{0})+(3K)^{-\frac{1}{4}}Mt,\quad\quad\quad\text{for
}\lambda\neq 1,$

and

$\displaystyle\ln(\tilde{C}_{0}\tau)$
$\displaystyle\leq\left[\ln(\tilde{C}_{0}\tau_{0})+(3K)^{-\frac{1}{4}}M(1+t)^{\alpha}t\right](1+t)^{-\alpha}$
$\displaystyle=(1+t)^{-\alpha}\ln(\tilde{C}_{0}\tau_{0})+(3K)^{-\frac{1}{4}}Mt$
$\displaystyle\leq\ln(\tilde{C}_{0}\tau_{0})+(3K)^{-\frac{1}{4}}Mt,\qquad\quad\text{for
}\lambda=1.$

Hence

$\displaystyle\tau\leq\tau_{0}e^{(3K)^{-\frac{1}{4}}Mt},\quad\text{for
}\lambda\geq 0,$

i.e.,

$\displaystyle\rho\geq$
$\displaystyle\tau_{0}^{-1}e^{-(3K)^{-\frac{1}{4}}Mt}\geq
C_{0}^{-1}e^{-\hat{M}t},\quad\text{for }\lambda\geq 0,$ (3.61)

where we have used (2.5) and $\hat{M}:=(3K)^{-\frac{1}{4}}M$.

Step 3. For $\gamma=3$, $\lambda<0$, By replacing $h(\tau)$ with $h_{1}(\tau)$
and using the same calculations in the Step 4 in subsection 3.1, we can obtain
(3.33), that is,

$(3K)^{1\over
4}\tau^{-1}\tau_{t}+{\alpha\over{(1+t)^{\lambda}}}(3K)^{\frac{1}{4}}\ln(\tilde{C}_{0}\tau)\leq
e^{-\frac{\lambda(1+t)^{1-\lambda}}{1-\lambda}}M,$

which implies that

$(\ln(\tilde{C}_{0}\tau))_{t}+{\alpha\over{(1+t)^{\lambda}}}\ln(\tilde{C}_{0}\tau)\leq(3K)^{-\frac{1}{4}}e^{-\frac{\lambda(1+t)^{1-\lambda}}{1-\lambda}}M.$

Thus

$\displaystyle\left(e^{\frac{\alpha(1+t)^{1-\lambda}}{1-\lambda}}\ln(\tilde{C}_{0}\tau)\right)_{t}\leq\left((3K)^{-\frac{1}{4}}e^{-\frac{\lambda(1+t)^{1-\lambda}}{1-\lambda}}M\right)e^{\frac{\alpha(1+t)^{1-\lambda}}{1-\lambda}}.$
(3.62)

Integrating (3.62) from $0$ to $t$ to see

$\displaystyle\ \ \ \
e^{\frac{\alpha(1+t)^{1-\lambda}}{1-\lambda}}\ln(\tilde{C}_{0}\tau)-e^{\frac{\alpha}{1-\lambda}}\ln(\tilde{C}_{0}\tau_{0})$
$\displaystyle\leq(3K)^{-\frac{1}{4}}M\int_{0}^{t}e^{\frac{(\alpha-\lambda)(1+s)^{1-\lambda}}{1-\lambda}}ds$
$\displaystyle=(3K)^{-\frac{1}{4}}Me^{\frac{(\alpha-\lambda)(1+t)^{1-\lambda}}{1-\lambda}}t.$

Therefore

$\displaystyle\ln(\tilde{C}_{0}\tau)$
$\displaystyle\leq\left[e^{\frac{\alpha}{1-\lambda}}\ln(\tilde{C}_{0}\tau_{0})+(3K)^{-\frac{1}{4}}Me^{\frac{(\alpha-\lambda)(1+t)^{1-\lambda}}{1-\lambda}}t\right]e^{-\frac{\alpha(1+t)^{1-\lambda}}{1-\lambda}}$
$\displaystyle=e^{\frac{\alpha(1-(1+t)^{1-\lambda})}{1-\lambda}}\ln(\tilde{C}_{0}\tau_{0})+(3K)^{-\frac{1}{4}}Me^{\frac{-\lambda(1+t)^{1-\lambda}}{1-\lambda}}t$
$\displaystyle\leq\ln\tilde{C}_{0}\tau_{0}+(3K)^{-\frac{1}{4}}Me^{\frac{-\lambda(1+t)^{1-\lambda}}{1-\lambda}}t$
$\displaystyle=\ln(\tilde{C}_{0}\tau_{0})+(3K)^{-\frac{1}{4}}Me^{\frac{-\lambda(1+t)^{1-\lambda}}{1-\lambda}}t.$

Hence

$\displaystyle\tau\leq\tau_{0}e^{(3K)^{-\frac{1}{4}}Mte^{\frac{-\lambda(1+t)^{1-\lambda}}{1-\lambda}}},\quad\text{for
}\lambda<0,$

i.e.,

$\displaystyle\rho\geq\tau_{0}^{-1}e^{-(3K)^{-\frac{1}{4}}Mte^{\frac{-\lambda(1+t)^{1-\lambda}}{1-\lambda}}}\geq
C_{0}^{-1}e^{-\hat{M}te^{\frac{-\lambda(1+t)^{1-\lambda}}{1-\lambda}}},\quad\text{for
}\lambda<0,$ (3.63)

where we have used (2.5) and $\hat{M}:=(3K)^{-\frac{1}{4}}M$..

Step 4. For any given $t_{0}>0$, it follows from (3.43), (3.50), (3.61) and
(3.63) that

$\displaystyle|a_{1}|\leq\bigg{[}\frac{\alpha}{2}+\alpha\ln\bar{M}(t_{0})\bigg{]}\max\\{(1+t_{0})^{-\lambda},1\\},$
(3.64)
$\displaystyle|a_{0}|\leq\bigg{[}\frac{|\lambda|\alpha}{2K_{c}}\ln\bar{M}(t_{0})+\frac{\alpha^{2}}{4K_{c}}\ln\bar{M}(t_{0})+\frac{\alpha^{2}}{4K_{c}}(\ln\bar{M}(t_{0}))^{2}\bigg{]}\max\\{(1+t_{0})^{-2\lambda},1\\},$

where for $\lambda\geq 0$,

$\bar{M}(t_{0}):=\max\big{\\{}\sqrt{3K}\tilde{C}_{0},(\sqrt{3K}\tilde{C}_{0})^{-1},\sqrt{3K}e^{-\hat{M}t_{0}},(\sqrt{3K}e^{-\hat{M}t_{0}})^{-1}\big{\\}}\geq
1,$ (3.65)

and for $\lambda<0$,

$\bar{M}(t_{0}):=\max\big{\\{}\sqrt{3K}\tilde{C}_{0},(\sqrt{3K}\tilde{C}_{0})^{-1},\sqrt{3K}e^{-\hat{M}t_{0}e^{\frac{-\lambda}{1-\lambda}(1+t_{0})^{-\lambda}}},(\sqrt{3K}e^{-\hat{M}t_{0}e^{\frac{-\lambda}{1-\lambda}(1+t_{0})^{-\lambda}}})^{-1}\big{\\}}\geq
1.$ (3.66)

It is directly to check that when $\bar{M}(t_{0})$ is non-decreasing with
respect to $t_{0}$.

Using (3.50), (3.64) and (3.65)-(3.66), there exists a function
$\tilde{M}(t_{0})$, which is non-decreasing with respect to $t_{0}$, such that

$\left|\frac{a_{1}}{a_{2}}\right|+\sqrt{\left|\frac{a_{1}}{a_{2}}\right|^{2}+2\left|\frac{a_{0}}{a_{2}}\right|}\leq\tilde{M}(t_{0}),\quad\text{for
any }0\leq t\leq t_{0}.$ (3.67)

If there exists one point $x_{0}$ such that
$y_{1}(x_{0},0)<-\max\big{\\{}\frac{2}{K_{c}t_{0}},\tilde{M}(t_{0})\big{\\}}$,
i.e., (2.19) holds, then using (3.67), one has

$-\frac{1}{2}a_{2}y_{1}(x_{0},0)^{2}-a_{1}(0)y_{1}(x_{0},0)+a_{0}(0)<0,$

and using (3.48), we obtain

$y_{1}^{\prime}|_{t=0}=-\frac{1}{2}a_{2}y_{1}(x_{0},0)^{2}+\left(-\frac{1}{2}a_{2}y_{1}(x_{0},0)^{2}-a_{1}(0)y_{1}(x_{0},0)+a_{0}(0)\right)<-\frac{1}{2}a_{2}y_{1}(x_{0},0)^{2}\leq
0.$

Thus along the forward characteristic line starting from $(x_{0},0)$, we have

$y_{1}(x_{+}(t),t)<y_{1}(x_{0},0)<-\max\big{\\{}\frac{2}{K_{c}t_{0}},\tilde{M}(t_{0})\big{\\}},\quad\text{for
all }0\leq t\leq t_{0}.$ (3.68)

Turning back to (3.48) and using (3.67) and (3.68), it holds

$y_{1}^{\prime}=-\frac{1}{2}a_{2}y_{1}^{2}+\left(-\frac{1}{2}a_{2}y^{2}-a_{1}y+a_{0}\right)<-\frac{1}{2}a_{2}y_{1}^{2},\quad\text{for
all }0\leq t\leq t_{0}.$ (3.69)

Integrating (3.69) with respect to $t$ along the forward characteristic line
starting from $(x_{0},0)$, we get

$0>y_{1}^{-1}(x_{+}(t),t)>y_{1}^{-1}(x_{0},0)+\int_{0}^{t}\frac{1}{2}a_{2}ds=y_{1}^{-1}(x_{0},0)+\frac{1}{2}K_{c}t,\quad\text{for
all }0\leq t\leq t_{0}.$ (3.70)

Since $y_{1}^{-1}(x_{0},0)\in(-\frac{1}{2}K_{c}t_{0},0)$, the right hand side
of (3.70) will goes to zero as $t\rightarrow t_{0}$, which implies that
$y(x_{+}(t),t)$ will blow up before $t_{0}$. By similar arguments, we obtain
that $q_{1}$ will blow up before $t_{0}$ if (2.20) holds.

Finally, using (3.61), (3.63) and (3.49), we see that $y_{1}$ or $q_{1}$ blows
up in finite time means $\|(\tau_{x},u_{x})\|_{L^{\infty}}$ blows up in finite
time. Therefore the proof of Theorem 2.6 is completed. $\hfill\square$

## References

  * [1] Y. Sui, H. Yu, Singularity formation for compressible Euler equations with time-dependent damping, Discrete & Continuous Dynamical Systems-A, 41 (2021) 4921–4941.
  * [2] C. Dafermos, Hyperbolic conservation laws in continuum physics, Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences], Springer-Verlag, Berlin, 325 (2000).
  * [3] G. Stokes, On a difficulty in the theory of sound, Philos. Magazine, 33 (1848) 349–356.
  * [4] P. Lax, Development of singularities of solutions of nonlinear hyperbolic partial differential equations, J. Math. Phys., 5 (1964) 611–614.
  * [5] F. John, Formation of Singularities in One-Dimensional Nonlinear Wave Propagation, Comm. Pure Appl. Math., 27 (1974) 377–405.
  * [6] T. Li, Y. Zhou, D. Kong, Global classical solutions for general quasilinear hyperbolic systems with decay initial data, Nonlinear Anal-Theor., 28 (1997) 1299–1332.
  * [7] T. Liu, Development of singularities in the nonlinear waves for quasilinear hyperbolic partial differential equations, J. Differ. Equations, 33 (1979) 92–111.
  * [8] X. Pan, Blow up of solutions to 1-d Euler equations with time-dependent damping, J. Math. Anal. Appl., 442 (2016) 435–445.
  * [9] K. Cheung, Blowup phenomena for the N-dimensional compressible Euler equations with damping, Z. Angew. Math. Phys., 68 (2017) 27.
  * [10] S. Chen, H. Li, J. Li, M. Mei, K. Zhang, Global and blow-up solutions for compressible Euler equations with time-dependent damping, J. Differ. Equations, 268 (2020) 5035–5077.
  * [11] X. Fang, H. Yu, Uniform boundedness in weak solutions to a specific dissipative system, J. Math. Anal. Appl., 461 (2018) 1153–1164.
  * [12] F. Huang, T. Li, H. Yu, D. Yuan, Large time behavior of entropy solutions to one-dimensional unipolar hydrodynamic model for semiconductor devices, Z. Angew. Math. Phys., 69 (2018) 69.
  * [13] Y. Sui, H. Yu, Vacuum and singularity formation problem for compressible Euler equations with general pressure law and time-dependent damping, Nonlinear Analysis: Real World Applications, 65 (2022) 103472.
  * [14] H. Yu, Asymptotic Behavior of the 3D Compressible Euler Equations with Nonlinear Damping and Slip Boundary Condition, J. App. Math., (2012) 1–16.
  * [15] H. Li, J. Li, M. Mei, K. Zhang, Convergence to nonlinear diffusion waves for solutions of p-system with time-dependent damping, J. Math. Anal. Appl., 456 (2017) 849–871.
  * [16] H. Li, J. Li, M. Mei, K. Zhang, Optimal convergence rate to nonlinear diffusion waves for Euler equations with critical overdamping, Appl. Math. Lett., 113 (2021) 106882.
  * [17] Y. Sugiyama, Singularity formation for the 1D compressible Euler equations with variable damping coefficient, Nonlinear Analysis, 170 (2018) 70–87.
  * [18] G. Chen, R. Pan, S. Zhu, Singularity formation for the compressible Euler equations, SIAM J. Math. Anal., 49 (2017) 2591–2614.
  * [19] G. Chen, Optimal time-dependent lower bound on density for classical solution of 1-D compressible Euler equations, Indiana Univ. Math. J., 66(3) (2017) 725–740.
  * [20] G. Chen, Optimal density lower bound on nonisentropic gas dynamics, J. Differential. Equations., 268 (2020) 4017–4028.


QUBO Algorithm for Eigenvector ProblemsB. Krakoff, S. Mniszewski, C. Negre
ex_supplement

# A QUBO Algorithm to Compute Eigenvectors of Symmetric Matrices ††thanks:
This article was submitted to the arXiv This research was supported in part by
an appointment with the National Science Foundation (NSF) Mathematical
Sciences Graduate Internship (MSGI) Program sponsored by the NSF Division of
Mathematical Sciences. This program is administered by the Oak Ridge Institute
for Science and Education (ORISE) through an interagency agreement between the
U.S. Department of Energy (DOE) and NSF. ORISE is managed for DOE by ORAU. All
opinions expressed in this paper are the author’s and do not necessarily
reflect the policies and views of NSF, ORAU/ORISE, or DOE. We acknowledge LANL
computational facilities. This document has been assigned the following
Unlimited Release Number (LAUR): ’LA-UR-21-21969’ dated Feb 26th, 2021.

Benjamin Krakoff University of Michigan Mathematics Department ().
<EMAIL_ADDRESS>Susan M. Mniszewski Computer, Computational, & Statistical
Sciences Division, Los Alamos National Laboratory<EMAIL_ADDRESS>Christian
F. A. Negre Theoretical Division, Los Alamos National Laboratory ().
<EMAIL_ADDRESS>
###### Abstract

We describe an algorithm to compute the extremal eigenvalues and corresponding
eigenvectors of a symmetric matrix by solving a sequence of Quadratic Binary
Optimization problems. This algorithm is robust across many different classes
of symmetric matrices, can compute the eigenvector/eigenvalue pair to
essentially arbitrary precision, and with minor modifications can also solve
the generalized eigenvalue problem. Performance is analyzed on small random
matrices and selected larger matrices from practical applications.

###### keywords:

QUBO, annealing, eigenvalue problems

##  1 Introduction

The problem of computing eigenvectors and eigenvalues to a desired precision
has many applications in science and mathematics, including web page ranking
[13], planar embeddings [8] and principal component analysis [14], among many
others. The recent development of new computing paradigms has led to the
production of various kinds of _annealers_ , which are specialized hardware
designed to solve Quadratic Binary Optimization Problems (QUBOs). Such
annealers include D-Wave’s quantumn annealers and Fujitsu’s Digital Annealer,
to name a few. This has lead to a corresponding interest in reformulating
computational tasks as QUBOs and solving them using these annealers. This
strategy has been applied to several problems including, graph partitioning
[16], solving polynomial equations [5] and vertex coloring [10]. Here we
compute eigenvectors of symmetric matrices by solving a sequence of QUBOs,
which allow the eigenvectors and eigenvalues to be found to any desired
precision. A mathematically similar approach to this problem is considered in
[15], but accuracy is increased by increasing the size of the associated QUBO.
In contrast, the proposed algorithm can compute eigenvectors to essentially
arbitrary precision without increasing the size of the QUBOs, which can have
as few as twice as many variables as the original eigenvalue problem. The
trade-off for using small QUBOs is that more iterations are required. A nearly
identical approach is considered in Appendix C of [1], although here the
effects of different parameters are more thoroughly studied, and the
presentation gives a very general optimization framework. The performance data
is collected using D-Wave’s Ocean simulated annealing (SA) package.

The paper is organized as follows. The relevant mathematical background for
symmetric matrices and use of QUBO solvers as a descent method is explained in
Section 2.1. The algorithm for computing the eigenvector/eigenvalue pair is
given in Section 2.3. Experimental results with various parameters and
matrices are presented in Section 3, followed by the conclusion in Section 4.

##  2 Methods

###  2.1 Mathematical Background

Let $A$ be a symmetric matrix. A well-known consequence of the spectral
theorem is that the smallest eigenvalue $\lambda$ and corresponding
eigenvector $\mathbf{v}$ are global minima for the Rayleigh quotient
$\mathbf{x}^{t}A\mathbf{x}/\mathbf{x}^{t}\mathbf{x}$

(1) $\displaystyle\lambda=\min_{||\mathbf{x}||=1}\mathbf{x}^{t}A\mathbf{x},\ \
\mathbf{v}=\underset{||\mathbf{x}||=1}{\text{argmin
}}\mathbf{x}^{t}A\mathbf{x}$

The proposed algorithm uses a QUBO formulation of the problem to both obtain a
good initial guess for the global minimum, and to implement an iterative
descent from the initial guess. Similar to classical descent methods such as
Newton Conjugate-Gradient and the BFGS algorithms [9], the algorithm requires
computing, but not inverting, a Hessian matrix at each descent step. We begin
with an overview of QUBOs and how they can be used to approximately solve
certain constrained quadratic optimization problems.
Let $\\{0,1\\}^{m}$ denote the set of binary vectors of length $m$, and let
$Q$ be a symmetric $m\times m$ matrix. The combinatorial optimization problem

$\underset{\mathbf{x}_{b}\in\\{0,1\\}^{m}}{\text{argmin
}}\mathbf{x}_{b}^{t}Q\mathbf{x}_{b}$

is called a _quantum unconstrained binary optimization_ problem, or QUBO, and
it is known to be NP-hard [3]. Interest in casting various problems as QUBOs
has increased due to the development of various kinds of annealers, which are
a class of hardware that use ideas from statistical mechanics to produce
approximate solutions to a QUBO. See for example [2] or [4].
To solve a real-variable optimization problem using a QUBO, we require a
method of approximating each real variable by $b$ binary variables. This
number $b$ will be a parameter referred to as the number of bits.
Let’s start with a few concrete examples of the arithmetic involved, beginning
with a demonstration of how to multiply two real numbers $x=-.5$ and $y=.5$
using 2 bits. Form the _precision vector_ $\mathbf{p}=(-1,.5)$ with
corresponding _precision matrix_
$P=\mathbf{p}^{t}\mathbf{p}=\begin{pmatrix}1&-.5\\\ -.5&.25\end{pmatrix}$. Set
$\mathbf{x}_{b}=\begin{pmatrix}1\\\ 1\end{pmatrix}$ and
$\mathbf{y}_{b}=\begin{pmatrix}0\\\ 1\end{pmatrix}$ so that

(2) $\displaystyle x=\mathbf{p}\cdot\mathbf{x}_{b},\ \
y=\mathbf{p}\cdot\mathbf{y}_{b}$ (3) $\displaystyle-.25=x\cdot
y=\mathbf{x}_{b}^{t}\mathbf{p}^{t}\mathbf{p}\mathbf{y}_{b}=(1,1)\begin{pmatrix}1&-.5\\\
-.5&.25\end{pmatrix}\begin{pmatrix}0\\\ 1\end{pmatrix}$

Now let $Q$ be a symmetric matrix, and we shall demonstrate how to compute the
quadratic form $(x,y,z)Q\begin{pmatrix}x\\\ y\\\ z\end{pmatrix}$ using binary
variables. As above, set $z=\mathbf{p}\cdot\mathbf{z}_{b}$ where
$\mathbf{z}_{b}$ is a binary vector and $z$ is the corresponding real number.
We can rewrite the quadratic form as

(4)
$\displaystyle(\mathbf{x}_{b}^{t},\mathbf{y}_{b}^{t},\mathbf{z}_{b}^{t})\begin{pmatrix}\mathbf{p}^{t}&&\\\
&\mathbf{p}^{t}&\\\
&&\mathbf{p}^{t}\end{pmatrix}Q\begin{pmatrix}\mathbf{p}&&\\\ &\mathbf{p}&\\\
&&\mathbf{p}\end{pmatrix}\begin{pmatrix}\mathbf{x}_{b}\\\ \mathbf{y}_{b}\\\
\mathbf{z}_{b}\end{pmatrix}$

The middle three terms can be written more succinctly as $(I_{3}\otimes
p^{t})Q(I_{3}\otimes p)=Q\otimes P$ where $I_{3}$ is the $3\times 3$ identity
matrix and $\otimes$ is the tensor product.
Now we describe the construction in full generality. Given a _precision
vector_
$\mathbf{p}=(-1,\frac{1}{2},\frac{1}{2^{2}},\frac{1}{2^{3}},\ldots\frac{1}{2^{b-1}})$
of length $b$, the set of integer multiples of $\frac{1}{2^{b-1}}$ in the
interval $[-1,1-\frac{1}{2^{b}}]$ is exactly the set

(5) $\displaystyle
C_{1,b}:=\\{\mathbf{p}\cdot\mathbf{x}_{b}|\mathbf{x}_{b}\in\\{0,1\\}^{b}\\}$

We use the sub-scripted $\mathbf{x}_{b}$ as a convention to emphasize that
$\mathbf{x}_{b}$ is a binary vector, i.e. the subscript does not refer to the
number of bits. More generally, the $n$-fold product of $C_{1,b}$ is the set

(6) $\displaystyle
C_{n,b}:=\\{(I_{n}\otimes\mathbf{p})\mathbf{x}_{b}|\mathbf{x}_{b}\in\\{0,1\\}^{nb}\\}$

where $I_{n}$ is the identity matrix, $\otimes$ is the tensor product. The set
$C_{n,b}$ will be referred to as a discretized cube. Let $Q$ be a symmetric
$n\times n$ matrix, $\mathbf{r}$ an $n$-vector and suppose we want to solve
the constrained quadratic programming problem

(7) $\displaystyle\underset{\mathbf{x}\in[-1,1]^{n}}{\text{argmin
}}\mathbf{r}^{t}\mathbf{x}+\mathbf{x}^{t}Q\mathbf{x}$

To get an approximate solution to 7 using a QUBO, first replace the unit cube
by a discretized unit cube to get the optimization problem.

(8) $\displaystyle\underset{\mathbf{x}\in C_{n,b}}{\text{argmin
}}\mathbf{r}^{t}\mathbf{x}+\mathbf{x}^{t}Q\mathbf{x}$

Setting $m=nb$, $\mathbf{x}=(I_{n}\otimes\mathbf{p})\mathbf{x}_{b}$, and
$P=\mathbf{p}^{t}\mathbf{p}$, this is equivalent to the QUBO

(9) $\displaystyle\underset{\mathbf{x}_{b}\in\\{0,1\\}^{m}}{\text{argmin }}\
\mathbf{r}^{t}(I_{n}\otimes\mathbf{p})\mathbf{x}_{b}$
$\displaystyle+\mathbf{x}_{b}^{t}(Q\otimes P)\mathbf{x}_{b}$ (10)
$\displaystyle=\underset{\mathbf{x}_{b}\in\\{0,1\\}^{m}}{\text{argmin }}\
\mathbf{x}_{b}^{t}\text{Diag}(\mathbf{r}^{t}I_{n}\otimes\mathbf{p})\mathbf{x}_{b}+\mathbf{x}_{b}^{t}(Q\otimes
P)\mathbf{x}_{b}$ (11)
$\displaystyle=\underset{\mathbf{x}_{b}\in\\{0,1\\}^{m}}{\text{argmin }}\
\mathbf{x}_{b}^{t}(\text{Diag}(\mathbf{r}^{t}I_{n}\otimes\mathbf{p})+Q\otimes
P)\mathbf{x}_{b}$

Here $\text{Diag}(\mathbf{v})$ refers to the diagonal matrix with entries from
$\mathbf{v}$. Going from lines 9 to 10 uses the following identity valid for
binary vectors:
$\mathbf{v}^{t}\mathbf{x}_{b}=\mathbf{x}_{b}^{t}\text{Diag}(\mathbf{v})\mathbf{x}_{b}$.
To summarize, given a number of bits $b$, this procedure approximates the real
optimization problem in $n$ variables 7 with the QUBO 11 of size $n\cdot b$.
We conclude by remarking that we are not restricted to the cube $[-1,1]^{n}$.
If we instead want to optimize over the cube $[-\delta,\delta]^{n}$, repeat
the same construction with the precision vector $\delta\cdot\mathbf{p}$.
It is worth briefly discussing the error introduced by replacing the real cube
with a discretized cube. Suppose we had an ideal annealer that always produces
a best solution to 8, call that solution $\mathbf{x}_{A}$. Let
$\mathbf{x}_{T}$ be a best solution to 7, and let $\hat{\mathbf{x}_{T}}$ be
the point in $C_{n,b}$ closest to $\mathbf{x}_{T}$. The gradient of
$\mathbf{r}^{t}\mathbf{x}+\mathbf{x}^{t}Q\mathbf{x}$ is
$2Q\mathbf{x}+\mathbf{r}$, and
$||\hat{\mathbf{x}_{T}}-\mathbf{x}_{T}||\leq\sqrt{\frac{n}{2^{b}}}$, leading
to the Lipschitz estimate

(12)
$\displaystyle|\mathbf{r}^{t}\hat{\mathbf{x}_{T}}+\hat{\mathbf{x}_{T}}^{t}Q\hat{\mathbf{x}_{T}}-(\mathbf{r}^{t}\mathbf{x}_{T}+\mathbf{x}_{T}^{t}Q\mathbf{x}_{T})|\leq\sqrt{\frac{n}{2^{b}}}\sup_{[-1,1]^{n}}||2Q\mathbf{x}+\mathbf{r}||$

combining with the inequality
$\mathbf{r}^{t}\mathbf{x}_{T}+\mathbf{x}_{T}^{t}Q\mathbf{x}_{T}\leq\mathbf{r}^{t}\mathbf{x}_{A}+\mathbf{x}_{A}^{t}Q\mathbf{x}_{A}\leq\mathbf{r}^{t}\hat{\mathbf{x}_{T}}+\hat{\mathbf{x}_{T}}^{t}Q\hat{\mathbf{x}_{T}}$
implies

(13)
$\displaystyle|\mathbf{r}^{t}\mathbf{x}_{A}+\mathbf{x}_{A}^{t}Q\mathbf{x}_{A}-(\mathbf{r}^{t}\mathbf{x}_{T}+\mathbf{x}_{T}^{t}Q\mathbf{x}_{T})|\leq\sqrt{\frac{n}{2^{b}}}\sup_{[-1,1]^{n}}||2Q\mathbf{x}+\mathbf{r}||$

This inequality makes a trade-off apparent. With more bits, the solution on
the discretized cube will better approximate the true solution of 7 but will
require solving a larger QUBO. Indeed, numerical experiments from subsequent
sections will show that $b=2$ generally requires more iterations than $b=8$,
indicating that the quality of the approximate solution at each step is worse,
although interestingly using $b=2$ takes less time overall since solving
smaller QUBOs is much faster. It is also worth noting that estimate 12 does
not control the actual distance between solutions
$||\mathbf{x}_{A}-\mathbf{x}_{T}||$. In the special case when $Q$ is positive
definite, this distance can be controlled, but it would be interesting to have
estimates in greater generality.
The annealers that one works with in practice are never ideal, and so will
rarely return the absolute best solution $\mathbf{x}_{A}$ but instead a
response consisting of many samples
$\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{l}$ of good solutions with
energies $E(\mathbf{x}_{i})\leq E(\mathbf{x}_{i+1})$. (Here _energy_ of a
solution $\mathbf{x}_{i}$ refers to the value of the objective function at
$\mathbf{x}_{i}$). An obvious approach is to treat the lowest energy solution
$\mathbf{x}_{0}$ as the best approximation of $\mathbf{x}_{T}$. A subtler
approach that can reap great benefits in practice is to take a linear
combination of the full response:

(14)
$\frac{1}{l}\sum_{i=1}^{l}e^{-\beta(E(\mathbf{x}_{i})-E(\mathbf{x}_{0}))}\mathbf{x}_{i}$

as an approximation of $\mathbf{x}_{T}$, where $\beta$ is a parameter.
Experiments in later sections were conducted either using the best response
$\mathbf{x}_{0}$ or the full response with $\beta=100$ and performance is
compared for several values of $n$ and $b$.
Another approach to get better approximations of $\mathbf{x}_{T}$ is to solve
a sequence of QUBOs with bias. More precisely, get an initial approximation of
$\mathbf{x}_{T}$ by following the previous procedure to produce
$\mathbf{x}_{T_{1}}$. Then modify the QUBO by adding a linear term
$-\alpha\mathbf{x}_{T_{1}}^{t}\mathbf{x}$ where $\alpha>0$ and find
approximate solutions to

(15) $\displaystyle\underset{C_{n,b}}{\text{argmin }}\
\mathbf{r}^{t}\mathbf{x}+\mathbf{x}^{t}Q\mathbf{x}-\alpha\mathbf{x}_{T_{1}}^{t}\mathbf{x}=\underset{C_{n,b}}{\text{argmin
}}\
(\mathbf{r}-\alpha\mathbf{x}_{T_{1}})^{t}\mathbf{x}+\mathbf{x}^{t}Q\mathbf{x}$

By Cauchy-Schwartz,
$\frac{\mathbf{x}_{T_{1}}}{||\mathbf{x}_{T_{1}}||}=\underset{||\mathbf{x}||=1}{\text{argmin
}}-\mathbf{x}_{T_{1}}^{t}\mathbf{x}$, thus in solving 15 the annealer is
encouraged to produce solutions in the direction of $\mathbf{x}_{T_{1}}$. The
annealer produces a new lower-energy solution $\mathbf{x}_{T_{2}}$ and this
process can be repeated until the new solution no longer has lower energy than
the previous. Experimental results in later sections contain data with
$\alpha=0$ and $.1$. Biasing is most helpful in the initial phase of the
algorithm when it is iteratively producing solutions close to previous
solutions. In later phases biasing is less useful, as will be evident in
Figures 4 and 5.

###  2.2 An Iterative Descent Algorithm

Algorithms that solve continuous optimization problems rely on a good initial
guess an an iterative descent rule. These tasks can be formulated as QUBO
problems when trying to minimize the Rayleigh quotient over the unit sphere.

####  2.2.1 Obtaining an Initial Guess

Let $\lambda_{1}\leq\lambda_{2}\leq\ldots\leq\lambda_{n}$ be the eigenvalues
of $A$. To get an initial approximation of $\lambda_{1}$, one can ask to solve

(16) $\displaystyle\underset{\mathbf{x}\in C_{n,b}}{\text{argmin }}\
\mathbf{x}^{t}A\mathbf{x}$

as an approximation of

(17) $\displaystyle\underset{||\mathbf{x}||=1}{\text{argmin
}}\mathbf{x}^{t}A\mathbf{x}$

An immediate problem is that if $A$ is positive definite, the solution to 16
is just $\mathbf{x}=0$. This can be remedied by replacing $A$ with $A-\lambda
I_{n}$, where $\lambda\in(\lambda_{i},\lambda_{i+1})$ for some $i$. The
eigenvectors are unaffected, the eigenvalues can be recovered from the new
matrix and the solutions to

(18) $\displaystyle\underset{\mathbf{x}\in C_{n,b}}{\text{argmin
}}\mathbf{x}^{t}(A-\lambda I_{n})\mathbf{x}$

tend to be long, nonzero vectors very close to the span of eigenvectors which
have negative eigenvalues for $A-\lambda I_{n}$, namely
$\mathbf{v}_{1},\ldots,\mathbf{v}_{i}$. A good initial choice is the average
of the eigenvalues $\lambda=\frac{tr(A)}{n}$ and as the algorithm progresses,
$\lambda$ will decrease towards $\lambda_{1}$. To converge in fewer
iterations, it’s better to choose $\lambda$ close to, but greater than
$\lambda_{1}$, as we will examine later.
These observations lead to the following iterative fixed-point method to
produce a good initial guess for the lowest eigenvector. Initially solve 18
with $\lambda=\frac{tr(A)}{n}$ to produce a guess $\mathbf{v}_{1}$. Update
$\lambda$ using the Rayleigh quotient
$\lambda=\frac{\mathbf{v}_{1}^{t}A\mathbf{v}_{1}}{||\mathbf{v}_{1}||^{2}}$ and
solve 18 again possibly using $\mathbf{v}_{1}$ as a bias vector to produce a
second guess $\mathbf{v}_{2}$. Repeat until $\lambda$ is no longer decreasing.
For small matrices, say $10\times 10$, this procedure often suffices to
produce the lowest eigenvalue with 2-3 digits of accuracy and the
corresponding eigenvector to within a distance of order $.1$ of the true
eigenvector. The descent stage of the algorithm increases the precision to
essentially arbitrary order.

####  2.2.2 Iterative Descent

Suppose we want to minimize a function
$f:\mathbb{R}^{m}\rightarrow\mathbb{R}$, and let $\nabla f$ and $H(f)$ denote
the gradient and Hessian of $f$, respectively. Starting with an initial guess
$\mathbf{x}_{0}$, a common strategy is to Taylor expand $f$ around
$\mathbf{x}_{0}$

(19) $\displaystyle f(\mathbf{x})$ $\displaystyle=f(\mathbf{x}_{0})+\nabla
f\cdot\delta+\delta^{t}\frac{H(f)}{2}\delta+o(||\delta||^{3})$ (20)
$\displaystyle\delta$ $\displaystyle:=\mathbf{x}-\mathbf{x}_{0}$

and choose $\delta$ to minimize $\nabla f\cdot\delta$, which is gradient
descent, or to minimize $\nabla f\cdot\delta+\delta^{t}\frac{H(f)}{2}\delta$,
which includes second order methods such as Newton’s method, BFGS, Newton
Conjugate-Gradient, etc. Once a better solution
$\mathbf{x}_{1}=\mathbf{x}_{0}+\delta$ has been found, Taylor expand around
$\mathbf{x}_{1}$ again and repeat. The proposed algorithm obtains a good
descent direction by using a QUBO to find good approximate solutions to

$\underset{\delta\in C_{n,b}}{\text{argmin }}\nabla
f\cdot\delta+\delta^{t}\frac{H(f)}{2}\delta$

Similar to Newton-CG and BFGS, this method requires computing, but not
inverting, the Hessian matrix, and benefits from a line search which possibly
increases the size of $\delta$. See [9] for more details on classical
optimization algorithms and the benefits of line search. Here the line search
step amounts to minimizing a quadratic, and so the optimal scaling can be
directly computed.
If $k^{th}$ approximate solution $\mathbf{x}_{k}$ is closer to the true
solution than any point in the discretized cube, one cannot expect minimizing
the QUBO to produce a better solution. A key part of the descent phase is
enforcing a minimum step size in addition to the line search so that the
candidate $k+1^{st}$ solution is possibly _worse_ than that $k^{th}$. If the
candidate solution is worse, the algorithm discards the candidate and replaces
the discretized unit cube $C_{n,b}$ by a scaled-down discretized cube $t\cdot
C_{n,b}$ where $t<<1$, which amounts to repeating the procedure outlined in
section 2 with the precision vector $t\cdot\mathbf{p}$. Once the discretized
cube has been scaled down, the algorithm continues running until it needs to
scale down the cube further, or exits having achieved the desired accuracy.

###  2.3 The Algorithm

With the key ingredients covered, we are in a position to present the
algorithm. As a reminder, at each step the objective function is of the form
$f(\mathbf{x})=\mathbf{x}^{t}(A-\lambda)\mathbf{x}$, whose gradient and
Hessian can be calculated as $\nabla f=2(A-\lambda)\mathbf{x}^{t}$,
$H(f)=2(A-\lambda)$, and these formulas are implicitly used in the descent
phase of the algorithm. At several stages, the algorithm solves optimization
problems of the form $\underset{\mathbf{x}\in t\cdot C_{n,b}}{\text{argmin}}\
\mathbf{r}^{t}\mathbf{x}+\mathbf{x}^{t}A\mathbf{x}$. These are turned in to
QUBOs as explained in section 2.1, and annealers are used to minimize the
QUBOs that appear, possibly using full responses or biasing. In subsequent
section the effects of biasing, full responses and other parameters will be
analyzed.

Algorithm 1 Controlled Precision QUBO-based Algorithm to Compute Eigenvectors
of Symmetric Matrices

Inputs: Symmetric $n\times n$ matrix $A$, bits for precision vector $b$,
desired precision $\epsilon_{tol}$

Outputs: Approx smallest evec, eval $\mathbf{v},\lambda$ within
$\epsilon_{tol}$ of true values

$\lambda\leftarrow\frac{tr(A)}{n}$

$H\leftarrow A-\lambda I_{n}$ // Enforcing Indefiniteness

$\mathbf{v}\leftarrow\underset{C_{n,b}}{\text{argmin}}\
\mathbf{x}^{t}H\mathbf{x}$ // Initial Guess Phase

$\mathbf{v}\leftarrow\frac{\mathbf{v}}{||\mathbf{v}||}$

while $\mathbf{v}^{t}A\mathbf{v}<\lambda$ do

$\lambda\leftarrow\frac{\mathbf{v}^{t}A\mathbf{v}}{||\mathbf{v}||^{2}}$

$H\leftarrow A-\lambda\cdot I_{n}$

$\mathbf{v}\leftarrow\underset{C_{n,b}}{\text{argmin
}}\mathbf{x}^{t}H\mathbf{x}$

$\mathbf{v}\leftarrow\frac{\mathbf{v}}{||\mathbf{v}||}$

end while

$precision\leftarrow.1$ // Descent Phase

while $precision>\epsilon_{tol}$ do

$H\leftarrow A-\lambda I_{n}$

$\delta\leftarrow\underset{precision\cdot C_{n,b}}{\text{argmin
}}2\mathbf{v}^{t}H\delta+\delta H\delta$ // Getting Descent Direction

$\delta\leftarrow\delta-\langle\mathbf{v},\delta\rangle\delta$ //
Orthogonalizing $\delta,\mathbf{v}$

$t_{min}=\max\\{-\mathbf{v}^{t}H\delta/(\delta^{t}H\delta),1\\}$ // Line
search step

$\delta\leftarrow t_{min}\cdot\delta$

if
$\frac{(\mathbf{v}+\delta)^{t}A(\mathbf{v}+\delta)}{||\mathbf{v}+\delta||^{2}}<\lambda$
then // Checking if solution improves

$\mathbf{v}\leftarrow\frac{\mathbf{v}+\delta}{||\mathbf{v}+\delta||}$

$\lambda\leftarrow\mathbf{v}^{t}A\mathbf{v}$

else

$precision\leftarrow.1\cdot precision$ // Increasing precision otherwise

end if

end while

return $\mathbf{v},\lambda$

Two steps merit a bit more explanation. Replacing $\delta$ by
$\delta-\langle\mathbf{v},\delta\rangle\delta$ forces $\mathbf{v},\delta$ to
be orthogonal. Since the Rayleigh quotient needs to be optimized over the
sphere, the update direction $\delta$ should be tangent to the sphere at
$\mathbf{v}$, and the tangent space of the sphere at $\mathbf{v}$ is precisely
the set of vectors orthogonal to $\mathbf{v}$. Second, the scaling $\delta$ is
computed as $t_{min}=\max\\{-\mathbf{v}^{t}H\delta/(\delta^{t}H\delta),1\\}$.
The expression $-\mathbf{v}^{t}H\delta/(\delta^{t}H\delta)$ is the line search
step coming from minimizing the quadratic
$(\mathbf{v}+t_{min}\delta)^{t}H(\mathbf{v}+t_{min}\delta)$. Strictly speaking
this quadratic only has a minimum when $\delta^{t}H\delta$ is positive, and in
practice when using this algorithm it almost always is, and if not, set
$t_{min}=1$. A minimum scaling $t_{min}\geq 1$ is enforced so that the
candidate update $\mathbf{v}+t_{min}\delta$ possibly overshoots the exact
solution, resulting in a worse estimate of the lowest eigenvector.
Overshooting is an indication that the discretized cube is no longer fine
enough to produce better solutions, and so the candidate update is discarded
and the discretized cube is scaled down. Intuitively the scaling at each step
should be about the order of $\frac{1}{2^{b-1}}$, and the numerical
experiments below all use $.1$.
Oftentimes in practice one wishes to solve a generalized eigenvalue problem of
the form $A\mathbf{v}=\lambda B\mathbf{v}$. In the case when $A,B$ are
symmetric and $B$ is strictly positive definite, the smallest generalized
eigenvalue minimizes the generalized Rayleigh quotient

(21)
$\displaystyle\lambda=\min_{||\mathbf{x}||=1}\frac{\mathbf{x}^{t}A\mathbf{x}}{\mathbf{x}^{t}B\mathbf{x}}\
\ \mathbf{v}=\underset{||\mathbf{x}||=1}{\text{argmin
}}\frac{\mathbf{x}^{t}A\mathbf{x}}{\mathbf{x}^{t}B\mathbf{x}}$

The following small changes solves the generalized eigenvalue problem, again
to essentially arbitrary precision. First, instead of initializing $\lambda$
as $\frac{tr(A)}{n}$, one can generate a random unit vector $\mathbf{w}$ (or
use a specified vector) and initialize
$\lambda=\frac{\mathbf{w}^{t}A\mathbf{w}}{\mathbf{w}^{t}B\mathbf{w}}$. Second,
replace every Rayleigh quotient with the corresponding generalized Rayleigh
quotient. Lastly, instead of updating $H$ as $H=A-\lambda I_{n}$, update as
$H=A-\lambda B$, as the latter preserves the $B$-eigenspectrum of $A$, while
the former does not.
We conclude by emphasizing that this algorithm reaches arbitrary precision
without increasing the size of the QUBOs, all of which involve $n\cdot b$
binary variables. Additionally, all quadratic problems are of the form
$\mathbf{r}^{t}\mathbf{x}+\mathbf{x}^{t}(A-\lambda)\mathbf{x}$ where $A$ is
fixed, implying that the potential non-zero coefficients of the QUBO do not
change (examine formula 11).

##  3 Experimental results

The algorithm and its variants are tested on a class of random matrices of
varying sizes. For each experiment, the algorithm ran until the Rayleigh
quotient of the approximate eigenvector was within $10^{-8}$ of the true
value.

###  3.1 Basic Performance

First, we demonstrate the convergence as a function of the number of
iterations using example matrices from the TAMU SuiteSparse collection [6].
Figure 1 shows performance on the breasttissue_10NN matrix, a weighted graph
adjacency matrix of size $106\times 106$ for $2,4,6$ and $8$ bits using best
response and no biasing.

Figure 1: Digits of accuracy plotted against number of iterations. _Total
annealing time_ is the total amount of time the algorithm spends on performing
SA, since SA dominates the cost of the method. _Evec error_ is the distance of
the computed vector from the true unit eigenvector. _Precision_ refers to the
scaling applied to the discretized cube at each iteration. The initial guess
phase corresponds to a precision of 1. Observe that whenever the error
increases, the algorithm responds by increasing the precision, which often
gives large accuracy gains within the subsequent 2-3 iterations.

Interestingly using fewer bits gives less time to reach desired accuracy
despite requiring more iterations. A $log-log$ regression on the MP matrices
(see next section) gives that the anneal time grows like $(n\cdot b)^{1.57}$
and the number of iterations grows like $n^{.44}b^{-.32}$ so the total time is
roughly $n^{2}b^{1.2}$. The algorithm works on even larger matrices, as is
demonstrated in Figure 2 using the spaceShuttleEntry_1 matrix, a $560\times
560$ control matrix.

Figure 2: Error plot for $560\times 560$ Space Shuttle Control Matrix

Figure 3 demonstrates the performance for the generalized eigenvalue problems
using mesh1em1 as the $A$ matrix and meshe1 as the $B$ matrix, two $48\times
48$ matrices from the SuiteSparse database.

Figure 3: Error plots for Generalized Eigenvalue Problem on two $48\times 48$
mesh matrices

In order to try to get algorithms that run as fast as possible, one might ask
if it is possible get the algorithm to work using 1 bit of precision. With the
current scheme this cannot be done. However, by reformulting the problem as an
Ising instead of QUBO, one can indeed use only one bit precision. Ising
problems are of the form

$\underset{\mathbf{x}\in\\{-1,1\\}^{n}}{\mathrm{argmin}}\mathbf{h}^{t}\mathbf{x}+\mathbf{x}^{t}Q\mathbf{x}$

the main distinction from QUBOs being the spin variables $\pm 1$. QUBOs or
Ising problems are mathematically equivalent, and most annealers are capable
of solving either.
Using the Ising formulation, its possible to mimic the same algorithm, which
works well on very small matrices. However, for larger matrices, such as for
the $106\times 106$ weighted adjacency matrix the single-bit version of the
algorithm takes longer than using two bits, taking $32.3$ seconds and
requiring over 400 iterations (compare with Figure 1). An educated guess for
why this might happen follows. Since the solutions produced by Ising problems
have coordinates that are all non-zero and of the same magnitude, if the
algorithm has already produced a solution whose $k^{th}$ coordinate is close
to the true value, the added solution from the Ising problem will force that
coordinate away from the optimal value. Using two bits is effective because
the solutions can have coordinates that are positive, negative or zero.

###  3.2 Analysis of Parameters

To demonstrate the effect of biasing and full response parameters, the
algorithm is tested on small matrices of sizes $3$, $10$ and $20$ with number
of bits $b=2,4,6$ and $8$. In the interest of not overwhelming the reader with
plots and data, only the data for matrices of size $10$ and $20$ is displayed.
We analyze the error at the end of the initial guess phase, and the average
number of iterations each method requires. For each choice of size, bits and
parameters, 10 Marchenko-Pasture matrices [11] with parameter $\lambda=.3$ are
generated and the average errors at the end of the descent phase is recorded.
Ideally this initial phase should end with the smallest possible error before
beginning the descent phase. Towards this end, taking full responses (equation
14) and biases (equation 15) can be very beneficial. However, this benefit
fades as the sizes of the QUBOs increase as one can see from Figures 4 and 5.

Figure 4: Eigenvector Error for MP matrices at the end of initial guess phase.
Observe that for small QUBOs the full response decreases the error
significantly.

Figure 5: Average number of iterations required for MP matrices for different
parameters.

The choice to initialize $\lambda$ as $\frac{tr(A)}{n}$ is motivated by a
desire to produce an initial guess which is close to, but greater than, the
true lowest eigenvector. To demonstrate this effect on $10\times 10$ and
$20\times 20$ matrices, we compare performance initializing as
$\frac{tr(A)}{n}$, which is the average of all the eigenvalues, against
initializing as the highest Gershgorin bound, which upper bounds the maximum
eigenvalue [7]. Choosing an initialization closer to the true eigenvalue often
leads to fewer iterations, although the difference is somewhat small and fades
as the number of bits increases, as seen in Figure 6.

Figure 6: Total Iterations required for different initializations of $\lambda$

###  3.3 Gap Size Analysis

Here we analyze the effect of the spacing between eigenvalues. In particular,
the gap $|\lambda_{1}-\lambda_{2}|$ can significantly affect the number of
iterations required to reach a given precision. For this experiment, given a
gap size $g=|\lambda_{1}-\lambda_{2}|$, an orthogonal matrix $U$ is chosen at
random with respect to the Haar measure using the SciPy implementation of
[12]. The algorithm is then analyzed on the matrix
$U^{t}\text{Diag}(0,g,1,\ldots,n-2)U$. As Figure 7 demonstrates, as the gap
size decreases the algorithm takes longer to achieve a given accuracy. The
exception is when the gap is 0, and the smallest eigenvalue appears with
multiplicity. In this case the algorithm actually requires fewer iterations.

Figure 7: Average number of iterations on 30 samples as a function of the gap
size $|\lambda_{1}-\lambda_{2}|$. The smaller the gap, the more iterations
required, especially when the number of bits is small. In the extreme case
where the gap is $0$ and the lowest eigenvalue appears with multiplicity, the
algorithm is actually faster in the sense that fewer iterations are needed for
the computed eigenvalue to approximate the true eigenvalue.

Figure 8 has two example error plots demonstrating the slower convergence.
Observe that the eigenvector error relative to both the precision and
eigenvalue error increases as the gap size decreases.

Figure 8: Sample plots for two $10\times 10$ matrices with gap size $.1$ and
$.01$. As the gap size decreases, the ratio of eigenvector error to precision
and eigenvalue error increases.

In the case when there is degeneracy, that is the gap $g$ is $0$, one might
want two eigenvectors that span the eigenspace. This can be accomplished by
running the algorithm once to get an approximate eigenvector $\mathbf{v}_{1}$,
replace the matrix $A$ with $A+\alpha\mathbf{v}_{1}\mathbf{v}_{1}^{t}$ where
$\alpha>0$, and run the algorithm again to get the eigenvector
$\mathbf{v}_{2}$. By the spectral theorem for the symmetic matrix
$A+\alpha\mathbf{v}_{1}\mathbf{v}_{1}^{t}$,
$\mathbf{v}_{1}^{t}\mathbf{v}_{2}=0$ implying that $\mathbf{v}_{2}$ is an
eigenvector for $A$. Replacing $A$ by
$A+\alpha\mathbf{v}_{1}\mathbf{v}_{1}^{t}$ is necessary for numeric purposes.
The gap $g$ is never numerically zero, so if the algorithm is run twice on the
matrix $A$ even with different randomization, it will often produce the same
vector. One can also try to take advantage of the first computation by
initializing the approximate eigenvalue to $\lambda_{n}+\epsilon$ in the
second run of the algorithm. The data shown below in Figure 9 was collected
for $\alpha=\frac{tr(A)}{n}-\lambda_{n}$ and $\epsilon=1$, and one can see a
slight boost in performance in the second run of the algorithm.

Figure 9: Average number of iterations running the algorithm twice on matrices
with complete degeneracy. One can see slightly better performance on the
second run, particularly on the $20\times 20$ matrices

##  4 Conclusion

We have proposed and tested an algorithm to find eigenvectors of symmetric
matrices by minimizing the corresponding Rayleigh quotient with an iterative
steepest-descent method. Initial guesses and subsequent descent directions are
found by looking for minima over discretized cubes of various sizes, encoded
as QUBO problem which is in turn solved with a SA method. The algorithm is
able to reach essentially arbitrary precision even for fairly large matrices.
We have performed a thorough study of the effect of the different parameters,
including, the eigenvalue spacing, initial guesses, and number of bits, and
the matrix size. We have explored the possibility of using a single bit
precision by reformulating the QUBO problem as an Ising problem. Finally, we
have introduced two novel approaches to accelerate the convergence such as
biasing and using a larger set of solution from the SA step. These two
approaches might be applicable to other QUBO based problems. We encourage the
reader to test these algorithms on other annealing devices.

## References

  * [1] S. a Rahman, R. Lewis, E. Mendicelli, and S. Powell, Su(2) lattice gauge theory on a quantum annealer, (2021), https://doi.org/10.13140/RG.2.2.17160.75522.
  * [2] M. Aramon, G. Rosenberg, E. Valiante, T. Miyazawa, H. Tamura, and H. G. Katzgraber, Physics-inspired optimization for quadratic unconstrained problems using a digital annealer, Frontiers in Physics, 7 (2019), p. 48, https://doi.org/10.3389/fphy.2019.00048, https://www.frontiersin.org/article/10.3389/fphy.2019.00048.
  * [3] F. Barahona, On the computational complexity of Ising spin glass models, Journal of Physics A Mathematical General, 15 (1982), pp. 3241–3253, https://doi.org/10.1088/0305-4470/15/10/028.
  * [4] S. Boixo, V. N. Smelyanskiy, A. Shabani, S. V. Isakov, M. Dykman, V. S. Denchev, M. H. Amin, A. Y. Smirnov, M. Mohseni, and H. Neven, Computational multiqubit tunnelling in programmable quantum annealers, Nature Communications, 7 (2016), p. 10327.
  * [5] C. C. Chang, A. Gambhir, T. S. Humble, and S. Sota, Quantum annealing for systems of polynomial equations, Scientific Reports, 9 (2019), p. 10258, https://doi.org/10.1038/s41598-019-46729-0, https://doi.org/10.1038/s41598-019-46729-0.
  * [6] T. A. Davis and Y. Hu, The university of florida sparse matrix collection, ACM Trans. Math. Softw., 38 (2011), https://doi.org/10.1145/2049662.2049663, https://doi.org/10.1145/2049662.2049663.
  * [7] S. Gershgorin, Über die Abgrenzung der Eigenwerte einer Matrix, Bull. Acad. Sci. URSS, 1931 (1931), pp. 749–754.
  * [8] K. M. Hall, An r-dimensional quadratic placement algorithm, Management Science, 17 (1970), pp. 219–229, https://doi.org/10.1287/mnsc.17.3.219, https://doi.org/10.1287/mnsc.17.3.219.
  * [9] R. B. J. Stoer, Introduction to Numerical Analysis, Springer, New York, NY, 1980.
  * [10] G. A. Kochenberger, F. Glover, B. Alidaee, and C. Rego, An unconstrained quadratic binary programming approach to the vertex coloring problem, Annals of Operations Research, 139 (2005), pp. 229–241, https://doi.org/10.1007/s10479-005-3449-7, https://doi.org/10.1007/s10479-005-3449-7.
  * [11] V. A. Marčenko and L. A. Pastur, Distribution of eigenvalues for some sets of random matrices, Mathematics of the USSR-Sbornik, 1 (1967), pp. 457–483, https://doi.org/10.1070/sm1967v001n04abeh001994, https://doi.org/10.1070/sm1967v001n04abeh001994.
  * [12] F. Mezzadri, How to generate random matrices from the classical compact groups, Notices of the American Mathematical Society, 54 (2006).
  * [13] L. Page, S. Brin, R. Motwani, and T. Winograd, The pagerank citation ranking: Bringing order to the web., Technical Report 1999-66, Stanford InfoLab, November 1999, http://ilpubs.stanford.edu:8090/422/.  Previous number = SIDL-WP-1999-0120.
  * [14] K. Pearson, LIII. On lines and planes of closest fit to systems of points in space, Nov. 1901, https://doi.org/10.1080/14786440109462720, https://doi.org/10.1080/14786440109462720.
  * [15] A. Teplukhin, B. K. Kendrick, S. Tretiak, and P. A. Dub, Electronic structure with direct diagonalization on a d-wave quantum annealer, Scientific Reports, 10 (2020), https://doi.org/10.1038/s41598-020-77315-4, http://dx.doi.org/10.1038/s41598-020-77315-4.
  * [16] H. Ushijima-Mwesigwa, C. F. A. Negre, and S. M. Mniszewski, Graph partitioning using quantum annealing on the d-wave system, in Proceedings of the Second International Workshop on Post Moores Era Supercomputing, PMES’17, New York, NY, USA, 2017, Association for Computing Machinery, p. 22–29, https://doi.org/10.1145/3149526.3149531, https://doi.org/10.1145/3149526.3149531.