# A Simple Boundary Condition Regularization Strategy for Image Velocimetry
Based Pressure Field Reconstruction

Connor Pryce University of Waterloo, Dept. of Mechanical and Mechatronics
Engineering, Waterloo, ON, Canada  Lanyu Li University of Waterloo, Dept. of
Mechanical and Mechatronics Engineering, Waterloo, ON, Canada  Jared P.
Whitehead Brigham Young University, Mathematical Dept., Provo, UT, USA  Zhao
Pan To whom correspondence may be addressed<EMAIL_ADDRESS>University
of Waterloo, Dept. of Mechanical and Mechatronics Engineering, Waterloo, ON,
Canada

###### Abstract

We propose a simple boundary condition regularization strategy to reduce error
propagation in pressure field reconstruction from corrupted image velocimetry
data. The core idea is to replace the canonical Neumann boundary conditions
with Dirichlet ones obtained by integrating the tangential part of the
pressure gradient along the boundaries. Rigorous analysis and numerical
experiments justify the effectiveness of this regularization.

##  1 Background and Statement of the Problem

Solving the Pressure Poisson Equation (PPE) is one popular approach to
pressure field reconstruction from image velocimetry data (Van Oudheusden,,
2013). Techniques like Particle Image Velocimetry (PIV) and Lagrangian
Particle Tracking (LPT) provide a non-intrusive estimate of the pressure
gradient $\nabla
p=\bm{g}(\bm{u})=-\bm{u}_{t}-(\bm{u}\cdot\nabla)\bm{u}+{\text{Re}}^{-1}\nabla^{2}\bm{u}$
in the domain and on the boundary via the momentum equations, where $\bm{g}$
is a function of the velocity field $\bm{u}$, and Re is the Reynolds number.

Applying the divergence to $\nabla p$ yields the PPE:

$\nabla^{2}p=f(\bm{u})=-\nabla\cdot\left(\frac{\partial\bm{u}}{\partial
t}+(\bm{u}\cdot\nabla)\bm{u}-\frac{1}{\text{Re}}\nabla^{2}\bm{u}\right),$ (1)

where $f$ is the data of the Poisson equation. To solve (1), boundary
conditions are required. Commonly used boundary conditions are i) Dirichlet
boundary conditions

$p=h\quad\text{on}~{}\partial\Omega,$ (2)

which define the value of the pressure on the boundaries and are often
measured by pressure transducers and/or estimated using Bernoulli’s principle
in irrotational regions of the domain (De Kat and Van Oudheusden,, 2012); as
well as ii) Neumann boundary conditions

$\bm{n}\cdot\nabla p=g_{n}\quad\text{on}~{}\partial\Omega,$ (3)

which define the pressure gradient normal to the boundary ($\bm{n}$ is the
outward pointing unit normal vector on the boundary $\partial\Omega$), and are
often derived from the momentum equations.

While the Poisson equation benefits from the well-posedness of elliptic
equations, superior numerical stability, ease of implementation, and high
computational efficiency, experimentation (Charonko et al., (2010); Sperotto
et al., (2022); Zhang et al., (2022)) has shown that pressure field
reconstruction on a domain with long Neumann boundaries may suffer from
excessive error propagation due to contaminated image velocimetry data. This
observation is expected, and supported analytically for two reasons: i)
locally, the error in the data near a Neumann boundary is amplified and
‘diffuses’ further into the interior of the domain (Faiella et al.,, 2021);
and ii) globally, long Neumann boundaries result in a larger Poincare constant
for the Laplacian operator and thus the error in the data may be further
intensified when compared to a domain with Dirichlet conditions of the same
length. For these reasons, and despite being available everywhere, Neumann
boundary conditions should be avoided if possible (Pan et al.,, 2016).
Alternatively, Dirichlet conditions are favourable for taming error
propagation, but they are not always easily accessible in practice.

This renders a fundamental dilemma that inhibits the performance of the PPE
for pressure field reconstruction from noisy velocimetry data. In an attempt
to break this bottleneck, one idea is to replace the Neumann boundary
conditions with derived Dirichlet ones. In this letter, we introduce a simple
regularization strategy of this type to address the boundary condition
dilemma, and provide a rigorous analysis and validation to show that this
specific regularization can improve the quality of the pressure field
reconstruction.

##  2 Boundary Condition Regularization

In addition to the natural choice of the canonical Neumann boundary conditions
(cNBC) in (3), the momentum equations also provide the pressure gradient
tangential to the boundary:

$\bm{\tau}\cdot\nabla p=g_{\tau}\quad\text{on}~{}\partial\Omega,$ (4)

where ${g}_{\tau}$ is the pressure gradient projected on $\bm{\tau}$, the unit
vector tangential to the boundary in the counter-clockwise direction.
Assigning a reference pressure $p_{0}$ at a location $\xi_{0}$ on the boundary
and integrating (4) directly gives the value of the pressure which is a
derived Dirichlet boundary condition (dDBC):

$p(\xi)=h_{d}=\int_{\xi_{0}}^{\xi}\bm{\tau}\cdot\nabla
pdS+p_{0}=\int_{\xi_{0}}^{\xi}g_{\tau}dS+p_{0}\quad\text{on}~{}\partial\Omega.$
(5)

Using this dDBC to replace the canonical Neumann boundary conditions for the
Poisson equation yields a Dirichlet problem.

##  3 Analysis and Error Estimate

We consider a simple example to see how the regularization strategy proposed
in Sect. 2 can improve the performance of the pressure field reconstruction
based on a PPE with only a single point Dirichlet condition provided on the
boundary. The error-contaminated problem for a domain with cNBC can be
described as

$\displaystyle\nabla^{2}\tilde{p}$ $\displaystyle=\tilde{f}\quad$
$\displaystyle\text{in}~{}\Omega$ (6a)
$\displaystyle\bm{n}\cdot\nabla\tilde{p}$ $\displaystyle=\tilde{g}_{n}\quad$
$\displaystyle\text{on}~{}\partial\Omega,$ (6b)

where $\tilde{p}=p+\epsilon_{p}$, $\tilde{f}=f+\epsilon_{f}$, and
$\tilde{g}_{n}=g_{n}+\epsilon_{g_{n}}$ are the contaminated pressure, data,
and cNBC, respectively; and $\epsilon_{p}$, $\epsilon_{f}$, and
$\epsilon_{g_{n}}$ are the corresponding error. Comparing (6) with (3) and (1)
leads to a new Poisson equation with respect to the error in the pressure
field:

$\displaystyle\nabla^{2}\epsilon_{p}$ $\displaystyle=\epsilon_{f}\quad$
$\displaystyle\text{in}~{}\Omega$ (7a)
$\displaystyle\bm{n}\cdot\nabla\epsilon_{p}$
$\displaystyle=\epsilon_{g_{n}}\quad$
$\displaystyle\text{on}~{}\partial\Omega,$ (7b)

For the error propagation problem of (7), if we measure $\epsilon_{p}$ with a
space-averaged error
($\|\epsilon_{p}\|_{L^{2}{\Omega}}=\sqrt{\int\epsilon_{p}^{2}d\Omega/|\Omega|}$),
the techniques in Pan et al., (2016) provide an estimate for the error in the
pressure field which we present here for convenience:

$\|\epsilon_{p}\|_{L^{2}(\Omega)}\leq
C_{N}\|\epsilon_{f}\|_{L^{2}(\Omega)}+D_{N}\|\epsilon_{g_{n}}\|_{L^{2}(\partial\Omega)},$
(8)

where $D_{N}=\sqrt{C_{N}C_{NB}{|\partial\Omega|}/{|\Omega|}}$ is the
coefficient corresponding to the error contribution from the boundary. $C_{N}$
and $C_{NB}$ are the Poincare constants for the problem, independent of the
image velocimetry technique, solver implementation and spatial resolution.
These constants are solely determined by the geometry (i.e., shape, size, and
dimension) of the domain and the setup of the boundary conditions.
$|\partial\Omega|$, $|\Omega|$ represent the size (length, area, or volume,
depending on the dimension of the problem) of the boundary and the domain.

If we apply the boundary regularization strategy described in Sect. 2, the
error contaminated dDBC and the corresponding error are:

$\displaystyle\tilde{p}$ $\displaystyle=\tilde{h}_{d}\quad$
$\displaystyle\text{on}~{}\partial\Omega$ (9a) $\displaystyle\epsilon_{p}$
$\displaystyle=\epsilon_{h_{d}}\quad$
$\displaystyle\text{on}~{}\partial\Omega,$ (9b)

where $\epsilon_{h_{d}}$ is the error in the dDBC. Replacing (6b) with (9a),
the Poisson problem (6a) is regularized, and the corresponding error
propagation can be modeled by replacing (7b) with (9b) in (7). Similar
techniques in Pan et al., (2016) provide an error estimate of the recovered
pressure for a Dirichlet problem:

$\|\epsilon_{p}\|_{L^{2}(\Omega)}\leq
C_{D}\|\epsilon_{f}\|_{L^{2}(\Omega)}+\|\epsilon_{h_{d}}\|_{L^{\infty}(\partial\Omega)},$
(10)

where $C_{D}$ is the corresponding Poincare constant. Noting that
$|\epsilon_{h_{d}}|\leq\int_{\xi_{0}}^{\xi}|\epsilon_{g_{\tau}}|dS\leq|\partial\Omega|\|\epsilon_{g_{\tau}}\|_{L^{2}(\partial\Omega)}$,
by the Cauchy-Schwarz inequality,
$\|\epsilon_{h_{d}}\|_{L^{\infty}(\partial\Omega)}$ in (10) can be estimated
as

$\|\epsilon_{h_{d}}\|_{L^{\infty}(\partial\Omega)}\leq|\partial\Omega|\|\epsilon_{g_{\tau}}\|_{L^{2}(\partial\Omega)}.$
(11)

Combining (10) and (11), we arrive at

$\|\epsilon_{p}\|_{L^{2}(\Omega)}\leq
C_{D}\|\epsilon_{f}\|_{L^{2}(\Omega)}+|\partial\Omega|\|\epsilon_{g_{\tau}}\|_{L^{2}(\partial\Omega)},$
(12)

which is an estimate of the error in the reconstructed pressure field when the
regularization described in (5) is practiced. For a domain with a given
geometry, the Poincare constant ($C_{D}$ in (12)) is almost always smaller
than that of the original Neumann problem ($C_{N}$ in (8). This suggests that
the proposed regularization in (5) can improve the accuracy of the pressure
reconstruction, especially for large domains.

Letting $\Omega=L\times L$ be a two-dimensional domain, a concrete example is
presented to demonstrate the effectiveness of the proposed regularization. For
the problem equipped with cNBC, (8) becomes:

$\|\epsilon_{p}\|_{L^{2}(\Omega)}\leq\frac{L^{2}}{\pi^{2}}\|\epsilon_{f}\|_{L^{2}(\Omega)}+\frac{4}{\pi^{3/2}}L\|\epsilon_{g_{n}}\|_{L^{2}(\partial\Omega)},$
(13)

where $C_{N}=L^{2}/\pi^{2}$, $C_{NB}=4L$, $|\Omega|=L^{2}$, and
$|\partial\Omega|=4L$. For the regularized problem on the same domain with
dDBC, (12) turns into

$\|\epsilon_{p}\|_{L^{2}(\Omega)}\leq\frac{L^{2}}{2\pi^{2}}\|\epsilon_{f}\|_{L^{2}(\Omega)}+4L\|\epsilon_{g_{\tau}}\|_{L^{2}(\partial\Omega)},$
(14)

where $C_{D}=L^{2}/2\pi^{2}$. Comparing (13) against (14), it is obvious that
the proposed regularization can produce a more accurate pressure
reconstruction. For a square domain, $C_{D}$ is half of $C_{N}$, and an
elongated rectangular domain reduces this factor again, further taming error
propagation. This is achieved at the potential cost of slightly increasing the
impact of the error in the pressure gradient on the boundary (the constant
$4L/\pi^{3/2}$ in front of $\|\epsilon_{g_{n}}\|_{L^{2}(\partial\Omega)}$ in
(13) is smaller than the constant $4L$ in front of
$\|\epsilon_{g_{\tau}}\|_{L^{2}(\partial\Omega)}$ in (14)). Thus, for
experiments with high uncertainty in a large domain, the regularization is
expected to be effective in reducing the error in the reconstructed pressure
field.

##  4 Curl-Free Regularization on a Closed Boundary

If the regularization proposed in (5) is employed on a closed boundary, it may
lead to a non-physical Dirichlet condition when the error on the boundary
($\epsilon_{g_{\tau}}$) is not curl-free. That is, after the integration
indicated in (5), the pressure value at $\xi=\xi_{0}$ may not be continuous
($p(\xi_{0}^{+})\neq p(\xi_{0}^{-})$). If integrated along another direction
starting from the same reference point and reference pressure, a different
pressure along the boundary is expected. This should not be the case since
$\nabla p$ is curl-free (including the path along the boundary), and
integration on $\nabla p$ should be path-independent. To resolve this issue,
we construct a continuous function $p(\xi)$ on the boundary such that $\nabla
p(\xi)$ is an interpolant or regression of $\tilde{g}_{\tau}$, subject to the
constraint $p(\xi_{0}^{+})=p(\xi_{0}^{-})=p_{0}$.

In 2D, one way to explicitly achieve such a Dirichlet boundary is by linear
interpolation:

$p(\xi)=\tilde{h}_{d}=\underbrace{(1-\theta)\int_{\xi_{0}}^{\xi}\bm{\tau}\cdot\nabla\tilde{p}dS}_{\text{counter~{}clockwise~{}int.}}+\underbrace{\theta\int_{\xi_{0}+|\partial\Omega|}^{\xi}\bm{\tau}\cdot\nabla\tilde{p}dS}_{\text{clockwise~{}int.}}+p_{0},$
(15)

where $\theta=(\xi-\xi_{0})/|\partial\Omega|$, and
$\xi_{0}\leq\xi\leq\xi_{0}+|\partial\Omega|$. Equation (15) has an apparent
physical meaning and can be thought of as a weighted average of integration of
$\nabla\tilde{p}$ along two different directions, with $\theta$ being a
weighting parameter. For $\theta=0~{}\text{or}~{}1$, the counter-clockwise or
the clockwise integral starting from $\xi_{0}$ has not yet accumulated any
error, and thus, $p(\xi_{0}^{+})=p(\xi_{0}^{-})=p_{0}$. When $\theta\in(0,1)$,
(15) is a linear interpolation between the two integrals, and it is easy to
verify that $p^{\prime}(\xi)=\bm{\tau}\cdot\nabla\tilde{p}=\tilde{g}_{\tau}$.
This modified regularization guarantees that the resulting $\tilde{h}_{d}$ is
continuous even if the contaminated pressure gradient on the boundary is not
curl-free, and thus this regularization can be considered as a curl-free
boundary correction. Note, (15) is one of the simplest solutions, but not
necessarily the only one or the best one; nevertheless, its advantage will be
demonstrated in Sect. 5.

##  5 Numerical Experiments and Validation

To validate the proposed boundary regulation strategies proposed in (5) and
(15), a synthetic flow based on the Taylor vortex contaminated by artificial
error is considered numerically. The pressure field of the Taylor vortex is

$p=-\frac{\rho H^{2}}{64\pi^{2}\nu t^{3}}\exp\left(-\frac{r^{2}}{2\nu
t}\right),$ (16)

where, $H$ represents the angular moment of the vortex, $\nu$ the kinematic
viscosity of the fluid, $\rho$ the density of the fluid, $t$ the time, and $r$
the distance from the origin in polar coordinates. We choose the parameters so
that the characteristic length scale of the vortex is $L_{0}=\sqrt{2\nu t}=1$
and the leading coefficient $\frac{\rho H^{2}}{64\pi^{2}\nu t^{3}}$ as well as
the amplitude of the characteristic pressure for this flow is unity
($P_{0}=1$). In this numerical experiment, the ground truth of the pressure
field, the pressure gradients (i.e., $g_{n}$ and $g_{\tau}$) and the pressure
Laplacian $f$ are generated from (16) on an $L\times L$ square domain. Bias
error (i.e., $\epsilon_{g_{n}}=1$, $\epsilon_{g_{\tau}}=1$ and
$\epsilon_{f}=-1$) were added to the true solution to generate the synthetic
data. This constant bias represents a global error with dominating low-
frequency components and is typically challenging (Faiella et al.,, 2021).
Using this synthetic data, the PPE was solved with the commonly used cNBC, the
proposed dDBC, and the improved curl-free dDBC while varying the size of the
domain. The reconstructed pressure was then compared to the ground truth to
evaluate the error associated with each numerical experiment and is
illustrated in Fig. 1.

Figure 1: Error in the reconstructed pressure field
($||\epsilon_{p}||_{L^{2}(\Omega)}/P_{0}$) scaling with the length scale of
the domain ($L/L_{0}$). Symbols represent the results from numerical
experiments; the red and blue solid lines illustrate the error estimates
indicated in (13) and (14), respectively.

As seen in Fig. 1, the results from the numerical experiments match the error
estimates derived in Sect. 3 well. For the cNBC setup, the error estimate in
(13) is relatively sharp, while the upper bound in (14) is slightly more
conservative for the dDBC setup, meaning that the actual performance of the
pressure solver with dDBC regularization is even better than expected. Both
error estimates reveal the scaling behavior of the error in the pressure with
the length scale of the domain $L$. For the two-dimensional cases, when the
domain is large, the error in the domain ($\epsilon_{f}$) dominates and
$\|\epsilon_{p}\|_{L^{2}(\Omega)}\sim L^{2}$; while for a small domain, error
on the boundary ($\epsilon_{g_{n}}$ or $\epsilon_{g_{\tau}}$) dominates and
$\|\epsilon_{p}\|_{L^{2}(\Omega)}\sim L^{1}$. The curl-free dDBC proposed in
Sect. 4 also performs as expected, resulting in even less error than the two
aforementioned methods, especially when the domain is small. This test
showcases the power of these simple regularization strategies and validates
the error estimates outlined in Sect. 3.

The proposed boundary regularization strategies are further tested using
synthetic data of a wake behind a cylinder with $Re=100$, based on free stream
velocity $U_{\infty}=1$. A cropped and down-sampled high-fidelity simulation
of the wake flow was used as the ground truth of the velocity and pressure
fields. The cylinder of diameter $D$ was centered at $(x/D,y/D)=(0,0)$, and
the cropped domain is $(x/D,y/D)\in[2,6]\times[-4,4]$. The non-dimensional
spatial and temporal resolution of the down-sampled data is $dx/D=dy/D=0.0625$
and $dtD/U_{\infty}=0.1$ respectively. Point-wise Gaussian noise with a
variance equal to $1\%$ of the magnitude of the free stream velocity was added
to the ground truth velocity to generate the synthetic data. We reconstruct
the pressure field using a second-order finite difference Poisson solver using
the artificially contaminated velocity data and one point reference pressure
at the bottom left corner of the domain $(x/D,y/D)=(2,-4)$ with cNBC, dDBC,
and curl-free dDBC. The reconstructed pressure and the error are normalized by
$P_{0}=\frac{1}{2}\rho U_{\infty}^{2}$. 500 independent tests were carried out
and depicted in Fig. 2.

Figure 2: Typical results and statistical tests for pressure reconstruction
with and without boundary regularization. (a) quiver plot of the velocity
field overlaid on the pressure field ground truth; (b–d) reconstructed
pressure field with cNBC, dDBC, and curl-free dDBC, respectively; (f–h) error
in the reconstructed pressure field by comparing (b–d) with (a). (e) box plot
of error in the pressure reconstruction from 500 independent tests with the
red, blue, and green boxes corresponding to the statistics of the error for
cNBC, dDBC, and curl-free dDBC cases, respectively. Horizontal bars in the
middle of the boxes show the median while the upper and lower edges of the box
indicate the 25 and 75 percentiles. The upper and lower whiskers bound the 95%
confidence intervals of the error while the symbols within the boxes mark
where the corresponding error shown in (f–h) lie within the data.

The statistical tests in Fig. 2 show the performance of the three different
boundary conditions with both the dDBCs performing much better than the cNBC
as expected. In particular, the median error in the dDBC is reduced by roughly
a factor of 10 when compared to the cNBC, and the curl-free dDBC reduces this
error by an additional factor of $\sim$2\. This is consistent with i) the
theory that an elongated rectangular domain can reduce error propagation even
further when compared to a square domain as indicated by (13) and (14); and
ii) that the error estimate in (12) is rather conservative. In addition to
greatly reducing the expected error in pressure (i.e., the mean of
$\|\epsilon_{p}\|_{L^{2}(\Omega)}/P_{0}$ in Fig. 2(e) are 1.83, 0.17, and 0.11
for cNBC, dDBC, and curl-free dDBC case, respectively), the variance of
$\|\epsilon_{p}\|_{L^{2}(\Omega)}/P_{0}$ is also reduced by the boundary
regularization (e.g., the standard deviation of
$\|\epsilon_{p}\|_{L^{2}(\Omega)}/P_{0}$ are 1.27, 0.10, and 0.05, cNBC, dDBC,
and curl-free dDBC cases, respectively). This means that the proposed boundary
regularization can improve both the accuracy and precision of pressure
reconstruction. The observation in this test can also be explained from the
numerical perspective: The resulting discretized systems from the PPE (i.e.,
$\mathcal{L}p=f$, where $\mathcal{L}$ is the discretized Laplacian with
associated boundary conditions) have 2-norm condition numbers of $1.1\times
10^{4}$ and $4.6\times 10^{5}$ with and without regularization, respectively,
meaning that the regularization improves the numerical conditioning of the
problem as well. The effectiveness of the boundary regularization can be
directly visualized by Fig. 2(b–d) or Fig. 2(f–h), which is a typical case out
of the 500 tests. Without boundary regularization, high bias in the
reconstructed pressure ($\tilde{p}$) occurs due to the application of long
Neumann boundaries, except for the bottom left corner where a transducer
provides an accurate reference pressure (see Fig. 2(b&f)). From Fig. 2(c&g),
it is obvious that replacing long cNBC with dDBC can effectively tame the
error propagation; however, error accumulating along the boundary due to the
practice of the regularization of (5) in the counter-clockwise direction is
propagated to the interior of the domain (i.e., the $\epsilon_{p}$ is higher
along the left edge of the domain than that of the bottom edge). This artifact
is resolved by the curl-free dDBC regularization (see Fig. 2(d&h)). This
suggests that even the curl-free correction on the boundary alone typically
results in further reduction in $\epsilon_{p}$ throughout the entire domain.

##  6 Concluding Remarks and Perspectives

This work portrays a simple yet overlooked idea to regularize the velocity-
based pressure field reconstruction: replacing the canonical Neumann boundary
conditions with the Dirichlet conditions derived from integrating the
tangential part of the pressure gradient on the boundary. This simple practice
can effectively improve the quality, in terms of accuracy and robustness, of
the pressure field reconstruction from corrupted image velocimetry data, which
is analytically proved and validated by numerical experiments. The
corresponding error estimates provided in this work are conservative i.e.
pessimistic, but still dictate the generic dynamics of error propagation.

The proposed regularization strategies are very simple, and presumably can
make Dirichlet conditions available everywhere with low computational cost.
This is attractive as the derived Dirichlet conditions provide a ‘stronger’
type (Faiella et al.,, 2021) of information on the boundary than cNBC alone.
This ‘additional’ information on the boundary provides great flexibility and
potential for more sophisticated algorithms to further improve the
reconstruction. Our analysis is independent of the dimension of the domain (2D
or 3D), and the continuous setting of the analysis indicates that the strategy
works for structured or unstructured data based on PIV or LPT. Beyond
establishing this family of regularization strategies, fulfilling the full
potential of this seemingly odd, yet effective idea, will be left for future
investigations.

##  7 Acknowledgement

We thank the discussion with Drs. Fernando Zigunov and Grady Wright. This work
is partially supported by the Natural Sciences and Engineering Research
Council of Canada (NSERC) Discovery Grant (RGPIN-2020-04486), and the
Undergraduate Research Assistantship (URA) program, University of Waterloo.
J.P.W. was partially supported by NSF grant DMS-2206762.

## References

  * Charonko et al., (2010) Charonko, J. J., King, C. V., Smith, B. L., and Vlachos, P. P. (2010).  Assessment of pressure field calculations from particle image velocimetry measurements.  Measurement Science and Technology, 21(10):105401.
  * De Kat and Van Oudheusden, (2012) De Kat, R. and Van Oudheusden, B. (2012).  Instantaneous planar pressure determination from PIV in turbulent flow.  Experiments in fluids, 52(5):1089–1106.
  * Faiella et al., (2021) Faiella, M., Macmillan, C. G. J., Whitehead, J. P., and Pan, Z. (2021).  Error propagation dynamics of velocimetry-based pressure field calculations (2): on the error profile.  Measurement Science and Technology, 32(8):084005.
  * Pan et al., (2016) Pan, Z., Whitehead, J., Thomson, S., and Truscott, T. (2016).  Error propagation dynamics of piv-based pressure field calculations: How well does the pressure poisson solver perform inherently?  Measurement Science and Technology, 27(8):084012.
  * Sperotto et al., (2022) Sperotto, P., Pieraccini, S., and Mendez, M. A. (2022).  A meshless method to compute pressure fields from image velocimetry.  Measurement Science and Technology, 33(9):094005.
  * Van Oudheusden, (2013) Van Oudheusden, B. (2013).  PIV-based pressure measurement.  Measurement Science and Technology, 24(3):032001.
  * Zhang et al., (2022) Zhang, J., Bhattacharya, S., and Vlachos, P. P. (2022).  Uncertainty of piv/ptv based eulerian pressure estimation using velocity uncertainty.  Measurement Science and Technology, 33(6):065303.


# Schwarzschild-Finsler-Randers spacetime:
Dynamical analysis, Geodesics and Deflection Angle

E. Kapsabelis<EMAIL_ADDRESS>Section of Astrophysics, Astronomy
and Mechanics, Department of Physics, National and Kapodistrian University of
Athens, Panepistimiopolis 15784, Athens, Greece  P. G. Kevrekidis
<EMAIL_ADDRESS>Department of Mathematics and Statistics, University of
Massachusetts, Amherst, MA 01003-9305, USA  P. C. Stavrinos
<EMAIL_ADDRESS>Department of Mathematics, National and Kapodistrian
University of Athens, Panepistimiopolis 15784, Athens, Greece  A.
Triantafyllopoulos<EMAIL_ADDRESS>Section of Astrophysics, Astronomy
and Mechanics, Department of Physics, National and Kapodistrian University of
Athens, Panepistimiopolis 15784, Athens, Greece

# Schwarzschild-Finsler-Randers spacetime:
Geodesics, Dynamical Analysis and Deflection Angle

E. Kapsabelis<EMAIL_ADDRESS>Section of Astrophysics, Astronomy
and Mechanics, Department of Physics, National and Kapodistrian University of
Athens, Panepistimiopolis 15784, Athens, Greece  P. G. Kevrekidis
<EMAIL_ADDRESS>Department of Mathematics and Statistics, University of
Massachusetts, Amherst, MA 01003-9305, USA  P. C. Stavrinos
<EMAIL_ADDRESS>Department of Mathematics, National and Kapodistrian
University of Athens, Panepistimiopolis 15784, Athens, Greece  A.
Triantafyllopoulos<EMAIL_ADDRESS>Section of Astrophysics, Astronomy
and Mechanics, Department of Physics, National and Kapodistrian University of
Athens, Panepistimiopolis 15784, Athens, Greece

###### Abstract

In this work, we extend the study of Schwarzschild-Finsler-Randers (SFR)
spacetime previously investigated by a subset of the present authors [1, 2].
We will examine the dynamical analysis of geodesics which provides the
derivation of the energy and the angular momentum of a particle moving along a
geodesic of SFR spacetime. This study allows us to compare our model with the
corresponding of general relativity (GR). In addition, the effective potential
of SFR model is examined and it is compared with the effective potential of
GR. The phase portraits generated by these effective potentials are also
compared. Finally we deal with the derivation of the deflection angle of the
SFR spacetime and we find that there is a small perturbation from the
deflection angle of GR. We also derive an interesting relation between the
deflection angles of SFR model and the corresponding result in the work of
[3]. These small differences are attributed to the anisotropic metric
structure of the model and especially to a Randers term which provides a small
deviation from the GR.

##  I Introduction

Einstein’s field equations in general relativity predict that the curvature is
produced not only by the distribution of mass-energy but also by its motion
[4]. Candidate metric geometries that can intrinsically describe the motion
are the Finsler and Finsler-like geometries which constitute metrical
generalizations of Riemannian geometry and depend on position and
velocity/momentum/scalar coordinates. These are dynamic geometries that can
describe locally anisotropic phenomena and Lorentz violations [5, 6, 7, 8, 9,
10, 11, 12, 13, 14, 15] as well as with field equations, FRW and Raychaudhuri
equations, geodesics, dark matter and dark energy effects [20, 19, 17, 18, 16,
21]. By considering this approach, the gravitational field is interpreted as
the metric of a generalized spacetime and constitutes a force-field which
contains the motion. This possibility reveals the Finslerian geometrical
character of spacetime.

In the framework of applications of Finsler geometry, many works in different
directions of geometrical and physical structures have contributed to the
extension of research for theoretical and observational approaches during the
last years. We cite some works from the literature of the applications of
Finsler geometry [22, 23, 7, 8, 24, 25, 26, 27, 28, 13, 29, 30, 31, 21, 32].

In the first period of development of applications of Finsler geometry to
Physics, especially to General Relativity, remarkable works were published by
G. Randers [33], J. I. Horváth [34] and A. Moór [35]. Later, Einstein’s field
equations were formulated in the Finslerian framework by the works of J. I.
Horváth [34, 35], Y. Takano [36] and S. Ikeda [37]. In these studies, the
field equations had been considered without calculus of variations. G. S.
Asanov [38] explored the Finslerian gravitational field by using Riemannian
osculating methods and derived Einstein field equations using the variational
principle. A class of Finsler spaces (FR standing for Finsler-Randers)
originated by G. Randers [33] who studied the physical properties of spacetime
with an asymmetrical metric which provides the uni-direction of time-like
intervals. This consideration gives a particular interest in a generalized
metric structure of the Riemannian spacetime. Based on this form of spacetime,
it is possible to investigate the gravitational field with more degrees of
freedom in the framework of a tangent/vector/scalar bundle [39, 40, 31]. The
FR cosmological model was first introduced in [41, 42]. It is of special
interest since the Friedmann equations include an extra geometrical term that
acts as a dark energy-fluid. The Finsler-Randers-type spacetime can be
considered as a direction-dependent motion of the Riemannian/FRW model.

The local anisotropic structure of spacetime affects the gravitational field
and leads to modified cosmological considerations. Based on Finsler or
Finsler-like cosmologies, the Friedmann equations include extra terms which
influence the cosmological evolution [17, 39, 31, 20]. When Lorentz symmetry
holds, the spacetime is isotropic in the sense that all directions and uniform
motions are equivalent. The introduction of a vector field in the structure of
spacetime causes relativity violations and local anisotropy which arise from
breaking the Lorentz symmetry and which affect the metric, curvature,
geodesics and null cone [43, 44, 45, 46, 47, 48, 49, 50].

In the framework of modified gravitational theories with Finsler-Randers type
structure, two fundamental theories of investigation for the gravitation and
cosmology can be developed. The first one is connected with the Friedmann-
Finsler-Randers cosmological model and the second one is related to the study
of SFR spacetime.

The introduction of a force field causes an asymmetry to a pseudo-Finslerian
metric. Asymmetrical and locally anisotropic models such as Finsler-Randers
spacetime can be connected with the chiral fields in Cosmology for
descriptions of the inflationary epoch and the present accelerated expansion
of the Universe [51].

An FR space has a metric function of the form

$F(x,y)=(-a_{\mu\nu}(x)y^{\mu}y^{\nu})^{1/2}+u_{\alpha}y^{\alpha}$ (1)

where $u_{\alpha}$ is a covector with $||u_{\alpha}||\ll 1$,
$y^{\alpha}=\frac{dx^{\alpha}}{d\tau}$ and $a_{\mu\nu}(x)$ is a pseudo-
Riemannian metric for which the Lorentzian signature $(-,+,+,+)$ has been
assumed and the indices $\mu,\nu,\alpha$ take the values $0,1,2,3$. The
geodesics of this space can be produced by (1) and the Euler-Lagrange
equations. If $u_{\alpha}$ denotes a force field $f_{\alpha}$ and $y^{\alpha}$
is substituted with $dx^{\alpha}$ then $f_{\alpha}dx^{\alpha}$ represents the
spacetime effective energy produced by the anisotropic force field
$f_{\alpha}$, therefore equation (1) is written as

$F(x,dx)=\left(-a_{\mu\nu}(x)dx^{\mu}dx^{\nu}\right)^{1/2}+f_{\alpha}dx^{\alpha}$
(2)

This form of metric provides a dynamical effective structure of spacetime. A
small differentiation is presented between GR and the FR gravitation model.
This is because of the work provided by the one-form $A_{\gamma}$ which gives
an external motion to the Riemannian spacetime. This motion is an internal
concept for the FR spacetime.

A cosmological model can be introduced by Eq. (2) if we assume the FRW
cosmological metric instead of the general type of the Riemannian one [41,
42]. In this case, we get a Friedmann-Finsler-Randers cosmological model in
the following form

$a_{\mu\nu}(x)=\mathrm{diag}\left[-1,\frac{a^{2}}{1-\kappa
r^{2}},a^{2}r^{2},a^{2}r^{2}\sin^{2}\theta\right]$ (3)

This model was also further studied later in [58, 59, 9, 60, 61, 52, 64, 62,
63, 65, 53, 54, 55, 66, 67, 68, 56, 21, 57].

In this work, we will follow the SFR model presented in [1]. The metric
$g_{\mu\nu}$ is the classic Schwarzschild one:

$\displaystyle
g_{\mu\nu}\mathrm{d}x^{\mu}\mathrm{d}x^{\nu}=-fdt^{2}+\frac{dr^{2}}{f}+r^{2}d\theta^{2}+r^{2}\sin^{2}\theta\,d\phi^{2}$
(4)

with $f=1-\frac{R_{s}}{r}$ and $R_{s}=2GM$ the Schwarzschild radius (we assume
units where the speed of light $c=1$).

Hereafter, we consider an $\alpha$-Randers type metric as the one in Eq. (1)
which is distinguished from the $\beta$-Randers type metric that is
investigated in the Standard Model Extension (SME) [7, 8, 46, 11].

The metric $v_{\alpha\beta}$ is derived from a metric function $F_{v}$ of the
$\alpha$-Randers type [1]:

$F_{v}=\sqrt{-g_{\alpha\beta}(x)y^{\alpha}y^{\beta}}+A_{\gamma}(x)y^{\gamma}$
(5)

where
$g_{\alpha\beta}=g_{\mu\nu}\tilde{\delta}^{\mu}_{\alpha}\tilde{\delta}^{\nu}_{\beta}$
is the Schwarzschild metric from Eq. (4) and $A_{\gamma}(x)$ is a covector
which expresses a deviation from general relativity, with $|A_{\gamma}(x)|\ll
1$, i.e., we assume that the deviation is small. In this work, we continue the
investigation of the Schwarzschild-Finsler-Randers spacetime (SFR) which has
been studied in previous works by a subset of the present authors [1, 2].

In this article we examine the influence of extra gravitational effect which
is imprinted in the geodesics of an SFR spacetime and we compare it with that
of GR case. We prove that the additional amount of energy included in the
equation of geodesics (see Eqs. (28)-(31)) stems from the geometry and the
form of SFR model.

This approach is also applied the Newtonian gravitational theory and we obtain
an interesting relation between the fundamental term $A_{0}(r)$ of SFR
gravitational theory and the Newtonian potential. This relation also provides
a physical interpretation of the extra term $A_{0}(r)$ in the Newtonian
framework.

In addition we calculate the deflection angle for the SFR model and we compare
its value with that of the GR case. We also show that there is a relation
between the deflection angle of SFR and of Shapiro et al. [3] observational
result of the deflection angle for very small angles ($\phi\approx 0$) near a
fiducial geodesic.

The present work is organized as follows. The structure of the model is given
in Section II. In this framework, the geodesics are studied and a dynamical
analysis is presented in Section III. We also compare our results with GR and
discuss the corresponding similarities and differences and we give an
application to the Newtonian framework. A dynamical analysis for the effective
potential of this spacetime is provided in the Section IV, where upon suitable
assumptions, the phase portraits of both models (SFR and GR) are presented. In
sec. V, we study the deflection angle of the SFR spacetime and we compare it
with very small values of the deflection angle of Shapiro et al. [3]. Finally,
the conclusions of our study and some possible directions for future
exploration are presented in Section VI.

##  II Basic structure of the model

In this section, we briefly present the underlying geometry of the SFR
gravitational model, as well as the field equations for the SFR metric. The
solution of these equations for this metric is presented at the end of the
section. An extended study of this model can be found in [1, 40]. The Lorentz
tangent bundle $TM$ is a fibered 8-dimensional manifold with local coordinates
$\\{x^{\mu},y^{\alpha}\\}$ where the indices of the $x$ variables are
$\kappa,\lambda,\mu,\nu,\ldots=0,\ldots,3$ and the indices of the $y$
variables are $\alpha,\beta,\ldots,\theta=4,\ldots,7$. The tangent space at a
point of $TM$ is spanned by the so-called adapted basis
$\\{E_{A}\\}=\,\\{\delta_{\mu},\dot{\partial}_{\alpha}\\}$ with

$\delta_{\mu}=\dfrac{\delta}{\delta x^{\mu}}=\dfrac{\partial}{\partial
x^{\mu}}-N^{\alpha}_{\mu}(x,y)\dfrac{\partial}{\partial y^{\alpha}}$ (6)

and

$\dot{\partial}_{\alpha}=\dfrac{\partial}{\partial y^{\alpha}}$ (7)

where $N^{\alpha}_{\mu}$ are the components of a nonlinear connection
$N=N^{\alpha}_{\mu}(x,y)\,\mathrm{d}x^{\mu}\otimes\dot{\partial}_{\alpha}$.

The nonlinear connection induces a split of the total space $TTM$ into a
horizontal distribution $T_{H}TM$ and a vertical distribution $T_{V}TM$. The
above-mentioned split is expressed with the Whitney sum:

$TTM=T_{H}TM\oplus T_{V}TM$ (8)

The anholonomy coefficients of the nonlinear connection are defined as

$\Omega^{\alpha}_{\nu\kappa}=\dfrac{\delta N^{\alpha}_{\nu}}{\delta
x^{\kappa}}-\dfrac{\delta N^{\alpha}_{\kappa}}{\delta x^{\nu}}$ (9)

A Sasaki-type metric [69, 70] $\mathcal{G}$ on $TM$ is:

$\mathcal{G}=g_{\mu\nu}(x,y)\,\mathrm{d}x^{\mu}\otimes\mathrm{d}x^{\nu}+v_{\alpha\beta}(x,y)\,\delta
y^{\alpha}\otimes\delta y^{\beta}$ (10)

where we have defined the metrics $g_{\mu\nu}$ and $v_{\alpha\beta}$ to be
pseudo-Finslerian.

A pseudo-Finslerian metric $f_{\alpha\beta}(x,y)$ is defined as one that has a
Lorentzian signature of $(-,+,+,+)$ and that also obeys the following form:

$\displaystyle
f_{\alpha\beta}(x,y)=\pm\frac{1}{2}\dfrac{\partial^{2}F^{2}}{\partial
y^{\alpha}\partial y^{\beta}}$ (11)

where the function $F$ satisfies the following conditions [69]:

  1. 1.

$F$ is continuous on $TM$ and smooth on $\widetilde{TM}\equiv
TM\setminus\\{0\\}$, i.e., the tangent bundle minus the null set $\\{(x,y)\in
TM|F(x,y)=0\\}$ .

  2. 2.

$F$ is positively homogeneous of first degree on its second argument:

$F(x^{\mu},ky^{\alpha})=kF(x^{\mu},y^{\alpha}),\qquad k>0$ (12)

  3. 3.

The form

$f_{\alpha\beta}(x,y)=\dfrac{1}{2}\dfrac{\partial^{2}F^{2}}{\partial
y^{\alpha}\partial y^{\beta}}$ (13)

defines a non-degenerate matrix:

$\det\left[f_{\alpha\beta}\right]\neq 0$ (14)

where the plus-minus sign in (11) is chosen so that the metric has the correct
signature.

In the following, we choose a non-linear connection with the following form:

$N^{\alpha}_{\mu}=\frac{1}{2}y^{\beta}g^{\alpha\gamma}\partial_{\mu}g_{\beta\gamma}$
(15)

The metric tensor $v_{\alpha\beta}$ is derived from (5) by using (11), after
omitting higher order terms $O(A^{2})$:

$v_{\alpha\beta}(x,y)=g_{\alpha\beta}(x)+h_{\alpha\beta}(x,y)$ (16)

where

$\displaystyle
h_{\alpha\beta}=\frac{1}{\tilde{a}}(A_{\beta}g_{\alpha\gamma}y^{\gamma}+A_{\gamma}g_{\alpha\beta}y^{\gamma}+A_{\alpha}g_{\beta\gamma}y^{\gamma})+\frac{1}{\tilde{a}^{3}}A_{\gamma}g_{\alpha\epsilon}g_{\beta\delta}y^{\gamma}y^{\delta}y^{\epsilon}$
(17)

with $\tilde{a}=\sqrt{-g_{\alpha\beta}y^{\alpha}y^{\beta}}$. The total metric
defined in the previous steps is called the Schwarzschild-Finsler-Randers
(SFR) metric. As we can see, the term $h_{\alpha\beta}(x,y)$ can be considered
as a perturbation of the Schwarzschild metric since $|A_{\gamma}(x)|\ll 1$.

The nonzero coefficients of a canonical and distinguished $d-$connection
$\mathcal{D}$ on $TM$ are:

$\displaystyle L^{\mu}_{\nu\kappa}$
$\displaystyle=\frac{1}{2}g^{\mu\rho}\left(\delta_{k}g_{\rho\nu}+\delta_{\nu}g_{\rho\kappa}-\delta_{\rho}g_{\nu\kappa}\right)$
(18) $\displaystyle L^{\alpha}_{\beta\kappa}$
$\displaystyle=\dot{\partial}_{\beta}N^{\alpha}_{\kappa}+\frac{1}{2}v^{\alpha\gamma}\left(\delta_{\kappa}v_{\beta\gamma}-v_{\delta\gamma}\,\dot{\partial}_{\beta}N^{\delta}_{\kappa}-v_{\beta\delta}\,\dot{\partial}_{\gamma}N^{\delta}_{\kappa}\right)$
(19) $\displaystyle C^{\mu}_{\nu\gamma}$
$\displaystyle=\frac{1}{2}g^{\mu\rho}\dot{\partial}_{\gamma}g_{\rho\nu}$ (20)
$\displaystyle C^{\alpha}_{\beta\gamma}$
$\displaystyle=\frac{1}{2}v^{\alpha\delta}\left(\dot{\partial}_{\gamma}v_{\delta\beta}+\dot{\partial}_{\beta}v_{\delta\gamma}-\dot{\partial}_{\delta}v_{\beta\gamma}\right)$
(21)

See Appendix A for more details.

The field equations for our model have been derived in previous works and can
be found in Appendix B. The solution of the field equations (98), (99) and
(100) to first order in $A_{\gamma}(x)$ in vacuum
($T_{\mu\nu}=Y_{\alpha\beta}=\mathcal{Z}^{\kappa}_{\alpha}=0$) is [1]:

$A_{\gamma}(x)=\left[\tilde{A}_{0}\left(1-\frac{R_{S}}{r}\right)^{1/2},0,0,0\right]=\left[\tilde{A}_{0}f^{1/2},0,0,0\right]$
(22)

with $\tilde{A}_{0}$ a constant. While this is an approximate solution, it
will be sufficient for our purposes given the assumption $|A_{\gamma}(x)|\ll
1$.

##  III Geodesics

In this section, we will study the geodesics of the SFR and perform a
dynamical analysis. We compare our results with the corresponding ones of GR.
From the definition of the metric function (5) we have:

$F(x,dx)=\left(-g_{\mu\nu}(x)dx^{\mu}dx^{\nu}\right)^{1/2}+A_{\gamma}(x)dx^{\gamma}$
(23)

where $g_{\mu\nu}(x)$ is the Schwarzschild metric and $A_{\gamma}(x)$ is a
one-form vector field with $|A_{\gamma}(x)|\ll 1$.
By using the rel. (4) the rel. (23) is written as:

$F(x,dx)=\Big{[}fdt^{2}-\frac{dr^{2}}{f}-r^{2}d\theta^{2}-r^{2}\sin^{2}\theta\,d\phi^{2}\Big{]}^{1/2}+A_{\gamma}(x)dx^{\gamma}$
(24)

We define the Lagrangian

$L(x,\dot{x})=F(x,\dot{x})=\Big{[}f\dot{t^{2}}-\frac{\dot{r^{2}}}{f}-r^{2}\dot{\theta}^{2}-r^{2}sin^{2}\theta\dot{\phi^{2}}\Big{]}^{1/2}+\tilde{A_{0}}f^{1/2}\dot{t}$
(25)

where we denote $\dot{x}=\frac{dx}{d\tau}$ and we have used Eqs. (24) and
(22). From the Euler-Lagrange equations

$\frac{d}{d\tau}\frac{\partial L}{\partial\dot{x}^{\mu}}=\frac{\partial
L}{\partial x^{\mu}}$ (26)

we find the equations for the geodesics:

$\ddot{x}^{\lambda}+\Gamma^{\lambda}_{\mu\nu}\dot{x}^{\mu}\dot{x}^{\nu}+g^{\kappa\lambda}\Phi_{\kappa\mu}\dot{x}^{\mu}=0$
(27)

where $\Gamma^{\lambda}_{\mu\nu}$ are the Christoffel symbols of Riemann
geometry, $\dot{x}^{\mu}=\frac{dx^{\mu}}{d\tau}$ and
$\Phi_{\kappa\mu}=\partial_{\kappa}A_{\mu}-\partial_{\mu}A_{\kappa}$ and
$A_{\mu}$ is the solution rel.(22). We notice that from the definition of
$\Phi_{\kappa\mu}$ we get a rotation form of geodesics. If $A_{\mu}$ is a
gradient of a scalar field, $A_{\mu}=\dfrac{\partial\Phi}{\partial x^{\mu}}$
then $\Phi_{\kappa\mu}=0$ and the geodesics of our model are identified with
the Riemannian ones.
The geodesics of our model can then be explicitly written in the form:

$\displaystyle\ddot{t}+\frac{1-f}{rf}\dot{r}\dot{t}=-\tilde{A}_{0}\dot{r}\frac{f^{-3/2}(1-f)}{2r}$
(28)
$\displaystyle\ddot{r}+\frac{f(1-f)}{2r}\dot{t}^{2}-\frac{1-f}{2rf}\dot{r}^{2}-rf\big{(}\dot{\theta}^{2}+\sin^{2}\theta\dot{\phi}^{2}\big{)}=-\tilde{A}_{0}\dot{t}\frac{f^{1/2}(1-f)}{2r}$
(29)
$\displaystyle\ddot{\theta}+\frac{2}{r}\dot{\theta}\dot{r}-\frac{1}{2}\sin
2\theta\,\dot{\phi}^{2}=0$ (30)
$\displaystyle\ddot{\phi}+\frac{2}{r}\dot{\phi}\dot{r}+2\cot\theta\,\dot{\theta}\dot{\phi}=0$
(31)

From the relations (28)-(31), we notice that the first two dynamical equations
involve a contribution of extra terms particular to the SFR spacetime while
the last two relations are the same as in GR. From the above mentioned
relations, we notice that the Riemannian geodesics are affected by a
contribution that involves the curl of the force field in the SFR spacetime
which can be interpreted as extra energy for the content of GR. We can notice
that the geodesics of SFR are destroyed at the singular value $r=0$ as in the
GR case, this can be seen from the relations (28)-(31). In our case, the
singularity is inherited from the Schwarzschild spacetime of GR. Moreover, it
is evident from their dynamical form that the equations are not meaningful for
$r\leq R_{S}$ in the context of the presently considered Schwarzschild metric,
hence we only consider radial displacements past this singular point.

We now make a key assumption regarding the angular dependence of the model.
Namely, by using $\theta=\frac{\pi}{2}$ we notice that Eq. (30) is satisfied
and equations (28), (29) and (31) can be written as:

$\displaystyle\ddot{t}+\frac{1-f}{rf}\dot{r}\dot{t}=-\tilde{A}_{0}\dot{r}\frac{f^{-3/2}(1-f)}{2r}$
(32)
$\displaystyle\ddot{r}+\frac{f(1-f)}{2r}\dot{t}^{2}-\frac{1-f}{2rf}\dot{r}^{2}-rf\dot{\phi}^{2}=-\tilde{A}_{0}\dot{t}\frac{f^{1/2}(1-f)}{2r}$
(33) $\displaystyle\ddot{\phi}+\frac{2}{r}\dot{\phi}\dot{r}=0$ (34)

From Eq.(34) we find:

$r^{2}\dot{\phi}=J=const.$ (35)

where $J$ is the angular momentum and the relevant equation represents its
conservation law. If we use the relation $f^{\prime}=\frac{1-f}{r}$ where
$f=1-\frac{2GM}{r}$ and the Leibniz chain-rule
$\frac{d}{d\tau}=\frac{dr}{d\tau}\frac{d}{dr}=\dot{r}\frac{d}{dr}$, then Eq.
(32) can be written as:

$f\ddot{t}+\frac{df}{d\tau}\dot{t}=-\tilde{A}_{0}\frac{df^{1/2}}{d\tau}$ (36)

which, in turn, gives us

$f\dot{t}+\tilde{A}_{0}f^{1/2}=\mathcal{E}_{R}=const.$ (37)

where $\mathcal{E}_{R}$ is the energy of the particle moving along the
geodesic. We notice that the first term constitutes the energy for a particle
moving along the geodesics in general relativity, $\mathcal{E}_{GR}=f\dot{t}$
and we can rewrite the relevant expression as:

$\mathcal{E}_{GR}+\tilde{A}_{0}f^{1/2}=\mathcal{E}_{R}.$ (38)

By using Eq. (33) with (35), and (37), we arrive at the (effectively one-
degree-of-freedom) radial equation:

$\ddot{r}+\frac{1-f}{2rf}(\mathcal{E}_{R}^{2}-\dot{r}^{2})-\frac{fJ^{2}}{r^{3}}=\tilde{A}_{0}\mathcal{E}_{R}\frac{f^{-1/2}(1-f)}{2r}$
(39)

where we omitted $O(\tilde{A}_{0}^{2})$ terms. As before, we use the relation
$f^{\prime}=\frac{1-f}{r}$ in (39) to bring it to the equivalent form:

$\ddot{r}+\frac{f^{\prime}}{2f}(\mathcal{E}_{R}^{2}-\dot{r}^{2})-\frac{fJ^{2}}{r^{3}}=\frac{\tilde{A}_{0}\mathcal{E}_{R}}{2}f^{-1/2}f^{\prime}$
(40)

We can further simplify the Eq. (40) by using the Leibniz chain-rule
$\frac{d}{d\tau}=\frac{dr}{d\tau}\frac{d}{dr}=\dot{r}\frac{d}{dr}$ and upon
deriving the first integral of the motion, we obtain:

$\dot{r}^{2}+f\left(\frac{J^{2}}{r^{2}}+\epsilon\right)+2\tilde{A}_{0}\mathcal{E}_{R}f^{1/2}=\mathcal{E}_{R}^{2}$
(41)

where $\epsilon$ is a constant and for $\epsilon=0$ we have null geodesics. It
is important to indicate here that for each value of $\epsilon$, we obtain a
different curve of this “first integral” of Eq. (41) and the inclusion of all
of the admissible ($r$,$\dot{r}$) curves will provide us with the phase
portraits presented below. The first two terms from (41) constitute the total
energy in general relativity (GR),
$\mathcal{E}^{2}_{GR}=\dot{r}^{2}+f(\frac{J^{2}}{r^{2}}+\epsilon)$ and the
third term emerges from the structure of SFR spacetime and its energetic
contribution. Therefore Eq. (41) can be written as:

$\mathcal{E}^{2}_{GR}+2\tilde{A}_{0}\mathcal{E}_{R}f^{1/2}=\mathcal{E}_{R}^{2}$
(42)

Eq. (42) shows that the term $A_{\gamma}(x)$ from (22) provides an additional
energy contribution to the system of GR.

For a particle moving along the geodesics in the SFR spacetime with metric
function $F_{v}(x,dx)$ (rel. 24), we have extra gravitational effect (energy)
compared to GR. An amount of energy comes from the gravitational field of
total space in the SFR model. Therefore, from the rel. (37) and (42), the
second term $\tilde{A}_{0}f^{1/2}$ can be interpreted as the difference of
energy between $\mathcal{E}_{R}$ and $\mathcal{E}_{GR}=f\dot{t}$.

Remark 1

The action of a vector field in the pseudo-Riemannian structure of space-time
can give a vector-dependent gravitational field of Finsler-Randers (FR) type
which can describe the asymmetry and the locally anisotropic form of spacetime
that the Riemannian geometry is unable to provide. In this approach, because
of broken Lorentz symmetry, particles and forces interact with this vector
field.

In the following, we give an application for the weak gravitational field
where we show that the term $A_{0}(r)$ in the Eq. (22) is connected to the
gravitational potential $\Phi$.

Application to the Newtonian limit.
In the Newtonian gravitational theory, we can consider the second term of
rel.(2) of the manuscript in the following form:

$f_{a}dx^{a}=\frac{GmM}{r^{2}}dr$ (43)

The gravitational potential is given by

$U=\frac{W}{m}=\frac{1}{m}\int
f_{a}dx^{a}=\frac{1}{m}\int^{r}_{\infty}\frac{GmM}{r^{2}}dr=-\frac{GM}{r}=\Phi$
(44)

where $W$ represents the gravitational potential energy to be done to bring a
unit mass m from infinity to a point and $r$ is the distance from a mass $M$.
This means that the work is converted to gravitational potential energy. In
addition, from the relations (4) and (22) we have

$g_{00}=-\left(1-\frac{R_{S}}{r}\right)=-f$ (45)

and

$A_{0}(r)=\tilde{A}_{0}\left(1-\frac{R_{S}}{r}\right)^{1/2}=\tilde{A}_{0}f^{1/2}$
(46)

From the relations (44), (45) and (46), we obtain

$A_{0}(r)=\tilde{A}_{0}\sqrt{-g_{00}}=\tilde{A}_{0}\left(1-\frac{\Phi}{2}\right)^{1/2}$
(47)

where we have used $R_{S}=2GM$. The rel. (47) gives a physical interpretation
in the geometrical term $A_{0}(r)$ because of its dependence on the
gravitational potential $\Phi$. The term $A_{0}(r)$ is included in the
relations (37), (38), (42) giving a physical meaning to these equations.

Remark 2

The relation (47) is also related to gravitational redshift of photons in the
SFR spacetime (ref. [2], paragraph 6). In that case we have proved that :

$\frac{E_{rec}}{E_{emit}}=\frac{\nu_{rec}}{\nu_{emit}}\approx\tilde{A_{0}}f^{1/2}$
(48)

where $\nu_{rec},\nu_{emit}$ denote the frequencies of receiver and emitter.

Below, we give the Figures 1(a), 1(b) and Fig. 2 for the geodesics of GR and
SFR we have obtained by solving the equations (28)-(31). The relevant ordinary
differential equations are solved via a standard solver within Mathematica and
$(r,\phi)$ are presented as a function of $\tau$, while Fig. 2 presents the
evolution in the original $(x,y)$ plane. In our case, we assume $R_{s}=2$ and
initial radial distance $r_{0}=3$, so the photons are found on the
photonsphere with $r_{ph}=\frac{3}{2}R_{s}=3$ in the GR case. The deviation
between the trajectories of the SFR and those of the GR is clearly discernible
in both figures.

(a) This is an $(\tau,r)$ graph for the geodesics of photons for angular
momentum $J=4$ and initial radial distance $r_{0}=3$. The red line shows the
SFR geodesics and the blue line the GR geodesics

(b) This is an $(\tau,\phi)$ graph for the geodesics of photons for angular
momentum $J=4$ and initial radial distance $r_{0}=3$. The red line shows the
SFR geodesics and the blue line the GR geodesics

Fig. 1: These are the graphs for the GR and SFR geodesics of photons for
angular momentum $J=4$ and initial radial distance $r_{0}=3$

From Fig. 1(a), we can see that the radial component in the SFR model takes
lower values compared to the GR one which remains constant. This difference
between the r-components of SFR and GR can be interpreted as the increase of
the radius of the photonsphere due to the one-form $A_{\gamma}$ as we have
shown in [2]. This leads the orbit of the photon to fall inside the event
horizon because the initial distance $r_{0}=3$ and energy are not sufficient
to allow circular orbits of the photonsphere. That means for an orbit with $r$
constant in the SFR model, the particle needs more energy compared to the GR
case. In Fig. 2, the geodesics of GR and SFR are depicted. In the case of GR,
the photons move in circular orbits around the black hole. In the SFR model,
the photons follow a spiral orbit and fall inside the event horizon.

It is important to remind the reader here that underlying these results is the
key assumption of $\theta=\frac{\pi}{2}$ which allows the reduction of the
model to an effective single degree-of-freedom system. It is important in
future work to consider how deviations from this equilibrium value (and the
corresponding incorporation of the full dynamical system) may affect the
conclusions presented above. However, as the latter is outside the scope of
the present study, we now focus on the further analysis of the effective
potential of the SFR model and its implications for the phase portrait of the
relevant system.

Fig. 2: This is an $x-y$ graph for the geodesics of photons for angular
momentum $J=4$ and initial radial distance $r_{0}=3$. The red line shows the
SFR geodesics and the blue line the GR geodesics

##  IV Effective Potential of SFR model

In this section, we will study the effective potential of the SFR model and
compare it with the effective potential of GR. The equation of the energy in
GR reads:

$\dot{r}^{2}+f\left(\frac{J^{2}}{r^{2}}+\epsilon\right)=\mathcal{E}_{GR}^{2}$
(49)

We see from (49) that the effective potential energy landscape is given by:

$V_{eff,GR}=\frac{1}{2}f\left(\frac{J^{2}}{r^{2}}+\epsilon\right)$ (50)

In Fig.3(b), we show the graph for the effective potential in GR for angular
momentum $J=3,J=4$ and $J=5$ to examine its variation for different values of
the angular momentum. Our effective potential (here and in what follows in
Fig. 3(a)) solely bears a maximum at a finite $r>0$. Based on the general
theory of dynamical systems for the conservative type of problems considered
herein, the relevant fixed point will be a saddle, as will be confirmed in the
corresponding phase portraits below (where the attracting direction of the
stable manifold and the repelling direction of the unstable manifold will be
evident).

We now recall the key difference (and associated additional contribution) to
the energetics of the SFR model. In particular, the energy equation for the
latter, derived from Eq. (40), is given as:

$\dot{r}^{2}+f\left(\frac{J^{2}}{r^{2}}+\epsilon\right)+2\tilde{A}_{0}\mathcal{E}_{R}f^{1/2}=\mathcal{E}_{R}^{2}$
(51)

In (51) the effective potential can be written in the form

$V_{eff,SFR}=\frac{1}{2}f\left(\frac{J^{2}}{r^{2}}+\epsilon\right)+\tilde{A}_{0}\mathcal{E}_{R}f^{1/2}$
(52)

The graph for the effective potential in the SFR model $(V_{eff},r)$ is
depicted in Fig. 3(a), in this case for different values of angular momentum.

In Fig. 3(c) and 3(d) we show the effective potentials of the SFR and GR
models comparing the two for $J=1$ and $J=5$. As we can see in Fig. 3(c), the
difference between GR and SFR is bigger than that of Fig. 3(d). Notably, when
the contribution of the angular momentum is weaker, the difference between the
two models is more substantial/clearly discernible. When the angular momentum
becomes large, the relevant difference is rather weak and the $V_{eff}$ of the
two models become proximal.

(a) Graph for the $V_{eff}(r)$ in the SFR model for
angular momentum $J=3$, $J=4$ and $J=5$

(b) Graph for the $V_{eff}(r)$ in the GR model for
angular momentum $J=3$, $J=4$ and $J=5$

(c) Graph for the $V_{eff}(r)$ in GR (blue line)
and SFR (red line) for angular momentum $J=1$

(d) This is a graph for the $V_{eff}(r)$ in GR (blue line) and SFR (red line)
for angular momentum $J=5$

Fig. 3: These are the graphs for the $V_{eff}(r)$ in the GR and SFR models for
various angular momenta

In Fig. 3 we observe the phase portraits associated with the effective
potentials depicted above. These phase portraits reflect the existence of an
energy barrier whose precise height depends on the value of the angular
momentum. Energies below this barrier height result in reflection from the
outside and trapping from the inside. On the other hand, energies higher than
those of the barrier result in reaching the Schwarzschild radius (if the
particle is coming from the outside) or reaching infinity (if the particle is
moving outward from the inside). The latter figure demonstrates the
differences between the two phase portraits which are quantitative but not
qualitative.

(a) This is a phase plot for the radial geodesics of
SFR, representing the trajectories in $(r,\dot{r})$ space

(b) This is a phase plot for the radial geodesics of GR representing the
trajectories in $(r,\dot{r})$ space

Fig. 4: These are phase plots for the radial geodesics of the GR and SFR
models Fig. 5: This is a comparison between the radial phase portraits of the
GR (blue) and SFR (yellow) models

##  V Deflection Angle

In this section, we will deal with the deflection angle of the SFR model and
we will compare our findings with the corresponding ones of the GR model. In
this consideration, we take into account photons that pass close to a central
mass M. From Eq. (41) for photons, we put $\epsilon=0$ and we get:

$\frac{\dot{r}^{2}}{J^{2}}+\frac{f}{r^{2}}+\frac{2\tilde{A_{0}}f^{1/2}}{Jb}=\frac{1}{b^{2}}$
(53)

where $b=J/\mathcal{E}_{R}$ is a composite constant formed by the ratio of the
angular momentum $J$ divided by the energy of the particle moving along the
geodesic $\mathcal{E}_{R}$.

By using the Leibniz chain-rule
$\dot{\phi}=\frac{d\phi}{d\tau}=\frac{d\phi}{dr}\frac{dr}{d\tau}=\frac{d\phi}{dr}\dot{r}$
with the relations (35) and (53) we have:

$\frac{\dot{r}^{2}}{\dot{\phi}^{2}}=r^{4}\left(\frac{1}{b^{2}}-\frac{f}{r^{2}}-\frac{2\tilde{A_{0}}f^{1/2}}{Jb}\right)$
(54)

After some rearrangements we find:

$\frac{d\phi}{dr}=\frac{1}{r^{2}}\left(\frac{1}{b^{2}}-\frac{f}{r^{2}}-\frac{2\tilde{A_{0}}f^{1/2}}{Jb}\right)^{-1/2}$
(55)

The deflection angle is calculated by the integration of (55):

$\Delta\phi_{SFR}=2\int_{r_{1}}^{\infty}\frac{dr}{r^{2}}\left[\frac{1}{b^{2}}-\frac{1}{r^{2}}\left(1-\frac{2GM}{r}\right)-\frac{2\tilde{A_{0}}}{Jb}\left(1-\frac{2GM}{r}\right)^{1/2}\right]^{-1/2}$
(56)

where we have used $f=1-\frac{2GM}{r}$.
We perform a change of variables in the integral of Eq. (56):

$\Delta\phi_{SFR}=2\int^{w_{1}}_{0}dw\left[1-w^{2}\left(1-\frac{2GM}{b}w\right)-2a\left(1-\frac{2GM}{b}w\right)^{1/2}\right]^{-1/2}$
(57)

where we have set $w=\frac{b}{r}$ and $a=\frac{\tilde{A_{0}}b}{J}$.
If we expand the integral in powers of $\frac{2GM}{b}$ and $a$ we find:

$\Delta\phi_{SFR}\approx
2\int_{0}^{w_{1}}dw\frac{1+\frac{GM}{b}w}{\left[\left(1-2a\right)+\frac{2GM}{b}w-w^{2}\right]^{1/2}}$
(58)

where we have omitted second order terms.
By evaluating the integral in Eq. (58) we get (see Appendix 3):

$\Delta\phi_{SFR}=\pi+\frac{4GM}{b}\frac{1-a}{\sqrt{1-2a}}$ (59)

The deflection angle $\delta\phi_{SFR}$ can be found as:

$\displaystyle\delta\phi_{SFR}=\Delta\phi_{SFR}-\pi\Rightarrow$
$\displaystyle\delta\phi_{SFR}=\frac{4GM}{b}\frac{1-a}{\sqrt{1-2a}}$ (60)

If we expand Eq. (60) in powers of $a=\frac{\tilde{A}_{0}b}{J}$ the deflection
angle can be written as:

$\delta\phi_{SFR}\approx\left(1+\frac{a^{2}}{2}\right)\frac{4GM}{b}$ (61)

The deflection angle $\delta\phi$ of GR [4, 71] is given by:

$\delta\phi_{GR}=\frac{4GM}{b}$ (62)

Therefore, we notice that the deflection angle of SFR includes a small
additional Randers contribution term $a$ which shows a small deviation from GR
because $|\tilde{A_{0}}|\ll 1$. We can see from Eq. (60) that:

$\lim_{\tilde{A_{0}}\rightarrow 0}\delta\phi_{SFR}=\delta\phi_{GR}$ (63)

The small difference of the deflection angle of the SFR model from the GR one
can plausibly be attributed to the Lorentz violations [7] or on the small
amount of energy which is added to the gravitational potential of SFR.

S. S. Shapiro et al. [3] have shown that the deflection angle for a light ray
is:

$\theta\simeq\frac{(1+\gamma)GM}{c^{2}b}(1+cos\phi)$ (64)

where M is the mass of the Sun, G is the gravitational constant, $\phi$ is the
angle between the source and the Sun for an observer on Earth and $\gamma$ is
the PPN parameter that characterizes the contribution of space curvature to
gravitational deflection and is estimated to be $\gamma_{Shap}=0.9998\pm
0.0004$. We compare this observational result of Shapiro et al. with Eq. (55)
of our model in order to obtain some constraints for the parameter $\alpha$.
If we do this we find:

$2\left(1+\frac{\alpha^{2}}{2}\right)\simeq 1+\gamma_{Shap}$ (65)

where we have taken $\phi\simeq 0$ for an infinitesimal neighbourhood of a
fiducial geodesic. By solving for $\alpha$ we get:

$\alpha\simeq\pm\sqrt{-1+\gamma_{Shap}}$ (66)

which gives the values for $\alpha\simeq\pm 0.0141421$. This result provides
an interesting relation between the deflection angle of a light ray in the SFR
model and the approach to the deflection angle used in [3].

Remark 3: By considering the following relation, we can connect the
geometrical concept of the curvature $\kappa_{\phi}=\frac{d\phi}{d\tau}$ of a
path with the deflection angle $\delta\phi$ in the following way:

$\displaystyle\dot{\phi}=\frac{d\phi}{d\tau}$
$\displaystyle=\frac{d\phi}{dr}\frac{dr}{d\tau}=\frac{d\phi}{dr}\dot{r}\Rightarrow$
(67) $\displaystyle\kappa_{\phi}=\frac{d\phi}{dr}\dot{r}\Rightarrow$ (68)
$\displaystyle\delta\phi=\int\frac{\kappa_{\phi}}{\dot{r}}dr$ (69)

This form of curvature can be called deflection curvature.

##  VI Conclusions & Future Challenges

In this article, we investigated the analytic form of the geodesics of the
model SFR which was introduced in previous works [1, 2]. A dynamical analysis
was presented based on the energy and angular momentum of a particle along of
geodesics (null or timelike) of the SFR spacetime. Comparisons between the SFR
and GR were provided. We found that there is a small deviation from the GR
model which is due to the dynamical term $A_{\gamma}(x)$. We also formulated
and studied an effective potential of our model and we compared the one of the
SFR case once again with the effective potential of GR attributing the small
but discernible differences to the specific structure of (and perturbation
incorporated within) the SFR spacetime. The relevant differences in the
trajectories were illustrated both in the evolution over the time-variable
$\tau$ and in the $(x,y)$ plane. In addition, we calculated the deflection
angle for the SFR spacetime and we compared with the corresponding one of GR.
The result is a small difference of the SFR model from GR, it is possibly
caused by Lorentz violations or by the small amount of energy which is added
to the gravitational potential of SFR spacetime.

In addition, we found an interesting relation between the observational data
of [3] and the SFR model’s prediction for the deflection angle of a null
geodesic. Finally, we presented an application which relates the covector
field $A_{\gamma}(x)$ of the SFR model with the Newtonian gravitational
potential.

It is important to note that this work opens a number of interesting
directions of further study for the future. On the one hand, the traditional
assumption of $\theta=\pi/2$ made over here is clearly a restrictive one that
simplifies the equations of motion automatically satisfying the dynamics for
the angular variable $\theta$ with the latter being at steady state. However,
more generally, one can straightforwardly envision scenarios where this
condition is no longer satisfied. It is then of interest to explore if one
starts in the vicinity of $\pi/2$ whether one stays in that neighborhood or
perhaps if one deviates away from this steady state and how the associated
dynamics of the full 4-degree-of-freedom space is accordingly explored.
Another aspect that is also worth further exploring is that of the small
amplitude covector deviation from the General Relativity standard model. Here,
we have limited our considerations to the realm of associated small amplitude
perturbations (where leading order expansions of the field would suffice).
However, it would also be of interest to explore the situation when one
gradually deviates from the realm of this approximation as well. In addition,
applications of geodesics of the SFR model can be pursued for more concrete
cosmological studies such as, e.g., for the case of the S2 stars orbiting the
black hole in Sagittarius A* in which the geodesics of the star are perturbed
from the classical Keplerian orbits because of the distribution of stellar
remnants. Indeed, our hope is that this work may pave the way towards testing
the Schwarzschild-Finsler-Randers gravitational model which incorporates
features going beyond the standard Riemannian geometry of spacetime. In this
vein, some of the above topics are presently under consideration and
associated results will be presented in future publications.

##  VII Acknowledgements

The authors would like to thank Prof. Gary Gibbons for his valuable
suggestions.

##  Appendix A Distinguished connection on $TM$.

In this work, we consider a distinguished connection ($d-$connection) ${D}$ on
$TM$ [69, 70]. This is a linear connection with coefficients
$\\{\Gamma^{A}_{BC}\\}=\\{L^{\mu}_{\nu\kappa},L^{\alpha}_{\beta\kappa},C^{\mu}_{\nu\gamma},C^{\alpha}_{\beta\gamma}\\}$
which preserves by parallelism the horizontal and vertical distributions:

$\displaystyle{D_{\delta_{\kappa}}\delta_{\nu}=L^{\mu}_{\nu\kappa}(x,y)\delta_{\mu}}\quad$
$\displaystyle,\quad
D_{\dot{\partial}_{\gamma}}\delta_{\nu}=C^{\mu}_{\nu\gamma}(x,y)\delta_{\mu}$
(70)
$\displaystyle{D_{\delta_{\kappa}}\dot{\partial}_{\beta}=L^{\alpha}_{\beta\kappa}(x,y)\dot{\partial}_{\alpha}}\quad$
$\displaystyle,\quad
D_{\dot{\partial}_{\gamma}}\dot{\partial}_{\beta}=C^{\alpha}_{\beta\gamma}(x,y)\dot{\partial}_{\alpha}$
(71)

From the above conditions, the definitions for partial covariant
differentiation follow immediately, e.g. for $X\in TTM$ the expression for the
covariant h-derivative is:

$X^{A}_{|\nu}\equiv D_{\nu}\,X^{A}\equiv\delta_{\nu}X^{A}+L^{A}_{B\nu}X^{B}$
(72)

and for the covariant v-derivative:

$X^{A}|_{\beta}\equiv
D_{\beta}\,X^{A}\equiv\dot{\partial}_{\beta}X^{A}+C^{A}_{B\beta}X^{B}$ (73)

The $d-$connection is metric-compatible when we have:

$D_{\kappa}\,g_{\mu\nu}=0,\quad D_{\kappa}\,v_{\alpha\beta}=0,\quad
D_{\gamma}\,g_{\mu\nu}=0,\quad D_{\gamma}\,v_{\alpha\beta}=0$ (74)

A $d-$connection can be uniquely defined when the following conditions are
satisfied:

  * •

The $d-$connection is metric compatible

  * •

Coefficients
$L^{\mu}_{\nu\kappa},L^{\alpha}_{\beta\kappa},C^{\mu}_{\nu\gamma},C^{\alpha}_{\beta\gamma}$
depend solely on the quantities $g_{\mu\nu}$, $v_{\alpha\beta}$ and
$N^{\alpha}_{\mu}$

  * •

Coefficients $L^{\mu}_{\kappa\nu}$ and $C^{\alpha}_{\beta\gamma}$ are
symmetric on the lower indices, i.e.
$L^{\mu}_{[\kappa\nu]}=C^{\alpha}_{[\beta\gamma]}=0$

We use the symbol $\mathcal{D}$ instead of $D$ for a connection satisfying
these conditions. We call $\mathcal{D}$ a canonical and distinguished
$d-$connection. The coefficients of this connection are

$\displaystyle L^{\mu}_{\nu\kappa}$
$\displaystyle=\frac{1}{2}g^{\mu\rho}\left(\delta_{k}g_{\rho\nu}+\delta_{\nu}g_{\rho\kappa}-\delta_{\rho}g_{\nu\kappa}\right)$
(75) $\displaystyle L^{\alpha}_{\beta\kappa}$
$\displaystyle=\dot{\partial}_{\beta}N^{\alpha}_{\kappa}+\frac{1}{2}v^{\alpha\gamma}\left(\delta_{\kappa}v_{\beta\gamma}-v_{\delta\gamma}\,\dot{\partial}_{\beta}N^{\delta}_{\kappa}-v_{\beta\delta}\,\dot{\partial}_{\gamma}N^{\delta}_{\kappa}\right)$
(76) $\displaystyle C^{\mu}_{\nu\gamma}$
$\displaystyle=\frac{1}{2}g^{\mu\rho}\dot{\partial}_{\gamma}g_{\rho\nu}$ (77)
$\displaystyle C^{\alpha}_{\beta\gamma}$
$\displaystyle=\frac{1}{2}v^{\alpha\delta}\left(\dot{\partial}_{\gamma}v_{\delta\beta}+\dot{\partial}_{\beta}v_{\delta\gamma}-\dot{\partial}_{\delta}v_{\beta\gamma}\right)$
(78)

Curvatures and torsions on $TM$ are defined by the linear maps:

$\mathcal{R}(X,Y)Z=[\mathcal{D}_{X},\mathcal{D}_{Y}]Z-\mathcal{D}_{[X,Y]}Z$
(79)

and

$\mathcal{T}(X,Y)=\mathcal{D}_{X}Y-\mathcal{D}_{Y}X-[X,Y]$ (80)

where $X,Y,Z\in TTM$. We use the following definitions for the curvature
components [69, 70]:

$\displaystyle\mathcal{R}(\delta_{\lambda},\delta_{\kappa})\delta_{\nu}=R^{\mu}_{\nu\kappa\lambda}\delta_{\mu}$
(81)
$\displaystyle\mathcal{R}(\delta_{\lambda},\delta_{\kappa})\dot{\partial}_{\beta}=R^{\alpha}_{\beta\kappa\lambda}\dot{\partial}_{\alpha}$
(82)
$\displaystyle\mathcal{R}(\dot{\partial}_{\gamma},\delta_{\kappa})\delta_{\nu}=P^{\mu}_{\nu\kappa\gamma}\delta_{\mu}$
(83)
$\displaystyle\mathcal{R}(\dot{\partial}_{\gamma},\delta_{\kappa})\dot{\partial}_{\beta}=P^{\alpha}_{\beta\kappa\gamma}\dot{\partial}_{\alpha}$
(84)
$\displaystyle\mathcal{R}(\dot{\partial}_{\delta},\dot{\partial}_{\gamma})\delta_{\nu}=S^{\mu}_{\nu\gamma\delta}\delta_{\mu}$
(85)
$\displaystyle\mathcal{R}(\dot{\partial}_{\delta},\dot{\partial}_{\gamma})\dot{\partial}_{\beta}=S^{\alpha}_{\beta\gamma\delta}\dot{\partial}_{\alpha}$
(86)

In addition, we use the following definitions for the torsion components:

$\displaystyle\mathcal{T}(\delta_{\kappa},\delta_{\nu})=$
$\displaystyle\mathcal{T}^{\mu}_{\nu\kappa}\delta_{\mu}+\mathcal{T}^{\alpha}_{\nu\kappa}\dot{\partial}_{\alpha}$
(87) $\displaystyle\mathcal{T}(\dot{\partial}_{\beta},\delta_{\nu})=$
$\displaystyle\mathcal{T}^{\mu}_{\nu\beta}\delta_{\mu}+\mathcal{T}^{\alpha}_{\nu\beta}\dot{\partial}_{\alpha}$
(88)
$\displaystyle\mathcal{T}(\dot{\partial}_{\gamma},\dot{\partial}_{\beta})=$
$\displaystyle\mathcal{T}^{\mu}_{\beta\gamma}\delta_{\mu}+\mathcal{T}^{\alpha}_{\beta\gamma}\dot{\partial}_{\alpha}$
(89)

From (81), the h-curvature tensor of the $d-$connection in the adapted basis
and the corresponding h-Ricci tensor read:

$\displaystyle
R^{\mu}_{\nu\kappa\lambda}=\delta_{\lambda}L^{\mu}_{\nu\kappa}-\delta_{\kappa}L^{\mu}_{\nu\lambda}+L^{\rho}_{\nu\kappa}L^{\mu}_{\rho\lambda}-L^{\rho}_{\nu\lambda}L^{\mu}_{\rho\kappa}+C^{\mu}_{\nu\alpha}\Omega^{\alpha}_{\kappa\lambda}$
(90) $\displaystyle
R_{\mu\nu}=R^{\kappa}_{\mu\nu\kappa}=\delta_{\kappa}L^{\kappa}_{\mu\nu}-\delta_{\nu}L^{\kappa}_{\mu\kappa}+L^{\rho}_{\mu\nu}L^{\kappa}_{\rho\kappa}-L^{\rho}_{\mu\kappa}L^{\kappa}_{\rho\nu}+C^{\kappa}_{\mu\alpha}\Omega^{\alpha}_{\nu\kappa}$
(91)

From (86), the v-curvature tensor of the $d-$connection in the adapted basis
and the corresponding v-Ricci tensor are:

$\displaystyle S^{\alpha}_{\beta\gamma\delta}$
$\displaystyle=\dot{\partial}_{\delta}C^{\alpha}_{\beta\gamma}-\dot{\partial}_{\gamma}C^{\alpha}_{\beta\delta}+C^{\epsilon}_{\beta\gamma}C^{\alpha}_{\epsilon\delta}-C^{\epsilon}_{\beta\delta}C^{\alpha}_{\epsilon\gamma}$
(92) $\displaystyle S_{\alpha\beta}$
$\displaystyle=S^{\gamma}_{\alpha\beta\gamma}=\dot{\partial}_{\gamma}C^{\gamma}_{\alpha\beta}-\dot{\partial}_{\beta}C^{\gamma}_{\alpha\gamma}+C^{\epsilon}_{\alpha\beta}C^{\gamma}_{\epsilon\gamma}-C^{\epsilon}_{\alpha\gamma}C^{\gamma}_{\epsilon\beta}$
(93)

The generalized Ricci scalar curvature in the adapted basis is:

$\mathcal{R}=g^{\mu\nu}R_{\mu\nu}+v^{\alpha\beta}S_{\alpha\beta}=R+S$ (94)

where

$\displaystyle R=g^{\mu\nu}R_{\mu\nu}\quad,\quad
S=v^{\alpha\beta}S_{\alpha\beta}$ (95)

##  Appendix B Field equations of the model.

A Hilbert-like action on $TM$ can be defined as

$K=\int_{\mathcal{N}}d^{8}\mathcal{U}\sqrt{|\mathcal{G}|}\,\mathcal{R}+2\kappa\int_{\mathcal{N}}d^{8}\mathcal{U}\sqrt{|\mathcal{G}|}\,\mathcal{L}_{M}$
(96)

for some closed subspace $\mathcal{N}\subset TM$, where $|\mathcal{G}|$ is the
absolute value of the metric determinant, $\mathcal{L}_{M}$ is the Lagrangian
of the matter fields, $\kappa$ is a constant and

$d^{8}\mathcal{U}=\mathrm{d}x^{0}\wedge\ldots\wedge\mathrm{d}x^{3}\wedge\mathrm{d}y^{4}\wedge\ldots\wedge\mathrm{d}y^{7}$
(97)

Variation with respect to $g_{\mu\nu}$, $v_{\alpha\beta}$ and
$N^{\alpha}_{\kappa}$ leads to the following field equations [40]:

$\displaystyle\overline{R}_{\mu\nu}-\frac{1}{2}({R}+{S})\,{g_{\mu\nu}}+\left(\delta^{(\lambda}_{\nu}\delta^{\kappa)}_{\mu}-g^{\kappa\lambda}g_{\mu\nu}\right)\left(\mathcal{D}_{\kappa}\mathcal{T}^{\beta}_{\lambda\beta}-\mathcal{T}^{\gamma}_{\kappa\gamma}\mathcal{T}^{\beta}_{\lambda\beta}\right)=T_{\mu\nu}$
(98) $\displaystyle
S_{\alpha\beta}-\frac{1}{2}({R}+{S})\,{v_{\alpha\beta}}+\left(v^{\gamma\delta}v_{\alpha\beta}-\delta^{(\gamma}_{\alpha}\delta^{\delta)}_{\beta}\right)\left(\mathcal{D}_{\gamma}C^{\mu}_{\mu\delta}-C^{\nu}_{\nu\gamma}C^{\mu}_{\mu\delta}\right)=Y_{\alpha\beta}$
(99) $\displaystyle
g^{\mu[\kappa}\dot{\partial}_{\alpha}L^{\nu]}_{\mu\nu}+2\mathcal{T}^{\beta}_{\mu\beta}g^{\mu[\kappa}C^{\lambda]}_{\lambda\alpha}=\mathcal{Z}^{\kappa}_{\alpha}$
(100)

where

$\mathcal{T}_{\nu\beta}^{\alpha}=\dot{\partial}_{\beta}N_{\nu}^{\alpha}-L_{\beta\nu}^{\alpha}$
(101)

are torsion components, where $L_{\beta\nu}^{\alpha}$ is defined in (19). From
the form of (10) it follows that $\sqrt{|\mathcal{G}|}=\sqrt{-g}\sqrt{-v}$,
with $g,v$ the determinants of the metrics $g_{\mu\nu},v_{\alpha\beta}$
respectively.

##  Appendix C Calculation of the deflection angle

We begin the calculation from Eq. (57)

$\Delta\phi_{SFR}=2\int^{w_{1}}_{0}dw\left[1-w^{2}\left(1-\frac{2GM}{b}w\right)-2a\left(1-\frac{2GM}{b}w\right)^{1/2}\right]^{-1/2}$
(102)
$\displaystyle\Delta\phi_{SFR}=2\int_{0}^{w_{1}}dw\left(1-\frac{2GM}{b}w\right)^{-1/2}\left[\left(1-\frac{2GM}{b}w\right)^{-1}-w^{2}-2a\left(1-\frac{2GM}{b}w\right)^{-1/2}\right]^{-1/2}\Rightarrow$
$\displaystyle\Delta\phi_{SFR}\approx
2\int_{0}^{w_{1}}dw\left(1+\frac{GM}{b}w\right)\left[\left(1+\frac{2GM}{b}w\right)-w^{2}-2a\right]^{-1/2}\Rightarrow$
$\displaystyle\Delta\phi_{SFR}\approx
2\int_{0}^{w_{1}}dw\frac{1+\frac{GM}{b}w}{\left[\left(1+\frac{2GM}{b}w\right)-w^{2}-2a\right]^{1/2}}\Rightarrow$
$\displaystyle\Delta\phi_{SFR}\approx
2\int_{0}^{w_{1}}dw\frac{1+\frac{GM}{b}w}{\left[\left(1-2a\right)+\frac{2GM}{b}w-w^{2}\right]^{1/2}}$
(103)

In order to find $w_{1}$ we solve the following equation from the denominator:

$\left(1-2a\right)+\frac{2GM}{b}w-w^{2}=0$ (104)

and we get:

$w_{1}=\frac{GM}{b}+\sqrt{\left(\frac{GM}{b}\right)^{2}+(1-2a)}$ (105)

which is the positive root of the denominator. The solution for the integral
in Eq. (103) is:

$\Delta\phi_{SFR}=\pi+\frac{1}{\sqrt{1-2a}}\frac{2GM}{b}+\frac{2GM}{b}\sqrt{1-2a}$
(106)

where we omit terms $O((\frac{2GM}{b})^{2})$, given their smallness. Hence, we
find:

$\Delta\phi_{SFR}=\pi+\frac{4GM}{b}\frac{1-a}{\sqrt{1-2a}}$ (107)

## References

  * [1] A. Triantafyllopoulos, S. Basilakos, E. Kapsabelis and P. C. Stavrinos, Schwarzschild-like solutions in Finsler–Randers gravity, Eur. Phys. J. C 80 (2020) no.12, 1200
  * [2] E. Kapsabelis, A. Triantafyllopoulos, S. Basilakos and P. C. Stavrinos, Applications of the Schwarzschild–Finsler–Randers model, Eur. Phys. J. C 81 (2021) no.11, 990
  * [3] S. S. Shapiro, J. L. Davis, D. E. Lebach, J. S. Gregory, Measurements of the solar gravitational deflection of radio waves using geodetic very-long-baseline interferometry data, 1979-1999, Phys. Rev. Lett. 92 (2004) no.12, 121101.
  * [4] J. Hartle, Gravity: An Introduction to Einstein’s General Relativity. Pearson Education Inc, Addison Wesley: San Francisco (2002)
  * [5] G. S. Asanov and P. C. Stavrinos, Finslerian deviations of Geodesics over tangent bundle, Reports on Mathematical Physics 30 (1991) no.1, 63–69
  * [6] P. C. Stavrinos and S. Ikeda, Some connections and variational principle to the Finslerian scalar-tensor theory of gravitation, Reports on Mathematical Physics 44 (1999) no.1–2, 221–230
  * [7] V. A. Kostelecký, Riemann-Finsler geometry and Lorentz-violating kinematics, Phys. Lett. B 701 (2011) 137–143
  * [8] V. Alan Kostelecký, N. Russell and R. Tso, Bipartite Riemann–Finsler geometry and Lorentz violation, Phys. Lett. B 716 (2012) 470–474
  * [9] P. Stavrinos, Weak Gravitational Field in Finsler-Randers Space and Raychaudhuri Equation, Gen. Rel. Grav. 44 (2012) 3029–3045
  * [10] E. Minguzzi, Raychaudhuri equation and singularity theorems in Finsler spacetimes, Class. Quant. Grav. 32 (2015) no.18, 185008
  * [11] J. Foster and R. Lehnert, Classical-physics applications for Finsler $b$ space, Phys. Lett. B 746 (2015) 164–170
  * [12] V. Antonelli, L. Miramonti and M. D. C. Torri, Neutrino oscillations and Lorentz Invariance Violation in a Finslerian Geometrical model, Eur. Phys. J. C 78 (2018) no.8, 667
  * [13] B. R. Edwards and V. A. Kostelecký, Riemann–Finsler geometry and Lorentz-violating scalar fields, Phys. Lett. B 786 (2018) 319–326
  * [14] S. Ikeda, E. N. Saridakis, P. C. Stavrinos and A. Triantafyllopoulos, Cosmology of Lorentz fiber-bundle induced scalar-tensor theories, Phys. Rev. D 100 (2019) no.12, 124035
  * [15] J. J. Relancio and S. Liberati, Constraints on the deformation scale of a geometry in the cotangent bundle, Phys. Rev. D 102 (2020) no.10, 104025
  * [16] P. Stavrinos and C. Savvopoulos, Dark Gravitational Field on Riemannian and Sasaki Spacetime, Universe 6 (2020) no.9, 138
  * [17] A. P. Kouretsis, M. Stathakopoulos and P. C. Stavrinos, The General Very Special Relativity in Finsler Cosmology, Phys. Rev. D 79 (2009) 104011
  * [18] N. E. Mavromatos, S. Sarkar and A. Vergou, Stringy Space-Time Foam, Finsler-like Metrics and Dark Matter Relics, Phys. Lett. B 696 (2011) 300–304
  * [19] P. Stavrinos and S. I. Vacaru, Broken Scale Invariance, Gravity Mass, and Dark Energy in Modified Einstein Gravity with Two Measure Finsler Like Variables, Universe 7 (2021) no.4, 89
  * [20] S. Konitopoulos, E. N. Saridakis, P. C. Stavrinos and A. Triantafyllopoulos, Dark gravitational sectors on a generalized scalar-tensor vector bundle model and cosmological applications, Phys. Rev. D 104 (2021) no.6, 064018
  * [21] R. Hama, T. Harko and S. V. Sabau, Dark energy and accelerating cosmological evolution from osculating Barthel–Kropina geometry, Eur. Phys. J. C 82 (2022) no.4, 385
  * [22] G. W. Gibbons, J. Gomis and C. N. Pope, General very special relativity is Finsler geometry, Phys. Rev. D 76 (2007) 081701
  * [23] J. Skakala and M. Visser, Bi-metric pseudo-Finslerian spacetimes, J. Geom. Phys. 61 (2011) 1396–1400
  * [24] R. Gallego Torrome, P. Piccione and H. Vitorio, On Fermat’s principle for causal curves in time oriented Finsler spacetimes, J. Math. Phys. 53 (2012), 123511
  * [25] P. Stavrinos, O. Vacaru and S. I. Vacaru, Modified Einstein and Finsler like theories on tangent Lorentz bundles, Int. J. Mod. Phys. D 23 (2014) no.11, 1450094
  * [26] A. Fuster and C. Pabst, Finsler pp-waves, Phys. Rev. D 94 (2016) no.10, 104072
  * [27] N. Voicu, Volume forms for time orientable Finsler spacetimes, J. Geom. Phys. 112 (2017), 85–94
  * [28] M. Hohmann, C. Pfeifer and N. Voicu, Finsler gravity action from variational completion, Phys. Rev. D 100 (2019) no.6, 064035
  * [29] D. Colladay and L. Law, Spontaneous CPT breaking and fermion propagation in the Schwarzschild geometry, Phys. Lett. B 795 (2019) 457–461
  * [30] E. Caponio and A. Masiello, On the analyticity of static solutions of a field equation in Finsler gravity, Universe 6 (2020) no.4, 59
  * [31] R. Hama, T. Harko, S. V. Sabau and S. Shahidi, Cosmological evolution and dark energy in osculating Barthel-Randers geometry, Eur. Phys. J. C 81 (2021) 742
  * [32] X. Li, X. Zhang and H. N. Lin, Probing a Finslerian Schwarzschild black hole with the orbital precession of Sagittarius A*, Phys. Rev. D 106 (2022) no.6, 064043
  * [33] G. Randers, On an Asymmetrical Metric in the Four-Space of General Relativity, Phys. Rev. 59 (1941) no.2, 195–199
  * [34] J. I. Horváth, A Geometrical Model for the Unified Theory of Physical Fields, Phys. Rev. 80 (1950) 901
  * [35] J. I. Horváth, A. Moór, Entwicklung einer einheitlichen Feldtheorie begründet auf die Finslersche Geometrie, Z. Physik 131 (1952) 544–570
  * [36] Y. Takano, Theory of Fields in Finsler Spaces. I, Prog. Theor. Phys. 40 (1968) no.5, 1159–-1180
  * [37] S. Ikeda, On the theory of fields in Finsler spaces, J. Math. Phys. 22 (1981) 1215
  * [38] G. S. Asanov, Gravitational field equations based on Finsler geometry, Foundations of Physics 13 (1983) no.5, 501–527
  * [39] A. Triantafyllopoulos and P. C. Stavrinos, Weak field equations and generalized FRW cosmology on the tangent Lorentz bundle, Class. Quant. Grav. 35 (2018) no.8, 085011
  * [40] A. Triantafyllopoulos, E. Kapsabelis and P. C. Stavrinos, Gravitational Field on the Lorentz Tangent Bundle: Generalized Paths and Field Equations, Eur. Phys. J. Plus 135 (2020) no.7, 557
  * [41] P. C. Stavrinos, Congruences of Fluids in a Finslerian Anisotropic Space-Time, Int. J. Theor. Phys. 44 (2005) 245-–254
  * [42] P. C. Stavrinos, A. P. Kouretsis and M. Stathakopoulos, Friedmann Robertson-Walker model in generalised metric space-time with weak anisotropy, Gen. Rel. Grav. 40 (2008) 1403–1425
  * [43] F. Girelli, S. Liberati and L. Sindoni, Planck-scale modified dispersion relations and Finsler geometry, Phys. Rev. D 75 (2007) 064015
  * [44] E. Minguzzi, Light cones in Finsler spacetime, Commun. Math. Phys. 334 (2015) no.3, 1529–1551
  * [45] M. A. Javaloyes and M. Sánchez, On the definition and examples of cones and Finsler spacetimes, RACSAM 114 (2019) no.1, 30
  * [46] J. E. G. Silva and C. A. S. Almeida, Kinematics and dynamics in a Bipartite-Finsler spacetime, Phys. Lett. B 731 (2014) 74–79
  * [47] V. A. Kostelecký and M. Mewes, Astrophysical Tests of Lorentz and CPT Violation with Photons, Astrophys. J. 689 (2008) L1
  * [48] S. I. Vacaru, Principles of Einstein-Finsler Gravity and Perspectives in Modern Cosmology, Int. J. Mod. Phys. D 21 (2012) 1250072
  * [49] P. C. Stavrinos and S. I. Vacaru, Cyclic and Ekpyrotic Universes in Modified Finsler Osculating Gravity on Tangent Lorentz Bundles, Class. Quant. Grav. 30 (2013) 055012
  * [50] M. Hohmann and C. Pfeifer, Geodesics and the magnitude-redshift relation on cosmologically symmetric Finsler spacetimes, Phys. Rev. D 95 (2017) no.10, 104021
  * [51] S. V. Chervon, Chiral Cosmological Models: Dark Sector Fields Description, Quant. Matt. 2 (2013) 71–82
  * [52] D. C. Brody, G. W. Gibbons and D. M. Meier, A Riemannian approach to Randers geodesics, J. Geom. Phys. 106 (2016) 98–101
  * [53] S. Chanda, G. W. Gibbons, P. Guha, P. Maraner and M. C. Werner, Jacobi-Maupertuis Randers-Finsler metric for curved spaces and the gravitational magnetoelectric effect, J. Math. Phys. 60 (2019) no.12, 122501
  * [54] S. Chanda and P. Guha, Eisenhart lift and Randers–Finsler formulation for scalar field theory, Eur. Phys. J. Plus 136 (2021) no.1, 66
  * [55] S. Heefer, C. Pfeifer and A. Fuster, Randers pp-waves, Phys. Rev. D 104 (2021) no.2, 024007
  * [56] H. Lou, J. Li, W. Yang, W. Feng, W. Liu, Q. Zhang, N. Zhang, Y. Qi and Y. Wu, Theoretical analysis on the Rényi holographic dark energy in the Finsler–Randers cosmology, Int. J. Mod. Phys. D 31 (2022) no.02, 2250002
  * [57] S. Angit, R. Raushan and R. Chaubey, Stability and bifurcation analysis of Finsler–Randers cosmological model, Pramana 96 (2022) no.3, 123. Indian Academy of Sciences
  * [58] P. C. Stavrinos and F. I. Diakogiannis, A Geometric anisotropic model of space-time based on Finslerian metric, Gravitation and Cosmology 10 (2004) 269–278
  * [59] Z. Chang and X. Li, Lorentz invariance violation and symmetry in Randers-Finsler spaces, Phys. Lett. B 663 (2008) 103
  * [60] S. Basilakos, A. P. Kouretsis, E. N. Saridakis and P. Stavrinos, Resembling dark energy and modified gravity with Finsler-Randers cosmology, Phys. Rev. D 88 (2013) 123510
  * [61] S. Basilakos and P. Stavrinos, Cosmological equivalence between the Finsler-Randers space-time and the DGP gravity model, Phys. Rev. D 87 (2013) no.4, 043506
  * [62] P. C. Stavrinos and M. Alexiou, Raychaudhuri equation in the Finsler–Randers space-time and generalized scalar-tensor theories, Int. J. Geom. Meth. Mod. Phys. 15 (2017) no.03, 1850039
  * [63] G. Papagiannopoulos, S. Basilakos, A. Paliathanasis, S. Savvidou and P. C. Stavrinos, Finsler–Randers cosmology: dynamical analysis and growth of matter perturbations, Class. Quant. Grav. 34 (2017) no.22, 225008
  * [64] J. E. G. Silva, R. V. Maluf and C. A. S. Almeida, A nonlinear dynamics for the scalar field in Randers spacetime, Phys. Lett. B 766 (2017) 263–267
  * [65] R. Chaubey, B. Tiwari, A. Shukla and M. Kumar, Finsler–Randers Cosmological Models in Modified Gravity Theories, Proc. Nat. Inst. Sci. India (Pt. A Phys. Sci. ) 89 (2019) no.4, 757–768
  * [66] R. Raushan and R. Chaubey, Finsler–Randers cosmology in the framework of a particle creation mechanism: a dynamical systems perspective, Eur. Phys. J. Plus 135 (2020) no.2, 228
  * [67] G. Papagiannopoulos, S. Basilakos, A. Paliathanasis, S. Pan and P. Stavrinos, Dynamics in varying vacuum Finsler–Randers cosmology, Eur. Phys. J. C 80 (2020) no.9, 816
  * [68] J. E. G. Silva, A field theory in Randers-Finsler spacetime, EPL 133 (2021) no.2, 21002
  * [69] R. Miron and M. Anastasiei, The Geometry of Lagrange Spaces: Theory and Applications, Fundam. Theor. Phys. , Springer Netherlands, Dordrecht (1994)
  * [70] S. Vacaru, P. C. Stavrinos, E. Gaburov and D. Gonta, Clifford and Riemann-Finsler structures in geometric mechanics and gravity, Differential Geometry – Dynamical Systems, Monograph 7, Geometry Balkan Press, Bucharest (2006)
  * [71] S. Carroll, Spacetime and Geometry, An Introduction to General Relativity, Pearson Education Inc, Addison Wesley (2004)