# The development of a split-tail heliosphere and the role of non-ideal
processes: a comparison of the BU and Moscow models

M. Kornbleuth11affiliationmark:  11affiliationmark:  Astronomy Department,
Boston University, Boston, MA 02215, USA<EMAIL_ADDRESS>M.
Opher11affiliationmark:  ,22affiliationmark:  22affiliationmark:  Radcliffe
Institute for Advanced Study at Harvard University, Cambridge, MA 02138, USA
I. Baliukin 33affiliationmark:  ,44affiliationmark:  ,55affiliationmark:
33affiliationmark:  Space Research Institute of Russian Academy of Sciences,
Profsoyuznaya Str. 84/32, Moscow, 117997, Russia  44affiliationmark:  Moscow
Center for Fundamental and Applied Mathematics, Lomonosov Moscow State
University, GSP-1, Leninskie Gory, Moscow, 119991, Russia  55affiliationmark:
HSE University, Moscow, Russia  M. Gkioulidou 66affiliationmark:
66affiliationmark:  Applied Physics Laboratory, Johns Hopkins University,
Laurel, MD 20723, USA  J. D. Richardson 77affiliationmark:  77affiliationmark:
Kavli Institute for Astrophysics and Space Research and Department of Physics,
Massachusetts Institute of Technology, Cambridge, MA, USA  G. P. Zank
88affiliationmark:  88affiliationmark:  Department of Space Science, The
University of Alabama in Huntsville 4 Huntsville, AL 35805, USA  A. T.
Michael66affiliationmark:  G. Tóth99affiliationmark:  99affiliationmark:
University of Michigan, Ann Arbor, MI 48109, USA  V.
Tenishev99affiliationmark:  V. Izmodenov33affiliationmark:
,44affiliationmark:  ,1010affiliationmark:  1010affiliationmark:  Institute
for Problems in Mechanics, Vernadskogo 101-1, Moscow, 119526, Russia  D.
Alexashov1010affiliationmark:  S. Fuselier1111affiliationmark:
1111affiliationmark:  Southwest Research Institute, San Antonio, TX, 78228,
USA  J. F. Drake1212affiliationmark:  1212affiliationmark:  Department of
Physics and the Institute for Physical Science and Technology, University of
Maryland, College Park, MD, USA  K. Dialynas 1313affiliationmark:
1313affiliationmark:  Office of Space Research and Technology, Academy of
Athens, 10679 Athens, Greece

###### Abstract

Global models of the heliosphere are critical tools used in the interpretation
of heliospheric observations. There are several three-dimensional
magnetohydrodynamic (MHD) heliospheric models that rely on different
strategies and assumptions. Until now only one paper has compared global
heliosphere models, but without magnetic field effects. We compare the results
of two different MHD models, the BU and Moscow models. Both models use
identical boundary conditions to compare how different numerical approaches
and physical assumptions contribute to the heliospheric solution. Based on the
different numerical treatments of discontinuities, the BU model allows for the
presence of magnetic reconnection, while the Moscow model does not. Both
models predict collimation of the solar outflow in the heliosheath by the
solar magnetic field and produce a split-tail where the solar magnetic field
confines the charged solar particles into distinct north and south columns
that become lobes. In the BU model, the ISM flows between the two lobes at
large distances due to MHD instabilities and reconnection. Reconnection in the
BU model at the port flank affects the draping of the interstellar magnetic
field in the immediate vicinity of the heliopause. Different draping in the
models cause different ISM pressures, yielding different heliosheath
thicknesses and boundary locations, with the largest effects at high
latitudes. The BU model heliosheath is 15% thinner and the heliopause is 7%
more inwards at the north pole relative to the Moscow model. These differences
in the two plasma solutions may manifest themselves in energetic neutral atom
measurements of the heliosphere.

ISM: atoms - magnetohydrodynamics (MHD) - solar wind - Sun: heliosphere

##  1 Introduction

The heliosphere is formed through the interaction between the solar wind and
the interstellar medium (ISM) as the Sun moves through its local environment.
Traditionally, the heliosphere has been regarded as having a long “comet-like”
tail extending for thousands of AU (Parker, 1961; Baranov et al., 1971;
Baranov & Malama, 1993; Pauls et al., 1995; Zank et al., 1996). This view is
being challenged by modeling and observations. Opher et al. (2015) suggest
that the heliosphere has a “croissant-like” shape instead of a “comet-like”
tail, with a short heliotail at low latitudes that only extends 300-400 AU
beyond the termination shock. A key feature of this model is that the tension
of the solar magnetic field resists stretching by the solar wind and confines
the slow solar wind in the inner heliosheath into two jets. Several papers
(Izmodenov & Alexashov, 2015, 2018; Pogorelov et al., 2015) see evidence for
this confinement, but argue for a heliosphere with a long tail that extends
for thousands of AU beyond the termination shock. One key difference between
the models of Izmodenov & Alexashov (2015) and Opher et al. (2015) is that the
latter includes non-ideal effects, such as magnetic reconnection (which the
numerical approach allows to develop). There are also MHD instabilities along
the axis of the solar magnetic field in the heliosheath in Opher et al. (2021)
which are absent in the works of Izmodenov & Alexashov (2015) and Izmodenov &
Alexashov (2020). Here we address how these differences give rise to changes
in the heliospheric structure.

With a number of different three-dimensional (3D) magnetohydrodynamic (MHD)
models being used in the community to simulate and understand the heliosphere,
it is critical to benchmark the differences between the models and to
understand how the model assumptions affect the results. Previously, Müller et
al. (2008) compared four hydrodynamic models to explore differences and
commonalities arising from different modeling strategies. They found model
predictions to be similar in the nose direction, with predicted termination
shock distances all within 7% in the nose and within 14% in the downwind
direction. Additionally, the nose heliopause location agreed within 5% for
each of the models. The main discrepancies arose with regard to the filtration
of neutral hydrogen (H) from the ISM and the strength of the interstellar bow
shock.

Although Müller et al. (2008) compared the boundary locations for different
hydrodynamic models, a comparison of the heliospheric boundaries using models
that include the solar and interstellar magnetic fields has not been conducted
until now. Opher et al. (2017) found that reconnection between the
interstellar and solar magnetic fields can lead to a twisting of the
interstellar magnetic field around the heliopause. In the model of Izmodenov &
Alexashov (2015) and Izmodenov & Alexashov (2020) (herein referred to as the
Moscow model), reconnection between the interstellar and solar magnetic fields
at the heliopause is not allowed by their numerical approach.

In this work, we compare the recently developed SHIELD model (Michael et al.,
2021) (herein referred to as the BU model), which extends the MHD model from
Opher et al. (2015) to include a kinetic treatment of neutrals, with the MHD
model from Izmodenov & Alexashov (2015) and Izmodenov & Alexashov (2020).
There are other 3D MHD models of the heliosphere (Florinski & Pogorelov, 2008;
Washimi et al., 2011; Zank et al., 2013; Pogorelov et al., 2015; Heerikhuisen
et al., 2019; Guo et al., 2019), but here we focus on the comparison of the BU
and Moscow models to understand the effects of non-ideal MHD processes on the
heliosphere. These effects are not explicitly suppressed in BU model but are
suppressed in the Moscow model. We therefore investigate the difference which
arise from the BU model including non-ideal MHD processes while the Moscow
model does not. For both models, we use the same boundary conditions. In
Section 2, we describe the two models and the boundary conditions. In Section
3, we compare the global structure of the heliosphere and the heliosheath
properties for the two models. In Section 4, we present a summary of our
results and conclusions.

##  2 Model

The two numerical models of the heliosphere used in this work are the Moscow
model (Izmodenov & Alexashov, 2015, 2020) and the BU model (Opher et al.,
2015; Michael et al., 2021). Both models are 3D kinetic-MHD models which treat
the cold solar wind plasma and hot pick-up ion plasma as a single fluid, while
treating neutral particles kinetically. Both models use the same inner
boundary conditions for the solar wind (Figure 1), an averaged solar cycle
condition from the years 1995 to 2017 and the same outer boundary conditions
for the interstellar neutral H and protons (Izmodenov & Alexashov, 2020).

###  2.1 Moscow Model

The Moscow model treats the partially ionized interstellar plasma as a two-
component gas consisting of neutral H atoms and a charged plasma consisting of
protons, electrons, and helium ions (He+). In the solar wind, the plasma
consists of protons, electrons, and alpha particles (He++).

The neutral atoms are treated kinetically, while the plasma is described via
the ideal MHD equations. Source terms for the MHD equations are calculated via
the kinetic treatment of neutrals, as integrals of the H-atom velocity
distribution are calculated using a Monte Carlo method (Malama, 1991) which
solves the kinetic equation (e.g. Izmodenov et al. (2001)). The velocity
distribution of the proton component is assumed to be locally Maxwellian,
though the method allows for a generalization to any isotropic distribution
function (see Malama et al. (2006)).

A global iteration method (Baranov & Malama, 1993) gives a self-consistent
steady-state solution. A 3D moving grid is used to fit discontinuities. A
fitting technique proposed by Godunov et al. (1979) allows for a fitting of
all the major discontinuities present in the simulation – the heliopause,
termination shock and bow shock (when applicable). A specific, non-regular
moving grid allows for an exact fitting of the termination shock and the
heliopause by increasing the resolution and the number of cells in the
vicinity of discontinuities (see Figure 2 of Izmodenov & Alexashov (2015)). In
the non-regular grid, one surface of the grid is chosen to be the heliopause,
i.e. tangential discontinuity. The mathematical conditions at the surface are
enforced to be exactly: 1) no mass flux through the surface, 2) $B_{n}$=0, and
3) balance of total pressure on the two sides of the surface. With the
termination shock and heliopause existing as grid surfaces, the Moscow model
is able to define different populations of neutral H atoms based on their
place of origin (see, e.g. Malama (1991), Baranov et al. (1998)). Optimization
of the critical weights for the Monte Carlo method are performed for each of
the four neutral H populations separately (Malama, 1990), allowing for an
increase in the efficiency of the method. Additionally, the nature of the grid
at the heliopause allows no direct communication (such as magnetic
reconnection) between the solar wind and interstellar medium, which is a key
difference between the Moscow and BU models. In this work, we interpolate the
results from the Moscow model onto a Cartesian grid with uniform resolution of
1 AU for our analysis. The cells of the original grid in the Moscow model
generally have 8 vertices, in which the gas-dynamic values are set. Each cell
is divided into 6 tetrahedrons, which densely fill the space. To determine the
values in an arbitrary point of space, linear interpolation is performed in
the corresponding tetrahedron. Using this procedure, the values for the 1 AU
Cartesian grid were obtained. The original grid for the Moscow model is shown
in Figure 2 of Izmodenov & Alexashov (2015).

The Moscow model solves the ideal MHD equations using a unipolar magnetic
field configuration. As noted in Izmodenov & Alexashov (2015), the sign of the
radial component of the solar magnetic field is not needed because the terms
in the ideal MHD equations which are responsible for the influence of the
magnetic field do not depend on its orientation. The solutions do not depend
on the orientations of the solar and interstellar magnetic fields along their
magnetic field lines. Once the solution is obtained, the polarity of the solar
magnetic field can be calculated based on the work of Barsky (1999).

The Moscow model uses the system of coordinates connected with the
interstellar flow and magnetic field vectors, as described in Izmodenov &
Alexashov (2015). The z-axis is directed toward the interstellar flow. The
x-axis is in the plane containing interstellar flow and magnetic field vectors
(BV—plane) and perpendicular to the z-axis. The direction of the x-axis is
chosen such that the projection of interstellar magnetic field to the x-axis
is negative. The y-axis completes the right-handed system of the coordinates.

###  2.2 BU Model

The BU model, also known as the Solar-wind with Hydrogen Ion Exchange and
Large-scale Dynamics (SHIELD) model, is a numerical model based on the work of
Michael et al. (2021). This model extends the MHD model of Opher et al. (2015)
to include a kinetic treatment of neutral H atoms instead of a multi-fluid
treatment of H atoms (Zank et al., 1996). The BU model uses the Space Weather
Modeling Framework (SWMF) (Tóth et al., 2005) and couples the Outer
Heliosphere (OH) and Particle Tracker (PT) components. The OH component is
based on the Block-Adaptive Tree Solar wind Roe-Type Upwind Scheme (BATS-R-US)
solver (Tóth et al., 2012), which is a three-dimensional (3D), block adaptive,
upwind finite-volume MHD code that is highly parallel. Opher et al. (2003)
adapted BATS-R-US to the outer heliosphere as the OH component, which is a
global 3D multi-fluid simulation. As a standalone component within SWMF, the
OH component is capable of treating multiple ion fluids in addition to
multiple neutral fluids (Opher et al., 2020), but here we use only a single-
ion fluid. In a multi-fluid treatment of the heliosphere, the ideal MHD
equations are solved for the ion species and Euler’s equations are solved
separately for the individual neutral H populations, with each neutral
population corresponding to a different region of the heliosphere (Zank et
al., 1996). Source terms from McNutt et al. (1998) connect the ion and neutral
fluids, where we treat the thermal solar wind ion and pick-up ions (PUIs) as a
single fluid population.

The PT component of SWMF is based on the Adaptive Mesh Particle Simulation
(AMPS) to treat neutral atoms kinetically. AMPS is a global, kinetic, 3D
kinetic particle code which solves the Boltzmann equation using a Direct
Simulation Monte Carlo method (Tenishev et al., 2021). The BU model couples
the PT component to the OH component, and in so doing the PT component is used
to solve the Boltzmann equation for neutral H atoms streaming through the
domain and only incorporates effects due to charge exchange. Based on the work
of Malama (1991), a statistical estimation of the momentum and energy source
terms is found by summing the change in these quantities from the individual
charge exchange events in a given cell over a specific time interval (Michael
et al., 2021). The neutral H atoms are the only modeled neutral population,
and are injected at the outer boundary with a Maxwell-Boltzmann distribution.

When running the BU model, the OH standalone component is used first with the
multi-fluid neutral approximation. This is done to relax the plasma to a
steady state solution in the MHD model (Opher et al., 2009; Opher et al.,
2015). Once steady state is achieved, the OH component is coupled with the PT
component and neutral atoms are modeled kinetically. The BU model cycles
between the OH and PT components as source terms are passed from the PT
component to the OH component. Approximately 148 million particles are modeled
within AMPS, and the source terms from charge exchange for the MHD equations
accrue for 5000 time steps in the PT component in order to obtain sufficient
statistics before being passed to the MHD solver to update the plasma solution
in the OH component. When a charge exchange event occurs in the PT component,
the MHD parameters from the OH component are interpolated to the location of
the particle in the event via linear interpolation. Source terms are
accumulated between couplings with the OH component.

The computational domain for both the OH and PT components is x=$\pm$1500 AU,
y=$\pm$2000 AU, and z=$\pm$2000\. The coordinate system is such that the
z-axis is parallel to the solar rotation axis and the x-axis is 5∘ above the
direction of interstellar flow, with y completing the right-handed coordinate
system. For AMPS in the PT component, the grid resolution within the
heliosphere is 4.7 AU from x=-280 AU to 560 AU, y=-500 AU to 500 AU, and
z=-380 AU to 380 AU. Beyond this region, the grid resolution gradually
decreases until it reaches 18.8 AU at the domain edges in the ISM. For BATS-R-
US in the OH component, the resolution is similar to that of Opher et al.
(2015) and Michael et al. (2021). We use a resolution of 1 AU (from x=-190 to
-110 AU, y=-100 to 100 AU, z=-85 to 85 AU) and 2 AU (from x=-190 to -110 AU,
y=-300 to 300 AU, z=-200 to 200 AU) at the nose of the heliosphere, and 4 AU
resolution (from x=-300 to 1000 AU, y=-410 to 410 AU, z=-450 to 450 AU) in the
heliosheath. Both the OH and PT components interpolate the results from the
other components when coupled using a bilinear interpolation technique and the
interpolated results are stored in cell centers.

As in the Moscow model, the BU model uses a unipolar solar magnetic field
configuration in both hemispheres, as in Opher et al. (2015). The unipolar
treatment of the solar magnetic field eliminates spurious numerical effects
due to numerical diffusion and reconnection of the solar magnetic field across
the heliospheric current sheet (Michael et al., 2018). The BU model uses a
multi-fluid description to generate the initial state, which requires defining
a criteria to identify the different regions where neutrals are created and
destroyed separated by boundaries like the heliopause and termination shock so
that one knows how to model the multiple neutral fluids. After the final
multi-fluid state is obtained, the kinetic description is used which is
independent of any criteria and the boundaries form naturally. Several
criteria, such as the iso-temperature criteria for the heliopause, are used
only in post-processing to try to best capture the iso-surface of the
heliopause. It is true that this iso-surface best identifies the forward
hemisphere of the heliosphere and is less effective in the tail region. While
preventing magnetic dissipation within the heliosheath, the unipolar magnetic
field allows reconnection to occur in the port side (defined as the left side
of the heliosphere as it moves through the ISM as seen within looking out) and
tail of the BU model. Opher et al. (2017, 2020) suggest that this reconnection
explains the draped magnetic field ahead of the heliopause as revealed by
Voyager 1 and 2 data. Using a dipolar solar magnetic field, Michael et al.
(2018) showed that reconnection across the current sheet with the interstellar
magnetic field at the heliopause leads to a draped interstellar magnetic field
that instead deviates from observations. In the BU model the solar and
interstellar magnetic field lines are able to reconnect, but in the Moscow
model no communication is allowed at the heliopause.

###  2.3 Boundary Conditions

To compare the BU and Moscow models, we use identical inner and outer boundary
conditions taken from Model 1 of Izmodenov & Alexashov (2020). A difference
between the two models is that the BU model does not include He ions or alpha
particles, but the Moscow model does. For this study we ran the Moscow model
with no He ions or alpha particles so the model comparison is valid. In the
interstellar medium, the proton density is assumed to be $n_{p,ISM}$ = 0.04
cm-3, while the neutral H atom density is $n_{H,ISM}$ = 0.14 cm-3. All of the
neutral and ionized populations in the interstellar medium are assumed to have
the same bulk velocity $v_{ISM}$ = 26.4 km/s (longitude = 75.4∘, latitude =
-5.2∘ in ecliptic J2000 coordinate system) and temperature $T_{ISM}$ = 6530 K
at the outer boundary, where there is pristine ISM not mediated by the
heliosphere. We use the interstellar magnetic field intensity and orientation
corresponding to $B_{ISM}$ = 3.75 $\mu$G and $\alpha$ = 60∘, where the
magnetic field is aligned with the hydrogen deflection plane (Lallement et
al., 2005) and $\alpha$ is the angle between the interstellar velocity and
magnetic field vectors.

For the inner boundary conditions, we use 22-year averaged solar cycle
conditions (1995-2017) as in Izmodenov & Alexashov (2020). Helio-latitudinal
variations of the solar wind density and speed are taken into account, and the
temperature is related to the speed via the Mach number. At Earth, a Mach
number of M=6.44 (corresponding to a solar wind temperature of $T_{SW}$ =
188500 K) is used. In the ecliptic plane, hourly-averaged solar wind data from
the OMNI 2 dataset is used for the density and speed. Heliolatitudinal
variations of the solar wind speed are based on analysis of interplanetary
scintillation (IPS) observations (Tokumaru et al., 2012) from 1990 to 2017.
For heliolatitudinal variations of the solar wind mass flux, SOHO/SWAN full-
sky maps of backscattered Lyman-alpha intensities are used (Quémerais et al.,
2006; Lallement et al., 2010; Katushkina et al., 2013, 2019). Data from
SOHO/SWAN are available from 1995 to the end of 2017, and an inversion
procedure is used to obtain solar wind mass flux as a function of time and
heliolatitude. For the solar magnetic field, a Parker solution is assumed for
both models, with the radial component of the magnetic field as $B_{SW}$ =
37.5 $\mu$G at 1 AU.

For the Moscow model, the inner boundary conditions are implemented at 1 AU,
while in the BU model the inner boundary conditions are implemented at 10 AU.
The boundary conditions are matched at 10 AU by extracting the solar wind
conditions from the Moscow model at that distance and implementing them in the
BU model. At 10 AU in the ecliptic plane, the plasma density is 0.067 cm-3 and
the plasma speed is 442 km/s. In Figure 1, we present the latitudinal profile
of the solar wind proton density and speed at 10 AU, which is assumed to be
independent of longitude, and the comparison of the solar wind conditions out
to 30 AU generated by the two models.

##  3 Results

###  3.1 Global Structure of the Heliosphere

One of the major questions that led to the code comparison documented in the
present paper is whether the solar magnetic field of the heliosphere is strong
enough to collimate the plasma flow into distinct lobes in the downtail
region. Both models do collimate the solar wind outflow into distinct lobes.
On the other hand, there are distinct differences in spite of both models
using the same boundary conditions. The low latitude region of the tail in the
BU model is dominated by ISM plasma while the heliopause in the Moscow model
prevents penetration of ISM into this region. The reasons for such differences
are discussed below.

Both models show collimation of the the heliosheath plasma by the solar
magnetic field (as shown by Opher et al. (2015) and Drake et al. (2015), see
also collimation in astrospheres by Golikov et al. (2017a, b) and Korolkov &
Izmodenov (2021)). Figure 2 shows a comparison of the two solutions in the
meridional plane (defined as the Y=0 AU plane in the BU model coordinate
system). Both models produce a split tail in the sense that the solar magnetic
field confines the charged solar particles into distinct north and south
columns that become lobes (Figure 3). In the BU model, the ISM flows between
the two lobes at distances larger than 300 AU from the Sun. There are
indications that instabilities that can develop along the axis of the solar
magnetic field in the BU model (Opher et al., 2021), which are not seen in the
Moscow model, are the key reason for the ISM flowing between the two lobes. In
the current paper, the BU model is in a quasi-steady state, determined by
evaluating whether the location of the termination shock and heliopause in the
nose and tail no longer change with time. The quasi-steady state model
captures the presence of time dependent features, such as instabilities, which
are absent in the Moscow model. The role of instabilities and the time
dependence of the solution will be studied in future work when the BU model is
extended to be fully time dependent.

The plasma density downstream of the termination shock both in the upwind
(nose) direction and downwind (tail) direction within 300 AU of the Sun are
similar for both the BU and Moscow models. Beyond 300 AU in the downwind
direction, the plasma densities are different in both models since the ISM
flows between the two lobes in the BU model (Figure 2). In both models, the
plasma is confined in two lobes in the tail; however, in the BU model
approximately 300 AU from the Sun, there is an enhancement of the plasma
density where solar magnetic field lines are open to the ISM. This is due to
instabilities (Opher et al., 2021) present in the BU model.

The comparison between the magnetic field intensity between the BU and Moscow
models is also shown in Figure 2. As with the plasma density, the magnetic
field in the Moscow model extends down the tail for thousands of AU
uninterrupted by the ISM. There is a clear evidence of a split in the
heliotail at high and low latitudes by the solar magnetic field in both
models. The neutral H solutions are similar for the two models.

In Figure 3, we compare the mass flux from the BU and Moscow models in the
meridional plane. Within 100 AU of the termination shock in the heliotail the
mass flux between the two models shows a similar profile. Black lines
representing the magnetic field intensity are overlaid on the mass flux, which
indicate that the mass flux is organized by the solar magnetic field, with
higher mass flux at high latitudes and lower mass flux at lower latitudes.
Beyond 100 AU from the termination shock, differences in the mass flux of the
two models arise in response to the deflection of the solar wind plasma in the
BU model due to the shortened heliotail.

###  3.2 Interstellar Medium Properties

Upstream of the heliopause in the ISM, the solution looks fairly similar in
both models, though there are differences. The differences in the ISM stem
from different treatments of the heliopause in the two models. The inclusion
of non-ideal MHD effects in the BU model causes the draping of the
interstellar magnetic field to be very different between the two models (Fig.
4). The red lines reflect reconnected magnetic field lines that originate at
the Sun and are open to the ISM, while green lines reflect non-reconnected
magnetic field lines. In the Moscow model (bottom panels) the draping of the
interstellar magnetic field reflects ideal draping around a surface with the
axis of symmetry along the direction of the undisturbed interstellar magnetic
field plane. For the BU model, the interstellar magnetic field ahead of the
heliopause twists to an azimuthal direction (top panels) close to the
heliopause. This twist is due to reconnection in the port flank of the
heliosphere (Washimi et al., 2015).

The BU model suppresses reconnection in the nose region while allowing it in
the flanks (Opher et al., 2017), consistent with recent ideas about
reconnection suppression from diamagnetic drifts. Most of the reconnection at
the nose occurs between the $B_{y}$ components of the solar and interstellar
magnetic fields. By choosing the same polarity of the ISM for the solar
magnetic field, the BU model suppresses the reconnection at the nose. The
stabilizing effect of diamagnetic drifts is an important kinetic effect missed
by MHD models, which can develop at boundaries such as the Earth’s
magnetospause or the heliopause (Swisdak et al., 2003, 2010). The stabilizing
influence of these drifts has been extensively documented with solar wind and
magnetospheric data (Phan et al., 2010, 2013). The diamagnetic drift
suppresses reconnection when the drift speed is larger than the Alfvén speed
in the reconnecting magnetic field. The jump in plasma $\beta$ (the ratio of
the plasma pressure to the magnetic pressure) across the nose of the
heliopause is much greater than in the flanks because the heliosheath $\beta$
is greater at the nose than in the flanks (Opher et al., 2017). Large-scale
reconnection is therefore suppressed in the nose but not at the flanks. The
twisting of the interstellar magnetic field to the azimuthal direction ahead
of the heliopause leads to an increase in the magnetic field intensity at
higher latitudes along the heliopause. In the BU model, there is a build-up of
the interstellar magnetic field at the nose of the heliosphere and at high
latitudes, as demonstrated in Figure 2.

Figure 5 shows the interstellar magnetic field intensity, the elevation angle
of the magnetic field ($\delta$), and the azimuthal angle of the magnetic
field ($\lambda$) in the Voyager 1 and 2 directions. The elevation angle is
given by ${\delta}=sin^{-1}(B_{N}/B)$ and the azimuthal angle is given by
${\lambda}=tan^{-1}(B_{T}/B_{R})$. We use the RTN coordinate system, where R
is radially outward from the Sun, T is parallel with the plane of the solar
equator and is positive in the direction of solar rotation, and N completes a
right-handed system. The termination shock spans approximately 2 grid cells
with 4 AU resolution in the BU model, therefore the numerical effects of
resolution makes the solar magnetic field intensity appear continuous at the
termination shock. At large distances from the heliopause, the profiles of the
interstellar magnetic field are very similar. They are not identical due to
the sub-alfvenic nature of the flow between the heliopause and the pristine
ISM. The differences at the heliopause propagate to larger distances ahead of
the heliopause. For the assumed ISM conditions, a fast bow shock is not
present in the solutions, likely attributable to neutrals mediating and
smoothing the bow shock through charge exchange with the plasma (Izmodenov et
al., 1999, 2009; Zank et al., 2013).

Figure 5 shows the magnetic field close to the heliopause for both models. For
the BU model, there is a smooth transition in the azimuthal angle ($\lambda$)
and the elevation angle ($\delta$) of the magnetic field. As described in
Opher et al. (2017), reconnection at the heliopause yields a magnetic field
configuration that appears “solar-like” due to the twisting of the magnetic
field around the heliopause (similar to the classic Parker spiral). For the
Moscow model, where magnetic reconnection is suppressed, there is a
discontinuity in both the azimuthal and elevation angles of the magnetic field
at the heliopause. The magnetic field change at the heliopause in both models
is consistent with Voyager 1 and 2 observations (Burlaga & Ness, 2014; Burlaga
et al., 2019).

The difference in the pressure at the heliopause affects the distance of the
heliospheric boundaries in the two models. As noted in Figure 6, the
heliopause is further from the Sun in the Moscow model than in the BU model.
The stronger magnetic pressure due to the twisting of the interstellar
magnetic field in the BU model leads to a higher overall pressure and
therefore more compression of the heliosphere. The twisting of the
interstellar magnetic field in the BU model that is absent in the Moscow model
primarily affects the heliosphere at high latitudes and in the heliotail.
Therefore, the heliopause in the BU model is closer to the Sun than in the
Moscow model at high latitudes (7% at the north pole). When comparing
heliopause distances here and throughout this work, we use streamlines to
trace the exact location of the boundary in the BU and Moscow models.

The magnetic pressure at high latitudes is higher outside the heliopause in
the BU model compared to the Moscow model (Figure 7). The total pressure
outside of the heliopause at high latitudes is also higher in the BU model,
giving rise to a more compressed heliosphere. Beyond 350 AU at low latitudes,
the ISM flows between the two lobes in the BU model, which leads to the lower
thermal pressure and higher magnetic pressure in the tail. Both models have
total pressures which agree with Voyager observations (Rankin et al., 2019;
Dialynas et al., 2020).

Toward the northern pole, the heliosheath is 13 AU thinner (15%) in the BU
model relative to the Moscow model even though the difference in termination
shock locations is only 1 AU. Figure 8 shows 1D cuts near the heliopause of
the pressure variation for both models in the northern pole of the
heliosphere. We shift the 1D cut for the BU model inwards by 14 AU to align
the two heliopause locations. Due to the lack of communication between the ISM
and solar wind in the Moscow model, there is a sharp jump in the thermal
pressure of the plasma and the magnetic pressure at the heliopause while the
magnetic field pressure replaces the plasma pressure so that the total
pressure is conserved. This effect is not present in the BU model. For the
Moscow model, there is a higher thermal pressure for the plasma just beyond
the heliopause as compared to the BU model, but 50 AU past the heliopause the
Moscow and BU models have the same thermal pressure. In contrast, the Moscow
model has a notably lower magnetic pressure just beyond the heliopause as
compared to the BU model, but 50 AU past the heliopause the BU model still
maintains a higher magnetic pressure by 17%. As mentioned previously, this
higher magnetic pressure present in the BU model is a consequence of the
twisting of the interstellar magnetic field due to reconnection present in the
model. This higher magnetic pressure outside of the heliosphere at high
latitudes gives rise to a higher total pressure by 13%, which leads to a
compression of the heliosphere relative to the Moscow model.

###  3.3 Heliosheath Properties

The compression of the heliosphere in the BU model relative to the Moscow
model primarily affects the heliopause (Fig. 6). In Figures 9 and 10, we
present meridional and equatorial slices, respectively, of the BU and Moscow
models showing the plasma density, speed, and temperature in the heliosheath.

Aside from the split of the heliotail present in the BU model, the properties
in the heliosheath for both models are very similar. Within 300 AU of the Sun,
the plasma density, speed, and temperature profiles in the heliosheath are
similar for both models. The splitting of the tail occurs at roughly 300 AU in
the meridional slice. The BU model flows stagnate in this region, which does
not occur in the Moscow model. This stagnation region reflects the change of
the plasma flow in response to the heliopause down the tail. Beyond 350 AU,
where the ISM splits the heliotail, the properties differ between the models.
The plasma density and temperature in the BU model reflects a mixing between
the heliosheath and ISM plasma in this region. This is also seen in the
equatorial speed profiles of the BU model in Figure 10, where turbulent mixing
is present beyond the stagnation region at 300 AU. This turbulence prevents
velocity streamline tracing from the ISM within the region of mixed ISM and
solar wind plasma; however, Figure 11 shows the intersection of reconnected
interstellar magnetic field lines penetrating the region between the lobes of
the BU model where the two ends of the field lines are connected to the ISM.

In Figure 12 we present a comparison of the plasma conditions and the magnetic
field along the Voyager 1 and 2 trajectories for the two models. While we do
not compare with Voyager data directly due to the constant solar wind boundary
conditions used in the simulations, the Voyager trajectories are valuable as
they are direct probes of heliospheric and ISM conditions by Voyager 1 and 2.
Near the poles the heliosheath thickness is notably thicker for the Moscow
model (14 AU difference in heliopause location at the northern pole); in
contrast the termination shock and heliopause in the models are at similar
locations giving similar heliosheath thicknesses (Table 1). Comparison of
these models along the Voyager lines-of-sight reveals interesting
characteristics of the heliosphere. The similar results in the Voyager 1 and 2
directions indicate that the shape of the heliotail does not have a
significant effect on the heliosheath flows in the direction of the nose. In
future work we will compare the BU and Moscow models with time dependent solar
cycle conditions corresponding to Voyager observations to compare with the
data.

##  4 Conclusions

Table 2 summarizes our results. We compared two MHD models using identical
inner and outer boundary conditions for the protons and neutral H atoms, and
investigated how the different assumptions present in each model influenced
the results. Both models produce a split-tail since the solar magnetic field
confines the solar particles into distinct north and south columns that become
lobes. In the BU model, the ISM flows between the two lobes at large
distances. MHD instabilities in the BU model (Opher et al., 2021), which are
absent in the Moscow model, affect the heliotail and create a short, split
tail in contrast with the long tail of the Moscow model. The BU model gives a
smaller heliosphere with a heliopause closer to the Sun than predicted in the
Moscow model in all directions. The communication between the solar and
interstellar magnetic fields at the heliopause in the BU model allows for
magnetic reconnection in the port flank (Opher et al., 2017) and leads to a
change in the magnetic topology, which does not occur in the Moscow model.

The BU model uses a polarity for the solar magnetic field that suppresses
reconnection at the nose. Reconnection occurs in the port flank, the main
location where reconnection should proceed since diamagnetic effects suppress
reconnection at the nose (Opher et al., 2017). Magnetic reconnection at the
heliopause in the BU model twists the interstellar magnetic field and leads to
an increase in the interstellar magnetic pressure relative to the Moscow
model.

In comparing the two models, we see a short low latitude tail in the BU model
and a long tail in the Moscow model. The question of the structure of the
heliotail is currently an active and vigorous topic of debate in the community
between several groups. This paper establishes that both the BU and Moscow
models predict magnetic confinement in the heliotail. This supports the
earlier conclusions of the split arrangement of the helitail made by Opher et
al. (2015), Izmodenov & Alexashov (2015), and Pogorelov et al. (2015), which
had still been debated. The main difference between the two models is the ISM
piercing the two lobes in the BU model, which is absent in the Moscow model.
The next key step is to pin down the driver for the splitting of the two lobes
by the ISM. The relevant driver may well be an instability (Opher et al.,
2021) related to the interaction of streaming ISM neutral H and the ionized
heliotail plasma that results in a loss of coherence of the heliospheric jets
and leads to the mixing of ISM and heliospheric tail material. This
instability is present in the BU model but not in the Moscow model (for
reasons that are not clear at the present time). Investigation of the
instability and its consequences such as reconnection and ISM-heliotail plasma
mixing are deferred to a future study.

Numerical solution strategies and grid resolutions might play an important
role in the different modeling outcomes. This is a topic that needs to be
addressed, as global heliosphere modeling suffers from the inability to model
phenomena on comparatively small scales, including reconnection, the
heliospheric current sheet, and instabilities. Therefore, additional
comparisons with other models both on the scale of the global heliosphere and
also at small scales are required to further address this.

In future work we will use time dependent MHD models to compare more realistic
solutions. The results are important since the different heliosphere
properties from these models, such as heliosheath thickness at high latitudes,
yield different ENA signatures at different energy bands. These differences
are potentially observable by the Interstellar Boundary Explorer. This will be
explored in a future study.

M.K., M.O., and J.R. were supported by NASA grant 18-DRIVE18$\\_$2-0029, Our
Heliospheric Shield, 80NSSC20K0603. The work of I.B., V.I., and D.A. was
conducted in the framework of a topic of the state assignment “Plasma” to the
Space Research Institute, Russian Academy of Sciences. The work at JHU/APL was
supported also by NASA under contracts NAS5 97271, NNX07AJ69G, and NNN06AA01C
and by subcontract at the Office for Space Research and Technology. Resources
supporting this work were provided by the NASA High-End Computing (HEC)
program through the NASA Advanced Supercomputing (NAS) Division at Ames
Research Center. The authors would like to thank the staff at NASA Ames
Research Center for the use of the Pleiades super-computer under the award
SMD-20-46872133.

## References

  * Baranov et al. (1971) Baranov, V. B., Krasnobaev, K. V., & Kulikovskii, A. G. 1971, Soviet Physics Doklady, 15, 791
  * Baranov & Malama (1993) Baranov, V. B., & Malama, Y. G. 1993, J. Geophys. Res., 98, 15
  * Baranov et al. (1998) Baranov, V. B., Izmodenov, V. V., & Malama, Y. G. 1998, J. Geophys. Res., 103, 9575. doi:10.1029/97JA03662
  * Barsky (1999) Barsky, E. V. 1999, Astronomy Letters, 25, 825
  * Burlaga & Ness (2014) Burlaga, L. F. & Ness, N. F. 2014, ApJ, 795, L19. doi:10.1088/2041-8205/795/1/L19
  * Burlaga et al. (2019) Burlaga, L. F., Ness, N. F., Berdichevsky, D. B., et al. 2019, Nature Astronomy, 3, 1007. doi:10.1038/s41550-019-0920-y
  * Dialynas et al. (2020) Dialynas, K., Galli, A., Dayeh, M. A., et al. 2020, ApJ, 905, L24. doi:10.3847/2041-8213/abcaaa
  * Drake et al. (2015) Drake, J. F., Swisdak, M., & Opher, M. 2015, ApJ, 808, L44. doi:10.1088/2041-8205/808/2/L44
  * Florinski & Pogorelov (2008) Florinski, V. & Pogorelov, N. V. 2008, Numerical Modeling of Space Plasma Flows, 385, 18
  * Florinski (2015) Florinski, V. 2015, ApJ, 813, 49. doi:10.1088/0004-637X/813/1/49
  * Gloeckler et al. (2004) Gloeckler, G., Möbius, E., Geiss, J., et al. 2004, A&A, 426, 845. doi:10.1051/0004-6361:20035768
  * Godunov et al. (1979) Godunov, S. K., Zabrodine, A. V., Ya, Ivanov M., et al. 1979, Resolution Numerique des Problems Multidimensionnels de la Dynamique des Gaz (Moscow: Editions Mir)
  * Golikov et al. (2017a) Golikov, E. A., Izmodenov, V. V., Alexashov, D. B., et al. 2017a, MNRAS, 464, 1065. doi:10.1093/mnras/stw2402
  * Golikov et al. (2017b) Golikov, E. A., Izmodenov, V. V., & Alexashov, D. B. 2017b, Journal of Physics Conference Series, 815, 012035. doi:10.1088/1742-6596/815/1/012035
  * Guo et al. (2019) Guo, X., Florinski, V., & Wang, C. 2019, ApJ, 879, 87. doi:10.3847/1538-4357/ab262b
  * Heerikhuisen et al. (2019) Heerikhuisen, J., Zirnstein, E. J., Pogorelov, N. V., et al. 2019, ApJ, 874, 76. doi:10.3847/1538-4357/ab05e3
  * Izmodenov et al. (1999) Izmodenov, V. V., Geiss, J., Lallement, R., et al. 1999, J. Geophys. Res., 104, 4731. doi:10.1029/1998JA900122
  * Izmodenov et al. (2001) Izmodenov, V. V., Gruntman, M., & Malama, Y. G. 2001, J. Geophys. Res., 106, 10681. doi:10.1029/2000JA000273
  * Izmodenov et al. (2003) Izmodenov, V., Malama, Y. G., Gloeckler, G., et al. 2003, ApJ, 594, L59. doi:10.1086/378387
  * Izmodenov et al. (2009) Izmodenov, V. V., Malama, Y. G., Ruderman, M. S., et al. 2009, Space Sci. Rev., 146, 329. doi:10.1007/s11214-009-9528-3
  * Izmodenov & Alexashov (2015) Izmodenov, V. V., & Alexashov, D. B. 2015, ApJS, 220, 32
  * Izmodenov & Alexashov (2018) Izmodenov, V. & Alexashov, D. 2018, Journal of Physics Conference Series, 1031, 012014. doi:10.1088/1742-6596/1031/1/012014
  * Izmodenov & Alexashov (2020) Izmodenov, V. V. & Alexashov, D. B. 2020, A&A, 633, L12. doi:10.1051/0004-6361/201937058
  * Katushkina et al. (2013) Katushkina, O. A., Izmodenov, V. V., Quemerais, E., et al. 2013, Journal of Geophysical Research (Space Physics), 118, 2800. doi:10.1002/jgra.50303
  * Katushkina et al. (2019) Katushkina, O., Izmodenov, V., Koutroumpa, D., et al. 2019, Sol. Phys., 294, 17. doi:10.1007/s11207-018-1391-5
  * Korolkov & Izmodenov (2021) Korolkov, S. & Izmodenov, V. 2021, MNRAS. doi:10.1093/mnras/stab1071
  * Lallement et al. (2005) Lallement, R., Quémerais, E., Bertaux, J. L., et al. 2005, Science, 307, 1447. doi:10.1126/science.1107953
  * Lallement et al. (2010) Lallement, R., Quémerais, E., Koutroumpa, D., et al. 2010, Twelfth International Solar Wind Conference, 1216, 555. doi:10.1063/1.3395925
  * Malama (1990) Malama, Yu. G. 1990, Preprint IPM RAS No, 450
  * Malama (1991) Malama, Y. G. 1991, Ap&SS, 176, 21. doi:10.1007/BF00643074
  * Malama et al. (2006) Malama, Y. G., Izmodenov, V. V., & Chalov, S. V. 2006, A&A, 445, 693. doi:10.1051/0004-6361:20053646
  * McNutt et al. (1998) McNutt, R. L., Lyon, J., & Goodrich, C. C. 1998, J. Geophys. Res., 103, 1905. doi:10.1029/97JA02143
  * Michael et al. (2018) Michael, A. T., Opher, M., & Tóth, G. 2018, ApJ, 860, 171
  * Michael et al. (2021) Michael, A. T., Opher, M., Tóth, G., et al. 2021, ApJ, 906, 37. doi:10.3847/1538-4357/abc953
  * Müller et al. (2008) Müller, H.-R., Florinski, V., Heerikhuisen, J., et al. 2008, A&A, 491, 43. doi:10.1051/0004-6361:20078708
  * Opher et al. (2003) Opher, M., Liewer, P. C., Gombosi, T. I., et al. 2003, ApJ, 591, L61. doi:10.1086/376960
  * Opher et al. (2009) Opher, M., Bibi, F. A., Toth, G., et al. 2009, Nature, 462, 1036. doi:10.1038/nature08567
  * Opher et al. (2015) Opher, M., Drake, J. F., Zieger, B., & Gombosi, T. I. 2015, ApJ, 800, L28
  * Opher et al. (2017) Opher, M., Drake, J. F., Swisdak, M., et al. 2017, ApJ, 839, L12. doi:10.3847/2041-8213/aa692f
  * Opher et al. (2020) Opher, M., Loeb, A., Drake, J., et al. 2020, Nature Astronomy, 4, 675. doi:10.1038/s41550-020-1036-0
  * Opher et al. (2021) Opher, M., Drake, J. F., Zank, G. P., et al. 2021, Submitted, https://www.researchsquare.com/article/rs-198925/v1
  * Parker (1958) Parker, E. N. 1958, ApJ, 128, 664
  * Parker (1961) Parker, E. N. 1961, ApJ, 134, 20
  * Pauls et al. (1995) Pauls, H. L., Zank, G. P., & Williams, L. L. 1995, J. Geophys. Res., 100, 21595. doi:10.1029/95JA02023
  * Phan et al. (2010) Phan, T. D., Gosling, J. T., Paschmann, G., et al. 2010, ApJ, 719, L199. doi:10.1088/2041-8205/719/2/L199
  * Phan et al. (2013) Phan, T. D., Paschmann, G., Gosling, J. T., et al. 2013, Geophys. Res. Lett., 40, 11. doi:10.1029/2012GL054528
  * Pogorelov et al. (2015) Pogorelov, N. V., Borovikov, S. N., Heerikhuisen, J., & Zhang, M. 2015, ApJ, 812, L6
  * Quémerais et al. (2006) QuéMerais, E., Lallement, R., Ferron, S., et al. 2006, Journal of Geophysical Research (Space Physics), 111, A09114. doi:10.1029/2006JA011711
  * Rankin et al. (2019) Rankin, J. S., McComas, D. J., Richardson, J. D., et al. 2019, ApJ, 883, 101. doi:10.3847/1538-4357/ab3d9d
  * Sokół et al. (2015) Sokół, J. M., Swaczyna, P., Bzowski, M., et al. 2015, Sol. Phys., 290, 2589. doi:10.1007/s11207-015-0800-2
  * Swisdak et al. (2003) Swisdak, M., Rogers, B. N., Drake, J. F., et al. 2003, Journal of Geophysical Research (Space Physics), 108, 1218. doi:10.1029/2002JA009726
  * Swisdak et al. (2010) Swisdak, M., Opher, M., Drake, J. F., et al. 2010, ApJ, 710, 1769. doi:10.1088/0004-637X/710/2/1769
  * Tenishev et al. (2021) Tenishev, V., Shou, Y., Borovikov, D., Lee, Y., Fougere, N., Michael, A., & Combi, M., J. Geophys. Res. Space Phys., 126, e2020JA028242. doi:10.1029/2020JA028242
  * Tokumaru et al. (2012) Tokumaru, M., Kojima, M., & Fujiki, K. 2012, J. Geophys. Res., 117, A06108
  * Tóth et al. (2005) Tóth, G., Sokolov, I. V., Gombosi, T. I., et al. 2005, Journal of Geophysical Research (Space Physics), 110, A12226. doi:10.1029/2005JA011126
  * Tóth et al. (2012) Tóth, G., van der Holst, B., Sokolov, I. V., et al. 2012, Journal of Computational Physics, 231, 870
  * Washimi et al. (2011) Washimi, H., Zank, G. P., Hu, Q., et al. 2011, MNRAS, 416, 1475. doi:10.1111/j.1365-2966.2011.19144.x
  * Washimi et al. (2015) Washimi, H., Zank, G. P., Hu, Q., et al. 2015, ApJ, 809, 16. doi:10.1088/0004-637X/809/1/16
  * Witte (2004) Witte, M. 2004, A&A, 426, 835. doi:10.1051/0004-6361:20035956
  * Zank et al. (1996) Zank, G. P., Pauls, H. L., Williams, L. L., et al. 1996, J. Geophys. Res., 101, 21639. doi:10.1029/96JA02127
  * Zank et al. (2013) Zank, G. P., Heerikhuisen, J., Wood, B. E., et al. 2013, ApJ, 763, 20. doi:10.1088/0004-637X/763/1/20

|  | BU (V1) |  | Moscow (V1) |  | BU (V2) |  | Moscow (V2)
---|---|---|---|---|---|---|---|---
$r_{TS}$ |  | 82 $\pm$ 4 AU |  | 80 $\pm$ 1 AU |  | 82 $\pm$ 4 AU |  | 80 $\pm$ 1 AU
$r_{HP}$ |  | 121 $\pm$ 1 AU |  | 123 $\pm$ 1 AU |  | 126 $\pm$ 1 AU |  | 124 $\pm$ 1 AU
$d_{HS}$ |  | 39 AU |  | 43 AU |  | 44 AU |  | 44 AU
Table 1: Comparison of termination shock ($r_{TS}$) and heliopause ($r_{HP}$) locations for the Voyager 1 and 2 directions between the BU and Moscow models, as well as the heliosheath thicknesses ($d_{HS}$). Error bars reflect the grid resolution at the location of the probed boundaries. For the Moscow model, we extract distances using a 1 AU interpolated grid. How does the size of the heliosphere change between the two models? |  The heliopause is closer to the Sun at high latitudes in the BU model than in the Moscow model, leading to a more compressed heliosphere at the poles.
---|---
How does the heliotail compare between the two models? |  Both models show collimation of the heliosheath plasma by the solar magnetic field. They also both produce a split tail since the solar magnetic field confines the solar particles into two distinct north and south columns that become lobes. In the BU model, reconnection allows the ISM to flow between the two lobes at large distances. The low latitude heliotail is shorter in the BU model while the Moscow model has a long, comet-like tail.
What assumptions have the largest impact on the heliospheric structure? |  Magnetic reconnection at the heliopause in the BU model leads to an increase in the magnetic pressure outside of the heliosphere, which compresses the heliosphere. Non-ideal MHD effects contribute to a split tail in the BU model.
How do the plasma properties in the heliosheath compare? |  Within 100 AU of the termination shock, the mass flux is similar between the two models, showing an organization by the solar magnetic field.
How do the plasma properties in the ISM compare? |  The ISM conditions are similar until roughly 50 AU ahead of the heliopause, where the interaction between the interstellar and solar magnetic field at the heliopause begins to have an effect on the plasma in the BU model.
Table 2: Summary of conclusions from this work. Figure 1: Left column:
Boundary conditions at 10 AU for the BU and Moscow models of the solar wind
density (black) and speed (red) as a function of heliolatitude. The inner
boundary conditions reflect 22-year averaged solar cycle data from 1995-2017.
Middle column: Solar wind plasma density out to 30 AU in the meridional plane
for the BU (top) and Moscow (bottom) models. Right column: Solar wind plasma
speed out to 30 AU in the meridional plane for the BU (top) and Moscow
(bottom) models. Figure 2: Comparison of BU (top) and Moscow (bottom)
solutions of the heliosphere in the meridional plane Included are color and
line contours of proton density [cm-3] (left), magnetic field intensity [nT]
(middle) and neutral H density [cm-3] (right).The white lines represent the
termination shock (inner) and heliopause (outer) for each model. Figure 3:
Meridional (top) and equatorial (bottom) slices of BU (left) and Moscow
(right) models showing the mass flux ($n_{p}V_{p}$) through the heliosheath.
Black lines represent the magnetic field intensity overlaid on the mass flux.
Figure 4: Draping of interstellar magnetic field lines around the heliopause
for the BU model (top) and Moscow model (bottom). Left: port-side view of the
heliosphere. Middle/Right: view from outside of the heliosphere along the
direction of ISM flow. Green lines reflect interstellar magnetic field lines
that have not undergone magnetic reconnection. Red lines reflect interstellar
magnetic that have undergone magnetic reconnection with solar magnetic field
lines. The yellow isosurfaces represent the heliopause, with the heliopause in
the BU model represented by the isosurface of lnT = 12.55. The top middle and
top right figures show reconnected interstellar magnetic field lines and non-
reconnected interstellar magnetic field lines in the BU model draped along the
nose of the heliosphere, respectively. Figures are not in the same scale.
Figure 5: 1D cuts along the Voyager 1 (top) and Voyager 2 (bottom)
trajectories for the BU (black) and Moscow (red) models. Included are the
magnetic field strength (left), azimuthal angle of the magnetic field
(middle), and the elevation angle of the magnetic field (right).The elevation
angle, $\delta$, is given by ${\delta}=sin^{-1}(B_{N}/B)$ and the azimuthal
angle, $\lambda$, is given by ${\lambda}=tan^{-1}(B_{T}/B_{R})$, where the RTN
coordinate system is the local Cartesian system centered at the spacecraft. R
is radially outward from the Sun, T is in parallel with the plane of the solar
equator and is positive in the direction of solar rotation, and N completes a
right-handed system. The heliopause is shifted to be at the same location (128
AU) for both models, represented by the dashed vertical blue line. Figure 6:
Comparison of heliospheric boundary locations in the BU and Moscow models.
Left: comparison of termination shock (inner line) and heliopause (outer line)
for the BU (black) and Moscow (red) models. The heliopause for the BU model is
defined here in post-processing as lnT = 12.55. Figure 7: Comparison of
pressures in the heliosphere and ISM for the BU (top) and Moscow (bottom)
models for the thermal pressure (left), magnetic pressure (middle), and the
combined magnetic and thermal pressures (right) in units of dynes/cm2. The
white lines represent the termination shock (inner) and heliopause (outer) for
each model.  Figure 8: 1D cuts of the thermal (left), magnetic (middle), and
the combined thermal and magnetic pressures (right) across the heliopause
(dashed black vertical line) for the BU (black) and Moscow (red) models at the
northern pole (heliolatitude = 90∘). Pressures for the BU model are shifted
inward by 14 AU to align the heliopause locations of the two models.  Figure
9: Meridional slices of plasma conditions within the heliosheath for the BU
model (top) and the Moscow model (bottom). Presented are color contours and
lines of the plasma density [cm-3] (left), plasma speed [km/s] (middle), and
plasma temperature [K] (right). The white lines represent the termination
shock (inner) and heliopause (outer) for each model. Figure 10: Same as Figure
9, but in the equatorial plane. Figure 11: Intersection of the meridional and
equatorial planes with contours of magnetic field strength for the BU model.
Grey lines reflect reconnected interstellar magnetic field lines passing
through a region of mixed ISM and solar wind plasma between the lobes of the
BU model. Left: port-view of the reconnected lines through the heliosphere.
Right: view from tail towards the nose of the reconnected lines through the
heliosphere. Figure 12: 1D cuts along the Voyager 1 (solid) and Voyager 2
(dashed) trajectory for the BU (black) and Moscow (red) models. Top row (from
left to right): plasma density, plasma speed, and plasma temperature. Bottom
row (from left to right): the radial, normal, and tangential components of the
plasma velocity.


∎

e1e-mail<EMAIL_ADDRESS>e2e-mail<EMAIL_ADDRESS>e3e-mail:
<EMAIL_ADDRESS>e4e-mail<EMAIL_ADDRESS>e5e-mail:
<EMAIL_ADDRESS>
11institutetext: Theoretical High Energy Physics Group, Department of Physics,
Institut Teknologi Bandung, Jl. Ganesha 10 Bandung, 40132, Indonesia
22institutetext: Indonesia Center for Theoretical and Mathematical Physics
(ICTMP), Institut Teknologi Bandung, Jl. Ganesha 10 Bandung, 40132, Indonesia

# Anisotropic Hartle-Thorne formalism in Rastall gravity

M. Lawrence Pattersonse1,addr1  Freddy Permana Zene2,addr1,addr2  Hadyan
Luthfan Prihadie3,addr1  Getbogi Hikmawane4,addr1  Agus Surosoe5,addr1

###### Abstract

In this work, we extend the Hartle-Thorne formalism for anisotropic pressure
condition in Rastall gravity. As an implication, Tolman-Oppeneimer-Volkoff
equation in anisotropic form is modified. The equations in anisotropic Hartle-
Thorne formalism are also modified due to the appearance of $\lambda$, i.e. a
constant that quantifies the deviation of Rastall gravity from general
relativity. However, our results do not cover $\kappa\lambda=1/4$ and
$\lambda=1/6$, since both would mathematically give singularity in the
calculation.

††journal:

##  1 Introduction

Many years ago, Hartle and Thorne Har1 ; Har2  have proposed an approximation
procedure to calculate mass correction and deformation for slowly rotating
relativistic stars. Rotation is an aspect that is really crucial to be
considered when one studies a relativistic star. In reality, celestial objects
are normally rotating due to the law of angular momentum conservation Siv .
According to Hartle-Thorne (HT) approach, rotation is treated as a
perturbation in non-rotating configuration of the stars. As a consequence, the
metric in HT formalism undergoes the multipole expansion, which is truncated
at quadrupole order Har1 ; Har2 . Over time, this formalism has been involved
in many studies Urb ; Gup ; Pin ; Bez , mostly concerning neutron stars Urb ;
Wen ; Ber ; Sil ; Pat ; Ara2 .

Pattersons and Sulaksono Pat  derived HT formalism for anisotropic neutron
stars, which we will call anisotropic HT (aHT) formalism. Anisotropic neutron
stars have anisotropic pressure, i.e. unequal value between radial and
tangential pressures in neutron star’s matter. Some authors included
anisotropic pressure in their works, see Refs. Sil ; Pat ; Sul ; Rah ; Est1 ;
Lop ; Sha2 ; Das2 ; Mor ; Mau3 ; Dev .

Both HT and aHT formalisms are based on general relativity. In general
relativity, there exists a convention that covariant divergence of energy-
momentum tensor vanishes, i.e.

$\nabla_{\mu}T^{\mu\nu}=0.$ (1)

Rastall Ras  questioned this convention, and proposed an idea that the
covariant divergence of energy-momentum tensor is proportional to covariant
divergent of Ricci scalar. Apparently, Rastall’s idea leads to violation of
the usual conservation law. This can be phenomenologically viewed as quantum
effects in the classical context Oli . Beside that, Rastall gravity has been
involved in many works on various objects, including black holes Had ; Sak ;
Spa ; Hey ; Nas ; Zou ; Kum , wormholes Bhar ; Mal ; Hal ; Mor2 ; Hey2 ; Bro ,
strange quark stars Aya , and gravastar Maj . It has also applied in the study
of accelerated expanding universe Cap . In the context of calculation, Rastall
gravity is relatively elegant, since the Rastall’s term in the field equation
could be part of matter, which is shown in sub-section 2.1.

According to Ref. Oli , neutron stars are promising laboratories for testing
Rastall gravity, since the essence of the Rastall gravity is associated to
high curvature environment. After all, we argue that it is important to extend
aHT formalism in Rastall gravity; considering that many works on HT formalism
are applied in neutron stars calculation, and aHT formalism in Ref. Pat  is
also applied to obtain mass correction and deformation of neutron stars.

In this study, we focus on the derivation of aHT formalism in Rastall gravity,
so we do not include numerical calculation for the equations we have obtained.
We leave it for our subsequent work. However, the results we obtained in this
work are important aspects in the study of relativistic stars, considering the
significant role of Rastall gravity as a modified theory of gravity. Moreover,
we found some mistakes in two of the main equations in Ref. Pat , and we fix
them in this work. It means that this work also covers the correct ordinary
aHT formalism, since these results will be back to ordinary aHT formalism once
we turn off all Rastall’s terms.

##  2 Formalism

We separate the formalism into two parts. In sub-section 2.1, we briefly
review the Rastall gravity, while our whole derivation of aHT formalism in
Rastall gravity is given in sub-section 2.2. In this paper we use geometrized
units $G=c=1$.

###  2.1 Brief review of Rastall gravity

In Rastall gravity, the covariant divergence of energy-momentum tensor can be
expressed as Sak

$\nabla_{\mu}T^{\mu\nu}=\lambda\nabla^{\nu}R,$ (2)

where $\lambda$ is a parameter that quantifies the deviation from general
relativity, and $R$ is Ricci scalar. The field equation is now can be written
as

$G^{\mu}_{\nu}+\kappa\lambda R\delta^{\mu}_{\nu}=\kappa T^{\mu}_{\nu}.$ (3)

From Eq. (3), one can obtain

$R=\frac{\kappa T}{4\kappa\lambda-1},$ (4)

where $T$ denotes the energy-momentum scalar.

From Eq. (4), $\kappa\lambda=1/4$ would mathematically give us singularity. It
leads to an irregular Ricci scalar, and is explained as unimodular
gravitational theory Sak , so the results in our derivation do not cover it.

For the sake of simplicity, we use the following expression of the field
equation

$G^{\mu}_{\nu}=\kappa\tilde{T}^{\mu}_{\nu},$ (5)

where

$\tilde{T}^{\mu}_{\nu}=T^{\mu}_{\nu}-\delta^{\mu}_{\nu}\left(\frac{\kappa\lambda}{4\kappa\lambda-1}T\right).$
(6)

Usually in general relativity, $\kappa=8\pi$. Authors in Refs. Oli ; Had ; Sak
; Aka  also used $\kappa=8\pi$ in their calculations in the context of Rastall
gravity. Even so, we have to keep in mind that the value is obtained by
imposing general relativity on Newtonian limit. Since now we are working with
Rastall gravity, we have to use the Rastall gravity version of $\kappa$. By
applying Rastall gravity in the weak-field limit, one can get the relation Ras
; Ma

$\frac{\kappa\left(6\kappa\lambda-1\right)}{4\kappa\lambda-1}=8\pi.$ (7)

From Eq. (7), we can clearly see that if we set $\lambda=0$, then
$\kappa=8\pi$, and our calculation is recovered to general relativity. In
principle, $\kappa\lambda=1/6$ is not allowed Ma .

###  2.2 Anisotropic Hartle-Thorne formalism in Rastall gravity

Due to Eq. (6), $\tilde{T}^{\mu}_{\nu}$ for anisotropic pressure condition can
be written as

$\tilde{T}^{\mu}_{\nu}=\begin{pmatrix}-\tilde{\varepsilon}&&0&&0&&0\\\
0&&\tilde{p}&&0&&0\\\ 0&&0&&\tilde{q}&&0\\\ 0&&0&&0&&\tilde{q}\end{pmatrix},$
(8)

where
$\tilde{\varepsilon}=\varepsilon+\left(\frac{\kappa\lambda}{4\kappa\lambda-1}\right)T$,
$\tilde{p}=p-\left(\frac{\kappa\lambda}{4\kappa\lambda-1}\right)T$,
$\tilde{q}=q-\left(\frac{\kappa\lambda}{4\kappa\lambda-1}\right)T$, and
$\tilde{p}-\tilde{q}\equiv\sigma$.

Due to Rastall gravity, TOV equation in anisotropic form is modifed. It can be
written as

$\frac{dp}{dr}=\frac{-\frac{\left(\varepsilon+p\right)\left(m+\frac{\kappa}{2}r^{3}\tilde{p}\right)}{r(r-2m)}-\frac{2\sigma}{r}+\frac{2\kappa\lambda}{4\kappa\lambda-1}\frac{d\sigma}{dr}}{1-\frac{\kappa\lambda}{4\kappa\lambda-1}\left(3-\frac{d\varepsilon}{dp}\right)},$
(9)

and

$\frac{dm}{dr}=\frac{\kappa}{2}r^{2}\tilde{\varepsilon}.$ (10)

Here $m$ is the mass of the star.

Now recall HT metric, i.e. Wen ; Pat

$\displaystyle ds^{2}=$
$\displaystyle-e^{2\varphi}[1+2(h_{0}+h_{2}P_{2}(\cos\theta))]dt^{2}+\left[1+\frac{2}{r}\left(m_{0}+m_{2}P_{2}(\cos\theta)\right)\left(1-\frac{2m}{r}\right)^{-1}\right]dr^{2}$
$\displaystyle+r^{2}\sin^{2}\theta[1+2(v_{2}-h_{2})P_{2}(\cos\theta)](d\phi-\omega
dt)^{2}+r^{2}[1+2(v_{2}-h_{2})P_{2}(\cos\theta)]d\theta^{2}.$ (11)

Here $\omega$ is the angular velocity of local inertial frame and it satisfies
$\omega=\Omega-\bar{\omega}$, where $\Omega$ is the angular velocity of the
star, and $\bar{\omega}$ is the angular velocity of the star relative to the
local inertial frame; $h_{0}$, $h_{2}$, $m_{0}$, $m_{2}$, and $v_{2}$ are
perturbation functions; $P_{2}(\cos\theta)$ is the second order of Legendre
polynomial; $e^{2\varphi}$ is a function that is constrained by

$\frac{d\varphi}{dr}=\frac{m+\frac{\kappa}{2}r^{3}\tilde{p}}{r(r-2m)},$ (12)

and $m$ is constrained by Eq. (10).

The angular velocity relative to the local inertial frame can be obtained from
the following differential equation,

$\frac{1}{r^{4}}\frac{d}{dr}\left(r^{4}j\frac{d\bar{\omega}}{dr}\right)+\frac{4}{r}\frac{dj}{dr}\left(1-\frac{\sigma}{\varepsilon+p}\right)=0,$
(13)

where

$\displaystyle j$ $\displaystyle=$ $\displaystyle
e^{-\varphi}\left(1-\frac{2m}{r}\right)^{1/2},$ (14)

and boundary conditions at $r=0$ are $\bar{\omega}=\bar{\omega}_{c}$ and
$\frac{d\bar{\omega}}{dr}=0$; $\bar{\omega}_{c}$ can be chosen arbitrarily.

By following HT procedure Har1  to calculate the corrections of Einstein
tensor, we obtain Rastall gravity version of the corrections of Einstein
tensor, i.e.

$\Delta G^{\mu}_{\nu}=\delta
G^{\mu}_{\nu}-\kappa\xi\frac{d}{dr}\tilde{T}^{\mu}_{\nu},$ (15)

where $\delta G^{\mu}_{\nu}$ denotes the expansion terms of $G^{\mu}_{\nu}$.
We limit the expansion terms of $\xi$ to $l=2$,

$\xi=\xi_{0}(r)+\xi_{2}(r)P_{2}(\cos\theta).$ (16)

The relations between $\xi$ and pressure perturbation factor $\mathscr{P}$ in
aHT formalism, which do not undergo modification under Rastall gravity, are
denoted by Pat

$\displaystyle\xi_{0}=-\mathscr{P}_{0}\left(\varepsilon+p\right)\left(\frac{dp}{dr}-\frac{d\sigma}{dr}\right)^{-1},$
(17)
$\displaystyle\xi_{2}=-\mathscr{P}_{2}\left(\varepsilon+p\right)\left(\frac{dp}{dr}-\frac{d\sigma}{dr}\right)^{-1}.$
(18)

According to aHT procedure Pat , the corrections of energy-momentum tensor are
given by

$\displaystyle(\Delta
T^{t}_{t})_{l=0}=-\frac{2}{3}\left(\varepsilon+p\right)e^{-2\varphi}r^{2}\bar{\omega}\left(1-\frac{\sigma}{\varepsilon+p}\right),$
(19) $\displaystyle(\Delta T^{r}_{r})_{l=0}=(\Delta
T^{\theta}_{\theta})_{l=0}=(\Delta T^{\phi}_{\phi})_{l=0}=0.$ (20)

For aHT formalism in Rastall gravity, we need the corrections of energy-
momentum scalar, so that we can calculate the corrections of
$\tilde{T}^{\mu}_{\nu}$. The corrections of energy-momentum scalar are denoted
by

$\displaystyle\left(\Delta T\right)_{l=0}$ $\displaystyle=$
$\displaystyle-\frac{2}{3}\left(\varepsilon+p\right)e^{-2\varphi}r^{2}\bar{\omega}\left(1-\frac{\sigma}{\varepsilon+p}\right),$
(21) $\displaystyle\left(\Delta T\right)_{l=2}$ $\displaystyle=$
$\displaystyle-2\left(\varepsilon+p\right)\left(1-\frac{\sigma}{\varepsilon+p}\right)\frac{d\varepsilon}{dp}$
(22)
$\displaystyle\times\left(1-\frac{d\sigma}{dp}\right)^{-1}\mathscr{P}_{2}.$

The next stage is imposing Rastall field equation on the tensors’ corrections,

$\Delta G^{\mu}_{\nu}=\kappa\Delta\tilde{T}^{\mu}_{\nu}=\kappa\left[\Delta
T^{\mu}_{\nu}-\delta^{\mu}_{\nu}\left(\frac{\kappa\lambda}{4\kappa\lambda-1}\Delta
T\right)\right].$ (23)

The ($tt$)-component of the $l=0$ order gives

$\displaystyle\frac{dm_{0}}{dr}$ $\displaystyle=$
$\displaystyle\frac{\kappa}{2}r^{2}\frac{d\tilde{\varepsilon}}{dp}\left(1-\frac{d\sigma}{dp}\right)^{-1}\left(\varepsilon+p\right)$
(24)
$\displaystyle\times\left(1-\frac{\sigma}{\varepsilon+p}\right)\mathscr{P}_{0}+\frac{1}{12}j^{2}r^{4}\left(\frac{d\bar{\omega}}{dr}\right)^{2}$
$\displaystyle-\frac{1}{3}r^{3}\left(\frac{dj^{2}}{dr}\right)\bar{\omega}^{2}\left(1-\frac{\sigma}{\varepsilon+p}\right)$
$\displaystyle\times\left(1-\frac{\kappa\lambda}{4\kappa\lambda-1}\right).$

From the ($rr$)-component of the $l=0$ order we obtain

$\displaystyle\frac{d\mathscr{P}_{0}}{dr}$ $\displaystyle=$
$\displaystyle-\frac{m_{0}\left(1+\kappa
r^{2}\tilde{p}\right)}{\left(r-2m\right)^{2}}-\frac{\kappa
r^{2}\left(\varepsilon+p\right)}{2\left(r-2m\right)}$ (25)
$\displaystyle\times\left(1-\frac{\sigma}{\varepsilon+p}\right)\left(1-\frac{d\sigma}{dp}\right)^{-1}\beta\mathscr{P}_{0}$
$\displaystyle+\frac{1}{12}\frac{r^{4}j^{2}}{r-2m}\left(\frac{d\bar{\omega}}{dr}\right)^{2}$
$\displaystyle+\frac{1}{3}\left[\frac{d}{dr}\left(\frac{r^{3}j^{2}\bar{\omega}^{2}}{r-2m}\right)-\gamma\right],$

where

$\displaystyle\beta$ $\displaystyle=$ $\displaystyle
1+\frac{\kappa\lambda}{4\kappa\lambda-1}\left(1-\frac{d\varepsilon}{dp}\right)+\frac{2\kappa\lambda}{4\kappa\lambda-1},$
(26) $\displaystyle\gamma$ $\displaystyle=$
$\displaystyle\frac{r^{3}}{r-2m}\frac{\kappa\lambda}{4\kappa\lambda-1}\left(\frac{dj^{2}}{dr}\right)\bar{\omega}^{2}\left(1-\frac{\sigma}{\varepsilon+p}\right).$
(27)

From the ($rr$)-component of the $l=2$ order, we obtain

$\displaystyle\frac{dh_{2}}{dr}=$
$\displaystyle\left\\{-2\frac{d\varphi}{dr}+\frac{r}{r-2m}\left(2\frac{d\varphi}{dr}\right)^{-1}\kappa(\varepsilon+p)\left(1-\frac{\sigma}{\varepsilon+p}\right)\left[\left(1-\frac{d\sigma}{dp}\right)^{-1}\left(1+\frac{2\kappa^{2}\lambda}{4\kappa-1}\frac{d\varepsilon}{dp}\right)\right]-\frac{4m}{r^{3}}\right\\}h_{2}$
$\displaystyle-\frac{4v_{2}}{r(r-2m)}\left(2\frac{d\varphi}{dr}\right)^{-1}+\frac{1}{6}r^{3}j^{2}\frac{d\bar{\omega}}{dr}^{2}\left[\frac{d\varphi}{dr}r-\frac{1}{r-2m}\left(2\frac{d\varphi}{dr}\right)^{-1}\right]$
$\displaystyle-\frac{1}{3}r^{2}\frac{dj^{2}}{dr}\bar{\omega}^{2}\left[\left(1-\frac{6\sigma}{\varepsilon+p}\right)\frac{d\varphi}{dr}r+\frac{1}{r-2m}\left(1-\frac{\sigma}{\varepsilon+p}\right)\left(2\frac{d\varphi}{dr}\right)^{-1}\left(1-\frac{d\sigma}{dp}\right)^{-1}\right]$
$\displaystyle-\frac{2}{3}\frac{\kappa\lambda}{4\kappa\lambda-1}\left(\frac{dj^{2}}{dr}\right)r\bar{\omega}^{2}\left(\frac{r}{r-2m}\right)\left(2\frac{d\varphi}{dr}\right)^{-1}\frac{d\varepsilon}{dp}\left(1-\frac{\sigma}{\varepsilon+p}\right)\left(1-\frac{d\sigma}{dp}\right)^{-1},$
(28)

The equation $R_{r\theta}=\kappa\tilde{T}_{r\theta}$ gives us

$\displaystyle\frac{dv_{2}}{dr}=-2\frac{d\varphi}{dr}h_{2}+\left(\frac{1}{r}+\frac{d\varphi}{dr}\right)r^{3}\frac{dj^{2}}{dr}\bar{\omega}^{2}\left(\frac{2\sigma}{\varepsilon+p}-\frac{1}{3}\right)+\frac{1}{6}j^{2}r^{4}\left(\frac{d\bar{\omega}}{dr}\right)^{2}\left(\frac{1}{r}+\frac{d\varphi}{dr}\right).$
(29)

For the differential equations consisting $m_{0}$, $\mathscr{P}_{0}$, $h_{2}$,
and $v_{2}$; the boundary conditions are $m_{0}(r=0)=0$,
$\mathscr{P}_{0}(r=0)=0$, $h_{2}(r=0)=0$, and $v_{2}(r=0)=0$.

However, we found some mistypes in Ref. Pat :

  1. 1.

the differentials in Eq. (40)–(43) should be total differentials instead of
partial differentials, see Eq. (74)–(77) in Ref. Har1  for comparison;

  2. 2.

the second term of the right-handed side in Eq. (62) should include the factor
$\left(1-\frac{d\sigma}{dp}\right)^{-1}$;

  3. 3.

the factor $\left(1-\frac{d\sigma}{dp}\right)$ in the last term of the right-
handed side of Eq. (71) should be $\left(1-\frac{d\sigma}{dp}\right)^{-1}$.

By setting $\lambda=0$ and applying the corrections of mistypes in Ref. Pat ,
we can see that our formalism is reduced to aHT formalism, which is based on
general relativity. We can easily find that in ordinary aHT formalism,
$\beta=1$ and $\gamma=0$. The expression of Eq. (29) is same with the one in
Ref. Pat , but actually it is different, since $j$ in the equation contains
$\varphi$, which is modified due to the appearance of $\lambda$. Moreover, by
setting $\lambda=0$ and $\sigma=0$, the ordinary HT formalism is recovered.

Mass correction of the stars is given by

$\delta M=m_{0}(r=R)+\frac{J^{2}}{R^{3}},$ (30)

where $J$ is angular momentum of the stars, and $R$ is radius of the star..

To calculate deformation of the stars, we have to find $\mathscr{P}_{2}$ first
by using the formula Har1

$\mathscr{P}_{2}=-\left(h_{2}+\frac{1}{3}e^{-2\varphi
r^{2}\bar{\omega}^{2}}\right).$ (31)

Therefore, we will have the values of $\mathscr{P}_{0}$ and $\mathscr{P}_{2}$
in our hands, so we can find $\xi_{0}$ and $\xi_{2}$ (see Eq. (17) and Eq.
(18)). Finally, the radius correction $\delta r$ of the stars can be
determined,

$\delta r=\xi_{0}(r=R)+\xi_{2}(r=R)P_{2}(\cos\theta).$ (32)

At the pole, $\theta=0$; while at the equator, $\theta=\pi/2$. According to
the mathematical expression of $\delta r$, we can see that the oblateness of
the stars is not related to the first term of Eq. (32), and it is only related
to its second term.

##  3 Concluding remarks

In this letter, we have extended the aHT formalism in Rastall gravity. As the
result, the equations in aHT formalism are modified due to the appearance of
$\lambda$, a constant that quantifies the deviation of Rastall gravity from
general relativity. The results can accommodate the calculation of mass
correction and deformation of slowly rotating anisotropic relativistic stars
in Rastall gravity framework. However, it is important to note that our
results do not cover $\kappa\lambda=1/4$ and $\kappa\lambda=1/6$. In this
work, we also have fixed the mistypes in some equations in Ref. Pat .

###### Acknowledgements.

MLP acknowledges Indonesia Endowment Fund for Education Agency (LPDP) for
financial support. FPZ and GH would like to thank Kemendikbudristek, Ministry
of Education, Culture, Research, and Technology, Republic of Indonesia for
partial financial support. HLP would like to thank Ganesha Talent
Assistantship (GTA) Institut Teknologi Bandung for financial support. MLP and
HLP would like to thank the members of the Theoretical Physics Groups of
Institut Teknologi Bandung for their hospitality. We also acknowledge Muhammad
Fitrah Alfian Rangga Sakti for precious early discussions.

## References

  * (1) J.B. Hartle, Astrophys. J. 150, 1005 (1967)
  * (2) J.B. Hartle, K.S. Thorne, Astrophys. J. 153, 807 (1968)
  * (3) C. Sivaram, K. Arun, Open Astron. J. 5, 7-11 (2012)
  * (4) M. Urbanec, J.C. Miller, Z. Stuchlík, Mon. Not. R. Astron. Soc. 433, 1903–1909 (2013)
  * (5) A. Gupta, S. Iyer, A.R. Prasanna, Class. Quantum Grav. 13, 2675 (1996)
  * (6) E.O. Pinho, C.C. Barros Jr., Eur. Phys. J. C 83, 745 (2023)
  * (7) B. Bezděková, J. Bičák, Phys. Rev. D 108, 084043 (2023)
  * (8) D.H. Wen, W. Chen, L.G. Liu, Chin. Phys. Lett. 22, 1604 (2005)
  * (9) E. Berti, F. White, A. Maniopoulou, M. Bruni, Mon. Not. R. Astron. Soc. 358, 923-938 (2005)
  * (10) H.O. Silva, C.F.B. Macedo, E. Berti, L.C.B. Crispino, Class. Quantum Grav. 32, 145008 (2015)
  * (11) M.L. Pattersons, A. Sulaksono, Eur. Phys. J. C 81, 698 (2021)
  * (12) E. Aranguren, J.A. Font, N. Sanchis-Gual, R. Vera, Phys. Rev. D 108, 104065 (2023)
  * (13) A. Sulaksono, Int. J. Mod. Phys. E 24, 1550007 (2015)
  * (14) A. Rahmansyah, A. Sulaksono, A.B. Wahidin, A.M. Setiawan, Eur. Phys. J. C 80, 769 (2020)
  * (15) J. Estevez-Delgado, A. Cleary-Balderas, A.T. Marguía, R. Soto-Espitia, Eur. Phys. J. Plus 134, 600 (2019)
  * (16) I. Lopes, G. Panotopoulos, A. Rincón, Eur. Phys. J. Plus 134, 454 (2019)
  * (17) M. Sharif, A. Majid, Eur. Phys. J. Plus 135, 558 (2020)
  * (18) H.C. Das, J.A. Pattnaik, S.K. Patra, Phys. Rev. D 107, 083007 (2023)
  * (19) P.H.R.S. Moraes, G. Panotopoulos, I. Lopes, Phys. Rev. D 103, 084023 (2021)
  * (20) S.K. Maurya, B. Mishra, S. Ray, R. Nag, Chin. Phys. C 46, 105105 (2022)
  * (21) K. Dev, M. Gleiser, Gen. Relativ. Gravit. 34, 1793–1818 (2002)
  * (22) P. Rastall, Phys. Rev. D 6, 3357 (1972)
  * (23) A.M. Oliveira, H.E.S. Velten, J.C. Fabris, L. Casarini, Phys. Rev. D 92, 044020 (2015)
  * (24) H.L. Prihadi, M.F.A.R. Sakti, G. Hikmawan, F.P. Zen, Int. J. Mod. Phys. D 29, 3 (2020)
  * (25) M.F.A.R. Sakti, A. Suroso, A. Sulaksono, F.P. Zen, Phys. Dark Univ. 35, 100974 (2022)
  * (26) E. Spallucci, A. Smailagic, Int. J. Mod. Phys. D 27, 02 (2018)
  * (27) Y. Heydarzade, H. Moradpour, F. Darabi, Can. J. Phys. 95, 12 (2017)
  * (28) G.G.L. Nashed, Universe 8, 10 (2022)
  * (29) D.C. Zou, M. Zhang, C. Wu, R.H. Yue, Adv. High Energy Phys. 2020, 4065254 (2020)
  * (30) R. Kumar, S.G. Ghosh, Eur. Phys. J. C 78, 750 (2018)
  * (31) P. Bhar, Chin. J. Phys. 87, 782 (2024)
  * (32) A. Malik, A. Ashraf, F. Mofarreh, A. Ali, M. Shoaib, Int. J. Geom. Methods Mod. Phys. 20, 09 (2023)
  * (33) S. Halder, S. Bhattacharya, S. Chakraborty, Mod. Phys Lett. A 34, 12 (2019)
  * (34) H. Moradpour, N. Sadeghnezhad, S.H. Hendi, Can. J. Phys. 95, 12 (2017)
  * (35) Y. Heydarzade, M. Ranjbar, Eur. Phys. J. Plus 138, 703 (2023)
  * (36) K.A. Bronnikov, V.A.G. Barcellos, L.P. de Carvalho, J.C. Fabris, Eur. Phys. J. C 81, 395 (2021)
  * (37) A. Banerjee, T. Tangphati, S. Hansraj, A. Pradhan, Ann. Phys. 451, 169267 (2023)
  * (38) K. Majeed, G. Abbas, J. Phys. Commun. 6, 045005 (2022)
  * (39) M. Capone, V.F. Cardone, M.L. Ruggiero, J. Phys. Conf. Ser. 222, 012012 (2010)
  * (40) Ö. Akarsu, N. Katırcı, S. Kumar, R.C. Nunes, B. Öztürk, S. Sharma, Eur. Phys. J. C. 80, 1050 (2020)
  * (41) M.S. Ma, R. Zhao, Eur. Phys. J. C. 77, 629 (2017)