[2]Krishnan Suresh 1]Computer Science, University of Wisconsin, Madison, 1210
W. Dayton Street, Madison, 53706, WI, USA [2]Mechanical Engineering,
University of Wisconsin, Madison, 1513 University Avenue, Madison, 53706, WI,
USA

# Computing a Sparse Approximate Inverse on Quantum Annealing Machines

Sanjay Suresh<EMAIL_ADDRESS><EMAIL_ADDRESS>[  *

###### Abstract

Many engineering problems involve solving large linear systems of equations.
Conjugate gradient (CG) is one of the most popular iterative methods for
solving such systems. However, CG typically requires a good preconditioner to
speed up convergence. One such preconditioner is the sparse approximate
inverse (SPAI).

In this paper, we explore the computation of an SPAI on quantum annealing
machines by solving a series of quadratic unconstrained binary optimization
(QUBO) problems. Numerical experiments are conducted using both well-
conditioned and poorly-conditioned linear systems arising from a 2D finite
difference formulation of the Poisson problem.

###### keywords:

QUBO; Linear system of equations; Quantum annealing; Conjugate gradient; Pre-
conditioner; sparse approximate inverse; D-WAVE; Quantum computing

##  1 Introduction

Many engineering problems result in linear systems of equations [1] of the
form:

$\mathbf{K}\mathbf{u}=\mathbf{f}$ (1)

where $\mathbf{K}$ is a symmetric, positive-definite, sparse $N\times N$
matrix, $\mathbf{u}$ is the unknown field (for example, the temperature
field), and $\mathbf{f}$ is the applied force (for example, the heat-flux).
Solving such linear systems for large $N$ is a computationally intensive task
[2], and can be time-consuming on classical computers. Quantum computers have
been proposed as an alternate since they can potentially accelerate the
computation; see [3, 4] for recent reviews on the potential role of quantum
computers in engineering.

In particular, the Harrow-Hassidim-Lloyd (HHL) algorithm is a landmark
strategy for solving linear systems of equations on _quantum-gate computers_.
In theory, it offers an exponential speed-up over classical algorithms [5],
and it has been further improved recently [6, 7, 8]. However, due to the
accumulation of errors in current noisy intermediate-scale quantum (NISQ)
computers [9], the HHL algorithm and its variants are limited, in practice, to
very small ($N<5$) systems. Other strategies such as Grover’s algorithm [10]
and quantum approximate optimization algorithm (QAOA) [11] have been proposed
to solve linear systems on quantum gate computers, but they suffer from the
same limitation. Hybrid solvers such as the variational quantum linear solver
[12] have also been proposed to mitigate some of these challenges.

In parallel, _quantum annealing machines_ , such as the D-wave systems with
several thousand qubits [13], have also been proposed for solving linear
systems since they are less susceptible to noise [14, 15]. The basic principle
is to pose the solution of a linear system of equations as a minimization
problem and then convert this into a series of quadratic unconstrained binary
optimization (QUBO) problems. For example, O’Malley and Vesselinov used a
least-squares formulation and a finite-precision qubit representation to pose
QUBO problems [16]. Borle and Lomonaco carried out a theoretical analysis of
this approach [17]. Park et. al. showed how QUBO problems can be simplified
using matrix congruence [18], while its application in solving 1D Poisson
problems was demonstrated in [19]. If the matrix is positive-definite, which
is often the case in engineering problems, it is much more efficient to use a
potential-energy formulation, rather than the least-squares formulation, to
pose QUBO problems. Using the potential-energy formulation, Srivastava et. al.
described a box algorithm to solve QUBO problems arising from finite element
analysis of one-dimensional differential equations [20].

However, despite these advances, only small-size ($N<100$) problems have been
demonstrated on current quantum annealing machines, while most engineering
problems of practical interest result in much larger-size ($N>10^{5}$)
matrices. To address this gap, we propose here an alternate paradigm where
_quantum annealing machines are used, not to directly solve such large linear
systems, but to compute sparse preconditioners (SPAI)_. Once again, the
strategy is to rely on the potential energy formulation to pose a series of
quadratic unconstrained binary optimization (QUBO) problems. However, the
sparsity can be fully exploited to significantly reduce the size of the QUBO
problems, making it more amenable to quantum computing. The computed SPAI can
then be used as a preconditioner to rapidly solve large systems of equations
using the conjugate gradient (CG) method on classical computers. The efficacy
of this approach is demonstrated by solving both well-conditioned and poorly-
conditioned linear systems of equations arising from finite difference
formulation of 2D Poisson problems. Furthermore, by exploiting the well-
structured nature of the finite-difference formulation, we show how one can
compute the SPAI in constant time, independent of the size ($N$) of the linear
system.

##  2 Proposed Methodology

###  2.1 Poisson Problem

To provide a context to this paper, we consider solving a Poisson problem in
2D, governed by [21]:

$k\bigg{(}\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial
y^{2}}\bigg{)}=-f$ (2)

where $u$ is the unknown (temperature) field, $f$ is the (heat) source, and
$k$ is the underlying material property (example: conductivity). We will
assume that the field takes a value of zero on the boundary.

A classic approach for solving Eq. 2 is the finite difference method [21],
where the geometry is sampled by a uniform grid, over which the field $u$ is
to be determined. Then, the partial derivatives are approximated as follows
[21]:

$\frac{\partial^{2}u}{\partial
x^{2}}\Big{\rvert}_{m,n}\approx\frac{u_{m+1,n}-2u_{m,n}+u_{m-1,n}}{\Delta
x^{2}}$ (3) $\frac{\partial^{2}u}{\partial
y^{2}}\Big{\rvert}_{m,n}\approx\frac{u_{m,n+1}-2u_{m,n}+u_{m,n-1}}{\Delta
y^{2}}$ (4)

Substituting the above approximations in Eq. 2, and with $\Delta x=\Delta
y=h$, results in:

$4ku_{m,n}-k\left(u_{m-1,n}+u_{m+1,n}+u_{m,n-1}+u_{m,n+1}\right)=f_{m,n}h^{2}$
(5)

When applied at all grid points, and after eliminating the rows and columns
corresponding to the boundary nodes, we arrive at the linear system in Eq.
(1). As one can observe, the resulting $\mathbf{K}$ matrix is symmetric.
Further, for this problem, the sparsity is $5$, i.e., $\mathbf{K}$ has at most
5 entries in any row or column. Once the rows and columns corresponding to the
boundary nodes are eliminated, and the resulting $\mathbf{K}$ matrix can be
shown to be positive-definite [21].

While the finite difference method applies to arbitrary domains, we will focus
on rectangular domains for simplicity; see Fig 1. Since the field is assumed
to take a value of zero on the boundary, we will only solve for the field in
the interior, i.e., $N$ denotes the number of interior nodes.

Figure 1: Finite-difference grid over a rectangle with zero Dirichlet boundary
conditions.

Finally, if the material property ($k$) is not a constant over the domain, Eq.
5 can be generalized to

$4k_{m,n}u_{m,n}-\left(k_{m,n}^{m-1,n}u_{m-1,n}+k_{m,n}^{m+1,n}u_{m+1,n}+...\right)=f_{m,n}h^{2}$
(6)

where $k_{m,n}$ is the average $k$ over the (4 or less) elements surrounding
the node $(m,n)$, while $k_{m,n}^{m-1,n}$ is the average $k$ over the (2 or 1)
elements adjacent to the edge joining node $(m,n)$ and $(m-1,n)$. The
resulting $\mathbf{K}$ matrix is still symmetric, sparse, and positive
definite, but can become poorly conditioned depending on the material
distribution (see Section 3).

###  2.2 Classic Linear Solvers

The two main strategies for solving Eq.(1) are: direct and iterative. Direct
solvers usually rely on LU-decomposition [22], and generally require a
significant amount of memory for large systems. Iterative solvers require less
memory but converge to a solution gradually. The rate of convergence of
iterative methods depends on several factors including the condition number of
$\mathbf{K}$, sparsity, etc [23]. One of the most popular iterative methods is
conjugate gradient whose run-time complexity is given by [24, 25]:

$CG\sim O\left(Ns\sqrt{\kappa}\log{\frac{1}{\epsilon}_{cg}}\right)$ (7)

where $N$ is the dimension of $\mathbf{K}$, $s$ is the sparsity ($s=5$ for the
problem described above), $\kappa$ is the condition number, and
$\epsilon_{cg}$ is the desired residual error. Thus, for poorly conditioned
systems (i.e. when $\kappa$ is large), CG does not perform well, and
preconditioners are essential. Several preconditioners are widely used today
[26, 27, 28]; these include Jacobi, incomplete Cholesky, sparse approximate
inverse, etc. In this paper, we will rely on the sparse approximate inverse of
$\mathbf{K}$, and we propose a simple algorithm to compute this preconditioner
on quantum annealing computers using a QUBO formulation.

###  2.3 Computing a Sparse Approximate Inverse

Note that $\mathbf{K}^{-1}$, the exact inverse of $\mathbf{K}$, can be
computed by solving

$\mathbf{K}\mathbf{y}_{j}=\mathbf{e}_{j},j=0,1...N-1$ (8)

where $\mathbf{e}_{j}$ is the unit vector corresponding to the $j^{th}$
dimension. In general, $\mathbf{K}^{-1}$ will be dense. Our objective is to
compute a _sparse approximate inverse_ (SPAI) $\mathbf{M}$. There are various
_a priori_ and adaptive techniques for forcing sparsity on $\mathbf{M}$; see
[26, 27]. We will use a well-known _a priori_ technique where the sparsity
pattern of $\mathbf{K}$ is imposed on $\mathbf{M}$ [27]. In other words, to
compute $\mathbf{M}$, we once again solve

$\mathbf{K}\mathbf{m}_{j}=\mathbf{e}_{j},j=0,1...N-1$ (9)

but with the constraint that $\mathbf{m}_{j}$ must have the same sparsity
pattern as the $j^{th}$ column of $\mathbf{K}$.

To compute $\mathbf{m}_{j}$, let the row-sparsity index of the $j^{th}$ column
of $\mathbf{K}$ be $\mathbf{s}$, i.e., $\mathbf{s}[i]$ stores the $i^{th}$
non-zero row of $\mathbf{K}[:,j]$, $0\leq i\leq s-1$, and $s\ll N$. We can
rearrange $\mathbf{m}_{j}$ as

$\mathbf{m}_{j}=\begin{bmatrix}\mathbf{\hat{m}}_{j}\\\
\mathbf{0}\end{bmatrix}$ (10)

where $\mathbf{\hat{m}}_{j}$ is of length $s$. We can also rearrange
$\mathbf{K}$ using the same reordering as:

$\mathbf{K}=\begin{bmatrix}\mathbf{A}_{j}&\mathbf{B}_{j}\\\
\mathbf{B}_{j}^{T}&\mathbf{C}_{j}\end{bmatrix}$ (11)

where $\mathbf{A}_{j}$ is $s\times s$, and $\mathbf{C}_{j}$ is
$(N-s)\times(N-s)$, and rearrange $\mathbf{e}_{j}$ as

$\mathbf{e}_{j}=\begin{bmatrix}\mathbf{\hat{e}}_{j}\\\
\mathbf{0}\end{bmatrix}$ (12)

Thus Eq. 9 reduces to:

$\begin{bmatrix}\mathbf{A}_{j}&\mathbf{B}_{j}\\\
\mathbf{B}_{j}^{T}&\mathbf{C}_{j}\end{bmatrix}\begin{bmatrix}\mathbf{\hat{m}}_{j}\\\
\mathbf{0}\end{bmatrix}=\begin{bmatrix}\mathbf{\hat{e}}_{j}\\\
\mathbf{0}\end{bmatrix}$ (13)

Discarding the last $(N-s)$ equations, we arrive at :

$\mathbf{A}_{j}\mathbf{\hat{m}}_{j}=\mathbf{\hat{e}}_{j}$ (14)

where $\mathbf{A}_{j}$ is a small $s\times s$ matrix constructed from the
sparsity pattern of the $j^{th}$ column. Eq. 14 must be posed and solved for
each column $j$ of $\mathbf{K}$. However, we show in Section 2.5 that we only
need to solve for a few columns by exploiting the structured nature of the
finite difference formulation.

###  2.4 QUBO Formulation

Observe that solving Eq. (14) is equivalent to minimizing the potential energy
(we have dropped the subscript $j$ to avoid clutter):

$\min_{\mathbf{\hat{m}}}\Pi=\frac{1}{2}\mathbf{\hat{m}}^{T}\mathbf{A}\mathbf{\hat{m}}-\mathbf{\hat{m}}^{T}\mathbf{\hat{e}}$
(15)

This is a quadratic unconstrained minimization problem involving real
variables $\mathbf{\hat{m}}$. To solve this on a quantum annealing machine, we
represent each real component ${\hat{m}}_{i}$ using qubit variables. A well
known strategy is the radix representation, also referred to here as the _box
representation_ [16, 29]:

$\mathbf{\hat{m}}=\mathbf{c}+L(-2\mathbf{q}_{1}+\mathbf{q}_{2})$ (16)

where $\mathbf{c}$ and $L$ are real value parameters that we will iteratively
improve, while $\mathbf{q}_{1}$ and $\mathbf{q}_{2}$ are qubit vectors of
length $s$ each, i.e., a total of $2s$ qubits is used to capture
$\mathbf{\hat{m}}$.

Since $\mathbf{\hat{m}}$ is linear in $\mathbf{q}_{1}$ and $\mathbf{q}_{2}$,
substituting Eq. (16) into Eq. (15) will lead to a quadratic unconstrained
_binary_ optimization (QUBO) problem:

$\min_{\mathbf{q}=\\{\mathbf{q}_{1},\mathbf{q}_{2}\\}}\Pi=\frac{1}{2}\mathbf{q}^{T}\mathbf{Q}^{\prime}\mathbf{q}+\mathbf{q}^{T}\mathbf{d}$
(17)

Further, since the qubit variables can only take the binary values $0$ or $1$,
the linear term can be absorbed into the quadratic term [30], resulting in the
standard form:

$\min_{\mathbf{q}=\\{\mathbf{q}_{1},\mathbf{q}_{2}\\}}\Pi=\frac{1}{2}\mathbf{q}^{T}\mathbf{Q}\mathbf{q}$
(18)

where $\mathbf{Q}$ is symmetric (but not positive definite).

The overall strategy therefore is as follows: For the first column of
$\mathbf{K}$, the parameters $\mathbf{c}$ and $L$ are initialized to
$\mathbf{c}=\mathbf{0}$ and $L=0$. The resulting QUBO problem in Eq. 18 is
solved on a quantum annealing machine. Then $\mathbf{c}$ and $L$ in Eq. 16 are
updated via the _sparse box algorithm_ , discussed in Section 2.6. The process
is repeated until convergence is reached (in typically $25\sim 35$ iterations;
see Section 3). Then, the next column of $\mathbf{K}$ is processed until the
entire $\mathbf{M}$ matrix is constructed.

###  2.5 Node Mapping for a Structured Grid

For a generic $\mathbf{K}$ matrix, one must explicitly process each of the $N$
columns. However, for the structured finite difference grid, one can
significantly reduce the computation. Observe that each column of $\mathbf{K}$
corresponds to a unique node in the grid. Now consider a typical node
highlighted using a square box in Fig. 2. The $5\times 5$ matrix $\mathbf{A}$
in Eq. 15 corresponding to this node (column), is entirely determined by the
rows and columns of $\mathbf{K}$ associated with this node and the 4
neighboring nodes (highlighted using circles in Fig. 2). From Eq. 6, we
observe that the diagonal entries of $\mathbf{A}$ depend only on the material
property ($k_{m,n}$) associated with the square node, while the non-diagonal
entries depend only on the material property (example: $k_{m,n}^{m-1,n}$)
associated with the edge joining the square node and a circle node.

Figure 2: A node and its 4 neighboring nodes.

Since this pattern repeats over the entire grid, the matrix $\mathbf{A}$
associated with the two square nodes in Fig. 3, for example, are identical.
Consequently, the solution vectors $\mathbf{\hat{m}}$ for the two nodes are
identical.

Figure 3: Two nodes with identical columns in $\mathbf{M}$.

With exceptions made at the corner and edge nodes, one can conclude that, for
a single material domain, it is sufficient to compute 9 independent columns of
$\mathbf{M}$, independent of the size $N$. A typical set of these 9
independent nodes is illustrated in Fig. 4. This can be further reduced to 3
independent nodes (one corner, one edge, and one interior node) with
appropriate transformations. Exploiting this _node mapping_ , one can compute
the SPAI matrix $\mathbf{M}$ in constant time, independent of $N$. To the best
of our knowledge, this has not been exploited previously to compute SPAI in
classical, or quantum settings.

Figure 4: For a single material, only 9 nodes need to be considered.

If the domain is composed of two materials as illustrated in Fig. 5, then one
can show that only 21 independent columns of $\mathbf{M}$ need to be computed,
independent of $N$. This can be further reduced to 10 independent columns
through appropriate transformations.

Figure 5: For this two-material configuration, only 21 nodes need to be
considered.

###  2.6 Proposed Algorithms

To summarize, the proposed algorithm to compute $\mathbf{M}$ is described in
Alg. 1 where

  1. 1.

The sparsity pattern of $\mathbf{K}$ is first copied over to $\mathbf{M}$.

  2. 2.

Then, for each column of $\mathbf{K}$, if that column (node) is mapped to a
previously computed column (node), we copy the previously computed solution

  3. 3.

Else we call the _sparse box algorithm_ (see below), and the computed solution
is pushed to $\mathbf{M}$.

  4. 4.

Finally, we force symmetry on $\mathbf{M}$ to address possible numerical
errors.

Algorithm 1 Sparse Approximate Inverse

1:procedure QUBOSparseApproximateInverse($\mathbf{K}$)

2: $\mathbf{M}\leftarrow\mathbf{K}$ $\triangleright$ M has same sparsity as K

3: $N\leftarrow\text{dim}\mathbf{K}$ $\triangleright$ number of rows in K

4: for $j\leftarrow$ $0$ to $N-1$ do

5: if nodeMap[j] $<$ j then $\triangleright$ If node is mapped to another node

6: $\mathbf{M}[:,j]\leftarrow\mathbf{M}[:,nodeMap[j]]$ $\triangleright$ copy
solution

7: else

8: $\mathbf{\hat{m}}\leftarrow$ SparseBoxAlg($\mathbf{K}$, $j$)
$\triangleright$ solve via QUBO

9: $\mathbf{s}\leftarrow\mathbf{K}[:,j]$ $\triangleright$ non-zero rows of
column j

10: $\mathbf{M}[\mathbf{s},j]\leftarrow\mathbf{\hat{m}}$ $\triangleright$ copy
solution

11: end if

12: end for

13: $\mathbf{M}\leftarrow(\mathbf{M}+\mathbf{M}^{T})/2$ $\triangleright$
impose symmetry

14:end procedure$\triangleright$ Output: solution $\mathbf{M}$

The above algorithm uses a sparse box algorithm, which is a generalization of
the box algorithm proposed in [20]. The original box algorithm solves small
dense linear systems using a QUBO formulation (see [20] for details). It is
modified here to solve the sparse problem in Eq. 14. A few observations
regarding the proposed sparse box algorithm (see Alg. 2) are:

  1. 1.

We have chosen an initial box size of $L=1$. This is an arbitrary choice; the
algorithm is robust for any reasonable value [20]; see numerical experiments.
Choosing a large initial value for $L$ will increase the number of contraction
steps, while a small initial value will increase the number of translation
steps.

  2. 2.

A total of 2s qubits are created, and the initial potential energy is zero.

  3. 3.

In the main iteration, the QUBO problem is constructed using the software
packages pyQUBO [30].

  4. 4.

The QUBO problem can be solved in exactly or through (3) quantum annealing;
see numerical experiments.

  5. 5.

If the computed potential energy is less than the current minimum, then we
have found a better solution. In this case, we move the box center
$\mathbf{c}$ to the new computed solution (translation) to improve the
accuracy. Otherwise, the current solution is optimal, and we shrink the box
size $L$ (contraction) to improve the precision.

  6. 6.

For termination, using a very small value, say, $\epsilon_{box}\approx
10^{-14}$, is not desirable since (1) it will increase the computation cost,
and (2) as the box size becomes small, the potential energy function $\Pi$
becomes relatively flat, and the computed potential energies from various
qubit configurations will be numerically equal; this will result in the
algorithm choosing the wrong step (i.e. translate versus contract). Further,
$\epsilon_{box}\approx 10^{-14}$ is not required since we are only solving for
an approximate inverse. A typical choice is $\epsilon_{box}\approx 10^{-6}$;
see numerical experiments.

Algorithm 2 Sparse Box Algorithm

1:procedure SparseBoxAlg($\mathbf{K}$, $j$)

2: $\mathbf{s}\leftarrow\mathbf{K}[:,j]$ $\triangleright$ non-zero rows of
column j

3: $s\leftarrow\text{dim}(\mathbf{s})$ $\triangleright$ sparsity of column j

4: $i\leftarrow\mathbf{(}s[i]==j)$ $\triangleright$ find the entry for column
j

5: $\mathbf{c}\leftarrow\mathbf{0}$ $\triangleright$ center of length $n$

6: $L\leftarrow 1$ $\triangleright$ initialize box size

7: $\mathbf{q}_{1},\mathbf{q}_{2}\leftarrow\text{Qubits}(s)$ $\triangleright$
create qubit arrays of length s

8: $iter=0$

9: $\Pi_{\text{min}}\leftarrow 0$ $\triangleright$ energy initialization

10: repeat$\triangleright$ until convergence

11:
$\mathbf{\hat{m}}\leftarrow\mathbf{c}+L(-2\mathbf{q}_{1}+\mathbf{q}_{2})$$\triangleright$
symbolic expression

12:
$\Pi\leftarrow\frac{1}{2}\mathbf{\hat{m}}^{\intercal}\mathbf{K}[\mathbf{s},\mathbf{s}]\mathbf{\hat{m}}-\mathbf{\hat{m}}[i]$
$\triangleright$ construct QUBO

13:
$\Pi^{*},\mathbf{q}_{1}^{*},\mathbf{q}_{2}^{*}\leftarrow\text{minimize}(H)$
$\triangleright$ solve QUBO

14: if $\Pi^{*}<\Pi_{\text{min}}$ then

15:
$\mathbf{c}\leftarrow\mathbf{c}+L(-2\mathbf{q}_{1}^{*}+\mathbf{q}_{2}^{*})$
$\triangleright$ translation

16: $\Pi_{\text{min}}\leftarrow\Pi^{*}$ $\triangleright$ new minimum

17: else

18: $L\leftarrow L/2$ $\triangleright$ contraction

19: end if

20: $iter=iter+1$

21: until ($L<\epsilon_{box}$) or ($iter>iter_{max}$) $\triangleright$
termination

22:end procedure$\triangleright$ Output: solution $\mathbf{c}$

###  2.7 D-WAVE Embedding

We now consider mapping the QUBO problems generated by Alg. 2 to the D-Wave
Advantage quantum annealing machine. The D-wave Advantage is equipped with
5000+ qubits, embedded in a Pegasus architecture. Each QUBO problem involves
at most 10 logical qubits with the connectivity graph illustrated in Fig.
6(a). Using the default embedding, these logical qubits are mapped to 18
physical qubits illustrated in Fig. 6(b). The default chain strength was found
to be sufficient for all numerical experiments.

((a))

((b))

Figure 6: (a) 10 logical qubits. (b) 18 physical qubits on Pegasus
architecture.

##  3 Numerical Experiments

In the following experiments, we consider a $(g_{x}+2)\times(g_{y}+2)$
rectangular finite-difference grid as shown earlier in Fig. 1. We assume that
the field on the boundary is zero, i.e. we only solve for the interior
$g_{x}\times g_{y}$ grid. The default values, unless otherwise noted, for all
the experiments are:

  * •

The grid size is $g_{x}=401$ and $g_{y}=301$.

  * •

The rectangle is composed of a single material with $k=1$.

  * •

The box-tolerance in Algorithm 2 is set to $\epsilon_{box}=10^{-6}$.

  * •

The box length in Algorithm 2 is initialized to $L=1$.

  * •

The maximum box iterations in Algorithm 2 is set to 100.

  * •

The conjugate gradient residual tolerance is set to $10^{-10}$

In the experiments, Q-PCG refers to the standard preconditioned conjugate
gradient method where the proposewd SPAI preconditioner is used. In each of
the following experiments, we graph the residual error against the number of
iterations for both Q-PCG and CG. The implementation is in Python, and uses
pyQUBO [30] to create the QUBO model, and pyAMG [31] to construct the
$\mathbf{K}$ matrix.

###  3.1 Performance Evaluation

In the first experiment, we compare the convergence of Q-PCG and CG using the
default values listed above; this results in $N=120,701$ (size of the
$\mathbf{K}$ matrix). The convergence plots using D-WAVE’s _dimod_ exact QUBO
solver are illustrated in Fig. 7(a). _Exactly the same Q-PCG convergence plot
was obtained when using the D-WAVE quantum annealing solver._ The SPAI does
not significantly improve the convergence of CG here since $\mathbf{K}$ is
relatively well-conditioned. The computed Poisson field is illustrated in Fig.
7(b).

((a))

((b))

Figure 7: (a) Convergence plot for $g_{x}=401$, $g_{y}=301$. (b) Computed
Poisson field.

The two solvers are compared in Table 1.

Table 1: Comparing the two solvers. | Exact | Quantum
---|---|---
Total QUBO solves | 294 | 294
Avg. time per QUBO solve | 0.9 milliseconds | 26 milliseconds
Time to compute $\mathbf{M}$ | 0.5 seconds | 209 seconds

Additional details on the QPU timing for a single QUBO solve (using the
default 100 samples) are provided in Table 2. As one can observe there is
significant overhead in D-WAVE QPU allocation, programming, access and post-
process.

Table 2: QPU timing per QUBO solve. Task | Milliseconds
---|---
Access | 26
Programming | 15
Sampling | 11
Readout | 7
Post-process | 2
Anneal | 2
Delay | 2
Access overhead | 1.6

###  3.2 Multiple Materials

The real advantage of the SPAI becomes evident when we consider two materials
with the $k_{1}$ on the left and $k_{2}$ on the right (see Fig. 5). For the
default grid size of $g_{x}=401$ and $g_{y}=301$, we observe in Fig. 8(a) that
Q-PCG provides a speed-up of 3.3 when $k_{1}=1$ and $k_{2}=10$, and a speed-up
of 11 when $k_{1}=1$ and $k_{2}=100$. Such multi-material problems are fairly
common in engineering, especially during topology optimization [32, 33].

((a))

((b))

Figure 8: Convergence plots for two-material domain: (a) $k_{1}=1,k_{2}=10$.
(b) $k_{1}=1,k_{2}=100$.

###  3.3 Impact of box tolerance

We now study the impact of the box tolerance on the convergence of Q-PCG. For
a single material problem, with default values, Table 3 summarizes the Q-PCG
iterations and total box iterations using the exact and quantum annealing
solvers (recall from Fig. 7(a) that regular CG converges in about 900
iterations). As the box tolerance is varied between $10^{-8}$ and $10^{-2}$,
one can observe in Table 3 that the number of Q-PCG iterations remains around
459, while the number of box iterations decreases, as expected. However, Q-PCG
did not converge to the desired tolerance, i.e., the SPAI matrix is not
effective when the box tolerance is too coarse ($10^{-1}$).

Table 3: Impact of box tolerance on Q-PCG iterations and (total box iterations). $\epsilon_{box}$ | Exact | Quantum
---|---|---
$10^{-8}$ | 458 (401) | 459 (401)
$10^{-6}$ | 459 (294) | 459 (294)
$10^{-4}$ | 459 (202) | 459 (202)
$10^{-2}$ | 464 (96) | 473 (96)
$10^{-1}$ | \- (50) | \- (50)

###  3.4 Impact of box length

Next we study the impact of the box length on the convergence of Q-PCG. For a
single material problem, with default values, Table 4 summarizes the Q-PCG
iterations and total box iterations using exact and quantum annealing solvers.
As the box length is varied between $10^{-2}$ and $10^{4}$, one can observe in
Table 4 that the number of Q-PCG iterations remains at 459, while the number
of box iterations is minium when $L=1$. A star (*) indicates that the desired
box tolerance was not achieved, and the box algorithm exited when the maximum
iteration was reached. However, even in this case, Q-PCG converged.

Table 4: Impact of box length on Q-PCG iterations and (total box iterations).  $L$ | Exact | Quantum
---|---|---
$10^{4}$ | 458 (437) | 459 (437)
$10^{2}$ | 458 (390) | 459 (365)
$10^{1}$ | 458 (315) | 458 (315)
$10^{0}$ | 458 (294) | 458 (294)
$10^{-1}$ | 458 (333) | 458 (333)
$10^{-2}$ | 691 (431*) | 691 (431*)

##  4 Conclusions

A hybrid classical-quantum strategy for solving large linear systems of
equations was proposed. The strategy relied on computing a sparse approximate
preconditioner (SPAI) on a quantum annealing machine and using this
preconditioner, along with an iterative solver, on a classical machine. Its
effectiveness was demonstrated on large ($N>100,000$) ill-conditioned matrices
arising from a finite-difference formulation of the Poisson problem.

There are many directions for continued research. (1) While we demonstrated
the method using quantum annealing machines, it can be extended to quantum
gate computers since QUBO problems can be solved (approximately) on such
machines via quantum approximate optimization algorithms [34, 35]. (2) The
finite difference formulation of the Poisson problem resulted in matrices with
small sparsity ($s=5$). Other formulations and other field problems would lead
to a larger $s$. For example, 3D structured-grid finite element analysis of
elasticity problems would result in $s=27\times 3=81$. The proposed strategy
and algorithm, in theory, extends to such scenarios. However, the performance
of SPAI and efficient embedding of the QUBO problems need to be investigated.
(3) We selected a simple _a priori_ sparsity pattern for the preconditioner
$M$; adaptive patterns and their impact on Q-PCG need to be explored. (4) One
can potentially exploit parallel annealing [36] for computing $M$. (5) We
limited the box representation to 2 qubits per real variable (see Eq. 16);
extending to multiple qubits will accelerate convergence but will require a
larger number of qubits.

## Compliance with ethical standards

The authors declare that they have no conflict of interest.

## Replication of Results

The Python code pertinent to this paper is available at https://github.com/UW-
ERSL/SPAI.

## Acknowledgments

We would like to thank the Graduate School of UW-Madison for the the Vilas
Associate grant.

## References

  *   * [1] O.C. Zienkiewicz, R.L. Taylor, J.Z. Zhu, _The finite element method: its basis and fundamentals_ (Elsevier, 2005)
  * [2] N.I. Gould, J.A. Scott, Y. Hu, A numerical evaluation of sparse direct solvers for the solution of large sparse symmetric linear systems of equations.  ACM Transactions on Mathematical Software (TOMS) 33(2), 10–es (2007)
  * [3] G. Tosti Balducci, B. Chen, M. Möller, M. Gerritsma, R. De Breuker, Review and perspectives in quantum computing for partial differential equations in structural mechanics.  Frontiers in Mechanical Engineering p. 75 (2022)
  * [4] Y. Wang, J.E. Kim, K. Suresh, Opportunities and challenges of quantum computing for engineering optimization.  Journal of Computing and Information Science in Engineering 23(6) (2023)
  * [5] A.W. Harrow, A. Hassidim, S. Lloyd, Quantum algorithm for linear systems of equations.  Physical review letters 103(15), 150502 (2009)
  * [6] A. Ambainis, Variable time amplitude amplification and a faster quantum algorithm for solving systems of linear equations.  arXiv preprint arXiv:1010.4458 (2010)
  * [7] A.M. Childs, R. Kothari, R.D. Somma, Quantum algorithm for systems of linear equations with exponentially improved dependence on precision.  SIAM Journal on Computing 46(6), 1920–1950 (2017)
  * [8] X. Liu, H. Xie, Z. Liu, C. Zhao, Survey on the improvement and application of HHL algorithm.  Journal of Physics: Conference Series 2333(1), 012023 (2022)
  * [9] J. Preskill, Quantum computing in the NISQ era and beyond.  Quantum 2, 79 (2018)
  * [10] K. Srinivasan, B.K. Behera, P.K. Panigrahi, Solving linear systems of equations by gaussian elimination method using grover’s search algorithm: an ibm quantum experience.  arXiv preprint arXiv:1801.00778 (2017)
  * [11] D. An, L. Lin, Quantum linear system solver based on time-optimal adiabatic quantum computing and quantum approximate optimization algorithm.  ACM Transactions on Quantum Computing 3(2), 1–28 (2022)
  * [12] C. Bravo-Prieto, R. LaRose, M. Cerezo, Y. Subasi, L. Cincio, P.J. Coles, Variational quantum linear solver.  arXiv preprint arXiv:1909.05820 (2019)
  * [13] S.W. Shin, G. Smith, J.A. Smolin, U. Vazirani, How quantum is the d-wave machine?  arXiv preprint arXiv:1401.7087 (2014)
  * [14] P. Hauke, H.G. Katzgraber, W. Lechner, H. Nishimori, W.D. Oliver, Perspectives of quantum annealing: Methods and implementations.  Reports on Progress in Physics 83(5), 054401 (2020)
  * [15] S. Yarkoni, E. Raponi, T. Bäck, S. Schmitt, Quantum annealing for industry applications: Introduction and review.  Reports on Progress in Physics (2022)
  * [16] D. O’Malley, V.V. Vesselinov, B.S. Alexandrov, L.B. Alexandrov, Nonnegative/binary matrix factorization with a d-wave quantum annealer.  PloS one 13(12), e0206653 (2018)
  * [17] A. Borle, S.J. Lomonaco, in _WALCOM: Algorithms and Computation: 13th International Conference, WALCOM 2019, Guwahati, India, February 27–March 2, 2019, Proceedings 13_ (Springer, 2019), pp. 289–301
  * [18] S.W. Park, H. Lee, B.C. Kim, Y. Woo, K. Jun, in _2021 International Conference on Information and Communication Technology Convergence (ICTC)_ (IEEE, 2021), pp. 1363–1367
  * [19] R. Conley, D. Choi, G. Medwig, E. Mroczko, D. Wan, P. Castillo, K. Yu, in _Quantum Computing, Communication, and Simulation III_ , vol. 12446 (SPIE, 2023), pp. 53–63
  * [20] S. Srivastava, V. Sundararaghavan, Box algorithm for the solution of differential equations on a quantum annealer.  Physical Review A 99(5), 052355 (2019)
  * [21] H.P. Langtangen, S. Linge, _Finite difference computing with PDEs: a modern software approach_ (Springer Nature, 2017)
  * [22] M. Bollhöfer, O. Schenk, R. Janalik, S. Hamm, K. Gullapalli, State-of-the-art sparse direct solvers.  Parallel algorithms in computational science and engineering pp. 3–33 (2020)
  * [23] O. Axelsson, in _Sparse Matrix Techniques: Copenhagen 1976 Advanced Course Held at the Technical University of Denmark Copenhagen, August 9–12, 1976_ (Springer, 2007), pp. 1–51
  * [24] J.R. Shewchuk, et al.  An introduction to the conjugate gradient method without the agonizing pain (1994)
  * [25] J.L. Nazareth, Conjugate gradient method.  Wiley Interdisciplinary Reviews: Computational Statistics 1(3), 348–353 (2009)
  * [26] E. Chow, A priori sparsity patterns for parallel sparse approximate inverse preconditioners.  SIAM Journal on Scientific Computing 21(5), 1804–1822 (2000)
  * [27] M. Benzi, Preconditioning techniques for large linear systems: a survey.  Journal of computational Physics 182(2), 418–477 (2002)
  * [28] A.J. Wathen, Preconditioning.  Acta Numerica 24, 329–376 (2015)
  * [29] M.L. Rogers, R.L. Singleton Jr, Floating-point calculations on a quantum annealer: Division and matrix inversion.  Frontiers in Physics 8, 265 (2020)
  * [30] M. Zaman, K. Tanahashi, S. Tanaka, Pyqubo: Python library for mapping combinatorial optimization problems to qubo form.  IEEE Transactions on Computers 71(4), 838–850 (2021)
  * [31] N. Bell, L.N. Olson, J. Schroder, B. Southworth, PyAMG: Algebraic multigrid solvers in python.  Journal of Open Source Software 8(87), 5495 (2023).  10.21105/joss.05495.  URL https://doi.org/10.21105/joss.05495
  * [32] W. Zuo, K. Saitou, Multi-material topology optimization using ordered simp interpolation.  Structural and Multidisciplinary Optimization 55, 477–491 (2017)
  * [33] K. Suresh, Efficient generation of large-scale pareto-optimal topologies.  Structural and Multidisciplinary Optimization 47(1), 49–61 (2013)
  * [34] E. Farhi, J. Goldstone, S. Gutmann, A quantum approximate optimization algorithm.  arXiv preprint arXiv:1411.4028 (2014)
  * [35] B.D. Clader, B.C. Jacobs, C.R. Sprouse, Preconditioned quantum linear system algorithm.  Physical review letters 110(25), 250504 (2013)
  * [36] E. Pelofske, G. Hahn, H.N. Djidjev, Parallel quantum annealing.  Scientific Reports 12(1), 4499 (2022)


# A family of stabilizer-free virtual elements on triangular meshes

Xuejun Xu  School of Mathematical Science, Tongji University, Shanghai,
200092, China  Institute of Computational Mathematics, AMSS, Chinese Academy
of Sciences, Beijing, 100190, China<EMAIL_ADDRESS>and  Shangyou Zhang
Department of Mathematical Sciences, University of Delaware, Newark, DE 19716,
USA<EMAIL_ADDRESS>
###### Abstract.

A family of stabilizer-free $P_{k}$ virtual elements are constructed on
triangular meshes. When choosing an accurate and proper interpolation, the
stabilizer of the virtual elements can be dropped while the quasi-optimality
is kept. The interpolating space here is the space of continuous $P_{k}$
polynomials on the Hsieh-Clough-Tocher macro-triangle, where the macro-
triangle is defined by connecting three vertices of a triangle with its
barycenter. We show that such an interpolation preserves $P_{k}$ polynomials
locally and enforces the coerciveness of the resulting bilinear form.
Consequently the stabilizer-free virtual element solutions converge at the
optimal order. Numerical tests are provided to confirm the theory and to be
compared with existing virtual elements.

Key words. virtual element, stabilizer free, elliptic equation, Hsieh-Clough-
Tocher macro-triangle, triangular mesh.

AMS subject classifications. 65N15, 65N30

##  1\. Introduction

In this work, we construct a family of stabilizer-free $P_{k}$ virtual
elements ([4, 5, 9, 10, 11, 12, 13, 15, 16, 20, 21, 22]) on triangular meshes.

For solving the following model equation,

(1.1) $\displaystyle\begin{aligned} -\Delta u&=f&&\hbox{\rm{in }}\ \Omega,\\\
u&=0&&\hbox{\rm{on }}\ \partial\Omega,\end{aligned}$

where $\Omega\subset\mathbb{R}^{2}$ is a bounded polygonal domain and $f\in
L^{2}(\Omega)$, the weak form reads: Find $u\in H^{1}_{0}(\Omega)$ such that

(1.2) $\displaystyle\begin{aligned} (\nabla u,\nabla v)&=(f,v)&&\forall v\in
H^{1}_{0}(\Omega),\end{aligned}$

where $(\cdot,\cdot)$ denotes the $L^{2}$ inner product on $\Omega$ and we
have $|v|_{1}^{2}=(\nabla v,\nabla v)$.

Let $\mathcal{T}_{h}=\\{K\\}$ be a quasi-uniform triangular mesh on $\Omega$
with $h$ as the maximum size of triangles $K$. Let $\mathcal{E}_{h}$ be the
set of edges $e$ in $\mathcal{T}_{h}$. For $k\geq 1$, the virtual element
space is defined as

(1.3) $\displaystyle\tilde{V}_{h}=\\{v\in
H^{1}_{0}(\Omega):\tilde{v}|_{\partial K}\in\mathbb{B}_{k}(\partial
K),\Delta\tilde{v}|_{K}\in P_{k-2}(K)\\},$

where $P_{-1}=\\{0\\}$ and $\mathbb{B}_{k}(\partial K)=\\{v\in C^{0}(\partial
K):v|_{e}\in P_{k}(e)\ \forall e\subset K\\}$. In computation, the standard
interpolated virtual finite element space on $\mathcal{T}_{h}$ is defined by

(1.4) $\displaystyle V_{h}=\\{v_{h}=\Pi_{h}^{\nabla}\tilde{v}\ :\
v_{h}|_{K}\in\mathbb{V}_{k}(K),\ K\in\mathcal{T}_{h};\
\tilde{v}\in\tilde{V}_{h}\\},$

where $\mathbb{V}_{k}(K)=P_{k}(K)$ for the standard virtual elements (and to
be defined below in (1.7) for the new method), and
$v_{h}=\Pi_{h}^{\nabla}\tilde{v}$ satisfies, for all
$w_{h}\in\mathbb{V}_{k}(K)$,

(1.5) $\displaystyle\langle v_{h}-\tilde{v},w_{h}\rangle_{\partial K}=0,\
\hbox{\rm{ and }}\ (\nabla(v_{h}-\tilde{v}),\nabla w_{h})_{K}=0.$

The stabilizer-free virtual element equation reads: Find
$u_{h}=\Pi_{h}^{\nabla}\tilde{u}\in V_{h}$ such that

(1.6) $\displaystyle(\nabla u_{h},\nabla
v_{h})_{h}=(f,v_{h})\quad\forall\tilde{v}\in\tilde{V}_{h},\
v_{h}=\Pi_{h}^{\nabla}\tilde{v},$

where $(\nabla u_{h},\nabla v_{h})_{h}=\sum_{K\in\mathcal{T}_{h}}(\nabla
u_{h},\nabla v_{h})_{K}$. But the dimension of $V_{h}$ is less than that of
$\tilde{V}_{h}$ unless $k=1$ and on triangular meshes. Thus the bilinear form
in (1.6) is not coercive and the equation does not have a unique solution. A
discrete stabilizer must be added to the equation (1.6).

In order to delete the stabilizer, [6] proposed to replace
$\mathbb{V}_{k}(K)=P_{k}(K)$ by $\mathbb{V}_{k}(K)=P_{k+l}(K)$ in the virtual
element space (1.4) for the case $k=1$ on polygonal meshes, where $l$ depends
on the maximum number of edges. Further numerical tests and comparisons are
given in [7]. Another stabilization-free method for $k=1$ is proposed in [8]
that $\mathbb{V}_{k}(K)=P_{k}(K)\cup H_{l}(K)$, where $H_{l}(K)$ is the set of
2D harmonic polynomials of degree $l$ or less, and $l$ depends on the maximum
number of edges. This is an excellent idea because the $H_{l}$ polynomials may
enforce corerciveness while not destroying the gradient approximation, as they
have vanishing Laplacian. The same idea has been implemented in some other
finite elements [1, 28, 29]. But the method [8] is shown not working for $k>3$
numerically, in the numerical comparison with the new method proposed in this
paper.

We propose to define the interpolation space $\mathbb{V}_{k}(K)$ in (1.4) as

(1.7) $\displaystyle\mathbb{V}_{k}(K)=\\{v_{h}\in C(K)\ :\ v_{h}|_{K_{i}}\in
P_{k}(K_{i}),\ K=\cup_{i=1}^{3}K_{i}\\},$

where $K$ is split in to three triangles by connecting its barycenter with
three vertices, cf. Figure 1. We call $K$ a Hsieh-Clough-Tocher macro-triangle
[14, 29, 30, 57, 58].

$K:$$K_{2}$$K_{1}$$K_{3}$ Figure 1. A Hsieh-Clough-Tocher macro-triangle
$K=\cup_{i=1}^{3}K_{i}$

The interpolation operator $\Pi^{\nabla}_{h}$ in (1.5) is naturally defined by
$v_{h}=\Pi^{\nabla}_{h}\tilde{v}\in\mathbb{V}_{k}(K)$ of (1.7) satisfying

(1.8) $\displaystyle\begin{aligned} v_{h}&=\tilde{v}&&\hbox{\rm{on }}\
\partial K,\\\ (\nabla v_{h},\nabla w_{h})_{K}&=(\nabla\tilde{v},\nabla
w_{h})_{K}&&\forall w_{h}\in H^{1}_{0}(K)\cap\mathbb{V}_{k}(K).\end{aligned}$

We note that a different interpolation space only changes the numerical
quadrature formula for computing $(\nabla u_{h},\nabla
v_{h})=(\nabla\Pi_{h}^{\nabla}\tilde{u},\nabla\Pi_{h}^{\nabla}\tilde{v})$ in
the virtual elements equation. An accurate calculation of local interpolation
does not increase the computational cost of a problem after the stiffness
matrix is generated.

Eliminating a stabilizer would not only reduce computational cost, but also
likely to improve the condition number of the resulting system. More
importantly, a stabilizer-free method would likely utilize fully every degree
of freedom in the discrete approximation. Thus it often leads to discoveries
of superconvergence. The stabilizer is eliminated first in the weak Galerkin
finite element method [2, 18, 19, 23, 31, 32, 37, 38, 39, 42, 43, 44], then in
the $H(\hbox{\rm{div}})$ finite element method [24, 41, 46], in the $C^{0}$ or
$C^{-1}$ finite element methods for the biharmonic equation [40, 48, 49] and
in the discontinuous Galerkin finite element method [17, 26, 35, 36, 40]. It
leads to two-order superconvergent WG finite elements [3, 33, 34, 47, 53, 56]
and two-order superconvergent DG finite elements [50, 51, 54] for second order
elliptic equations, one or two-order superconvergent WG finite elements for
the Stokes equations [25, 45, 52], four-order superconvergent WG finite
elements [55] and four-order superconvergent DG finite elements [51, 56] for
the biharmonic equation. That is, for an example, a $P_{3}$ discontinuous
finite element method, with the stabilizer-free technique, produces the same
order approximate solution as a $C^{1}$-$P_{7}$ finite element method does, in
solving a 2D biharmonic equation.

In this paper, we show that with the new interpolation (1.8), the stabilizer-
free virtual element equation (1.6) has a unique and quasi-optimal solution,
on triangular meshes. Numerical tests on the new stabilizer-free virtual
elements are performed. Numerical comparisons are presented, with the other
stabilizer-free virtual elements and with the standard virtual elements.

##  2\. The well-posedness

We show in this section that the stabilizer-free virtual element equation has
a unique solution.

######  Lemma 2.1.

The interpolation operator $\Pi_{h}^{\nabla}$ is well defined in (1.8) and it
preserves $P_{k}$ polynomials,

(2.1) $\displaystyle\Pi_{h}^{\nabla}\tilde{v}=\tilde{v}\quad\ \hbox{\rm{if }}\
\tilde{v}\in P_{k}(K).$

###### Proof.

Because $\tilde{v}|_{\partial K}\in\mathbb{B}_{k}(\partial K)$, $v_{h}$ can
assume the boundary condition $v_{h}=\tilde{v}$ exactly on $\partial K$. The
linear system of equations in (1.8) is a finite dimensional square system. The
existence is implied by the uniqueness. To show the uniqueness, we let
$\tilde{v}=0$ in (1.8). Letting $w_{h}=v_{h}$ in (1.8), we get

$\displaystyle\nabla v_{h}={\bf 0}\quad\hbox{\rm{ on }}\ K.$

Thus $v_{h}=c$ is a constant on $K$. As $v_{h}$ is continuous on edges,
$v_{h}=c$ is a global constant on the whole domain. By the boundary condition,
we get $0=\tilde{v}|_{\partial\Omega}=v_{h}|_{\partial\Omega}=c$. Hence
$v_{h}=0$ and (1.8) has a unique solution.

If $\tilde{v}\in P_{k}(K)\subset\mathbb{V}_{k}(K)$, defined in (1.4), then the
solution of (1.8) says, letting $w_{h}=v_{h}-\tilde{v}$,

$\displaystyle\nabla(v_{h}-\tilde{v})={\bf 0}.$

Thus $v_{h}-\tilde{v}$ is a global constant which must be zero as it vanishes
at all $\partial K$. (2.1) is proved. ∎

######  Lemma 2.2.

The stabilizer-free virtual element equation (1.6) has a unique solution,
where the interpolation $\Pi_{h}^{\nabla}$ is defined in (1.8).

###### Proof.

As both $\tilde{u},\tilde{v}\in\tilde{V}_{h}$, (1.6) is a finite square system
of linear equations. The uniqueness of solution implies the existence. To show
the uniqueness, we let $f=0$ and $\tilde{v}=\tilde{u}$ in (1.6). It follows
that

$\displaystyle|\Pi_{h}^{\nabla}\tilde{u}|_{1,h}=0.$

Thus $\Pi_{h}^{\nabla}\tilde{u}=c$ is constant on each $K$. But
$\Pi_{h}^{\nabla}\tilde{u}$ is continuous on the whole domain. By the boundary
condition, we get $0=\Pi_{h}^{\nabla}\tilde{u}|_{\partial\Omega}=c$. That is,

(2.2) $\displaystyle\Pi_{h}^{\nabla}\tilde{u}=0.$

On one triangle $K=\cup_{i=1}^{3}K_{i}$,

(2.3) $\displaystyle\tilde{u}=\Pi_{h}^{\nabla}\tilde{u}=0\quad\ \hbox{\rm{on
}}\ \partial K.$

Inside the triangle, by (1.8), (2.2), (2.3) and integration by parts, we have

(2.4) $\displaystyle(-\Delta\tilde{u},w_{h})=(\nabla\tilde{u},\nabla
w_{h})=0\quad\forall w_{h}\in H^{1}_{0}(K)\cap\mathbb{V}_{k}(K).$

By the space $\tilde{V}_{h}$ definition (1.3), we denote

(2.5) $\displaystyle p_{k-2}=-\Delta\tilde{u}\in P_{k-2}(K).$

Let the $w_{h}$ in (2.4) be

(2.6) $\displaystyle w_{h}=p_{k-2}\phi_{{\bf x}_{0}}({\bf x})\in
H^{1}_{0}(K)\cap\mathbb{V}_{k}(K),$

where $\phi_{{\bf x}_{0}}({\bf x})$ is the $P_{1}$ Lagrange basis function at
node ${\bf x}_{0}$, cf. Figure 2. That is, $\phi_{{\bf x}_{0}}({\bf
x})|_{K_{i}}=\lambda_{i,0}$ is the barycentric coordinate at $x_{0}$ on
triangle $K_{i}$, i.e., a linear function which assumes value 1 at ${\bf
x}_{0}$ and vanishing on the line $e_{i}$, cf. Figure 2.

$K:$$K_{2}$$K_{1}$$K_{3}$${\bf x}_{0}$$e_{2}$$e_{1}$$e_{3}$ Figure 2. A Hsieh-
Clough-Tocher macro-triangle $K=\cup_{i=1}^{3}K_{i}$

With the $w_{k}$ in (2.6), we get from (2.4) and (2.5) that

$\displaystyle\int_{K}p_{k-2}^{2}\phi_{{\bf x}_{0}}({\bf x})d{\bf x}=0.$

As $\phi_{{\bf x}_{0}}({\bf x})>0$ inside $K$, it follows that

$\displaystyle p_{k-2}^{2}=0\ \hbox{\rm{ and }}\ p_{k-2}=0\ \hbox{\rm{ on }}\
K.$

By (2.3) and (2.5), $\Delta\tilde{u}=0$ in $K$ and $\tilde{u}=0$ on $\partial
K$. Thus, by the unique solution of the Laplace equation, $\tilde{u}=0$. The
lemma is proved. ∎

##  3\. Convergence

We prove the optimal order convergence of the stabilizer-free virtual element
solutions in this section.

######  Theorem 3.1.

Let $u\in H^{k+1}\cap H^{1}_{0}(\Omega)$ be the exact solution of (1.2). Let
$u_{h}$ be the stabilizer-free virtual element solution of (1.6). It holds
that

(3.1) $\displaystyle|u-u_{h}|_{1}\leq Ch^{k}|u|_{k+1}.$

###### Proof.

Since $w_{h}\in V_{h}\subset H^{1}_{0}(\Omega)$, we subtract (1.6) from (1.2)
to get

$\displaystyle(\nabla(u-u_{h}),\nabla w_{h})_{h}=0\quad\forall w_{h}\in
V_{h}.$

Applying the Schwarz inequality, it follows that

$\displaystyle|u-u_{h}|_{1}^{2}$
$\displaystyle=(\nabla(u-u_{h}),\nabla(u-I_{h}u))$
$\displaystyle\leq|u-u_{h}|_{1}|u-I_{h}u|_{1}\leq
Ch^{k}|u|_{k+1}|u-u_{h}|_{1},$

where $I_{h}u$ is the Scott-Zhang interpolation on quasi-uniform triangulation
$\mathcal{T}_{h}$ [27]. The proof is complete. ∎

To get the optimal order $L^{2}$ error bound, we need a full regularity of the
dual equation that the solution of the equation,

$\displaystyle-\Delta w$ $\displaystyle=g\quad\hbox{\rm{ in }}\ \Omega,$
$\displaystyle w$ $\displaystyle=0\quad\hbox{\rm{ on }}\ \partial\Omega,$

satisfies

(3.2) $\displaystyle|w|_{2}\leq C\|g\|_{0}.$

######  Theorem 3.2.

Let $u\in H^{k+1}\cap H^{1}_{0}(\Omega)$ be the exact solution of (1.2). Let
$u_{h}$ be the stabilizer-free virtual element solution of (1.6). Assuming
(3.2), it holds that

$\displaystyle\|u-u_{h}\|_{0}\leq Ch^{k+1}|u|_{k+1}.$

###### Proof.

Let $w\in H^{2}(\Omega)\cap H^{1}_{0}(\Omega)$ be the dual solution,

(3.3) $\displaystyle(\nabla w,\nabla v)$ $\displaystyle=(u-u_{h},v),\ \forall
v\in H^{1}_{0}(\Omega).$

Thus, by (3.3), (3.2) and (3.1), we get

$\displaystyle\|u-u_{h}\|_{0}^{2}$ $\displaystyle=(\nabla
w,\nabla(u-u_{h}))=(\nabla(w-w_{h}),\nabla(u-u_{h}))$ $\displaystyle\leq
Ch|w|_{2}h^{k}|u|_{k+1}\leq Ch^{k+1}|u|_{k+1}\|u-u_{h}\|_{0},$

where $w_{h}$ is the virtual element solution to the equation (3.3). We obtain
the $L^{2}$ error bound. ∎

##  4\. Numerical test

We solve numerically the Poisson equation (1.1) on domain
$\Omega=(0,1)\times(0,1)$, where an exact solution is chosen as

(4.1) $\displaystyle u(x,y)=\sin(\pi x)\sin(\pi y).$

Grid 1:Grid 2: Figure 3. The first two levels of grids in the computation of
Tables 1–4.

The computation is done first on a family of slightly irregular triangular
meshes shown in Figure 3.

In Table 1, we first list the errors and the computed order of convergence for
the $P_{1}$ and $P_{2}$ stabilizer-free virtual elements (1.4) with the Hsieh-
Clough-Tocher macro-triangle interpolation space $\mathbb{V}_{k}$ in (1.7).
Optimal orders are achieved for both elements and in both $L^{2}$ and $H^{1}$
norms. Here we use $\Pi_{h}^{\nabla}u$ instead of $u$ to check the error so
that we can detect possible superconvergence. At the bottom of Table 1, we
test a method of [8] where the interpolation space is $P_{2}\cup H_{3}$, i.e.,
enriching the $P_{2}(K)$ space by two harmonic $P_{3}$ polynomials. The method
is not proved yet. The numerical test shows it works well, producing optimal
order errors.

Table 1. Error profile on Figure 3 meshes for (4.1).  Grid | $\|\Pi^{\nabla}_{h}u-u_{h}\|_{0}$ | $O(h^{r})$ | $|\Pi^{\nabla}_{h}u-u_{h}|_{1}$ | $O(h^{r})$
---|---|---|---|---
| By the $P_{1}$ virtual element with HTC interpolation.
7 | 0.8208E-04 | 2.00 | 0.7621E-02 | 1.00
8 | 0.2052E-04 | 2.00 | 0.3817E-02 | 1.00
9 | 0.5129E-05 | 2.00 | 0.1911E-02 | 1.00
| By the $P_{2}$ virtual element with HTC interpolation.
7 | 0.7537E-07 | 3.00 | 0.6668E-04 | 1.99
8 | 0.9423E-08 | 3.00 | 0.1670E-04 | 2.00
9 | 0.1178E-08 | 3.00 | 0.4179E-05 | 2.00
| By the $P_{2}$ virtual element with 2 h.p. [8].
7 | 0.8627E-07 | 3.00 | 0.7861E-04 | 2.00
8 | 0.1079E-07 | 3.00 | 0.1968E-04 | 2.00
9 | 0.1350E-08 | 3.00 | 0.4925E-05 | 2.00

In Table 2, we first list the errors and the computed order of convergence for
the $P_{3}$ stabilizer-free virtual elements, $k=3$ in (1.4) with the Hsieh-
Clough-Tocher macro-triangle interpolation space $\mathbb{V}_{k}$ in (1.7).
Optimal orders are achieved for the element in both $L^{2}$ and $H^{1}$ norms.
At the bottom of Table 2, we test a method of [8] where the interpolation
space is $P_{3}\cup H_{5}$, i.e., enriching the $P_{3}(K)$ polynomial space by
four harmonic $P_{4}$ and $P_{5}$ polynomials. The method is not proved yet.
The numerical test shows it works well, producing optimal order errors.

Table 2. Error profile on Figure 3 meshes for (4.1).  Grid | $\|\Pi^{\nabla}_{h}u-u_{h}\|_{0}$ | $O(h^{r})$ | $|\Pi^{\nabla}_{h}u-u_{h}|_{1}$ | $O(h^{r})$
---|---|---|---|---
| By the $P_{3}$ virtual element with HTC interpolation.
6 | 0.4664E-08 | 4.00 | 0.3068E-05 | 2.99
7 | 0.2916E-09 | 4.00 | 0.3844E-06 | 3.00
8 | 0.1824E-10 | 4.00 | 0.4811E-07 | 3.00
| By the $P_{3}$ virtual element with 8 h.p. [8].
6 | 0.5676E-08 | 4.01 | 0.3392E-05 | 2.99
7 | 0.3542E-09 | 4.00 | 0.4257E-06 | 2.99
8 | 0.2222E-10 | 3.99 | 0.5372E-07 | 2.99

In Table 3, we first list the errors and the computed orders of convergence
for the $P_{4}$ stabilizer-free virtual elements, $k=4$ in (1.4) with the
Hsieh-Clough-Tocher macro-triangle interpolation space $\mathbb{V}_{k}$ in
(1.7). Optimal orders are achieved for the element in both $L^{2}$ and $H^{1}$
norms. Here the computer accuracy is exhausted when computing the last grid
solution. At the bottom of Table 3, we test a method of [8] where the
interpolation space is $P_{4}\cup H_{10}$, i.e., enriching the $P_{4}(K)$
polynomial space by 12 harmonic $P_{5}$, $P-6$, $P_{7}$, $P_{8}$, $P_{9}$ and
$P_{10}$ polynomials. Since the method does not work, we tested by adding more
harmonic polynomials until the error can not be reduced anymore.

Table 3. Error profile on Figure 3 meshes for (4.1).  Grid | $\|\Pi^{\nabla}_{h}u-u_{h}\|_{0}$ | $O(h^{r})$ | $|\Pi^{\nabla}_{h}u-u_{h}|_{1}$ | $O(h^{r})$
---|---|---|---|---
| By the $P_{4}$ virtual element with HTC interpolation.
5 | 0.7612E-09 | 5.00 | 0.3811E-06 | 3.99
6 | 0.2374E-10 | 5.00 | 0.2389E-07 | 4.00
7 | 0.7616E-12 | — | 0.1495E-08 | 4.00
| By the $P_{4}$ virtual element with 12 h.p. [8].
4 | 0.4589E-07 | 5.01 | 0.7336E-05 | 3.94
5 | 0.1525E-08 | 4.91 | 0.5830E-06 | 3.65
6 | 0.9493E-10 | 4.01 | 0.9239E-07 | 2.66

In Table 4, we first list the errors and the computed orders of convergence
for the $P_{5}$ stabilizer-free virtual elements, $k=5$ in (1.4) with the
Hsieh-Clough-Tocher macro-triangle interpolation space $\mathbb{V}_{k}$ in
(1.7). Optimal orders are achieved for the element in both $L^{2}$ and $H^{1}$
norms. At the middle of Table 4, we test a method of [8] where the
interpolation space is $P_{5}\cup H_{10}$, i.e., enriching the $P_{5}(K)$
polynomial space by 10 harmonic $P_{6}$, $P_{7}$, $P_{8}$, $P_{9}$ and
$P_{10}$ polynomials. Since the method does not work, we tested by adding more
harmonic polynomials until the error can not be reduced anymore. Comparing to
the last case, the convergent order deteriorates a lot. At the bottom of Table
4, we list the errors and the computed orders of convergence for the $P_{6}$
stabilizer-free virtual elements, $k=6$ in (1.4) with the Hsieh-Clough-Tocher
macro-triangle interpolation space $\mathbb{V}_{k}$ in (1.7). Optimal orders
are achieved for the element in both $L^{2}$ and $H^{1}$ norms.

Table 4. Error profile on Figure 3 meshes for (4.1).  Grid | $\|\Pi^{\nabla}_{h}u-u_{h}\|_{0}$ | $O(h^{r})$ | $|\Pi^{\nabla}_{h}u-u_{h}|_{1}$ | $O(h^{r})$
---|---|---|---|---
| By the $P_{5}$ virtual element with HTC interpolation.
3 | 0.3339E-07 | 6.01 | 0.5275E-05 | 4.96
4 | 0.5170E-09 | 6.01 | 0.1665E-06 | 4.99
5 | 0.8033E-11 | 6.01 | 0.5223E-08 | 4.99
| By the $P_{5}$ virtual element with 10 h.p. [8].
2 | 0.4262E-05 | 5.18 | 0.2158E-03 | 4.24
3 | 0.1566E-06 | 4.77 | 0.2319E-04 | 3.22
4 | 0.1819E-07 | 3.11 | 0.5702E-05 | 2.02
| By the $P_{6}$ virtual element with HTC interpolation.
2 | 0.1696E-06 | 7.30 | 0.1576E-04 | 6.19
3 | 0.1320E-08 | 7.01 | 0.2515E-06 | 5.97
4 | 0.1025E-10 | 7.01 | 0.3956E-08 | 5.99

The second part of computation is done on the uniform triangular meshes, shown
as in Figure 4. This is mainly for detecting possible superconvergence.

Grid 1:Grid 2:Grid 3: Figure 4. The first three levels of grids for the
computation in Tables 5–7.

In Table 5, we first list the errors and the computed orders of convergence
for the $P_{1}$ stabilizer-free virtual elements, $k=1$ in (1.4) with the
Hsieh-Clough-Tocher macro-triangle interpolation space $\mathbb{V}_{k}$ in
(1.7). Optimal orders are achieved for the element in both $L^{2}$ and $H^{1}$
norms. In fact, we have one-order superconvergence in $H^{1}$ semi-norm. At
the bottom of Table 5, we list the errors and the computed orders of
convergence for the $P_{1}$ standard virtual elements, $k=1$ in (1.4). Optimal
orders are achieved for the element in both $L^{2}$ and $H^{1}$ norms. Again,
we have one-order $H^{1}$ superconvergence for this element. Comparing the
errors, the new method is slightly better which is understandable as their
interpolation space $P_{1}(K)$ is a subspace of our interpolation space
$\mathbb{V}_{1}=C(K)\cap\cup_{i=1}^{3}P_{1}(K_{i})$.

Table 5. Error profile on Figure 4 meshes for (4.1).  Grid | $\|\Pi^{\nabla}_{h}u-u_{h}\|_{0}$ | $O(h^{r})$ | $|\Pi^{\nabla}_{h}u-u_{h}|_{1}$ | $O(h^{r})$
---|---|---|---|---
| By the $P_{1}$ SF virtual element with HTC interpolation.
7 | 0.3032E-03 | 2.00 | 0.1382E-02 | 2.00
8 | 0.7586E-04 | 2.00 | 0.3457E-03 | 2.00
9 | 0.1897E-04 | 2.00 | 0.8644E-04 | 2.00
| By the standard $P_{1}$ virtual element [4].
6 | 0.1210E-02 | 1.98 | 0.5518E-02 | 1.99
7 | 0.3032E-03 | 2.00 | 0.1382E-02 | 2.00
8 | 0.7586E-04 | 2.00 | 0.3457E-03 | 2.00

In Table 6, we first list the errors and the computed orders of convergence
for the $P_{2}$ stabilizer-free virtual elements, $k=2$ in (1.4) with the
Hsieh-Clough-Tocher macro-triangle interpolation space $\mathbb{V}_{k}$ in
(1.7). Optimal orders are achieved for the element in both $L^{2}$ and $H^{1}$
norms. Unlike the traditional $P_{2}$ finite element, we do not have any
superconvergence. Comparing to the traditional $P_{2}$ finite element, we
compute a solution in a larger vector space but get a worse result. This is
because the added Hsieh-Clough-Tocher macro-bubbles destroy the symmetry of
$P_{2}$ finite element equations on uniform triangular meshes. At the bottom
of Table 6, we list the errors and the computed orders of convergence for the
standard $P_{2}$ virtual elements, $k=2$ in (1.4). Optimal orders are achieved
for the element in both $L^{2}$ and $H^{1}$ norms. Comparing the two errors,
the new method is much better which is understandable as there is no
stabilizer here.

Table 6. Error profile on Figure 4 meshes for (4.1).  Grid | $\|\Pi^{\nabla}_{h}u-u_{h}\|_{0}$ | $O(h^{r})$ | $|\Pi^{\nabla}_{h}u-u_{h}|_{1}$ | $O(h^{r})$
---|---|---|---|---
| By the $P_{2}$ SF virtual element with HTC interpolation.
7 | 0.5604E-07 | 3.72 | 0.2302E-04 | 2.48
8 | 0.5261E-08 | 3.41 | 0.4990E-05 | 2.21
9 | 0.5901E-09 | 3.16 | 0.1194E-05 | 2.06
| By the standard $P_{2}$ virtual element [4].
6 | 0.1885E-05 | 3.15 | 0.3225E-03 | 2.08
7 | 0.2285E-06 | 3.04 | 0.7958E-04 | 2.02
8 | 0.2833E-07 | 3.01 | 0.1983E-04 | 2.00

In Table 7, we list the errors and the computed orders of convergence for the
$P_{3}$, $P_{4}$, $P_{5}$ and $P_{6}$ stabilizer-free virtual elements,
$k=3,4,5,$ or $6$ in (1.4) with the Hsieh-Clough-Tocher macro-triangle
interpolation space $\mathbb{V}_{k}$ in (1.7). Optimal orders are achieved for
all the elements in both $L^{2}$ and $H^{1}$ norms. Comparing the errors of
same virtual elements, the uniform triangular meshes are much better than the
triangular meshes shown in Figure 3.

Table 7. Error profile on Figure 3 meshes for (4.1).  Grid | $\|\Pi^{\nabla}_{h}u-u_{h}\|_{0}$ | $O(h^{r})$ | $|\Pi^{\nabla}_{h}u-u_{h}|_{1}$ | $O(h^{r})$
---|---|---|---|---
| By the $P_{3}$ SF virtual element with HTC interpolation.
6 | 0.4906E-07 | 3.98 | 0.1628E-04 | 2.98
7 | 0.3087E-08 | 3.99 | 0.2048E-05 | 2.99
8 | 0.1935E-09 | 4.00 | 0.2567E-06 | 3.00
| By the $P_{4}$ SF virtual element with HTC interpolation.
5 | 0.1038E-07 | 5.03 | 0.3216E-05 | 4.00
6 | 0.3216E-09 | 5.01 | 0.2011E-06 | 4.00
7 | 0.1001E-10 | 5.01 | 0.1258E-07 | 4.00
| By the $P_{5}$ SF virtual element with HTC interpolation.
5 | 0.2035E-09 | 6.00 | 0.7749E-07 | 4.98
6 | 0.3173E-11 | 6.00 | 0.2433E-08 | 4.99
7 | 0.5018E-13 | 5.98 | 0.7612E-10 | 5.00
| By the $P_{6}$ SF virtual element with HTC interpolation.
3 | 0.5536E-07 | 6.99 | 0.5949E-05 | 5.94
4 | 0.4294E-09 | 7.01 | 0.9376E-07 | 5.99
5 | 0.3333E-11 | 7.01 | 0.1468E-08 | 6.00

##  5\. Ethical Statement

###  5.1. Compliance with Ethical Standards

The submitted work is original and is not published elsewhere in any form or
language.

###  5.2. Funding

Xuejun Xu was supported by National Natural Science Foundation of China (Grant
Nos. 12071350), Shanghai Municipal Science and Technology Major Project No.
2021SHZDZX0100, and Science and Technology Commission of Shanghai
Municipality.

###  5.3. Conflict of Interest

There is no potential conflict of interest .

###  5.4. Ethical approval

This article does not contain any studies involving animals. This article does
not contain any studies involving human participants.

###  5.5. Informed consent

This research does not have any human participant.

###  5.6. Availability of supporting data

This research does not use any external or author-collected data.

###  5.7. Authors’ contributions

All authors made equal contribution.

###  5.8. Acknowledgments

None.

## References

  * [1] A. Al-Taweel, Y. Dong, S. Hussain and X. Wang, A weak Galerkin harmonic finite element method for Laplace equation, Commun. Appl. Math. Comput. 3 (2021), no. 3, 527–543.
  * [2] A. Al-Twaeel, S. Hussian and X. Wang, A stabilizer free weak Galerkin finite element method for parabolic equation, J. Comput. Appl. Math., 392 (2021), 113373.
  * [3] A. Al-Twaeel, X. Wang, X. Ye and S. Zhang, A stabilizer free weak Galerkin element method with supercloseness of order two, Numer. Methods Partial Differential Equations 37 (2021), no. 2, 1012–1029.
  * [4] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini and A. Russo, Basic principles of virtual element methods, Math. Models Methods Appl. Sci. 23 (2013) 199–214.
  * [5] L. Beirão da Veiga, F. Brezzi, L. D. Marini and A. Russo, H(div) and H(curl)-conforming virtual element methods, Numer. Math. 133 (2016) 303–332.
  * [6] S. Berrone, A. Borio, F. Marcon and G. Teora, A first-order stabilization-free virtual element method, Appl. Math. Lett. 142 (2023), Paper No. 108641, 6 pp.
  * [7] S. Berrone, A. Borio and F. Marcon, Comparison of standard and stabilization free Virtual Elements on anisotropic elliptic problems, 2022, arXiv:2202.08571v1.
  * [8] S. Berrone, A. Borio and F. Marcon, Lowest order stabilization free Virtual Element Method for the Poisson equation, 2023, arXiv:2103.16896v1.
  * [9] S. Cao, L. Chen and R. Guo, A virtual finite element method for two-dimensional Maxwell interface problems with a background unfitted mesh, Math. Models Methods Appl. Sci. 31 (2021), no. 14, 2907–2936.
  * [10] S. Cao, L. Chen and R. Guo, Immersed virtual element methods for electromagnetic interface problems in three dimensions, Math. Models Methods Appl. Sci. 33 (2023), no. 3, 455–503.
  * [11] L. Chen, H. Wei and M. Wen, An interface-fitted mesh generator and virtual element methods for elliptic interface problems, J. Comput. Phys. 334 (2017), 327–348.
  * [12] L. Chen and J. Huang, Some error analysis on virtual element methods, Calcolo 55 (2018), no. 1, Paper No. 5, 23 pp.
  * [13] L. Chen and X. Huang, Nonconforming virtual element method for 2mth order partial differential equations in Rn, Math. Comp. 89 (2020), no. 324, 1711–1744.
  * [14] R.W. Clough and J.L. Tocher, Finite element stiffness matrices for analysis of plates in bending, in: Proceedings of the Conference on Matrix Methods in Structural Mechanics, Wright Patterson A.F.B. Ohio, 1965.
  * [15] F. Feng, W. Han and J. Huang, Virtual element methods for elliptic variational inequalities of the second kind, J. Sci. Comput. 80 (2019), no. 1, 60–80.
  * [16] F. Feng, J. Huang and Y. Yu, A non-consistent virtual element method for reaction diffusion equations, East Asian J. Appl. Math. 10 (2020), no. 4, 786–799.
  * [17] Y. Feng, Y. Liu, R. Wang and S. Zhang, A conforming discontinuous Galerkin finite element method on rectangular partitions, Electron. Res. Arch. 29 (2021), no. 3, 2375–2389.
  * [18] Y. Feng, Y. Liu, R. Wang and S. Zhang, A stabilizer-free weak Galerkin finite element method for the Stokes equations, Adv. Appl. Math. Mech. 14 (2022), no. 1, 181–201.
  * [19] F. Gao, X. Ye and S. Zhang, A discontinuous Galerkin finite element method without interior penalty terms, Adv. Appl. Math. Mech. 14 (2022), no. 2, 299–314.
  * [20] J. Huang and S. Lin, A C0P2 time-stepping virtual element method for linear wave equations on polygonal meshes, Electron. Res. Arch. 28 (2020), no. 2, 911–933.
  * [21] J. Huang and Y. Yu, Some estimates for virtual element methods in three dimensions, Comput. Methods Appl. Math. 23 (2023), no. 1, 177–187.
  * [22] J. Huang and Y. Yu, A medius error analysis for nonconforming virtual element methods for Poisson and biharmonic equations, J. Comput. Appl. Math. 386 (2021), Paper No. 113229, 20 pp.
  * [23] L. Mu, X. Ye and S. Zhang, A stabilizer free, pressure robust and superconvergence weak Galerkin finite element method for the Stokes Equations on polytopal mesh, SIAM J. Sci, Comput., 43 (2021), A2614–A2637.
  * [24] L. Mu, X. Ye and S. Zhang, Development of pressure-robust discontinuous Galerkin finite element methods for the Stokes problem, J. Sci. Comput. 89 (2021), no. 1, Paper No. 26, 25 pp.
  * [25] L. Mu, X. Ye and S. Zhang, A stabilizer free, pressure robust, and superconvergence weak Galerkin finite element method for the Stokes Equations on polytopal mesh, SIAM J. Sci. Comput. 43 (2021), no. 4, A2614–A2637.
  * [26] L. Mu, X. Ye, S. Zhang and P. Zhu, A DG method for the Stokes equations on tensor product meshes with $[P_{k}]^{d}-P_{k-1}$ element, Communications on Applied Mathematics and Computation, 2023, doi:10.1007/s42967-022-00243-9.
  * [27] L. R. Scott and S. Zhang, Finite element interpolation of nonsmooth functions satisfying boundary conditions, Math. Comp. 54 (1990), no. 190, 483–493.
  * [28] T. Sorokina and S. Zhang, Conforming and nonconforming harmonic finite elements, Appl. Anal. 99 (2020), no. 4, 569–584.
  * [29] T. Sorokina and S. Zhang, Conforming harmonic finite elements on the Hsieh-Clough-Tocher split of a triangle, Int. J. Numer. Anal. Model. 17 (2020), no. 1, 54–67.
  * [30] X. Xu and S. Zhang, A new divergence-free interpolation operator with applications to the Darcy-Stokes-Brinkman equations, SIAM J. Sci. Comput. 32 (2010), no. 2, 855–874.
  * [31] J. Wang, X. Ye and S. Zhang, Numerical investigation on weak Galerkin finite elements, Int. J. Numer. Anal. Model. 17 (2020) no. 4, 517–531.
  * [32] J. Wang, X. Ye and S. Zhang, Weak Galerkin finite element methods with or without stabilizers, Numer. Algorithms 88 (2021), no. 3, 1361–1381.
  * [33] J. Wang, X. Wang, X. Ye, S. Zhang and Peng Zhu, Two-order superconvergence for a weak Galerkin method on rectangular and cuboid grids, Numer. Methods Partial Differential Equations 39 (2023), no. 1, 744–758.
  * [34] J. Wang, X. Wang, X. Ye, S. Zhang and Peng Zhu, On the superconvergence of a WG method for elliptic problem with variable coefficients, Science China, 2023, doi:10.1007/s11425-022-2097-8.
  * [35] X. Ye and S. Zhang, A conforming discontinuous Galerkin finite element method, Int. J. Numer. Anal. Model. 17 (2020) no. 1, 110–117.
  * [36] X. Ye and S. Zhang, A conforming discontinuous Galerkin finite element method: Part II, Int. J. Numer. Anal. Model. 17 (2020), no. 2, 281–296.
  * [37] X. Ye and S. Zhang, A stabilizer-free weak Galerkin finite element method on polytopal meshes. J. Comput. Appl. Math. 371 (2020), 112699, 9 pp.
  * [38] X. Ye and S. Zhang, A Stabilizer Free Weak Galerkin Method for the Biharmonic Equation on Polytopal Meshes, SIAM J. Numer. Anal. 58 (2020), no. 5, 2572–2588.
  * [39] X. Ye, S. Zhang and Y. Zhu, Stabilizer-free weak Galerkin methods for monotone quasilinear elliptic PDEs, Results Appl. Math. 8 (2020), Paper No. 100097, 10 pp.
  * [40] X. Ye, S. Zhang and Y. Zhu, A conforming discontinuous Galerkin finite element method: Part III, Int. J. Numer. Anal. Model. 17 (2020), no. 6, 794–805.
  * [41] X. Ye and S. Zhang, A stabilizer-free pressure-robust finite element method for the Stokes equations, Adv. Comput. Math. 47 (2021), no. 2, Paper No. 28, 17 pp.
  * [42] X. Ye and S. Zhang, A stabilizer free weak Galerkin finite element method on polytopal mesh: Part II, J. Comput. Appl. Math. 394 (2021), Paper No. 113525, 11 pp.
  * [43] X. Ye and S. Zhang, A stabilizer free weak Galerkin finite element method on polytopal mesh: Part III, J. Comput. Appl. Math. 394 (2021), Paper No. 113538, 9 pp.
  * [44] X. Ye and S. Zhang, A new weak gradient for the stabilizer free weak Galerkin method with polynomial reduction, Discrete Contin. Dyn. Syst. Ser. B 26 (2021), no. 8, 4131–4145.
  * [45] X. Ye and S. Zhang, A stabilizer free WG Method for the Stokes Equations with order two superconvergence on polytopal mesh, Electron. Res. Arch. 29 (2021), no. 6, 3609–3627.
  * [46] X. Ye and S. Zhang, A numerical scheme with divergence free H-div triangular finite element for the Stokes equations, Appl. Numer. Math. 167 (2021), 211–217.
  * [47] X. Ye and S. Zhang, A $P_{k+2}$ polynomial lifting operator on polygons and polyhedrons, Appl. Math. Lett. 116 (2021), Paper No. 107033, 6 pp.
  * [48] X. Ye and S. Zhang, A weak divergence CDG method for the biharmonic equation on triangular and tetrahedral meshes, Appl. Numer. Math. 178 (2022), 155–165.
  * [49] X. Ye and S. Zhang, A $C^{0}$-conforming DG finite element method for biharmonic equations on triangle/tetrahedron, J. Numer. Math. 30 (2022), no. 3, 163–172.
  * [50] X. Ye and S. Zhang, Achieving superconvergence by one-dimensional discontinuous finite elements: the CDG method, East Asian J. Appl. Math. 12 (2022), no. 4, 781–790.
  * [51] X. Ye and S. Zhang, Order two superconvergence of the CDG method for the Stokes equations on triangle/tetrahedron, Journal of Applied Analysis and Computation, 12 (2022), no. 6, 2578–2592.
  * [52] X. Ye and S. Zhang, Achieving superconvergence by one-dimensional discontinuous finite elements: weak Galerkin method, East Asian J. Appl. Math. 12 (2022), no. 3, 590–598.
  * [53] X. Ye and S. Zhang, Constructing order two superconvergent WG finite elements on rectangular meshes, Numer. Math. Theory Methods Appl. 16 (2023), no. 1, 230–241.
  * [54] X. Ye and S. Zhang, Order two superconvergence of the CDG finite elements on triangular and tetrahedral meshes, CSIAM Trans. Appl. Math. 4 (2023), no. 2, 256–274.
  * [55] X. Ye and S. Zhang, Four-order superconvergent weak Galerkin methods for the biharmonic equation on triangular meshes, Communications on Applied Mathematics and Computation, Doi:10.1007/s42967-022-00201-5, 2023.
  * [56] X. Ye and S. Zhang, Four-order superconvergent CDG finite elements for the biharmonic equation on triangular meshes, J. Comput. Appl. Math., 2023, doi: 10.1016/j.cam.2023.115516.
  * [57] S. Zhang, An optimal order multigrid method for biharmonic, $C^{1}$ finite-element equations, Numer. Math. 56 (1989), 613–624.
  * [58] S. Zhang, A new family of stable mixed finite elements for the 3D Stokes equations. Math. Comp. 74 (2005), no. 250, 543–554.


# Radial velocity study of the CP star $\epsilon$ Ursae Majoris

N. A. Sokolov1,2
1Central Astronomical Observatory at Pulkovo, St. Petersburg 196140, Russia
2Isaac Newton Institute of Chile, Branch at St. Petersburg E-mail:
<EMAIL_ADDRESS>on observations collected with the ELODIE spectrograph
on the 193-cm telescope at the Observatoire de Haute-Provence (CNRS), France

(Accepted 2007 November 16. Received 2007 November 16; in original form 2007
July 27)

###### Abstract

In this Letter, the radial velocity variability of the chemically peculiar
star $\epsilon$ Ursae Majoris ($\epsilon$ UMa) from the sharp cores of the
hydrogen lines is investigated. This study is based on the ELODIE archival
data obtained at different phases of the rotational cycle. The star exhibits
low-amplitude radial velocity variations with a period of $P$ = 5.0887 d. The
best Keplerian solution yields an eccentricity $e$ = 0.503 and a minimum mass
$\sim$14.7$M_{\rm Jup}$ on the hypothesis that the rotational axis of
$\epsilon$ UMa is perpendicular to the orbital plane. This result indicate
that the companion is the brown-dwarf with the projected semi-amplitude
variation of the radial velocity $K_{\rm 2}$ = 135.9 km ${\rm s}^{\rm-1}$ and
the sine of inclination times semi-major axis $a_{2}$$\sin$(i) = 0.055 au.

###### keywords:

binaries: spectroscopic – stars: chemically peculiar – stars: individual:
$\epsilon$ Ursae Majoris.

††pagerange: Radial velocity study of the CP star $\epsilon$ Ursae
Majoris–References††pubyear: 2008

##  1 Introduction

Epsilon Ursae Majoris ($\epsilon$ UMa, HD 112185, HR 4905) is the brightest
($V$ = 1.77 mag) chemically peculiar (CP) star and has been extensively
studied during the last century. Guthnick (1934) established a period of
5.d0887 d from variations in the intensity of the Ca II K line and also
noticed a periodic splitting of some lines. Struve & Hiltner (1943)
subsequently reported doubling of lines of Cr II, Fe II, V II and other
elements at certain phases. Since the overall widths of the double lines are
the same and not all lines double, they ruled out the orbital motion and
instead suggested that the phenomenon is related to rotation of the star.
Provin (1953) measured a double wave light variation with the same 5.d0887
period. The star is brightest when the Ca II K line intensity is near its
minimum and other elements are near their maximum strength.

Perhaps the most interesting and controversial aspect of $\epsilon$ UMa is its
the radial velocity variations. Abt & Snowden (1973) carried out the radial
velocity analysis of the 62 brightest northern CP stars for spectroscopic
binaries. They were used special observing and measuring techniques of the
hydrogen lines and did not find the radial velocity variations for $\epsilon$
UMa. Not surprisingly, they could not to measure the radial velocity
variations because the mean internal error per spectrum of 1.2 km ${\rm
s}^{-1}$ is near the limit at which the photographic technique is useful. On
the other hand, Morgan et al. (1978) detected the duplicity for the system
$\epsilon$ UMa with a separation of 0.053 arcsec, using speckle interferometry
technique. Although, the authors noted that the observed separation
corresponds to the diffraction limit of the 2.5-m Isaac Newton telescope.

Woszczyk & Jasiński (1980) measured the radial velocity variations of many
lines in the star and found sinusoidal variations in the radial velocities of
Fe, Cr and Ti lines with amplitudes of about 20 km ${\rm s}^{-1}$ and
attributed this to the existence of a spots of enhanced abundance of these
elements. Several surface abundance Doppler images of $\epsilon$ UMa have been
produced (Wehlau et al., 1982; Rice & Wehlau, 1990; Hatzes, 1991; Rice et al.,
1997; Holmgren & Rice, 2000). Most of the published Doppler images relate to
Fe, O, Ca, or Cr abundance. Recently, Lueftinger et al. (2003) determined for
$\epsilon$ UMa for the first time the abundance distributions of Mn, Ti, Sr,
and Mg. Attempt to determine orbital motion using these spectral lines would
be impossible. However, the hydrogen lines do not show significant rotational
effect in their radial velocity. They would be preferable for the measurement
of binary motion if the radial velocity of the hydrogen lines could be
measured with precision.

In this Letter, we present the radial velocity measurement of the hydrogen
lines in the spectrum of CP star $\epsilon$ UMa. The ELODIE observations of
this star and data reduction are described in Sect. 2. The orbital solution
derived for the satellite candidate are presented in Sect. 3. Section 4
reports discussion of our results for $\epsilon$ UMa and conclusions are
presented in Sect. 5.

##  2 Observational material and data reduction

Table 1: Measured radial velocities of $\epsilon$ UMa. All data are relative to the Solar system barycentre. File name | Exposure time | MJD | Phase | Radial velocity (km ${\rm s}^{-1}$) |  |
---|---|---|---|---|---|---
| (s) | (2,450,000+) |  | $H_{\delta}$ | $H_{\gamma}$ | $H_{\beta}$ | $H_{\alpha}$ | RVMean | $\sigma_{Mean}$
19960204 0022 | 168.73 | 0117.9814 | 0.821 | -9.56 | -9.19 | -9.93 | -9.47 | -9.54 | 0.31
19970218 0032 | 200.81 | 0497.9778 | 0.496 | -9.06 | -8.72 | -9.67 | -9.27 | -9.18 | 0.40
19970218 0033 | 200.37 | 0497.9819 | 0.497 | -8.90 | -8.77 | -9.51 | -9.19 | -9.09 | 0.33
20000124 0030 | 60.58 | 1568.1285 | 0.795 | -9.60 | -9.07 | -9.73 | -9.48 | -9.47 | 0.29
20000124 0031 | 90.57 | 1568.1440 | 0.798 | -9.50 | -9.16 | -9.71 | -9.63 | -9.50 | 0.24
20000126 0040 | 101.56 | 1570.1447 | 0.191 | -8.11 | -8.24 | -8.99 | -8.94 | -8.57 | 0.46
20000126 0041 | 50.58 | 1570.1499 | 0.192 | -8.29 | -8.28 | -9.07 | -8.99 | -8.66 | 0.43
20000519 0012 | 120.79 | 1683.9584 | 0.557 | -9.04 | -8.82 | -9.17 | -8.99 | -9.01 | 0.14
20000520 0009 | 120.76 | 1684.8706 | 0.737 | -9.56 | -9.02 | -9.32 | -9.13 | -9.26 | 0.24
20000521 0006 | 600.78 | 1685.8694 | 0.933 | -10.14 | -9.49 | -9.71 | -9.36 | -9.68 | 0.34
20000521 0007 | 600.44 | 1685.8784 | 0.935 | -9.85 | -9.44 | -9.63 | -9.44 | -9.59 | 0.20
20000606 0017 | 30.54 | 1701.8490 | 0.073 | -8.05 | -7.89 | -9.04 | -8.55 | -8.38 | 0.52
20000607 0007 | 30.27 | 1702.8632 | 0.272 | -8.44 | -8.33 | -8.73 | -8.44 | -8.48 | 0.17
20000608 0006 | 30.53 | 1703.8467 | 0.466 | -9.09 | -8.80 | -9.24 | -8.87 | -9.00 | 0.20
20000609 0003 | 30.97 | 1704.8443 | 0.662 | -9.34 | -8.64 | -9.35 | -8.68 | -9.00 | 0.40
20000615 1417 | 300.17 | 1710.8469 | 0.841 | -9.51 | -8.94 | -9.76 | -9.08 | -9.32 | 0.38
20000615 1420 | 50.47 | 1710.8505 | 0.842 | -9.46 | -8.74 | -9.79 | -9.30 | -9.32 | 0.44
20000615 1423 | 30.76 | 1710.8608 | 0.844 | -9.61 | -8.63 | -9.77 | -9.05 | -9.27 | 0.52
20000615 1425 | 30.26 | 1710.8634 | 0.844 | -9.68 | -8.84 | -9.67 | -9.17 | -9.34 | 0.41
20000615 1426 | 30.26 | 1710.8661 | 0.845 | -9.67 | -8.99 | -10.04 | -9.13 | -9.46 | 0.49
20000616 1418 | 50.92 | 1711.8591 | 0.040 | -8.35 | -7.88 | -9.11 | -8.85 | -8.55 | 0.55
20000616 1421 | 20.71 | 1711.8630 | 0.041 | -8.17 | -7.81 | -9.17 | -8.94 | -8.52 | 0.64
20000616 1424 | 20.65 | 1711.8670 | 0.042 | -8.13 | -8.00 | -9.33 | -8.94 | -8.60 | 0.64
20000618 0006 | 30.59 | 1713.8460 | 0.431 | -8.52 | -8.45 | -9.04 | -8.80 | -8.70 | 0.27
20000618 0007 | 30.71 | 1713.8506 | 0.432 | -8.52 | -8.40 | -8.76 | -8.74 | -8.61 | 0.17
20000618 0008 | 30.97 | 1713.8542 | 0.432 | -8.43 | -8.42 | -8.84 | -8.80 | -8.62 | 0.23
20000619 0024 | 30.49 | 1714.8368 | 0.625 | -8.85 | -8.47 | -9.47 | -8.77 | -8.89 | 0.42
20000619 0025 | 30.54 | 1714.8715 | 0.632 | -9.37 | -8.68 | -9.22 | -8.76 | -9.01 | 0.34
20000619 0026 | 30.04 | 1714.8742 | 0.633 | -9.23 | -8.96 | -9.17 | -8.83 | -9.05 | 0.19
Table 2: Best Keplerian orbital solution derived for $\epsilon$ UMa Parameter |  | Value | Error
---|---|---|---
P (fixed) | (d) | 5.0887 |
T (periastron) | (JD-2450000) | 1743.932 | 0.161
e |  | 0.503 | 0.063
${V}_{\rm 0}$ | (km ${\rm s}^{-1}$) | -8.920 | 0.030
$\omega$ | ($\degr$) | 260.53 | 12.91
a $\sin$(i) | (${\rm 10}^{-3}$ au) | 0.256 | 0.022
K | (km ${\rm s}^{-1}$) | 0.634 | 0.067
f(m) | (${\rm 10}^{-7}$ ${\rm M}_{\odot}$) | 0.868 | 0.227

The spectra of $\epsilon$ UMa were retrieved from the ELODIE archive (Moultaka
et al., 2004). This archive contains the complete collection of high-
resolution echelle spectra using the ELODIE fiber-fed echelle spectrograph
(Baranne et al., 1996) mounted on the 1.93-m telescope at the Haute-Provence
Observatory (France). The spectra have a resolution
($\lambda$/$\triangle$$\lambda$) of about 42000. The archived signal-to-noise
ratio was between $\sim$200 and $\sim$400 in the spectral region near
$\lambda$ 5550 Å. In addition, eight spectra obtained 15 and 16 June 2000 at
the same telescope were retrieved from the Hypercat Fits Archive
(http://leda.univ-lyon1.fr/11/spectrophotometry.html). In Table 1, each
spectrum is presented by its file name, exposure time, Julian date of the
observations and its corresponding phase. Note that the phases were computed
with the ephemeris from Table 2 (see below).

Special technique of radial velocity measurement of the hydrogen lines was
used and should be explained. Well known, that in slowly rotating late B-type
stars the hydrogen lines have broad wings and sharp cores. For $\epsilon$ UMa
the value of $v$ sin $i$ is equal 35 km ${\rm s}^{-1}$ (Lueftinger et al.,
2003). Unfortunately, the wings of the hydrogen lines can be affected by
different spectral lines. For example, in spectrum of $\epsilon$ UMa the wing
of $H_{\beta}$ line affected by many Cr lines (Žižnovský & Zverko, 1995). On
the other hand, the sharp cores of the hydrogen lines are not affected by
spectral lines and the flux is formed in the upper layers of the atmosphere.
Thus our technique was to use only the sharp cores of the hydrogen lines in
spectra of $\epsilon$ UMa. The spectra were processed using the spectral
reduction software SPE developed by S. Sergeev at Crimean Astrophysical
Observatory (CrAO). The program allows detecting the variations of the centre
gravity of the sharp cores of the hydrogen lines. After processing, all the
spectra were corrected for the motion of the Earth around the Sun. For
example, the demonstration of the positional variability of the core
$H_{\alpha}$ line is presented in Fig. 1. In Table 1, each spectrum is
presented by its radial velocities computed from the hydrogen lines, the mean
radial velocity and errors of the mean radial velocity.

Figure 1: The observed intensity profiles of the cores of $H_{\alpha}$ line
obtained on May 21, 2000 (phase = 0.93) and on June 6, 2000 (phase = 0.07)
marked by the solid and dashed lines, respectively. The velocity scale is
given with respect to the $\lambda$ = 6562.797 Å.

##  3 Orbital parameters

Orbital elements have been determined by a non-linear least-squares fitting of
the mean radial velocities from Table 1 using the program BINARY writing D.H.
Gudehus from Georgia State University
(http://www.chara.gsu.edu/$\sim$gudehus/binary.html). The solution for a
single-lined binary is modelled by up to six parameters:

  1. 1.

P period

  2. 2.

T time of periastron passage

  3. 3.

e eccentricity

  4. 4.

${V}_{\rm 0}$ system radial velocity

  5. 5.

$\omega$ longitude of periastron

  6. 6.

a $\sin$(i) sine of inclination times semi-major axis

The expected radial velocities are

$RV=K[\cos{(\theta+\omega)}+e\cos{\omega}]$ (1)

where $\theta$ is the angular position of the star measured from the centre of
mass at a given instant. The program also calculate the projected semi-
amplitude variation of the radial velocity:

$K=\frac{2\pi a\sin{i}}{P\sqrt{1-e^{2}}},$ (2)

though this is never used as a parameter in the solution and the mass
function:

$f(m)=\frac{M^{3}_{2}{\sin}^{3}i}{(M_{1}+M_{2})^{2}}.$ (3)

The program BINARY gives the estimated standard deviations of the orbital
parameters as well. The orbital solution of Table 2 was obtained by fitting a
Keplerian orbit to the 29 ELODIE radial velocity measurements. Note that most
of observations were obtained between JD=2451568 (January 2000) and JD=2451714
(June 2000) (see Table 1).

Experience shows that the best Keplerian fit to the data with the fixed period
P = 5.0887 d, the eccentricity e = 0.503 and the semi-amplitude K = 0.634 km
${\rm s}^{-1}$. The parameters of the best Keplerian orbital solution for
$\epsilon$ UMa are presented in Table 2. In close binary system with Bp-Ap
stars, there is evidence for a tendency toward synchronization between the
rotational and orbital motions. This effect is thought to be produced by the
tidal forces (Gerbaldi et al., 1985). The radial velocity curve is displayed
in Fig. 2 with the residuals around solution. A linear trend is not observed
in the residuals around the orbital solution that can be explained by the
absent of a second companion in a longer-period orbit. Although, the weighted
r.m.s. around the best Keplerian solution ($\sigma$(O-C)) is equal to 0.131 km
${\rm s}^{-1}$. This value is a bit large compared to the typical radial
velocity measurements from the ELODIE cross-correlation function (Naef et al.,
2003).

Figure 2: $Top:$ Phase diagram of the radial velocity measurements and
Keplerian orbital solution for $\epsilon$ UMa. $Bottom:$ Residuals around the
solution.

##  4 Discussion

The presence of the spots on the stellar surface of $\epsilon$ UMa can change
the observed spectral line profiles and induces a periodic radial velocity
signal similar to the one expected from the presence of a satellite. The
hydrogen lines analysis is one of the best tools to discriminate between
radial velocity variations due to changes in the spectral line shapes and
variations due to the real orbital motion of the star. It is obviously of
interest to compare the phase diagram of the radial velocity of $\epsilon$ UMa
derived above with the phase diagram of the radial velocity computed from
metallic lines. In this way, we selected the spectral region of the very
prominent unblended Cr II line at $\lambda$ 4558 Å.

Figure 3: Positional variations of strong photosphere lines in the spectral
region of the Cr II line at $\lambda$ 4558 Å obtained on June 15, 2000 (phase
= 0.84) and on January 26, 2000 (phase = 0.19) marked by the solid and dashed
lines, respectively. The velocity scale is given with respect to the $\lambda$
= 4558.65 Å.

Figure 3 shows the strong photosphere lines in the spectral region of the Cr
II line at $\lambda$ 4558 Å obtained at phases before and after the epoch of
periastron passage. The radial velocity shift of the Cr II $\lambda\lambda$
4554, 4558 and 4565 ÅÅ lines and the Fe II $\lambda$ 4555 Å line at different
phases is clearly seen. The Ti II line at $\lambda$ 4563 Å shows the line
doubling that appears at phase 0.84 and disappears at phase 0.19. The maximum
splitting of this line is at phase 0.04. The splitting of the lines was first
observed by Struve & Hiltner (1943) in the spectrum of $\epsilon$ UMa.
However, the authors concluded that the doubling is not caused by orbital
motion but may be due to a combination of the physical effects with Doppler
effect in a rotating star.

Wade (1997) estimated the masses for 10 magnetic CP stars using the position
of stars in the log$R_{\sun}$ – log$T_{\rm eff}$ plane. Position of $\epsilon$
UMa in the log$R_{\sun}$ – log$T_{\rm eff}$ plane seems to be among the most
evolved CP stars and gives the value of $M_{\epsilon{\rm UMa}}$ =
3.0$\pm$0.4$M_{\sun}$. According to Lueftinger et al. (2003) the radius of
$\epsilon$ UMa is equal to 4.2$\pm$0.2$R_{\sun}$ corresponding to their choice
of $v\sin i$ = 35 km ${\rm s}^{\rm-1}$ and $i$ = ${\rm 45}^{\circ}$ using the
trigonometric parallax measured by Hipparcos, $\pi$ = 40.30 mas (ESO, 1997)
and an angular diameter of 1.561 mas. The effective temperature of $\epsilon$
UMa was taken from the paper by Sokolov (1998) and is equal to 9340$\pm$530 K.
On the hypothesis that the rotational axis of $\epsilon$ UMa is perpendicular
to the orbital plane we can estimate the mass of the secondary star. For the
value of $M_{\rm 1}$ = 3.0$M_{\sun}$ and $i$ = ${\rm 45}^{\circ}$ Eq.(3) gives
$M_{\rm 2}$ = 0.014$M_{\sun}$. This result gives a value $\sim$14.7$M_{\rm
Jup}$, strongly suggesting that the companion is in the typical brown-dwarf
regime.

If we know the value $M_{\rm 1}$/$M_{\rm 2}$ then it is possible to estimate
the projected semi-amplitude variation of the radial velocity for the
companion according to formula:

$K_{2}=K_{1}\frac{M_{1}}{M_{2}},$ (4)

where $K_{1}$ is the projected semi-amplitude variation of the radial velocity
for $\epsilon$ UMa taken from Table 2. For the value of $K_{\rm 1}$ = 0.634 km
${\rm s}^{\rm-1}$ Eq.(4) gives $K_{\rm 2}$ = 135.9 km ${\rm s}^{\rm-1}$. Thus,
we can estimate the sine of inclination times semi-major axis using Eq.(2).
Computation gives the value of $a_{2}$$\sin$(i) = 0.055 au. These estimates
show that the proposed brown-dwarf is quite close to the surface of $\epsilon$
UMa at periastron. But, such close orbits are not new. For example, the
subgiant star HD 118203 have the planet with eccentric orbit ($e$ = 0.31), the
period of $P$ = 6.1335 d, and is close to its parent star ($a$ = 0.06 au) (Da
Silva et al., 2005).

Another way to interpret the radial velocity variations is the radial
pulsation of $\epsilon$ UMa. Retter et al. (2004) are analysed observations of
$\epsilon$ UMa obtained with the star tracker on the Wide Field Infrared
Explorer satellite. The authors observed that a light curve has about 2 per
cent amplitude of photometric variability. On the other hand, Molnar (1975)
has presented ultraviolet light curves for $\epsilon$ UMa from the OAO-2
satellite which indicate that the photometric variations of this star are due
to variable ultraviolet absorption effects which redistribute flux into the
visible region (see his Fig. 5). Note that no colour index changes (Provin,
1953). Certainly the radial pulsation appears unlikely given that the
rotational period is synchronized with the orbital period of $\epsilon$ UMa.

##  5 Conclusions

The archival ELODIE high-resolution echelle spectra of $\epsilon$ UMa permit
us to analyse the radial velocity variations of the sharp cores of the
hydrogen lines. This allowed determining the orbital elements of binary system
for the CP star $\epsilon$ UMa. The best Keplerian fit to the data shown that
the rotational period is synchronized with the orbital period. We are
estimated the mass of the secondary star which is equal $\sim$14.7$M_{\rm
Jup}$. This result indicate that the companion is the brown-dwarf with the
projected semi-amplitude variation of the radial velocity $K_{\rm 2}$ = 135.9
km ${\rm s}^{\rm-1}$ and the sine of inclination times semi-major axis
$a_{2}$$\sin$(i) = 0.055 au.

## Acknowledgements

The author would like to thank the referee Dr. J.B. Rice of this Letter for
his extremely helpful comments.

## References

  * Abt & Snowden (1973) Abt H.A., Snowden M.S., 1973, ApJSS, 25, 137
  * Baranne et al. (1996) Baranne A., Queloz D., Mayor M., Adrianzyk G., Knispel G., Kohler D., Lacroix D., Meunier J.-P., et al, 1996, A&AS, 119, 373
  * Da Silva et al. (2005) Da Silva R., Udry S., Bouchy F., Mayor M., Moutou C., Pont F., Queloz D., Santos N.C., et al, 2005, (astro-ph/0510048)
  * Gerbaldi et al. (1985) Gerbaldi M., Floquet M., Hauck B., 1985, A&A, 146, 341
  * Guthnick (1934) Guthnick P., 1934, Sitz. Preuss. Akad. Wiss. Berlin, 30, 506
  * Hatzes (1991) Hatzes A.P., 1991, MNRAS, 253, 89
  * Holmgren & Rice (2000) Holmgren D.E., Rice J.B., 2000, A&A, 364, 660
  * Lueftinger et al. (2003) Lueftinger T., Kuschnig R., Piskunov N.E., Weiss W.W., 2003, A&A, 406, 1033
  * Molnar (1975) Molnar M.R., 1975, AJ, 80, 137
  * Morgan et al. (1978) Morgan B.L., Beddoes D.R., Scaddan R.J., Dainty J.C., 1978, MNRAS, 183, 701
  * Moultaka et al. (2004) Moultaka J., Ilovaisky S.A., Prugniel P., Soubiran C., 2004, PASP, 116, 693
  * Naef et al. (2003) Naef D., Mayor M., Beuzit J.L., Perrier C., Queloz D., Sivan J.P., Udry S., 2003, (astro-ph/0310261)
  * Provin (1953) Provin S.S., 1953, ApJ, 118, 489
  * Retter et al. (2004) Retter A., Bedding T,R., Buzasi D.L., Kjeldsen H., Kiss L.L., 2004, ApJ, 601, L95
  * Rice & Wehlau (1990) Rice J.B., Wehlau W.H., 1990, A&A, 233, 503
  * Rice et al. (1997) Rice J.B., Wehlau W.H., Holmgren D.E., 1997, A&A, 326, 988
  * Sokolov (1998) Sokolov N.A., 1998, A&AS, 130, 215
  * Struve & Hiltner (1943) Struve O., Hiltner W.A, 1943, ApJ, 98, 225
  * Stibbs (1950) Stibbs D.W.N., 1950, MNRAS, 110,395
  * Wade (1997) Wade G.A., 1997, A&A, 325,1063
  * Wehlau et al. (1982) Wehlau W., Rice J., Piskunov N., Khokhlova V., 1982, Pis’ma Astron. Zh., 8, 30
  * Woszczyk & Jasiński (1980) Woszczyk A., Jasiński M., 1980, Acta Astron., 30, 331
  * Žižnovský & Zverko (1995) Žižňovský J., Zverko J., 1995, Contrib. Astron. Obs. Skalnaté Pleso, 25, 39