# Investigating the significance of zero-point motion in small molecular
clusters of sulphuric acid and water

Jake L. Stinson  Email<EMAIL_ADDRESS>Department of Physics and
Astronomy and London Centre for Nanotechnology, University College London,
Gower Street, London, WC1E 6BT, United Kingdom  Shawn M. Kathmann  Physical
Sciences Division, Pacific Northwest National Laboratory, Richland, Washington
99352, United States  Ian J. Ford  Department of Physics and Astronomy and
London Centre for Nanotechnology, University College London, Gower Street,
London, WC1E 6BT, United Kingdom

###### Abstract

The nucleation of particles from trace gases in the atmosphere is an important
source of cloud condensation nuclei (CCN), and these are vital for the
formation of clouds in view of the high supersaturations required for
homogeneous water droplet nucleation. The methods of quantum chemistry have
increasingly been employed to model nucleation due to their high accuracy and
efficiency in calculating configurational energies; and nucleation rates can
be obtained from the associated free energies of particle formation. However,
even in such advanced approaches, it is typically assumed that the nuclei have
a classical nature, which is questionable for some systems. The importance of
zero-point motion (also known as quantum nuclear dynamics) in modelling small
clusters of sulphuric acid and water is tested here using the path integral
molecular dynamics (PIMD) method at the density functional theory (DFT) level
of theory. The general effect of zero-point motion is to distort the mean
structure slightly, and to promote the extent of proton transfer with respect
to classical behaviour. In a particular configuration of one sulphuric acid
molecule with three waters, the range of positions explored by a proton
between a sulphuric acid and a water molecule at 300 K (a broad range in
contrast to the confinement suggested by geometry optimisation at 0 K) is
clearly affected by the inclusion of zero point motion, and similar effects
are observed for other configurations.

##  I Introduction

The role of sulphuric acid in the formation of cloud condensation nuclei (CCN)
is believed to be significant(Roy, 2009; Zhang  _et al._ , 2012), on account
of its low vapour pressure, relatively high atmospheric concentration and its
affinity to water. However, simple attempts to understand the binary
nucleation of sulphuric acid and water in detail have proved problematic. It
is clear that classical nucleation theory (CNT) is insufficient for describing
this process, since the critical cluster size suggested from experimental data
appears to be small, and consequently several extensions and alternatives have
been studied (Ford, 2004; Vehkamäki, 2006; Napari  _et al._ , 2010).

One approach, the use of atomistic models that explicitly treat individual
molecules or atoms within numerical simulations, has proliferated as a
consequence of increasing computational power; especially based on quantum
chemistry methods which treat the electronic interactions explicitly. Popular
quantum chemistry methods include electronic density functional theory (DFT)
(Re  _et al._ , 1999; Bandy and Ianni, 1998; Ianni and Bandy, 1999; Arstila
_et al._ , 1998; Larson  _et al._ , 2000; Ding  _et al._ , 2003; A. Natsheh
_et al._ , 2004; Arrouvel  _et al._ , 2005; Ding and Laasonen, 2004; Nadykto
_et al._ , 2008; Kurtén  _et al._ , 2009) and Møller-Plesset perturbation
theory (MPn where n refers to the order of the perturbation) (Larson  _et al._
, 2000; Ding and Laasonen, 2004; Kurtén  _et al._ , 2009; Temelso  _et al._ ,
2012). The usual strategy is to identify the lowest energy molecular
configuration and then to use the rigid-rotor-harmonic-approximation (RRHO) to
compute free energies, and thereby investigate cluster stability and
nucleation through specific growth and decay routes.

The Born-Oppenheimer approximation (Martin, 2004) is employed by both DFT and
MPn. It involves the separation of the wavefunctions of electrons and nuclei
followed by a classical treatment of the dynamics of the nuclei. The DFT
approach has been used to describe sulphuric acid and water clusters (Choe
_et al._ , 2007; Anderson  _et al._ , 2008; Hammerich  _et al._ , 2008). In
simulations based on such approaches, the sulphuric acid and water system has
been observed to exhibit proton transfers. Such events are of particular
importance in this system and a challenge to the modelling. A question that
arises is whether we can account for such processes correctly while
representing the nuclei as classical particles. Might a quantum treatment of
the proton dynamics be more accurate? Perhaps the additional uncertainty in
proton position can alter the delicate balance between neutral and ionised
structures? In this paper we employ Path Integral Molecular Dynamics (PIMD) to
study the quantum nuclear degrees of freedom (also known as zero-point motion)
of sulphuric acid/water molecular clusters to address this question. A
particular issue for consideration is the level of hydration of a single
sulphuric acid molecule that is required for proton transfer to occur, a
matter that can be addressed either by zero temperature calculations or
dynamics performed at finite temperature. It has been suggested that the
threshold is around three or more water molecules(Arrouvel  _et al._ , 2005).
Transfer of the second proton was studied by Ding and Laasonen (Ding and
Laasonen, 2004) who concluded that it is likely for a level of hydration of
around eight or nine water molecules.

PIMD emulates the quantum behaviour of a particle by using a classical
quasiparticle or bead description, a detailed derivation of which is given by
Tuckerman (2010). The PIMD method has been shown to have a significant effect
on the properties in some hydrogen bonded systems (Li  _et al._ , 2011; Walker
and Michaelides, 2010). PIMD has been employed previously together with a
parametrised version of the PM6 model (Stewart, 2007) (a semi-empirical model
of electronic structure) to study sulphuric acid and water clusters (Kakizaki
_et al._ , 2009; Sugawara  _et al._ , 2011). Kakizaki  _et al._ (2009)
concluded that the PIMD technique (using the normal mode transformation
(Tuckerman, 2010)) increased thermal fluctuations and produced more liquid-
like behaviour in systems at a temperature of ${\rm 250\>K}$ (Kakizaki  _et
al._ , 2009). Sugawara  _et al._ (2011) studied the degree of hydration
required for the first and second ionisation events for the sulphuric acid
molecule, and concluded that the first ionisation takes place when four water
molecules are present in the cluster in agreement with earlier work (Arrouvel
_et al._ , 2005). The second ionisation event occurred in the presence of
${\rm 10-12}$ water molecules in contrast with the study by Ding and Laasonen
(Ding and Laasonen, 2004) though the latter was based on geometry optimisation
techniques rather than on molecular dynamics. As the purpose of this paper was
to gauge the importance of zero-point motion in the sulphuric acid and water
system as accurately as possible, it was decided to use DFT rather than the
semi-empirical PM6 model developed by Kakizaki  _et al._ (2009) .

We study the importance of zero-point motion in a small cluster of sulphuric
acid and water using PIMD (Feynman and Hibbs, 2010; Tuckerman, 2010) as
implemented in the CASTEP code(Clark  _et al._ , 2005). Section II describes
the theory used, section III details our results, and section IV concludes our
study where we comment on the significance of zero-point motion in the
sulphuric acid and water system.

Figure 1: A 16 bead representation of a system containing one sulphuric acid
and four water molecules: the distribution of bead positions conveys the
quantum uncertainty.

##  II Methods

According to the PIMD technique each particle (nucleus) is represented by a
set of quasiparticles (known as beads) connected by harmonic springs. The
following Hamiltonian describing the bead dynamics can be derived using the
Trotter approximation (Tuckerman, 2010):

$\mathcal{H}(x_{k},p_{k})=\sum_{k=1}^{P}\left[\frac{p_{k}^{2}}{2m^{\prime}}+\frac{1}{2}m\omega_{P}^{2}(x_{k+1}-x_{k})^{2}+\frac{1}{P}U(x_{k})\right]$

under the condition $x_{P+1}=x_{1}$, where $P$ is the number of beads and
$x_{k}$ and $p_{k}$ are respectively the position and momentum of bead $k$.
$\omega_{P}$ is the harmonic frequency of the inter-bead springs and is given
by $\sqrt{P}/\beta\hbar$ where $\beta=(k_{B}T)^{-1}$ and $T$ and $k_{B}$ are
the system temperature and the Boltzmann constant respectively. While $m$
denotes the mass of the particle, the mass of the beads is represented by
$m^{\prime}$, and $U(x_{k})$ is the classical potential in which the particle
moves. The quantum nuclear behaviour is reflected in both the position and
momentum of the beads under the influence of this Hamiltonian, which is
controlled by the stiffness of the inter-bead springs. Since the latter is
proportional to the mass of the particle, the hydrogen nucleus is expected to
be the most susceptible to zero-point effects.

| Time step [${\rm fs}$] | Simulation time [${\rm ps}$]
---|---|---
SATH ${\rm 1}$ bead | ${\rm 1.00}$ | ${\rm 11.000}$
SATH ${\rm 4}$ bead | ${\rm 0.50}$ | $\,\,\,{\rm 1.500}$
SATH ${\rm 8}$ bead | ${\rm 0.25}$ | $\,\,\,{\rm 0.875}$
SATH ${\rm 16}$ bead | ${\rm 0.50}$ | ${\rm 10.673}$
SATH ${\rm 32}$ bead | ${\rm 0.50}$ | $\,\,\,{\rm 0.512}$
SAQH ${\rm 1}$ bead | ${\rm 1.00}$ | $\,\,\,{\rm 1.000}$
SAQH ${\rm 4}$ bead | ${\rm 0.50}$ | $\,\,\,{\rm 1.500}$
SAQH ${\rm 8}$ bead | ${\rm 0.50}$ | $\,\,\,{\rm 1.500}$
SAQH ${\rm 16}$ bead | ${\rm 0.50}$ | $\,\,\,{\rm 1.500}$
config H ${\rm 1}$ bead | ${\rm 1.00}$ | ${\rm 10.900}$
config H ${\rm 16}$ bead | ${\rm 1.00}$ | ${\rm 10.647}$
Table 1: Compilation of the simulation length and time step for the MD runs
performed. SATH refers to sulphuric acid trihydrate and SAQH refers to
sulphuric acid tetrahydrate, structures that assume typical configurations
shown in Figure 2a and 2b respectively. Config H refers to the trihydrate
configuration shown in Figure 4a. Note that the longest simulations were
performed for SATH and config H.

Figure 1 is a snapshot from a 16 bead simulation representing the behaviour of
a cluster of one sulphuric acid and four water molecules. The spatial
separation of the beads clearly illustrates the greater positional uncertainty
of the hydrogen nuclei compared to that of the oxygen and the sulphur nuclei.

(a)  (b)  (c)

Figure 2: Geometry optimised configurations for tri- and tetrahydrated (SATH
and SAQH) clusters are shown in (a) and (b) respectively. The labelling of
various hydrogen bonds is referred to in Section III.2. (c) shows the binding
energies of configurations (a) and (b) as a function of the system box size,
converging to values obtained by Temelso  _et al._ (2012) at the MP2 level.

Molecular dynamics simulations at $\mathrm{300\>K}$ incorporating both
classical nuclear dynamics and PIMD were performed using the CASTEP (Clark
_et al._ , 2005) (version ${\rm 5.5}$) code. The standard on-the-fly ultrasoft
pseudopotential provided internally by the CASTEP code was employed for all
calculations. The Perdew-Burke-Ernzerhof (Perdew  _et al._ , 1996) (PBE)
functional was used with a plane wave basis set. The PBE functional has been
found to perform well for hydrogen bonded systems(Thanthiriwatte  _et al._ ,
2011; Ireta  _et al._ , 2004). A cut off energy of ${\rm 550\>eV}$ was found
to converge the plane wave basis set sufficiently for all systems studied. A
time step of $\mathrm{1\>fs}$ was used for classical (single bead) simulations
and a time step of $\mathrm{0.5\>fs}$ or shorter was used for the PIMD
simulations due to the stiffness of the inter-bead springs. CASTEP utilizes
the Born-Oppenheimer version of ab initio MD and the Langevin thermostat with
a friction constant of ${\rm 0.01\>fs^{-1}}$ was used in all simulations. The
equilibration period was judged by observing when the running mean energy of
the system had relaxed (usually requiring less than $\mathrm{0.5\>ps}$) and
also by monitoring the distribution of cluster ’temperature’ (or kinetic
energy in the centre of mass frame), which ought to be approximately Gaussian
(Haile, 1997) with a standard deviation ($\sigma$) obeying
$\sigma/\left\langle T\right\rangle\sim N^{-1/2}$. A typical temperature
histogram satisfying this requirement is shown in Figure 3. The inset in
Figure 3 shows that the typical relaxation time was in the order of
$\mathrm{0.5\>ps}$ relaxation time. Initial configurations of sulphuric acid
and water identified from the literature were constructed under a classical
potential (MMFF94s) using the Avogadro (Hanwell  _et al._ , 2012) (version
${\rm 1.0.3}$) package. The choices of time step and simulation time for
various cases are given in Table 1.

Figure 3: Histogram of cluster kinetic energy (represented as a temperature)
from the equilibrated simulation referred to as SATH ${\rm 1}$ bead in Table
1. The inset shows the system’s energy as a function of time for the first
$\mathrm{1\>ps}$ of the same simulation.

Configurations of ${\rm[H_{2}SO_{4}][H_{2}O]_{n=3-4}}$ were studied at a
target temperature of ${\rm 300\>K}$. PIMD uses a certain number of beads to
approximate the zero-point motion, and cases with $P=1$, ${\rm 4}$, ${\rm 8}$,
${\rm 16}$ and ${\rm 32}$ beads were tested in this study. The staging
transformation(Tuckerman, 2010) was used for all PIMD simulations. The $P=1$
case represents the classical limit of the PIMD technique and corresponds to
the complete neglect of zero-point motion. The box size of the system was
optimised against MP${\rm 2}$ level data (Temelso  _et al._ , 2012) as shown
in Figure 2c(c). The binding energies, at zero temperature, of the two
configurations in Figures 2a and 2b are compared against MP${\rm 2}$ level
data. A box size of ${\rm 15\>\hbox{\AA}}$ was chosen as a compromise between
accuracy and computational demand.

##  III Results

###  III.1 DFT without zero-point motion

Molecular configurations likely to feature a dissociated sulphuric acid
molecule were identified from the literature and investigated. One such
configuration was labelled III-i-1 by Re  _et al._ (1999) and is illustrated
here in Figure 4a and denoted config H. Our single bead simulations at
$\mathrm{300\>K}$ show that the proton labelled H1 moves with considerable
freedom between oxygens O1 and O5. Furthermore, Figure 4b demonstrates an
anticorrelation between the length $R_{{\rm c}}$ of the dissociating bond
O1-H1 and the sum of the lengths of the neighbouring hydrogen bonds, labelled
O3-H7 and O4-H6 in Figure 4a, and denoted $R_{{\rm hy}}$. The formation of the
‘ionised’ state due to the switch to the O5-H1 bond (such that the value of
$R_{{\rm c}}$ is large) is seen to depend upon the prior existence of both the
neighbouring hydrogen bonds (namely a low value of $R_{{\rm hy}}$). If either
neighbouring hydrogen bond is broken the system remains ‘neutral’ (with a low
value of $R_{{\rm c}}$), which is not surprising since the configuration is
then similar to the SATH structure shown in Figure 2. This is an important
corollary to conclusions acquired from consideration of geometry optimisation
at $\mathrm{0\>K}$, where config H has been shown to ionise (Re  _et al._ ,
1999). At ${\rm 300\>K}$ the behaviour can most certainly not be represented
by harmonic fluctuations about an ionised mean structure and a free energy
based on the rigid-rotor-harmonic-approximation for this configuration would
fail due to significant anharmonic contributions. We shall return to this
system in the next section.

(a)  (b)

Figure 4: The configuration denoted config H is shown in (a) with labels that
identify certain O-H pairs. Plot (b) illustrates the probability density
(given in arbitrary units) as a function of two structural features labelled
$R_{{\rm c}}$ (the length of the covalent bond O1-H1) and $R_{{\rm hy}}$ (the
sum of the lengths of prospective hydrogen bonds O4-H6 and O3-H7), obtained at
DFT level, equivalent to using a single bead in PIMD. The associated potential
of mean force takes the form of a broad, shallow well where the ionisation of
the configuration is correlated with the status of the adjacent hydrogen
bonds, as denoted by $R_{{\rm hy}}$.

###  III.2 PIMD

A PIMD study was performed first for two low energy configurations (denoted
SATH and SAQH) identified in the literature (Re  _et al._ , 1999; Bandy and
Ianni, 1998) and shown in Figures 2a and 2b. It is envisaged that hydrogen
bonds, in particular those associated with the sulphuric acid, would be the
most susceptible to zero-point effects due to the inherent tendency of
sulphuric acid to dissociate. Figure 5 shows the average oxygen-oxygen
distance (${\rm d_{OO}}$) of specific hydrogen bonds as a function of the
number of beads representing atoms in the system. The bonds labelled hb1 and
hb2 in the SATH structure contract in length by around ${\rm 2-5\%}$ with
respect to the outcome of classical dynamics while the situation for hb3 is
less clear. Note that the longest simulations were performed for the single
bead and 16 bead representations of the SATH structure, as indicated in Table
1. For other cases shorter studies were performed to illustrate the trends,
though the accuracy of the results is lower.

Figure 5: The average oxygen-oxygen separation ${\rm d_{OO}}$ of specific
hydrogen bonds as a function of the number of beads used in the simulation.
Labels hb1 and hb2 refer to Figure 2a and hb3 is shown in Figure 2b. The error
bars were determined by the standard blocking procedure (Flyvbjerg and
Petersen, 1989; Frenkel and Smit, 2001) and a blocking length of
$\mathrm{0.256\>ps}$ was found to give independent sampling. The calculations
correspond to the cases listed in Table 1.

Next we examine in detail how the behaviour of the hydrogen atom in hydrogen
bond hb2 is affected by PIMD. This is explored by constructing a potential of
mean force (PMF) for the hydrogen, defined by:

$W(R,\beta)=-k_{B}T\ln g(R,\beta)$

where $R$ and $\beta$ are geometric parameters illustrated in Figure 6a and
$g(R,\beta)$ is the proportion of simulation snapshots with the hydrogen
located within the region defined by $R\rightarrow R+dR$ and
$\beta\rightarrow\beta+d\beta$ divided by the equivalent proportion for
noninteracting particles. For the PIMD simulations the centroid of the beads
representing the hydrogen atom was used to produce the PMF. The method is
described extensively by Kumar  _et al._ (2007). Figure 6b and 6c show the
PMFs acquired using classical MD and PIMD, respectively, for hydrogen bond
hb2.

The PMF plots in Figure 6 visualise the differences between the dynamics of
the hb2 bond in Figure 2a under classical MD and the PIMD schemes. Such a
comparison is limited by the computationally expensive techniques employed.
However it does offer an insight into the importance of zero-point effects in
small clusters of sulphuric acid and water. The main effect is a shift in the
minimum of the PMF of hydrogen bond length $R$ by about
$\mathrm{0.2\>}\textrm{\AA}$ going from the DFT to the PIMD result indicating
that the zero-point motion has a mean configurational influence on this bond.
Figure 6d is a one dimensional version of Figures 6b and 6b obtained by
integrating over the $\beta$ parameter.

(a)  (b)  (c)  (d)

Figure 6:  Contour plots of the potential of mean force $W(R,\beta)$ in units
of $k_{B}T$ for the hydrogen in the bond labelled hb2 in Figure 2a. The green
dashed lines indicate contour levels of $\mathrm{0}$, $\mathrm{-2}$,
$\mathrm{-4}$ and $\mathrm{-6}$. The coordinates for the PMF are defined by
sketch (a) and the method follows the approach described by Kumar  _et al._
(2007). Plot (b) shows results from standard DFT molecular dynamics and plot
(c) arises from PIMD using 16 beads. The simulation times are given in Table
1. Plot (d) shows a 1D version of plots (b) and (c) obtained by integrating
over the $\beta$ parameter. Figure 7: Plot of the O1-H1 bond length (denoted
as $\mathrm{R_{O1H1}}$) against time in config H in Figure 4a from the
$\mathrm{1}$ bead simulation as detailed in Table 1. The plot clearly
illustrates the motion between the neutral and ionised states. The horizontal
line drawn at ${\rm 1.22}\>{\rm\hbox{\AA}}$ provides a simple threshold
between covalent and hydrogen bond-like behaviour of the O1-H1 bond. Under
quantum nuclear dynamics the fraction of time spent above this threshold
increases in the order of 10%.

The effects of zero-point motion are clearly rather subtle. To explore this
further, we return to the delicate switching behaviour of the O1-H1-O5 bonds
discussed in section III.1 and contrast the classical and quantum nuclear
dynamics. Figure 7 illustrates the motion of the proton between the neutral
and ionised positions, discussed earlier, in terms of the O1-H1 bond length.
Which of the nuclei O1 or O5 was the nearest neighbour to the H1 nucleus (see
Figure 4a) was monitored to quantify this hopping behaviour. It was found that
in the classical case H1 was closer to O1 for $\mathrm{21.5\%}$ percent of the
simulation with standard error $\sigma_{SE}=3.2\%$ whereas in the 16 bead PIMD
simulation this figure dropped to $\mathrm{14.8\%}$ with $\sigma_{SE}=2.7\%$.
. This property was further investigated by defining a threshold for the O1-H1
bond length below which the system is considered neutral, and beyond which it
is better described as ionised. We define a ${\rm 1.22\>}\hbox{\AA}$ distance
to separate the two regimes, and this is shown as a horizontal line in Figure
7. For the classical dynamics, the percentage of time the system remains
neutral according to this criterion is ${\rm 20.1\%}$ with
${\rm\sigma_{SE}=2.9\%}$. An analysis of the PIMD simulation with 16 beads
yields a corresponding percentage of neutral residence time of ${\rm 12.5\%}$
with $\sigma_{SE}={\rm 2.4\%}$. These results are consistent with those
determined from the nearest neighbour criterion. The proportion of time spent
in the ionised configuration rises from $\mathrm{79.1\%}$ to
$\mathrm{87.5\%}$. This suggests that the inclusion of zero-point motion
promotes the formation of the ionised state; quantum uncertainty favours
proton transfer.

##  IV Conclusions

As a consequence of the computational expense of the PIMD technique,
especially when using many beads, the simulations presented are limited in
duration to around $\mathrm{10\>ps}$ for some configurations, and rather less
for others. The statistics on the structural and dynamical behaviour are
therefore preliminary. However, it is possible to extract some important
features from these simulations that correspond to intuitive expectation, and
which can be explored further with more extensive calculations.

Our study of small clusters of water and sulphuric acid molecules leads us to
two main conclusions. Firstly, we have demonstrated that molecular dynamics
can reveal features that are not available from knowledge of the geometry
optimised structure at zero temperature. The prime example of this is the
complex behaviour of cluster configuration III-i-1 identified by Re  _et al._
(Re  _et al._ , 1999) and here denoted config H. This configuration has been
regarded as the most stable ionised configuration for the trihydrated
sulphuric acid molecule (Bandy and Ianni, 1998; Re  _et al._ , 1999; Temelso
_et al._ , 2012), but our results indicate that the structure exhibits both
neutral and ionised characteristics at ${\rm 300\>K}$. This conclusion
highlights the limitations of the RRHO approximation(Kathmann  _et al._ ,
2007a, b) for free energy estimation using a single optimised structure.

Secondly, the inclusion of zero-point motion through PIMD simulations of the
sulphuric acid-water system has been shown to produce small but clear
structural distortion at ${\rm 300\>K}$ in a selected number of configurations
when compared with classical dynamics. The mean oxygen-oxygen separation of
hydrogen bonds hb1 and hb2 in the structure shown in Figure 2a is reduced by
${\rm 2-5\%}$. We observe a mild shortening of the hb2 hydrogen bond length,
shown by constructing potentials of mean force for the classical and PIMD
schemes, as illustrated in Figure 6. Furthermore, our results indicate that
zero-point motion brings about a greater propensity for proton transfer in the
O1-H1-O5 substructure of the configuration shown in Figure 4a at ${\rm
300\>K}$.

This conclusion is consistent with the paper by Li  _et al._ (2011) where
quantum nuclear effects on the hydrogen bond are studied, specifically Figure
3 in reference (23) where the OO length is compared with the length of the
projection of the covalent OH bond on the OO vector. The implication is that
the projected covalent bond length is increased by quantum effects when the
hydrogen bond is considered to be strong, as judged by a shift in vibrational
frequency of the covalent OH bond due to the presence of the hydrogen bond.

Our research supports the view that the zero-point effect is most significant
in configurations where proton transfer is intrinsically likely. Classical and
PIMD simulations of the cluster shown in Figure 4a have demonstrated frequent
proton transfer. Using an O-H separation of ${\rm 1.22}\>{\rm\hbox{\AA}}$ as a
threshold for distinguishing the ionised from the neutral state, the cluster
is found to remain neutral ${\rm 20.1\%}$ of the time
(${\rm\sigma_{SE}=2.9\%}$) with classical MD and ${\rm 12.5\%}$
($\sigma_{SE}={\rm 2.4\%}$) according to PIMD. It is possible to infer that
quantum effects have increased the degree of proton transfer. It is expected
that simulations at lower temperatures would increase the significance of the
zero-point effects, making this an avenue for future research. In addition,
since substances such as ammonia and amines are increasingly thought to be
relevant to atmospheric nucleation (Kurtén  _et al._ , 2008; Kirkby  _et al._
, 2011), assessing the importance of zero-point motion in these systems would
also be of interest.

In summary, zero-point motion does affect the structure of small clusters of
sulphuric acid and water, particularly the lengths of hydrogen bonds. At
$\mathrm{300\;K}$, the contribution appears to be most significant for cases
that are intrinsically susceptible to proton transfer.

##  V Acknowledgements

We thank Prof. Angelos Michaelides and his group at UCL for practical advice
and helpful discussions and this work benefited from interactions within the
Thomas Young Centre. SMK was supported by the U.S. Department of Energy,
Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences,
and Biosciences. JLS was supported by the IMPACT scheme at UCL and by the U.S.
Department of Energy, Office of Basic Energy Sciences, Division of Chemical
Sciences, Geosciences, and Biosciences. We are grateful for use of the UCL
Legion High Performance Computing Facility and the resources of the National
Energy Research Scientific Computing Center (NERSC), which is supported by the
U.S. Department of Energy, Office of Science under Contract No. DE-
AC02-05CH11231.

## References

  * Roy (2009) _Geoengineering the climate: science, governance and uncertainty_ (The Royal Society, 2009).
  * Zhang  _et al._ (2012) R. Zhang, A. Khalizov, L. Wang, M. Hu, and W. Xu, Chem. Rev. 112, 1957 (2012).
  * Ford (2004) I. J. Ford, Proc. Inst. Mech. Eng., Part C: J. Mech. Eng. Sci. 218, 883 (2004).
  * Vehkamäki (2006) H. Vehkamäki,  _Classical Nucleation Theory in Multicomponent Systems_ (Springer, 2006).
  * Napari  _et al._ (2010) I. Napari, J. Julin, and H. Vehkamäki, J. Chem. Phys. 133, 154503 (2010).
  * Re  _et al._ (1999) S. Re, Y. Osamura, and K. Morokuma, J. Phys. Chem. A 103, 3535 (1999).
  * Bandy and Ianni (1998) A. R. Bandy and J. C. Ianni, J. Phys. Chem. A 102, 6533 (1998).
  * Ianni and Bandy (1999) J. C. Ianni and A. R. Bandy, J. Mol. Struc-THEOCHEM 497, 19 (1999).
  * Arstila  _et al._ (1998) H. Arstila, K. Laasonen, and A. Laaksonen, J. Phys. Chem. A 108, 1031 (1998).
  * Larson  _et al._ (2000) L. J. Larson, M. Kuno, and F.-M. Tao, J. Chem. Phys. 112, 8830 (2000).
  * Ding  _et al._ (2003) C.-G. Ding, K. Laasonen, and A. Laaksonen, J. Phys. Chem. A 107, 8648 (2003).
  * A. Natsheh  _et al._ (2004) A. A. Natsheh, A. B. Nadykto, K. V. Mikkelsen, F. Yu, and J. Ruuskanen, J. Phys. Chem. A 108, 8914 (2004).
  * Arrouvel  _et al._ (2005) C. Arrouvel, V. Viossat, and C. Minot, J. Mol. Struc-THEOCHEM 718, 71 (2005).
  * Ding and Laasonen (2004) C.-G. Ding and K. Laasonen, Chem. Phys. Lett. 390, 307 (2004).
  * Nadykto  _et al._ (2008) A. B. Nadykto, F. Yu, and J. Herb, Phys. Chem. Chem. Phys. 10, 7073 (2008).
  * Kurtén  _et al._ (2009) T. Kurtén, I. K. Ortega, and H. Vehkamäki, J. Mol. Struc-THEOCHEM 901, 169 (2009).
  * Temelso  _et al._ (2012) B. Temelso, T. E. Morrell, R. M. Shields, M. A. Allodi, E. K. Wood, K. N. Kirschner, T. C. Castonguay, K. A. Archer, and G. C. Shields, J. Phys. Chem. A 116, 2209 (2012).
  * Martin (2004) R. M. Martin,  _Electronic structure: Basic Theory and practical methods_ (Cambridge University Press, 2004).
  * Choe  _et al._ (2007) Y.-K. Choe, E. Tsuchida, and T. Ikeshoji, J. Chem. Phys. 126, 154510 (2007).
  * Anderson  _et al._ (2008) K. E. Anderson, J. I. Siepmann, P. H. McMurry, and J. VandeVondele, J. Am. Chem. Soc. 130, 14144 (2008).
  * Hammerich  _et al._ (2008) A. Hammerich, V. Buch, and F. Mohamed, Chem. Phys. Lett. 460, 423 (2008).
  * Tuckerman (2010) M. E. Tuckerman,  _Statistical Mechanics: Theory and Molecular Simulations_ (Oxford University Press, 2010).
  * Li  _et al._ (2011) X.-Z. Li, B. Walker, and A. Michaelides, Proceedings of the National Academy of Sciences 108, 6369 (2011).
  * Walker and Michaelides (2010) B. Walker and A. Michaelides, J. Chem. Phys. 133, 174306 (2010).
  * Stewart (2007) J. J. P. Stewart, J. Mol. Model. 13, 1173 (2007).
  * Kakizaki  _et al._ (2009) A. Kakizaki, H. Motegi, T. Yoshikawa, T. Takayanagi, M. Shiga, and M. Tachikawa, J. Mol. Struc-THEOCHEM 901, 1 (2009).
  * Sugawara  _et al._ (2011) S. Sugawara, T. Yoshikawa, T. Takayanagi, M. Shiga, and M. Tachikawa, J. Phys. Chem. A 115, 11486 (2011).
  * Feynman and Hibbs (2010) R. P. Feynman and A. R. Hibbs,  _Quantum Mechanics and Path Integrals: Emended Edition (Dover Books on Physics)_ (Dover Publications, 2010).
  * Clark  _et al._ (2005) S. J. Clark, M. D. Segall, and C. J. Pickard, Z. Kristallogr. 220, 567 (2005).
  * Perdew  _et al._ (1996) J. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Thanthiriwatte  _et al._ (2011) K. S. Thanthiriwatte, E. G. Hohenstein, L. A. Burns, and C. D. Sherrill, J. Chem. Theory Compt. 7, 88 (2011).
  * Ireta  _et al._ (2004) J. Ireta, J. Neugebauer, and M. Scheffler, J. Chem. Phys. A 108, 5692 (2004).
  * Haile (1997) J. M. Haile,  _Molecular Dynamics Simulation: Elementary Methods (Wiley Professional)_ (Wiley-Interscience, 1997) p. 512.
  * Hanwell  _et al._ (2012) M. D. Hanwell, D. E. Curtis, D. C. Lonie, T. Vandermeersch, E. Zurek, and G. R. Hutchison, J. Cheminf. 4, 17 (2012).
  * Flyvbjerg and Petersen (1989) H. Flyvbjerg and H. G. Petersen, J. Chem. Phys. 91, 461 (1989).
  * Frenkel and Smit (2001) D. Frenkel and B. Smit,  _Understanding Molecular Simulations: from Algorithms to Application_ , 2nd ed. (Elsevier, 2001).
  * Kumar  _et al._ (2007) R. Kumar, J. R. Schmidt, and J. L. Skinner, J. Chem. Phys. 126, 204107 (2007).
  * Kathmann  _et al._ (2007a) S. Kathmann, G. Schenter, and B. Garrett, Phys. Rev. Lett. 98, 109603 (2007a).
  * Kathmann  _et al._ (2007b) S. Kathmann, G. Schenter, and B. Garrett, J. Phys. Chem. C 111, 4977 (2007b).
  * Kurtén  _et al._ (2008) T. Kurtén, V. Loukonen, H. Vehkamäki, and M. Kulmala, Atmospheric Chemistry and Physics Discussions 8, 7455 (2008).
  * Kirkby  _et al._ (2011) J. Kirkby, J. Curtius, J. Almeida, E. Dunne, J. Duplissy, S. Ehrhart, A. Franchin, S. Gagné, L. Ickes, A. Kürten, A. Kupc, A. Metzger, F. Riccobono, L. Rondo, S. Schobesberger, G. Tsagkogeorgas, D. Wimmer, A. Amorim, F. Bianchi, M. Breitenlechner, A. David, J. Dommen, A. Downard, M. Ehn, R. C. Flagan, S. Haider, A. Hansel, D. Hauser, W. Jud, H. Junninen, F. Kreissl, A. Kvashin, A. Laaksonen, K. Lehtipalo, J. Lima, E. R. Lovejoy, V. Makhmutov, S. Mathot, J. Mikkilä, P. Minginette, S. Mogo, T. Nieminen, A. Onnela, P. Pereira, T. Petäjä, R. Schnitzhofer, J. H. Seinfeld, M. Sipilä, Y. Stozhkov, F. Stratmann, A. Tomé, J. Vanhanen, Y. Viisanen, A. Vrtala, P. E. Wagner, H. Walther, E. Weingartner, H. Wex, P. M. Winkler, K. S. Carslaw, D. R. Worsnop, U. Baltensperger, and M. Kulmala, Nature 476, 429 (2011).


# Orbital Motion of a test particle around a Schwarzschild’s Black Hole in
STVG gravity

Devansh Shukla
Sardar Vallabhbhai National Institute of Technology, Surat
<EMAIL_ADDRESS>
&Abhay Menon A
Sardar Vallabhbhai National Institute of Technology, Surat
<EMAIL_ADDRESS>
&Kamlesh Pathak
Sardar Vallabhbhai National Institute of Technology, Surat
<EMAIL_ADDRESS>

###### Abstract

In this article, we have examined the existence of a static spherically
symmetric solution in the Scalar Tensor Vector Gravity (STVG) and investigated
its horizon distances to develop boundary limitations for our test particle.
We have computed the Kretschmann invariant of the metric to study the
singularities and verify that it reduces to general relativity’s Kretschmann
invariant as $\alpha\rightarrow 0$. Further, we investigated the orbital
motion of a time-like and light-like test particle around the static solution
by developing an effective potential and the radius of the innermost stable
circular orbit(ISCO).

##  1 Introduction

####

The Scalar Tensor Vector Gravity theory, or STVG for short, is a recent
modified theory of gravity developed by John W. Moffat in 2005. The theory
allows the gravitational constant $G$, a vector field coupling $\omega$ and
the vector field mass $\mu$ to vary with spacetime by being scalar fields[1].
The effect is that there is a Yukawa-like modification of the force,
essentially predicting that the gravity is attractive and stronger than the
Newtonian prediction at large distances and repulsive at small distances.

Although Einstein’s General Relativity(GR) has been quite successful, it falls
short in explaining the cosmic acceleration of the universe, galaxy rotation
curves and galaxy velocities without the inclusion of the mysterious "dark"
matter, which has to be about six times the baryonic matter to account for the
disparity[2, 3]. These arguments provide sufficient motivation to investigate
modified gravity theories. In this article, we will study one such modified
theory called STVG, which has been able to explain the solar system
observations[1], the rotation curves of galaxies[4, 5] and the dynamics of
galactic clusters[6] without the use of "dark" content.

Throughout the years, there has been excitement in the field of black hole
physics in figuring out geodesics for test particles, whether they be neutral
or charged, massive or massless, they help to study the nature of the
spacetime in question. In order to investigate the geodesic motion, we utilize
the parameter called innermost stable circular orbit(ISCO). The ISCO can give
us knowledge of the inner edge of the accretion disk, the spacetime geometry
and the gravitational waves from the binary system[7, 8].

As J.W. Moffat states in his paper[9], SVTG, nonsymmetric gravity theory(NGT),
metric-skew-tensor gravity (MSTG) have an essential feature that the modified
acceleration law for weak gravitational fields has a repulsive Yukawa force,
causing effective mass and coupling to matter to vary with distance rather
than being attractive everywhere. He computed the line elements for
Schwarzschild, Reissner-Nordström, and the Kerr solution. In this article, we
will utilize the line elements to compute and study the nature of the
innermost stable circular orbits in the Schwarzschild’s solution.

The paper is organized as follows. In section 2, we provide a brief
explanation of the Scalar-Tensor-Vector gravity and its field equations. In
section 3, we develop the static spherically symmetric solution’s line element
and examine its horizon distance(s). Next, we turn our attention to the
innermost stable circular orbit for time-like and light-like particles in
section 4. Finally, in section 5, we summarize the results.

We use the units $c=G_{N}=k_{B}=\hbar=1$ and the sign convention $(-,+,+,+)$
throughout the paper. The metric in STVG assumes that $G=G_{N}(1+\alpha)$,
where $\alpha$ is a real non-negative constant and $G_{N}$ is the Newtonian
gravitational constant.

##  2 The ansatz model

####

In order to obtain the ansatz model for STVG, we assume $G(x)$, $\mu(x)$ and
$\omega(x)$ to be three scalar fields, with associated potential functions
$V(G)$, $V(\mu)$ and $V(\omega)$ respectively.

Consider an anti-symmetric field $\tensor{\mathcal{B}}{{}_{\mu}{}_{\nu}}$ is
formed out of the vector field $\phi_{\mu}$ such that:

$\tensor{\mathcal{B}}{{}_{\mu}{}_{\nu}}=\partial_{\mu}\phi_{\mu}-\partial_{\nu}\phi_{\mu}$
(1)

The total action, with $A_{M}$ as the matter action, reads:

$A=A_{g}+A_{\phi}+A_{S}+A_{M}$ (2)

where

$A_{g}=\frac{1}{16\pi}\int d^{4}x\sqrt{-g}\left[\dfrac{R}{G}+2\Lambda\right]$
(3) $A_{\phi}=-\int
d^{4}x\sqrt{-g}\omega\left[\dfrac{1}{4}\tensor{\mathcal{B}}{{}_{\mu}{}_{\nu}}\tensor{\mathcal{B}}{{}^{\mu}{}^{\nu}}+V(\phi)\right]$
(4) $\displaystyle\begin{aligned} A_{S}=&\int
d^{4}x\sqrt{-g}\left[\dfrac{1}{G^{3}}\left(\frac{1}{2}\tensor{g}{{}^{\mu}{}^{\nu}}\tensor{\nabla}{{}_{\mu}}G\tensor{\nabla}{{}_{\nu}}G-V(G)\right)+\dfrac{1}{G}\left(\frac{1}{2}\tensor{g}{{}^{\mu}{}^{\nu}}\tensor{\nabla}{{}_{\mu}}\omega\tensor{\nabla}{{}_{\nu}}\omega-V(\omega)\right)\right.\\\
&\left.+\dfrac{1}{\mu^{2}G}\left(\frac{1}{2}\tensor{g}{{}^{\mu}{}^{\nu}}\tensor{\nabla}{{}_{\mu}}\mu\tensor{\nabla}{{}_{\nu}}\mu-V(\mu)\right)\right]\end{aligned}$
(5)

Here $\tensor{\nabla}{{}_{\mu}}$ denotes the covariant derivate with respect
to the metric $\tensor{g}{{}_{\mu}{}_{\nu}}$

The total energy-momentum tensor is given by:

$\tensor{T}{{}_{\mu}{}_{\nu}}=\tensor{T}{{}_{M}{}_{\mu}{}_{\nu}}+\tensor{T}{{}_{\phi}{}_{\mu}{}_{\nu}}+\tensor{T}{{}_{S}{}_{\mu}{}_{\nu}}$
(6)

where $\tensor{T}{{}_{M}{}_{\mu}{}_{\nu}}$ denotes the energy-momentum
contribution by ordinary matter; $\tensor{T}{{}_{\phi}{}_{\mu}{}_{\nu}}$ is
the contribution from the vector field and
$\tensor{T}{{}_{S}{}_{\mu}{}_{\nu}}$ is by the scalar fields $G$, $\mu$ and
$\omega$.

$\tensor{T}{{}_{M}{}_{\mu}{}_{\nu}}=-\dfrac{2}{\sqrt{-g}}\dfrac{\delta
A_{M}}{\delta\tensor{g}{{}^{\mu}{}^{\nu}}};\quad\tensor{T}{{}_{\phi}{}_{\mu}{}_{\nu}}=-\dfrac{2}{\sqrt{-g}}\dfrac{\delta
A_{\phi}}{\delta\tensor{g}{{}^{\mu}{}^{\nu}}};\quad\tensor{T}{{}_{S}{}_{\mu}{}_{\nu}}=-\dfrac{2}{\sqrt{-g}}\dfrac{\delta
A_{S}}{\delta\tensor{g}{{}^{\mu}{}^{\nu}}}$ (7)

We use the variational principle to compute the field equation by varying the
Eq. (2) against the metric $\tensor{g}{{}^{\mu}{}^{\nu}}$ to compute the field
equation

$\tensor{\mathcal{G}}{{}_{\mu}{}_{\nu}}-\tensor{g}{{}_{\mu}{}_{\nu}}\Lambda+Q_{\mu\nu}=8\pi
G\tensor{T}{{}_{\mu}{}_{\nu}}$ (8)

where
$Q_{\mu\nu}=G\left(\tensor{g}{{}_{\mu}{}_{\nu}}\tensor{\nabla}{{}^{\alpha}}\tensor{\nabla}{{}_{\alpha}}\dfrac{1}{G(x)}-\tensor{\nabla}{{}_{\mu}}\tensor{\nabla}{{}_{\nu}}\dfrac{1}{G(x)}\right)$

This term is identical to the boundary term in the Brans-Dicke theory and STVG
reduces to it under suitable conditions[10].

##  3 Static spherically symmetrical solution

####

In order to obtain the Schwarzschild’s line element, we take an ansatz
solution for the most general spherically symmetric metric:

$ds^{2}=-Z(r)dt^{2}+Z^{-1}(r)dr^{2}+r^{2}d\Omega^{2}$ (9)

J.W. Moffat[11] computed the value of $Z(r_{h})$ by assuming $\Lambda$ and
$V(\phi)$ vanish while keeping the coupling constant $\omega=1$:

$Z(r)=1-\dfrac{2GM}{r}+\dfrac{\alpha G_{N}GM^{2}}{r^{2}}$

Since, $G_{N}=G(1+\alpha)$, $Z(r)$ becomes:

$Z(r)=1-\dfrac{2(1+\alpha)M}{r}+\dfrac{\alpha(1+\alpha)M^{2}}{r^{2}}$ (10)

The line element reads:

$ds^{2}=-\left(1-\dfrac{2(1+\alpha)M}{r}+\dfrac{\alpha(1+\alpha)M^{2}}{r^{2}}\right)dt^{2}+\dfrac{dr^{2}}{\left(1-\dfrac{2(1+\alpha)M}{r}+\dfrac{\alpha(1+\alpha)M^{2}}{r^{2}}\right)}+r^{2}d\Omega^{2}$
(11)

The Kretschmann invariant, defined as the contraction of the Riemannian
tensor,
$K=\tensor{R}{{}_{\alpha}{}_{\beta}{}_{\mu}{}_{\nu}}\tensor{R}{{}^{\alpha}{}^{\beta}{}^{\mu}{}^{\nu}}$
reads:

$K=\dfrac{48M^{2}(1+\alpha)^{2}r^{2}-96M^{3}\alpha(1+\alpha)^{2}r+56M^{4}\alpha^{2}(1+\alpha)^{2}}{r^{8}}$
(12)

The Kretschmann invariant reduces to general relativity’s case for
$\alpha\rightarrow 0$, and as $r\rightarrow 0$, the invariant becomes
indeterminant, implying a singularity at $r=0$.

The modified line element Eq. (11) will describe a black hole if there is a
locatable horizon present where $Z(r)$ vanishes and $Z^{\prime}(r_{h})>0$ for
positive surface gravity[12]:

$r_{\pm}=M\left\\{\left(1+\alpha\right)\pm\sqrt{1+\alpha}\right\\}$ (13)

where $r_{+}$ and $r_{-}$ represents the outer(event horizon) and inner
horizon(Cauchy horizon) respectively.

Let’s consider the case where $\alpha=0$, the horizons exists at $r_{+}=2M$
and $r_{-}=0$, which are identical to the general relativity’s predictions.

(a) $MZ^{\prime}(r_{+})$ vs $\alpha$ (b) $MZ^{\prime}(r_{-})$ vs $\alpha$

Figure 1: The plots for $Z^{\prime}(r)$ with respect to $\alpha$ for both
solutions $r_{\pm}$

$MZ^{\prime}(r_{+})$ remains positive against $\alpha$, indicating positive
surface gravity

The static spherically symmetric solution, given by Eq. (11), describes a
physical black hole because the root $r_{+}$ is real and positive while
$MZ^{\prime}(r_{+})$ remains positive against $\alpha$ (Fig. 1a).

##  4 Test particle trajectory

Let us assume a test particle of a non-negative mass $\mu$ in the equatorial
frame ($\theta=\pi/2$). Its line element reads:

$ds^{2}=\tensor{g}{{}_{\alpha}{}_{\beta}}dx^{\alpha}dx^{\beta}=-Z(r)dt^{2}+Z^{-1}(r)dr^{2}+r^{2}d\phi^{2}=\epsilon$
(14)

where $\epsilon=-1,0,1$ for time-like, null-like, and space-like trajectories
respectively

The Lagrangian of a free particle in a potential-less spacetime reads:

$\mathcal{L}=\mu\sqrt{\tensor{g}{{}_{\alpha}{}_{\beta}}\dot{x}^{\alpha}\dot{x}^{\beta}}$
(15)

Since Eq. (15) is independent of $t$ and $\phi$, the first integrals of motion
obtained using the Euler-Lagrange equation reads:

$\mu Z(r)\dot{t}=E;\quad\mu r^{2}\dot{\phi}=L$ (16)

where $E$ and $L$ represent the energy and the angular momentum of the test
particle respectively.

Substituting Eq. (16) into the spacetime interval Eq. (14):

$-\dfrac{E^{2}}{\mu^{2}Z(r)}+\dfrac{\dot{r}^{2}}{Z(r)}+\dfrac{L^{2}}{\mu^{2}r^{2}}=\epsilon$
(17) $\dot{r}^{2}+V^{2}_{\text{eff}}=\dfrac{E^{2}}{\mu^{2}}$ (18)

with the effective potential represented by
$V_{\text{eff}}=\sqrt{\left(\dfrac{L^{2}}{\mu^{2}r^{2}}-\epsilon\right)\left(1-\dfrac{2GM}{r}+\dfrac{\alpha
G_{N}M^{2}}{r^{2}}\right)}$

The extremum values of $V_{\text{eff}}$ determine the stable and unstable
orbits of the test particle, with a circular orbit defined by $\dot{r}=0$ and
$V_{\text{eff},\;r}=0$.

Solving for $\tensor{\partial}{{}_{r}}V_{\text{eff}}=0$, we obtain:

$\dfrac{L^{2}}{\mu^{2}}=\frac{Mr^{2}\epsilon(1+\alpha)(\alpha
M-r)}{r^{2}-3(1+\alpha)Mr+2\alpha(1+\alpha)M^{2}}$ (19)

Restricting ourselves to only circular orbits ($\dot{r}=0$), the energy
equation reads:

$\dfrac{E^{2}}{\mu^{2}}=-\frac{\epsilon\left(r^{2}-2(1+\alpha)Mr+\alpha(1+\alpha)M^{2}\right)^{2}}{r^{2}\left(r^{2}-3(1+\alpha)Mr+2\alpha(1+\alpha)M^{2}\right)}$
(20)

###  4.1 Time-like orbits: $\epsilon=-1$

####

Choosing $\epsilon=-1$, we restrict ourselves to time-like orbits, with
angular momentum and energy defined by Eq. (19) and Eq. (20).

The existence of a real-valued angular momentum implies two additional
conditions for the existence of motion.

$r>\alpha M=r^{*};\quad r^{2}-3(1+\alpha)Mr+2\alpha(1+\alpha)M^{2}>0$ (21)

Solving Eq. (21), we obtain two additional roots:

$\widetilde{r}_{\pm}=\frac{M}{2}\left(3(1+\alpha)\pm\sqrt{\alpha^{2}+10\alpha+9}\right)$
(22)

Implying from Fig. 2, we can see that the $r$’s follow linear functions with
respect to $\alpha$, with $r_{-}$ and $\widetilde{r}_{+}$ as the extremum
values:

$r_{-}<\widetilde{r}_{-}<r^{*}<r_{+}<\widetilde{r}_{+}$ (23)

In the case of time-like circular orbits, by choosing the extremum value, the
motion occurs in the region $r>\widetilde{r}_{+}$.

To obtain the radius of the circular orbit, we substitute the value of the
angular momentum Eq. (19) into the effective potential and get its extremum
values by solving for $\partial_{r}V_{\text{eff}}(r)=0$:

$\left(r^{2}+GM(\alpha
G_{N}M-2r)\right)\left(r^{3}-GM\left(4\alpha^{2}G_{N}^{2}M^{2}-9\alpha
G_{N}Mr+6r^{2}\right)\right)=0$ (24) $\displaystyle\begin{aligned}
\dfrac{r_{isco}}{M}=2(1+\alpha)&+\left({\alpha\left(\alpha+1\right)\sqrt{(\alpha+5)}+\left(\alpha^{3}+8\alpha^{2}+15\alpha+8\right)}\right)^{1/3}\\\
&+\frac{\alpha^{2}+5\alpha+4}{\left({\alpha\left(\alpha+1\right)\sqrt{(\alpha+5)}+\left(\alpha^{3}+8\alpha^{2}+15\alpha+8\right)}\right)^{1/3}}\end{aligned}$
(25)

Evaluating Eq. (25) at $\alpha=0$, we obtain $r_{isco}=6M$, which agrees with
the predictions of general relativity and as expected
$r_{isco}>\widetilde{r}_{+}$, as shown in Fig. 2.

Figure 2: Plot of $\widetilde{r}_{\pm}$, $r_{\pm}$, $r^{*}$ and $r_{isco}$ per
unit $M$ vs $\alpha$

At $\alpha=0$, $r_{isco}=6M$, $\widetilde{r}_{+}=3M$, $r_{+}=2M$ and
$\widetilde{r}_{-}=r_{-}=r^{*}=0$

(a) Energy per unit mass as a function of $\dfrac{r}{M}$(for $\dfrac{r}{M}>1$)

(b) Angular momentum per unit mass as a function of $\dfrac{r}{M}$

Figure 3: Plot of (a) $\dfrac{E}{\mu}$ and (b) $\dfrac{L}{\mu M}$ as a
function of $\dfrac{r}{M}$, for $\alpha=0,0.5,1$

###  4.2 Light-like orbits: $\epsilon=0$

Choosing $\epsilon=0$, we restrict ourselves to light-like particles(null
trajectories).

Differentiating the effective potential and computing its extremum values
while assuming a non-zero angular momentum, we obtain:

$r^{2}+(\alpha+1)M(2\alpha M-3r)=0$ (26)

The real and positive roots of Eq. (26) will define the innermost stable
orbits for light-like particles:

$r_{isco}^{\pm}=\frac{M}{2}\left(3(\alpha+1)\pm\sqrt{\alpha^{2}+10\alpha+9}\right)$
(27)

Evaluating Eq. (27) at $\alpha=0$, we obtain the $r^{+}_{isco}=3M$ and
$r^{-}_{isco}=0$, we can safely neglect the $r^{-}_{isco}$ since no motion is
not possible below $\widetilde{r}_{+}$. Solely focusing on
$r_{isco}=r^{+}_{isco}=3M$, we can see it is just as predicted by general
relativity.

Figure 4: Plot of light-like $r_{+}$, $r_{-}$ and $r_{isco}$ per unit $M$ as a
function of $\alpha$

At $\alpha=0$, $r_{isco}=3M$, $r_{+}=2M$, $r_{-}=0$

As shown in Fig. 4, the $r_{isco}$, $r_{+}$ and $r_{-}$ are all proportional
to $\alpha$ and $r_{isco}>r_{+}>r_{-}$

##  5 Conclusion

####

In this article, we have developed the Schwarzschild’s line element in the
Scalar-Tensor-Vector gravity and examined it’s Kretschmann’s invariant, given
by Eq. (12). The Kretschmann’s invariant is indeterminable at $r=0$ implying a
singularity. Analyzing the horizon distances, we obtain two values: the event
horizon ($r_{+}$) and the Cauchy horizon($r_{-}$). As expected, the Cauchy
horizon is unstable, similar to the case of the Reissner-Nordström solution.

The orbit equation of a neutral particle around the black hole solution is
given by Eq. (18) with the angular momentum defined by Eq. (19) and the energy
for the circular orbit by Eq. (20). We study the nature of energy momentum and
energy for time-like particles in Fig. 3 for multiple values of $\alpha$.

We see that for $\epsilon=0$(null-like) and $\epsilon=-1$(time-like), the
innermost stable orbit is a strictly increasing function of $\alpha$. This
means that, in both cases, we see an outward shift of $r_{isco}$ as the
strength of the gravitational field increases. This is expected since stronger
gravitational fields cause greater perturbations closer to the source. We
notice that as $\alpha\rightarrow 0$, the $r_{isco}$ reaches the GR
predictions of $6M$ and $3M$ for time-like and null-like trajectories,
respectively.

## References

  * [1] J. W. Moffat.  Scalar and vector field constraints, deflection of light and lensing in modified gravity (mog), 2014.
  * [2] César Augusto Zen Vasconcellos.  Topics on Strong Gravity.  WORLD SCIENTIFIC, jan 2020.
  * [3] Wasley S. Krogdahl.  A critique of general relativity, 2007.
  * [4] J. W. Moffat and S. Rahvar.  The MOG weak field approximation and observational test of galaxy rotation curves.  Monthly Notices of the Royal Astronomical Society, 436(2):1439–1451, sep 2013.
  * [5] J. W. Moffat and S. Rahvar.  The MOG weak field approximation – II. observational test of chandra x-ray clusters.  Monthly Notices of the Royal Astronomical Society, 441(4):3724–3732, jun 2014.
  * [6] J. W. Moffat and V. T. Toth.  Rotational velocity curves in the milky way as a test of modified gravity.  Physical Review D, 91(4), feb 2015.
  * [7] Don N. Page and Kip S. Thorne.  Disk-accretion onto a black hole. time-averaged structure of accretion disk.  The Astrophysical Journal, 191:499–506, July 1974.
  * [8] Robert M Wald.  General relativity.  University of Chicago press, 2010.
  * [9] J W Moffat.  Scalar–tensor–vector gravity theory.  Journal of Cosmology and Astroparticle Physics, 2006(03):004–004, mar 2006.
  * [10] C. Brans and R. H. Dicke.  Mach’s principle and a relativistic theory of gravitation.  Phys. Rev., 124:925–935, Nov 1961.
  * [11] J. W. Moffat.  Black holes in modified gravity (MOG).  The European Physical Journal C, 75(4), apr 2015.
  * [12] Marco Calzà, Massimiliano Rinaldi, and Lorenzo Sebastiani.  A special class of solutions in f(r)-gravity.  The European Physical Journal C, 78(3), mar 2018.
  * [13] Timothy Clifton, Pedro G. Ferreira, Antonio Padilla, and Constantinos Skordis.  Modified gravity and cosmology.  Physics Reports, 513(1-3):1–189, mar 2012.
  * [14] Luciano Rezzolla and Alexander Zhidenko.  New parametrization for spherically symmetric black holes in metric theories of gravity.  Physical Review D, 90(8), oct 2014.
  * [15] M. Sharif and Misbah Shahzadi.  Neutral particle motion around a schwarzschild black hole in modified gravity.  Journal of Experimental and Theoretical Physics, 127:491–502, 9 2018.
  * [16] Abhay Menon A.  Particle motion under schwarzschild gravity: A meta-study of different theories of gravity.  Master’s thesis, Sardar Vallabhbhai National Institute of Technology, Surat, 2020.
  * [17] Priti Mishra and Tejinder P. Singh.  Fourth order gravity, scalar-tensor-vector gravity, and galaxy rotation curves.  Physical Review D, 88(10), nov 2013.
  * [18] Xin-Chang Cai and Yan-Gang Miao.  High-dimensional schwarzschild black holes in scalar–tensor–vector gravity theory.  The European Physical Journal C, 81(6), jun 2021.
  * [19] Hyun-Chul Lee and Yong-Jin Han.  Innermost stable circular orbit of kerr-mog black hole, Sep 2017.