text:
\section{Introduction}
\label{intro}

Accreting black holes (BH) in active galactic nuclei (AGN) and BH binaries as well as neutron star (NS) binaries show strong variability of the X-ray flux over a very broad (several decades) frequency range: the time scale extends down to milliseconds in binary systems and minutes in AGN \citep{2000A&A...363.1013R,2004MNRAS.348..783M}.
The broad band component of the power density spectrum (PDS) is generally well described by a twice-broken power-law. In the case of some BH binaries, the PDS displays a narrow feature peaking in the range $\sim\,0.1-10\,{\rm Hz}$s indicated as quasi-periodic oscillations (QPOs). Light curves both in BH and NS systems show a log-normal flux distribution and a linear relationship between the root-mean square (rms) variability and the luminosity of the source \citep{2001MNRAS.323L..26U,2002PASJ...54L..69N,2004MNRAS.347L..61U,2004ApJ...612L..21G,2005MNRAS.359..345U}.

The light curves in different energy bands are correlated with each other, but the hard X-ray variability usually lags the soft X-ray variability \citep{1979ApJ...233..350P,1981ApJ...246..494N}. The magnitude of the time-lag depends on frequency, but typically it is of order of 1 per cent of the variability time scale.
Some sources show negative time-lags, with the soft-band variability lagging the hard-band variability \citep{2007MNRAS.382..985M,2011MNRAS.416L..94E,2013MNRAS.431.2441D}. This type of behavior has been detected in both supermassive and stellar mass BH.
The cross-correlation function between any two energy bands peaks at zero lag \citep{2000ApJ...537L.107M,2001AIPC..599..310P},
and emission from different energy bands is coherent at low frequencies and incoherent at relatively high frequencies \citep{1999ApJ...510..874N}. It is interesting that the same features are observed in the case of accreting non-magnetized white dwarfs \citep{2013MNRAS.431.2535S} and even young stellar objects \citep{2016MNRAS.463.2265S}.

The observed variability is naturally explained by the propagating fluctuations model \citep{1997MNRAS.292..679L}, in which  fluctuations arise at different radial coordinates of the accretion disc and then propagate towards the central object. Different time scales are injected into the accretion flow at different distances from the central object  \citep{1997MNRAS.292..679L,2001MNRAS.321..759C,2015ebha.confE...8I}. More specifically, the initial perturbations of viscosity lead to perturbations of the local mass accretion rate and surface density. The perturbations propagate through the accretion disc and cause local variability of the mass accretion rate at all radii. The variability of the local mass accretion rate leads to photon flux variability observed from different parts of the accretion disc. The propagating fluctuations model naturally explains the lagging of hard X-ray variability \citep{2001MNRAS.327..799K,2013MNRAS.434.1476I}.
Simplified Green functions were used for propagating fluctuations modeling of aperiodic X-ray variability by \cite{1997MNRAS.292..679L} and \cite{2001MNRAS.327..799K}, where the fluctuations were assumed to be additive. Multiplicative fluctuations in the accretion disc were considered by \cite{2013MNRAS.434.1476I} and \cite{2011MNRAS.415.2323I} but in all cases only inward propagation was considered.

In the case of accretion onto magnetized NSs the accretion disc is disrupted at the magnetospheric boundary, which is much larger than the NS radius. From the magnetospheric radius, matter follows the magnetic field lines and reaches the NS surface in small regions located at the magnetic poles, where most of the kinetic energy is released in the form of X-rays. The mass accretion rate variability produced in the accretion disc results in mass accretion rate variability at the magnetospheric radius producing variability of the final X-ray flux \citep{2009A&A...507.1211R}. Thus, the basic physical processes of generating mass accretion rate fluctuations and their evolution in the accretion disc that can be diagnosed by the X-ray flux variability are similar for BH and NS X-ray binaries (XRB).

Accretion onto weakly magnetized NSs results in a boundary layer between the accretion disc and NS surface, where a significant fraction of total X-ray flux is produced \citep{2010AstL...36..848I}. The variability in this case is strongly affected by processes within the boundary layer, which is beyond the scope of this paper.

{The propagating mass accretion rate fluctuations are determined by the geometry of the accretion disc and the physical conditions in the disc (in particular local viscosity and geometrical thickness). On the other hand, the propagating fluctuations can influence the general features of accretion disc such us, under certain conditions, its stability (see e.g. \citealt{2012A&A...540A.114J}).}

In order to explain the fast variability of accreting BHs and NSs it is necessary to precisely describe the processes of generation and subsequent propagation of fluctuations in the accretion disc. In this paper we focus on the propagation of fluctuations by viscous diffusion. The propagation of mass accretion rate fluctuations can be described by the diffusion equation written under specific assumptions about the accretion disc structure. In some cases the solution of the diffusion equation can be simplified and obtained in the form of Green functions \citep{1974MNRAS.168..603L}. In general the diffusion equation has to be solved numerically  \citep{2015arXiv151102356M}.

The key principles of the construction of a Green function for the diffusion equation of disc accretion was formulated by \cite{1952Lust}. The specific Green functions obtained so far in the literature were obtained for the case of a Newtonian potential and various sets of initial assumptions encoded in boundary conditions. (i) The Green function for the case of an infinite accretion disc with the inner radius located at $R_{\rm in}=0$ was derived by \cite{1974MNRAS.168..603L}. \cite{1991MNRAS.248..754P} discussed disks with a central source of angular momentum. (ii) The Green functions for the case of nonzero inner radius but infinite outer radius $R_{\rm out}$ were found by \cite{2011MNRAS.410.1007T}. This particular Green function can be useful for the description of accretion discs around magnetized sources, where the inner disc radius is defined to be the size of magnetosphere. (iii) The evolution of accretion discs with fixed outer radius was discussed by \cite{1998MNRAS.293L..42K} and \cite{2015ApJ...804...87L}. The same problem was solved numerically by \cite{2009MNRAS.399.1633Z} and a particular case was discussed by \cite{2001ApJ...563..246W}.

In this paper we discuss how Green functions can be included in the model in order to describe propagation of fluctuations of mass accretion rate through the accretion disc, {which is assumed to be radiatively effective (i.e. most of the viscously dissipated energy is radiated locally)}.
{For the first time we consider the exact description of the diffusion process in application to propagating mass accretion rate fluctuations. This description requires propagation of the fluctuations not only inwards but also outwards, which has not been considered before.}
We describe how the resulting aperiodic X-ray flux variability in different energy bands is affected by diffusion processes.
{We reproduce the basic features of propagating mass accretion rate fluctuations discussed earlier in literature and get the basic relations as a limiting cases of our general approach. At the same time we show that the accurate description of the diffusion  process implies features that have not been expected before: we show that high frequency mass accretion rate variability at larger radii lags the variability at smaller radii under certain conditions. It leads to important conclusions on the nature of soft lags observed in supermassive and stellar mass BHs.}
{For the first time} we consider the possibility of interactions between fluctuations. Our numerical results use the Linden-Bell/Pringle Green function \citep{1974MNRAS.168..603L}, but the developed formalism allows the use of any appropriate Green function. Though the Green functions we use are constructed for flat space time, it is possible to get basic features arising from a detailed description of propagation fluctuation process. The results are mainly addressed to BH XRBs, but the technique we apply is applicable to accreting magnetized NS binaries.

\section{Viscous diffusion in an accretion disc}

\subsection{The equation of viscous diffusion}

We consider a thin Keplerian accretion disc with geometrical thickness $H$ much smaller than the radial scale $R$: $H\ll R$. The gravitational potential is given by
\be\label{eq:NewtonianPot}
\phi_{\rm N}=-\frac{GM}{R},
\ee
where $M$ is a mass of a central object. The Keplerian angular velocity is $\Omega_{\rm K}=\sqrt{GM/R^3}$. Differential rotation and viscous stress lead to an angular momentum exchange between adjacent rings. As a result, angular momentum is transported outwards.

The viscous force is dissipative and the work it does on adjacent rings goes into heat.
The viscous stress $t_{r\varphi}$ is defined as the viscous force per unit area. It is proportional to the shear and can be written as $t_{r\varphi}=\eta R \d\Omega_{\rm K}/\d R$, where the coefficient of proportionality $\eta$ is called dynamic viscosity. The kinematic viscosity is defines as $\nu=\eta/\rho$, where $\rho$ is local mass density.

The evolution of the disc surface density due to viscous processes is described by the diffusion equation, which can be obtained from the laws of conservation of energy and momentum (e.g. \citealt{2002apa..book.....F}):
\be\label{eq:DifEqGen}
\frac{\partial \Sigma(R,t)}{\partial t}=\frac{1}{R}\frac{\partial}{\partial R}\left[R^{1/2}\frac{\partial}{\partial R}\left(3\nu \Sigma R^{1/2}\right)\right],
\ee
where $\Sigma$ is the surface density and $t$ is time. The kinematic viscosity $\nu$ has the meaning of a diffusion coefficient.
It has to be mentioned that the evolution of the surface density is expected to be described by diffusion process on sufficiently long time scales only. On very short time scales the variability is defined by other processes including propagation of sound and magnetohydrodynamic (MHD) waves.

According to the $\alpha$-prescription \citep{1972AZh....49..921S,1973A&A....24..337S} the viscous stress is described by $t_{r\varphi}=\alpha P$, {where $P$ is the gas pressure in the disc and $\alpha<1$ is a numerical factor.}
\footnote{{The total pressure in the accretion disc is composed of gas pressure and radiation pressure. If the viscous stress is parameterised by total pressure dominated by radiation pressure, the disc is locally unstable \citep{1974ApJ...187L...1L}. However, if the stress is due to magnetic turbulence it scales with gas pressure alone \citep{1989ApJ...342...49S}.}}
It is equivalent to the prescription
\be\label{eq:AlphaPrescrNu}
\nu=\frac{2}{3}\alpha c_{\rm s}H=\frac{2}{3}\frac{\alpha c^2_{\rm s}}{\Omega_{\rm K}},
\ee
where $c_{\rm s}=(P/\rho)^{1/2}$ is the local isothermal sound speed.

The diffusion equation (\ref{eq:DifEqGen}) can be simplified by using specific angular momentum $h=(GMR)^{1/2}$ as a variable instead of $R$. Then the equation is reduced to \citep{1974MNRAS.168..603L}
\be\label{eq:DifEqN}
\frac{\partial \Sigma}{\partial t}=D_{\rm N}(h,M)\frac{\partial^2 \left(h\nu\Sigma\right)}{\partial h^2},\quad
D_{\rm N}(h,M)=\frac{3}{4}\frac{(GM)^2}{h^3}.
\ee

\subsection{Green function}

In general the diffusion equation is not linear because the kinematic viscosity $\nu$ can be a function of local physical conditions (for example, surface density or local temperature). But if this is not the case and the kinematic viscosity is a function of radius only, the diffusion equations (\ref{eq:DifEqGen}) and (\ref{eq:DifEqN}) are linear and can be solved by using the Green function method \citep{1952Lust}. Then the solution of the diffusion equation is given by
\be
\Sigma(R,t)=\int\limits_{R_{\rm in}}^{R_{\rm out}}G(R,R',t-t_0)\Sigma(R',t_0)\d R' ,
\ee
where $G(R,R',t)$ is a Green function, $\Sigma(R,t_0)$ is the surface density distribution over the radial coordinate $R$ at $t=t_0$, and $R_{\rm in}$ and $R_{\rm out}$ are inner and outer disc radius respectively. The Green function is defined by the kinematic viscosity $\nu(R)$, the inner and outer disc radii and the boundary conditions there.

The local mass accretion rate is defined by the surface density and radial velocity of accretion
$v_{\rm R}$: $\dot{M}(R,t)=2\pi R\Sigma(R,t)v_{\rm R}$. It can be also expressed as a function of surface density and viscosity:
\beq\label{eq:mdot00}
\dot{M}(R,t)=6\pi R^{1/2}\frac{\partial}{\partial R}\left(\nu\Sigma(R,t)R^{1/2}\right).
\eeq
If the kinematic viscosity is given by a power law $\nu\propto R^n$ the mass accretion rate $\dot{M}(R,t)\propto R^{1/2} \frac{\partial}{\partial R}\left(\Sigma(R,t)R^{n+1/2}\right)$ and can be represented by Green functions as well:
\beq\label{eq:mdot01}
\dot{M}(R,t)=\int\limits_{R_{\rm in}}^{R_{\rm out}}G_{\dot{M}}(R,R',t-t_0)\Sigma(R',t_0)\d R',
\eeq
where $G_{\dot{M}}(R,R',t)$ is a Green function for the mass accretion rate.
If the initial surface density fluctuations on top of the diffusion process are described by $\Sigma^*(R,t)$, then the resulting fluctuations of the mass accretion rate are given by a convolution \citep{1997MNRAS.292..679L}:
\beq\label{eq:mdot01_}
\dot{M}(R,t)=\int\limits_{R_{\rm in}}^{R_{\rm out}}G_{\dot{M}}(R,R',t)\otimes_t \frac{\partial \Sigma^*(R',t)}{\partial t}\d R'
\eeq
($\otimes_x$ is sing of convolution in $x$-variable).
According to (\ref{eq:mdot00}) the Green function for the mass accretion rate $G_{\dot{M}}(R,R',t)$ is
\be
\label{eq:mdot02}
G_{\dot{M}}(R,R',t)=6\pi R^{1/2}\frac{\partial}{\partial R}\left(\nu G(R,R',t)R^{1/2}\right).
\ee

If the kinematic viscosity does not depend on time and is given by
\be\label{eq:KinVic00}
\nu=\nu_0 (R/R_0)^n,
\ee
where $n\in (0;2)$ is a model parameter, the solution of the differential equation simplifies significantly. The corresponding Green functions for the surface density were found for a few particular cases (see \citealt{1974MNRAS.168..603L,2011MNRAS.410.1007T,2015ApJ...804...87L}).
The typical time scale of processes in the accretion disc which arise out of viscous diffusion is the viscous time scale
\be
t_{\rm v}(R)\approx \frac{2R^2}{3 \nu(R)}.
\ee
If the kinematic viscosity is described by equation (\ref{eq:KinVic00}) the viscous time scale $t_{\rm v}\propto R^{2-n}$.
It is useful to introduce the corresponding viscous frequency
\be\label{eq:viscous_frequency}
f_{\rm v}(R)\equiv \frac{1}{t_{\rm v}(R)}.
\ee
In case of viscosity prescription given by (\ref{eq:KinVic00}) the viscous frequency $f_{\rm v}\propto R^{n-2}$.

If the viscosity in the accretion disc is described by constant $\alpha$-parameter, it puts constraints on the power-law index $n$ in the description of kinematic viscosity (\ref{eq:KinVic00}). According to equation (\ref{eq:AlphaPrescrNu}) the kinematic viscosity can be represented as $\nu= \alpha\Omega_{\rm K}R^2(H/R)^2$. If the gas pressure and Kramer opacity dominate in the accretion disc, then $(H/R)\propto R^{1/8}$ \citep{2007ARep...51..549S}. As a result, the kinematic viscosity $\nu\propto R^{3/4}$ and therefore $n=3/4$ for the case of constant $\alpha$-parameter (although we note that there is no a priori reason to assume constant $\alpha$).
The viscous frequency is related to the Keplerian frequency $f_{\rm K}=\Omega_{\rm K}/(2\pi)$ as
\be\label{eq:frequencies}
f_{\rm v}(R)=3\pi\alpha f_{\rm K}(R)\left(\frac{H(R)}{R}\right)^2.
\ee

\subsection{Lynden-Bell/Pringle Green functions}
\label{sec:LB}

\begin{figure*}
\centering
\includegraphics[width=8cm]{pic1a.pdf}
\includegraphics[width=8cm]{pic1b.pdf}
\caption{Evolution of (a) surface density $\Sigma$ and (b) local mass accretion rate $\dot{M}$. The inner radius of the disc is taken to be $R_{\rm in}=0$, the surface density at $t=0$ is given by $\delta$-function at $R'=10$. Different curves are given for various time intervals measured in viscous time at radial coordinate $R=1$ ($\tau\equiv t/t_{\nu}(R=1)$): $\tau=0.002,\,0.01$ and $0.05$. It is clearly seen that fluctuations of the surface density and mass accretion rate propagate inwards and outwards. Parameter $n=1$.}
\label{pic:GfTanaka}
\end{figure*}

\begin{figure*}
\centering
\includegraphics[width=8cm]{pic2a.pdf}
\includegraphics[width=8cm]{pic2b.pdf}
\caption{The mass accretion rate caused by viscous diffusion of infinitely thin ring, which surface density at $t=0$ is given by $\delta$-function located at $R'=10$, as a function of dimensionless time (a) at radius $R=1$ and (b) at radius $R=50$. The inner disc radius is taken at $R_{\rm in}=0$. At radius $R=50$, which outside  radius $R'=10$, the mass accretion rate is negative initially, but it becomes positive after sufficiently large time. Parameter $n=1$.}
\label{pic:Mdot(t)}
\end{figure*}

For the case of $R_{\rm in}=0$ and $R_{\rm out}=\infty$ the Green function for the surface density was found by \cite{1974MNRAS.168..603L}. The Green function is given by
\beq
G(R,R',\tau)=\frac{1}{2(2-n)}\tau^{-1}R_c^{-1}\left(\frac{R'}{R_c}\right)^{1-n}\left(\frac{R'}{R}\right)^{n+1/4} \nonumber \\
\times\exp\left[-\frac{2l^2}{\tau}\left( \left(\frac{R'}{R_{\rm c}}\right)^{2-n}+ \left(\frac{R}{R_{\rm c}}\right)^{2-n}  \right)\right] \nonumber \\
\times I_{l}\left(\frac{4l^2}{\tau}\left(\frac{RR'}{R^2_{\rm c}}\right)^{\frac{2-n}{2}}\right),
\eeq
where $\tau=t/t_{\rm v}$ is time measured in units of viscous time at fiducial radius $R_{\rm c}$, $I_{\nu}$ is a modified Bessel function of the first kind and $l=(4-2n)^{-1}$.

According to equation (\ref{eq:mdot02}), the Green function for the mass accretion rate is
\beq
\label{eq:Gf_MdotLB}
G_{\dot{M}}(R,R',\tau)=
6\pi\nu_0 l^2 \tau^{-2} R^{\frac{3-2n}{4}} R ^{\prime \frac{9-2n}{4}} R_{\rm c}^{2n-2} \nonumber \\
\times \exp\left[-\frac{2l^2}{\tau} \left( \left(\frac{R'}{R_{\rm c}}\right)^{2-n}+ \left(\frac{R}{R_{\rm c}}\right)^{2-n}  \right) \right] \nonumber \\
\times \left\{
I_{l-1}\left(\psi \right)-I_{l}\left( \psi \right)\left(\frac{R}{R'}\right)^{1-n/2}
\right\},
\eeq
where
$$\psi=\frac{4l^2}{\tau}\left(\frac{RR'}{R^2_{\rm c}}\right)^{\frac{2-n}{2}}.$$
It is interesting to note a symmetry property of the mass accretion rate Green function:
\be\label{eq:dotMsym}
G_{\dot{M}}(R,R',\tau)=a^{1-n}G_{\dot{M}}(aR,aR',a^{2-n}\tau),
\ee
where $a$ is an arbitrary constant.

If the argument of the Bessel functions in (\ref{eq:Gf_MdotLB}) is small enough, which is the case at sufficiently large $\tau$:
$$\tau\gg 4l^2\left(\frac{\max(R,R')}{R_{\rm c}}\right)^{2-n}= 4l^2\frac{\max(t_{\rm v}(R),t_{\rm v}(R'))}{t_{\rm v}(R_{\rm c})},$$
one can use the approximation:
\beq
G_{\dot{M}}(R,R',\tau)\approx
6\pi\nu_0 l^2 R^{\frac{3-2n}{4}} R^{\prime\frac{9-2n}{4}} R_{\rm c}^{2n-2} \nonumber \\
\times\frac{\tau^{-(l+1)}}{\Gamma(l)}
\left[2l^2 \left(\frac{R R'}{R^2_{\rm c}}\right)^{1-n/2}\right]^{l-1},
\eeq
where $\Gamma(x)$ is the gamma function.

According to (\ref{eq:Gf_MdotLB}) the mass accretion rate caused by surface density concentrated at $R'$ at $t=0$ is positive in the inner disc region ($R<R'$). In the outer regions ($R>R'$)  the Green function (\ref{eq:Gf_MdotLB}) gives a negative mass accretion rate at first,
\footnote{Negative mass accretion rate {(i.e. mass transfer from smaller to larger radii)} is to be understood as the disturbance on top of a, typically larger, positive average mass accretion rate $\dot{M}_0$, i.e., one resulting in a smaller, but still positive, mass accretion rate. }
which takes away the angular momentum of the inner parts of the disc. However, after a sufficiently long time the mass accretion rate at a given radius $R>R'$ becomes positive. The change of sign of the mass accretion rate happens at approximately the viscous time scale corresponding to the radius $R$ (see Appendix \ref{App:Mdot}).

\subsection{Green function for truncated disc with given inner radius}

In a more general case with inner radius $R_{\rm in}$ \citep{2011MNRAS.410.1007T}, zero torque at $R_{\rm in}$ and $R_{\rm out}=\infty$ the Green functions are given by
\beq\label{eq:GfTanaka}
& G(R,R',t)=\left(1-\frac{n}{2}\right)R^{-n-1/4}R'^{5/4}R_{\rm in}^{n-2} & \\
&
\times\int\limits_{0}^{\infty}[J_l(kx)Y_l(k)-Y_l(kx)J_l(k)][J_l(kx')Y_l(k)-Y_l(kx')J_l(k)]& \nonumber \\
&
\times(J_l^2(k)+Y_l^2(k))^{-1}
\nonumber \\
&
\times\exp\left[-2\left(1-\frac{n}{2}\right)^2 k^2\frac{t}{t_{{\rm v,\,in}} }\right]k\d k,
& \nonumber
\eeq
where $x=(R/R_{\rm in})^{1-n/2}$, $t_{\rm v,\, in}=\frac{2}{3}R_{\rm in}^2/\nu(R_{\rm in})$ is local viscous time at the inner disc radius and $J_{l}(x)$ and $Y_l(x)$ are the Bessel functions of the first and second kind. The corresponding Green function for the mass accretion rate can be found according to equation (\ref{eq:mdot02}).

The Green functions for the case of $R_{\rm in}=0$ and $R_{\rm out}<\infty$ were found by \cite{2015ApJ...804...87L}.

The exact boundary conditions affect the typical time scales in the accretion disc (see Appendix\,\ref{App:Mdot}). The influence of advanced Green functions (with boundary conditions given at arbitrary radial coordinates and derived under the assumption of non-Newtonian gravitational potential) will be explored in further works.

\section{Exploring properties of the Green function}

In this section, we explore the properties of the Lynden-Bell \& Pringle Green function introduced in the previous Section. The properties of a more general Green function (for the case of given inner or outer radius) are qualitatively similar.

\subsection{Time domain}

The Green function of the surface density in the time domain $G(R,R',t)$ describes the evolution of the surface density given initially by a $\delta$-function at radial coordinate $R'$: $\Sigma(R,t=0)=\delta(R-R')$. The viscous force produces an angular momentum exchange between adjacent rings, such that the fast inner rings pass their angular momentum to slower outer rings. As a result, a part of the material accretes inwards losing angular momentum while another part moves outwards carrying out angular momentum. This process leads to spreading of an initially narrow ring (see Fig.\,\ref{pic:GfTanaka}a). The corresponding evolution of the mass accretion rate is given by the Green function for the mass accretion rate $G_{\dot{M}}(R,R',t)$ (see Fig.\,\ref{pic:GfTanaka}b). Spreading of the initial ring is manifested initially by a positive mass accretion rate at $R<R'$ (see Fig.\,\ref{pic:Mdot(t)}a) and a negative mass accretion rate at $R>R'$ (see Fig.\,\ref{pic:Mdot(t)}b). The radial coordinate where the mass accretion rate changes its sign moves outwards with time, which means that the angular momentum is picked up finally by a small fraction of the material, while most of it accretes onto the central object (see Fig.\,\ref{pic:Mdot(t)_}).

\begin{figure}
\centering
\includegraphics[width=8cm]{pic3.pdf}
\caption{The mass accretion rate caused by surface density perturbation at $R'=10$ at $t=0$ as a function of dimensionless time at radius $R=50$. The inner disc radius is taken at $R_{\rm in}=0$. One can see that almost all material accretes inside after sufficient amount of time. Parameter $n=1$.}
\label{pic:Mdot(t)_}
\end{figure}

We apply Green functions in order to describe the evolution of fluctiations arising on top of the average accretion flow. {Therefore, the negative mass accretion rates predicted by the Green functions results only in lower mass accretion rate at certain regions of the accretion disc, rather than genuinely negative accretion rates.
\footnote{Note that we consider the case of linear equation of viscous diffusion.}
The diffusion of a perturbation requires transfer of angular momentum from the inner parts of the disc to its outer parts. The extra angular momentum supplied from the inner parts of the disc slightly slows down the local accretion process and results in lower mass accretion rate.}

\subsection{Frequency domain}

The evolution of mass accretion rate in the time domain is described by convolution of the Green function and the initial fluctuations of the surface density (\ref{eq:mdot01_}).
The convolution turns to production in the frequency domain, which simplifies the problem. It makes the properties of Green functions in the frequency domain essentially important.

The Green function in the frequency domain is obtained by the Fourtier transform
\be
\overline{G}_{\dot{M}}(R,R',f)=\int\limits_{-\infty}^{\infty}\d x\,G_{\dot{M}}(R,R',x)\,e^{-2\pi i f x}.
\ee
{The obtained Green function shapes variability that originated from radius $R'$, at given radius $R$.}

The qualitative behavior of the Green function in the frequency domain is different for perturbations propagating initially inwards ($R<R'$) and outwards ($R>R'$).

In the case of inwards propagation, the absolute value of the Green function is $|\overline{G}_{\dot{M}}|\approx 1$ for frequencies much lower that the local viscous frequency $f_{\rm v}(R)$, and is suppressed for frequencies above $f_{\rm v}(R)$ (see Fig.\,\ref{pic:PS13} \textit{left}). It results in suppression of high frequency variability originating from distant outer radii. The suppression at high frequencies depends on the separation between radial coordinates of initial and final fluctuations ($R'$ and $R$): the larger the separation, the stronger Green function suppression. The suppression is different for different viscosity prescriptions given by the parameter $n$: the lower $n$ is, the larger the range of viscous time scales over the radial coordinates in the disc is, the stronger the suppression of the Green function (see Fig.\,\ref{pic:PS18} \textit{left}).

In the case of fluctuations propagating outwards ($R>R'$) the absolute value of the Green function $|\overline{G}_{\dot{M}}|$ is suppressed at high frequencies and at low frequencies, reaching its maximum at frequency close to local viscous frequency $f_{\rm v}(R')$ (see Fig.\,\ref{pic:PS13} \textit{right}). This property is caused by the fact that initially negative mass accretion rate at $R>R'$ turns to zero after a while and becomes positive on a sufficiently long timescale. The typical timescale of this process is close to the local viscous time scale (see Appendix\,\ref{App:Mdot}), which results in maximum of the $|\overline{G}_{\dot{M}}|$ being at $f_{\rm v}(R')$ (see Fig,\,\ref{pic:PS18} \textit{right}).

The absolute value of the Green function for outwards propagation at its maximum can be comparable to that for inward propagation (see Fig.\,\ref{pic:14_LB}). In this case the influence of the inner parts of the accretion disc on the local mass accretion rate variability can dominate over the influence of the outer parts of the disc. However, the influence of inner disc applies only to the closest radial coordinates and strongly is suppressed for distant outer regions (see Fig.\,\ref{pic:PS13} \textit{right}).

In Fig.\,\ref{pic:GF_phase}, we plot the phase angle in radians of the Green function for two values of $n$. For low frequencies ($f<<f_v$), we see that $\dot{M}(R<R',t)$ lags behind $\dot{M}(R=R',t)$, which simply represents the inwards propagation of the density perturbation at $R'$. Conversely, $\dot{M}(R>R',t)$ \textit{leads} $\dot{M}(R=R',t)$ (note that a phase angle $>\pi$ and $<2\pi$ is indistinguishable from a phase angle $<0$ and $>-\pi$), because the accretion rate upstream of the perturbation is negative (e.g. Fig.\,\ref{pic:Mdot(t)}b). For higher frequencies, we see that phase wrapping causes oscillations in the phase angle.

According to (\ref{eq:dotMsym}) the Lynden-Bell Green function (see Section\,\ref{sec:LB}) has the property:
\beq
\overline{G}_{\dot{M}}(R,R',f)=a^{-1}\overline{G}_{\dot{M}}\left(aR,aR',\frac{f}{a^{2-n}}\right),
\eeq
where $a>0$ is an arbitrary constant. Therefore the Green function in this particular case can be considered as a function of two variables instead of three (i.e. no need to consider a different $R'$ value to what we already plotted in Fig.\,\ref{pic:14_LB}).

\begin{figure*}
\centering
\includegraphics[width=16cm]{pic4.pdf}
\caption{The absolute value of Green function of the mass accretion rate in the frequency domain. The inner disc radius is located at $R_{\rm in}=0$, the initial perturbation of the surface density is at $R'=10$. The radial coordinates of final perturbations are taken to be $R=1,\,2,\,4,\,8$ (solid, dashed, dotted and dashed-dotted lines respectively), i.e. perturbations propagate inwards (left) and $R=11,\,20,\,40,\,80$ (solid, dashed, dotted and dashed-dotted lines respectively), i.e. perturbations propagate outwards (right). The parameter $n=1$.}
\label{pic:PS13}
\end{figure*}

\begin{figure*}
\centering
\includegraphics[width=16cm]{pic5.pdf}
\caption{The absolute value of the Green function of the mass accretion rate in the frequency domain. The inner disc radius is fixed at $R_{\rm in}=0$, the initial perturbations are at $R'=10$. The radial coordinate of final variability is at $R=1$ (left) and $R=40$ (right). Black solid, red dashed and blue dotted lines are given for $n=1.5,\,1$ and $0.5$ respectively. Suppression of variability at high frequencies is stronger for lower parameter $n$. Suppression of variability at low frequencies for outwards propagation is stronger for higher parameter $n$.}
\label{pic:PS18}
\end{figure*}

\begin{figure*}
\centering
\includegraphics[width=18cm]{pic6.pdf}
\caption{The absolute value of the Green  function for mass accretion rate in the frequency domain as a function of final radial coordinate $R$ and dimensionless frequency  $f/f_{\rm v}(R=1)$. The results are given for parameters $n=0.5$ (left) and $n=1.5$ (right). The initial perturbation is given by $\delta$-function and located at radial coordinate $R'=25$.}
\label{pic:14_LB}
\end{figure*}

\begin{figure*}
\centering
\includegraphics[width=18cm]{pic7.pdf}
\caption{The phase angle of the Green function for the mass accretion rate in the frequency domain as a function of radial coordinate $R$ and dimensionless frequency $f/f_{\rm v}(R=1)$. The results are given for for parameters $n=0.5$ (left) and $n=1.5$ (right). The initial perturbation is given by $\delta$-function and located at $R'=25$.}
\label{pic:GF_phase}
\end{figure*}

\section{The propagating fluctuations model using Green functions}

In the propagating fluctuation model, the initial variability of mass accretion rate is assumed to be stirred up through out the disc, i.e. fluctuations are introduced at a \textit{range} of $R'$. These fluctuations propagate to other parts of the accretion disc and, in general, interact (multiplicatively) with one another. In previous treatments of the model, only inwards propagation has been considered and/or a very simple Green function has been assumed to describe propagation \citep{1997MNRAS.292..679L,2001MNRAS.327..799K,2013MNRAS.434.1476I}.
Here, we consider both inwards and outwards propagation of fluctuations. In the simplest case, we may consider the fluctuations to simply pass one another. In this case, linearity is preserved. In general though, fluctuations interact with one another as they propagate: e.g. inward propagating fluctuations from a large $R'$ will modulate and be modulated by outward propagating fluctuations from a smaller $R'$. This general case is non-linear.

\subsection{Mass accretion rate fluctuations: the general case}

The total mass accretion rate consists of the average mass accretion rate $\dot{M}_0$ and fluctuations $\dot{m}(R,t)$ on top of the average value:
\be
\dot{M}(R,t)=\dot{M}_0+\dot{m}(R,t).
\ee
We assume that local perturbations of surface density $A(R,t)$ caused by MHD are $A(R,t)\propto \dot{M}(R,t)/R$.

Indeed, the evolution of surface density is given by
\footnote{{Using the following designation for the partial derivative: $\partial_x\equiv \partial/\partial x$.}}
$$R\partial_t\Sigma=\partial_R\dot{M}.$$
The local mass accretion rate can be represented as $\dot{M}=2\pi R v_R\Sigma$, where $v_R$ is local radial velocity. The radial velocity is defined by the diffusion process and can be disturbed by local interaction of the accretion flow with magnetic field. Then we can represent the radial velocity as
$$
v_r(R,t)=v_r^{\rm (d)}(R,t)[1+\vartheta(R,t)],
$$
where $v_r^{\rm (d)}(R,t)$ is the radial velocity defined by the diffusion process and $\vartheta(R,t)\ll 1$ represents fluctuations of radial velocity caused by local MHD processes. Then the local surface density evolution is defined by the relation {
\beq
R\partial_t \Sigma&=&\partial_R\left\{R\Sigma v_r^{\rm (d)}(R,t)[1+\vartheta(R,t)]\right\}\nonumber\\
&=&\partial_R\left\{\dot{M}^{\rm (d)}(R,t)[1+\vartheta(R,t)]\right\},\nonumber
\eeq
where $\dot{M}^{\rm (d)}(R,t)$ is the mass accretion rate caused by viscous diffusion only. Differentiation over the radial coordinate gives
\beq
R\partial_t \Sigma &=&\partial_{R}\dot{M}^{\rm (d)}(R,t)\left[1+\vartheta(R,t)\right]+\dot{M}^{\rm (d)}(R,t)\partial_R\vartheta(R,t) \nonumber \\
&\approx&\partial_{R}\dot{M}^{\rm (d)}(R,t)+\dot{M}^{\rm (d)}(R,t)\partial_R\vartheta(R,t).
\eeq
The first term in the right-hand side (rhs) of relation defines variations of the surface density due to the diffusion process, while the second term defines new perturbations due to local fluctuations of viscosity. Hence, the initial perturbations of the radial density $R\Sigma(R)$ are proportional to local mass accretion rate: $\partial_t\Sigma^*\propto \dot{M}/R$. Therefore, the perturbations of the surface density can be represented as follows:}
\be\label{eq:LocPert01}
A(R,t)=\frac{a(R,t)}{R}\frac{\dot{M}(R,t)}{\dot{M}_0}=\frac{a(R,t)}{R}\left(1+\frac{\dot{m}(R,t)}{\dot{M}_0}\right),
\ee
where $a(R,t)$ is defined by local physics of interaction between accretion flow and magnetic field and does not depend on conditions in other parts of accretion disc.
Then we get a non-linear equation for the perturbations of the mass accretion rate in the time domain:
\be
\label{eq:mdot_tdomain}
\dot{m}(R,t)=\int\limits_{R_{\rm in}}^{R_{\rm out}}\d R'\,G_{\dot{M}}(R,R',t)\otimes_t
 \frac{a(R',t)}{R'}\left(1+\frac{\dot{m}(R',t)}{\dot{M}_0}\right)
\ee

The local perturbations of the surface density in the frequency domain are given by
\be\label{eq:LocalPertA}
\overline{A}(R,f)=
\frac{\overline{a}(R,f)}{R}\otimes_f
\left(\delta(f)+\overline{\frac{\dot{m}(R,f)}{\dot{M}_0}}\right),
\ee
where $\otimes_f$ denotes a convolution in frequency.

The local perturbations of surface density propagate through the accretion disc according to the solution of the diffusion equation (\ref{eq:DifEqGen}) and define the final perturbations of the mass accretion rate. If the final mass accretion rates are given by convolution according to equation (\ref{eq:mdot01}) the perturbations of mass accretion rate in the frequency domain are given by
\beq\label{eq:NoiceMult0}
\overline{\dot{m}}(R,f)=
\int\limits_{R_{\rm in}}^{R_{\rm out}}\d R'\,\overline{G}_{\dot{M}}(R,R',f)\overline{A}(R',f),
\eeq
where $\overline{G}_{\dot{M}}(R,R',f)$ is the corresponding Green function in the frequency domain and $R_{\rm in}$ and $R_{\rm out}$ are inner and outer radii of accretion disc.
\footnote{The inner and outer radii in equation (\ref{eq:NoiceMult}) are defined by the accretion disc geometry and the physical conditions in the disc. The outer radius $R_{\rm out}$ cannot exceed the radius where the effective temperature of accretion disc is lower than $\simeq 6500\,{\rm K}$. Otherwise the accretion disc becomes thermally unstable because of hydrogen recombination and a dramatic change of opacity \citep{2001NewAR..45..449L}. Moreover, the initial perturbations of the mass accretion rate due to MHD processes are expected to be weaker in case of low temperature and weakly ionized plasma.}
Substitution of expression (\ref{eq:LocalPertA}) into (\ref{eq:NoiceMult0}) gives
\beq\label{eq:NoiceMult}
&\overline{\dot{m}}(R,f)&=\int\limits_{R_{\rm in}}^{R_{\rm out}}\frac{\d R'}{R'}\,\overline{G}_{\dot{M}}(R,R',f)\overline{a}(R',f)+ \\
&&\frac{1}{\dot{M}_0}\int\limits_{R_{\rm in}}^{R_{\rm out}}\frac{\d R'}{R'}\,
\overline{G}_{\dot{M}}(R,R',f)
\left[\overline{a}(R',f)\otimes_f\overline{\dot{m}}(R',f)\right].  \nonumber
\eeq
If the perturbations of the mass accretion rate are significantly smaller than the average mass accretion rate ($\dot{m}(R,t)\ll\dot{M}_0$) then the perturbations of mass accretion rate are defined by the first term on the rhs of equation (\ref{eq:NoiceMult}). In this case they are represented by a superposition of fluctuations originating at different radii and modified by the diffusion process.
If perturbations of the mass accretion rate are comparable with the average mass accretion rate then the final result is affected by both terms of the rhs of equation (\ref{eq:NoiceMult}).

\subsection{The equation in discrete form: particular case of Green function given by $\delta$-function}

In order to directly compare our formalism with that of \cite{2013MNRAS.434.1476I}
let us consider how the formalism transforms if the flow is split into $N$ rings of radii $R_N<R_{N-1}<...<R_2<R_1$. Then the equation (\ref{eq:mdot_tdomain}) takes form
\be
\dot{m}(R_i,t)\simeq\sum \limits_{j>i} \Delta R_j\,G(R_i,R_j,t)\otimes_t
 \frac{a(R_j,t)}{R_j}\left(1+\frac{\dot{m}(R_j,t)}{\dot{M}_0}\right).\nonumber
\ee
If the Green function is given by $\delta$-function (the case of arbitrary Green function is discussed in Appendix \ref{Sec:EqDiscForm2}): $G(R_i,R_j,t)=\delta(t-\Delta t_{ij})$, where $\Delta t_{ij}$ is the propagation time from $R_j$ to $R_i$, we get
\be
\dot{m}(R_i,t)\simeq\sum \limits_{j>i} \Delta R_j\,\frac{a(R_j,t-\Delta t_{ij})}{R_j}\left(1+\frac{\dot{m}(R_j,t-\Delta t_{ij})}{\dot{M}_0}\right).\nonumber
\ee
Then the total mass accretion rate at the ring $i$ can be represented as follows:
\beq
& \dot{M}(R_i,t)=\dot{M}(R_{i-1},t-\Delta t_{i,(i-1)})&   \nonumber  \\
&+ \Delta R_{i-1}\,\frac{a(R_{i-1},t-\Delta t_{i,(i-1)})}{R_{i-1}}\left(1+\frac{\dot{m}(R_{i-1},t-\Delta t_{i,(i-1)})}{\dot{M}_0}\right)&
\nonumber
\eeq
or
\beq
&\dot{M}(R_i,t)=\dot{M}(R_{i-1},t-\Delta t_{i,(i-1)})& \nonumber\\
 &\times\mathcal{A}_{n-1}(t-\Delta t_{i,(i-1)}),&
\eeq
where
\beq
\mathcal{A}_{n-1}(t-\Delta t_{i,(i-1)})=
1+ \frac{\Delta R_{i-1}}{R_{i-1}}\,\frac{a(R_{i-1},t-\Delta t_{i,(i-1)})}{\dot{M}_0}. \nonumber
\eeq
Therefore considering this particular case, the total mass accretion rate in the $i$-th ring can be obtained from the mass accretion rate in the $(i+1)$th ring by multiplying on the stochastic function with a mean of unity.
Then the total mass accretion rate in the frequency domain can be calculated:
\beq
\overline{\dot{M}}(R_n,f)=\overline{\dot{M}}(R_1)\coprod\limits_{i=1}^{n-1}\overline{\mathcal{A}}_{i}(f)e^{2\pi i f\Delta t_{i,(i-1)}},
\eeq
where $\overline{\mathcal{A}}_{i}(f)$ is the Fourier transform of the function $\mathcal{A}_{i}(t)$. If fluctuations at different radii are not correlated with each other the power of the mass accretion rate is
\beq
|\overline{\dot{M}}(R_n,f)|^2=\overline{\dot{M}}^2(R_1)\coprod\limits_{i=1}^{n-1}|\overline{\mathcal{A}}_{i}(f)|^2,
\eeq
This result agrees with the formalism, which was developed earlier by \cite{2013MNRAS.434.1476I}.

\section{Power spectra of initial perturbations}

We assume that the initial fluctuations of the surface density are stochastic and have a random phase but well-defined power spectrum. This power spectrum depends on the physical process generating variability.

The initial perturbations might be caused by MHD turbulence, which serves as the source of viscosity in the accretion flow, in a form of fluctuations in magnetic stress \citep{1991ApJ...376..214B,1995ApJ...440..742H,1995ApJ...446..741B}.

The typical time-scale can be defined by the dynamo process \citep{2004MNRAS.348..111K}, which is generally shorter than the viscous time scale. However, the time scale for producing magnetic field sufficiently coherent to affect the accretion rate can be longer and comparable with viscous time-scale \citep{2003ApJ...593..184L}.

An accurate description of the emergence of new fluctuations of the mass accretion rate requires the introduction of their cross-spectrum (see Appendix \ref{App:PSD}), while the PDS of the initial fluctuations at each radial coordinate can be used for approximate description of real process.

Our formalism allows us to use any power spectra for the new fluctuations. Here we will focus only on some possibilities.

\subsection{White noise}

The power spectrum of white noise is given by a constant. Though this kind of process is not physical, it can be used as an approximation in cases when the break frequency of the initial perturbations is well above local viscous frequency.

\subsection{Lorentzian profile}

The power spectrum at each radius can be represented by a zero-centred Lorentzian breaking at frequency $f_{\rm br}(R)$:
\be\label{eq:LorentzProf}
S_a(f)=\left|a(f)\right|^2 = \frac{2 P(R)}{\pi}\frac{f_{\rm br}(R)}{(f_{\rm br}(R))^2+f^2},
\ee
where $P=\int_0^\infty\d f S_a(f)$ is the total power of the local process.
The frequency $f_{\rm br}$ is a model parameter and has to be provided by the physical model of the initial fluctuations in the accretion disc. The breaking frequency is likely to be within the range between viscous $f_{\rm v}$ and Keplerian frequency $f_{\rm K}$ \citep{2015arXiv151205350H}: $f_{\rm v}<f_{\rm br}<f_{\rm K}$, which are separated by about 2 orders of magnitude (see Eq.\,(\ref{eq:frequencies})).
The integrated power of initial perturbations can depend on the radial coordinate in the accretion disc. \textit{In this paper, however, we take the integrated power to be constant over the accretion disc.}

The Lorentzian profile arises from the idea that the initial fluctuations in the time domain decay exponentially. Indeed, the Fourier transform of the function
$g(t)=\sin(2\pi f_0 t)e^{-t/T}$, where $f_0$ and $T$ are constants, is given by
$$
\overline{g}(f)\propto \frac{1}{((f-f_0)^2+(1/T)^2)^{1/2}}.
$$
The corresponding power is
$$
|\overline{g}(f)|^2\propto \frac{1}{(f-f_0)^2+(1/T)^2}.
$$
In that sense the Lorentzian profile given by equation (\ref{eq:LorentzProf}) implies that the initial fluctuations decay on timescale $\sim 1/\Delta f_{\rm br}$. In this paper we use zero-centured Lorentzian profiles, but if the accretion disc at some radial coordinate is disturbed by some process of characteristic frequency $f^*$, it would be reasonable to use the Lorentzian profile centured at frequency $f^*$.

\section{Fourier properties of the local mass accretion rate variability}

\subsection{Power spectrum}

The perturbations of the mass accretion rate originate from initial stochastic variability and are, therefore, also stochastic. They cannot be well described by any particular realisation of the mass accretion rate process in either time or frequency.  However, they can be described in terms of the average power spectrum and cross-spectrum. In the most general case when the mass accretion rate perturbations are given by (\ref{eq:mdot_tdomain}), PDS of the mass accretion rate at $R$ is defined by equation (see Appendix \ref{App:PSD}):
\beq
\label{eq:S_mdot}
S_{\dot{m}}(R,f)=\int\limits_{R_{\rm in}}^{R_{\rm out}}\frac{\d R'}{(R')^2}\Delta R(R')
|\overline{G}_{\dot{M}}(R,R',f)|^2 \nonumber \\
\times
S_{a}(R',f)\otimes_f \left[\delta(f) + \frac{S_{\dot{m}}(R',f)}{\dot{M}^2_0}\right],
\eeq
where $S_{a}(R,f)$ is the power spectrum of initial perturbations at radius $R$, $\Delta R(R)$ is the radial scale on which the initial perturbations can be considered as coherent with one another. This scale can be different at different radii in the accretion disc. It is likely that this coherence scale is $\Delta R\sim H$, where $H$ is the geometrical thickness of the accretion disc at given radius \citep{2015arXiv151205350H}. However, in the case of accretion onto highly magnetized objects the scale can be larger.

In the linear approach, when only the first term in the square brackets in the rhs of equation (\ref{eq:S_mdot}) is taken into account, the final power spectrum of the mass accretion rate is defined by the initial variability at every radius given by $S_{a}(R',f)$, and the transfer function power $|\overline{G}(R,R',f)|^2$. Because the transfer function suppresses variability at frequencies higher than local viscous frequency $f_{\rm v}(R')$, every region in the accretion disc effectively produces fluctuations at frequencies below the corresponding $f_{\rm v}(R')$.  As a result, there should be a break in the power spectrum, which corresponds to the lowest viscous frequency in that part of the accretion disc, where initial fluctuations are produced. The lowest frequency corresponds to the outer radius $R_{\rm out}$ and the break arises even in the case when the initial perturbations are given by white noise (see Fig.\,\ref{pic:17_2}). If the initial fluctuations are produced below some local frequency $f_{\rm br}(R)$ only, then there should be another break in power spectrum, which locates at $\sim f_{\rm br}(R)$ (see red dashed line in Fig.\,\ref{pic:17_2}). At different radii the location of the second break is also different in general. Above this frequency the variability is not produced effectively by local processes and additionally suppressed by the transfer function. Thus, one would expect double broken power spectrum with breaks corresponding to the inner and outer disc radii.

The inclination of the power spectrum between two slopes and above the second slope depends on the power spectra of initial perturbations. We describe the power spectrum of initial perturbations by zero-centred Lorentzians breaking at frequency $f_{\rm br}(R)$. If, for example, the power spectrum drops exponentially at frequencies above $f_{\rm br}(R)$ (like in the model of \citealt{1997MNRAS.292..679L}), the final power spectrum is steeper. In this case it is, indeed, $\propto 1/f$ between two slopes as found by \citealt{1997MNRAS.292..679L}.

The second term in equation (\ref{eq:S_mdot}) describes the response of the newly arising fluctuations on local variability of the mass accretion rate (the amplitude of new fluctuations is proportional to the local mass accretion rate, which is variable). This term includes all non-linear interactions between inward, outward and local fluctuations. The second term is negligible if the mass accretion rate variability amplitude is much smaller than the average mass accretion rate.

Because the initial fluctuations of the surface density scale with the mass accretion rate, the PDS of initial perturbations is $S_{a}(R,f)\propto \dot{M}^2$.
The rms of local fluctuations of the mass accretion rate is
$$({\rm rms})^2 \propto \int\limits_{f>0} \d f\, S_{\dot{m}}(f)$$
Therefore, if the nonlinear effects are not taken into account and there is only the first term in the integrand in equation (\ref{eq:S_mdot}), the rms of local mass accretion rate is proportional to the mass accretion rate:
\be
{\rm rms}\propto \dot{M}.
\ee
This provides a neat explanation for the linear rms-flux relation ubiquitously observed from accreting objects.

Interestingly, the nonlinear term in equation (\ref{eq:S_mdot}) can cause a deviation from this linear rms-flux relation.

\begin{figure}
\centering
\includegraphics[width=8.5cm]{pic8a.pdf}
\includegraphics[width=8.5cm]{pic8b.pdf}
\includegraphics[width=8.5cm]{pic8c.pdf}
\caption{The power spectra of the mass accretion rate at dimensionless radius $R=50$ for different parameters $n$. Black solid (red dashed) line corresponds to the initial fluctuations given by white noise (zero-centred Lorentzian breaking at local viscous frequency). The accretion disc inner radius is at $R_{\rm in}=0$. The initial fluctuations are produced at radii $R'\in(0,200)$. $\Delta R=1$.}
\label{pic:17_2}
\end{figure}

\subsection{Cross-spectrum}

The cross-spectrum of mass accretion rate variability for radii $R_1$ and $R_2$ is (see Appendix \ref{App:Cross})
\beq
\label{eq:C_mdot}
C_{\dot{m}}(R_1,R_2\,|\,f)&=\int\limits_{R_{\rm in}}^{R_{\rm out}}\frac{\d R'}{(R')^2}\Delta R(R')\overline{G}_{\dot{M}}(R_1,R',f)\overline{G}^*_{\dot{M}}(R_2,R',f)& \nonumber \\
& \times S_{a}(R',f)\otimes_f
\left[\delta(f)+\frac{S_{\dot{m}}(R',f)}{\dot{M}^2_0}\right],&
\eeq
where $S_{a}(R,f)$ is the PDS of initial perturbations, $S_{\dot{m}}(R,f)$ is PDS of the mass accretion rate defined by (\ref{eq:S_mdot}), $\Delta R$ is local radial scale, on which the initial fluctuations are correlated, and `*' denotes complex conjugation. The second term in the rhs of (\ref{eq:C_mdot}) is negligible if the perturbations of the mass accretion rate are much smaller than the average mass accretion rate. It is obvious that $C_{\dot{m}}(R_2,R_1\,|\,f)=C^*_{\dot{m}}(R_1,R_2\,|\,f)$.
It is also obvious that the cross-spectrum gives the power spectrum when $R_1=R_2$: $C_{\dot{m}}(R,R\,|\,f)=S_{\dot{m}}(R,f)$.

In order to measure the coherence between two radii $R_1$ and $R_2$ it is convenient to use coherence function:
\be
{\rm Coh}(R_1,R_2\,|\,f) \equiv\frac{|C_{\dot{m}}(R_1,R_2\,|\,f)|^2}{S_{\dot{m}}(R_1,f)S_{\dot{m}}(R_2,f)},
\ee
which takes the value within $[0,1]$.

Because the high frequency variability is suppressed by the transfer functions, coherent high frequency variability is taking place only at radii close to the origin of the initial perturbations (see Fig.\,\ref{pic:23}a). At low frequencies the variability is coherent at a broader range of radii. It is interesting that at low frequencies the mass accretion rate variability at smaller radii lags variability at larger radii, however, at high frequency (see Fig.\,\ref{pic:23}b) the situation might be opposite: the variability at larger radii lags relative to the variability at smaller radii. This happens for two reasons: (i) the mass accretion rate fluctuations can propagate outwards and (ii) the viscous frequency drops with radial coordinate and at sufficiently high frequency, variability from outer parts of a disc is more suppressed then variability from inner parts of accretion disc, i.e. fluctuations from outside are damped out whereas those propagating outwards from smaller radii are not.

\begin{figure}
\centering
\includegraphics[width=8.5cm]{pic9a.pdf}
\includegraphics[width=8.5cm]{pic9b.pdf}
\caption{{(a) The coherence function ${\rm Coh}(15,R\,|\,f)$ for variability at radial coordinate $R_1=15$ and other radii in the accretion disc. (b) The corresponding phase angle of the cross-spectrum. At low Fourier frequencies the variability at outer radii leads variability at $R=15$ (red color) and variability at inner radii lags variability at $R=15$. At sufficiently high frequencies the situation is opposite: the variability at smaller radii lead variability at larger radii. Parameters: $n=1$, $R_{\rm in}=0$, $R_{\rm out}=100$, $\Delta R=1$, $f_{\rm br}=100 f_{\rm v}$.}}
\label{pic:23}
\end{figure}

\section{Fourier properties of the output flux}

In this section we discuss the fourier properties of the output flux mainly in context of BHs, for which the flux is coming entirely from the accretion disc.
We assume that the total flux at the surface of the accretion disc drops with the radial coordinate as $F_{\rm sur}\propto R^{-3}$.

\subsection{Power spectra for a given energy band}

The total luminosity available to be radiated at a given region of the accretion disc is proportional to the local mass accretion rate. If local variability of the mass accretion rate is small in comparison with the average mass accretion rate, the fluctuations of flux in some energy band are $\propto \dot{m}(R,t)$. If the energy dependence stays constant in time, the flux in a given energy band is
\be\label{eq:Flux}
F_{\rm h}(t)\simeq\int\limits_{R_{\rm in}}^{R_{\rm out}}\d R'\,\frac{h(R')}{R'^2}\dot{m}(R',t),
\ee
where $h(R)$ is a weighting function which corresponds to the energy band. Then the PDS of the flux variability in the energy band is
\be
S_{F_{\rm h}}(f)=\int \limits_{R_{\rm in}}^{R_{\rm out}}\d R_1\,\int \limits_{R_{\rm in}}^{R_{\rm out}}\d R_2\, \frac{h(R_1)}{R_1^2}\frac{h(R_2)}{R_2^2} C_{\dot{m}}(R_1,R_2\,|\,f),
\ee
where $S_{\dot{m}}(R,f)$ is PDS of local mass accretion rate given by equation (\ref{eq:S_mdot}).

For simplicity the weighting function $h(R)$ for the hard energy band is taken to be either a step function:
\begin {equation}
h(R) \! =\!
\left\{  \begin{array}{ll}
\strut\displaystyle \! \!
1, & \mbox{for}\ R<3R_{\rm in}, \\
\! \!  \frac{1}{3},                                             & \mbox{for}\ R\geq 3 R_{\rm in},
 \end{array} \right. \nonumber
\end{equation}
{which corresponds to abrupt hardening of the spectrum within the innermost region or exponent: $h(R)\propto R^{-1}$, which corresponds to gradual hardening of the spectrum}. The weighting function for the soft energy band is taken to be $s(R)=1$ all over the accretion disc.

\begin{figure}
\centering
\includegraphics[width=8.7cm]{pic10.pdf}
\caption{{The power spectrum of total flux variability. The solid lines represent the result based on accurate Linden-Bell Green functions, while the dashed lines show the result based on Green functions given by $\delta$-function in the time domain. The real diffusion process suppress significantly the variability at high frequencies. Parameters: $R_{\rm in}=1$, $R_{\rm out}=100$, rms$_{\rm ini}=0.01$. The initial mass accretion rate perturbations are described by zero-centured Lorentzian profile with $f_{\rm br}=f_{\rm v}$. The hardness function is given by the step function.}}
\label{pic:PowerFlux}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=8.7cm]{pic11.pdf}
\caption{{The power spectra of total flux variability given for different distribution of power of initial variability over the accretion disc. The initial mass accretion rate perturbations are described by a zero-centured Lorentzian profile with $f_{\rm br}=100\,f_{\rm v}$. Parameters: $R_{\rm in}=1$, $R_{\rm out}=100$, $\Delta R=1$. The hardness function is given by the step function.}}
\label{pic:PowerFlux2}
\end{figure}

The power spectrum of the photon flux variability imprints information about the inner and outer radii of the accretion disc and, like the power spectrum of the mass accretion rate, has at least two breaks at about the viscous frequencies at inner and outer radii (see Fig.\,\ref{pic:PowerFlux}). The shape of the power spectrum strongly depends on the distribution of power of initial mass accretion rate fluctuations over the radial coordinate (see Fig.\,\ref{pic:PowerFlux2}).

\begin{figure}
\centering
\includegraphics[width=9cm]{pic12.pdf}
\caption{{(a) Power spectra of variability in hard (solid black line) and soft energy bands (dashed red line), (b) the coherence function between variability at hard and soft energy bands and (c) the phase lags between variability in hard and soft energy bands. The initial mass accretion rate perturbations are described by zero-centured Lorentzian profile with $f_{\rm br}=100\,f_{\rm v}$. Parameters: $R_{\rm in}=1$, $R_{\rm out}=100$, $\Delta R=1$. The hardness function is taken to be $\propto R^{-1}$.}}
\label{pic:PowerFlux3}
\end{figure}

\subsection{Cross-spectrum between two energy bands: hard and soft lags}

One can consider two energy bands with different emissivities given by weighting functions $h(R)$ (hard band) and $s(R)$ (soft band), which result in different power spectra in hard and soft energy bands (see Fig.\,\ref{pic:PowerFlux3}a). Then the cross-spectrum is calculated as follows:
\be\label{eq:CrossEnergyBand}
C_{\rm h,s}(f)= \int \limits_{R_{\rm in}}^{R_{\rm out}}\d R_1 \int \limits_{R_{\rm in}}^{R_{\rm out}}\d R_2\,
\frac{h(R_1)}{R^2_1} \frac{s(R_2)}{R^2_2} C_{\dot{m}}(R_1,R_2\,|\,f),
\ee
where $C_{\dot{m}}(R_1,R_2\,|\,f)$ is the cross-spectrum for the mass accretion rate given by equation (\ref{eq:C_mdot}).
{The corresponding coherence function (see Fig.\,\ref{pic:PowerFlux3}b) is given by
\be
{\rm Coh}_{\rm h,s}(f)=\frac{\left| C_{\rm h,s}(f)  \right|^2}{S_{F_{\rm h}}(f) S_{F_{\rm s}}(f)},
\ee
}
while the phase lag $\Phi(f)$ between the two energy bands (see Fig.\,\ref{pic:PowerFlux3}c) is given by
\be
\tan\Phi (f)=\frac{{\rm Im}[C_{\rm h,s}(f)]}{{\rm Re}[C_{\rm h,s}(f)]},
\ee
where ${\rm Im}[C_{\rm h,s}(f)]$ and ${\rm Re}[C_{\rm h,s}(f)]$ are imaginary and real parts of the cross-spectrum at frequency $f$. The corresponding time lag is
$t_{\rm lag}(f)=\Phi (f)/(2\pi f)$.

It is natural that the variability in the hard energy band lags the variability in the soft energy band. In this case the time lags are defined by the typical time of propagation of fluctuations from outer radii, where soft flux is produced, to inner radii, where the most energy flux is produced in the hard energy band. However, negative phase lags (when the soft energy band variability lags hard energy variability) are also possible at high frequencies (see Fig.\,\ref{pic:PhaseLag}). The negative lags are caused by the fact that the variability of mass accretion rate at the inner radii can affect the variability at the outer radii. It becomes important at high frequencies where the Green functions in the frequency domain suppress variability from the outer radii. The frequency where the soft variability lags hard variability is comparable to the viscous frequency at the inner disc radius. {It is important to note that the appearance of the negative phase lags is strongly dependent on specific conditions in the accretion disc. Moreover, negative lags arise at very high Fourier frequencies, where the variability can be principally affected by sound waves instead of viscous diffusion.}

The predicted negative phase lags are comparable to the observed negative lags in stellar mass black hole systems \citep{2011MNRAS.414L..60U,2015ApJ...814...50D} and AGN \citep{2010MNRAS.401.2419Z,2013ApJ...777L..23W,2014MNRAS.439.1548A}. Therefore, viscous diffusion in accretion disc can in principle be a key process responsible for negative time lags, but we note that Fe K lags in AGNs are solid evidence of reverberation.

The coherence is typically smaller at higher Fourier frequency. It is a result of suppression of mass accretion rate variability due to viscous diffusion, when the variability at high frequency can be coherent only at close radii. However, because of the possibility of outward propagation (note that outward propagation effectively transfers variability of a particular characteristic frequency, suppressing variability at higher and lower frequencies), there is a range of Fourier frequencies, for which the coherence becomes higher under the additional influence of outward propagation (see the high frequency bump in Fig.\,\ref{pic:PowerFlux3}b). Depending on the displacement of regions producing hard and soft photon energy flux, the coherence becomes higher or even reaches its local maximum at sufficiently high frequency. It is remarkable, that the increasing of coherence at high frequencies is expected also from simple reverberation models \citep{2010MNRAS.401.2419Z,2014A&ARv..22...72U}.

\begin{figure}
\centering
\includegraphics[width=8.7cm]{pic13.pdf}
\caption{The phase lags as a function of Fourier frequency. Hard X-ray variability lags soft X-ray variability at low frequencies. However, high frequencies soft variability lags hard variability. It is caused by principal possibility of diffusion from inner to outer radii. Parameters: $R_{\rm in}=1$, $R_{\rm out}=100$. The initial perturbations are given by zero-centured Lorentzian profile. The hardness function is given by step function.}
\label{pic:PhaseLag}
\end{figure}

\subsection{Influence of non-linearity}

New local fluctuations of the surface density arise on top of the average mass accretion rate at each radial coordinate . Thus, new fluctuations can be affected by existing variability in the accretion flow, which gives rise to non-linear effects. The non-linear effects are described by the second terms in the rhs of equation\,(\ref{eq:S_mdot}) and equation\,(\ref{eq:C_mdot}).
The detailed study of non-linear effects is beyond the scope of the present work. However, a few features have to be marked: (i) the contribution of non-linear effects influence the shape of the PDS shape and Fourier frequency (see Fig.\,\ref{pic:NonLin1}), (ii) the higher the relative rms, the stronger the non-linear effects (see Fig.\,\ref{pic:NonLin2}).

\begin{figure}
\centering
\includegraphics[width=8.7cm]{pic14.pdf}
\caption{{The ratio of PDS with non-linear effects taken into account to PDS without non-linearity. Non-linear effects increase PDS.
 The initial mass accretion rate perturbations are described by zero-centured Lorentzian profile with $f_{\rm br}=100\,f_{\rm v}$. Parameters: $R_{\rm in}=1$, $R_{\rm out}=100$, $\Delta R=1$. The hardness function is taken to be $\propto R^{-1}$.}}
\label{pic:NonLin1}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=8.7cm]{pic15.pdf}
\caption{{The rms of total flux variability with effects of non-linearity taken into account as a function of rms, where non-linear effects are neglected (rms$_{0}$, dashed line represents rms=rms$_{0}$ corresponding to linear case). The stronger the variability, the stronger non-linear effects. The initial mass accretion rate perturbations are described by zero-centured Lorentzian profile with $f_{\rm br}=100\,f_{\rm v}$. Parameters: $R_{\rm in}=1$, $R_{\rm out}=100$, $\Delta R=1$. The hardness function is taken to be $\propto R^{-1}$.}}
\label{pic:NonLin2}
\end{figure}

\subsection{Accreting neutron stars}

In the case of highly magnetized accreting NSs the accretion disc is interrupted at the magnetospheric radius, where the accretion disc temperature is too low to produce X-ray flux (unless the mass accretion rate is extremely high, as it is for the case of pulsating ultra luminous X-ray sources powered by NSs). Most of the flux is produced at the NS surface, where the matter channeled by the magnetic field looses its kinetic energy (however, at extremely high mass accretion rate the flux can be affected by the matter at the magnetosphere, see e.g. \citealt{2017MNRAS.tmp..143M}). As a result, the observed fast aperiodic variability of the X-ray flux in the case of highly magnetized accreting NS is caused by mass accretion rate variability at the thin inner ring of the accretion disc. This feature causes observed variability of the PDS with the luminosity: the higher the mass accretion rate, the smaller the inner disc radius, the higher the Fourier frequency of a break corresponding to the inner disc radius \citep{2009A&A...507.1211R}. The PDS of observed X-ray variability, therefore, imprints information about the distant regions of accretion disc \citep{2005astro.ph..1215G}, which is useful for diagnostics of the whole system. It is interesting that accreting highly magnetized NSs provide an opportunity to study timing properties of cold accretion discs, where the temperature drops below $\sim 6500\,{\rm K}$ and hydrogen is weakly ionized and viscous properties of the disc can be significantly different \citep{2017arXiv170304528T}. The same condition of an accretion disc can be realized in BH XRBs, but only at much lower mass accretion rates.

\section{Summary}
\label{discussion}

We have developed a formalism that describes propagation of mass accretion rate fluctuations in an accretion disc using an exact solution of the diffusion equation \citep{1974MNRAS.168..603L}. The solution of the equation is given by Green functions, which can be constructed for various boundary conditions. The Green functions in the frequency domain play the role of transfer functions and describe how initial variability at a given radius results in variability at other radii. The fluctuations of the surface density arising at some radius due to local processes result in mass accretion fluctuations propagating both inwards and outwards. The developed formalism provides the basis for the construction of accurate models of fast aperiodic variability of mass accretion rate in accretion discs around BHs and NSs.

It is shown that the Green functions suppress variability at frequencies higher than the local viscous frequency due to diffusive spreading. The suppression depends on the Fourier frequency of variability and is stronger for lower parameter $n$, which describes dependence of kinematic viscosity on the radial coordinate $\nu\propto R^n$ (see Fig.\,\ref{pic:PS13}a and Fig.\,\ref{pic:PS18}a). As a result, the outer parts of the accretion disc effectively produce variability only below their viscous frequency even if high frequency variability is produced locally. High frequency variability can propagate only to the neighboring radii without significant suppression. The suppression at high Fourier frequencies encoded in Green functions is in agreement with recent numerical simulations based on the first principles \citep{2014ApJ...791..126C,2015arXiv151205350H}. This property of viscous diffusion means that different parts of the accretion disc effectively produce mass accretion rate fluctuations only at timescales longer than the local viscous time scale.

It is interesting that according to solutions of the diffusion equation, the mass accretion rate fluctuations propagate outwards as well, which has not been explored before in the context of propagating mass accretion rate fluctuations. In these conditions, the outer parts of accretion disc ``feel" processes which occur in the inner part of the disc. However, in the case of outward propagation of fluctuations, the transfer functions (see Fig.\,\ref{pic:PS13}b and Fig.\,\ref{pic:PS18}b) differ from the transfer function describing inward propagation. The high frequency variability is strongly suppressed similarly to the case of inward propagation, but the transfer functions have a maximum at a frequency close to local viscous frequency. It means that the low frequency variability is also suppressed for fluctuations propagating outwards. In this case the suppression is stronger for lower frequencies and higher parameter $n$.

The power spectrum of the mass accretion rate variability at any radial coordinate is affected by properties of initial fluctuations all over the disc (distribution of variability power over the radial coordinate and break frequency $f_{\rm br}(R)$ of initial variability) and by viscous properties of the accretion disc, which are described by Green functions. The observed variability of the photon energy flux in a given energy band is defined by mass accretion rate variability at the inner parts of accretion disc. Thus, it is in principal possible to decode the properties of initial variability using observed photon energy flux variability and appropriate Green functions.

Our model assumes that there is an effective radial scale $\Delta R$ in the accretion disc where the initial fluctuations are correlated. This assumption is in agreement with recent numerical simulations \citep{2015arXiv151205350H}, where the radial scale is associated with the typical scale at which the dynamo process develops \citep{2004MNRAS.348..111K}. The influence of initial fluctuations on the final variability depends on the local radial scale $\Delta R$. {Previous, discretized, treatments of the propagating fluctuations model (e.g. \citealt{2013MNRAS.434.1476I}) also implicitly make this assumption.}

The developed model naturally predicts two breaks in the PDS: the lower frequency break corresponds to the outer disc radius and settles at about the viscous frequency of the outer radius, the higher frequency break corresponds to the inner disc radius, its location depends on viscous frequency at the inner radius and PDS of initial perturbation there. Between the breaks the PDS can be described by power laws except the highest frequencies, where the PDS can drop faster than it is described by power law.

In previous, simpler treatments of propagating fluctuations, the predicted power spectrum has always consisted of flat top noise between a low and high frequency  break if the total variability power is assumed to vary smoothly with radius. Observed power spectra, on the other hand, routinely show a series of discrete humps, that have been modeled by injecting extra variability at characteristic radii \citep{2016MNRAS.462.4078R,2017arXiv170407705R,2017MNRAS.469.2011R}. Here, however, we show our more sophisticated treatment of the problem naturally leads to discrete humps in the predicted power spectrum, even for input variability varying smoothly with radius.
{\cite{2017MNRAS.469.2011R} recently showed that observations of XTE J1550-564 pose problems for the propagating fluctuations model, since the power spectra in two energy bands are almost identical, but there are still phase lags between the two bands. This is very challenging to explain, because the lags are assumed to result from the different bands having different radial emissivity laws, which inevitably leads to different power spectral shapes. This is still the case in our more sophisticated treatment of the model.
{The unique properties of XTE J1550-564 were reported earlier by \citealt{2016A&A...591A..77S} who reported on the evidence for non-linear dynamics in the lightcurves of the source.}
However, the more nuanced predictions explored here leave more scope for the XTE J1550-564 data to be explained. In particular, if the property seen in the observations of \cite{2017MNRAS.469.2011R} is fairly rare, it is reasonable to fit it with rather ``fine tuned" parameters. We also note that the same observation is a problem for all lag models in the literature (e.g. \citealt{2015A&A...584A.109R}). }

The model in the linear regime, i.e., when the initial fluctuations are not affected by mass accretion rate variability, reproduces the observed linear relation between rms and luminosity \citep{2005MNRAS.359..345U}. However, it does not exclude the possibility of deviation from the linear relation. A deviation can be caused by nonlinear interaction between the fluctuations of mass accretion rate.
We predict deviation from the linear rms-flux relation for very high levels of rms variability. Discovery of such a property in observational data would allow us to break degeneracies associated with propagating fluctuation models.

The variability at smaller radii usually lags variability at larger radii, but we show that at high frequency the situation might be opposite and variability at larger radii lags variability from the inner parts of accretion disc. It leads to the possibility of negative time-lags in BH XRBs: the variability in soft X-rays lags variability in hard X-rays at sufficiently high Fourier frequencies. The typical value of the corresponding negative phase lags is $\sim 0.05\,{\rm rad}$, which is close to the lags observed in AGNs \citep{2010MNRAS.401.2419Z,2013ApJ...764L...9C,2014MNRAS.439.1548A} and stellar-mass black holes \citep{2011MNRAS.414L..60U,2015ApJ...814...50D}. The negative lags are commonly explained by Comptonization time delays \citep{1997ApJ...480..735K} or by Compton-scattering reverberation from material located at some distance from the primary X-ray source \citep{2009Natur.459..540F,2012ApJ...760...73L}.
Thus, we conclude that viscous diffusion can potentially contribute significantly to the observed lags of soft X-ray radiation.

\section*{Acknowledgments}

This research was supported by the Netherlands Organization for Scientific Research (AAM) and Veni grant (AI). Partial support comes from the EU COST Action MP1304 ``NewCompStar". The authors are grateful to Galina Lipunova and Andrey Semena for useful discussions.

{


text:
\section{Introduction}
\label{intro}

Turbulence is associated with large fluctuations of velocity gradients which
appear preferentially at the smallest scales of the flow. The amplitudes of
the fluctuations exceed the mean values by orders of magnitude when the Reynolds
number of the flow is sufficiently large. This behaviour is known as
small-scale intermittency \cite{Sreenivasan1999}. It has been discussed recently that
a deeper understanding of fluid turbulence as a whole requires a detailed resolution
of the intermittent dynamics at the small-scale end of the inertial range \cite{Sreenivasan2004}.
This will include scales that are smaller than the (mean) Kolmogorov dissipation length
$\eta_K=\nu^{3/4}/\langle\epsilon\rangle^{1/4}$ with the kinematic viscosity $\nu$ and the mean
energy dissipation rate $\langle\epsilon\rangle$. Although significant progress in
measurement techniques has been made, the finest structures remain still spatially
unresolved \cite{Zeff2003} or the flows are limited to very low Reynolds numbers
\cite{Mullin2006}. Direct numerical simulations can reach today sufficiently high Reynolds
numbers while resolving the gradients at the small-scale end of the inertial range
cascade properly \cite{Schumacher2007,Donzis2008,Horiuti2008}.

The small-scale structure and statistics of turbulence has been mostly studied for the case
of homogeneous, isotropic and statistically stationary turbulence. A wall-bounded shear flow consists
of a boundary layer which is dominated by coherent streamwise structures and a central bulk region
in which they are basically absent. The statistics in the wall-normal direction is inhomogeneous and
requires a height-dependent analysis on account of the varying strength of the mean shear.
As a consequence, there is significantly less work on this aspect for wall-bounded flows than for homogeneous isotropic
turbulence
\cite{Kim1993,Saddoughi1994,Blackburn1996,Jimenez2006}.
The coupling between large-scale shear and intermittency in the inertial cascade range has been
studied in channel flow simulations \cite{Toschi1999,Toschi2000}. Further studies in the same
flow configuration demonstrated the violation Kolmogorov's refined similarity hypothesis (K62) \cite{Kolmogorov1962}
when approaching the walls \cite{Benzi1999,Casciola2001}. It was therefore suggested to replace the
third-order increment moment in  K62 by a second-order one which is directly coupled to the large-scale
shear in the corresponding balance equations. As a consequence, shear-improved eddy viscosities
in Smagorinsky subgrid-scale models can improve the representation of Reynolds stresses as shown in Ref.
\cite{Leveque2007}.

In comparison with homogeneous, isotropic and statistically stationary turbulence,
the presence of an inhomogeneous direction increases
the level of flow complexity significantly. With regard to the
velocity gradient statistics one can identify the following two aspects.
First, it is more challenging to measure all nine derivatives of the velocity
gradient tensor with the sufficient resolution in an open shear-flow-setup (see e.g. \cite{Adrian2007}).
Second,  it is frequently believed that shear flow turbulence is in a state of local isotropy at
the small-scale end for larger Reynolds numbers \cite{Saddoughi1994}. More recent
experimental and numerical studies demonstrated however that significant deviations persist,
in particular when higher-order moments are discussed \cite{Shen2000,Schumacher2003,Biferale2005}.
It can thus be expected that some fundamental differences to the homogeneous isotropic case persist
to very high Reynolds numbers. All this suggests to our view a systematic study of the height-dependence
of the statistics of the velocity gradient fields in turbulent shear flows.

In the present work, we want to make a first step in this direction and conduct an analysis
of the small-scale statistics of the velocity gradient in a wall-bounded shear flow. It is based
on three-dimensional direct numerical simulations of a turbulent channel flow. Our study will be focussed to two
fields, the local enstrophy field and the energy dissipation rate field. The gradient tensor of
the velocity fluctuations, $m_{ij}$, can be decomposed into a symmetric and an antisymmetric part
\begin{equation}
m_{ij}=\frac{\partial v_i^{\prime}}{\partial x_j}=a_{ij}+s_{ij}\,.
\end{equation}
The local enstrophy field,
${\bm \omega}^2=({\bm\nabla\times v^{\prime}})^2$, probes the magnitude of the antisymmetric part
$a_{ij}=-\epsilon_{ijk}\omega_k/2$ of the velocity fluctuation gradient tensor.
Here $\epsilon_{ijk}$ is the fully antisymmetric Levi-Civita tensor. The energy dissipation rate field,
$\epsilon=2\nu s_{ij}s_{ji}$, measures basically the magnitude of the symmetric part $s_{ij}$ of
$m_{ij}$. The ensemble average of both fields are connected by
\begin{equation}
\langle\epsilon\rangle=\nu \langle {\bm \omega}^2\rangle
\label{epsom2}
\end{equation}
for homogeneous isotropic turbulence. This does however not imply that the local statistical fluctuations
are synchronized \cite{Zeff2003,Chen1997}. In addition to homogeneous isotropic turbulence, we have
to conduct the statistical analysis for both fields as a function of distance from the wall and to separate
the turbulent fluctuations from the mean flow.

The outline of the paper is as follows. We discuss the equations of motion and the numerical
model in the next section. Afterwards we investigate the dependence of the small-scale statistics
on the spectral resolution which sets the stage for the analysis. The subsequent section
investigates the single-field and joint statistics of the local enstrophy and dissipation rate
fields as a function of the distance from the wall. We summarize our findings and give a brief
outlook at the end.

\begin{figure}
\begin{center}
\includegraphics[angle=0,scale=0.5,draft=false]{fig0_sketch_channel}
\caption{Sketch of the channel flow geometry and the coordinates. The simulation
domain is $L_x \times L_y\times L_z= 4\pi \times 2\pi \times 2$ in units of the
half-width $L$.}
\label{fig1}
\end{center}
\end{figure}
\section{Equations and numerical model}
The numerical method solves  the Navier-Stokes equations
for an incompressible flow which are given in dimensionless form by
\begin{equation}
\label{nsviscous}
\frac{\partial \bm{v}}{\partial t} -
\bm{v}\times\bm{\Omega}
= -\nabla \left(p+v^2/2\right)
+ \frac{1}{{\rm Re}} \nabla^2 \bm{v},\hskip4mm
\nabla \cdot\bm{v}  =  0,
\end{equation}
with $\bm{\Omega}=\nabla\times\bm{v}$ and boundary conditions
\begin{equation}
\label{bctopviscous}
\bm{v}=0
\hskip3mm\mbox{at $z=\pm1$}.
\end{equation}
Figure \ref{fig1} illustrates the geometry and the coordinates of the channel flow.
The applied pseudospectral method is based on Fourier series in the horizontal
directions $x$ and $y$ and a Chebyshev polynomial expansion in the
vertical $z$--direction \cite{Canuto1977,Gottlieb1977}. The solenoidal velocity field is
represented in terms of two scalar quantities  $\psi$ and $\varphi$ using the
poloidal-toroidal decomposition
\begin{equation}
\bm{v}({\bm x},t)=
\bm{\nabla} \times \left(\bm{\nabla} \times \bm{e}_z\varphi({\bm x},t) \right) +
   \bm{\nabla} \times \bm{e}_z\psi({\bm x},t).
\end{equation}
Equations for $\varphi$ and $\psi$ are derived by taking the curl and
twice the curl of the momentum equation and projection onto the vertical
direction. We obtain two equations for the vertical velocity component
$v_z=-\Delta_h \varphi$ and the vertical vorticity component
$\Omega_z= -\Delta_h \psi$, where $\Delta_h=\partial_x^2+\partial_y^2$
denotes the horizontal Laplace operator. The quantities $v_z$ and $\Omega_z$
determine the velocity field up to a mean flow $U(z,t) \bm{e}_x +V(z,t) \bm{e}_y$.
Equations for $U$ and $V$ are obtained by averaging the momentum
equation over horizontal cross-sections of the periodicity domain.
Upon introducing the definitions
\begin{equation}
\zeta=\Omega_z,\hskip4mm \eta=\nabla^2 v_z,\hskip4mm \xi=v_z
\end{equation}
for ease of notation, the evolution equations based on the poloidal-toroidal
representation take the form
\begin{eqnarray}
&&\nabla^2 \zeta - {\rm Re}
\, \partial_t \zeta   = -{\rm Re} \,\bm{e}_z \cdot\left[  \bm{\nabla}
\times \left(\bm{v} \times  \bm{\Omega}\right)\right],\\
\label{eta-eq}
&& \nabla^2 \eta - {\rm Re} \, \partial_t \eta= {\rm Re}\left[\partial_z  \bm{\nabla} \cdot
\left(\bm{v} \times  \bm{\Omega}\right)-
\bm{e}_z \cdot  \nabla^2  \left(\bm{v} \times  \bm{\Omega}\right)\right],\\
\label{xi-eq}
&& \nabla^2 \xi = \eta,\\
\label{meanflow-x-eq}
&& \partial_z^2 U-{\rm Re}\, \partial_t U =
{\rm Re}\,\left[ \partial_z \left\langle v_x v_z\right\rangle_A + \partial_x p \right], \\
&& \partial_z^2 V-{\rm Re}\, \partial_t V=
{\rm Re}\,\left[ \partial_z \left\langle v_y v_z\right\rangle_A  \right].
\end{eqnarray}
The angular brackets $\langle\cdot\rangle_A$ denote horizontal averages. The
mean flow equations also contain the  mean driving pressure gradient  $\partial_x p$.
The appropriate boundary conditions are readily derived using the incompressibility
constraint $\nabla \cdot\bm{v}=0$. We obtain
\begin{equation}
\xi=\zeta=\partial_z\xi=U= V=0 \hskip1cm \mbox{at} \hskip4mm z=\pm 1\,.
\end{equation}
For the time discretization of the evolution equations we use a
method of second-order accuracy. If we write such an equation
symbolically as $\partial_t f = {\cal L} f + {\cal N}(f)$ where $\cal L$ denotes
a linear operator and $\cal N$ the remaining terms, our time-stepping scheme reads
\begin{equation}
\frac{3 f^{n+1}-4 f^n + f^{n-1}}{2 \Delta t}= {\cal L} f^{n+1}
+ 2{\cal N}(f^n)-
{\cal N}(f^{n-1})\,,
\end{equation}
where $\Delta t$ is the time step. The left hand side approximates the time
derivative $\partial_t f$ at the time level $n+1$ using the previous two time
levels. The linear term ${\cal L}f$ is treated implicitely, and $\cal N$
explicitly through the second-order Adams-Bashforth method, where the prefactors
correspond to a linear extrapolation to the time level $n+1$ \cite{Krasnov2004}.
The solution procedure for the coupled eqns. (\ref{eta-eq}) and (\ref{xi-eq}) requires auxiliary
functions to satisfy the boundary conditions $\partial_z \xi=0$ for eq. (\ref{eta-eq})
as explained in \cite{Krasnov2004}. The numerical code implements distributed-memory
parallelization through MPI.

\begin{figure}
\begin{center}
\includegraphics[angle=0,scale=0.5,draft=false]{fig1_tau13}
\includegraphics[angle=0,scale=0.5,draft=false]{fig2_diss}
\caption{Vertical profiles of the Reynolds shear stress $\tau_{13}$ (left)
and the energy dissipation rate $\epsilon$ (right) averaged in planes of constant $z$.
The wall-normal coordinate $z$ is given in wall units. The results of
small- ($128^3$ points, blue curves), medium- ($256^3$ points, green curves) and
high-resolution ($512^3$ points, red curves) simulations are shown. In each case,
the quantities are made nondimensional with
the wall friction velocity $u_\tau$ from the simulation with $512^3$ points. }
\label{fig2}
\end{center}
\end{figure}
\section{Single quantity statistics}
\subsection{Spectral resolution effects}
Our statistical analysis is focussed to the small-scale properties of the
channel flow turbulence. In order to set the stage for this investigation,
we first compare  statistical properties at three different spectral resolutions and the
same set of parameters. The first case with $128^3$ (or more precisely $128^8\times 129$) grid
points corresponds with the classical simulation\footnote{More precisely, the grid resolution in
\cite{Kim1987} was $N_x\times N_y\times N_z=190\times 129\times 160$ (with $y$ being the
wall-normal coordinate) for a channel with the same sidelengths as in the present work.}
of Kim, Moin and Moser \cite{Kim1987} in terms of the vertical resolution. The second case is
two times better resolved in each space direction with $256^3$ grid points. The highest spectral
resolution is  four times better resolved than the first  one and contains $512^3$ (or more precisely
$512^2\times 513$) collocation points. The three runs  were performed with identical fixed volume flux
such that the resulting friction Reynolds number approaches the grid-independent
limit  $Re_{\tau}=180$.  The definition
\begin{equation}
Re_{\tau}=\frac{u_{\tau}L}{\nu}\,.
\end{equation}
is based on the wall friction velocity $u_{\tau}$, the half-channel height
$L$, and the kinematic viscosity by $\nu$. We follow the usual convention and
denote the turbulent velocity fluctuations by
\begin{equation}
{\bm v}^{\prime}({\bm x},t)={\bm v}({\bm x},t)-U(z){\bm e}_x\,.
\end{equation}
Figure \ref{fig2} (left) compares the Reynolds shear stress $\tau_{13}=\langle
v_x^{\prime} v_z^{\prime}\rangle_{A,t}$ as a function of $z^+=z \nu/u_{\tau}$.
The three vertical profiles of $\tau_{13}$ coincide quite well,
except for slight deviations around the maximum in case of the lowest resolution.
The additional subscript $t$ in $\langle \cdot\rangle_{A,t}$ indicates time averaging.
We repeat the analysis for the vertical profile of the
energy dissipation rate averaged
over time and in planes at constant distance from the wall. The quantity is defined as
\begin{equation}
\epsilon_{avr}(z)=\langle\epsilon({\bm x},t)\rangle_{A,t}=2\nu\langle s_{ij}s_{ji}\rangle_{A,t}=\frac{\nu}{2}
\left\langle\left(\frac{\partial v_i^{\prime}}{\partial x_j}+\frac{\partial v_j^{\prime}}{\partial x_i}
\right)^2\right\rangle_{A,t}\,.
\end{equation}
The plot yields almost perfectly collapsing curves for most of the channel. Deviations for the
smallest resolution  exist for the viscous layer at $z^+\le 20$. We can conclude
that the lowest resolution is sufficient when first-order moments are considered,
such as shear stresses. It is therefore suitable for many purposes, such as validations
of closures of the Reynolds averaged equations \cite{Wilcox}. We can however see already,
that the mean dissipation rate profile deviates significantly in the viscous layer.

A full statistical information requires the knowledge of all moments or the
probability density function (PDF).  We proceed therefore by comparing the PDFs of the
energy dissipation rate, $\epsilon({\bm x},t)$, and the local enstrophy,
${\bm \omega}^2({\bm x},t)=({\bm\nabla\times v^{\prime}})^2$, for both resolutions.
\begin{figure}
\begin{center}
\includegraphics[angle=0,scale=0.53,draft=false]{fig3_pdf_128e}
\includegraphics[angle=0,scale=0.53,draft=false]{fig3_pdf_512e}
\includegraphics[angle=0,scale=0.53,draft=false]{fig3_pdf_128o}
\includegraphics[angle=0,scale=0.53,draft=false]{fig3_pdf_512o}
\caption{Comparison of the probability density functions of the energy dissipation
rate (upper row) and local enstrophy (lower row) for the two spectral resolutions.
Data are obtained in planes at four different distances from the wall which are indicated
in the legends. The left column stands for the lowest resolution case, the right one for the
highest resolution. Both quantities are given in units of their statistical mean in the
corresponding plane. We denote $\epsilon_{avr}=\langle\epsilon\rangle_{A,t}$ and
$\omega^2_{avr}=\langle\omega^2\rangle_{A,t}$.}
\label{fig3}
\end{center}
\end{figure}
Four different distances have been chosen: $z^+=10$ is in the viscous buffer layer,
$z^+=30$ at the begining of the loga\-rithmic layer where the Reynolds shear stress
has its maximum (see Fig. \ref{fig2} (left)), $z^+=72$ in the loga\-rithmic region,
and $z^+=180$ in the midplane of the channel. While the first-order moments are
basically unchanged at the higher resolution, we can detect clear differences
in the far tails of the PDFs for all cases. The far tails are exactly the parts of the PDF
that become relevant for higher-order moments. Their sufficient resolution requires
the efforts made in the present work (see also \cite{Schumacher2007} for a detailed discussion).

\begin{table}[h]
\caption{Moments of the normalized energy dissipation rate as a function of height and resolution. }

\begin{center}
\begin{tabular}{@{}|c|ccc|ccc|ccc|@{}}
\hline
$~~~$ & \multicolumn{3}{|c|}{$<\epsilon^2>$} & \multicolumn{3}{|c|}{$<\epsilon^3>$} & \multicolumn{3}{|c|}{$<\epsilon^4> / <\epsilon^2>^2$}\\
\hline
$z^+$ & $128$ & $256$ & $512$ & $128$ & $256$ & $512$ & $128$ & $256$ & $512$\\
\hline
$ 10$ & $2.12$ & $2.15$ & $2.14$ & $ 7.73$ & $ 7.93$ & $ 7.83$ & $ 9.59$ & $ 9.45$ & $ 9.34$\\
$ 30$ & $2.25$ & $2.93$ & $3.14$ & $ 9.92$ & $17.95$ & $21.83$ & $15.45$ & $20.40$ & $25.44$\\
$ 72$ & $2.46$ & $2.96$ & $3.15$ & $13.00$ & $20.43$ & $24.36$ & $21.34$ & $32.53$ & $36.77$\\
$180$ & $2.55$ & $2.81$ & $3.06$ & $14.56$ & $18.67$ & $23.07$ & $25.24$ & $30.47$ & $35.94$\\
\hline
\end{tabular}

\end{center}
\label{tab:epsilon}
\end{table}

\subsection{Height dependence of statistics}
Let us turn now to the high-resolution case of both fields in the right column
of Fig. \ref{fig3}. Note first that the amplitudes are always rescaled by the corresponding mean
values for the given plane at fixed $z^+$. Close to the wall, the turbulent velocity fluctuations are small.
The mean dissipation profile indicates however that this is not the case for the magnitude of the
derivatives. As we can see in Fig. \ref{fig2}, the mean energy dissipation rate is largest
close to the wall and drops by more than two orders of magnitude toward the channel centre plane.
The tails of both distributions are the sparsest ones in the vicinity of the wall. While the means are large
the variations about the means remain small. The extension of the tails varies with increasing
distance from the wall as the data for the three other values of $z^+$ indicate. Since the mean values
decrease monotonically towards the midplane, high-amplitude events are shifted ever further into the tail.
The fattest tails are observed at $z^+=72$, i.e. in the middle of the logarithmic region.
Particularly here, high amplitudes can be expected further away from the wall
due to bursting events. They are for example caused by coherent hairpin vortices that can grow
out into the logarithmic layer. Such a large-scale hairpin vortex packet formation has been recently
suggested by Adrian and co-workers on the basis of experiments and direct numerical simulations
\cite{Adrian2007}. Our statistical analysis is consistent with this observation.

The high resolution case reproduces another feature which is well-known from the homogeneous
isotropic case. The tail of the local enstrophy is fatter than that of the energy
dissipation rate \cite{Zeff2003,Donzis2008}. Simulations of isotropic box turbulence by Donzis
{\it et al.} \cite{Donzis2008} suggest that this difference in the tails decreases with increasing
Reynolds number, a point that certainly has to be studied in the shear flow case as well.

Height and resolution dependences of higher moments of the dissipation rate
are given in table \ref{tab:epsilon}. The fatter tails of the PDFs at higher resolution account
for the growth in the second and third moment with the resolution
for all $z^+$ values except in the viscous region at $z^+=10$. Likewise, the flatness
$\langle \epsilon^4\rangle/\langle \epsilon^2\rangle^2$ is growing with the
resolution. This quantity is also indicative of the degree of
intermittency at the smallest lengthscales.
Remarkably, it is not constant, but increases significantly from
$z^+=10$ to $z^+=72$, i.e. the level of small-scale intermittency depends
on the wall distance outside the viscous region.
 However, it changes only slightly from $z^+=72$ to the midplane.

Regarding the relation between planar averages of enstrophy and energy
dissipation we remark that they are
connected by the relation
\begin{equation}
\langle \epsilon\rangle_A=\nu \langle{\bm\omega}^2\rangle_A +\nu \partial_z^2\langle
v_z^{\prime\, 2}\rangle_A\,
\label{epsom3}
\end{equation}
in a wall-bounded shear flow.
In contrast to Eq. (\ref{epsom2}), this relation contains  an additional derivative term
arising from the averaging which is restricted to planes of fixed $z$. We found that the second term on the right
hand side of (\ref{epsom3}) is at least two orders of magnitude smaller than $\langle{\bm\omega}^2\rangle_A$.
We have also verified that relation (\ref{epsom3}) is satisfied in our numerical simulations.
Furthermore, the statistics of the energy dissipation rate field deviated from the log-normal distribution
for all four planes, particularly in the far tails at the small and large amplitudes. This property is similar to
isotropic turbulence.

\section{Joint statistics of energy dissipation and local enstrophy}
In view of the differences in the PDFs for enstrophy and energy dissipation it is
interesting  to study how both fields are statistically correlated.
For this reason we plot in Fig. \ref{fig4} the joint PDF
of $\epsilon$ and ${\bm \omega}^2$ normalized by both single quantity PDFs
\begin{equation}
\Pi(\epsilon, {\bm \omega}^2)=\frac{p(\epsilon, {\bm \omega}^2)}{p(\epsilon)\,p({\bm \omega}^2)}
\label{jointpdf}
\end{equation}
\begin{figure}
\begin{center}
\includegraphics[angle=0,scale=0.38,draft=false,clip=]{fig4_pdf_512_0256}
\includegraphics[angle=0,scale=0.38,draft=false,clip=]{fig4_pdf_512_0361}
\includegraphics[angle=0,scale=0.38,draft=false,clip=]{fig4_pdf_512_0417}
\includegraphics[angle=0,scale=0.38,draft=false,clip=]{fig4_pdf_512_0457}
\caption{Joint and normalized PDFs for the energy dissipation rate
$\epsilon$ and the square of vorticity $\omega^2$ as given by Eq. (\ref{jointpdf}).
Data are taken again at four $z^+$-positions. Upper left: $z^+ = 180$. Upper right:
$z^+ = 72$. Lower left: $z^+ = 30$. Lower right: $z^+ = 10$. Results of DNS with
$512^3$ grid points are shown, the joint PDF fields are normalized by the product of
separate PDFs for $\epsilon / \epsilon_{avr}$ and $\omega^2 / \omega^2_{avr}$. The
white diagonal line in the lower right figure corresponds with $\epsilon=\nu{\bm\omega}^2$.
The colour coding is the same in all four figures and given in units of decadic logarithm.
It varies from $10^0$ (red) to $10^{-4}$ (blue).}
\label{fig4}
\end{center}
\end{figure}
High amplitudes of $\Pi$ indicate a strong statistical correlation between both fields.
The PDF of the data in the midplane has the broadest support. Very high-amplitude events of
dissipation and vorticity are strongly correlated as indicated by the local maximum which is
stretched out into the upper right corner of the plane. This maximum is even more pronounced in the
logarithmic layer as can be seen in the upper right panel of Fig. \ref{fig4}. The local
maximum of $\Pi$ in the outer right corner becomes narrower with decreasing distance from the wall.
It is in line with an overall decrease of the support for highly correlated events. The observed
trend is consistent with our observations for the single quantity statistics. Again, we can detect
a decrease of the fluctations of both fields about the mean. As a guide to the eye, we added the line
$\epsilon=\nu {\bm\omega}^2$ (in white) to the lower left panel, i.e. to the joint-PDF plot closest
to the wall. It can be seen that  a broad range of amplitudes of both fields is basically concentrated now
around this line, which connects the means of both fields (see eqns. (\ref{epsom2}) or
(\ref{epsom3})). To conclude, at all
distances from the wall there is a strong statistical correlation between the local enstrophy and the
dissipation. Close to the wall high-amplitude events of both fields seem to become synchronized.
The support of the statistically correlated events varies significantly across the channel.

\begin{figure}
\begin{center}
\includegraphics[angle=0,scale=0.65,draft=false]{fig5_shot_eps}
\vspace{0.1cm}
\includegraphics[angle=0,scale=0.65,draft=false]{fig5_shot_om2}
\caption{Isosurfaces of the energy dissipation rate (upper figure) and the local enstrophy
(lower figure). The isolevels are taken at three times the mean dissipation and the mean
squared vorticity. }
\label{fig5}
\end{center}
\end{figure}

Figure \ref{fig5} displays the structures of both fields. The largest amplitudes
of both fields are clearly concentrated close to the walls.  We also observe the hairpin-like
vortices that burst into the channel center (see lower panel). Isosurfaces of  the energy
dissipation rate display more sheet-like structures (see upper panel). They reach also out into the logarithmic
region. Such shear layers can be formed when pairs of streamwise streaks are lifted up
\cite{Schoppa2000} or when they become linearly unstable. The inspection of the structures
in both plots supports thus our statistical results obtained before. Maxima of both fields are
found close together in space.  In particular, we expect, that exactly the bursting events will
contribute significantly to the local maximum of the joint PDF in the upper right panel of Fig. \ref{fig4}.
A firm proof would require a separation of these coherent structures in the data record, which
will be part of our future work on this subject.

\section{Summary and conclusions}
We have presented a study of the statistics of the energy dissipation rate and
local enstrophy in a turbulent channel flow. We showed that a full statistical
analysis of both fields requires a spectral resolution which is larger than the standard
one. This was one objective of our study, namely to make this case in a turbulent wall-bounded
shear flow. It clearly limits the range of accessible Reynolds numbers. Based on this result,
we investigated the height-dependence of the statistics of both fields in a first example.
We showed that the variations of both fields about their mean values are strongest in the
logarithmic layer. They decrease slightly towards the channel centre and significantly towards
the wall. The
statistics of the velocity gradients
is thus strongly height-dependent. It is clear that a more comprehensive parametric study has to
follow, e.g., in order to analyze these trends at different Reynolds numbers.

There is a further interesting aspect onto which our studies can shed new light.
Recently, new ideas for multifractal subgrid-scale models have been developed which are based on exactly
such statistical fluctuations of velocity derivatives which have been studied
here \cite{Burton2005}. The model of \cite{Burton2005} has been successfully validated in homogeneous
isotropic turbulence and it seems to be a hopeful candidate for  application to shear flow turbulence.
Investigations in this direction are under way.

{\em Acknowledgements.} Some parts of this work grew out of discussions with
K. R. Sreenivasan and V. Yakhot. We wish to thank both for their suggestions.
Our work is supported by the Emmy-Noether-Program (DK and TB) and the
Heisenberg-Program (JS) of the Deutsche Forschungsgemeinschaft (DFG). The
supercomputing ressources were provided by the J\"ulich Supercomputing Centre (Germany) under grant HIL01.

redpajama_set_name: RedPajamaArXiv