# Gravitational Waves in $\mathbf{f(R)}$ Gravity Power Law Model

Dhruba Jyoti Gogoi1 and Umananda Dev Goswami1 $1.$ Department of Physics,
Dibrugarh University, Dibrugarh 786004, Assam, India

###### Abstract

We’ve investigated the different polarization modes of Gravitational Waves
(GWs) in $f(R)$ gravity power law model in de Sitter space. It is seen that
the massive scalar field polarization mode exists in this model. The mass of
the scalar field depends highly on the background curvature and the power term
$n$. However, we found that the model doesnot exhibit a massive scalar mode
for $n=2$ and instead it shows a breathing mode in addition to the tensor plus
and cross modes. Thus mass of the scalar field is found to vary with $n$
within the range $1\leq n\leq 2$.

Modified Gravity; Power Law Model; Gravitational Waves

###### pacs:

04.30.Tv, 04.50.Kd

##  I Introduction

Inspired by scientific curiosity, modification of Einstein’s gravity was at
first attempted by Weyl in 1919 weyl_neue_1919  and thereafter by Eddington in
1923 eddington_mathematical_1924  by including higher order curvature
invariants in gravitational action. It is seen that these modifications had no
experimental motivations and they were almost ignored up to 1962, when it was
for the first time realized that modification of gravitational action had some
merits. Gravity from General Relativity (GR) is not renormalizable and thus it
is not possible to quantize it according to conventional methods. In 1962,
Utiyama and DeWitt showed that renormalization at one loop demands that the
Einstein-Hilbert action be modified by higher order curvature terms
utiyama_renormalization_1962 . Later, a number of drawbacks of GR have been
observed including the inability to explain cosmic acceleration. To overcome
these drawbacks of GR a number of new theories have been proposed suggesting
the further modification of GR. Among them, $f(R)$ gravity is a prominent one
keeping itself standing against the challenges and problems till now
sotiriou_fr_2010 ; nojiri_2011 ; nojiri_2007 . Another important reason to
choose $f(R)$ gravity is its less complexity compared to the other modified
theories of gravity nojiri_2017 .

After around 100 years of Einstein’s prediction of Gravitational Waves (GWs),
the Laser Interferometer Gravitational Wave Observatory (LIGO) Scientific
Collaboration and Virgo Collaboration announced the detection of GWs for the
first time in September 14, 2015, which was named as GW150914
abbott_observation_2016 . This observation was later followed by four other
events, viz. GW151226 abbott_gw151226_2016 , GW170104 abbott_gw170104_2017 ,
GW170814 abbott_gw170814:_2017  and GW170817 abbott_gw170817:_2017 . These
events opened new directions in testing GR and modified theories of gravity.
Therefore the study of GWs in modified gravity will play a very important role
in modifying GR and on the predictions of GR. The experimental results
obtained from the sector of GWs can effectively constrain the modified gravity
models to a fair extent and can even contribute to the bottom-up approach in
forming or proposing modified gravity theories and models in an efficient way.

In the metric formalism of $f(R)$ gravity, the GWs have other modes of
polarizations besides the tensor modes of polarization found in GR. Production
and potential detection of GWs in $f(R)$ gravity were studied by C. Corda
using Starobinsky model in corda01  and using a model of the form
$f(R)=R+R^{-1}$ in corda02 . These studies showed that massive mode of GW
radiation corda03  could exist in $f(R)$ theory of gravity. In a recent study,
the propagating degrees of freedom of GWs in $f(R)$ gravity are found to be 3
myung_propagating_2016 , and in another recent study on the Starobinsky model
liang_polarizations_2017  showed that there exists a mixed state of massless
breathing mode and massive scalar mode of polarization besides the usual
tensor modes in the model. Thus the Newmann Penrose (NP) formalism cannot be
applied in that study due to the existence of massive scalar modes of GW
polarization. In our study, we’ll restrict ourselves to the power law model of
$f(R)$ gravity given by $f(R)=\alpha R^{n}$, where $\alpha$ is an arbitrary
constant and $n$ is a real number. This model has been pursued in
jaime_about_2013  and where it was concluded that the model can produce a late
time acceleration in the Universe or a transient matter domination era, but
not both. Same results were also obtained in the Refs. amendola_power-
laws_2007 ; amendola_are_2007 ; amendola_conditions_2007  previously. In the
Refs. capozziello_quintessence_2003 ; carloni_cosmological_2005 ;
goswami_cosmological_2013  this model was discussed for the different special
cases along with $n=-1$ and $n=3/2$, in which $n=3/2$ case is conformally
equivalent to Liouville field theory. In another paper faraoni_r_2011 , the
model was used for invoking the chameleon mechanism to work within the solar
system and in which it was shown that for this to happen the value of $n$ is
required to be very close to $1$. The result rendered from the work implied
that the model is a poor candidate for a realistic alternative to dark energy.
The cosmological dynamics of the model also has been studied in Ref.
carloni_cosmological_2005 . It is to be noted that although the model used in
our case is already ruled out in solar system tests zakharov , we have noticed
that in several studies faraoni_r_2011 ; capozz ; borka ; zakharov2  the model
has been used as a toy model. These studies have motivated us to use this
model for its simple form and more importantly because of the existence of
massless scalar mode of polarization of GWs even in de sitter space. Also we
think that a wide range of studies with this sort of new result will
definitely help to decisively rule out or accept the model from future
astrophysical/cosmological observations.

In this work, we’d like to test the model in the GW regime to find out the
possibility of its viability and the drawbacks. Another important point of
this model is that it is the only model in $f(R)$ theories in the metric
formalism which can give pure massless breathing mode as a special case. Here
we’ll study the model with the special case of $n=2$ along with the stability
and polarization modes.

The paper is organized as follows. In the next section the variation of mass
of the scalar field in de Sitter space with the power term $n$ is studied. In
the third section stability of the model is studied in de Sitter space. The
fourth section contains a study towards the equivalence of the theory with the
Scalar-Tensor theory and the required conditions for such equivalence. The
tensor and the scalar polarization modes are studied with the help of geodesic
deviation in the fifth section and with the help of Newman-Penrose formalism
in the sixth section. In section seven, we have discussed a possible way to
check the validity of the model experimentally. In the last section, we
conclude the paper with a very brief discussion of the results and the future
aspects of the model in such type of studies.

##  II Scalar and tensor fields from the model

In this study of GWs we use de Sitter vacuum space. The action of $f(R)$
gravity in this space is given by

$S=\dfrac{1}{2\kappa}\int d^{4}x\sqrt{-g}\,f(R),$ (1)

where $f(R)$ is a function of Ricci curvature $R$. The vacuum field equations
obtained from this action are given by

$f^{\prime}(R)R_{\mu\nu}-\dfrac{1}{2}f(R)g_{\mu\nu}-\nabla_{\mu}\nabla_{\nu}f^{\prime}(R)+g_{\mu\nu}\,\square
f^{\prime}(R)=0,$ (2)

where $\square=g^{\mu\nu}\,\nabla_{\mu}\nabla_{\nu}$ and the prime over $f(R)$
denotes the derivative with respect to $R$. Taking trace of vacuum field Eq.
(2), we obtain

$f^{\prime}(R)R+3\,\square f^{\prime}(R)-2f(R)=0.$ (3)

Now, we assume that there is a propagating GWs in the spacetime, which will
perturb the metric around its background value. Since the perturbation is
usually very small and hence if we consider the background metric as
$\bar{g}_{\mu\nu}$ then to the first order of perturbation value $h_{\mu\nu}$,
we may express the spacetime metric as

$g_{\mu\nu}=\bar{g}_{\mu\nu}+h_{\mu\nu},\;\;\mbox{where}\;|h_{\mu\nu}|<<|\bar{g}_{\mu\nu}|.$
(4)

In view of this perturbation, expanding the Ricci tensor and the Ricci scalar
upto the first order of $h_{\mu\nu}$, we may write

$\displaystyle R_{\mu\nu}$ $\displaystyle\simeq\bar{R}_{\mu\nu}+\delta
R_{\mu\nu}+O(h^{2})$
$\displaystyle=\bar{R}_{\mu\nu}-\dfrac{1}{2}(\nabla_{\mu}\nabla_{\nu}h-\nabla_{\mu}\nabla^{\lambda}h_{\lambda\nu}-\nabla_{\nu}\nabla^{\lambda}h_{\mu\lambda}$
$\displaystyle+\square h_{\mu\nu})+O(h^{2})$ (5)

and

$\displaystyle R$ $\displaystyle\simeq\bar{R}+\delta R+O(h^{2})$
$\displaystyle=\bar{R}-\square
h+\nabla^{\mu}\nabla^{\nu}h_{\mu\nu}-\bar{R}_{\mu\nu}h^{\mu\nu}+O(h^{2}),$ (6)

where $\bar{R}$ is the de Sitter curvature. Thus, due to this perturbation the
trace Eq. (3) can be rewritten as

$3f^{\prime\prime}(\bar{R})\square\delta
R+\left[f^{\prime\prime}(\bar{R})\bar{R}-f^{\prime}(\bar{R})\right]\delta
R=0.$ (7)

Using our model

$f(R)=\alpha R^{n}$ (8)

the Eq. (7) can be expressed in the following compact form:

$(\square-m^{2})\delta R=0,$ (9)

where

$m^{2}=\dfrac{(2-n)}{3(n-1)}\bar{R}.$ (10)

The Eq. (9) is the conventional form of the scalar field equation (Klein-
Gordon) with $\delta R$ as the scalar field and hence here, $m$ can be
identified as the mass of this scalar field. Moreover, since this scalar field
corresponding to $\delta R$ is due to ripple in spacetime, so the field can be
considered as the quantized massive gravitational field.

From the expression (10) we see that the mass term of the scalar field depends
highly on the exponent term $n$ and the background curvature $\bar{R}$. This
indicates that near the massive objects like neutron stars etc., where the
background curvatures are very large, the mass of the scalar form of
gravitational field will increase and consequently this will slow down the
propagation speed of the massive mode of GWs in such regions according to this
particular model. Further, for a massive scalar field or massive mode of GWs
$n$ should lie in between $1$ and $2$ beyond which it will result tachyonic
instabilities. When $n=1$ the model will give the Einstein’s case. When $n$
starts to increase from $1$ to $2$ mass of the field decreases monotonically
and finally, for $n=2$ the massive scalar mode of polarization will vanish.
i.e. for the value $n=2$ model will give a massless scalar field. The
variation of the scalar field mass with respect to $n$ for the unit value of
the background curvature is shown in the Fig.1. One interesting observation
from the mass term is that when $n$ approaches to $1$ the mass of the scalar
field increases very rapidly and for $n=1$ the mass becomes infinite. Thus for
$n=1$ the scalar mode vanishes and only the massless spin-2 modes propagate.

Figure 1: Variation of mass $m$ of the scalar field with respect to $n$ for
the unit value of the background curvature $\bar{R}$.

The tensor field Eq. (2) in vacuum for the power law model (8) takes the from:

$nR^{n-1}R_{\mu\nu}-\dfrac{1}{2}g_{\mu\nu}R^{n}-n\nabla_{\mu}\nabla_{\nu}R^{n-1}+n\,g_{\mu\nu}\square
R^{n-1}=0.$ (11)

Perturbing this field equation around the de Sitter curvature $\bar{R}$ (see
the Eq. (5)) we may write,

$\displaystyle n\bar{R}\bar{R}_{\mu\nu}+n(n-1)\bar{R}_{\mu\nu}\delta
R-\dfrac{1}{2}\bar{g}_{\mu\nu}\bar{R}^{2}-\dfrac{1}{2}n\,\bar{g}_{\mu\nu}\bar{R}\,\delta
R$ $\displaystyle+n\bar{R}\,\delta
R_{\mu\nu}-\,n(n-1)\nabla_{\mu}\nabla_{\nu}(\delta
R)+n(n-1)\bar{g}_{\mu\nu}\square\delta R=0.$ (12)

In the ideal case, i.e. when $R=\bar{R}$, the de Sitter curvature should
satisfy the above field equation. Thus for the ideal case the above equation
gives,

$n\bar{R}_{\mu\nu}=\dfrac{1}{2}\bar{g}_{\mu\nu}\bar{R}.$ (13)

Trace of this equation gives

$n\bar{R}=2\bar{R},$

which can be written as

$n\bar{R}^{2}-2\bar{R}^{2}=0.$ (14)

Again the trace of Eq. (II) is,

$n\bar{R}^{2}-2\bar{R}^{2}+n(n-2)\bar{R}\delta R+3n(n-1)\square\delta R=0.$

Using Eq. (14) in the above equation we find,

$(\square-m^{2})\delta R=0,$ (15)

where $m^{2}=\dfrac{(2-n)}{3(n-1)}\bar{R}.$ This equation is identical to the
Eq. (9). For $n=2$ above equation reduces to,

$\square\delta R=0.$ (16)

This is the equation of massless scalar field for the case of $n=2$. This is
an interesting result, which shows that if the de Sitter curvature has to
satisfy the field Eq. (II) then the model (8) has to take the value $n=2$,
consequently, the massive scalar mode of polarization vanishes and a pure
massless scalar mode of polarization (also known as breathing mode of
polarization) is obtained.

Thus the Eq. (16) can be treated as a stability point for the model in de
Sitter space which identically satisfies the stability condition discussed in
the following section without modifying the Hubble constant in de Sitter
space.

Again, since $n=2$ in the model (8) (we’ll refer this model with $n=2$ as pure
$R^{2}$ model) makes the scalar field massless and independent of the the
background curvature $\bar{R}$, so for simplicity we take this advantage to
choose $\bar{R}=0$ and perturb the tensor field equation for pure $R^{2}$
model (see the Eq. (2)) to give

$R\delta
R_{\mu\nu}-\dfrac{1}{4}R^{2}\bar{g}_{\mu\nu}-\partial_{\mu}\partial_{\nu}R+\bar{g}_{\mu\nu}\square
R=0.$ (17)

Now we define a parameter,

$h_{\mu\nu}=\bar{h}_{\mu\nu}-\dfrac{1}{2}\bar{g}_{\mu\nu}\bar{h}-\bar{g}_{\mu\nu}\delta
R,$ (18)

where taking the trace we find $\bar{h}$ as

$\bar{h}=-\,h-4\delta R,$ (19)

and substituting $\bar{h}$ in the definition (18), $\bar{h}_{\mu\nu}$ can be
expressed as

$\bar{h}_{\mu\nu}=h_{\mu\nu}-\dfrac{1}{2}\bar{g}_{\mu\nu}h-\bar{g}_{\mu\nu}\delta
R.$ (20)

Choosing transverse traceless gauge condition:
$\partial^{\mu}\bar{h}_{\mu\nu}=0$, $\bar{h}=0$ and using the definition (18)
in Eq. (5), we get

$\delta
R_{\mu\nu}=\dfrac{1}{2}\left[-\,\square\bar{h}_{\mu\nu}+2\partial_{\mu}\partial_{\nu}(\delta
R)+\bar{g}_{\mu\nu}\square\delta R\right].$ (21)

Finally, using this Eq. (21) and Eq. (16) in the Eq. (17), we may write

$\square\bar{h}_{\mu\nu}=0,$ (22)

which is the usual equation that gives the tensor modes, i.e. plus and cross
modes of polarization of GWs. The solution of this equation is

$\bar{h}_{\mu\nu}=e_{\mu\nu}\exp(iq_{\mu}x^{\mu})+c.c.,$ (23)

where $\bar{g}_{\mu\nu}q^{\mu}q^{\mu}=0$ and $q^{\mu}e_{\mu\nu}=0.$ Thus, this
clearly shows that the pure $R^{2}$ gravity model gives three polarization
modes of GWs, all massless, two of them are tensor modes while the other is
the scalar or the pure breathing mode as disccussed above.

##  III Stability of the Model

According to John D. Barrow and A. C. Ottewill barrow_stability_1983  the
condition for the existence of de Sitter solution for a constant curvature
$\bar{R}$ is

$\bar{R}f^{\prime}(\bar{R})=2f(\bar{R}).$ (24)

In Einstein’s case $\bar{R}=4\Lambda$ barrow_stability_1983 , where $\Lambda$
is the cosmological constant. For our model (8), we get from the above
condition that

$n\bar{R}^{n}=2\bar{R}^{n}.$ (25)

This condition is similar to the condition obtained in Eq. (13), giving $n=2$
and which transforms the Eq. (9) into the Eq. (16) with mass for the scalar
field $m=0$ as mentioned in the previous section. However, there is another
possible solution for this Eq. (18) is $\bar{R}=0$, which will lead to
Minkowski spacetime having zero background curvature.

In faraoni_stability_2005  it was shown that the model (8) produce power law
inflation $a\propto t^{\bar{\alpha}}$, with

$\bar{\alpha}=\dfrac{(n-1)(1-2n)}{(n-2)}$

for generic values of $n\neq 0,1/2,1$ and as we have already seen that this
model does not allow or admit de-Sitter solution for $n\neq 2$. This will be a
case only when a cosmological constant is added to the action corresponding to
a term with $n=0$ barrow_stability_1983 ; faraoni_stability_2005 . But for any
positive value of $n$ Minkowski space is a solution without any cosmological
constant faraoni_stability_2005 . The stability condition for this model from
faraoni_stability_2005  reads

$\bar{R}\dfrac{2-n}{n(n-1)}\geq 0.$ (26)

This condition gives, $1<n\leq 2$ and $n<0$. From Eq. (10) it can be easily
seen that the first stability condition i.e. $1<n\leq 2$ agrees well with the
non-tachyonic range of the scalar field. $n>2$ and $n<1$ in Eq. (10) will make
$m^{2}$ negative resulting tachyonic instabilities. However, the second
stability condition obtained from Eq. (26) i.e. $n<0$ in Eq. (10) can also
result tachyonic instabilities. To avoid such tachyonic instabilities in our
toy model we would not consider the situation satisfied by the later
condition. Therefore, the model is considered to be stable in the range
$1<n\leq 2$ only. For $n=2$ it is identically satisfied without imposing any
constraints on the Hubble parameter in de Sitter space and also shows that for
any values of $n$ the Minkowski space is stable faraoni_stability_2005 . Thus
our previous results are supported by these inferences.

##  IV Equivalence with the Scalar Tensor Theory

A straight forward way to study $f(R)$ gravity is to see its equivalence with
the Scalar Tensor Theory (STT) chiba_1/r_2003 . This is reasonable because as
like $f(R)$ theories, STTs also have two types of polarization modes, tensor
modes and scalar mode. By introducing a scalar field $\phi\equiv R$, the
general form of the action (1) can be given as

$S=\dfrac{1}{2\kappa}\int d^{4}x\sqrt{-g}\,\big{[}\psi(\phi)R-V(\phi)\big{]},$
(27)

where $\psi(\phi)=f^{\prime}(\phi)$ and $v(\phi)=\phi
f^{\prime}(\phi)-f(\phi).$ Now, differentiating Eq. (27) with respect to
$\phi$, we find

$f^{\prime\prime}(R)(R-\phi)=0.$ (28)

This shows that when $f^{\prime\prime}(R)\neq 0$, we must have $R=\phi$. Thus
$f^{\prime\prime}(R)\neq 0$ is a very important condition required to be
satisfied in order to compare the $f(R)$ theory with the STT. The pure $R^{2}$
model satisfies this condition. In metric formalism, the STT with Brans-Dicke
parameter $\omega=0$ is equivalent to the $f(R)$ gravity.

Again the action (27) can be written as

$S=\dfrac{1}{2\kappa}\int
d^{4}x\sqrt{-g}\big{[}f(\phi)+f^{\prime}(\phi)(R-\phi)\big{]}.$ (29)

The field equation obtained from the above action is given as

$\displaystyle
G_{\mu\nu}=\dfrac{1}{f^{\prime}(\phi)}\big{[}\nabla_{\mu}\nabla_{\nu}f^{\prime}(\phi)-g_{\mu\nu}\,\square
f^{\prime}(\phi)+\dfrac{1}{2}\,g_{\mu\nu}\\{f(\phi)$ $\displaystyle-\phi
f^{\prime}(\phi)\\}\big{]}.$ (30)

Trace of the above equation is

$\square f^{\prime}=\dfrac{2}{3}f(\phi)-\dfrac{1}{3}\phi f^{\prime}(\phi).$
(31)

For the pure $R^{2}$ model the Eq. (IV) gives

$G_{\mu\nu}=\phi^{\,-1}\big{(}\partial_{\mu}\partial_{\nu}\phi-
g_{\mu\nu}\,\square\,\phi-\dfrac{1}{4}g_{\mu\nu}\,\phi^{2}\big{)}.$ (32)

And from the trace of this equation we get,

$\square\,\phi=0.$ (33)

This equation is similar to the Eq. (16) which gives a massless breathing mode
of polarization. In Minkowski space $\square\,\phi\equiv\square\,\delta\phi$,
since $\square\,\bar{\phi}=0$ in this space, so the solution of this equation
can be written as,

$\delta\phi=\phi_{0}\exp(ip_{\mu}x^{\mu})+c.c.,$ (34)

where $\bar{g}_{\mu\nu}p^{\mu}p^{\nu}=0.$ In case of non zero background, a
damping factor can appear in the solution. But the effect of this damping
factor is very negligible and we can safely neglect the term damping_sol .
Hence, this equation is valid in both de Sitter spacetime and Minkowski
spacetime. Now, combining the solutions for the massless tensor modes and the
scalar mode, and considering that the GW is travelling along the z-axis with
the speed of light $c=1$, we get the solution for $h_{\mu\nu}$ as

$h_{\mu\nu}=\bar{h}_{\mu\nu}(t-z)-\bar{g}_{\mu\nu}\,\delta\phi(t-z).$ (35)

Hence, for this solution we may take $q_{\mu}=\Omega\,(1,0,0,1)$ and
$p_{\mu}=\omega\,(1,0,0,1)$. It should be noted that for $n\neq 2$ and $n\neq
1$ the above solution will contain a mixed state of massive scalar mode and
breathing mode besides the tensor modes of polarization
liang_polarizations_2017 .

##  V Geodesic Deviation and Polarization Modes

The de Sitter spacetime can be expressed as a hypersurface in the host pseudo-
Euclidean space with metric $\eta_{AB}=(+1,-1,-1,-1,-1)$ with points in
Cartesian coordinates $\chi^{A}$ satisfying
$\eta_{AB}\chi^{A}\chi^{B}=-\,l^{2}$, here $l$ is the de Sitter length
parameter. The four-dimensional stereographic coordinates ${x^{\mu}}$ are
obtained through a stereographic projection from the de Sitter hypersurface
into a target Minkowski spacetime geo00 . Taking into account the local
transitivity properities of spacetime, we consider a family of modified de
Sitter geodesics that are able to connect any two points in spacetime and they
are satisfied by the equation geo01 ; geo02

$\dfrac{dU_{\rho}}{ds}-\Lambda^{\gamma}_{\rho\sigma}U_{\gamma}u^{\sigma}=0,$
(36)

where $s$ is an affine parameter, $u^{\sigma}$ is the usual $4$-velocity of
particles and $\Lambda^{\gamma}_{\rho\sigma}$ is the Christoffel connection.
$U^{\mu}=\xi^{\mu}_{\alpha}u^{\alpha}$ is the anholonomic $4$-velocity, which
is related with the usual $4$-velocity of particles via the Killing vector
$\xi^{\mu}_{\alpha}$. The geodesic deviation equation obtained from the Eq.
(36) can be written as

$\dfrac{D^{2}V^{\mu}}{Ds^{2}}\delta\lambda=R^{\mu}_{\nu\sigma\rho}U^{\nu}u^{\sigma}\eta^{\sigma}+\dfrac{D}{Ds}\left[\left(\eta^{\alpha}u^{\beta}-\eta^{\beta}u^{\alpha}\right)\triangledown_{\beta}\xi^{\mu}_{\alpha}\right],$
(37)

where $\dfrac{D}{Ds}=u^{\sigma}\triangledown_{\sigma}$ represents the
covariant derivative with respect to $s$,
$V^{\mu}=\xi^{\mu}_{\alpha}v^{\alpha}$, $\eta^{\mu}=\dfrac{\partial
x^{\mu}}{\partial\lambda}\delta\lambda=v^{\mu}\delta\lambda$ and $\lambda$ is
some parameter such that for each $\lambda=$ constant,
$x^{\mu}=x^{\mu}(s,\lambda)$ will represent the equation of a geodesic in Eq.
(36). This equation shows that geodesic deviation in de Sitter spacetime
contains translational (first term on RHS) as well as rotational (second term
on RHS) terms geo03 . Since, the second term on RHS is model independent as it
doesnot contain Riemannian tensor term, the rotational effect is a generic
property of the spacetime only. That is why, in Newtonian limit with
$l\rightarrow\infty$ we can safely neglect the second term on RHS and thus can
have an approximated and simplified form of geodesic deviation equation with
only first term on the RHS. Now to see the geodesic deviation in our model, we
first substitute the solution (34) in the Eq. (32), which gives,

$R_{\mu\nu}=-\,\dfrac{1}{4}\,\bar{g}_{\mu\nu}\,\delta\phi-p_{\mu}p_{\nu}.$
(38)

For the wave travelling along the $z$ direction, the nonzero components of
$R_{\mu\nu}$ are $R_{tt}$, $R_{tz}$ and $R_{zz}$. Now to linear order we may
write,

$R_{\mu\nu\alpha\beta}\approx\dfrac{1}{2}(h_{\nu\alpha,\mu\beta}+h_{\mu\beta,\nu\alpha}-h_{\mu\alpha,\nu\beta}-h_{\nu\beta,\mu\alpha}).$

This can be further simplified as

$R_{itjt}=-\,\dfrac{1}{2}(h_{ij,tt}+h_{tt,ij}).$ (39)

Thus, using only the scalar part of the Eq. (35) we can write the simplified
and approximated geodesic deviation equation from the Eq. (37) in the limit
$l\rightarrow\infty$ as

$\ddot{x}^{i}=\dfrac{1}{2}\,\Big{(}\delta^{i}_{j}\delta\ddot{\phi}-\delta\phi_{\,,j}^{\,,i}\Big{)}x^{j},$
(40)

which gives,

$\displaystyle\ddot{x}$ $\displaystyle=\frac{1}{2}\,\delta\ddot{\phi}x,$ (41)
$\displaystyle\ddot{y}$ $\displaystyle=\frac{1}{2}\,\delta\ddot{\phi}y,$ (42)
$\displaystyle\ddot{z}$ $\displaystyle=0.$ (43)

Thus the geodesic equations show that there is no longitudinal component of
the GW polarization for the pure $R^{2}$ model. So, for this model besides the
tensor modes there exists only a scalar mode which is massless and pure in
nature, known as the breathing mode. Again, from the Eq. (35) considering only
the scalar part we may write,

$\delta\ddot{\phi}=-\,\omega^{2}\delta\phi.$ (44)

Using this equation in the Eq. (41) we have,

$\ddot{x}+\frac{1}{2}\,\omega^{2}\delta\phi\,x=0,$ (45)

which for a GW propagating along $z$ axis takes the form:

$\ddot{x}+\frac{1}{2}\,\phi_{0}\,\omega^{2}\cos\omega(t-z)x=0.$ (46)

The solution of this Eq. (46) gives the time variation of deviation of $x$ at
some fixed $z$. This equation is a special form of the well known Mathieu’s
equation canosa_numerical_1971 ; gheorghe_charged_0000  and the solution is
graphically shown in Fig.2 for two different set of parameters. It is seen
that the time period of geodesic deviation depends highly on $\omega$. The
results are identical along the $y$ direction also as clear from the Eq. (42).

Figure 2: Geodesic deviation along $x$ for a fixed $z$. The top plot is for
$\omega=1.5$, $\phi_{0}=0.02$ and $z=1$, and the bottom is for $\omega=1$,
$\phi_{0}=0.2$ and $z=1$. All parameters are in arbitrary units.

##  VI Polarization Modes With Newman-Penrose Formalism

The Newman-Penrose (NP) formalism newman_approach_1962  can be used to find
out the different polarization modes of GWs in a model. However, one major
drawback of NP formalism is that it is only applicable to null waves. In
$f(R)$ theory metric formalism, usually GWs have massive longitudinal mode of
polarization due to which the NP formalism fails and shows deviated results
liang_polarizations_2017 . But in the particular case of pure $R^{2}$ model,
as we have already seen that the scalar field is massless and thus this allows
us to use the NP formalism in the study of polarization modes of GWs in the
theory. In NP formalism a GW is described with the help of six amplitudes
$\\{\psi_{2},\psi_{3},\psi_{4},\phi_{22}\\}$ representing six polarization
modes in a particular coordinate system or frame newman_approach_1962 . All
these amplitudes are defined as eardley_gravitational-wave_1973 :

$\displaystyle\psi_{2}$ $\displaystyle=-\,\dfrac{1}{6}R_{ztzt},$
$\displaystyle\psi_{3}$
$\displaystyle=-\,\dfrac{1}{2}R_{xtzt}+\dfrac{1}{2}iR_{ytzt},$
$\displaystyle\psi_{4}$ $\displaystyle=-\,R_{xtxt}+R_{ytyt}+2iR_{xtyt},$
$\displaystyle\phi_{22}$ $\displaystyle=-\,R_{xtxt}-R_{ytyt}.$

Each of the complex amplitudes $\psi_{3}$ and $\psi_{4}$ are actually
equivalent to two real amplitudes eardley_gravitational-wave_1973 . In Brans-
Dicke theory, the massless scalar field appears as the breathing mode showing
$\phi_{22}=-\,R_{xtxt}-\,R_{ytyt}\neq 0.$ For our model these amplitudes are
found as

$\displaystyle\psi_{2}$ $\displaystyle=0,$ $\displaystyle\psi_{3}$
$\displaystyle=0,$ $\displaystyle\psi_{4}$
$\displaystyle=\ddot{\bar{h}}_{xx}-\ddot{\bar{h}}_{yy}+i\,(2\,\ddot{\bar{h}}_{xy}),$
$\displaystyle\phi_{22}$ $\displaystyle=-\,2\,\delta\ddot{\phi}.$

These results show that our model, i.e. pure $R^{2}$ model gives three modes
of massless polarization: two for tensor modes and one for the breathing mode
as found earlier, and according to Lorentz - invariant $E(2)$ classification
of plane waves, the model results GWs of class $N_{3}$ and all modes are
independent of the observer eardley_gravitational-wave_1973 . It should be
mentioned that a similar study was done in Ref. kausar_gravitational_2016  for
the same toy model in Palatini formalism. In this case, using the NP formalism
one can obtain,

$\displaystyle\psi_{2}$ $\displaystyle=0,$ $\displaystyle\psi_{3}$
$\displaystyle=0,$ $\displaystyle\psi_{4}$
$\displaystyle=\ddot{\bar{h}}_{xx}-\ddot{\bar{h}}_{yy}+i\,(2\,\ddot{\bar{h}}_{xy}),$
$\displaystyle\phi_{22}$ $\displaystyle=0.$

This shows that in Palatini formalism, the model $f(R)=\alpha R^{2}$ gives
only tensor modes of polarization of GWs and there is no scalar mode of
polarization or scalar degrees of freedom.

##  VII Pulsar Timing Arrays and Correlation of Polarization Modes

Pulsar timing arrays (PTAs) can be used to detect the polarization modes of
GWs. A good number of works have been going on in this field and PTAs are
found to be effective in the detection of extra polarization modes to be
present in GWs. Hence they can be used as a tool to test modified theories of
gravity lee_pulsar_2008 ; mingarelli_characterizing_2013 ; lee_detecting_2010
.

Presence of GWs disturb the null geodesic of the signals from pulsars. Due to
this reason the time of arrival of pulsar signal changes. So, by tracing the
changes in time of arrival of the radio signals from the pulsars, it is
possible to detect the GWs. For the mathematical treatment of the PTA
procedure let us consider a PTA detecting radio signals in the regime of GWs.
In case of GWs from the pure $R^{2}$ model, the information about the source
is carried by the three polarization amplitudes, viz., $h_{+}(t)$,
$h_{\times}(t)$ and $h_{\phi}(t)$, where the first two stands for the Einstein
or tensor modes and the third one for scalar (breathing) mode of polarization.
Thus the GW signal from a source for the pure $R^{2}$ model can be wriiten as

$h_{ij}(t,\hat{\Omega})=e^{+}_{ij}(\hat{\Omega})h_{+}(t,\hat{\Omega})+e^{\times}_{ij}(\hat{\Omega})h_{\times}(t,\hat{\Omega})+e^{\phi}_{ij}h_{\phi}(t,\hat{\Omega}).$
(47)

Here $e^{+}_{ij}$ and $e^{\times}_{ij}$ are the polarization tensors as given
by,

$\displaystyle
e^{+}_{ij}(\hat{\Omega})=\hat{m}_{i}\hat{m}_{j}-\hat{n}_{i}\hat{n}_{j},$
$\displaystyle
e^{\times}_{ij}(\hat{\Omega})=\hat{m}_{i}\hat{n}_{j}+\hat{n}_{i}\hat{m}_{j}.$

Hence, these polarization tensors are uniquely defined once the unit vectors
$\hat{m}$ and $\hat{n}$ used to describe the principal axies of wave are
specified. $\hat{\Omega}$ is the direction of GW propagation given by
$\hat{m}\times\hat{n}$. As per our previous convention, we set
$\hat{\Omega}=\hat{z}$, and consider a null vector $S^{\mu}$ that points to
the Solar system barycentre from the pulsar in Minkowski spacetime. In
purturbed spacetime vector $S^{\mu}$ will change to another vector
$\sigma^{\mu}$ as given by,

$\sigma^{\mu}=S^{\mu}-\dfrac{1}{2}\eta^{\mu\nu}h_{\mu\nu}S^{\nu}+\delta\phi\,\eta^{\mu\nu}b_{\mu\nu}S^{\nu},$
(48)

where the second part of this equation is due to tensor modes of purturbation
and the third part is due to the scalar mode of purturnation in spacetime with
$b_{\mu\nu}$ is a unit breathing mode matrix having two non-zero unit
components $b_{11}$ and $b_{22}$, obtained by the application of transverse
condition to this mode of GWs. It is to be noted that for this mode we can not
apply the traceless condition as the application of this condition to this
mode will retain only the tensor modes or GR modes kausar_gravitational_2016 .
Now, If we define $S^{\mu}$ as $S^{\mu}=\nu(1,-\alpha,-\beta,-\gamma)$, where
$\alpha,\beta$ and $\gamma$ are direction cosines, and $\nu$ is the frequencey
of the radio pulses from the source, then from Eq. (48), we may write,

$\displaystyle\sigma^{t}$ $\displaystyle=\nu\,$ $\displaystyle\sigma^{x}$
$\displaystyle=-\,\nu\,\big{[}\alpha\,\big{(}1-\dfrac{1}{2}h_{+}-\delta\phi\big{)}-\dfrac{\beta}{2}h_{\times}\big{]},$
$\displaystyle\sigma^{y}$
$\displaystyle=-\,\nu\,\big{[}\beta\big{(}1+\dfrac{1}{2}h_{+}+\delta\phi\big{)}-\dfrac{\alpha}{2}h_{\times}\big{]},$
$\displaystyle\mbox{and}\;\;\sigma^{z}$ $\displaystyle=-\,\nu\,\gamma.$

The radio pulses from the pulsar follow a null geodesic through spacetime. The
geodesic equation of the pulses with the affine parameter $\lambda$ is

$\dfrac{d\nu}{d\lambda}=-\,\nu^{2}\alpha\,\beta\,\dot{h}_{\times}+\dfrac{1}{2}\,\nu^{2}\big{[}\beta^{2}(\dot{h}_{+}-\delta\dot{\phi})-2\alpha^{2}(\dot{h}_{+}+\delta\dot{\phi})\big{]}.$
(49)

Using this equation in the geodesic deviation equation we can have the
frequency shift at Solar system barycentre as given by,

$\displaystyle z(t,\hat{\Omega})$
$\displaystyle=\;\mid\dfrac{\nu(t)-\nu_{0}}{\nu_{0}}\mid$ (50)
$\displaystyle=\dfrac{1}{2(1+\gamma)}\Big{[}\alpha^{2}(\Delta
h_{+}+\Delta\delta\phi)-\beta^{2}(\Delta h_{+}-\Delta\delta\phi)\Big{]}$
$\displaystyle+\dfrac{\alpha\beta}{(1+\gamma)}\Delta h_{\times},$

where $\nu(t)$ is the frequency observed at Solar system barycentre. Thus it
is seen that the presence of GWs will give rise to a frequency shift and these
shifts can be observed with the help of PTAs. It is also clear that the
presence of massless breathing mode results in a contribution of a shift in
the frequency besides a contribution from the usual tensor (GR) modes.
Following lee_pulsar_2008 , for stochastic GW background with normalized
frequency $f$, the correlation function for GR modes and breathing mode of
polarization are calculated. This function for tensor modes is found as

$C^{+,\times}(\theta)=\xi^{GR}(\theta)\int_{0}^{\infty}\dfrac{|h_{c}^{+,\times}|^{2}}{24\pi^{2}f^{3}}\,df,$
(51)

where

$\xi^{GR}(\theta)=\dfrac{3\,(1-\cos\theta)}{4}\log\,(\dfrac{1-\cos\theta}{2})+\dfrac{1}{2}-\dfrac{1-\cos\theta}{8}+\dfrac{\delta(\theta)}{2},$

and $\theta$ is the angular separation between two pulsars. For the scalar
modes it is

$C^{b}(\theta)=\xi^{b}(\theta)\int_{0}^{\infty}\dfrac{|h_{c}^{b}|^{2}}{12\pi^{2}f^{3}}\,df,$
(52)

where

$\xi^{b}(\theta)=\dfrac{1}{8}\big{[}\cos\theta+3+4\,\delta(\theta)\big{]}.$

The correlation coefficients $\xi^{GR}$ and $\xi^{b}$ as a function of
$\theta$ are plotted in the Fig.3 It is seen that provided the cases in
lee_pulsar_2008  hold good, PTAs can effectively distinguish between the
breathing mode and tensor modes present in GWs. It needs to mention that the
correlation versus angular separation curve for the tensor modes of
polarizations was first obtained by Hellings and Downs in 1983 and hence it is
usually known as Hellings-Downs (HD) curve hellings_upper_1983 ;
jenet_understanding_2014 .

Figure 3: Variation of correlation coefficient $\xi(\theta)$ of different
polarization modes with respect to angular separation between two pulsars.

The calculations of correlation functions are model independent, but for
different polarization modes there are different correlation functions
lee_detecting_2010 . Since our toy model $f(R)=\alpha R^{2}$ gives only tensor
plus, cross and massless breathing modes of polarization of GWs, we have
included the correlation functions for tensor modes and massless breathing
mode only. This provides an experimental way to test the viability of the toy
model. However, in other $f(R)$ gravity models (in metric formalism), there
are four polarization modes viz. tensor plus, cross, massless breathing and
massive longitudinal modes of polarizations kausar_gravitational_2016 ;
liang_polarizations_2017 . Therefore, in those models, the experimental
viability can be checked by considering three correlation functions for tensor
modes, massless transverse breathing mode and massive longitudinal mode. In
the near future, with increased sensitivity and sufficient data, the
experimentally obtained correlation functions from pulsar timing array data
can hopefully help us to distinguish between different polarization modes
which will provide us a way to check the viability of the toy model.

##  VIII Summary and Conclusions

In this work, we have used the $f(R)$ gravity power law model to study the
polarization modes of Gravitational Waves in de Sitter spacetime. We have seen
that the field equation in de Sitter spacetime gives a massive scalar mode
with a mass term $m^{2}=\dfrac{(2-n)\bar{R}}{(n-1)}.$ The mass term varies
widely with the exponent term $n$ and the background curvature or the de
Sitter curvature, and the mass of the scalar field becomes zero when
background curvature is zero or $n=2$. Later using the stability condition for
the theory in de Sitter spacetime we have seen that for a constant curvature
$\bar{R}$ the theory has stable solutions for $n=2$. It has been observed that
for this particular case of the model, the massive longitudinal mode of
polarization vanishes. Thus, this is the only case in $f(R)$ theory in metric
formalism treatment where only three polarization modes are obtained and the
third mode is a pure breathing mode. To validate this result we’ve studied the
geodesic deviations for the pure $R^{2}$ model explicitly and later we used NP
formalism to confirm the validity of our result. The absence of the massive
scalar field in this model allows us to use NP formalism, which is a powerful
tool to check the polarization contents of null GWs. The results from the
analysis show absence of massive polarization mode of the GWs in the theory,
which establish that the polarization modes and the mass of the scalar field
in $f(R)$ gravity are model dependent. Another important result obtained from
this work is that the scalar field is independent of the background curvature
in pure $R^{2}$ model of $f(R)$ gravity and evidently chameleonic behaviour is
not observed here. The absence of massive longitudinal mode makes this pure
$R^{2}$ model different from the other $f(R)$ gravity models and hence a more
detail study in this model is required.

The pure $R^{2}$ model can be extended by adding other terms to the action
which can account for the missing part giving rise to mixed polarization
states of GWs besides tensor modes of polarization kausar_gravitational_2016 ;
liang_polarizations_2017 . Those extended models can be constrained with the
experimental results obtained so far. It is expected that future experiments
can provide better constraints to the $f(R)$ gravity models using which the
existence of extra polarization modes can be hopefully confirmed, which in
turn give options to test the reliability of $f(R)$ gravity in the
modifications and extensions of GR. Extensions like nonminimal matter field
coupling and other modified gravity theories can be included as the future
aspects of this type of works, which can be tested with future experiments on
GWs.

## Acknowledgments

UGD is thankful to the Inter-University Centre for Astronomy and Astrophysics
(IUCAA), Pune, India for hospitality during his visits as a Visiting
Associate. The authors thank anonymous reviewers for their valuable
suggestions and comments which helped to improve the manuscript.

## References

  * (1) H. Weyl, “Eine neue Erweiterung der Relativitätstheorie”, Annalen der Physik 364, 101 (1919).
  * (2) A. S. Eddington, The mathematical theory of relativity, Cambridge Univ. Press, Cambridge (1924).
  * (3) R. Utiyama and B. S. DeWitt,“Renormalization of a classical gravitational field interacting with quantized matter”, J. Math. Phys. 3, 608 (1962).
  * (4) T.P. Sotiriou and V. Faraoni, “f (R) theories of gravity”, Rev. Mod. Phys. 82, 451 (2010) [arXiv:0805.1726].
  * (5) S. Nojiri and S. D. Odintsov, “Unified cosmic history in modified gravity: from F(R) theory to Lorentz non-invariant models”, Phys. Rept. 505, 59 (2011) [arXiv:1011.0544].
  * (6) S. Nojiri and S. D. Odintsov, “Introduction to Modified Gravity and Gravitational Alternative for Dark Energy”, Int. J. Geom. Meth. Mod. Phys. 4, 115 (2007) [arXiv:hep-th/0601213].
  * (7) S. Nojiri, S. D. Odintsov and V. K. Oikonomou, “Modified Gravity Theories on a Nutshell: Inflation, Bounce and Late-time Evolution”, Phys. Rept. 692, 1 (2017) [arXiv:1705.11098].
  * (8) B. P. Abbott and R. Abbott et al., “Observation of Gravitational Waves from a Binary Black Hole Merger”, Phys. Rev. Lett. 116, 061102 (2016) [doi: 10.1103/PhysRevLett.116.061102].
  * (9) B. P. Abbott and R. Abbott et al., “GW151226: Observation of Gravitational Waves from a 22-Solar-Mass Binary Black Hole Coalescence”, Phys. Rev. Lett. 116, 241103 (2016) [doi: 10.1103/PhysRevLett.116.241103].
  * (10) B. P. Abbott and R. Abbott et al., “GW170104: Observation of a 50-Solar-Mass Binary Black Hole Coalescence at Redshift 0.2”, Phys. Rev. Lett. 118, 221101 (2017) [doi: 10.1103/PhysRevLett.118.221101].
  * (11) B. P. Abbott and R. Abbott et al., “GW170814: A Three-Detector Observation of Gravitational Waves from a Binary Black Hole Coalescence”, Phys. Rev. Lett. 119, 141101 (2017) [doi: 10.1103/PhysRevLett.119.141101].
  * (12) B. P. Abbott and R. Abbott et al., “GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral”, Phys. Rev. Lett. 119, 161101 (2017) [doi:10.1103/PhysRevLett.119.161101].
  * (13) C. Corda, “Massive gravitational waves from the $R^{2}$ theory of gravity: Production and response of interferometers”, Int. J. Mod. Phys. A 23 (10), 1521 (2008).
  * (14) C. Corda, “The production of matter from curvature in a particular linearized high order theory of gravity and the longitudinal response function of interferometers”, JCAP 0704, 009 (2007).
  * (15) S. Capozziello, C. Corda and M. F. De Laurentis, “Massive gravitational waves from $f(R)$ theories of gravity:Potential detection with LISA”, Physics Letters B 669, 255 (2008).
  * (16) Y. S. Myung, “Propagating degrees of freedom in f(R) gravity”, arXiv:1608.01764 (2016).
  * (17) D. Liang, Y. Gong, S. Hou and Y. Liu, “Polarizations of gravitational waves in $f(R)$ gravity”, Phys. Rev. D 95, 104034 (2017) [arXiv:1701.05998].
  * (18) L. G. Jaime, L. Patio and M. Salgado, “About matter and dark-energy domination eras in $R^{n}$ gravity or lack thereof” Phys. Rev. D 87, 024029 (2013) [arXiv:1212.2604].
  * (19) L. Amendola, D. Polarski and S. Tsujikawa, “Power-laws f(R) theories are cosmologically unacceptable”, IJMP D 16, 1555 (2007) [arXiv:astro-ph/0605384].
  * (20) L. Amendola, D. Polarski and S. Tsujikawa, “Are f(R) dark energy models cosmologically viable ?”, Phys. Rev. Lett. 98 131302 (2007) [arXiv:astro-ph/0603703].
  * (21) L. Amendola, R. Gannouji, D. Polarski and S. Tsujikawa, “Conditions for the cosmological viability of f(R) dark energy models”, Phys. Rev. D 75, 083504 (2007) [arXiv:gr-qc/0612180].
  * (22) S. Capozziello, S. Carloni and A. Troisi, “Quintessence without scalar fields”, arXiv:astro-ph/0303041 (2003).
  * (23) S. Carloni, P. K. S. Dunsby, S. Capozziello and A. Troisi, “Cosmological dynamics of $R^{n}$ gravity”, Class. Quant. Grav. 22, 4839 (2005) [arXiv:gr-qc/0410046].
  * (24) U. D. Goswami and K. Deka, “Cosmological dynamics of f(R) gravity scalar degree of freedom in Einstein frame”, IJMP D 22, 1350083 (2013) [arXiv:1303.5868].
  * (25) V. Faraoni, “$R^{n}$ gravity and the chameleon”, Phys. Rev. D 83, 124044 (2011).
  * (26) A. F. Zakharov et al., “Solar system constraints on $R^{n}$ gravity”, Phys. Rev. D 74, 107101 (2006)[arXiv:astro-ph/0611051].
  * (27) S. Capozziello, “Gravitational lensing in fourth order gravity”, Phys. Rev. D 73, 104019 (2006) [arXiv:astro-ph/0604435].
  * (28) D. Borka et al., “Constraints on $R^{n}$ gravity from precession of orbits of S2-like stars”, Phys. Rev. D 85, 124004 (2012) [arXiv:1206.0851].
  * (29) A. F. Zakharov et al., “Constraints on $R^{n}$ gravity from precession of orbits of S2-like stars: a case of a bulk distribution of mass”, Adv. in Space Research 54, 1108 (2014) [arXiv:1407.0366].
  * (30) J. D. Barrow and A. C. Ottewill, “The Stability of General Relativistic Cosmological Theory”, J. Phys. A: Math. Gen. 16, 2757 (1983).
  * (31) V. Faraoni and S. Nadeau, “The stability of modified gravity models”, Phys. Rev. D 72, 124005 (2005) [arXiv:gr-qc/0511094].
  * (32) T. Chiba, “$\dfrac{1}{R}$ gravity and scalar-tensor gravity”, Phys. Lett. B 575, 1 (2003).
  * (33) Louis Yang et al., “Gravitational waves in viable $f(R)$ models”, JCAP 08, 029 (2011).
  * (34) A. Araujo and J.G. Pereira, “Entropy in locally-de Sitter spacetimes”, Int. J. Mod. Phys. D. 24, 1550099 (2015).
  * (35) S. W. Hawking and G. F. R. Ellis, “The Large Scale Structure of Space-Time”, Cambridge University Press (1973).
  * (36) J.G. Pereira and A.C. Sampson, “de Sitter geodesics: reappraising the notion of motion”, Gen Relativ Gravit. 44, 1299 (2012).
  * (37) J. Podolský and R. Švarc, “Interpreting spacetimes of any dimension using geodesic deviation”, Phys. Rev. D. 85, 044057 (2012).
  * (38) J. Canosa, “Numerical solution of Mathieu’s equation”, J. Comput. Phys. 7, 255 (1971).
  * (39) V. N. Gheorghe and G. Werth, “Charged particle traps: Physics and techniques of charged particle field confinement”, (Physica-Verlag NetLibrary, Inc. distributor, Heidelberg; Boulder, 2005), ISBN 9783540265764.
  * (40) E. Newman and R. Penrose, “An approach to gravitational radiation by a method of spin coefficients”, J. Math. Phys. 3, 566 (1962).
  * (41) D. M. Eardley, D. L. Lee and A. P. Lightman, “Gravitational-wave observations as a tool for testing relativistic gravity”, Phys. Rev. D 8, 3308 (1973).
  * (42) H. R. Kausar, L. Philippoz and P. Jetzer, “Gravitational wave polarization modes in f (R) theories”, Phys. Rev. D 93, 124071 (2016) [arXiv:1606.07000].
  * (43) K. J. Lee, F. A. Jenet and R. H. Price, “Pulsar timing as a probe of non‐Einsteinian polarizations of gravitational waves”, The Astrophysical Journal 685, 1304 (2008).
  * (44) C. M. F. Mingarelli, T. Sidery, I. Mandel and A. Vecchio, “Characterizing gravitational wave stochastic background anisotropy with pulsar timing arrays”, Phys. Rev. D 88, 062005 (2013) [arXiv:1306.5394].
  * (45) K. Lee, F. A. Jenet, R. H. Price, N. Wex and M. Kramer, “Detecting massive gravitons using pulsar timing arrays”, The Astrophysical Journal 722, 1589 (2010) [arXiv:1008.2561].
  * (46) R. W. Hellings and G. S. Downs, “Upper limits on the isotropic gravitational radiation background from pulsar timing analysis”, The Astrophysical Journal 265, L39 (1983).
  * (47) F. A. Jenet and J. D. Romano, “Understanding the gravitational-wave Hellings and Downs curve for pulsar timing arrays in terms of sound and electromagnetic waves”, arXiv:1412.1142 (2014).


# Radial Kohn-Sham problem via integral-equation approach

Jānis Užulis  Andris Gulans  Department of Physics, University of Latvia,
Jelgavas iela 3, Riga, LV-1058 Latvia

###### Abstract

We present a numerical tool for solving the non-relativistic Kohn-Sham problem
for spherically-symmetric atoms. It treats the Schrödinger equation as an
integral equation relying heavily on convolutions. The solver supports
different types of exchange-correlation functionals including screened and
long-range corrected hybrids. We implement a new method for treating range
separation based on the complementary error function kernel. The present tool
is applied in non-relativistic total energy calculations of atoms. A
comparison with ultra-precise reference data[Cinal, JOMC 58, 1571 (2020)]
shows a 14-digit agreement for Hartree-Fock results. We provide further
benchmark data obtained with 5 different exchange-correlation functionals.

##  I Introduction

Particle in a spherically symmetric potential is a classic problem that enters
quantum mechanics textbooks. This problem has a fundamental importance in
modern computational materials science, as it has to be solved in electronic-
structure codes for a variety of purposes. One comes across it when generating
pseudopotentials [1], constructing basis functions in several formalisms such
as numerical atomic orbitals[2], linearized augmented plane waves (LAPW) and
linearized muffin-tin orbitals[3, 4]. This problem also appears in a
calculation of core orbitals in LAPW.

At present, the majority of electronic-structure calculations employ
(semi)local density-functional theory (DFT), as evident from data-centered
high-throughput initiatives [5, 6, 7]. In this case, the exchange and
correlation are described by means of the local-density approximation (LDA) or
the generalized gradient approximation (GGA). These models lead to the Kohn-
Sham (KS) equation

$-\frac{\nabla^{2}}{2}\psi(\mathbf{r})+v(\mathbf{r})\psi(\mathbf{r})=\varepsilon\psi(\mathbf{r})$
(1)

with a local effective KS potential $v(\mathbf{r})$. Assuming the spherical
symmetry of the potential, this equation reduces to a one-dimensional problem
that can be solved by an outward integration as an initial value problem. This
task is typically performed using the Numerov’s method[8], although other
options such as the predictor-corrector and Bulirsch-Stoer[9] method are also
applied in electronic-structure codes [10, 11]. If Eq. 1 is intended as an
eigenvalue problem, the outward integration is coupled with the shooting
method. In other words, the procedure of solving the radial problem is already
well-established in the case of a local potential.

The (semi)local approximations in DFT perform with limitations for a range of
applications that require accurate atomization (reaction) energies [12] or
good estimates of the band structure. A frequently applied solution to this
problem is hybrid exchange-correlation functionals where the GGA and/or LDA
are combined with the Fock, screened [13] or the long-range [14] exchange.
These three exchange models introduce a non-local potential
$\hat{v}^{\mathrm{nl}}$ that makes the usual approach to solving the radial
problem not applicable. Two recent studies addressed this problem with an aim
to implement an efficient and precise solver [15, 16]. Cinal introduced a
compact basis and used it with the Hartree-Fock (HF) equations [15]. The
obtained total energies of atoms were computed in double and quadruple
precision ensuring an extreme precision level. In the other study, Lehtola
implemented a radial solver based on the finite-element method and, yet again,
employed a compact basis with a small number of support points [16]. Both, HF
and DFT calculations, were performed in this study, and the reported total
energies were given with the decimal signs up to $10^{-6}$ Ha. These two
implementations have a drawback that they require a specific radial grid which
is not compatible with the applications in electronic-structure codes
explained above.

In this study, we explore an alternative approach for solving the spherically-
symmetric problem. It is suitable for calculations with the types of grids
that are employed in electronic-structure codes and can be integrated into
them. Our method uses the idea that the Schrödinger equation can be rewritten
as an integral equation following the ideas first published by Kalos in
1962[17]. This approach is commonly used in calculations with wavelet basis
sets in Ref. [18, 19, 20, 21, 22], because it allows one to avoid computation
of derivatives. This method was also applied in problems with model potentials
[23] and for optimising orbitals in small molecules within the HF theory [24].
However, we are not aware of previous studies applying this method in radial
solvers.

This paper is structured as follows. We give an introduction to the integral-
equation method in Sec. II and provide expressions for calculating the total
energy in Sec. III. The formalism described in these sections relies heavily
on convolutions and radial integrals discussed in Secs. IV and II,
respectively. The considered types of convolutions involve the Coulomb and
Yukawa kernels as well as complimentary error function kernel (erfc) kernel
(required in hybrid functionals with screened exchange) for which we propose a
new method. Our implementation employs a large number of integral evaluations
requiring attention to how they are performed. We show an efficient high-order
approach for computing the integrals on arbitrary grids. In Sec. VI, we
evaluate the performance of logarithmic and inverse polynomial radial grids.
Finally, total-energy calculations of closed-shell atoms are presented and
analysed in Sec. VII.

##  II Integral-equation method

We consider the non-relativistic KS equation shown in Eq. 1. We assume that
the potential consists of local and non-local contributions
$\hat{v}=v^{\mathrm{L}}(\mathbf{r})+\hat{v}^{\mathrm{NL}}$. The local part
comprises the electron-nuclear, Hartree and the LDA or GGA exchange-
correlation terms expressed as

$v_{\mathrm{n}}(\mathbf{r})=-Z/r$ (2)

with the nuclear charge $Z$,

$v_{\mathrm{H}}(\mathbf{r})=\int\frac{\rho(\mathbf{r}^{\prime})}{|\mathbf{r}^{\prime}-\mathbf{r}|}d\mathbf{r}^{\prime}$
(3)

with the electron density $\rho(\mathbf{r})$, and

$v^{\mathrm{L}}_{\mathrm{xc}}(\mathbf{r})=(1-\alpha)v^{\mathrm{GGA}}_{\mathrm{x}}(\mathbf{r})+v^{\mathrm{GGA}}_{\mathrm{c}}(\mathbf{r}),$
(4)

respectively. The non-local part is defined by its action on a trial
wavefunction $\psi_{n}(\mathbf{r})$:

$\hat{v}^{\mathrm{NL}}_{\mathrm{x}}\psi_{n}(\mathbf{r})=\alpha\int
v^{\mathrm{NL}}_{\mathrm{x}}(\mathbf{r},\mathbf{r}^{\prime})\psi_{n}(\mathbf{r}^{\prime})d\mathbf{r}^{\prime}$
(5)

containing the (screened) Fock exchange

$v^{\mathrm{NL}}_{\mathrm{x}}(\mathbf{r},\mathbf{r}^{\prime})=\sum\limits_{n}f_{n}\psi_{n}(\mathbf{r})V(|\mathbf{r}^{\prime}-\mathbf{r}|)\psi^{\ast}_{n}(\mathbf{r}^{\prime}),$
(6)

where the sum runs over all considered orbitals, $f_{n}$ is the occupation
number and $V(|\mathbf{r}|)$ is the interaction kernel (Coloumb, Yukawa, erfc
or erf). Finally, the parameter $\alpha$ is the weight of the non-local
exchange in the potential. We evaluate
$v^{\mathrm{L}}_{\mathrm{xc}}(\mathbf{r})$ using the libxc library[25]. The
details on calculating $v_{\mathrm{H}}(\mathbf{r})$ and
$\hat{v}^{\mathrm{NL}}_{\mathrm{x}}\psi_{n}(\mathbf{r})$ are given in Sec. IV.

The effective potential depends on the electron density, and, therefore, the
KS problem has to be solved self-consistently. There are numerous methods for
ensuring convergence to the ground-state solution (see the review in Ref.
[26]), but a method as straightforward as linear mixing of potentials works
sufficiently well for spherical atoms. Therefore, we focus on solving the KS
equation for a fixed potential.

A formal rearrangement of terms in Eq. 1 leads to

$\psi(\mathbf{r})=2(\nabla^{2}-\lambda^{2})^{-1}[\hat{v}\psi(\mathbf{r})]$ (7)

with $\lambda^{2}=-2\varepsilon$. An expression with
$\psi(\mathbf{r})=(\nabla^{2}-\lambda^{2})^{-1}f(\mathbf{r})$ has the meaning
that $\psi(\mathbf{r})$ is the solution of the screened Poisson equation

$(\nabla^{2}-\lambda^{2})\psi(\mathbf{r})=-f(\mathbf{r}),$ (8)

where the usual prefactor $4\pi$ on the right-hand side is omitted. If this
equation is applied in the context of electrostatics, the functions
$f(\mathbf{r})$ and $\psi(\mathbf{r})$ have the meaning of the charge and the
resulting screened potential, respectively.

We solve the KS equation for bound states meaning that $\varepsilon<0$ and
$\psi(\mathbf{r})$ decays to 0 as $r\rightarrow\infty$. Then, it is
appropriate to express $(\nabla^{2}-\lambda^{2})^{-1}$ in Eq. 7 via its
Green’s function in the following manner:

$\psi(\mathbf{r})=2\int\frac{e^{-\lambda|\mathbf{r}-\mathbf{r}^{\prime}|}}{4\pi|\mathbf{r}-\mathbf{r}^{\prime}|}\hat{v}\psi(\mathbf{r}^{\prime})d\mathbf{r}^{\prime}.$
(9)

Thus, the KS problem given in Eq. 1 is transformed into the integral equation
that does not contain any differential operators. In spherically symmetric
atoms, a wavefunction can be expressed according to its quantum numbers $n$,
$\ell$ and $m$ as

$\psi_{n\ell m}(\mathbf{r})=u_{n\ell}(r)Y_{\ell m}(\hat{r}),$ (10)

where $u_{n\ell}(r)$ is the radial part, and only one spherical harmonic
$Y_{\ell m}(\hat{r})$ enters the expression. A calculation of the convolution
in Eq. 9 reduces to an evaluation of a one-dimensional integral as discussed
below in Sec. IV.

Following Ref. [17], Eq. 9 is solved iteratively by inserting a trial function
$\psi_{n\ell m}^{(i)}(\mathbf{r})$ to the right-hand side and obtain an
updated function $\psi_{n\ell m}^{(i+1)}(\mathbf{r})$, where $i$ is the step
number. The parameter $\lambda$ is updated along with the wavefunction by
evaluating

$\varepsilon^{(i)}=\langle\psi_{n\ell m}^{(i)}|\hat{H}|\psi_{n\ell
m}^{(i)}\rangle/\langle\psi_{n\ell m}^{(i)}|\psi_{n\ell m}^{(i)}\rangle.$ (11)

The described procedure converges to the lowest-energy solution for given
quantum numbers $\ell$ and $m$.

Now follows a description of a procedure that we use for solving for a few
orbitals with a given $\ell$ and $m$ with the lowest Kohn-Sham energies.
Suppose that, after $i$ iterations, estimates for the KS energies
$\varepsilon^{(i)}_{n\ell}$ and wavefunctions $\psi_{n\ell
m}^{(i)}(\mathbf{r})$ have been obtained. In step 1, we construct basis
functions

$\chi_{n}^{(i)}(\mathbf{r})=2(\nabla^{2}+\varepsilon_{n\ell}^{(i)})^{-1}\psi_{n\ell
m}^{(i)}(\mathbf{r}).$ (12)

In step 2, we calculate the matrix elements

$H_{nn^{\prime}}=\langle\chi_{n}^{(i)}|-\nabla^{2}/2+\hat{v}|\chi_{n^{\prime}}^{(i)}\rangle$
(13)

and

$S_{nn^{\prime}}=\langle\chi_{n}^{(i)}|\chi_{n^{\prime}}^{(i)}\rangle.$ (14)

In step 3, we solve the matrix eigenproblem

$Hz=\sigma Sz,$ (15)

where $\sigma$ and $z$ are an eigenvalue and an eigenvector, respectively. In
step 4, we update the estimates

$\varepsilon_{n\ell}^{(i+1)}=\sigma_{n}$ (16)

and

$\psi_{n\ell
m}^{(i+1)}(\mathbf{r})=\sum_{n^{\prime}}z_{nn^{\prime}}\chi_{n^{\prime}}(\mathbf{r}).$
(17)

If $|\varepsilon_{n\ell}^{(i+1)}-\varepsilon_{n\ell}^{(i)}|$ is greater than a
predefined threshold value, the procedure is repeated from step 1. Once this
process converges, $\varepsilon_{n\ell}^{(i)}$ and $\psi_{n\ell m}^{(i)}$ are
the eigenpairs of the KS equation for a given potential.

Figure 1:  Convergence of inner loop (difference between each eigenvalue and
its converged value depending on iteration number) for Rn ($\ell=0$ channel)
in first external cycle iteration.

Taking into account that the complete problem has to be solved self-
consistently, we find it more practical to use $\psi_{n\ell
m}^{(0)}(\mathbf{r})$ in Eq. 12 in all steps. With such an approach, only the
eigenvalue estimates $\varepsilon_{n\ell}^{(i)}$ are updated, and the obtained
result (an energy and an orbital) for a fixed potential is not exactly an
eigenpair. However, this simplification provides a good enough approximation
for $\varepsilon_{n\ell}$ and $\psi_{n\ell m}(\mathbf{r})$ such that the self-
consistency procedure always converges to the solution of the KS problem with
an excellent precision. The only exception is when KS energies of an atom are
close to 0 (typically $\sim-10^{-2}$ Ha). In such a case, however, the full
procedure with the update of the wavefunction fails to converge too.

To illustrate the performance of the method, we solve the KS equation for the
Rn atom. The initial guess $\psi_{n\ell m}^{(0)}(\mathbf{r})$ is set to the
wavefunctions in hydrogen-like ion with the nuclear charge of Rn. Fig. 1 shows
the performance of the described algorithm for the $s$-orbitals. The KS
energies converge to the eigenvalue estimates within $10^{-8}$ Ha in 12-22 and
52 steps in the case of the core and valence orbitals, respectively. The
quality of the initial guess improves with every step of the self-consistency
step, and the number of the integral-equation iterations reduces dramatically.

##  III Total energy

To calculate the total energy, we consider it as a sum of the five components:

$E_{\mathrm{tot}}=E_{\mathrm{kin}}+E_{\mathrm{n}}+E_{\mathrm{H}}+E_{\mathrm{xc}}^{\mathrm{L}}+E_{\mathrm{x}}^{\mathrm{NL}}.$
(18)

where $E_{\mathrm{kin}}$ is the kinetic energy and the remaining terms are
labelled consistently with the contribution to the effective potential in Eqs.
2–5. These components are evaluated using radial integrals as follows:

$E_{\mathrm{kin}}=\frac{1}{2}\sum_{n\ell}f_{n\ell}\int\limits_{0}^{\infty}\biggl{(}(ru^{\prime}_{n\ell}(r))^{2}+\ell(\ell+1)u^{2}_{n\ell}(r)\biggr{)}dr,$
(19)

where $f_{n\ell}$ is the occupation number of the electron shells,

$E_{\mathrm{ext}}=\int_{0}^{\infty}-Z\rho(r)rdr,$ (20)
$E_{\mathrm{H}}=\frac{1}{2}\int_{0}^{\infty}v_{\mathrm{H}}(r)\rho(r)r^{2}dr,$
(21)
$E_{\mathrm{xc}}^{L}=\int_{0}^{\infty}\epsilon_{\mathrm{xc}}^{L}(r)\rho(r)r^{2}dr,$
(22)

where $\epsilon_{\mathrm{xc}}^{L}(r)$ is the (semi)local exchange-correlation
energy density per particle, and

$E_{\mathrm{x}}^{\mathrm{NL}}=\frac{1}{2}\sum_{n\ell}f_{n\ell}\int_{0}^{\infty}\psi_{n\ell}(r)[\hat{v}_{\mathrm{x}}^{\mathrm{NL}}\psi_{n\ell}(r)]r^{2}dr.$
(23)

##  IV  Convolutions with bare and screened Coulomb kernels

The algorithm described above employs convolutions with types of kernels: (i)
the Coulomb kernel $V_{\mathrm{C}}(\mathbf{r})=1/r$, (ii) the Yukawa kernel
$V_{\mathrm{Y}}(\mathbf{r})=e^{-\lambda r}/r$ and (iii) the erfc kernel
$V_{\mathrm{SR}}(\mathbf{r})=\mathrm{erfc}(\mu r)/r$.

First, we consider the Coulomb kernel $V(\mathbf{r})=1/r$. It is applied in a
calculation of the Hartree potential $v_{\mathrm{H}}(\mathbf{r})$ in Eq. 3 and
the non-local Fock exchange. The standard approach is to apply the Laplace
expansion:

${}\frac{1}{|\mathbf{r}-\mathbf{r}^{\prime}|}=\sum_{\ell=0}^{\infty}\frac{4\pi}{2\ell+1}\frac{r^{\ell}_{<}}{r^{\ell+1}_{>}}Y^{\ast}_{\ell
m}(\hat{r})Y_{\ell m}(\hat{r}^{\prime}),$ (24)

where $r_{<}=\mathrm{min}(r,r^{\prime})$ and
$r_{>}=\mathrm{max}(r,r^{\prime})$. This equation allows us to express the
result of the convolution or simply the potential $v(\mathbf{r})$ due to the
density $\rho_{\ell m}(r)Y_{\ell m}(\hat{r})$ as

$\displaystyle{}v(\mathbf{r})=\frac{4\pi Y_{\ell m}(\hat{r})}{2\ell+1}$
$\displaystyle\left[\frac{1}{r^{\ell+1}}\int\limits_{0}^{r}r^{\prime\ell+2}\rho_{\ell
m}(r^{\prime})dr^{\prime}\right.$ (25)
$\displaystyle+\left.r^{\ell}\int\limits_{r}^{\infty}\frac{1}{r^{\prime\ell-1}}\rho_{\ell
m}(r^{\prime})dr^{\prime}\right].$

This equation decouples $r$ from $r^{\prime}$ and, thus, makes the evaluation
of the potential efficient.

In a calculation of the Hartree potential, we assume spherically symmetric
density $\rho(\mathbf{r})=\rho_{00}(r)Y_{00}(\hat{r})$. Also in the case of
the Fock exchange, the spherical symmetry is assumed, and the result of the
non-local operator acting on a trial wavefunction can be expressed as
$\hat{v}_{x}^{\mathrm{NL}}\chi_{n\ell
m}(\mathbf{r})=\hat{v}_{x}^{\mathrm{NL}}\left[\tilde{u}_{n\ell}(r)Y_{\ell
m}(\hat{r})\right]=h(r)Y_{\ell m}(\hat{r})$. Using Eqs. 5, 6 and 24, the
radial part of the result is obtained as

$\displaystyle h(r)=$
$\displaystyle-\frac{1}{2}\sum_{n^{\prime}\ell^{\prime}}\sum\limits_{\ell^{\prime\prime}=|\ell-\ell^{\prime}|}^{\ell+\ell^{\prime}}{\\!\\!\\!\\!\\!\\!\\!}{{}^{\prime}}f_{n^{\prime}\ell^{\prime}}\begin{pmatrix}\ell&\ell^{\prime}&\ell^{\prime\prime}\\\
0&0&0\end{pmatrix}^{2}u_{n^{\prime}\ell^{\prime}}(r)$ (26)
$\displaystyle\cdot\left(\int_{0}^{\infty}dr^{\prime}\frac{r_{<}^{\ell^{{}^{\prime\prime}}}}{r_{<}^{\ell^{{}^{\prime\prime}}+1}}r^{\prime
2}u_{n^{\prime}\ell^{\prime}}(r^{\prime})\tilde{u}_{n\ell}(r^{\prime})\right),$

where the Wigner 3-j symbol is applied and $\sum^{\prime}$ has the meaning of
of the summations with step 2 (see Ref. [15, 27]).

Similarly to $V_{\mathrm{C}}(\mathbf{r}-\mathbf{r}^{\prime})$, the Yukawa
kernel can be written in a separable form:

$\frac{e^{-\lambda|\mathbf{r}-\mathbf{r}^{\prime}|}}{|\mathbf{r}-\mathbf{r}^{\prime}|}=4\pi\lambda\sum_{\ell=0}^{\infty}i_{\ell}(\lambda
r_{<})k_{\ell}(\lambda r_{>})\sum_{m=-\ell}^{\ell}Y^{*}_{\ell
m}(\hat{r})Y_{\ell^{\prime}m^{\prime}}(\hat{r}^{\prime}),$ (27)

where $\lambda$ is the screening parameter, $i_{\ell}$ and $k_{\ell}$ are the
modified spherical Bessel functions of the first and second kinds,
respectively. This equation allows us to factorize the kernel into parts that
depend on either on $r$ or $r^{\prime}$ similarly to how it is done in the
case of the Coulomb kernel in Eq. 25.

The third considered convolution kernel $V_{\mathrm{SR}}(\mathbf{r})$ appears
in hybrid exchange-correlation functionals where the range separation is
introduced. These methods decompose the Coulomb interaction into the short-
and long-range contributions (SR and LR, respectively) as follows:

$\frac{1}{|\mathbf{r}-\mathbf{r}^{\prime}|}=\underbrace{\frac{\mathrm{erfc}(\mu|\mathbf{r}-\mathbf{r}^{\prime}|)}{|\mathbf{r}-\mathbf{r}^{\prime}|}}_{SR}+\underbrace{\frac{\mathrm{erf}(\mu|\mathbf{r}-\mathbf{r}^{\prime}|)}{|\mathbf{r}-\mathbf{r}^{\prime}|}}_{LR},$
(28)

where $\mu$ is known as the screening or range-separation parameter. Only one
of these two terms remains in the non-local exchange in Eq. 6, whereas the
other one is fully omitted. Hybrid functionals designed for solid-state
applications keep only the SR part and thus mimic screening in bulk materials.
In calculations of molecules, it is more common to require the correct long-
range behaviour of the exchange functional and then the LR term is kept.

Table 1: Fitting parameters $A$ and $a$ for the complementary error function. The definition is given in Eq. 29.  $\Re(A)$ | $\Im(A)$ | $\Im(a)$
---|---|---
$N_{\mathrm{f}}=3$ $\Re(a)$=3.62435558273902
3.34075293484632 | 0.00000000000000 | 0.00000000000000
-1.14940001089277 | 1.50422252025891 | 1.43494484832716
-0.02638283853366 | -0.32974525523286 | 2.93170108446798
0.00540643159457 | 0.01222658029598 | 4.59735152907432
$N_{\mathrm{f}}=4$ $\Re(a)$=4.03651233059082
5.88649164984020 | 0.00000000000000 | 0.00000000000000
-2.42999364511904 | 2.89879784656921 | 1.27567126568312
-0.06523181639778 | -0.98148171395886 | 2.58394634761278
0.05351129281419 | 0.07580726842369 | 3.97319627725743
-0.00153165601674 | -0.00161854017762 | 5.52839978282953
$N_{\mathrm{f}}=5$ $\Re(a)$=4.44903651057752
11.68016225012443 | 0.00000000000000 | 0.00000000000000
-5.58814846386175 | 5.91087741873711 | 1.16136522308687
0.04005009651638 | -2.71301646794361 | 2.34291495991632
0.22742059158799 | 0.34113912678865 | 3.57030720633364
-0.01968628999525 | -0.01379539810330 | 4.88590055700420
0.00028294069148 | 0.00022566427890 | 6.37901581267627
$N_{\mathrm{f}}=6$ $\Re(a)$=4.42897206653722
10.43888617723020 | 0.00000000000000 | 0.00000000000000
-4.15230200875813 | 6.36698813144879 | 1.07811959900154
-1.09637578338379 | -2.73158403646281 | 2.16055514404309
0.58534033454693 | 0.16941238316725 | 3.25051246017243
-0.05722146054555 | 0.03586257692680 | 4.35355433640684
0.00111445947001 | -0.00356355734717 | 5.50798572653769
0.00000137005560 | 0.00007066173359 | 6.80680166786776
$N_{\mathrm{f}}=7$ $\Re(a)$=4.56174027039906
12.86098703034591 | 0.00000000000000 | 0.00000000000000
-4.90510587402478 | 8.42049068837062 | 1.01602163523968
-1.93579633844306 | -3.74443704209627 | 2.04133345120767
1.01179259239745 | 0.15201429359842 | 3.08051399458118
-0.09901358250763 | 0.09599744431866 | 4.11644818289646
-0.00338136152879 | -0.01174395914925 | 5.06467929848477
0.00103305167262 | 0.00015238534910 | 5.97856237235781
-0.00002200273877 | 0.00000584510008 | 7.20275534500376
$N_{\mathrm{f}}=8$ $\Re(a)$=4.55393944225591
11.79355692851662 | 0.00000000000000 | 0.00000000000000
-3.39975426063274 | 8.58842130087254 | 0.94528885879869
-3.19998047591390 | -3.10247031448467 | 1.89068584187872
1.24319452630443 | -0.50090970998006 | 2.83585047371902
-0.01592303613588 | 0.25412899494655 | 3.78108079372661
-0.02633000853185 | -0.01569161687470 | 4.72860249427096
0.00201070487884 | -0.00122269463134 | 5.67051703001310
0.00000598719210 | 0.00012240497422 | 6.61841395394285
-0.00000190141931 | -0.00000217573991 | 7.70094251177546

As the LR term can be expressed as
$V_{\mathrm{LR}}(\mathbf{r})=V_{\mathrm{C}}(\mathbf{r})-V_{\mathrm{SR}}(\mathbf{r})$,
we discuss in detail only the SR term (the erfc kernel). As shown in previous
studies, $V_{\mathrm{SR}}(\mathbf{r}-\mathbf{r}^{\prime})$ is not separable in
the same way as the Coulomb and Yukawa kernels, and the standard approach to
performing such a convolution involves a two-dimensional quadrature [28, 16].
In this work, we propose an alternative method that reduces the calculation to
evaluating a few one-dimensional integrals. The $\mathrm{erfc}$ function can
be efficiently represented in the following form:

$\mathrm{erfc}(r)\approx\sum_{j=-N_{\mathrm{f}}}^{N_{\mathrm{f}}}A_{j}e^{a_{j}r},$
(29)

where $A_{j}$ and $a_{j}$ are complex parameters determined in a fitting
procedure. To ensure that $\mathrm{erfc}(r)$ is a real function, we impose
that $A_{j}=A^{\ast}_{-j}$, $a_{j}=a^{\ast}_{-j}$ as well as $A_{0}$ and
$a_{0}$ are real. The obtained fitting parameters are given Tab. 1.

Figure 2:  Unsigned difference between complementary error function
($\mathrm{erfc}$) and its fit with $2N_{f}+1$ functions in Eq. 29.

Fig. 2 shows the quality of the fit depending on the number of the functions.
Increasing the number of the exponents in Eq. 29 reduces the error of the
representation. A linear combination of 17 functions ($N_{\mathrm{f}}=8$)
reproduces the $\mathrm{erfc}(r)$ values within $\sim 10^{-12}$ for all
allowed values of $r$ and within $\sim 10^{-13}$ for all $r>0.1$.

Using the fit, we express $V_{\mathrm{SR}}(\mathbf{r}-\mathbf{r}^{\prime})$ as
a linear combination of complex Yukawa potentials:

$V_{\mathrm{SR}}(\mathbf{r}-\mathbf{r}^{\prime})=\sum_{j=-N_{f}}^{N_{f}}\frac{A_{j}e^{a_{j}\mu|\mathbf{r}-\mathbf{r}^{\prime}|}}{|\mathbf{r}-\mathbf{r}^{\prime}|}.$
(30)

Furthermore, due to the restrictions imposed on $A_{j}$ and $a_{j}$, Eq. 30
reduces to

$V_{\mathrm{SR}}(\mathbf{r}-\mathbf{r}^{\prime})=\frac{A_{0}e^{a_{0}\mu|\mathbf{r}-\mathbf{r}^{\prime}|}}{|\mathbf{r}-\mathbf{r}^{\prime}|}+2\Re\left[\sum_{j=1}^{N_{f}}\frac{A_{j}e^{a_{j}\mu|\mathbf{r}-\mathbf{r}^{\prime}|}}{|\mathbf{r}-\mathbf{r}^{\prime}|}\right].$
(31)

A convolution with the $V_{\mathrm{SR}}(\mathbf{r}-\mathbf{r}^{\prime})$
reduces to a sum of $N_{\mathbf{f}}+1$ convolutions with
$V_{\mathrm{Y}}(\mathbf{r}-\mathbf{r}^{\prime})$ that has the separable form
shown above.

Figure 3:  Errors in LC-BLYP total energies due to the fitted erfc kernel in
the non-local exchange. The fit employs $2N_{\mathbf{f}}+1$ functions. The
energy reference was obtained with $N_{\mathrm{f}}=8$.

To test the introduced approach for the range-separated kernels, we perform a
total-energy calculation of the Ar and Rn atoms using the LC-BLYP exchange-
correlation functional [14]. Its non-local part of the exchange consists
entirely of the LR term. Fig. 3 shows the convergence of the total energy with
respect to $N_{\mathrm{f}}$. The errors reduce consistently with an increase
of the number of fitting functions, and $N_{\mathrm{f}}=5$ is sufficient to
guarantee the precision of $10^{-7}$ Ha for both considered atoms.

##  V Numerical integration and derivatives

All radial functions are represented on a predefined grid $r_{i}$ that spans
from $r_{\mathrm{min}}$ to $r_{\mathrm{max}}$, and the integrals that appear
in this method require a numerical evaluation. First, we consider an integral
of the kind $\int_{r_{i}}^{r_{i+1}}f(r)dr$. Assuming $f(r)$ is a smooth
function, we interpolate it using a $(2k-1)$th-degree Lagrange polynomial

$L(r)=\sum\limits_{j=i-k+1}^{i+k}f(r_{j})\ell_{j}(r),$ (32)

where

$\ell_{j}(r)=\prod_{\begin{subarray}{c}j^{\prime}\neq
j\end{subarray}}\frac{r-r_{j^{\prime}}}{r_{j}-r_{j^{\prime}}}.$ (33)

$L(r)$ is chosen such that $L(r_{j})=f(r_{j})$ for $j=i-k+1,\dots,i+k$. The
required integral is approximated then as $\int_{r_{i}}^{r_{i+1}}L(r)dr$.
Evaluating the coefficients of this polynomial to compute the integral is
impractical. Instead, we apply a Newton-Cotes formula with $Q+1$ points for
numerical integration using equally spaced abscissas [29] and obtain

$\int\limits_{r_{i}}^{r_{i+1}}L(r)dr=\frac{\delta
r_{i}}{Q}\sum\limits_{q=0}^{Q}W_{q}L\left(r_{i}+\frac{q}{Q}\delta
r_{i}\right),$ (34)

where $W_{q}$ are the integration weights and $\delta r_{i}=r_{i+1}-r_{i}$.
Note that the Newton-Cotes formulas yield an exact answer for a polynomial of
a sufficiently low degree. Combining Eqs. 32 and 34 yields

$\int\limits_{r_{i}}^{r_{i+1}}L(r)dr=\frac{\delta
r_{i}}{Q}\sum\limits_{j=i-k+1}^{i+k}\sum\limits_{q=0}^{Q}W_{q}\ell_{j}\left(r_{i}+\frac{q}{Q}\delta
r_{i}\right)f(r_{j}).$ (35)

The result can be written in the compact form

$\int\limits_{r_{i}}^{r_{i+1}}f(r)dr\approx\sum\limits_{j=i-k+1}^{i+k}w_{ij}f(r_{j}),$
(36)

where

$w_{ij}=\frac{\delta
r_{i}}{Q}\sum\limits_{q=0}^{Q}W_{q}\ell_{j}\left(r_{i}+\frac{q}{Q}\delta
r_{i}\right).$ (37)

The weights $w_{ij}$ depend only on a grid and have to be computed only once.
If the integration boundaries are $r_{\mathrm{min}}$ and $r_{\mathrm{max}}$,
the integral can be expressed in a similar form

$\int\limits_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}f(r)dr\approx\sum\limits_{i}\tilde{w}_{i}f(r_{i}).$
(38)

In our further calculations, we use the Bode’s rule in Eq. 34 and the Lagrange
polynomial with $p=9$. For $i=1,\dots,4$, it is not possible to select the
support points for interpolation the same way as suggested above. In this
case, we interpolate $f(r)$ using $j=0,\dots,9$ assuming that $f(r_{0})=0$,
and it allows us also to include an integral over the range $0<r<r_{1}$ in the
calculation. In comparison to neglecting this tiny region, taking it into
account makes it possible to define $r_{1}$ by 2-3 orders of magnitude larger.

Calculation of the kinetic energy requires a derivative, and we use the
Lagrange interpolation again. An approximation to the derivative is expressed
as

$f^{\prime}(r)=\sum_{j=0}^{k}f_{j}\ell^{(1)}_{j}(r)$ (39)

with the weights

$\ell^{(1)}_{j}(r)=\sum_{\begin{subarray}{c}i=0\\\ i\neq
j\end{subarray}}^{k}\Biggl{(}\frac{1}{r_{j}-r_{i}}\prod_{\begin{subarray}{c}m=0\\\
m\neq j\\\ m\neq i\end{subarray}}^{k}\frac{r-r_{m}}{r_{j}-r_{m}}\Biggr{)}.$
(40)

##  VI  Radial grids

Figure 4:  Error in total HF energy of Ar calculated using a radial grid with
$N$ points. The reference energy corresponds to a calculation with $N=2000$.
The grid types shown in the legend are defined in Eqs. 42 and 43. Figure 5:
Error in total HF energy of Rn calculated using a radial grid with $N$ points.
The reference energy corresponds to a calculation with $N=2000$. The grid
types shown in the legend are defined in Eqs. 42 and 43.

The logarithmic grid is the de facto standard choice electronic-structure
codes. It is defined as

$r_{i}=r_{\mathrm{min}}\left(\frac{r_{\mathrm{max}}}{r_{\mathrm{min}}}\right)^{\frac{i-1}{N-1}},$
(41)

where $r_{\mathrm{min}}$ and $r_{\mathrm{max}}$ correspond to the innermost
and the outermost points of the grid, respectively. $N$ is the number of the
grid points. It was argued in Ref. [30] that an inverse cubic grid is more
appropriate in atomic calculations than the equidistant or logarithmic ones.
This idea was tested in the LAPW exciting code [11] where the inverse cubic
grid is currently the default option.

In this work, we introduce order-$p$ polynomial grids with the following
definition:

$r_{i}=ir_{\mathrm{min}}+\left(\frac{i-1}{N-1}\right)^{p}(r_{\mathrm{max}}-Nr_{\mathrm{min}}).$
(42)

The linear term $ir_{\mathrm{min}}$ ensures that the distance between $r_{1}$
and $r_{2}$ is reasonable and roughly is equal to $r_{\mathrm{min}}$. Too
small $r_{2}-r_{1}$ leads to a numerical noise that appears due subtraction of
nearly equal numbers during the Lagrange interpolation. We use the same
approach to adjust the logarithmic grid and obtain

$r_{i}=(i-1)r_{\mathrm{min}}+r_{\mathrm{min}}\left(\frac{r_{\mathrm{max}}}{r_{\mathrm{min}}}N-1\right)^{\frac{i-1}{N-1}}.$
(43)

To assess the described grid types, we compare their performance. Figs. 4 and
5 show errors in the total energies of the Ar and Rn atoms with respect to
$N$. The inverse polynomial grids with $p=3,5$ and $7$, perform similarly for
Ar, but there is a much larger distinction between them for Rn where the
setting $p=7$ performs the best. In neither case, the logarithmic grid shows
the optimal performance.

##  VII Atomic energies

Table 2: Non-relativistic total energies (in Hartrees) of closed-shell atoms.  atom | HF | VWN | PBE | PBE0 | B3LYP | LC-BLYP
---|---|---|---|---|---|---
He | -2.861679996 | -2.834835624 | -2.892934867 | -2.895178376 | -2.915218663 | -2.866810561
Be | -14.573023168 | -14.447209474 | -14.629947716 | -14.636641425 | -14.673328176 | -14.584722714
Ne | -128.547098109 | -128.233481269 | -128.866427745 | -128.871759474 | -128.980973238 | -128.816627071
Mg | -199.614636425 | -199.139406315 | -199.955115169 | -199.970695270 | -200.103549936 | -199.907035649
Ar | -526.817512803 | -525.946194919 | -527.346128774 | -527.388217197 | -527.567834997 | -527.321256048
Ca | -676.758185925 | -675.742282614 | -677.348819102 | -677.392363604 | -677.595272840 | -677.329114053
Zn | -1777.848116191 | -1776.573849681 | -1779.182796711 | -1779.191450269 | -1779.503834064 | -1779.205663017
Kr | -2752.054977346 | -2750.147940421 | -2753.416108936 | -2753.512330119 | -2753.851525806 | -2753.494137871
Sr | -3131.545686439 | -3129.453161377 | -3132.948659228 | -3133.055470684 | -3133.414387383 | -3133.036729104
Pd | -4937.921024070 | -4935.368405699 | -4939.793447386 | -4939.923012608 | -4940.317844921 | -4939.893211695
Cd | -5465.133142530 | -5462.390982009 | -5467.051582940 | -5467.198445049 | -5467.607200062 | -5467.162476372
Xe | -7232.138363872 | -7228.856106486 | -7234.233211984 | -7234.433366465 | -7234.867433857 | -7234.417540215
Ba | -7883.543827330 | -7880.111578015 | -7885.731166873 | -7885.930820635 | -7886.384144735 | -7885.925529978
Yb | -13391.456193118 | -13388.048594318 | -13395.288895142 | -13395.345840849 | -13395.890072201 | -13395.505337325
Hg | -18408.991494945 | -18404.274219990 | -18412.743896725 | -18413.020005924 | -18413.551671176 | -18412.999474421
Rn | -21866.772240873 | -21861.346868935 | -21870.576884534 | -21870.947943763 | -21871.481295618 | -21870.907206607
Ra | -23094.303666425 | -23088.688083054 | -23098.174709109 | -23098.554874219 | -23099.101083655 | -23098.513007616

We apply the radial solver for calculating non-relativistic total energies of
atoms using the following approximations for exchange and correlation: (i) the
HF approximation, (ii) the LDA parametrization by Vosko, Wilk and Nusair
(VWN5) [31], (iii) the GGA parametrization by Perdew, Burke and Ernzerhof
(PBE) [32], (iv) the PBE0 hybrid [33], (v) the B3LYP hybrid [34] and (vi) the
LC-BLYP hybrid with the long-range part ($\mu=0.3$) of the exchange [14]. The
same grid parameters are chosen for all atoms in all calculations, namely,
$r_{\mathrm{min}}=10^{-6}$ bohr, $r_{\mathrm{max}}=30$ bohr $N=800$. Based on
the findings shown in Sec. VI, we employ the inverse polynomial grid with
$p=7$.

The total energies obtained for atoms with closed subshells are given in Tab.
2. The complete set of inputs and outputs (total energies, Kohn-Sham
eigenenergies and wavefunctions) from our non-relativistic spin-restricted
spherically-symmetric calculations for all atoms from H through U is available
in an open-access repository [35]. With the chosen parameters, our HF data
agree with the results published in Ref. [15] to all signs given in Tab. 2.
Moreover, for all these elements, we find a 14-digit agreement approaching the
double-precision limit in the floating-point representation.

Our VWN5 energies agree perfectly with the data in Ref. [16], where these
calculations we performed using a finite-element basis. The total energies
were given with signs up to 1 $\mu$Ha, and our results are consistent with all
given signs in that work. In comparison to another study [36] with VWN5
calculations, we find discrepancies up to 3 $\mu$Ha in the case of the Rn
atom.

The agreement between our HF energies and those in Ref. [15] as well as the
convergence tests performed in Sec. VI make us confident that all digits given
in the VWN5 results are significant. The remaining four exchange-correlation
functionals employ the GGA, and, in this case, it is common to introduce
density thresholds below which $v_{\mathrm{xc}}^{\mathrm{L}}(\mathbf{r})$ is
not computed. Such a threshold is also introduced in the libxc library and is
the main precision-limiting factor in these calculations. Nevertheless, we
anticipate that this feature does not lead to errors exceeding a few nHa.

Our PBE and PBE0 energies for He, Be, Ne, Mg and Ar agree with the multi-
resolution analysis data well within 1 $\mu$Ha. Finally, also the LC-BLYP
total energies for He, Be, Ne and Mg agree perfectly (within all signs given
in the reference) with the data from the finite-element calculation reported
in Ref. [16], and only for Ar with find a difference of 1 $\mu$Ha. The
obtained agreement in the LC-BLYP data is noteworthy, because the two methods
employed different approaches for calculating the integrals with the long-
range part of the exchange.

##  VIII Conclusions

We implemented a highly precise tool for solving the Kohn-Sham problem for
spherically-symmetric atoms. This tool solves the Schrödinger equation as an
integral equation rather than a differential one and heavily relies one-
dimensional convolutions. For this reason, we implemented efficient
quadratures based on Lagrange interpolation and investigated different grid
types. The frequently used logarithmic grids are far from being the optimal
choice. The proposed alternative, inverse polynomial grids, are superior to
the logarithmic ones (especially with the polynomial degree 7) and allow for a
significant reduction of the number of points. These technical tricks are
potentially useful in electronic-structure codes, especially, those employing
linearized augmented plane waves or linearized muffin-tin orbitals.

To support hybrid functionals with the screened or long-range exchange in the
presented radial solver, we introduced a new method for calculating
convolutions with the complementary error function kernel. This method
presents an alternative to the standard approach that applies bi-variate
quadrature. Our approach can be transferred to full-potential all-electron
codes where it will serve as a reference method for hybrid functionals with
screened exchange.

Finally, we applied the radial solver in calculations of closed-shell atoms
and observed a remarkable agreement with highly precise Hartree-Fock data in
literature. This result implies that our further calculations performed with
five local and hybrid exchange-correlation functionals set a reliable
benchmark.

###### Acknowledgements.

This work was funded by the Latvian Research Council via the project Precise
Fock Exchange (PREFEX) with the grant agreement No. lzp-2020/2-0251.

## References

  * Troullier and Martins [1991] N. Troullier and J. L. Martins, Efficient pseudopotentials for plane-wave calculations, Phys. Rev. B 43, 1993 (1991).
  * Blum  _et al._ [2009] V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter, and M. Scheffler, Ab initio molecular simulations with numeric atom-centered orbitals, Computer Physics Communications 180, 2175 (2009).
  * Slater [1937] J. C. Slater, Wave functions in a periodic potential, Phys. Rev. 51, 846 (1937).
  * Andersen [1975] O. K. Andersen, Linear methods in band theory, Phys. Rev. B 12, 3060 (1975).
  * Cur [2012] Aflow: An automatic framework for high-throughput materials discovery, Computational Materials Science 58, 218 (2012).
  * Jain  _et al._ [2013] A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, and K. A. Persson, Commentary: The materials project: A materials genome approach to accelerating materials innovation, APL Materials 1, 011002 (2013).
  * Kirklin  _et al._ [2015] S. Kirklin, J. E. Saal, B. Meredig, A. Thompson, J. W. Doak, M. Aykol, S. Rühl, and C. Wolverton, The open quantum materials database (oqmd): assessing the accuracy of dft formation energies, npj Computational Materials 1, 15010 (2015).
  * Noumerov [1924] B. V. Noumerov, A Method of Extrapolation of Perturbations, Monthly Notices of the Royal Astronomical Society 84, 592 (1924), https://academic.oup.com/mnras/article-pdf/84/8/592/3661174/mnras84-0592.pdf .
  * Bul [1966] Numerical treatment of ordinary differential equations by extrapolation methods, Numerische Mathematik 8, 1 (1966).
  * [10] The Elk Code, http://elk.sourceforge.net/.
  * Gulans  _et al._ [2014] A. Gulans, S. Kontur, C. Meisenbichler, D. Nabok, P. Pavone, S. Rigamonti, S. Sagmeister, U. Werner, and C. Draxl, exciting: a full-potential all-electron package implementing density-functional theory and many-body perturbation theory, Journal of Physics: Condensed Matter 26, 363202 (2014).
  * Peverati and Truhlar [2014] R. Peverati and D. G. Truhlar, Quest for a universal density functional: the accuracy of density functionals across a broad spectrum of databases in chemistry and physics, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 372, 20120476 (2014).
  * Heyd  _et al._ [2003] J. Heyd, G. E. Scuseria, and M. Ernzerhof, Hybrid functionals based on a screened coulomb potential, The Journal of Chemical Physics 118, 8207 (2003), https://doi.org/10.1063/1.1564060 .
  * Tawada  _et al._ [2004] Y. Tawada, T. Tsuneda, S. Yanagisawa, T. Yanai, and K. Hirao, A long-range-corrected time-dependent density functional theory, The Journal of Chemical Physics 120, 8425 (2004), https://doi.org/10.1063/1.1688752 .
  * Cinal [2020] M. Cinal, Highly accurate numerical solution of hartree–fock equation with pseudospectral method for closed-shell atoms, Journal of Mathematical Chemistry 58, 1571 (2020).
  * Lehtola [2020] S. Lehtola, Fully numerical calculations on atoms with fractional occupations and range-separated exchange functionals, Phys. Rev. A 101, 012516 (2020).
  * Kalos [1962] M. H. Kalos, Monte carlo calculations of the ground state of three- and four-body nuclei, Physical Review 128, 1791 (1962).
  * Harrison  _et al._ [2004] R. J. Harrison, G. I. Fann, T. Yanai, Z. Gan, and G. Beylkin, Multiresolution quantum chemistry: Basic theory and initial applications, The Journal of Chemical Physics 121, 11587 (2004).
  * Bischoff and Valeev [2011] F. A. Bischoff and E. F. Valeev, Low-order tensor approximations for electronic wave functions: Hartree–fock method with guaranteed precision, The Journal of Chemical Physics 134, 104104 (2011).
  * Bischoff  _et al._ [2012] F. A. Bischoff, R. J. Harrison, and E. F. Valeev, Computing many-body wave functions with guaranteed precision: The first-order møller-plesset wave function for the ground state of helium atom, The Journal of Chemical Physics 137, 104103 (2012).
  * Jensen  _et al._ [2017] S. R. Jensen, S. Saha, J. A. Flores-Livas, W. Huhn, V. Blum, S. Goedecker, and L. Frediani, The elephant in the room of density functional theory calculations, The Journal of Physical Chemistry Letters 8, 1449 (2017), pMID: 28291362, https://doi.org/10.1021/acs.jpclett.7b00255 .
  * Ratcliff  _et al._ [2020] L. E. Ratcliff, W. Dawson, G. Fisicaro, D. Caliste, S. Mohr, A. Degomme, B. Videau, V. Cristiglio, M. Stella, M. D’Alessandro, S. Goedecker, T. Nakajima, T. Deutsch, and L. Genovese, Flexibilities of wavelets as a computational basis set for large-scale electronic structure calculations, The Journal of Chemical Physics 152, 194110 (2020), https://doi.org/10.1063/5.0004792 .
  * Hu  _et al._ [2000] X.-G. Hu, T.-S. Ho, and H. Rabitz, Solving the bound-state schrödinger equation by reproducing kernel interpolation, Physical Review E 61, 2074 (2000).
  * Solala  _et al._ [2017] E. Solala, S. A. Losilla, D. Sundholm, W. Xu, and P. Parkkinen, Optimization of numerical orbitals using the helmholtz kernel, The Journal of Chemical Physics 146, 084102 (2017).
  * Lehtola  _et al._ [2018] S. Lehtola, C. Steigemann, M. J. Oliveira, and M. A. Marques, Recent developments in libxc — a comprehensive library of functionals for density functional theory, SoftwareX 7, 1 (2018).
  * Woods  _et al._ [2019] N. D. Woods, M. C. Payne, and P. J. Hasnip, Computing the self-consistent field in kohn–sham density functional theory, Journal of Physics: Condensed Matter 31, 453001 (2019).
  * Johnson [2007] W. R. Johnson,  _Atomic Structure Theory_ (Springer Berlin Heidelberg, 2007).
  * Ángyán  _et al._ [2006] J. G. Ángyán, I. Gerber, and M. Marsman, Spherical harmonic expansion of short-range screened coulomb interactions, Journal of Physics A: Mathematical and General 39, 8613 (2006).
  * Abramowitz [1965] M. Abramowitz,  _Handbook of mathematical functions, with formulas, graphs, and mathematical tables_ (Dover Publications, New York, 1965).
  * Havlová and Smrčka [1984] H. Havlová and L. Smrčka, Numerical solution of the radial schrödinger equation in the inverse cubic scale, Czechoslovak Journal of Physics 34, 961 (1984).
  * Vosko  _et al._ [1980] S. H. Vosko, L. Wilk, and M. Nusair, Accurate spin-dependent electron liquid correlation energies for local spin density calculations: a critical analysis, Canadian Journal of Physics 58, 1200 (1980), https://doi.org/10.1139/p80-159 .
  * Perdew  _et al._ [1996a] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996a).
  * Adamo and Barone [1999] C. Adamo and V. Barone, Toward reliable density functional methods without adjustable parameters: The pbe0 model, The Journal of Chemical Physics 110, 6158 (1999), https://doi.org/10.1063/1.478522 .
  * Stephens  _et al._ [1994] P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J. Frisch, Ab initio calculation of vibrational absorption and circular dichroism spectra using density functional force fields, The Journal of Physical Chemistry 98, 11623 (1994), https://doi.org/10.1021/j100096a001 .
  * Uzulis and Gulans [2022] J. Uzulis and A. Gulans, Total energies of atoms from integral-equation radial solver (dataset under review, preliminary doi given), https://doi.org/10.24435/materialscloud:2w-ev (2022).
  * Kraisler  _et al._ [2010] E. Kraisler, G. Makov, and I. Kelson, Ensemble $v$-representable ab initio density-functional calculation of energy and spin in atoms: A test of exchange-correlation approximations, Phys. Rev. A 82, 042516 (2010).