# Three-component Stäckel model of the Galaxy based
on the rotation curve from maser data

A. O.Gromov, I.I. Nikiforov e-mail<EMAIL_ADDRESS>
(St. Petersburg State University, St. Petersburg, 198504 Russia)

###### Abstract

A three-component Stäckel model of the Galaxy, including the bulge, disk, and
halo, is constructed. Parameter estimates of the potential are obtained as a
result of fitting the model rotation curve to azimuthal velocities found from
data on trigonometric parallaxes and spatial velocities of masers. The fitting
method takes into account the measurement and natural dispersions of azimuthal
velocities and uses an algorithm for excluding objects with excessive
residuals. In order to obtain more uniform samples, the objects were divided
into two groups: masers associated with high-mass star forming regions and
masers of other types. A significant kinematic inhomogeneity of these groups
was identified and taken into account: the azimuthal velocity dispersion is
$\sigma_{0,1}=4.3\pm 0.4$ km s-1, in the first group and $\sigma_{0,2}=15.2\pm
1.3$ km s-1 in the second. After constructing the model of the Galactic-plane
potential, it was generalized to the entire space under the assumption of the
existence of a third quadratic integral of motion. When reconstructing the
Galactic rotation curve in detail, the used algorithm gives an analytical
expression for the Stäckel potential, which significantly simplifies the task
of constructing the Galaxy’s phase density model in the Stäckel approximation.
In order to make the Stäckel model more realistic, one needs to develop
methods of direct account of data on the vertical distribution of density in
the Galaxy.

Keywords: methods: analytical — methods: data analysis — Galaxy: structure —
Galaxy: kinematics and dynamics

##  1 INTRODUCTION

Stäckel models are widely discussed in papers related in one way or another to
constructing phase models of star systems, including those of our Galaxy. This
is due to the possibility of the existence in such models of a third integral
of motion, quadratic in velocity (in addition to the two classical integrals
of energy and areas)

$I_{3}=(Rv_{z}-zv_{R})^{2}+z^{2}v_{\lambda}^{2}+z_{0}^{2}\left(v_{z}^{2}-2\Phi^{*}\right)\,,$
(1)

where

$\displaystyle z_{0}^{2}\frac{\partial\Phi^{*}}{\partial
R}=z^{2}\frac{\partial\Phi}{\partial R}-Rz\frac{\partial\Phi}{\partial
z}\,,\qquad$ $\displaystyle z_{0}^{2}\frac{\partial\Phi^{*}}{\partial
z}=\left(R^{2}+z_{0}^{2}\right)\frac{\partial\Phi}{\partial
z}-Rz\frac{\partial\Phi}{\partial R}\,.$

Here $R$, $\lambda$, $z$ are the cylindrical coordinates, $\Phi$ is the
potential of the model considered, and $z_{0}$ is some constant of the
dimension of length.

Such an integral is widely applied in stellar dynamics, for example, in the
works of Kuzmin (1952), Kuzmin and Malasidze (1987), Osipkov (1975) and many
other authors.

The velocity ellipsoid for most subsystems in the Galaxy, determined from the
proper motions and radial velocities of the stars, is triaxial (see, e.g., the
data on the principal velocity dispersions of main sequence stars with various
$B-V$ color indices in a monograph by Binney and Merrifield (1998)). One of
the ellipsoid axes coincides with the $v_{\lambda}$ axis. If the phase density
depends on three integrals of motion, then this fact can be explained within
the framework of the stationary Galaxy theory. Additionally, for models with
three integrals of motion, the problem of computing orbits can be solved in
quadratures. Since the radial velocity dispersion exceeds significantly the
vertical one not only for nearby stars but also for objects representing a
large interval of galactocentric distances (e.g., Bond et al., 2010;
Rastorguev et al., 2017), hereinafter we shall assume for the entire Galaxy
the existence of a global third integral as a phase density argument. Note
that the choice of parameter $z_{0}$ (see Section 4), the possible variations
of which represent the difference in the local integrals for individual stars,
can be considered as a selection of some average $z_{0}$ for the latter, which
even in this case allows us to talk about some common quasi integral of
motion.

For practical purposes, one should consider only the single-valued integrals
of motion: if non-isolating integrals become arguments of phase density, then
it becomes infinitely multivalued, which lacks physical sense (see, e.g., the
monograph by Ogorodnikov, 1958). Hénon and Heiles (1964) have shown
numerically that when energy reaches some critical value, ergodic areas are
formed where the third integral becomes non-isolating. Rodionov (1974)
determined that in order to avoid such «ergodic layers» a sixfold continuous
differentiability of the potential at the center of the system is needed.

The possibility of the existence of such third integrals has also been
considered for other galaxies (e.g., Binney et al., 1990; Merrifield, 1991).

The condition of the existence of the third quadratic integral of motion (1)
is written in the form of

$\displaystyle 3\left(z\frac{\partial\Phi}{\partial
R}-R\frac{\partial\Phi}{\partial
z}\right)-\left(R^{2}+z_{0}^{2}-z^{2}\right)\frac{\partial^{2}\Phi}{\partial
R\partial z}+Rz\left(\frac{\partial^{2}\Phi}{\partial
R^{2}}-\frac{\partial^{2}\Phi}{\partial z^{2}}\right)=0\,.$ (2)

For further computations it is convenient to switch to elliptical coordinates:
$\xi_{1}\in[1;\infty)$, $\xi_{2}\in[-1;1]$,

$R=z_{0}\sqrt{\left(\xi^{2}_{1}-1\right)\left(1-\xi^{2}_{2}\right)}\,,\quad
z=z_{0}\xi_{1}\xi_{2}.$ (3)

The limitation (2) can then be written as

$\frac{\partial^{2}}{\partial\xi_{1}\partial\xi_{2}}\left[\left(\xi_{1}^{2}-\xi_{2}^{2}\right)\Phi\right]=0\,,$
(4)

and the potential, in turn, as follows:

$\Phi=\frac{\varphi(\xi_{1})-\varphi(\xi_{2})}{\xi^{2}_{1}-\xi^{2}_{2}}\,,$
(5)

where $\varphi(\xi)$ is an arbitrary function.

Expression (4) is the condition of separating the variables, which means that
such potentials are separable, i.e., they allow one to separate the variables
in the Hamilton–Jacobi equation. Stäckel (1890) was the first to introduce in
his works a class of potentials that satisfy (5). They later became known as
Stäckel potentials and were introduced into stellar dynamics by Eddington
(1915).

Another advantage of Stäckel potentials in phase modeling is the ability to
describe models in angle–action variables. Multiple expressions for phase
density with such variables are encountered in the literature (e.g., Dehnen,
1999; Posti et al., 2015). The algorithm for determining the actions and
angles in separable potentials (Stäckel potentials are one type of such
potentials) was developed by J. Binney and his team. This approach replaced
the method of constructing tori (Kaasalainen and Binney, 1994), which turned
out to be inconvenient since it gives a dependence of the phase variables on
actions and angles and not vice versa, as is required, and the method of
adiabatic invariants by Binney (2010), which was satisfactory only for stars
close to the equatorial plane.

To introduce the Stäckel potential into models, Sanders (2012), Binney (2012),
Sanders and Binney (2016) use the «Stäckel fudge» algorithm, where the
$\varphi(\xi)$ function is an interpolation based on some number of orbital
points in the assumption that the non-Stäckel potential in those points has
the properties of the Stäckel one. However, this approach is approximate.
Additionally, the $z_{0}$ parameter changes from orbit to orbit in such an
approach. Due to this, one can only talk of some local third integral, since
the constancy requirement for $z_{0}$ is key to keeping integral $I_{3}$.
Which integral is the phase density argument in the case of a non-constant
$z_{0}$ is not quite clear.

Undoubtedly, studies within the framework of the Stäckel fudge have provided a
large contribution to the development of star system phase modeling methods.
However, since global integrals should, strictly speaking, be phase density
arguments, we suggest that an alternative approach based on the $I_{3}$
integral with constant $z_{0}$, as another variant of fitting a more complex
reality, is also of interest. In contrast to the numerical multi-step Stäckel
fudge algorithm, such an approach is analytical and more simple. This paper is
dedicated to the attempt to implement this analytical approach. The future
will show which method will prove to be more accurate and convenient.

The first step to construct a phase model is obtaining the model of the
potential that is consistent with observations. Unfortunately, there are only
a few studies where Stäckel potentials are constructed based on observations.
One of the first was the work of Satoh and Miyamoto (1976), where a one-
component model of the Galaxy was determined by the data on Galactic rotation
and by the density in the solar vicinity. However, that study uses only 18
objects with distances up to 10 kpc. Additionally, the paper by Satoh and
Miyamoto (1976) adopts the value $\rho_{\odot}=0.148$ $M_{\odot}\,pc^{-3}$,
which does not correspond to the contemporary results (Bland-Hawthorn and
Gerhard, 2016; Loktin and Marsakov, 2009).

In a later work, Famaey and Dejonghe (2003) constructed Stäckel models of the
Galaxy based on the derivatives, mainly local, of the dynamical
characteristics (local circular velocity, flat rotation curve, Oort constants
and others). Such models are only representative; in future, as Famaey and
Dejonghe (2003) themselves note, more extensive kinematic data should be used,
but such work has not been undertaken.

Formally, the models presented in Famaey and Dejonghe (2003) are three-
component, however, the third component only appears in them due to an
addition, mainly for demonstrative purposes, of a thick disk component to the
spheroid and (thin) disk, the dynamic contribution of which is small compared
to the thin disk (Bland-Hawthorn and Gerhard, 2016).

We should also note the paper by Binney and Wong (2017), where the «Stäckel
fudge» algorithm is applied to data on Galactic globular clusters and a phase
model is constructed for that cluster system. However, as was mentioned above,
using the «Stäckel fudge» algorithm gives only an approximation of the Stäckel
model.

This work, where we implement an analytical approach to constructing the
Stäckel model for the Galaxy, is the first in a series of papers which as the
end result aims to construct a phase model of the Galaxy which agrees with the
large amount of information obtained from observations. Based on kinematic
data on masers with trigonometric parallaxes we find the optimal parameters
for the model potential in the Galactic plane, which we then generalize to
three dimensions in the assumption of the existence of a third quadratic
integral of motion. As a result we construct a three-component (halo, disk,
bulge) Stäckel model of the Galaxy based on the current rotation data; we then
discuss the feasibility of this model.

##  2 METHOD

Stäckel models of the Galaxy constructed in the papers mentioned above are
merely approximate (representative or roughly Stäckel). The problem may be
solved by a method proposed by Rodionov (1974), which allows one to generalize
the equatorial plane potential to the entire space in the Stäckel way,
obtaining an analytical expression as a result. Since no limitations are
applied in this case to the equatorial plane potential, any model may be used
to construct the Stäckel model, for example, one obtained by approximating the
data on the rotation of the «cool» Galactic subsystem encompassing a large
interval of distances $R$. In particular, one can use maser data, and, as a
possibility, the Gaia catalog.

If the potential in set in the equatorial plane, $\varphi(\xi)$ is determined
as (Rodionov, 1974):

$\varphi(\xi)=\xi^{2}\;\Phi\left(R=z_{0}\sqrt{|\xi^{2}-1|},z=0\right)\,.$ (6)

Such a choice of function $\varphi(\xi)$ facilitates the fulfillment of the
following conditions:

  1. 1)

$\varphi(0)=0$;

  2. 2)

$\varphi(1)=\Phi_{0}$, where $\Phi_{0}$ — is the potential in the center of
the model;

  3. 3)

when $\xi\to\infty$
$\displaystyle\frac{\varphi(\xi)}{\xi}\to\frac{GM}{z_{0}}$, where $G$ — is the
gravitational constant and, $M$ — is the total mass.

The author of the method (Rodionov, 1974) also proposed ways to estimate the
function $\varphi(\xi)$, if the potential is defined on the $z$ and in an
arbitrary column $(R_{*},z)$, where $R_{*}=\text{const}$.

We used this method earlier to construct one- and two-component Stäckel test
models of the Galaxy (Gromov et al., 2015, 2016). We used the data on neutral
hydrogen rotation as well as those on masers. A good agreement between model
rotation curves and observed data was obtained, and also between model values
and estimates of density in the solar vicinity, mass in a sphere of 50 kpc
radius and other dynamic characteristics.

In this work we use the same method to construct a more realistic three-
component Stäckel model of the Galaxy with the classical composition: halo,
disk, and bulge. Since the goal is to determine the applicability of
Rodionov’s method in a separate representation of the main dynamic components
of the Galaxy in the model, further refinement of the model (for instance,
adding a thick disk) seems excessive at this stage.

We chose a quasi-isothermal potential proposed in Kuzmin et al. (1986) to
describe the halo:

$\Phi_{1}(R,0)=\Phi_{0,1}\ln\left(1+\frac{\beta}{w(R)}\right),$ (7)

where the $w(R)$ function is determined as:

$w^{2}(R)=1+\kappa_{1}^{2}\,R^{2}\,.$

The disk is presented as a general isochronous potential (Kuzmin and
Malasidze, 1969):

$\Phi_{2}(R,0)=\Phi_{0,2}\,\frac{\alpha}{\left(\alpha-1\right)+\sqrt{1+\kappa_{2}^{2}R^{2}}}\,.$
(8)

To describe the central bulge we used the Hernquist potential (Hernquist,
1990):

$\Phi_{3}(R,0)=\Phi_{0,3}\,\frac{1}{R+\kappa_{3}}\,.$ (9)

According to (6) the function $\varphi_{i}(\xi)$ for individual components has
the following form:

$\displaystyle\varphi_{1}(\xi)=\xi^{2}\Phi_{0,1}\ln\left(1+\frac{\beta}{\sqrt{1+\kappa_{1}^{2}z_{0}^{2}|\xi^{2}-1|}}\right),$
(10)
$\displaystyle\varphi_{2}(\xi)=\xi^{2}\Phi_{0,2}\,\frac{\alpha}{\left(\alpha-1\right)+\sqrt{1+\kappa_{2}^{2}z_{0}^{2}|\xi^{2}-1|}},$
(11)
$\displaystyle\varphi_{3}(\xi)=\xi^{2}\Phi_{0,3}\,\frac{1}{\sqrt{z_{0}^{2}|\xi^{2}-1|}+\kappa_{3}}.$
(12)

As a result, the final expression for the Stäckel potential function

$\varphi(\xi)=\varphi_{1}(\xi)+\varphi_{2}(\xi)+\varphi_{3}(\xi)$ (13)

is analytical, which allows us to find the function value and, as a
consequence, the values of the action variables in the entire region of
defined phase coordinates, and not in isolated points. The values of
$\varphi(\xi)$ are determined exactly here, and the accuracy of determining
the action variables that are expressed through complex integrals is limited
only by the accuracy of the numerical methods used to find the latter.

##  3 OBSERVED DATA

In this work we use data on masers located in star forming regions. VLBI
observations allow one to obtain accurate estimates of trigonometric
parallaxes and proper motions, including even for maser sources that are
rather distant from the Sun ($r<10$ kpc) (Nikiforov and Veselova, 2018).

We used the catalog from Rastorguev et al. (2017) as the main maser database,
which includes 103 high-mass star forming regions (HMSFRs) from the
homogeneous catalog of Reid et al. (2014) and 38 additional maser sources of
which, as determined by Nikiforov and Veselova (2018), nine belong to the
HMSFR class (see Table 1), and the remaining ones are of other types
(hereafter non-HMSFR masers). Furthermore, we added another 3 masers with full
data: Sh2-76EMM1, Sh2-76EMM2 (G040.44+02.45) and AFGL 5142 (G174.20-00.07);
they all belong to the non-HMSFR group. The total sample thus consists of 144
objects

Table 1: List of additional maser sources from the catalog of Rastorguev et al. (2017), which belong to high-mass starforming regions (HMSFRs) according to Nikiforov and Veselova (2018) G170.66${-}$00.25 | G059.83+00.67
---|---
G213.70${-}$12.60 | G071.52${-}$00.38
G054.10${-}$00.08 | G305.200+0.019
G058.77+00.64 | G305.202+0.208
G059.47${-}$00.18 |

As was shown by Nikiforov and Veselova (2018), the measurement uncertainty for
non-HMSFR maser parallaxes is on average higher than for the HMSFR masers.
Additionally, unlike the HMSFR sources belonging to the same class of objects,
non-HMSFR masers belong to 14 different classes (Nikiforov and Veselova,
2018). This fact allows us to assume that the sample uniting these objects is
kinematically inhomogeneous. We therefore differentiated between HMSFR (112
objects) and non-HMSFR (32 objects) masers when processing data in the general
case. In particular, dynamic modeling was carried out separately for the more
homogeneous HMSFR group and for the joint (HMSFR + non-HMSFR) sample. We did
not consider the non-HMSFR sample separately, since its volume is insufficient
for completing the outlined task. As we show below, HMSFR and non-HMSFR masers
are indeed kinematically inhomogeneous with respect to each other.

##  4 PARAMETER OPTIMIZATION

The parameters of the model potential described in Section 2 were estimated
using the nonlinear least squares method. The residual for an individual
object was the difference between the model circular velocity
$\theta_{C}(R_{i})$, where $R_{i}$ is the galactoaxial distance of the said
object, and azimuthal velocity $\theta_{i}$, computed by the measured position
and three velocity components of the object (see details in Gromov et al.
(2016)111In item 7 of the Appendix in Gromov et al. (2016) the correct formula
is $\sin\beta=\frac{r\cos b}{R}\sin l$.). Since the natural azimuthal velocity
dispersion of masers was considered as a-priori unknown, at the initial stage
we minimized the objective function

$L^{2}=\sum_{i=1}^{N}p_{i}\left[\theta_{i}-\theta_{C}(R_{i})\right]^{2},$ (14)

where

$\displaystyle p_{i}=\frac{1}{\sigma^{2}_{i}}$ (15)

is the weight of the $i$th object that accounts only for the measurement
uncertainty $\sigma_{i}$ of velocity $\theta_{i}$ (see Gromov et al. (2016));
$\displaystyle\theta_{C}^{2}=-R\,\frac{d\Phi}{dR}$,
$\Phi(R,0)=\Phi_{1}(R,0)+\Phi_{2}(R,0)+\Phi_{3}(R,0)$; $N$ is the volume of
the sample considered.

The parameter estimates for HMSFR masers and for those of the joint sample,
obtained by minimizing (14) are presented in Table 2, the corresponding model
rotation curves are fitted to the observed data in Fig 1. Since some
computation versions in this work have limiting parameter $\beta$ values —
$\beta\rightarrow\infty$ (in this case the quasi-isothermal potential becomes
the Schuster–Plummer potential),— the associated parameter
$q=\displaystyle\frac{\beta}{\beta+1}\in[0;1)\,$ was estimated instead.

Table 2: Solution without account for the natural azimuthal velocity dispersion Characteristic | HMSFRs | HMSFRs +
---|---|---
|  | non-HMSFRs
$N$ | 112 | 144
$q$ | $0.6748\pm 0.0073$ | $0.8421\pm 0.0022$
$\kappa_{1}$, kpc-1 | $0.0611\pm 0.0013$ | $0.2888\pm 0.0035$
$\Phi_{0,1}$, km2 s-2 | $305.5\pm 6.5$ | $262.8\pm 1.1$
$\alpha$ | $1.214\pm 0.023$ | $1.240\pm 0.019$
$\kappa_{2}$, kpc-1 | $0.533\pm 0.015$ | $0.10102\pm 0.00029$
$\Phi_{0,2}$, km2 s-2 | $400.6\pm 5.0$ | $310.83\pm 0.97$
$\Phi_{0,3}$, km2 s-2 | $230.4\pm 14.3$ | $233.7\pm 3.7$
$\kappa_{3}$, kpc | $0.25\pm 0.12$ | $0.70\pm 0.24$
$L^{2}/N_{\text{free}}$ | $236$ | $713$

Figure 1: Model circular velocity curves, obtained without account of the
natural dispersion of azimuthal velocities: from HMSFR masers (a), from a
joint HMSFR + non- HMSFR sample (b). The solid curve is the three-component
model of the Galaxy; the dashed curve shows the quasi-isothermal component
representing the halo; the dotted line is the general isochronous component
that describes the disk; the dashed-and-dotted curve is the Hernquist
component representing the bulge. Also shown are the maser azimuthal
velocities found from their measured parallaxes, radial velocities, and proper
motions. The circles show the HMSFR masers from the catalog of Reid et al.
(2014); rhombi show the HMSFR masers from the additional list in Rastorguev et
al. (2017); small squares show the non-HMSFR masers from the list of
Rastorguev et al. (2017) (they have little effect on the solution). The bars
show the initial data measurement errors. The absence of a bar means that it
is smaller than the size of the symbol representing the object. The horizontal
axis shows the distance of the Sun from the center of the Galaxy adopted in
this work, $R_{0}=8.34$ kpc Reid et al. (2014).

When minimizing (14) with weights (15) the values of $L^{2}/N_{\text{free}}$
(reduced statistics $\chi^{2}$), where $N_{\text{free}}=N-M$ is the number of
degrees of freedom, $M$ is the number of parameters to be determined, turned
out to be much larger than unity (Table 2). This implies that the deviations
of azimuthal velocities $\theta_{i}$ of masers from the model cannot be
explained only by measurement errors, and therefore, the natural velocity
dispersion has to be taken into account even for such the «cool» disk
subsystem as masers (in accordance with the results of Rastorguev et al.,
2017).

Therefore, at the second stage, the natural dispersion $\sigma_{0}$ of maser
azimuthal velocities was introduced into the objective function as an unknown
parameter by writing the weight coefficients as:

$p_{i}=\left(\sigma_{i}^{2}+\sigma_{0}^{2}\right)^{-1}\,,$ (16)

where $\sigma_{i}$ — is the measurement uncertainty $\theta_{i}$.

The value $\sigma_{0}$ was estimated in an iterative procedure, each step of
which had fixed parameters of the potential determined in the previous step,
and the natural dispersion was determined from the following equation:

$L^{2}\left(\sigma_{0}^{2}\right)=N_{\text{free}}\,.$ (17)

Using the $\sigma_{0}$ obtained at this stage, the parameters of the potential
were determined anew as a result of minimizing function (14) with the system
of weights (16). The found parameter estimates were used in the next step to
find a new approximation of $\sigma_{0}$. In all cases the procedure converged
in three steps. As an initial approximation we took the result obtained with
the system of weights (15).

To verify and account for the kinematic inhomogeneity of the non-HMSFR maser
group with respect to the HMSFR group, in addition to the natural dispersion,
common for the joint HMSFR + non-HMSFR sample (let us call this variant
«approach 1»), we also estimated this dispersion for the HMSFR group,
$\sigma_{0,1}$, and the non-HMSFR group, $\sigma_{0,2}$, individually (let us
call this variant «approach 2»). In the latter case we used the weight
coefficients $p_{i,1}=\left(\sigma_{i}^{2}+\sigma_{0,1}^{2}\right)^{-1}$ for
the HMSFR objects, and
$p_{i,2}=\left(\sigma_{i}^{2}+\sigma_{0,2}^{2}\right)^{-1}$ for non-HMSFR
objects.

In the case of two natural dispersions equation (17) can be written as:

$L^{2}\left(\sigma_{0,1}^{2},\sigma_{0,2}^{2}\right)=\sum_{i=1}^{N_{1}}\frac{\left[\theta_{i}-\theta_{C}(R_{i})\right]^{2}}{\sigma_{i}^{2}+\sigma_{0,1}^{2}}+\sum_{i=N_{1}+1}^{N}\frac{\left[\theta_{i}-\theta_{C}(R_{i})\right]^{2}}{\sigma_{i}^{2}+\sigma_{0,2}^{2}}=N_{\text{free}}\,,$

where the first sum is taken by the HMSFR objects, and the second — by the
non-HMSFR objects. If the dispersions in the denominators of the two sums are
correct, the contribution of each of them into the total result is
asymptotically equal to the fraction of each group of objects in the joint
sample. This gives the equations for determining $\sigma_{0,1}$ and
$\sigma_{0,2}$:

$\sum_{i=1}^{N_{1}}p_{i,1}\left[\theta_{i}-\theta_{C}(R_{i})\right]^{2}=\frac{N_{1}}{N}\,N_{\text{free}}\,,$
(18)
$\sum_{i=N_{1}+1}^{N}p_{i,2}\left[\theta_{i}-\theta_{C}(R_{i})\right]^{2}=\frac{N_{2}}{N}\,N_{\text{free}}\,,$
(19)

where $N_{1}$ and $N_{2}$ are the numbers of masers in HMSFR and non-HMSFR
groups respectively, $N=N_{1}+N_{2}$; for the total sample $N_{1}=112$,
$N_{2}=32$, $N=144$. The $\sigma_{0,1}$ and $\sigma_{0,2}$ values were
obtained as a result of an iterative procedure similar to the one described
above.

The results obtained at the second stage are presented in Table 3 and Fig. 2.

Table 3: Solution with account for the natural azimuthal velocity dispersions Characteristic | HMSFRs | HMSFRs + non-HMSFRs | HMSFRs + non-HMSFRs
---|---|---|---
|  | (approach 1) | (approach 2)
$N$ | 112 | 144 | 144
$\sigma_{0}$ or $\sigma_{0,1}$, km s-1 | $7.67\pm 0.74$ | $12.1\pm 1.0$ | $7.59\pm 0.65$
$\sigma_{0,2}$, km s-1 |  |  | $20.9\pm 1.8$
$q$ | $1_{-0.035}$ | $1_{-0.053}$ | $1_{-0.035}$
$\kappa_{1}$, kpc-1 | $0.0773\pm 0.0026$ | $0.0734\pm 0.0031$ | $0.0956\pm 0.0032$
$\Phi_{0,1}$, km2 s-2 | $265.2\pm 9.4$ | $264.9\pm 15.1$ | $265.0\pm 10.2$
$\alpha$ | $0.1437\pm 0.0041$ | $0.1148\pm 0.0048$ | $0.1651\pm 0.0062$
$\kappa_{2}$, kpc-1 | $0.0533\pm 0.0013$ | $0.04686\pm 0.0016$ | $0.0504\pm 0.0012$
$\Phi_{0,2}$, km2 s-2 | $311.2\pm 3.2$ | $318.4\pm 4.1$ | $317.4\pm 3.2$
$\Phi_{0,3}$, km2 s-2 | $224.5\pm 15.6$ | $223.4\pm 20.1$ | $223.8\pm 15.7$
$\kappa_{3}$, kpc | $0.71\pm 0.32$ | $0.72\pm 0.34$ | $0.76\pm 0.31$

Figure 2: Same as in Fig. 1, but with account for the natural dispersion of
azimuthal velocities: for the HMSFR sample (a), for the joint HMSFR + non-
HMSFR sample with the same (b), and with different (c) natural dispersions for
the HMSFR and non-HMSFR types. The top left parts of the panels show (in
vertical axis scale) the bars of root-meansquare deviations corresponding to
the derived natural dispersions.

The solution without account for the natural azimuthal velocity dispersion was
unstable with respect to the sample composition (see the results for two
samples in Table 2 and in Fig. 1), obviously, due to the exceedingly strong
influence of objects with a high ratio of the absolute residual value to the
(very small) measurement error (examples are evident in Fig. 1). Introducing
the natural dispersion has, for the most part, eliminated this effect and
stabilized the solution: all three results turned out to be close (Table 3,
Fig. 2). Only the natural dispersions differ, and significantly: adding non-
HMSFR masers to the sample increases $\sigma_{0}$ by a factor of 1.5; direct
dispersion measurements for the two maser groups show that $\sigma_{0,2}$ is
multiple times larger than $\sigma_{0,1}$. These facts imply a kinematic
inhomogeneity of the joint sample and justify the introduction of two
dispersions.

An analysis of the residuals of the derived solutions has shown that even
after taking into account the natural dispersion, the sample still contains
objects with data outliers, which can substantially shift the solution. To
eliminate objects with excessive residuals at the third and final stage, we
used the following algorithm with a flexible boundary for criterion statistics
(Nikiforov, 2012). We searched for objects for which

$\frac{|\theta_{i}-\theta_{C}(R_{i})|}{\sqrt{\sigma_{i}^{2}+\sigma_{0}^{2}}}>k,$
(20)

where $k$ is determined from the equation

$\left[1-\psi(k)\right]N_{\text{free}}=1,$

$\psi(z)=\displaystyle\sqrt{\frac{2}{\pi}}\int^{z}_{0}e^{-\displaystyle\frac{1}{2}\,t^{2}}dt$
is the probability integral and $\sigma_{0}$ is the corresponding natural
dispersion, i.e. $\sigma_{0}$, $\sigma_{0,1}$ or $\sigma_{0,2}$ depending on
the approach and object. Then out of $L$ objects that satisfied the condition
(20), $L-L^{\prime}$ object were excluded, where $L^{\prime}=3$ at the stage
of determining $\sigma_{0}$, $\sigma_{0,1}$, $\sigma_{0,2}$ or $L^{\prime}=1$
at the final stage of estimating model parameters with the derived (fixed)
natural dispersions (see Nikiforov, 2012). Objects with absolute residual
values exceeding $k_{0.05}$, i.e., the root of equation

$\left[1-\psi(k_{0.05})\right]N_{\text{free}}=0.05,$

were also excluded, if found among the remaining objects.

After each iteration of exclusions, the problems of parameter optimization,
natural dispersion determination, and finding possible outliers in the data
were solved anew until no object was excluded during the next iteration.

The object G110.19+02.47 with a negative residual $<3\sigma$, which is not
formally excessive according the criteria listed above, was also forcefully
excluded within the framework of approach 2. However, it is isolated from the
array of other non-HMSFR masers in the histogram of the distribution of
relative residuals. This is the only object the class of which could not be
determined from the original work of Chibueze et al. (2014). These authors
also note the low azimuthal velocity of this object, linking it with the local
kinematic anomaly, which does not characterize the Galaxy as a whole, of the
Perseus arm in the region of the rotation curve $R\sim 9$ kpc. Additionally,
including this object in the sample noticeably increases $\sigma_{0,2}$ from
$15.2\pm 1.3$ (final estimate) to $19.6\pm 1.7$ km/s. Maser G110.19+02.47 was
therefore excluded as a possibly anomalous object.

The list of excluded objects for the samples and approaches considered is
shown in Table 4. The final model parameter optimization results from the
maser data are presented in Table 5 and in Fig. 3. Model component masses
within a 50 kpc radius sphere are as follows: $M_{\text{bulge}}=1.1\times
10^{10}\leavevmode\nobreak\ M_{\odot}$, $M_{\text{disk}}=8.7\times
10^{10}\leavevmode\nobreak\ M_{\odot}$, $M_{\text{halo}}=7.2\times
10^{11}\leavevmode\nobreak\ M_{\odot}$.

Table 4: Lists of objects excluded from consideration in order of their exclusion during the final estimation of model parameters ($L^{\prime}=1$). The critical values of $k$ and $k_{0.05}$ are shown for the final sample volumes Object | $\rule{0.0pt}{16.99998pt}\displaystyle\frac{|\theta_{i}-\theta_{C}(R_{i})|}{\sqrt{\sigma_{i}^{2}+\sigma_{0,i}^{2}}}$
---|---
HMSFRs
($k=2.59$, $k_{0.05}=3.49$)
G348.70-01.04 | $5.53$
G023.44-00.18 | $4.01$
G213.70-12.60 | $4.12$
HMSFRs+non-HMSFRs
Approach 1
($k=2.67$, $k_{0.05}=3.55$)
G348.70-01.04 | $3.75$
G110.10+02.47 | $4.87$
G353.27+00.64 | $4.28$
G173.72-02.70 | $4.35$
G023.44-01.18 | $4.14$
G045.37-00.22 | $4.24$
G339.88-01.26 | $3.80$
Approach 2
($k=2.67$, $k_{0.05}=3.56$)
G348.70-01.04 | $5.83$
G110.10+02.47 | —
G353.27+00.64 | $3.06$
G213.70-12.60 | $4.00$
G023.44-01.18 | $4.11$
G045.37-00.22 | $4.13$
Table 5: Solution with account for the natural azimuthal velocity dispersion after excluding from consideration the masers with data outliers and the anomalous object Characteristic | HMSFRs | HMSFRs + non-HMSFRs | HMSFRs + non-HMSFRs
---|---|---|---
|  | (approach 1) | (approach 2)
$N$ | 109 | 137 | 138
$\sigma_{0}$ or $\sigma_{0,1}$, km s-1 | $3.85\pm 0.38$ | $6.56\pm 0.57$ | $4.34\pm 0.38$
$\sigma_{0,2}$, km s-1 |  |  | $15.2\pm 1.3$
$q$ | $1^{+0}_{-0.018}$ | $0.943\pm 0.023$ | $1^{+0}_{-0.014}$
$\kappa_{1}$, kpc-1 | $0.1173\pm 0.0023$ | $0.1012\pm 0.0025$ | $0.0830\pm 0.0020$
$\Phi_{0,1}$, km2 s-2 | $264.6\pm 6.0$ | $264.5\pm 7.2$ | $264.6\pm 5.3$
$\alpha$ | $0.2495\pm 0.0046$ | $0.1651\pm 0.0046$ | $0.1068\pm 0.0031$
$\kappa_{2}$, kpc-1 | $0.05455\pm 0.00085$ | $0.0575\pm 0.0011$ | $0.05645\pm 0.00055$
$\Phi_{0,2}$, km2 s-2 | $321.2\pm 1.9$ | $318.9\pm 3.1$ | $315.6\pm 1.7$
$\Phi_{0,3}$, km2 s-2 | $226.5\pm 9.7$ | $223.6\pm 13.4$ | $224.3\pm 7.4$
$\kappa_{3}$, kpc | $0.75\pm 0.12$ | $0.66\pm 0.25$ | $0.79\pm 0.12$

Figure 3: Same as in Fig. 2, but after excluding from consideration the masers
with data outliers and the anomalous object. The open symbols show the
discarded objects.

Parameter $z_{0}$ of the Stäckel potential was determined from the following
relation:

$\displaystyle z_{0}^{2}(R)=\frac{\displaystyle
3\,\frac{\partial\Phi}{\partial
R}(R,0)+R\left[\displaystyle\frac{\partial^{2}\Phi}{\partial
R^{2}}(R,0)-4\,\frac{\partial^{2}\Phi}{\partial
z^{2}}(R,0)\right]}{\displaystyle\frac{\partial^{3}\Phi}{\partial
z^{2}\partial R}(R,0)}-R^{2}\,.$ (21)

Similar functions for $z_{0}$ were constructed by Einasto and Rummel (1970),
and also Osipkov (1975). We used the three-component potential from Gardner et
al. (2011) as potential $\Phi(R,z)$ to find $z_{0}$ since data on the vertical
structure of the Galaxy were used in its construction. The derived value
$z_{0}=5.3$ kpc agrees with the estimates presented in other papers:
$z_{0}=2.2\div 7$ kpc (Hori, 1962; Kuzmin, 1953, 1956; Malasidze, 1973). We
should note that some of these studies use distances to the center of the
Galaxy $R_{0}$, that are lower than the current estimate. For example,
$R_{0}=7.0$ kpc in Kuzmin (1956) ($z_{0}=3.1$ kpc), $R_{0}=7.2$ kpc in Kuzmin
(1953) ($z_{0}=3.6$ kpc). Additionally, all the mentioned papers considered
spherical models, thus decreasing $z_{0}$. The value of $z_{0}$ derived by us
is close to $z_{0}=4.8$ kpc, reported by Malasidze (1973).

##  5 DISCUSSION

The options for the procedure of determining model parameters for the
potential of the Galaxy from masers considered in this work show that taking
into account even a small natural azimuthal velocity dispersion for masers is
necessary, given the modern high-accuracy data on these objects. The derived
azimuthal dispersion for the joint sample $\sigma_{0}=6.56\pm 0.57$ km s-1 is
intermediary between the radial ($\sigma_{U,0}=9.4\pm 0.9$ km s-1) and
vertical ($\sigma_{W,0}=5.9\pm 0.8$ km s-1) natural dispersions found by
Rastorguev et al. (2017) using almost the same sample, i.e. it agrees with
them. However the azimuthal dispersion for the homogeneous HMSFR sample is
significantly lower: close values were derived by two methods
$\sigma_{0,1}=3.9\pm 0.4$ and $4.3\pm 0.4$ km s-1 (Table 5). On the other
hand, the dispersion for non-HMSFR masers $\sigma_{0,2}=15.2\pm 1.3$ km s-1
(even after excluding the anomalous object) is about 3.5–3.9 times higher than
for HMSFRs, despite the expected general similarity of these types, which not
only confirms the kinematic inhomogeneity of non-HMSFRs compared to HMSFRs,
but implies a significant internal inhomogeneity of the non-HMSFR group. The
small volume of the non-HMSFR sample ($N=32$) and the large number of object
types there (14 types) do not allow us to determine the dispersion for
individual classes.

That said, introducing two natural dispersions for the joint sample lead to a
solution without significant shifts of point parameter estimates compared to
the results for the homogeneous HMSFR group, moreover, with smaller
statistical parameter errors on the whole. We therefore chose the model with
two natural dispersions as the final version for the Stäckel generalization,
derived after the exclusion of masers with excessive residuals and the
anomalous object (Table 5, approach 2; Fig. 3c).

The derived model rotation curves (Fig. 3) are close to models constructed
from masers at $R>3$–$4$ kpc in Rastorguev et al. (2017); Reid et al. (2019).
Our model reconstructs more accurately the decline of the rotation curve in
the outer part of the Galaxy than the one in Reid et al. (2019), since the
latter a priori used a smoother function form for the model. The model in
Rastorguev et al. (2017) reveals more details since they used a polynomial
fourth-order model for fitting. A good agreement is achieved in the common $R$
interval and with the rotation curves based on classical Cepheids in the
recent works of Ablimit et al. (2020); Mroz et al. (2019), which also show a
slow but steady decreasing velocity trend at, $R>6$ kpc. The same is also true
for the rotation curve from red giants (Eilers et al., 2019). On the whole,
modern rotation curves of the flat subsystems coincide rather well outside of
the central region of the Galaxy (see also Bland-Hawthorn and Gerhard (2016)).
The halo circular velocity curves derived from Gaia data on RR Lyrae stars
(Wegg et al., 2019) are close to our model, especially in the outer part of
the galaxy. At $R<3$ kpc our model cannot be reliable, if only due to the fact
that there it is based on only a few objects, to say nothing of the bar
influence (Chemin et al., 2015), but it is still sufficiently close to the
rotation curve constructed in the survey by Bland-Hawthorn and Gerhard (2016).

Dynamic modeling in this work was done in the assumption that the average
rotation velocity of masers is close to the circular velocity. This assumption
agrees with the low azimuthal velocity dispersion which we derived for masers,
especially for HMSFRs. Direct estimates of the asymmetric shift for masers are
represented by $\overline{V_{s}}$ in Reid et al. (2019): they are located in
the $-2\div-9$ km s -1 interval for different versions of analysis, but even
the most accurate of the estimates, $\overline{V_{s}}=-3.1\pm 2.2$ km s -1,
does not differ significantly from zero. Similar results were obtained also in
Reid et al. (2014). Thus, the asymmetric shift can be neglected for masers, at
least with the current level of accuracy.

We used Poisson’s equation to obtain an analytical expression for spatial
density for the constructed Stäckel potential [(formulas (5), (10)–(13)] The
corresponding equidensites (equal density curves) are presented in Fig. 4 for
a density of $\rho=0.1$ $M_{\odot}$ pc-3. As is evident from the figure, the
halo in the constructed model is non-spherical. The model also does not
adequately agree with the modern estimates of the vertical scale of the thin
($300\pm 50$ pc) and thick disk ($900\pm 180$ pc)(Bland-Hawthorn and Gerhard,
2016).

Figure 4: Model equidensites for $\rho=0.1$ $M_{\odot}$ pc-3. The solid curve
is the three-component model; the dashed curve shows the halo; the dotted
curve represents the disk, and the dashed-and-dotted curve is the bulge.

The derived non-too-realistic equidensites raise the question of the necessity
to additionally introduce into the model a vertical distribution of density in
one form or another when using Rodionov’s method. I.e., despite the
interesting results of Stäckel modeling in a series of theoretical studies
(e.g., Kuzmin (1952); Kuzmin and Malasidze (1987); Osipkov (1975)), it turns
out that the practical use of such methods in general does not lead to a
construction of models of acceptable adequacy. We should note that besides
some test work carried out by our group (Gromov et al. (2015, 2016)), there
were no earlier attempts to use the Stäckel generalization in practice
(applicable to real data). Such a result was therefore not obvious.
Considering the unrealistic equidensites obtained for disks in Binney and Wong
(2017); Famaey and Dejonghe (2003), where other methods of Stäckel modeling
were used (including the «Stäckel fudge»), this problem is of a more general
nature and not merely a shortcoming of the Stäckel generalization method.

Note that preliminary modeling based on the new catalog of Reid et al. (2019)
gives similar results, i.e. the nature of the vertical distribution is data
insensitive.

Thus, the vertical structure of the model obtained in this work implies that
using a direct Stäckel generalization method is insufficient for constructing
realistic Stäckel models. This would require solving the question of
accounting for the known properties of the vertical mass distribution in the
Galaxy in Stäckel modeling. We considered several approaches and their
combinations as possibilities: redefining the function $\varphi(\xi)$ with
account for the vertical distribution; introducing several additional
conditions when estimating parameters (for example, we were able to obtain an
adequate disk thickness in the assumption of a barometric vertical density
distribution in preliminary calculations); varying the parameter $z_{0}$ in
order to make the shape of the halo more spherical. We propose to study these
approaches in our next paper.

For the constructed model the solar neighborhood density ($R_{0}=8.34$ kpc
[32]) $\rho_{\odot}=0.082$ $M_{\odot}$ pc-3 is in good agreement with the
current estimates: $\rho_{\odot}=0.08$–$0.11$ $M_{\odot}$ pc-3(Loktin and
Marsakov, 2009), $\rho_{\odot}=0.097\pm 0.013$ $M_{\odot}$ pc-3(Bland-Hawthorn
and Gerhard, 2016).

The model mass in a sphere of 50 kpc radius is $M_{50}=0.82\times
10^{12}\leavevmode\nobreak\ M_{\odot}$, which is somewhat higher than in the
other studies. For example, Deason et al. (2012) found $M_{50}=(0.42\pm
0.04)\times 10^{12}\leavevmode\nobreak\ M_{\odot}$, and Williams and Evans
(2015) — $M_{50}=(0.45\pm 0.15)\times 10^{12}\leavevmode\nobreak\ M_{\odot}$.
The model mass in a 20 kpc sphere turned out equal to $0.35\times
10^{12}\leavevmode\nobreak\ M_{\odot}$. The overestimated masses are probably
also linked to omitting the data on the vertical density distribution.

##  6 CONCLUSIONS

We constructed an analytical three component Stäckel model of the Galaxy using
the data on masers with trigonometric parallaxes, proper motions, and radial
velocities.

The maser data were used to estimate the parameters of the components — halo,
disk, and bulge — for a model representing the potential in the Galactic
plane. We show that non-HMSFR masers, i.e., masers not related to high-mass
star forming regions (HMSFRs), are kinematically inhomogeneous with respect to
HMSFR masers: the natural (non-instrumental) dispersion of their azimuthal
velocity components, $\sigma_{0,2}=15.2\pm 1.3$ km s-1 in 3.5–3.9 times higher
than the similar dispersion $\sigma_{0,1}=(3.9\div 4.3)\pm 0.4$ km s-1 for
HMSFRs. This fact should be taken into account when using data on the
kinematics of inhomogeneous maser types. The final model was obtained with
account for the difference in natural dispersions for these two groups —
HMSFRs and non-HMSFRs. The model circular velocity curve is in good agreement
with the observed data.

The derived model of the Galactic plane potential served as a basis for
constructing a Stäckel model of the Galaxy using Rodionov’s method (Rodionov,
1974), which allows us to determine the values of the potential in the entire
space. The model gives close to real values of the density in the solar
neighborhood of $\rho_{\odot}=0.082$ $M_{\odot}$ pc-3 and the total mass of
$M_{50}=0.82\times 10^{12}\leavevmode\nobreak\ M_{\odot}$ (accurate to an
order of magnitude). However, the nature of equidensites in the meridian plane
of the model shows that direct use of the Stäckel generalization method alone
is not enough to construct realistic Stäckel models. To that end, one needs to
develop methods of directly accounting for data on the vertical density
distribution in the Galaxy when using Rodionov’s method.

Despite the uncovered but solvable complication, the approach on the whole has
significant potential: the proposed simple and, especially, analytical method
of computing function $\varphi(\xi)$ of the Stäckel potential can be used when
finding actions and may significantly simplify the corresponding algorithms of
constructing the phase model of the Galaxy.

The authors are grateful to V.V. Bobylev for providing the information on
three masers with full data and A.V. Veselova for the help with data
reduction, and especially for the complex task of identifying object types.
The authors are grateful to the anonymous referee for useful comments.

##  7 FUNDING

The reported study was funded by RFBR, project number 19-32-90144.

##  8 CONFLICT OF INTEREST

The authors declare no conflict of interest.

## References

  * 1. I. Ablimit, G. Zhao, C. Flynn, and S. A. Bird, Astrophys. J. 895 (1), L12 (2020).
  * 2. J. Binney, Monthly Notices Royal Astron. Soc. 401 (4), 2318 (2010).
  * 3. J. Binney, Monthly Notices Royal Astron. Soc. 426 (2), 1324 (2012).
  * 4. J. Binney and M. Merrifield, Galactic Astronomy (Princeton University Press, Princeton, NJ, 1998).
  * 5. J. Binney and L. K. Wong, Monthly Notices Royal Astron. Soc. 467 (2), 2446 (2017).
  * 6. J. Binney, R. L. Davies, and G. D. Illingworth, Astrophys. J. 361, 78 (1990).
  * 7. J. Bland-Hawthorn and O. Gerhard, Annual Rev. Astron. Astrophys. 54, 529 (2016).
  * 8. N. A. Bond, Z. Ivezic, B. Sesar, et al., Astrophys. J. 716 (1), 1 (2010).
  * 9. R. A. Burns, T. Handa, H. Imai, et al., Monthly Notices Royal Astron. Soc. 467 (2), 2367 (2017).
  * 10. L. Chemin, F. Renaud, and C. Soubiran, Astron. and Astrophys. 578, A14 (2015).
  * 11. J. O. Chibueze, H. Hamabata, T. Nagayama, et al., Monthly Notices Royal Astron. Soc. 466 (4), 4530 (2017).
  * 12. J. O. Chibueze, H. Sakanoue, T. Nagayama, et al., Publ. Astron. Soc. Japan 66 (6), 104 (2014).
  * 13. A. J. Deason, V. Belokurov, N. W. Evans, and J. An, Monthly Notices Royal Astron. Soc. 424 (1), L44 (2012).
  * 14. W. Dehnen, Astron. J. 118 (3), 1201 (1999).
  * 15. A. S. Eddington, Monthly Notices Royal Astron. Soc. 76, 37 (1915).
  * 16. A.-C. Eilers, D. W. Hogg, H.-W. Rix, and M. K. Ness, Astrophys. J. 871 (1), 120 (2019).
  * 17. J. Einasto and U. Rummel, Astrophysics 6 (2), 120 (1970).
  * 18. B. Famaey and H. Dejonghe, Monthly Notices Royal Astron. Soc. 340 (3), 752 (2003).
  * 19. E. Gardner, P. Nurmi, C. Flynn, and S. Mikkola, Monthly Notices Royal Astron. Soc. 411 (2), 947 (2011).
  * 20. A. O. Gromov, I. I. Nikiforov, and L. P. Ossipkov, Baltic Astronomy 24, 150 (2015).
  * 21. A. O. Gromov, I. I. Nikiforov, and L. P. Ossipkov, Baltic Astronomy 25, 53 (2016).
  * 22. M. Henon and C. Heiles, Astron. J. 69, 73 (1964).
  * 23. L. Hernquist, Astrophys. J. 356, 359 (1990).
  * 24. G. Hori, Publ. Astron. Soc. Japan 14, 353 (1962).
  * 25. M. Kaasalainen and J. Binney, Monthly Notices Royal Astron. Soc. 268, 1033 (1994).
  * 26. G. G. Kuzmin, Publ. Tartu Astrofiz. Obs. 32, 332 (1952).
  * 27. G. G. Kuzmin, Izvestiya Akad. Nauk ESSR 2, 3 (1953).
  * 28. G. G. Kuzmin, Astron. Zh. 33, 27 (1956).
  * 29. G. G. Kuzmin and G. A. Malasidze, Akademiia Nauk Gruzinskoj SSR 54, 565 (1969).
  * 30. G. G. Kuzmin and G. A. Malasidze, Publ. Tartu Astrofiz. Obs. 52, 48 (1987).
  * 31. G. G. Kuzmin, U. I. K. Veltmann, and P. L. Tenjes, Publ. Tartu Astrofiz. Obs. 51, 232 (1986).
  * 32. A. V. Loktin and Marsakov V. A., Lecscii po zvezdnoi astronomii (SFedU, Rostov-na-Dony, 2009).
  * 33. G. A. Malasidze, Proc. conf. on Dynamics of Galaxies and Star Clusters, pp. 93–98 (1973).
  * 34. M. R. Merrifield, Astron. J. 102, 1335 (1991).
  * 35. P. Mroz, A. Udalski, D. M. Skowron, et al., Astrophys. J. 870 (1), L10 (2019).
  * 36. I. I. Nikiforov, Astron. and Astrophys. Transact. 27 (3), 537 (2012).
  * 37. I. I. Nikiforov and A. V. Veselova, Astr. Lett. 44 (11), 699 (2018).
  * 38. K. F. Ogorodnikov, Dinamika zvezdnykh sistem (Gosydarstvennoe izdatel’stvo fiziko-matematichskoi literatury, Moskva, 1958) [in Russian].
  * 39. L. P. Osipkov, in Vestnik Leningradskogo Universiteta. Seriya 1. Matematika. Mekhanika. Astronomiia, pp. 151–158 (1975) [in Russian].
  * 40. L. Posti, J. Binney, C. Nipoti, and L. Ciotti, Monthly Notices Royal Astron. Soc. 447 (4), 3060 (2015).
  * 41. A. S. Rastorguev, N. D. Utkin, M. V. Zabolotskikh, et al., Astrophysical Bulletin 72 (2), 122 (2017).
  * 42. M. J. Reid, K. M. Menten, A. Brunthaler, et al., Astrophys. J. 885 (2), 131 (2019).
  * 43. V. I. Rodionov, in Vestnik Leningradskogo Universiteta. Seriya 1. Matematika. Mekhanika. Astronomiia, pp. 142–148 (1974) [in Russian].
  * 44. J. Sanders, Monthly Notices Royal Astron. Soc. 426 (1), 128 (2012).
  * 45. J. Sanders and J. Binney, Monthly Notices Royal Astron. Soc. 457 (2), 2107 (2016).
  * 46. M. J. Reid, K. M. Menten, A. Brunthaler, et al., Astrophys. J. 793 (2), 72 (2014).
  * 47. C. Satoh and M. Miyamoto, Publ. Astron. Soc. Japan 28 (4), 599 (1976).
  * 48. P. Stäckel, Math. Annal. 35, 91 (1890).
  * 49. C. Wegg, O. Gerhard, and M. Bieth, Monthly Notices Royal Astron. Soc. 485 (3), 3296 (2019).
  * 50. A. A. Williams and N. W. Evans, Monthly Notices Royal Astron. Soc. 454 (1), 698 (2015).


# Character of Doped Holes in Nd1-xSrxNiO2

Tharathep Plienbumrung  Institute for Functional Matter and Quantum
Technologies, University of Stuttgart, Pfaffenwaldring 57, D-70550 Stuttgart,
Germany Center for Integrated Quantum Science and Technology, University of
Stuttgart, Pfaffenwaldring 57, D-70550 Stuttgart, Germany Michael Schmid
Waseda Research Institute for Science and Engineering, Waseda University,
Okubo, Shinjuku, Tokyo, 169-8555, Japan Maria Daghofer  Institute for
Functional Matter and Quantum Technologies, University of Stuttgart,
Pfaffenwaldring 57, D-70550 Stuttgart, Germany Center for Integrated Quantum
Science and Technology, University of Stuttgart, Pfaffenwaldring 57, D-70550
Stuttgart, Germany Andrzej M. Oleś<EMAIL_ADDRESS>Max Planck Institute
for Solid State Research, Heisenbergstrasse 1, D-70569 Stuttgart, Germany
Institute of Theoretical Physics, Jagiellonian University, Profesora
Stanisława Łojasiewicza 11, PL-30348 Kraków, Poland

(29 July, 2021)

###### Abstract

We investigate charge distribution in the recently discovered high-$T_{c}$
superconductors, layered nickelates. With increasing value of charge-transfer
energy we observe the expected crossover from the cuprate to the local triplet
regime upon hole doping. We find that the $d-p$ Coulomb interaction $U_{dp}$
plays a role and makes Zhang-Rice singlets less favorable, while the amplitude
of local triplets is enhanced. By investigating the effective two-band model
with orbitals of $x^{2}-y^{2}$ and $s$ symmetries we show that
antiferromagnetic interactions dominate for electron doping. The screened
interactions for the $s$ band suggest the importance of rare-earth atoms in
superconducting nickelates.

##  I Introduction: Superconducting infinite-layered nickelates NdNiO2

The discovery of Bednorz and Muller Bed86  started intense search for novel
superconductors with high values of the critical temperature $T_{c}$. But in
spite of a tremendous effort in the theory, the mechanism responsible for the
pairing in cuprates is still unknown Kei15 . Yet, this is one of the
fundamental open problems in modern physics.

Perhaps less spectacular was the recent discovery of superconductivity in
infinite-layered NdNiO2 doped by Sr Li19  as the values of $T_{c}$ are ”only”
close to 15 K Li20 . Nevertheless, it gave a new impulse to the theory of
high-$T_{c}$ superconductivity at large. To some extent, the nickelate
superconductor family is rather similar to cuprate superconductors Bot20  as
once again two-dimensional (2D) planes of transition metal and oxygen ions
play a central role here Bia19 . With Ni+ ions one has again $d^{9}$
electronic configuration and similar lattice structure but the apical oxygens
are absent. By following the same analysis for a NiO2 layer as one performed
for a CuO2 layer, we expect Ni1+ to exhibit antiferromagnetic (AFM) order. One
finds that the charge-transfer gap in nickelate is larger than that of cuprate
Jia20 ; Tha21 . The doped holes reside on oxygen sites in cuprates forming the
Zhang-Rice singlet Zha88 . On the contrary, the doped holes will likely reside
on Ni sites in doped Nd1-xSrxNiO2.

##  II Charge-transfer model: NiO2 model revisited

We introduced the multiband $d-p$ Hamiltonian for a NiO2 plane Tha21  starting
from a 2D Ni4O8 $2\times 2$ cluster with periodic boundary conditions (PBCs).
The basis set includes four orbitals per NiO2 unit cell: two $e_{g}$ orbitals
$\\{3z^{2}-r^{2},x^{2}-y^{2}\\}\equiv\\{z,\bar{z}\\}$ at each Ni+ ion and one
bonding $2p_{\sigma}$ orbital (either $2p_{x}$ or $2py$) at each oxygen ion in
the 2D plane,

${\cal H}=H_{dp}+H_{pp}+H_{\rm diag}+H_{\rm int}^{d}+H_{\rm int}^{p}.$ (1)

Here the first two terms in the Hamiltonian (1) stand for the kinetic energy:
$H_{dp}$ includes the $d-p$ hybridization $\propto t_{pd}$ and $H_{pp}$
includes the interoxygen $p-p$ hopping $\propto t_{pp}$,

$\displaystyle H_{dp}$ $\displaystyle=$
$\displaystyle\sum_{\\{m\alpha;j\nu\\},\sigma}\left(t_{pd}\hat{d}^{\dagger}_{m\alpha\sigma}\hat{p}_{j\nu\sigma}+{\rm
H.c.}\right),$ (2) $\displaystyle H_{pp}$ $\displaystyle=$
$\displaystyle\sum_{\\{i\mu;j\nu\\},\sigma}\left(t_{pp}\hat{p}^{\dagger}_{i\mu\sigma}\hat{p}_{j\nu\sigma}+{\rm
H.c.}\right),$ (3)

where $\hat{d}_{m\alpha\sigma}^{\dagger}$ ($\hat{p}_{j\nu\sigma}^{\dagger}$)
is the creation operator of an electron at nickel site $m$ (oxygen site $i$)
in an orbital $\alpha\in\\{z,\bar{z}\\}$ ($\mu\in\\{x,y\\}$). Here $z$ and
$\bar{z}$ stand for $3z^{2}-r^{2}$ and $x^{2}-y^{2}$ orbitals, while
$\\{x,y\\}$ stand for $p_{x}$ and $p_{y}$ orbital. The elements
$\\{t_{pd},t_{pp}\\}$ are accompanied by the phase factors which follow from
orbital phases Dag09 .

The one-particle (level) energies are included in $H_{\rm diag}$, where we
introduce the charge-transfer energy between $d$ and $p$ orbitals,

$\Delta=\varepsilon_{d}-\varepsilon_{p}.$ (4)

The model is completed by the Coulomb interactions in the $d$ and $p$
orbitals. For the $d$ electrons,

$\displaystyle H_{\rm int}^{d}\\!$ $\displaystyle=$
$\displaystyle\sum_{m\alpha}U_{\alpha}n_{m\alpha\uparrow}n_{m\alpha\downarrow}+\left(U^{\prime}-\frac{1}{2}J_{H}\right)\sum_{i}n_{m1}n_{m2}$
(5) $\displaystyle-$
$\displaystyle\\!2J_{H}\\!\sum_{i}\\!\vec{S}_{m1}\\!\cdot\\!\vec{S}_{m2}+\\!J_{H}\\!\sum_{i}\\!d^{\dagger}_{m1\uparrow}d^{\dagger}_{m1\downarrow}d_{m2\uparrow}d_{m2\downarrow},$

where $\alpha=z,\bar{z}$ stands for $z^{2}$ and $x^{2}-y^{2}$ orbitals.
$d^{\dagger}_{m\alpha\sigma}$ ($d_{m\alpha\sigma}$) is electron creation
(annihilation) operator at site $m$ in orbital $\alpha$ with spin $\sigma$.
$n_{m\alpha\sigma}$ is a number operator. $\vec{S}_{m\alpha}$ stands for spin
operator in $\alpha$ orbital at site $m$. The two Kanamori parameters to
describe the interactions between $3d$ electrons are $\\{U,J_{H}\\}$.
$U_{\alpha}$ is Coulomb repulsion element for $\alpha$ orbital. $J_{H}$ and
$U^{\prime}=U-2J_{H}$ are Hund’s exchange and interorbital Coulomb interaction
Ole83 . Similar interactions with $\\{U^{p},J_{H}^{p}\\}$ are written for $p$
electrons.

Table 1: Parameters of the NiO2 charge-transfer model (all in eV) used in Tha21 . $t_{pd}$ | $t_{pp}$ | $\Delta$ | $U_{z}=U_{\bar{z}}$ | $J_{H}$ | $U_{p}$ | $J_{H}^{p}$
---|---|---|---|---|---|---
1.30 | 0.55 | $7.0$ | 8.4 | 1.2 | 4.4 | 0.8
Figure 1: Doping tendency of NiO2 model as a function of crystal-field
splitting $\Delta$: (a) $U_{dp}=0$; (b) $U_{dp}=1$ eV. (Blue) undoped,
$N_{\uparrow}=2$ and $N_{\downarrow}=2$; (Orange) one hole doped,
$N_{\uparrow}=3$, $N_{\downarrow}=2$; (Green) two holes doped,
$N_{\uparrow}=3$, $N_{\downarrow}=3$ in Ni4O8 (Cu4O8) cluster. The nickelate
(cuprate) regime is highlighted in blue (orange). The model parameters are
given in Table 1.

Figure 1 shows different weight distributions for hole occupations at Ni and O
sites for the parameters of Table 1. The Ni-O charge-transfer model has been
studied via impurity Jia20  as well as lattice approach Tha21 . It has been
established that the holes in undoped compounds remain within
$d_{x^{2}-y^{2}}$ orbitals. Thereby we have assumed that Nd does not
contribute to the electronic structure and the system without Sr is a Mott
insulator. Doping by Sr gives a doped hole which tends to reside at nickel
sites rather than at oxygen sites. The nickelate (cuprate) regime Jia20  is
highlighted in blue (orange) in Fig. 1. In the nickelate regime the holes
reside predominantly at Ni sites. This is the essential difference with
cuprates where a doped hole (for hole doping) resides predominantly at oxygen
and forms a Zhang-Rice singlet Zha88 .

Including intersite Coulomb repulsion $U_{dp}$ favors the hole occupancy at Ni
sites Tha21  and shifts the doping crossover to lower values of the charge
transfer energy $\Delta$. All the on-site energies of the Ni($3d$) orbitals
have been included in the Ni-O hybridization terms $t_{pd}$ Jia20 . Similar
results were obtained for finite $e_{g}$ orbital splitting, where
$\Delta_{z}=1$ eV should be considered the upper limit.

The NiO2 compound is a Mott-Hubbard insulator. It is then possible to replace
the charge-transfer model by the $d$-only Hubbard model. The effective Ni-Ni
hoppings can be derived from second-order perturbation theory Hu19 . The next
question is on which $d$-orbitals the doped holes are?

The asymmetric distribution of holes suggests that one could replace the Ni-O
model (1) with Ni $d$-only model as oxygen $p$-orbitals become unimportant.
The DFT calculation shows the band structure of NdNiO2 that two bands are
crossing Fermi level. In the orbital-resolved band structure, the lower band
has $d_{x^{2}-y^{2}}$ character and the upper band contains both Nd and Ni
contributions. The large charge transfer energy as well as the presence of
electron pocket at $\Gamma$ are the two striking features of the nickelate
compound. The empty 5d states of Nd are responsible for providing electron
pocket. The empty $5d$ states are below the Fermi level, in other words, these
states provide the hole states into Ni band by the so-called ’self-doped’
effect ZYZ20 . Furthermore, the $5d$ states was shown to be hybridized with Ni
apical orbitals i.e., $3d_{z^{2}}$ and 4$s$ orbitals Adh20  lead us to
construct the effective two-band model consisting of Ni in-plane orbital,
$d_{x^{2}-y^{2}}$, and Ni off-plane orbital, the modified $s$ orbital. In this
work, we present the character of doped holes in the realistic two-band model
of NdNiO2 compound.

##  III Electronic structure calculations

Figure 2: DFT band structure and DOS of the Wannier Hamiltonian: (a) DFT band
structure (black solid lines) compared with Wannier band structure (red dashed
lines), and (b) Wannier DOS for the Ni(${x^{2}-y^{2}}$) (red) and apical
Ni($s$) orbital (blue).

The ab initio electronic structure calculations using density functional
theory (DFT) were performed with the Quantum Espresso code QE1 ; QE2 ; QE3
using a plane-wave pseudopotential method, combining a projector augmented
wave method Lej16  and a specific choice of pseudopotentials pot . Within this
work we chose an energy cutoff of $600$ eV and a $\Gamma$-centered Brillouin
zone mesh of $16\times 16\times 16$. For all calculations we used the same
crystal structure published in Ref. Gu20 . The two-band model is derived by
performing a Wannier projection onto the DFT band structure as implemented
within the Wannier90 interface Wan90 , including onsite energies and hopping
parameters of each Wannier orbital. The projection was performed onto a
Ni(${x^{2}-y^{2}}$) and an apical Ni($s$) Wannier orbital within an energy
window ranging from $-2.0$ to $4.0$ eV. The DFT band structure, represented
here by the projected Wannier bands, and its density of states (DOS), are
shown in Fig. 2. Panel Fig. 2a shows the DFT band structure (black lines) and
Wannier bands (red dashed lines). The right panel Fig. 2b contains the
resulting Wannier DOS around the Fermi energy for the bands of two symmetries:
${x^{2}-y^{2}}$ and $s$. Our result qualitatively matches previous studies
Gu20 ; Adh20 .

Table 2: Hoppings parameters for the two-band model (all in eV) used in ED calculations. $(x,y,z)$ | $t^{\alpha\beta}_{ijk}$
---|---
(0,0,0) | $\begin{pmatrix}0.2&0\\\ 0&1.2\end{pmatrix}$
(1,0,0) | $\begin{pmatrix}-0.380&-0.050\\\ -0.050&-0.031\end{pmatrix}$
(0,1,0) | $\begin{pmatrix}-0.380&0.050\\\ 0.050&-0.031\end{pmatrix}$
(0,0,1) | $\begin{pmatrix}-0.039&0\\\ 0&-0.076\end{pmatrix}$
(1,1,0) | $\begin{pmatrix}0.088&0\\\ 0&-0.111\end{pmatrix}$
(1,0,1) | $\begin{pmatrix}0.001&-0.009\\\ -0.009&-0.252\end{pmatrix}$
(0,1,1) | $\begin{pmatrix}0.001&0.009\\\ 0.009&-0.252\end{pmatrix}$
(1,1,1) | $\begin{pmatrix}0.015&0\\\ 0&0.056\end{pmatrix}$

The Wannier Hamiltonian is given within the real space notation and is of the
form,

$\displaystyle\mathcal{H}(\mathbf{R})$
$\displaystyle=\sum_{ij\alpha\sigma}t_{ij\alpha}(\mathbf{R})c_{i\alpha\sigma}^{\dagger}c_{j\alpha\sigma}.$
(6)

The basis is ordered following the convention $\\{d_{x^{2}-y^{2}},s\\}$.
Explicit numerical values for the hopping parameters are given in Table 2 for
the terms with leading contributions, i.e., terms larger than $0.001$ (the
other terms were neglected). Note that the relation
$\mathcal{H}(-\mathbf{R})=\mathcal{H}^{T}(\mathbf{R})$ holds for Wannier
models. As a consequence, both terms (for distances $\mathbf{R}$ and
$-\mathbf{R}$) need to be included in further calculations. From
$\mathcal{H}(\mathbf{R})$ a tight-binding Hamiltonian can be constructed by
applying a Fourier transformation of the form,

$\displaystyle\mathcal{H}_{ij}(\mathbf{k})=\sum_{\mathbf{R}}\mathrm{e}^{\mathrm{i}\mathbf{k}\mathbf{R}}\mathcal{H}_{ij}(\mathbf{R}),$
(7)

where $\mathbf{R}$ describes the distances of the Wannier orbitals $|i-j|$ and
is typically represented in terms of the lattice vectors.

##  IV Effective two-band Hamiltonian

Following the idea of neglecting oxygen orbitals, we construct an effective
model containing only Ni($3d$) orbitals. The band structure calculation shows
that only two bands contribute at Fermi level. Therefore, the two-band model
of nickelates is capable of reproducing the physics of nickelates. It consists
of $x^{2}-y^{2}$ orbital and the $s$ orbital which includes rare-earth $5d$
states and Ni apical orbitals.

We consider the two-band Hamiltonian $H=H_{\rm kin}+H_{\rm int}$. The kinetic
part is

$H_{\rm
kin}=\sum_{i\alpha\sigma}\epsilon_{\alpha}a^{\dagger}_{i\alpha\sigma}a_{i\alpha\sigma}+\sum_{ij\alpha\beta\sigma}t_{ij}^{\alpha\beta}a^{\dagger}_{i\alpha\sigma}a_{j\beta\sigma}\,,$
(8)

while the interactions are given in a similar way to Eq. (5). The orbitals
${x^{2}-y^{2}}$ and $s$ have the diagonal energies $\epsilon_{\alpha}$, being
0 and $\epsilon$. The bands are constructed following Ref. Adh20 . The oxygen
$2p$ orbitals are included implicitly in ${x^{2}-y^{2}}$ ones. Note that one
should not confuse here $x^{2}-y^{2}$ symmetry with Ni($x^{2}-y^{2}$) orbital.
The former is an effective orbital including in-plane oxygen while the latter
is solely a Ni($3d$) orbital. The $s$ orbital contains of contributions from
Nd($5d$), Ni($3d_{z^{2}}$) and Ni($4s$) orbitals. The symmetry of this hybrid
orbital is the same as that of atomic $s$ orbital.

Instead of Nd($4f$) electronic states, the empty Nd($5d$) orbitals are
responsible for the striking electron pocket at the $\Gamma$ point in the band
structure Bot20 . The Nd atoms are originally located in the off-plane
direction of Ni-O plane. In this effective model the rare-earth atoms are
included into Ni atoms via $s$ orbital, indicating the model is three-
dimensional and the Coulomb interaction as well as Hund’s coupling need to be
screened. For instance, the two electrons sitting within $s$ orbital could be
located either at Nd or at Ni atom. Therefore, the Coulomb repulsion between
these two electrons within the $s$ orbital is reduced. We introduce a
parameter $\alpha\in[0,1]$ to represent the reduced Coulomb interaction and
Hund’s coupling, i.e., $U_{2}=\alpha U_{1}$ and $J=\alpha J_{H}$.
$\alpha=1(0)$ stands for weak (strong) screening effect from rare-earth atoms,
and we consider two parameter sets $A$ and $B$, given in Table 3. We study the
effective two-band model via Lanczos algorithm Koch  on a $2\times 2\times 2$
unit cell.

Table 3: Parameters of the two-band model (all in eV) used in exact
diagonalization calculations.
The reference energies for the two bands of $x^{2}-y^{2}$ and $s$ symmetry are 0 and $\epsilon$, respectively. set | $\epsilon$ | $t$ | $U_{1}$ | $J_{H}$
---|---|---|---|---
$A$ | 1.21 | 0.38 | 8.0 | 1.2
$B$ | 1.21 | 0.38 | 4.0 | 0.6
Figure 3: Orbital-resolved electron densities as obtained for: (a) set $A$,
and (b) set $B$, see Table 3. Strongly anisotropic electron distribution over
the $d_{x^{2}-y^{2}}$ and $s$ orbitals is favored when the interactions are
screened by $\alpha>0.5$, as shown for: (c) set $A$, and (d) set $B$.

The undoped nickelate correponds to quarter-filling, i.e., 8 electrons. The
stoichimetric compound has $d^{9}$ configuration, where Ni($x^{2}-y^{2}$)
orbital is half-filled, and the weak hybridization of Ni and Nd causes
Ni($x^{2}-y^{2}$) orbital to be away from half-filling and creates the self-
doping effect ZYZ20 ; Lec21 . The Ni-Ni hopping integrals are obtained via
fitting Wannier functions on DFT band calculation. We aim to address the
question: where the doped holes are located in the two-band model?

In the absence of electron hoppings Hamiltonian, the ground state is a trivial
antiferromagnet where $x^{2}-y^{2}$ orbital is half-filled and the two
canonical AFM phases, $C$-AFM and $G$-AFM, are degenerate. Including hopping
elements to further neighbors leads to metallic behavior. In one-band Hubbard
model the metal-insulator transition occurs when $U=2zt$, where $z=4$ is the
number of neighbors in the 2D plane. Similarly, the two-band description can
lead to partial orbital-selective Mott transition where one band is insulating
and the other one remains metallic. In the one-band version, the two relavant
configurations are singly occupyied or form double occupation within a single
site. The quarter filling two-band version, however, contains serveral
possible configurations where the lowest energy is still a single occupancy
followed by a local triplet state with energy $(U-3J_{H})$ Ole83 .

##  V Results and Discussion

We begin with electronic density distribution on $2\times 2\times 2$ clusters
obtained by exact diagonalization with twisted boundary condition (TBC). The
hopping parameters used here are given in Table 2. With PBC, one requires
hopping integrals $t^{\alpha\beta}_{ij}$ to be scaled by the factor of 1/2 due
to double-counting at the boundary. To avoid this additional factor in the
hopping terms, we replace the PBC with TBC. Instead of having a constant phase
$t_{N+i}=t_{i}$ as in PBCs, the hopping terms at the boundary are modified by
twisted angles $(\phi_{x},\phi_{y},\phi_{z})$, giving
$t_{N+i}=e^{i\vec{\phi}\cdot\vec{r}_{i}}t_{i}$; for more details see Refs.
Poi91 ; Shi97 . The PBC corresponds to $(0,0,0)$, while $(\pi,\pi,\pi)$ is
obtained for the anti-periodic boundary condition. In what follows, the
observables are obtained by averaging over several twisted angles.

Figure 4: Spin structure factor (9) of the three-dimensional two-band model,
as obtained for: (a) set $A$, and (b) set $B$, see Table 3.

In Figs. 3a-3b, we show the undoped orbital-resolved densities, i.e., the
electron occupancy within each orbital: $x^{2}-y^{2}$ (predominantly occupied)
and $s$ (it usually has decent amount of electrons). Finite occupancy of $s$
orbital is expected due to self-doping effect. We further see that increasing
of the parameter $\alpha$ causes the occupation on $x^{2}-y^{2}$ orbital to
approach half-filling. Once the $x^{2}-y^{2}$ orbitals are almost half-filled,
we show below that the $G$-AFM ground state is the only dominant magnetic
ground state. Including screening effect into $x^{2}-y^{2}$ orbital by
lowering its Coulomb interaction gives this orbital as slightly more
favorable, see also Fig. 3a-3b.

By adding one hole, we see that both orbitals can be occupied by holes
depending on how strong the screening effects of $s$ orbital are, see Figs.
3c-3d. The screening effect on $x^{2}-y^{2}$ orbitals makes little changes in
the one-hole case. An interesting scenerio arises when two holes are added
into the undoped Ni-O plane. When taking the screening effect into account,
holes are effectively occupying the $x^{2}-y^{2}$ orbital regardless of the
parameter $\alpha$.

Theoretical studies based on Ni($e_{g}$) bands Lec21 ; Tha21 ; ZYZ20  have
suggested that high-spin $S=1$ state Ni2+ is favorable for hole doping. In
contrast, RIXS measurements Ros20  show that two holes are residing mainly
within Ni($x^{2}-y^{2}$). Due to the limitation of the measurement, we cannot
determine the occupation of the rare-earth $5d$ states. The difference between
the predictions of the theoretical model and the experiment arises from the
rare-earth atoms. In the Ni($e_{g}$) with $d^{8}$ configuration, the
$x^{2}-y^{2}$ is half-filled and forming a high-spin state together with
another $3d$ orbital with energy $(U-3J_{H})$ Tha21 . This energy is
significantly smaller than the energy $U$ for a low-spin $d^{8}$ state. On the
contrary, the $s$ orbital, in the hole picture, has lower energy than
$x^{2}-y^{2}$ orbital and the quarter-filling of electrons corresponds to
$\frac{3}{4}$-filling by holes via particle-hole transformation. In this case,
the $s$ orbital is filled by holes and leaves half-filled $x^{2}-y^{2}$
orbital—then holes reside mainly on $x^{2}-y^{2}$ orbital. This intuitive
picture of hole configuration requires the $x^{2}-y^{2}$ to be nearly half-
filled, in other words, it suggests that the ground state is a strong
antiferromagnet as in cuprates Zha88 ; Ani99 . However, NMR experiments Hay03
; Hay99  report no observation of long-range AFM order in the $R$NiO2
($R$=Nd,La) down to 2 K.

Figure 5: Competition between low-spin and high-spin states as obtained for
the parameters of: (a,c) set $A$, and (b,d) set $B$, see Table 3. Panels (a)
and (b) show local double occupancy $D_{\alpha}$ in each
$\alpha=x^{2}-y^{2},s$ orbital, while panels (c) and (d) show local triplet
states $T$.

According to DFT+sicDMFT approach Lec21 , the paramagnetic ground state has
the lowest energy followed by $C$-AFM with energy difference about 20 meV/atom
and $G$-AFM with energy 105 meV/atom. The $C$-AF state has parallel spins
along the $c$ axis, while $G$-AFM phase has antiparallel spin alignment to its
nearest neighbor in all directions. To address the magnetic order in this two-
band model, the spin structure factor of undoped profile,

$\langle
S(k)S(-k)\rangle=\sum_{ij}e^{i\vec{k}(\vec{r}_{i}-\vec{r}_{j})}\left\langle\vec{S}_{i}\cdot\vec{S}_{j}\right\rangle,$
(9)

is shown in Fig. 4. With this small cluster size we can observe only six
relevant ${\bf k}$-points in the Brillouin zone. The spin structure factor at
$(\pi,\pi,\pi)$ is enhanced as the $x^{2}-y^{2}$ is close to half-filled band,
indicating that $G$-AFM order is favored, see Fig. 4. At small $\alpha$, the
spin structure factor (9) shows a competition between several magnetic ordered
states, dominated by $(\pi,\pi,0)$ and $(\pi,\pi,\pi)$. In the recent $x$-ray
scattering experiment Lu213 , the existence of AFM correlations was confirmed
but the small cluster size prevents us from concluding whether AFM long-range
order could be stable in this regime.

The screening effect from rare-earth atoms Sak20 ; ZYZ20  could be responsible
for the competition of magnetic orders. At small $\alpha$ where the screening
effect from rare-earth is strong, the electrons tend to form on-site triplets
with energy $\epsilon+\alpha(U_{1}-3J_{H})$ Tha21 , coexisting with singly
occupied $x^{2}-y^{2}$ orbitals. Then increasing $\alpha$ enhances the
interactions on $s$ orbitals and the on-site triplet energy surpasses the
bandwidth of $s$ orbital, i.e., $\epsilon+\alpha(U_{1}-3J_{H})>2zt$ ($z=4$).
The largest hopping elements are the in-plane hoppings along $x$ and $y$
direction of $x^{2}-y^{2}$ orbital. Another transition then occurs when the
on-site triplet energy overcome the bandwidth of $x^{2}-y^{2}$ orbital,
becoming a Mott insulator.

Figure 5 shows the competition between low-spin and high-spin states in NiO2
planes of the novel Ni-layered superconductors for decreasing number of holes,
i.e., for electron doping. Here we use
$D_{\alpha}=\frac{1}{N}\sum_{i}\left\langle
n_{i\alpha\uparrow}n_{i\alpha\downarrow}\right\rangle$ and
$T=\frac{1}{N}\sum_{i}\left\langle\frac{3}{4}+\vec{S}_{i1}\cdot\vec{S}_{i2}\right\rangle$,
where $N$ is the number of lattice sites. We show here that the double
occupancy $D_{\alpha}$ is reduced with decreasing number of holes in the
plane. Similarly, the amplitude of high-spin states $T$ at Ni sites is reduced
in this regime. So we conclude that electron doped materials have low-spin
configuration. On the contrary, hole doping may favor locally high-spin
$(S=1)$ states instead of singlet $(S=0)$ double occupancies of $x^{2}-y^{2}$
orbitals Tha21 .

##  VI Summary and conclusions

In summary, we have replaced the half-filled $d-p$ charge-transfer model by
the effective two-band model at quarter-filling. The characters of each band
are given by ${x^{2}-y^{2}}$ and $s$ symmetry. The contributions of Nd and Ni
atoms to $s$ orbital lead to the reduction of Hund’s exchange and Coulomb
repulsion. Its strength is scaled by the parameter $\alpha$. The model shows
the three competing phases: metal, orbital-selective and Mott insulator. The
Mott insulator is realized when the splitting between on-site triplet and
singly occupied state is larger than the size of the bandwidth, similar to
one-band Hubbard model. It is followed by the orbital-selective Mott insulator
where the orbital splitting separates the $s$ and $x^{2}-y^{2}$ bandwidth.

The nonmagnetic ordering is unlikely within this small unit cell and only AFM
configuration can be realized. While the $G$-AFM phase is clearly dominating
at large $\alpha$, the competition between $C$-AFM and $G$-AFM is found at low
$\alpha$ where the on-site triplet competes with AFM ground state, indicating
the tendency toward AFM ordering. The holes are doped differently among the
three phases. In metallic phase, the itinerant $s$ orbital is favorable for
holes. On the other hand, in the orbital-selective phase, the situation is
slightly complicated while adding one hole favors $s$ orbital but adding two
holes can effectively lead again to $x^{2}-y^{2}$ occupation, as in the Mott
insulating phase.

When the screening effect on $x^{2}-y^{2}$ orbital is included, the stability
of orbital-selective phase is enhanced and metallic phase vanishes. The two-
band model at quarter-filling, therefore, connects the two controversial
scenerios of which orbitals preferred by doped-holes via the parameter
$\alpha$. The parameter somehow represents the screening effect from rare-
earth orbital, showing the importance of rare-earth atoms for the electronic
structure of superconducting nickelate.

Authors’ Contributions: Tharathep Plienbumrung perfomed numerical analysis of
the two-band model. The two-band model was derived from the $d-p$ charge-
transfer model by Michael Schmid. All authors selected the relevant
information, analyzed the predictions of the theory versus experimental data,
developed the interpretation of the numerical results, and wrote the
manuscript.

Funding: T. P. acknowledges Development and Promotion of Science and
Technology Talents Project (DPST). A. M. Oleś kindly acknowledges Narodowe
Centrum Nauki (NCN, National Science Centre, Poland) Project No.
2016/23/B/ST3/00839.

Acknowledgments: We would like to thank Andres Greco, Peter Horsch, Krzysztof
Rościszewski, and George A. Sawatzky for many insightful discussions. A. M.
Oleś is grateful for the Alexander von Humboldt Foundation Fellowship
(Humboldt-Forschungspreis).

Conflicts of Interest: We declare no conflicts of interest.

## References

  * (1) Bednorz, J.G.; Muller, K.A. Possible high-$T_{c}$ superconductivity in the Ba-La-Cu-O system. Z. Phys. B 1986, 64, 189-193.
  * (2) Keimer, B.; Kivelson, S.A.; Norman, M.R.; Uchida, S.; Zaanen, J. From quantum matter to high-temperature superconductivity in copper oxides. Nature 2011, 518, 179-186.
  * (3) Li, D.; Lee, K.; Wang, B.Y.; Osada, M.; Crossley, S.; Lee, H.R.; Cui, Yi; Hikita, Y.; Hwang, H.Y. Superconductivity in an infinite-layer nickelate. Nature 2019, 572, 624-627.
  * (4) Li, D.; Wang, B.Y.; Lee, K.; Harvey, S.P.; Osada, M.; Goodge, B.H.; Kourkoutis, L.F.; Hwang, H.Y. Superconductiviting Dome in Nd1-xSrxNiO2 Infinite Layer Films. Phys. Rev. Lett. 2020, 125, 027001.
  * (5) Botana, A.S.; Norman, M.R. Simitarities and differences between LaNiO2 and CaCuO2 and implications for superconductivity. Phys. Rev. X 2020, 10, 011024.
  * (6) Jarlborg, T.; Bianconi, A. Multiple Electronic Components and Lifshitz Transitions by Oxygen Wires Formation in Layered Cuprates and Nickelates. Condensed Matter 2019, 4, 15.
  * (7) Jiang, M.; Berciu, M.; Sawatzky, G.A. Critical Nature of the Ni Spin State in Doped NdNiO2. Phys. Rev. Lett. 2020, 124, 207004.
  * (8) Plienbumrung, T.; Daghofer, M.; Oleś, A.M. Interplay between Zhang-Rice singlet and high-spin states in a model for doped NiO2 planes. Phys. Rev. B 2021, 103, 104513.
  * (9) Zhang, F.C.; Rice, T.M. Effective Hamiltonian for the superconducting Cu Oxides. Phys. Rev. B 1988, 37, 3759-3761.
  * (10) Arrigoni, E.; Aichhorn, M.; Daghofer, M.; Hanke, W. Phase diagram and single-particle spectrum of CuO2 high-$T_{c}$ layers: Variational cluster approach to the three-band Hubbard model. New J. Phys. 2009, 11, 055066.
  * (11) Oleś, A.M. Antiferromagnetism and correlation of electrons in transition metals. Phys. Rev. B 1983, 28, 327-339.
  * (12) Hu, L.-H.; Wu, C. Two-band model for magnetism and superconductivity in nickelates. Phys. Rev. Research 2019, 1, 032046.
  * (13) Zhang, G.-M.; Yang, Y.F.; Zhang, F.C. Self-doped Mott insulator for parent compounds of nickelate superconductors. Phys. Rev. B 2020, 101, 020501(R).
  * (14) Adhikary, P.; Bandyopadhyay, S.; Das, T.; Dasgupta, I.; Saha-Dasgupta, T. Orbital selective superconductivity in a two-band model of infinite-Layer nickelates. Phys. Rev. B 2020, 102, 100501(R).
  * (15) Giannozzi, O.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G.L.; Cococcioni, M.; Daboet, I.; et al. QUANTUM ESPRESSO: A modular and open-source software project for quantum simulations of materials. J. Phys.: Condensed Matter 2009, 21, 395502.
  * (16) Giannozzi, O.; Andreussi, O.; Brumme, T.; Bunau, O.; Buongiorno Nardelli, M.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Cococcioni, M.; et al. Advanced capabilities for materials modelling with Quantum ESPRESSO. J. Phys.: Condensed Matter 2017, 29, 465901.
  * (17) Giannozzi, O.; Baseggio, O.; Bonfà, P.; Brunato, D.; Car, R.; Carnimeo, I.; Cavazzoni, C.; de Cironcoli, S.; Delugas, P.; Ruffino, F.F.; et al. Quantum ESPRESSO toward the exascale. J. Chem. Phys. 2020, 152, 154105.
  * (18) Lejaeghere, K.; Bihlmayer, G.; Björkman, T.; Blaha, P.; Blügel, S.; Blum, V.; Caliste, D.; Castelli, I.E.; Clark, S.J.; Dal Corso, A.; et al. Reproducibility in density functional theory calculations of solids. Science 2016, 351, aad3000.
  * (19) Dal Corso, A. Pseudopotentials periodic table: From H to Pu. Comp. Mat. Science 2014, 95, 337-350.
  * (20) Gu, Y.; Zhu, S.; Wang, X.; Hu, J.; Chen, H. A substantial hybridization between correlated Ni-$d$ orbital and itinerant electrons in infinite-layer nickelates. Communications Phys. 2020, 3, 84.
  * (21) Pizzi, G.; Zhu. S.; Wang, X.; Hu, J.; Chen, H. Wannier90 as a community code: new features and applications. J. Phys.: Condensed Matter 2020, 32, 165902.
  * (22) Koch, E. The Lanczos Method, in: The LDA+DMFT approach to strongly correlated materials. edited by E. Pavarini, E. Koch, D. Vollhardt, and A. Lichtenstein. Forschungzentrum Jülich, 2011.
  * (23) Lechermann, F. Doping-dependent character and possible magnetic ordering of NdNiO2. Phys. Rev. Materials 2021, 5, 044803.
  * (24) Shiroshi, M.; Wadati, M. Integrable boundary conditions for the one-dimensional Hubbard model. J. Phys. Soc. Jpn. 2012, 66, 2288-2301.
  * (25) Poilblanc, D. Twisted boundary conditions in cluster calculations of the optical conductivity in the two-dimensional lattice models. Phys. Rev. B 1991, 44, 9562-9581.
  * (26) Rossi, M.; Lu, H.; Nag, A.; Li, D.; Osada, M.; Lee, K.; Wang, B.Y.; Agrestini, S.; Garcia-Fernandez, M.; Chuang, Y.-D.; et al. Orbital and Spin Character of Doped Carriers in Infinite-Layer Nickelates. ArXiv:2011.00595 2020.
  * (27) Anisimov, V. I.; Bukhvalov, D.; Rice, T. M. Electronic structure of possible nickelate analogs to the cuprates. Phys. Rev. B 1999, 59, 7901-7906.
  * (28) Hayward, M. A.; Green, M. A.; Rosseinsky, M. J.; Sloan, J. Sodium Hydride as a Powerful Reducing Agent for Topotactic Oxide Deinteraction: Synthesis and Characterization of the Nickel(I) Oxide LaNiO2. J. Am. Chem. Soc. 1999, 121, 8843-8854.
  * (29) Hayward, M. A.; Rosseinsky, M. J. Synthesis of the infinite layer Ni(I) phase NdNiO2+x by low temperature reduction of NdNiO3 with Sodium Hydride. Solid State Sciences 2003, 5, 839-850.
  * (30) Lu, H.; Rossi, M.; Nag, A.; Osada, M.; Li, D.F.; Lee, K.; Wang, B.Y.; Garcia-Fernandez, M.; Agrestini, S.; Shen, Z.X.; et al. Magnetic excitations in infinite-layer nickelates. Science 2021, 373, 213-216.
  * (31) Sakakibara, H.; Usui, H.; Suzuki, K.; Kotani, T.; Aoki, H.; Kuroki K. Model Construction and a Possibility of Cupratelike Pairing in a New $d^{9}$ Nickelate Superconductor (Nd,Sr)NiO2. Phys. Rev. Lett. 2020, 125, 077003.