# Effect of hydrostatic strain on the mechanical properties and topological
phase transition of bi-alkali pnictogen NaLi2Bi

Seyed Mohammad bagher Malek Hosseini  Department of Physics, Vali-e-Asr
University of Rafsanjan, Rafsanjan, Iran  Shahram Yalameha  Faculty of
Physics, University of Isfahan, 81746-73441, Isfahan, Iran.

###### Abstract

The bi-alkali pnictogens have attracted significant attention for
optoelectronic and photocathodic device applications. However, in most of the
compounds belonging to this family, there has been less effort put into
investigating the mechanical properties and topological phase transitions
(TPT) of the compounds. Here, in the framework of density functional theory,
the mechanical properties and topological phase transition of NaLi2Bi under
hydrostatic pressures are investigated. Elastic constants and phonon
calculations have shown the mechanical and dynamical stability of this
compound under hydrostatic tension and compression. The analysis of the
elastic constants show that the NaLi2Bi in the equilibrium state is an auxetic
material with a negative Poisson’s ratio of -0.285, which changes to a
material with a positive Poisson’s ratio under hydrostatic tension. Meanwhile,
Poisson’s ratio and Pugh ratio indicate that this compound has brittle
behavior and maintains it under hydrostatic pressures. The calculated results
of the band structure within the generalized gradient approximation (GGA)
(Tran-Blaha modified Becke-Johnson exchange potential approximation (TB-mBJ))
show that NaLi2Bi is a nontrivial topological material (trivial topological
material). It was found that hydrostatic compression (tension) in the GGA (TB-
mBJ) approach leads to a transition from a nontrivial (trivial) to a trivial
(nontrivial) topological phase for this compound. Moreover, the calculated
Wannier charge centers confirm the TPT. Identifying the mechanisms controlling
the auxetic behavior and TPT of this compound offers a valuable feature for
designing and developing high-performance nanoscale electromechanical and
spintronic devices.

Mechanical properties; Bi-alkali pnictogen; Topological phase transition;
auxetic material; Hydrostatic pressure

##  I I. INTRODUCTION

In general, bi-alkali pnictogens have the formula XY2Z, in which X and Y
represent alkali metals and Z represents group 7 elements. These compounds
comprise an interesting class of semiconducting compounds that have
technological importance in optoelectronics, sensing, and photocathodic
devices, due to their favorable electronic properties [1]. It is expected that
these compounds will provide high-intensity X-ray electron Lasers in the
future [2, 3, 4]. In addition, these compounds are used in electron radiative
devices because of their feasible work function [5]. Recently, some of these
compounds have attracted the attention of researchers due to their narrow band
gap and strong spin-orbit coupling (SOC) [6, 7].

From the experimental point of view, bi-alkali pnictogens have also been
investigated. Using the in-situ X-ray diffraction method, Ruiz et al., [8]
investigated the crystalline structure of multi-alkali photocathodes. Also,
Cultrera et al., [9] have synthesized and characterized top-ordered bi and
multi-alkali pnictogen photocathodes and involved the Rb element in their
chemical structure. Feng et al., [10] studied the effect of surface roughness
on the emission of alkali penictogen cathodes in electron guns and developed
K-Cs-Sb cathodes using co-deposition. Ding et al. [11] showed that co-
evaporation is the best growth method for KCs2Sb photocathodes because it
produces films with a smooth surface, high quantum efficiency, and proper
crystal structures. Schubert et al., [12] investigated the bi-alkali pnictogen
photocathodes for high-brightness accelerators.

From the theoretical perspective, the major attention on the bi-alkali
pnictogens is to investigate their electronic and structural properties.
However, there are also reports on the optical, elastic, and topological
properties of a number of bi-alkali pnictogen compounds. Khan et al., [13]
reported the structural and optoelectronic properties of bi-alkali pnictogen
compounds XY2Z (X and Y = Li, Na, K, Rb, Cs; Z = N, P, As, Sb, Bi). Amador et
al., [14] described the electronic structure and optical properties of cubic
Na2KSb and hexagonal NaK2Sb compounds. Murtaza et al., [15] investigated the
structural, elastic, electronic, and optical properties of some bi-alkali
antimonides. Cocchi et al., [16] studied the electronic and optical properties
of the CsK2Sb compound. Kalarasse et al., [17] reported the electronic and
elastic properties of the alkali pnictide compounds Li3Sb, Li3Bi, Li2NaSb, and
Li2NaBi in the equilibrium state. Jin et al., [18] investigated the electronic
structure, doping effect, and topological signature in hexagonal compounds
Li3-xNaxM (x = 3, 2, 1, 0; M = N, P, As, Sb, Bi). Sklyadneva et al., [7]
reported that the KNa2Bi compound is nontrivial topological material. Yalameha
et al., [19] investigated that the crystal structure of KNa2Sb was also a
topological insulator (TI) under hydrostatic pressure. Furthermore, they have
recently shown that Rb(Na, K)2Bi and Cs(Na, K)2Bi compounds can take a diverse
set of topological states by strain-engineering [20, 21]. They also
investigated the mechanical and anisotropic elastic properties of Cs(Na, K)2Bi
compounds under hydrostatic pressure [22]. Very recently, Song et al. [23]
comprehensively investigated the thermal transport and thermoelectric
properties of the bi-alkali bismuthide compound NaLi2Bi.

So far, no comprehensive study has been devoted to examining the mechanical,
anisotropic elastic, and topological properties of the cubic bi-alkali
pnictogen compound NaLi2Bi. In this work, the effect of hydrostatic tension
and compression on mechanical properties (e.g., Young’s modulus, Poisson’s
ratio, Pugh ratio, Cauchy pressure, etc.) and the TPT of this compound have
been investigated. The mechanical and dynamic stability of the compound under
the studied pressures has also been studied. It should be noted that in this
composition, the TPT point in PBE-GGA (from nontrivial to trivial state) and
TB-mBJ (from trivial to nontrivial state) approaches and their critical
pressures have been observed and investigated. Even though this compound does
not have any experimental reports, it is hoped that, based on these results,
experimental research on this compound can be conducted.

##  II II. Computational details

In this paper, the density functional theory (DFT)-based Wien2k package [24]
has been used to explore various physical properties of the titled compound.
This package implements the full potential linearized augmented plane wave
method (FP-LAPW) [25]. In this work, the generalized gradient approximation
(GGA) of the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional is
used [26]. The Tran-Blaha modified Becke-Johnson exchange potential
approximation (TB-mBJ) [27] calculations were also conducted in order to
obtain improved band gaps. The important parameters Gmax (i.e., the potential
and charge density Fourier expansion parameter) and RMT×Kmax, where RMT is the
smallest muffin-tin (MT) sphere radii and Kmax is the largest reciprocal
lattice vector used in the plane wave expansion, are set to 15.0 Bohr-1 and
9.5 (Ry)1/2, respectively. In the calculations of this compound, a 13×13×13
k-mesh of points is used to achieve self-consistency. The phonon dispersion is
calculated using the Phonopy package [28]. In Phonopy, first-principles phonon
calculations are performed using the finite displacement method [29]. So, we
used a supercell model and the finite displacement approach with a 2×2×2
conventional unit cell supercell and atomic displacement distance of 0.01 Å.
Irelast [30] and ElATools [31] codes are used for elastic constants
calculations and their analysis, respectively. Wannier90 [32] and Z2pack [33]
packages have been used to calculate Wannier charge centers (WCC). It is
noteworthy that the hydrostatic pressures (tension- and compression-type) are
evaluated using the volume ratio V/V0, where V0 is the equilibrium volume of
the primitive unit cell. In all calculations presented in this study, spin-
orbit coupling (with full relativistic effect) has been considered unless
otherwise stated.

##  III III. Results and discussion

###  III.1 III.A Structural properties, stability conditions

Figure 1: (a) Conventional crystal structure of bi-alkali pnictide NaLi2Bi
with space group Fm-3m (No. 225). The atomic positions in a unit cell are: Na
atoms at (0, 0, 0), Bi atoms at (0.5, 0, 0) and Li(1)/Li(2) atoms at (0.75,
0.25, 0.25)/(0.25, 0.75, 0.75). Phonon dispersion curves along
K-$\Gamma$-L-W-X high-symmetry points at (b) V/V0 = 0. 930, (c) V/V0 = 0.940,
(d) V/V0 = 1.0, (e) V/V0 = 1.060 and (f) V/V0 = 1.070. The phonon dispersion
exhibits no imaginary frequency (negative frequency), confirming the dynamic
stability of these compounds. (g) The first 3D Brillouin zone and the high-
symmetry points. Note that here the meaning of the G (Gamma) point is the same
as $\Gamma$ point.

The crystal structure of the NaLi2Bi compound is cubic with space group Fm-3m
(No. 225), and the schematic conventional crystal structure of this compound
is shown in Fig. 1(a). The atomic positions in a unit cell for Li(1)/Li(2),
Na, and Bi atoms are (0.75, 0.25, 0.25)/(0.25, 0.75, 0.75), (0, 0, 0), and
(0.5, 0, 0), respectively. A list of the results of first-principles
calculations of the structural and elastic properties of this compound along
with available theoretical lattice constants, is provided in Table 1. By
fitting the Murnaghan equation of state to the total energy as a function of
volume, the equilibrium lattice parameter (a0) can be determined [34] (see
Fig. 1S of Supplementary Information). The calculated equilibrium lattice
constant (7.021 Å) is in good agreement with theoretical predictions (see
Table 1).

There are only three independent elastic constants for materials with cubic
symmetry: C11, C12, and C44. As can be seen in Table 1, the calculated elastic
constant values at the equilibrium state (V/V0 = 1.0) in this work are
consistent with other theoretical results. In order for a cubic crystal system
to be mechanically stable, it needs to satisfy the following Born-Huang
criteria [35],

$C_{44}>0,C_{11}-C_{12}>0,C_{11}+2C_{12}>0;$ (1)

According to Table 1 and Eq. 1, all elastic constants at the studied pressures
of this compound satisfy the criteria of mechanical stability. In the
equilibrium state and for hydrostatic tension (V/V0 = 1.06 and 1.07) and
compression (V/V0 = 0.930 and 0.940) values, phonon dispersion curves for
NaLi2Bi along the high-symmetry points in the Brillouin zone (see Fig. 1(g))
are plotted in Figs. 1(b-f). It reveals that the NaLi2Bi compound is
dynamically stable at the studied pressures and equilibrium state since there
are no modes with negative frequency values. To find the thermodynamical
stability of the NaLi2Bi compound, the Formation energies (Ef) of volume
ratios V/V0 = 0.93, 0.94, 1.0, 1.06, and 1.07 are calculated. Table 1S shows
the calculated values of the Ef per formula unit. In this table, the
calculated values of Ef per formula unit are negative, indicating this
compound in the studied pressures is thermodynamically stable. The thermal
stability of this compound has been recently investigated by Song et al. [23].
Thus, the NaLi2Bi compound is stable and can be synthesized in the laboratory.

Table 1: Calculated lattice constants (a0), elastic constants (Cij), Cauchy
pressure (CP), Poisson’s ratio ($\nu$), Pugh ratio (B/G), machinability index
($\mu_{M}$) and universal anisotropy index (AU) relative to volume ratios
(V/V0) for NaLi2Bi compound. It should be noted that the equivalent pressures
are specified with each of the volume ratios.

V/V0 |  | Pressure (GPa)
---
| a0 (Å)
---
| C11
---
(GPa)
| C12
---
(GPa)
| C44
---
(GPa)
| CP
---
(GPa)
$\nu$ | B/G | $\mu_{M}$ | AU | AE
0.940 | 1.7 | 6.887 | 50.48 | 19.89 | 42.35 | -22.46 | 0.1433 | 1.6527 | 1.2434 | 1.3562 | 2.20
0.955 | 1.0 | 6.922 | 47.89 | 19.14 | 41.52 | -22.37 | 0.1417 | 1.5054 | 1.0570 | 1.4815 | 2.27
0.970 | 0.8 | 6.958 | 40.65 | 18.00 | 38.07 | -20.07 | 0.1500 | 1.3864 | 0.9831 | 1.9913 | 2.58
0.985 | 0.4 | 6.993 | 39.96 | 17.11 | 36.10 | -18.98 | 0.1475 | 1.3440 | 0.8880 | 1.7721 | 2.45
1.0 | 0 |  | 7.021,
---
6.97 a,
7.016 b
| 36.03,
---
41.70 b
41.73 c
| 15.97,
---
17.30 b
17.23 c
| 33.96,
---
30.61 b
25.99 c
-17.99 | 0.1473 | 1.0853 | 0.6672 | 2.0192 | 2.60
1.015 | -0.35 | 7.063 | 31.91 | 15.26 | 23.43 | -8.17 | 0.2026 | 1.0849 | 0.6851 | 1.4036 | 2.32
1.030 | -0.65 | 7.098 | 32.21 | 14.55 | 20.79 | -6.24 | 0.2096 | 1.0883 | 0.6710 | 0.9356 | 2.01
1.045 | -1.3 | 7.133 | 29.83 | 14.47 | 18.53 | -4.06 | 0.2284 | 1.0579 | 0.6919 | 0.9937 | 2.07
1.060 | -1.5 | 7.168 | 32.05 | 15.44 | 16.87 | -1.43 | 0.2482 | 1.0683 | 0.7104 | 0.6287 | 1.81

aRef. [13], bRef. [17], cRef. [23].

###  III.2 III.B Mechanical properties

The elastic constants (Cij) are key parameters for the mechanical behaviour of
materials, including mechanical stability, anisotropic response, stiffness,
and hardness. The computed Cij and their related mechanical properties, such
as Poisson’s ratio ($\nu$), Pugh ratio (B/G), machinability index ($\mu_{M}$),
and the universal anisotropic index (AU) under hydrostatic pressures, are
reported in Table 1. It is worth mentioning that the studied hydrostatic
pressures are of two types, tension and compression, which are addressed by
the volume ratios (V/V0) in this table. However, the range of pressures
related to volume ratios is also included. The volume ratios in this study
have been chosen in a way that allows us to observe changes in mechanical
behaviour (e.g., negative Poisson’s ratio) as well as topological phase
transition (TPT).

Among the three constants C11, C12, and C44, the C44 constant represents
resistance to shear deformation, while the C11 relates to unidirectional
strain along the principal crystallographic directions. The value of C44 in
Table 1 is lower than C11 in all pressures, which means that the cubic cell is
more easily deformed by a shear in comparison to the unidirectional stress. A
material’s Cauchy pressure (CP) can also be a useful mechanical parameter, and
it is defined as CP = (C12 \- C44) [31]. Angular characteristics of atomic
bonding as well as the brittleness and ductility of a material can be
described by Cauchy pressure [31, 35]. The Cauchy pressure of a ductile
material is positive, while the Cauchy pressure of a brittle material is
negative. On the other hand, ionic and covalent bonding in a material is
correlated with positive and negative Cauchy pressures, respectively.
According to Pettifor’s rule [36], a material with a large positive CP will
possess significant metallic bonds and exhibit high levels of ductility, while
a material that has a negative CP will possess more angular bonds, resulting
in greater brittleness and covalent bonding. Hence, the negative value of CP
predicts that the NaLi2Bi compound is brittle in nature, and some covalent
bonding is present in this compound. In addition, these results demonstrate
that in tension (compression) hydrostatic pressures, this material approaches
a high level of ductility (brittleness).

Figure 2: 3D representation surfaces of Young’s modulus (E) of NaLi2Bi under
hydrostatic tension and compression. The more localized the colours and the
more complex the colour distribution pattern, the higher the level of
anisotropy. The purple arrow in V/V0 = 1.060 shows that the localization in
the color’s distribution is reduced.

Using the Voigt-Reuss-Hill (VRH) approximation, the Poisson’s ratio ($\nu$)
and Pugh ratio (B/G) are estimated. Details of this approximation are
addressed in Ref. [31]. Poisson’s ratio and Pugh ratio are two factors that
determine a system’s brittle and ductile nature in addition to Cauchy
pressure. The Poisson’s ratio measures deformation (i.e., expansion or
contraction) along the perpendicular to the loading direction. Frantsevich et
al. [37] have separated the brittleness and ductility of a material based on
Poisson’s ratio ($\nu$), defining the border between the two as $\nu$ = 0.26.
So, materials with Poisson’s ratio values greater and smaller than 0.26 are
considered ductile and brittle, respectively. Moreover, in 1954, Pugh has
proposed that the bulk modulus to shear modulus ratio (or vice versa)
represents a good criterion to identify the brittleness and ductility of
materials [38]. The critical value that separates brittle and ductile
materials is about 1.75. Therefore, materials with a Pugh’s ratio greater than
1.75 are predicted to be ductile; otherwise, they are expected to be brittle.
Based on these two criteria, it can be seen that the NaLi2Bi compound has
brittle in nature in the equilibrium state (V/V0 = 1.0), and this
characteristic is preserved in the studied hydrostatic pressures. It is
noteworthy that a significant agreement exists between the two conditions
related to $\nu$ and B/G, and the results associated with CP. Poisson’s ratio
is more sensitive to the detection of physical characteristics than two CP and
B/G parameters. The Poisson’s ratio also provides information about the
bonding nature of the materials. Materials that are covalently bonded tend to
have Poisson’s ratios around 0.1, whereas materials that are ionic tend to
have Poisson’s ratios around 0.25 [39]. Based on Poisson’s ratio values, it is
predicted that most of the chemical bonds in this composition are covalent,
and the covalent bonds become more intense at compression hydrostatic
pressures. Also, in tensile hydrostatic pressures, ionic chemical bonds play a
prominent role compared to covalent bonds. A Poisson’s ratio can provide
insight into the nature of interatomic forces in solids [40, 41]. Materials
that have Poisson’s ratios within 0.25 and 0.50 are called central force
solids, whereas non-central force solids have Poisson’s ratios either less
than 0.25 or greater than 0.50 [41]. In the NaLi2Bi compound, therefore, non-
central forces of atoms prevail, while in tension hydrostatic pressures,
central forces become more important (see Table 1).

Another useful mechanical performance index is the machinability index
($\mu_{M}$), and it is defined as $\mu_{M}$ = B/C44 (B is bulk modulus). This
parameter measures a material’s machinability, or its ability to be machined
using cutting or shaping tools. A material’s degree of machinability is
becoming increasingly important to the industry today because it determines
factors such as how much power is needed for cutting, what temperature to use,
and how much plastic strain to apply. Furthermore, it can be used to determine
a solid’s plasticity and dry lubricating ability [42, 43]. Having a high
$\mu_{M}$ value indicates excellent lubricating properties, which results in
lower friction. The value of $\mu_{M}$ (= 0.6672) indicates that the NaLi2Bi
compound has a good level of machinability. Under hydrostatic compression
(V/V0 ¡ 1), however, this characteristic increases, indicating it could be an
attractive feature for machine tool applications.

Figure 3: 3D representation surfaces of Poisson’s ratio for NaLi2Bi in (a)
V/V0 = 0.940 , (b) V/V0 = 1.0, and (c) V/V0 = 1.060 states. The maximum
positive, minimum positive, and negative values of Poisson’s ratio are
$\nu_{max}$, $\nu_{min}$, and $\nu_{negative}$, respectively. The maximum and
minimum values of Poisson’s ratio are located at [110]([101] or [011]) and
[100]([001] or [010]) directions, respectively.

Plastic deformation in crystals, the formation of microscale cracks in
ceramics, and plastic relaxation in thin films are all influenced by elastic
anisotropy. Hence, evaluating elastic anisotropy factors is essential for
predicting solid behaviour under different external stress conditions. The
universal anisotropy index (AU) and the anisotropy of Young’s modulus (AE) of
NaLi2Bi are listed in Table 1 (the calculation of this two indices is detailed
in Ref. [31]). In V/V0 = 1.0, there is the highest degree of anisotropy (AU =
2.0192). According to the trend of changes in this index under hydrostatic
tension and compression, it can be found that the degree of anisotropy (i.e.,
AU) decreases slowly in tension-type of pressures and more rapidly in
compression-type of pressures. Young’s modulus (E) can also be a good
indicator of the level of anisotropy. The 3D representation surfaces of
Young’s modulus are calculated by ElaTools and included in Fig. 2. The more
localizing colors and the more complex the pattern of color distribution in
this figure, the higher the level of anisotropy is. According to AE in Table 1
and Fig. 2, it can be seen that NaLi2Bi has high anisotropy. With increasing
tension-type pressure, anisotropy decreases, and its minimum value corresponds
to a state with V/V0 = 1.060 and AE = 1.81. This is due to the reduced
localization in the colors distribution pattern (see purple arrow).

Figure 4: 2D polar representation of Poisson’s ratio in the (100)-plane for
NaLi2Bi in equilibrium state (V/V0 =1.0 ), and hydrostatic tension (V/V0 ¿ 1
), and compression (V/V0 ¡ 1). Blue (green) and red colors represent the
maximum (minimum) positive and negative values of Poisson’s raito,
respectively. With the increase of hydrostatic tension (V/V0 ¿ 1), the value
of NPR decreases until the critical point V/V0 = 1.060, while in hydrostatic
compression (V/V0 ¡ 1), the value of NPR remains almost constant.

Materials belonging to the auxetic family possess negative Poisson’s ratio
(NPR), causing expansion under tension rather than contraction. So, this
enables the development of innovative technologies for impact-resistant
composites, ultra-precise sensors, tougher ceramics, and high-performance
armor due to auxetic materials’ atypical elastic behaviour [44, 45]. This
feature is investigated using 3D surface representations (see Fig. 3) and 2D
polar plots (see Fig. 4) of Poisson’s ratio of NaLi2Bi in (100)-plane under
hydrostatic tension and compression. Based on Fig. 3, the NaLi2Bi compound has
a NPR of $\nu_{negative}$ = -0.285 along the [110]([101] or [011])-direction
in the equilibrium state. As a result, this compound can be regarded as an
auxetic material. The results show that with increasing hydrostatic tension,
the negative Poisson’s rate decreases and reaches a critical value of -0.017
(see Fig. 3(a)). From this point on, the compound loses its auxetic behaviour
(see Fig. 2S). Therefore, a type of mechanical phase transition (MPT) from the
auxetic state to the non-auxetic state has occurred for this compound. Fig. 4
displays the 2D polar representation of Poisson’s ratio in the (100)-plane for
the studied pressures. The blue (green) represents the maximum (minimum) value
and the red represents the negative value. In this figure, it can be clearly
seen that with the increase of hydrostatic tension (V/V0 ¿ 1), the value of
NPR decreases to the critical point of V/V0 = 1.060. In hydrostatic
compression (V/V0 ¡ 1), the value of NPR remains almost constant. It is
noteworthy that the maximum and minimum values of Poisson’s ratio are located
at [110]( [101] or [011]) and [100]([001] or [010]) directions, respectively.

###  III.3 III.C Topological properties

Figure 5:  (a) The band structures of the NaLi2Bi compound under hydrostatic
tension (V/V0 ¿ 1 ), compression (V/V0 ¡ 1) , and equilibrium states (V/V0 =
1.0) using the PBE-GGA approach. The $\Gamma_{6}$ (Bi-s-liked) and
$\Gamma_{8}$ (Bi-p-liked) states are marked with orange and green arrows,
respectively. The black arrow shows the critical point of the topological
phase transition in the V/V0 = 0.955 state. The orbital projected band
structures of the NaLi2Bi in the (b) V/V0 = 0.940, (c) equilibrium state (V/V0
= 1.0), and (d) V/V0 = 1.060 along high symmetry line L- $\Gamma$-W. The
evolution lines of the Wannier charge centers along $k_{y}$ in the $k_{x}$ = 0
and $k_{x}$ = 0.5 planes in the (e) V/V0 = 0.940, (f) V/V0 = 1.0, and (g) V/V0
= 1.060 states. The arbitrary reference line (green line) is shown for $k_{x}$
= 0 and $k_{x}$ = 0.5 planes. Figure 6:  (a) The band structures of the
NaLi2Bi compound under hydrostatic tension (V/V0 ¿ 1 ), compression (V/V0 ¡ 1
) , and equilibrium states (V/V0 = 1.0) using the TB-mBJ approach. The
$\Gamma_{6}$ (Bi-s-liked ) and $\Gamma_{8}$ (Bi-p-liked) states are marked
with orange and green arrows, respectively. The black arrow shows the critical
point of the topological phase transition in the V/V0 = 1.06 state. The
orbital projected band structures of the NaLi2Bi in the (b) V/V0 = 0.930, (c)
equilibrium state (V/V0 = 1.0), and (d) V/V0 = 1.070 along high symmetry line
L- $\Gamma$-W. The evolution lines of the Wannier charge centers along $k_{y}$
in the $k_{x}$ = 0 and $k_{x}$= 0.5 planes in the (e) V/V0 = 0.930, (f) V/V0 =
1.0, and (g) V/V0= 1.070 states. The arbitrary reference line (green line) is
shown for $k_{x}$= 0 and $k_{x}$ = 0.5 planes.

Topological insulators (TIs) with novel surface states have recently drawn
considerable theoretical and experimental attention [21, 46, 47]. The spin-
orbit coupling (SOC) and band inversion at the Fermi energy give TI remarkable
electronic properties. This makes TI ideal for applications in electronics,
thermoelectrics, optoelectronics, superconductors, and magnetics [47, 48, 49,
50]. It has been found that it is possible to control topological band order
using electromagnetic fields, alloying, pressure, strain, and temperature
effects [19, 20, 51, 52, 53]. Here, the mechanism of hydrostatic strain is
employed in this study in order to investigate the TPT of the NaLi2Bi
compound. Here, the mechanism of hydrostatic strain is employed in this study
in order to investigate the TPT of the NaLi2Bi compound using PBE-GGA and TB-
mBJ approaches.

• PPE-GGA approach: To investigate topological band order, we first examine
the equilibrium state of the NaLi2Bi compound. Fig. 3S illustrates the orbital
projected band structure of NaLi2Bi without SOC at an equilibrium state.
Compared to Ref. [13], the bandgap of this compound along high symmetry line
L-$\Gamma$-W is 0.41 eV, which is in good agreement with this reference (0.4
eV). In this figure, based on the orbital projected bands (fat bands), it is
given that the Bi-s-liked ($\Gamma_{1}$) is above Bi-p-liked ($\Gamma_{5}$) at
$\Gamma$ point. It can be concluded that in the absence of SOC, the band
inversion did not occur and thus the compound is a semiconductor with an
trivial bandgap. This configuration is similar to HgTe compound in the absence
of SOC [54]. When the SOC is included, the bandgap is closed and Bi-s-liked
orbitals are exchanged with Bi-p-liked (Bi-s-liked below Bi-p-liked band), and
a band inversion of sp-type occurs (see Fig. 5(a and c)). In this case (V/V0 =
1.0), the compound is reduced to a zero-gap semimetal with an inverted band
order. This is the result of familiarity, which is created due to the presence
of strong SOC in TIs and topological semimetals [7, 54, 55]. In cubic crystals
and in the presence of SOC, the $\Gamma$ point around Fermi level can be
traced back to the s-like doublet $\Gamma$6 (j = 1/2) state, the p-like
quartet $\Gamma$8 (j = 3/2) state, and the split-off doublet $\Gamma$7 state
(j = 1/2) (($\Gamma_{7}$ state is indicated in Fig. 4S)) [56]. Here ”j” is the
total angular momentum. In fact, SOC splits the $\Gamma_{5}$ state (see Fig.
3S in the absence of SOC) into $\Gamma_{8}$ (j = 3/2, denoted as p3/2) and
$\Gamma_{7}$ (j= 1/2, denoted as p1/2) states [54, 55]. Based on the energy
difference between the $\Gamma_{8}$ and $\Gamma_{6}$ states in the band
structure, a parameter can be defined (i.e., $\delta$ = [$\Gamma_{6}$
-$\Gamma$8]) that shows the band inversion behaviour, which is known as the
band inversion strength (BIS). Therefore, for nontrivial topological phases,
this value is negative ($\delta$ ¡ 0) and for trivial topological phases, this
value is positive ($\delta$ ¿ 0). Also, this value becomes zero for the TPT
point. In order to investigate TPT, the evolution of the band structure at the
$\Gamma$ point under hydrostatic tension (V/V0 ¿ 1) ,and compression (V/V0 ¡
1) has been evaluated in Fig. 5(a). In addition, in this figure, $\Gamma_{6}$
and $\Gamma_{8}$ states are marked with orange and green arrows, respectively.
It can be found that when hydrostatic pressures of the tension type are
applied, the nontrivial topological phase is preserved and the —$\delta$—
($\delta$ ¡ 0) also increases. For example, the orbital projected band
structure in the V/V0 = 1.060 state from Fig. 5(d) shows that there is the sp-
band inversion and the only difference with the equilibrium state is in the
value of the $\delta$. Meanwhile, with the increase of the hydrostatic
pressure of the compression type, the $\delta$ decreases and reaches zero
value in V/V0 = 0.955 state. This state shows the critical point of the
topological phase transition, which is indicated by the black arrow. The value
of critical pressure, in this case, is estimated as 1.0 GPa, according to
Table 1. It is worth noting that the TPT is from the nontrivial to the trivial
state. When we cross the critical point of the TPT, at V/V0 = 0.940, the Bi-p-
liked ($\Gamma_{8}$) is below the Bi-s-liked ($\Gamma_{6}$) and the compound
is converted into a semiconductor with a trivial bandgap of $\approx$ 0.1 eV
(see Figs. 5(a and b)).

Figure 7:  (a) Band inversion strength ($\delta$) of NiLi2Bi compound in two
approaches, PBE-GGA and TB-mBJ, as a function of V/V0. The black arrows show
the topological phase transitions (TPT) points. (b) Schematic evolution of the
atomic energy levels at $\Gamma$ point under hydrostatic pressure with
(denoted as W SOC) and without SOC (denoted as W/O SOC) within PBE-GGA and TB-
mBJ approaches. The symbols ”+” and ”-” label the parities of the bands at the
$\Gamma$ point. Orange, black, and green rectangles represent nontrivial
states, TPT points, and trivial states, respectively.

To confirm that the topological phase transition occurs in this structure
under hydrostatic tension and compression, we calculated the topological
$\mathbb{Z}_{2}$-index by the evolution lines of Wannier Charge Centres (WCC)
along $k_{y}$ in the two time-reversal invariant planes $k_{x}$= 0 and $\pi$
in the Brillouin zone, as shown in Figs. 5(e-g). For V/V0 ¿ 0.955 states, an
odd number of points in the $k_{x}$=0 plane are crossed by the evolution curve
in any arbitrary line horizontal axis (green line), while an even number of
points in the $k_{x}$= 0.5 plane are crossed (e.g., Figs. 5(f and g)); thus
yielding (z2)0 = 1, (z2)0.5 = 0 and $\mathbb{Z}_{2}$ = [((z2)0 \+ (z2)0.5) mod
2] = 1. Accordingly, NaLi2Bi is a nontrivial topological zero-gap semimetal
with a topological index $\mathbb{Z}_{2}$ = 1 for V/V0 ¿ 0.955. In the same
way, for the states that are after the TPT point V/V0 ¡ 0.955 (e.g., V/V0 =
0.940 state in Fig. 5(e)), the topological index becomes $\mathbb{Z}_{2}$ = 0
((z2)0 = 0, (z2)0.5 = 0 and $\mathbb{Z}_{2}$ = [((z2)0 \+ (z2)0.5) mod 2] =
0), which shows that the compound has transitioned to the trivial topological
phase.

• TB-mBJ approach: It is well known that TB-mBJ gives a very well estimate of
the band gap of solids as is the case for the GW method and the hybrid
functional, but with a less cost [57, 27]. In Fig. 6, we re-performed the
calculations related to Fig. 5 (PBE-GGA approach) in the TB-mBJ approach to
better understand the topological properties of the studied compound. In the
equilibrium state (V/V0 = 1.0), the compound has a trivial bandgap (0.62 eV)
in the presence of SOC, as shown in Fig. 6(a). As can be seen in Fig. 6(c),
this compound is a conventional (trivial topological semiconductor)
semiconductor since the band inversion has not occurred (Bi-p-liked below Bi-
s-liked band). The calculations of the topological $\mathbb{Z}_{2}$-index
($\mathbb{Z}_{2}$ = 0) performed by the evolution lines of WCC are also a
confirmation of this issue (see Fig. 6(f)). With the TB-mBJ approach, in
contrast to the PBE-GGA approach, where the TPT occurs under hydrostatic
compression, the trivial bandgap increases during hydrostatic compression and
the compound maintains its trivial topological phase as a result (see Fig.
6(a)). For example, Figs. 6(b and e) show the orbital projected band
structures and WCC at V/V0 = 0.93, which indicates that the band inversion did
not occur and that the topological $\mathbb{Z}_{2}$-index is zero
($\mathbb{Z}_{2}$ = 0). As can be seen in Fig. 6(a), with the increase of the
hydrostatic pressure of the tension-type, the bandgap decreases and in the
V/V0 = 1.06 state, we reach the critical point of the TPT (see black arrow)
where the BIS is zero. The hydrostatic pressure of this critical point is -3.4
GPa (in this literature tensile hydrostatic pressure is indicated by the
negative sign). It is at this critical point that the TPT from a trivial to a
nontrivial state occurs. A band inversion of sp-type has occurred above the
critical pressure (i.e., V/V0 = 1.07 state), as exhibited in Fig. 6(a and d),
and the BIS has become negative. The compound in this case results in a
nontrivial topological zero-gap semimetal. As seen in Fig. 6(g), the evolution
lines of WCC also show that the topological $\mathbb{Z}_{2}$-index is nonzero
( $\mathbb{Z}_{2}$ = 1).

In Fig. 7(a), we can compare the TPT of PBE-GGA and TB-mBJ approaches based on
how the band inversion strength ($\delta$) changes as a function of V/V0. The
yellow and blue areas represent $\delta$ ¡ 0 (nontrivial state) and $\delta$ ¿
0 (trivial state), respectively. As seen, in PBE-GGA approach, TPT point
corresponds to hydrostatic compression (V/V0 ¡ 1.0), while in TB-mBJ approach,
it corresponds to hydrostatic tension (V/V0 ¿ 1.0) (see black arrows). It
should be noted that the TPT point in the PBE-GGA approach is from a
nontrivial state to a trivial state, while it is reversed in the TB-mBJ
approach. It is possible to define $\delta_{SOC}$ (= [$\Gamma_{8}$ \-
$\Gamma_{7}$]) as another parameter similar to $\delta$, which indicates the
strength of SOC. The changes of this parameter within PBE-GGA and TB-mBJ
approaches are shown in Fig. 4S. The results obtained in two approaches show
that this parameter increases (decreases) with the increase of hydrostatic
pressure of the compression (tensile) type.

To explain the physical mechanism underlying the band inversion and the parity
exchange of this compound, the schematically evolution of the energy levels of
the atomic levels at $\Gamma$ point under hydrostatic pressure (compression
and tension types), in the presence and absence of SOC for the two approaches
studied, is shown in Fig. 7(b). As per the normal band flling order near Fermi
energy, the valence band consists of Bi-p orbitals, while the conduction band
is composed of s-type orbitals (see green rectangle), which, due to the
formation of chemical bonds in this compound, hybridize and split into bonding
and antibonding states in the absence of SOC (i.e., p+ (s+) and p- (s-) with
different parities) [19, 58]. The system is therefore an conventional
semiconductor with trivial bandgap as shown in Fig. 7(b) in both PBE-GGA and
TB-mBJ approaches. In the PBE-GGA approach, when the SOC is taken into account
under hydrostatic pressures, V/V0 ¿ 0.955, (see orange rectangle of stage (I))
the p- state splits into the —p-, $\pm$ j⟩ with a total angular momentum of j
= 3/2 and j = 1/2. On the other hand, the energy band gap reaches zero and the
compound becomes a topological semimetal system. In this case, the parity of s
and p±3/2 states near the Fermi energy have been exchanged and band inversion
occurs (s+ state below p${}^{-}_{\pm 3/2}$ state). At the critical point (V/V0
= 0.955) from stage (II) (see black rectangle), these two states are close
together and at the V/V0 ¡ 0.955, their parity is exchanged with respect to
the stage (I), and the system transitions from the nontrivial state to the
trivial state (see blue rectangle of stage (III)). There is a similar
mechanism in the TB-mBJ approach, but the only difference is that during stage
(I), the compound has a trivial bandgap, and after reaching the TPT point, the
parity of s and p±3/2 states is exchanged and the compound becomes a semimetal
material with nontrivial topological phase (from the trivial state to the
nontrivial state).

##  IV I.V. Summary and outlook

In summary, the mechanical properties and topological phase transition of bi-
alkali pnictogen NaLi2Bi are investigated in the PBE-GGA approach from the
first-principles calculation. Calculations based on phonon modes and elastic
constants indicate that this compound is stable the studied hydrostatic
tension and compression. The results show that applying hydrostatic tension
leads to a kind of mechanical phase transition in the compound so that the
auxetic behaviour of the material changes to non-auxetic behaviour after the
critical pressure of -1.5 GPa (V/V0 = 1.060). According to the Cauchy
pressure, Poisson’s ratio, and Pugh ratio, it has been determined that this
compound has a brittle behaviour in the equilibrium state, and this behaviour
is maintained under the hydrostatic pressures studied in this work. Based on
the study of the electronic band structure in the equilibrium state, it can be
concluded that the NaLi2Bi compound has a zero-gap semimetal in the PBE-GGA
approach, in which the sp-band inversion has occurred, and its topological
index is $\mathbb{Z}_{2}$ = 1. On the other hand, these calculations show
that, in tension-type hydrostatic pressures, the topological phase transition
occurs from the non-trivial phase to the trivial phase, and the critical
pressure is $\approx$1 GPa. Additionally, the TB-mBJ approach indicates that
this compound is a semiconductor with a trivial bandgap. It is possible to
close the bandgap by increasing the tension-type hydrostatic pressure by
$\approx$-3.4 GPa, and after passing through this point the compound
transforms from a trivial phase into a topological nontrivial phase and
becomes a zero-gap semimetal with an inverted band structure (sp-band
inversion occurs). By comparing the two studied approaches in the equilibrium
state of the NaLi2Bi compound, it can be concluded that conventional
approaches like GGA may provide false results in the detection of topological
phases and band inversion in materials. The final point we suggest in this
paper is that uniaxial strain may induce other topological classes suitable
for practical application in this compound, which deserves further study.

##  V acknowledgement

##  VI References

## References

  * Schmeißer  _et al._ [2015a] M. Schmeißer, A. Jankowiak, T. Kamps, J. Kühn, _et al._ , in  _ERL_ , Vol. 15 (2015) p. 97.
  * Cultrera  _et al._ [2015] L. Cultrera, S. Karkare, H. Lee, X. Liu, I. Bazarov, and B. Dunham, Physical Review Special Topics-Accelerators and Beams 18, 113401 (2015).
  * Gulliford  _et al._ [2013] C. Gulliford, A. Bartnik, I. Bazarov, L. Cultrera, J. Dobbins, B. Dunham, F. Gonzalez, S. Karkare, H. Lee, H. Li, _et al._ , Physical Review Special Topics-Accelerators and Beams 16, 073401 (2013).
  * Gulliford  _et al._ [2015] C. Gulliford, A. Bartnik, I. Bazarov, B. Dunham, and L. Cultrera, Applied Physics Letters 106, 094101 (2015).
  * Spicer [1958] W. E. Spicer, Physical review 112, 114 (1958).
  * Yalameha  _et al._ [2021a] S. Yalameha, Z. Nourbakhsh, A. Ramazani, and D. Vashaee, Materials Science and Engineering: B 273, 115430 (2021a).
  * Sklyadneva  _et al._ [2016] I. Y. Sklyadneva, I. P. Rusinov, R. Heid, K.-P. Bohnen, P. M. Echenique, and E. V. Chulkov, Scientific reports 6, 1 (2016).
  * Ruiz-Osés  _et al._ [2014] M. Ruiz-Osés, S. Schubert, K. Attenkofer, I. Ben-Zvi, X. Liang, E. Muller, H. Padmore, T. Rao, T. Vecchione, J. Wong, _et al._ , APL materials 2, 121101 (2014).
  * Cultrera  _et al._ [2017] L. Cultrera, C. Gulliford, A. Bartnik, H. Lee, and I. Bazarov, Journal of Applied Physics 121, 055306 (2017).
  * Feng  _et al._ [2017] J. Feng, S. Karkare, J. Nasiatka, S. Schubert, J. Smedley, and H. Padmore, Journal of Applied Physics 121, 044904 (2017).
  * Ding  _et al._ [2017] Z. Ding, M. Gaowei, J. Sinsheimer, J. Xie, S. Schubert, H. Padmore, E. Muller, and J. Smedley, Journal of Applied Physics 121, 055305 (2017).
  * Schubert  _et al._ [2013] S. Schubert, M. Ruiz-Osés, I. Ben-Zvi, T. Kamps, X. Liang, E. Muller, K. Müller, H. Padmore, T. Rao, X. Tong, _et al._ , Apl Materials 1, 032119 (2013).
  * Khan  _et al._ [2021] Z. Khan, G. Murtaza, A. A. Khan, A. Laref, N. A. Kattan, and M. Haneef, International Journal of Energy Research 45, 7703 (2021).
  * Amador  _et al._ [2021] R. Amador, H.-D. Saßnick, and C. Cocchi, Journal of Physics: Condensed Matter 33, 365502 (2021).
  * Murtaza  _et al._ [2016] G. Murtaza, M. Ullah, N. Ullah, M. Rani, M. Muzammil, R. Khenata, S. M. Ramay, and U. Khan, Bulletin of Materials Science 39, 1581 (2016).
  * Cocchi  _et al._ [2018] C. Cocchi, S. Mistry, M. Schmeißer, J. Kühn, and T. Kamps, Journal of Physics: Condensed Matter 31, 014002 (2018).
  * Kalarasse  _et al._ [2011] L. Kalarasse, B. Bennecer, and F. Kalarasse, Computational materials science 50, 2880 (2011).
  * Jin  _et al._ [2020] L. Jin, X. Zhang, T. He, W. Meng, X. Dai, and G. Liu, Physical Chemistry Chemical Physics 22, 5847 (2020).
  * Yalameha  _et al._ [2018] S. Yalameha, Z. Nourbakhsh, and A. Vaez, Journal of Magnetism and Magnetic Materials 468, 279 (2018).
  * Yalameha  _et al._ [2021b] S. Yalameha, Z. Nourbakhsh, A. Ramazani, and D. Vashaee, Materials Science and Engineering: B 273, 115430 (2021b).
  * Yalameha  _et al._ [2021c] S. Yalameha, Z. Nourbakhsh, and D. Vashaee, Journal of Physics: Condensed Matter 34, 105702 (2021c).
  * Yalameha  _et al._ [2021d] S. Yalameha, Z. Nourbakhsh, A. Ramazani, and D. Vashaee, Nanomaterials 11, 2739 (2021d).
  * Song  _et al._ [2023] X. Song, Y. Zhao, J. Ni, S. Meng, and Z. Dai, Materials Today Physics , 100990 (2023).
  * Blaha  _et al._ [2020] P. Blaha, K. Schwarz, F. Tran, R. Laskowski, G. K. Madsen, and L. D. Marks, The Journal of Chemical Physics 152, 074101 (2020).
  * Petersen  _et al._ [2000] M. Petersen, F. Wagner, L. Hufnagel, M. Scheffler, P. Blaha, and K. Schwarz, Computer Physics Communications 126, 294 (2000).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Physical review letters 77, 3865 (1996).
  * Jiang [2013] H. Jiang, The Journal of chemical physics 138, 134115 (2013).
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, Scripta Materialia 108, 1 (2015).
  * Kresse  _et al._ [1995] G. Kresse, J. Furthmüller, and J. Hafner, EPL (Europhysics Letters) 32, 729 (1995).
  * Jamal  _et al._ [2018] M. Jamal, M. Bilal, I. Ahmad, and S. Jalali-Asadabadi, Journal of Alloys and Compounds 735, 569 (2018).
  * Yalameha  _et al._ [2022] S. Yalameha, Z. Nourbakhsh, and D. Vashaee, Computer Physics Communications 271, 108195 (2022).
  * Mostofi  _et al._ [2008] A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Computer physics communications 178, 685 (2008).
  * Gresch  _et al._ [2017] D. Gresch, G. Autes, O. V. Yazyev, M. Troyer, D. Vanderbilt, B. A. Bernevig, and A. A. Soluyanov, Physical Review B 95, 075146 (2017).
  * Murnaghan [1944] F. D. Murnaghan, Proceedings of the National Academy of Sciences 30, 244 (1944).
  * Saeidi  _et al._ [2019] P. Saeidi, S. Yalameha, _et al._ , Computational Condensed Matter 18, e00358 (2019).
  * Pettifor [1992] D. Pettifor, Materials science and technology 8, 345 (1992).
  * Frantsevich [1982] I. N. Frantsevich, Reference book (1982).
  * Pugh [1954] S. Pugh, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 45, 823 (1954).
  * Savin  _et al._ [1992] A. Savin, H.-J. Flad, J. Flad, H. Preuss, and H. G. von Schnering, Angewandte Chemie International Edition in English 31, 185 (1992).
  * Anderson and Demarest Jr [1971] O. L. Anderson and H. H. Demarest Jr, Journal of Geophysical Research 76, 1349 (1971).
  * Hadi  _et al._ [2017] M. Hadi, M. Roknuzzaman, A. Chroneos, S. Naqib, A. Islam, R. Vovk, and K. Ostrikov, Computational Materials Science 137, 318 (2017).
  * Vitos  _et al._ [2003] L. Vitos, P. A. Korzhavyi, and B. Johansson, Nature Materials 2, 25 (2003).
  * Puttlitz and Stalter [2004] K. J. Puttlitz and K. A. Stalter,  _Handbook of lead-free solder technology for microelectronic assemblies_ (CRC Press, 2004).
  * Greaves  _et al._ [2011] G. N. Greaves, A. L. Greer, R. S. Lakes, and T. Rouxel, Nature materials 10, 823 (2011).
  * Evans and Alderson [2000] K. E. Evans and A. Alderson, Advanced materials 12, 617 (2000).
  * Polash  _et al._ [2021] M. M. H. Polash, S. Yalameha, H. Zhou, K. Ahadi, Z. Nourbakhsh, and D. Vashaee, Materials Science and Engineering: R: Reports 145, 100620 (2021).
  * Veldhorst  _et al._ [2012] M. Veldhorst, C. Molenaar, X. Wang, H. Hilgenkamp, and A. Brinkman, Applied physics letters 100, 072602 (2012).
  * Maassen and Lundstrom [2013] J. Maassen and M. Lundstrom, Applied Physics Letters 102, 093103 (2013).
  * Yalameha and Nourbakhsh [2020] S. Yalameha and Z. Nourbakhsh, Journal of Physics: Condensed Matter 32, 295502 (2020).
  * Harrison  _et al._ [2014] S. Harrison, L. Collins-McIntyre, S. Li, A. Baker, L. Shelford, Y. Huo, A. Pushp, S. Parkin, J. Harris, E. Arenholz, _et al._ , Journal of Applied Physics 115, 023904 (2014).
  * Kim  _et al._ [2012] M. Kim, C. H. Kim, H.-S. Kim, and J. Ihm, Proceedings of the National Academy of Sciences 109, 671 (2012).
  * Monserrat and Vanderbilt [2016] B. Monserrat and D. Vanderbilt, Physical review letters 117, 226801 (2016).
  * Saeidi and Nourbakhsh [2018] P. Saeidi and Z. Nourbakhsh, Journal of Magnetism and Magnetic Materials 451, 681 (2018).
  * Zhu  _et al._ [2012] Z. Zhu, Y. Cheng, and U. Schwingenschlögl, Physical Review B 85, 235401 (2012).
  * Zhang  _et al._ [2016] X. Zhang, Z. Hou, Y. Wang, G. Xu, C. Shi, E. Liu, X. Xi, W. Wang, G. Wu, and X.-x. Zhang, Scientific reports 6, 23172 (2016).
  * Schmeißer  _et al._ [2015b] M. Schmeißer, A. Jankowiak, T. Kamps, J. Kühn, _et al._ , in  _ERL_ , Vol. 15 (2015) p. 97.
  * Perdew  _et al._ [2008] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Physical review letters 100, 136406 (2008).
  * Yalameha and Nourbakhsh [2022] S. Yalameha and Z. Nourbakhsh, Materials Science and Engineering: B 281, 115742 (2022).


# Construction of explicit symplectic integrators in general relativity. III.
Reissner-Nordström-(anti)-de Sitter black holes

Ying Wang1,2, Wei Sun1, Fuyao Liu1, Xin Wu1,2,3,† 1\. School of Mathematics,
Physics and Statistics, Shanghai University of Engineering Science, Shanghai
201620, China
2\. Center of Application and Research of Computational Physics, Shanghai
University of Engineering Science, Shanghai 201620, China
3\. Guangxi Key Laboratory for Relativistic Astrophysics, Guangxi University,
Nanning 530004, China  Emails<EMAIL_ADDRESS>(Y. W.),
<EMAIL_ADDRESS>(W. S<EMAIL_ADDRESS>(F. L.); ${\dagger}$
Corresponding Author<EMAIL_ADDRESS>(X. W.)

###### Abstract

We give a possible splitting method to a Hamiltonian for the description of
charged particles moving around the Reissner-Nordström-(anti)-de Sitter black
hole with an external magnetic field. This Hamiltonian can be separated into
six analytical solvable pieces, whose solutions are explicit functions of
proper time. In this case, second- and fourth-order explicit symplectic
integrators are easily available. They exhibit excellent long-term behavior in
maintaining the boundness of Hamiltonian errors regardless of ordered or
chaotic orbits if appropriate step-sizes are chosen. Under some circumstances,
an increase of positive cosmological constant gives rise to strengthening the
extent of chaos from the global phase space; namely, chaos of charged
particles occurs easily for the accelerated expansion of the universe.
However, an increase of the magnitude of negative cosmological constant does
not. The different contributions on chaos are because the cosmological
constant acts as a repulsive force in the Reissner-Nordström-de Sitter black
hole, but an attractive force in the Reissner-Nordström-anti-de Sitter black
hole.

_Unified Astronomy Thesaurus concepts_ : Black hole physics (159);
Computational methods (1965); Computational astronomy (293); Celestial
mechanics (211)

##  1 Introduction

Many exact solutions on black holes can be derived from Einstein’s general
relativistic field equations. The first statical spherically symmetric
solution is the Schwarzschild metric (Schwarzschild 1916). When the mass
source in the origin has an electric charge, the Schwarzschild metric is the
famous Reissner-Nordström (RN) solution (Reissner 1916; Nordström 1918). When
a positive constant $\Lambda$ for describing the static cosmological model
introduced by Einstein is included in the Schwarzschild metric, the
Schwarzschild-de Sitter metric is provided. With the inclusion of cosmological
constant $\Lambda$, the Reissner-Nordström-(anti)-de Sitter [RN-(A)dS]
solution (Reissner 1916; Kottler 1918; Nordström 1918) is also obtained.

A small positive cosmological constant acts as a gravitationally repulsive
energy component, and models dark energy for dominating the energy budget of
the universe leads to an accelerated expansion of the universe (Ashtekar
2017). This result is strongly supported by observational evidences (Riess et
al. 1998; Perlmutter et al. 1999; Ade et al. 2016). In fact, the cosmological
constant describes a weak field gravity for the dynamics of galaxies, galaxy
groups and clusters, and is supported by the parameters of galactic halos and
of a sample of galaxy groups (Gurzadyan 2019; Gurzadyan $\&$ Stepanian 2019).
Surprisingly, the theory of gravitational waves from isolated systems links up
with the presence of a positive cosmological constant (Ashtekar et al. 2016;
Ashtekar 2017). Although the known observational results do not support the
existence of a negative cosmological constant, the negative cosmological
constant becomes useful in the study of the black hole thermodynamics under
the existence of quantum fields in gravitational fields. Appropriate boundary
conditions should be imposed in black holes so that the black holes remain
thermally stable. The AdS boundary, as a reflecting wall, can cause
asymptotically AdS black holes to remain thermally stable. The phase
transition between Schwarzschild AdS black holes and thermal AdS space was
discovered by Hawking $\&$ Page (1983). Finite temperature configurations in
the decoupled field theory are associated to black hole configurations in AdS
spacetimes. The black holes yield the Hawking radiation in the AdS spacetime
(Maldacena 1998). On the other hand, the cosmological constant affects some
properties of accretion discs orbiting black holes in quasars and active
galactic nuclei (Stuchlík 2005). These accretion discs are developed through
the accretion process of dust and gases around compact objects. The mass of
the accreting body increases during the accretion process, but it could also
decrease (Jamil et al. 2008; Jamil 2009). Karkowski $\&$ Malec (2013) found
that the damping effect on the mass accretion rate is weaker for positive
values of the cosmological constant than for negative values of the
cosmological constant, and it is related to the ratio of the dark energy
density to the fluid density. The authors mainly focused on the question of
how the cosmological constant exerts influences on the possible accretion
rate. Ficek (2015) considered relativistic Bondi-type accretion in the
RN-(A)dS spacetime.

The Schwarzschild, RN and RN-(A)dS spacetimes are highly nonlinear, but they
are integrable and non-chaotic. Magnetic fields can destroy the integrability
of these spacetimes. In fact, a strong magnetic field with hundred Gauss is
indeed present in the vicinity of the supermassive black hole at the Galactic
center (Eatough et al. 2013). The presence of an external magnetic field was
shown through high-frequency quasiperiodic oscillations and relativistic jets
observed in microquasars (Stuchlík et al. 2013; Kološ et al. 2017). It can
influence oscillations in the accretion disk and the creation of jets. The
accretion of charged particles from the accretion disk goes towards the black
hole due to a radiation-reaction force. The external magnetic field that is
usually considered is so weak that it does not perturb a spacetime metric;
namely, the external magnetic field does not appear in the spacetime metric
and is not a solution of the field equations of general relativity. However,
such weak magnetic field will essentially affect motions of electrons or
protons when the ratios of charges of electrons or protons to masses of
electrons or protons are large (Tursunov et al. 2018). Of course, some
magnetic fields can perturb spacetime metrics, such as the magnetized Ernst-
Schwarzschild spacetime geometry (Ernst 1976). The chaotic motions of charged
particles around the Schwarzschild and Kerr black holes were confirmed in a
series of papers (Nakamura $\&$ Ishizuka 1993; Takahashi $\&$ Koyama 2009;
Kopáček et al. 2010; Kopáček $\&$ Karas 2014; Stuchlík $\&$ Kološ 2016;
Kopáček $\&$ Karas 2018; Pánis et al. 2019; Stuchlík et al. 2020; Yi $\&$ Wu
2020). In particular, the chaotic dynamics of neutral particles in the
magnetized Ernst-Schwarzschild spacetim is possible due to the magnetic field
acting as a gravitational effect although no electromagnetic forces exist
(Karas $\&$ Vokroulflický 1992; Li $\&$ Wu 2019). Besides external magnetic
fields, spin of particles is an important source for inducing chaos. Chaos can
occur when a spinning test particle moves in the spacetime background of a
Kerr black hole (Lukes-Gerakopoulos et al. 2016). Resonances and chaos for
spinning test particles in the Schwarzschild and Kerr backgrounds were shown
by Brink et al. (2015a, b) and Zelenka et al. (2020). Some spacetime
backgrounds, e.g., the Manko, Sanabria-Gómez, Manko (MSM) metric can be
chaotic (Dubeibe et al. 2007; Han 2008; Seyrich $\&$ Lukes-Gerakopoulos 2012).
The word “chaos” is a nonlinear behavior of a dynamical system and describes
the system with exponentially sensitive dependence on initial conditions. The
chaotic or regular feature of orbital motions can be reflected on the
gravitational waves (Zelenka et al. 2019).

Detecting the chaotical behavior requires that a numerical integration scheme
have good stability and precision, i.e., provide reliable results to the
trajectories. The Schwarzschild, RN and RN-(A)dS spacetimes are integrable,
but their solutions are unknown. Thus, numerical schemes are still necessary
to solve these systems. It is well known that symplectic integrators are the
most appropriate solvers for long-term evolution of Hamiltonian systems.
Clearly, Hamiltonians of the Schwarzschild, RN and RN-(A)dS spacetimes are not
separations of variables. Perhaps each of the Hamiltonians may be split into
two integrable parts, but no explicit functions of proper time can be given to
analytical solutions of the two parts in the general case. Explicit symplectic
integrators (Swope et al. 1982; Wisdom 1982; Ruth 1983; Forest $\&$ Ruth 1990;
Wisdom $\&$ Holman 1991111The Wisdom-Holman method is not strictly explicit
because an implicit iterative method is required to calculate the eccentric
anomaly in the solution of the Kepler main part.) that have been developed in
the solar system and share good long-term numerical performance in the
preservation of geometric structures and constants of motion222The
preservation of constants of motion (such as energy integral) in symplectic
integrators does not mean that the constants are exactly conserved, but means
that errors of the constants are bounded. The energy-preserving integrators
(Bacchini et al. 2018a, b; Hu et al. 2019, 2021) and the projection methods
(Fukushima 2003a,b,c, 2004; Ma et al. 2008; Wang et al. 2016, 2018; Deng et
al. 2020) are not symplectic, whereas they do exactly conserve energy. become
useless in this case. However, completely implicit or explicit and implicit
combined symplectic integrators should be available without doubt (Brown 2006;
Preto $\&$ Saha 2009; Kopáček et al. 2010; Lubich et al. 2010; Zhong et al.
2010; Mei et al. 2013a, 2013b; Tsang et al. 2015). In most cases, explicit
algorithms have advantage over same order implicit schemes at the expense of
computational time. Considering this point, we proposed second- and fourth-
order explicit symplectic integration algorithms for the Schwarzschild
spacetime by separating the Hamiltonian of this spacetime into four parts with
analytical solutions as explicit functions of proper time in Paper I (Wang et
al. 2021a). When five parts similar to the four separable pieces in the
Hamiltonian of the Schwarzschild metric are given to a splitting form of the
Hamiltonian of the RN black hole with an external magnetic field, explicit
symplectic integrators were also set up in Paper II (Wang et al. 2021b). In
other words, a fact concluded from Papers I and II is that the construction of
explicit symplectic integrations is possible for a curved spacetime when a
Hamiltonian of the spacetime is separated into more parts with analytical
solutions as explicit functions of proper time. How to split the Hamiltonian
depends on what this spacetime is.

Following Papers I and II, we want to know in this paper how explicit
symplectic integrators are designed for the RN-(A)dS spacetime with an
external magnetic field. We are particularly interested in investigating how
the cosmological constant effects the regular and chaotic dynamical behavior
of charged particles moving around the black hole. For the sake of these
purposes, we introduce the relativistic gravitational field in Sect. 2. In
Sect. 3, we give new explicit symplectic integrators and check their numerical
performance. Then, a new explicit symplectic method with an appropriate step
size is applied to survey the effect of the cosmological constant on the
dynamics of orbits of charged particles. Finally, the main results are
concluded in Sect. 4.

##  2 RN black hole with cosmological constant and external magnetic field

In dimensionless spherical-like coordinates $(t,r,\theta,\phi)$, the RN black
hole with electric charge $Q$ and cosmological constant $\Lambda$ describing
the expansion of the universe is given in the following spacetime (Ahmed et
al. 2016)

$\displaystyle-\tau^{2}$ $\displaystyle=$ $\displaystyle
ds^{2}=g_{\alpha\beta}dx^{\alpha}dx^{\beta}$ (1) $\displaystyle=$
$\displaystyle-(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}}-\frac{\Lambda}{3}r^{2})dt^{2}$
$\displaystyle+(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}}-\frac{\Lambda}{3}r^{2})^{-1}dr^{2}$
$\displaystyle+r^{2}d\theta^{2}+r^{2}\sin^{2}\theta d\phi^{2}.$

Here, geometrized units $c=G=M=1$ are given to the speed of light $c$, the
constant of gravity $G$, and the mass of black hole $M$. Dimensionless
operations act on proper time $\tau$, coordinate time $t$, separation $r$,
charge $Q$, and cosmological constant $\Lambda$. Namely, $\tau\rightarrow
M\tau$, $t\rightarrow Mt$, $r\rightarrow Mr$, $Q\rightarrow MQ$ and
$\Lambda\rightarrow\Lambda/M^{2}$. The four-velocity
$\mathbf{U}=U^{\mu}=\dot{x}^{\mu}$ is also dimensionless and satisfies the
constraint

$\mathbf{U}\bullet\mathbf{U}=U^{\alpha}U_{\alpha}=g_{\mu\nu}\dot{x}^{\mu}\dot{x}^{\nu}=-1.$
(2)

Unlike the Euclid inner product, the Riemann inner product $\bullet$ in Eq.
(2) can be negative.

$\Lambda=0$ indicates that Equation (1) is an asymptotically flat metric. A
non-vanishing cosmological constant $\Lambda$ represents a weak field gravity
proper for the dynamics of galaxies, galaxy groups and clusters. $\Lambda>0$
and $\Lambda<0$ correspond to asymptotically de Sitter (dS) and anti-de Sitter
(AdS) spacetimes, respectively. Without loss of generality, we consider a
relation between realistic and dimensionless values of positive cosmological
constant. The realistic value of cosmological constant obtained by Planck data
is $\Lambda_{PL}=1.11\times 10^{-52}m^{-2}=1.11\times 10^{-46}km^{-2}$ (Ade et
al. 2016). Various realistic values of cosmological constant for dark matter
considering the data on galaxy systems from galaxy pairs to galaxy clusters
were given in the range $\Lambda=10^{-52}m^{-2}\sim 10^{-50}m^{-2}$ (Gurzadyan
$\&$ Stepanian 2019). Suppose that the value $\Lambda_{0}=10^{-46}km^{-2}$
considered in a paper of Xu $\&$ Wang (2017) is regarded as a realistic value
of cosmological constant, a dimensionless value of cosmological constant can
be estimated in terms of $\Lambda=\Lambda_{0}M^{2}G^{2}/c^{4}$. Obviously, it
depends on the mass of a black hole. For the supermassive black hole candidate
in the center of the giant elliptical galaxy M87 having mass $M=(6.5\pm
0.7)\times 10^{9}M_{\odot}$ with $M_{\odot}$ being the Sun’s mass (EHT
Collaboration et al. 2019a, b, c), we obtain dimensionless cosmological
constant $\Lambda=\Lambda_{0}M^{2}G^{2}/c^{4}=10^{-46}km^{-2}\times[(6.5\pm
0.7)\times 10^{9}M_{\odot}G/c^{2}]^{2}=10^{-46}km^{-2}\times(6.5\pm
0.7)^{2}\times 10^{18}\times 1.475^{2}km^{2}=7.3\times 10^{-27}\sim 1.1\times
10^{-26}$. Such a small dimensionless value does not change the spacetime
geometry of black hole if the separation $r$ is not large enough. The effect
of cosmological constant on the spacetime should also be negligible for the
$150M_{\odot}$ binary black hole merger GW190521 (Abbott et al. 2020a, 2020b).
The dimensionless value of cosmological constant increases with an increase of
the black hole’s mass. If the black hole’s mass is $M=(1/1.475)\times
10^{20}M_{\odot}$, we have dimensionless cosmological constant
$\Lambda=10^{-6}$. For $M=(1/1.475)\times 10^{21}M_{\odot}$,
$\Lambda=10^{-4}$. Given $M=2.14391705774127\times 10^{22}M_{\odot}$,
$\Lambda=0.1$. When $M=(1/1.475)\times 10^{23}M_{\odot}$, $\Lambda=1$. These
larger dimensionless values should exert explicit influences on the spacetime
geometries. The larger dimensionless values of cosmological constant (e.g.
$\Lambda=-0.075,-0.1,0.1$) in the references (Ahmed et al. 2016; Shahzad $\&$
Jawad 2019; Shafiq et al. 2020) may be taken for the existence of extremely
supermassive black holes. In what follows, various values of cosmological
constant are dimensionless. Equation (1) is called the RN-(A)dS spacetime. It
has at most three event horizons (Ahmed et al. 2016; Shahzad $\&$ Jawad 2019).

Seen from the expression $g_{00}=-1-2\varphi$ ($\varphi$ is the gravitational
potential) in Eq. (1), the potential $-1/r$ plays a role in an attractive
force proportional to $-1/r^{2}$. Similarly, the cosmological constant makes
the potential $-\Lambda r^{2}/6$ yield an attractive force proportional to
$\Lambda r/3$ for $\Lambda<0$, but a repulsive force for $\Lambda>0$. The
force is independent of the black hole’s mass. However, the black hole’s
charge leads to the potential $Q^{2}/(2r^{2})$ giving a repulsive force
proportional to $Q^{2}/r^{3}$. As the distance increases, the attraction from
the heavy body and the repulsive force from the charge decrease; the latter is
smaller than the former. However, the attraction or repulsiveness from the
cosmological constant grows. It eventually dominates the other two forces when
the distance is large enough.

The metric corresponds to a Lagrangian system

$\mathcal{L}=\frac{1}{2}(\frac{ds}{d\tau})^{2}=\frac{1}{2}g_{\mu\nu}\dot{x}^{\mu}\dot{x}^{\nu}.$
(3)

It defines a covariant generalized momentum

$p_{\mu}=\frac{\partial\mathcal{L}}{\partial\dot{x}^{\mu}}=g_{\mu\nu}\dot{x}^{\nu}.$
(4)

Since the Lagrangian does not explicitly depend on $t$ and $\phi$, the Euler-
Lagrangian equations show that two constant momentum components with respect
to proper time $\tau$ are

$\displaystyle p_{t}$ $\displaystyle=$
$\displaystyle-(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}}-\frac{\Lambda}{3}r^{2})\dot{t}=-\mathbb{E},$
(5) $\displaystyle p_{\phi}$ $\displaystyle=$ $\displaystyle
r^{2}\sin^{2}\theta\dot{\phi}=\mathbb{L}.$ (6)

$\mathbb{E}$ is the energy of a test particle moving in the gravitational
field, and $\mathbb{L}$ denotes the angular momentum of the particle.

Transform the Lagrangian into a Hamiltonian

$\displaystyle H$ $\displaystyle=$
$\displaystyle\mathbf{U}\bullet\mathbf{p}-\mathcal{L}=\frac{1}{2}g^{\mu\nu}p_{\mu}p_{\nu}$
(7) $\displaystyle=$
$\displaystyle-\frac{1}{2}(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}}-\frac{\Lambda}{3}r^{2})^{-1}\mathbb{E}^{2}$
$\displaystyle+\frac{1}{2}(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}}-\frac{\Lambda}{3}r^{2})p^{2}_{r}+\frac{1}{2}\frac{p^{2}_{\theta}}{r^{2}}$
$\displaystyle+\frac{1}{2}\frac{\mathbb{L}^{2}}{r^{2}\sin^{2}\theta}.$

Due to Eq. (2), this Hamiltonian is always identical to a given constant

$H=-\frac{1}{2}.$ (8)

In fact, a second integral (Carter 1968) excluding $E$ and $\mathbb{L}$ can be
found by the variables separated in the Hamilton-Jacobi equation. Thus, this
system is integrable.

Let the black hole suffer from the interaction of an external magnetic field.
The electromagnetic field is described by a four-vector potential with two
nonzero covariant components (Felice $\&$ Sorge 2003; Wang $\&$ Zhang 2020)

$\displaystyle A_{t}$ $\displaystyle=$ $\displaystyle-\frac{Q}{r},$ (9)
$\displaystyle A_{\phi}$ $\displaystyle=$
$\displaystyle\frac{B}{2}g_{\phi\phi}=\frac{B}{2}r^{2}\sin^{2}\theta.$ (10)

$A_{t}$ is the Coulomb part of the electromagnetic potential. It does not
appear in the metric $g_{\alpha\beta}$, but it should be given because a
charged particle and the black hole’s charge $Q$ yields a Coulomb force. This
Coulomb force is compensated by $A_{t}$. $B$ denotes the strength of the
magnetic field parallel to the $z$ axis. As claimed in the Introduction, $B$
is so weak that it does not affect the metric but affects the motion of
charged particles. The motion of particles with charge $q$ around the
magnetized black hole is determined by the Hamiltonian

$\displaystyle K$ $\displaystyle=$
$\displaystyle\frac{1}{2}g^{\mu\nu}(p_{\mu}-qA_{\mu})(p_{\nu}-qA_{\nu})$ (11)
$\displaystyle=$
$\displaystyle-\frac{1}{2}(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}}-\frac{\Lambda}{3}r^{2})^{-1}(E-\frac{\tilde{Q}}{r})^{2}$
$\displaystyle+\frac{1}{2}(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}}-\frac{\Lambda}{3}r^{2})p^{2}_{r}+\frac{1}{2}\frac{p^{2}_{\theta}}{r^{2}}$
$\displaystyle+\frac{L^{2}}{2r^{2}\sin^{2}\theta}+\frac{1}{8}\beta^{2}r^{2}\sin^{2}\theta-\frac{\beta
L}{2},$

where $Q^{*}=qQ$ and $\beta=qB$. Dimensionless operations (see also Yi $\&$ Wu
2020) are given as follows: $K\rightarrow Km^{2}_{p}$, $E\rightarrow Em_{p}$,
$p_{r}\rightarrow p_{r}m_{p}$, $p_{\theta}\rightarrow p_{\theta}Mm_{p}$,
$L\rightarrow LMm_{p}$, $B\rightarrow B/M$ and $q\rightarrow qm_{p}$, where
$m_{p}$ is the particle’s mass. Unlike $\mathbb{E}$ and $\mathbb{L}$ in Eqs.
(5) and (6), the energy $E$ and angular momentum $L$ are written as

$\displaystyle E$ $\displaystyle=$
$\displaystyle(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}}-\frac{\Lambda}{3}r^{2})\dot{t}+\frac{Q^{*}}{r},$
(12) $\displaystyle L$ $\displaystyle=$ $\displaystyle
r^{2}\sin^{2}\theta\dot{\phi}+\frac{1}{2}\beta r^{2}\sin^{2}\theta.$ (13)

The magnetic field causes the radial potential $\beta^{2}r^{2}/8$ to yield an
attractive force proportional to $-\beta^{2}r/4$. The attraction grows with
the distance. However, the radial potential $L^{2}/(2r^{2})$ from the angular
momentum gives a repulsive force proportional to $L^{2}/r^{3}$. In fact, this
repulsive force is a centrifugal force of inertia. As Wang et al. (2021b)
claimed, the Coulomb term gives the particle a repulsive force effect for
$Q^{*}>0$, but a gravitational force effect for $Q^{*}<0$. It is more
important than the repulsive force effect from the black hole’s charge.
Although the potential from the black hole’s charge and the potential from the
angular momentum yield the repulsive forces, the latter repulsive forces vary
their directions when the charged particles’ motions are not restricted to the
radial motions. Similarly, the attractive forces from the magnetic field
unlike those from the negative cosmological constant (or those from the heavy
body) have also different directions in the general case.

The Hamiltonian $K$ is still equal to the given constant

$K=-\frac{1}{2}.$ (14)

However, no second integral excluding the two integrals $E$ and $L$ exists in
the case. Thus, this Hamiltonian is non-integrable. Its dynamics is
complicated and mainly employs numerical techniques to be investigated.

##  3 Numerical investigations

First, new explicit symplectic integrators are specifically designed for the
Hamiltonian system (11). Then, the performance of the new algorithms is
evaluated. Finally, the dynamics of order and chaos of orbits in the
Hamiltonian is surveyed, and the relation between the chaotic behavior and the
cosmological constant is mainly discussed.

###  3.1 Construction of explicit symplectic integrators

In Papers I and II (Wang et al. 2021a, b), the explicit symplectic integrators
were proposed for the Hamiltonian of Schwarzschild’s metric with four
analytical integrable splitting parts, and the Hamiltonian of RN’s metric with
five analytical integrable separable parts. These successful examples in
general relativity tell us that the key of the construction of an explicit
symplectic integrator lies in splitting the considered Hamiltonian in an
appropriate way. Not only the splitting parts in this Hamiltonian should be
analytically solvable, but also all analytical solutions in the separable
parts should be explicit functions of proper time $\tau$. Noticing this idea,
we begin our work.

Split the Hamiltonian (11) into six parts

$K=K_{1}+K_{2}+K_{3}+K_{4}+K_{5}+K_{6}.$ (15)

The sub-Hamiltonians are expressed as

$\displaystyle K_{1}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}}-\frac{\Lambda}{3}r^{2})^{-1}(E-\frac{Q^{*}}{r})^{2}$
(16)
$\displaystyle+\frac{L^{2}}{2r^{2}\sin^{2}\theta}+\frac{1}{8}\beta^{2}r^{2}\sin^{2}\theta-\frac{\beta
L}{2},$ $\displaystyle K_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}p^{2}_{r},$ (17) $\displaystyle K_{3}$
$\displaystyle=$ $\displaystyle-\frac{1}{r}p^{2}_{r},$ (18) $\displaystyle
K_{4}$ $\displaystyle=$ $\displaystyle\frac{p^{2}_{\theta}}{2r^{2}},$ (19)
$\displaystyle K_{5}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{Q^{2}}{r^{2}}p^{2}_{r},$ (20) $\displaystyle
K_{6}$ $\displaystyle=$ $\displaystyle-\frac{\Lambda}{6}r^{2}p^{2}_{r}.$ (21)

$K_{2}$, $K_{3}$, $K_{4}$ and $K_{5}$ are consistent with those in the
decomposing method of the RN black hole in Paper II (Wang et al. 2021b).

$K_{1}$, $\cdots$, $K_{6}$ correspond to differential operators $\mathcal{A}$,
$\mathcal{B}$, $\mathcal{C}$, $\mathcal{D}$, $\mathcal{E}$ and $\mathcal{F}$,
respectively. These operators are written in the forms

$\displaystyle\mathcal{A}$ $\displaystyle=$ $\displaystyle-\frac{\partial
K_{1}}{\partial r}\frac{\partial}{\partial p_{r}}-\frac{\partial
K_{1}}{\partial\theta}\frac{\partial}{\partial p_{\theta}}$ (22)
$\displaystyle=$
$\displaystyle[-(E-\frac{Q^{*}}{r})^{2}(\frac{1}{r^{2}}-\frac{Q^{2}}{r^{3}}-\frac{\Lambda}{3}r)$
$\displaystyle\cdot(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}}-\frac{\Lambda}{3}r^{2})^{-2}$
$\displaystyle+\frac{Q^{*}}{r^{2}}(E-\frac{Q^{*}}{r})(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}}-\frac{\Lambda}{3}r^{2})^{-1}$
$\displaystyle+\frac{L^{2}}{r^{3}\sin^{2}\theta}-\frac{\beta^{2}}{4}r\sin^{2}\theta]\frac{\partial}{\partial
p_{r}}$
$\displaystyle+(\frac{L^{2}\cos\theta}{r^{2}\sin^{3}\theta}-\frac{1}{8}\beta^{2}r^{2}\sin(2\theta))\frac{\partial}{\partial
p_{\theta}}$ $\displaystyle=$ $\displaystyle
f_{1}(r,\theta)\frac{\partial}{\partial
p_{r}}+f_{2}(r,\theta)\frac{\partial}{\partial p_{\theta}},$
$\displaystyle\mathcal{B}$ $\displaystyle=$ $\displaystyle
p_{r}\frac{\partial}{\partial r},$ (23) $\displaystyle\mathcal{C}$
$\displaystyle=$ $\displaystyle-\frac{2}{r}p_{r}\frac{\partial}{\partial
r}-\frac{p^{2}_{r}}{r^{2}}\frac{\partial}{\partial p_{r}},$ (24)
$\displaystyle\mathcal{D}$ $\displaystyle=$
$\displaystyle\frac{p_{\theta}}{r^{2}}\frac{\partial}{\partial\theta}+\frac{p^{2}_{\theta}}{r^{3}}\frac{\partial}{\partial
p_{r}},$ (25) $\displaystyle\mathcal{E}$ $\displaystyle=$
$\displaystyle\frac{Q^{2}}{r^{2}}p_{r}\frac{\partial}{\partial
r}+Q^{2}\frac{p^{2}_{r}}{r^{3}}\frac{\partial}{\partial p_{r}},$ (26)
$\displaystyle\mathcal{F}$ $\displaystyle=$
$\displaystyle-\frac{\Lambda}{3}r^{2}p_{r}\frac{\partial}{\partial
r}+\frac{\Lambda}{3}rp^{2}_{r}\frac{\partial}{\partial p_{r}}.$ (27)

Using them to act on the solutions
$\textbf{z}(0)=(r_{0},\theta_{0},p_{r0},p_{\theta 0})$ at the beginning of
proper time step $h$, we can obtain their corresponding analytical solutions
$\textbf{z}=(r,\theta,p_{r},p_{\theta})$ over a step size $h$. All the
analytical solutions are indeed explicit functions of the proper time step
$h$. For example, the analytical solutions z for the sub-Hamiltonian $K_{1}$
are expressed in an exponential operator $e^{h\mathcal{A}}$ as
$\textbf{z}=(p_{r},p_{\theta})=e^{h\mathcal{A}}\textbf{z}(0)$. That is,
$p_{r}=p_{r0}+hf_{1}(r_{0},\theta_{0})$ and $p_{\theta}=p_{\theta
0}+hf_{2}(r_{0},\theta_{0})$. These exponential operators symmetrically
compose a second-order explicit symplectic integrator

$\displaystyle S^{K}_{2}(h)$ $\displaystyle=$ $\displaystyle
e^{\frac{h}{2}\mathcal{F}}\circ e^{\frac{h}{2}\mathcal{E}}\circ
e^{\frac{h}{2}\mathcal{D}}\circ e^{\frac{h}{2}\mathcal{C}}\circ
e^{\frac{h}{2}\mathcal{B}}\circ e^{h\mathcal{A}}$ (28) $\displaystyle\circ
e^{\frac{h}{2}\mathcal{B}}\circ e^{\frac{h}{2}\mathcal{C}}\circ
e^{\frac{h}{2}\mathcal{D}}\circ e^{\frac{h}{2}\mathcal{E}}\circ
e^{\frac{h}{2}\mathcal{F}}.$

This explicit symplectic method is specifically designed for the Hamiltonian
(11). It is easily used to yield a fourth-order symplectic scheme of Yoshida
(1990)

$S^{K}_{4}(h)=S^{K}_{2}(\gamma h)\circ S^{K}_{2}(\delta h)\circ
S^{K}_{2}(\gamma h),$ (29)

where $\delta=1-2\gamma$ and $\gamma=1/(2-\sqrt[3]{2})$.

Here are three points illustrated. (i) the Hamiltonian (11) becomes the RN
Hamiltonian (7) for the case of $\beta=0$. It is natural that the newly
established algorithms are appropriate for the RN problem. (ii) The splitting
technique to the Hamiltonian (11) may not be unique. Perhaps there may be
other splitting methods that satisfy the requirement for the construction of
such similar explicit symplectic integration algorithms. (iii) The present
operator splitting method can be applied to an extended version of the
Hamiltonian (11) as follows:

$\displaystyle\Gamma$ $\displaystyle=$ $\displaystyle
f_{1}(r,\theta)+g^{rr}p^{2}_{r}+g^{\theta\theta}p^{2}_{\theta\theta},$ (30)
$\displaystyle g^{rr}$ $\displaystyle=$
$\displaystyle(a_{0}+a_{1}r+\cdots+a_{n}r^{n}+\frac{b_{1}}{r}+\frac{b_{2}}{r^{2}}$
$\displaystyle+\cdots+\frac{b_{m}}{r^{m}})f_{2}(\theta),$ $\displaystyle
g^{\theta\theta}$ $\displaystyle=$ $\displaystyle f_{3}(r),$

where $f_{1}$, $f_{2}$ and $f_{3}$ are arbitrary continuous differentiable
functions when $r$ is restricted to $r>2$, and $a_{0}$, $\cdots$, $a_{n}$,
$b_{1}$, $\cdots$, $b_{m}$ are constant parameters. For the Kerr metric, the
operator splitting method is not available because the denominators of
$g^{rr}$ and $g^{\theta\theta}$ are $r^{2}+a^{2}\cos^{2}\theta$. Similarly,
the method becomes useless for the magnetized Ernst-Schwarzschild spacetime
(Ernst 1976).

###  3.2 Numerical evaluations

Now, let us apply the second-order explicit symplectic method S2 or the
fourth-order method S4 to work out the magnetized RN-dS spacetime (11). We
take the parameters $E=0.975$, $L=4.2$, $\beta=8\times 10^{-3}$, $Q=0.3$,
$Q^{*}=1\times 10^{-6}$, and $\Lambda=2\times 10^{-6}$. A test orbit has
initial conditions $r=25$, $p_{r}=0$ and $\theta=\pi/2$. The starting value of
$p_{\theta}>0$ is determined by Eq. (14). Proper time step $h=1$ is adopted.
It can be seen clearly from Fig. 1 (a) that algorithm S2 performs bounded
Hamiltonian errors $\Delta K=-1/2-K$ with an order of $\mathcal{O}(10^{-6})$
in a long numerical integration of $10^{9}$ steps. The errors of S4 are
smaller in 3 orders than those of S2, but have a slightly secular drift due to
roundoff errors. Fortunately, this secular error growth is absent when a
larger proper time step $h=4$ is used in the fourth-order method S4*. The
fourth-order algorithm with the larger time step $h=4$ is the same as the
second-order method with the smaller time step $h=1$ in the order of
Hamiltonian errors. The former efficiency is 4/3 times faster than the latter
one.

In fact, the test orbit is a regular orbit colored red in Fig. 1 (b). This
regular orbit seems to consist of 7 loops on Poincaré section at the plane
$\theta=\pi/2$ and $p_{\theta}>0$ (note that the 7 loops are unclearly
visible). The order of Hamiltonian errors for each algorithm is not altered
when the red regular orbit is replaced with the green regular orbit, blue
regular orbit, or black weak chaotic orbit in Fig. 1 (b). In other words, the
algorithms’ accuracies are independent of dynamical behavior of orbits. In
what follows, we use the second-order method S2 with the smaller time step
$h=1$ to trace the dynamics of order and chaos of orbits.

###  3.3 Dynamical transition with a variation of the cosmological constant

To show cosmological constant $\Lambda$ how to affect the orbital dynamics of
the system (11), we consider $\Lambda$ to increase. The red orbit is regular
for $\Lambda=2\times 10^{-6}$ in Fig. 1 (b), but chaotic for
$\Lambda=2.1\times 10^{-6}$ in Fig. 2 (a). The extent of chaos in Figs. 2 (b)
and (c) is further strengthened with an increase of positive cosmological
constant $\Lambda$. In fact, strong chaos occurs for $\Lambda=1\times
10^{-5}$. However, no chaos exists for $\Lambda=0$ in Fig. 2 (d). In
particular, the red orbit is clearly composed of 7 small islands and is
approximately a periodic orbit in this case. These results seem to show that
chaos gets stronger as a positive cosmological constant increases. For
negative cosmological constants, such as $\Lambda=-5\times 10^{-6}$ and
$\Lambda=-1\times 10^{-5}$ in Figs. 2 (e) and (f), chaos is absent, either.

The magnetic parameter $\beta=8\times 10^{-3}$ is taken in Figs. 1 and 2. What
about the dynamics of the system (11) when $\beta$ gets slightly larger? For
$\beta=9\times 10^{-3}$ with $\Lambda=0$ in Fig. 3 (a), two strong chaotic
orbits exist, but do not for $\beta=8\times 10^{-3}$ with $\Lambda=0$ in Fig.
2 (d). As the positive cosmological constant further increases, chaos is
easier to occur. For $\Lambda=2\times 10^{-6}$ in Fig. 3 (b), three chaotic
orbits are present. For $\Lambda=1\times 10^{-5}$ in Fig. 3 (c), the four
orbits are chaotic. However, an increase of the magnitude of negative
cosmological constant weakens the extent of chaos. For example, two weak
chaotic orbits exist for $\Lambda=-1.2\times 10^{-6}$ in Fig. 3 (e). When
$\Lambda=-1.5\times 10^{-6}$ in Fig. 3 (f), no orbits can be chaotic.

The aforementioned demonstrations are based on the choice of a small value of
the Coulomb parameter $Q^{*}$. Now, let us focus on the orbital dynamical
transition when larger values of $Q^{*}$ are considered. Comparisons between
Figs. 4 (a) and (b), Figs. 4 (c) and (d), Figs. 4 (e) and (f), Figs. 5 (g) and
(h), and Figs. 5 (i) and (j) still consistently support the result that chaos
becomes easier with an increase of positive cosmological constant $\Lambda$.
In particular, this result is irrespective of whether $Q$ and $Q^{*}$ are
larger or smaller, and $Q^{*}$ is positive or negative. However, an increase
of the magnitude of negative cosmological constant leads to suppressing the
occurrence of chaos or weakening the extent of chaos, as shown by the
comparison between Figs. 5 (k) and (l). The chaos weakened is qualitatively
observed from the Poincaré section. The blue orbit is chaotic for
$\Lambda=-2\times 10^{-6}$ in panel (k), but it is not for $\Lambda=-1\times
10^{-5}$ in panel (l). In addition, other results are visible in Figs. 4 and
5. Figs. 1 (b), 2 (c), and 4 (a) and (b) describe that chaos is gradually
weakened as a positive Coulomb parameter $Q^{*}$ increases. And chaos is much
easily induced when the magnitude of negative Coulomb parameter increases. See
Figs. 4 (e) and (f), and 5 (g) and (h) for more details. These results are
consistent with those of Paper II (Wang et al. 2021b). As another one of the
main results in Paper II, an increase of $Q$ does not exert a typical
influence on the occurrence of chaos but twists the shape of orbits (such as
the green orbit). This result is also confirmed by the comparisons among Figs.
4 (a)-(d), and 5 (g)-(j).

It can be concluded from Figs. 1 (b), 2, 3, 4 and 5 that an increase of
positive cosmological constant (or an increase of the magnitude of negative
Coulomb parameter $Q^{*}$) leads to strengthening the extent of chaos from the
global phase space. However, an increase of the magnitude of negative
cosmological constant (or an increase of positive Coulomb parameter $Q^{*}$)
does not. These results can be explained similarly with the aid of slightly
modified Equation (36) in Paper II:

$\displaystyle K_{1}$ $\displaystyle\approx$ $\displaystyle-\frac{1}{2}(\beta
L+E^{2})-\frac{\Lambda}{6}Q^{*2}-\frac{1}{6}\Lambda E^{2}r^{2}$ (31)
$\displaystyle+\frac{1}{3}\Lambda
EQ^{*}r-\frac{E^{2}}{r}+\frac{\beta^{2}}{8}r^{2}\sin^{2}\theta$
$\displaystyle+\frac{EQ^{*}}{r}+\frac{Q^{2}E^{2}}{2r^{2}}+\frac{L^{2}}{2r^{2}\sin^{2}\theta}$
$\displaystyle+\frac{Q^{*}}{2r^{2}}(4E-Q^{*})+\cdots.$

The third and fourth terms in the potential (31) correspond to the forces
associated with the cosmological constant. The third term denominates the
fourth term. It acts as a repulsive force for the RN-dS black hole. The
repulsive force reduces the attractive force given by the potential $-E^{2}/r$
of the black hole. However, the third term yields an attractive force for the
RN-AdS black hole. The attractive force enhances the attractive force from the
black hole. On the other hand, the potential $EQ^{*}/r$ corresponds to the
Coulomb force from the Coulomb part of the magnetic field potential. For
$Q^{*}>0$, the Coulomb force is a repulsive force, which weakens the Lorentz
force as an attractive force determined by the potential
$\beta^{2}r^{2}\sin^{2}\theta/8$ of the magnetic field. However, the Coulomb
force is an attractive force for $Q^{*}<0$. This leads to enhancing the
attractive force from the magnetic field. Clearly, the Coulomb force from the
Coulomb part of the magnetic field potential is larger than the repulsive
force governed by the potential $Q^{2}E^{2}/(2r^{2})$ from the black hole’s
charge. Only when the gravitational attractions of charged particles from the
black hole and the magnetic field are approximately balanced, may chaos occur.
Therefore, chaos is easily caused in some circumstances when any one of the
positive cosmological constant, the magnitude of negative Coulomb parameter,
and the magnetic field parameter increases.

##  4 Conclusion

The Hamiltonian for the description of charged particles moving around the
RN-(A)dS black hole with an external magnetic field can be split into six
parts, which have analytical solutions as explicit functions of proper time.
In this case, second- and fourth-order explicit symplectic integration methods
are easily designed for this Hamiltonian system.

It is shown via numerical simulations that the newly proposed algorithms with
appropriate choices of step-sizes perform good long-term performance in
maintaining the boundness of Hamiltonian errors. The result should always be
the same, regardless of whether a test orbit is regular or chaotic.

The cosmological constant in the RN-dS black hole acts as a repulsive force,
which weakens the gravitational force of charged particle from the black hole.
As a result, the extent of chaos is strengthened from the global phase space
when a positive cosmological constant increases. In other words, the
accelerated expansion of the universe can easily induce chaos of charged
particles in some cases. However, chaos is weakened when the magnitude of
negative cosmological constant increases. This is because the cosmological
constant in the RN-AdS black hole acts as an attractive force, which causes
the gravitational effect from the black hole to dominate the attractive force
from the magnetic field. For the presence of chaos in the present metric
background, the effect of the cosmological constant on the accretion rate or
the gravitational waves will be worth considering.

## Acknowledgments

The authors are very grateful to a referee for useful suggestions. This
research has been supported by the National Natural Science Foundation of
China [Grant Nos. 11973020 (C0035736), 11533004, 11803020, 41807437,
U2031145], and the Natural Science Foundation of Guangxi (Grant Nos.
2018GXNSFGA281007 and 2019JJD110006).

## References

  * Abbott et al. (2020a) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020a, ApJL, 900, L13
  * Abbott et al. (2020b) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020b, Phys. Rev. Lett., 125, 101102
  * Ade et al. (2016) Ade, P. A. R., Aghanim, N., Arnaud, M., et al. 2016, A$\&$A, 594, A13
  * Ahmed et al. (2016) Ahmed, A. K., Camci, U., $\&$ Jamil, M. 2016, Class. Quant. Grav., 33, 215012
  * Ashtekar (2017) Ashtekar, A. 2017, Rep. Prog. Phys., 80, 102901
  * Ashtekar et al. (2016) Ashtekar, A., Bonga, B., $\&$ Kesavan, A. 2016, Phys. Rev. Lett., 116, 051101
  * Bacchini et al. (2018a) Bacchini, F., Ripperda, B., Chen, A. Y., $\&$ Sironi, L. 2018a, ApJS, 237, 6
  * Bacchini et al. (2018b) Bacchini, F., Ripperda, B., Chen, A. Y., $\&$ Sironi, L. 2018b, ApJS, 240, 40
  * Brink et al. (2015a) Brink, J., Geyer, M., $\&$ Hinderer, T. 2015a, Phys. Rev. D 91, 083001
  * Brink et al. (2015b) Brink, J., Geyer, M., $\&$ Hinderer, T. 2015b, Phys. Rev. Lett. 114, 081102
  * Brown (2006) Brown, J. D. 2006, Phys. Rev. D, 73, 024001
  * Carter (1968) Carter, B. 1968, Phy. Rev., 174, 1559
  * Deng et al. (2020) Deng, C., Wu, X., $\&$ Liang, E. 2020, MNRAS, 496, 2946
  * Dubeibe et al. (2007) Dubeibe, F. L., Pachón, L. A., $\&$ Sanabria-Gómez, J. D. 2007, Phys. Rev. D, 75, 023008
  * Eatough et al. (2013) Eatough, R. P., Falcke, H., Karuppusamy, R., et al. 2013, Nature , 501, 391
  * EHT et al. (2019a) EHT Collaboration, et al. 2019a, ApJL, 875, L1 (Paper I)
  * EHT et al. (2019b) EHT Collaboration et al. 2019b, ApJL, 875, L4 (Paper IV)
  * EHT et al. (2019c) EHT Collaboration et al. 2019c, ApJL, 875, L6 (Paper VI)
  * Ernst (1976) Ernst, F. J. 1976, J. Math. Phys. 17, 54
  * Ficek (2015) Ficek, F. 2015, Class. Quantum Grav., 32, 235008
  * Felice & Sorge (2003) Felice, D. d, $\&$ Sorge, F. 2003, Class. Quantum Grav., 20, 469
  * Forest & Ruth (1990) Forest, E., $\&$ Ruth, R. D. 1990, Physica D, 43, 105
  * Fukushima (2003a) Fukushima, T. 2003a, AJ, 126, 1097
  * Fukushima (2003b) Fukushima, T. 2003b, AJ, 126, 2567
  * Fukushima (2003c) Fukushima, T. 2003c, AJ, 126, 3138
  * Fukushima (2004) Fukushima, T. 2004, AJ, 127, 3638
  * Gurzadyan (2019) Gurzadyan, V. G. 2019, Eur. Phys. J. Plus, 134, 14
  * Gurzadyan & Stepanian (2019) Gurzadyan, V. G., $\&$ Stepanian, A. 2019, Eur. Phys. J. C, 79, 169
  * Han (2008) Han, W. B. 2008, Phys. Rev. D, 77, 123007
  * Hawking & Page (1983) Hawking, S. W., $\&$ Page, D. N. 1983, Commun. Math. Phys., 87, 577
  * Hu et al. (2019) Hu S., Wu X., Huang G., $\&$ Liang E. 2019, ApJ, 887, 191
  * Hu et al. (2021) Hu S., Wu X., $\&$ Liang E. 2021, ApJS, accepted
  * Jamil (2009) Jamil, M. 2009, Eur. Phys. J. C, 62, 609
  * Jamil et al. (2008) Jamil, M, Rashid, M. A., $\&$ Qadir, A. 2008, Eur. Phys. J. C, 58, 325
  * Karkowski1 & Malec (2014) Karkowski1, J., $\&$ Malec, E. 2013, Phys. Rev. D, 87, 044007
  * Karas & Vokroulflický (1992) Karas, V., $\&$ Vokroulflický, D. 1992, Gen. Relativ. Gravit. 24, 729
  * Kološ et al. (2017) Kološ, M., Tursunov, A., $\&$ Stuchlík, Z. 2017, EPJC, 77, 860
  * Kopáček et al. (2010) Kopáček, O., Karas, V., Kovář, J., $\&$ Stuchlík, Z. 2010, ApJ, 722, 1240
  * Kopáček & Karas (2014) Kopáček, O., $\&$ Karas, V. 2014, ApJ, 787, 117
  * Kopáček & Karas (2018) Kopáček, O., $\&$ Karas, V. 2018, ApJ, 853, 53
  * Kottler (1918) Kottler, F. 1918, Ann. Phys., 56, 401
  * Li & Wu (2019) Li, D., $\&$ Wu, X. 2019, Eur. Phys. J. Plus, 134, 96
  * Lubich et al. (2010) Lubich, C., Walther, B., $\&$ Brügmann, B. 2010, Phys. Rev. D, 81, 104025
  * Lukes-Gerakopoulos et al. (2016) Lukes-Gerakopoulos, G., Katsanikas, M., Patsis, P. A., $\&$ Seyrich, J. 2016, Phys. Rev. D, 94, 024024
  * Ma et al. (2008) Ma, D. Z., Wu, X., $\&$ Zhu, J. F. 2008, New Astrom., 13, 216
  * Maldacena (1998) Maldacena, J. 1998, Adv. Theor. Math. Phys., 2, 231
  * Mei et al. (2013a) Mei, L., Ju, M., Wu, X., $\&$ Liu, S. 2013a, Mon. Not. R. Astron. Soc., 435, 2246
  * Mei et al. (2013b) Mei, L., Wu, X., $\&$ Liu, F. 2013b, Eur. Phys. J. C, 73, 2413
  * Nakamura & Ishizuka (1993) Nakamura, Y., $\&$ Ishizuka, T. 1993, Astrophys. Space Sci. 210, 105
  * Nordström (1918) Nordström, G. 1918, Proc. Kon. Ned. Akad. Wet., 20, 1238
  * Pánis et al. (2019) Pánis, R., Kološ, M., Stuchlík, Z. 2019, Eur. Phys. J. C, 79, 479
  * Perlmutter (1999) Perlmutter, S., et al. 1999, ApJ, 517, 565
  * Preto & Saha (2009) Preto, M., $\&$ Saha, P. 2009, ApJ, 703, 1743
  * Riess et al. (1998) Riess, A. G., et al. 1998, AJ, 116, 1009
  * Reissner (1916) Reissner, H. 1916, Ann. Phys., 50, 106
  * Ruth (1983) Ruth, R. D, 1983, IEEE Trans. Nucl. Sci. NS 30, 2669
  * Jonathan & Lukes-Gerakopoulos (2012) Jonathan Seyrich, J., $\&$ Lukes-Gerakopoulos, G. 2012, Phys. Rev. D, 86, 124013
  * Schwarzschild (1916) Schwarzschild, K. 1916, Stizber. Deut. Akad. Wiss., Berlin, K1. Math.-Phys. Tech. s., 189
  * Shafiq (2020) Shafiq, S., Hussain, S., Ozair, M., Aslam, A., $\&$ Hussain, T. 2020, Eur. Phys. J. C, 80, 744
  * Shahzad & Jawad (2019) Shahzad, M. U., $\&$ Jawad, A. 2019, Can. J. Phys., 97, 742
  * Stuchlík (2005) Stuchlík, Z. 2005, Modern Physics Letters A, 20, 561
  * Stuchlík & Kološ (2016) Stuchlík, Z., $\&$ Kološ, M. 2016, Eur. Phys. J. C, 76, 32
  * Stuchlík et al. (2020) Stuchlík, Z., Kološ, M., Kovář, J., $\&$ Tursunov, A. 2020, Universe, 6, 26
  * Stuchlík et al. (2013) Stuchlík, Z., Kotrlová, A., $\&$ Török, G. 2013, A$\&$A, 552, A10
  * Swope et al. (1982) Swope, W. C., Andersen, H. C., Berens, P. H., $\&$ Wilson, K. R. 1982, J. Chem. Phys. 76, 637
  * Takahashi & Koyama (2009) Takahashi, M., $\&$ Koyama, H. 2009, ApJ, 693, 472
  * Tsang et al. (2015) Tsang, D., Galley, C. R., Stein, L. C., $\&$ Turner, A. 2015, ApJL, 809, L9
  * Tursunov1 et al. (2018) Tursunov1, A., Kološ, M., Stuchlík, Z., $\&$ Gal’tsov, D. V. 2018, ApJ, 861, 2
  * Wang et al. (2018) Wang, S. C., Huang, G. Q., $\&$ Wu, X. 2018, AJ, 155, 67
  * Wang & Jiang (2020) Wang, X., $\&$ Jiang, J. 2020, JCAP, 07, 052
  * Wang et al. (2021a) Wang Y., Sun W., Liu F., Wu X. 2021a, ApJ, 907, 66 (Paper I)
  * Wang et al. (2021b) Wang Y., Sun W., Liu F., Wu X., 2021b, ApJ, 909, 22 (Paper II)
  * Wang et al. (2016) Wang, S. C., Wu, X., $\&$ Liu, F. Y. 2016, MNRAS, 463, 1352
  * Wisdom (1982) Wisdom, J. 1982, AJ, 87, 577
  * Wisdom & Holman (1982) Wisdom, J., $\&$ Holman, M. 1991, AJ, 102, 1528
  * Xu & Wang (2017) Xu. Z., $\&$ Wang, J. 2017, Phys. Rev. D, 95, 064015
  * Yi & Wu (2020) Yi, M., $\&$ Wu, X. 2020, Phys. Scr., Phys. Scr., 95, 085008
  * Yoshida (1990) Yoshida, H. 1990, Phys. Lett. A, 150, 262
  * Zelenka et al. (2019) Zelenka, O., Lukes-Gerakopoulos, G., $\&$ Witzany, V. 2019, arXiv: 1903. 00360 [gr-qc]
  * Zelenka et al. (2020) Zelenka, O., Lukes-Gerakopoulos, G., Witzany, V., $\&$ Kopáček, O. 2020, Phys. Rev. D, 101, 024037
  * Zhong et al. (2010) Zhong, S. Y., Wu, X., Liu, S. Q., $\&$ Deng, X. F. 2010, Phys. Rev. D, 82, 124040

Figure 1: (a) Hamiltonian errors $\Delta K=-1/2-K$ for the new explicit
symplectic algorithms solving the magnetized RN-dS spacetime (11). The
parameters $E=0.975$, $L=4.2$, $\Lambda=2\times 10^{-6}$, $\beta=8\times
10^{-3}$, $Q=0.3$ and $Q^{*}=1\times 10^{-6}$ are used, and the initial
conditions of a test orbit are $r=25$, $p_{r}=0$ and $\theta=\pi/2$ with the
starting value of $p_{\theta}>0$ determined by Eq. (14) are chosen. The
second-order method S2 and the fourth-order method S4 take proper time step
$h=1$. The realistic errors for S4 are 1000 times smaller than the plotted
errors. The errors for S2 are stabilized in an order of
$\mathcal{O}(10^{-6})$, whereas do not remain bounded for S4 due to roundoff
errors. Fortunately, such a secular drift in Hamiltonian errors lose when a
larger time step $h=4$ is used in method S4*. (b) Poincaré sections on the
plane $\theta=\pi/2$ and $p_{\theta}>0$, given by algorithm S2 with proper
time step $h=1$. The test orbit in panel (a) is a regular orbit colored red in
panel (b). The other three orbits are regular tori colored green with initial
separation $r=10$ and blue with initial separation $r=18$, and a weak chaotic
orbit colored black with initial separation $r=30$.

Figure 2: Same as Fig. 1 (b) but for different values of cosmological constant
$\Lambda$.

Figure 3: Same as Fig. 2 but magnetic parameter $\beta=9\times 10^{-3}$ is
fixed.

Figure 4: Same as Fig. 2 but different combinations of parameters $Q$, $Q^{*}$
and $\Lambda$ are given.

Figure 5: Fig. 4 continued.