# Converging tetrahedron method calculations for the nondissipative parts of
spectral functions

Minsu Ghim<EMAIL_ADDRESS>Center for Correlated Electron
Systems, Institute for Basic Science, Seoul 08826, Korea  Department of
Physics and Astronomy, Seoul National University, Seoul 08826, Korea  Center
for Theoretical Physics, Seoul National University, Seoul 08826, Korea  Cheol-
Hwan Park<EMAIL_ADDRESS>Center for Correlated Electron Systems,
Institute for Basic Science, Seoul 08826, Korea  Department of Physics and
Astronomy, Seoul National University, Seoul 08826, Korea  Center for
Theoretical Physics, Seoul National University, Seoul 08826, Korea  Donostia
International Physics Center, 20018 San Sebastián, Spain  Centro de Física de
Materiales, Universidad del País Vasco, EHU, 20018 San Sebastián, Spain

###### Abstract

Many physical quantities in solid-state physics are calculated from $k$-space
summation. For spectral functions, the frequency-dependent factor can be
decomposed into the energy-conserving delta function part and the
nondissipative principal value part. A very useful scheme for this $k$-space
summation is the tetrahedron method. Tetrahedron methods have been widely used
to calculate the summation of the energy-conserving delta function part such
as the imaginary part of the dielectric function. On the other hand, the
corresponding tetrahedron method for the nondissipative part such as the real
part of the dielectric function has not been used much. In this paper, we
address the technical difficulties in the tetrahedron method for the
nondissipative part and present an easy-to-implement, stable method to
overcome those difficulties. We demonstrate our method by calculating the
static and dynamical spin Hall conductivity of platinum. Our method can be
widely applied to calculate linear static or dynamical conductivity, self-
energy of an electron, and electric polarizability, to name a few.

##  I INTRODUCTION

In the field of electronic structure calculations, many physical quantities of
a periodic system are obtained from a Brillouin zone (BZ) integral of the
following form:

$\sum_{n,m\neq
n}\int_{\textrm{BZ}}\frac{d^{3}k}{(2\pi)^{3}}f_{{n\mathbf{k}}}\frac{F_{nm\textbf{k}}}{\hbar\omega-\left(\epsilon_{m{\bf
k}}-\epsilon_{n{\bf k}}\right)+i\eta}\,.$ (1)

Here $f_{{n\mathbf{k}}}$ and $\epsilon_{n{\bf k}}$ are the Fermi-Dirac
occupation factor and energy eigenvalue for an electronic state with band
index $n$ and Bloch wavevector k, respectively, $F_{nm\textbf{k}}$, a complex
function of two band indices $n$ and $m$ and Bloch wavevector k, is a product
of proper matrix elements, and $\eta=0^{+}$ is an infinitesimal positive
number. Physical quantities such as correlation functions from linear response
theory (conductivity, susceptibility, dielectric functions, etc .) or any
integration of Green’s functions in $k$ space (e. g. , self-energy) can be
obtained by calculating this integral.

The Brillouin-zone integration has been mainly performed from two different
schemes: the smearing method and the tetrahedron method. The smearing method
assumes that $\eta$ is small but finite. The integral can be replaced by a sum
over a special set of $k$ points. Next, the tetrahedron method is the
analytical integration method assuming that $F_{nm\textbf{k}}$,
$\epsilon_{m{\bf k}}$, and $\epsilon_{n{\bf k}}$ vary linearly inside a
tetrahedron in $k$ space. Six tetrahedra are constructed for each
parallelepiped in $k$ space whose vertices belong to the regular $k$-point
grid: described in Appendix A. In a single tetrahedron, the matrix elements
and band energies are linearized.

The integral gives two terms, the principal value part and the delta function
part, thanks to the following decomposition:

$\begin{split}\frac{1}{\hbar\omega-\left(\epsilon_{m{\bf k}}-\epsilon_{n{\bf
k}}\right)+i\eta}&\,=\\\ {\rm
P}\,\frac{1}{\hbar\omega-\left(\epsilon_{{m\mathbf{k}}}-\epsilon_{{n\mathbf{k}}}\right)}&-i\pi\delta\left(\hbar\omega-(\epsilon_{{m\mathbf{k}}}-\epsilon_{{n\mathbf{k}}})\right)\,.\end{split}$
(2)

While being used to calculate the energy-conserving delta function part [1],
the tetrahedron method has not been used much for the nondissipative principal
value part. In principle, the nondissipative part can be obtained from the
delta function part, since the two are connected by the Kramers-Kronig
relation. However, the Kramers-Kronig transform may require a significant
amount of computational time: The Kramers-Kronig relation is written as an
integral with respect to the frequency from zero to infinity. It is thus
required to calculate the delta function part at a dense frequency grid
covering a wide frequency range, even if a single frequency component of the
nondissipative part is needed.

Previous studies suggested the tetrahedron method for the nondissipative part
[2, 3]. The method can also be applied to metals with a proper division of
tetrahedra at the Fermi surface at zero temperature [4]. However, this method
has some technical difficulties. First, the explicit formula for the
integration is in seven different complicated forms, depending on the energy
eigenvalues at the four vertices of a tetrahedron. In more recent studies [5,
6], the multiple-case method was similarly adopted. Second, the criteria to
distinguish between these cases are not well-defined numerically: They depend
on whether some variables are exactly zero or not. Third, most importantly,
the round-off error occurring during the evaluation of the logarithmic
functions is not controlled. To implement the tetrahedron method for the
nondissipative part, these difficulties should be handled appropriately.

In this paper, we present an easy-to-implement, numerically stable tetrahedron
method for the nondissipative part to overcome these technical difficulties.
After a detailed explanation of the method, we compare the convergence and the
computational time of the tetrahedron method with those of the adaptive
smearing technique [7] for the intrinsic static and dynamical spin Hall
conductivity (SHC) of fcc platinum using the interpolation scheme based on
maximally locallized Wannier functions (MLWFs) [8, 9]. We show that, to reach
the same level of convergence, our improved tetrahedron method requires orders
of magnitude shorter computation times than the adaptive smearing method.
Finally, we discuss how the convergence can be achieved efficiently.

##  II METHOD

###  II.1 Tetrahedron Method for the Kubo Formula

Consider the following three types of integrals:

$\displaystyle I_{1,nm}(\hbar\omega)$
$\displaystyle=\int_{\textrm{tet.}}d^{3}k\,{\rm
P}\,\frac{F_{nm\textbf{k}}}{\Delta_{nm\textbf{k}}+\hbar\omega}\,,$ (3)
$\displaystyle I_{2,nm}(\hbar\omega)$
$\displaystyle=\int_{\textrm{tet.,}\Delta_{nm\textbf{k}}=\hbar\omega}d^{2}S_{k}\frac{F_{nm\textbf{k}}}{|\nabla_{\textbf{k}}\Delta_{nm\textbf{k}}|}\,,$
(4) $\displaystyle I_{3,nm}$
$\displaystyle=\int_{\textrm{tet.}}d^{3}k\frac{F_{nm\textbf{k}}}{\Delta_{nm\textbf{k}}^{2}}\,,$
(5)

where the integration in momentum space is performed over a particular
tetrahedron, and
$\Delta_{nm\textbf{k}}=\epsilon_{{n\mathbf{k}}}-\epsilon_{{m\mathbf{k}}}$.
Especially, $I_{1,nm}$ and $I_{3,nm}$ is obtained by evaluating the
nondissipative parts. The formulas to calculate Green’s-function-like or
susceptibility-like quantities such as anomalous hall conductivity [10, 11],
Fan-Migdal self-energy [12, 13, 14], and spin Hall conductivity [15, 16], can
be transformed and decomposed to terms of the forms in Eqs. (3) – (5). To
allow analytical integration in Eqs. (3), (4), and (5), we used the Padé
approximant, the best approximation by a rational function. The Padé
approximant is obtained by using

$\displaystyle F_{nm\textbf{k}}=$
$\displaystyle\,F_{0}+F_{x}k_{x}+F_{y}k_{y}+F_{z}k_{z}$ (6)
$\displaystyle\Delta_{nm\mathbf{k}}=$
$\displaystyle\,\Delta_{0}+\Delta_{x}k_{x}+\Delta_{y}k_{y}+\Delta_{z}k_{z}$
(7)

where $F_{0},\,\cdots,F_{z}$, and $\Delta_{0},\,\cdots,\Delta_{z}$ are the
parameters reproducing the values at each of the four vertices of a
tetrahedron. As a result,

$\displaystyle I_{1,nm}(\hbar\omega)$
$\displaystyle=\int_{\textrm{tet.}}d^{3}k\,{\rm
P}\,\frac{F_{0}+F_{x}k_{x}+F_{y}k_{y}+F_{z}k_{z}}{\Delta_{0}+\hbar\omega+\Delta_{x}k_{x}+\Delta_{y}k_{y}+\Delta_{z}k_{z}}$
(8) $\displaystyle I_{3,nm}$
$\displaystyle=\int_{\textrm{tet.}}d^{3}k\frac{F_{0}+F_{x}k_{x}+F_{y}k_{y}+F_{z}k_{z}}{[\Delta_{0}+\Delta_{x}k_{x}+\Delta_{y}k_{y}+\Delta_{z}k_{z}]^{2}}\,,$
(9)

and these can be integrated by hand since we know that the integration of
rational functions contains other rational functions with logarithms. For
example, the one-dimensional case for Eq. (8) is equivalent to

$\displaystyle\begin{split}\int_{k_{i}}^{k_{f}}dk_{x}&{\rm
P}\,\frac{F_{0}+F_{x}k_{x}}{\Delta_{0}+\hbar\omega+\Delta_{x}k_{x}}\\\
=&\,\frac{F_{0}\Delta_{x}-F_{x}(\Delta_{0}+\hbar\omega)}{\Delta_{x}^{2}}{\rm
log}\left|\frac{\Delta_{x}k_{f}+\Delta_{0}+\hbar\omega}{\Delta_{x}k_{i}+\Delta_{0}+\hbar\omega}\right|\\\
&+\frac{F_{x}}{\Delta_{x}}(k_{f}-k_{i})\,.\end{split}$ (10)

The closed-form expressions for the three-dimensional case are presented in
Appendix B. We note that the equivalent results for $I_{1,nm}$ [Eq. (3)] are
in Ref. [3] and those for $I_{2,nm}$ [Eq. (4)] are in Ref. [1]. The results
for $I_{3,nm}$ [Eq. (5)] are presented here for the first time.

A specific component of the intrinsic spin Hall conductivity, for example, is
expressed as [15, 16]

$\displaystyle\begin{split}\sigma^{z}_{xy}(\omega)={}\frac{e}{\hbar}&\frac{\Omega}{(2\pi)^{3}}\sum_{n,m\neq
n}\int_{\textrm{BZ}}d^{3}k(f_{{n\mathbf{k}}}-f_{{m\mathbf{k}}})\times\\\
&\frac{{\rm
Im}[\matrixelement{u_{{n\mathbf{k}}}}{\hat{j}^{z}_{x\textbf{k}}}{u_{{m\mathbf{k}}}}\matrixelement{u_{{m\mathbf{k}}}}{\hat{v}_{y\textbf{k}}}{u_{{n\mathbf{k}}}}]}{\Delta_{nm\textbf{k}}^{2}-(\hbar\omega+i\eta)^{2}}\end{split}$
(11)

where $\Omega$ is the volume of a unit cell and $\ket{u_{{n\mathbf{k}}}}$ is
the periodic part of the Bloch state satisfying
$H_{\mathbf{k}}\ket{u_{{n\mathbf{k}}}}=\epsilon_{{n\mathbf{k}}}\ket{u_{{n\mathbf{k}}}}$
where
$H_{\mathbf{k}}=e^{-i\mathbf{k}\cdot\mathbf{r}}He^{i\mathbf{k}\cdot\mathbf{r}}$,
and $H$ is a lattice-periodic Hamiltonian.
$\hat{j}^{z}_{x\textbf{k}}=\frac{\hbar}{4}\\{\sigma^{z},v_{x\textbf{k}}\\}$,
and $\hat{v}_{y\textbf{k}}=\frac{1}{\hbar}\partial_{y}H_{\textbf{k}}$ are the
spin current and velocity operators, respectively. After dividing the
Brillouin zone according to the scheme presented in Appendix A, we can
decompose Eq. (11) into integrations over a single tetrahedron of the forms in
Eqs. (3)–(5). We obtain the dynamical conductivity of the form

$\displaystyle\begin{split}&\frac{1}{2\hbar\omega}\sum_{n,m\neq
n}\left[\left\\{I_{1,nm}(-\hbar\omega)-I_{1,nm}(\hbar\omega)\right\\}\right.\\\
&\left.+i\pi\left\\{I_{2,nm}(-\hbar\omega)+I_{2,nm}(\hbar\omega)\right\\}\right]\end{split}$
(12)

and the static conductivity of the form

$\displaystyle\sum_{n,m\neq n}I_{3,nm}\,.$ (13)

###  II.2 Numerical Problem and Its Solution

If we perform the integrations in Eqs. (3) and (5) analytically, we obtain the
following results:

$\displaystyle\begin{split}I_{1,nm}(\hbar\omega)={}&-\frac{(1+x_{1})(1+x_{2})(1+x_{3})}{6x_{1}^{2}x_{2}^{2}x_{3}^{2}(x_{1}-x_{2})^{2}(x_{2}-x_{3})^{2}(x_{3}-x_{1})^{2}}\\\
&\times\frac{\textrm{det}(\textbf{t})}{\Delta_{4}+\hbar\omega}\sum_{i=1}^{4}F_{i}\left(\sum_{j=1}^{3}C_{ij}^{(1)}\xi_{j}+B_{i}^{(1)}\right)\end{split}$
(14)
$\displaystyle\begin{split}I_{3,nm}={}&\frac{(1+x_{1})(1+x_{2})(1+x_{3})}{2x_{1}^{2}x_{2}^{2}x_{3}^{2}(x_{1}-x_{2})^{2}(x_{2}-x_{3})^{2}(x_{3}-x_{1})^{2}}\\\
&\times\frac{\textrm{det}(\textbf{t})}{\Delta_{4}^{2}}\sum_{i=1}^{4}F_{i}\left(\sum_{j=1}^{3}C_{ij}^{(3)}\xi_{j}+B_{i}^{(3)}\right)\,,\end{split}$
(15)

where $F_{i}$ and $\Delta_{i}$ ($i\in\\{1,2,3,4\\}$) are the values of
$F_{nm\textbf{k}}$ and $\Delta_{nm\textbf{k}}$ at the four vertices of a
single tetrahedron, respectively. $x_{j}$ and $\xi_{j}$ ($j\in\\{1,2,3\\}$)
are defined as:

$\displaystyle x_{j}=$
$\displaystyle\frac{\Delta_{4}-\Delta_{j}}{\Delta_{j}+\hbar\omega},\,\xi_{j}=\textrm{log}|1+x_{j}|\,.$
(16)

The expressions for $C_{ij}$ and $B_{i}$ in terms of $x_{i}$’s are given in
Appendix B.

However, there remains a numerical problem of calculating a very small number
divided by another very small number with floating-point arithmetic, which
occurs when the values of $\Delta_{nm\textbf{k}}$ of two or more vertices are
very close to each other. According to the definition above, $x_{1}$, $x_{2}$,
$x_{3}$, $x_{1}-x_{2}$, $x_{2}-x_{3}$, and $x_{3}-x_{1}$ are measures of how
close $\Delta_{i}$ and $\Delta_{j}$ are for pairs $(i,j)$ =
$(1,4),(2,4),(3,4),(1,2),(2,3)$, and $(3,1)$, respectively. As two or more
$\Delta_{i}$’s get closer, the elements in a subset
$X\subseteq\\{x_{1},x_{2},x_{3},x_{1}-x_{2},x_{2}-x_{3},x_{3}-x_{1}\\}$,
become closer to zero. Then, let us define $I_{1,nm}|_{X\rightarrow 0}$ and
$I_{3,nm}|_{X\rightarrow 0}$ as the limits of $I_{1,nm}$ and $I_{3,nm}$ if all
the elements in $X$ approach zero. (The expressions for
$I_{1,nm}|_{X\rightarrow 0}$ are presented in Ref. [3]). However, errors are
generated in evaluating the logarithmic functions because of rounding of
floating-point arithmetic and are amplified due to the small denominator in
the prefactors of Eqs. (14) and (15), resulting in non-convergence.

Let us define $\epsilon$ as the maximum value which a computer cannot store
when it is added to $1$ due to rounding of floating-point arithmetic. Roughly
speaking, $\epsilon$ $\approx$ $10^{-8}$, $10^{-16}$, and $10^{-32}$ for
single-, double-, and quadruple-precision data types, respectively.
Explicitly, we may write this as $1\doteq 1+O(\epsilon)$, where $\doteq$ means
“equality in machines.” For example, if a machine cannot distinguish two real
numbers $A$ and $B$, $A\doteq B$. For a real number $x$, $x\doteq
x+O(x\epsilon)$.

A round-off error comes from calculating
$\xi_{j}=\textrm{log}|1+x_{j}|=O(x_{j})$ (in case $|x_{j}|\ll 1$) or
$\xi_{k}=\textrm{log}|1+x_{k}|=\textrm{log}|1+x_{j}+(x_{k}-x_{j})|=\xi_{j}+O(x_{j}-x_{k})$
(in case $|x_{j}-x_{k}|\ll|1+x_{j}|$). Since the default computer algorithm to
compute ${\rm log}$ uses the Taylor expansion, ${\rm log}\,(1+x)$ is
calculated by $(1+x-1)-(1+x-1)^{2}/2+(1+x-1)^{3}/3+\cdots$, not by
$x-x^{2}/2+x^{3}/3+\cdots$. The value of $1+x-1$ is not the same as that of
$x$, but $1+x-1\doteq x+O(\epsilon)$ if $|x|\ll 1$. Hence,

$\textrm{log}|1+x|\doteq O(x)+O(\epsilon)$ (17)

for small $|x|$ ($\ll 1$). For a similar reason,
log$|1+x|=\textrm{log}|1+y+(x-y)|\doteq\textrm{log}|1+y|+O(x-y)+O(\epsilon)$
for small $|x-y|$ ($\ll 1$).

This round-off error can be problematic if $x_{j}\rightarrow 0$ or
$x_{j}-x_{k}\rightarrow 0$. In the first case, a term such as

$\displaystyle\begin{split}&\frac{\sum_{l}C_{il}\xi_{l}+B_{i}}{x_{1}^{2}x_{2}^{2}x_{3}^{2}(x_{1}-x_{2})^{2}(x_{2}-x_{3})^{2}(x_{3}-x_{1})^{2}}\qquad\quad\\\
&\doteq O(1)\times\frac{\sum_{l\neq
j}C_{il}\xi_{l}+C_{ij}(O(x_{j})+O(\epsilon))+B_{i}}{x_{j}^{2}}\\\
&=O(1)\times\frac{(\sum_{l\neq
j}C_{il}\xi_{l}+C_{ij}O(x_{j})+B_{i})+C_{ij}O(\epsilon)}{x_{j}^{2}}\\\
&=\frac{O(x_{j}^{2})+O(\epsilon)}{x_{j}^{2}}\end{split}$ (18)

should be finite as $x_{j}\rightarrow 0$, but is divergent due to
$O(\epsilon)$. In the last equality, we have used the fact that the
coefficients of $O(x_{j})$ should vanish in order to have a non-divergent
value in the limit $x_{j}\to 0$ if $\epsilon=0$. We can also check this
equality using Eqs. (B7) and (B8) in Appendix B. In the second case,

$\displaystyle\begin{split}&\frac{\sum_{l}C_{il}\xi_{l}+B_{i}}{x_{1}^{2}x_{2}^{2}x_{3}^{2}(x_{1}-x_{2})^{2}(x_{2}-x_{3})^{2}(x_{3}-x_{1})^{2}}\qquad\quad\\\
&\doteq O(1)\times\Bigg{[}\frac{\sum_{l\neq
k}C_{il}\xi_{l}+B_{i}}{(x_{j}-x_{k})^{2}}\\\
&\qquad\qquad\quad+\frac{C_{ik}(\xi_{j}+O(x_{j}-x_{k})+O(\epsilon))}{(x_{j}-x_{k})^{2}}\Bigg{]}\\\
&=O(1)\times\Bigg{[}\frac{\sum_{l\neq
k}C_{il}\xi_{l}+C_{ik}(\xi_{j}+O(x_{j}-x_{k}))+B_{i}}{(x_{j}-x_{k})^{2}}\\\
&\qquad\qquad\quad+\frac{C_{ik}O(\epsilon)}{(x_{j}-x_{k})^{2}}\Bigg{]}\\\
&=\frac{O((x_{j}-x_{k})^{2})+O(\epsilon)}{(x_{j}-x_{k})^{2}}\end{split}$ (19)

should also be finite as $x_{j}-x_{k}\rightarrow 0$, but is divergent due to
$O(\epsilon)$. In the last equality, the coefficients of $O(x_{j}-x_{k})$
should vanish similarly as Eq. (18). The same phenomenon happens when two or
more $x_{j},x_{j}-x_{k}$ are close to zero.

It is also well known that log1P is a more accurate way to compute
$\textrm{log}(1+x)$ in case $|x|\ll 1$ [17], but the round-off error does not
completely disappear here either. It is expressed as
$\textrm{log1P}(x)=x\,\textrm{log}(1+x)\,/\,(1+x-1)$. The error is
$O(x\epsilon)$, not $O(\epsilon)$, because

$\displaystyle\begin{split}\textrm{log1P}(x)=&\,x\,\frac{\textrm{log}(1+x)}{1+x-1}\\\
=&\,x\,\frac{(1+x-1)-(1+x-1)^{2}/2+\cdots}{1+x-1}\\\
=&\,x\,\left[1-(1+x-1)/2+\cdots\right]\\\ \doteq&\,x\,(O(1)+O(\epsilon))\\\
=&\,O(x)+O(x\epsilon)\,.\end{split}$ (20)

Therefore, the numerical problem is not yet in general resolved even with
log1P function because of the remaining error term.

To avoid this difficulty, we may first consider using the formulas for
$I_{1,nm}|_{X\rightarrow 0}$ in Ref. [3] and finding and using the formulas
for $I_{3,nm}|_{X\rightarrow 0}$. However, this method has two problems in
practice. First, since there are six parameters in total, many different
limiting cases should be considered depending on which of the six parameters
are much smaller in magnitude than one, making the implementation complicated.
Second, since they are real numbers which are not exactly zero, it is not easy
to define the criteria to determine whether one of the limiting-case formulas
should be applied or not.

By introducing a criterion for determining whether the elements of $X$ are
close to zero, we can solve the numerical problem in a simple way, without
implementing many cases. Let us redefine $x_{1}$, $x_{2}$, $x_{3}$,
$x_{1}-x_{2}$, $x_{2}-x_{3}$, and $x_{3}-x_{1}$ symmetrically:

$x_{ij}=\frac{\Delta_{i}-\Delta_{j}}{(\Delta_{i}+\Delta_{j})/2+\hbar\omega}\,.$
(21)

If one of $x_{ij}$ gets close to zero, the numerical integration will have the
error shown in Eqs. (18) and (19). In order to prevent this error, we enforce
$|x_{ij}|$ to be not smaller than a specified value $\epsilon_{\textrm{tet.}}$
(cutoff) by modifying $\Delta_{i}$ and $\Delta_{j}$ as

$\Delta_{i}\rightarrow\frac{\Delta_{i}+\Delta_{j}}{2}\pm
0.5\left|\frac{\Delta_{i}+\Delta_{j}}{2}+\hbar\omega\right|\epsilon_{\textrm{tet.}}$
(22)

and

$\Delta_{j}\rightarrow\frac{\Delta_{i}+\Delta_{j}}{2}\mp
0.5\left|\frac{\Delta_{i}+\Delta_{j}}{2}+\hbar\omega\right|\epsilon_{\textrm{tet.}}\,,$
(23)

where the upper signs are selected if $\Delta_{i}\geq\Delta_{j}$ and the lower
signs are selected if $\Delta_{i}<\Delta_{j}$. Accordingly, the updated
$|x_{ij}|=\epsilon_{\textrm{tet.}}$, and the round-off error is not amplified
by the small denominator [see Eqs. (18) and (19)].

The cutoff should make the error term in the denominator $O(\epsilon)$ small
enough compared to the true value term, $O(x_{j}^{2})$ in Eq. (18) or
$O((x_{j}-x_{k})^{2})$ in Eq. (19). $\epsilon_{\textrm{tet.}}$ prevents the
error from becoming comparable to the true value with the condition,

$O(\epsilon)\ll O(\epsilon_{\textrm{tet.}}^{2})\,.$ (24)

Therefore, $\epsilon_{\textrm{tet.}}$ depends on the precision, or,
equivalently, the data type. If a higher-precision data type is used, the
required $\epsilon_{\textrm{tet.}}$ is smaller. We adopted
$\epsilon_{\textrm{tet.}}$ = $10^{-3}$ for double precision numbers, and we
suggest using this value for other materials. Since the material-dependent
variables, $F_{i}$ (velocity and spin-current matrix elements) and $x_{ij}$
(determined by band structure), are not included in the condition above, the
same value would work well for other materials.

###  II.3 Higher-Order Correction

In the ordinary tetrahedron method, the matrix elements and the energy
eigenvalues are linearly interpolated using the values at the four vertices.
The integrand can thus deviate from the true value due to the neglected
curvature in $k$ space of the numerator and the denominator. To correct this
error, Blöchl suggested an improved tetrahedron method [18], which has been
widely implemented in many computer programs. The method can take the
curvature effect into account using only the matrix elements and energy
eigenvalues at the four vertices without higher-order interpolation using the
Gauss theorem. Quantities such as the total energy or the charge (spin)
density can be calculated using the Blöchl’s correction. However, this method
is limited to smoothly varying functions which are interpolated by
polynomials, not by rational functions. Near the zero of the denominator in
$k$ space where $\Delta_{nm\textbf{k}}^{2}-(\hbar\omega+i\eta)^{2}=0$, the
integrand is rapidly varying. The gradient and the curvature inside a
tetrahedron are not well represented, resulting in the failure of the method.

Kawamura et al. introduced a novel, improved tetrahedron method applicable to
response functions [19]. The work is about finding a good linear interpolation
function from the third-order correction. The method uses additional 16 $k$
points surrounding the original tetrahedron together with the four original
vertices, and finds the 20 coefficients to construct third-order interpolating
polynomials for the matrix element and for the energy eigenvalue. Finally, one
fits a linear function using the least square method. As a result, the
formulas for the linear fitting functions are automatically calculated from
the matrix elements and the energy eigenvalues at the 20 $k$ points. This
scheme can also be employed with our solution to the numerical problem (round-
off errors) at the same time.

##  III COMPUTATIONAL DETAILS

The intrinsic contribution to both the static and dynamical SHC of platinum
was investigated by the simple summation with the adaptive smearing scheme [7]
and the tetrahedron method proposed here. The entire procedure from the self-
consistent calculations to the Wannier-interpolation is the same as that in
Ref. [9]. The norm-conserving pseudopotential and the pbesol exchange-
correlation energy functional [20] were used. The kinetic energy cutoff for
the wave functions was set to 60 Ry. The energy eigenvalues, the corresponding
Fermi-Dirac occupancy factor, and the matrix elements for the SHC were
computed via self-consistent and non-self-consistent calculations on $12\times
12\times 12$ and $8\times 8\times 8$ $k$-point grids, respectively, using the
plane-wave pseudopotential code PWSCF from the Quantum ESPRESSO package [21,
22], via the wannier90 package, and via the pw2wannier90 which is the
connecting program between PWSCF and wannier90  [23]. The initial guess at
MLWFs were $s$, $p$, and $d$ atom-centered nodeless orbitals, and the inner
energy window is 15 eV’s wide from the bottom of the valence-band minimum up
to 4 eV above the Fermi level.

Regarding the adaptive smearing for the simple summation, we followed the
scheme of Ref. [7]: We set $\eta$ in Eq. (11) as
$\eta=a|\frac{\partial\epsilon_{{n\mathbf{k}}}}{\partial\textbf{k}}-\frac{\partial\epsilon_{{m\mathbf{k}}}}{\partial\textbf{k}}|\Delta
k$ where $a$ is a dimensionless constant of order one (we set $a=\sqrt{2}$),
and $\Delta k$ is the distance between the nearest-neighbor $k$ points in the
interpolation grid. To avoid the aforementioned numerical instability for the
tetrahedron method, we set $\epsilon_{\textrm{tet.}}$ to be $10^{-3}$.

##  IV RESULTS

Figure 1: The combined electronic band structure of fcc platinum. The black
dots indicate the band structure obtained from ab initio calculations, and the
red curves were obtained from the Wannier interpolation. The Fermi level of
undoped platinum is located on the top of the valence band, as displayed by a
solid horizontal line. The black dashed horizontal line is the ceiling of the
inner window. Figure 2: The dynamical SHC of platinum calculated by the simple
summation with the adaptive smearing method and by the tetrahedron method with
higher-order correction. The $k$-point grid used for the numerical integration
is specified.

The band structure of platinum is shown in Fig. 1. 18 MLWFs successfully
reproduced the band structure within the inner window ($<4$ eV).

Figure 2 shows the comparison of the two methods for the dynamical SHC. We
used unshifted Monkhorst-Pack grids. The adaptive smearing scheme requires a
finer $k$ mesh to reach convergence. The tetrahedron method achieved
convergence at a relatively coarse $40\crossproduct 40\crossproduct 40$ grid.
On the other hand, the peak value at 0.57 eV of the real part obtained on
$100\crossproduct 100\crossproduct 100$ grids, using the adaptive smearing,
does not reach half of that obtained on $40\crossproduct 40\crossproduct 40$
grids using the tetrahedron method.

The imaginary part is zero below 0.48 eV, which is the threshold for a direct
band-to-band transition. The tetrahedron method is useful to see this feature,
since it reproduces the delta function better than the smearing method. This
result agrees with previous studies on the original tetrahedron method for the
integrals over the Fermi surface, such as in the case of density of states.

Figure 3: (a) $d_{\rm dyn.}$ [Eq. (25)] as a function of
$\epsilon_{\textrm{tet.}}$. (b) A comparison of computational time among the
three tetrahedron methods and the adaptive smearing method. (c) and (d)
Similar quantities as (a) and (b), respectively, for $d_{\rm stat.}$ [Eq.
(26)]. The vertical dashed line in (a) and (c) denotes
$\epsilon_{\textrm{tet.}}=10^{-3}$ selected for other calculations [the
results shown in Fig. 2 and Figs. 3(b) and 3(d)].

To quantify convergence, we defined the standard deviation of the dynamical
SHC as follows,

$d_{\rm
dyn.}=\sqrt{\frac{\displaystyle\int_{\omega_{1}}^{\omega_{2}}d\omega|\sigma^{z}_{xy}(\omega)-\sigma^{z}_{xy,{\rm
conv.}}(\omega)|^{2}}{\textstyle\omega_{2}-\omega_{1}}}\,,$ (25)

where $\omega$ is the frequency, $\omega_{1}=0$ eV, and $\omega_{2}=7$ eV. For
the static SHC, the Fermi-energy shift $E_{\rm F}$ is substituted for
$\omega$,

$d_{\rm stat.}=\sqrt{\frac{\displaystyle\int_{E_{{\rm F}1}}^{E_{{\rm
F}2}}dE_{\rm F}|\sigma^{z}_{xy}(0)-\sigma^{z}_{xy,{\rm
conv.}}(0)|^{2}}{\textstyle E_{{\rm F}2}-E_{{\rm F}1}}}\,,$ (26)

where $E_{{\rm F}1}=-7.1$ eV and $E_{{\rm F}2}=+3$ eV. The integrals were
evaluated by a summation using $\Delta\omega=0.1$ eV and $\Delta E_{\rm
F}=0.1$ eV. Using the tetrahedron method, the SHC calculated on $140\times
140\times 140$ grids for the dynamical and static SHC was considered the
converged one ($\sigma^{z}_{xy,{\rm conv.}}$).

To see the cutoff dependence of our tetrahedron method, we obtained $d_{\rm
dyn.}$ and $d_{\rm stat.}$ with five different $k$ meshes as a function of
$\epsilon_{\textrm{tet.}}$ values (Fig. 3). In general, the smaller
$\epsilon_{\textrm{tet.}}$ is, the larger the round-off error is. On the other
hand, the larger $\epsilon_{\textrm{tet.}}$ is, the larger the deviation of
the energy eigenvalues at the four vertices of a tetrahedron [Eqs. (22) and
(23)] is. Therefore, as Figs. 3(a) and 3(c) show, there is an optimal range
for $\epsilon_{\textrm{tet.}}$ to converge SHC, with which the standard
deviation of the error ($d_{\rm stat.}$ or $d_{\rm dyn.}$) is small enough
compared to the total SHC value $\sigma^{z}_{xy}\sim
O(10^{3})(\hbar/e){\rm(S/cm)}$. The dashed vertical lines indicate
$\epsilon_{\textrm{tet.}}=10^{-3}$ located in the middle of this optimal
range. This value was adopted in obtaining the results shown in Fig. 2 and
Figs. 3(b) and 3(d).

We compared the computational time to see which method is the most effective.
Figures 3(b) and 3(d) show, respectively, that $d_{\rm dyn.}$ and $d_{\rm
stat.}$ decrease as the $k$-point grid for integration becomes dense.
Convergence is not reached at all if tetrahedron method without the cutoff
treatment is used, whereas the tetrahedron method with the cutoff treatment
achieves convergence faster than the adaptive smearing method. The higher-
order correction proposed in Ref. [19] improves the convergence further.

Figure 4: The difference (absolute value) between the value of static SHC from
a single tetrahedron calculated using the cutoff treatment and the ‘exact’
value versus $x=x_{12}=x_{23}=x_{34}$ [Eq. (21)]. The ‘exact’ value was
obtained using quadruple precision numbers producing an ignorable round-off
error within the shown $x$ range. The black curve indicates the absolute value
of the difference between the exact value at a given $x$ and the value when
$x$ is strictly zero.

Figure 4 illustrates how much the SHC values for a single tetrahedron with the
cutoff treatment deviate from the ‘exact’ value as a function of
$x\,\,(=x_{12}=x_{23}=x_{34})$ [Eq. (21)]. With the treatment, the integration
values are fixed to be that of $|x|=\epsilon_{\textrm{tet.}}$; hence, the
nearly horizontal lines emerge when $|x|<\epsilon_{\textrm{tet.}}$. Too small
or too large $\epsilon_{\textrm{tet.}}$ leads to sizable errors, and
$\epsilon_{\textrm{tet.}}=10^{-3}$ displays the smallest error. For example,
$\epsilon_{\textrm{tet.}}=10^{-5}$ is too small, since the round-off error is
not yet corrected by the treatment until this error gets very large. On the
other hand, $\epsilon_{\textrm{tet.}}=10^{0}$ is too large because the energy
eigenvalues at the four vertices are changed too much by the treatment. These
facts are consistent with the results shown in Figs. 3(a) and 3(c). Without
the treatment, the error increases with decreasing $x$ due to the round-off
term $O(\epsilon)$. The slope of the diagonally scattered data in Fig. 4 is
$-3$, from which we can infer that the round-off error is proportional to
$1/x^{3}$. If all $|x_{ij}|$’s decrease at the same time, the error is not
proportional to $1/x^{2}$ as Eqs. (18) and (19). The reason for the slope
being -3 is explained in Appendix C.

Let $f(x)$ the formula evaluating the contribution to the integration from a
single tetrahedron with a given $x$. Then the slope of the black curve at
small $x$ approaches $+1$ because $f(x)-f(0)=f^{\prime}(0)x+O(x^{2})$, where
$f(0)$ is the limiting formula. Also, $f^{\prime}(0)$ can be found from the
explicit formulas of the tetrahedron method, such as Eqs. (B4) and (B6). Thus,
the slope of the logarithmic plot is +1 when $x$ is small enough.

If we adopt a cutoff, $\epsilon_{\textrm{tet.}}$, the error follows the black
curve when $x<\epsilon_{\textrm{tet.}}$. and the red scattered points when
$x>\epsilon_{\textrm{tet.}}$. The optimal cutoff is around $2\crossproduct
10^{-4}$ from the figure above since the error is minimized both when
$x<\epsilon_{\textrm{tet.}}$ and $x>\epsilon_{\textrm{tet.}}$ at the same
time. Even if we adopt this cutoff, however, the error is comparable to or
even larger than the error from our proposed simpler method with
$\epsilon_{\textrm{tet.}}=0.001$ (the green curve).

The Fermi-Dirac distribution we used is the step function, which cannot
describe the smooth Fermi-Dirac distribution for finite temperatures. By
scanning multiple Fermi levels and conducting a convolution of the integral of
the non-dissipative part and the Fermi-Dirac distribution, we can take finite
temperatures into account. This idea has been implemented in the WannierBerri
code [24].

Our method can be applied to evaluate the integration of much more complicated
functions with a different pole structure. From other correlation functions to
Green’s function or various self-energy, (for example, see Ref. [25]) there
are many quantities which can be expressed in the form of Eqs. (3)–(5) after
partial fraction decomposition.

The proposed tetrahedron method for the Kubo formula can also be
straightforwardly extended to deal with the ${n\mathbf{k}}$-dependent, finite
imaginary part of the quasiparticle self-energy. The arguments of logarithms
in Eqs. (14) and (15) are now complex numbers instead of absolute values of
real numbers. For this extension, when calculating ${\rm log}A-{\rm log}B$,
one should use the same branch cut for both log terms. For example, if
$A=e^{ia}$ and $B=e^{ib}$ with $a,b\in(-\pi,\pi]$, ${\rm log}A-{\rm
log}B=i(a-b)$, while ${\rm log}(A/B)=i(a-b+2n\pi)$ where $n$ is 0, +1, or -1
depending on $a$ and $b$. Therefore, one should use ${\rm log}A-{\rm log}B$
and not ${\rm log}(A/B)$ in deriving the results such as Eqs. (14) and (15).

To summarize, we developed a new tetrahedron method to accurately evaluate the
dissipative part of spectral functions. By discovering and solving a numerical
problem due to rounding of floating-point arithmetic, we achieved improved
convergence than the adaptive smearing method. A higher-order correction was
also combined with our method. We demonstrated the advantage of our method by
calculating the dynamical and static spin Hall conductivity of platinum. Our
method can be applied to a wide range of other physical quantities. We will
make the implementation of our proposed method publicly available through
wannier90  [23] soon.

##  V ACKNOWLEDGEMENTS

###### Acknowledgements.

We thank Ji Hoon Ryoo for providing a code and directions to calculate SHC
using MLWFs [9]. We also thank Sophie Beck, Marco Gibertini, Jerome Jackson,
Jason Kaye, Jae-Mo Lihm, and Jonathan R. Yates for useful discussions. This
work was supported by the Institute for Basic Science (No. IBSR009-D1) and by
the Creative-Pioneering Research Program through Seoul National University.
Computational resources were provided by KISTI Supercomputing Center (Grant
No. KSC-2020-INO-0078).

Figure 5: A schematic showing how six tetrahedra are generated from a given
parallelepiped (cube in the figure): 1248 (red lines), 1628, 1568, 1758, 1378,
and 1438. Here 18 is adopted to be the main diagonal (the red, dash-dotted
line), and six edges (12, 13, 15, 84, 86, and 87) and five diagonals
traversing a surface (dashed lines; 14, 17, 28, 38, 48) become edges of
tetrahedra.

##  Appendix A STANDARD WAY TO CONSTRUCT TETRAHEDRA

The original way to construct tetrahedra is as follows [18]. In the first
Brillouin zone, $N$ $k$ points are chosen uniformly, and the $N$
parallelepipeds whose vertices are the neighboring $k$ points are built. Each
parallelepiped is cut into six tetrahedra as depicted in Fig. 5. As a result,
the number of the simplexes becomes $6N$. The choice of the main diagonal
($18$ in the figure) is not unique, but it is advisable to select the shortest
one.

When the Fermi surface or the surface $\hbar\omega-\Delta_{nm\textbf{k}}=0$
passes through a tetrahedron, a surface of constant energy
$\epsilon_{{m\mathbf{k}}}=\epsilon_{\textrm{F}}$ or constant energy difference
$\Delta_{nm\textbf{k}}=\hbar\omega$ cuts it into smaller polyhedra. Figures
6(a)–6(c) depict how the surface cuts a tetrahedron depending on the value of
$\epsilon_{\textrm{F}}$ or $\hbar\omega$. Since $\epsilon_{{n\mathbf{k}}}$ and
$\epsilon_{{m\mathbf{k}}}$ are linearly interpolated inside a single
tetrahedron, each surface cuts a tetrahedron by a plane. The pieces cut by the
surface are polyhedrons that can be constructed from smaller tetrahedra [Fig.
6(d)].

##  Appendix B ANALYTIC TETRAHEDRON METHOD

Figure 6: Illustrations of how the Fermi surface or the surface $\hbar\omega-\Delta_{nm\textbf{k}}=0$ cuts a tetrahedron in the Brillouin zone. It is assumed that the indices of the vertices are in order of increasing $\epsilon_{{m\mathbf{k}}}$ or $\Delta_{nm\textbf{k}}$, and they are denoted as $\epsilon_{i}$ or $\Delta_{i}$ for $i=1,2,3,4$. (a) $\epsilon_{1}<\epsilon_{\textrm{F}}<\epsilon_{2}$ or $\Delta_{1}<\hbar\omega<\Delta_{2}$. (b) $\epsilon_{2}<\epsilon_{\textrm{F}}<\epsilon_{3}$ or $\Delta_{2}<\hbar\omega<\Delta_{3}$; the tetragonal surface is cut into two triangular pieces, A and B. (c) $\epsilon_{3}<\epsilon_{\textrm{F}}<\epsilon_{4}$ or $\Delta_{3}<\hbar\omega<\Delta_{4}$. (d) A schematic showing that the polyhedron in panel (b) is constructed from three tetrahedra — 1285 plus 5267 plus 5268. Range | Parametrization, $\partial(s_{1},s_{2},s_{3})/\partial(u,v)$
---|---
$\Delta_{1}<\hbar\omega<\Delta_{2}$ |  $\begin{pmatrix}-y_{1}&-y_{1}\\\ y_{2}&0\\\ 0&y_{3}\end{pmatrix}$, $y_{1}=\frac{\hbar\omega-\Delta_{1}}{\Delta_{2}-\Delta_{1}},y_{2}=\frac{\hbar\omega-\Delta_{1}}{\Delta_{3}-\Delta_{1}},y_{3}=\frac{\hbar\omega-\Delta_{1}}{\Delta_{4}-\Delta_{1}}$
$\Delta_{2}<\hbar\omega<\Delta_{3}$ |  $\begin{pmatrix}y_{1}&0\\\ 0&y_{2}\\\ 1-y_{1}-y_{3}&-y_{3}\end{pmatrix}$ for A, $\begin{pmatrix}y_{1}-y_{4}&-y_{4}\\\ y_{4}-1&y_{4}-1+y_{2}\\\ 1-y_{1}&0\end{pmatrix}$ for B, $y_{1}=\frac{\hbar\omega-\Delta_{4}}{\Delta_{2}-\Delta_{4}},y_{2}=\frac{\hbar\omega-\Delta_{1}}{\Delta_{3}-\Delta_{1}},y_{3}=\frac{\hbar\omega-\Delta_{1}}{\Delta_{4}-\Delta_{1}},y_{4}=\frac{\hbar\omega-\Delta_{3}}{\Delta_{2}-\Delta_{3}}$
$\Delta_{3}<\hbar\omega<\Delta_{4}$ |  $\begin{pmatrix}y_{1}&0\\\ 0&y_{2}\\\ 1-y_{1}-y_{3}&1-y_{2}-y_{3}\end{pmatrix}$, $y_{1}=\frac{\hbar\omega-\Delta_{4}}{\Delta_{2}-\Delta_{4}},y_{2}=\frac{\hbar\omega-\Delta_{4}}{\Delta_{3}-\Delta_{4}},y_{3}=\frac{\hbar\omega-\Delta_{4}}{\Delta_{4}-\Delta_{1}}$
Table 1: Parametrizations of the surface integrals [Eq. (4)], assuming
$\Delta_{1}\leq\Delta_{2}\leq\Delta_{3}\leq\Delta_{4}$. For
$\Delta_{2}<\hbar\omega<\Delta_{3}$, both triangles A and B shown in Fig. 6(b)
are parametrized.

Equations (3) – (5) can be integrated over each tetrahedron using linear
interpolations. Equations (3) and (5) are volume integrals, while Eq. (4) is a
surface integral. The integrands of the volume integrals are rational
functions in $k$ space, so the results are associated with logarithmic
functions. The volume integral is a triple integral of a linear function
divided by another linear function over a polyhedron, whose shape depends on
$\epsilon_{F}$ [Fig. 6]. The surface integral is a double integral on a plane
cut by a surface with a constant energy difference
$\Delta_{nm\textbf{k}}=\hbar\omega$ [Fig. 6].

All the tetrahedra and triangles can be parametrized in new coordinate
systems, $(s_{1},s_{2},s_{3})$ and $(u,v)$, respectively, defined by

$\displaystyle k_{i}=$ $\displaystyle\sum_{j}t_{ij}s_{j}$ (B1)
$\displaystyle\textbf{t}_{i}=$ $\displaystyle(t_{1i},t_{2i},t_{3i})\,,$ (B2)

and

$\begin{pmatrix}s_{1}\\\ s_{2}\\\
s_{3}\end{pmatrix}=\frac{\partial(s_{1},s_{2},s_{3})}{\partial(u,v)}\begin{pmatrix}u\\\
v\end{pmatrix},$ (B3)

where $\textbf{t}_{i}$ is the vector from vertex 1 to vertex $(i+1)$ [Fig. 6].
The Jacobian matrix in Eq. (B3) is given in Tab. 1. Under this change of
variables, one vertex is set to be the origin, while the other two vertices
are defined as $(u,v)=$(1,0), (0,1) for the triangular surface integral, and
the other three vertices are defined as $(s_{1},s_{2},s_{3})=$(1,0,0),
(0,1,0), (0,0,1) for the tetrahedral volume integral.

The expressions for Eqs. (3) - (5) are:

$\displaystyle\begin{split}I_{1,nm}(\hbar\omega)={}&-\frac{(1+x_{1})(1+x_{2})(1+x_{3})}{6x_{1}^{2}x_{2}^{2}x_{3}^{2}(x_{1}-x_{2})^{2}(x_{2}-x_{3})^{2}(x_{3}-x_{1})^{2}}\\\
&\times\frac{\textrm{det}(\textbf{t})}{\Delta_{4}+\hbar\omega}\sum_{i=1}^{4}F_{i}\left(\sum_{j=1}^{3}C_{ij}^{(1)}\xi_{j}+B_{i}^{(1)}\right)\end{split}$
(B4)
$\displaystyle\begin{split}I_{2,nm}(\hbar\omega)={}&\frac{|\frac{\partial{\textbf{k}}}{\partial{u}}\crossproduct\frac{\partial{\textbf{k}}}{\partial{v}}|}{|\nabla_{\textbf{k}}\Delta_{nm\textbf{k}}|}\crossproduct\\\
&\int_{0}^{1}\int_{0}^{1-u}du\,dv\,F_{nm{\bf k}(u,v)}\end{split}$ (B5)
$\displaystyle\begin{split}I_{3,nm}={}&\frac{(1+x_{1})(1+x_{2})(1+x_{3})}{2x_{1}^{2}x_{2}^{2}x_{3}^{2}(x_{1}-x_{2})^{2}(x_{2}-x_{3})^{2}(x_{3}-x_{1})^{2}}\\\
&\times\frac{\textrm{det}(\textbf{t})}{\Delta_{4}^{2}}\sum_{i=1}^{4}F_{i}\left(\sum_{j=1}^{3}C_{ij}^{(3)}\xi_{j}+B_{i}^{(3)}\right)\,.\end{split}$
(B6)

The expressions for $C_{ij}$, $B_{i}$,
$|\nabla_{\textbf{k}}\Delta_{nm\textbf{k}}|$, and $F_{nm{\bf k}(u,v)}$ are
given below. The definitions of the other used variables are presented in the
main text.

For $a,b,c\in\\{1,2,3\\}$, define $b\equiv a+1$ (mod 3), and $c\equiv a+2$
(mod 3). In other words, they are cyclic.

$\displaystyle\begin{split}C_{aa}^{(1)}={}&-(x_{b}-x_{c})^{2}x_{b}^{2}x_{c}^{2}\\\
&\times(3x_{a}^{2}-2(x_{b}+x_{c})x_{a}+x_{b}x_{c})\\\
C_{ba}^{(1)}={}&-(x_{b}-x_{c})^{2}x_{b}^{2}x_{c}^{2}x_{a}(1+x_{b})(x_{c}-x_{a})\\\
C_{ca}^{(1)}={}&+(x_{b}-x_{c})^{2}x_{b}^{2}x_{c}^{2}x_{a}(1+x_{c})(x_{a}-x_{b})\\\
C_{4a}^{(1)}={}&-(x_{b}-x_{c})^{2}x_{b}^{2}x_{c}^{2}(x_{a}-x_{b})(x_{c}-x_{a})\\\
B_{a}^{(1)}={}&x_{a}C_{4a}^{(1)}\\\
B_{4}^{(1)}={}&-x_{1}x_{2}x_{3}(x_{1}-x_{2})^{2}(x_{2}-x_{3})^{2}(x_{3}-x_{1})^{2}\end{split}$
(B7)
$\displaystyle\begin{split}C_{aa}^{(3)}={}&-(x_{b}-x_{c})^{2}x_{b}^{2}x_{c}^{2}(1+x_{a})\\\
&\times(2x_{a}^{3}+(3-x_{b}-x_{c})x_{a}^{2}-2(x_{b}+x_{c})x_{a}+x_{b}x_{c})\\\
C_{ba}^{(3)}={}&(1+x_{a})C_{ba}^{(1)}\\\
C_{ca}^{(3)}={}&(1+x_{a})C_{ca}^{(1)}\\\
C_{4a}^{(3)}={}&(1+x_{a})C_{4a}^{(1)}\\\ B_{a}^{(3)}={}&x_{a}C_{4a}^{(3)}\\\
B_{4}^{(3)}={}&B_{4}^{(1)}\end{split}$ (B8)

The parametrizations for the surface integral $I_{2,nm}$ are written in terms
of the values of energy differences. The detailed formulas are in Tab. 1. The
gradient of $\Delta_{nm\textbf{k}}$ is constant on a single surface:

$|\nabla_{\textbf{k}}\Delta_{nm\textbf{k}}|=\sqrt{|\sum_{i,j,k}t^{-1}_{ik}t^{-1}_{jk}\Delta_{i}\Delta_{j}|}\,.$
(B9)

Using the coefficients of $F_{nm\textbf{k}}$ in the newly parametrized
coordinate, $\tilde{F_{0}},\tilde{F_{1}},\tilde{F_{2}}$,

$\displaystyle\begin{split}F_{nm{\bf
k}(u,v)}=&\,F_{1}+(F_{2}-F_{1})\,s_{1}(u,v)\\\
&+(F_{3}-F_{1})\,s_{2}(u,v)+(F_{4}-F_{1})\,s_{3}(u,v)\\\
\equiv&\,\tilde{F_{0}}+\tilde{F_{1}}\,u+\tilde{F_{2}}\,v\,,\end{split}$ (B10)

we find that the surface integral is given by

$\int_{0}^{1}\int_{0}^{1-u}du\,dv\,F_{nm{\bf
k}(u,v)}=\frac{\tilde{F_{0}}}{2}+\frac{\tilde{F_{1}}+\tilde{F_{2}}}{6}\,.$
(B11)

Also,
$|\frac{\partial{\textbf{k}}}{\partial{u}}\crossproduct\frac{\partial{\textbf{k}}}{\partial{v}}|$
in Eq. (B5) is obtained using Eqs. (B1) – (B3) and the parametrization in Tab.
1 between the two coordinate systems $(s_{1},s_{2},s_{3})$ and $(u,v)$. The
integrals are hence calculable by collecting all the terms in Eqs. (B4) –
(B6).

##  Appendix C ANALYSIS OF THE ROUND-OFF ERROR SHOWN IN FIG. 4

According to Eqs. (18) and (19), the round-off error is proportional to
$O(\epsilon)/x^{2}$ if only one of $|x_{j}|$ or $|x_{j}-x_{k}|$ approaches 0.
In Fig. 4, however, the slope in logarithmic scale is $-3$, implying the error
is proportional to $O(\epsilon)/x^{3}$.

In order to understand this behavior, suppose that $x_{1}$, $x_{2}$, $x_{3}$,
$x_{1}-x_{2}$, $x_{2}-x_{3}$, and $x_{3}-x_{1}$ are all proportional to
$O(x)$. We are using log$(x)\doteq O(x)+O(\epsilon)$ or log1P$(x)\doteq
O(x)+O(x\epsilon)$, so $s=0$ or $1$ for $\xi_{l}=\textrm{log}(1+x_{l})\doteq
O(x_{l})+O(x_{l}^{s}\epsilon)$, respectively.

$\displaystyle\begin{split}&\frac{\sum_{l}C_{il}\xi_{l}+B_{i}}{x_{1}^{2}x_{2}^{2}x_{3}^{2}(x_{1}-x_{2})^{2}(x_{2}-x_{3})^{2}(x_{3}-x_{1})^{2}}\qquad\quad\\\
&\doteq\frac{O(1)}{x^{12}}[C_{ii}(O(x_{i})+O(x_{i}^{s}\epsilon))+C_{ij}(O(x_{j})+O(x_{j}^{s}\epsilon))\\\
&\qquad\qquad+C_{ik}(O(x_{k})+O(x_{k}^{s}\epsilon))+B_{i}]\\\
&=\frac{1}{x^{12}}[O(x^{8})(O(x)+O(x^{s}\epsilon))+O(x^{8})(O(x)+O(x^{s}\epsilon))\\\
&\qquad\qquad+O(x^{8})(O(x)+O(x^{s}\epsilon))+O(x^{9})]\\\
&=\frac{1}{x^{12}}[O(x^{12})+O(x^{8})O(x^{s}\epsilon)]\\\
&=\frac{O(x^{4-s})+O(\epsilon)}{x^{4-s}}\,.\end{split}$ (C1)

In the second to last equality in Eq. (C1), we have used the fact that the
coefficients of $O(x^{9})$, $O(x^{10})$, and $O(x^{11})$ should vanish in
order to have a non-divergent value in the limit $x\to 0$ if $\epsilon=0$.
Hence, the round-off error term is

$\frac{O(\epsilon)}{x^{4-s}}\,.$ (C2)

The error is proportional to $1/x^{3}$ as shown in Fig. 4 since log1P was used
in our implementation, i. e. , $s=1$.

## References

  * MacDonald  _et al._ [1979] A. H. MacDonald, S. H. Vosko, and P. T. Coleridge, Extensions of the tetrahedron method for evaluating spectral properties of solids, J. Phys. C 12, 2991 (1979).
  * Gilat and Bharatiya [1975] G. Gilat and N. R. Bharatiya, Tetrahedron method of zone integration: Inclusion of matrix elements, Phys. Rev. B 12, 3479 (1975).
  * Brener  _et al._ [1981] N. E. Brener, J. L. Fry, and R. A. Johnson, Matrix elements in the analytic tetrahedron method, Int. J. Quantum Chem. 20, 459 (1981).
  * Lehmann and Taut [1972] G. Lehmann and M. Taut, On the numerical calculation of the density of states and related properties, Phys. Status Solidi B 54, 469 (1972).
  * Kaprzyk and Mijnarends [1986] S. Kaprzyk and P. E. Mijnarends, A simple linear analytic method for brillouin zone integration of spectral functions in the complex energy plane, J. Phys. C 19, 1283 (1986).
  * Kaprzyk [2012] S. Kaprzyk, An algebraic formulation and implementation of the tetrahedron linear method for the brillouin zone integration of spectral functions, Comp. Phys. Commun. 183, 347 (2012).
  * Yates  _et al._ [2007] J. R. Yates, X. Wang, D. Vanderbilt, and I. Souza, Spectral and Fermi surface properties from Wannier interpolation, Phys. Rev. B 75, 195121 (2007).
  * Souza  _et al._ [2001] I. Souza, N. Marzari, and D. Vanderbilt, Maximally localized Wannier functions for entangled energy bands, Phys. Rev. B 65, 035109 (2001).
  * Ryoo  _et al._ [2019] J. H. Ryoo, C.-H. Park, and I. Souza, Computation of intrinsic spin Hall conductivities from first principles using maximally localized Wannier functions, Phys. Rev. B 99, 235113 (2019).
  * Mainkar  _et al._ [1996] N. Mainkar, D. A. Browne, and J. Callaway, First-principles lcgo calculation of the magneto-optical properties of nickel and iron, Phys. Rev. B 53, 3692 (1996).
  * Nagaosa  _et al._ [2010] N. Nagaosa, J. Sinova, S. Onoda, A. H. MacDonald, and N. P. Ong, Anomalous hall effect, Rev. Mod. Phys. 82, 1539 (2010).
  * Giustino [2017] F. Giustino, Electron-phonon interactions from first principles, Rev. Mod. Phys. 89, 015003 (2017).
  * Fan [1951] H. Y. Fan, Temperature dependence of the energy gap in semiconductors, Phys. Rev. 82, 900 (1951).
  * Migdal [1958] A. B. Migdal, Interaction between electrons and lattice vibrations in a normal metal, Zh. Eksp. Teor. Fiz. [Sov. Phys. JETP] 34, 996 (1958).
  * Guo  _et al._ [2005] G. Y. Guo, Y. Yao, and Q. Niu, Ab initio calculation of the intrinsic spin Hall effect in semiconductors, Phys. Rev. Lett. 94, 226601 (2005).
  * Guo  _et al._ [2008] G. Y. Guo, S. Murakami, T.-W. Chen, and N. Nagaosa, Intrinsic spin hall effect in platinum: First-principles calculations, Phys. Rev. Lett. 100, 096401 (2008).
  * Beebe [2017] N. H. F. Beebe,  _The mathematical-function computation handbook_ (Springer, Salt Lake City, 2017).
  * Blöchl  _et al._ [1994] P. E. Blöchl, O. Jepsen, and O. K. Andersen, Improved tetrahedron method for Brillouin-zone integrations, Phys. Rev. B 49, 16223 (1994).
  * Kawamura  _et al._ [2014] M. Kawamura, Y. Gohda, and S. Tsuneyuki, Improved tetrahedron method for the Brillouin-zone integration applicable to response functions, Phys. Rev. B 89, 094515 (2014).
  * Perdew  _et al._ [2008] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Restoring the density-gradient expansion for exchange in solids and surfaces, Phys. Rev. Lett. 100, 136406 (2008).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, _et al._ , Quantum ESPRESSO : a modular and open-source software project for quantum simulations of materials, J. Phys. Cond. Matt. 21, 395502 (2009).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, _et al._ , Advanced capabilities for materials modelling with Quantum ESPRESSO , J. Phys. Cond. Matt. 29, 465901 (2017).
  * Pizzi  _et al._ [2020] G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, _et al._ , wannier90 as a community code: new features and applications, J. Phys. Cond. Matt. 32, 165902 (2020).
  * Tsirkin [2021] S. S. Tsirkin, High performance Wannier interpolation of Berry curvature and related quantities with WannierBERRI code, npj Comput. Mater. 7, 1 (2021).
  * Lihm and Park [2020] J.-M. Lihm and C.-H. Park, Phonon-induced renormalization of electron wave functions, Phys. Rev. B 101, 121102 (2020).


# Modeling the response of a standard accretion disc to stochastic viscous
fluctuations

Naveel Ahmad<EMAIL_ADDRESS>Ranjeev Misra  Naseer Iqbal  Bari Maqbool
Mubashir Hamid  Department Of Physics, University of Kashmir, Srinagar-190006,
India  Inter-University Center for Astronomy and Astrophysics, Post Bag 4,
Ganeshkhind, Pune-411007, India

###### Abstract

The observed variability of X-ray binaries over a wide range of time-scales
can be understood in the framework of a stochastic propagation model, where
viscous fluctuations at different radii induce accretion rate variability that
propagate inwards to the X-ray producing region. The scenario successfully
explains the power spectra, the linear rms-flux relation as well as the time-
lag between different energy photons. The predictions of this model have been
obtained using approximate analytical solutions or empirically motivated
models which take into account the effect of these propagating variability on
the radiative process of complex accretion flows. Here, we study the variation
of the accretion rate due to such viscous fluctuations using a hydro-dynamical
code for the standard geometrically thin, gas pressure dominated $\alpha$-disc
with a zero torque boundary condition. Our results confirm earlier findings
that the time-lag between a perturbation and the resultant inner accretion
rate variation depends on the frequency (or time-period) of the perturbation.
Here we have quantified that the time-lag $t_{lag}\propto f^{-0.54}$, for
time-periods less than the viscous time-scale of the perturbation radius and
is nearly constant otherwise. This, coupled with radiative process would
produce the observed frequency dependent time-lag between different energy
bands. We also confirm that if there are random Gaussian fluctuations of the
$\alpha$-parameter at different radii, the resultant inner accretion rate has
a power spectrum which is a power-law.

###### keywords:

accretion, accretion discs; instabilities; hydrodynamics; methods: numerical;
X-rays: binaries

††journal: New Astronomy

##  1 Introduction

Black hole X-ray binaries are powered by an accretion disc. The X-rays are
emitted from the inner regions of the disc close to the black hole where the
characteristic time-scales are expected to be of the order of less than
seconds. Thus, it was rather surprising that the X-ray emission from these
systems vary over a wide range of time-scales. It is now believed that this
wide range occurs because viscous fluctuations with long characteristic time-
scales occur in the outer parts of the disc and subsequently these variations
propagate to the inner parts (Lyubarskii, 1997). Such a model qualitatively
explains the power-law form of the power spectra of the X-ray emission.
Moreover, these long time-scale fluctuations arising from large radii on
propagating inwards gets superimposed on the smaller time-scale fluctuations
arising there. Hence, this explains the linear positive relation between the
$r.m.s$ variability and the X-ray flux, that are observed in many X-ray
binaries (Uttley & McHardy, 2001; Gaskell, 2004; Uttley, 2004; Uttley et. al.,
2005; Heil & Vaughan, 2010; Heil et al., 2011, 2012). Thus, the model provides
a natural connection between the power spectra and the rms-flux relation which
has been studied by several authors (e.g. Ingram & Done, 2011; Heil et al.,
2012). Moreover, the same process should be acting on accreting systems such
as Active Galactic Nuclei and hence in general can be used to obtain the size
scale of the system (Kelly et al., 2009, 2011). By modeling the observed
fluctuations, it may also be possible to obtain direct constraints on the
hydrodynamic structure of the accretion disc (Titarchuk et al., 2007).

Another critical aspect of the temporal behaviour of black hole X-ray binaries
is that typically, the systems show hard lags such that the high energy
photons are delayed relative to the soft ones. The time-lags increase
logarithmically with energy and are surprisingly frequency dependent,
$t_{lag}\propto f^{-0.7}$ (Miyamoto et al., 1988; Nowak et al., 1999; Wilms
et. al., 2004). The magnitude of these lags are too large for them to be due
to light crossing time-scales or due to Comptonization unless the corona is
unphysically larger than $10^{9}$ cms (Hua et. al., 1999; Kroon & Becker,
2016). Instead, the propagation model provides an interpretation where if
softer photons are produced at a larger radii compared to the high energy
ones, there will be a time-lag between them corresponding to the time take for
the fluctuation to propagate from the outer radius to the inner regions
(Misra, 2000; Kotov et al., 2001). The energy and frequency dependence of
these time-lags can also be explained within such a framework (Misra, 2000;
Kotov et al., 2001). For many years, it has been puzzling why high accretion
rate systems do not seem to be affected strongly by the radiation pressure
instability (Lightman, 1974; Lightman & Eardley, 1974). The instability should
have made bright black hole X-ray binaries in their soft state (when they are
supposed to have a standard accretion disc) to have X-ray flux variations of
several orders of magnitude over time-scales of minutes. Interestingly, the
absence of such variation could be due viscous fluctuations in the disc
(Janiuk & Misra, 2012).

Given the success and the potential of such a propagation model to explain
many of the timing properties of black hole X-ray binaries, there have been
attempts to integrate the basic idea of the model with one which can also
explain the time averaged spectrum and other timing features such as Quasi-
periodic Oscillations (QPOs). A popular model to explain the hard X-ray
emission from these sources during the hard spectral state, is that the
standard accretion disc is truncated and the inner region is a hot flow
(Shapiro et. al., 1976; Done et. al., 2007). For such an accretion disc
geometry, there have been a series of works (Ingram & Done, 2011, 2012; Ingram
& van der Klis, 2013; Rapisarda et al., 2016), which incorporate the basic
effects of stochastic propagation model. Such an integrated approach, can
explain the time-averaged spectra and the temporal properties such as rms-flux
and time-lags arising due to propagation. The model can also explain the $\sim
1$ Hz QPO observed in these sources as relativistic precession of the hot
flow. There is no doubt that in the future, there will be more such attempts
to provide a comprehensive picture using more complex accretion geometries
with the propagation model being an integral part of such an approach.

Since such generic models are complex it is natural that these initial
endeavours that attempt to quantitatively explain the temporal features,
include stochastic propagation in an approximate manner where it is assumed
that the accretion rate varies at different radii and that the time taken for
the propagation to move inwards is the local viscous time-scale. On the other
hand, there are approximate analytical techniques which have been used to
describe the diffusion process in the accretion flow and quantify the basic
prediction of the propagating fluctuations (Lynden-Bell & Pringle, 1974;
Lyubarskii, 1997; Misra, 2000; Kotov et al., 2001). For example, while Misra
(2000) describes the propagation in terms of cylindrical sound waves, Kotov et
al. (2001) have used the approximate Green’s function formalism of Lynden-Bell
& Pringle (1974) and Lyubarskii (1997). Solving the disc equations in a
linearized form Psaltis & Norman (2000) have shown that any transition radius
in the disc will act like a low-pass filter to such fluctuations and have
argued that resonant features may cause Quasi-periodic behaviour. It is
possible, however, to solve the time dependent hydro-dynamical equations with
a stochastic viscosity to obtain the full global non-linear behaviour of the
disc. Cowperthwaite & Reynolds (2014) simulated such a disc with global
stochastic fluctuations in the viscosity and confirmed the non-linear
behaviour of the $r.m.s$ being proportional to the flux. They found that the
accretion rate variations at frequencies lower than the local viscous
frequency is coherent for different radii and hence revealing the basic idea
of propagation fluctuation. They computed the expected time-lags and found
that their magnitudes are smaller than the viscous time-scale and that they
are frequency dependent. Further, Hogg & Reynolds (2016) performed global
magneto-hydrodynamic simulation of accretion disc and confirmed that the
Maxwell stresses in the low frequency regime, indeed produce fluctuations
which propagate inwards. These simulation results, especially that the time-
lags are a function of frequency need to be incorporated into future models
which would strive to quantitatively explain the temporal behaviour of
accreting systems.

In these previous works, the fluctuations have been introduced in the
viscosity as inspired by magneto-hydrodynamic simulations (Hogg & Reynolds,
2016). However, in the standard $\alpha$-disc model (Shakura & Sunyaev, 1973),
the viscosity is parametrised by $\alpha$ such that the viscous stress is
$\alpha P$, where $P$ is the pressure. Viscous variations can be then
represented by variation in the $\alpha$ parameter. Such simulations have been
undertaken for, e.g., to study the effect of these variations on unstable
radiation pressure dominated discs (Janiuk & Misra, 2012) and to understand
the effect of X-ray irradiation on the time dependent properties of the outer
regions of an accretion disc (Maqbool et al., 2015). Since the nature
especially the time-dependent behaviour of the turbulent viscosity is largely
unknown, it is prudent to test different mechanisms by which the viscosity may
fluctuate.

In this work we study the response of a gas pressure dominated standard
accretion disc to viscus fluctuations at different radii and time-periods,
concentrating on the temporal behaviour of the accretion rate at different
radii rather than on the emergent spectrum. One of our motivation is to
confirm whether the salient features, especially the frequency dependent time-
lag are similar when the fluctuations are introduced in $\alpha$ rather than
the viscous stress. More importantly, the aim is to introduce simple empirical
functions which capture the radial and frequency dependence of the variability
and time-lags which may then be used in models that quantitatively explain
temporal features. The results obtained may perhaps be fairly generic despite
the simplifying assumption of a standard gas pressure dominated disc. The
expectation here is that empirical scaling laws obtained, may form the basis
for more complex studies involving more complex accretion geometry and the
coupling of the flow to the local radiative process.

In the next section, we list the basic time-dependent equations that describe
the standard gas pressure dominated accretion disc. In §3,we study the
response of such a disc to sinusoidal perturbations at different radii and
with different time-periods. In §5, we introduce Gaussian stochastic
perturbations at each radii and study the temporal behaviour of the inner
accretion rate. In the last section, we summarise and discusses the important
results of the work.

##  2 Time-dependent disc structure equations

Since we are interested in the long term temporal evolution of an accretion
disc on viscous time-scales we assume as in the standard disc, vertical hydro-
dynamical equilibrium leads to

$P=\frac{GM\Sigma}{r^{3}}\frac{h}{2}$ (1)

where $h$ is the half thickness of the disc, $\Sigma=2\rho h$ is the surface
density and the pressure $P$ is assumed to be due to gas pressure only i.e.
$P={\Sigma}kT/{2\mu hm_{H}}$. Here $T$ is the mid-plane temperature, $m_{H}$
is the mass of hydrogen atom and $\mu\sim 0.5$ is the mean molecular weight.
The disc is also assumed to be in thermal equilibrium such that the local
viscous heat dissipated is radiated. The viscous dissipation rate per unit
area is taken in the standard form to be

$Q_{vis}^{+}=\frac{3{\dot{M}}}{8\pi
r^{2}}\frac{GM}{r}\left(1-{\left(\frac{R}{r}\right)}^{0.5}\right)$ (2)

where $M$ is the mass of the black hole, $R$ is the inner radius of the disc,
$\dot{M}$ is the accretion rate and $r$ is the radius of disc. The disc is
assumed to be radiating as a black body, hence the energy flux per unit area
is given by (e.g. Shapiro & Teukolsky, 1986)

$F(r)\simeq\frac{\sigma T^{4}}{4\tau}=Q_{vis}^{+}$ (3)

where $T$ is the mid-plane temperature of the disc and $\sigma$ is the Stefan-
Boltzmann constant. The optical depth, $\tau$ is considered both due to
electron-scattering, ${\bar{\kappa}_{es}}$ as well as due to free-free
interaction, ${\bar{\kappa}_{ff}}$ and is given as

$\tau=\Sigma({\bar{\kappa}_{ff}}+{\bar{\kappa}_{es}})$ (4)

Opacity due to free-free interaction, ${\bar{\kappa}_{ff}}$ is given by
Kramer’s opacity law written as

$\bar{\kappa}_{ff}=0.64\times{10^{23}}\rho T^{-7/2}cm^{2}g^{-1}$

and the opacity due to electron scattering, ${\bar{\kappa}_{es}}=0.41$ cm2
g-1.

In the time dependent case, the surface density, $\Sigma$ changes with time if
the accretion-rate, $\dot{M}$ varies with radius and in such a situation the
disc behaviour is governed by the continuity equation;

$\frac{\partial{\Sigma}}{\partial t}+{\frac{1}{2\pi
r}}{\frac{\partial{\dot{M}}}{\partial{r}}}=0$ (5)

and the angular-momentum conservation equation needs to be retained in its
differential form,

$\frac{\partial}{\partial t}\left({\Sigma 2\pi
r}\sqrt{GMr}\right)+\frac{\partial}{\partial
r}\left(\dot{M}\sqrt{GMr}\right)=\frac{\partial}{\partial r}\left(\alpha P4\pi
r^{2}h\right)$

which can be simplified to

$\dot{M}\frac{\partial}{\partial
r}\left(\sqrt{GMr}\right)=\frac{\partial}{\partial r}\left(\alpha P4\pi
r^{2}h\right)$ (6)

The partial differential Equation (5) is solved numerically to obtain the time
dependent behaviour of an accretion disc using the standard finite
differencing scheme with the courant condition imposed as described in Maqbool
et al. (2015). At an outer boundary radius, the disc values are fixed to the
steady state ones while for the inner one we maintain the viscous stress free
condition.

We find that after a transient period the disc settles to a steady state with
a constant accretion rate $\dot{M_{0}}$ and having the steady state surface
density profile $\Sigma_{0}(r)$. A pertinent time-scale for the system at a
given radius is the viscous time-scale which is given by

$T_{vis}=\frac{2\pi}{\dot{M_{0}}}\Sigma_{0}{r^{2}}$ (7)

After this transient period which is roughly equal to few times the viscous
time-scale of the outer radius, time varying perturbations are introduced to
the viscous parameter $\alpha$ as described in the following sections, and the
time dependent properties of the disc are studied. For all the results
presented in this work, we adopt $\alpha=0.1$ and the steady state value of
the accretion rate to be $\dot{M_{0}}=10^{17}$ g s-1.

##  3 Sinusoidal perturbation at different radii

We first study the response of the disc to a sinusoidal perturbation of the
viscous parameter $\alpha$

${\alpha(r_{p},t)}=\alpha+\triangle{\alpha}sin\left(\frac{2\pi
t}{T_{p}}\right)$ (8)

at a particular perturbation radius, $r_{p}$ and with a time-period $T_{p}$.
The $\alpha$ variation causes the accretion rate to vary at $r_{p}$ and the
perturbation moves inwards to smaller radii. This is illustrated in Figure 1
which shows snapshots in time for a perturbation at $r_{p}=100r_{g}$ with
$\triangle{\alpha}=0.01$. While the amplitude of the accretion rate variation
is large at $r_{p}$, it decreases as the perturbation travels inwards. This is
further illustrated in Figure 2, where the root mean square ($r.m.s$) of the
accretion rate variability is plotted as a function of radii. Here the time-
period for the perturbation is 0.2 times the viscous time-scale at the
perturbation radius and the Figure shows the behaviour for three different
perturbation radii. The variability decreases rapidly and then saturates at
smaller radii. Such a behaviour is expected as the perturbation gets damped
when its time-period is smaller than the local viscous time-scale and is
unabated when the time-period is larger.

Figure 1: Snapshots of accretion rate versus radii at different times. The
perturbation has been introduced at the radius of 100$r_{g}$ with a time-
period of 50 s. Figure 2: The root mean square of the accretion rate
variability as a function of radius. The $r.m.s$ have been normalised to that
of the accretion rate at the perturbation radius. The time-periods for the
curves marked A, B and C are 0.2 times the viscous time-scale of the
perturbation radius of 100, 150 and 200 $r_{g}$, respectively. The amplitude
of the perturbation in all these cases decreases initially, but then saturates
with radius.

Figure 3: The figures show the variation of ratio of the $r.m.s$ of accretion
rate at the inner most radius (6$r_{g}$) normalised to the $r.m.s$ of
accretion rate at the perturbation radius vs time-period, $T_{p}$. The top
panel is for the case, when the perturbation radius, $r_{p}$ is at 100$r_{g}$,
where the viscous time-scale, $T_{vis}$ is 48 seconds, which is marked by a
vertical line. The lower panel is for the case, when the perturbation radius
is at 200$r_{g}$, where the viscous time-scale is 134 seconds. Both curves
have been fitted using a function of the form
$\frac{1}{1+\left({AT_{vis}}/{T_{p}}\right)^{p}}$ (shown by dotted lines). The
values of $A$ and $p$ are 0.21 (0.22) and 1.45 (1.42), respectively for
100$r_{g}$ (200$r_{g}$).

The saturation of the $r.m.s$ at small radii, allows us to quantify the effect
a perturbation has, by computing the ratio of the saturated $r.m.s$ (computed
at the inner edge) with that at the perturbation radius
$R=(r.m.s)_{in}/(r.m.s)_{r_{p}}$. This is illustrated in Figure 3, where the
ratio is plotted as a function of the time-period, $T_{p}$, when the
perturbation radius are $r_{p}=100$ (Top Panel) and $200$ $r_{g}$ (Bottom
panel). The corresponding viscous time-scales for these radii are $\sim 48$
and $\sim 134$ seconds, respectively and they are marked as vertical lines in
the Figure. The amplitude ratio remains nearly constant at unity when the
oscillation time-period is larger than the viscous time-scales and decreases
for smaller values. A phenomenological fit to the dependence of the ratio to
the time-period can be

$R(T_{p})=\frac{(r.m.s)_{in}}{(r.m.s)_{r_{p}}}=\frac{1}{1+\left(\frac{AT_{vis}}{T_{p}}\right)^{p}}$
(9)

where $T_{vis}$ is the viscous time-scale at the perturbation radius, $r_{p}$.
Comparison with the numerical values show that $A=$ 0.21 and 0.22 and $p=$
1.45 and 1.42 for $r_{p}=100$ and $200$ $r_{g}$ respectively, and are
represented by dotted lines in Figure 3. Thus, this functional form captures
the dependence of the inner disc variation to an outer disc fluctuation.

##  4 Time-lag for sinusoidal perturbations

It is expected that a perturbation at a radius $r_{p}$ would propagate inwards
and reach the inner radii after a time-lag equal to the viscous time-scale at
$r_{p}$. There are a couple of reasons this expectation might not be met.
First, the viscous time-scale computed using equation (7) is after all an
approximation and hence the actual time-lag associated with the perturbation
may differ by a factor of few. More importantly, in a cylindrical geometry,
one can expect that the time-lag should have a dependence on the time-period
(or frequency) of the perturbation. Unlike planar waves, the speed of
cylindrical sound waves depends on its frequency, which has been exploited to
explain frequency dependent time-lags in Cygnus X-1 (Misra, 2000). The
diffusion equation that describe the propagation here should also have this
property and we test whether this is the case, by computing the time-lag
between the perturbation and its response in the inner disc region.

Figure 4 shows the computed time-lag, $t_{lag}$ divided by the viscous time-
scale, $T_{vis}$ of the perturbation radius versus the time-period, $T_{p}$ of
the oscillation which is also divided by the same viscous time-scale. The
curves marked A, B and C correspond to the perturbation radii 100, 150 and 200
$r_{g}$ respectively. In this representation, all three curves line up along
each other and show that in general, for time-periods longer than the viscous
time-scale, the time-lags are nearly independent of the time-period and are
roughly 0.4 times the viscous time-scale. For shorter time-periods the time-
lags decrease proportionally with the time-period. A functional form for the
time-lag for time-periods less than viscous time-scale is

$\frac{t_{lag}}{T_{vis}}\sim 0.38\left(\frac{T_{p}}{T_{vis}}\right)^{0.54}$
(10)

We note this result that the time-lag increases with the time-period is
consistent with those found by Cowperthwaite & Reynolds (2014).

Figure 4: The figure shows the variation of the time-lag, $t_{lag}$ with time-
period, $T_{p}$, both normalised to the viscous time-scale, $T_{vis}$,
corresponding to the radius where sinusoidal perturbation of different time-
periods are given. The curves marked A, B and C correspond to the perturbation
radii of 100, 150 and 200 $r_{g}$, respectively. All these curves are
represented by $\frac{t_{lag}}{T_{vis}}\sim
0.38\left(\frac{T_{p}}{T_{vis}}\right)^{0.54}$ shown in the figure by curve D.
This function was fitted only for $\frac{T_{p}}{T_{vis}}<1$.

##  5 Multi-Gaussian perturbation at different radii

In an accretion disc, it is expected that there will be stochastic
perturbations at all radii which propagate inwards, causing the inner
accretion rate to vary in a wide range of frequencies. To get insight into the
nature of such a system we conduct the following numerical toy experiment.

We consider a disc with an outer radius of 35 $r_{g}$ where the corresponding
viscous time-scale is 10 seconds. We divide the disc in 64 bins which are
logarithmically spaced in radii. At each radial bin, the viscous parameter
$\alpha$ is perturbed by $\Delta\alpha$ where $\Delta\alpha$ is derived from a
Gaussian distribution having a width of 0.0005. A narrow distribution has been
chosen because for a larger width, there are occasionally large $\alpha$
variations. For these large variations, the numerical code requires very small
time steps and hence leads to very time consuming computations. The
perturbation is given after a time equal to the local viscous time-scale.
Thus, $\alpha$ at each radial bin, is forced to vary in a step like fashion
where the change occurs on a time equal to the local viscous time-scale. This
step like variation in $\alpha$ causes a spike in the local accretion rate
whenever $\alpha$ is suddenly changed.

Near the outer radii (i.e. at two radial bins away from the last one), this
causes the magnitude of the change in the accretion rate to have a sharp rise
and fall as shown in the top panel of Figure 5. The accretion rate variations
at this radius, correspond to the local $\alpha$ variation at that bin as well
as the variation induced in the outer two bins. The time-scales for these
spikes are approximately ten seconds i.e. the local viscous time-scale at 35
$r_{g}$.

Figure 5: The figure shows the variation of accretion-rate, $\dot{M}$, when
random Gaussian perturbations are added to the viscosity parameter, $\alpha$
at different radial points of the disc. The curves on the top, middle and
lower panel show the variation of accretion rates, $\dot{M}$ with time near
the outer radius ($34r_{g}$), in between radius ($21r_{g}$) and at the inner
radius ($7r_{g}$), respectively.

For smaller radii, the accretion rate varies due to the local step-like
variation of $\alpha$ and due to the propagation of the fluctuations from the
outer regions. The crucial element here is that the two causes for the local
accretion rate variation are multiplicative rather than being additive. Thus a
significant enhancement of the variability occurs at smaller radii as compared
to the larger ones. The accretion rate variation in an intermediate
(21$r_{g}$) and at an inner radius (7$r_{g}$) are shown in the middle and
bottom panels of Figure 5. As expected the variability increases as the radii
decreases. This is quantified further in Figure 6 where the power spectra of
the accretion rate variability of the three radii considered in Figure 5 are
shown and marked as A, B, and C. While the power spectra of the outer radius
(Curve A at 34$r_{g}$) is weaker, there is enhanced power in the inner radii
(Curve B and C at 21 $r_{g}$ and 7$r_{g}$). The inner radii power spectra have
a power-law like behaviour beyond 0.1 Hz with an index of $\sim 1.5$. We have
verified that changing the number of annular bins does not cause any
qualitative change in the power spectra.

Figure 6: The curves A, B and C shows the power spectrum corresponding to the
radii, whose accretion-rate variations with time are shown in the upper,
middle and lower panels, respectively in Fig 5. The spectral index for curve A
is $\sim 0.6$, while for curves B and C, is of the order of $\sim 1.5$. The
vertical line is the inverse of the viscous time-scale corresponding to
$35r_{g}$.

We caution against over interpretation of the results of this numerical toy
experiment, since for computational purposes the region of disc is rather
small and the perturbations are given in a rather ad-hoc manner. Nevertheless,
some of the essential elements of the propagation fluctuation model seem to be
captured, i.e. the power spectra of the accretion rate in the inner regions
has a power-law form and since the effect is multiplicative, the variability
is stronger in the inner region than in the outer.

##  6 Summary and discussion

We solve the time-dependent equations for a standard gas pressure dominated
accretion disc with perturbations on the local viscous parameter $\alpha$ to
understand the response of the accretion rate at different locations of the
disc.

For a sinusoidal perturbation at an outer radius $r_{p}$, we find that the
ratio of the amplitude of the accretion rate variation at the inner disc to
that at the perturbation radius, can be empirically represented by

$R(T_{p})=\frac{(r.m.s)_{in}}{(r.m.s)_{r_{p}}}=\frac{1}{1+\left(\frac{0.2T_{vis}}{T_{p}}\right)^{1.4}}$

where $T_{p}$ is the time-period of the sinusoidal perturbation and $T_{vis}$
is the viscous time-scale at the perturbation radius. We find that the time-
lag between the inner accretion rate variability with respect to that of the
perturbation radius, depends on the time-period of the oscillation and can be
empirically represented as

$\frac{t_{lag}}{T_{vis}}\sim 0.38\left(\frac{T_{p}}{T_{vis}}\right)^{0.54}$

for $T_{P}<T_{vis}$ and is nearly a constant for larger values.

We further show by having all along the disc, step-like fluctuations of the
viscous parameter $\alpha$ on the local viscous time-scale, that the power
spectrum of resultant inner accretion rate variability has a power-law shape
in frequency with index $\sim 1.5$.

The results obtained are consistent with those found by Cowperthwaite &
Reynolds (2014) which show that the qualitative response of the disc is same
when the fluctuations are introduced in the $\alpha$ or directly to the
viscosity which implies that perhaps they do not depend on the details of the
viscosity prescription. The empirical functions obtained here, especially the
dependence of the time-lag on frequency, will have utility in making models
that predict the power spectra and time-lags of X-ray binaries within the
framework of propagating fluctuations such as those described in Ingram & van
der Klis (2013); Rapisarda et al. (2016). In particular, the dampening of the
fluctuation can be expressed using the empirical relation above and the time-
lag between the accretion rate variation can be taken to be time-period (or
frequency) dependent. These will no doubt have an impact on both the
qualitative and quantitative interpretations that such models make.

However, it should also be emphasized that the results here are based on
simple assumptions which may not be valid in real accretion discs in X-ray
binaries. In X-ray binaries, the observed variability is often seen and
quantified in the X-ray emission arising from a corona on top of a disc or
from a hot inner region where the standard disc is truncated. In both cases
one needs to consider a more generalized disc than the standard gas pressure
dominated one assumed here. A perhaps important omission is that the pressure
support in the radial direction has been neglected, which while true for low
luminosity discs, would be violated in high luminosity ones and if there is a
hot inner disc. Radial pressure support will give rise to radially moving
sound waves which may significantly alter the results obtained here.

In conclusion, while we verify by solving the time-dependent disc equations
the salient features of the propagation fluctuation model, we also provide an
empirical relation for the time-period dependent time-lag between the
accretion rate fluctuations. Clearly, there is a scope to do more advanced
computation using more detailed physics with predictive powers. Such analysis
will provide insight into the geometry and nature of these systems, when
compared with the existing data and future ones expected from the newly
launched AstroSat.

##  7 Acknowledgments

The authors are grateful to ISRO-RESPOND Program for providing the financial
assistance under grant No. ISRO/RES/2/370 for carrying out this work. Four of
us (NA, NI, BM and MH) are highly thankful to IUCAA, Pune for providing the
computing and library facilities in completing this manuscript.

##  8 References

## References

  * Cowperthwaite & Reynolds (2014) Cowperthwaite P. S. & Reynolds C. S., 2014, ApJ, 791, 126.
  * Done et. al. (2007) Done C., Gierlinski M. & Kubota A., 2007 A&ARv, 15, 1.
  * Gaskell (2004) Gaskell C. M., 2004, ApJ, 612, L21.
  * Heil & Vaughan (2010) Heil L. M. & Vaughan S., 2010, MNRAS, 405, L86.
  * Heil et al. (2011) Heil L. M., Vaughan S. & Uttley P., 2011, MNRAS, 411, L66.
  * Heil et al. (2012) Heil L. M., Vaughan S. & Uttley P., 2012, MNRAS, 422, 2620.
  * Hogg & Reynolds (2016) Hogg J. D. & Reynolds C. S., 2016, ApJ, 826, 40.
  * Hua et. al. (1999) Hua X. M., Kazanas D. & Cui W., 1999, ApJ, 512, 793.
  * Ingram & Done (2011) Ingram A. & Done C., 2011, MNRAS, 415, 2323.
  * Ingram & Done (2012) Ingram A. & Done C., 2012, MNRAS, 419, 2369.
  * Ingram & van der Klis (2013) Ingram A. & van der Klis M., 2013, MNRAS, 434, 1476.
  * Janiuk & Misra (2012) Janiuk A. & Misra R., 2012, A&A, 540, A114.
  * Kelly et al. (2009) Kelly B. C., Bechtold J. & Siemiginowska A., 2009, ApJ, 698, 895.
  * Kelly et al. (2011) Kelly B. C., Sobolewska M. & Siemiginowska A., 2011, ApJ, 730, 52.
  * Kotov et al. (2001) Kotov O., Churazov E. & Gilfanov M., 2001, MNRAS, 327, 799.
  * Kroon & Becker (2016) Kroon J. J. & Becker P. A., 2016, ApJ, arXiv, 160301728K.
  * Lightman (1974) Lightman A. P., 1974, ApJ, 194, 419.
  * Lightman & Eardley (1974) Lightman A. P. & Eardley D. M., 1974, ApJ, 187, L1.
  * Lynden-Bell & Pringle (1974) Lynden-Bell D. & Pringle J. E., 1974, MNRAS, 168, 603.
  * Lyubarskii (1997) Lyubarskii Yu. E., 1997, MNRAS, 292, 679.
  * Maqbool et al. (2015) Maqbool B., Misra R., Iqbal N. & Ahmad N., 2015, MNRAS, 448, 3242.
  * Misra (2000) Misra R., 2000, ApJ, 529, 95.
  * Miyamoto et al. (1988) Miyamoto S., Kitamoto S., Mitsuda K. & Dotani T., 1988, Nature, 336, 450-452.
  * Nowak et al. (1999) Nowak M. A., Vaughan B. A., Wilms J., Dove J. B. & Begelman M. C., 1999, ApJ, 510, 874.
  * Psaltis & Norman (2000) Psaltis D. & Norman C., 2000, astro-ph/0001391.
  * Rapisarda et al. (2016) Rapisarda S., Ingram A., Kalamkar M., & van der Klis M., 2016, arXiv:1607.08178.
  * Shakura & Sunyaev (1973) Shakura N. I. & Sunyaev R. A., 1973, A&A, 24, 337.
  * Shapiro et. al. (1976) Shapiro S. L., Lightman A. P. & Eardley D. M., 1976, ApJ, 204, 187.
  * Shapiro & Teukolsky (1986) Shapiro S. L. & Teukolsky S. A., 1986, Black Holes, White Dwarfs and Neutron Stars: The Physics of Compact Objects, by Stuart L. Shapiro, Saul A. Teukolsky, pp. 672. ISBN 0-471-87316-0. Wiley-VCH , June 1986, 672.
  * Titarchuk et al. (2007) Titarchuk L., Shaposhnikov N. & Arefiev V., 2007, ApJ, 660, 556.
  * Uttley (2004) Uttley P., 2004, MNRAS, 347, L61.
  * Uttley et. al. (2005) Uttley P., McHardy I. M. & Vaughan S., 2005, MNRAS, 359, 345.
  * Uttley & McHardy (2001) Uttley P. & McHardy I. M., 2001, MNRAS, 323, 26.
  * Wilms et. al. (2004) Wilms J., Pottschmidt K., Nowak M. A., Gleissner T., Pooley G. G., Remillard R., Staubert R., Heindl W. A., Uttley P. & Fender R. P., 2004, NuPhS, 132, 420.