# An extended two-dimensional vocal tract model for fast acoustic simulation
of single-axis symmetric three-dimensional tubes

###### Abstract

The simulation of two-dimensional (2D) wave propagation is an affordable
computational task and its use can potentially improve time performance in
vocal tracts’ acoustic analysis. Several models have been designed that rely
on 2D wave solvers and include 2D representations of three-dimensional (3D)
vocal tract-like geometries. However, until now, only the acoustics of
straight 3D tubes with circular cross-sections have been successfully
replicated with this approach. Furthermore, the simulation of the resulting 2D
shapes requires extremely high spatio-temporal resolutions, dramatically
reducing the speed boost deriving from the usage of a 2D wave solver. In this
paper, we introduce an in-progress novel vocal tract model that extends the 2D
Finite-Difference Time-Domain wave solver (2.5D FDTD) by adding tube depth,
derived from the area functions, to the acoustic solver. The model combines
the speed of a light 2D numerical scheme with the ability to natively simulate
3D tubes that are symmetric in one dimension, hence relaxing previous
resolution requirements. An implementation of the 2.5D FDTD is presented,
along with evaluation of its performance in the case of static vowel modeling.
The paper discusses the current features and limits of the approach, and the
potential impact on computational acoustics applications.

Index Terms: computational acoustics, vocal tract, FDTD, articulatory speech
synthesis

##  1 Introduction

The shape of the human upper vocal tract is quite intricate and the study of
its acoustics proves a challenging research topic. The cavity that connects
the larynx to the mouth opening forms a curved tube with varying irregular
cross-sectional shape, in turn connected to several side branches, e.g., the
nasal tract and the piriform fossae. Since every detail of this complex
structure has an effect on voice production, the most typical approach to the
analysis of vocal tract acoustics is the employment of three-dimensional (3D)
simulations [1, 2, 3]. Such simulations capitalize on the computation of
pressure wave propagation within precise 3D geometrical reconstructions of
real vocal tracts, and are capable of producing highly accurate results.
However, the computational power required to solve wave equations in three
dimensions is remarkably high, yielding to simulation times that range from
several minutes to more than a day for just a few milliseconds of audio [1,
4].

Several research groups explored two-dimensional (2D) simulations as an
alternative to a precise but costly 3D approach [5, 6, 7, 4]. The reason for
this interest is two-fold. Compared to the 3D case, 2D wave equation solvers
are characterized by much lower computational requirements and can potentially
achieve real-time or quasi real-time performance [8, 9]. Furthermore, this
class of simulations appears to be capable of preserving most of the
geometrical details of the vocal tract that is being analyzed. In 2D a vocal
tract is represented by a flat contour, typically extracted from the mid-
sagittal slice of a 3D tube built from its area function [7, 4]; it can
include the vocal tract’s original curvature as well as 2D equivalents of its
side-branches [6]. In other words, the only missing information regards the
varying cross-sectional shape of the 3D vocal tract, which is replaced by a
set of circumferences with the same area progression.

However, in a practical scenario the advantages of 2D simulations prove quite
limited. As thoroughly discussed by Arnela and Guasch for the case of vowel
synthesis [4], a direct mid-sagittal representation of the tube built from the
vocal tract’s area function produces erroneous formant locations and
bandwidths, i.e., not matching the results obtained from the 3D acoustic
simulation of the same 3D shape. To fix these discrepancies, the extracted 2D
contours have to undergo a non-linear deformation process that leads to a
significant downsizing of the vocal tract’s constrictions, in some cases up to
an order of magnitude [9]. A direct consequence of this modification regards
simulation times. To model such narrow constrictions, 2D simulations need to
run with extremely high spatio-temporal resolutions, in turn increasing the
computational load of the solver and producing waiting times far from real-
time performance. Moreover, the methodology proposed by Arnela and Guasch is
applicable only to the case of straight tubes, and it is still unclear how to
obtain 2D contours that acoustically match curved and/or branching 3D
geometries.

In this paper, we present a novel approach to vocal tract modeling, that stems
from the 2D rationale [9] and improves upon it. It capitalizes on an extended
2D Finite-Difference Time-Domain solver (2.5D FDTD), capable of simulating how
pressure propagates in 3D geometries that are symmetric at least in one
dimension, typically along the $z$ axis (with $x$ and $y$ being the dimensions
of the starting 2D scheme). At the core of this approach there is the
inclusion within the model of extra impedance terms, that derive from the
tube’s depth (i.e., its continuous extension along $z$) and that are sampled
in every point of the scheme; the result is a depth map, that combined with
the 2D mid-sagittal contour of the original 3D geometry allows for a fast
simulation of its acoustics, with computational requirements comparable to the
case of standard 2D numerical schemes. Depth maps can be retrieved from area
functions as well as from full 3D models of real vocal tracts, like Magnetic
Resonance Imaging scans; furthermore, the overall 2.5D representation of the
analyzed geometry leaves its original proportions intact, discarding the need
for the extremely high spatio-temporal resolutions that characterizes 2D
simulations. However, as discussed in the next section, the development of the
2.5D model is still in progress and at its current stage of development it can
be employed only for the acoustic simulation of straight tubes with circular
cross-sectional shape.

##  2 Methods

###  2.1 Acoustics Equation

Let us consider a sound wave propagating in air, and let us call $\bm{\xi}$
the displacement vector of the medium caused by the motion of the wave and $V$
the volume of a generic element of the medium (particle). Now suppose that the
medium is contained between two rigid surfaces, that extend towards infinity
along both the $x$ and the $y$ axis without ever intersecting and that are one
the reflected image of the other (being $x/y$ the plane of symmetry). In the
described scenario, the wave moves in a space that is 3D but constrained on
one axis (i.e., $z$); such a space can be described by the scalar field
$D(x,y)$, which contains the euclidean distance between the points of the two
surfaces at the same $x,y$ coordinates. In other words, $D$ describes the
space in terms of its depth along $z$.

If we choose to focus on a particle of volume $V=dxdyD$, with $dx$ and $dy$
relatively small, $\bm{\xi}$ will only have $x$ and $y$ components due to the
presence of the enclosing surfaces. Consequently, during the passage of the
sound wave, the particle’s volume change can be expressed as:

$\displaystyle dV=(\frac{\partial D\xi_{x}}{\partial x}\ dx)\
dy+(\frac{\partial D\xi_{y}}{\partial y}\ dy)\ dx$ (1)

if $\bm{\xi}$ is quite small (reasonable assumption in a realistic scenario)
and $D$ is slowly varying in space. The two products in parentheses represent
the first order Taylor Series approximations of the increment of $\bm{\xi}$
and $D$, on $x$ and $y$ respectively. Hence, by choosing $dx=dy=ds$, the
fractional volume change can be written as:

$\displaystyle\frac{dV}{V}\ =\ \frac{1}{D}\ (\ \frac{\partial
D\xi_{x}}{\partial x}\ +\ \frac{\partial D\xi_{y}}{\partial y}\ )$ (2)

If we describe the motion of the particle via Newton’s equation and we combine
it with the equation of state $p\ =-K\frac{dV}{V}$ (where $K$ is the bulk
modulus of air), it is possible to prove the following relationship between
the displacement and the fractional volume change just obtained:

$\displaystyle\frac{\partial^{2}\bm{\xi}}{\partial t^{2}}\ =\
\frac{c^{2}}{D}\bm{\nabla}(\ \frac{\partial D\xi_{x}}{\partial x}\ +\
\frac{\partial D\xi_{y}}{\partial y}\ )$ (3)

with $c$ defined as the speed of sound in air. Finally, by following Fletcher
and Rossing [10] and considering $D$ independent of time, we can obtain from
Equation 3 the following acoustic wave equation:

$\displaystyle\frac{\partial^{2}p}{\partial t^{2}}\ =\ \frac{c^{2}}{D}(\
\frac{\partial}{\partial x}(D\frac{\partial p}{\partial x})\ +\
\frac{\partial}{\partial y}(D\frac{\partial p}{\partial y})\ )$ (4)

where $p$ is sound pressure. The equation describes air wave motion in a 3D
space constrained in one dimension; given the similarity with the standard 2D
acoustic wave equation, this 3D space can be considered as an extended 2D
(2.5D) space. By symmetrically constraining the 2.5D space along the $y$
dimension too (hence turning $D$ into a circular area), and by assuming plane-
wave propagation only (due to large distance from the source), Equation 4
reduces to Webster’s horn equation [11].

###  2.2 2.5D FDTD Wave Solver

The procedure followed for the design of the 2.5D wave solver shares many
similarities with the case of the 2D vocal tract model described in [9]. We
decomposed the 2.5D acoustic wave equation into extended versions of the 2D
equation of continuity and 2D equation of motion:

$\displaystyle\frac{\partial p}{\partial t}\ =\ -\frac{\rho c^{2}}{D}(\
\frac{\partial Dv_{x}}{\partial x}\ +\ \frac{\partial Dv_{y}}{\partial y}\ )$
(5) $\displaystyle\beta\frac{\partial\bm{v}}{\partial t}\ +\
\left(1-\beta\right)\bm{v}\ =\ -\beta^{2}\frac{\bm{\nabla}p}{\rho}\ +\
\left(1-\beta\right)\bm{v}_{b}$ (6)

where $\bm{v}$ is the 2D acoustic particle velocity. Equation 6 does not
include any explicit dependence on the depth of the considered space (via
$D$). This is due to the fact that the equation of motion derives from
Newton’s second law, whose 2D and 2.5D forms coincide. However, similarly to
what proposed by Allen and Raghuvanshi [12], we augmented the obtained
standard 2D equation with the scalar field $\beta(x,y,t)$. By varying between
1 and 0, this term allows for the transition between the momentum equation and
the enforcement of a prescribed velocity $\bm{v_{b}}$. Through this mechanism
it is possible to enforce boundary conditions (more details in the next
paragraph) and to simulate dynamic geometries [9].

We applied a standard 2D Yee scheme, where each grid point consists of a
squared 2D cell [12, 13]. Per each time step $n$, pressure values $p^{(n)}$
are sampled across the whole domain at the center of the cells, while the
$v_{x}^{(n)}$ and $v_{y}^{(n)}$ components of the velocity vectors are sampled
on the edges shared with the right and top neighbor cell respectively. Hence,
in a generic cell at discrete coordinates $(\hat{x},\ \hat{y})$, $p$ will be
sampled at $(\hat{x},\ \hat{y})$, while $v_{x}$ at $(\hat{x}+1/2,\ \hat{y})$
and $v_{y}$ at $(\hat{x},\ \hat{y}+1/2)$. This leads to the following discrete
update rules of Equations 5 and 6, where we denote with
$\bm{\widetilde{\nabla}}$ the standard discrete spatial derivatives as
performed in 2D FDTD:

$\displaystyle p^{(n+1)}\ =\ \frac{\bar{D}p^{(n)}-\rho c^{2}\Delta t\
\bm{\widetilde{\nabla}}\cdot\bm{V}^{(n)}}{\bar{D}}$ (7)
$\displaystyle\bm{v}^{(n+1)}\ =\ \frac{\beta\bm{v}^{(n)}-\beta^{2}\Delta t\
\widetilde{\nabla}p^{(n+1)}/\rho\ +\ \Delta t(1-\beta)\bm{v_{b}}}{\beta+\Delta
t(1-\beta)}$ (8)

with:

$\displaystyle\bm{V}\ =\ (\ D_{(x)}v_{x},\ D_{(y)}v_{y}\ )$ (9)

The terms $\bar{D}$, $D_{(x)}$ and $D_{(y)}$ are the components of the depth
map of the 2.5D space and they effectively act as extra impedance terms in two
dimensions. Their relation with the field $D$ will be explained in detail in
the next subsection.

In line with the work of Takemoto and Mokthari [1], the vocal tract’s walls
are simulated by adapting the local reactive boundary approach originally
proposed by Yokota et al. [14]. This is done by means of setting $\beta=0$ and
$\bm{v_{b}}=\rho c\ \mu\ p_{w}\bm{n}$, with $p_{w}$ equal to the pressure
value sampled from the cell in front of the wall, $\bm{n}$ the unit vector
normal to it and directed towards the wall itself and $\mu$ the boundary
admittance. Similarly, arbitrary excitation can be injected into the domain
via boundary cells (i.e, $\beta=0$), by setting $\bm{v_{b}}$ equal to the
output velocity of a glottal model; this can be coupled with the vocal tract
by feeding back the pressure value of the cell in front of the excitation
velocity stream. Alike its 2D version [9], the 2.5D vocal tract model allows
for the coupling with both 1D [15, 16] and 2D glottal models [17].

###  2.3 Depth Map

The depth map is a discretized version of the continuous field $D$, processed
according to the described 2D scheme. Since every individual FDTD cell holds
three acoustic parameters sampled in three different locations (center, right
and top), the discretization of $D$ will likewise produce three depth values
per each cell, namely $\bar{D}$, $D_{(x)}$ and $D_{(y)}$. As a result, in the
FDTD grid the pressure values are aligned with $\bar{D}$ values, while the
velocity vector components are aligned with $D_{(x)}$ and $D_{(y)}$ values.

In a practical scenario, the field $D$ can be replaced by the 3D model of the
geometry whose acoustics is being analyzed. Even if typically composed of
discrete elements (e.g., vertices), the model can be treated as a continuous
geometry via interpolation. However, $D$ always represents a symmetric domain;
this is not necessarily true for a 3D model, hence additional processing is
required to derive the depth map from it. In the case of a generic 3D tube,
the depth map is computed as follows:

  1. 1.

a 2D contour is extracted as the intersection between the tube and its mid-
sagittal plane (i.e., 2D vocal tract geometry). The contour itself represents
the vocal tract’s walls, while the mid-sagittal plane becomes the $x/y$ plane
of the 2.5D space;

  2. 2.

in every wall cell, $\bar{D}$, $D_{(x)}$ and $D_{(y)}$ are set to zero;

  3. 3.

in every cell outside of the the tube’s walls, $\bar{D}$, $D_{(x)}$ and
$D_{(y)}$ are assigned a fixed open space depth;

  4. 4.

for every cell with discrete coordinates $(\hat{x},\ \hat{y})$ inside the
tube’s walls, the two lines perpendicular to the $x/y$ plane and passing by
$(\hat{x}+1/2,\ \hat{y})$ and $(\hat{x},\ \hat{y}+1/2)$ are intersected with
the model (Figure 1). The length of the two resulting segments is assigned to
$D_{(x)}(\hat{x},\ \hat{y})$ and $D_{(y)}(\hat{x},\ \hat{y})$ respectively;

  5. 5.

for every cell with discrete coordinates $(\hat{x},\ \hat{y})$ inside the
tube’s walls, $D_{(x)}(\hat{x},\ \hat{y})$ and $D_{(y)}(\hat{x},\ \hat{y})$
are interpolated as $(\ D_{(x)}(\hat{x},\ \hat{y})+D_{x}(\hat{x}+1,\ \hat{y})\
)/2$ and $(\ D_{(y)}(\hat{x},\ \hat{y})+D_{(y)}(\hat{x},\ \hat{y}+1)\ )/2$
respectively;

  6. 6.

for every cell with discrete coordinates $(\hat{x},\ \hat{y})$ inside the
tube’s walls, the depth $\bar{D}(\hat{x},\ \hat{y})$ is obtained as $(\
D_{(x)}(\hat{x},\ \hat{y})+D_{(x)}(\hat{x}-1,\ \hat{y})+D_{(y)}(\hat{x},\
\hat{y})+D_{(y)}(\hat{x},\ \hat{y}-1)\ )/4$;

  7. 7.

a minimum depth threshold is defined and applied in every cell (depth values
are clamped). A typical threshold is at least an order of magnitude lower than
the smallest non-zero depth value found across the domain.

The result of this process is a single-axis symmetric equivalent of the
original 3D model. The cross-sections of such a 2.5D shape have same areas and
widths (i.e., extension along $y$) as the 3D model’s sections; moreover, their
shapes can be irregular, but always symmetric along $z$, as showed in Figure
1.

The interpolations at steps 5 and 6 are introduced to satisfy the assumption
that $D$ is slowly varying in space, required to obtain Equation 1 (Section
2.1). This enhances the overall stability of the solver too.

Figure 1: on the left, depth map extraction ($D_{(x)}$ solid line, $D_{(y)}$
dashed line) from a slice $C$ of an irregular tube, at coordinate $x=\hat{x}$
(for simplicity we assume the same shape at $\hat{x}+1/2$). On the right,
equivalent single-axis symmetric cross-section $C^{\prime}$, as represented in
the 2.5D FDTD.

###  2.4 Model Applicability

The development of the 2.5D vocal tract model is still in progress, yet some
of its features can be analyzed already.

The use of the depth map imposes to slightly change the cross-sectional shapes
of the original geometry by means of single-axis symmetry. As a result, the
acoustics of the 3D and of the 2.5D vocal tracts will be somewhat different in
the portion of the spectrum that is beyond 5 kHz [3]. However, the 2.5D
representation is capable of preserving curvature on the $x/y$ plane, area
progression as well as several details of the varying irregular cross-
sectional shapes of real vocal tract geometries. These are features beyond the
capabilities of 1D and 2D models, and whose combined acoustics effects are yet
to be studied.

At the current stage of development, using the 2.5D model to analyze the
acoustics of realistic vocal tracts is impractical though. The model
capitalizes on the solution of a lossless system (Equations 5, 6). The solver
adds wall impedance by means of the 2D walls boundary condition (Section 2.2),
but the losses happening on the surface described by the depth map (i.e., 2.5D
boundaries) are largely ignored, as the extra depth-related impedance terms
$\bar{D}$, $D_{(x)}$ and $D_{(y)}$ model only the effects of the surface’s
spatial derivative111This scenario is analogous to the case of Webster’s horn
equation, where the absolute size of the modeled tube does not affect
propagation.. In other words, there is no consistent way of matching the non-
planar modes of a generic 3D shape.

This issue can be easily overcome for the modeling of 3D tubes with circular
cross-sections, which can also include bending. By extending the calculation
of modes in plates proposed by Fletcher and Rossing [10], it is possible to
prove that non-planar modes in each section of a lossless 2.5D tube are
identical to the ones found in its 2D mid-sagittal cross-section. As a result,
we can apply the methodology proposed by Arnela and Guasch [4] to match the
first non-planar mode of a circular cross-section, by scaling the radii of
every section of the tube by the constant factor $\frac{0.5\pi}{1.84}$ and by
modifying the depth map accordingly.

##  3 Experiments and Results

###  3.1 Model Validation

A comparative study was carried out to validate the model’s precision, by
using as reference a high-accuracy 3D Finite Element Method (FEM) [18]. To
this end, we computed the transfer functions of a set of vocal tract-like
geometries, we extracted the corresponding formants’ positions and compared
them with the results obtained with the 3D FEM. The study was carried out for
the following static vowels: /a/, /i/ and /u/. We chose to construct the 3D
tubes from Story’s area functions [19]. This dataset defines a standard in
vocal tracts’ acoustic analysis and is ideal to work with circular cross-
sections (see Section 2.4).

As further validation, we compared the model’s time performance with two
different versions of the same 2D FDTD: a standard serial implementation
(2DS), and a highly optimized parallel implementation (2DP). Both the models
were timed while computing the transfer function of Story’s /u/. The target of
such comparisons was to estimate the extra computational cost introduced by
the depth terms, as well as the time boost made possible by a parallel
pipeline.

Figure 2: 2.5D tube discretization (mid-sagittal contour and depth map) for
vowel /a/.

###  3.2 Experimental Setup

The 2.5D FDTD acoustic wave solver was implemented in MATLAB. The resolution
of the FDTD grid was set to $\Delta s=0.74$ mm; this is the largest value that
preserves the geometrical details of the tubes and was determined empirically.
Any resolution above this threshold produces less accurate results. Figure 2
illustrates the resulting 2.5D tube geometry for vowel /a/. We set the size of
the grid to $270\times 45$ cells, to obtain a domain that could fit any of the
three tubes. The temporal resolution $\Delta t$ of the simulation is
restricted by the Courant-Friedrichs-Lewy condition in 2D, $\Delta t\leq\Delta
s/\sqrt{2}c$, where $c$ is the speed of sound ($c=350$ m/s within the vocal
tract). We therefore set $\Delta t=1.51\times 10^{-6}$ s, equal to a
simulation rate of 661,500 Hz.

Air density was set to $\rho=1.14$ kg/m2s and all wall cells were assigned a
standard boundary admittance $\mu=0.005$. Since the 2.5D FDTD does not yet
include a way to precisely model the radiation losses, we implemented an open-
end termination at the mouth opening by imposing a Dirichlet boundary
condition in the same fashion of what done in [4]. Following the procedure
described in [9], the impulse responses were obtained by exciting the model
with a band-passed velocity pulse injected from the glottal end; the pressure
variation was recorded via a microphone placed 3 mm inside the mouth opening,
keeping track of computational times. Per each vocal tract, the MATLAB solver
serially iterated across the full grid to simulate a total of 50 ms of audio.
The application ran on a workstation equipped with an Intel Core i7-8700K
processor.

The same parameters were maintained also for the 2DS and the 2DP, but all the
depth-related impedance terms were removed from their solvers222These settings
produce incorrect 2D formants. Yet, the 2D models were used for time
comparisons only and their output was ignored.. Furthermore, while the 2DS was
implemented in MATLAB, the 2DP was implemented as a shader as described in [9]
and ran on a Nvidia GTX 1080 graphics card.

###  3.3 Results

We applied the Fast Fourier Transform to the impulse responses to obtain the
vocal tracts’ transfer functions for the three vowels. The positions of the
first 8 formants were extracted and compared with the results of the 3D FEM
approach. The exact formants’ positions for the 3D FEM can be found in
Arnela’s PhD dissertation [20]. Table 1 shows the formants’ positional
differences and the percentage error between the 2D and the 3D simulation
results for /a/, /i/ and /u/ respectively. For each vowel, the 2.5D run-time
for producing 50 ms of audio was 15 minutes and 21 seconds; for the 2DS it was
14 minutes and 41 seconds, while for the 2DP 1.4 seconds.

Formants |  /a/ |  /i/ |  /u/
---|---|---|---
F1 |  4Hz 0.57% |  -3Hz -1.14% |  1Hz 0.38%
F2 |  -28Hz -2.62% |  49Hz 2.32% |  -37Hz -4.88%
F3 |  -11Hz -0.36% |  50Hz 1.66% |  36Hz 1.59%
F4 |  56Hz 1.37% |  82Hz 1.98% |  17Hz 0.47%
F5 |  -59Hz -1.17% |  1Hz 0.01% |  107Hz 2.56%
F6 |  -23Hz -0.40% |  -13Hz -0.22% |  -102Hz -2.01%
F7 |  -5Hz -0.07% |  -95Hz -1.44% |  -46Hz -0.75%
F8 |  25Hz 0.32% |  35Hz 0.45% |  26Hz 0.39%
Table 1: positional errors of the first 8 formants in 2.5D, computed for vowel
/a/, /i/, and /u/ with respect to 3D FEM values.

##  4 Discussion and Conclusion

The acoustic analysis of vowels’ area functions is a standard validation
methodology that allows to test the fundamental properties of a wave
propagation model. In the case of the 2.5 FDTD, all the three vowels are
characterized by positional errors that tend to stay below 2%, and in several
cases do not go above 1%. To achieve comparable results, a standard 2D FDTD
requires a spatial resolution that is almost three times higher than the value
used in this test [9]. Such a level of accuracy suggests that 2.5D equations
are a valid means to model how waves travel in a constrained space, and that
the proposed solver is correct. Furthermore, the affordable spatial resolution
remarkably decreases the computational load of the simulation compared to 2D
models’ massive grids [4, 7, 9].

The current serial MATLAB implementation of the 2.5D FDTD proves to be
extremely lightweight, even when compared with the performance of the equally-
sized 2DS. The time difference between the two systems is 20 s (2%); this
value is quite small and derives from the fact that the instruction sets of
the 2.5D and 2DS solvers differ by four multiplications only (see [9]). Given
the similar computational load of the two systems, it is safe to assume that
an optimized parallel implementation of the 2.5D running on a modern graphics
card will reach quasi-real-time performance, close to what achieved with the
2DP.

These results gain particular interest in the context of the future
development of the 2.5D vocal tract model. We are currently working on the
inclusion of free-radiation effects (via 2.5D perfect matching layers) and on
the validation of the modeling of the vocal tract’s curvature. The inclusion
of these features is relatively straightforward and will make the system
capable of simulating geometrical and acoustic details still unavailable in
one-dimensional and 2D models. The following step will consist of the
introduction in the acoustic equations of extra loss terms, to model 2.5D
boundaries. By doing so, it will be possible to simulate the effects of
varying irregular cross-sections, thus fully exploring the potential of the
2.5D model.

##  5 Acknowledgements

This work is supported by the Natural Sciences and Engineering Research
Council (NSERC) of Canada and Canadian Institutes for Health Research (CIHR).

## References

  * [1] H. Takemoto, P. Mokhtari, and T. Kitamura, “Acoustic analysis of the vocal tract during vowel production by finite-difference time-domain method,” _The Journal of the Acoustical Society of America_ , vol. 128, no. 6, pp. 3724–3738, 2010.
  * [2] T. Vampola, J. Horáček, A.-M. Laukkanen, and J. G. Švec, “Human vocal tract resonances and the corresponding mode shapes investigated by three-dimensional finite-element modelling based on ct measurement,” _Logopedics Phoniatrics Vocology_ , vol. 40, no. 1, pp. 14–23, 2015.
  * [3] M. Arnela, S. Dabbaghchian, R. Blandin, O. Guasch, O. Engwall, A. Van Hirtum, and X. Pelorson, “Influence of vocal tract geometry simplifications on the numerical simulation of vowel sounds,” _The Journal of the Acoustical Society of America_ , vol. 140, no. 3, pp. 1707–1718, 2016.
  * [4] M. Arnela and O. Guasch, “Two-dimensional vocal tracts with three-dimensional behavior in the numerical generation of vowels,” _The Journal of the Acoustical Society of America_ , vol. 135, no. 1, pp. 369–379, 2014.
  * [5] J. Mullen, D. M. Howard, and D. T. Murphy, “Waveguide physical modeling of vocal tract acoustics: flexible formant bandwidth control from increased model dimensionality,” _IEEE Transactions on Audio, Speech, and Language Processing_ , vol. 14, no. 3, pp. 964–971, 2006.
  * [6] M. Speed, D. Murphy, and D. M. Howard, “Characteristics of two-dimensional finite difference techniques for vocal tract analysis and voice synthesis,” in _Tenth Annual Conference of the International Speech Communication Association_ , 2009.
  * [7] Y. Wang, H. Wang, J. Wei, and J. Dang, “Mandarin vowel synthesis based on 2d and 3d vocal tract model by finite-difference time-domain method,” in _Proceedings of The 2012 Asia Pacific Signal and Information Processing Association Annual Summit and Conference_. IEEE, 2012, pp. 1–4.
  * [8] J. Mullen, D. M. Howard, and D. T. Murphy, “Real-time dynamic articulations in the 2-d waveguide mesh vocal tract model,” _IEEE Transactions on Audio, Speech, and Language Processing_ , vol. 15, no. 2, pp. 577–585, 2007.
  * [9] V. Zappi, A. Vasuvedan, A. Allen, N. Raghuvanshi, and S. Fels, “Towards real-time two-dimensional wave propagation for articulatory speech synthesis,” in _Proceedings of Meetings on Acoustics 171ASA_ , vol. 26, no. 1. ASA, 2016, p. 045005.
  * [10] N. H. Fletcher and T. D. Rossing, _The physics of musical instruments_. Springer Science & Business Media, 2012.
  * [11] G. Ungeheuer, _Elemente einer akustischen Theorie der Vokalartikulation_. Springer-Verlag, 1962\.
  * [12] A. Allen and N. Raghuvanshi, “Aerophones in flatland: Interactive wave simulation of wind instruments,” _ACM Transactions on Graphics (TOG)_ , vol. 34, no. 4, p. 134, 2015.
  * [13] K. Yee, “Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media,” _IEEE Transactions on antennas and propagation_ , vol. 14, no. 3, pp. 302–307, 1966.
  * [14] T. Yokota, S. Sakamoto, and H. Tachibana, “Visualization of sound propagation and scattering in rooms,” _Acoustical science and technology_ , vol. 23, no. 1, pp. 40–46, 2002.
  * [15] K. Ishizaka and J. L. Flanagan, “Synthesis of voiced sounds from a two-mass model of the vocal cords,” _Bell system technical journal_ , vol. 51, no. 6, pp. 1233–1268, 1972.
  * [16] B. H. Story and I. R. Titze, “Voice simulation with a body-cover model of the vocal folds,” _The Journal of the Acoustical Society of America_ , vol. 97, no. 2, pp. 1249–1260, 1995.
  * [17] A. Vasudevan, V. Zappi, P. Anderson, and S. Fels, “A fast robust 1d flow model for a self-oscillating coupled 2d fem vocal fold simulation.” in _INTERSPEECH_ , 2017, pp. 3482–3486.
  * [18] M. Arnela and O. Guasch, “Finite element computation of elliptical vocal tract impedances using the two-microphone transfer function method,” _The Journal of the Acoustical Society of America_ , vol. 133, no. 6, pp. 4197–4209, 2013.
  * [19] B. H. Story, “Comparison of magnetic resonance imaging-based vocal tract area functions obtained from the same speaker in 1994 and 2002,” _The Journal of the Acoustical Society of America_ , vol. 123, no. 1, pp. 327–335, 2008\.
  * [20] M. Arnela Coll, “Numerical production of vowels and diphthongs using finite element methods,” Ph.D. dissertation, Universitat Ramon Llull, 2015.


# Timelike orbits around accelerating black holes

Mohammad Bagher Jahani Poshteh<EMAIL_ADDRESS>School of Physics, Institute for
Research in Fundamental Sciences (IPM), P.O. Box 19395-5531, Tehran, Iran

###### Abstract

We study the geodesics of massive particles around accelerating Schwarzschild
black hole. We show that the radius of the innermost stable circular orbit
increases and the angular momentum of particle at this orbit decreases by
increasing the acceleration. Apart from quantitative influence, the
acceleration qualitatively changes the physics. We show that in accelerating
black hole spacetime there exist an outermost stable circular orbit in flat,
de Sitter, and anti de Sitter background. Investigations of radial geodesics
show that the acceleration acts like a repulsive force in the sense that test
particles around accelerating black holes can move radially outward, unless
there exist a large negative value of cosmological constant in the background
to compensate the repulsive force. We also investigate the precession of
perihelion of orbits around accelerating black holes. The precession would be
larger compared to the non-accelerating case. It is also shown that the
precession in anti de Sitter background could be opposite to the particle
motion.

##  I Introduction

The study of relativistic motion of particles in gravitational field began
short after the general theory of relativity was introduced [1, 2, 3]. In
light of the investigations of lightlike geodesics we have found many
interesting phenomena, e.g. gravitational lensing [4, 5, 6, 7, 8, 9] and black
hole shadow [10, 11, 12]. These features are in common among black holes and
naked singularities [13, 14] and ultra-compact objects [15, 16].

The investigation of timelike geodesics around massive objects could help us
understand the nature of the objects. In this paper we are interested in
timelike geodesics in accelerating black hole spacetime. The motion of massive
test particles around non-accelerating black holes in general relativity has
been thoroughly studied in [17]. In [18, 19, 20] fundamental frequencies of
particles’ motion around Kerr black holes have been investigated. Circular
motions on equatorial plane around Kerr–Newman black holes and Kerr–Newman
naked singularities are compared in [21]. In [22] null and timelike geodesics
on equatorial plane of a distorted Schwarzschild black hole are studied. In
[23], for a class of naked singularity spacetimes, it has been shown that the
perihelion precession can be in the opposite direction of particle’s motion.
Precession of timelike bound orbits in Kerr spacetime has also been studied
and it is found that the precession is positive [24].

Black holes can be pair produced on cosmic strings [25, 26, 27] as well as in
a de Sitter [28, 29, 30] or magnetic field [31, 32, 33] background (see [34]
for the black hole production rate on cosmic string in a de Sitter space with
background magnetic field). On the other hand, we might have primordial black
holes got attached to cosmic strings in the early Universe [35]. All of these
black holes will be accelerating (due to the tension of the cosmic string [35]
and/or the force exerted by the positive cosmological constant and/or the
magnetic field [34]).

Accelerating black holes could evolve to supermassive black holes [36] (see
also [37]). However, if these black holes have played a role in structure
formation, their velocity should be small [36]; which means that the
acceleration should be small. It has been speculated, from observational point
of view, that the black hole at the center of our Galaxy is connected to
cosmic string [38].

Study of null geodesics around accelerating black hole has recently attracted
some attention [39, 40, 41]. The shadow of the accelerating black hole has
been studied in [42, 43]. On the other hand, thermodynamics of these black
holes have been investigated in [44, 45, 46]. Near horizon symmetries of
accelerating black holes have also been studied [47].

Accelerating black holes are described by C metric. In [48], for the particles
co-accelerating with the black hole at constant distance, null and timelike
geodesics are found in the standard $\\{x,y\\}$ coordinates [49] as well as
Weyl coordinates [50] and the coordinates adapted to the boost-rotation
symmetry [51]. Generalization to AdS C metric [52] and spinning C metric [53]
have also been done. In [54], by studying the effective potential of test
particles around the black hole, accessible regions for null/timelike
geodesics are presented. Stability of geodesics is also discussed.

The aim of this paper is to study the motion of massive particles in
accelerating black hole spacetime with cosmological constant. We would like to
study radial geodesics as well as the innermost stable circular orbit (ISCO)
of particles around accelerating supermassive black holes. We would also like
to investigate the perihelion precession of such orbits. In our examples we
consider a non-rotating black hole with the same mass as the black hole at the
center of Milky Way Galaxy, Sgr A*. We also use the orbital data of S2 star
around Sgr A*. Considering the upper bound on the acceleration of supermassive
black holes [36], we show that we can take S2 to be nearly on equatorial plane
of accelerating black hole during one complete period.

Study of radial motion of test particles in accelerating black hole spacetime
shows that these particles can move radially outward. This means that
accelerating black hole exert some sort of repulsive force on test particles.
We also find the interesting result that there exist an upper bound on the
radius of stable circular orbits around the accelerating black hole. These
features are in common among accelerating black holes and black holes in de
Sitter spacetime.

We also investigate the precession of orbits around accelerating black holes.
We show that the precession of perihelion is larger for larger values of black
hole acceleration. However, for black holes in anti de Sitter background, the
precession can be negative (the orbit precesses in the opposite direction of
the motion).

The outline of our paper is as follows. In the next section we present the
Lagrangian of motion on a plane perpendicular to the direction of acceleration
(equatorial plane) and study the radial motion on this plane. Stable circular
orbits on the equatorial plane is studied in Sec. III. In Sec. IV we study the
precession of the orbits. We conclude our paper in Sec. V. We work in
geometric units where $G=c=1$ and use mostly-positive signature for the
spacetime metric.

##  II Radial motions on equatorial plane

It has been shown in [36] that the velocity of supermassive black holes
connected to cosmic string should be less than $100\,{\rm km}/{\rm s}$ so that
they could be captured by galaxies during the structure formation. This, in
turn, constrain tension of the cosmic string to $\mu\lesssim
10^{-19}\sqrt{M/M_{\astrosun}}$, where $M$ is the mass of supermassive black
hole connected to the cosmic string. For the black hole at the center of our
Galaxy, with $M_{{\rm SgrA}^{*}}\simeq 4.23\times 10^{6}M_{\astrosun}$ [55],
the tension of the cosmic string attached to it would be $\mu\lesssim
2.06\times 10^{-16}$. Therefore we find from the Newton second law,
$\mu=\alpha M_{{\rm SgrA}^{*}}$, that the acceleration is $\alpha\lesssim
3.29\times 10^{-26}{\rm m}^{-1}$.

The orbital period of S2 star around Sgr A* is about 16 years [56, 57]. We
take the acceleration of Sgr A* to be $\alpha=10^{-26}{\rm m}^{-1}$ — this is
about the largest value allowed by [36]. Assuming the initial velocity of the
black hole to be zero, the displacement of the black hole during the period of
S2 is $2.29\times 10^{8}{\rm m}$. This is much smaller than the Schwarzschild
radius of the black hole which is of the order of $10^{10}{\rm m}$. Therefore
we can assume that, if S2 starts its period near the equatorial plane of the
black hole, it will remain near this plane during one period of the orbit
around the slowly accelerating black hole.

Spacetime around uniformly accelerating black holes is described by C metric.
111This spacetime is of algebraic type D. It has been proved in [58] that
vacuum type D spacetimes admit Killing and Killing-Yano tensors if (and only
if) they are without acceleration. (A)dS C metric has been studied in, for
example, [29, 59, 60, 61, 62]. For a black hole of mass $m$ and uniform
acceleration $\alpha$, in a background spacetime with cosmological constant
$\Lambda$, the metric can be written in the following form [63]

$\displaystyle ds^{2}$ $\displaystyle=$ $\displaystyle\frac{1}{(1+\alpha
r\cos\theta)^{2}}$ $\displaystyle\times$
$\displaystyle\left[-Q(r)dt^{2}+\frac{dr^{2}}{Q(r)}+\frac{r^{2}d\theta^{2}}{P(\theta)}+P(\theta)r^{2}\sin^{2}\theta
d\phi^{2}\right],$

where

$\displaystyle Q(r)$ $\displaystyle=$
$\displaystyle(1-\alpha^{2}r^{2})\left(1-\frac{2m}{r}\right)-\frac{\Lambda}{3}r^{2},$
$\displaystyle P(\theta)$ $\displaystyle=$ $\displaystyle 1+2\alpha
m\cos\theta.$ (2)

The coordinate $t$ ranges over all of ${\rm I\\!R}$ and $0\leq\theta\leq\pi$.
The periodicity of $\phi$ is $-C_{0}\pi\leq\phi\leq C_{0}\pi$. One can choose
$C_{0}$ so as to eliminate the conical singularity at one of the poles. We
take $C_{0}=(1+2\alpha m)^{-1}$ to eliminate the singularity at $\theta=0$
[54]. 222Setting $C_{0}=(1-2\alpha m)^{-1}$, on the other hand, eliminates the
singularity at $\theta=\pi$ [54]. For other choices of $C_{0}$ neither the
singularity at $\theta=0$ nor the one at $\theta=\pi$ would be eliminated and
one has cosmic strings at both poles.

We would like to study the geodesics on the equatorial plane of the
accelerating black holes. We assume that the direction of acceleration is
perpendicular to this plane. On equatorial plane the metric
(LABEL:eqn:metric:griffiths) reduces to

$ds^{2}=-Qdt^{2}+\frac{dr^{2}}{Q}+r^{2}d\theta^{2}+r^{2}d\phi^{2},$ (3)

with the metric function given by Eq. (2).

The metric function $Q$ has been plotted in Fig. 1. We have taken some large
values of acceleration and cosmological constant so that we could better
illustrate the behavior of the metric function in different cases. For the
case of flat and de Sitter background the metric function has two zeros for
positive $r$ and has a maximum between them. The behavior of the metric
function in the case of anti de Sitter background is more interesting. For
$\Lambda>-3\alpha^{2}$ there are two zeros for positive $r$, however if
$\Lambda<-3\alpha^{2}$ there would be only one zero and the metric function
goes to $+\infty$ as $r$ increases. 333For $\Lambda>-3\alpha^{2}$ there is a
maximum between two zeros of the metric function. As
$\Lambda\rightarrow(-3\alpha^{2})^{+}$ the point at which the maximum appears
goes to $+\infty$.

The region between the two zeros of $Q$ — one associated to event horizon of
the black hole, the other to acceleration horizon — is the domain of outer
communication in which $Q>0$. For $\Lambda<-3\alpha^{2}$ this domain contains
any radius larger than the radius of black hole’s event horizon. The motion of
test body around the black hole is in the region of outer communication.
Therefore we only consider the range of $r$ which is in the domain of outer
communication.

Figure 1: Top: The metric function of an accelerating black hole in flat
(dashed blue plot) and de Sitter with $m^{2}\Lambda=10^{-4}$ (solid red plot)
background. Bottom: The metric function of an accelerating black hole in anti
de Sitter background with $m^{2}\Lambda=-10^{-3}$ (dashed blue plot) and
$m^{2}\Lambda=-10^{-6}$ (solid red plot). We have taken $m\alpha=10^{-2}$.

The Lagrangian governing the geodesic motion is
$2\mathcal{L}=g_{\gamma\sigma}\dot{x}^{\gamma}\dot{x}^{\sigma}=-\mu^{2}$,
where $\mu$ is the rest mass of the infalling particle and dot, for timelike
geodesics, represents differentiation with respect to the proper time (for a
detailed study of timelike geodesics in Schwarzschild and Kerr background see
[17]). Using the metric (3) the Lagrangian on equatorial plane would be

$\mathcal{L}=\frac{1}{2}\left(-Q\dot{t}^{2}+\frac{\dot{r}^{2}}{Q}+r^{2}\dot{\phi}^{2}\right).$
(4)

The canonical momenta $p_{t}=-\frac{\partial\mathcal{L}}{\partial\dot{t}}$ and
$p_{\phi}=\frac{\partial\mathcal{L}}{\partial\dot{\phi}}$ are conserved
because of the symmetry of the metric, i.e.
$\dot{p}_{t}=-\frac{\partial\mathcal{L}}{\partial t}=0$ and
$\dot{p}_{\phi}=\frac{\partial\mathcal{L}}{\partial\phi}=0$. They are the
energy and angular momentum of the infalling particle and can be found to be

$E=p_{t}=Q\dot{t},\qquad L_{z}=p_{\phi}=r^{2}\dot{\phi}.$ (5)

We note here that if the acceleration has a component parallel to the
equatorial plane which is of the same order or larger than the perpendicular
component, then the third component of the angular momentum is no longer a
constant of motion. In this case geodesic equations cannot be analytically
integrated. By using these quantities we can rewrite Eq. (4) as

$\dot{r}^{2}=\tilde{E}^{2}-Q\left[1+\frac{\tilde{L}_{z}^{2}}{r^{2}}\right],$
(6)

where $\tilde{E}$ and $\tilde{L}_{z}$ are respectively the energy and angular
momentum per unit rest mass and the dot now represents differentiation with
respect to proper time per unit rest mass. Note also that the second term on
the right hand side of Eq. (6),

$\tilde{V}=Q\left[1+\frac{\tilde{L}_{z}^{2}}{r^{2}}\right],$ (7)

is the effective potential.

Now we consider radial motion in accelerating black hole spacetime. 444For a
black hole accelerating in $z$ direction, we consider radial motion in $x-y$
plane. For radial geodesics we have $\tilde{L}_{z}=0$. Therefore Eq. (6)
reduces to

$\dot{r}^{2}=\tilde{E}^{2}-Q.$ (8)

If a test particle is falling from rest at $r_{0}$ we find
$\tilde{E}^{2}=Q(r_{0})$.

For non-accelerating Schwarzschild black hole in anti de Sitter background it
has been shown that the test particle always plunges into the black hole [64].
This is due to the attractive force generated by the negative cosmological
constant. However, for non-accelerating Schwarzschild black hole in de Sitter
background, it is known that the test particle can move radially
outward/inward if $r_{0}$ is larger/smaller than a specific value [65].

The acceleration of the test particle which is at rest at $r_{0}$ is given by
[65]

$\ddot{r}=\frac{Q^{\prime}(r)\left[\dot{r}^{2}-Q(r_{0})\right]}{2Q(r)}.$ (9)

For $\Lambda>-3\alpha^{2}$ the metric function $Q$ has two zeros at positive
$r$ (see Fig. 1). Let us denote the smaller root by $r_{+}$ and the larger one
by $r_{++}$. Between these two roots the metric function has a maximum at
$r_{M}$. Suppose a test particle starts its motion at $r_{0}$. For
$r_{+}<r_{0}<r_{M}$ the slope of the metric function is positive, therefore
the initial acceleration $\ddot{r}=-Q^{\prime}(r_{0})/2$ is negative. This
means that the test particle plunges into the black hole. For
$r_{M}<r_{0}<r_{++}$, however, the acceleration is positive and the test
particle moves radially outward. On the other hand, in anti de Sitter
background with $\Lambda<-3\alpha^{2}$ the test particle always plunges into
the event horizon.

Therefore we find that acceleration acts like a positive cosmological
constant. Similarity between a positive cosmological constant and cosmic
strings has already been claimed in the context of black hole pair production
in the author’s earlier work [34]. We also see that a negative cosmological
constant can compensate the acceleration ($\Lambda=-3\alpha^{2}$ completely
compensate the effect of $\alpha$).

##  III Stable circular orbits

Here we consider the circular geodesics for which $\tilde{L}_{z}\neq 0$. We
have plotted the effective potential in Fig. 2 for different values of the
angular momentum in the region near the event horizon of the black hole (we
take approximately the observed value of the cosmological constant,
$\Lambda\simeq 10^{-52}{\rm m}^{-2}$ [66]). In this region, we see that for
$\tilde{L}_{z}$ slightly above $\tilde{L}_{z,{\rm ISCO}}$ there are a minimum
and a maximum in the plot of the potential which are associated to stable and
unstable circular orbits, respectively. The minimum gets closer to the black
hole as $\tilde{L}_{z}$ decreases. For $\tilde{L}_{z}=\tilde{L}_{z,{\rm
ISCO}}$ we have the ISCO which is at the inflection point of the effective
potential.

Figure 2: The effective potential for different values of angular momentum in
the region near the black hole horizon. We have taken $m\alpha=10^{-26}$,
$m^{2}\Lambda=10^{-52}$, and $\tilde{L}_{z}=3m$ (dotted black plot),
$\tilde{L}_{z}=\tilde{L}_{z,{\rm ISCO}}\approx 3.4641m$ (dashed red plot), and
$\tilde{L}_{z}=4m$ (solid blue plot). The maximum and minimum of the solid
blue plot are respectively at $r\approx 4m$ and $r\approx 12m$. The inflection
point of the dashed red plot is at $r\approx 6m$.

We find that the radius of ISCO and the angular momentum of particle at this
orbit depends on the acceleration. In Fig. 3, using numerical techniques, we
have plotted the radius of ISCO and the angular momentum at ISCO as a function
of the acceleration. We see that the radius of ISCO increases and the angular
momentum of particle at ISCO decreases by increasing the acceleration. This
means that if the black hole is accelerating, the infalling particle can orbit
on a stable circular orbit around the hole with a larger radius and smaller
angular momentum.

Figure 3: Top: The radius of ISCO as a function of the acceleration of the
black hole. Bottom: The angular momentum of the particle on ISCO. We have
taken $m^{2}\Lambda=10^{-52}$.

For non-accelerating black holes, it is known that there exists an outermost
stable circular orbit (OSCO) in a background spacetime with a positive
cosmological constant [67, 68]. We have find that, for accelerating black
holes, OSCO exists for flat and AdS background as well as in dS spacetime.
This point has been presented in Fig. 4. Let us call the region in which the
OSCO takes place as the far region. We see that for $\tilde{L}_{z}$ slightly
below $\tilde{L}_{z,{\rm OSCO}}$ there are a minimum and a maximum in the plot
of the potential which are, respectively, associated to stable and unstable
circular orbits. The minimum goes away from the black hole as $\tilde{L}_{z}$
increases. For $\tilde{L}_{z}=\tilde{L}_{z,{\rm OSCO}}$ we have the OSCO which
is at the inflection point of the effective potential.

In Fig. 5 we see that the radius of OSCO and the angular momentum of particle
at this orbit decrease by increasing the acceleration. We also note that there
is a critical value of acceleration for which OSCO and ISCO would be on top of
each other. (For the value of the cosmological constant taken in Figs. 3 and 5
the critical acceleration is $m\alpha\simeq 1.61466\times 10^{-2}$.)

Figure 4: The effective potential for different values of angular momentum in
the region far from the black hole horizon. We have taken $m\alpha=10^{-3}$,
$\Lambda=0$, and $\tilde{L}_{z}=7.05800m$ (dotted black plot),
$\tilde{L}_{z}=\tilde{L}_{z,{\rm OSCO}}\approx 7.06232m$ (dashed red plot),
and $\tilde{L}_{z}=7.06700m$ (solid blue plot).

Figure 5: Top: The radius of OSCO as a function of the acceleration of the
black hole. Bottom: The angular momentum of the particle on OSCO. We have
taken $m^{2}\Lambda=10^{-52}$.

Behavior of the potential in both the near and the far region is presented in
Fig. 6. For the case in which $\tilde{L}_{z}<\tilde{L}_{z,{\rm ISCO}}$
($\tilde{L}_{z}>\tilde{L}_{z,{\rm OSCO}}$) the potential has a maximum in the
far (near) region. If the energy of the particle coming from infinity is
larger than the maximum value of the potential it will plunge into the black
hole. The case in which $\tilde{L}_{z,{\rm
ISCO}}<\tilde{L}_{z}<\tilde{L}_{z,{\rm OSCO}}$ is more interesting and has
been presented in the middle (dashed red) plot of Fig. 6. In this case there
are two maxima at $r_{{\rm max}_{1}}$ and $r_{{\rm max}_{2}}$ (assuming
$r_{{\rm max}_{1}}<r_{{\rm max}_{2}}$ and $\tilde{V}(r_{{\rm
max}_{1}})>\tilde{V}(r_{{\rm max}_{2}})$) and one minimum between them at
$r_{{\rm min}}$.

Ignoring the particles confining to $r<r_{{\rm max}_{1}}$, depending on the
energy of the particle, the following cases could happen:

  * •

For a particle approaching the black hole, if the energy satisfies
$\tilde{E}^{2}>\tilde{V}(r_{{\rm max}_{1}})$, the particle would plunge into
the black hole.

  * •

If $\tilde{V}(r_{{\rm max}_{2}})<\tilde{E}^{2}<\tilde{V}(r_{{\rm max}_{1}})$,
the particle reaches a minimum radius in the range $r_{{\rm
max}_{1}}<r<r_{{\rm min}}$ and then escape to infinity.

  * •

If $\tilde{V}(r_{{\rm min}})<\tilde{E}^{2}<\tilde{V}(r_{{\rm max}_{2}})$, and
if the particle is coming from infinity it reaches a minimum radius of
$r>r_{{\rm max}_{2}}$. The other possibility is that the particle rotates
around the black hole in a bound orbit with radius $r_{{\rm
max}_{1}}<r<r_{{\rm max}_{2}}$.

  * •

If $\tilde{E}^{2}<\tilde{V}(r_{{\rm min}})$ the particles coming from infinity
would reach a minimum at $r>r_{{\rm max}_{2}}$ (here $r$ is larger than the
case in which $\tilde{V}(r_{{\rm min}})<\tilde{E}^{2}<\tilde{V}(r_{{\rm
max}_{2}})$) and then escape to infinity.

Also if $\tilde{E}^{2}=\tilde{V}(r_{{\rm max}_{1}})$ or
$\tilde{E}^{2}=\tilde{V}(r_{{\rm max}_{2}})$ the particle will rotate the
black hole on an unstable circular orbit, and if
$\tilde{E}^{2}=\tilde{V}(r_{{\rm min}})$ the particle rotate the black hole on
a stable circular orbit.

Figure 6: The effective potential for different values of angular momentum in
the near and far region from the black hole horizon. We have taken
$m\alpha=10^{-2}$, $\Lambda=0$, and $\tilde{L}_{z}=3.4m<\tilde{L}_{z,{\rm
ISCO}}$ (dotted black plot), $\tilde{L}_{z,{\rm
ISCO}}<\tilde{L}_{z}=3.6m<\tilde{L}_{z,{\rm OSCO}}$ (dashed red plot), and
$\tilde{L}_{z}=3.7m>\tilde{L}_{z,{\rm OSCO}}$ (solid blue plot). Qualitatively
similar behaviors take place in the de Sitter and anti de Sitter background as
well.

We also note that for $\Lambda<-3\alpha^{2}$ the potential goes to $+\infty$
as $r$ increases. However for $\Lambda>-3\alpha^{2}$ (whether in de Sitter,
flat, or anti de Sitter background) the potential goes to $-\infty$ as $r$
increases. These point have been illustrated in Fig. 7.

Figure 7: The effective potential for $m\alpha=10^{-6}$ and
$m^{2}\Lambda=-10^{-3}$ (solid blue plot) and $m\alpha=10^{-2}$ and
$m^{2}\Lambda=-10^{-10}$ (dashed red plot). The solid blue plot satisfies
$\Lambda<-3\alpha^{2}$, while for the dashed red plot $\Lambda>-3\alpha^{2}$.
We have taken $\tilde{L}_{z}=3m$ in both plots.

##  IV Precession of perihelion

In this section we investigate the precession of perihelion for orbits around
accelerating black holes. Eq. (4) can be written as

$\left(\frac{dr}{d\phi}\right)^{2}=\frac{r^{4}\tilde{E}^{2}}{\tilde{L}_{z}^{2}}-\frac{Qr^{4}}{\tilde{L}_{z}^{2}}-Qr^{2}.$
(10)

At the perihelion and aphelion we have $dr/d\phi=0$. If we denote the radius
of the perihelion and aphelion respectively by $r_{p}$ and $r_{a}$, we can
write the energy and angular momentum per unit rest mass as

$\displaystyle\tilde{E}^{2}$ $\displaystyle=$
$\displaystyle\frac{Q(r_{a})Q(r_{p})\left(r_{p}^{2}-r_{a}^{2}\right)}{Q(r_{a})r_{p}^{2}-Q(r_{p})r_{a}^{2}},$
(11) $\displaystyle\tilde{L}_{z}^{2}$ $\displaystyle=$
$\displaystyle\frac{r_{a}^{2}r_{p}^{2}\left[Q(r_{p})-Q(r_{a})\right]}{Q(r_{a})r_{p}^{2}-Q(r_{p})r_{a}^{2}}.$
(12)

As we mentioned earlier the periodicity of $\phi$ is $2C_{0}\pi$, with
$C_{0}=(1+2\alpha m)^{-1}$.Therefore, the precession can be found, by using
Eq. (10), as

$\displaystyle\Delta\phi$ $\displaystyle=$ $\displaystyle
2\left[\phi(r_{a})-\phi(r_{p})\right]-2(1+2\alpha m)^{-1}\pi$ $\displaystyle=$
$\displaystyle
2\int_{r_{p}}^{r_{a}}\frac{dr}{\sqrt{\frac{r^{4}\tilde{E}^{2}}{\tilde{L}_{z}^{2}}-\frac{Qr^{4}}{\tilde{L}_{z}^{2}}-Qr^{2}}}-2(1+2\alpha
m)^{-1}\pi,$

where $\tilde{E}^{2}$ and $\tilde{L}_{z}^{2}$ are given by eqs. (11) and (12).

The radius of the perihelion and aphelion are related to the eccentricity $e$
and semi-major axis $a$ of the orbit through $r_{p}=(1-e)a$ and
$r_{a}=(1+e)a$. As an example we take the semi-major axis and eccentricity of
S2 in its orbit around Sgr A*. They are $a=1.543\times 10^{14}{\rm m}$ and
$e=0.88$ [69].

In the top panel of Fig. 8 we have plotted the precession
(LABEL:eqn:precession) as a function of acceleration for S2. We see that the
precession angle increases by increasing the acceleration of the central black
hole. Also in the bottom panel of Fig. 8 we see that the precession increase
by increasing the cosmological constant. It is very interesting that the
precession angle is negative for $m^{2}\Lambda\lesssim-1.51363\times
10^{-16}$.

Figure 8: Top: Precession of perihelion as a function of the acceleration. We
have taken $\Lambda=0$. Bottom: Precession of perihelion as a function of the
cosmological constant. We have taken $\alpha=0$.

##  V Concluding remarks

We have studied the geodesics of massive particles on equatorial plane of
accelerating black holes. We have only considered the component of the
acceleration which is perpendicular to the plane of the particle motion. If
the component parallel to this plane is considerably large, then the angular
momentum of the particle is no longer a constant of motion and we would only
have one constant of motion (the energy). In that case the geodesic equations
cannot be solved analytically.

We have found several new and interesting results. The study of radial
geodesics and stable circular orbits around an accelerating black hole show
that there exists some sort of similarity between accelerating black hole and
black holes in de Sitter background. Radially outward motion of test particles
in accelerating black hole spacetime points out to the fact that acceleration,
like a positive cosmological constant, generate a repulsive force. (However,
in the case of de Sitter space, the “cosmic acceleration” is spherically
symmetric, whereas here in acceleration is in $z$ direction.) Another
similarity is the existence of OSCO in accelerating black hole spacetime which
also exists for particles orbiting a black hole in de Sitter background.

We have found that the precession of perihelion is larger around an
accelerating black hole. It is also shown that in anti de Sitter background
the precession can be negative. We postpone a detailed study of the precession
in anti de Sitter background to a future work. It would also be very
interesting to study the timelike geodesics around rotating accelerating black
holes.

##  VI Acknowledgments

We are grateful to M.M. Sheikh-Jabbari for reading several pre-publication
versions of this paper and providing us with useful comments. We also
acknowledge the support of Iran Science Elites Federation and the hospitality
of the University of Guilan. We would also like to thank B. Turimov and O.
Yunusov for pointing out the mistake in the plots on the top panels of Fig. 3
and Fig. 5 in the earlier versions of this paper.

## References

  * Einstein [1916] A. Einstein, The foundation of the general theory of relativity, Annalen Phys. 49, 769 (1916).
  * Dyson  _et al._ [1920] F. W. Dyson, A. S. Eddington, and C. Davidson, A Determination of the Deflection of Light by the Sun’s Gravitational Field, from Observations Made at the Total Eclipse of May 29, 1919, Phil. Trans. Roy. Soc. 220A, 291 (1920).
  * Einstein [1936] A. Einstein, Lens-Like Action of a Star by the Deviation of Light in the Gravitational Field, Science 84, 506 (1936).
  * Schneider  _et al._ [1992] P. Schneider, J. Ehlers, and E. E. Falco,  _Gravitational lenses_ (Springer, 1992).
  * Darwin [1959] C. Darwin, The Gravity Field of a Particle, Proc. R. Soc. A 249, 180 (1959).
  * Darwin [1961] C. Darwin, The Gravity Field of a Particle. II, Proc. R. Soc. A 263, 39 (1961).
  * Jahani Poshteh and Mann [2019] M. B. Jahani Poshteh and R. B. Mann, Gravitational Lensing by Black Holes in Einsteinian Cubic Gravity, Phys. Rev. D 99, 024035 (2019), arXiv:1810.10657 [gr-qc] .
  * Virbhadra and Ellis [2000] K. S. Virbhadra and G. F. R. Ellis, Schwarzschild black hole lensing, Phys. Rev. D 62, 084003 (2000), arXiv:astro-ph/9904193 .
  * Virbhadra [2009] K. S. Virbhadra, Relativistic images of Schwarzschild black hole lensing, Phys. Rev. D 79, 083004 (2009), arXiv:0810.2109 [gr-qc] .
  * Synge [1966] J. L. Synge, The Escape of Photons from Gravitationally Intense Stars, Mon. Not. Roy. Astron. Soc. 131, 463 (1966).
  * Hennigar  _et al._ [2018] R. A. Hennigar, M. B. Jahani Poshteh, and R. B. Mann, Shadows, Signals, and Stability in Einsteinian Cubic Gravity, Phys. Rev. D 97, 064041 (2018).
  * Akiyama  _et al._ [2019] K. Akiyama, A. Alberdi, W. Alef, K. Asada, R. Azulay, A.-K. Baczko, D. Ball, M. Baloković, J. Barrett, D. Bintley, _et al._ , First m87 event horizon telescope results. iv. imaging the central supermassive black hole, The Astrophysical Journal Letters 875, L4 (2019).
  * Virbhadra and Ellis [2002] K. S. Virbhadra and G. F. Ellis, Gravitational lensing by naked singularities, Physical Review D 65, 103004 (2002).
  * Joshi  _et al._ [2020] A. B. Joshi, D. Dey, P. S. Joshi, and P. Bambhaniya, Shadow of a naked singularity without photon sphere, Physical Review D 102, 024022 (2020).
  * Cunha and Herdeiro [2018] P. V. Cunha and C. A. Herdeiro, Shadows and strong gravitational lensing: a brief review, General Relativity and Gravitation 50, 1 (2018).
  * Shaikh  _et al._ [2019] R. Shaikh, P. Banerjee, S. Paul, and T. Sarkar, Analytical approach to strong gravitational lensing from ultracompact objects, Physical Review D 99, 104040 (2019).
  * Chandrasekhar [1998] S. Chandrasekhar,  _The mathematical theory of black holes_ , Vol. 69 (Oxford university press, 1998).
  * Wilkins [1972] D. C. Wilkins, Bound geodesics in the kerr metric, Physical Review D 5, 814 (1972).
  * Glampedakis and Kennefick [2002] K. Glampedakis and D. Kennefick, Zoom and whirl: Eccentric equatorial orbits around spinning black holes and their evolution under gravitational radiation reaction, Physical Review D 66, 044002 (2002).
  * Fujita and Hikida [2009] R. Fujita and W. Hikida, Analytical solutions of bound timelike geodesic orbits in kerr spacetime, Classical and Quantum Gravity 26, 135002 (2009).
  * Pugliese  _et al._ [2013] D. Pugliese, H. Quevedo, and R. Ruffini, Equatorial circular orbits of neutral test particles in the kerr-newman spacetime, Physical Review D 88, 024042 (2013).
  * Shoom  _et al._ [2016] A. A. Shoom, C. Walsh, and I. Booth, Geodesic motion around a distorted static black hole, Phys. Rev. D 93, 064019 (2016), arXiv:1512.01566 [gr-qc] .
  * Bambhaniya  _et al._ [2019] P. Bambhaniya, A. B. Joshi, D. Dey, and P. S. Joshi, Timelike geodesics in Naked Singularity and Black Hole Spacetimes, Phys. Rev. D 100, 124020 (2019), arXiv:1908.07171 [gr-qc] .
  * Bambhaniya  _et al._ [2021] P. Bambhaniya, D. N. Solanki, D. Dey, A. B. Joshi, P. S. Joshi, and V. Patel, Precession of timelike bound orbits in Kerr spacetime, Eur. Phys. J. C 81, 205 (2021), arXiv:2007.12086 [gr-qc] .
  * Hawking and Ross [1995] S. W. Hawking and S. F. Ross, Pair production of black holes on cosmic strings, Phys. Rev. Lett. 75, 3382 (1995), arXiv:gr-qc/9506020 .
  * Eardley  _et al._ [1995] D. M. Eardley, G. T. Horowitz, D. A. Kastor, and J. H. Traschen, Breaking cosmic strings without monopoles, Phys. Rev. Lett. 75, 3390 (1995), arXiv:gr-qc/9506041 .
  * Ashoorioon and Mann [2014] A. Ashoorioon and R. B. Mann, Black Holes as Beads on Cosmic Strings, Class. Quant. Grav. 31, 225009 (2014), arXiv:1402.2072 [hep-th] .
  * Mellor and Moss [1989] F. Mellor and I. Moss, Black Holes and Quantum Wormholes, Phys. Lett. B 222, 361 (1989).
  * Mann and Ross [1995] R. B. Mann and S. F. Ross, Cosmological production of charged black hole pairs, Phys. Rev. D 52, 2254 (1995), arXiv:gr-qc/9504015 .
  * Dias and Lemos [2004] O. J. C. Dias and J. P. S. Lemos, Pair creation of de Sitter black holes on a cosmic string background, Phys. Rev. D 69, 084006 (2004), arXiv:hep-th/0310068 .
  * Garfinkle and Strominger [1991] D. Garfinkle and A. Strominger, Semiclassical Wheeler wormhole production, Phys. Lett. B 256, 146 (1991).
  * Hawking  _et al._ [1995] S. W. Hawking, G. T. Horowitz, and S. F. Ross, Entropy, Area, and black hole pairs, Phys. Rev. D 51, 4302 (1995), arXiv:gr-qc/9409013 .
  * Dowker  _et al._ [1994] F. Dowker, J. P. Gauntlett, S. B. Giddings, and G. T. Horowitz, On pair creation of extremal black holes and Kaluza-Klein monopoles, Phys. Rev. D 50, 2662 (1994), arXiv:hep-th/9312172 .
  * Ashoorioon and Jahani Poshteh [2021] A. Ashoorioon and M. B. Jahani Poshteh, Black hole pair production on cosmic strings in the presence of a background magnetic field, Phys. Lett. B 816, 136224 (2021), arXiv:2006.16983 [hep-th] .
  * Vilenkin [1985] A. Vilenkin, Cosmic Strings and Domain Walls, Phys. Rept. 121, 263 (1985).
  * Vilenkin  _et al._ [2018] A. Vilenkin, Y. Levin, and A. Gruzinov, Cosmic strings and primordial black holes, JCAP 11, 008, arXiv:1808.00670 [astro-ph.CO] .
  * Gußmann [2021] A. Gußmann, Polarimetric signatures of the photon ring of a black hole that is pierced by a cosmic axion string, (2021), arXiv:2105.06659 [astro-ph.HE] .
  * Morris  _et al._ [2017] M. R. Morris, J.-H. Zhao, and W. Goss, A nonthermal radio filament connected to the galactic black hole?, The Astrophysical Journal Letters 850, L23 (2017).
  * Lim [2021] Y.-K. Lim, Null geodesics in the $C$ metric with a cosmological constant, Phys. Rev. D 103, 024007 (2021), arXiv:2011.05503 [gr-qc] .
  * Frost and Perlick [2021] T. C. Frost and V. Perlick, Lightlike Geodesics and Gravitational Lensing in the Spacetime of an Accelerating Black Hole, Class. Quant. Grav. 38, 085016 (2021), arXiv:2010.11908 [gr-qc] .
  * Ashoorioon  _et al._ [2021] A. Ashoorioon, M. B. Jahani Poshteh, and R. B. Mann, Distinguishing a Slowly Accelerating Black Hole by Differential Time Delays of Images, (2021), arXiv:2110.13132 [gr-qc] .
  * Grenzebach  _et al._ [2015] A. Grenzebach, V. Perlick, and C. Lämmerzahl, Photon Regions and Shadows of Accelerated Black Holes, Int. J. Mod. Phys. D 24, 1542024 (2015), arXiv:1503.03036 [gr-qc] .
  * Zhang and Jiang [2021] M. Zhang and J. Jiang, Shadows of accelerating black holes, Phys. Rev. D 103, 025005 (2021), arXiv:2010.12194 [gr-qc] .
  * Appels  _et al._ [2016] M. Appels, R. Gregory, and D. Kubiznak, Thermodynamics of Accelerating Black Holes, Phys. Rev. Lett. 117, 131303 (2016), arXiv:1604.08812 [hep-th] .
  * Anabalón  _et al._ [2018] A. Anabalón, M. Appels, R. Gregory, D. Kubizňák, R. B. Mann, and A. Ovgün, Holographic Thermodynamics of Accelerating Black Holes, Phys. Rev. D 98, 104038 (2018), arXiv:1805.02687 [hep-th] .
  * Anabalón  _et al._ [2019] A. Anabalón, F. Gray, R. Gregory, D. Kubizňák, and R. B. Mann, Thermodynamics of Charged, Rotating, and Accelerating Black Holes, JHEP 04, 096, arXiv:1811.04936 [hep-th] .
  * Anabalón  _et al._ [2021] A. Anabalón, S. Brenner, G. Giribet, and L. Montecchio, Closer look at black hole pair creation, Phys. Rev. D 104, 024044 (2021), arXiv:2103.05782 [hep-th] .
  * Pravda and Pravdova [2001] V. Pravda and A. Pravdova, Coaccelerated particles in the C metric, Class. Quant. Grav. 18, 1205 (2001), arXiv:gr-qc/0010051 .
  * Kinnersley and Walker [1970] W. Kinnersley and M. Walker, Uniformly accelerating charged mass in general relativity, Phys. Rev. D 2, 1359 (1970).
  * Bonnor [1983] W. Bonnor, The sources of the vacuum c-metric, General relativity and Gravitation 15, 535 (1983).
  * Bičák and Schmidt [1989] J. Bičák and B. Schmidt, Asymptotically flat radiative space-times with boost-rotation symmetry: The general structure, Physical Review D 40, 1827 (1989).
  * Chamblin [2001] A. Chamblin, Capture of bulk geodesics by brane world black holes, Class. Quant. Grav. 18, L17 (2001), arXiv:hep-th/0011128 .
  * Pravda and Pravdová [2002] V. Pravda and A. Pravdová, On the spinning c-metric, in  _Gravitation: Following The Prague Inspiration: A Volume in Celebration of the 60th Birthday of Jirí Bičák_ (World Scientific, 2002) pp. 247–262.
  * Lim [2014] Y.-K. Lim, Geodesic motion in the vacuum C-metric, Phys. Rev. D 89, 104016 (2014), arXiv:1405.2611 [gr-qc] .
  * Johannsen  _et al._ [2016] T. Johannsen, A. E. Broderick, P. M. Plewa, S. Chatzopoulos, S. S. Doeleman, F. Eisenhauer, V. L. Fish, R. Genzel, O. Gerhard, and M. D. Johnson, Testing general relativity with the shadow size of sgr ${\mathrm{a}}^{*}$, Phys. Rev. Lett. 116, 031101 (2016).
  * Boehle  _et al._ [2016] A. Boehle, A. Ghez, R. Schödel, L. Meyer, S. Yelda, S. Albers, G. Martinez, E. Becklin, T. Do, J. Lu, _et al._ , An improved distance and mass estimate for sgr a* from a multistar orbit analysis, The Astrophysical Journal 830, 17 (2016).
  * Hees  _et al._ [2017] A. Hees _et al._ , Testing General Relativity with stellar orbits around the supermassive black hole in our Galactic center, Phys. Rev. Lett. 118, 211101 (2017), arXiv:1705.07902 [astro-ph.GA] .
  * Demianski and Francaviglia [1980] M. Demianski and M. Francaviglia, Separability structures and killing-yano tensors in vacuum type-d space-times without acceleration, International Journal of Theoretical Physics 19, 675 (1980).
  * Dias and Lemos [2003] O. J. Dias and J. P. Lemos, Pair of accelerated black holes in an anti–de sitter background: The ads c metric, Physical Review D 67, 064001 (2003).
  * Podolsky and Griffiths [2001] J. Podolsky and J. B. Griffiths, Uniformly accelerating black holes in a de Sitter universe, Phys. Rev. D 63, 024006 (2001), arXiv:gr-qc/0010109 .
  * Podolsky [2002] J. Podolsky, Accelerating black holes in anti-de Sitter universe, Czech. J. Phys. 52, 1 (2002), arXiv:gr-qc/0202033 .
  * Chen  _et al._ [2015] Y. Chen, Y.-K. Lim, and E. Teo, New form of the c metric with cosmological constant, Physical Review D 91, 064014 (2015).
  * Griffiths and Podolskỳ [2006] J. B. Griffiths and J. Podolskỳ, A new look at the plebański–demiański family of solutions, International Journal of Modern Physics D 15, 335 (2006).
  * Cruz  _et al._ [2005] N. Cruz, M. Olivares, and J. R. Villanueva, The Geodesic structure of the Schwarzschild anti-de Sitter black hole, Class. Quant. Grav. 22, 1167 (2005), arXiv:gr-qc/0408016 .
  * Jaklitsch  _et al._ [1989] M. Jaklitsch, C. Hellaby, and D. Matravers, Particle motion in the spherically symmetric vacuum solution with positive cosmological constant, General relativity and gravitation 21, 941 (1989).
  * Padmanabhan [2003] T. Padmanabhan, Cosmological constant: The Weight of the vacuum, Phys. Rept. 380, 235 (2003), arXiv:hep-th/0212290 .
  * Howes [1979] R. Howes, Existence and stability of circular orbits in a schwarzschild field with nonvanishing cosmological constant, Australian Journal of Physics 32, 293 (1979).
  * Boonserm  _et al._ [2020] P. Boonserm, T. Ngampitipan, A. Simpson, and M. Visser, Innermost and outermost stable circular orbits in the presence of a positive cosmological constant, Phys. Rev. D 101, 024050 (2020), arXiv:1909.06755 [gr-qc] .
  * Eisenhauer  _et al._ [2003] F. Eisenhauer, R. Schödel, R. Genzel, T. Ott, M. Tecza, R. Abuter, A. Eckart, and T. Alexander, A geometric determination of the distance to the galactic center, The Astrophysical Journal Letters 597, L121 (2003).