# Metric tensors for the interpolation error and its gradient in $L^{p}$ norm

Hehu Xie111LSEC, ICMSEC, Academy of Mathematics and Systems Science, CAS,
Beijing 100080, China email<EMAIL_ADDRESS>Xiaobo Yin222Department of
Mathematics, Central China Normal University, Wuhan 430079, China email:
<EMAIL_ADDRESS>
> Abstract. A uniform strategy to derive metric tensors in two spatial
> dimension for interpolation errors and their gradients in $L^{p}$ norm is
> presented. It generates anisotropic adaptive meshes as quasi-uniform ones in
> corresponding metric space, with the metric tensor being computed based on a
> posteriori error estimates in different norms. Numerical results show that
> the corresponding convergence rates are always optimal.
>
> Keywords. metric tensor, interpolation; gradient; anisotropic.
>
> AMS subject classification. 65N30, 65N50

##  1 Introduction

Generation of adaptive meshes is now the standard option in most software
packages. Traditionally, isotropic mesh adaptation has received much
attention, where regular mesh elements are only adjusted in size based on an
error estimate. However, for problems with anisotropic solutions (with, say,
sharp boundary or internal layers), the shape of elements can be further
optimized and an equidistribution of a scalar error density is not sufficient
to ensure that a mesh is optimally efficient [14]. Indeed anisotropic meshes
have been used successfully in many areas, for example in singular
perturbation and flow problems [4, 5, 6, 21, 22, 35, 43] and in adaptive
procedures [2, 7, 8, 10, 11, 23, 35, 36]. For anisotropic mesh adaptation, the
common practice is to generate the needed anisotropic mesh as a quasi-uniform
one in the metric space determined by a tensor (or a matrix-valued function),
always called monitor function or metric tensor. Both the monitor function
(denoted by the letter $M$) and metric tensor (denoted by the calligraphy
letter $\mathcal{M}$) play the same role in mesh generation, i.e., they are
used to specify the size, shape, and orientation of mesh elements throughout
the physical domain. The only difference lies in the way they specify the size
of elements. Indeed, the former specifies the element size through the
equidistribution condition, while the latter determines the element size
through the unitary volume requirement. Readers could regard the metric tensor
as normalization for the monitor function. Examples of anisotropic meshing
strategies include blue refinement [29, 30], directional refinement [36],
Delaunay-type triangulation method [7, 8, 11, 35], advancing front method
[19], bubble packing method [41], local refinement and modification [21, 37],
variational methods [9, 17, 24, 27, 28, 31], and so on. Readers are referred
to [18] and [34] for an overview.

Among these meshing strategies, the definition of the metric tensor (or
monitor function) based on the Hessian of the solution seems widespread in the
meshing community[1, 11, 12, 13, 14, 15, 20, 21, 23, 25, 26, 24, 37, 39].
Especially, Huang and Russell [26] propose the monitor function

$\displaystyle
M=\det{\Big{(}}I+\frac{1}{\alpha}|H(u)|{\Big{)}}^{-\frac{1}{d+p(2-m)}}{\Big{\|}}I+\frac{1}{\alpha}|H(u)|{\Big{\|}}^{\frac{mp}{d+p(2-m)}}{\Big{[}}I+\frac{1}{\alpha}|H(u)|{\Big{]}},$
(1.1)

for the interpolation error in $W^{m,p}$ norm ($m=0,1$, $p\in[1,+\infty)$),
where $d$ stands for the spatial dimension. Set
$\mathcal{H}=I+\frac{1}{\alpha}|H(u)|$, when $d=2$,

$\displaystyle
M_{m,p}=\det(\mathcal{H})^{-\frac{1}{2+p(2-m)}}\|\mathcal{H}\|^{\frac{mp}{2+p(2-m)}}\mathcal{H}.$
(1.2)

Separately, it becomes

$\displaystyle M_{0,p}=\det(\mathcal{H})^{-\frac{1}{2(p+1)}}\mathcal{H},$
(1.3)

for the interpolation error in $L^{p}$ norm and

$\displaystyle
M_{1,p}=\det(\mathcal{H})^{-\frac{1}{p+2}}\|\mathcal{H}\|^{\frac{p}{p+2}}\mathcal{H}.$
(1.4)

for the gradient of interpolation error in $L^{p}$ norm.

The objective of this paper is to give a unified strategy deriving metric
tensors in two spatial dimension for interpolation error and its gradients in
$L^{p}$ norm. The development begin with the error estimates [32] for $L^{2}$
norm and our recent work [42] for $H^{1}$ norm on linear interpolation for
quadratic functions on triangles. These estimates are anisotropic in the sense
that they allow a full control of the shape of elements when used within a
mesh generation strategy. Using the relationship between different norms, a
posterior error estimates for other norms ($W^{m,p},m=0,1$, $p\neq 2$) can be
gained. We will apply these error estimates to formulate corresponding metric
tensors in a unified way. The procedure is based on two considerations: on the
one hand the anisotropic mesh is generated as a quasi-uniform mesh in the
metric tensor. On the other hand, the anisotropic mesh is required to minimize
the error for a given number of triangles. To compare with those existing
methods, we list our main results using monitor function style, that is

$\displaystyle M_{0,p}^{n}({\bf
x})=\det(\mathcal{H})^{-\frac{1}{2(p+1)}}\mathcal{H},$ (1.5)

for interpolation errors in $L^{p}$ norm and

$\displaystyle M_{1,p}^{n}({\bf
x})=\det(\mathcal{H})^{-\frac{1}{p+2}}\mbox{tr}(\mathcal{H})^{\frac{p}{p+2}}\mathcal{H},$
(1.6)

for gradient of interpolation errors in $L^{p}$ norm. To sum up, the metric
tensor can be expressed by

$\displaystyle M_{m,p}^{n}({\bf
x})=\det(\mathcal{H})^{-\frac{1}{2+p(2-m)}}\mbox{tr}(\mathcal{H})^{\frac{mp}{2+p(2-m)}}\mathcal{H},$
(1.7)

for the $W^{m,p}$ norm ($m=0,1$, $p\in(0,+\infty]$) of the interpolation
error.

The paper is organized as follows. In Section 2, we describe the anisotropic
error estimates on linear interpolation for quadratic functions on triangles
obtained in our recent work [42]. The formulation of the monitor function and
metric tensor is developed in Section 3. Numerical results are presented in
Section 4 to illustrating our analysis. Finally, conclusions are drawn in
Section 5.

##  2 Estimates for interpolation error and its gradient

As we know, the interpolation error depends on the solution, the size and
shape of the elements in the mesh. Understanding this relation is crucial for
the generating efficient meshes for the finite element method. In the mesh
generation community, this relation is studied more closely for the model
problem of interpolating quadratic functions. This treatment yields a reliable
and efficient estimator of the interpolation error for general functions
provided a saturation assumption is valid [3, 16]. For instance, Nadler [32]
derived an exact expression for the $L^{2}$-norm of the linear interpolation
error in terms of the three sides ${\bf\ell}_{1}$, ${\bf\ell}_{2}$, and
${\bf\ell}_{3}$ of the triangle $K$,

$\displaystyle\|u-u_{I}\|^{2}_{L^{2}(K)}=\frac{|K|}{180}{\Big{[}}{\Big{(}}d_{1}+d_{2}+d_{3}{\Big{)}}^{2}+d_{1}d_{2}+d_{2}d_{3}+d_{1}d_{3}{\Big{]}},$
(2.1)

where $|K|$ is the area of the triangle, $d_{i}={\bf\ell}_{i}\cdot
H{\bf\ell}_{i}$ with $H$ being the Hessian of $u$. Assuming
$u=\lambda_{1}x^{2}+\lambda_{2}y^{2}$, D’Azevedo and Simpson [13] derived the
exact formula for the maximum norm of the interpolation error

$\displaystyle\|(u-u_{I})\|^{2}_{L^{\infty}(K)}=\frac{D_{12}D_{23}D_{31}}{16\lambda_{1}\lambda_{2}|K|^{2}},$
(2.2)

where
$D_{ij}={\bf\ell}_{i}\cdot\mbox{diag}(\lambda_{1},\lambda_{2}){\bf\ell}_{j}$.
Based on the geometric interpretation of this formula, they proved that for a
fixed area the optimal triangle, which produces the smallest maximum
interpolation error, is the one obtained by compressing an equilateral
triangle by factors $\sqrt{\lambda_{1}}$ and $\sqrt{\lambda_{2}}$ along the
two eigenvectors of the Hessian of $u$. Furthermore, the optimal incidence for
a given set of interpolation points is the Delaunay triangulation based on the
stretching map (by factors $\sqrt{\lambda_{1}}$ and $\sqrt{\lambda_{2}}$ along
the two eigenvector directions) of the grid points. Rippa [38] showed that the
mesh obtained in this way is also optimal for the $L^{p}$-norm of the error
for any $1\leq p\leq\infty$.

The element-wise error estimates in the following theorem are developed in
[42] using the theory of interpolation and proper numerical quadrature
formula.

######  Theorem 2.1.

Let $u$ be a quadratic function and $u_{I}$ is the Lagrangian linear finite
element interpolation of $u$. The following relationship holds:

$\displaystyle\|\nabla(u-u_{I})\|^{2}_{L^{2}(K)}=\frac{1}{48|K|}\sum_{i=1}^{3}({\bf\ell}_{i+1}\cdot
H{\bf\ell}_{i+2})^{2}|{\bf\ell}_{i}|^{2},$ (2.3)

where we prescribe $i+3=i,i-3=i$.

To get the a posteriori error estimate of the interpolation error in $L^{p}$
and $W^{1,p}$ norms for $p\neq 2$, we need some lemmas below.

######  Lemma 2.1.

For any $d$ positive numbers $a_{1},\cdots,a_{d}$, the inequalities

$\displaystyle{\Big{(}}\sum\limits_{j=1}^{d}a_{j}^{2}{\Big{)}}^{\frac{1}{2}}\leq{\Big{(}}\sum\limits_{j=1}^{d}a_{j}^{p}{\Big{)}}^{\frac{1}{p}}\leq
d^{\frac{1}{p}-\frac{1}{2}}{\Big{(}}\sum\limits_{j=1}^{d}a_{j}^{2}{\Big{)}}^{\frac{1}{2}},$
(2.4)

and

$\displaystyle
d^{\frac{1}{p}-\frac{1}{2}}{\Big{(}}\sum\limits_{j=1}^{d}a_{j}^{2}{\Big{)}}^{\frac{1}{2}}\leq{\Big{(}}\sum\limits_{j=1}^{d}a_{j}^{p}{\Big{)}}^{\frac{1}{p}}\leq{\Big{(}}\sum\limits_{j=1}^{d}a_{j}^{2}{\Big{)}}^{\frac{1}{2}}$
(2.5)

hold for numbers $0<p<2$ and $p>2$, respectively.

###### Proof.

We just give the proof for the case $0<p<2$, it is similar for the case $p>2$.

For any number $0<p<2$,

$\displaystyle{\Big{(}}\sum\limits_{j=1}^{d}a_{j}^{2}{\Big{)}}^{\frac{1}{2}}\leq{\Big{(}}\sum\limits_{j=1}^{d}a_{j}^{p}{\Big{)}}^{\frac{1}{p}}$

holds due to the Jensen’s inequality. From the generalized arithmetic-mean
geometric-mean inequality, for any positive numbers $a_{1},\cdots,a_{d}$,

$\displaystyle{\Big{(}}\sum\limits_{j=1}^{d}\frac{1}{d}a_{j}^{p}{\Big{)}}^{\frac{1}{p}}\leq{\Big{(}}\sum\limits_{j=1}^{d}\frac{1}{d}a_{j}^{2}{\Big{)}}^{\frac{1}{2}}.$

Then

$\displaystyle{\Big{(}}\sum\limits_{j=1}^{d}a_{j}^{p}{\Big{)}}^{\frac{1}{p}}\leq
d^{\frac{1}{p}-\frac{1}{2}}{\Big{(}}\sum\limits_{j=1}^{d}a_{j}^{2}{\Big{)}}^{\frac{1}{2}}.$

∎

To sum up, for any $d$ positive numbers $a_{1},\cdots,a_{d}$, the inequalities

$\displaystyle\underline{C}_{p}{\Big{(}}\sum\limits_{j=1}^{d}a_{j}^{2}{\Big{)}}^{\frac{1}{2}}\leq{\Big{(}}\sum\limits_{j=1}^{d}a_{j}^{p}{\Big{)}}^{\frac{1}{p}}\leq\overline{C}_{p}{\Big{(}}\sum\limits_{j=1}^{d}a_{j}^{2}{\Big{)}}^{\frac{1}{2}}$
(2.6)

holds for any numbers $p>0$, where $\underline{C}_{p}=1$ for $0<p<2$ and
$d^{\frac{1}{p}-\frac{1}{2}}$ for $p>2$,
$\overline{C}_{p}=d^{\frac{1}{p}-\frac{1}{2}}$ for $0<p<2$ and $1$ for $p>2$.

######  Lemma 2.2.

[3] For any $p\in(0,+\infty]$ and any non-negative $v\in P_{2}(K)$ it holds

$\displaystyle
C_{1/p}^{-\frac{1}{p}}|K|^{\frac{1}{p}-1}\|v\|_{L^{1}(K)}\leq\|v\|_{L^{p}(K)}\leq
C_{p}|K|^{\frac{1}{p}-1}\|v\|_{L^{1}(K)}$ (2.7)

with

$\left\\{\begin{array}[]{lll}C_{p}=1&\mbox{if}\,\,0<p\leq 1,&\\\
C_{p}=(d+1)(d+2)(d!)^{\frac{1}{p}}{\Big{(}}\prod\limits_{j=1}^{d}(p+j){\Big{)}}^{-\frac{1}{p}}&\mbox{if}\,\,1<p<+\infty,&\\\
C_{\infty}=\lim\limits_{p\rightarrow+\infty}C_{p}=(d+1)(d+2),\\\
C_{1/\infty}=\lim\limits_{p\rightarrow+\infty}C_{1/p}=1.\end{array}\right.$

###  2.1 Estimates for interpolation errors in $L^{p}$ norm

We consider the error of linear interpolation $e=u-u_{I}$ for a quadratic
function $u$ on $K$. Since the function $e$ is quadratic on $K$, we can apply
Lemma 2.2 to obtain

$\displaystyle
C_{1/p}^{-1/p}|K|^{\frac{1}{p}-1}\|e\|_{L^{1}(K)}\leq\|e\|_{L^{p}(K)}\leq
C_{p}|K|^{\frac{1}{p}-1}\|e\|_{L^{1}(K)}.$ (2.8)

Set $p=2$,

$\displaystyle|K|^{-\frac{1}{2}}\|e\|_{L^{1}(K)}\leq\|e\|_{L^{2}(K)}\leq
C_{2}|K|^{-\frac{1}{2}}\|e\|_{L^{1}(K)},$

or

$\displaystyle
C_{2}^{-1}|K|^{\frac{1}{2}}\|e\|_{L^{2}(K)}\leq\|e\|_{L^{1}(K)}\leq|K|^{\frac{1}{2}}\|e\|_{L^{2}(K)}.$
(2.9)

Combine (2.8) and (2.9), we get

$\displaystyle
C_{1/p}^{-1/p}C_{2}^{-1}|K|^{\frac{1}{p}-\frac{1}{2}}\|e\|_{L^{2}(K)}\leq\|e\|_{L^{p}(K)}\leq
C_{p}|K|^{\frac{1}{p}-\frac{1}{2}}\|e\|_{L^{2}(K)}.$ (2.10)

In this article, $A\sim B$ stands for that there exist two constants
$\underline{C}$ and $\overline{C}$ such that

$\displaystyle\underline{C}A\leq B\leq\overline{C}A,$

where the two constants $\underline{C}$ and $\overline{C}$ may depend on the
prescribed error, the index $p$, the dimension $d$, and the numbers of
elements $N$, however are independent of function at hand. So (2.10) can be
rewritten as

$\displaystyle\|e\|_{L^{p}(K)}\sim|K|^{\frac{1}{p}-\frac{1}{2}}\|e\|_{L^{2}(K)}.$

Together with the expression (2.1) for the $L^{2}$ norm of the linear
interpolation error derived by Nadler[32], we have the a posteriori error
estimate in $L^{p}$ norms as follows:

$\displaystyle\|e\|^{2}_{L^{p}(K)}\sim|K|^{\frac{2}{p}-1}\|e\|^{2}_{L^{2}(K)}$
$\displaystyle=$
$\displaystyle\frac{|K|^{\frac{2}{p}}}{180}{\Big{[}}{\Big{(}}\sum_{i=1}^{3}d_{i}{\Big{)}}^{2}+d_{1}d_{2}+d_{2}d_{3}+d_{1}d_{3}{\Big{]}}.$
(2.11)

###  2.2 Estimates for gradient of interpolation errors in $L^{p}$ norm

Now we consider the gradient of linear interpolation error $\nabla
e=\nabla(u-u_{I})$ for a quadratic function $u$. Since the function

$\displaystyle v_{j}(x)={\Big{(}}\frac{\partial e}{\partial
x_{j}}{\Big{)}}^{2}$

is quadratic on $K$, we can apply Lemma 2.2 to obtain

$\displaystyle\|v_{j}\|_{L^{p/2}(K)}^{1/2}\geq
C_{2/p}^{-1/p}|K|^{\frac{1}{p}-\frac{1}{2}}\|v_{j}\|_{L^{1}(K)}^{1/2}=C_{2/p}^{-1/p}|K|^{\frac{1}{p}-\frac{1}{2}}{\Big{\|}}\frac{\partial
e}{\partial x_{j}}{\Big{\|}}_{L^{2}(K)},$ (2.12)

and

$\displaystyle\|v_{j}\|_{L^{p/2}(K)}^{1/2}\leq
C_{p/2}^{1/2}|K|^{\frac{1}{p}-\frac{1}{2}}\|v_{j}\|_{L^{1}(K)}^{1/2}=C_{p/2}^{1/2}|K|^{\frac{1}{p}-\frac{1}{2}}{\Big{\|}}\frac{\partial
e}{\partial x_{j}}{\Big{\|}}_{L^{2}(K)}.$ (2.13)

Since

$\displaystyle\|\nabla
e\|_{L^{p}(K)}^{p}=\sum\limits_{j=1}^{d}\int_{K}{\Big{|}}\frac{\partial
e}{\partial x_{j}}{\Big{|}}^{p}d{\bf
x}=\sum\limits_{j=1}^{d}\|v_{j}\|_{L^{p/2}(K)}^{p/2},$

then together with (2.12) and (2.13), we have

$\displaystyle
C_{2/p}^{-1/p}|K|^{\frac{1}{p}-\frac{1}{2}}{\Big{(}}\sum\limits_{j=1}^{d}{\Big{\|}}\frac{\partial
e}{\partial x_{j}}{\Big{\|}}_{L^{2}(K)}^{p}{\Big{)}}^{\frac{1}{p}}\leq\|\nabla
e\|_{L^{p}(K)}\leq
C_{p/2}^{1/2}|K|^{\frac{1}{p}-\frac{1}{2}}{\Big{(}}\sum\limits_{j=1}^{d}{\Big{\|}}\frac{\partial
e}{\partial x_{j}}{\Big{\|}}_{L^{2}(K)}^{p}{\Big{)}}^{\frac{1}{p}}.$

From (2.6), the inequality

$\displaystyle\underline{C}_{p}C_{2/p}^{-1/p}|K|^{\frac{1}{p}-\frac{1}{2}}\|\nabla
e\|_{L^{2}(K)}\leq\|\nabla
e\|_{L^{p}(K)}\leq\overline{C}_{p}C_{p/2}^{1/2}|K|^{\frac{1}{p}-\frac{1}{2}}\|\nabla
e\|_{L^{2}(K)},$

holds, or simply

$\displaystyle\|\nabla e\|_{L^{p}(K)}\sim|K|^{\frac{1}{p}-\frac{1}{2}}\|\nabla
e\|_{L^{2}(K)}.$

Together with the a posteriori error estimate (2.3) of the interpolation error
in $H^{1}$($=W^{1,2}$) norm, we have the a posteriori error estimate in
$W^{1,p}$ norms as follows:

$\displaystyle\|\nabla e\|^{2}_{L^{p}(K)}$ $\displaystyle\sim$
$\displaystyle|K|^{\frac{2}{p}-1}\|\nabla e\|_{L^{2}(K)}^{2}$ (2.14)
$\displaystyle=$
$\displaystyle|K|^{\frac{2}{p}-1}\frac{1}{48|K|}\sum_{i=1}^{3}({\bf\ell}_{i+1}\cdot
H_{K}{\bf\ell}_{i+2})^{2}|{\bf\ell}_{i}|^{2}$ $\displaystyle=$
$\displaystyle\frac{|K|^{\frac{2}{p}-2}}{48}\sum_{i=1}^{3}({\bf\ell}_{i+1}\cdot
H_{K}{\bf\ell}_{i+2})^{2}|{\bf\ell}_{i}|^{2}.$

##  3 Metric tensors for anisotropic mesh adaptation

We now use the results of Section 2 to develop metric tensors for
interpolation errors and their gradients in $L^{p}$ norm in a unified way. As
a common practice in anisotropic mesh generation, the metric tensor,
$\mathcal{M}({\bf x})$, is used in a meshing strategy in such a way that an
anisotropic mesh is generated as a quasi-uniform mesh in the metric space
determined by $\mathcal{M}({\bf x})$. Mathematically, this can be interpreted
as the shape, size and equidistribution requirements as follows.

The shape requirement. The elements of the new mesh, $\mathcal{T}_{h}$, are
(or are close to being) equilateral in the metric.

The size requirement. The elements of the new mesh $\mathcal{T}_{h}$ have a
unitary volume in the metric, i.e.,

$\displaystyle\int_{K}\sqrt{\det(\mathcal{M}({\bf x}))}d{\bf x}=1,\quad\forall
K\in\mathcal{T}_{h}.$ (3.1)

The equidistribution requirement.  The anisotropic mesh is required to
minimize the error for a given number of mesh points (or equidistribute the
error on every element).

Notice that to derive the monitor function, we just need the shape and
equidistribution requirements.

###  3.1 Metric tensors for gradients of interpolation errors in $L^{p}$ norm

$\mathcal{F}_{K}$$\bf\ell_{1}$$\bf\ell_{2}$$\bf\ell_{3}$$\theta_{1}$$\theta_{2}$$\theta_{3}$$\hat{\theta}_{1}$$\hat{\theta}_{2}$$\hat{\theta}_{3}$$\hat{\ell}_{1}$$\hat{\ell}_{2}$$\hat{\ell}_{3}$

Figure 1: Affine map ${\bf\hat{x}}=\mathcal{F}_{K}{\bf x}$ from $K$ to the
reference triangle $\hat{K}$.

We derive the monitor function $M({\bf x})$ first. Assume $H(u)$ be a
symmetric positive definite matrix on every point ${\bf x}$, this assumption
will be dropped later. Set $M({\bf x})=C({\bf x})H(u)$. Consider the $L^{2}$
projection of $H(u)$ on $K$, denoted by $H_{K}$, then so does $M_{K}$. Since
$H_{K}$ is a symmetric positive definite matrix, we consider the singular
value decomposition $H_{K}=R^{T}\Lambda R$, where
$\Lambda=\mbox{diag}(\lambda_{1},\lambda_{2})$ is the diagonal matrix of the
corresponding eigenvalues ($\lambda_{1},\lambda_{2}>0$) and $R$ is the
orthogonal matrix having as rows the eigenvectors of $H_{K}$. Denote by
$F_{K}$ and ${\bf t}_{K}$ the matrix and the vector defining the invertible
affine map $\hat{\bf x}=\mathcal{F}_{K}({\bf x})=F_{K}{\bf x}+{\bf t}_{K}$
from the generic element $K$ to the reference triangle $\hat{K}$ (see Figure
1).

Obviously, $M_{K}=C_{K}H_{K}$. Let $M_{K}=F_{K}^{T}F_{K}$, then
$F_{K}=C_{K}^{\frac{1}{2}}\Lambda^{\frac{1}{2}}R$. Mathematically, the shape
requirement can be expressed as

$\displaystyle|\hat{\ell}_{i}|=L\,\,\mbox{and}\,\,\cos\hat{\theta}_{i}=\frac{\hat{\ell}_{i+1}\cdot\hat{\ell}_{i+2}}{L^{2}}=\frac{1}{2},\,i=1,2,3,$
(3.2)

where $L$ is a constant for every element $K$. Enforcing the shape
requirement, we get

$\displaystyle\|\nabla e\|^{2}_{L^{p}(K)}$ $\displaystyle\sim$
$\displaystyle\frac{|K|^{\frac{2}{p}-2}}{48}\sum_{i=1}^{3}({\bf\ell}_{i+1}\cdot
H_{K}{\bf\ell}_{i+2})^{2}|{\bf\ell}_{i}|^{2}$ $\displaystyle=$
$\displaystyle\frac{|K|^{\frac{2}{p}-2}}{48C_{K}^{2}}\sum_{i=1}^{3}({\bf\ell}_{i+1}\cdot
M_{K}{\bf\ell}_{i+2})^{2}|{\bf\ell}_{i}|^{2}$ $\displaystyle=$
$\displaystyle\frac{L^{4}|K|^{\frac{2}{p}-2}}{48C_{K}^{2}}\sum_{i=1}^{3}(\cos\hat{\theta}_{i})^{2}|{\bf\ell}_{i}|^{2}=\frac{L^{4}|K|^{\frac{2}{p}-2}}{192C_{K}^{2}}\sum_{i=1}^{3}|{\bf\ell}_{i}|^{2}.$

Notice that,

$\displaystyle|K|=\frac{|\hat{K}|}{C_{K}\sqrt{\det(H_{K})}},$

we have

$\displaystyle\|\nabla e\|^{2}_{L^{p}(K)}$ $\displaystyle\sim$
$\displaystyle\frac{L^{4}|\hat{K}|^{\frac{2}{p}-2}C_{K}^{2-\frac{2}{p}}\det(H_{K})^{1-\frac{1}{p}}}{192C_{K}^{2}}\sum_{i=1}^{3}{\Big{|}}C_{K}^{-\frac{1}{2}}R^{-1}\Lambda^{-\frac{1}{2}}\hat{\ell}_{i}{\Big{|}}^{2}$
$\displaystyle=$
$\displaystyle\frac{L^{4}|\hat{K}|^{\frac{2}{p}-2}\det(H_{K})^{1-\frac{1}{p}}}{192C_{K}^{1+\frac{2}{p}}}\sum_{i=1}^{3}{\Big{|}}\Lambda^{-\frac{1}{2}}\hat{\ell}_{i}{\Big{|}}^{2}$
$\displaystyle=$
$\displaystyle\frac{L^{4}|\hat{K}|^{\frac{2}{p}-2}\det(H_{K})^{1-\frac{1}{p}}}{192C_{K}^{1+\frac{2}{p}}}\frac{\mbox{tr}(H_{K})}{\det(H_{K})}$
$\displaystyle=$
$\displaystyle\frac{L^{4}|\hat{K}|^{\frac{2}{p}-2}\det(H_{K})^{-\frac{1}{p}}\mbox{tr}(H_{K})}{192C_{K}^{1+\frac{2}{p}}}$
$\displaystyle\sim$
$\displaystyle\frac{\det(H_{K})^{-\frac{1}{p}}\mbox{tr}(H_{K})}{C_{K}^{1+\frac{2}{p}}},$

then

$\displaystyle\|\nabla e\|^{p}_{L^{p}(K)}=(\|\nabla
e\|^{2}_{L^{p}(K)})^{p/2}\sim\frac{\det(H_{K})^{-\frac{1}{2}}\mbox{tr}(H_{K})^{\frac{p}{2}}}{C_{K}^{1+\frac{p}{2}}}.$

To satisfy the equidistribution requirement, let

$\displaystyle\|\nabla
e\|^{p}_{L^{p}(K)}={\Big{(}}\sum\limits_{K\in\mathcal{T}_{h}}e_{K}^{p}{\Big{)}}/N=\epsilon^{p}/N,$

where $N$ is the number of elements of $\mathcal{T}_{h}$. Then

$\displaystyle
C_{K}\sim\det(H_{K})^{-\frac{1}{p+2}}\mbox{tr}(H_{K})^{\frac{p}{p+2}}.$

So $M({\bf x})$ could be the form

$\displaystyle M({\bf
x})=\det(H)^{-\frac{1}{p+2}}\mbox{tr}(H)^{\frac{p}{p+2}}H(u),$

since $M({\bf x})$ can be modified by multiplying a constant. Since it
corresponds the gradient of interpolation errors in $L^{p}$ norm, we denote it
by $M_{1,p}^{n}({\bf x})$.

To establish the metric tensor $\mathcal{M}_{1,p}^{n}({\bf x})$, set
$\mathcal{M}_{1,p}^{n}({\bf x})=\theta_{1,p}M_{1,p}^{n}({\bf x})$, at this
time, the size requirement (3.1) should be used, which leads to

$\displaystyle\theta_{1,p}\int_{K}\rho_{1,p}({\bf x})d{\bf x}=1,$

where

$\displaystyle\rho_{1,p}({\bf x})=\sqrt{\det(M_{1,p}^{n}({\bf x}))}.$

Summing the above equation over all the elements of $\mathcal{T}_{h}$, one
gets

$\displaystyle\theta_{1,p}\sigma_{1,p}=N,$

where

$\displaystyle\sigma_{1,p}=\int_{\Omega}\rho_{1,p}({\bf x})d{\bf x}.$

Thus, we get

$\displaystyle\theta_{1,p}=\frac{N}{\sigma_{1,p}},$

and as a consequence,

$\displaystyle\mathcal{M}_{1,p}^{n}({\bf
x})=\frac{N}{\sigma_{1,p}}\det(H)^{-\frac{1}{p+2}}\mbox{tr}(H)^{\frac{p}{p+2}}H(u).$

###  3.2 Metric tensor for the interpolation errors in $L^{p}$ norm

Using the error estimates (2.11) for interpolation errors in $L^{p}$ norm and
the shape requirement (3.2), we have

$\displaystyle\|e\|^{2}_{L^{p}(K)}$ $\displaystyle\sim$
$\displaystyle\frac{|K|^{\frac{2}{p}}}{180}{\Big{[}}{\Big{(}}\sum_{i=1}^{3}d_{i}{\Big{)}}^{2}+d_{1}d_{2}+d_{2}d_{3}+d_{1}d_{3}{\Big{]}}$
$\displaystyle=$
$\displaystyle\frac{|K|^{\frac{2}{p}}}{180C_{K}^{2}}{\Big{[}}{\Big{(}}\sum_{i=1}^{3}|\hat{\ell}_{i}|^{2}{\Big{)}}^{2}+\sum_{i=1}^{3}{\Big{(}}|\hat{\ell}_{i+1}||\hat{\ell}_{i+2}|{\Big{)}}^{2}{\Big{]}}$
$\displaystyle=$
$\displaystyle\frac{L^{4}|K|^{\frac{2}{p}}}{15C_{K}^{2}}=\frac{L^{4}|\hat{K}|^{\frac{2}{p}}}{15C_{K}^{2+\frac{2}{p}}\det(H)^{\frac{1}{p}}}\sim\frac{1}{C_{K}^{2+\frac{2}{p}}\det(H)^{\frac{1}{p}}}.$

Then,

$\displaystyle\|e\|^{p}_{L^{p}(K)}=(\|e\|^{2}_{L^{p}(K)})^{p/2}\sim\frac{1}{C_{K}^{p+1}\det(H)^{\frac{1}{2}}}.$

To satisfy the equidistribution requirement, let

$\displaystyle\|e\|^{p}_{L^{p}(K)}={\Big{(}}\sum\limits_{K\in\mathcal{T}_{h}}e_{K}^{p}{\Big{)}}/N=\epsilon^{p}/N.$

Using similar argument in last subsection, we easily get monitor functions

$\displaystyle M_{0,p}^{n}({\bf x})=\det(H)^{-\frac{1}{2(p+1)}}H(u),$

and metric tensors

$\displaystyle\mathcal{M}_{0,p}^{n}({\bf
x})=\frac{N}{\sigma_{0,p}}\det(H)^{-\frac{1}{2(p+1)}}H(u),$

for the interpolation errors in $L^{p}$ norm.

###  3.3 Practice use of metric tensor

So far we assume that $H(u)$ is a symmetric positive definite matrix at every
point. However this assumption doesn’t hold in many cases. In order to obtain
a symmetric positive definite matrix, the following procedure are often
implemented. First, the Hessian $H$ is modified into
$|H|=R^{T}\,\,\mbox{diag}(|\lambda_{1}|,|\lambda_{2}|)R$ by taking the
absolute value of its eigenvalues ([22]). Since $|H|$ is only semi-positive
definite, $\mathcal{M}_{m,p}^{n}$ cannot be directly applied to generate the
anisotropic meshes. To avoid this difficulty, we regularize the expression
with the flooring parameter $\alpha_{m,p}>0$ (see, e.g., [24]). Replacing
$|H|$ with

$\displaystyle\mathcal{H}=\alpha_{m,p}I+|H|,$

we get the modified metric tensors, also denoted by $\mathcal{M}_{m,p}^{n}$,
that is

$\displaystyle\mathcal{M}_{m,p}^{n}({\bf
x})=\frac{N}{\sigma_{m,p}}\det(\mathcal{H})^{-\frac{1}{2+p(2-m)}}\mbox{tr}(\mathcal{H})^{\frac{mp}{2+p(2-m)}}\mathcal{H},$
(3.3)

which are suitable for practical mesh generation.

###  3.4 Comparison with existing methods using monitor function style

When $m=0$, the new monitor function $M_{0,p}^{n}$ (1.5) is in fact the same
with (1.3) in [26, 25]. Chen, Sun and Xu [12] proved that under suitable
conditions, the error estimate

$\displaystyle\|u-u_{I}\|_{L^{p}(\Omega)}\leq CN^{-\frac{2}{d}}\|\sqrt[d]{\det
H}\|_{L^{\frac{pd}{2p+d}}(\Omega)},1\leq p\leq\infty,$

holds on the quasi-uniform mesh determined by the metric
$(\det{H})^{-\frac{1}{2p+d}}H$, where $H$ is a majorant of the Hessian matrix,
$N$ is the number of elements in the triangulation and the constant $C$ does
not depend on $u$ and $N$. This estimate is optimal in the sense that it is a
lower bound if $u$ is strictly convex or concave. Note that $\mathcal{H}$ can
be chosen as a majorant of the Hessian matrix.

When $m=1$, the new monitor function $M_{1,p}^{n}$ (1.6) is different with
(1.4) [26] that the former refers to $\mbox{tr}(\mathcal{H})$ and the latter
involves $\|\mathcal{H}\|$. In some cases, the two monitor functions are
pretty much alike. However, in other cases, the effect of the former is
superior to the latter for mesh generation. Numerical results in [40] have
shown our approach’s superiority for the error in $H^{1}$ norm.

Figure 2: Example 1: Interpolation error and its gradient in $L^{p}$ norm

##  4 Numerical experiments

In this section, we present some numerical results for three problems with
given analytical solutions. The numerical results are performed by using the
BAMG software [23]. Given a background mesh and an approximation solution,
BAMG generates the mesh according to the metric tensor. The code allows the
user to supply his/her own metric tensor defined on a background mesh. In our
computation, the background mesh has been taken as the most recent mesh
available.

Denote by $nbt$ the number of triangles in the current mesh. The number of
triangles is adjusted when necessary by trial and errors through the
modification of the multiplicative coefficient of the metric tensors.

Figure 3: Example 1: plots of the solution (a) and corresponding mesh (b)
using $\mathcal{M}_{1,1}$

Example 1 This example is to generate adaptive meshes for

$\displaystyle u({\bf
x})=\frac{1}{1+e^{-200(\sqrt{x_{1}^{2}+x_{2}^{2}}-0.8)}},\quad{\bf
x}\in(0.1,1)\times(0.1,1).$ (4.1)

This function is anisotropic along the quarter circle
$x_{1}^{2}+x_{2}^{2}=0.8^{2}$ and changes sharply in the direction normal to
this curve. A similar example was presented in [33] where the region is
$(0,1)\times(0,1)$. In the current computation, each run is stopped after 15
iterations to guarantee that the adaptive procedure tends towards stability.
We show in Figure 2 the $L^{p}$ norms of the interpolation error and its
gradient using corresponding metric tensors, for $p=1,2,4,\infty$. For
example, the curve $p=2$ in (a) stands for the interpolation error using the
metric tensor $\mathcal{M}_{0,2}$, while $p=\infty$ in (b) stands for the
gradient interpolation error using the metric tensor $\mathcal{M}_{1,\infty}$.
We see that the convergence rates for the interpolation error and its gradient
are always nearly optimal, i.e. $\|e\|_{L^{p}}\sim N^{-1}$ and $\|\nabla
e\|_{L^{p}}\sim N^{-0.5}$. We also show in Figure 3 plots of the solution and
corresponding mesh using the metric tensor $\mathcal{M}_{1,1}$.

Figure 4: Example 2: Interpolation error and its gradient in $L^{p}$ norm

Example 2 This example is to generate adaptive meshes for

$\displaystyle u({\bf
x})=x_{1}^{2}x_{2}+x_{2}^{3}+\tanh(10(\sin(5x_{2})-2x_{1})),\quad{\bf
x}\in(-1,1)\times(-1,1).$ (4.2)

This function is anisotropic along the zigzag curve $\sin(5x_{2})-2x_{1}=0$
and changes sharply in the direction normal to this curve (taken from [3]). In
the current computation, each run is stopped after 20 iterations to guarantee
that the adaptive procedure tends towards stability. We show in Figure 4 the
$L^{p}$ norms of the interpolation error and its gradient using corresponding
metric tensors, for $p=1,2,4,\infty$. As in Example 1, the convergence rates
for the interpolation error and its gradient here are always nearly optimal.
In Figure 5 we select 6 meshes with 4000 triangles generated by corresponding
metric tensors. We can learn that the optimal meshes in different norms are
different. For example, the mesh generated by the metric tensor
$\mathcal{M}_{1,\infty}$ concentrates more triangles and nodes along the
zigzag line.

Figure 5: Example 2: Meshes generated by the metric tensor $\mathcal{M}_{m,p}$
for (a)$m=0,p=1$, (b)$m=0,p=2$, (c)$m=0,p=\infty$, (d)$m=1,p=1$, (e)$m=1,p=2$,
(f)$m=1,p=\infty$.

Example 3 (Taken from [40]) This example is to solve the boundary value
problem of Poisson’s equation

$\displaystyle-\triangle u$ $\displaystyle=$ $\displaystyle f,\quad{\bf
x}\in\Omega\equiv(-1.2,1.2)\times(-1.2,1.2),$ (4.3)

with the Dirichlet boundary condition and the right-hand side term being
chosen such that the exact solution is given by

$\displaystyle u({\bf
x})=\sum_{i=1}^{5}\big{[}(1+e^{\frac{x+y-c_{i}}{2\epsilon}})^{-1}+(1+e^{\frac{x-y-
d_{i}}{2\epsilon}})^{-1}\big{]},$ (4.4)

where $c_{i}=0,-0.6,0.6,-1.2,1.2;\,d_{i}=0,-0.6,0.6,-1.2,1.2.$ The solution
exhibits ten sharp layers on lines $x+y-c_{i}=0$ and $x-y-d_{i}=0$,
$i=1,2,\cdots,5$, when $\epsilon$ is small. In our computations, $\epsilon$ is
taken as 0.01. Numerical results in [40] have shown that our approach’s
superiority for the error in $H^{1}$ norm. In the current computation, each
run is stopped after 20 iterations to guarantee that the adaptive procedure
tends towards stability, except that governed by $\mathcal{M}_{1,\infty}$,
which need 30 iterations. We show in Figure 6 the $L^{p}$ norms of the
interpolation error and its gradient using corresponding metric tensors, for
$p=1,2,4,\infty$. As in Example 1 and Example 2, the convergence rates for the
interpolation error and its gradient here are always nearly optimal. Another
purpose to select this example is to describe the difference of finding layers
using different norms. In Figure 7 we list meshes in different stage during
one selected run governed by corresponding metric tensors. While in Figure 8
convergence history is shown. From the three figures we can learn that most of
the metric tensors can quickly find the layers except the metric tensor
$\mathcal{M}_{1,\infty}$ when dealing with the complex problems, e.g., with
multiple layers.

Figure 6: Example 3: Interpolation error and its gradient in $L^{p}$ norm

Figure 7: Example 3: Meshes generated by the metric tensors
$\mathcal{M}_{1,2}$ after (a) 5 step, (b) 10 step, (c) 15 step, and
$\mathcal{M}_{1,\infty}$ after (d) 5 step, (e) 10 step, (f) 15 step, (g) 20
step.

Figure 8: Example 3: Convergence history versus number of triangles by the
metric tensor $\mathcal{M}_{0,p}$ and $\mathcal{M}_{1,p}$ for
$p=1,2,4,\infty$.

##  5 Conclusions

In the previous sections we have developed a uniform strategy to derive metric
tensors in two spatial dimension for interpolation errors and their gradients
in $L^{p}$ norm. The metric tensor $\mathcal{M}_{0,p}^{n}$ for the $L^{p}$
norm of the interpolation error is similar to some existing methods. However,
the metric tensor $\mathcal{M}_{1,p}^{n}$ is essentially different with those
metric tensors existed. There is a fine distinction between the new metric
tensor $\mathcal{M}_{1,p}^{n}$ and $\mathcal{M}_{1,p}^{h}$ proposed by Huang
and Russell[26] that the former refers to $\mbox{tr}(\mathcal{H})$ and the
latter involves $\|\mathcal{H}\|$. In some cases, the two metric tensors are
pretty much alike. However, when dealing with the complex problems, e.g., with
multiple layers, the effect of the former is superior to the latter for mesh
generation. Numerical results show that the corresponding convergent rates are
always almost optimal.

## References

  * [1] A. Agouzal, K. Lipnikov, Y. Vassilevski, Adaptive generation of quasi-optimal tetrahedral meshes, East-West J. Numer. Math. 7 (1999) 223 C244.
  * [2] A. Agouzal, K. Lipnikov, Y. Vassilevski, Hessian-free metric-based mesh adaptation via geometry of interpolation error, Comp. Math. Math. Phys. 50 (1) (2010) 124-138.
  * [3] A. Agouzal, Y. Vassilevski, Minimization of gradient errors of piecewise linear interpolation on simplicial meshes, Comput. Meth. Appl. Mech. Eng. 199 (2010) 2195-2203.
  * [4] D. Ait-Ali-Yahia, W. Habashi, A. Tam, M.-G. Vallet, M. Fortin, A directionally adaptive methodology using an edge-based error estimate on quadrilateral grids, Int. J. Numer. Methods Fluids 23 (1996) 673-690.
  * [5] T. Apel, G. Lube, Anisotropic mesh refinement in stabilized Galerkin methods, Numer. Math. 74(3) (1996) 261-282.
  * [6] R. Becker, An adaptive finite element method for the incompressible Navier-stokes equations on time-dependent domains, Ph.D. thesis, Ruprecht-Karls-Universit$\ddot{a}$t Heidelberg, 1995.
  * [7] H. Borouchaki, P.L. George, F. Hecht, P. Laug and E. Saltel, Delaunay mesh generation governed by metric specifications Part I. Algorithms, finite elem. anal. des. 25 (1997) 61-83.
  * [8] H. Borouchaki, P.L. George, B. Mohammadi, Delaunay mesh generation governed by metric specifications Part II. Applications, finite elem. anal. des. 25 (1997) 85-109.
  * [9] J.U. Brackbill, J.S. Saltzman, Adaptive zoning for singular problems in two dimensions, J. Comput. Phys. 46 (1982) 342-368.
  * [10] G. Buscaglia, E. Dari, Anisotropic Mesh Optimization and its Application in Adaptivity, Int. J. Numer. Meth. Eng. 40(22) (1997) 4119-4136.
  * [11] M.J. Castro-D$\acute{\imath}$az, F. Hecht, B. Mohammadi, O. Pironneau, Anisotropic unstructured mesh adaption for flow simulations, Internat. J. Numer. Methods Fluids 25(4) (1997) 475-491.
  * [12] L. Chen, P. Sun, J. Xu, Optimal anisotropic meshes for minimizing interpolation errors in $L^{p}$-norm, Math. Comp. 76(257) (2007) 179-204.
  * [13] E.F. D’Azevedo, R.B. Simpson, On optimal interpolation triangle incidences, SIAM J. Sci. Statist. Comput. 10 (1989) 1063-1075.
  * [14] E.F. D’Azevedo, Optimal triangular mesh generation by coordinate transformation, SIAM J. Sci. Stat. Comput. 12 (1991) 755-786.
  * [15] E.F. D’Azevedo, R.B. Simpson, On optimal triangular meshes for minimizing the gradient error, Numer. Math. 59 (1991) 321-348.
  * [16] W. Dorfler, R. Nochetto, Small data oscillation implies the saturation assumption, Numer. Math. 91 (2002) 1-12.
  * [17] A.S. Dvinsky, Adaptive grid generation from harmonic maps on Riemannian manifolds, J. Comput. Phys. 95 (1991) 450-476.
  * [18] P. Frey, P.L. George, Mesh Generation: Application to Finite Elements, Hermes Science, Oxford and Paris, 2000.
  * [19] R.V. Garimella, M.S. Shephard, Boundary layer meshing for viscous flows in complex domain. in: Proceedings of the 7th International Meshing Roundtable, Sandia National Laboratories, Albuquerque, NM, 1998, 107-118.
  * [20] P.L. George, F. Hecht. Nonisotropic grids, in: J.F. Thompson, B.K. Soni, N.P. Weatherill, (Eds.), Handbook of Grid Generation, CRC Press, Boca Raton, 1999 20.1-20.29.
  * [21] W.G. Habashi, J. Dompierre, Y. Bourgault, D. Ait-Ali-Yahia, M. Fortin, M.-G. Vallet, Anisotropic mesh adaptation: towards user-indepedent, mesh-independent and solver-independent CFD. Part I: general principles, Int. J. Numer. Meth. Fluids 32 (2000) 725-744.
  * [22] W.G. Habashi, M. Fortin, J. Dompierre, M.-G. Vallet, Y. Bourgault, Anisotropic mesh adaptation: a step towards a mesh-independent and user-independent CFD, Barriers and challenges in computational fluid dynamics (Hampton, VA, 1996), 99-117, Kluwer Acad. Publ., Dordrecht, 1998\.
  * [23] F. Hecht, Bidimensional anisotropic mesh generator, Technical Report, INRIA, Rocquencourt, 1997.
  * [24] W. Huang. Measuring mesh qualities and application to variational mesh adaptation. SIAM J. Sci. Comput. 26(5) (2005) 1643-1666.
  * [25] W. Huang, Metric tensors for anisotropic mesh generation, J. Comput. Phys. 204(2) (2005) 633-665.
  * [26] W. Huang and R.D. Russell, Adaptive Moving Mesh Methods, Series in Applied Mathematical Sciences, Springer, 2011.
  * [27] O.P. Jacquotte, A mechanical model for a new grid generation method in computational fluid dynamics, Comput. Meth. Appl. Mech. Eng. 66 (1988) 323-338.
  * [28] P. Knupp, L. Margolin, M. Shashkov, Reference jacobian optimization-based rezone strategies for arbitrary lagrangian eulerian methods, J. Comput. Phys. 176 (2002) 93-128.
  * [29] R. Kornhuber, R. Roitzsch, On adaptive grid refinement in the presence of internal or boundary layers, IMPACT Comput. Sci. Eng. 2 (1990) 40-72.
  * [30] J. Lang, An adaptive finite element method for convection-diffusion problems by interpolation techniques, Technical Report TR 91-4, Konrad-Zuse-Zentrum Berlin, 1991.
  * [31] R. Li, T. Tang, and P. Zhang, Moving mesh methods in multiple dimensions based on harmonic maps, J. Comput. Phys. 170(2) (2001) 562-588.
  * [32] E.J. Nadler, Piecewise linear approximation on triangulations of a planar region, Ph.D. Thesis, Division of Applied Mathematics, Brown University, Providence, RI, 1985.
  * [33] H. Nguyen, M. Gunzburger, L. Ju, J. Burkardt, Adaptive anisotropic meshing for steady convection-dominated problems, Comput. Meth. Appl. Mech. Eng. 198 (2009) 2964-2981.
  * [34] S. Owen, Meshing software survey, 1998, http://www.andrew.cmu.edu/user/sowen/softsurv.html.
  * [35] J. Peraire, M. Vahdati, K. Morgan, O.C. Zienkiewicz, Adaptive remeshing for compressible flow computation, J. Comp. Phys. 72(2) (1987) 449-466.
  * [36] W. Rachowicz, An anisotropic h-adaptive finite element method for compressible Navier-Stokes equations, Comput. Meth. Appl. Mech. Eng. 146 (1997) 231-252.
  * [37] J. Remacle, X. Li, M.S. Shephard, and J.E. Flaherty, Anisotropic adaptive simulation of transient flows using discontinuous Galerkin methods, Int. J. Numer. Meth. Eng., 62(7) (2005) 899-923.
  * [38] S. Rippa, Long and thin triangles can be good for linear interpolation, SIAM J. Numer. Anal. 29 (1992) 257-270.
  * [39] Y. Vassilevski, K. Lipnikov, Adaptive algorithm for generation of quasi-optimal meshes, Comp. Math. Math. Phys. 39 (1999) 1532 C1551.
  * [40] H. Xie, X. Yin, A strategy to derive metric tensors for anisotropic mesh generation, to appear.
  * [41] S. Yamakawa and K. Shimada, High quality anisotropic tetrahedral mesh generation via ellipsoidal bubble packing. in: Proceedings of the 9th International Meshing Roundtable, Sandia National Laboratories, Albuquerque, NM, 2000. Sandia Report 2000-2207, 263-273.
  * [42] X. Yin, H. Xie, A-posteriori error estimators suitable for moving mesh methods under anisotropic meshes, to appear.
  * [43] O.C. Zienkiewicz, J. Wu, Automatic directional refinement in adaptive analysis of compressible flows, Int. J. Numer. Meth. Eng. 37 (1994) 2189-2210.


###### Abstract

Scalarized black holes (BH) have been shown to form dynamically in extended-
scalar-tensor theories, either through spontaneous scalarization – when the BH
is unstable against linear perturbations – or through a non-linear
scalarization. In the latter, linearly stable BHs can ignite scalarization
when sufficiently perturbed. These phenomena are, however, not incompatible
and mixed scalarization is also possible. The objective of this work is
twofold: first, study mixed scalarization on a family of Einstein-Maxwell-
scalar models; and second, study the effect of the counter scalarization that
occurs when one of the coupling parameters has a sign opposite to the one that
generates scalarization. Both objectives are addressed by constructing and
examining the mixed scalarization’s domain of existence. An overall dominance
of the spontaneous scalarization over the non-linear scalarization is
observed. Thermodynamically, an entropical preference for mixed over the
standard scalarization (spontaneous or non-linear) exists. In the presence of
counter scalarization, a quench of the scalarization occurs, mimicking the
effect of a scalar particle’s mass/positive self-interaction term.

Mixed scalarization of charged black holes: from spontaneous to non-linear
scalarization

Zakaria Belkhadria1,2,3,4 and Alexandre M. Pombo5

1Département de Physique Théorique, Université de Genève, 24 quai Ernest
Ansermet, CH-1211 Geneva 4, Switzerland

2Gravitational Wave Science Center (GWSC), Université de Genève, CH-1211
Geneva, Switzerland

3Dipartimento di Matematica, Università di Cagliari, via Ospedale 72, 09124
Cagliari, Italy

4INFN, Sezione di Cagliari, Cittadella Universitaria, 09042 Monserrato, Italy

5CEICO, Institute of Physics of the Czech Academy of Sciences, Na Slovance 2,
182 21 Praha 8, Czechia

##  1 Introduction

With the recent observational data from the LIGO-VIRGO collaboration (e.g. [1,
2]) and the direct imaging by the Event Horizon Telescope [3, 4], a heightened
interest in alternative black hole (BH) solutions has emerged.

Of particular interest are hairy black hole solutions (see [5] for a review)
resulting from a scalarization process (aka scalarized BHs).

A possible mechanism of scalarization occurs in extended-scalar-tensor
theories (eST) [6]. In these, the scalar field non-minimally couples to the
model’s invariant and scalar perturbations of the vacuum BH can ignite the
growth of scalar hair around the BH. The resulting objects may present
significant deviations when compared with the standard vacuum General
relativity (GR) solutions, which can give a better insight into the nature of
gravity and particle physics.

A family of eST theories that have undergone extensive analysis is the
Einstein-Maxwell-Scalar (EMS) model [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
18, 19, 20, 21]. In the latter, scalarization is triggered by a non-minimal
coupling between a real scalar field, $\phi$, and the Maxwell invariant,
$\mathcal{I}=F_{\mu\nu}F^{\mu\nu}$, through a coupling function $f(\phi)$

$\mathcal{S}=\int
d^{4}x\sqrt{-g}\left[R-2\,\partial_{\mu}\phi\partial^{\mu}\phi-f(\phi)\mathcal{I}\right]\
,$ (1)

and minimally coupled with the Ricci scalar, $R$, associated with the metric
ansatz $g_{\mu\nu}$. Scalarization is triggered by a sufficiently large
charge-to-mass ratio.

Although BHs within this model are considered of less astrophysical
significance111In a dynamical astrophysical environment, the presence of
plasmas around the BH leads to prompt discharge. Alternatively, the
neutralization can occur through Hawking charge evaporation [22]. compared to
those in other eST models such as the Scalar-Gauss-Bonnet [23, 24, 25], their
relative computational simplicity has proven valuable to developing a better
insight into the scalarization phenomenon [26, 27, 28, 29, 30, 31, 32, 33, 34,
35, 36, 37, 38, 39, 40, 41].

In particular, the simplicity associated with the EMS model has allowed the
study of several coupling functions [8, 42, 43], from which two classes of
solutions were found222The same classification also exists for Scalar-Gauss-
Bonnet models [44, 45]: class I (or dilatonic-type) and class II (aka
scalarization type). While in the former $\phi=0$ does not solve the field
equation 333For the EMS model, there is an exceptional case: if $Q=P$,
$\phi=0$ solves this class so that the dyonic, equal charges Reissner-
Nordstrom BH is a solution., in the latter, the trivial scalar field does
solve the field equation and scalarized solutions can result from
perturbations of the vacuum BH solution 444A similar phenomenon was observed
for a vector field instead of a scalar field [46, 47], however, these seem to
be prone to ghost instabilities [48]. This demands that

$\frac{df(\phi)}{d\phi}\bigg{|}_{\phi=0}=0\ .$ (2)

This condition is naturally implemented, for instance, if one requires the
model to be $\mathbb{Z}_{2}$-invariant under $\phi\rightarrow-\phi$.
Scalarized solutions can be further divided into two sub-classes.

When vacuum BH are unstable against linear perturbations, a tachyonic
instability arises when

$\frac{d^{2}f(\phi)}{d\phi^{2}}\bigg{|}_{\phi=0}\neq 0\ ,$ (3)

and with the opposite sign of $\mathcal{I}$. The scalar hair around the BH
grows spontaneously from a perturbed vacuum Reissner-Nordstrom BH (RN BH). In
this case, the scalarized solutions bifurcate from the vacuum RN BH solutions.
These are known as class II.A or spontaneous/normal scalarization.

A second family of scalarized solutions, class II.B or non-linear
scalarization, is also possible when the BH is stable against linear
perturbations,

$\frac{d^{2}f(\phi)}{d\phi^{2}}\bigg{|}_{\phi=0}=0\ ,$ (4)

but unstable against non-linear perturbations. In this case, a sufficiently
large non-linear perturbation can trigger the growth of the scalar hair of the
vacuum BH.

These two scalarization types are, however, not incompatible and mixed
scalarization is also possible [44]. A coupling function compatible with both
spontaneous and non-linear scalarization is

$f(\phi)=e^{-\alpha\phi^{2}-\beta\phi^{4}}\ ,$ (5)

with $\alpha<0$ and $\beta<0$. The sign of the coupling parameters $(\alpha,\
\beta)$ was chosen in accordance with the literature. Pure – aka non-mixed –
spontaneous (non-linear) scalarization being recovered when $\beta=0$
($\alpha=0$).

The existence of scalarization is highly dependent on the sign of the
spontaneous/non-linear function parameters ($\alpha,\,\beta$). While for pure
scalarization the sign of either parameter is well defined, the presence of an
additional scalarization mechanism, say spontaneous scalarization ($\alpha$),
allows $\beta$ to have the “wrong” sign. We call this counter-scalarization.

The objective of this work is twofold: first, to study the interplay between
spontaneous and non-linear scalarization in the mixed scalarization scenario;
second, to examine the effects of incorporating a counter-scalarization term
into the coupling function.

The paper is organized as follows: Sec. 2 introduces the basics of the EMS
model, including a description of the equations of motion, boundary
conditions, and relevant relations. Sec. 3 is dedicated to the numerical
results. These include the computation of the domain of existence, Sec. 3.2,
and both the entropical, Sec. 3.2, and perturbative, Sec. 3.3, stabilities.
The paper concludes with some final remarks in Sec. 4.

Throughout this paper, we set $4\pi G=1=4\pi\epsilon_{0}$ for convenience. The
spacetime signature is chosen to be $(-,+,+,+)$. We focus exclusively on
spherically symmetric solutions, which implies that the metric and matter
functions depend solely on the radial coordinate. For simplicity in notation,
once a function is introduced with its radial dependency, such as $X(r)$, we
will subsequently denote it by $X$ with the understanding that it is a
function of $r$. Derivatives with respect to the radial coordinate $r$ and the
scalar field $\phi$ are represented by $X^{\prime}\equiv\frac{dX}{dr}$ and
$X_{,\phi}\equiv\frac{dX}{d\phi}$, respectively.

##  2 The EMS model

As already stated in the Introduction (Sec. 1), in this work we are going to
restrict ourselves to spherically symmetric EMS models described by action
(1). For the line element, let us consider a standard metric ansatz that is
compatible with spherical symmetry and has two unknown functions

$ds^{2}=-N(r)e^{-2\delta(r)}dt^{2}+\frac{dr^{2}}{N(r)}+r^{2}\big{(}d\theta^{2}+\sin^{2}\theta\,d\varphi^{2}\big{)}\
,~{}~{}{\rm with}~{}~{}N(r)\equiv 1-\frac{2\,m(r)}{r},$ (6)

where $m(r)$ is the Misner-Sharp mass function [49] and $\delta(r)$ is an
unknown metric function. Spherical symmetry, in the absence of a magnetic
charge555A magnetic charge would also be compatible with spherical symmetry,
but that shall not be considered here - see [43] for magnetically charged BHs
in this context., imposes an electrostatic $4-$vector potential,
$A(r)=V(r)\,dt$, and a scalar field solely radial dependent,
$\phi(t,r,\theta,\varphi)\equiv\phi(r)$. The absence of angular dependence
allows one to obtain the effective Lagrangian,

$\mathcal{L}_{\rm
eff}=e^{-\delta}m^{\prime}-\frac{1}{2}r^{2}e^{-\delta}N\phi^{\prime
2}+\frac{1}{2}f(\phi)e^{\delta}r^{2}V^{\prime 2}\ .$ (7)

Variation of the effective Lagrangian with respect to the metric and matter
functions yields the field equations:

$\displaystyle m^{\prime}=\frac{r^{2}N\phi^{\prime
2}}{2}+\frac{Q^{2}}{2r^{2}f(\phi)}\ ,\qquad\delta^{\prime}+r\phi^{\prime 2}=0\
,\qquad V^{\prime}=-\frac{Qe^{-\delta}}{f(\phi)r^{2}}\ ,$ (8)
$\displaystyle\phi^{\prime\prime}+\frac{1+N}{rN}\phi^{\prime}-\frac{Q^{2}}{r^{3}Nf(\phi)}\left(\phi^{\prime}-\frac{f_{,\phi}(\phi)}{2rf(\phi)}\right)=0\
.$

Where the electrostatic potential $V$ is under a first integral which was used
to simplify the remaining equations. The constant of integration is
interpreted as the electric charge, $Q$.

To solve the set of four coupled ordinary differential equations (8), one must
implement the appropriate boundary conditions. At the horizon, $r=r_{H}$, the
field equations can be approximated by a power series expansion in $r-r_{H}$
as

$\displaystyle
m=\frac{r_{H}}{2}+\frac{Q^{2}}{2r_{H}^{2}f(\phi_{0})}(r-r_{H})+\cdots\ \
,\qquad\delta=\delta_{0}-r_{H}\bigg{(}\frac{Q^{2}}{2r_{H}\big{(}Q^{2}-r_{H}^{2}f(\phi_{0})\big{)}}\frac{f_{,\phi}(\phi_{0})}{f(\phi_{0})}\bigg{)}^{2}\,(r-r_{H})+\cdots\
,$ (9)
$\displaystyle\phi=\phi_{0}+\frac{Q^{2}f_{,\phi}(\phi_{0})}{2r_{H}f(\phi_{0})\big{(}Q^{2}-r_{H}^{2}f(\phi_{0})\big{)}}(r-r_{H})+\cdots\
,\qquad\qquad\ \ \ \ \
V=-\frac{e^{-\delta_{0}}Q}{r_{H}^{2}f(\phi_{0})}(r-r_{H})+\cdots\ ,$

in terms of the two essential parameters $\phi_{0}$ and $\delta_{0}$, where
the subscript 0 denotes functions evaluated at the horizon $r_{H}$. At spatial
infinity, asymptotic flatness is ensured by a power series expansion in $1/r$.

$\displaystyle m(r)=M-\frac{Q^{2}+Q_{s}^{2}}{2r}+\cdots\
,\qquad\qquad\qquad\delta(r)\approx\frac{Q_{s}^{2}}{2\,r^{2}}+\cdots\ ,$ (10)
$\displaystyle\phi(r)=\frac{Q_{s}}{r}+\frac{MQ_{s}}{r^{2}}+\cdots\
,\qquad\qquad\qquad\quad V(r)=\psi_{e}+\frac{Q}{r}+\cdots\ ,$

with $M$ representing the ADM mass, $Q$ the BH’s electric charge, and $Q_{s}$
the scalar “charge”666The term scalar “charge” is used due to the similar
radial decay to a true electric charge, not because of an associated conserved
Noether current.777The so-called “scalar” hair associated with the EMS
scalarization model is of a secondary nature and does not add any additional
degree of freedom., while $\psi_{e}$ is the electrostatic potential at
infinity.

Equation (8), together with the boundary conditions arising from the power
series expansion (9) and (10), constitute a Dirichlet boundary condition
problem that must be numerically integrated (see Sec. 3).

###  2.1 Identities and physical quantities of interest

Scalarized solutions are physically characterised by the dimensionless
quantities: charge-to-mass ratio, $q$, reduced horizon area, $a_{H}$, and
reduced horizon temperature, $t_{H}$,

$q\equiv\frac{Q}{M}\ ,\qquad a_{H}\equiv\frac{A_{H}}{16\pi
M^{2}}=\frac{r_{H}^{2}}{4M^{2}}\ ,\qquad t_{H}\equiv 8\pi
MT_{H}=2MN^{\prime}(r_{H})e^{-\delta_{0}}\ \ ,$ (11)

where $A_{H}=4\pi r_{H}^{2}$ and $T_{H}=N^{\prime}(r_{H})e^{-\delta_{0}}/4\pi$
are the area and temperature of the BH’s horizon, respectively. Regularity of
the solutions is guaranteed by the Ricci scalar, $R$, and the Kretschmann
scalar, $K\equiv R_{\mu\nu\delta\lambda}R^{\mu\nu\delta\lambda}$,

$\displaystyle R$
$\displaystyle=\frac{N^{\prime}}{r}\big{(}3\,r\,\delta^{\prime}-4\big{)}+\frac{2}{r^{2}}\Big{\\{}1+N\Big{[}r^{2}\delta^{\prime\prime}-\big{(}1-r\delta^{\prime}\big{)}^{2}\Big{]}\Big{\\}}-N^{\prime\prime}\
,$ (12) $\displaystyle K$
$\displaystyle=\frac{4}{r^{4}}\big{(}1-N\big{)}^{2}+\frac{2}{r^{2}}\left[N^{\prime\,2}+\big{(}N^{\prime}-2N\delta^{\prime}\big{)}^{2}\right]+\left[N^{\prime\prime}-3\,\delta^{\prime}N^{\prime}+2N\big{(}\delta^{\prime\,2}-\delta^{\prime\prime}\big{)}\right]^{2}\
.$

Physical accuracy is ensured through the so-called virial identity, Smarr law
and non-linear Smarr relation. The virial identity is obtained through a
Derrick-like scaling argument [50, 51, 52, 53] and is given by

$\int_{r_{H}}^{\infty}dr\left\\{e^{-\delta}\,r^{2}\,\phi^{\prime\,2}\left[1+\frac{2\,r_{H}}{r}\bigg{(}\frac{m}{r}-1\bigg{)}\right]\right\\}=\int_{r_{H}}^{\infty}dr\left[e^{-\delta}\left(1-\frac{2\,r_{H}}{r}\right)\frac{1}{r^{2}}\frac{Q^{2}}{f(\phi)}\right]\
,$ (13)

which is independent of the equations of motion and displays that
scalarization can occur only in the presence of an electric charge $Q\neq 0$
in an EMS model 888Since $1+\frac{2r_{H}}{r}\Big{(}\frac{m}{r}-1\Big{)}>0$,
the left-hand side of the equation is strictly positive and can only be
counterbalanced by a non-zero electric charge in the right-hand side..

The Smarr law for this family of solutions is found not to be explicitly
dependent on the scalar field:

$M=\frac{1}{2}T_{H}A_{H}+\psi_{e}Q\ .$ (14)

The first law of black hole thermodynamics for EMS black holes is expressed as

$dM=\frac{1}{4}T_{H}\,dA_{H}+\psi_{e}\,dQ\ .$ (15)

At last, it can be shown that scalarized solutions also obey the so-called
non-linear Smarr formula [8],

$M^{2}+Q_{s}^{2}=Q^{2}+\frac{1}{4}A_{H}^{2}T_{H}^{2}\ .$ (16)

###  2.2 Coupling function

As previously mentioned, the choice of the coupling function is crucial in
determining whether the RN BH is susceptible to scalarization, and hence one
should be careful in designing the coupling function. Let us recap the
conditions for scalarization. The first requirement is that the GR BH solution
should also be a solution within the EMS model. Analysis of the field
equations (8) shows that this can be secured by imposing (2). This condition
is naturally implemented if one requires the model to be
$\mathbb{Z}_{2}-$invariant under the transformation $\phi\rightarrow-\phi$.
The type of scalarization, on the other hand, is controlled by the second
derivative of $f(\phi)$. For this, it is important to recall that the scalar
field is described by the Klein-Gordon equation

$\Box\phi=\frac{\mathcal{I}}{4}f_{,\phi}.$ (17)

Let us now consider a small-$\phi$ expansion of the coupling function

$f(\phi)=f(0)+\frac{1}{2}\frac{d^{2}f}{d\phi^{2}}\bigg{|}_{\phi=0}\phi^{2}+\cdots\
,$ (18)

the Klein-Gordon equation (17) linearized for small-$\phi$ reads:

$\big{(}\Box-\mu_{\rm eff}^{2}\big{)}\phi=0\ ,\qquad{\rm where}\qquad\mu_{\rm
eff}^{2}=\mathcal{I}\,\frac{d^{2}f(\phi)}{d\phi^{2}}\bigg{|}_{\phi=0}\ .$ (19)

The instability arises if the scalar field’s effective mass $\mu_{\rm
eff}^{2}<0$ (aka tachyonic mass), which, in particular, requires
$f_{,\phi\phi}$ to obey (3) and with the opposite sign of $\mathcal{I}$. This
constitutes the set of solutions known as normal/spontaneous scalarization
associated with class II.A, of which an exemplary function is

$f(\phi)=e^{-\alpha\phi^{2}}\ ,$ (20)

Class II.B on the other hand, occurs when solutions are always linearly stable
and no tachyonic instability exists, i.e. for (4). This condition is easily
satisfied by considering higher order terms in the expansion of $f(\phi)$ such
as

$f(\phi)=e^{-\beta\phi^{4}}\ .$ (21)

An exemplary function that exhibits, simultaneously, the tachyonic and non-
linear instabilities is the previously introduced (5),
$f(\phi)=e^{-\alpha\phi^{2}-\beta\phi^{4}}$. Observe that, due to the presence
of non-linear (linear) scalarization in the coupling function, scalarized BHs
in the EMS model can incorporate a counteracting linear (non-linear) term.

In particular, it is feasible to have scalarized solutions with $\alpha>0$
supported by the non-linear term $\beta<0$; or a scalarized solution with
$\beta>0$ supported by the tachyonic (linear) instability $\alpha<0$.

Analysis of (19) reveals that these counteracting terms do not contribute to
scalarization. Instead, they can be interpreted as a mass term for the scalar
field – in the case of the linear term $\alpha>0$ –, or as a positive quartic
self-interaction – for $\beta>0$ –, both proportional to the BH’s electric
charge due to the non-minimal coupling between the scalar field and the
Maxwell invariant. Investigation of the latter’s impact on the scalarization
phenomena constitutes the second objective of this work.

##  3 Numerical results

The set of four coupled ODEs (8), with the proper boundary conditions
(9)-(10), are solved numerically. The chosen routines automatically impose the
proper boundary conditions through a shooting method on the two unknown
parameters $\phi_{0}$ and $\delta_{0}$, with the maximum integration error and
boundary conditions automatically ensured to be less than $10^{-15}$.

Figure 1: Radial profile of the scalar field, $\phi$, metric, $(m,\,\delta)$,
and energy density, $\rho$, of a mixed scalarized BH solutions with
$r_{H}=0.78$, $\phi_{0}=0.40$ and $q=0.90$ and coupling parameters
$\alpha=-10.0$ and $\beta=-10.0\,$.

Physical accuracy is guaranteed through the virial identity, with a relative
error of $10^{-6}$, and Smarr and non-linear Smarr relations, with a relative
error of $10^{-7}$. The resulting numerical solution’s profile can be seen in
Fig. 1 for an exemplary solution with $r_{H}=0.78$, $\phi_{0}=0.40$ and
$q=0.90$ and coupling parameters $\alpha=-10.0$ and $\beta=-10.0\,$.

All obtained solutions are everywhere regular at and outside the horizon,
$r_{H}$. Each scalarized BH solution is uniquely defined by
$(r_{H},\,\alpha,\,\beta,\,q)$.

###  3.1 Domain of existence

Variation of the horizon radius, $r_{H}$, and coupling parameters,
$(\alpha,\,\beta)$, for a fixed electric charge $Q$, form the $3-$dimensional
domain of existence that characterizes the EMS model under analysis. To avoid
the additional complexity associated with the $3-$dimensional domain of
existence, let us fix one of the coupling parameters (say $\alpha$) and vary
the other ($\beta$) – see Fig. 2.

In Fig. 2 is graphically represented the projection of the domain of existence
of a scalarized BH as a function of $\alpha$ for three values of
$\beta=\\{-10.0,\,0,\,+10.0\\}$ (left panel); and as a function of $\beta$ for
four values of $\alpha=\\{-10.0,\,-1.0,\,0,\,+1.0\\}$ (right panel).

Figure 2: Domain of existence of scalarized BHs in EMS models (shaded blue
regions) with coupling function $f(\phi)=e^{-\alpha\phi^{2}-\beta\phi^{4}}$ as
a function of $\alpha$ with $\beta=\\{-10.0,\,0,\,+10.0\\}$ (left panel); or
$\beta$ for $\alpha=\\{-10.0,\,-1.0,\,0,\,+1.0\\}$ (right panel). The domain
of existence is bounded from below by the existence/turning line (solid/dot-
dashed line) and from above by the critical line (dashed line).

A common feature of both domains of existence is a region with $q\leqslant 1$
– where degeneracy occurs– and a region with overcharged solutions $q>1$. The
former region is limited from below by the bifurcation/turning line and from
above by the extremal line $q=1$, and represents a region of the parameter
space where, at least, two BH configurations with the same $q$ coexist: a bald
RN BH and one or more scalarized BH solution. These solutions, however,
possess different horizon radii (Fig. 3), temperatures and entropy (see Sec.
3.2).

The overcharged region goes from the extremal line, $q=1$, and is bounded from
above by the critical line, which is highly coupling parameter’s dependent:
$q_{crit.}\equiv q_{crit.}(\alpha,\beta)$. At the critical line, the
solution’s horizon radius tends to zero, $r_{H}\to 0$. In this region, no bald
RN BH exist and no degeneracy between scalarized BHs was observed. After the
critical line, no solutions exist.

In addition, Fig. 2 demonstrates a dominance of the tachyonic instability over
the non-linear scalarization: for the same value of $\beta$ and $\alpha$, the
domain of existence is bounded from below by the bifurcation line and has a
$q_{crit.}$ comparable to the pure spontaneous scalarization, $\beta=0.0\,$.
Such is expected since the maximum of the scalar field amplitude occurs at the
horizon (see Fig. 1) and is smaller than unity $max(phi)=\phi_{0}<1$,
resulting in a decreased contribution of the higher powers of the scalar
field999A similar behaviour was observed for the scalar-Gauss-Bonnet model
[44].

The lower bound, however, depends on the relative value of $\alpha$ in
relation to $\beta$. If $\beta\lesssim\alpha$, the scalarization is dominated
by the tachyonic instability and the domain of existence is bounded from below
by the bifurcation line; when $\beta\gg\alpha$, non-linear scalarization
dominates and the domain of existence is bounded from below by the turning
line (see Fig. 2). In this case, the tachyonic instability is not strong
enough to sustain a negligible amount of scalar field (see Fig. 3), and a
minimum value of $\phi_{0}$ exists. With the increase of $\alpha$, the size of
the initial $\phi_{0}$ jump decreases until $\phi_{0}\to 0$ and the tachyonic
instability dominates.

Let us now analyse the individual domains of existence. Starting with a fixed
$\beta=\\{-10.0,0.0,+10.0\\}$ and vary $\alpha$ – see Fig. 2 (left panel). In
the majority of the domain of existence, solutions bifurcate from the RN BH at
$q\leqslant 1$ – where scalarized BHs with a negligibly amount of scalar field
can exist – and stops at the critical line, $r_{H}\to 0$. In between a
continuously monotonic increase of the scalar field amplitude at the horizon
and a reduction of the latter occurs.

As expected, the bifurcation line is insensitive to the non-linear term since
the bifurcation is only dependent on the linear terms of $f_{,\phi}(\phi)$
that enter into the rhs of the Klein-Gordon equation (17).

The critical line, on the other hand, is highly dependent on the coupling
parameters. The maximum charge-to-mass ratio for which the critical solutions
exist, $q_{crit.}$, increases/decreases with the addition of a
negative/positive $\beta$ term. While for $\beta<0$, a non-linear instability
that amplifies the scalarization exists, a positive value has the opposite
effect. The $\beta>0$ decreases the width of the domain of existence due to a
decrease in $q_{crit.}$, resulting in a quench of the scalarization similar to
the one observed for scalarization with a positive self-interaction [54].

Observe now the case with fixed $\alpha=\\{-10,-1,0,+1\\}$ and varying $\beta$
– Fig. 2 (right panel). In this case, for all the chosen $\alpha\neq-10$,
solutions start at the extremal line $q=1$ for which a minimal, non-
negligible, amount of scalar field exists around the BH, $\phi_{0}\neq 0$. The
BH’s scalar field’s associated jump high decreases (increases) with the
addition of the $\alpha<0$ ($\alpha>0$) term (see also Fig. 3).

Figure 3: Scalar field horizon amplitude, $\phi_{0}$, as a function of the
charge-to-mass ratio, $q$, for a mixed scalarized BH for six different
coupling parameters configurations $(\alpha,\beta)$: (dot-dashed purple) mixed
scalarization $(-10.0,-10.0)$; (dotted green) pure spontaneous scalarization
$(-10.0,0.0)$; (dotted red) spontaneous scalarization $(-10.0,+10.0)$; (dashed
blue) mixed scalarization with $(-1.0,-10.0)$; (dot-dashed orange) pure non-
linear scalarization $(0.0,-10.0)$; mixed scalarization with $(+1.0,-10.0)$.

Following the initial jump in the scalar field amplitude at the horizon,
solutions observe a simultaneous increase of $\phi_{0}$ and $r_{H}$ until a
maximum mass, $r_{H}$ is reached – Fig. 3 $\beta=-10$ lines. After this point,
a second branch with decreasing $r_{H}$ and increasing $\phi_{0}$ exists until
the critical solution is achieved. The latter is known as the hot branch and
is known to be stable for $\alpha=0$ [55, 56]; the former is known as the cold
branch and is unstable (the denomination will be clear in the thermodynamics
section 3.2).

As mentioned before, while the effect of a positive $\alpha$ is similar to
scalarization by a massive scalar field, the positive $\beta$ mimics a
positive self-interacting (attractive) potential. The present results follow
the same pattern as the one presented for a massive/self-interacting scalar
field in EMS [54, 57] and scalar-Gauss-Bonnet models [58, 59, 45]. In
particular, a quench of the scalarization phenomena due to the decrease in the
domain of existence width (see Fig. 2). However, due to the non-minimal
coupling between the ”mass” or ”self-interaction” terms to the Maxwell
invariant, the impact of the latter is proportional to the electric charge
$Q$.

###  3.2 Thermodynamics

A solution is said to be stable if it is simultaneously entropically
preferable and stable against radial perturbations. The latter will be studied
in Sec. 3.3. In EMS models, entropy is given by the Bekeinstein-Hawking
formula [60, 61, 62] and reduces to the analysis of the reduced horizon area,
$a_{H}$, observe Fig. 4 (left panel).

Figure 4: Reduced area, $a_{H}$, (left panel) and reduced horizon temperature,
$t_{H}$, (right panel) as a function of the charge-to-mass ratio $q$ for a
mixed scalarized BH with six different coupling parameter’s configurations
$(\alpha,\beta)$: (dot-dashed purple) mixed scalarization $(-10.0,-10.0)$;
(dotted green) pure spontaneous scalarization $(-10.0,0.0)$; (dotted red)
spontaneous scalarization $(-10.0,+10.0)$; (dashed blue) mixed scalarization
with $(-1.0,-10.0)$; (dot-dashed orange) pure non-linear scalarization
$(0.0,-10.0)$; mixed scalarization with $(+1.0,-10.0)$.

Analysis of the horizon area shows that solutions dominated by spontaneous
scalarization are always entropically preferable when compared with electro-
vacuum GR; while non-linearly dominated solutions have a first branch that is
everywhere entropically unfavourable (aka cold branch), and a second branch
which contains a set of solutions that are entropically preferable (aka hot
branch). The second branch’s entropically non-preferable region decreases
(increases) with the addition of $\alpha<0$ ($\alpha>0$). The terminology
hot/cold comes from the second branch having a higher temperature than the
first, see Fig. 4 (right panel).

###  3.3 Comments on stability

Stability against radial perturbations is studied through a standard strategy
that considers spherically symmetric, linear perturbations of the equilibrium
solutions while keeping the metric ansatz (6), but allowing the functions
$N,\,\delta,\,\phi,\,V$ time, $t$, dependent besides $r$:

$ds^{2}=-\tilde{N}(r,t)\,e^{-2\,\tilde{\delta}(r,t)}dt^{2}+\frac{dr^{2}}{\tilde{N}(r,t)}+r^{2}\left(d\theta^{2}+\sin^{2}\theta\,d\varphi^{2}\right),\quad
A=\tilde{V}(r,t)\,dt,\quad\phi=\tilde{\phi}(r,t)\ .$ (22)

Each function can be further expanded into the equilibrium solution (aka
“bare”) that where obtained previously (Sec. 3-3.2), plus a perturbation term.
The latter contains the time-dependence through a Fourier mode with frequency
$\Omega$,

$\displaystyle\tilde{N}(r,t)=N(r)+\epsilon N_{p}(r)e^{-i\Omega t}\
,\quad\tilde{\delta}(r,t)=\delta(r)+\epsilon\delta_{p}(r)e^{-i\Omega t}\ ,$
(23) $\displaystyle\tilde{\phi}(r,t)=\phi(r)+\epsilon\phi_{p}(r)e^{-i\Omega
t}\ ,\quad\tilde{V}(r,t)=V(r)+\epsilon V_{p}(r)e^{-i\Omega t}\ ,$

where the subscript p denotes perturbations of the equilibrium solutions. The
linearized field equations around the background solution yield the metric and
$V_{p}(r)$ perturbations expressed in terms of the scalar field perturbation,
$\phi_{p}(r)$,

$N_{p}=-2\,r\,N\phi^{\prime}\,\phi_{p}\ ,\qquad\delta_{p}=-2\int
dr\,r\,\phi^{\prime}\phi_{p}^{\prime}\ ,\qquad
V_{p}^{\prime}=-V^{\prime}\bigg{[}\delta_{p}+\phi_{p}\frac{f_{,\phi}(\phi)}{f(\phi)}\bigg{]}\
,$ (24)

This leads to a single perturbation equation for $\phi_{p}(r)$ which can be
written in a Schrödinger-like form by redefining $\Psi(r)=r\,\phi_{p}$ and
inserting the ’tortoise’ coordinate $x$ defined by $dx/dr=e^{\delta}/N$:

$-\frac{d^{2}\Psi}{dx^{2}}+U_{\Omega}\Psi=\Omega^{2}\,\Psi\ ,$ (25)

with the perturbation potential $U_{\Omega}$ defined as:

$U_{\Omega}\equiv\frac{e^{-2\delta}N}{r^{2}}\left\\{1-N-2\,r^{2}\phi^{\prime\,2}-\frac{Q^{2}}{2r^{2}}\left[\frac{2}{f(\phi)}\left(1-2\,r^{2}\phi^{\prime
2}\right)-\frac{2f_{,\phi}^{2}(\phi)}{f^{3}(\phi)}+\frac{1}{f^{2}(\phi)}\Big{(}f_{,\phi\phi}(\phi)+4\,r\,\phi^{\prime}f_{,\phi}(\phi)\Big{)}\right]\right\\}$
(26)

The potential $U_{\Omega}$, which is regular across the entire domain
$-\infty<x<\infty$, diminishes to zero at both the black hole (BH) event
horizon and infinity. A mode is classified as unstable if $\Omega^{2}<0$,
which, within the asymptotic boundary conditions of our framework, indicates a
bound state. However, a standard quantum mechanical result dictates that
equation (25) will not exhibit bound states if $U_{\Omega}$ consistently
exceeds its minimal asymptotic value, implying it must be positive (see e.g.
[63]). Therefore, a uniformly positive effective potential serves as evidence
of mode stability against spherical perturbations.

Figure 5: Effective potential for a set of scalarized RN BH solutions. (Left
panel) Spontaneous scalarized dominated solutions with $\alpha=-10.0$ and
(dashed red) $\beta=-10.0$ mixed scalarization, (solid green) $\beta=0.0$ pure
spontaneous scalarization, and (dot-dashed blue) $\beta=+10.0$ counter-non-
linear. (Right panel) non-linear scalarization dominated solutions with
$\beta=-10.0$ for both hot (solid) and cold (dashed) branches with (red)
$\alpha=-1.0$ mixed scalarization, (black) $\alpha=0$ pure non-linear
scalarization, and (blue) $\alpha=+1.0$ counter-spontaneous scalarization.

To discern the interplay between the parameters $\alpha$ and $\beta$ on the
effective potential, we analyzed the profile of $U_{\Omega}$ for several
scalarized solutions, holding $\alpha$ constant while varying $\beta$ – Fig. 5
(left panel)–, and vice-versa – Fig. 5 (right panel).

The effective potential for radial spherical perturbations of spontaneously
scalarized dominated solutions corresponding to a fixed $\alpha=-10.0$ – Fig.
5 (left panel) –, exhibit an everywhere positive effective potential,
$U_{\Omega}>0$, suggesting the absence of instabilities. Additionally, one
observes an increase (decrease) of the maximum value of $U_{\Omega}$ with the
addition of the non-linear (counter-non-linear) parameter $\beta$.

In the case of non-linearly dominated scalarized solutions with fixed
$\beta=-10.0$ – Fig. 5 (right panel) –, for both the hot and cold branches, a
region where $U_{\Omega}$ changes sign exists. Which, while not indicating
instabilities, doesn’t guarantee stability.

The addition of a tachyonic term, $\alpha=-1.0$, to the non-linearly dominated
scalarized solution reduces the amplitude of the negative $U_{\Omega}$ region,
bringing it closer to the stable spontaneously lead scalarization, while a
counter-scalarized term $\alpha=+1.0$ deepens it.

At last, it is important to note that a non-uniformly positive potential for
the non-linear scalarized solutions is not a guarantee of instabilities. To
further investigate these solutions, the application of more intricate methods
like the S-deformation method [64, 65, 66] is required. Such is beyond the
scope of this paper (see [55] for a similar study of the pure non-linear
scalarization, and [54, 57] for the impact of self-interacting/mass term).

##  4 Conclusions

In this work, we have investigated the interplay between the spontaneous and
non-linear scalarization of charged black holes within the Einstein-Maxwell-
Scalar model. The resulting mixed scalarized solutions possess both properties
of pure non-linear and spontaneous scalarization.

The interplay between the two types of scalarization shows a dominance of the
spontaneus properties over the non-linear ones. In particular, the domain of
existence for comparable coupling parameters possesses the same structure of
the spontaneous scalarization with a slight increase in the width from the
non-linear scalarization. The influence of non-linear scalarization becomes
apparent only when its coupling parameter is significantly larger than that of
spontaneous scalarization.

The tachyonic instability associated with the spontaneous scalarization in the
mixed coupling makes the black hole more susceptible to scalarization. As a
result, it requires a weaker coupling of the scalar field to the Maxwell
invariant to achieve the same level of scalarization.

The presence of a mixed scalarization also allowed the study of a “counter-
scalarization´´ term, wherein one of the scalarization parameters has the
“wrong´´ sign and, instead of supporting/intensifying the scalarizaton,
suppresses it. From the linearized KG equation, one observes that these terms
possess the properties of a scalar field’s mass (for the linear parameter
$\alpha>0$) or of a positive self-interaction (for the non-linear $\beta>0$).
The resulting quench associated with these parameters mimics the effect of a
scalar field’s mass/self-interaction observed in previous studies [54, 57].
However, due to the coupling to the Maxwell invariant, the effect is not
constant.

Thermodynamically, mixed scalarized solutions with higher negative values of
both $\alpha$ and $\beta$ are favourable. Solutions dominated by spontaneous
scalarization are entropically preferable over their General Relativity
counterparts, while non-linear dominated solutions show mixed thermodynamic
behaviour, i.e. an entropically favourable and unfavourable regions.

Perturbative stability analysis against radial perturbations indicated
stability for dominant spontaneously scalarized black holes. In contrast, for
non-linear dominated scalarization, such a conclusion is not possible to be
made. Further studies into these must be made.

However, it seems like the addition of the spontaneous scalarization parameter
to the non-linear scalarization makes the resulting mixed scalarization tend
to a more stable configuration.

Future research directions include applying the S-Deformation method for an
in-depth analysis of the radial stability of non-linearly dominated solutions.
A comprehensive study of quasi-normal modes in mixed scalarized solutions
could provide further insights extending the work of [67, 56, 68, 69].
Additionally, incorporating rotation into these mixed scalarization
configurations presents an intriguing avenue for exploration.

## Acknowledgments

We would like to express our sincere gratitude to Eugen Radu for reading and
commenting on the manuscript and for the insightful discussions. Our thanks
also go to Pedro G.S. Fernandes and Nuno M. Santos for their valuable
discussions. Z.B extends special thanks to the Gravitation Group of Aveiro
(CIDMA) and to Carlos Herdeiro for their hospitality during the initial phase
of studying EMS models. Z.B also gratefully acknowledges the networking
support provided by the COST Action CA18108. A. M. Pombo is supported by the
Czech Grant Agency (GAĈR) under grant number 21-16583M.

## References

  * [1] B. P. Abbott, R. Abbott, T. Abbott, M. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. Adhikari, et al., “Observation of gravitational waves from a binary black hole merger,” Physical review letters, vol. 116, no. 6, p. 061102, 2016.
  * [2] R. Abbott, T. Abbott, S. Abraham, F. Acernese, K. Ackley, A. Adams, C. Adams, R. Adhikari, V. Adya, C. Affeldt, et al., “Gwtc-2: compact binary coalescences observed by ligo and virgo during the first half of the third observing run,” Physical Review X, vol. 11, no. 2, p. 021053, 2021.
  * [3] D. Ball, C.-k. Chan, P. Christian, B. T. Jannuzi, J. Kim, D. P. Marrone, L. Medeiros, F. Ozel, D. Psaltis, M. Rose, et al., “First m87 event horizon telescope results. i. the shadow of the supermassive black hole,” IOP PUBLISHING LTD, 2019.
  * [4] K. Akiyama, A. Alberdi, W. Alef, J. C. Algaba, R. Anantua, K. Asada, R. Azulay, U. Bach, A.-K. Baczko, D. Ball, et al., “First sagittarius a* event horizon telescope results. i. the shadow of the supermassive black hole in the center of the milky way,” The Astrophysical Journal Letters, vol. 930, no. 2, p. L12, 2022.
  * [5] C. A. Herdeiro and E. Radu, “Asymptotically flat black holes with scalar hair: a review,” International Journal of Modern Physics D, vol. 24, no. 09, p. 1542014, 2015.
  * [6] T. Damour and G. Esposito-Farese, “Nonperturbative strong-field effects in tensor-scalar theories of gravitation,” Physical Review Letters, vol. 70, no. 15, p. 2220, 1993.
  * [7] Y. S. Myung and D.-C. Zou, “Instability of reissner–nordström black hole in einstein-maxwell-scalar theory,” The European Physical Journal C, vol. 79, pp. 1–11, 2019.
  * [8] C. A. Herdeiro, E. Radu, N. Sanchis-Gual, and J. A. Font, “Spontaneous scalarization of charged black holes,” Physical review letters, vol. 121, no. 10, p. 101102, 2018.
  * [9] R. Konoplya and A. Zhidenko, “Analytical representation for metrics of scalarized einstein-maxwell black holes and their shadows,” Physical Review D, vol. 100, no. 4, p. 044015, 2019.
  * [10] Q. Gan, P. Wang, H. Wu, and H. Yang, “Photon ring and observational appearance of a hairy black hole,” Physical Review D, vol. 104, no. 4, p. 044049, 2021\.
  * [11] Q. Gan, P. Wang, H. Wu, and H. Yang, “Photon spheres and spherical accretion image of a hairy black hole,” Physical Review D, vol. 104, no. 2, p. 024003, 2021.
  * [12] G. Guo, P. Wang, H. Wu, and H. Yang, “Scalarized einstein–maxwell-scalar black holes in anti-de sitter spacetime,” The European Physical Journal C, vol. 81, pp. 1–14, 2021.
  * [13] Y. Peng, “Scalarization of horizonless reflecting stars: neutral scalar fields non-minimally coupled to maxwell fields,” Physics Letters B, vol. 804, p. 135372, 2020.
  * [14] Y. S. Myung and D.-C. Zou, “Scalarized black holes in the einstein-maxwell-scalar theory with a quasitopological term,” Physical Review D, vol. 103, no. 2, p. 024010, 2021.
  * [15] M. Khalil, N. Sennett, J. Steinhoff, and A. Buonanno, “Theory-agnostic framework for dynamical scalarization of compact binaries,” Physical Review D, vol. 100, no. 12, p. 124013, 2019.
  * [16] C.-Y. Zhang, Q. Chen, Y. Liu, W.-K. Luo, Y. Tian, and B. Wang, “Critical phenomena in dynamical scalarization of charged black holes,” Physical Review Letters, vol. 128, no. 16, p. 161105, 2022.
  * [17] G. Guo, P. Wang, H. Wu, and H. Yang, “Thermodynamics and phase structure of an einstein-maxwell-scalar model in extended phase space,” Physical Review D, vol. 105, no. 6, p. 064069, 2022.
  * [18] S. Hod, “Spontaneous scalarization of charged reissner-nordström black holes: Analytic treatment along the existence line,” Physics Letters B, vol. 798, p. 135025, 2019.
  * [19] S. Hod, “Analytic treatment of near-extremal charged black holes supporting non-minimally coupled massless scalar clouds,” The European Physical Journal C, vol. 80, pp. 1–7, 2020.
  * [20] S. Hod, “Spin-charge induced scalarization of kerr-newman black-hole spacetimes,” Journal of High Energy Physics, vol. 2022, no. 8, pp. 1–12, 2022.
  * [21] J. L. Blázquez-Salcedo, S. Kahlen, and J. Kunz, “Critical solutions of scalarized black holes,” Symmetry, vol. 12, no. 12, p. 2057, 2020.
  * [22] G. W. Gibbons, “Vacuum polarization and the spontaneous loss of charge by black holes,” Communications in Mathematical Physics, vol. 44, pp. 245–264, 1975.
  * [23] D. D. Doneva and S. S. Yazadjiev, “New gauss-bonnet black holes with curvature-induced scalarization in extended scalar-tensor theories,” Physical review letters, vol. 120, no. 13, p. 131103, 2018.
  * [24] H. O. Silva, J. Sakstein, L. Gualtieri, T. P. Sotiriou, and E. Berti, “Spontaneous scalarization of black holes and compact stars from a gauss-bonnet coupling,” Physical review letters, vol. 120, no. 13, p. 131104, 2018.
  * [25] G. Antoniou, A. Bakopoulos, and P. Kanti, “Evasion of no-hair theorems and novel black-hole solutions in gauss-bonnet theories,” Physical review letters, vol. 120, no. 13, p. 131102, 2018.
  * [26] C. A. Herdeiro, T. Ikeda, M. Minamitsuji, T. Nakamura, and E. Radu, “Spontaneous scalarization of a conducting sphere in maxwell-scalar models,” Physical Review D, vol. 103, no. 4, p. 044019, 2021.
  * [27] C.-Y. Zhang, P. Liu, Y.-Q. Liu, C. Niu, and B. Wang, “Dynamical charged black hole spontaneous scalarization in anti–de sitter spacetimes,” Physical Review D, vol. 104, no. 8, p. 084089, 2021.
  * [28] C.-Y. Zhang, Q. Chen, Y. Liu, W.-K. Luo, Y. Tian, and B. Wang, “Dynamical transitions in scalarization and descalarization through black hole accretion,” Physical Review D, vol. 106, no. 6, p. L061501, 2022.
  * [29] F. Yao, “Scalarized einstein–maxwell-scalar black holes in a cavity,” The European Physical Journal C, vol. 81, no. 11, p. 1009, 2021.
  * [30] G. Antoniou, A. Papageorgiou, and P. Kanti, “Probing modified gravity theories with scalar fields using black-hole images,” Universe, vol. 9, no. 3, p. 147, 2023.
  * [31] J.-Y. Jiang, Q. Chen, Y. Liu, Y. Tian, W. Xiong, C.-Y. Zhang, and B. Wang, “Type i critical dynamical scalarization and descalarization in einstein-maxwell-scalar theory,” arXiv preprint arXiv:2306.10371, 2023\.
  * [32] W.-K. Luo, C.-Y. Zhang, P. Liu, C. Niu, and B. Wang, “Dynamical spontaneous scalarization in einstein-maxwell-scalar models in anti–de sitter spacetime,” Physical Review D, vol. 106, no. 6, p. 064036, 2022.
  * [33] C. Niu, W. Xiong, P. Liu, C.-Y. Zhang, and B. Wang, “Dynamical descalarization in einstein-maxwell-scalar theory,” arXiv preprint arXiv:2209.12117, 2022\.
  * [34] M.-Y. Lai, Y. S. Myung, R.-H. Yue, and D.-C. Zou, “Spin-charge induced spontaneous scalarization of kerr-newman black holes,” Physical Review D, vol. 106, no. 8, p. 084043, 2022.
  * [35] W. Xiong, P. Liu, C. Niu, C.-Y. Zhang, and B. Wang, “Dynamical spontaneous scalarization in einstein-maxwell-scalar theory,” Chinese Physics C, vol. 46, no. 9, p. 095103, 2022.
  * [36] Q. Chen, Z. Ning, Y. Tian, X. Wu, C.-Y. Zhang, and H. Zhang, “Time evolution of einstein-maxwell-scalar black holes after a thermal quench,” Journal of High Energy Physics, vol. 2023, no. 10, pp. 1–26, 2023.
  * [37] C. Promsiri, T. Tangphati, E. Hirunsirisawat, and S. Ponglertsakul, “Scalarization of planar anti de sitter charged black holes in einstein-maxwell-scalar theory,” arXiv preprint arXiv:2302.04654, 2023\.
  * [38] F. Corelli, T. Ikeda, and P. Pani, “Challenging cosmic censorship in einstein-maxwell-scalar theory with numerically simulated gedanken experiments,” Physical Review D, vol. 104, no. 8, p. 084069, 2021.
  * [39] Q. Chen, Z. Ning, Y. Tian, B. Wang, and C.-Y. Zhang, “Nonlinear dynamics of hot, cold, and bald einstein-maxwell-scalar black holes in ads spacetime,” Physical Review D, vol. 108, no. 8, p. 084016, 2023.
  * [40] H. Xu and S.-J. Zhang, “Tachyonic instability of reissner-nordström-melvin black holes in einstein-maxwell-scalar theory,” Nuclear Physics B, vol. 987, p. 116110, 2023.
  * [41] G. Guo, P. Wang, H. Wu, and H. Yang, “Scalarized kerr-newman black holes,” Journal of High Energy Physics, vol. 2023, no. 10, pp. 1–20, 2023.
  * [42] P. G. Fernandes, C. A. Herdeiro, A. M. Pombo, E. Radu, and N. Sanchis-Gual, “Spontaneous scalarisation of charged black holes: coupling dependence and dynamical features,” Classical and Quantum Gravity, vol. 36, no. 13, p. 134002, 2019.
  * [43] D. Astefanesei, C. Herdeiro, A. Pombo, and E. Radu, “Einstein-maxwell-scalar black holes: classes of solutions, dyons and extremality,” Journal of High Energy Physics, vol. 2019, no. 10, pp. 1–27, 2019.
  * [44] D. D. Doneva and S. S. Yazadjiev, “Beyond the spontaneous scalarization: New fully nonlinear mechanism for the formation of scalarized black holes and its dynamical development,” Physical Review D, vol. 105, no. 4, p. L041502, 2022.
  * [45] A. M. Pombo and D. D. Doneva, “Effects of mass and self-interaction on nonlinear scalarization of scalar-gauss-bonnet black holes,” arXiv preprint arXiv:2310.08638, 2023.
  * [46] J. M. Oliveira and A. M. Pombo, “Spontaneous vectorization of electrically charged black holes,” Physical Review D, vol. 103, no. 4, p. 044004, 2021\.
  * [47] S. Barton, B. Hartmann, B. Kleihaus, and J. Kunz, “Spontaneously vectorized einstein-gauss-bonnet black holes,” Physics Letters B, vol. 817, p. 136336, 2021.
  * [48] L. Pizzuti and A. M. Pombo, “The spooky ghost of vectorization,” arXiv preprint arXiv:2310.18399, 2023.
  * [49] C. W. Misner and D. H. Sharp, “Relativistic equations for adiabatic, spherically symmetric gravitational collapse,” Physical Review, vol. 136, no. 2B, p. B571, 1964.
  * [50] C. A. Herdeiro, J. M. Oliveira, A. M. Pombo, and E. Radu, “Deconstructing scaling virial identities in general relativity: spherical symmetry and beyond,” Physical Review D, vol. 106, no. 2, p. 024054, 2022.
  * [51] J. M. Oliveira and A. M. Pombo, “A convenient gauge for virial identities in axial symmetry,” Physics Letters B, vol. 837, p. 137646, 2023.
  * [52] C. A. Herdeiro, J. M. Oliveira, A. M. Pombo, and E. Radu, “Virial identities in relativistic gravity: 1d effective actions and the role of boundary terms,” Physical Review D, vol. 104, no. 10, p. 104051, 2021.
  * [53] G. Derrick, “Comments on nonlinear wave equations as models for elementary particles,” Journal of Mathematical Physics, vol. 5, no. 9, pp. 1252–1254, 1964.
  * [54] P. G. Fernandes, “Einstein–maxwell-scalar black holes with massive and self-interacting scalar hair,” Physics of the Dark Universe, vol. 30, p. 100716, 2020.
  * [55] J. L. Blázquez-Salcedo, C. A. Herdeiro, J. Kunz, A. M. Pombo, and E. Radu, “Einstein-maxwell-scalar black holes: the hot, the cold and the bald,” Physics Letters B, vol. 806, p. 135493, 2020.
  * [56] J. L. Blázquez-Salcedo, C. A. Herdeiro, S. Kahlen, J. Kunz, A. M. Pombo, and E. Radu, “Quasinormal modes of hot, cold and bald einstein–maxwell-scalar black holes,” The European Physical Journal C, vol. 81, pp. 1–16, 2021.
  * [57] D.-C. Zou and Y. S. Myung, “Scalarized charged black holes with scalar mass term,” Physical Review D, vol. 100, no. 12, p. 124055, 2019.
  * [58] C. F. Macedo, J. Sakstein, E. Berti, L. Gualtieri, H. O. Silva, and T. P. Sotiriou, “Self-interactions and spontaneous black hole scalarization,” Physical Review D, vol. 99, no. 10, p. 104041, 2019.
  * [59] D. D. Doneva, K. V. Staykov, and S. S. Yazadjiev, “Gauss-bonnet black holes with a massive scalar field,” Physical Review D, vol. 99, no. 10, p. 104045, 2019.
  * [60] S. W. Hawking, “Particle creation by black holes,” Communications in mathematical physics, vol. 43, no. 3, pp. 199–220, 1975.
  * [61] J. D. Bekenstein, “Black holes and the second law,” in JACOB BEKENSTEIN: The Conservative Revolutionary, pp. 303–306, World Scientific, 2020.
  * [62] P. Majumdar, “Black hole entropy and quantum gravity,” arXiv preprint gr-qc/9807045, 1998.
  * [63] A. Messiah, Quantum Mechanics.  No. Chapter III2 in 1, (North Holland Publishing Company, 1961.
  * [64] M. Kimura, “A simple test for the stability of a black hole by s-deformation,” Classical and Quantum Gravity, vol. 34, no. 23, p. 235007, 2017.
  * [65] M. Kimura and T. Tanaka, “Robustness of the s-deformation method for black hole stability analysis,” Classical and Quantum Gravity, vol. 35, no. 19, p. 195008, 2018.
  * [66] M. Kimura and T. Tanaka, “Stability analysis of black holes by the s-deformation method for coupled systems,” Classical and Quantum Gravity, vol. 36, no. 5, p. 055005, 2019.
  * [67] Y. S. Myung and D.-C. Zou, “Quasinormal modes of scalarized black holes in the einstein–maxwell–scalar theory,” Physics Letters B, vol. 790, pp. 400–407, 2019.
  * [68] Y. S. Myung and D.-C. Zou, “Stability of scalarized charged black holes in the einstein–maxwell–scalar theory,” The European Physical Journal C, vol. 79, pp. 1–11, 2019.
  * [69] G. Guo, P. Wang, H. Wu, and H. Yang, “Quasinormal modes of black holes with multiple photon spheres,” Journal of High Energy Physics, vol. 2022, no. 6, pp. 1–24, 2022.