# Higher-Order Corrections to
Optimisers based on Newton’s Method

Stephen Brooks<EMAIL_ADDRESS>
###### Abstract

The Newton, Gauss–Newton and Levenberg–Marquardt methods all use the first
derivative of a vector function (the Jacobian) to minimise its sum of squares.
When the Jacobian matrix is ill-conditioned, the function varies much faster
in some directions than others and the space of possible improvement in sum of
squares becomes a long narrow ellipsoid in the linear model. This means that
even a small amount of nonlinearity in the problem parameters can cause a
proposed point far down the long axis of the ellipsoid to fall outside of the
actual curved valley of improved values, even though it is quite nearby. This
paper presents a differential equation that ‘follows’ these valleys, based on
the technique of geodesic acceleration, which itself provides a 2nd order
improvement to the Levenberg–Marquardt iteration step. Higher derivatives of
this equation are computed that allow $n^{\mathrm{th}}$ order improvements to
the optimisation methods to be derived. These higher-order accelerated methods
up to 4th order are tested numerically and shown to provide substantial
reduction of both number of steps and computation time.

††journal: Journal of Computational Physics

[SJB]organization=Collider–Accelerator Department,addressline=Brookhaven
National Laboratory, city=Upton, state=NY, postcode=11973, country=USA

##  1 Definitions and Introduction

Consider finding the value of a vector $\mathbf{x}$ such that the vector-
valued function $\mathbf{f}(\mathbf{x})=\mathbf{0}$, noting the input and
output of $\mathbf{f}$ might have different dimensions.

Newton’s method solves
$J(\mathbf{x})(\mathbf{x}_{\mathrm{new}}-\mathbf{x})=-\mathbf{f}(\mathbf{x})$
where $J$ is the Jacobian matrix. The Gauss–Newton algorithm generalises this
to rectangular $J$ using pseudo-inverses that may be calculated using Singular
Value Decomposition (SVD). The Levenberg–Marquardt algorithm [1, 2] introduces
a damping factor into this pseudo-inverse, which allows progress along
‘easier’ directions without having to go far in ‘difficult’ directions that
may exhibit nonlinearity.

The remainder of this paper will be written for the simpler Newton method,
where $J^{-1}$ is the inverse of the square Jacobian matrix. However, the
algorithms derived also work for the Gauss-Newton pseudo-inverse
$[J^{-1}]_{GN}=(J^{T}J)^{-1}J^{T}$ and the damped Levenberg–Marquardt version
$[J^{-1}]_{LM(\lambda)}=(J^{T}J+\lambda I)^{-1}J^{T}$. The latter is used in
numerical tests.

One common source of slow convergence is that $J$ is ill-conditioned, so the
optimisation valley is much narrower in some directions than others, while the
problem contains some nonlinearity, which may seem small but is amplified once
you change into coordinates where $J$ is well-conditioned. This is because
even a small amount of nonlinearity can make the long narrow valley stop
overlapping with its approximation in the linear model. This reduced range of
validity of the linear model means many small steps have to be taken.

##  2 Natural Optimisation Pathway

The goal of the Newton step is to reduce the error vector $\mathbf{f}$,
ideally to zero. For a nonlinear function, the optimisation follows a curved
pathway [3, 4] and one natural such pathway is $\mathbf{x}(t)$ defined
implicitly by

$\mathbf{f}(\mathbf{x}(t))=(1-t)\mathbf{f}(\mathbf{x}(0))$

for $t\in[0,1]$. This scales down all components of the error equally and at
$t=1$ it reaches the true solution.

Taking the first derivative of this equation gives

$\sum_{i}\partial_{i}\mathbf{f}(\mathbf{x}(t))\dot{x}_{i}(t)=J(\mathbf{x}(t))\dot{\mathbf{x}}(t)=-\mathbf{f}(\mathbf{x}(0))=-\frac{\mathbf{f}(\mathbf{x}(t))}{1-t},$

which at $t=0$ makes $\dot{\mathbf{x}}$ equal to the Newton step and to a
scaling of it for all $0<t<1$. So this pathway is always tangent to the Newton
step direction and corresponds to the limit of a Newton algorithm run with
steps scaled down to be infinitesimally small.

##  3 Higher-Order Derivatives

If the pathway curves, one may wonder if longer steps can be taken if the
curvature is taken into account. The second and higher derivatives of the
equation defining the natural pathway have the form

$\frac{\mathrm{d}^{n}}{\mathrm{d}t^{n}}\mathbf{f}(\mathbf{x}(t))=\mathbf{0}$

for $n\geq 2$. Multiple derivatives of a function composition
($\mathbf{f}\circ\mathbf{x}$ here) are given by Faà di Bruno’s formula [7, 8,
9]

$\frac{\mathrm{d}^{n}}{\mathrm{d}t^{n}}\mathbf{f}(\mathbf{x}(t))=\sum_{\pi\in\Pi_{n}}\mathbf{f}^{(|\pi|)}(\mathbf{x}(t))\bigotimes_{p\in\pi}\mathbf{x}^{(|p|)}(t),$

where $\Pi_{n}$ is the set of all partitions of $\\{1,2,...,n\\}$. The
$d^{\mathrm{th}}$ derivative of the vector function $\mathbf{f}$ is a tensor
that takes $d$ vectors as input and outputs a vector, with elements defined by

$f^{(d)}(\mathbf{x})^{i}_{j_{1}j_{2}...j_{d}}=\frac{\partial^{d}f_{i}(\mathbf{x})}{\partial
x_{j_{1}}\partial x_{j_{2}}...\partial x_{j_{d}}}.$

Note that $\mathbf{f}^{(1)}=J$. This paper will adopt compact notation where
tensor products of vectors $\mathbf{u}\otimes\mathbf{v}\otimes\mathbf{w}$ will
be written $\mathbf{u}\mathbf{v}\mathbf{w}$ so that
$(\mathbf{u}\mathbf{v}\mathbf{w})_{ijk}=u_{i}v_{j}w_{k}$. These may be
contracted with the derivative tensor to give a vector written in the form
$\mathbf{f}^{(3)}\mathbf{u}\mathbf{v}\mathbf{w}$, where
$(\mathbf{f}^{(3)}\mathbf{u}\mathbf{v}\mathbf{w})_{n}=\sum_{i,j,k}f^{(3)n}_{ijk}u_{i}v_{j}w_{k}$.

###  3.1 Second Order

For $n=2$, $\Pi_{2}=\\{\\{\\{1\\},\\{2\\}\\},\\{\\{1,2\\}\\}\\}$ and

$\frac{\mathrm{d}^{2}}{\mathrm{d}t^{2}}\mathbf{f}(\mathbf{x}(t))=\mathbf{f}^{(2)}(\mathbf{x}(t))\dot{\mathbf{x}}(t)\dot{\mathbf{x}}(t)+\mathbf{f}^{(1)}(\mathbf{x}(t))\ddot{\mathbf{x}}(t)=\mathbf{0}$
$\Rightarrow\qquad\ddot{\mathbf{x}}(t)=-J^{-1}(\mathbf{x}(t))\mathbf{f}^{(2)}(\mathbf{x}(t))\dot{\mathbf{x}}(t)\dot{\mathbf{x}}(t).$

This agrees with the well-known [3, 4, 5, 6] quadratic acceleration term for
Levenberg–Marquardt if the $J^{-1}$ is replaced by a damped pseudo-inverse
$[J^{-1}]_{LM(\lambda)}$.

###  3.2 Third Order

For conciseness, $\mathbf{x}$ and its derivatives will be evaluated at $t=0$
unless otherwise stated and $\mathbf{f}$ and its derivatives at $\mathbf{x}$.
For $n=3$,

$\Pi_{3}=\\{\\{\\{1\\},\\{2\\},\\{3\\}\\},\\{\\{1\\},\\{2,3\\}\\},\\{\\{2\\},\\{1,3\\}\\},\\{\\{3\\},\\{1,2\\}\\},\\{\\{1,2,3\\}\\}\\}$

and

$\frac{\mathrm{d}^{3}}{\mathrm{d}t^{3}}\mathbf{f}=\mathbf{f}^{(3)}\dot{\mathbf{x}}\dot{\mathbf{x}}\dot{\mathbf{x}}+3\mathbf{f}^{(2)}\dot{\mathbf{x}}\ddot{\mathbf{x}}+\mathbf{f}^{(1)}\mathbf{x}^{(3)}=\mathbf{0}.$

This gives the third derivative of $\mathbf{x}$ as

$\mathbf{x}^{(3)}=-J^{-1}(\mathbf{f}^{(3)}\dot{\mathbf{x}}\dot{\mathbf{x}}\dot{\mathbf{x}}+3\mathbf{f}^{(2)}\dot{\mathbf{x}}\ddot{\mathbf{x}}).$

###  3.3 Fourth Order, Recurrence and General Case

Higher-order expressions can be obtained either from the set partitions
$\Pi_{n}$ or the equivalent differentiation chain and product rules that
obtain $\frac{\mathrm{d}^{n+1}}{\mathrm{d}t^{n+1}}\mathbf{f}$ from
$\frac{\mathrm{d}^{n}}{\mathrm{d}t^{n}}\mathbf{f}$. The formulae

$\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{f}^{(n)}=\mathbf{f}^{(n+1)}\dot{\mathbf{x}}\qquad\mathrm{and}\qquad\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{x}^{(n)}=\mathbf{x}^{(n+1)}$

together with the product rule are enough to generate the full sequence.
Starting from $n=2$,

$\displaystyle\mathbf{f}^{(2)}\dot{\mathbf{x}}\dot{\mathbf{x}}+\mathbf{f}^{(1)}\ddot{\mathbf{x}}$
$\displaystyle=$ $\displaystyle\mathbf{0}$
$\displaystyle\mathbf{f}^{(3)}\dot{\mathbf{x}}\dot{\mathbf{x}}\dot{\mathbf{x}}+3\mathbf{f}^{(2)}\dot{\mathbf{x}}\ddot{\mathbf{x}}+\mathbf{f}^{(1)}\mathbf{x}^{(3)}$
$\displaystyle=$ $\displaystyle\mathbf{0}$
$\displaystyle\mathbf{f}^{(4)}\dot{\mathbf{x}}\dot{\mathbf{x}}\dot{\mathbf{x}}\dot{\mathbf{x}}+6\mathbf{f}^{(3)}\dot{\mathbf{x}}\dot{\mathbf{x}}\ddot{\mathbf{x}}+4\mathbf{f}^{(2)}\dot{\mathbf{x}}\mathbf{x}^{(3)}+3\mathbf{f}^{(2)}\ddot{\mathbf{x}}\ddot{\mathbf{x}}+\mathbf{f}^{(1)}\mathbf{x}^{(4)}$
$\displaystyle=$ $\displaystyle\mathbf{0}$

and so on. A computer algebra system can generate these terms based on a rule
like

$\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{f}^{(n)}\mathbf{x}^{(a)}\mathbf{x}^{(b)}\mathbf{x}^{(c)}=$
$\mathbf{f}^{(n+1)}\mathbf{x}^{(1)}\mathbf{x}^{(a)}\mathbf{x}^{(b)}\mathbf{x}^{(c)}+\mathbf{f}^{(n)}\mathbf{x}^{(a+1)}\mathbf{x}^{(b)}\mathbf{x}^{(c)}+\mathbf{f}^{(n)}\mathbf{x}^{(a)}\mathbf{x}^{(b+1)}\mathbf{x}^{(c)}+\mathbf{f}^{(n)}\mathbf{x}^{(a)}\mathbf{x}^{(b)}\mathbf{x}^{(c+1)}$

and collecting like terms, for example by sorting the $\mathbf{x}$ derivatives
in increasing order.

The highest derivative $\mathbf{x}^{(n)}$ may be moved to the other side to
get a formula like

$\mathbf{x}^{(4)}=-J^{-1}(\mathbf{f}^{(4)}\dot{\mathbf{x}}\dot{\mathbf{x}}\dot{\mathbf{x}}\dot{\mathbf{x}}+6\mathbf{f}^{(3)}\dot{\mathbf{x}}\dot{\mathbf{x}}\ddot{\mathbf{x}}+4\mathbf{f}^{(2)}\dot{\mathbf{x}}\mathbf{x}^{(3)}+3\mathbf{f}^{(2)}\ddot{\mathbf{x}}\ddot{\mathbf{x}}),$

shown for the $n=4$ case, which expresses it in terms of lower derivatives of
$\mathbf{x}$.

##  4 Taking Finite Steps

The derivatives $\mathbf{x}^{(n)}$ calculated above can produce a corrected
higher-order step using the Taylor series of $\mathbf{x}$ around $t=0$

$\mathbf{x}(\epsilon)=\sum_{n=0}^{\infty}\frac{1}{n!}\epsilon^{n}\mathbf{x}^{(n)},$

where the step is thought of as stopping at a time $t=\epsilon$ in the
parameterisation of the natural pathway. This unknown $\epsilon$ may seem like
a problem but it can be made to cancel. Define the correction at order $n$ to
be the $n^{\mathrm{th}}$ term of the Taylor series:

$\mathbf{c}_{n}=\frac{1}{n!}\epsilon^{n}\mathbf{x}^{(n)}.$

The step begins at $\mathbf{c}_{0}=\mathbf{x}$ and the first order uncorrected
step ends at $\mathbf{c}_{0}+\mathbf{c}_{1}$, so has length $\mathbf{c}_{1}$.
Now recall that for $n\geq 2$, the derivatives of $\mathbf{f}\circ\mathbf{x}$
are zero and use Faà di Bruno’s formula as before:

$\frac{\mathrm{d}^{n}}{\mathrm{d}t^{n}}\mathbf{f}(\mathbf{x}(t))=\sum_{\pi\in\Pi_{n}}\mathbf{f}^{(|\pi|)}(\mathbf{x}(t))\bigotimes_{p\in\pi}\mathbf{x}^{(|p|)}(t)=\mathbf{0}.$

Multiplying both sides by $\epsilon^{n}$ gives

$\sum_{\pi\in\Pi_{n}}\mathbf{f}^{(|\pi|)}(\mathbf{x}(t))\bigotimes_{p\in\pi}\epsilon^{|p|}\mathbf{x}^{(|p|)}(t)=\mathbf{0},$

using the fact that $\pi$ is a partition of $\\{1,2,...,n\\}$, so the sum of
sizes $|p|$ of all its elements is $n$. Noting that
$\epsilon^{n}\mathbf{x}^{(n)}=n!\mathbf{c}_{n}$ and evaluating at $t=0$ gives

$\sum_{\pi\in\Pi_{n}}\mathbf{f}^{(|\pi|)}\bigotimes_{p\in\pi}|p|!\mathbf{c}_{|p|}=\mathbf{0}.$

This formula is the basis for calculating corrections $\mathbf{c}_{n}$ for
finite steps in the following sections.

###  4.1 The Meaning of $\epsilon$

Observant readers might have noticed that
$\mathbf{c}_{1}=\epsilon\dot{\mathbf{x}}$ and in an earlier section,
$\dot{\mathbf{x}}=-J^{-1}\mathbf{f}$, so taking a full Newton step would imply
$\epsilon=1$. This paper treats $\epsilon$ as a small value because when
experiencing slow convergence from the ‘narrow curving valleys’ problem, the
area of validity for the local linear model (the trust region) is much smaller
than what is required to go all the way to the model minimum. This means the
steps taken that succeed in reducing the function sum of squares would only be
a fraction of the Newton step, for example a Levenberg–Marquardt step with
$\lambda$ chosen large enough to damp away the longest-range movement axes of
the exact Newton scheme.

##  5 Finite Difference Schemes

The higher-order corrections $\mathbf{c}_{n}$ are expressible in terms of
multiple directional derivatives of $\mathbf{f}$. For a numerical method,
these derivatives must be calculated from function values, or at most, the
Jacobian used by the algorithm. In this paper finite difference schemes are
used, some of which have their ‘stencils’ of sampled points spread in multiple
axes to give mixed derivatives.

###  5.1 Second Order

For $n=2$, the general formula gives

$\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{1}+\mathbf{f}^{(1)}2\mathbf{c}_{2}=\mathbf{0}$
$\Rightarrow\qquad\mathbf{c}_{2}=-\tfrac{1}{2}J^{-1}\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{1}.$

Taylor expansion of $\mathbf{f}$ in the direction $\mathbf{c}_{1}$ of the
original uncorrected step gives

$\mathbf{f}(\mathbf{x}+\mathbf{c}_{1})=\mathbf{f}+J\mathbf{c}_{1}+\tfrac{1}{2}\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{1}+O(\epsilon^{3})$
$\Rightarrow\qquad\tfrac{1}{2}\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{1}=\mathbf{f}(\mathbf{x}+\mathbf{c}_{1})-(\mathbf{f}+J\mathbf{c}_{1})+O(\epsilon^{3}).$

In other words, the difference between $\mathbf{f}(\mathbf{x}+\mathbf{c}_{1})$
and a linear estimate using the $\mathbf{f}$ and $J$ already calculated at
$\mathbf{x}(0)$, is to leading order a second derivative term similar to the
one required for calculating $\mathbf{c}_{2}$. Thus,

$\mathbf{c}_{2}=-J^{-1}(\mathbf{f}(\mathbf{x}+\mathbf{c}_{1})-(\mathbf{f}+J\mathbf{c}_{1}))+O(\epsilon^{3}).$

The evaluations required for this calculation are shown in Figure 1. In this
case, only one other point besides the evaluations of $\mathbf{f}$ and $J$ at
$\mathbf{x}$ is needed.

Figure 1: Finite difference stencil for calculating the second order
correction $\mathbf{c}_{2}$. Points represent evaluations of the function
$\mathbf{f}$ and rings represent evaluations of its Jacobian.

###  5.2 Third Order

For $n=3$, the general formula gives

$\mathbf{f}^{(3)}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{1}+3\mathbf{f}^{(2)}\mathbf{c}_{1}2\mathbf{c}_{2}+\mathbf{f}^{(1)}6\mathbf{c}_{3}=\mathbf{0}.$
$\Rightarrow\qquad\mathbf{c}_{3}=-\tfrac{1}{6}J^{-1}(\mathbf{f}^{(3)}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{1}+6\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{2}).$

There are a few differences from the second order case:

  * 1.

There is a third order derivative
$\mathbf{f}^{(3)}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{1}$, which will
require an additional stencil point in the direction of $\mathbf{c}_{1}$.

  * 2.

Errors will now have to be $O(\epsilon^{4})$ as the main terms have size
$O(\epsilon^{3})$.

  * 3.

There is a mixed derivative $\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{2}$,
requiring a two dimensional stencil pattern.

The mixed derivative requires knowledge of the direction $\mathbf{c}_{2}$,
which must be evaluated first. The second order stencil for $\mathbf{c}_{2}$
had error $O(\epsilon^{3})$ and now $O(\epsilon^{4})$ is needed, so even the
lower-order derivative $\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{1}$ will
have to be evaluated using a third order stencil. Fortunately, this stencil is
also needed for evaluating
$\mathbf{f}^{(3)}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{1}$, so all
coefficients of a cubic approximation to $\mathbf{f}$ in this direction can be
known.

####  5.2.1 Phase One: Calculating $\mathbf{c}_{2}$

The additional stencil point in the $\mathbf{c}_{1}$ direction is chosen to be
$\mathbf{x}+\frac{1}{2}\mathbf{c}_{1}$ here, although other choices are
possible. To third order,

$\mathbf{f}(\mathbf{x}+\tfrac{1}{2}\mathbf{c}_{1})=\mathbf{f}+\tfrac{1}{2}J\mathbf{c}_{1}+\tfrac{1}{8}\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{1}+\tfrac{1}{48}\mathbf{f}^{(3)}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{1}+O(\epsilon^{4})$
$\mathbf{f}(\mathbf{x}+\mathbf{c}_{1})=\mathbf{f}+J\mathbf{c}_{1}+\tfrac{1}{2}\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{1}+\tfrac{1}{6}\mathbf{f}^{(3)}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{1}+O(\epsilon^{4}).$

Writing the nonlinear part of $\mathbf{f}$ as
$\mathbf{f}_{nl}(\mathbf{x}+\mathbf{a})=\mathbf{f}(\mathbf{x}+\mathbf{a})-(\mathbf{f}+J\mathbf{a})$,
the derivatives in the $\mathbf{c}_{1}$ direction can be expressed

$\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{1}=16\mathbf{f}_{nl}(\mathbf{x}+\tfrac{1}{2}\mathbf{c}_{1})-2\mathbf{f}_{nl}(\mathbf{x}+\mathbf{c}_{1})+O(\epsilon^{4})$
$\mathbf{f}^{(3)}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{1}=12\mathbf{f}_{nl}(\mathbf{x}+\mathbf{c}_{1})-48\mathbf{f}_{nl}(\mathbf{x}+\tfrac{1}{2}\mathbf{c}_{1})+O(\epsilon^{4})$

and $\mathbf{c}_{2}$ calculated from the formula
$\mathbf{c}_{2}=-\frac{1}{2}J^{-1}\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{1}$
with $O(\epsilon^{4})$ error.

####  5.2.2 Phase Two: Calculating $\mathbf{c}_{3}$

This step requires the mixed derivative
$\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{2}$. Expressions of the form
$\mathbf{f}^{(3)}\mathbf{u}\mathbf{v}\mathbf{w}=(\mathbf{u}\cdot\nabla)(\mathbf{v}\cdot\nabla)(\mathbf{w}\cdot\nabla)\mathbf{f}$
are iterated directional derivatives. Each directional derivative can be
approximated to leading order as

$\displaystyle(\mathbf{u}\cdot\nabla)\mathbf{f}(\mathbf{x})$ $\displaystyle=$
$\displaystyle\frac{\mathbf{f}(\mathbf{x}+\epsilon\mathbf{u})-\mathbf{f}(\mathbf{x})}{\epsilon}+O(\epsilon)$
$\displaystyle\Rightarrow\qquad(\epsilon\mathbf{u}\cdot\nabla)\mathbf{f}(\mathbf{x})$
$\displaystyle=$
$\displaystyle\mathbf{f}(\mathbf{x}+\epsilon\mathbf{u})-\mathbf{f}(\mathbf{x})+O(\epsilon^{2})$
$\displaystyle\Rightarrow\qquad(\epsilon^{n}\mathbf{u}\cdot\nabla)\mathbf{f}(\mathbf{x})$
$\displaystyle=$
$\displaystyle\mathbf{f}(\mathbf{x}+\epsilon^{n}\mathbf{u})-\mathbf{f}(\mathbf{x})+O(\epsilon^{2n}).$

In the last formula above, $\epsilon^{n}\mathbf{u}$ represents an
$O(\epsilon^{n})$ sized term such as $\mathbf{c}_{n}$. Using this multiple
times allows mixed derivatives to be expressed to leading order as
combinations of function evaluations at different points (i.e. finite
difference stencils). For example,

$\displaystyle\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{2}$ $\displaystyle=$
$\displaystyle(\mathbf{c}_{1}\cdot\nabla)(\mathbf{c}_{2}\cdot\nabla)\mathbf{f}$
$\displaystyle\simeq$
$\displaystyle(\mathbf{c}_{1}\cdot\nabla)(\mathbf{f}(\mathbf{x}+\mathbf{c}_{2})-\mathbf{f}(\mathbf{x}))$
$\displaystyle\simeq$
$\displaystyle\mathbf{f}(\mathbf{x}+\mathbf{c}_{2}+\mathbf{c}_{1})-\mathbf{f}(\mathbf{x}+\mathbf{c}_{1})-(\mathbf{f}(\mathbf{x}+\mathbf{c}_{2})-\mathbf{f}(\mathbf{x})).$

Here, all terms have size $O(\epsilon^{3})$ and all approximations are leading
order accurate meaning the error is no worse than $O(\epsilon^{4})$, as
required.

In general, a $d^{\mathrm{th}}$ derivative of different directions would
require $2^{d}$ evaluations. An $n$ times repeated derivative in the same
direction only requires $n+1$ as some of the evaluation points are coincident.
A mixture like $\mathbf{f}^{(a+b+c)}\mathbf{u}^{\otimes a}\mathbf{v}^{\otimes
b}\mathbf{w}^{\otimes c}$ would require evaluation at $(a+1)(b+1)(c+1)$
points.

Some efficiencies may be gained from coincident evaluation points and the fact
that the full first derivative $J$ is usually evaluated at $\mathbf{x}$
already. This was used in the previous ‘second order’ section, which only
required one additional evaluation point for
$\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{1}$ rather than three.

Now everything is in place to evaluate
$\mathbf{c}_{3}=-\frac{1}{6}J^{-1}(\mathbf{f}^{(3)}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{1}+6\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{2})$.
The full step will be corrected from $\mathbf{c}_{1}$ to
$\mathbf{c}_{1}+\mathbf{c}_{2}+\mathbf{c}_{3}$. The evaluations required are
shown in Figure 2.

Figure 2: Finite difference stencil for calculating the third order correction
$\mathbf{c}_{3}$ along with the second order correction $\mathbf{c}_{2}$ that
is also required.

###  5.3 Fourth Order

For $n=4$, the general formula gives

$\mathbf{f}^{(4)}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{1}+6\mathbf{f}^{(3)}\mathbf{c}_{1}\mathbf{c}_{1}2\mathbf{c}_{2}+4\mathbf{f}^{(2)}\mathbf{c}_{1}6\mathbf{c}_{3}+3\mathbf{f}^{(2)}2\mathbf{c}_{2}2\mathbf{c}_{2}+\mathbf{f}^{(1)}24\mathbf{c}_{4}=\mathbf{0}$
$\Rightarrow\qquad\mathbf{c}_{4}=-\tfrac{1}{24}J^{-1}(\mathbf{f}^{(4)}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{1}+12\mathbf{f}^{(3)}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{2}+24\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{3}+12\mathbf{f}^{(2)}\mathbf{c}_{2}\mathbf{c}_{2}).$

As expected, there are more higher-order and mixed derivatives. The double
derivative $\mathbf{f}^{(2)}\mathbf{c}_{2}\mathbf{c}_{2}$ can take advantage
of the Jacobian to eliminate a point from the stencil, just as previous
unidirectional derivatives did. The direction $\mathbf{c}_{3}$ is now involved
in the derivatives, so three evaluation phases are required. All errors have
to be $O(\epsilon^{5})$ including those of $\mathbf{c}_{2}$ and
$\mathbf{c}_{3}$.

####  5.3.1 Phase One: Calculating $\mathbf{c}_{2}$

An additional stencil point $\mathbf{x}+\frac{3}{2}\mathbf{c}_{1}$ will be
added to increase the order of accuracy in the $\mathbf{c}_{1}$ direction.
Defining
$\mathbf{f}_{nl}(\mathbf{x}+\mathbf{a})=\mathbf{f}(\mathbf{x}+\mathbf{a})-(\mathbf{f}+J\mathbf{a})$
as before,

$\displaystyle\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{1}$
$\displaystyle\simeq$ $\displaystyle
24\mathbf{f}_{nl}(\mathbf{x}+\tfrac{1}{2}\mathbf{c}_{1})-6\mathbf{f}_{nl}(\mathbf{x}+\mathbf{c}_{1})+\tfrac{8}{9}\mathbf{f}(\mathbf{x}+\tfrac{3}{2}\mathbf{c}_{1})$
$\displaystyle\mathbf{f}^{(3)}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{1}$
$\displaystyle\simeq$
$\displaystyle-120\mathbf{f}_{nl}(\mathbf{x}+\tfrac{1}{2}\mathbf{c}_{1})+48\mathbf{f}_{nl}(\mathbf{x}+\mathbf{c}_{1})-8\mathbf{f}(\mathbf{x}+\tfrac{3}{2}\mathbf{c}_{1})$
$\displaystyle\mathbf{f}^{(4)}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{1}$
$\displaystyle\simeq$ $\displaystyle
192\mathbf{f}_{nl}(\mathbf{x}+\tfrac{1}{2}\mathbf{c}_{1})-96\mathbf{f}_{nl}(\mathbf{x}+\mathbf{c}_{1})+\tfrac{64}{3}\mathbf{f}(\mathbf{x}+\tfrac{3}{2}\mathbf{c}_{1}),$

all with errors $O(\epsilon^{5})$. These formulae came from writing out the
Taylor expansions and inverting the system of equations, which can also be
done by inverting a matrix as the equations are linear in $\mathbf{f}$ and its
derivatives. Calculate
$\mathbf{c}_{2}=-\frac{1}{2}J^{-1}\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{1}$
using the first formula above.

####  5.3.2 Phase Two: Calculating $\mathbf{c}_{3}$

The mixed derivative
$\mathbf{f}^{(3)}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{2}$ needs a grid with
three points in the $\mathbf{c}_{1}$ direction and two in the $\mathbf{c}_{2}$
direction for a total of six. Many previous points can be re-used, with the
only new point for fourth order in this phase being
$\mathbf{x}+\frac{1}{2}\mathbf{c}_{1}+\mathbf{c}_{2}$.

$\displaystyle\mathbf{f}^{(3)}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{2}$
$\displaystyle\simeq$
$\displaystyle(4\mathbf{f}(\mathbf{x}+\mathbf{c}_{2})-8\mathbf{f}(\mathbf{x}+\tfrac{1}{2}\mathbf{c}_{1}+\mathbf{c}_{2})+4\mathbf{f}(\mathbf{x}+\mathbf{c}_{1}+\mathbf{c}_{2}))-$
$\displaystyle(4\mathbf{f}-8\mathbf{f}(\mathbf{x}+\tfrac{1}{2}\mathbf{c}_{1})+4\mathbf{f}(\mathbf{x}+\mathbf{c}_{1}))$
$\displaystyle\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{2}$
$\displaystyle\simeq$
$\displaystyle(-3\mathbf{f}(\mathbf{x}+\mathbf{c}_{2})+4\mathbf{f}(\mathbf{x}+\tfrac{1}{2}\mathbf{c}_{1}+\mathbf{c}_{2})-\mathbf{f}(\mathbf{x}+\mathbf{c}_{1}+\mathbf{c}_{2}))-$
$\displaystyle(-3\mathbf{f}+4\mathbf{f}(\mathbf{x}+\tfrac{1}{2}\mathbf{c}_{1})-\mathbf{f}(\mathbf{x}+\mathbf{c}_{1}))$
$\displaystyle\mathbf{f}^{(2)}\mathbf{c}_{2}\mathbf{c}_{2}$
$\displaystyle\simeq$ $\displaystyle
2\mathbf{f}_{nl}(\mathbf{x}+\mathbf{c}_{2}).$

Again, all errors are $O(\epsilon^{5})$ and
$\mathbf{c}_{3}=-\frac{1}{6}J^{-1}(\mathbf{f}^{(3)}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{1}+6\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{2})$
can now be calculated to the required accuracy.

####  5.3.3 Phase Three: Calculating $\mathbf{c}_{4}$

This phase requires the single mixed derivative
$\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{3}$, which can be handled
analogously to when $\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{2}$ first
appeared, by extending the stencil in the $\mathbf{c}_{3}$ direction with the
two points $\mathbf{x}+\mathbf{c}_{3}$ and
$\mathbf{x}+\mathbf{c}_{1}+\mathbf{c}_{3}$.

$\displaystyle\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{3}$
$\displaystyle\simeq$
$\displaystyle\mathbf{f}(\mathbf{x}+\mathbf{c}_{1}+\mathbf{c}_{3})-\mathbf{f}(\mathbf{x}+\mathbf{c}_{3})-(\mathbf{f}(\mathbf{x}+\mathbf{c}_{1})-\mathbf{f}(\mathbf{x})).$

Now all values are available to evaluate the fourth order correction
$\mathbf{c}_{4}=-\tfrac{1}{24}J^{-1}(\mathbf{f}^{(4)}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{1}+12\mathbf{f}^{(3)}\mathbf{c}_{1}\mathbf{c}_{1}\mathbf{c}_{2}+24\mathbf{f}^{(2)}\mathbf{c}_{1}\mathbf{c}_{3}+12\mathbf{f}^{(2)}\mathbf{c}_{2}\mathbf{c}_{2})$.
The full step will be corrected from $\mathbf{c}_{1}$ to
$\mathbf{c}_{1}+\mathbf{c}_{2}+\mathbf{c}_{3}+\mathbf{c}_{4}$ and the
evaluations required are shown in Figure 3.

Figure 3: Finite difference stencil for calculating the second, third and
fourth order corrections $\mathbf{c}_{2}+\mathbf{c}_{3}+\mathbf{c}_{4}$ at
fourth order accuracy.

It is clear that this process could be continued to even higher orders,
although the stencils would require more and more points. Practically,
automated selection of points and calculations of the stencil coefficients
would also be required.

##  6 Numerical Test Problem

The higher-order algorithms were tested on a simple function to verify their
performance. This function had to exhibit the ‘narrow curving valleys’ problem
in its sum of squares, so a very anisotropic function $(x,y)\mapsto(x,Ky)$ for
$K\gg 1$ was chosen and then some nonlinearity added in parameter space. The
resulting function is

$\mathbf{f}(x,y)=(x+y^{2},K(y-x^{2})).$

Typical values of $K=10^{6}$ were used and the iteration was started from the
arbitrary point $(x,y)=(\pi,e)$, moving towards the minimum sum of squares at
$\mathbf{f}(0,0)=(0,0)$.

Figure 4 shows the improved performance of the higher order corrected methods
on the test problem with $K=10^{6}$. The error norm in the plot is the value
of $|\mathbf{f}(x,y)|$ after each iteration. There is a slow convergence
region for $0.1\leq|\mathbf{f}|\leq 10$ after which the algorithm converges
rapidly. When the full Newton step using the linear model becomes a valid,
convergence should be quadratic with the error norm roughly squaring on each
iteration. This rapid convergence appears as the near-vertical descending
lines on the graph for $|\mathbf{f}|<0.1$.

Figure 4: Performance of higher-order corrected Levenburg–Marquardt methods on
a test problem.

###  6.1 Varying Valley Anisotropy

Varying $K$ should show the relative performance of the different order
algorithms as the valley gets narrower, while the curvature is kept constant.
Table 1 shows the number of iterations required to converge as $K$ is varied
between $1$ and $10^{12}$.

Table 1: Test problem convergence times for different order methods as the anisotropy factor $K$ is varied. Anisotropy $K$ | 1st order | 2nd order | 3rd order | 4th order
---|---|---|---|---
1 | 8 | 6 | 5 | 5
10 | 15 | 8 | 6 | 5
100 | 47 | 16 | 9 | 8
1000 | 196 | 30 | 18 | 11
10000 | 880 | 68 | 24 | 18
100000 | 4041 | 162 | 50 | 27
$10^{6}$ | 18733 | 397 | 88 | 43
$10^{7}$ |  $>$20000 | 971 | 166 | 70
$10^{8}$ |  $>$20000 | 2432 | 312 | 110
$10^{9}$ |  $>$20000 | 5828 | 631 | 243
$10^{10}$ |  $>$20000 |  $>$20000 | 2876 | 968
$10^{11}$ |  $>$20000 |  $>$20000 | 10886 | 2706
$10^{12}$ |  $>$20000 |  $>$20000 |  $>$20000 | 9159

A simplified model is that the valley has width $1/K$ while the
$n^{\mathrm{th}}$ order method has error term $O(\epsilon^{n+1})$, so this
error would push the proposed step out of the valley when
$O(\epsilon^{n+1})=1/K$ or $\epsilon=O(K^{-\frac{1}{n+1}})$. This would give
$O(1/\epsilon)=O(K^{\frac{1}{n+1}})$ steps to convergence.

Figure 5: Test problem convergence times for different order methods as the
anisotropy factor $K$ is varied.

Plotting the data on a log-log plot in Figure 5 reveals straight lines in
parts of the data that suggest a power law relationship. For $K\geq 10^{9}$
there is an additional increase in convergence time, which may be from the
limits of double precision used for the calculation. Taking the gradient
through the last three available points with $K\leq 10^{8}$ gives power law
exponents of $0.660,0.392,0.265,0.203$ for the first through fourth order
methods. This is somewhat similar to the simplified model’s
$\frac{1}{2},\frac{1}{3},\frac{1}{4},\frac{1}{5}$ although lower-order methods
seem slower, particularly the first order.

###  6.2 Integration with Levenburg–Marquardt Method

These numerical experiments were performed with a Levenburg–Marquardt method
enhanced with the higher-order corrections. Step length $\epsilon$ is
controlled by choosing the damping parameter $\lambda\geq 0$ in the pseudo-
inverse $[J^{-1}]_{LM(\lambda)}=(J^{T}J+\lambda I)^{-1}J^{T}$. The note [10]
shows that $\lambda=0$ gives the full Gauss–Newton step, while
$\lambda\rightarrow\infty$ produces infinitesimal steepest gradient steps,
with the values of $\lambda$ in between producing optimal reductions in the
linear model for a given step size.

A good choice of $\lambda$ should be chosen at each step. In this study the
values

$\lambda_{n}=\lambda_{\mathrm{old}}10000^{(n/10)^{3}}\qquad\mbox{for}\quad-10\leq
n\leq 10,$

where $\lambda_{\mathrm{old}}$ is the value from the previous step, are run in
parallel and the one that produces the best reduction in $|\mathbf{f}|$
chosen. The initial step uses $\lambda_{\mathrm{old}}=1$.

The fact these steps are run in parallel means that for the higher order
methods, many initial Levenburg–Marquardt steps $\mathbf{c}_{1}(\lambda_{n})$
are calculated, each of which has higher order corrections
$\mathbf{c}_{i}(\lambda_{n})$. The function $\mathbf{f}$ is evaluated at all
the corrected step points
$\mathbf{x}_{\mathrm{out}}(\lambda_{n})=\sum_{i=0}^{order}\mathbf{c}_{i}(\lambda_{n})$
and the one with lowest $|\mathbf{f}|$ and its corresponding $\lambda_{n}$ is
chosen.

The higher order methods enable longer steps and thus smaller values of
$\lambda$ to be used. The scheme above is somewhat wasteful by trying 21
values of $\lambda$ each step, but on modern computers these can be
parallelised, unlike the slow progress along the narrow optimisation valley,
which is a serial calculation.

##  7 Performance on a Practical Problem

These algorithms have also been used on a more complex optimisation problem
(which motivated their development). This problem has 180 parameters and 300
output variables and features levels of successively more difficult narrow
curved valleys in its optimisation space.

The full details of this problem are not the point of this paper but a brief
summary will be given here. An initial distribution of 100 Ca+ ions is
accelerated through a potential of 1 kV and given a $\pm$2% energy chirp. It
is transported through a curved electrostatic channel, where the electric
field is produced by 15 rings of 12 configurable electrodes. Each electrode is
modelled as a point charge and these 180 charges are the optimisation
variables. The output vector whose norm should be minimised contains the
$(x,y,z)$ position coordinates of the 100 ions on exiting the channel (so 300
entries in all), with the bunch centroid subtracted. The ions do not interact
in this model and their trajectories are calculated by the 4th order
Runge–Kutta method [11, 12] with a timestep $\delta t=10^{-7}$ seconds.

This problem is interesting because focussing the ions to a point in the
linear dynamics approximation can be done with standard optics, but optical
aberrations at higher order will remain, for example spherical aberration from
large angle effects. These higher order aberrations can also be corrected by
careful choice of the electrode voltages, although this gets more difficult
the smaller the focal point becomes and the more aberrations have to be
cancelled simultaneously.

The figure of merit is the RMS focal size of the ion bunch, which is
$\frac{1}{\sqrt{100}}$ of the norm $|\mathbf{f}|$. Figure 6 shows how this is
reduced by the Levenburg–Marquardt optimisation method with various levels of
higher order correction.

Figure 6: Performance of higher-order corrected Levenburg–Marquardt methods on
a physics problem with 180 parameters and 300 output variables.

Figure 6 shows that higher order methods signficantly accelerate the
optimisation progress while the ion bunch focal size is greater than $10^{-7}$
metres, with comparative performance analogous to that on the test problem in
Figure 4. Once the focal size reduces below $10^{-7}$ metres, progress in the
optimisation slows down greatly for all orders of method reaching this level.
The reason for this slow-down is yet to be determined. Lack of numerical
precision would produce a similar ‘noise floor’, although in this study care
was taken to calculate the Jacobian $J$ with automatic differentiation rather
than the noisier finite difference schemes.

Figure 7: Calculation time vs. performance for higher-order methods on a
physics problem with 180 parameters.

One potential concern with these higher-order methods is that the additional
evaluations of $\mathbf{f}$ in the stencil will cost too much time to make the
method worth using. To measure this effect, the optimised focal size is
plotted as a function of calculation time in Figure 7. The differences in
total execution time can be seen at the right-hand end of each line, which is
after 1000 iterations. The fourth order method takes 44% longer per step but
still manages to pull ahead of the other methods in real time before the
‘floor’ is reached. Table 2 gives the amount of time for each method to reach
an RMS focal size of under 1 micron.

Table 2: Real time taken to reach a focal size of $<10^{-6}$ m in the ion focussing problem with 180 parameters. Method Order | Iterations | Calculation Time (s)
---|---|---
1 | 532 | 3604.863
2 | 94 | 665.812
3 | 62 | 495.010
4 | 46 | 418.301

##  8 Conclusion

Methods such as Levenburg–Marquardt are not only for curve fitting, but also
powerful optimisers where the function to be minimised is a sum of squares.
Like many optimisers, they can get stuck for long periods of time in ‘curved
narrow valleys’. This paper derives higher-order corrections beyond the
already known second order [3, 4] that further accelerate the optimiser
performance in these situations. A general formula for deriving the
$n^{\mathrm{th}}$ order correction is given, with suggestions on how to build
finite difference stencils to evaluate it.

These successful methods have been derived using the concept of a ‘natural
pathway’ for the optimisation, which is an ordinary differential equation
(ODE) that is meant to follow the valleys. The form of this ODE is chosen
somewhat arbitrarily here but it appears to work well, perhaps because it is a
continuous version of the optimiser’s path in the limit where step size
$\epsilon\rightarrow 0$. Using a higher-order step method on this ODE thus
makes the optimiser behave ‘as if’ it had done a large number of very small
steps.

This link with ODEs also suggests potential future work in applying well-known
higher-order methods for ODEs such as Runge–Kutta [11, 12] or Bulirsch–Stoer
[13] to difficult optimisation problems, as well as methods for stiff ODEs. In
this paper the RK4 method was not preferred because it would have required
four evaluations of the Jacobian, which is still much more expensive than the
eight additional function points in the stencil of the fourth order method.

## References

  * [1] “A Method for the Solution of Certain Non-Linear Problems in Least Squares”, Kenneth Levenberg, Quarterly of Applied Mathematics 2 (2), pp.164–168 (1944) doi:10.1090/qam/10666
  * [2] “An Algorithm for Least-Squares Estimation of Nonlinear Parameters”, Donald Marquardt, SIAM Journal on Applied Mathematics 11 (2), pp.431–441 (1963) doi:10.1137/0111030
  * [3] “Why are Nonlinear Fits to Data so Challenging?”, Mark K. Transtrum, Benjamin B. Machta and James P. Sethna, Phys. Rev. Lett. 104, 060201 (2010) doi:10.1103/PhysRevLett.104.060201
  * [4] “Geometry of nonlinear least squares with applications to sloppy models and optimization”, Mark K. Transtrum, Benjamin B. Machta, and James P. Sethna, Phys. Rev. E 83, 036701 (2011) doi:10.1103/PhysRevE.83.036701
  * [5] “Improvements to the Levenberg-Marquardt algorithm for nonlinear least-squares minimization”, Mark K. Transtruma and James P. Sethna, arXiv:1201.5885v1 [physics.data-an] (2012).
  * [6] “Geodesic acceleration and the small-curvature approximation for nonlinear least squares”, Mark K. Transtruma and James P. Sethna, arXiv:1207.4999v1 [math.OC] (2012).
  * [7] “Du calcul des derivations” [On the calculus of derivatives] (in French), L.F.A. Arbogast, Strasbourg: Levrault, pp. xxiii+404 (1800). Entirely freely available from Google books https://books.google.com/books?id=YoPq8uCy5Y8C
  * [8] “Sullo sviluppo delle funzioni” [On the development of the functions] (in Italian), F. Faà di Bruno, Annali di Scienze Matematiche e Fisiche, 6: 479–480, LCCN 06036680 (1855). Entirely freely available from Google books https://books.google.com/books?id=ddE3AAAAMAAJ&pg=PA479
  * [9] “Note sur une nouvelle formule de calcul differentiel” [On a new formula of differential calculus] (in French), F. Faà di Bruno, The Quarterly Journal of Pure and Applied Mathematics, 1: 359–360 (1857). Entirely freely available from Google books https://books.google.com/books?id=7BELAAAAYAAJ&pg=PA359
  * [10] “Bounded Approximate Solutions of Linear Systems using SVD”, Stephen Brooks, tech note BNL-223624-2022-TECH, available from https://technotes.bnl.gov/ or https://stephenbrooks.org/ap/report/2015-3/svdboundedsolve.pdf (2015).
  * [11] Runge, C., _Math. Ann._ 46, p.167 (1895); and Kutta, M.W.Z., _für Math. u. Phys._ 46, p.435 (1901).
  * [12] _Numerical Recipes in C_ , W.H. Press, B.P. Flannery, S.A. Teukolsky and W.T. Vettering, chapter 15.1: Integration of Ordinary Differential Equations - Runge–Kutta Method.
  * [13] _Introduction to Numerical Analysis_ , J. Stoer and R. Bulirsch, section 7.2.14, New York: Springer-Verlag (1980).


# Dynamical scalarization in Einstein-Maxwell-dilaton theory

Cheng-Yong Zhang1, Peng Liu1, Yunqi Liu2, Chao Niu1, Bin Wang2,3
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We study the process of fully nonlinear dynamical scalarization starting from
a charged black hole or a naked singularity in asymptotically flat spacetime
in the Einstein-Maxwell-dilaton theory. Initially the dilaton field is
negligible compared to the gravitational and the Maxwell field. Then the
dilaton field experiences an immediate growth, later it oscillates with
damping amplitude and finally settles down to a finite value. For a hairy
black hole develops from an original Reissner-Nordström black hole, since the
dilaton oscillation and decay are almost independent of the coupling
parameter, unlike the Anti-de Sitter spacetime it is not easy to distinguish
the resulting hairy black hole from the original asymptotically flat charged
hole. For a hairy black hole evolves from an original naked singularity, the
resulting hairy black hole has rich structures. In the scalarization process,
the naked singularity is soon enveloped by one outer horizon, then another
horizon is developed and in the end a stable hairy black hole forms and two
horizons degenerate into one to protect the singularity. The hairy black hole
mass saturates exponentially in the scalarization.

1\. Department of Physics and Siyuan Laboratory, Jinan University, Guangzhou
510632, China
2\. Center for Gravitation and Cosmology, College of Physical Science and
Technology, Yangzhou University, Yangzhou 225009, China
3\. School of Aeronautics and Astronautics, Shanghai Jiao Tong University,
Shanghai, China

##  1 Introduction

The Einstein-Maxwell-dilaton (EMD) theory origins from the Kaluza-Klein
compactification [1] and also appears as the low energy limit of string theory
and is ubiquitous in supergravity [2, 3]. In Jordan frame, dropping all the
fields except the metric $\tilde{g}_{\mu\nu}$, the dilaton $\phi$ and the
Maxwell field $F_{\mu\nu}$, the action reads

$S=\frac{1}{16\pi}\int
d^{4}x\sqrt{-\tilde{g}}e^{-2\phi}\left[\tilde{R}+4\tilde{\nabla}_{\mu}\phi\tilde{\nabla}^{\mu}\phi-
F_{\mu\nu}F^{\mu\nu}\right],$ (1)

in which $\tilde{R}$ is the Ricci scalar of the metric $\tilde{g}_{\mu\nu}$.
By a conformal transformation $g_{\mu\nu}=e^{-2\phi}\tilde{g}_{\mu\nu}$, we
get another convenient representation of EMD theory in the Einstein frame, in
which the action now reads

$S=\frac{1}{16\pi}\int
d^{4}x\sqrt{-g}\left[R-2\nabla_{\mu}\phi\nabla^{\mu}\phi-e^{-b\phi}F_{\mu\nu}F^{\mu\nu}\right].$
(2)

Here $R$ is the Ricci scalar of the metric $g_{\mu\nu}$ and we have introduced
a dimensionless constant $b$ to parameterize among theories. When $b=0$, the
theory reduces to the Einstein-Maxwell theory minimally coupled to a free
scalar. When $b=1$, the theory describes the low energy limit of superstring
theory. When $b=\sqrt{3}$, the theory gives the four dimensional reduction of
Kaluza-Klein theory.

The EMD theory has been widely studied in holography due to the rich phase
structures and dynamics of charged black holes in asymptotic Anti-de Sitter
(AdS) spacetimes [4, 5, 6, 7, 8, 9, 10]. In the asymptotic flat spacetime, it
has also attracted heavy attentions since it admits hairy black hole solutions
with the scalar, vector radiative modes, in addition to the tensor channels.
It is a well-motivated theoretical laboratory to explore the impact of new
degrees of freedom in the context of testing the no-hair conjecture. The exact
static charged dilaton black hole solutions of the action (2) were found in
[11, 12]. These black holes always carry scalar hair and their charge to mass
ratio could exceed unity. The linear studies show that these solutions are
stable for generic values of the dilaton coupling and the black hole charge
[13, 14, 15, 16, 17, 18]. The existence of the dilaton breaks the
isospectrality between the axial and polar sectors of the linear
perturbations. At the nonlinear level, the dynamical evolution of an
individual black hole and the binary black hole merger were studied
numerically [19]. The binary system was also analyzed by post-Newtonian
approximation [20, 21]. The results show that these black hole systems are
difficult to distinguish from their analogs within general relativity when the
black hole charge is small, dramatic changes occur only for nearly-extremal
charged black holes on very compact orbits.

The scalarization of black hole has attracted a lot of attention due to the
discovery of spontaneous scalarization in the Einstein-scalar-Gauss-Bonnet
(EsGB) theory [22, 23, 24, 25, 26, 27, 28] and Einstein-Maxwell-scalar theory
recently [29]. Spontaneous scalarization endows black hole with scalar hair
without altering the predictions from general relativity in the weak field
limit. This mechanism was first studied in scalar-tensor theories [30, 31, 32,
33, 34, 35] and has been found in many other theories [36, 37, 38, 39, 40,
41]. Most of the studies focus on the static properties or the linear
stability of the hairy black holes [42, 43, 44, 45, 46, 47, 48, 49, 50, 51,
52, 53, 54, 55, 56, 57]. Several works focused on the nonlinear dynamics of
individual black hole scalarization [58, 59, 60] and binary black hole merger
in EsGB theory [61, 62, 63, 64, 65, 66]. The nonlinear equations of motion in
EsGB theory may not be well posed, while the Einstein-Maxwell-scalar (EMS) or
EMD models have no higher curvature corrections and allow a technical
simplification for the nonlinear studies [19, 29, 67, 68, 70, 71].

In this paper we focus on the nonlinear dynamical scalarization in spherically
symmetric spacetime in the EMD theory. When the initial dilaton is small
enough, we observe that the resulting configuration is very close to the
Reissner-Nordström (RN) black hole solution. Interestingly, we find that our
numerical method works well not only for simulations starting from initial
black holes, but also from naked singularities. We show that a naked
singularity can be protected after the scalarization and a hairy black hole
can be formed finally. But for a hairy black hole evolves from the original RN
black hole through scalarization, we find difficulty in distinguishing them
since the spectrum of the dilaton oscillation is almost independent of the
coupling parameter. This is different from the observation in the AdS
spacetimes [70].

This paper is organized as follows. In section 2, we present the equations of
motion in EMD model. In section 3, we describe our numerical method and show
the numerical results. In section 4 we summaries the results.

##  2 Equations of motion

Varying the action (2) with respect to $g_{\mu\nu},\phi$ and $A_{\mu}$, we get
the equations of motion for the metric, dilaton and gauge field, respectively.

$\displaystyle R_{\mu\nu}-\frac{1}{2}Rg_{\mu\nu}=$ $\displaystyle
2\left[\partial_{\mu}\phi\partial_{\nu}\phi-\frac{1}{2}g_{\mu\nu}\nabla_{\rho}\phi\nabla^{\rho}\phi+e^{-b\phi}\left(F_{\mu\rho}F_{\nu}^{\
\rho}-\frac{1}{4}g_{\mu\nu}F_{\rho\sigma}F^{\rho\sigma}\right)\right],$ (3)
$\displaystyle\nabla_{\mu}\nabla^{\mu}\phi=$
$\displaystyle-\frac{b}{4}e^{-b\phi}F_{\mu\nu}F^{\mu\nu},$ (4)
$\displaystyle\nabla_{\mu}\left(e^{-b\phi}F^{\mu\nu}\right)=$ $\displaystyle
0.$ (5)

The theory has a symmetry $(b,\phi)\leftrightarrow-(b,\phi)$. In the following
we consider the cases with $b<0$. To simulate the dynamic scalarization in
spherically symmetric spacetime, we adopt the Painlevé-Gullstrand (PG)-like
coordinates ansatz

$ds^{2}=-\left(1-\zeta^{2}\right)\alpha^{2}dt^{2}+2\zeta\alpha
dtdr+dr^{2}+r^{2}(d\theta^{2}+\sin^{2}\theta d\phi^{2}).$ (6)

Here $\zeta,\alpha$ are functions of $t,r$. The apparent horizon locates at
$\zeta=1$. The PG coordinates are regular on the apparent horizon, and have
been used to study the black hole formation both analytically and numerically
[72, 73, 74, 75, 60, 59]. For RN black hole, the metric functions read
$\alpha=1,\zeta=\sqrt{\frac{2M}{r}-\frac{Q^{2}}{r^{2}}}.$ Here $M$ is the
black hole mass and $Q$ the black hole charge. Note that the Arnowitt-Deser-
Misner mass in PG coordinates always evaluates to be zero, and does not
capture the correct physical mass of the spacetime [76]. Hence, we use the
Misner-Sharp mass which is defined as

$m_{MS}(t,r)=\frac{r}{2}\left(1-g^{\mu\nu}\partial_{\mu}r\partial_{\nu}r\right)=\frac{r}{2}\zeta(t,r)^{2}.$
(7)

The Misner-Sharp mass can be thought as the radially integrated energy density
of the stress energy tensor and the spacetime mass
$M=\lim_{r\to\infty}m_{MS}(t,r)$ is evaluated at the spacial infinity. The
black hole irreducible mass $M_{h}=\sqrt{\frac{A}{4\pi}}$ in which $A$ is the
area of the apparent horizon.

We take the gauge potential $A_{\mu}dx^{\mu}=A(t,r)dt$ and the dilaton field
$\phi(t,r)$. Introducing auxiliary variables

$\Phi=\partial_{r}\phi,\ \ \ P=\frac{1}{\alpha}\partial_{t}\phi-\zeta\Phi,\ \
\ E=\frac{1}{\alpha}\partial_{r}A,$ (8)

the Maxwell equations give

$E=\frac{Qe^{b\phi}}{r^{2}},$ (9)

in which $Q$ is a constant interpreted as the electric charge. The field
strength of the Maxwell field is always singular at the center. This prevents
the study of gravitational collapse starting from regular initial spherically
symmetric condition in the whole space in this model. The Einstein equations
give

$\displaystyle\alpha^{\prime}=$ $\displaystyle-\frac{rP\Phi\alpha}{\zeta},$
(10) $\displaystyle\zeta^{\prime}=$
$\displaystyle\frac{r}{2\zeta}\left(\Phi^{2}+P^{2}+\frac{Q^{2}e^{b\phi}}{r^{4}}\right)-\frac{\zeta}{2r}+rP\Phi,$
(11) $\displaystyle\partial_{t}\zeta=$
$\displaystyle\frac{r\alpha}{\zeta}\left(P+\Phi\zeta\right)\left(P\zeta+\Phi\right).$
(12)

The scalar equation becomes

$\displaystyle\partial_{t}\phi$
$\displaystyle=\alpha\left(P+\Phi\zeta\right),$ (13)
$\displaystyle\partial_{t}P$
$\displaystyle=\frac{\left(\left(P\zeta+\Phi\right)\alpha
r^{2}\right)^{\prime}}{r^{2}}-\frac{b\alpha}{2}\frac{Q^{2}e^{b\phi}}{r^{4}}.$
(14)

##  3 The numerical results

We focus on the dynamic evolution of the black hole irreducible mass $M_{h}$
and the value of the dilaton on the apparent horizon $\phi_{h}$ in this work.
We can simulate the nonlinear evolution starting from spacetimes with black
holes or naked singularities at the center. The results are reliable both from
the physics and the convergence of the numerical method.

###  3.1 The numerical setup

Let us consider the boundary conditions at first. Due to the auxillary freedom
in the metric ansatz (6), we could always fix

$\alpha|_{r\to\infty}=1,$ (15)

by rescaling the time coordinate. From (11) we see that
$\zeta\to\sqrt{\frac{2M}{r}}$ when $r\to\infty$. Here $M$ is the Misner-Sharp
mass at infinity and is a constant. We thus replace $\zeta$ by a new variable
$s=\sqrt{r}\zeta$ in the numerical simulation and set the boundary condition
for $s$ as

$s|_{r\to\infty}=\sqrt{2M}.$ (16)

We take the initial condition

$\phi=\kappa e^{-(\frac{r-6M}{M})^{2}},\ \ P=0,$ (17)

in which $\kappa$ is of order $10^{-10}$ such that the initial energy of the
dilaton can be neglected compared to the whole initial spacetime. Other types
of initial condition would not change the results qualitatively.

Given the above boundary and initial conditions, we can work out the initial
$\Phi$ from (8), $\zeta$ from (11) and $\alpha$ from (10). The initial metric
functions $\alpha,\zeta$ are very close to those of the RN black hole solution
due to the small $\phi$. Thus our simulation can be considered as a
perturbation to the RN black hole and then study its evolution in some sense,
though the RN black hole is not the exact solution of the EMD theory. This
strategy is often adopted in the dynamic simulation, such as [59, 69].

Using (12,13,14), we can derive the values of $\zeta,\phi,P$ on the next time
slice. Then from (8,10) we obtain the corresponding $\Phi,\alpha$. Iterating
this procedure, we can write the metric functions $\alpha,\zeta$ and dilaton
field $\phi,\Phi,P$ on all the following time slices. The nonlinear equation
(11) is used only once to solve the initial $\zeta$ and will not be used
again, since it is expensive to solve it.

The radial computational region ranges from $r_{0}$ to $\infty$. Here
$r_{0}\simeq 0.8M$ typically in our simulation. In fact, the initial apparent
horizon locates at $r_{h}\simeq M+\sqrt{M^{2}-Q^{2}}.$ So $r_{0}$ always lies
in the apparent horizon and the information there would not affect the outside
world in principle. We compactify the space by a coordinate transformation
$z=\frac{r}{r+M}$ and the computational region becomes $z\in(z_{0},1)$. We use
the finite difference method in the radial direction and the fourth order
Runge-Kutta method in the time direction. The radial direction is discretized
by uniform grid with $2^{11}\sim 2^{12}$ points. The Kreiss-Oliger dissipation
is employed to stabilize numerical evolution. For the first step, the equation
(11) is solved by Newton-Raphson method.

####  3.1.1 The evolution starting from a naked singularity

|
---|---
Figure 1: Left: The evolution of the metric component $-g_{tt}$ at early times
starting from a black hole (left) or a naked singularity configuration
(right). The time step between adjacent lines is $\Delta t\simeq 0.38M$.

It is known that the charge to mass ratio $Q/M$ of the hairy black hole in EMD
theory can exceed unity. To have a hairy black hole with $Q/M>1$ in the end,
we can only start the simulation from a spacetime with naked singularity,
since the equation (9) tells that the charge does not change in the evolution.
It is surprising that our numerical code works well in this case. In Fig.1, we
show the evolution of the metric component
$-g_{tt}=\left(1-\zeta^{2}\right)\alpha^{2}$ at early times. The right panel
tells that an apparent horizon forms at $t=0.76M$ (the yellow line). The
initial naked singularity is enveloped by a single horizon at first. With the
further evolution of the nonlinear perturbation, we observe some interesting
phenomena for the resulting hairy black hole in the EMD theory. The hairy
black hole can gradually develop two horizons to protect the central
singularity. However this configuration is not stable, after further
scalarization the hairy black hole grows and finally settles down to a hole
with only one horizon left surrounding the singularity. Our simulation is
reliable since the apparent horizon forms fast. There is no time for the
information on the inner cutoff $r_{0}$ to affect the region outside the
horizon. In fact one can estimate the region that can be affected. From the
metric (6) we see that the light propagates with
$\frac{dr}{dt}=\frac{1-\zeta^{2}}{2\zeta}\alpha$, which decreases with both
$t$ and $r$ near the cutoff. It turns out that the largest value of
$\frac{dr}{dt}$ is about 0.04. So the information on $r_{0}$ propagates at
most $\Delta r\simeq 0.03$ when the apparent horizon forms. This region is
surrounded by the horizon so that our simulation is safe. Furthermore, we can
check the convergence and accuracy of our numerical code, as shown in the
following subsection.

####  3.1.2 Convergence

|
---|---
Figure 2: Left: the convergence and truncation error estimate for the the
irreducible mass of the black hole $M_{h}$ (right) and the the dilaton field
on the apparent horizon $\phi_{h}$ (right). The results show that both
$\frac{M_{2N}-M_{N}}{M_{4N}-M_{2N}}$ and
$\frac{\phi_{2N}-\phi_{N}}{\phi_{4N}-\phi_{2N}}$ are approximately equal to
$2^{4}$. Here $M_{N},\phi_{N}$ stand for the corresponding results obtained by
using $N$ grid points. Note that there is no value of $M_{h},\phi_{h}$ for the
simulation with $Q=1.2$ when $t<0.35M$ since the apparent horizon has not been
formed.

For finite difference method, one often uses
$\frac{V_{2N}-V_{N}}{V_{4N}-V_{2N}}=2^{p}+O(\frac{1}{N})$ to estimate the
convergence order $p$. Here $V_{N}$ is the quantity that obtained with $N$
grid points. In Fig.2, we show the convergence of our simulation. The
evolution of both the irreducible mass $M_{h}$ and dilaton $\phi_{h}$ show
that $p\simeq 4$. This is expected since our numerical method is of the fourth
order. The accuracy of the simulation with $N=1024$ descends when $t\apprge
150M$, at where the outgoing wave have reached the far region and the uniform
gird is insufficient for high resolution. Using more dense grid can improve
the resolution. However, it is impossible and unnecessary to simulate the
system forever, since we are interested only in the phenomenon in the near
horizon region and the outgoing wave will not affect the inner region again in
the asymptotic flat spacetime. It is good enough to use $2^{11}\sim 2^{12}$
grid points in our work. Note that the accuracy is also good enough for the
evolution starting from a naked singularity.

|
---|---
Figure 3: Left: the final value $\phi_{f}$ of the dilaton on the apparent
horizon and $M_{f}$ of the black hole irreducible mass with various $b$.
Right: the distribution of $\phi$ and the Misner-Sharp mass in the radial
direction with various $b$. Here we fix $Q=0.9M$.

###  3.2 Effects of $b$ on the scalarization

In the left panel of Fig.3, we show the final value of the black hole
irreducible mass $M_{f}$ and the final value $\phi_{f}$ of the dilaton on the
apparent horizon with respect to $-b$.

The final $\phi_{f}$ increases with $-b$ at first and then decreases with
$-b$. This can be understood by two competing factors. On one hand, the
dilaton field absorbs the energy from the Maxwell field through the nonlinear
coupling term $e^{-b\phi}F_{\mu\nu}F^{\mu\nu}$ in the action. $\phi_{f}$
increases with $-b$ as the coupling becomes stronger. This argument is
supported by equation (9) in which we see that the field strength of the
Maxwell field becomes weaker for stronger coupling. On the other hand,
resembling the static Schrödinger equation, there is an effective repulsive
potential near the black hole which can be derived from the perturbation
equation for the dilaton field
$\nabla_{\mu}\nabla^{\mu}\delta\phi=V\delta\phi,$ in which
$V=\frac{b^{2}}{4}e^{-b\phi}F_{\mu\nu}F^{\mu\nu}$. For bigger $-b$, the
potential barrier becomes bigger which leads to the weaker dilaton field near
the horizon. The left panel of Fig.3 reflects the results of such competition
with the increase of the coupling parameter $-b$. In the region away from the
black hole boundary, the factor $e^{-b\phi}$ decreases rapidly with $r$ and
forms a steep effective potential away from the horizon so that the dilaton is
driven away from the black hole. In the right panel of Fig.3 we observe that
in the large $r$ region, $\phi(r)$ decreases monotonically and sharply away
from the horizon. Dilaton field with $b=-25$ has a bit higher value than that
of $b=-4$ at large $r$, but the difference is tiny.

The above argument can also shed light on understanding the distribution of
the Misner-Sharp mass. When the dilaton absorbs more energy from the Maxwell
field for stronger coupling, more energetic dilaton field is further swallowed
by the black hole, so that $M_{f}$ increases with $-b$ monotonically. For
small coupling $-b$, the dilaton can accumulate around the black hole so that
the Misner-Sharp mass increases rapidly near the black hole. For large $-b$,
the dilaton distributes more smoothly in a wider region so that the Misner-
Sharp mass increases slower with $r$.

|
---|---
Figure 4:  The evolution of the black hole irreducible mass $M_{h}$ and the
value of the dilaton on the apparent horizon $\phi_{h}$ when $Q=0.9M$. The
right panel shows the evolution of $\ln|M_{h}(t)-M_{f}|$ and
$\ln|\phi_{h}(t)-\phi_{f}|$ in which $M_{f}$ is the saturation value of
$M_{h}$ at late times.

The evolution of the value $\phi_{h}$ of the dilaton on the apparent horizon
can be divided into two stages roughly, as shown in Fig.4. At early times, the
dilaton grows abruptly and then oscillates with damping amplitude, which can
be estimated by $\phi_{h}(t)\propto\phi_{f}+e^{-\omega_{I}t}\sin\omega_{R}t$.
Here $\omega_{I}$ is the damping rate and $\omega_{R}$ is the oscillating
frequency. This behavior resembles the quasinormal mode. At late times, it
converges to $\phi_{f}$ exponentially and can be fitted as
$\phi_{h}(t)\propto\phi_{f}-e^{-\eta t}$, in which $\eta$ is a constant.
Fitting the curves shows that $\omega_{I}\simeq 0.12M$ and $\eta\simeq
0.022M$. Note that $\omega_{R}$ is sensitive to the parameter $b$, while
$\omega_{I}$ and $\eta$ are insensitive to $b$. The perturbation analysis also
reveals that $\omega_{I}$ is almost independent of $b$ in EMD theory [15, 16,
17]. These behaviors is qualitatively different from those in asymptotic AdS
spacetime in which $\omega_{I}$ is sensitive to $b$ [70]. This is because in
asymptotically flat spacetime there is no potential wall like the AdS boundary
to bounce back the perturbation and magnify the difference caused by the
coupling $b$.

The evolution of black hole irreducible mass $M_{h}$ can also be divided into
two stages. At early time, it increases exponentially in a step like form
which coincides with the pulse of the dilaton growth. The growth of the
irreducible mass can be fitted by $M_{h}(t)\propto M_{f}-e^{-\gamma_{i}t}$. At
late times, it saturates smoothly to $M_{f}$ with $M_{h}(t)\propto
M_{f}-e^{-\gamma_{f}t}$. Here both $\gamma_{i,f}$ are constants that are
insensitive to $b$. Furthermore, there are relations

$\gamma_{i}=2\omega_{I},\ \ \ \gamma_{f}=2\eta.$ (18)

---
Figure 5: The time needed to form the apparent horizon starting from a naked
singularity when $Q=1.2M$ with respect to $b$. When $-b\lesssim 3.7$, the
numerical codes does not work well.

This can be seen explicitly from the inset of Fig.4. Though these relations
are obtained from fully nonlinear evolution, it can be understood from the
viewpoint of the linear perturbation analysis. The perturbation of the scalar
field invokes the back-reaction of the metric only at the second order [78],
or from the Einstein equation (3), there is $\delta\zeta\sim\delta\phi^{2}$.
The black hole irreducible mass $M_{h}$ is directly related to the radius of
the apparent horizon $r_{h}$, which is the zero point of $\zeta-1=0$. So we
can expect that $\delta r_{h}\sim\delta\phi^{2}$ during the evolution which
implies (18). This argument is independent of the special theory and should be
hold in general. Indeed, these relations were also observed in asymptotic AdS
spacetime [70, 71].

We also studied the effect of $b$ on the dynamic scalarization starting from
spacetime with naked singularity. For large $-b$, the apparent horizon forms
quickly so that the naked singularity can be protected immediately. The
following evolution is qualitatively similar to that starting from spacetime
with charged black holes. For small $-b$, the coupling between the Maxwell
field and the dilaton is weak, it needs longer period of time to form the
apparent horizon, as shown in Fig.5. It should be noted that $-b$ cannot be
too small, otherwise the spacetime geometry changes too slowly to form the
apparent horizon. The information from the inner cutoff would affect the
outside world and this makes our numerical simulation crash.

##  4 Summary and discussion

We have studied the fully nonlinear dynamic scalarization in spherically
symmetric spacetime in EMD theory, starting from initial charged black holes
or naked singularities at the center. The energy of the initial dilaton is
very small compared to those of the gravitational and Maxwell field. Due to
the coupling between the dilaton and the Maxwell field, part of the energy is
transferred from the Maxwell field to the dilaton, resulting in a nontrivial
dilaton field at the end. The black holes absorb some of the energy from the
dilaton, such that their irreducible mass also increases. On the other hand,
the coupling between the dilaton and the Maxwell field provides an effective
repulsive force to the dilaton and drives the dilaton away from the black
hole. When the initial configuration is a naked singularity, apparent horizon
can be developed immediately to protect the singularity. The naked singularity
evolves into a hairy hole with one horizon, further scalarization can
accommodate a black hole with two horizons. Finally the hairy black hole
evolves into a stable configuration with only one horizon left. The following
evolutions are qualitatively the same as those starting from spacetimes with
ordinary black holes.

At early times, the dilaton grows abruptly. Then it oscillates with damping
amplitude. At late times, it converges to a nonvanishing final value
exponentially. Both the damping rate and the converging rate are almost
independent of the coupling parameter and the black hole charge to mass ratio.
The irreducible mass of the black hole also grows exponentially at early times
and then saturates to the final value at late times. The saturating rates are
just twice the damping rate and the convergence rate of the dilaton. Though
these behaviors are found from fully nonlinear evolutions, they can be
understood intuitively from the viewpoint of the linear perturbation. Our
study showed that it is hard to distinguish the charged dilaton black hole
from the RN black hole by studying the wave spectrum.

This paper focused on the configuration with negligibly small initial matter
perturbation. As an extension, one can generalize the discussion with large
initial perturbation. The initial black hole at the center could be taken as a
small seed black hole. The model can then be considered as a dynamical process
of gravitational collapse of the matter. Since the PG coordinate are horizon
penetrating, we can evolve the system until the black hole becomes stable and
confirm whether a hairy black hole forms or not in the end. This could
compensate the disadvantages of the method used in [58]. These works are in
progress.

## Acknowledgments

Peng Liu would like to thank Yun-Ha Zha for her kind encouragement during this
work. The authors thank Yu Tian, Qian Chen and Peng-Cheng Li for the helpful
discussions. This research is supported by National Key R&D Program of China
under Grant No.2020YFC2201400, and the Natural Science Foundation of China
under Grant Nos.11690021, 11947067, 12005077, 11847055, 11905083, 11805083 and
Guangdong Basic and Applied Basic Research Foundation under Grant No.
2021A1515012374.

## References

  * [1] T. Kaluza, On the Problem of Unity in Physics, Sitzungsber. Preuss. Akad. Wiss. Berlin (Math. Phys.) 1921 (1921) 966972\.
  * [2] E. Cremmer and B. Julia, “The N=8 Supergravity Theory. 1. The Lagrangian,” Phys. Lett. 80B, 48 (1978).
  * [3] J. Polchinski, String Theory, Volume 2: Superstring Theory and Beyond, Cambridge, 1998.
  * [4] O. DeWolfe, S. S. Gubser and C. Rosen, “Dynamic critical phenomena at a holographic critical point,” Phys. Rev. D 84 (2011), 126014 [arXiv:1108.2029 [hep-th]].
  * [5] O. DeWolfe, S. S. Gubser and C. Rosen, “A holographic critical point,” Phys. Rev. D 83 (2011), 086005 [arXiv:1012.1864 [hep-th]].
  * [6] Y. Liu and Y. W. Sun, “Holographic Superconductors from Einstein-Maxwell-Dilaton Gravity,” JHEP 07 (2010), 099 [arXiv:1006.2726 [hep-th]].
  * [7] D. Giataganas, U. Gürsoy and J. F. Pedraza, “Strongly-coupled anisotropic gauge theories and holography,” Phys. Rev. Lett. 121 (2018) no.12, 121601 [arXiv:1708.05691 [hep-th]].
  * [8] J. Knaute, R. Yaresko and B. Kämpfer, “Holographic QCD phase diagram with critical point from Einstein–Maxwell-dilaton dynamics,” Phys. Lett. B 778 (2018), 419-425 [arXiv:1702.06731 [hep-ph]].
  * [9] A. Ballon-Bayona, H. Boschi-Filho, E. F. Capossoli and D. M. Rodrigues, “Criticality from Einstein-Maxwell-dilaton holography at finite temperature and density,” Phys. Rev. D 102 (2020) no.12, 126003 [arXiv:2006.08810 [hep-th]].
  * [10] J. X. Mo and S. Q. Lan, “Dynamic phase transition of charged dilaton black holes,” [arXiv:2105.00868 [gr-qc]].
  * [11] G. W. Gibbons and K. i. Maeda, “Black Holes and Membranes in Higher Dimensional Theories with Dilaton Fields,” Nucl. Phys. B 298 (1988), 741-775
  * [12] D. Garfinkle, G. T. Horowitz and A. Strominger, “Charged black holes in string theory,” Phys. Rev. D 43 (1991), 3140 [erratum: Phys. Rev. D 45 (1992), 3888]
  * [13] S. Fernando and K. Arnold, “Scalar perturbations of charged dilaton black holes,” Gen. Rel. Grav. 36 (2004), 1805-1819 [arXiv:hep-th/0312041 [hep-th]].
  * [14] V. Ferrari, M. Pauri and F. Piazza, “Quasinormal modes of charged, dilaton black holes,” Phys. Rev. D 63 (2001), 064009 [arXiv:gr-qc/0005125 [gr-qc]].
  * [15] R. A. Konoplya, “Quasinormal modes of the electrically charged dilaton black hole,” Gen. Rel. Grav. 34 (2002), 329-335 [arXiv:gr-qc/0109096 [gr-qc]].
  * [16] C. Y. Zhang, S. J. Zhang and B. Wang, “Charged scalar perturbations around Garfinkle–Horowitz–Strominger black holes,” Nucl. Phys. B 899 (2015), 37-54 [arXiv:1501.03260 [hep-th]].
  * [17] R. Brito and C. Pacilio, “Quasinormal modes of weakly charged Einstein-Maxwell-dilaton black holes,” Phys. Rev. D 98 (2018) no.10, 104042 [arXiv:1807.09081 [gr-qc]].
  * [18] J. L. Blázquez-Salcedo, S. Kahlen and J. Kunz, “Quasinormal modes of dilatonic Reissner–Nordström black holes,” Eur. Phys. J. C 79 (2019) no.12, 1021 [arXiv:1911.01943 [gr-qc]].
  * [19] E. W. Hirschmann, L. Lehner, S. L. Liebling and C. Palenzuela, “Black Hole Dynamics in Einstein-Maxwell-Dilaton Theory,” Phys. Rev. D 97 (2018) no.6, 064032 [arXiv:1706.09875 [gr-qc]].
  * [20] F. L. Julié, “On the motion of hairy black holes in Einstein-Maxwell-dilaton theories,” JCAP 01 (2018), 026 [arXiv:1711.10769 [gr-qc]].
  * [21] M. Khalil, N. Sennett, J. Steinhoff, J. Vines and A. Buonanno, “Hairy binary black holes in Einstein-Maxwell-dilaton theory and their effective-one-body description,” Phys. Rev. D 98 (2018) no.10, 104010 [arXiv:1809.03109 [gr-qc]].
  * [22] D. D. Doneva and S. S. Yazadjiev, New Gauss-Bonnet Black Holes with Curvature-Induced Scalarization in Extended Scalar-Tensor Theories, Phys. Rev. Lett. 120, no.13, 131103 (2018) [arXiv:1711.01187 [gr-qc]].
  * [23] 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, Phys. Rev. Lett. 120, no.13, 131104 (2018) [arXiv:1711.02080 [gr-qc]].
  * [24] G. Antoniou, A. Bakopoulos and P. Kanti, Evasion of No-Hair Theorems and Novel Black-Hole Solutions in Gauss-Bonnet Theories, Phys. Rev. Lett. 120, no.13, 131102 (2018) [arXiv:1711.03390 [hep-th]].
  * [25] P. V. Cunha, C. A. Herdeiro and E. Radu, Spontaneously Scalarized Kerr Black Holes in Extended Scalar-Tensor-Gauss-Bonnet Gravity, Phys. Rev. Lett. 123, no.1, 011101 (2019) [arXiv:1904.09997 [gr-qc]].
  * [26] A. Dima, E. Barausse, N. Franchini and T. P. Sotiriou, “Spin-induced black hole spontaneous scalarization,” Phys. Rev. Lett. 125 (2020) no.23, 231101. [arXiv:2006.03095 [gr-qc]].
  * [27] C. A. R. Herdeiro, E. Radu, H. O. Silva, T. P. Sotiriou and N. Yunes, Spin-induced scalarized black holes, Phys.Rev.Lett. 126 (2021) 1, 011103. [arXiv:2009.03904 [gr-qc]].
  * [28] E. Berti, L. G. Collodel, B. Kleihaus and J. Kunz, Spin-induced black-hole scalarization in Einsteinscalar-Gauss-Bonnet theory, Phys.Rev.Lett. 126 (2021) 1, 011104. [arXiv:2009.03905 [gr-qc]].
  * [29] C. A. R. Herdeiro, E. Radu, N. Sanchis-Gual and J. A. Font, “Spontaneous Scalarization of Charged Black Holes,” Phys. Rev. Lett. 121, no. 10, 101102 (2018). [arXiv:1806.05190 [gr-qc]].
  * [30] T. Damour and G. Esposito-Farese, Nonperturbative strong field effects in tensor-scalar theories of gravitation, Phys. Rev. Lett., vol. 70, pp. 2220-2223, 1993.
  * [31] T. Damour and G. Esposito-Farese, “Tensor - scalar gravity and binary pulsar experiments,” Phys. Rev. D 54 (1996), 1474-1491 [arXiv:gr-qc/9602056 [gr-qc]].
  * [32] T. Harada, “Stability analysis of spherically symmetric star in scalar-tensor theories of gravity,” Prog. Theor. Phys. 98 (1997), 359-379 [arXiv:gr-qc/9706014 [gr-qc]].
  * [33] V. Cardoso, I. P. Carucci, P. Pani and T. P. Sotiriou, “Black holes with surrounding matter in scalar-tensor theories,” Phys. Rev. Lett. 111 (2013), 111101 [arXiv:1308.6587 [gr-qc]].
  * [34] V. Cardoso, I. P. Carucci, P. Pani and T. P. Sotiriou, “Matter around Kerr black holes in scalar-tensor theories: scalarization and superradiant instability,” Phys. Rev. D 88 (2013), 044056 [arXiv:1305.6936 [gr-qc]].
  * [35] C. Y. Zhang, S. J. Zhang and B. Wang, “Superradiant instability of Kerr-de Sitter black holes in scalar-tensor theory,” JHEP 08 (2014), 011 [arXiv:1405.3811 [hep-th]].
  * [36] C. A. R. Herdeiro and E. Radu, “Black hole scalarization from the breakdown of scale invariance,” Phys. Rev. D 99 (2019) no.8, 084039 [arXiv:1901.02953 [gr-qc]].
  * [37] Y. Brihaye, C. Herdeiro and E. Radu, “The scalarised Schwarzschild-NUT spacetime,” Phys. Lett. B 788, 295-301 (2019) [arXiv:1810.09560 [gr-qc]].
  * [38] J. M. S. Oliveira and A. M. Pombo, “Spontaneous vectorization of electrically charged black holes,” Phys. Rev. D 103 (2021) no.4, 044004 [arXiv:2012.07869 [gr-qc]].
  * [39] Y. Brihaye and Y. Verbin, “Scalarized dyonic black holes in vector-tensor Horndeski gravity,” Phys. Rev. D 104 (2021) no.2, 024047 [arXiv:2105.11402 [gr-qc]].
  * [40] G. Antoniou, A. Lehébel, G. Ventagli and T. P. Sotiriou, “Black hole scalarization with Gauss-Bonnet and Ricci scalar couplings,” Phys. Rev. D 104 (2021) no.4, 044002 [arXiv:2105.04479 [gr-qc]].
  * [41] Y. Brihaye and B. Hartmann, “Spontaneous scalarization of boson stars,” JHEP 09 (2019), 049. [arXiv:1903.10471 [gr-qc]].
  * [42] C. F. B. Macedo, J. Sakstein, E. Berti, L. Gualtieri, H. O. Silva and T. P. Sotiriou, “Self-interactions and Spontaneous Black Hole Scalarization,” Phys. Rev. D 99 (2019) no.10, 104041 [arXiv:1903.06784 [gr-qc]].
  * [43] Y. Brihaye, C. Herdeiro and E. Radu, “Black Hole Spontaneous Scalarisation with a Positive Cosmological Constant,” Phys. Lett. B 802, 135269 (2020). [arXiv:1910.05286 [gr-qc]].
  * [44] H. Guo, S. Kiorpelidi, X. M. Kuang, E. Papantonopoulos, B. Wang and J. P. Wu, “Spontaneous holographic scalarization of black holes in Einstein-scalar-Gauss-Bonnet theories,” Phys. Rev. D 102 (2020) no.8, 084029. [arXiv:2006.10659 [hep-th]].
  * [45] D. Astefanesei, C. Herdeiro, A. Pombo and E. Radu, “Einstein-Maxwell-scalar black holes: classes of solutions, dyons and extremality,” JHEP 10 (2019), 078. [arXiv:1905.08304 [hep-th]].
  * [46] K. Lin, S. Zhang, C. Zhang, X. Zhao, B. Wang and A. Wang, “No static regular black holes in Einstein-complex-scalar-Gauss-Bonnet gravity,” Phys. Rev. D 102 (2020) no.2, 024034. [arXiv:2004.04773 [gr-qc]].
  * [47] H. S. Liu, H. Lu, Z. Y. Tang and B. Wang, “Black hole scalarization in Gauss-Bonnet extended Starobinsky gravity,” Phys. Rev. D 103 (2021) no.8, 084043 [arXiv:2004.14395 [gr-qc]].
  * [48] H. Guo, X. M. Kuang, E. Papantonopoulos and B. Wang, “Horizon curvature and spacetime structure influences on black hole scalarization,” Eur. Phys. J. C 81 (2021) no.9, 842 [arXiv:2012.11844 [gr-qc]].
  * [49] J. L. Blázquez-Salcedo, D. D. Doneva, J. Kunz and S. S. Yazadjiev, “Radial perturbations of the scalarized Einstein-Gauss-Bonnet black holes,” Phys. Rev. D 98 (2018) no.8, 084011 [arXiv:1805.05755 [gr-qc]].
  * [50] J. L. Blázquez-Salcedo, D. D. Doneva, S. Kahlen, J. Kunz, P. Nedkova and S. S. Yazadjiev, “Axial perturbations of the scalarized Einstein-Gauss-Bonnet black holes,” Phys. Rev. D 101 (2020) no.10, 104006 [arXiv:2003.02862 [gr-qc]].
  * [51] J. L. Blázquez-Salcedo, D. D. Doneva, S. Kahlen, J. Kunz, P. Nedkova and S. S. Yazadjiev, “Polar quasinormal modes of the scalarized Einstein-Gauss-Bonnet black holes,” Phys. Rev. D 102 (2020) no.2, 024086 [arXiv:2006.06006 [gr-qc]].
  * [52] H. O. Silva, C. F. B. Macedo, T. P. Sotiriou, L. Gualtieri, J. Sakstein and E. Berti, “Stability of scalarized black hole solutions in scalar-Gauss-Bonnet gravity,” Phys. Rev. D 99 (2019) no.6, 064011 [arXiv:1812.05590 [gr-qc]].
  * [53] S. J. Zhang, B. Wang, A. Wang and J. F. Saavedra, “Object picture of scalar field perturbation on Kerr black hole in scalar-Einstein-Gauss-Bonnet theory,” Phys. Rev. D 102 (2020) no.12, 124056 [arXiv:2010.05092 [gr-qc]].
  * [54] D. C. Zou and Y. S. Myung, “Radial perturbations of the scalarized black holes in Einstein-Maxwell-conformally coupled scalar theory,” Phys. Rev. D 102 (2020) no.6, 064011 [arXiv:2005.06677 [gr-qc]].
  * [55] Y. S. Myung and D. C. Zou, “Stability of scalarized charged black holes in the Einstein–Maxwell–Scalar theory,” Eur. Phys. J. C 79 (2019) no.8, 641 [arXiv:1904.09864 [gr-qc]].
  * [56] Y. S. Myung and D. C. Zou, “Instability of Reissner–Nordström black hole in Einstein-Maxwell-scalar theory,” Eur. Phys. J. C 79 (2019) no.3, 273 [arXiv:1808.02609 [gr-qc]].
  * [57] D. D. Doneva and S. S. Yazadjiev, “Dynamics of the nonrotating and rotating black hole scalarization,” Phys. Rev. D 103 (2021) no.6, 064024 [arXiv:2101.03514 [gr-qc]].
  * [58] J. L. Ripley and F. Pretorius, “Gravitational collapse in Einstein dilaton-Gauss–Bonnet gravity,” Class. Quant. Grav. 36 (2019) no.13, 134001. [arXiv:1903.07543 [gr-qc]].
  * [59] J. L. Ripley and F. Pretorius, “Scalarized Black Hole dynamics in Einstein dilaton Gauss-Bonnet Gravity,” Phys. Rev. D 101 (2020) no.4, 044015. [arXiv:1911.11027 [gr-qc]].
  * [60] J. L. Ripley and F. Pretorius, “Dynamics of a $\mathbb{Z}_{2}$ symmetric EdGB gravity in spherical symmetry,” Class. Quant. Grav. 37 (2020) no.15, 155003. [arXiv:2005.05417 [gr-qc]].
  * [61] H. Witek, L. Gualtieri, P. Pani and T. P. Sotiriou, “Black holes and binary mergers in scalar Gauss-Bonnet gravity: scalar field dynamics,” Phys. Rev. D 99 (2019) no.6, 064035 [arXiv:1810.05177 [gr-qc]].
  * [62] F. L. Julié and E. Berti, “Post-Newtonian dynamics and black hole thermodynamics in Einstein-scalar-Gauss-Bonnet gravity,” Phys. Rev. D 100 (2019) no.10, 104061 [arXiv:1909.05258 [gr-qc]].
  * [63] W. E. East and J. L. Ripley, “Evolution of Einstein-scalar-Gauss-Bonnet gravity using a modified harmonic formulation,” Phys. Rev. D 103 (2021) no.4, 044040 [arXiv:2011.03547 [gr-qc]].
  * [64] W. E. East and J. L. Ripley, “Dynamics of Spontaneous Black Hole Scalarization and Mergers in Einstein-Scalar-Gauss-Bonnet Gravity,” Phys. Rev. Lett. 127 (2021) no.10, 101102 [arXiv:2105.08571 [gr-qc]].
  * [65] H. O. Silva, H. Witek, M. Elley and N. Yunes, “Dynamical Descalarization in Binary Black Hole Mergers,” Phys. Rev. Lett. 127 (2021) no.3, 031101 [arXiv:2012.10436 [gr-qc]].
  * [66] H. J. Kuan, D. D. Doneva and S. S. Yazadjiev, “Dynamical Formation of Scalarized Black Holes and Neutron Stars through Stellar Core Collapse,” Phys. Rev. Lett. 127 (2021) no.16, 161103 [arXiv:2103.11999 [gr-qc]].
  * [67] P. G. S. Fernandes, C. A. R. Herdeiro, A. M. Pombo, E. Radu and N. Sanchis-Gual, “Spontaneous Scalarisation of Charged Black Holes: Coupling Dependence and Dynamical Features,” Class. Quant. Grav. 36 (2019) no.13, 134002 [erratum: Class. Quant. Grav. 37 (2020) no.4, 049501]. [arXiv:1902.05079 [gr-qc]].
  * [68] P. G. S. Fernandes, C. A. R. Herdeiro, A. M. Pombo, E. Radu and N. Sanchis-Gual, “Charged black holes with axionic-type couplings: Classes of solutions and dynamical scalarization,” Phys. Rev. D 100 (2019) no.8, 084045. [arXiv:1908.00037 [gr-qc]].
  * [69] D. Astefanesei, J. L. Blázquez-Salcedo, C. Herdeiro, E. Radu and N. Sanchis-Gual, “Dynamically and thermodynamically stable black holes in Einstein-Maxwell-dilaton gravity,” JHEP 07 (2020), 063 [arXiv:1912.02192 [gr-qc]].
  * [70] C. Y. Zhang, P. Liu, Y. Liu, C. Niu and B. Wang, “Evolution of Anti-de Sitter black holes in Einstein-Maxwell-dilaton theory,” [arXiv:2104.07281 [gr-qc]].
  * [71] C. Y. Zhang, P. Liu, Y. Liu, C. Niu and B. Wang, “Dynamical charged black hole spontaneous scalarization in anti–de Sitter spacetimes,” Phys. Rev. D 104 (2021) no.8, 084089 [arXiv:2103.13599 [gr-qc]].
  * [72] R. J. Adler, J. D. Bjorken, P. Chen and J. S. Liu, “Simple analytic models of gravitational collapse,” Am. J. Phys. 73 (2005), 1148-1159 [arXiv:gr-qc/0502040 [gr-qc]].
  * [73] J. Ziprick and G. Kunstatter, “Numerical study of black-hole formation in Painleve-Gullstrand coordinates,” Phys. Rev. D 79 (2009), 101503 [arXiv:0812.0993 [gr-qc]].
  * [74] Y. Kanai, M. Siino and A. Hosoya, “Gravitational collapse in Painleve-Gullstrand coordinates,” Prog. Theor. Phys. 125 (2011), 1053-1065 [arXiv:1008.0470 [gr-qc]].
  * [75] J. L. Ripley, “Excision and avoiding the use of boundary conditions in numerical relativity,” Class. Quant. Grav. 36 (2019) no.23, 237001 [arXiv:1908.04234 [gr-qc]].
  * [76] M. Shibata, Numerical Relativity (100 Years of General Relativity) - World Scientific Publishing (2015).
  * [77] C. W. Misner and D. H. Sharp, Relativistic Equations for Adiabatic, Spherically Symmetric Gravitational Collapse, Phys. Rev. 136, B571 (1964).
  * [78] R. Brito, V. Cardoso and P. Pani, “Superradiance: New Frontiers in Black Hole Physics,” Lect. Notes Phys. 906 (2015), pp.1-237 [arXiv:1501.06570 [gr-qc]].