# Black hole perturbations in spatially covariant gravity with just two
tensorial degrees of freedom

Jin Saito<EMAIL_ADDRESS>Department of Physics, Rikkyo University,
Toshima, Tokyo 171-8501, Japan  Tsutomu Kobayashi<EMAIL_ADDRESS>Department of Physics, Rikkyo University, Toshima, Tokyo 171-8501, Japan

###### Abstract

We study linear perturbations around a static and spherically symmetric black
hole solution in spatially covariant gravity with just two tensorial degrees
of freedom. In this theory, gravity modification is characterized by a single
time-dependent function that appears in the coefficient of $K^{2}$ in the
action, where $K$ is the trace of the extrinsic curvature. The background
black hole solution is given by the Schwarzschild solution foliated by the
maximal slices and has a universal horizon at which the lapse function
vanishes. We show that the quadratic action for the odd-parity perturbations
is identical to that in general relativity upon performing an appropriate
coordinate transformation. This in particular implies that the odd-parity
perturbations propagate at the speed of light, with the inner boundary being
the usual event horizon. We also derive the quadratic action for even-parity
perturbations. In the even-parity sector, one of the two tensorial degrees of
freedom is mixed with an instantaneous scalar mode, rendering the system
distinct from that in general relativity. We find that monopole and dipole
perturbations, which are composed solely of the instantaneous scalar mode,
have no solutions regular both at the universal horizon and infinity (except
for the trivial one corresponding to the constant shift of the mass
parameter). We also consider stationary perturbations with higher multipoles.
By carefully treating the locations of the inner boundary, we show that also
in this case there are no solutions regular both at the inner boundary and
infinity. Thus, the black hole solution we consider is shown to be
perturbatively unique.

††preprint: RUP-23-14

##  I Introduction

Lovelock’s theorem [1, 2] states that a diffeomorphism-invariant theory
constructed only from the metric tensor in four dimensions leads uniquely to
the Einstein equations in general relativity (GR). To modify gravity, one must
therefore break at least one of the postulates of the theorem, A typical way
of modifying gravity is to add a new dynamical degree of freedom (DOF). The
simplest example is scalar-tensor gravity having one scalar and two tensorial
DOFs, which has been studied extensively, with a particular focus over the
past decade on the Horndeski theory [3, 4, 5] and its extensions called
degenerate higher-order scalar-tensor (DHOST) theories [6, 7, 8, 9, 10, 11]
(see Ref. [12] for a review). An apparently different approach is abandoning
full diffeomorphism invariance. For example, one can consider spatially
covariant theories respecting only spatial diffeomorphism invariance under the
transformations of spatial coordinates, $x^{i}\to\tilde{x}^{i}(t,x^{j})$ [13,
14]. However, this approach is basically equivalent to adding new dynamical
degrees of freedom and spatially covariant gravity is regarded as a gauge-
fixed version of a fully covariant scalar-tensor theory. ; Indeed, starting
from a spatially covariant theory, one can introduce a Stückelberg scalar
field to restore full diffeomorphism invariance and write an action for a
corresponding fully covariant scalar-tensor theory [15].

In scalar-tensor theories, the dynamical scalar DOF (say $\phi$) obeys a wave-
like equation with some propagation speed $c_{s}$, the concrete form of which
depends on the Lagrangian. An interesting twist is a case where $\phi$ is an
instantaneous mode, $c_{s}=\infty$, which occurs when the form of the
Lagrangian is chosen appropriately. In such theories, only tensorial DOFs
propagate, while the configuration of the instantaneous scalar DOF is
determined completely from boundary conditions. Within the so-called
$P(\phi,X)$ theory (where $X=-(\partial\phi)^{2}/2$), the cuscuton theory,
$P=\mu(\phi)\sqrt{X}-V(\phi)$, gives rise to such a nonpropagating scalar
field [16], leading to interesting cosmology [17]. The cuscuton theory was
extended in Ref. [18] by demanding that $c_{s}=\infty$ in the Horndeski and
slightly more general theories. Performing a more rigorous Hamiltonian
analysis, one can determine the subset of spatially covariant theories of
gravity having two tensorial DOFs only. This was done in the case where the
action is linear in the lapse function in Ref. [19] and then the general
conditions to render the scalar mode nondynamical were derived in Ref. [20].
Spatially covariant gravity satisfying these conditions is dubbed as a TTDOF
theory. Since the dynamical DOFs in TTDOF theories are the same as those in
GR, TTDOF gravity may be thought of as minimal modification of GR. Determining
the most general form of the action for TTDOF theories has turned out to be
very difficult, and the authors of Ref. [20] managed to obtain a particular
example of a family of TTDOF theories assuming that the action is quadratic in
the extrinsic curvature and linear in the intrinsic curvature as in the
Arnowitt-Deser-Misner (ADM) expression for the Einstein-Hilbert term. The
resultant action contains that for the extended cuscuton theory as a special
case.

A family of TTDOF theories obtained in Ref. [20] has several time-dependent
parameters, which can in principle be freely chosen from a purely theoretical
point of view. Some of them can however be fixed by requiring that the speed
of gravitational waves is equal to that of light and the usual GR behavior is
restored in the weak-gravity regime in the Solar System [21]. The TTDOF theory
with the thus reduced parameter space can still mimic the background evolution
of the standard $\Lambda$CDM model [22]. Interestingly, Ref. [22] studied the
CMB constraints on the TTDOF cosmological model and reported a $\sim$
4$\sigma$ deviation from the $\Lambda$CDM model based on GR.

Ignoring the terms that are relevant only on cosmological scales, we have a
TTDOF theory characterized by a single time-dependent parameter, the impacts
of which are not seen in the propagation of gravitational waves and in the
weak-gravity regime in the Solar System. The theory admits the Schwarzschild
solution foliated by the maximal slices [21]. This can easily be seen by
noting that the modified gravity parameter enters only in the coefficient of
the $K^{2}$ term in the TTDOF action we are considering, where $K$ is the
trace of the extrinsic curvature. The situation here is quite similar to that
of Einstein-Aether theory and the low-energy limit of Horava gravity [23, 24].
Such a black hole solution in modified gravity is sometimes called stealth, as
the geometry is described by the Schwarzschild solution in GR even though the
Stückelberg scalar exhibits a nontrivial configuration. In this paper, we
study linear perturbations around this stealth black hole solution in the
TTDOF theory, paying particular attention to the behavior of the instantaneous
scalar mode whose inner boundary conditions are imposed at the universal
horizon (the location at which the lapse function vanishes) rather than the
usual event horizon. One of our goals is to see whether the stealth
Schwarzschild solution is unique or not. The question was partially addressed
in Ref. [21], but there only static monopole perturbations were considered. In
the present paper, we start with deriving quadratic actions for all (time-
dependent) perturbations with higher multipoles, following closely the
formulation of black hole perturbation theory in Horndeski [25, 26] and DHOST
gravity [27, 28, 29] (see also Refs. [30, 31, 32, 33, 34]). This result itself
has wider applications than merely investigating perturbative uniqueness of
the background black hole solution. See also Ref. [35] for a related analysis
of the behavior of the instantaneous scalar field around a black hole with a
universal horizon.

This paper is organized as follows. In the next section, we briefly review the
TTDOF theory and its static and spherically symmetric black hole solution. In
Sec. III, we study the odd-parity perturbations. Then, we consider the even-
parity perturbations, presenting the main results separately for monopole,
dipole, and higher-multipole modes in Sec. IV. Finally, we draw our
conclusions in Sec. V.

##  II A black hole solution in TTDOF theory

In this section, we introduce the TTDOF theory [20] and its static and
spherically symmetric black hole solution [21].

###  II.1 TTDOF theory

We consider modified gravity with just two tensorial degrees of freedom
(TTDOF) developed in Ref. [20]. The symmetry of the theory is the spatial
diffeomorphism invariance under $x^{i}\to\tilde{x}^{i}(t,x^{j})$, and
therefore we use the ADM variables to write the action. Specifically, the
action that we consider in this paper is given by

$\displaystyle S=\frac{1}{2}\int{\rm d}t{\rm
d}^{3}x\sqrt{\gamma}N\left[K_{ij}K^{ij}-\frac{1}{3}\left(\frac{2N}{\beta+N}+1\right)K^{2}+R\right],$
(1)

where $N$ is the lapse function, $R$ is the three-dimensional Ricci scalar
calculated from the spatial metric $\gamma_{ij}$, and $K_{ij}$ is the
extrinsic curvature of constant time hypersurfaces,

$\displaystyle
K_{ij}:=\frac{1}{2N}\left(\partial_{t}\gamma_{ij}-D_{i}N_{j}-D_{j}N_{i}\right),$
(2)

with $N_{i}$ being the shift vector and $D_{i}$ the covariant derivative
operator associated with $\gamma_{ij}$. Here, $\beta=\beta(t)$ is an arbitrary
function of time characterizing a modification of GR, and by setting $\beta=0$
the action (1) reduces to the ADM expression of the Einstein-Hilbert action.
Among a large family of TTDOF theories introduced in Ref. [20], we focus on
its particular subset described by the action (1), because standard gravity is
reproduced in the solar system (in the sense that the parameterized post-
Newtonian parameter $\gamma=1$) and gravitational waves propagate at the speed
of light [21], while rendering the cosmology nontrivial in this theory [22].

One can recover the full four-dimensional diffeomorphism invariance by
introducing the Stückelberg scalar field. The resultant covariant expression
of the action (1) belongs to the so-called U-degenerate theory [21], which is
a higher-order scalar-tensor theory satisfying the degeneracy conditions only
when one takes the unitary gauge [36]. However, the apparent scalar degree of
freedom is in fact an instantaneous mode having infinite propagation speed and
obeying an elliptic equation. Therefore, the scalar field does not propagate
and its behavior is determined completely by boundary conditions.

###  II.2 A black hole solution

The ADM variables for a static and spherically symmetric solution are of the
form

$\displaystyle N=N(r),\quad N_{i}{\rm d}x^{i}=B(r)F(r){\rm
d}r,\quad\gamma_{ij}{\rm d}x^{i}{\rm d}x^{j}=F^{2}{\rm
d}r^{2}+r^{2}\sigma_{ab}{\rm d}x^{a}{\rm d}x^{b},$ (3)

with $\sigma_{ab}{\rm d}x^{a}{\rm d}x^{b}={\rm d}\theta^{2}+\sin^{2}\theta{\rm
d}\varphi^{2}$ ($a,b=\theta,\varphi$). Note at this point that we cannot set
$B=0$ in general when working in the ADM action (i.e., in the unitary gauge),
because we no longer have the freedom to perform a coordinate transformation
$t\to\tilde{t}(t,r)$.

Substituting Eq. (3) to the action (1) and varying it with respect to $N$,
$B$, and $F$, we obtain a set of the field equations for these variables. A
general solution has not been derived, but one can at least find the following
particular solution:

$\displaystyle N=N_{0}\sqrt{f(r)},\quad F(r)=\frac{1}{\sqrt{f(r)}},\quad
B(r)=\frac{N_{0}b_{0}}{r^{2}},$ (4)

with

$\displaystyle f(r):=1-\frac{\mu_{0}}{r}+\frac{b_{0}^{2}}{r^{4}},$ (5)

where $N_{0}$, $\mu_{0}$, and $b_{0}$ are integration constants. The solution
represents a foliation of the Schwarzschild geometry by maximal slices
($K=0$). Indeed, by performing the coordinate transformation

$\displaystyle\tau=N_{0}t-\int^{r}\frac{b_{0}/r^{2}}{\sqrt{f}(1-\mu_{0}/r)}{\rm
d}r,$ (6)

and introducing the Stückelberg field $\phi=t(\tau,r)$ accordingly, one
obtains ${\rm d}s^{2}=-(1-\mu_{0}/r){\rm d}\tau^{2}+(1-\mu_{0}/r)^{-1}{\rm
d}r^{2}+r^{2}\sigma_{ab}{\rm d}x^{a}{\rm d}x^{b}$, which is nothing but the
standard Schwarzschild metric. From this observation, we see that $\mu_{0}$ is
the mass parameter, $b_{0}$ characterizes the foliation, and $N_{0}$ just
corresponds to the rescaling of the time coordinate. It should be emphasized
that the TTDOF theory with the action (1) admits the above solution even if
$\beta$ is an arbitrary time-dependent function.

For $b_{0}\leq b_{0,c}:=3\sqrt{3}\mu_{0}^{2}/16$, $N(r)$ vanishes at some
$r\,(>0)$. The location at which $N(r)=0$ is called the universal horizon,
which is the causal boundary for the scalar mode with infinite propagation
speed. Only for $b_{0}=b_{0,c}$ the universal horizon is regular, while $f(r)$
changes its sign at the universal horizon for $b_{0}<b_{0,c}$. In the rest of
the paper, we basically consider the black hole solution with a regular
universal horizon (though the perturbation equations can be derived without
regard to the value of $b_{0}$). For $b_{0}=b_{0,c}$, one has

$\displaystyle
f(r)=\left(1-\frac{3\mu_{0}}{4r}\right)^{2}\left(1+\frac{\mu_{0}}{2r}+\frac{3\mu_{0}^{2}}{16r^{2}}\right)\geq
0,$ (7)

and the universal horizon is located at $r=3\mu_{0}/4$. Constant $t$
hypersurfaces and the regular universal horizon are depicted in the Penrose
diagram of the Schwarzschild geometry in Fig. 1. For $b_{0}>b_{0,c}$, $f(r)$
is positive everywhere, and therefore the instantaneous mode is accessible to
the singularity at $r=0$. We do not consider this case in this paper.

Figure 1: Colored lines represent $t=\,$const hypersurfaces, and the black
dashed line shows the location of the universal horizon, $r=3\mu_{0}/4$.

Here it should be noted again that gravitational waves propagate at the speed
of light in our TTDOF theory. While the causal boundary for the instantaneous
scalar mode is given by the universal horizon at $r=3\mu_{0}/4$, that for
gravitational waves is still the usual event horizon at $r=\mu_{0}$.

In the rest of the paper, we set $\mu_{0}=1$ to simplify the expressions.

###  II.3 A test scalar field on the black hole background

We are going to discuss the metric perturbations on the black hole background
introduced above, which are composed of the usual tensorial degrees of freedom
and the instantaneous scalar mode. Since the system of the perturbation
equations for the metric is highly involved, it is instructive to consider
here a test scalar field $\pi(t,x^{i})$ obeying an elliptic equation on
constant $t$ hypersurfaces as the simplest model of an instantaneous mode:

$\displaystyle S_{\pi}:=-\frac{1}{2}\int{\rm d}t{\rm
d}^{3}x\sqrt{\gamma}N\nabla_{i}\pi\nabla^{i}\pi.$ (8)

This example helps us to see what will be done in the next two sections.

For the above black hole background, the equation of motion for $\pi$ reads

$\displaystyle
f\pi^{\prime\prime}+f^{\prime}\pi^{\prime}-\left[\frac{\ell(\ell+1)}{r^{2}}+\frac{f^{\prime}}{r}\right]\pi=0,$
(9)

where a prime stands for differentiation with respect to $r$ and the angular
part of the Laplacian was replaced by the eigenvalue $-\ell(\ell+1)/r^{2}$.
The natural inner boundary is seen to be the universal horizon at which $f=0$.
In the case of $\ell=0$, it is easy to find the analytic solution

$\displaystyle\pi(t,r)=r\left[c_{1}(t)+c_{2}(t)\int_{r}^{\infty}\frac{{\rm
d}r}{r^{2}f}\right],$ (10)

where $c_{1}$ and $c_{2}$ are integration functions. For large $r$ we have

$\displaystyle\pi\simeq c_{1}r+c_{2}\left(1+\frac{1}{2r}\right),$ (11)

while in the vicinity of the universal horizon, we have

$\displaystyle\pi\simeq\frac{3c_{2}}{8(r-3/4)}-\frac{c_{2}}{3}\ln(r-3/4).$
(12)

In order for $\pi$ not to diverge at infinity, we require that $c_{1}=0$.
However, if $c_{2}\neq 0$ then $\pi$ diverges at the universal horizon.
Therefore, the only possible regular solution is given by $\pi=0$.

We will perform a similar analysis for metric perturbations. In the case of
metric perturbations, multiple components of the metric are coupled, leading
to a more complicated system. Moreover, the system contains higher spatial
derivatives (for even-parity perturbations with $\ell\geq 1$), allowing for
more linearly independent solutions.

##  III Odd-parity perturbations

In this section and the next, we study linear perturbations around the black
hole solution introduced in the previous section. We first consider odd-parity
perturbations in this section. Only the gravitational-wave degrees of freedom
participate in the odd-parity sector. Noting that the TTDOF theory (1) is
built so that gravitational waves propagate in the same way as in GR, we
expect that the odd-parity sector is identical to that in GR, except that the
background metric is written in the nonstandard coordinate system. Indeed, a
modification from GR appears only in the $K^{2}$ term in action, and it is
easy to see that this term does not include any contribution from odd-parity
perturbations. Below we will see this more explicitly.

Using the spherical harmonics $Y_{\ell m}(\theta,\varphi)$, we expand the odd-
parity perturbations as

$\displaystyle\delta N_{a}$ $\displaystyle=\sum_{\ell,m}h_{0}^{(\ell
m)}(t,r)\epsilon^{b}_{~{}a}\partial_{b}Y_{\ell m},$ (13)
$\displaystyle\delta\gamma_{ra}$ $\displaystyle=\sum_{\ell,m}h_{1}^{(\ell
m)}(t,r)\epsilon^{b}_{~{}a}\partial_{b}Y_{\ell m},$ (14)
$\displaystyle\delta\gamma_{ab}$
$\displaystyle=\frac{1}{2}\sum_{\ell,m}h_{2}^{(\ell
m)}(t,r)\left[\epsilon_{a}^{~{}c}\nabla_{c}\nabla_{b}+\epsilon_{b}^{~{}c}\nabla_{c}\nabla_{a}\right]Y_{\ell
m},$ (15)

where $\nabla_{a}$ is the covariant derivative defined with respect to
$\sigma_{ab}$ and $\epsilon_{ab}$ is the Levi-Civita tensor with
$\epsilon_{\theta\varphi}=\sin\theta$ and
$\epsilon_{\theta}^{~{}\varphi}=\epsilon_{\theta
a}\sigma^{a\varphi}=1/\sin\theta$. Note that $h_{2}^{(\ell m)}=0$ for
$\ell=1$. Under a gauge transformation

$\displaystyle x^{a}\to x^{a}+\xi^{a},\quad\xi_{a}:=\sum_{\ell,m}\xi^{(\ell
m)}(t,r)\epsilon^{b}_{~{}a}\partial_{b}Y_{\ell m},$ (16)

$h_{2}^{(\ell m)}$ transforms as

$\displaystyle h_{2}^{(\ell m)}\to h_{2}^{(\ell m)}+2r^{2}\xi^{(\ell m)}.$
(17)

We thus choose to remove $h_{2}^{(\ell m)}$ for $\ell\geq 2$.

The quadratic action for the odd-parity perturbations can be written in the
form

$\displaystyle S^{\textrm{(odd)}}=\sum_{\ell,m}\int{\rm d}t{\rm
d}r\mathcal{L}_{\ell m}^{\textrm{(odd)}}.$ (18)

Let us consider the quadratic Lagrangian with $\ell\geq 2$ (the $\ell=1$
sector must be treated separately). Thanks to the symmetry, it is sufficient
to study the Lagrangian with $m=0$:

$\displaystyle\mathcal{L}_{\ell\geq 2,m=0}^{\textrm{(odd)}}$
$\displaystyle=\frac{1}{4N_{0}r^{2}}\left[2j^{2}+\frac{c_{\ell}}{f}\right]h_{0}^{2}-\frac{c_{\ell}N_{0}}{4r^{2}}\left(1-\frac{1}{r}\right)h_{1}^{2}+\frac{j^{2}}{4N_{0}}\left[\left(\dot{h}_{1}-h_{0}^{\prime}\right)^{2}+\frac{4}{r}h_{0}\dot{h}_{1}\right]-\frac{c_{\ell}}{2}\frac{b_{0}}{r^{4}\sqrt{f}}h_{0}h_{1},$
(19)

where $j^{2}:=\ell(\ell+1)$ and $c_{\ell}:=(\ell-1)\ell(\ell+1)(\ell+2)$.
Here, a dot and a prime denote derivatives with respect to $t$ and $r$,
respectively. To simplify the expression, we omit the labels $(\ell m)$ from
$h_{0}$ and $h_{1}$. We find no terms that depend on $\beta$. Therefore, the
Lagrangian (19) must coincide with that of the odd-parity perturbations of a
Schwarzschild black hole in GR expressed in the nonstandard coordinate system.
It is obvious that one can arrive at Eq. (19) even if $\beta$ is a function of
time.

To go back to the standard Schwarzschild coordinates, we use the time
coordinate $\tau$ defined in Eq. (6) and define

$\displaystyle\widetilde{h}_{0}:=\frac{h_{0}}{N_{0}},\quad\widetilde{h}_{1}:=h_{1}+\frac{b_{0}/r^{2}}{\sqrt{f}(1-1/r)}\frac{h_{0}}{N_{0}}.$
(20)

Then, the quadratic action for $\ell\geq 2$ and $m=0$ can be written as

$\displaystyle S^{(\textrm{odd})}_{\ell\geq 2,m=0}=\frac{1}{4}\int{\rm
d}\tau{\rm
d}r\left\\{\frac{1}{r^{2}}\left[2j^{2}+\frac{c_{\ell}}{1-1/r}\right]\widetilde{h}_{0}^{2}-\frac{c_{\ell}}{r^{2}}\left(1-\frac{1}{r}\right)\widetilde{h}_{1}^{2}+j^{2}\left[\left(\partial_{\tau}\widetilde{h}_{1}-\widetilde{h}_{0}^{\prime}\right)^{2}+\frac{4}{r}\widetilde{h}_{0}\partial_{\tau}\widetilde{h}_{1}\right]\right\\}.$
(21)

This coincides with the quadratic action for the odd-parity perturbations of
the Schwarzschild solution in GR under an appropriate identification of
variables. Obviously, the same conclusion should hold for the dipole
perturbation, which corresponds to a slow rotation of the black hole.

A comment is now in order. In Ref. [37], quasi-normal frequencies of
gravitational perturbations of a similar black hole in Einstein-Aether theory
were calculated imposing the inner boundary conditions at the universal
horizon. We argue, however, that the odd-parity metric perturbations in the
present case obey the same equation as in GR, i.e. the Regge-Wheeler equation,
and propagate at the speed of light, which indicates that the inner boundary
must be set at the horizon for photons, i.e. the usual event horizon at $r=1$.

##  IV Even-parity perturbations

Let us move to discuss the even-parity perturbations. In usual scalar-tensor
theories in which the scalar field is dynamical, both the gravitational-wave
and scalar degrees of freedom take part in the odd-parity sector. Therefore,
even-parity perturbations with $\ell\geq 2$ are composed of a mixture of
gravitational waves and the scalar, while monopole and dipole perturbations
are composed solely of the scalar degree of freedom. A similar mixing occurs
for $\ell\geq 2$ in a tricky way in the present case where the scalar is an
instantaneous mode obeying an elliptic equation. Monopole and dipole
perturbations in the present case are nondynamical and their configurations
are determined entirely from boundary conditions. Below we will see these
points in more detail.

We write the even-parity perturbations of the ADM variables as

$\displaystyle\delta N$ $\displaystyle=N(r)\sum_{\ell,m}\mathsf{H}_{0}^{(\ell
m)}(t,r)Y_{\ell m},$ (22) $\displaystyle\delta N_{r}$
$\displaystyle=\sum_{\ell,m}\mathsf{H}_{1}^{(\ell m)}(t,r)Y_{\ell m},$ (23)
$\displaystyle\delta\gamma_{rr}$
$\displaystyle=F^{2}(r)\sum_{\ell,m}\mathsf{H}_{2}^{(\ell m)}(t,r)Y_{\ell m},$
(24) $\displaystyle\delta N_{a}$ $\displaystyle=\sum_{\ell,m}\mathsf{b}^{(\ell
m)}(t,r)\partial_{a}Y_{\ell m},$ (25) $\displaystyle\delta\gamma_{ra}$
$\displaystyle=\sum_{\ell,m}\mathsf{a}^{(\ell m)}(t,r)\partial_{a}Y_{\ell m},$
(26) $\displaystyle\delta\gamma_{ab}$
$\displaystyle=r^{2}\sum_{\ell,m}\mathsf{K}^{(\ell m)}(t,r)\sigma_{ab}Y_{\ell
m}+r^{2}\sum_{\ell,m}\mathsf{G}^{(\ell m)}(t,r)\nabla_{a}\nabla_{b}Y_{\ell
m}.$ (27)

Note that $\mathsf{a}^{(00)}=\mathsf{b}^{(00)}=\mathsf{G}^{(00)}=0$. Under
infinitesimal coordinate transformations

$\displaystyle r\to r+\sum_{\ell,m}\xi_{r}^{(\ell m)}(t,r)Y_{\ell m},\quad
x^{a}\to x^{a}+\sum_{\ell,m}\xi_{\Omega}^{(\ell m)}(t,r)\nabla^{a}Y_{\ell m},$
(28)

the perturbation variables transform as

$\displaystyle\mathsf{H}_{0}\to\mathsf{H}_{0}-\frac{N^{\prime}}{N}\xi_{r},\quad\mathsf{H}_{1}\to\mathsf{H}_{1}-(BF\xi_{r})^{\prime}-F^{2}\dot{\xi}_{r},\quad\mathsf{H}_{2}\to\mathsf{H}_{2}-2\xi_{r}^{\prime}-\frac{2F^{\prime}}{F}\xi_{r},$
$\displaystyle\mathsf{b}\to\mathsf{b}-r^{2}\dot{\xi}_{\Omega}-BF\xi_{r},\quad\mathsf{a}\to\mathsf{a}-r^{2}\xi_{\Omega}^{\prime}-F^{2}\xi_{r},\quad\mathsf{K}\to\mathsf{K}-\frac{2}{r}\xi_{r},\quad\mathsf{G}\to\mathsf{G}-2\xi_{\Omega},$
(29)

where the labels $(\ell m)$ were omitted. In the following analysis, we will
use these gauge degrees of freedom to remove some of the variables. Note that
we are working in the unitary gauge and hence we do not have the freedom to
change the time coordinate.

As in the analysis of the odd-parity perturbations, we calculate the quadratic
action,

$\displaystyle S^{\textrm{(even)}}=\sum_{\ell,m}\int{\rm d}t{\rm
d}r\mathcal{L}_{\ell m}^{\textrm{(even)}},$ (30)

and study the perturbations with $\ell=0$, $\ell=1$, and $\ell\geq 2$
separately.

###  IV.1 $\ell=0$

For the monopole perturbations, we are left with
$\mathsf{H}_{0},\mathsf{H}_{1},\mathsf{H}_{2}$, and $\mathsf{K}$, where we
omit the labels $(00)$. From the transformation rules (29) it can be seen that
we can impose the gauge condition $\mathsf{K}=0$. The quadratic Lagrangian for
$\ell=0$ is then given by

$\displaystyle{\cal L}_{\ell=0}^{(\textrm{even})}$ $\displaystyle=\frac{\beta
r^{2}}{12N_{0}\left(\beta+N_{0}\sqrt{f}\right)}\left[\frac{\dot{\mathsf{H}}_{2}}{\sqrt{f}}-\frac{2N_{0}b_{0}}{r^{2}}\widetilde{\mathsf{H}}_{0}^{\prime}-\frac{2}{r^{2}}\left(r^{2}\widetilde{\mathsf{H}}_{1}\right)^{\prime}\right]^{2}-\frac{2b_{0}}{r}\widetilde{\mathsf{H}}_{0}^{\prime}\widetilde{\mathsf{H}}_{1}-N_{0}rf\widetilde{\mathsf{H}}_{0}^{\prime}\mathsf{H}_{2}-\frac{r}{N_{0}\sqrt{f}}\mathsf{H}_{2}\dot{\widetilde{\mathsf{H}}}_{1},$
(31)

where we defined the convenient combinations of the variables as

$\displaystyle\widetilde{\mathsf{H}}_{0}$
$\displaystyle:=\mathsf{H}_{0}+\frac{1}{2}\mathsf{H}_{2},$ (32)
$\displaystyle\widetilde{\mathsf{H}}_{1}$
$\displaystyle:=\sqrt{f}\mathsf{H}_{1}-\frac{N_{0}b_{0}}{r^{2}}(\mathsf{H}_{0}+\mathsf{H}_{2}).$
(33)

The Euler-Lagrange equations for
$\widetilde{\mathsf{H}}_{0},\widetilde{\mathsf{H}}_{1}$, and $\mathsf{H}_{2}$
read

$\displaystyle\left(\frac{2b_{0}}{r}\widetilde{\mathsf{H}}_{1}+N_{0}rf\mathsf{H}_{2}+\frac{2N_{0}b_{0}}{r^{2}}\mathsf{A}\right)^{\prime}$
$\displaystyle=0,$ (34)
$\displaystyle-\frac{2b_{0}}{r}\widetilde{\mathsf{H}}_{0}^{\prime}+\frac{r}{N_{0}\sqrt{f}}\dot{\mathsf{H}}_{2}+2r^{2}\left(\frac{\mathsf{A}}{r^{2}}\right)^{\prime}$
$\displaystyle=0,$ (35) $\displaystyle
N_{0}rf\widetilde{\mathsf{H}}_{0}^{\prime}+\frac{r}{N_{0}\sqrt{f}}\dot{\widetilde{\mathsf{H}}}_{1}+\frac{\dot{\mathsf{A}}}{\sqrt{f}}$
$\displaystyle=0,$ (36)

where $\mathsf{A}$ is defined as

$\displaystyle\mathsf{A}:=\frac{\beta
r^{2}}{6N_{0}(\beta+N_{0}\sqrt{f})}\left[\frac{\dot{\mathsf{H}}_{2}}{\sqrt{f}}-\frac{2N_{0}b_{0}}{r^{2}}\widetilde{\mathsf{H}}_{0}^{\prime}-\frac{2}{r^{2}}\left(r^{2}\widetilde{\mathsf{H}}_{1}\right)^{\prime}\right].$
(37)

Equations (34)–(36) can be solved exactly (for any $\beta=\beta(t)$) as
follows. First, Eq. (34) can be integrated to give

$\displaystyle\frac{2b_{0}}{r}\widetilde{\mathsf{H}}_{1}+N_{0}rf\mathsf{H}_{2}+\frac{2N_{0}b_{0}}{r^{2}}\mathsf{A}=C_{1}(t),$
(38)

where $C_{1}(t)$ is a time-dependent integration function. Using Eqs. (35),
(36), and (38), we obtain

$\displaystyle
2N_{0}^{2}r^{2}f^{3/2}\left(\frac{\mathsf{A}}{r^{2}}\right)^{\prime}=-\dot{C}_{1},$
(39)

which can be integrated to give

$\displaystyle\mathsf{A}=\frac{\dot{C}_{1}r^{2}}{2N_{0}^{2}}\int^{\infty}_{r}\frac{{\rm
d}r}{r^{2}f^{3/2}}+C_{2}(t)r^{2},$ (40)

where $C_{2}(t)$ is another integration function. We then use Eqs. (35) and
(37) to get

$\displaystyle(r^{2}\widetilde{\mathsf{H}}_{1})^{\prime}=-N_{0}r^{3}\left(\frac{\mathsf{A}}{r^{2}}\right)^{\prime}-\frac{3N_{0}}{\beta}(\beta+N_{0}\sqrt{f})\mathsf{A}.$
(41)

We impose the boundary conditions
$\mathsf{H}_{0},\mathsf{H}_{1},\mathsf{H}_{2}\to 0$ at infinity. Then, noting
that $\mathsf{A}\approx(\dot{C}_{1}/2N_{0}^{2})r+C_{2}r^{2}$ for large $r$, we
see from Eq. (41) that $\dot{C}_{1}=C_{2}=0$, i.e., $\mathsf{A}=0$, though
$C_{1}$ can still be a nonvanishing, time-independent constant. Thus, we
obtain

$\displaystyle\frac{\widetilde{\mathsf{H}}_{1}}{\sqrt{4\pi}}=\frac{N_{0}C_{0}(t)}{r^{2}},$
(42)

where $C_{0}(t)$ is an integration function. (The factor $1/\sqrt{4\pi}$ comes
from $Y_{00}$.) It follows from Eq. (38) that

$\displaystyle\frac{\mathsf{H}_{2}}{\sqrt{4\pi}}=\frac{1}{f}\left[\frac{\delta\mu}{r}-\frac{2b_{0}C_{0}(t)}{r^{4}}\right],$
(43)

where now we write $C_{1}=\sqrt{4\pi}N_{0}\delta\mu$ with $\delta\mu$ being a
time-independent constant. Finally, we have

$\displaystyle\frac{\widetilde{\mathsf{H}}_{0}}{\sqrt{4\pi}}=\frac{\dot{C}_{0}}{N_{0}}\int^{\infty}_{r}\frac{{\rm
d}r}{r^{2}f^{3/2}}.$ (44)

All the perturbation variables are regular at the usual event horizon, $r=1$.
Notice, however, that $\widetilde{\mathsf{H}}_{0}$ diverges as
$\widetilde{\mathsf{H}}_{0}\sim(r-3/4)^{-2}$ at the universal horizon. To
avoid this, we are forced to set $\dot{C}_{0}=0$. Now it is easy to see that
the above solution for the monopole perturbations can be reproduced by
perturbing the parameters of the background solution as
$\mu_{0}\to\mu_{0}+\delta\mu$ and $b_{0}\to b_{0}+C_{0}$, which shows that
$\delta\mu$ simply corresponds to a shift of the mass parameter and
nonvanishing $C_{0}$ renders the universal horizon singular by enforcing
$b_{0}\neq b_{0,c}$. Therefore, we conclude that no nontrivial regular
solution for the monopole perturbations exists.

Let us look at the above solution from the viewpoint of the Stückelberg field.
All the metric perturbations with the coefficients $C_{0}$ and $\dot{C}_{0}$
can be eliminated by performing the coordinate transformation

$\displaystyle t\rightarrow
T=t-\frac{C_{0}(t)}{N_{0}}\int^{\infty}_{r}\frac{{\rm d}r}{r^{2}f^{3/2}}.$
(45)

We thus see that the geometry remains Schwarzschild (with the mass parameter
$1+\delta\mu$), while one has the fluctuation of the Stückelberg field,
$\phi=T+\delta\phi$, where $\delta\phi$ is the minus of the second term in Eq.
(45). However, $\phi$ is singular at $r=3/4$ unless $C_{0}(t)=0$, and the same
conclusion follows.

Let us give a comment on the difference between the present analysis and that
in Ref. [21]. In Ref. [21], the authors investigated static monopole
deformations of the same black hole solution in the same theory without paying
particular attention to the regularity of the inner boundary, while in the
present paper we have started from generic time-dependent monopole
perturbations and showed explicitly how any time-dependent deformation is
prohibited by the boundary condtions. The result obtained here is not trivial
as the fluctuation of the instantaneous scalar could in principle be time-
dependent.

###  IV.2 $\ell=1$

In the quadratic Lagrangian for the dipole perturbations, $\mathsf{K}^{(1m)}$
and $\mathsf{G}^{(1m)}$ appears only through the combination
$\mathsf{K}^{(1m)}-\mathsf{G}^{(1m)}$. We may therefore impose the gauge
condition $\mathsf{K}^{(1m)}-\mathsf{G}^{(1m)}=0$ along with
$\mathsf{H}_{0}^{(1m)}=0$ by choosing $\xi_{r}^{(1m)}$ and
$\xi_{\Omega}^{(1m)}$ appropriately. Focusing on the $m=0$ part without loss
of generality, the quadratic Lagrangian for $\ell=1$ is found to be

$\displaystyle{\cal L}_{\ell=1,m=0}^{(\textrm{even})}$
$\displaystyle=\frac{\beta
r^{2}}{12N_{0}(\beta+N_{0}\sqrt{f})}\left[\frac{\dot{\mathsf{H}}_{2}}{\sqrt{f}}-\frac{N_{0}b_{0}}{r^{2}}\mathsf{H}_{2}^{\prime}-\frac{2}{r^{2}}\left(r^{2}\widetilde{\mathsf{H}}_{1}\right)^{\prime}+\frac{4\widetilde{\mathsf{b}}}{r}+\frac{6(\beta+N_{0}\sqrt{f})}{\beta
r}\left(\widetilde{\mathsf{H}}_{1}-\widetilde{\mathsf{b}}\right)\right]^{2}$
$\displaystyle\quad+\frac{r^{2}}{2N_{0}}\left[\dot{\widetilde{\mathsf{a}}}-rN_{0}b_{0}\left(\frac{\sqrt{f}}{r^{3}}\widetilde{\mathsf{a}}\right)^{\prime}-r\left(\frac{\sqrt{f}}{r}\widetilde{\mathsf{b}}\right)^{\prime}-\frac{\widetilde{\mathsf{H}}_{1}}{r\sqrt{f}}\right]^{2}-\frac{3\sqrt{f}}{\beta}\left(\widetilde{\mathsf{H}}_{1}-\widetilde{\mathsf{b}}\right)^{2}+\frac{6b_{0}}{r^{2}}\widetilde{\mathsf{H}}_{1}\widetilde{\mathsf{a}}-\frac{3b_{0}}{r^{2}}\mathsf{H}_{2}\widetilde{\mathsf{b}}$
$\displaystyle\quad+N_{0}f\widetilde{\mathsf{a}}^{2}-\frac{N_{0}}{2}\left[f+(rf)^{\prime}\right]\widetilde{\mathsf{a}}\mathsf{H}_{2}+\frac{N_{0}}{4}(rf)^{\prime}\mathsf{H}_{2}^{2},$
(46)

where to simplify the expression we defined

$\displaystyle\widetilde{\mathsf{a}}:=\frac{\mathsf{a}}{r},\quad\widetilde{\mathsf{b}}:=\frac{\mathsf{b}}{r\sqrt{f}}-\frac{N_{0}b_{0}}{r^{3}}\mathsf{a},\quad\widetilde{\mathsf{H}}_{1}:=\sqrt{f}\mathsf{H}_{1}-\frac{N_{0}b_{0}}{r^{2}}\mathsf{H}_{2},$
(47)

and omitted the labels $(10)$.

We then introduce new variables $\mathsf{P}$ and $\mathsf{Q}$ and define the
Lagrangian $\widetilde{\mathcal{L}}^{(2)}_{\ell=1,m=0}$ as

$\displaystyle\widetilde{\mathcal{L}}^{(2)}_{\ell=1,m=0}$
$\displaystyle=\mathcal{L}^{(2)}_{\ell=1,m=0}$ $\displaystyle\quad-\frac{\beta
r^{2}}{12N_{0}(\beta+N_{0}\sqrt{f})}\left[\frac{\dot{\mathsf{H}}_{2}}{\sqrt{f}}-\frac{N_{0}b_{0}}{r^{2}}\mathsf{H}_{2}^{\prime}-\frac{2}{r^{2}}\left(r^{2}\widetilde{\mathsf{H}}_{1}\right)^{\prime}+\frac{4\widetilde{\mathsf{b}}}{r}+\frac{6(\beta+N_{0}\sqrt{f})}{\beta
r}\left(\widetilde{\mathsf{H}}_{1}-\widetilde{\mathsf{b}}\right)-\frac{\beta+N_{0}\sqrt{f}}{\beta}\mathsf{Q}\right]^{2}$
$\displaystyle\quad+\frac{r^{2}}{2N_{0}}\left[\dot{\widetilde{\mathsf{a}}}-rN_{0}b_{0}\left(\frac{\sqrt{f}}{r^{3}}\widetilde{\mathsf{a}}\right)^{\prime}-r\left(\frac{\sqrt{f}}{r}\widetilde{\mathsf{b}}\right)^{\prime}-\frac{\widetilde{\mathsf{H}}_{1}}{r\sqrt{f}}-\mathsf{P}\right]^{2}.$
(48)

This Lagrangian is equivalent to the original one
$\mathcal{L}^{(2)}_{\ell=1,m=0}$ because the additional parts vanish upon
substituting the solutions to the equations of motion for $\mathsf{P}$ and
$\mathsf{Q}$. However, the new Lagrangian
$\widetilde{\mathcal{L}}^{(2)}_{\ell=1,m=0}$ is more useful for our anaysis.
Varying $\widetilde{\mathcal{L}}^{(2)}_{\ell=1,m=0}$ with respect to
$\widetilde{\mathsf{H}}_{1}$, $\mathsf{H}_{2}$, $\widetilde{\mathsf{a}}$, and
$\widetilde{\mathsf{b}}$, we obtain the equations of motion for these
variables, which turn out to be the constraint equations. One can easily solve
them to express $\widetilde{\mathsf{H}}_{1}$, $\mathsf{H}_{2}$,
$\widetilde{\mathsf{a}}$, and $\widetilde{\mathsf{b}}$ in terms of
$\mathsf{P}$, $\mathsf{Q}$, and their first derivatives:
$\widetilde{\mathsf{H}}_{1}=(\dots)\mathsf{P}+(\dots)\mathsf{Q}+(\dots)\dot{\mathsf{P}}+(\dots)\dot{\mathsf{Q}}+(\dots)\mathsf{P}^{\prime}+(\dots)\mathsf{Q}^{\prime}$,
$\mathsf{H}_{2},\widetilde{\mathsf{a}},\widetilde{\mathsf{b}}=\dots$, where
the explicit expressions are messy. Substituting these back to
$\widetilde{\mathcal{L}}^{(2)}_{\ell=1,m=0}$, we can express it in terms of
$\mathsf{P}$, $\mathsf{Q}$, and their derivatives. We then introduce the new
variable $\chi:=\mathsf{Q}+3\sqrt{f}\mathsf{P}$ and replace $\mathsf{Q}$ in
the Lagrangian by $\chi$. By doing so we can remove all the derivatives acting
on $\mathsf{P}$. The equation of motion for $\mathsf{P}$ can therefore be
solved to express it in terms of $\chi$ and its derivatives. Substituting the
solution back to the Lagrangian to remove $\mathsf{P}$ and performing
integration by parts, we finally arrive at the reduced Lagrangian for a single
master variable $\chi$. Despite lengthy expressions at each intermediate step,
the final form of the reduced Lagrangian is rather simple:

$\displaystyle\widetilde{\mathcal{L}}^{(2)}_{\ell=1,m=0}=\frac{\beta(t)}{81N_{0}^{2}}\left[d_{1}(r)(\chi^{\prime\prime})^{2}+d_{2}(r)(\chi^{\prime})^{2}+d_{3}(r)\chi^{2}\right],$
(49)

where the coefficients $d_{1}$, $d_{2}$, and $d_{3}$ are independent of
$\beta(t)$ and are given explicitly by

$\displaystyle d_{1}(r)$ $\displaystyle=\frac{r^{8}f^{3/2}}{3},$ (50)
$\displaystyle d_{2}(r)$
$\displaystyle=-2r^{6}f^{1/2}\left(2-\frac{1}{r}-\frac{5b_{0}^{2}}{r^{4}}\right),$
(51) $\displaystyle d_{3}(r)$
$\displaystyle=\frac{r^{4}}{f^{1/2}}\left(\frac{2}{r}-\frac{20b_{0}^{2}}{r^{4}}+\frac{21b_{0}^{2}}{r^{5}}\right).$
(52)

(We are primarily interested in the case where $b_{0}=b_{0,c}$, but the above
expressions are valid for any $b_{0}$.) Note that there are no time
derivatives of $\chi$ in the Lagrangian, which means that there is no
propagating dipole mode.

The equation of motion for $\chi$ is given by the fourth-order differential
equation with respect to $r$,

$\displaystyle\left(d_{1}\chi^{\prime\prime}\right)^{\prime\prime}-\left(d_{2}\chi^{\prime}\right)^{\prime}+d_{3}\chi=0,$
(53)

and the configuration of $\chi$ is determined once one specifies the boundary
conditions. As implied by Eqs. (50)–(52), the inner boundary conditions must
be imposed at the location at which $f$ vanishes, i.e., the universal horizon.
All the other variables can then be obtained straightforwardly from $\chi$.

We have not found an analytic solution to Eq. (53). However, we can obtain
solutions valid near the boundaries. Let us first look for a solution valid
for large $r$ in the form

$\displaystyle\chi=\mathcal{C}_{0}(t)+\frac{\mathcal{C}_{1}(t)}{r}+\frac{\mathcal{C}_{2}(t)}{r^{2}}+\ldots+\log
r\left[\mathcal{D}_{0}(t)+\frac{\mathcal{D}_{1}(t)}{r}+\frac{\mathcal{D}_{2}(t)}{r^{2}}+\dots\right].$
(54)

Substituting this into Eq. (53), one can derive the algebraic relations among
the coefficients. One may take $\mathcal{C}_{0}$, $\mathcal{C}_{2}$,
$\mathcal{C}_{3}$, and $\mathcal{C}_{5}$ as independent coefficients and
express all the other coefficients in terms of them. One can then obtain the
metric perturbations $\mathsf{H}_{1}$, $\mathsf{H}_{2}$, $\mathsf{a}$, and
$\mathsf{b}$ in a series expansion similar to Eq. (54). In order for the
metric perturbations to be vanishing at infinity, it turns out that
$\mathcal{C}_{0}=\mathcal{C}_{2}=\mathcal{C}_{3}=0=\ddot{\mathcal{C}}_{5}$
must be imposed. Then, $\mathcal{C}_{5}$ may have the form
$\mathcal{C}_{5}=v_{1}t+v_{0}$, where $v_{1}$ and $v_{0}$ are constants,
showing that $\chi$ diverges as $t\to\infty$ unless
$\dot{\mathcal{C}}_{5}=v_{1}=0$. We thus see that only a constant
$\mathcal{C}_{5}$ is allowed to be nonvanishing from the boundary conditions
at infinity and $\chi$ is given by

$\displaystyle\chi=\frac{\mathcal{C}_{5}}{r^{5}}\left(1+\frac{3}{4r}+\frac{9}{16r^{2}}+\frac{7}{16r^{3}}+\dots\right).$
(55)

In this way, $\chi$ is uniquely determined up to an overall constant. The
dipole metric perturbations for large $r$ are given accordingly by

$\displaystyle\widetilde{\mathsf{H}}_{1}$
$\displaystyle=\frac{\mathcal{C}_{5}}{r^{3}}\left(\frac{2}{9}+\frac{1}{12}\frac{1}{r}+\frac{1}{24}\frac{1}{r^{2}}+\dots\right),$
(56) $\displaystyle\mathsf{H}_{2}$
$\displaystyle=\frac{b_{0}}{N_{0}\mu_{0}}\frac{\mathcal{C}_{5}}{r^{4}}\left(\frac{8}{3}+\frac{4}{3}\frac{1}{r}+\frac{5}{6}\frac{1}{r^{2}}+\dots\right),$
(57) $\displaystyle\widetilde{\mathsf{a}}$
$\displaystyle=\frac{b_{0}}{N_{0}}\frac{\mathcal{C}_{5}}{r^{4}}\left(\frac{4}{3}+\frac{1}{3}\frac{1}{r}+\dots\right),$
(58) $\displaystyle\widetilde{\mathsf{b}}$
$\displaystyle=-\frac{\mathcal{C}_{5}}{r^{3}}\left(\frac{1}{9}+\frac{1}{12}\frac{1}{r}+\frac{1}{16}\frac{1}{r^{2}}+\dots\right),$
(59)

where $\beta$ does not appear in the series expansion.

Next, let us examine the behavior of the dipole perturbations near the
(regular) universal horizon at $r=3/4$, setting now $b_{0}=b_{0,c}$. Near the
universal horizon, we have $d_{1}\simeq(81/256\sqrt{2})\cdot(r-3/4)^{3}$,
$d_{2}\simeq(81/256\sqrt{2})\cdot 3(r-3/4)$, and
$d_{3}\simeq(81/256\sqrt{2})\cdot 4(r-3/4)^{-1}$. It then follows that, near
the universal horizon, $\chi$ is given by a linear combination of the four
independent solutions

$\displaystyle(r-3/4)^{\sqrt{2}},\quad(r-3/4)^{\sqrt{2}}\ln(r-3/4),\quad(r-3/4)^{-\sqrt{2}},\quad(r-3/4)^{-\sqrt{2}}\ln(r-3/4).$
(60)

By numerically solving Eq. (53) from some large $r$ toward the universal
horizon, we find that the coefficients of the latter two are nonvanishing and
$\chi$ diverges, as can be seen from Fig. 2. This forces the metric
perturbations reconstructed from $\chi$ to diverge as

$\displaystyle(BF)^{-1}\mathsf{H}_{1}\sim\mathsf{H}_{2}\sim\mathsf{a}\sim\mathsf{b}\sim(r-3/4)^{-1}\chi.$
(61)

One can also see that the linear perturbation of the three-dimensional Ricci
scalar diverges as $\delta R\sim(r-3/4)^{-1}\chi$. We, therefore, conclude
that no regular dipole perturbations are allowed in the present setup.

Figure 2: The behavior of $\chi$ near the universal horizon, obtained by
solving Eq. (53) numerically with the boundary condition (55) (with
$\mathcal{C}_{5}=1$).

###  IV.3 $\ell\geq 2$

Finally, let us consider the even-parity perturbations with $\ell\geq 2$,
which contains both gravitational waves and instantaneous mode. By choosing
appropriately $\xi_{r}^{(\ell m)}$ and $\xi_{\Omega}^{(\ell m)}$, one can set
$\mathsf{K}^{(\ell m)}=\mathsf{G}^{(\ell m)}=0$. The remaining variables are
$\mathsf{H}_{0}^{(\ell m)}$, $\mathsf{H}_{1}^{(\ell m)}$,
$\mathsf{H}_{2}^{(\ell m)}$, $\mathsf{a}^{(\ell m)}$, and $\mathsf{b}^{(\ell
m)}$. Focusing again on the perturbations with $m=0$ and omitting the labels
$(\ell m)$, the quadratic Lagrangian is given by

$\displaystyle{\cal L}_{\ell\geq 2,m=0}^{(\textrm{even})}$
$\displaystyle=\frac{\beta
r^{2}}{12N_{0}(\beta+N_{0}\sqrt{f})}\biggl{[}\frac{\dot{\mathsf{H}}_{2}}{\sqrt{f}}-\frac{N_{0}b_{0}}{r^{2}}\mathsf{H}_{2}^{\prime}-\frac{2}{r^{2}}\left(r^{2}\widetilde{\mathsf{H}}_{1}\right)^{\prime}+\frac{6(\beta+N_{0}\sqrt{f})}{\beta
r}\widetilde{\mathsf{H}}_{1}$ $\displaystyle\quad-\frac{6N_{0}b_{0}}{\beta
r^{3}}(\beta+N_{0}\sqrt{f})\mathsf{H}_{0}-\frac{j^{2}(\beta+3N_{0}\sqrt{f})}{\beta
r}\widetilde{\mathsf{b}}\biggr{]}^{2}+\frac{j^{2}r^{2}}{4N_{0}}\left[\dot{\widetilde{\mathsf{a}}}-rN_{0}b_{0}\left(\frac{\sqrt{f}}{r^{3}}\widetilde{\mathsf{a}}\right)^{\prime}-r\left(\frac{\sqrt{f}}{r}\widetilde{\mathsf{b}}\right)^{\prime}-\frac{\widetilde{\mathsf{H}}_{1}}{r\sqrt{f}}\right]^{2}$
$\displaystyle\quad-\frac{3\sqrt{f}}{\beta}\left(\widetilde{\mathsf{H}}_{1}-\frac{N_{0}b_{0}}{r^{2}}\mathsf{H}_{0}-\frac{j^{2}}{2}\widetilde{\mathsf{b}}\right)^{2}+\frac{3j^{2}b_{0}}{r^{2}}\widetilde{\mathsf{H}}_{1}\widetilde{\mathsf{a}}-\frac{3j^{2}b_{0}}{2r^{2}}\mathsf{H}_{2}\widetilde{\mathsf{b}}+\frac{j^{2}}{2}N_{0}f\widetilde{\mathsf{a}}^{2}-\frac{j^{2}N_{0}}{4}\left[f+(rf)^{\prime}\right]\widetilde{\mathsf{a}}\mathsf{H}_{2}$
$\displaystyle\quad+\frac{N_{0}}{4}(rf)^{\prime}\mathsf{H}_{2}^{2}+\frac{(\ell-1)\ell(\ell+1)(\ell+2)}{4N_{0}}\widetilde{\mathsf{b}}^{2}+N_{0}f\left[r\mathsf{H}_{2}^{\prime}-j^{2}(r\widetilde{\mathsf{a}})^{\prime}\right]\mathsf{H}_{0}$
$\displaystyle\quad+N_{0}\left[\frac{j^{2}}{2}+(rf)^{\prime}\right]\mathsf{H}_{0}\mathsf{H}_{2}-\frac{j^{2}}{2}N_{0}[f+(rf)^{\prime}]\widetilde{\mathsf{a}}\mathsf{H}_{0},$
(62)

where $\widetilde{\mathsf{H}}_{1}$, $\widetilde{\mathsf{a}}$, and
$\widetilde{\mathsf{b}}$ are defined in the same way as in the case of
$\ell=1$ (see Eq. (47)) and $j^{2}:=\ell(\ell+1)$ as before.

Since we only have one of the two tensorial modes as a dynamical degree of
freedom in the even-parity sector with $\ell\geq 2$ and there is no dynamical
scalar degree of freedom, we expect that by integrating out nondynamical
variables in the Lagrangian (62) we would end up with a reduced Lagrangian
written in terms of a single master variable. Unfortunately, however, we have
not been able to do so. Therefore, in what follows we focus on stationary
perturbations by dropping all time derivatives and solve directly the
equations of motion derived from the Lagrangian (62) without trying to rewrite
the system in terms of a single variable. (In doing so we also assume that
$\beta=\,$const.) This amounts to considering stationary deformations of the
black hole while discarding propagating gravitational waves. Varying the
action, one can straightforwardly obtain the equations of motion for
$\mathsf{H}_{0}$, $\widetilde{\mathsf{H}}_{1}$, $\mathsf{H}_{2}$,
$\widetilde{\mathsf{a}}$, and $\widetilde{\mathsf{b}}$:

$\displaystyle\frac{\delta S^{\textrm{(even)}}}{\delta\mathsf{H}_{0}^{(\ell
0)}}=\frac{\delta S^{\textrm{(even)}}}{\delta\widetilde{\mathsf{H}}_{1}^{(\ell
0)}}=\frac{\delta S^{\textrm{(even)}}}{\delta\mathsf{H}_{2}^{(\ell
0)}}=\frac{\delta S^{\textrm{(even)}}}{\delta\widetilde{\mathsf{a}}^{(\ell
0)}}=\frac{\delta S^{\textrm{(even)}}}{\delta\widetilde{\mathsf{b}}^{(\ell
0)}}=0.$ (63)

First, we determine the behavior of the metric perturbations at large $r$,
assuming the series expansion form

$\displaystyle\mathsf{H}_{0}$
$\displaystyle=r^{-\ell}\left(\mathcal{C}_{0}^{(0)}+\frac{\mathcal{C}_{1}^{(0)}}{r}+\frac{\mathcal{C}_{2}^{(0)}}{r^{2}}+\dots\right),\quad\widetilde{\mathsf{H}}_{1}=r^{-\ell}\left(\mathcal{C}_{0}^{(1)}+\frac{\mathcal{C}_{1}^{(1)}}{r}+\frac{\mathcal{C}_{2}^{(1)}}{r^{2}}+\dots\right),$
$\displaystyle\mathsf{H}_{2}$
$\displaystyle=r^{-\ell}\left(\mathcal{C}_{0}^{(2)}+\frac{\mathcal{C}_{1}^{(2)}}{r}+\frac{\mathcal{C}_{2}^{(2)}}{r^{2}}+\dots\right),\quad\widetilde{\mathsf{a}}=r^{-\ell}\left(\mathcal{C}_{0}^{(a)}+\frac{\mathcal{C}_{1}^{(a)}}{r}+\frac{\mathcal{C}_{2}^{(2)}}{r^{2}}+\dots\right),$
$\displaystyle\widetilde{\mathsf{b}}$
$\displaystyle=r^{-\ell}\left(\mathcal{C}_{0}^{(b)}+\frac{\mathcal{C}_{1}^{(b)}}{r}+\frac{\mathcal{C}_{2}^{(b)}}{r^{2}}+\dots\right).$
(64)

Substituting these to the equations of motion (63), we can derive the
algebraic relations among the coefficients
$\mathcal{C}_{0}^{(0)},\mathcal{C}_{0}^{(1)},\dots$. By inspecting the
relations, we find that only two of the coefficients are free and independent,
and all the other coefficients are expressed using the two. Explicitly, we
have

$\displaystyle\mathsf{H}_{0}$
$\displaystyle=\frac{\mathcal{B}_{0}}{r^{\ell+4}}\frac{b_{0}}{N_{0}}\left[1+\frac{2\ell^{2}+5\ell+4}{4(\ell+1)}\frac{1}{r}+\frac{2\ell^{3}+11\ell^{2}+21\ell+16}{16(\ell+1)}\frac{1}{r^{2}}+\dots\right]$
$\displaystyle\quad-\frac{1}{2(\ell+1)}\frac{\mathcal{B}_{2}}{r^{\ell+1}}\left[1+\frac{\ell+2}{2}\frac{1}{r}+\dots\right],$
(65) $\displaystyle\widetilde{\mathsf{H}}_{1}$
$\displaystyle=\frac{\mathcal{B}_{0}}{r^{\ell+2}}\left[1+\frac{\ell(2\ell+1)}{4(\ell+1)}\frac{1}{r}+\frac{\ell(2\ell+1)}{16}\frac{1}{r^{2}}+\dots\right]$
$\displaystyle\quad-\frac{\mathcal{B}_{2}}{r^{\ell+3}}N_{0}b_{0}\left[\frac{\ell+2}{4(\ell+1)}+\frac{4\ell^{4}+35\ell^{3}+85\ell^{2}+98\ell+48}{16(\ell+1)(\ell+2)(2\ell+3)}\frac{1}{r}+\dots\right],$
(66) $\displaystyle\mathsf{H}_{2}$
$\displaystyle=-\frac{2\mathcal{B}_{0}}{r^{\ell+3}}\frac{b_{0}}{N_{0}}\left[1+\frac{2\ell^{2}+5\ell+4}{4(\ell+1)}\frac{1}{r}+\frac{2\ell^{3}+11\ell^{2}+21\ell+16}{16(\ell+1)}\frac{1}{r^{2}}+\dots\right]$
$\displaystyle\quad+\frac{\mathcal{B}_{2}}{r^{\ell+1}}\left[1+\frac{\ell^{3}+5\ell^{2}+7\ell+4}{2(\ell+1)(\ell+2)}\frac{1}{r}+\dots\right],$
(67) $\displaystyle\widetilde{\mathsf{a}}$
$\displaystyle=\frac{1}{\ell+1}\frac{\mathcal{B}_{0}}{r^{\ell+4}}\frac{b_{0}}{N_{0}}\left[1+\frac{2\ell+5}{4}\frac{1}{r}+\frac{2\ell^{2}+13\ell+22}{16}\frac{1}{r^{2}}+\dots\right]$
$\displaystyle\quad-\frac{\mathcal{B}_{2}}{2r^{\ell+1}}\left[\frac{1}{\ell+1}+\frac{\ell+4}{2(\ell+2)}+\dots\right],$
(68) $\displaystyle\widetilde{\mathsf{b}}$
$\displaystyle=\frac{1}{\ell+1}\frac{\mathcal{B}_{0}}{r^{\ell+2}}\left[1+\frac{2\ell+1}{4}\frac{1}{r}+\frac{(2\ell+1)(\ell+2)}{16}\frac{1}{r^{2}}+\dots\right]$
$\displaystyle\quad-\frac{\mathcal{B}_{2}}{r^{\ell+3}}N_{0}b_{0}\left[\frac{1}{4(\ell+1)}+\frac{4\ell^{2}+13\ell+18}{16(\ell+1)(2\ell+3)}\frac{1}{r}+\dots\right],$
(69)

where we redefined the two independent coefficients and introduced the
constants $\mathcal{B}_{0}$ and $\mathcal{B}_{2}$. The above expressions are
valid for any $b_{0}$.

To see the physical interpretation of the two constants $\mathcal{B}_{0}$ and
$\mathcal{B}_{2}$, let us perform similar calculations for the Schwarzschild
solution, ${\rm d}s^{2}=-(1-1/r){\rm d}t^{2}+(1-1/r)^{-1}{\rm
d}r^{2}+r^{2}\sigma_{ab}{\rm d}x^{a}{\rm d}x^{b}$, in GR. The quadratic
Lagrangian for the even-parity perturbations of the Schwarzschild solution in
GR can be reproduced by setting $N_{0}=1$, $b_{0}=0$, and $\beta=0$ in the
Lagrangian (62) and removing $\mathsf{b}$ with an infinitesimal coordinate
transformation of $t$:

$\displaystyle{\cal L}_{\ell\geq 2,m=0}^{(\textrm{even, GR})}$
$\displaystyle=r\left[\frac{\dot{\mathsf{H}}_{2}}{\sqrt{f}}-\frac{2}{r^{2}}\left(r^{2}\widetilde{\mathsf{H}}_{1}\right)^{\prime}\right]\widetilde{\mathsf{H}}_{1}+\frac{j^{2}r^{2}}{4}\left[\dot{\widetilde{\mathsf{a}}}-\frac{\widetilde{\mathsf{H}}_{1}}{r\sqrt{f}}\right]^{2}+\frac{j^{2}}{2}f\widetilde{\mathsf{a}}^{2}-\frac{j^{2}}{4}\left[f+(rf)^{\prime}\right]\widetilde{\mathsf{a}}\mathsf{H}_{2}+\frac{1}{4}(rf)^{\prime}\mathsf{H}_{2}^{2}$
$\displaystyle\quad+f\left[r\mathsf{H}_{2}^{\prime}-j^{2}(r\widetilde{\mathsf{a}})^{\prime}\right]\mathsf{H}_{0}+\left[\frac{j^{2}}{2}+(rf)^{\prime}\right]\mathsf{H}_{0}\mathsf{H}_{2}-\frac{j^{2}}{2}[f+(rf)^{\prime}]\widetilde{\mathsf{a}}\mathsf{H}_{0},$
(70)

with $f=1-1/r$. Similar manipulations as described above lead to

$\displaystyle\mathsf{H}_{0}$
$\displaystyle=-\frac{1}{2(\ell+1)}\frac{\bar{\mathcal{B}}_{2}}{r^{\ell+1}}\left[1+\frac{\ell+2}{2}\frac{1}{r}+\dots\right],$
(71) $\displaystyle\widetilde{\mathsf{H}}_{1}$ $\displaystyle=0,$ (72)
$\displaystyle\mathsf{H}_{2}$
$\displaystyle=\frac{\bar{\mathcal{B}}_{2}}{r^{\ell+1}}\left[1+\frac{\ell^{3}+5\ell^{2}+7\ell+4}{2(\ell+1)(\ell+2)}\frac{1}{r}+\dots\right],$
(73) $\displaystyle\widetilde{\mathsf{a}}$
$\displaystyle=-\frac{\bar{\mathcal{B}}_{2}}{2r^{\ell+1}}\left[\frac{1}{\ell+1}+\frac{\ell+4}{2(\ell+2)}\frac{1}{r}+\dots\right],$
(74)

where it is found that we are allowed to have only one integration constant
$\bar{\mathcal{B}}_{2}$ in the case of GR. Equations (71)–(74) agree with Eqs.
(65)–(69) with $b_{0}=0=\mathcal{B}_{0}$, if one identifies $\mathcal{B}_{2}$
as $\bar{\mathcal{B}}_{2}$. This observation clearly shows that the terms
proportional to $\mathcal{B}_{2}$ come from the tensorial degrees of freedom,
while those with $\mathcal{B}_{0}$ originate from the instantaneous scalar
mode.

Having determined the solution for large $r$, let us investigate the behavior
of the perturbations near the inner boundaries, focusing on the case of the
regular universal horizon, $b_{0}=b_{0,c}$. If one sets $\mathcal{B}_{2}=0$ to
kill the tensorial degree of freedom and integrates the equations of motion
inwards, the metric perturbations are found to diverge at the universal
horizon, $r=3/4$, as presented in Fig. 4. This result shows that the terms
with the coefficient $\mathcal{B}_{0}$ are indeed induced by the instantaneous
mode, whose causal boundary is the universal horizon. By inspecting the
equations of motion in the vicinity of $r=3/4$, we find that the metric
perturbations diverge as

$\displaystyle\widetilde{\mathsf{H}}_{1}$
$\displaystyle\simeq\varepsilon^{1-k_{\ell}}(c+c^{\prime}\ln\varepsilon),\quad\mathsf{H}_{2}\simeq-\frac{2\sqrt{3}}{N_{0}}\widetilde{\mathsf{H}}_{1},\quad\mathsf{H}_{0}\simeq\frac{\sqrt{3}}{N_{0}}\varepsilon^{1-k_{\ell}}\left[c-\frac{4\sqrt{2}\beta}{N_{0}}\frac{(2-k_{\ell})(3-k_{\ell})}{\ell(\ell+1)}c^{\prime}+c^{\prime}\ln\varepsilon\right],$
$\displaystyle\widetilde{\mathsf{a}}$
$\displaystyle\simeq-\frac{3\sqrt{3}}{8N_{0}(1-k_{\ell})}\varepsilon^{-k_{\ell}}\left(c-\frac{c^{\prime}}{1-k_{\ell}}+c^{\prime}\ln\varepsilon\right),$
$\displaystyle\widetilde{\mathsf{b}}$
$\displaystyle\simeq\frac{4}{1-k_{\ell}}\varepsilon^{2-k_{\ell}}\left[c-c^{\prime}\left(\frac{1}{1-k_{\ell}}+\frac{2\sqrt{2}\beta}{N_{0}}\right)+\frac{2(1-k_{\ell})(3-k_{\ell})}{\ell(\ell+1)}c^{\prime}\ln\varepsilon\right],$
(75)

where $\varepsilon:=r-3/4$, $k_{\ell}:=2+\sqrt{1+\ell(\ell+1)/2}$, and $c$ and
$c^{\prime}$ are constants.

If one sets $\mathcal{B}_{2}\neq 0$ and integrates the equations of motion
inwards, one finds that the metric perturbations diverge at the usual event
horizon, $r=1$, before reaching the universal horizon, as shown in Fig. 4.
Recalling that the tensorial degrees of freedom propagate at the speed of
light in the present theory, it is reasonable that the inner boundary in this
case is given by the usual event horizon. In this case, one can study the
behavior of the metric perturbations near the event horizon by expanding the
equations of motion around $r=1$ to obtain

$\displaystyle\widetilde{\mathsf{H}}_{1}$
$\displaystyle\simeq\bar{c}(r-1)^{-1},\quad\mathsf{H}_{2}\simeq\frac{32\bar{c}}{3\sqrt{3}N_{0}}(r-1)^{-1},\quad\mathsf{H}_{0}\simeq-\frac{1184\bar{c}}{81\sqrt{3}N_{0}}\ln{(r-1)},$
$\displaystyle\widetilde{\mathsf{a}}$
$\displaystyle\simeq\frac{32\bar{c}}{3\sqrt{3}\ell(\ell+1)N_{0}}(r-1)^{-1},\quad\widetilde{\mathsf{b}}\simeq\frac{2\bar{c}}{\ell(\ell+1)}(r-1)^{-1},$
(76)

where $\bar{c}$ is a constant. In any case the stationary perturbations
diverge at the inner boundaries, leading to the conclusion that no regular
perturbations with $\ell\geq 2$ are allowed.

Figure 3: The behavior of the numerical solution for $\mathsf{H}_{2}$ near the
universal horizon. The boundary condition at large $r$ is given by
$\mathcal{B}_{0}=1$ and $\mathcal{B}_{2}=0$. The modified gravity parameter is
given by $\beta=1/5$.

Figure 4: The behavior of the numerical solution for $\mathsf{H}_{2}$ near the
event horizon (blue solid line). The boundary condition at large $r$ is given
by $\mathcal{B}_{0}=0$ and $\mathcal{B}_{2}=1$. The modified gravity parameter
is given by $\beta=1/5$. The analytic result is also shown as the orange
dotted line, with $\bar{c}=1$.

##  V Conclusions

We have considered a spatially covariant theory of gravity having just two
tensorial degrees of freedom and a non-propagating (instantaneous) scalar
mode. There is no propagating scalar mode, and therefore the number of
propagating degrees of freedom is the same as in general relativity (GR). In a
particular subset of such theories [20], the standard Newtonian behavior of
gravity is reproduced and the propagation of gravitational waves in a
cosmological background obeys the same equation as in GR [21]. Moreover, in
the same subset of theories of [20], the Schwarzschild solution foliated by
the maximal slices is a solution to the field equations, as in the case of
Einstein-Aether theory [38]. The solution forms a universal horizon, which is
the causal boundary for the scalar mode with infinite propagation speed [21].
In this paper, we have studied linear perturbations of this black hole
solution.

First, we have studied the odd-parity perturbations. The odd-parity sector
contains only one of the two tensorial modes but is devoid of contributions
from the instantaneous scalar mode. We have presented the quadratic Lagrangian
for the odd-parity perturbations, which turned out to be identical to that in
GR after an appropriate identification of the variables. Since the tensorial
modes propagate at the speed of light, the causal boundary of the odd-parity
perturbations is given by the usual even horizon rather than the universal
horizon inside of it.

Next, we have considered the even-parity sector of black hole perturbations,
in which the instantaneous scalar mode and one of the two tensorial modes are
mixed. We have derived quadratic actions for monopole, dipole, and higher
multipole modes ($\ell\geq 2$), and studied them separately. The tensorial
mode does not contribute to the monopole and dipole perturbations. As a
result, we have found no radiative behavior in these perturbations. We have
solved the set of field equations for the monopole perturbations and found no
solution that is regular both at infinity and at the inner boundary, i.e. the
universal horizon, except for the trivial one corresponding to a shift of the
mass parameter. For the dipole perturbations, we have derived a single master
equation which is a fourth-order differential equation with respect to the
radial coordinate. Again, we have found no solution that is regular both at
infinity and at the universal horizon.

The situation gets more complicated for the even-parity perturbations with
$\ell\geq 2$. We have derived the general quadratic action. Although there is
only one propagating degree of freedom in the even-parity sector with
$\ell\geq 2$, we have not been able to reduce the system to a single master
equation because of the complexity of the equations stemming from the mixing
with the instantaneous scalar mode. We therefore focused on stationary
perturbations and investigated their properties. We have found that the
stationary perturbations with $\ell\geq 2$ at large $r$ are characterized by
two integration constants, $\mathcal{B}_{2}$ and $\mathcal{B}_{0}$, governing
respectively the contributions from the tensorial degree of freedom and the
instantaneous scalar mode. We carefully identified the locations of the
appropriate inner boundaries and showed that the perturbations diverge at the
inner boundaries in any case unless $\mathcal{B}_{2}=\mathcal{B}_{0}=0$.

To conclude, we have developed black hole perturbation theory for the
Schwarzschild solution foliated by the maximal slices in spatially covariant
gravity with just two tensorial degrees of freedom and established its
perturbative uniqueness. As a future direction, it would be interesting to
study the wave-like behavior of the even-parity perturbations with $\ell\geq
2$ to see the impact of mixing with the instantaneous scalar mode on
gravitational waves, which is technically more difficult and is left for
further study.

###### Acknowledgements.

The work of JS was supported by the Rikkyo University Special Fund for
Research. The work of TK was supported by JSPS KAKENHI Grant No. JP20K03936
and MEXT-JSPS Grant-in-Aid for Transformative Research Areas (A) “Extreme
Universe”, No. JP21H05182 and No. JP21H05189.

## References

  * [1] D. Lovelock, _The Einstein tensor and its generalizations_ , _J. Math. Phys._ 12 (1971) 498.
  * [2] D. Lovelock, _The four-dimensionality of space and the einstein tensor_ , _J. Math. Phys._ 13 (1972) 874.
  * [3] G. W. Horndeski, _Second-order scalar-tensor field equations in a four-dimensional space_ , _Int. J. Theor. Phys._ 10 (1974) 363.
  * [4] C. Deffayet, X. Gao, D. A. Steer and G. Zahariade, _From k-essence to generalised Galileons_ , _Phys. Rev. D_ 84 (2011) 064039 [1103.3260].
  * [5] T. Kobayashi, M. Yamaguchi and J. Yokoyama, _Generalized G-inflation: Inflation with the most general second-order field equations_ , _Prog. Theor. Phys._ 126 (2011) 511 [1105.5723].
  * [6] M. Zumalacárregui and J. García-Bellido, _Transforming gravity: from derivative couplings to matter to second-order scalar-tensor theories beyond the Horndeski Lagrangian_ , _Phys. Rev. D_ 89 (2014) 064046 [1308.4685].
  * [7] J. Gleyzes, D. Langlois, F. Piazza and F. Vernizzi, _Healthy theories beyond Horndeski_ , _Phys. Rev. Lett._ 114 (2015) 211101 [1404.6495].
  * [8] D. Langlois and K. Noui, _Degenerate higher derivative theories beyond Horndeski: evading the Ostrogradski instability_ , _JCAP_ 02 (2016) 034 [1510.06930].
  * [9] M. Crisostomi, K. Koyama and G. Tasinato, _Extended Scalar-Tensor Theories of Gravity_ , _JCAP_ 04 (2016) 044 [1602.03119].
  * [10] J. Ben Achour, M. Crisostomi, K. Koyama, D. Langlois, K. Noui and G. Tasinato, _Degenerate higher order scalar-tensor theories beyond Horndeski up to cubic order_ , _JHEP_ 12 (2016) 100 [1608.08135].
  * [11] J. Ben Achour, D. Langlois and K. Noui, _Degenerate higher order scalar-tensor theories beyond Horndeski and disformal transformations_ , _Phys. Rev. D_ 93 (2016) 124005 [1602.08398].
  * [12] T. Kobayashi, _Horndeski theory and beyond: a review_ , _Rept. Prog. Phys._ 82 (2019) 086901 [1901.07183].
  * [13] X. Gao, _Unifying framework for scalar-tensor theories of gravity_ , _Phys. Rev. D_ 90 (2014) 081501 [1406.0822].
  * [14] X. Gao, _Hamiltonian analysis of spatially covariant gravity_ , _Phys. Rev. D_ 90 (2014) 104033 [1409.6708].
  * [15] X. Gao and Y.-M. Hu, _Higher derivative scalar-tensor theory and spatially covariant gravity: the correspondence_ , _Phys. Rev. D_ 102 (2020) 084006 [2004.07752].
  * [16] N. Afshordi, D. J. H. Chung and G. Geshnizjani, _Cuscuton: A Causal Field Theory with an Infinite Speed of Sound_ , _Phys. Rev. D_ 75 (2007) 083513 [hep-th/0609150].
  * [17] N. Afshordi, D. J. H. Chung, M. Doran and G. Geshnizjani, _Cuscuton Cosmology: Dark Energy meets Modified Gravity_ , _Phys. Rev. D_ 75 (2007) 123509 [astro-ph/0702002].
  * [18] A. Iyonaga, K. Takahashi and T. Kobayashi, _Extended Cuscuton: Formulation_ , _JCAP_ 12 (2018) 002 [1809.10935].
  * [19] C. Lin and S. Mukohyama, _A Class of Minimally Modified Gravity Theories_ , _JCAP_ 10 (2017) 033 [1708.03757].
  * [20] X. Gao and Z.-B. Yao, _Spatially covariant gravity theories with two tensorial degrees of freedom: the formalism_ , _Phys. Rev. D_ 101 (2020) 064018 [1910.13995].
  * [21] A. Iyonaga and T. Kobayashi, _Distinguishing modified gravity with just two tensorial degrees of freedom from general relativity: Black holes, cosmology, and matter coupling_ , _Phys. Rev. D_ 104 (2021) 124020 [2109.10615].
  * [22] T. Hiramatsu and T. Kobayashi, _Testing gravity with the cosmic microwave background: constraints on modified gravity with two tensorial degrees of freedom_ , _JCAP_ 07 (2022) 040 [2205.04688].
  * [23] E. Barausse, T. Jacobson and T. P. Sotiriou, _Black holes in Einstein-aether and Horava-Lifshitz gravity_ , _Phys. Rev. D_ 83 (2011) 124043 [1104.2889].
  * [24] D. Blas and S. Sibiryakov, _Horava gravity versus thermodynamics: The Black hole case_ , _Phys. Rev. D_ 84 (2011) 124043 [1110.2195].
  * [25] T. Kobayashi, H. Motohashi and T. Suyama, _Black hole perturbation in the most general scalar-tensor theory with second-order field equations I: the odd-parity sector_ , _Phys. Rev. D_ 85 (2012) 084025 [1202.4893].
  * [26] T. Kobayashi, H. Motohashi and T. Suyama, _Black hole perturbation in the most general scalar-tensor theory with second-order field equations II: the even-parity sector_ , _Phys. Rev. D_ 89 (2014) 084042 [1402.6740].
  * [27] K. Takahashi, H. Motohashi and M. Minamitsuji, _Linear stability analysis of hairy black holes in quadratic degenerate higher-order scalar-tensor theories: Odd-parity perturbations_ , _Phys. Rev. D_ 100 (2019) 024041 [1904.03554].
  * [28] K. Tomikawa and T. Kobayashi, _Perturbations and quasinormal modes of black holes with time-dependent scalar hair in shift-symmetric scalar-tensor theories_ , _Phys. Rev. D_ 103 (2021) 084041 [2101.03790].
  * [29] K. Takahashi and H. Motohashi, _Black hole perturbations in DHOST theories: master variables, gradient instability, and strong coupling_ , _JCAP_ 08 (2021) 013 [2106.07128].
  * [30] C. de Rham and J. Zhang, _Perturbations of stealth black holes in degenerate higher-order scalar-tensor theories_ , _Phys. Rev. D_ 100 (2019) 124023 [1907.00699].
  * [31] J. Khoury, M. Trodden and S. S. C. Wong, _Existence and instability of hairy black holes in shift-symmetric Horndeski theories_ , _JCAP_ 11 (2020) 044 [2007.01320].
  * [32] D. Langlois, K. Noui and H. Roussille, _Black hole perturbations in modified gravity_ , _Phys. Rev. D_ 104 (2021) 124044 [2103.14750].
  * [33] D. Langlois, K. Noui and H. Roussille, _On the effective metric of axial black hole perturbations in DHOST gravity_ , _JCAP_ 08 (2022) 040 [2205.07746].
  * [34] K. Nakashi, M. Kimura, H. Motohashi and K. Takahashi, _Black hole perturbations in higher-order scalar–tensor theories: initial value problem and dynamical stability_ , _Class. Quant. Grav._ 39 (2022) 175003 [2204.05054].
  * [35] K. Lin, S. Mukohyama, A. Wang and T. Zhu, _No static black hole hairs in gravitational theories with broken Lorentz invariance_ , _Phys. Rev. D_ 95 (2017) 124053 [1704.02990].
  * [36] A. De Felice, D. Langlois, S. Mukohyama, K. Noui and A. Wang, _Generalized instantaneous modes in higher-order scalar-tensor theories_ , _Phys. Rev. D_ 98 (2018) 084024 [1803.06241].
  * [37] C. Zhang, A. Wang and T. Zhu, _Universal horizons and black hole spectroscopy in gravitational theories with broken Lorentz symmetry_ , 2209.04735.
  * [38] P. Berglund, J. Bhattacharyya and D. Mattingly, _Mechanics of universal horizons_ , _Phys. Rev. D_ 85 (2012) 124019 [1202.4497].


# Fast robust peg-in-hole insertion
with continuous visual servoing

Rasmus Laurvig Haugaard
SDU Robotics
University of Southern Denmark Denmark
<EMAIL_ADDRESS>
&Jeppe Langaa
SDU Robotics
University of Southern Denmark Denmark
<EMAIL_ADDRESS>
&Christoffer Sloth
SDU Robotics
University of Southern Denmark Denmark
<EMAIL_ADDRESS>
&Anders Glent Buch
SDU Robotics
University of Southern Denmark Denmark
<EMAIL_ADDRESS>

###### Abstract

This paper demonstrates a visual servoing method which is robust towards
uncertainties related to system calibration and grasping, while significantly
reducing the peg-in-hole time compared to classical methods and recent
attempts based on deep learning. The proposed visual servoing method is based
on peg and hole point estimates from a deep neural network in a multi-cam
setup, where the model is trained on purely synthetic data. Empirical results
show that the learnt model generalizes to the real world, allowing for higher
success rates and lower cycle times than existing approaches.

> Keywords: Peg-In-Hole, Visual Servoing

##  1 Introduction

Uncertainties in robotic systems accumulate from various system-parts,
including robot kinematics, tools, object grasping, camera calibrations, etc.
In high-precision tasks, some uncertainties between target objects may have to
be reduced. One approach to reduce the uncertainty is better global
calibration, but that is likely to be in conflict with other desired system
features, like passive grippers, passive compliance, or modular tool systems
which enable more agile systems, but also introduce further uncertainties.
Another approach to reduce uncertainties which remedies the above limitations
is to handle the uncertainty live and locally, e.g. with force and/or visual
feedback, depending on the task.

While many high-precision tasks exist, we focus on a common industrial task,
inserting pegs into holes (peg-in-hole). Compliant insertion based on force
feedback is effective, when the uncertainty is low enough to enable rich peg-
hole interaction. If the uncertainty is slightly larger, classical methods
like random search and spiral search can effectively reduce the uncertainty to
enable compliant insertion. When the uncertainty becomes larger, the classical
search methods become slow, as they blindly search the uncertainty region.
More importantly, they assume that the surface is flat around the hole within
the uncertainty region, which is only a good assumption if the region is
small, and are likely to fail otherwise. In addition, these methods can be
harmful for the surface quality of both the peg and the hole-object, and they
might result in dropping the peg, if the grasp-force is inadequate, which is
not unlikely with passive grippers.

Visual feedback has the potential to effectively reduce larger uncertainties.
Only detecting one of the target objects though, only reduces some of the
uncertainties. E.g. in the peg-in-hole case, mounting in-hand cameras on the
peg-robot to detect the hole in the robot’s tool-frame will ideally eliminate
uncertainties from robot kinematics and the hole position relative to the peg-
robot base, but it will not reduce the uncertainties from tooling and grasping
but instead introduce new uncertainties from camera calibration. Classical
vision methods tend to either be sensitive to lighting and background clutter
or require visual markers.

We propose learnt visual servoing for the initial peg-hole alignment which is
robust to system calibration, grasping uncertainties, surface geometry and
lighting conditions, while being significantly faster than spiral search,
random search and recent methods based on deep learning. The servoing is based
on peg and hole point estimates from a deep neural network in a multi-cam
setup, where the model is trained on purely synthetic data.

Figure 1 shows our setup consisting of two cameras mounted on one robot
(camera-robot), a peg mounted on another robot (peg-robot), and four holes.
Our method also allows the cameras to be attached to the peg-robot or be fixed
to the table, but having the cameras on a separate robot provides high
flexibility, enabling visual servoing for multiple tools on multiple robots.
To test generalization, the setup includes two peg types and four hole types,
covering multiple materials, scales and tolerances, as well as different peg
and hole surface geometries.

Figure 1:  Setup consisting of two UR5 robots where two Intel Realsense D435
are mounted on one robot (camera-robot) and a peg is mounted on the other
robot (peg-robot). Four hole-types, metal, plastic, wide and cap are mounted
to the table. For the cap hole, an M4 bolt is used as the peg. Note that we
only use the color sensors, and thus the size of the camera tool can be
greatly reduced.

##  2 Related Work

A vast amount of robotic applications rely mainly on predefined and static
movements, which do not account for dynamical changes of the environment or
task. In contrast, visual servoing incorporates a visual feedback into the
control loop, which enables the robot to deal with large deviations. A variety
of different visual servoing methods exists including image-based visual
servoing [1], pose-based visual servoing [2] and hybrid approaches [3] that
use parts of both. Visual servoing opens up for a lot of possibilities by
providing feedback without tactility, e.g. tracking and grasping objects with
a robot [4, 5]. Classical visual servoing requires that the robot is
accurately calibrated, while more recent methods have tried to address this
problem with uncalibrated visual servoing methods [6, 7, 8], but these methods
needs more time to converge, due to the omitted models of the kinematics and
camera calibrations. Instead, we propose to make an initial calibration of the
system, but to be robust towards uncertainties in the calibration, similar to
[9]. This ensures that small changes to the setup do not affect the
performance of the method, while still exploiting the advantages of a
calibrated system.

In deep learning, domain randomization has enabled training models on
synthetic data and successfully transferring them to the real world without
further training in a variety of tasks, including object detection [10], using
a robotic hand to manipulate a cube with reinforcement learning [11], car
detection and pose estimation [12] and affordance learning [13]. [10, 12, 13]
use synthetic 3D distractor objects to improve their model’s generalization.
Instead, we overlay natural images on our synthetic images in order to improve
generalization.

More closely related to our proposed method is [14], which applies visual
servoing for peg-in-hole alignment assuming a full pose- estimation from
classical methods, where peg and hole are clearly marked. Although the full
pose enables the method to handle the rotational alignment of the peg, the
requirement of peg and hole markers make such methods undesirable in most use
cases. Our method does not require visual markers.

Recently, [15] proposed image-based visual servoing for peg and hole alignment
prior to insertion, which is similar to our approach. They fix two in-hand
cameras to the peg-robot in a top-view configuration and train a deep neural
network to classify if the robot should move in one of four directions, or
whether they should start the insertion. Their visual servoing is step-wise,
not continuous. While they outperform spiral search in their results, forming
the problem as classification and using step-wise visual servoing makes the
method too slow for practical use with insertion durations in the range of
20-70 seconds. Also, since they only detect the hole, their method does not
handle grasping and tool uncertainties. In contrast, we estimate the
positional error of the peg-robot based on both peg and hole point estimates,
which handles grasping and tool uncertainties, and reduce the error
continuously, enabling significantly faster insertions.

##  3 Methods

This section presents an overview of the alignment methods examined in this
work, including the proposed visual servoing. Our visual servoing alignment is
divided into two tasks: First, the hole and peg center points are estimated in
images using a deep neural network. Secondly, based on the point estimates,
the peg is aligned to the hole in the plane that is perpendicular to the
insertion direction. After alignment, the peg is inserted with compliance
using force-feedback.

###  3.1 Overview of peg-in-hole alignment

We examine three methods for peg-hole positional alignment: random search,
spiral search and our proposed visual servoing alignment. The methods are
illustrated in Figure 2. Note that all three methods only align the position,
not the orientation. In random search, a point within the uncertainty boundary
is sampled and an insertion is attempted at that point. This is repeated until
the peg reaches below the hole-surface or a time limit is exceeded. In spiral
search, the peg moves outwards in a spiral while pressing against the hole-
surface. The peg is partially inserted in the hole when it moves over the hole
within a success region based on the peg- and hole tolerances. Spiral search
will continue until a force limit is reached, indicating that the peg is in
contact with the hole, or until the peg exceeds the uncertainty boundary. In
our proposed visual servoing, both peg and hole center points are continuously
estimated in two or more cameras. An error vector is estimated based on the
estimated points, and the error is attempted minimized in a servoing loop.
Visual servoing will continue either until convergence, until a time limit is
reached or until the peg exceeds the uncertainty boundary. At convergence, the
peg is moved downward along the insertion direction until force-feedback
indicates contact. To further increase robustness during alignment, visual
servoing can be followed by spiral search or random search with a smaller
uncertainty region where appropriate.

After successful positional alignment, there may still be small position
uncertainties as well as orientation uncertainties. To remedy this, the peg is
inserted with compliance using force-feedback.

random searchvisual servoingspiral searchhole outlinepeg start outlinepeg
searchuncertainty boundarypeg start positionpeg success region Figure 2:
Illustration of the three examined methods for peg-in-hole alignment. The
paths between the points in random search are left out for clarity. Note that
the peg is in contact with the hole-surface at the search points in random
search and throughout spiral search while the peg moves above the surface in
visual servoing.

###  3.2 Point estimation

We make point estimates based on heatmaps from deep neural networks, similar
to [16, 17] which estimate human pose keypoints, but here adapted to pegs and
holes. As shown by [18], outputting heatmaps allows for better spatial
generalization than directly regressing point coordinates from a global,
latent representation. Let $p=(x,y)^{T}$ be the pixel coordinates of a point
in an image, $I$. Then a desired heatmap, $\Phi$, from a desired point,
$p^{*}$, can be defined using a Gaussian kernel:

$\Phi_{p^{*}}=\exp\left(-\frac{||p-p^{*}||^{2}}{2\sigma^{2}}\right),$ (1)

where $\sigma$ is a hyper parameter controlling the size of the active area of
the heatmap. We set $\sigma=3\text{px}$ and have separate heatmaps for the peg
and hole points. Figure 4.f shows an example of target heatmaps.

A deep neural network, $f$, with trainable parameters, $\phi$, estimates the
heatmap, $\hat{\Phi}=f_{\phi}(I)$, with the same resolution as the input
image. We choose a U-Net [19] architecture with an ImageNet pre-trained
ResNet18 [20] backbone for simplicity and to leverage transfer learning while
keeping the inference time low. As will be apparent from the results in
Section 4, this architecture is sufficient to model the problem at hand.

At inference, the point with the largest corresponding value in the estimated
heatmap is chosen as the point estimate, $\hat{p}$:

$\hat{p}=\operatorname*{\arg\\!\max}_{p}\hat{\Phi}_{p}$ (2)

####  3.2.1 Data generation

Our proposed model used in our real-world experiments is trained purely on
synthetic data with no manual annotations, inspired by [10, 11, 12, 13]. The
virtual scene consists of a camera, a cylinder and a surface with a hole in
it. For each image, peg- and camera transforms are sampled, a high dynamic
range image (HDRI) environment map is sampled, and a procedural material is
sampled. We do not apply randomization to the geometry.

When synthesizing an image, the distance from the camera to the hole is
sampled uniformly between 12 cm and 15 cm. The elevation (the angle between
the hole plane and the optical axis) is sampled between 35° and 45°. The roll
angle (the rotation around the optical axis) is sampled between -5°and 5°. The
peg position is sampled from a disc centered around the hole with a radius of
15 mm and a height relative to the hole between 5 mm and 15 mm. The peg
orientation error is sampled as an axis angle vector, where the axis is
sampled from the surface of a unit sphere, and the angle is sampled between 0°
and 5°.

We use an in-house developed framework [21] for synthetic image generation in
Blender111https://www.blender.org/ with procedural materials and HDRI
environment maps from HDRI Haven222https://hdrihaven.com/. Environment maps
provide a large variety of natural lighting and reflections. We hypothesize
that realistic, diverse lighting and reflections enable the model to learn
features that generalize well to the real world. The sampled, procedural
material is applied to both the hole-surface and the peg to discourage the
model from learning to distinguish between different samples of the procedural
material. We render 1000 images in about 10 minutes on an NVIDIA GeForce RTX
2080. Examples of the synthetic renders are shown in Figure 3.

Figure 3:  Examples of synthetic renders with procedural materials and HDRI
environment maps.

If the images shown to the model during training only contain the peg and the
hole-surface, the model will likely be sensitive to background clutter. One
approach to remedy this is to use natural images as backgrounds as in [22],
but this encourages the model to discriminate between the natural- and
synthetic domain. [22] alleviates the problem by using somewhat realistic
models and freezing a pre-trained backbone. [10, 12, 13] overcome the problem
by using 3D distractor objects, such that both falses and positives are from
the synthetic domain. Since it is not obvious how to obtain a diverse set of
3D distractor objects, we instead introduce distractors from natural images.
We sample a natural image as an overlay image and another natural image,
converted to grayscale, as the alpha channel of the overlay image. We sample
the natural images from the MS COCO dataset [23]. The augmentation is
illustrated in Figure 4(a-d). Since the overlay from the natural domain is
imposed on the whole image, the model will not be able to simply distinguish
between the synthetic- and natural domain. Also, natural overlays on peg and
hole will effectively serve as extra randomization on peg and hole appearance,
which can be hard to simulate. Since compositing is computationally efficient,
the augmentation is done live during training with different overlays at each
epoch, effectively reducing the amount of required renders.

Figure 4:  Example of a synthetic training image. a) Synthetic render. b)
Overlay. c) Mask. d) Composite image using a-c. e) Composite image with random
crop, random horizontal flip, and random blur. f) Target heatmaps (hole: blue
channel, peg: red channel).  Figure 5: Examples from the three training
datasets, from top- to bottom row: synth, plastic, metal.

To evaluate the model trained on synthetic images, we obtain natural datasets
for the metal and plastic holes. To obtain the datasets, a rigid camera mount
is attached to the camera-robot, the whole setup is calibrated, and robot
movements between consecutive sampled peg and camera transforms are ensured to
be safe. The peg and camera transforms are sampled like for the synthetic
images. We do not attempt to add any variations to the natural environment,
e.g. turning on and off lights or capturing images from multiple rooms, since
it is impractical.

####  3.2.2 Training details

All images are resized to 224x224 before they are handed to the model. The
loss is defined as the mean squared error between the desired and output
heatmaps. Each of our datasets consists of 1000 images, where 100 of them are
reserved for validation. We also trained on 10x more images on the synthetic
domain while keeping the total amount of parameter updates, but with little to
no effect. We conclude that the data augmentation provides enough variation
within the domain. We train our models with the one-cycle [24] policy for 30
epochs with a maximum learning rate of 0.001 and weight decay set to 0.0001.
During training, both synthetic and natural images are augmented with random
crop, random horizontal flip and a 50% chance of box blur with a kernel size
of either 3 or 5 px. The augmentation is illustrated for a synthetic image in
Figure 4. One training takes about 8 minutes on an NVIDIA GeForce RTX 2080,
and we obtain single-image inference at 150 Hz.

###  3.3 Visual servoing

Our method requires regions of interest (ROIs) in the images, providing crops
as seen in Figure 5, as well as approximated depths (distance from cameras to
peg and hole). ROIs and depths can either be calculated based on an
approximate hole frame or by manually marking the ROIs in the full images once
for a given setup. Our visual servoing alignment method is outlined in
Algorithm 1.

Input: Insertion direction, $l$. Estimated depths, $z_{i}$, in $n$ cameras.
Target convergence threshold, $\phi_{t}$. Filter parameters $\alpha_{\tau}$,
$\alpha_{\gamma}$ and $\alpha_{\phi}$.

$\tau\leftarrow$ initial peg-robot position

// filtered target position

$\phi\leftarrow 10\phi_{t}$

// filtered error magnitude

while  _$\phi >\phi_{t}$_ do

$q\leftarrow$get position of peg-robot

for  _$i=1,...,n$_ do

$I_{i}\leftarrow$ get image from camera $i$

$p_{i}^{*},h_{i}^{*}\leftarrow$ estimate peg and hole points in $I_{i}$

$p_{i},h_{i}\leftarrow$ estimate 3D points based on $p_{i}^{*},h_{i}^{*}$,
$z_{i}$ and the system state- and calibration

$c_{i}\leftarrow$ get position of camera $i$

$v_{i}\leftarrow(p_{i}+h_{i})/2-c_{i}$

// view direction

$u_{i}\leftarrow\dfrac{v_{i}\times l}{||v_{i}\times l||}$

// direction of error from this view

$b_{i}\leftarrow u_{i}\cdot(h_{i}-p_{i})$

// magnitude of error from this view

end for

$(A,b)\leftarrow\left(\begin{bmatrix}u_{1x}&u_{1y}&u_{1z}\\\
u_{2x}&u_{2y}&u_{2z}\\\ &\vdots\end{bmatrix},\begin{bmatrix}b_{1}\\\ b_{2}\\\
\vdots\end{bmatrix}\right)$

$\hat{e}\leftarrow$ solve $(A,b)$ by least squares

// estimated error vector

$t\leftarrow q+\hat{e}$

// unfiltered target position

$\tau\leftarrow\alpha_{\tau}\tau+(1-\alpha_{\tau})t$

// update filtered target position

set peg-robot reference to $\tau$

// move robot towards $\tau$

$\gamma\leftarrow\alpha_{\gamma}\gamma+(1-\alpha_{\gamma})\hat{e}$ if $\gamma$
is defined, else $\hat{e}$

// update filtered error vector

$\phi\leftarrow\alpha_{\phi}\phi+(1-\alpha_{\phi})||\gamma||$

// update filtered error magnitude

end while

Algorithm 1 Visual servo positional alignment

While performing the alignment with visual servoing, the robot moves in a
plane perpendicular to a provided insertion direction, $l$. Note that the
insertion direction is the known direction in which the peg-robot will move
when the positional alignment along the plane is done. It does thus not need
to be exactly parallel to the axis of the peg or the hole. Because of depth-
ambiguity, points from the images of the $i$’th camera provide information
only along the direction, $u_{i}$, in the movement plane that is perpendicular
to the view direction, $v_{i}$. The point estimates from the image are
projected to 3D based on the system state, system calibration and constant
depth estimates, $z_{i}$, used for both the peg- and hole point,
$p_{i},h_{i}$. Peg and hole estimates from multiple cameras then allow us to
estimate the peg-robot error vector $\hat{e}$.

A target position, $t$, for the peg-robot is found by adding the estimated
error vector, $\hat{e}$, to the current peg-robot position, similar to visual
servoing algorithms that use the estimated image Jacobian for controlling the
robot. The target position and the error magnitude are filtered as shown in
Algorithm 1 with coefficients
$\alpha_{\tau}=\alpha_{\gamma}=\alpha_{\phi}=0.9$. Finally, we use the built-
in servoing api from the UR5 robot to move towards the target position, and
terminate the visual servoing when the filtered error magnitude is smaller
than the threshold $\phi_{t}$, which is set to be 1/20 of the peg-hole
diameter.

##  4 Experiments

In this section, we empirically validate different aspects of our visual
servoing system. We start with an accuracy test of the visual part, the point
estimation, and then we go on to test the full system, including both
alignment and insertion.

###  4.1 Point estimation accuracy

Three models are trained on images of plastic, images of metal (see Figure 1),
and the synthetic images, respectively (900 images per dataset). Each model is
evaluated on all validation datasets (100 images per dataset). The cross-
performance is illustrated in Figure 6, showing detection success rates versus
increasing distance tolerances. The models trained on the natural domains fail
on the synthetic domain, and more importantly, the model trained on the metal
domain fails to detect the plastic hole consistently. In contrast, the models
trained on the synthetic domain perform consistently well in all three
domains, indicating that the synthetic images force the models to learn more
general features.

Figure 6:  Point detection success rates at different pixel thresholds. 10
models are trained for each of the three training datasets, and each model is
evaluated on all of the validation datasets. Both mean performance per
training dataset and the individual model performances are plotted. The
individual performances are plotted in a lighter color for clarity. Note that
the peg- and hole diameter is approximately 50 px in the images.

###  4.2 Peg-in-hole performance

Using models trained on all three domains, we compare our visual servoing to
random search and spiral search. The results are shown in Table 1 See
supplementary material for a video of the methods.

The initial peg positions are sampled from a disc centered around the hole
with a diameter of three times the peg diameter and a height relative to the
hole between 5 mm and 15 mm (3 mm and 5 mm for cap). The peg-hole orientation
error is sampled up to 2°. To show robustness towards calibration, camera
extrinsic errors are sampled with up to 2° in orientation and 10 mm in
position. Note that these are large calibration errors.

The peg used for bearing, plastic and wide is the 10 mm shaft with an h7
tolerance shown in Figure 1. The bearing hole is 10 mm with an H7 tolerance
while the plastic and wide holes are 3D-printed with a measured diameter of
10.6 mm and 10.4 mm, respectively. The cap hole is 4.4 mm and the peg is an M4
bolt with a measured diameter of 3.9 mm.

Random search relies on detecting the peg reaching below the hole-surface. For
this reason we do not apply randomization in peg height. It also means that
random search requires a rather large tolerance compared to the uncertainty,
which is why random search is not considered for metal with H7/h7 tolerance.
Visual servoing is followed by spiral search as mentioned in Section 3.1 for
metal, plastic and wide, but is avoided for cap to avoid damaging the hole
surface with the M4 bolt. For the same reason, we do not attempt force-based
insertion for cap.

In the plastic* experiment, a light source was placed pointing into the
cameras, causing flares and hard reflections on the hole-surface (see Figure
7). All attempts were successful. The peg for the cap hole is threaded, and is
thus far from any of the three examined domains. The models trained on the
natural domains completely fail to generalize, while the model trained on the
synthetic domain succeeds in all the attempts.

The attempts with our proposed visual servoing are consistently both faster
and more robust than random search and spiral search. Furthermore, the model
trained on the synthetic domain generalize to all the examined cases, in
contrast to the models trained on natural domains. [15] also performs peg-in-
hole visual servoing based on deep learning on synthetic images. They use a 10
mm peg and a 10.4 mm hole, similar to wide, and a similar but only positional
error in the start positions. They achieve insertion times between 20 s and 70
s, while our average insertion time is 2.3 s.

| classical search methods | visual servoing, trained on: |
---|---|---|---
| random | spiral | plastic | metal | synth (proposed) | optimal
metal | - | 10% (7.1 s) |  99% (4.1 s) |  100% (4.0 s) |  100% (3.6 s) | 3.2 s
plastic | 0% (-) | 27% (10 s) |  100% (2.4s) |  100% (2.7 s) |  100% (2.2 s) | 1.0 s
plastic* | - | - | 93% (2.6 s) | 11% (3.7 s) |  100% (2.3 s) | -
wide | 4% (15 s) |  100% (36 s) |  99% (2.4 s) |  100% (2.2 s) |  100% (2.3 s) | 1.0 s
cap | 34% (16 s) | - | 1% (1.7 s) | 0% (-) |  100% (1.6 s) | 0.4 s
Table 1:  Success rates and mean time for successful attempts. Timings
indicate the total time from the initial position to reaching an insertion
depth of 10 mm (5 mm for cap). Each entry is based on 100 runs. Highest
success rates are marked in bold. If an entry is not significantly lower than
the highest success rate with a p-value of 5%, it is also marked in bold.
Optimal: peg is aligned to the hole based on calibration (no search).

Figure 7:  Example images from the cameras during visual servoing of the peg-
in-hole cases. From left to right: metal, plastic, plastic*, wide and cap.

##  5 Conclusion

We studied the use of deep learning based point estimation and continuous
visual servoing for peg-in-hole positional alignment. We demonstrated that our
method can significantly increase speed and robustness compared to classical
methods like random search and spiral search. We also showed that our method
is significantly faster than a previous attempt based on deep learning. Since
visual servoing does not require peg-hole contact, our method is robust to
hole surface geometry and gentle to the peg and hole surfaces. We trained
point estimation models on both natural and synthetic domains. We demonstrated
that synthetic domain randomization and using distractors from natural images
introduces enough variation to enable the model to generalize to all the
examined peg-in-hole cases. While compliant force insertion is likely not
enough to account for larger angle errors in tight tolerances, integrating
axis alignment in visual servoing could remedy this problem and is an area for
future work. Finally, our method is focused on peg-in-hole tasks but we
hypothesize that deep learning based visual servoing can be applied
successfully in many assembly sub-tasks.

#### Acknowledgments

The authors gratefully acknowledge the economic support from Innovation Fund
Denmark through the project MADE FAST.

## References

  * Weiss et al. [1987] L. Weiss, A. Sanderson, and C. Neuman.  Dynamic sensor-based control of robots with visual feedback.  _IEEE Journal on Robotics and Automation_ , 3(5):404–417, 1987.
  * Thuilot et al. [2002] B. Thuilot, P. Martinet, L. Cordesses, and J. Gallice.  Position based visual servoing: keeping the object in the field of vision.  In _Proceedings 2002 IEEE International Conference on Robotics and Automation (Cat. No.02CH37292)_ , volume 2, pages 1624–1629 vol.2, 2002.
  * Malis et al. [1999] E. Malis, F. Chaumette, and S. Boudet.  2 1/2 d visual servoing.  _IEEE Transactions on Robotics and Automation_ , 15(2):238–250, 1999.
  * Allen et al. [1993] P. K. Allen, A. Timcenko, B. Yoshimi, and P. Michelman.  Automated tracking and grasping of a moving object with a robotic hand-eye system.  _IEEE Transactions on Robotics and Automation_ , 9(2):152–165, 1993.
  * Wu et al. [2017] H. Wu, T. T. Andersen, N. A. Andersen, and O. Ravn.  Application of visual servoing for grasping and placing operation in slaughterhouse.  In _2017 3rd International Conference on Control, Automation and Robotics (ICCAR)_ , pages 457–462, 2017.
  * Piepmeier et al. [2004] J. A. Piepmeier, G. V. McMurray, and H. Lipkin.  Uncalibrated dynamic visual servoing.  _IEEE Transactions on Robotics and Automation_ , 20(1):143–147, 2004.
  * Qiu et al. [2019] Z. Qiu, S. Hu, and X. Liang.  Model predictive control for constrained image-based visual servoing in uncalibrated environments.  _Asian Journal of Control_ , 21(2):783–799, 2019\.  doi:10.1002/asjc.1756.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/asjc.1756.
  * Chang et al. [2020] Y. Chang, L. Li, Y. Wang, and K. You.  Toward fast convergence and calibration-free visual servoing control: A new image based uncalibrated finite time control scheme.  _IEEE Access_ , 8:88333–88347, 2020.
  * Taylor and Ostrowski [2000] C. J. Taylor and J. P. Ostrowski.  Robust vision-based pose control.  In _Proceedings 2000 ICRA. Millennium Conference. IEEE International Conference on Robotics and Automation. Symposia Proceedings (Cat. No.00CH37065)_ , volume 3, pages 2734–2740 vol.3, 2000.
  * Hinterstoisser et al. [2019] S. Hinterstoisser, O. Pauly, H. Heibel, M. Marek, and M. Bokeloh.  An annotation saved is an annotation earned: Using fully synthetic training for object instance detection.  _CoRR_ , abs/1902.09967, 2019.  URL http://arxiv.org/abs/1902.09967.
  * OpenAI et al. [2018] OpenAI, M. Andrychowicz, B. Baker, M. Chociej, R. Jozefowicz, B. McGrew, J. Pachocki, A. Petron, M. Plappert, G. Powell, A. Ray, J. Schneider, S. Sidor, J. Tobin, P. Welinder, L. Weng, and W. Zaremba.  Learning dexterous in-hand manipulation, 2018.
  * Khirodkar et al. [2018] R. Khirodkar, D. Yoo, and K. M. Kitani.  Domain randomization for scene-specific car detection and pose estimation, 2018.
  * Hämäläinen et al. [2019] A. Hämäläinen, K. Arndt, A. Ghadirzadeh, and V. Kyrki.  Affordance learning for end-to-end visuomotor robot control.  _CoRR_ , abs/1903.04053, 2019.  URL http://arxiv.org/abs/1903.04053.
  * Huang et al. [2013] S. Huang, K. Murakami, Y. Yamakawa, T. Senoo, and M. Ishikawa.  Fast peg-and-hole alignment using visual compliance.  In _2013 IEEE/RSJ international conference on intelligent robots and systems_ , pages 286–292. IEEE, 2013.
  * Triyonoputro et al. [2019] J. C. Triyonoputro, W. Wan, and K. Harada.  Quickly inserting pegs into uncertain holes using multi-view images and deep network trained on synthetic data.  _CoRR_ , abs/1902.09157, 2019.  URL http://arxiv.org/abs/1902.09157.
  * Xiao et al. [2018] B. Xiao, H. Wu, and Y. Wei.  Simple baselines for human pose estimation and tracking.  _CoRR_ , abs/1804.06208, 2018.  URL http://arxiv.org/abs/1804.06208.
  * Wang et al. [2020] J. Wang, K. Sun, T. Cheng, B. Jiang, C. Deng, Y. Zhao, D. Liu, Y. Mu, M. Tan, X. Wang, W. Liu, and B. Xiao.  Deep high-resolution representation learning for visual recognition.  _IEEE transactions on pattern analysis and machine intelligence_ , 2020.
  * Nibali et al. [2018] A. Nibali, Z. He, S. Morgan, and L. Prendergast.  Numerical coordinate regression with convolutional neural networks.  _CoRR_ , abs/1801.07372, 2018.  URL http://arxiv.org/abs/1801.07372.
  * Ronneberger et al. [2015] O. Ronneberger, P. Fischer, and T. Brox.  U-net: Convolutional networks for biomedical image segmentation.  _CoRR_ , abs/1505.04597, 2015.  URL http://arxiv.org/abs/1505.04597.
  * He et al. [2015] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  _CoRR_ , abs/1512.03385, 2015.  URL http://arxiv.org/abs/1512.03385.
  * Rasmus Laurvig Haugaard [2020] Rasmus Laurvig Haugaard.  synth-ml, 2020.  URL https://gitlab.com/sdurobotics/vision/synth-ml.
  * Hinterstoisser et al. [2017] S. Hinterstoisser, V. Lepetit, P. Wohlhart, and K. Konolige.  On pre-trained image features and synthetic images for deep learning.  _CoRR_ , abs/1710.10710, 2017.  URL http://arxiv.org/abs/1710.10710.
  * Lin et al. [2014] T. Lin, M. Maire, S. J. Belongie, L. D. Bourdev, R. B. Girshick, J. Hays, P. Perona, D. Ramanan, P. Dollár, and C. L. Zitnick.  Microsoft COCO: common objects in context.  _CoRR_ , abs/1405.0312, 2014.  URL http://arxiv.org/abs/1405.0312.
  * Smith [2018] L. N. Smith.  A disciplined approach to neural network hyper-parameters: Part 1 - learning rate, batch size, momentum, and weight decay.  _CoRR_ , abs/1803.09820, 2018.  URL http://arxiv.org/abs/1803.09820.